Site response measurements and implications to soil liquefaction potential using microtremor H/V in Greater Metro Manila, Philippines

Abstract

This research explores the use of microtremor horizontal-to-vertical spectral ratio (H/V) in obtaining site response characteristics and investigating its relationship with soil liquefaction potential in Greater Metro Manila. We performed single station microtremor measurements in 61 sites along with in situ geotechnical techniques to verify liquefaction potential. The resulting 238 spectral curves were classified according to dominant features and subsequently grouped with the calculated liquefaction potential index (LPI) of the soil. Based on a robust comparison of obtained primary parameters, it is revealed that the shape of the H/V curve, its predominant period and relative amplitude are fundamentally linked to the spatial variability and the shear strength of soils. Therefore, areas of high seismic demand can also have high liquefaction potential, and vice versa. We then correlated the predominant period with the LPI of the soil and extracted a boundary using simple statistical techniques to classify high and low potential for liquefaction subsequently validating its use as a complementary tool for rapid site-specific liquefaction assessment. Such findings are a novel contribution to liquefaction studies employing rapid techniques since the application of microtremors to liquefaction in the Philippines has not been practiced extensively.

Keywords:

• Site response
• microtremor H/V
• predominant period
• liquefaction
• LPI

Introduction

Understanding how soils respond to strong ground shaking is a significant concern in seismic hazard estimation across all geological environments, most specially in metro regions where infrastructure is continuously being developed. Various factors control the response of soils to strong earthquakes including, source proximity, near-surface soil properties, acoustic impedance contrasts and basin configuration, to name a few. The profound effects of amplified ground shaking in soft soils are historically demonstrated in the 1985 Ms 8.1 Mexico City earthquake and the 1989 Ms 7.2 Loma Prieta earthquake, where buildings sustained considerable damages despite being hundreds of kilometers away from the earthquake epicenter. In the Philippines, a comparable amplified ground shaking phenomenon was observed during the 1968 Ms 7.3 Casiguran earthquake and the 1990 Ms 7.9 Luzon earthquake where buildings in the capital city of Manila suffered severe damages despite being located about 225 and 120 kilometers away from the earthquake epicenters, respectively (Rodolfo 2014; PHIVOLCS n.d.-a, n.d.-b).

In addition to damage on buildings due to amplified ground shaking, damages associated with soil liquefaction is also a quintessential consequence of earthquakes. Liquefaction occurs when the soil material below the water table temporarily loses strength and behaves as a viscous liquid rather than a solid (Earthquake Engineering Research Institute 1994), rendering structures founded on liquefiable soil to sustain damages. Comparatively, a hallmark of both phenomena is the amplification and prolonging of the duration of seismic energy contributing to the severity of the damage at the surface (Kramer 1996; Towhata 2008; Beroya and Aydin 2010; Mase et al. 2022). Although the mechanism for seismic amplification varies per location, site response characteristics can be investigated in advance of an earthquake using microtremors or ambient noise (Nakamura 1989, 2000, 2019; Mase et al. 2023).

In recent years, global applications of microtremor horizontal-to-vertical spectral ratio (H/V) for liquefaction studies utilized the seismic vulnerability index of Nakamura (1996, 1997) (and its variations) that basically constitute the product of the site response parameters - predominant period of the soil and its relative amplitude to quantify liquefaction prediction (Huang and Tseng 2002; Hardesty et al. 2010; Choobbasti et al. 2015; Sathyaseelan et al. 2017; Singh et al. 2017; Herrera et al. 2018; Ramos et al. 2019; Meneisy et al. 2020; Arango-Serna et al. 2021; Kang et al. 2021). The relevant downside of this technique is that it is based on multiple and complex assumptions as listed in Herrera et al. (2018) and Arango-Serna et al. (2021). Moreover, some of the listed works have a lack in integrating a design earthquake model which can trigger liquefaction. If present, there is a lack in the robust comparison of the more basic parameters which are the predominant period and its relative amplitude with significant subsurface information that can contribute to liquefaction occurrence. With a large set of subsurface data, the relationship of the site response parameters with liquefaction potential can be better understood.

In this paper, the microtremor H/V was used to study potential site response on seismic ground motion in the Greater Metro Manila area (GMMA), a densely populated region that hosts an actively seismic zone. This study also presents the first application of the method in terms of liquefaction studies in the study area. We focused on sites that are potentially liquefiable by virtue of liquefaction potential index (LPI) calculated from in-situ geotechnical data acquired in conjunction with the microtremor H/V. This is to facilitate the possible implications of site response with liquefaction potential. Predominant period values, relative amplitudes, and soil thicknesses are then compared with geological and geotechnical data, particularly the soil shear-wave velocities, soil N-values, and the LPI of the soil deposit. In addition, a boundary was extracted separating high and low potentials of liquefaction with the predominant period using simple statistical analysis which can be used as a first-pass indication for liquefaction occurrence.

Geological setting and earthquake scenario

GMMA is a large region encompassing the provinces of Bulacan, Cavite and Metro Manila. It lies at the southernmost extension of the Central Luzon Basin (Figure 1) and is bounded to the west by the Manila Bay and to the east by the mountains of Sierra Madre.

Figure 1. Experimental sites and geological map of Greater Metro Manila area (GMMA), Philippines.

