# Inversion of the horizontal-to-vertical spectral ratio in presence of strong lateral heterogeneity

Inversion of the horizontal-to-vertical spectral ratio in presence of strong lateral heterogeneity Abstract The horizontal-to-vertical spectral ratio (HVSR) method has recently gained a resurgence of interest because of its new interpretation through the diffuse field approach. According to this theory, the HVSR measurement is linked to the Green’s function retrieval. While the HVSR method is traditionally used to evaluate simple and shallow velocity model, we demonstrate that complex Earth’ structure (that includes composite geological layering) can be assessed by choosing appropriate data processing and inversion scheme. In particular, we explore the effects of a strong lateral heterogeneity on the HVSR by using measurements along a vertical crater wall. However, wave propagation computation with lateral heterogeneity is inefficient to conduct an inversion. Then, to accelerate the forward calculation, we identified a simplified model consisting of an unbounded 2-D multilayer medium that allows describing the observed measurements with the use of a correction factor. The 2-D character comes from a non-standard noise illumination related to the wind blowing perpendicularly to the crater walls. It is supported by several theoretical and experimental arguments. We also show that the H and V components can be retrieved independently, which may improve the sensitivity of the inversion at high frequency. To reduce the non-uniqueness of the problem, we consider the joint inversion of the H and HVSR’s at several positions where the thicknesses of the layers were evaluated thanks to a large geological outcrop. The good agreement between the recovered velocity structure and the geology shows that the method is not limited to unbounded simple shallow structures, which opens the domains of application of the method to exploration geophysics or to hazard assessment. Joint inversion, Waveform inversion, Seismic noise, Site effects, Volcano seismology 1 INTRODUCTION The horizontal-to-vertical spectral ratio (the so-called HVSR or H/V) is a widely used method to determine the natural frequency of a site and to estimate a simple velocity model of the subsurface; that is, typically one layer over a half-space (e.g. Nakamura 1989; Lermo & Chávez-García 1994; Bard 1999). Such a simple interpretation comes primarily from the fact that the technique was developed empirically without clear theoretical basis. Actually, the fundamental peak of the HVSR curve was interpreted based on the assumption of predominance of body or surface waves in the ambient seismic field [see review of Lunedei & Malischewsky (2015); Piña-Flores et al. (2016)]. However, it is now well known that ambient seismic field is composed of both body and surface waves (e.g. Bonnefoy-Claudet et al. 2008) and therefore, an HVSR theory must account for it to explain the entire HVSR curve (i.e. not only the main peak). In this contribution, we base our argumentation on a recent theoretical developments that account for all the seismic phases, being then more robust to assess with details a velocity structure (Sánchez-Sesma et al. 2011a; Lontsi et al. 2015; García-Jerez et al. 2016; Piña-Flores et al. 2016; Spica et al.2017). So far, two different theories have been proposed for the modeling of the whole HVSR curve: the distributed surface sources (e.g. Lunedei & Albarello 2010; Albarello & Lunedei 2011; Lunedei & Malischewsky 2015) and the diffuse field approach (DFA, Sánchez-Sesma et al. 2011a,b; Piña-Flores et al. 2016). Even if these two branches of investigations have apparent solid divergence regarding their theoretical assumptions, other studies have already demonstrated their compliance through the equivalence of the reciprocity theorem and of the diffuse field theory (Weaver & Lobkis 2004; Perton et al. 2016). The DFA underlines the connection between the directional energy densities (DED), and the Green’s function (GF, Perton et al. 2009). This connection allows conducting the inversion of the HVSR measurement for physical properties evaluation. In this view, the data must be processed consistently with the theory in order to enhance the diffuse character of the ambient seismic field and to retrieve the balanced amplitude for any phases. Several data processing’s, such as the energy normalization or the spectral whitening have already been proposed for the special case of the HVSR (Spica et al. 2015; García-Jerez et al. 2016) but their reliability still need to be attested. From the theoretical point of view, the forward modeling generally considers a uniform layered medium in 3-D (e.g. Sánchez-Sesma et al. 2011a; Piña-Flores et al. 2016). Matsushima et al. (2014) also focused on the 2-D modeling by considering a directional dependency of the HVSR in the presence of smooth lateral heterogeneities. Following this approach, we investigate the use of directional HVSR at the edge of a volcanic crater; that is, when the heterogeneities are strong and apparent. We discuss the use of a suitable data processing and show the effect of the crater wall on the noise illumination by comparing with theoretical models in the 2-D and 3-D cases. We demonstrate that an appropriately projected HVSR is approximately represented when considering a 2-D model of unbounded layered media. The use of the discrete wave number (DWN) method for the forward modeling is discussed. The directional HVSR inversion process is finally discussed and performed. Because of the strong non-uniqueness of the problem (Piña-Flores et al. 2016), the inversion is realized by considering jointly a subset of measurements. The vertical outcrop of the crater wall additionally allows direct evaluation of the thicknesses parameters of the geological layers and then acknowledges for the reliability of the results. 2 KAWAH IJEN CRATER FEATURES Since the goal is to explore new domains of application of the HVSR method rather than to discuss the volcano geology, we only briefly introduce the geological and structural features relevant for this study. Additional information on the Kawah Ijen volcano geological features can be found in van Hinsberg et al. (2010), Caudron et al. (2015a), Spica et al. (2015) or in the geological map of Sujanto et al. (1988). Kawah Ijen is the youngest and only active volcano inside the large Ijen caldera. The Caldera of about 210 km2 is located on the eastern part of Java, Indonesia (see Fig. 1a). The volcano hosts the worldwide largest hyperacidic crater-lake [0.03 km3, pH ≈ 0.0; Caudron et al. (2017)] and presents an important and continuous degassing (Delmelle & Bernard 1994). The crater of this basaltic–andesitic composite volcano (Delmelle et al. 2000) has a regular elliptical shape of about 1.2 km by 1 km with steep vertical cliffs (∼20° to 30°) rising up to 250 m over the lake level (van Hinsberg et al. 2010). On the western side, a break in the crater rim extends almost to lake level and acts as an open door for the wind inside the crater (Caudron et al. 2015b). The eastern wall of the crater has an elevation of ∼180 m over the lake level and shows a subhorizontal layered stratigraphy (Figs 1b and 2). It consists of thick lava flows covered by layered phreatomagmatic pyroclasts and air fall deposits interbedded with thin lava flows (Fig. 2; Spica et al. 2015). Silica lacustrine superficial sediments cover an important part of the crater walls and attest for the former higher elevation of the lake. A thin surface of unconsolidated sulphur-bearing mud, derived from recent fumarolic activities, covers the entire crater rim (van Hinsberg et al. 2010). In order to have a full picture of the crater geology, the reader can refer to the digital elevation model of Gunawan et al. (2016). Figure 1. View largeDownload slide (a) Kawah Ijen volcano topography and locations of the seismic stations. Grey hexagons symbolize the permanent seismic stations of the Kawah Ijen network. Black inverted triangles refer to the positions of all the HVSR measurements. The location of the Ijen caldera on Java Island is reported by the white dot on the map in the upper left corner. (b) Zoom on the eastern part of the crater. Stations used for inversion are shown in black. Station codes are reported on the right-hand side. The tangent (T) to the crater at a station is given by the two contiguous stations on both side and is marked as a grey line. It is used to apply the rotation of the north–south (N) and east–west (E) directions and to obtain the radial direction (R). Figure 1. View largeDownload slide (a) Kawah Ijen volcano topography and locations of the seismic stations. Grey hexagons symbolize the permanent seismic stations of the Kawah Ijen network. Black inverted triangles refer to the positions of all the HVSR measurements. The location of the Ijen caldera on Java Island is reported by the white dot on the map in the upper left corner. (b) Zoom on the eastern part of the crater. Stations used for inversion are shown in black. Station codes are reported on the right-hand side. The tangent (T) to the crater at a station is given by the two contiguous stations on both side and is marked as a grey line. It is used to apply the rotation of the north–south (N) and east–west (E) directions and to obtain the radial direction (R). Figure 2. View largeDownload slide Picture of the northeastern part of the crater wall shown in Fig. 1(b) and sketch of the main geological layers. The numbers at the bottom of the figure refer to the station number. The bullet points at the right correspond to the different geological units as follows: (1) airfall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows [covered by silica lacustrine sediments (5)]. Red dots depict the depth interfaces obtained from the inversion as described in Section 5.2. Figure 2. View largeDownload slide Picture of the northeastern part of the crater wall shown in Fig. 1(b) and sketch of the main geological layers. The numbers at the bottom of the figure refer to the station number. The bullet points at the right correspond to the different geological units as follows: (1) airfall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows [covered by silica lacustrine sediments (5)]. Red dots depict the depth interfaces obtained from the inversion as described in Section 5.2. 3 INSTRUMENT AND DATA Ambient seismic field measurements were performed by moving a single seismic station every 100 m all over the Ijen crater between 2010 October 8 and 24. The station was composed of three-component Lennartz LE-3-D seismometer with natural period of 5 s and a CITYSHARK digital acquisition system sampling each channel at 100 Hz. When the sensor was installed on soft ground, a 1.2-cm buried-concrete plate was used in order to optimize the soil–sensor coupling (Bard & SESAME-team 2004). All recordings were performed during 20 min, a duration verified to be sufficient to ensure the stability of the DED measurements at Kawah Ijen volcano and for the frequency band considered here (Spica et al. 2015). In this research, we focus on four measurements located on the eastern part of the crater (Figs 1 and 2). On this area, the geology varies smoothly from one measurement to another, and the same number of layers can be considered, allowing a joint inversion of all the measurements (see Section 5.1). HVSR inversion results for other stations were already presented in Spica et al. (2015), and the shallow geology of western part of the crater is too complex for the theoretical assumptions made in this study (see Section 4.2). The ambient seismic field contains most of its energy below 6 Hz (see Fig. 3). Since ambient seismic field recorded at stations in the eastern part of the crater are not influenced by anthropogenic sources (Spica et al. 2015), the origin of this ‘noise’ is most likely due to the wind hitting the crater walls (Caudron et al. 2015a). Caudron (2013) showed that unless during volcanic crises, the spectral contents, polarization and wavefield compositions do not significantly change at Kawah Ijen volcano. Based on Caudron (2013), we know that HVSR measurements were obtained when volcano-related seismic activity was at its quiet background level (Spica et al. 2015). Figure 3. View largeDownload slide Power spectra of the ambient seismic field component (Z, E, N) recorded at station s6. Figure 3. View largeDownload slide Power spectra of the ambient seismic field component (Z, E, N) recorded at station s6. 4 HV CALCULATION According to its definition, the H/V spectral ratio corresponds to the square root of the spectral energies ratio of the horizontal directions (with indices 1, and 2) over the vertical direction (index 3, Arai & Tokimatsu 2004):   $$\frac{H}{V}(\mathbf {x},\omega ) = \sqrt{\frac{E_1(\mathbf {x},\omega )+E_2(\mathbf {x},\omega )}{E_3(\mathbf {x},\omega )}}.$$ (1)At an observed location, the spectral energies, the so-called DED in Perton et al. (2009), are proportional to the average autocorrelations of the diffused wavefield components, which in turn are proportional to the imaginary parts of the GF (ImGF):   $$E_i(\mathbf {x},\omega ) = \rho \omega ^2\left\langle \big | u_i(\mathbf {x},\omega )\big |^2 \right\rangle \propto -\omega ^D\mbox{Im} \big [ \mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega ) \big ]$$ (2)where ω = 2πf is the angular frequency, ui(x, ω) is the displacement field in the i direction at a point x. Im[ ] stands for the imaginary part operator, $$\mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )$$ is the displacement GF 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 exponent D is equal to 2 in 2-D and to 1 in 3-D (Sánchez-Sesma & Campillo 2006). The quantity $$\big | u_i(\mathbf {x},\omega )\big |^2 = u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )$$, with * the conjugate operator, corresponds in frequency domain to the autocorrelation of the displacement field in time domain. In what follows, the dependence in ω and x is implicit. The ∝ means that the expressions are proportional by a factor C that is independent from ω and x. The brackets 〈 〉 represent the average operator. Theoretically, it is realized over all the possible modes of a equipartition field. In fact, in the context of the DFA, this equation is only verified when the seismic wave field is equipartitioned, that is, all the incident waves of a same state (e.g. P, S and surface waves) have the same energies (Perton et al. 2016). However, there is no way to experimentally verify the equipartition hypothesis since we only have access to the whole field, where all the states and incident, transmitted, reflected and diffracted are indiscernible. Then the average is practically realized on sliced time windows that attempt to isolate each contribution associated to the incidence of a single mode. ρ is the mass density of a medium very far away where are originated the incident waves. Its value, as well as the value of the proportional factor C may not be easily fixed with real data but they can be determined for the theoretical cases (Perton et al. 2016). 4.1 Experimental HV Although equipartition is an assumption difficult to verify and surely not true in the majority of ambient seismic field measurements, several signal processing’s exist to reinforce equipartition of the seismic wavefield (e.g. Bensen et al. 2007). Here, we only apply the operation of spectral whitening as in Spica et al. (2015). Because several sources can act in different frequency bands, the whitening consists in normalizing the signals by the source energies computed in each time window (i.e. source deconvolution) and in several frequency bands as $$\tilde{v}_i(\mathbf {x},\omega ) = v_i(\mathbf {x},\omega )/\sqrt{\sum _{i = 1}^3|v_i(\mathbf {x},\Delta \omega )|^2}$$. Δω is a frequency band of about 2.5 Hz width and centred on ω. This bandwidth is larger than the width of the peaks in Fig. 3 (since they are related to the GF) and narrow enough to remove the spectral envelope due to the seismic sources. We use the particle velocity vj(x, ω) = iωuj(x, ω) to keep the link with the energy. Experimentally, the H/V ratio is computed in terms of the wavefield autocorrelations as:   \begin{eqnarray} \frac{H}{V}^{\text{3-D}} &=& \sqrt{\frac{\langle \big | \tilde{v}_1\big |^2 \rangle +\langle \big | \tilde{v}_2\big |^2 \rangle }{\langle \big | \tilde{v}_3\big |^2 \rangle }}. \end{eqnarray} (3)The average is computed on the spectra obtained over several (≈80) time windows of 20 s duration with an overlap of 75 per cent. By denoting the square root horizontal energy in a single time window as $$\tilde{H}_W = \sqrt{ (|\tilde{v}_1(\mathbf {x}, \omega )|^2 + |\tilde{v}_2(\mathbf {x}, \omega )|^2)}$$ and the square root vertical energy as $$\tilde{V}_W = \sqrt{|\tilde{v}_3(\mathbf {x}, \omega )|^2}$$, eq. (3) is then related to the ratio of the averages: $$\sqrt{{\langle \tilde{H}_W\rangle }/ { \langle \tilde{V}_W\rangle }}$$. In that sense, it is different from the usual calculation of the HVSR (Nakamura 1989) that is related to the average of the ratios: $$\sqrt{\langle { H_W}/{V_W}\rangle }$$ (without spectral whitening). Evaluating the ratio first and then the average (e.g. Nakamura 1989) allows the source deconvolution in each time window; however, it removes the link with the GF, which is shown here to be of interest. Furthermore, the deconvolution has the disadvantage to create instabilities that are generally overcome by applying a smoothing on VW (e.g. Konno & Ohmachi 1998). Besides, in our approach, the effect of the source deconvolution is averaged over many windows and is statistically more stable. We present in Fig. 4 the velocity-averaged autocorrelations in the three directions of a local polar frame; where er and eθ directions are the radial and tangential directions in the horizontal plane as shown in Fig. 1(b). Because of the presence of the crater, the tangent and perpendicular components differ. By computing the autocorrelations along all the directions of the horizontal plane by steps of 5°, we observed that the autocorrelation in the eθ direction corresponds almost to the minimum of the energy and that autocorrelation in the er direction is almost maximum. In fact, the minimum is 20° apart from the tangent direction and its value is only 6 per cent lower. Additionally, we add the upper and lower bounds of the autocorrelations (Fig. 4). They are the maximum and minimum values at each frequency from all the autocorrelations computed with half of the windows. These results justify that the 20 min duration of recordings allows the convergence of the DED. Since the velocities autocorrelation are proportional to the ImGF by the same unknown factor −CωD/ρ (see eq. 2), we present them normalized with respect to the maximum of the three curves. This allows comparing relatively the autocorrelation with the ImGF, and determining the factor D according to the observed frequency trend (see Section 4.2); that is, determining the space dimension in which the field is diffused. This identification of D is impossible when observing only the HVSR because the factor of proportionality simplifies in the HVSR (see eq. 3). Figure 4. View largeDownload slide Autocorrelations at station s6 for the tangential (dotted black), radial (dashed grey) and vertical (continuous black) velocity components with their respective upper and lower bounds calculated with half the windows number (thin curves). Figure 4. View largeDownload slide Autocorrelations at station s6 for the tangential (dotted black), radial (dashed grey) and vertical (continuous black) velocity components with their respective upper and lower bounds calculated with half the windows number (thin curves). 4.2 Theoretical HV According to eq. (2), the H/V ratio can be computed in terms of the ImGF (Sánchez-Sesma et al. 2011a):   $$\frac{H}{V} = \sqrt{\frac{\mbox{Im}(\mathcal {G}_{11}+\mathcal {G}_{22})}{\mbox{Im}(\mathcal {G}_{33})}}.$$ (4)In this expression, the ImGF components are associated to a specific geometry and to material properties that explain adequately the experimental data. This configuration is detailed in what follows. 4.2.1 Identification of a simplified geometry for the forward problem The soil properties assessment requires simulating the propagation of elastic waves in a geometrical configuration representative of the real one; that is, by considering irregular layers, the presence of the crater wall and the heterogeneities inside the layers. However, such complex simulation (i.e. computationally achievable with heavy finite-element method, FEM) would involve too many parameters and would take too much time to provide a valid solution. Therefore, several simplifications are proposed, following four important assumptions. First, the heterogeneities inside the layers are discarded supposing that the local wavelengths are larger than the typical size of the heterogeneities by choosing an appropriate frequency band. This assumption is justified according to the assessed velocities (Section 5.2) and according to the size of the heterogeneities that are observable in Fig. 2. Second, we assume that the layers are all homogeneous, elastic, isotrope and horizontally flat. This assumption is made locally, according to the small variation of the thicknesses in horizontal directions in function of the wavelength (Fig. 2). The top surface is a free surface and the last layer is the half-space. Third, we investigate the effect of the wall over the ambient seismic source directionality by conducting numerical simulations in 3-D and 2-D and we observe that the 2-D model better represents the observed data. Fourth, we show from a theoretical approach that the HVSR and H obtained for a 2-D unbounded layered medium (i.e. without the crater wall) and when introducing a correction factor, reproduce approximately the HVSR and the H obtained when considering the complex geometry with the crater. In order to validate the third assumption, we consider a multilayered medium with the presence of a crater, respectively in 3-D and 2-D. The crater is (axi-)symmetric and surrounded by a multilayer medium made of four layers (with respective thicknesses of 18, 10, 9 and 54 m) on top of a half-space. The associated velocities for the longitudinal waves are α = 443, 6350, 647, 2200 and 4000 m s−1 and for the shear waves β = 213, 3450, 353, 1060 and 2230 m s−1. These values result from average wave velocities in rocks (Table 1) described in Fig. 2. Since they are also close to the values we obtained after inversion (Table 1), the simulations are pertinent here. The crater wall has a slope of about 30°, as for Ijen. The three imaginary part of the Gii components are calculated for superimposed source and receiver located at d = 15 m from the crater wall, as in the real case. A section of the geometry is presented in Fig. 5. Figure 5. View largeDownload slide Axisymmetric section of the geometry considered for the modelization of the H/V near a crater wall (case a). Figure 5. View largeDownload slide Axisymmetric section of the geometry considered for the modelization of the H/V near a crater wall (case a). Table 1. Comparison between average observed shear wave velocities in rocks and inverted shear velocities at Ijen. The range of velocities in the second column is based on Christensen (1982), Christensen et al. (1980), Christensen & Stanley (2003), Mavko (n.d.), Mavko et al. (2009) and Jolly et al. (2012). Layer number (from top)  Rock types  Average shear velocity (m s−1)  Inverted shear velocity (m s−1)  1 and 3  Aifall/scoria  100–600  200–400  2 and 5  Lava flows  2800–3500  2700–3400  3  Phreatomagmatic  300–800  850  6  Unknown/altered rock  .  200–400  Layer number (from top)  Rock types  Average shear velocity (m s−1)  Inverted shear velocity (m s−1)  1 and 3  Aifall/scoria  100–600  200–400  2 and 5  Lava flows  2800–3500  2700–3400  3  Phreatomagmatic  300–800  850  6  Unknown/altered rock  .  200–400  View Large The computation of the ImGF as of the HVSR is performed using the recent formulation of the indirect boundary element method (IBEM, Perton et al. 2016) that allows simulating effectively the wave propagation in piecewise irregular regions and unbounded layered media. This formulation discretizes only the crater where the virtual sources of the IBEM formulation emit directly in the unbounded stratified media. This limits drastically in 3-D the size of the model and avoids the presence of spurious wave reflections that may occur at the edges of the geometrical model in FEM or traditional IBEM or BEM. These reflections, even attenuated by appropriate absorbing boundaries may slow down the computation or affect strongly the spectra and would compromise the computation of the HVSR. In addition, to avoid memory saturation at high frequency, the radius of the crater is reduced to 250 m comparing to real configuration. 3-D configuration: the HVSR (eq. 4), as well as the square root energy densities $$H^{\text{3-D}} = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{11}+\mathcal {G}_{22})}$$ and $$V^{\text{3-D}} = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{33})}$$ are shown as continuous lines in the Fig. 6. Additionally, we superimposed the HVSR and the H obtained when the square root horizontal energy density is restrained to the separate contributions in the radial $$H^{\text{3-D}}_r = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{rr})}$$ (dotted lines) and tangential $$H^{\text{3-D}}_{\theta } = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{\theta \theta })}$$ (dashed lines) directions. Figure 6. View largeDownload slide HVSR, H and V curves for case a at d = 15 m for the 3-D geometry. Full 3-D is in continuous black lines, and curves when horizontal energy is restrained to eθ and er are shown, respectively, in dashed and dot black lines. H and V curves are normalized with respect to the maximum of V, since the absolute amplitude is of no importance here. Figure 6. View largeDownload slide HVSR, H and V curves for case a at d = 15 m for the 3-D geometry. Full 3-D is in continuous black lines, and curves when horizontal energy is restrained to eθ and er are shown, respectively, in dashed and dot black lines. H and V curves are normalized with respect to the maximum of V, since the absolute amplitude is of no importance here. The H/V presents some usual characteristics: a large peak followed by a minimum that are well described in Piña-Flores et al. (2016). Generally, little attention is paid to the independent H and V curves or even to the two H components $$H^{\text{3-D}}_r$$ and $$H^{\text{3-D}}_{\theta }$$, and here they disagree strongly with the observations. The amplitude of the main peak of $$H^{\text{3-D}}_{\theta }/V$$ is about 1.5 times the one of $$H^{\text{3-D}}_r/V$$, in contradiction with the autocorrelation amplitude shown in Fig. 4, where $$\langle | \tilde{v}_r |^2 \rangle$$ is generally much larger than $$\langle | \tilde{v}_{\theta } |^2 \rangle$$. Furthermore, the curves H and V show a trend proportional to the frequency, rather the curves in Fig. 4 have a trend almost flat. Considering that the wind is the main source of noise in that frequency band, we reasonably assume that it corresponds to a coherent source all along the wall, mainly applying a force in the direction er over a large surface. Comparatively, the forces along ez and eθ should be much lower on the wall due to symmetry except at the top surface for Fz where symmetry breaks. This assumption is supported by the results shown in Fig. 4 and by the fact that it exists a direction that minimizes strongly one of the horizontal energy. Another argument comes from the wind direction measurements realized in 2010. A WindSonic anemometer connected to a Campbell CR200 data logger was programmed to record the mean and maximum wind speeds and mean speed direction. The station was installed close to the dam, which is a narrow corridor of wind (Fig. 1) and the only entry of wind inside the crater. In average, the wind speed was about 2 m s−1 with an angle of 72$${^{\circ}_{.}}$$6 from north with a standard deviation of 51$${^{\circ}_{.}}$$8. Supposing that the wind kept its direction inside the crater until reaching the stations, the stations were almost aligned perpendicularly to the wind direction. More details about the wind measurements can be find in Caudron (2013). We then conclude that the 3-D configuration is not suitable for the interpretation of the measurements. 2-D configuration: we consider now the wave propagation in the 2-D plane (er, ez). The square root energy densities in 2-D and 3-D configurations differ of a factor ω (eq. 2), so that they are now evaluated as $$H^{\text{2-D}} = \sqrt{-\omega ^2 \,\text{Im}(\mathcal {G}_{rr})}$$ and $$V^{\text{2-D}} = \sqrt{-\omega ^2 \,\text{Im}(\mathcal {G}_{33})}$$. Note that $$\mathcal {G}_{rr}$$ is a 2-D GF and that the index r indicates only the radial direction at the source–receiver position. The HVSR simplifies in $${H/V}^{\text{2-D}} = \sqrt{{\mbox{Im}(\mathcal {G}_{rr})}/{\mbox{Im}(\mathcal {G}_{33})}}$$. These quantities are shown as black lines in the Fig. 7 for two distances d = 5 m (dashed) and d = 15 m (continuous). H and V are normalized with respect to their respective maximum because the factor C is unknown. The H/V2-D differs from its 3-D counterpart, with a step around 2 Hz. Comparing to the 3-D case, the H and V present sharper peaks with maxima located at different frequencies and a different frequency trend, which agrees well with the velocity autocorrelation (Fig. 4). These differences are due to different diffraction regimes in 2-D and 3-D that modify the bulk to surface wave energy ratio (Perton et al. 2009, 2016). This configuration is then considered more adapted to represent the data and confirm the 2-D characteristic of the diffuse field. Figure 7. View largeDownload slide Comparison of the HVSR, H and V curves for case a at d = 5 m (continuous black line) and d = 15 m (dashed black line) and for case b (red lines) in 2-D. H and V have been normalized to their respective maximum. Figure 7. View largeDownload slide Comparison of the HVSR, H and V curves for case a at d = 5 m (continuous black line) and d = 15 m (dashed black line) and for case b (red lines) in 2-D. H and V have been normalized to their respective maximum. As discussed previously, the wave propagation simulation in geometries involving lateral heterogeneities, even using the IBEM method in a 2-D configuration, is not conceivable for an inverse problem. Given the positions of the source and receiver and the frequency range, we observed on simulation that the waves of interest for the HVSR only undergo reflections on one side of the crater wall and do not even probe the horizontal floor. In order to validate the fourth assumption, we consider the simulation of wave propagation in a multilayered medium only; that is, without the crater, and we add a correction factor in order to account from the wave reflection on a vertical wall (we neglect the slope of the crater wall). This configuration is denoted as case b in opposition of the previous configuration denoted as case a. The calculation for case b is realized using the DWN method (Bouchon 2003). This method is described in detail for the 3-D case in Sánchez-Sesma et al. (2011a) and for the 2-D case in the Section 4.2.2. The radiation pattern, as the sensitivity to reflected waves is different for H and V components. For H, that is, when the force is horizontal and when the displacement measured is also horizontal, the measured reflected waves are mainly longitudinal waves. Since the distance d is small with respect to the wavelength in case a, we can assume almost plane wave behaviour. In this context, the wave reflection imposes that the displacements in case a are twice the one in case b, when d = 0. However, for d ≠ 0, the incident and reflected waves are no more superimposed. We use then a correction factor of 2 cos (kd); with k = α/ω being the wave number associated to the longitudinal wave of velocity α. When introducing this factor to the H spectra in case b, the HVSR and H for d = 5 and 15 m are almost identical to their respective curve in case a. No correction was found for the V component, mainly because in this case the wall produces the interference of incident and reflected surface waves that cannot be approximated with a simple factor. Discrepancies only occur for V at d = 15 m above 5.5 Hz and do not affect the HVSR. The 2-D unbounded model is then suitable to reproduce the projected experimental data and is used in what follows as the forward modeling. 4.2.2 Forward problem method Several methods have already been proposed for the particular modelization of HVSR in multilayers and under the DFA. The DWN method (Bouchon 2003) was first proposed in Sánchez-Sesma et al. (2011b). Sophisticated methods such the one proposed in García-Jerez et al. (2016) or Perton & Sánchez-Sesma (2016) allow an efficient calculation of wave propagation in layered medium by being able to differentiate the integration of bulk and surface wave. However, the DWN method being popular, we show here that it is as efficient as the aforementioned methods in the special case when source and receiver coincide. We only summarize the method in the 2-D case with in-plane displacement, since the method was already presented for the 3-D case in Sánchez-Sesma et al. (2011b). In this method, the wavefield is decomposed as a sum of plane waves. The decomposition is realized according to the horizontal component of the wave vector kx because this component is conserved by reflection and transmission. A compact form of the displacement component $$u_i^n (\mathbf {x},\mathbf {x}_s )$$ inside the layer n, at receiver located in x and for a point source in xs and oriented in direction j is:   \begin{eqnarray} u_i^n(\mathbf {x},\mathbf {x}_s)& = &\delta (z,z_s)G_{ij}(\mathbf {x},\mathbf {x}_s) + \int _{-\infty }^{+\infty } \left( \begin{array}{c}S^{\nearrow Pn} u_i^{\nearrow Pn} e^{ik_z^{Pn} z}+ S^{\nearrow SVn} u_i^{\nearrow SVn} e^{ik_z^{SVn} z}+\\ S^{\searrow Pn} u_i^{\searrow Pn} e^{-ik_z^{Pn} z} +S^{\searrow SVn} u_i^{\searrow SVn} e^{-ik_z^{SVn} z} \end{array}\right) e^{-ik_x (x-x_s)}dk_x \end{eqnarray} (5)For an extended description of the integrand, the interested reader can consult the textbook of Aki & Richards (2002). Gij(x, xs) is the incident field as in an unbounded homogeneous medium that has to be taken into account only if the source and the receiver belong to the same layer (δ(z, zs)). $$k_z^{Pn} = \sqrt{{k^{Pn}}^2-k_x^2}$$ (respectively, $$k_z^{SVn}$$) is the vertical component of the wavenumber kPn = ω/αn (respectively, kSVn = ω/βn) associated with the P (respectively, SV) wave of the layer n with a negative sign of its imaginary part. uPn and uSVn are unitary vector of polarization of the waves. Here αn and βn are the longitudinal and shear velocities of the layer n. Then, the wave amplitudes S↗Pn, S$${\searrow}$$Pn, S↗SVn, S$${\searrow}$$SVn are solved for each discretized values of kx by considering the boundary equations: displacement and traction continuities. This equation simplifies when source and receiver coincide at the top of the layered media as:   $$u_i^1 (\mathbf {0},\mathbf {0}) = G_{ij} (\mathbf {0},\mathbf {0})+\int _{-\infty }^{+\infty } \left(S^{\nearrow P} u_i^{\nearrow P}+S^{\searrow P} u_i^{\searrow P}+S^{\nearrow SV} u_i^{\nearrow SV}+S^{\searrow SV} u_i^{\searrow SV}\right) dk_x$$ (6)The integrand presents poles associated to the surface waves. These pole contributions are usually very complicated to evaluate numerically and require using a fine discretization of kx and to replace the real frequency by a complex one: ω = ωR + iωI (Bouchon 2003). Real and imaginary parts of the integrand of G33(0, 0; ω = 6π) are presented in Fig. 8. Figure 8. View largeDownload slide Real and imaginary parts of the integrand amplitude of the component Green’s function G33(0, 0; ω) in function of horizontal wavenumber (expressed in km−1). The amplitudes are compressed using the formulae indicated in the legend. This modification preserved the sign of the integrand. Figure 8. View largeDownload slide Real and imaginary parts of the integrand amplitude of the component Green’s function G33(0, 0; ω) in function of horizontal wavenumber (expressed in km−1). The amplitudes are compressed using the formulae indicated in the legend. This modification preserved the sign of the integrand. The geometry considered is a layer (α1 = 2000 m s−1, β1 = 1000 m s−1, ρ1 = 1000 kg m−3 and h1 = 1 km) on top of a half-space (αHS = 4000 m s−1, βHS = 2000 m s−1 and ρHS = 1000 kg m−3). In order to better visualize the integrand, we use the logarithm function of the amplitude and a compression factor. However, the sign of the amplitudes is preserved. According to eq. (2), the ImGF has to be negative in sign. Here, the direct source contribution, that is, Gij(0, 0), is not considered and only the integrand is presented. Therefore, no restrictions on the sign of the integrand are assumed. The contribution of the bulk waves to the GF is delimited by the value kx = ω/β1 ≈ 18.85 km−1 and corresponds to the smooth amplitude variation below the pole peaks. Symmetries allow simplifying the calculation for positive value of kx only. Contrary to the imaginary part, the real part presents a change of sign around the poles. When evaluating numerically the continuous integral in eq. (6) by a discrete sum, this sudden change of sign makes the numerical integration unstable. Also, the integration limit has to be very large for the real part but not for the imaginary part. Same properties are observed for G11(0, 0; ω). Besides, the non-diagonal components of the GF show the opposite behaviour: the real part of the integrand are easy to integrate and not the imaginary parts. The evaluation of the imaginary part is then much easier than of the real part. This is an important observation for the HVSR and ImGF calculation and as far as we observe, it is not dependent to the layer properties. This is only the case when source and receiver coincide. Otherwise, the term $$e^{-ik_x (x-x_s ) }$$ in eq. (5) makes the previous real and imaginary parts mixed, so that difficulties encountered for real part integration also arise for the imaginary part integration. With an appropriate choice of ωI, the discretization over kx can be coarse and the integral converges quickly. The consequences of the introduction of ωI in the integration is corrected by calculating the time response through inverse Fourier transform multiplied by $$e^{-\omega _I t}$$ and then by calculating back the associated spectrum through direct Fourier transforms. Therefore, the DWN method offers a fast way to compute the HVSR following eq. (4). 5 DATA INVERSIONS 5.1 HVSR and H inversion The HVSR inversion has recently been extensively discussed (Lontsi et al. 2015; García-Jerez et al. 2016; Piña-Flores et al. 2016; Sánchez-Sesma 2017). As shown in Piña-Flores et al. (2016), a single consideration of the HVSR is insufficient for characterizing soil properties since a scale factor for velocities and thicknesses leads to the same HVSR. The additional consideration of dispersion curves helps to reduce the number of solutions by constraining the velocity model but they are not assessable with our data set. We rather explored the joint inversion of four HVSR at nearby positions, with different layer thicknesses and by considering the same elastic properties inside the layers. The thicknesses are not known with precision but they can be evaluated from the outcrop in Fig. 2. Although the thicknesses may slightly vary from the outcrop because the seismic stations are 15 m away from the crater wall, it represents an important a priori information to reduce the range of thicknesses in our inversion. The joint inversion that considers same velocities for different thicknesses avoids the multiple solutions related to the scale factor. Since the HVSR computation under DFA theory also implies the separated computation of H and V components (see Sections 4.1 and 4.2.1), they also can be used during inversion. As shown for the 3-D case in Piña-Flores et al. (2016) and Sánchez-Sesma et al. (2017), the GFs or the H and V associated to a layered media display large oscillations at high frequency related to the body waves over a trend proportional to frequency. Since the oscillations amplitudes are small regarding to the trend, they appear flatten in the HVSR, and the information that can be retrieved is mainly due to the surface waves at low frequency; however surface waves are not linearly sensitive to the thicknesses and to the body wave velocities. Then, the HVSR inversion leads to non-unique solutions. In the opposite, the use of H and V during the inversion may improve the sensitivity directly to bulk wave velocities, diminishing the number of solutions and may also be able to better characterize the shallow structure due to information collected at high frequency. This is also true in 2-D; however, since the frequency trend in GF is flat [see eq. 2 and Sánchez-Sesma et al. (2017)], the oscillations related to bulk wave are more pronounced in the HVSR and HVSR are then sensitive to bulk wave velocities as well. Although the amplitude retrieval is still a controversial topic in ambient seismic field correlation (e.g. Prieto et al. 2011; Weaver et al. 2011), several studies confirmed the good retrieval of the phase through the retrieval of surface or bulk wave velocities (e.g. Shapiro et al. 2005; Roux et al. 2005). This proves at least that the amplitude associated to individual phases might be retrieved correctly in narrow frequency bands. On larger frequency band, they might be modulated by a function, related to source envelopes, with variations slower than those of the GF (as discussed in Section 4.1). The results presented in Piña-Flores et al. (2016), such as the amplitude inversion of the HVSR, suggest that the ratio of the ImGF reduces the influence of these envelopes. It might be because of the division of the potentially same envelopes or because of the special configuration of the superimposed source and receiver that reduce the sensitivity to non-isotropic illumination or to attenuation. Here, we assume that the spectral whitening (using narrow frequency band as discussed in Section 4.1), removes any envelope. Consequently, the inversion is realized by considering jointly the H component and the HVSR. The V component is discarded because it has been shown in Fig. 7 that the one obtained in the simplified unbounded layered media might be different at high frequency to the one obtained in the configuration that includes the wall; however, if we used the more sophisticated model that integrates the crater wall for the forward model, the V should also be considered in the inversion. Then, the soil properties and the layer thicknesses are assessed by inverting jointly the HVSR and the H at several positions. The objective function (ε) is defined as the root mean square difference between the experimental HVSR and the H spectra obtained by considering only the horizontal displacement in the radial direction and the HVSR and the H spectra computed with the DWN for the 2-D case:   $$\epsilon = \sqrt{\frac{1}{N_x N_{\omega } }\sum _{x_i}\sum _{\omega } \big ( \big(H_r^{\text{exp}}\big/V^{\text{exp}}-H^{2\text{-D}}\big/V^{\text{2-D}}\big)^2+\big(H_r^{\text{exp}}-H^{\text{2-D}}\big)^2 \big ) }.$$ (7)In this expression, the theoretical and experimental HVSR spectra are normalized according to the maximum amplitude in the experimental spectra in order to balance the error weight with the H spectra, which are also normalized as discussed in Section 4.1 since the constant C is not known. Nx = 4 is the number of positions xi considered in the joint inversion and Nω = 30 is the number of points taken linearly separated in the spectra from 1.5 to 6 Hz. The lower frequency bound is taken according to the results observed in the HVSR in Fig. 7. Curves at higher frequency are associated to small and local heterogeneities, which are of minor importance here. Rather than using the simulated annealing or genetic algorithm methods discussed in García-Jerez et al. (2016) and Piña-Flores et al. (2016), we use a ‘Pattern-Search’ method (Audet & Dennis 2002) because of its efficiency for highly non-linear problem and its ability to integrate constraints. Beside the algorithm has the drawback to fail to converge to the global minimum when the solution space is very large, the non-uniqueness of the solution is strongly reduced thanks to the large number of constraints considered. The elastic velocities β and the layer thicknesses at each position are the free parameters. The thickness of each layer is allowed to vary by 20 per cent only from its original a priori value. No variation restriction is considered for the shear wave velocities but they are taken constant all along the geological layers for the different positions. During the first inversions, α was also taken as a free parameter but results showed that it might be equally taken as α = 1.9β for all the layers. This value is consistent with values reported in the references given in Table 1. This simplification allowed reducing the number of free parameters and obtaining better fit. The density is related to α through the linear relationship: ρ = 0.32α + 0.77 with α expressed in km s−1 and ρ in kg m−3 (Berteussen 1977). 5.2 Results and discussion Results of the inversion for the subset of stations s1, s5, s6 and s8 are shown in Fig. 9. The H and V spectra are shown normalized. An important feature of these results is that both the amplitude and the entire shape of the HVSR and H curves (i.e. not only the main peak) are well retrieved for all the stations and for frequencies ranging from 1.5 to 6 Hz. As expected for the V curves, the frequency peak positions agree well until 4 Hz, but clear discrepancies are visible at higher frequencies. For the HVSR curves, the agreement is still acceptable for higher frequencies, but it is due to the absence of local peak. For the H curves, a peak around 4 Hz differs from the experimental ones, probably meaning that a superficial and thin layer might be missing in the model. The top layers could then be refined according to the high number of sublayers observed in the outcrop but their thickness evaluation is uncertain. Figure 9. View largeDownload slide Results of the joint HVSR inversion for the subset of stations s1, s5, s6 and s8 (shown in Fig. 2). Black and red dotted lines depict the target and best fit (ε0) modeled HVSR, H and V curves, respectively. The velocity models associated to the best fit are plotted with thick dotted red lines in the last column. Other profiles with higher misfit are also shown in continuous thin lines. The scale is from blue (≈ε0) to red lines (≈1.5ε0). The mean lake level is depicted as the cyan area. Figure 9. View largeDownload slide Results of the joint HVSR inversion for the subset of stations s1, s5, s6 and s8 (shown in Fig. 2). Black and red dotted lines depict the target and best fit (ε0) modeled HVSR, H and V curves, respectively. The velocity models associated to the best fit are plotted with thick dotted red lines in the last column. Other profiles with higher misfit are also shown in continuous thin lines. The scale is from blue (≈ε0) to red lines (≈1.5ε0). The mean lake level is depicted as the cyan area. The models associated to the best fit are superimposed in Fig. 2 and we can observe that they agree well with the local geology observed in the outcrop. As discussed in Spica et al. (2015), the outcrop of the north wall is composed by the following succession of layers (from top to bottom): (1) air-fall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows. These geological units are highlighted by numbered points in Fig 2. All these geological units are well represented in the velocity model and their thicknesses change in accordance with the outcrop characteristics. As an example, a strong impedance contrast is observed at the interface between the phreatic/superficial deposits (850 m s−1) lying on the thick lava flows (3400 m s−1; layer 4 over layer 5). The thin lava flow with an important velocity of 2700 m s−1 is particularly well observed in the models and in the outcrop. Although this layer is particularly thin, it appears remarkably well on the velocity model, suggesting the model is well assessed. The retrieved velocities are comparable to those obtained in a previous study of the Ijen crater (Spica et al. 2015) and charts of seismic velocities for these geological units (Table 1). Interestingly the half-space, referred by (5) in Fig. 2, corresponds to a layer in contact with the lake with a velocity lower than the above thick lava flow referred by the number (4). A possible interpretation is that it corresponds to the same geological unit but that has been strongly altered or corroded by the acidity of the lake making it more porous. Other profiles with higher misfit error are also depicted in Fig. 9. The four first layers have their thicknesses and velocities retrieved with confidence. The deepest layers are retrieved with large incertitudes, as expected, since they depend strongly to low frequency and then to the surface waves. These results validate the 2-D character of the ambient seismic source, the appropriate data processing and the importance of considering independently the H and V components. Here, the V component has been discarded only because of the presence of the strong lateral heterogeneity but might be considered elsewhere. The H and V components show important peak at higher frequency than the HVSR and allow then a better identification of the shallow layer properties. 6 CONCLUSIONS We explored here the effects of a strong lateral heterogeneity on the HVSR by using measurements along a vertical crater wall. We identified that the horizontal square root energy density H restrained to the direction normal to the crater wall is very close to the H in a 2-D layered medium. This is also the case for the V component but in a more limited frequency bandwidth. This identification is realized from theoretical simulations using IBEM when considering the crater embedded in a layered medium and the DWN method when considering only the same layered medium and with similar elastic properties as the one obtained after inversion. We demonstrated that the DWN method is efficient for the calculation of the ImGF when source and receiver are superimposed and may then be used for the forward modeling. The inversion was realized jointly by considering the HVSR curves and the independent H at several positions and with initial thicknesses roughly evaluated. The V energy density might also be used in other circumstances. The experimental H and V have been processed according to the theoretical assumption of a diffuse field, that is, by applying a spectral whitening and by conducting the mean autocorrelation as for GF retrieval. The good agreement of the results with the local geology validates the assumptions made such as the 2-D character of the ambient seismic source and the appropriate data processing. This paper also underlines that the H and V components could also be considered independently to the HVSR. It demonstrates then the potential of the method in complex geometrical configuration. This methodology paves the way for identifying layer elastic velocities at outcrop and then to follow their thickness variation far inside. Acknowledgements We warmly thank Thierry Camelbeeck and Thomas Lecocq from the Royal Observatory of Belgium for their contribution in the early stage of this research and the development of the Kawah Ijen monitoring project. The field data acquisition was funded by the Royal Observatory of Belgium. We also thank Cipta Anthanasius, Sarane Sterckx and Arnaud Watlet for their help on the field and data acquisition. We are grateful to the Indonesian Centre for Volcanology and Geological Hazard Mitigation support on the field, and particularly to the observers of Kawah Ijen observatory: Bambang Heri Purwanto and Suparjan. We appreciate the valuable suggestions of Thomas Bodin, Agostiny Marrios Lontsi and two anonymous reviewers. REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , University Science Book. Albarello D., Lunedei E., 2011. Structure of an ambient vibration wavefield in the frequency range of engineering interest ([0.5, 20] hz): insights from numerical modelling, Near Surf. Geophys. , 9( 6), 543– 559. Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am. , 94( 1), 53– 63. Google Scholar CrossRef Search ADS   Audet C., Dennis J.E., Jr, 2002. Analysis of generalized pattern searches, SIAM J. Optim. , 13( 3), 889– 903. Google Scholar CrossRef Search ADS   Bard P.Y., 1999. Microtremor measurements: a tool for site effect estimation, in The Effects of Surface Geology on Seismic Motion , 3, pp. 1251– 1279, Balkema, Rotterdam. Bard P.-Y., 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: ftp://delta.geo.uib.no/pub/seismo/REPORTS/SESAME/USER-GUIDELINES/SESAME-HV-User-Guidelines.doc, last accessed 1st November 2017. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. Google Scholar CrossRef Search ADS   Berteussen K.A., 1977. Moho depth determinations based on spectral-ratio analysis of NORSAR long-period P waves, Phys. Earth planet. Inter. , 15( 1), 13– 27. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Köhler A., Cornou C., Wathelet M., Bard P.-Y., 2008. Effects of Love waves on microtremor h/v ratio, Bull. seism. Soc. Am. , 98( 1), 288– 300. Google Scholar CrossRef Search ADS   Bouchon M., 2003. A review of the discrete wavenumber method, Pure appl. Geophys. , 160( 3–4), 445– 465. Google Scholar CrossRef Search ADS   Caudron C., 2013. Multi-disciplinary continuous monitoring of Kawah Ijen volcano, East Java, Indonesia, PhD thesis , Université Libre de Bruxelles, Brussels, Belgium. Caudron C., Lecocq T., Syahbana D.K., McCausland W., Watlet A., Camelbeeck T., Bernard A., Surono S., 2015a. Stress and mass changes at a ‘wet’ volcano: example during the 2011–2012 volcanic unrest at Kawah Ijen volcano (Indonesia), J. geophys. Res. , 120( 7), 5117– 5134. Google Scholar CrossRef Search ADS   Caudron C.et al.  , 2015b. Kawah Ijen volcanic activity: a review, Bull. Volcanol. , 77( 3), 1– 39. Google Scholar CrossRef Search ADS   Caudron C., Campion R., Rouwet D., Lecocq T., Capaccioni B., Syahbana D., Purwanto B.H., Bernard A., 2017. Stratification at the Earth’s largest hyperacidic lake and its consequences, Earth planet. Sci. Lett. , 459, 28– 35. Google Scholar CrossRef Search ADS   Christensen N.I., 1982. Seismic velocities, in Handbook of Physical Properties of Rocks , Vol. 2, pp. 1– 228. Christensen N.I., Stanley D., 2003. Seismic velocities and densities of rocks, Int. Geophys. , 81, 1587– 1594. Google Scholar CrossRef Search ADS   Christensen N., Wilkens R., Blair S., Carlson R., 1980. 10. Seismic velocities, densities, and elastic constants of volcanic breccias and basalt from deep sea drilling project LEG 59, in Initial Reports of the Deep Sea Drilling Project , 59, pp. 515– 517, U.S. Govt. Printing Office, Washington, DC. Delmelle P., Bernard A., 1994. Geochemistry, mineralogy, and chemical modeling of the acid crater lake of Kawah Ijen volcano, Indonesia, Geochim. Cosmo. Acta , 58( 11), 2445– 2460. Google Scholar CrossRef Search ADS   Delmelle P., Bernard A., Kusakabe M., Fischer T.P., Takano B., 2000. Geochemistry of the magmatic–hydrothermal system of Kawah Ijen volcano, east Java, Indonesia, J. Volcanol. Geotherm. Res. , 97( 1–4), 31– 53. Google Scholar CrossRef Search ADS   García-Jerez A., Piña-Flores J., Sánchez-Sesma F.J., Luzón F., Perton M., 2016. A computer code for forward calculation and inversion of the H/V spectral ratio under the diffuse field assumption, Comput. Geosci. , 97, 67– 78. Google Scholar CrossRef Search ADS   Gunawan H. et al.  , 2016. New insights into Kawah Ijen’s volcanic system from the wet volcano workshop experiment, Geol. Soc. Lond. Spec. Publ. , 437, doi:10.1144/SP437.7. Jolly A.D., Chardot L., Neuberg J., Fournier N., Scott B.J., Sherburn S., 2012. High impact mass drops from helicopter: a new active seismic source method applied in an active volcanic setting, Geophys. Res. Lett. , 39( 12), L12306, doi:10.1029/2012GL051880. Konno K., Ohmachi T., 1998. Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor, Bull. seism. Soc. Am. , 88( 1), 228– 241. Lermo J., Chávez-García F.J., 1994. Are microtremors useful in site response evaluation?, Bull. seism. Soc. Am. , 84( 5), 1350– 1364. Lontsi A.M., Sánchez-Sesma F.J., Molina-Villegas J.C., Ohrnberger M., Krüger F., 2015. Full microtremor H/V(z, f) inversion for shallow subsurface characterization, Geophys. J. Int. , 202( 1), 298– 312. Google Scholar CrossRef Search ADS   Lunedei E., Albarello D., 2010. Theoretical HVSR curves from full wavefield modelling of ambient vibrations in a weakly dissipative layered Earth, Geophys. J. Int. , 181( 2), 1093– 1108. Lunedei E., Malischewsky P., 2015. A review and some new issues on the theory of the H/V technique for ambient vibrations, in Perspectives on European Earthquake Engineering and Seismology , pp. 371– 394, Springer International Publishing. Google Scholar CrossRef Search ADS   Matsushima S., Hirokawa T., Martin F.D., Kawase H., Sánchez-Sesma F.J., 2014. The effect of lateral heterogeneity on horizontal to vertical spectral ratio of microtremors inferred from observation and synthetics, Bull. seism. Soc. Am. , 104( 1), 381– 393. Google Scholar CrossRef Search ADS   Mavko G., n.d. Conceptual Overview of Rock and Fluid Factors that Impact Seismic Velocity and Impedance , Stanford School of Earth Sciences . Available at: https://pangea.stanford.edu/courses/gp262/Notes/8.SeismicVelocity.pdf, last accessed 1st November 2017. Mavko G., Mukerji T., Dvorkin J., 2009. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media , Cambridge University Press. Google Scholar CrossRef Search ADS   Nakamura Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, Q. Rep. Railw. Tech. Res. Inst. , 30( 1), 25– 33. Perton M., Contreras-Zazueta M.A., Sánchez-Sesma F.J., 2016. Indirect boundary element method to simulate elastic wave propagation in piecewise irregular and flat regions, Geophys. J. Int. , 205( 3), 1832– 1842. Google Scholar CrossRef Search ADS   Perton M., Sánchez-Sesma F.J., 2016. Green’s function calculation from equipartition theorem, J. acoust. Soc. Am. , 140( 2), 1309– 1318. Google Scholar CrossRef Search ADS PubMed  Perton M., Sánchez-Sesma F.J., Rodríguez-Castellanos A., Campillo M., Weaver R.L., 2009. Two perspectives on equipartition in diffuse elastic fields in three dimensions, J. acoust. Soc. Am. , 126( 3), 1125– 1130. Google Scholar CrossRef Search ADS PubMed  Piña-Flores J., Perton M., García-Jerez A., Carmona E., Luzón F., Molina-Villegas J.C., Sánchez-Sesma F.J., 2016. The inversion of spectral ratio H/V in a layered system using the diffuse field assumption (DFA), Geophys. J. Int.  208( 1), 577– 588. Google Scholar CrossRef Search ADS   Prieto G.A., Denolle M., Lawrence J.F., Beroza G.C., 2011. On amplitude information carried by the ambient seismic field, C. R. Geosci. , 343( 8), 600– 614. Google Scholar CrossRef Search ADS   Roux P., Sabra K.G., Gerstoft P., Kuperman W.A., Fehler M.C., 2005. P-waves from cross-correlation of seismic noise, Geophys. Res. Lett. , 32( 19), L19303, doi:10.1029/2005GL023803. Sánchez-Sesma F.J., 2017. Modeling and inversion of the microtremor H/V spectral ratio: physical basis behind the diffuse field approach, Earth Planets Space , 69, 92, doi:10.1186/s40623-017-0667-6. Sánchez-Sesma F.J., Campillo M., 2006. Retrieval of the Green’s function from cross correlation: the canonical elastic problem, Bull. seism. Soc. Am. , 96( 3), 1182– 1191. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J.et al.  , 2011a. A theory for microtremor H/V spectral ratio: application for a layered medium, Geophys. J. Int. , 186( 1), 221– 225. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Weaver R.L., Kawase H., Matsushima S., Luzón F., Campillo M., 2011b. Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid, Bull. seism. Soc. Am. , 101( 4), 1704– 1709. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Iturrarán-Viveros U., Perton M., 2017. Some properties of Green’s functions for diffuse field interpretation, Math. Methods Appl. Sci. , 40( 9), 3348– 3354. Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307( 5715), 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Spica Z.et al.  , 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, J. Volc. Geotherm. Res. , 302, 173– 189. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Nakata N., Liu X., Beroza G., 2017. Site characterization at Groningen gas field area through joint surface-borehole H/V analysis, Geophys. J. Int. , doi:10.1093/gji/ggx426. Sujanto M., Syarifuddin M.Z., Sitorus K., 1988. Geological Map of The Ijen Caldera Complex, East Java, Volcanol. Surv. Indones . van Hinsberg V., Berlo K., van Bergen M., Williams-Jones A., 2010. Extreme alteration by hyperacidic brines at Kawah Ijen volcano, east Java, Indonesia: I. textural and mineralogical imprint, J. Volc. Geotherm. Res. , 198( 1–2), 253– 263. Google Scholar CrossRef Search ADS   Weaver R.L., Lobkis O.I., 2004. Diffuse fields in open systems and the emergence of the Green’s function (l), J. acoust. Soc. Am. , 116( 5), 2731– 2734. Google Scholar CrossRef Search ADS   Weaver R.L., Hadziioannou C., Larose E., Campillo M., 2011. On the precision of noise correlation interferometry, Geophys. J. Int. , 185( 3), 1384– 1392. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Inversion of the horizontal-to-vertical spectral ratio in presence of strong lateral heterogeneity