The Quaternary Alluvium is characterized by soft and unconsolidated sediments of dominantly sand and clay interlayers of varying thickness. The alluvium is extensively distributed along the coastal lowlands and terminates to the east as it thins out towards the central plateau. Based on the accounts of Gervasio (1968), the formation of the coastal lowlands started during the Pliocene, followed by its abrupt seaward expansion at the advent of the Pleistocene. The expansion was supplied by volcanic sediments associated with the calderagenic episodes of the Taal Caldera or Laguna de Bay. The episodes deposited thick tuffaceous rocks and bedded tuff in the east corresponding to the Guadalupe Formation. Bounding the plateau to the east is the West Valley Fault (WVF) of the Marikina Valley Fault System (MVFS). The WVF is an approximately 100-km north to south trending right lateral strike slip fault with surface manifestation that runs from the mountains of Bulacan in the north to the highlands of Tagaytay in the south.

Based on paleoseismic investigation of Nelson et al. (2000), the WVF has generated 4 surface-rupturing earthquakes over a period of <1300 years with a preferred recurrence interval of 400-600 years. Its last major earthquake was documented in 1658, based on accounts of historians and priests (Masó 1910). Considering the lower bound of the recurrence interval, the WVF has a possibility for a near-future seismic event. The Metro Manila Earthquake Impact Reduction Study (2004) and Rimando and Knuepfer (2006) estimated a magnitude 7.2 and 7.3 earthquake, respectively, should the WVF rupture. These magnitude estimates translate to an Intensity VIII in the PHIVOLCS Earthquake Intensity Scale (PEIS) (Table 1).

Table 1. PHIVOLCS Earthquake Intensity Scale (PEIS) with modified mercalli Intensity (MMI) equivalent.

PEIS	Shaking	MMI
I	Scarcely Perceptible	I
II	Slightly Felt	II
III	Weak	III
IV	Moderately Strong	IV
V	Strong	V
VI	Very Strong	VI
VII	Destructive	VII
VIII	Very Destructive	VIII, IX
IX	Devastating	X, XI
X	Completely Devastating	XII

Materials and methods

For the microtremor H/V analysis, single station microtremor measurements were employed to sample raw ambient noise data. Essentially, it estimates the amplification characteristics at a certain period of motion of the ground by only calculating the ratio of the horizontal with the vertical Fourier amplitude spectra of only surface measurements (Nakamura 1989). The technique follows the formula, (1) $$H/V=\frac{\sqrt{(N^2+E^2)/2}}{V}$$(1) where $N$ and $E$ are the Fourier amplitude spectra of the horizontal components and $V$ is the vertical component. The resulting spectral peak corresponds to the predominant period and its relative amplitude consistent to the major impedance contrast in the subsurface profile. The microtremor H/V differs from the standard spectral ratio such that transfer function is calculated by using only the recorded ground motions at a single soft soil site, eliminating the need for a reference site (Nakamura 1989; Lermo and Chávez-García 1994).

For the geotechnical analysis and liquefaction potential investigations, Refraction Microtremor (REMI) surveys and Screw Driving Sounding (SDS) tests were used. The REMI survey was developed by Louie et al. (2001) and is used in this study to obtain soil shear-wave velocity profiles of the subsurface. It is a non-invasive technique that analyzes and extracts dispersion characteristics of rayleigh waves determined from ambient noise recordings. Louie et al. (2001), Mase et al. (2018) and (Daag et al. 2022) provide details on further surface wave analysis using this technique. On the other hand, the SDS test is a relatively new penetration test developed in Japan and it is used in this study to obtain soil N-values comparable to the Standard Penetration Test (SPT) (Orense et al. 2019; Daag et al. 2023). It is essentially an improved and automated version of the Swedish Weight Sounding (SWS) test, capable of classifying sands from clays by percent fines estimation with high probability. Basically, we semi-quantitatively compared general and readily observable primary output of the methods (for e.g. predominant period and relative amplitude against soil shear-wave velocity and N-value) to facilitate in the understanding of the possible relationship between site response and liquefaction potential. Supplementary depth to groundwater table data necessary for liquefaction investigation was determined by the Ground Penetrating Radar (GPR).

Data acquisition and processing

Majority of the sites for microtremor H/V measurements were along the coastal lowlands of GMMA, where the soil deposits are susceptible to liquefaction. Additional sites were acquired along the central plateau for comparison. We acquired single station microtremor data at 61 temporary sites using the OYO McSEIS MT-NEO. It is an all-in-one (battery and sensor) triaxial accelerometer with a natural frequency response of 0.1 to 200 Hz. Sampling interval was set at 10 milliseconds. The recording duration for each station was set at 20 min, with some reaching over an hour, depending on the presence of transients. All the recordings followed the guidelines set by the SESAME Project (SESAME 2004). Data recorded on the instrument were stored in an SD card and transferred to the computer. Raw motion files of each station were converted to multi-column ASCII format to be analyzed in Geopsy, an open-source software capable of processing and visualizing multiple signals at once (Wathelet et al. 2020). Processing parameters were kept as consistent as possible with modifications only applied as recommended by the SESAME Project.

REMI survey and SDS test were also conducted in the proximity of the microtremor measurements to obtain the shear-wave velocity and soil N-values of the shallow subsurface, respectively for validation of the spectral curves. Both methods were also used for the calculation of the liquefaction potential index (LPI), together with water table depth estimates from the Ground Penetrating Radar (GPR). For the REMI survey, twelve 4.5 Hz geophones were installed on a linear array with about 4- to 8-meter spacing, depending upon site conditions. A sledgehammer was used to induce vibration and microtremor data for analysis were collected using DAQLink II seismograph and Vscope software. Data processing was carried out using REMI Vspect and Disper software to retrieve the shear wave velocity profile of the particular soil deposit. For the SDS test, it involves the penetration of a rotating rod into the soil at 7 incremental loading steps (0.25 to 1.0 kN) until reaching 25 centimeters of penetration. At every 25 centimeters of penetration, several parameters are measured such as torque, load, penetration speed, and rod friction. These parameters are then sent to the proprietary Geoweb server where it automatically evaluates the equivalent N-values. Further details on test operations and parameter estimations are in Orense et al. (2019) and (Daag et al. 2023). For the GPR survey, we employed the MALA ProEx GPR using antenna center frequency of 500 MHz to resolve the shallow water table depth. Raw data were processed in RadExplorer to obtain the radargrams.

Results and discussion

Site period distribution

Out of 61 sites, we collected 238 H/V curves. Most sites contained more than one H/V measurement and thus, Figure 2 shows the representative site period or simply the average of the predominant period of all measurements within the site.

Figure 2. Site period distribution and relative amplitude across the study area (amplitude values from 2 to 8 scales with the dot size). for sites with multiple H/V measurements, the predominant period and their relative amplitudes were simply averaged to obtain the representative values. Labelled sites are used as an example in this section and latter sections of the paper.

The sites were classified based on predominant period from the H/V analysis. Long period sites have predominant periods > 0.6 s; moderate period sites have predominant periods in the range of 0.2 s to 0.6 s; short period sites have predominant periods in the range of 0.05 s to 0.2 s; and no peak sites have no identifiable H/V peaks. The two peaks case was not observed much in this study suggesting that there is only one major impedance contrast present in the study area. Majority (43%) of the experimental sites are classified as moderate period sites and are distributed throughout the coastal lowlands. Long period sites (25%) are observed mostly in the northern coast where the environment is dominated by tidal plains and sand bars while there are a few that are near the reclaimed area in Metro Manila indicating that soil deposits are generally thicker and softer. On the other hand, short period sites (15%) and sites with no peaks (18%) are found sparsely throughout the study area but are nearing or atop the central plateau where the soft soil deposit is relatively thin. The period distribution is consistent with the results of Narag et al. (2000) specifically within the coastal area of Metro Manila. Examples of H/V curves are presented in Figure 3 showing the distinguishable H/V curve shapes across the period classifications.

Figure 3. Examples of the (a) long period sites (predominant period > 0.6 s), (b) moderate period sites (0.2 s to 0.6 s), (c) short period sites (0.05 s to 0.2 s) and (d) sites with no peaks.

The relationship between the predominant period of the soil, $T,$ and the overall soil thickness, $h,$ down to the major impedance contrast is given by the quarter wavelength formula, (2) $$T=\frac{4h}{v_s}$$(2) where $v_s$ is the local average shear-wave velocity. Using the shear-wave velocity obtained from the REMI survey as a constraint, overall soil thickness across the survey area is shown in Figure 4. It is important to note that the shear-wave velocity of the major impedance contrast in and around Metro Manila is ∼600 to 760 m/s equivalent to the engineering bedrock (Metro Manila Earthquake Impact Reduction Study 2004; Grutas and Yamanaka 2012).

Figure 4. Spatially interpolated distribution of soil thickness across the study area calculated from the predominant period and local average shear-wave velocity.

Throughout the coastal lowlands, soil thickness varies from ∼2.5 meters near the rocky coast of Cavite to 56 meters in Bulacan. A general increase in soil thickness is also observed from east to west which is typical in the study area following the gentle dipping topography of the central plateau towards Manila Bay.

In general, based on the distribution of the predominant period and amplitude as well as the overall soil thickness across the study area, the potential for amplification of seismic ground motion is expected to be greatest towards the north in Bulacan, surpassing earlier microtremor measurements of Abeki et al. (1996) and Narag et al. (2000) in Metro Manila where the period and amplitude distribution is consistent with their results. These findings in Bulacan correlate well with its topographical location being proximal to the center of the Central Luzon Basin, thereby having thicker soils as compared to Metro Manila and Cavite. Overall soil thickness decreases towards the south and east as unconsolidated sediments thin out or increase in stiffness. In Cavite, no amplification is expected along its rocky coastline, however, there are a few sites that exhibit some degree of amplification in sites along the spit feature and nearby river systems which are likely to be underlain by reworked sediments due to coastal or fluvial activity, respectively.

Liquefaction potential index (LPI) distribution

In this section, we determine the liquefaction potential distribution of the soil deposit nearby microtremor H/V measurements. Soil profiles produced from the REMI surveys and the SDS tests reveal shallow soil layers with varying shear strengths and thicknesses (Figure 5). These profiles are grouped according to period classification determined from the previous section. Groundwater table depths are also reflected in the soil profiles showing the near-surface saturation of the soil in almost all sites.

Figure 5. Example of soil profiles showing depth vs. N-value from the SDS tests and depth vs. shear-wave velocity from the REMI surveys grouped according to period classification. a1-a6 are profiles in long period sites, b1-b7 are profiles in moderate period sites, c1-c6 are profiles in short period sites, and d1-d5 are profiles in sites with no peaks.

Current practice of site-specific liquefaction assessment revolves around modifications of the ‘simplified procedure’ of Seed and Idriss (1971, 1982) using different invasive (Cetin et al. 2004; Moss et al. 2006; Boulanger and Idriss 2012, 2016) and non-invasive (Andrus and Stokoe 2000; Kayen et al. 2013) techniques. This simplified procedure involves empirical stress-based calculations where the liquefaction resistance of soil is compared to a triggering mechanism acting on the soil to induce liquefaction (Seed and Idriss 1971; Youd and Perkins 1978; Youd et al. 2001). However, in this paper, we followed Maeda et al. (2015) in estimating liquefaction resistance for SDS test data based on the ‘Specifications for Highway Bridges’ of Japan Road Association (JRA) (Japan Road Association 1996). This is to duplicate their calculations since the SDS test is a relatively new technique and it is not an objective of this paper to explore different methodologies. We also followed the procedure proposed by Kayen et al. (2013) based on the shear-wave velocity obtained from the REMI survey since their work is comprehensive and is built on the work of Andrus and Stokoe (2000). Both methods are independent in estimating liquefaction risk at depth.

Using the JRA approach, the liquefaction potential of the soil can be computed through the liquefaction resistance factor ($F_L$), which is the ratio between the dynamic shear stress ($R$) and seismic shear stress ($L$) as seen in Equation (3)(3) $$F_L=\frac{R}{L}$$(3) . (3) $$F_L=\frac{R}{L}$$(3) $F_L$ specifies the risk of liquefaction in depth, where values less than 1 will indicate the possibility of the hazard. The seismic shear stress ($L$) is calculated using Equation (4)(4) $$L={\gamma }_d{k}_{\mathit{hgL}}\frac{{\sigma }_v}{{\sigma \prime }_v}\hspace{1em}\mathrm{ }{\gamma }_d=1.0-0.015x$$(4) , (4) $$L={\gamma }_d{k}_{\mathit{hgL}}\frac{{\sigma }_v}{{\sigma \prime }_v}\hspace{1em}\mathrm{ }{\gamma }_d=1.0-0.015x$$(4) where $k_{\mathit{hgL}}$ is the peak ground acceleration in g, ${\gamma }_d$ is the reduction coefficient, $x$ is the depth from the surface in m, and ${\sigma }_v$ and ${\sigma \prime }_v$ as the total and effective overburden pressure, respectively, at depth $x.$ On the other hand, the dynamic shear stress ($R$) is calculated using Equation (5)(5) $$R=C_w{R}_L$$(5) , (5) $$R=C_w{R}_L$$(5) where $C_w$ is the coefficient of seismic motion and $R_L$ is the cyclic triaxial strength ratio calculated from the SDS N-value and fines content. Since the SDS test estimates parameters every 25 centimeters, $F_L$ can be calculated for every 25 centimeters.

Using the formula developed by Kayen et al. (2013), the potential for liquefaction can also be evaluated by the factor of safety against liquefaction ($F_S$), which is obtained by comparing the earthquake-induced seismic loading or cyclic stress ratio on a soil layer ($\mathit{CSR}$) with the capacity of the soil to resist liquefaction with 15% probability or cyclic resistance ratio (${\mathit{CRR}}_{PL(15\%)}$) as seen in Equation (6)(6) $$F_S=\frac{{\mathit{CRR}}_{PL\left(15\mathrm{\%}\right)}}{\mathit{CSR}}$$(6) . (6) $$F_S=\frac{{\mathit{CRR}}_{PL\left(15\mathrm{\%}\right)}}{\mathit{CSR}}$$(6)

Similar to $F_L,$ $F_S$ specifies the risk of liquefaction in depth, where values less than 1 will indicate the possibility of liquefaction. The cyclic stress ratio ($\mathit{CSR}$) is expressed by Equation (7)(7) $$\mathit{CSR}=0.65\left(\frac{a_{\max }}{g}\right)\left(\frac{{\sigma }_v}{{{\sigma }_v}_0}\right)r_d$$(7) , (7) $$\mathit{CSR}=0.65\left(\frac{a_{\max }}{g}\right)\left(\frac{{\sigma }_v}{{{\sigma }_v}_0}\right)r_d$$(7) $$\begin{array}{c}{r}_d=1-0.00765z;\mathrm{ }\mathit{for}\mathrm{ }z<9.2\mathrm{ }m\\ r_d=1.174-0.0267z;\mathrm{ }\mathit{for}\mathrm{ }z\ge 9.2\mathrm{ }m\end{array}$$ where $a_{\max }$ is the peak ground acceleration in g, g is the acceleration of gravity, $r_d$ is the reduction coefficient, z is the depth from the surface in m, and ${\sigma }_v$ and ${{\sigma }_v}_0$ as the total and effective overburden pressure, respectively, at depth $z.$ On the other hand, cyclic stress ratio $(\mathit{CRR})$ is calculated using Equation (8)(8) $$\mathit{CRR}=\exp \left\{\frac{\left[{\left(0.0073v_{s1}\right)}^{2.8011}-2.6168ln\left(M_w\right)-0.0099ln\left({{{\sigma }^{\prime }}_v}_0\right)+0.0028FC-0.4809{\varphi }^{-1}\left(P_L\right)\right]}{1.946}\right\}$$(8) , (8) $$\mathit{CRR}=\exp \left\{\frac{\left[{\left(0.0073v_{s1}\right)}^{2.8011}-2.6168ln\left(M_w\right)-0.0099ln\left({{{\sigma }^{\prime }}_v}_0\right)+0.0028FC-0.4809{\varphi }^{-1}\left(P_L\right)\right]}{1.946}\right\}$$(8) where $v_{s1}$ is the effective stress-normalized shear-wave velocity, $M_w$ is the simulated moment magnitude used in this study, ${{{\sigma }^{\prime }}_v}_0$ is the reference vertical overburden stress, $FC$ is the fines content (obtained from the SDS test). $v_{s1}$ is given in Equation (10)(10) $$\mathit{LPI}=\underset{0}{\overset{20}{\int }}{F}_{L}or{F}_s\text{*}\left(10-0.5z\right)dz$$(10) where $P^a$ is a reference vertical overburden stress of 100 kPa. (9) $$v_{s1}=v_s{\left(\frac{P^a}{{{{\sigma }^{\prime }}_v}_0}\right)}^{0.25}$$(9)

In contrast to the resolution of calculating $F_L$ in the SDS test, $F_S$ can only be calculated per layer present in the shear-wave velocity profile. LPI was then calculated using the formula of Iwasaki et al. (1984) given by, (10) $$\mathit{LPI}=\underset{0}{\overset{20}{\int }}{F}_{L}or{F}_s\text{*}\left(10-0.5z\right)dz$$(10) where $z$ is the soil depth. We also follow the criteria of Iwasaki et al. (1984) in liquefaction risk (Table 2).

Table 2. Liquefaction potential index (LPI) classification and description (after Iwasaki et al. 1984).

LPI Value	Description
0	Very low
0 to 5	Low
5 to 15	High
>15	Very high

The descriptive statistics of the resulting LPI in terms of predominant period values are summarized in Table 3 and their statistical distribution is presented in Figure 6 in a box and whisker plot.

Figure 6. Box and whisker plot of the sites in terms of predominant period and LPI.

Table 3. Summary of the descriptive statistics of the liquefaction potential classification in terms of predominant period.

	Very high	High	Low	Very low
No. of values	183	19	13	55
Minimum	0.12	0.11	0	0
25% percentile	0.45	0.19	0	0
Median	0.60	0.22	0.07	0.10
75% percentile	0.91	0.32	0.21	0.29
Maximum	1.31	0.78	0.21	0.29
Mean	0.66	0.27	0.09	0.09
Std. dev.	0.28	0.16	0.08	0.08
Std. error	0.02	0.04	0.02	0.01

There is a super majority of the sites that have very high liquefaction potential. These sites correspond to a minimum and maximum period range of 0.12 s to 1.31 s, respectively, but most sites are found within 0.45 s to 0.91 s. Comparatively, sites with high liquefaction potential are those found at a lower minimum and maximum period range of 0.11 s to 0.78 s, but most are found within 0.19 s to 0.32 s. It is interesting to note that Figure 6 shows high liquefaction potential can still occur at sites with a predominant period at about 0.12 s. In contrast, those sites with low to very low liquefaction potential contain sites that have no peaks up to periods of about 0.29 s. Another interesting find is that in terms of the predominant period, general trends can be observed separating the period classes. All (but one) long period sites have very high LPI. Those situated at moderate period sites have dominantly high to very high LPI with a few very low LPI and one low LPI. Short period sites have dominantly very low to low LPI with a few high to very high LPI. All sites with no peaks have evidently very low LPI with a few low LPI.

Discussion

Comparison with soil profiles

To understand the possible link between site response and liquefaction potential, we detail the comparison of the H/V curves in Figure 3 against the soil shear-wave velocity and N-value from Figure 5 which is generalized but contains necessary information to understand characteristic trends. The H/V curves contain diagnostic key features of the subsurface structure (Molnar et al. 2018; Molnar et al. 2022) and it is evident from the spectral curves that there are distinct contrasts in their shapes and peaks across the period classes. We already know that the predominant period is directly related to the overall thickness of the soil cover from Equation (2)(2) $$T=\frac{4h}{v_s}$$(2) whereas the amplitude scales with the impedance contrast between the soft soil and the rigid bedrock (Oubaiche et al. 2012; Uebayashi et al. 2012; Castellaro and Mulargia 2014; Likitlersuang et al. 2020) thereby providing information as to the stiffness of the soil compared to the bedrock. As an aid in interpreting H/V curve shapes, we used the papers of LeBrun et al. (2004) and Beroya et al. (2009) as well as the SESAME guidelines (SESAME 2004).

Sites a1, a2, a3 and a5 (Figure 3a) show predominant periods of 1.31 s, 1.06 s, 1.19 s and 1.03 s, respectively with amplitudes ∼4. These sites can have soil thicknesses of ∼38 to 56 meters. In comparison with the soil profiles (Figure 5, a1-a3 and a5), N-values of the upper 20 meters in these sites commonly range from 0.5 to 5 corresponding to very loose sands or soft to very soft clays whereas shear-wave velocities in the same depth range show about 92 m/s to 188 m/s. In fact, it is common for these low values to extend more than 30 meters at depth. In some areas such as in sites a4 and a6, although the predominant period is shorter (therefore less thick deposit), the amplitudes equate to about ∼7, indicative of a larger contrast in the acoustic impedance of the soft soil layer above the bedrock. As mentioned earlier, the engineering bedrock corresponding to the Guadalupe Formation has about ∼600 m/s to 760 m/s shear-wave velocities with the difference possibly attributed to the state of weathering. The apparent consistency of the low shear strength values and the considerable thickness of the deposit demonstrated a clear singular high amplitude peak qualitatively indicating a simple soil structure such as a homogenous single soft soil layer over a rigid bedrock with high impedance contrast.

Moderate period sites show more variability in the shape of the H/V curves. Sites b1, b3, b4, b5 and b7 (Figure 3b) show singular peaks which is the case for most of the H/V curves in this classification, however site b2 has a slightly broader maximum and site b6 has two peaks. Broad peaks and two peaks cases are characteristic of a more complex soil structure and are most likely attributed to strong lateral heterogeneity or irregularities (Guillier et al. 2006; Moisidi et al. 2012). Soil profiles corresponding to site b2 (Figure 5, b2) show the presence of a ∼2-meter-thick layer at very shallow depths with somewhat higher N-values (17 to 25) and shear-wave velocity (450 m/s) than the layers above and below it. Another REMI survey conducted about 45 meters away from where b2 was conducted showed thinning of this high velocity layer from 2 meters to just 90 centimeters. On the other hand, the soil profile in site b6, particularly the N-value profile (Figure 5, b6) shows a very dense sand (N-value reaching 50) present at about 5 meters then lessens at depth. Shear-wave velocity profiles taken nearby site b6 did not resolve this layer. It is questionable if the presence of these shallow high velocity layers is sufficient to complicate H/V curve interpretation and it was not explored in more detail. Therefore, these sites lean towards an inconclusive analysis of the amplitudes, although the predominant period is reliable and usable. Other reasons supplementary to the sufficiency of this layer to complicate H/V curve interpretation is also possible such as the presence of subsurface sloping interfaces between hard and soft layers (SESAME 2004) but was not explored in this paper. Nevertheless, the sites in this period classification are underlain by ∼13 to 30 meters of soil deposit having N-values commonly ranging from 1 to 15 corresponding to the increased presence of stiff and medium dense soils at the shallow subsurface. The shear-wave velocities also reflect increased values ranging from ∼105 m/s to 250 m/s, with some layers reaching 300 m/s at depth. LPI values in the long period and moderate period sites are dominated by high to very high LPI with only a few low and very low LPI, as previously stated (Figure 6).

At short period sites, the location of the engineering bedrock is nearer to the surface, commonly less than 15 meters as shown in the shear-wave velocity profiles in sites c1, c3, c5 and c6 (Figure 5) coincident with the thickness calculations ranging from ∼8 to 14 meters (Figure 4). These thicknesses translate to predominant periods of 0.19 s, 0.11 s, 0.13 s and 0.10 s, respectively (Figure 3c). Soil N-value ranges are generally higher from ∼3 to 25 corresponding to soft to very stiff clays and loose to dense sands. Similarly, shear-wave velocities at these sites are commonly ranging from ∼200 m/s to 360 m/s, with some layers reaching 450 m/s near the surface. These values of high shear strength at relatively lesser thickness of the soil deposit reflect a low amplitude peak equating to 2 to 2.5 at shorter periods from the H/V curves (Figure 3c). Contrary to the direct relationship of the predominant period with the overall soil thickness, short periods can also result from a low impedance contrast. This is shown in site c2 where the soil layer has an average shear-wave velocity of 335 m/s (Figure 5, c2). This is due to the strong tradeoff of the soil thickness and average shear-wave velocity represented in Equation (2)(2) $$T=\frac{4h}{v_s}$$(2) . In this period classification, the sites are dominated by low to very low LPI (Figure 6) but there are still sufficiently low shear strength layers with considerable thicknesses that could generate a high to very high LPI, manifested by clear, albeit low amplitude peaks in the short period range in Figure 3c.

Sites d1 to d5 have H/V curves with no observable peaks (Figure 3d). This simply means that measurements were placed on hard rock areas or there is an insufficient thickness of low shear strength material on top of rigid bedrock to produce a prominent peak in the H/V curve. These sites were mostly observed on measurements placed on or near the tuffaceous rocks of the Guadalupe Formation. As mentioned earlier, LPI values in this period classification are mainly very low LPI with only a few low LPI. It is interesting to note that with regards to those with low LPI classification (Figure 5, d3 and d4), a ∼2.5- meter-thick sand layer with N-values at about 1 to 5 is enough to trigger liquefaction. In terms of shear-wave velocity, this relatively thin layer was not resolved in much detail. It is very common in these sites that the SDS test has only penetrated a few tens of centimeters with a select few reaching a couple of meters. A unique limitation of the test is that its probing termination is controlled primarily by the shallow presence of very stiff or dense material (N-values > 15) (Orense et al. 2019; Daag et al. 2023). Early termination of the test already supplements the geology of the area indicative of a rock site.

Adding to the understanding of liquefaction at depth as well as taking advantage of the data on the depth of the termination of the probing of the SDS test, $h_{\mathit{SDS}}$ (assuming that it is the approximate expression of the bedrock below), a power-law regression equation was generated (Ibs-von Seht and Wohlenberg 1999; Delgado et al. 2000; Parolai et al. 2002) between predominant period, $T,$ and depth in the form of, (11) $$h_{\mathit{SDS}}=a{T}^b$$(11) where $a$ and $b$ are empirical parameters (Figure 7).

Figure 7. Relationship between predominant period and SDS test penetration depth.

There exists a strong correlation between the predominant period and the depth of penetration of the SDS test. This relationship was found and useful in terms of a first approximation of the thickness of the liquefiable layer since the termination of the SDS test at encountering N-values > 15 may suggest that material above may likely liquefy, and the material below may be the starting depth to where the material becomes very stiff or dense assuming that the stiffness of the soil deposit increases linearly with depth. Some papers regard the depth of the soil material corresponding to an N-value > 15 as an important target material for anchoring in shallow foundation studies (Kishida 1966; Towhata 2008; Mittal et al. 2013).

Comparison with LPI

The correlation of the predominant period with the LPI was investigated since it is practice to obtain the strength of the relationship by performing linear regression (Figure 8). This also complements the discussion in the previous section. Figure 8 shows that there is a moderately strong correlation between LPI and predominant period. Scattering in the predominant period can be observed across all LPI classifications specially in the long period sites suggesting a simple linear relationship is not achievable, although there is somewhat a linear trend that occurs in the short period to moderate period range at 0.2 s to 0.6 s. Overlaps in the predominant period are also observed, especially in the 0.1 s to 0.3 s period range where all LPI classifications exist, similar to Figure 6, wherein long and moderate periods occupy very high and high LPI respectively.

Figure 8. Predominant period against LPI showing moderately strong correlation. Red circles are sites with very high LPI, purple circles are sites with high LPI, yellow circles are sites with low LPI and green circles are sites with very low LPI.

The reason behind the scattering is assumed to be a function of the distance of long period sites to the West Valley Fault (WVF). From Equations (5)(5) $$R=C_w{R}_L$$(5) and (8)(8) $$\mathit{CRR}=\exp \left\{\frac{\left[{\left(0.0073v_{s1}\right)}^{2.8011}-2.6168ln\left(M_w\right)-0.0099ln\left({{{\sigma }^{\prime }}_v}_0\right)+0.0028FC-0.4809{\varphi }^{-1}\left(P_L\right)\right]}{1.946}\right\}$$(8) , the maximum ground acceleration at each site was factored in. Peak accelerations are greatly controlled by the distance to the source fault. The farther the site is to the fault, the lower the value of acceleration and the higher the value of the factor of safety. Probabilistic peak accelerations in the study area are taken from The Philippine Earthquake Model of PHIVOLCS (PHIVOLCS 2017) and range from 0.23 to 0.60 g. A number of long period sites are located in Bulacan in the north, far from the WVF thus resulting in lower acceleration values. In contrast, short period and moderate period sites are located more proximal to the WVF therefore having higher acceleration values. Nevertheless, a minimum LPI value of 16 is observed in the long period scattered region already corresponding to a very high LPI.

Despite the moderately strong correlation, it appears that we can possibly extract a boundary of the predominant period between high and low LPIs. We assume, based on the interquartile ranges from Figure 6, that we can group very high and high LPIs and low and very low LPIs resulting in just two classifications – high and low potential. An LPI equal to 5 serves as the boundary between the two classes. In searching for the equivalent boundary in terms of predominant period, a confusion matrix was employed to explore the value with highest accuracy to separate the two classes (Figure 9).

Figure 9. Confusion matrix-derived relationship showing peak at 0.2 s with 94% accuracy. Periods 0 and 1 are equivalent to < 0.2 s and > 0.2 s, respectively. LPIs 0 and 1 are equivalent to <5 and >5, respectively.

A confusion matrix summarizes the classification performance of a classifier with respect to some test data (Ting 2010). It can be as simple as a 2 × 2 matrix, listed in one dimension by the actual class (predominant period) and in the other by the class that the classifier assigns (LPI). We employed a Boolean logic to further simplify analysis. Values of 0 and 1 were assigned corresponding to low and high potential to liquefaction, respectively, to both predominant period and LPI. True positives (long period correctly predicting high LPI) are where 1s agree and true negatives (short period correctly predicting low LPI) are where 0s agree. We arrived at 0.2 s with 94% accuracy suggesting that in 94% of the observations, 0.2 s correctly classified high from low LPI. This value of the predominant period may be used as a first pass indicator that liquefaction potential is high in the study area. Microtremor measurements can thereby serve as a tool to fill data gaps and complement site-specific liquefaction delineation where in situ geotechnical data is limited or not readily available, considering the cost of a detailed geotechnical analysis.

As mentioned earlier, a hallmark to strong ground shaking and liquefaction is the amplified and prolonged duration of seismic energy. From the earlier comparisons, we have observed somewhat clear direct relationships in terms of the general period classifications with the liquefaction potential. It is evident that longer predominant periods with relatively high amplitudes can have higher potential for liquefaction and the opposite is true that shorter predominant periods with relatively lower amplitudes can have lower potential for liquefaction. The apparent reason for this is the link between the predominant period with the soil N-values and shear-wave velocities. Varying thickness and shear strength of the soil material is reflected upon the shape of the H/V curve and its predominant period and relative amplitude.

Considerations on topography and groundwater table

It is clear that there is a possible link between site response and liquefaction potential from the previous discussion. This section emphasizes the role of topography and groundwater table depth in facilitating possible misinterpretations of the microtremor H/V. Figure 10 shows the distribution of the sites in terms of topography and groundwater table depth and Figure 11 shows the distribution of the groundwater table across the study area.

Figure 10. Plot showing summary of the distribution of the experimental sites in terms of elevation and groundwater table depth.

Figure 11. Distribution of water table level across the study area obtained from ground Penetrating Radar (GPR).

The characteristics of thick and soft soil deposits influence seismic energy at the near surface making the sites more prone to amplified ground shaking and liquefaction (Qodri et al. 2021). A scenario that would oppose this statement is that if the site is located in a topographical high (e.g. plateau, hill, mountain) but still produces a prominent peak in the H/V curve. Such topographical position would already contradict the possibility of liquefaction. Similarly, one can still obtain long periods even if the subsurface material is clay as observed in Figure 5 or if the depth of the groundwater table is significantly deep. This suggests that an H/V curve with a prominent peak at a given period range does not necessarily equate to a high LPI. Caution should therefore be exercised in using the predominant period at 0.2 s in separating high and low liquefaction potentials in a different geological environment, and it is important that information on geology, geomorphology and groundwater table depth must be considered foremost prior to interpreting the H/V results when applied to liquefaction studies.

Conclusion

This research demonstrated the application of the microtremor H/V in obtaining site response characteristics in coastal GMMA and explored the relationship of site response parameters to soil liquefaction potential. Listed below are the main findings of the research:

• The distribution of the predominant period and its relative amplitude in coastal GMMA ranges from no peaks observed in the rocky coast of Cavite in the south as well as in the central plateau of Metro Manila in the east to 1.31 s with high amplitudes in the coastal lowlands of Bulacan in the north. This translates to a general increasing trend in overall soil thickness from east to west and south to north. This also augments earlier microtremor data of Abeki et al. (1996) and Narag et al. (2000) which are limited to Metro Manila.

• The distribution of the LPI grouped with the predominant period classification reveals that long period sites (>0.6 s) have generally very high LPI, moderate period sites (0.2 s to 0.6 s) have dominantly high to very high LPI with a few very low LPI, short period sites (0.05 s to 0.2 s) have dominantly very low to low LPI with a few high to very high LPI, an sites with no peaks have generally very low LPI.

• The relationship of site response characteristics to soil liquefaction potential is the fundamental link between shear strength properties of the soil. Variations in the soil N-values and shear-wave velocities are generally reflected in the microtremor H/V curves.

• Correlation of the predominant period against the LPI yielding a moderately strong relationship (r² = 0.51) as well as the statistical extraction of a predominant period boundary (0.2 s) from high and low liquefaction potential suggests its potential use as a first-pass tool in liquefaction assessment. This boundary can be useful for future studies on liquefaction assessment in coastal GMMA, specially where subsurface data is limited or not readily available.

The microtremor H/V is a powerful technique for obtaining semi-quantitative information of the subsurface. However, it is noted that it lacks the quantitative layer-per-layer analysis that other in situ geotechnical techniques provide. This is particularly important in more dynamic environments where the soil deposits exhibit lateral heterogeneity in properties and therefore complicating microtremor H/V interpretation. A dense array of microtremor H/V measurements may be recommended to deliver the necessary level of accuracy in those settings. Hence, the microtremor H/V technique cannot replace conventional techniques in the assessment of liquefaction, but it can serve as a useful complementary tool. We also point out that the results of the liquefaction assessment are limited to coastal GMMA, and we recommend that a similar kind of analysis should be done in other geological environments.

Author’s contributions

A. D. conceived this study and handled the overall project supervision, administration, and funding acquisition. L. E. A. and O. L. carried out the data acquisition, processing and analysis, and manuscript preparation. A. D., R. G. and R. S. Jr. provided resource and review of the manuscript.

Acknowledgements

The authors are most grateful to the staff of the Liquefaction Project namely, A. T. Serrano, M. J. V. Reyes, O. P. C. Halasan, K. S. Sochayseng, A. A. T. Magnaye, M. C. Dela Cruz, A. O. Amandy, and E. J. M. Arnoco, for their invaluable insights and expertise for the holistic improvement of this paper. The authors also thank M. I. T. Abigania, M. P. Dizon, D. J. L. Buhay, and E. D. Mitiam for their assistance throughout the entire project duration. The authors also express their gratitude to the SpecificEQ Project for the introduction of the methodology and the software. The authors are also thankful to the Department of Education (DepEd), particularly to the Educational Facilities Division for the administrative assistance and comments that helped in various stages of the project. Special thanks are given to the school principals and Division Engineers for the hospitality and assistance during the fieldwork campaigns.

Data availability statement

The data to support the findings of this study are available from the corresponding author, A. D. upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was funded by the Department of Science and Technology – Philippine Council for Industry, Energy, and Emerging Technology Research and Development (DOST-PCIEERD).