, Volume 212 (2) – Feb 1, 2018
12 pages

/lp/ou_press/inversion-of-the-horizontal-to-vertical-spectral-ratio-in-presence-of-A0GIQ7beHR
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx458
Publisher site
See Article on Publisher Site

### Abstract

Abstract The horizontal-to-vertical spectral ratio (HVSR) method has recently gained a resurgence of interest because of its new interpretation through the diffuse field approach. According to this theory, the HVSR measurement is linked to the Green’s function retrieval. While the HVSR method is traditionally used to evaluate simple and shallow velocity model, we demonstrate that complex Earth’ structure (that includes composite geological layering) can be assessed by choosing appropriate data processing and inversion scheme. In particular, we explore the effects of a strong lateral heterogeneity on the HVSR by using measurements along a vertical crater wall. However, wave propagation computation with lateral heterogeneity is inefficient to conduct an inversion. Then, to accelerate the forward calculation, we identified a simplified model consisting of an unbounded 2-D multilayer medium that allows describing the observed measurements with the use of a correction factor. The 2-D character comes from a non-standard noise illumination related to the wind blowing perpendicularly to the crater walls. It is supported by several theoretical and experimental arguments. We also show that the H and V components can be retrieved independently, which may improve the sensitivity of the inversion at high frequency. To reduce the non-uniqueness of the problem, we consider the joint inversion of the H and HVSR’s at several positions where the thicknesses of the layers were evaluated thanks to a large geological outcrop. The good agreement between the recovered velocity structure and the geology shows that the method is not limited to unbounded simple shallow structures, which opens the domains of application of the method to exploration geophysics or to hazard assessment. Joint inversion, Waveform inversion, Seismic noise, Site effects, Volcano seismology 1 INTRODUCTION The horizontal-to-vertical spectral ratio (the so-called HVSR or H/V) is a widely used method to determine the natural frequency of a site and to estimate a simple velocity model of the subsurface; that is, typically one layer over a half-space (e.g. Nakamura 1989; Lermo & Chávez-García 1994; Bard 1999). Such a simple interpretation comes primarily from the fact that the technique was developed empirically without clear theoretical basis. Actually, the fundamental peak of the HVSR curve was interpreted based on the assumption of predominance of body or surface waves in the ambient seismic field [see review of Lunedei & Malischewsky (2015); Piña-Flores et al. (2016)]. However, it is now well known that ambient seismic field is composed of both body and surface waves (e.g. Bonnefoy-Claudet et al. 2008) and therefore, an HVSR theory must account for it to explain the entire HVSR curve (i.e. not only the main peak). In this contribution, we base our argumentation on a recent theoretical developments that account for all the seismic phases, being then more robust to assess with details a velocity structure (Sánchez-Sesma et al. 2011a; Lontsi et al. 2015; García-Jerez et al. 2016; Piña-Flores et al. 2016; Spica et al.2017). So far, two different theories have been proposed for the modeling of the whole HVSR curve: the distributed surface sources (e.g. Lunedei & Albarello 2010; Albarello & Lunedei 2011; Lunedei & Malischewsky 2015) and the diffuse field approach (DFA, Sánchez-Sesma et al. 2011a,b; Piña-Flores et al. 2016). Even if these two branches of investigations have apparent solid divergence regarding their theoretical assumptions, other studies have already demonstrated their compliance through the equivalence of the reciprocity theorem and of the diffuse field theory (Weaver & Lobkis 2004; Perton et al. 2016). The DFA underlines the connection between the directional energy densities (DED), and the Green’s function (GF, Perton et al. 2009). This connection allows conducting the inversion of the HVSR measurement for physical properties evaluation. In this view, the data must be processed consistently with the theory in order to enhance the diffuse character of the ambient seismic field and to retrieve the balanced amplitude for any phases. Several data processing’s, such as the energy normalization or the spectral whitening have already been proposed for the special case of the HVSR (Spica et al. 2015; García-Jerez et al. 2016) but their reliability still need to be attested. From the theoretical point of view, the forward modeling generally considers a uniform layered medium in 3-D (e.g. Sánchez-Sesma et al. 2011a; Piña-Flores et al. 2016). Matsushima et al. (2014) also focused on the 2-D modeling by considering a directional dependency of the HVSR in the presence of smooth lateral heterogeneities. Following this approach, we investigate the use of directional HVSR at the edge of a volcanic crater; that is, when the heterogeneities are strong and apparent. We discuss the use of a suitable data processing and show the effect of the crater wall on the noise illumination by comparing with theoretical models in the 2-D and 3-D cases. We demonstrate that an appropriately projected HVSR is approximately represented when considering a 2-D model of unbounded layered media. The use of the discrete wave number (DWN) method for the forward modeling is discussed. The directional HVSR inversion process is finally discussed and performed. Because of the strong non-uniqueness of the problem (Piña-Flores et al. 2016), the inversion is realized by considering jointly a subset of measurements. The vertical outcrop of the crater wall additionally allows direct evaluation of the thicknesses parameters of the geological layers and then acknowledges for the reliability of the results. 2 KAWAH IJEN CRATER FEATURES Since the goal is to explore new domains of application of the HVSR method rather than to discuss the volcano geology, we only briefly introduce the geological and structural features relevant for this study. Additional information on the Kawah Ijen volcano geological features can be found in van Hinsberg et al. (2010), Caudron et al. (2015a), Spica et al. (2015) or in the geological map of Sujanto et al. (1988). Kawah Ijen is the youngest and only active volcano inside the large Ijen caldera. The Caldera of about 210 km2 is located on the eastern part of Java, Indonesia (see Fig. 1a). The volcano hosts the worldwide largest hyperacidic crater-lake [0.03 km3, pH ≈ 0.0; Caudron et al. (2017)] and presents an important and continuous degassing (Delmelle & Bernard 1994). The crater of this basaltic–andesitic composite volcano (Delmelle et al. 2000) has a regular elliptical shape of about 1.2 km by 1 km with steep vertical cliffs (∼20° to 30°) rising up to 250 m over the lake level (van Hinsberg et al. 2010). On the western side, a break in the crater rim extends almost to lake level and acts as an open door for the wind inside the crater (Caudron et al. 2015b). The eastern wall of the crater has an elevation of ∼180 m over the lake level and shows a subhorizontal layered stratigraphy (Figs 1b and 2). It consists of thick lava flows covered by layered phreatomagmatic pyroclasts and air fall deposits interbedded with thin lava flows (Fig. 2; Spica et al. 2015). Silica lacustrine superficial sediments cover an important part of the crater walls and attest for the former higher elevation of the lake. A thin surface of unconsolidated sulphur-bearing mud, derived from recent fumarolic activities, covers the entire crater rim (van Hinsberg et al. 2010). In order to have a full picture of the crater geology, the reader can refer to the digital elevation model of Gunawan et al. (2016). Figure 1. View largeDownload slide (a) Kawah Ijen volcano topography and locations of the seismic stations. Grey hexagons symbolize the permanent seismic stations of the Kawah Ijen network. Black inverted triangles refer to the positions of all the HVSR measurements. The location of the Ijen caldera on Java Island is reported by the white dot on the map in the upper left corner. (b) Zoom on the eastern part of the crater. Stations used for inversion are shown in black. Station codes are reported on the right-hand side. The tangent (T) to the crater at a station is given by the two contiguous stations on both side and is marked as a grey line. It is used to apply the rotation of the north–south (N) and east–west (E) directions and to obtain the radial direction (R). Figure 1. View largeDownload slide (a) Kawah Ijen volcano topography and locations of the seismic stations. Grey hexagons symbolize the permanent seismic stations of the Kawah Ijen network. Black inverted triangles refer to the positions of all the HVSR measurements. The location of the Ijen caldera on Java Island is reported by the white dot on the map in the upper left corner. (b) Zoom on the eastern part of the crater. Stations used for inversion are shown in black. Station codes are reported on the right-hand side. The tangent (T) to the crater at a station is given by the two contiguous stations on both side and is marked as a grey line. It is used to apply the rotation of the north–south (N) and east–west (E) directions and to obtain the radial direction (R). Figure 2. View largeDownload slide Picture of the northeastern part of the crater wall shown in Fig. 1(b) and sketch of the main geological layers. The numbers at the bottom of the figure refer to the station number. The bullet points at the right correspond to the different geological units as follows: (1) airfall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows [covered by silica lacustrine sediments (5)]. Red dots depict the depth interfaces obtained from the inversion as described in Section 5.2. Figure 2. View largeDownload slide Picture of the northeastern part of the crater wall shown in Fig. 1(b) and sketch of the main geological layers. The numbers at the bottom of the figure refer to the station number. The bullet points at the right correspond to the different geological units as follows: (1) airfall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows [covered by silica lacustrine sediments (5)]. Red dots depict the depth interfaces obtained from the inversion as described in Section 5.2. 3 INSTRUMENT AND DATA Ambient seismic field measurements were performed by moving a single seismic station every 100 m all over the Ijen crater between 2010 October 8 and 24. The station was composed of three-component Lennartz LE-3-D seismometer with natural period of 5 s and a CITYSHARK digital acquisition system sampling each channel at 100 Hz. When the sensor was installed on soft ground, a 1.2-cm buried-concrete plate was used in order to optimize the soil–sensor coupling (Bard & SESAME-team 2004). All recordings were performed during 20 min, a duration verified to be sufficient to ensure the stability of the DED measurements at Kawah Ijen volcano and for the frequency band considered here (Spica et al. 2015). In this research, we focus on four measurements located on the eastern part of the crater (Figs 1 and 2). On this area, the geology varies smoothly from one measurement to another, and the same number of layers can be considered, allowing a joint inversion of all the measurements (see Section 5.1). HVSR inversion results for other stations were already presented in Spica et al. (2015), and the shallow geology of western part of the crater is too complex for the theoretical assumptions made in this study (see Section 4.2). The ambient seismic field contains most of its energy below 6 Hz (see Fig. 3). Since ambient seismic field recorded at stations in the eastern part of the crater are not influenced by anthropogenic sources (Spica et al. 2015), the origin of this ‘noise’ is most likely due to the wind hitting the crater walls (Caudron et al. 2015a). Caudron (2013) showed that unless during volcanic crises, the spectral contents, polarization and wavefield compositions do not significantly change at Kawah Ijen volcano. Based on Caudron (2013), we know that HVSR measurements were obtained when volcano-related seismic activity was at its quiet background level (Spica et al. 2015). Figure 3. View largeDownload slide Power spectra of the ambient seismic field component (Z, E, N) recorded at station s6. Figure 3. View largeDownload slide Power spectra of the ambient seismic field component (Z, E, N) recorded at station s6. 4 HV CALCULATION According to its definition, the H/V spectral ratio corresponds to the square root of the spectral energies ratio of the horizontal directions (with indices 1, and 2) over the vertical direction (index 3, Arai & Tokimatsu 2004):   $$\frac{H}{V}(\mathbf {x},\omega ) = \sqrt{\frac{E_1(\mathbf {x},\omega )+E_2(\mathbf {x},\omega )}{E_3(\mathbf {x},\omega )}}.$$ (1)At an observed location, the spectral energies, the so-called DED in Perton et al. (2009), are proportional to the average autocorrelations of the diffused wavefield components, which in turn are proportional to the imaginary parts of the GF (ImGF):   $$E_i(\mathbf {x},\omega ) = \rho \omega ^2\left\langle \big | u_i(\mathbf {x},\omega )\big |^2 \right\rangle \propto -\omega ^D\mbox{Im} \big [ \mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega ) \big ]$$ (2)where ω = 2πf is the angular frequency, ui(x, ω) is the displacement field in the i direction at a point x. Im[ ] stands for the imaginary part operator, $$\mathcal {G}_{ii}(\mathbf {x},\mathbf {x},\omega )$$ is the displacement GF 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 exponent D is equal to 2 in 2-D and to 1 in 3-D (Sánchez-Sesma & Campillo 2006). The quantity $$\big | u_i(\mathbf {x},\omega )\big |^2 = u_i(\mathbf {x},\omega )u_i^*(\mathbf {x},\omega )$$, with * the conjugate operator, corresponds in frequency domain to the autocorrelation of the displacement field in time domain. In what follows, the dependence in ω and x is implicit. The ∝ means that the expressions are proportional by a factor C that is independent from ω and x. The brackets 〈 〉 represent the average operator. Theoretically, it is realized over all the possible modes of a equipartition field. In fact, in the context of the DFA, this equation is only verified when the seismic wave field is equipartitioned, that is, all the incident waves of a same state (e.g. P, S and surface waves) have the same energies (Perton et al. 2016). However, there is no way to experimentally verify the equipartition hypothesis since we only have access to the whole field, where all the states and incident, transmitted, reflected and diffracted are indiscernible. Then the average is practically realized on sliced time windows that attempt to isolate each contribution associated to the incidence of a single mode. ρ is the mass density of a medium very far away where are originated the incident waves. Its value, as well as the value of the proportional factor C may not be easily fixed with real data but they can be determined for the theoretical cases (Perton et al. 2016). 4.1 Experimental HV Although equipartition is an assumption difficult to verify and surely not true in the majority of ambient seismic field measurements, several signal processing’s exist to reinforce equipartition of the seismic wavefield (e.g. Bensen et al. 2007). Here, we only apply the operation of spectral whitening as in Spica et al. (2015). Because several sources can act in different frequency bands, the whitening consists in normalizing the signals by the source energies computed in each time window (i.e. source deconvolution) and in several frequency bands as $$\tilde{v}_i(\mathbf {x},\omega ) = v_i(\mathbf {x},\omega )/\sqrt{\sum _{i = 1}^3|v_i(\mathbf {x},\Delta \omega )|^2}$$. Δω is a frequency band of about 2.5 Hz width and centred on ω. This bandwidth is larger than the width of the peaks in Fig. 3 (since they are related to the GF) and narrow enough to remove the spectral envelope due to the seismic sources. We use the particle velocity vj(x, ω) = iωuj(x, ω) to keep the link with the energy. Experimentally, the H/V ratio is computed in terms of the wavefield autocorrelations as:   \begin{eqnarray} \frac{H}{V}^{\text{3-D}} &=& \sqrt{\frac{\langle \big | \tilde{v}_1\big |^2 \rangle +\langle \big | \tilde{v}_2\big |^2 \rangle }{\langle \big | \tilde{v}_3\big |^2 \rangle }}. \end{eqnarray} (3)The average is computed on the spectra obtained over several (≈80) time windows of 20 s duration with an overlap of 75 per cent. By denoting the square root horizontal energy in a single time window as $$\tilde{H}_W = \sqrt{ (|\tilde{v}_1(\mathbf {x}, \omega )|^2 + |\tilde{v}_2(\mathbf {x}, \omega )|^2)}$$ and the square root vertical energy as $$\tilde{V}_W = \sqrt{|\tilde{v}_3(\mathbf {x}, \omega )|^2}$$, eq. (3) is then related to the ratio of the averages: $$\sqrt{{\langle \tilde{H}_W\rangle }/ { \langle \tilde{V}_W\rangle }}$$. In that sense, it is different from the usual calculation of the HVSR (Nakamura 1989) that is related to the average of the ratios: $$\sqrt{\langle { H_W}/{V_W}\rangle }$$ (without spectral whitening). Evaluating the ratio first and then the average (e.g. Nakamura 1989) allows the source deconvolution in each time window; however, it removes the link with the GF, which is shown here to be of interest. Furthermore, the deconvolution has the disadvantage to create instabilities that are generally overcome by applying a smoothing on VW (e.g. Konno & Ohmachi 1998). Besides, in our approach, the effect of the source deconvolution is averaged over many windows and is statistically more stable. We present in Fig. 4 the velocity-averaged autocorrelations in the three directions of a local polar frame; where er and eθ directions are the radial and tangential directions in the horizontal plane as shown in Fig. 1(b). Because of the presence of the crater, the tangent and perpendicular components differ. By computing the autocorrelations along all the directions of the horizontal plane by steps of 5°, we observed that the autocorrelation in the eθ direction corresponds almost to the minimum of the energy and that autocorrelation in the er direction is almost maximum. In fact, the minimum is 20° apart from the tangent direction and its value is only 6 per cent lower. Additionally, we add the upper and lower bounds of the autocorrelations (Fig. 4). They are the maximum and minimum values at each frequency from all the autocorrelations computed with half of the windows. These results justify that the 20 min duration of recordings allows the convergence of the DED. Since the velocities autocorrelation are proportional to the ImGF by the same unknown factor −CωD/ρ (see eq. 2), we present them normalized with respect to the maximum of the three curves. This allows comparing relatively the autocorrelation with the ImGF, and determining the factor D according to the observed frequency trend (see Section 4.2); that is, determining the space dimension in which the field is diffused. This identification of D is impossible when observing only the HVSR because the factor of proportionality simplifies in the HVSR (see eq. 3). Figure 4. View largeDownload slide Autocorrelations at station s6 for the tangential (dotted black), radial (dashed grey) and vertical (continuous black) velocity components with their respective upper and lower bounds calculated with half the windows number (thin curves). Figure 4. View largeDownload slide Autocorrelations at station s6 for the tangential (dotted black), radial (dashed grey) and vertical (continuous black) velocity components with their respective upper and lower bounds calculated with half the windows number (thin curves). 4.2 Theoretical HV According to eq. (2), the H/V ratio can be computed in terms of the ImGF (Sánchez-Sesma et al. 2011a):   $$\frac{H}{V} = \sqrt{\frac{\mbox{Im}(\mathcal {G}_{11}+\mathcal {G}_{22})}{\mbox{Im}(\mathcal {G}_{33})}}.$$ (4)In this expression, the ImGF components are associated to a specific geometry and to material properties that explain adequately the experimental data. This configuration is detailed in what follows. 4.2.1 Identification of a simplified geometry for the forward problem The soil properties assessment requires simulating the propagation of elastic waves in a geometrical configuration representative of the real one; that is, by considering irregular layers, the presence of the crater wall and the heterogeneities inside the layers. However, such complex simulation (i.e. computationally achievable with heavy finite-element method, FEM) would involve too many parameters and would take too much time to provide a valid solution. Therefore, several simplifications are proposed, following four important assumptions. First, the heterogeneities inside the layers are discarded supposing that the local wavelengths are larger than the typical size of the heterogeneities by choosing an appropriate frequency band. This assumption is justified according to the assessed velocities (Section 5.2) and according to the size of the heterogeneities that are observable in Fig. 2. Second, we assume that the layers are all homogeneous, elastic, isotrope and horizontally flat. This assumption is made locally, according to the small variation of the thicknesses in horizontal directions in function of the wavelength (Fig. 2). The top surface is a free surface and the last layer is the half-space. Third, we investigate the effect of the wall over the ambient seismic source directionality by conducting numerical simulations in 3-D and 2-D and we observe that the 2-D model better represents the observed data. Fourth, we show from a theoretical approach that the HVSR and H obtained for a 2-D unbounded layered medium (i.e. without the crater wall) and when introducing a correction factor, reproduce approximately the HVSR and the H obtained when considering the complex geometry with the crater. In order to validate the third assumption, we consider a multilayered medium with the presence of a crater, respectively in 3-D and 2-D. The crater is (axi-)symmetric and surrounded by a multilayer medium made of four layers (with respective thicknesses of 18, 10, 9 and 54 m) on top of a half-space. The associated velocities for the longitudinal waves are α = 443, 6350, 647, 2200 and 4000 m s−1 and for the shear waves β = 213, 3450, 353, 1060 and 2230 m s−1. These values result from average wave velocities in rocks (Table 1) described in Fig. 2. Since they are also close to the values we obtained after inversion (Table 1), the simulations are pertinent here. The crater wall has a slope of about 30°, as for Ijen. The three imaginary part of the Gii components are calculated for superimposed source and receiver located at d = 15 m from the crater wall, as in the real case. A section of the geometry is presented in Fig. 5. Figure 5. View largeDownload slide Axisymmetric section of the geometry considered for the modelization of the H/V near a crater wall (case a). Figure 5. View largeDownload slide Axisymmetric section of the geometry considered for the modelization of the H/V near a crater wall (case a). Table 1. Comparison between average observed shear wave velocities in rocks and inverted shear velocities at Ijen. The range of velocities in the second column is based on Christensen (1982), Christensen et al. (1980), Christensen & Stanley (2003), Mavko (n.d.), Mavko et al. (2009) and Jolly et al. (2012). Layer number (from top)  Rock types  Average shear velocity (m s−1)  Inverted shear velocity (m s−1)  1 and 3  Aifall/scoria  100–600  200–400  2 and 5  Lava flows  2800–3500  2700–3400  3  Phreatomagmatic  300–800  850  6  Unknown/altered rock  .  200–400  Layer number (from top)  Rock types  Average shear velocity (m s−1)  Inverted shear velocity (m s−1)  1 and 3  Aifall/scoria  100–600  200–400  2 and 5  Lava flows  2800–3500  2700–3400  3  Phreatomagmatic  300–800  850  6  Unknown/altered rock  .  200–400  View Large The computation of the ImGF as of the HVSR is performed using the recent formulation of the indirect boundary element method (IBEM, Perton et al. 2016) that allows simulating effectively the wave propagation in piecewise irregular regions and unbounded layered media. This formulation discretizes only the crater where the virtual sources of the IBEM formulation emit directly in the unbounded stratified media. This limits drastically in 3-D the size of the model and avoids the presence of spurious wave reflections that may occur at the edges of the geometrical model in FEM or traditional IBEM or BEM. These reflections, even attenuated by appropriate absorbing boundaries may slow down the computation or affect strongly the spectra and would compromise the computation of the HVSR. In addition, to avoid memory saturation at high frequency, the radius of the crater is reduced to 250 m comparing to real configuration. 3-D configuration: the HVSR (eq. 4), as well as the square root energy densities $$H^{\text{3-D}} = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{11}+\mathcal {G}_{22})}$$ and $$V^{\text{3-D}} = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{33})}$$ are shown as continuous lines in the Fig. 6. Additionally, we superimposed the HVSR and the H obtained when the square root horizontal energy density is restrained to the separate contributions in the radial $$H^{\text{3-D}}_r = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{rr})}$$ (dotted lines) and tangential $$H^{\text{3-D}}_{\theta } = \sqrt{-\omega \,\text{Im}(\mathcal {G}_{\theta \theta })}$$ (dashed lines) directions. Figure 6. View largeDownload slide HVSR, H and V curves for case a at d = 15 m for the 3-D geometry. Full 3-D is in continuous black lines, and curves when horizontal energy is restrained to eθ and er are shown, respectively, in dashed and dot black lines. H and V curves are normalized with respect to the maximum of V, since the absolute amplitude is of no importance here. Figure 6. View largeDownload slide HVSR, H and V curves for case a at d = 15 m for the 3-D geometry. Full 3-D is in continuous black lines, and curves when horizontal energy is restrained to eθ and er are shown, respectively, in dashed and dot black lines. H and V curves are normalized with respect to the maximum of V, since the absolute amplitude is of no importance here. The H/V presents some usual characteristics: a large peak followed by a minimum that are well described in Piña-Flores et al. (2016). Generally, little attention is paid to the independent H and V curves or even to the two H components $$H^{\text{3-D}}_r$$ and $$H^{\text{3-D}}_{\theta }$$, and here they disagree strongly with the observations. The amplitude of the main peak of $$H^{\text{3-D}}_{\theta }/V$$ is about 1.5 times the one of $$H^{\text{3-D}}_r/V$$, in contradiction with the autocorrelation amplitude shown in Fig. 4, where $$\langle | \tilde{v}_r |^2 \rangle$$ is generally much larger than $$\langle | \tilde{v}_{\theta } |^2 \rangle$$. Furthermore, the curves H and V show a trend proportional to the frequency, rather the curves in Fig. 4 have a trend almost flat. Considering that the wind is the main source of noise in that frequency band, we reasonably assume that it corresponds to a coherent source all along the wall, mainly applying a force in the direction er over a large surface. Comparatively, the forces along ez and eθ should be much lower on the wall due to symmetry except at the top surface for Fz where symmetry breaks. This assumption is supported by the results shown in Fig. 4 and by the fact that it exists a direction that minimizes strongly one of the horizontal energy. Another argument comes from the wind direction measurements realized in 2010. A WindSonic anemometer connected to a Campbell CR200 data logger was programmed to record the mean and maximum wind speeds and mean speed direction. The station was installed close to the dam, which is a narrow corridor of wind (Fig. 1) and the only entry of wind inside the crater. In average, the wind speed was about 2 m s−1 with an angle of 72$${^{\circ}_{.}}$$6 from north with a standard deviation of 51$${^{\circ}_{.}}$$8. Supposing that the wind kept its direction inside the crater until reaching the stations, the stations were almost aligned perpendicularly to the wind direction. More details about the wind measurements can be find in Caudron (2013). We then conclude that the 3-D configuration is not suitable for the interpretation of the measurements. 2-D configuration: we consider now the wave propagation in the 2-D plane (er, ez). The square root energy densities in 2-D and 3-D configurations differ of a factor ω (eq. 2), so that they are now evaluated as $$H^{\text{2-D}} = \sqrt{-\omega ^2 \,\text{Im}(\mathcal {G}_{rr})}$$ and $$V^{\text{2-D}} = \sqrt{-\omega ^2 \,\text{Im}(\mathcal {G}_{33})}$$. Note that $$\mathcal {G}_{rr}$$ is a 2-D GF and that the index r indicates only the radial direction at the source–receiver position. The HVSR simplifies in $${H/V}^{\text{2-D}} = \sqrt{{\mbox{Im}(\mathcal {G}_{rr})}/{\mbox{Im}(\mathcal {G}_{33})}}$$. These quantities are shown as black lines in the Fig. 7 for two distances d = 5 m (dashed) and d = 15 m (continuous). H and V are normalized with respect to their respective maximum because the factor C is unknown. The H/V2-D differs from its 3-D counterpart, with a step around 2 Hz. Comparing to the 3-D case, the H and V present sharper peaks with maxima located at different frequencies and a different frequency trend, which agrees well with the velocity autocorrelation (Fig. 4). These differences are due to different diffraction regimes in 2-D and 3-D that modify the bulk to surface wave energy ratio (Perton et al. 2009, 2016). This configuration is then considered more adapted to represent the data and confirm the 2-D characteristic of the diffuse field. Figure 7. View largeDownload slide Comparison of the HVSR, H and V curves for case a at d = 5 m (continuous black line) and d = 15 m (dashed black line) and for case b (red lines) in 2-D. H and V have been normalized to their respective maximum. Figure 7. View largeDownload slide Comparison of the HVSR, H and V curves for case a at d = 5 m (continuous black line) and d = 15 m (dashed black line) and for case b (red lines) in 2-D. H and V have been normalized to their respective maximum. As discussed previously, the wave propagation simulation in geometries involving lateral heterogeneities, even using the IBEM method in a 2-D configuration, is not conceivable for an inverse problem. Given the positions of the source and receiver and the frequency range, we observed on simulation that the waves of interest for the HVSR only undergo reflections on one side of the crater wall and do not even probe the horizontal floor. In order to validate the fourth assumption, we consider the simulation of wave propagation in a multilayered medium only; that is, without the crater, and we add a correction factor in order to account from the wave reflection on a vertical wall (we neglect the slope of the crater wall). This configuration is denoted as case b in opposition of the previous configuration denoted as case a. The calculation for case b is realized using the DWN method (Bouchon 2003). This method is described in detail for the 3-D case in Sánchez-Sesma et al. (2011a) and for the 2-D case in the Section 4.2.2. The radiation pattern, as the sensitivity to reflected waves is different for H and V components. For H, that is, when the force is horizontal and when the displacement measured is also horizontal, the measured reflected waves are mainly longitudinal waves. Since the distance d is small with respect to the wavelength in case a, we can assume almost plane wave behaviour. In this context, the wave reflection imposes that the displacements in case a are twice the one in case b, when d = 0. However, for d ≠ 0, the incident and reflected waves are no more superimposed. We use then a correction factor of 2 cos (kd); with k = α/ω being the wave number associated to the longitudinal wave of velocity α. When introducing this factor to the H spectra in case b, the HVSR and H for d = 5 and 15 m are almost identical to their respective curve in case a. No correction was found for the V component, mainly because in this case the wall produces the interference of incident and reflected surface waves that cannot be approximated with a simple factor. Discrepancies only occur for V at d = 15 m above 5.5 Hz and do not affect the HVSR. The 2-D unbounded model is then suitable to reproduce the projected experimental data and is used in what follows as the forward modeling. 4.2.2 Forward problem method Several methods have already been proposed for the particular modelization of HVSR in multilayers and under the DFA. The DWN method (Bouchon 2003) was first proposed in Sánchez-Sesma et al. (2011b). Sophisticated methods such the one proposed in García-Jerez et al. (2016) or Perton & Sánchez-Sesma (2016) allow an efficient calculation of wave propagation in layered medium by being able to differentiate the integration of bulk and surface wave. However, the DWN method being popular, we show here that it is as efficient as the aforementioned methods in the special case when source and receiver coincide. We only summarize the method in the 2-D case with in-plane displacement, since the method was already presented for the 3-D case in Sánchez-Sesma et al. (2011b). In this method, the wavefield is decomposed as a sum of plane waves. The decomposition is realized according to the horizontal component of the wave vector kx because this component is conserved by reflection and transmission. A compact form of the displacement component $$u_i^n (\mathbf {x},\mathbf {x}_s )$$ inside the layer n, at receiver located in x and for a point source in xs and oriented in direction j is:   \begin{eqnarray} u_i^n(\mathbf {x},\mathbf {x}_s)& = &\delta (z,z_s)G_{ij}(\mathbf {x},\mathbf {x}_s) + \int _{-\infty }^{+\infty } \left( \begin{array}{c}S^{\nearrow Pn} u_i^{\nearrow Pn} e^{ik_z^{Pn} z}+ S^{\nearrow SVn} u_i^{\nearrow SVn} e^{ik_z^{SVn} z}+\\ S^{\searrow Pn} u_i^{\searrow Pn} e^{-ik_z^{Pn} z} +S^{\searrow SVn} u_i^{\searrow SVn} e^{-ik_z^{SVn} z} \end{array}\right) e^{-ik_x (x-x_s)}dk_x \end{eqnarray} (5)For an extended description of the integrand, the interested reader can consult the textbook of Aki & Richards (2002). Gij(x, xs) is the incident field as in an unbounded homogeneous medium that has to be taken into account only if the source and the receiver belong to the same layer (δ(z, zs)). $$k_z^{Pn} = \sqrt{{k^{Pn}}^2-k_x^2}$$ (respectively, $$k_z^{SVn}$$) is the vertical component of the wavenumber kPn = ω/αn (respectively, kSVn = ω/βn) associated with the P (respectively, SV) wave of the layer n with a negative sign of its imaginary part. uPn and uSVn are unitary vector of polarization of the waves. Here αn and βn are the longitudinal and shear velocities of the layer n. Then, the wave amplitudes S↗Pn, S$${\searrow}$$Pn, S↗SVn, S$${\searrow}$$SVn are solved for each discretized values of kx by considering the boundary equations: displacement and traction continuities. This equation simplifies when source and receiver coincide at the top of the layered media as:   $$u_i^1 (\mathbf {0},\mathbf {0}) = G_{ij} (\mathbf {0},\mathbf {0})+\int _{-\infty }^{+\infty } \left(S^{\nearrow P} u_i^{\nearrow P}+S^{\searrow P} u_i^{\searrow P}+S^{\nearrow SV} u_i^{\nearrow SV}+S^{\searrow SV} u_i^{\searrow SV}\right) dk_x$$ (6)The integrand presents poles associated to the surface waves. These pole contributions are usually very complicated to evaluate numerically and require using a fine discretization of kx and to replace the real frequency by a complex one: ω = ωR + iωI (Bouchon 2003). Real and imaginary parts of the integrand of G33(0, 0; ω = 6π) are presented in Fig. 8. Figure 8. View largeDownload slide Real and imaginary parts of the integrand amplitude of the component Green’s function G33(0, 0; ω) in function of horizontal wavenumber (expressed in km−1). The amplitudes are compressed using the formulae indicated in the legend. This modification preserved the sign of the integrand. Figure 8. View largeDownload slide Real and imaginary parts of the integrand amplitude of the component Green’s function G33(0, 0; ω) in function of horizontal wavenumber (expressed in km−1). The amplitudes are compressed using the formulae indicated in the legend. This modification preserved the sign of the integrand. The geometry considered is a layer (α1 = 2000 m s−1, β1 = 1000 m s−1, ρ1 = 1000 kg m−3 and h1 = 1 km) on top of a half-space (αHS = 4000 m s−1, βHS = 2000 m s−1 and ρHS = 1000 kg m−3). In order to better visualize the integrand, we use the logarithm function of the amplitude and a compression factor. However, the sign of the amplitudes is preserved. According to eq. (2), the ImGF has to be negative in sign. Here, the direct source contribution, that is, Gij(0, 0), is not considered and only the integrand is presented. Therefore, no restrictions on the sign of the integrand are assumed. The contribution of the bulk waves to the GF is delimited by the value kx = ω/β1 ≈ 18.85 km−1 and corresponds to the smooth amplitude variation below the pole peaks. Symmetries allow simplifying the calculation for positive value of kx only. Contrary to the imaginary part, the real part presents a change of sign around the poles. When evaluating numerically the continuous integral in eq. (6) by a discrete sum, this sudden change of sign makes the numerical integration unstable. Also, the integration limit has to be very large for the real part but not for the imaginary part. Same properties are observed for G11(0, 0; ω). Besides, the non-diagonal components of the GF show the opposite behaviour: the real part of the integrand are easy to integrate and not the imaginary parts. The evaluation of the imaginary part is then much easier than of the real part. This is an important observation for the HVSR and ImGF calculation and as far as we observe, it is not dependent to the layer properties. This is only the case when source and receiver coincide. Otherwise, the term $$e^{-ik_x (x-x_s ) }$$ in eq. (5) makes the previous real and imaginary parts mixed, so that difficulties encountered for real part integration also arise for the imaginary part integration. With an appropriate choice of ωI, the discretization over kx can be coarse and the integral converges quickly. The consequences of the introduction of ωI in the integration is corrected by calculating the time response through inverse Fourier transform multiplied by $$e^{-\omega _I t}$$ and then by calculating back the associated spectrum through direct Fourier transforms. Therefore, the DWN method offers a fast way to compute the HVSR following eq. (4). 5 DATA INVERSIONS 5.1 HVSR and H inversion The HVSR inversion has recently been extensively discussed (Lontsi et al. 2015; García-Jerez et al. 2016; Piña-Flores et al. 2016; Sánchez-Sesma 2017). As shown in Piña-Flores et al. (2016), a single consideration of the HVSR is insufficient for characterizing soil properties since a scale factor for velocities and thicknesses leads to the same HVSR. The additional consideration of dispersion curves helps to reduce the number of solutions by constraining the velocity model but they are not assessable with our data set. We rather explored the joint inversion of four HVSR at nearby positions, with different layer thicknesses and by considering the same elastic properties inside the layers. The thicknesses are not known with precision but they can be evaluated from the outcrop in Fig. 2. Although the thicknesses may slightly vary from the outcrop because the seismic stations are 15 m away from the crater wall, it represents an important a priori information to reduce the range of thicknesses in our inversion. The joint inversion that considers same velocities for different thicknesses avoids the multiple solutions related to the scale factor. Since the HVSR computation under DFA theory also implies the separated computation of H and V components (see Sections 4.1 and 4.2.1), they also can be used during inversion. As shown for the 3-D case in Piña-Flores et al. (2016) and Sánchez-Sesma et al. (2017), the GFs or the H and V associated to a layered media display large oscillations at high frequency related to the body waves over a trend proportional to frequency. Since the oscillations amplitudes are small regarding to the trend, they appear flatten in the HVSR, and the information that can be retrieved is mainly due to the surface waves at low frequency; however surface waves are not linearly sensitive to the thicknesses and to the body wave velocities. Then, the HVSR inversion leads to non-unique solutions. In the opposite, the use of H and V during the inversion may improve the sensitivity directly to bulk wave velocities, diminishing the number of solutions and may also be able to better characterize the shallow structure due to information collected at high frequency. This is also true in 2-D; however, since the frequency trend in GF is flat [see eq. 2 and Sánchez-Sesma et al. (2017)], the oscillations related to bulk wave are more pronounced in the HVSR and HVSR are then sensitive to bulk wave velocities as well. Although the amplitude retrieval is still a controversial topic in ambient seismic field correlation (e.g. Prieto et al. 2011; Weaver et al. 2011), several studies confirmed the good retrieval of the phase through the retrieval of surface or bulk wave velocities (e.g. Shapiro et al. 2005; Roux et al. 2005). This proves at least that the amplitude associated to individual phases might be retrieved correctly in narrow frequency bands. On larger frequency band, they might be modulated by a function, related to source envelopes, with variations slower than those of the GF (as discussed in Section 4.1). The results presented in Piña-Flores et al. (2016), such as the amplitude inversion of the HVSR, suggest that the ratio of the ImGF reduces the influence of these envelopes. It might be because of the division of the potentially same envelopes or because of the special configuration of the superimposed source and receiver that reduce the sensitivity to non-isotropic illumination or to attenuation. Here, we assume that the spectral whitening (using narrow frequency band as discussed in Section 4.1), removes any envelope. Consequently, the inversion is realized by considering jointly the H component and the HVSR. The V component is discarded because it has been shown in Fig. 7 that the one obtained in the simplified unbounded layered media might be different at high frequency to the one obtained in the configuration that includes the wall; however, if we used the more sophisticated model that integrates the crater wall for the forward model, the V should also be considered in the inversion. Then, the soil properties and the layer thicknesses are assessed by inverting jointly the HVSR and the H at several positions. The objective function (ε) is defined as the root mean square difference between the experimental HVSR and the H spectra obtained by considering only the horizontal displacement in the radial direction and the HVSR and the H spectra computed with the DWN for the 2-D case:   $$\epsilon = \sqrt{\frac{1}{N_x N_{\omega } }\sum _{x_i}\sum _{\omega } \big ( \big(H_r^{\text{exp}}\big/V^{\text{exp}}-H^{2\text{-D}}\big/V^{\text{2-D}}\big)^2+\big(H_r^{\text{exp}}-H^{\text{2-D}}\big)^2 \big ) }.$$ (7)In this expression, the theoretical and experimental HVSR spectra are normalized according to the maximum amplitude in the experimental spectra in order to balance the error weight with the H spectra, which are also normalized as discussed in Section 4.1 since the constant C is not known. Nx = 4 is the number of positions xi considered in the joint inversion and Nω = 30 is the number of points taken linearly separated in the spectra from 1.5 to 6 Hz. The lower frequency bound is taken according to the results observed in the HVSR in Fig. 7. Curves at higher frequency are associated to small and local heterogeneities, which are of minor importance here. Rather than using the simulated annealing or genetic algorithm methods discussed in García-Jerez et al. (2016) and Piña-Flores et al. (2016), we use a ‘Pattern-Search’ method (Audet & Dennis 2002) because of its efficiency for highly non-linear problem and its ability to integrate constraints. Beside the algorithm has the drawback to fail to converge to the global minimum when the solution space is very large, the non-uniqueness of the solution is strongly reduced thanks to the large number of constraints considered. The elastic velocities β and the layer thicknesses at each position are the free parameters. The thickness of each layer is allowed to vary by 20 per cent only from its original a priori value. No variation restriction is considered for the shear wave velocities but they are taken constant all along the geological layers for the different positions. During the first inversions, α was also taken as a free parameter but results showed that it might be equally taken as α = 1.9β for all the layers. This value is consistent with values reported in the references given in Table 1. This simplification allowed reducing the number of free parameters and obtaining better fit. The density is related to α through the linear relationship: ρ = 0.32α + 0.77 with α expressed in km s−1 and ρ in kg m−3 (Berteussen 1977). 5.2 Results and discussion Results of the inversion for the subset of stations s1, s5, s6 and s8 are shown in Fig. 9. The H and V spectra are shown normalized. An important feature of these results is that both the amplitude and the entire shape of the HVSR and H curves (i.e. not only the main peak) are well retrieved for all the stations and for frequencies ranging from 1.5 to 6 Hz. As expected for the V curves, the frequency peak positions agree well until 4 Hz, but clear discrepancies are visible at higher frequencies. For the HVSR curves, the agreement is still acceptable for higher frequencies, but it is due to the absence of local peak. For the H curves, a peak around 4 Hz differs from the experimental ones, probably meaning that a superficial and thin layer might be missing in the model. The top layers could then be refined according to the high number of sublayers observed in the outcrop but their thickness evaluation is uncertain. Figure 9. View largeDownload slide Results of the joint HVSR inversion for the subset of stations s1, s5, s6 and s8 (shown in Fig. 2). Black and red dotted lines depict the target and best fit (ε0) modeled HVSR, H and V curves, respectively. The velocity models associated to the best fit are plotted with thick dotted red lines in the last column. Other profiles with higher misfit are also shown in continuous thin lines. The scale is from blue (≈ε0) to red lines (≈1.5ε0). The mean lake level is depicted as the cyan area. Figure 9. View largeDownload slide Results of the joint HVSR inversion for the subset of stations s1, s5, s6 and s8 (shown in Fig. 2). Black and red dotted lines depict the target and best fit (ε0) modeled HVSR, H and V curves, respectively. The velocity models associated to the best fit are plotted with thick dotted red lines in the last column. Other profiles with higher misfit are also shown in continuous thin lines. The scale is from blue (≈ε0) to red lines (≈1.5ε0). The mean lake level is depicted as the cyan area. The models associated to the best fit are superimposed in Fig. 2 and we can observe that they agree well with the local geology observed in the outcrop. As discussed in Spica et al. (2015), the outcrop of the north wall is composed by the following succession of layers (from top to bottom): (1) air-fall deposits and scorias covered by a veneer of consolidated sulphur-bearing mud; (2) interbedded thin lava flow; (3) unconsolidated phreatomagmatic deposits and (4) thick lava flows. These geological units are highlighted by numbered points in Fig 2. All these geological units are well represented in the velocity model and their thicknesses change in accordance with the outcrop characteristics. As an example, a strong impedance contrast is observed at the interface between the phreatic/superficial deposits (850 m s−1) lying on the thick lava flows (3400 m s−1; layer 4 over layer 5). The thin lava flow with an important velocity of 2700 m s−1 is particularly well observed in the models and in the outcrop. Although this layer is particularly thin, it appears remarkably well on the velocity model, suggesting the model is well assessed. The retrieved velocities are comparable to those obtained in a previous study of the Ijen crater (Spica et al. 2015) and charts of seismic velocities for these geological units (Table 1). Interestingly the half-space, referred by (5) in Fig. 2, corresponds to a layer in contact with the lake with a velocity lower than the above thick lava flow referred by the number (4). A possible interpretation is that it corresponds to the same geological unit but that has been strongly altered or corroded by the acidity of the lake making it more porous. Other profiles with higher misfit error are also depicted in Fig. 9. The four first layers have their thicknesses and velocities retrieved with confidence. The deepest layers are retrieved with large incertitudes, as expected, since they depend strongly to low frequency and then to the surface waves. These results validate the 2-D character of the ambient seismic source, the appropriate data processing and the importance of considering independently the H and V components. Here, the V component has been discarded only because of the presence of the strong lateral heterogeneity but might be considered elsewhere. The H and V components show important peak at higher frequency than the HVSR and allow then a better identification of the shallow layer properties. 6 CONCLUSIONS We explored here the effects of a strong lateral heterogeneity on the HVSR by using measurements along a vertical crater wall. We identified that the horizontal square root energy density H restrained to the direction normal to the crater wall is very close to the H in a 2-D layered medium. This is also the case for the V component but in a more limited frequency bandwidth. This identification is realized from theoretical simulations using IBEM when considering the crater embedded in a layered medium and the DWN method when considering only the same layered medium and with similar elastic properties as the one obtained after inversion. We demonstrated that the DWN method is efficient for the calculation of the ImGF when source and receiver are superimposed and may then be used for the forward modeling. The inversion was realized jointly by considering the HVSR curves and the independent H at several positions and with initial thicknesses roughly evaluated. The V energy density might also be used in other circumstances. The experimental H and V have been processed according to the theoretical assumption of a diffuse field, that is, by applying a spectral whitening and by conducting the mean autocorrelation as for GF retrieval. The good agreement of the results with the local geology validates the assumptions made such as the 2-D character of the ambient seismic source and the appropriate data processing. This paper also underlines that the H and V components could also be considered independently to the HVSR. It demonstrates then the potential of the method in complex geometrical configuration. This methodology paves the way for identifying layer elastic velocities at outcrop and then to follow their thickness variation far inside. Acknowledgements We warmly thank Thierry Camelbeeck and Thomas Lecocq from the Royal Observatory of Belgium for their contribution in the early stage of this research and the development of the Kawah Ijen monitoring project. The field data acquisition was funded by the Royal Observatory of Belgium. We also thank Cipta Anthanasius, Sarane Sterckx and Arnaud Watlet for their help on the field and data acquisition. We are grateful to the Indonesian Centre for Volcanology and Geological Hazard Mitigation support on the field, and particularly to the observers of Kawah Ijen observatory: Bambang Heri Purwanto and Suparjan. We appreciate the valuable suggestions of Thomas Bodin, Agostiny Marrios Lontsi and two anonymous reviewers. REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , University Science Book. Albarello D., Lunedei E., 2011. Structure of an ambient vibration wavefield in the frequency range of engineering interest ([0.5, 20] hz): insights from numerical modelling, Near Surf. Geophys. , 9( 6), 543– 559. Arai H., Tokimatsu K., 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum, Bull. seism. Soc. Am. , 94( 1), 53– 63. Google Scholar CrossRef Search ADS   Audet C., Dennis J.E., Jr, 2002. Analysis of generalized pattern searches, SIAM J. Optim. , 13( 3), 889– 903. Google Scholar CrossRef Search ADS   Bard P.Y., 1999. Microtremor measurements: a tool for site effect estimation, in The Effects of Surface Geology on Seismic Motion , 3, pp. 1251– 1279, Balkema, Rotterdam. Bard P.-Y., 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: ftp://delta.geo.uib.no/pub/seismo/REPORTS/SESAME/USER-GUIDELINES/SESAME-HV-User-Guidelines.doc, last accessed 1st November 2017. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. Google Scholar CrossRef Search ADS   Berteussen K.A., 1977. Moho depth determinations based on spectral-ratio analysis of NORSAR long-period P waves, Phys. Earth planet. Inter. , 15( 1), 13– 27. Google Scholar CrossRef Search ADS   Bonnefoy-Claudet S., Köhler A., Cornou C., Wathelet M., Bard P.-Y., 2008. Effects of Love waves on microtremor h/v ratio, Bull. seism. Soc. Am. , 98( 1), 288– 300. Google Scholar CrossRef Search ADS   Bouchon M., 2003. A review of the discrete wavenumber method, Pure appl. Geophys. , 160( 3–4), 445– 465. Google Scholar CrossRef Search ADS   Caudron C., 2013. Multi-disciplinary continuous monitoring of Kawah Ijen volcano, East Java, Indonesia, PhD thesis , Université Libre de Bruxelles, Brussels, Belgium. Caudron C., Lecocq T., Syahbana D.K., McCausland W., Watlet A., Camelbeeck T., Bernard A., Surono S., 2015a. Stress and mass changes at a ‘wet’ volcano: example during the 2011–2012 volcanic unrest at Kawah Ijen volcano (Indonesia), J. geophys. Res. , 120( 7), 5117– 5134. Google Scholar CrossRef Search ADS   Caudron C.et al.  , 2015b. Kawah Ijen volcanic activity: a review, Bull. Volcanol. , 77( 3), 1– 39. Google Scholar CrossRef Search ADS   Caudron C., Campion R., Rouwet D., Lecocq T., Capaccioni B., Syahbana D., Purwanto B.H., Bernard A., 2017. Stratification at the Earth’s largest hyperacidic lake and its consequences, Earth planet. Sci. Lett. , 459, 28– 35. Google Scholar CrossRef Search ADS   Christensen N.I., 1982. Seismic velocities, in Handbook of Physical Properties of Rocks , Vol. 2, pp. 1– 228. Christensen N.I., Stanley D., 2003. Seismic velocities and densities of rocks, Int. Geophys. , 81, 1587– 1594. Google Scholar CrossRef Search ADS   Christensen N., Wilkens R., Blair S., Carlson R., 1980. 10. Seismic velocities, densities, and elastic constants of volcanic breccias and basalt from deep sea drilling project LEG 59, in Initial Reports of the Deep Sea Drilling Project , 59, pp. 515– 517, U.S. Govt. Printing Office, Washington, DC. Delmelle P., Bernard A., 1994. Geochemistry, mineralogy, and chemical modeling of the acid crater lake of Kawah Ijen volcano, Indonesia, Geochim. Cosmo. Acta , 58( 11), 2445– 2460. Google Scholar CrossRef Search ADS   Delmelle P., Bernard A., Kusakabe M., Fischer T.P., Takano B., 2000. Geochemistry of the magmatic–hydrothermal system of Kawah Ijen volcano, east Java, Indonesia, J. Volcanol. Geotherm. Res. , 97( 1–4), 31– 53. Google Scholar CrossRef Search ADS   García-Jerez A., Piña-Flores J., Sánchez-Sesma F.J., Luzón F., Perton M., 2016. A computer code for forward calculation and inversion of the H/V spectral ratio under the diffuse field assumption, Comput. Geosci. , 97, 67– 78. Google Scholar CrossRef Search ADS   Gunawan H. et al.  , 2016. New insights into Kawah Ijen’s volcanic system from the wet volcano workshop experiment, Geol. Soc. Lond. Spec. Publ. , 437, doi:10.1144/SP437.7. Jolly A.D., Chardot L., Neuberg J., Fournier N., Scott B.J., Sherburn S., 2012. High impact mass drops from helicopter: a new active seismic source method applied in an active volcanic setting, Geophys. Res. Lett. , 39( 12), L12306, doi:10.1029/2012GL051880. Konno K., Ohmachi T., 1998. Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor, Bull. seism. Soc. Am. , 88( 1), 228– 241. Lermo J., Chávez-García F.J., 1994. Are microtremors useful in site response evaluation?, Bull. seism. Soc. Am. , 84( 5), 1350– 1364. Lontsi A.M., Sánchez-Sesma F.J., Molina-Villegas J.C., Ohrnberger M., Krüger F., 2015. Full microtremor H/V(z, f) inversion for shallow subsurface characterization, Geophys. J. Int. , 202( 1), 298– 312. Google Scholar CrossRef Search ADS   Lunedei E., Albarello D., 2010. Theoretical HVSR curves from full wavefield modelling of ambient vibrations in a weakly dissipative layered Earth, Geophys. J. Int. , 181( 2), 1093– 1108. Lunedei E., Malischewsky P., 2015. A review and some new issues on the theory of the H/V technique for ambient vibrations, in Perspectives on European Earthquake Engineering and Seismology , pp. 371– 394, Springer International Publishing. Google Scholar CrossRef Search ADS   Matsushima S., Hirokawa T., Martin F.D., Kawase H., Sánchez-Sesma F.J., 2014. The effect of lateral heterogeneity on horizontal to vertical spectral ratio of microtremors inferred from observation and synthetics, Bull. seism. Soc. Am. , 104( 1), 381– 393. Google Scholar CrossRef Search ADS   Mavko G., n.d. Conceptual Overview of Rock and Fluid Factors that Impact Seismic Velocity and Impedance , Stanford School of Earth Sciences . Available at: https://pangea.stanford.edu/courses/gp262/Notes/8.SeismicVelocity.pdf, last accessed 1st November 2017. Mavko G., Mukerji T., Dvorkin J., 2009. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media , Cambridge University Press. Google Scholar CrossRef Search ADS   Nakamura Y., 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, Q. Rep. Railw. Tech. Res. Inst. , 30( 1), 25– 33. Perton M., Contreras-Zazueta M.A., Sánchez-Sesma F.J., 2016. Indirect boundary element method to simulate elastic wave propagation in piecewise irregular and flat regions, Geophys. J. Int. , 205( 3), 1832– 1842. Google Scholar CrossRef Search ADS   Perton M., Sánchez-Sesma F.J., 2016. Green’s function calculation from equipartition theorem, J. acoust. Soc. Am. , 140( 2), 1309– 1318. Google Scholar CrossRef Search ADS PubMed  Perton M., Sánchez-Sesma F.J., Rodríguez-Castellanos A., Campillo M., Weaver R.L., 2009. Two perspectives on equipartition in diffuse elastic fields in three dimensions, J. acoust. Soc. Am. , 126( 3), 1125– 1130. Google Scholar CrossRef Search ADS PubMed  Piña-Flores J., Perton M., García-Jerez A., Carmona E., Luzón F., Molina-Villegas J.C., Sánchez-Sesma F.J., 2016. The inversion of spectral ratio H/V in a layered system using the diffuse field assumption (DFA), Geophys. J. Int.  208( 1), 577– 588. Google Scholar CrossRef Search ADS   Prieto G.A., Denolle M., Lawrence J.F., Beroza G.C., 2011. On amplitude information carried by the ambient seismic field, C. R. Geosci. , 343( 8), 600– 614. Google Scholar CrossRef Search ADS   Roux P., Sabra K.G., Gerstoft P., Kuperman W.A., Fehler M.C., 2005. P-waves from cross-correlation of seismic noise, Geophys. Res. Lett. , 32( 19), L19303, doi:10.1029/2005GL023803. Sánchez-Sesma F.J., 2017. Modeling and inversion of the microtremor H/V spectral ratio: physical basis behind the diffuse field approach, Earth Planets Space , 69, 92, doi:10.1186/s40623-017-0667-6. Sánchez-Sesma F.J., Campillo M., 2006. Retrieval of the Green’s function from cross correlation: the canonical elastic problem, Bull. seism. Soc. Am. , 96( 3), 1182– 1191. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J.et al.  , 2011a. A theory for microtremor H/V spectral ratio: application for a layered medium, Geophys. J. Int. , 186( 1), 221– 225. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Weaver R.L., Kawase H., Matsushima S., Luzón F., Campillo M., 2011b. Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid, Bull. seism. Soc. Am. , 101( 4), 1704– 1709. Google Scholar CrossRef Search ADS   Sánchez-Sesma F.J., Iturrarán-Viveros U., Perton M., 2017. Some properties of Green’s functions for diffuse field interpretation, Math. Methods Appl. Sci. , 40( 9), 3348– 3354. Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307( 5715), 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Spica Z.et al.  , 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, J. Volc. Geotherm. Res. , 302, 173– 189. Google Scholar CrossRef Search ADS   Spica Z., Perton M., Nakata N., Liu X., Beroza G., 2017. Site characterization at Groningen gas field area through joint surface-borehole H/V analysis, Geophys. J. Int. , doi:10.1093/gji/ggx426. Sujanto M., Syarifuddin M.Z., Sitorus K., 1988. Geological Map of The Ijen Caldera Complex, East Java, Volcanol. Surv. Indones . van Hinsberg V., Berlo K., van Bergen M., Williams-Jones A., 2010. Extreme alteration by hyperacidic brines at Kawah Ijen volcano, east Java, Indonesia: I. textural and mineralogical imprint, J. Volc. Geotherm. Res. , 198( 1–2), 253– 263. Google Scholar CrossRef Search ADS   Weaver R.L., Lobkis O.I., 2004. Diffuse fields in open systems and the emergence of the Green’s function (l), J. acoust. Soc. Am. , 116( 5), 2731– 2734. Google Scholar CrossRef Search ADS   Weaver R.L., Hadziioannou C., Larose E., Campillo M., 2011. On the precision of noise correlation interferometry, Geophys. J. Int. , 185( 3), 1384– 1392. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial