http://journal.uir.ac.id/index.php/JGEET E-ISSN : 2541-5794 P-ISSN : 2503-216X Journal of Geoscience, Engineering, Environment, and Technology Vol 08 No 02-2 2023 Special Edition Special Issue from “The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022” Amirudin et al./ JGEET Vol 08 No 02-2 2023 23 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Seismic Vulnerability Analysis Using the Horizontal to Vertical Spectral Ratio (HVSR) Method on the West Palu Bay Coastline Amirudin1, Iktri Madrinovella1*, Sofian2 1 Geophysical Engineering, Pertamina University, DKI Jakarta 12220, Indonesia 2 Palu Class I Geophysical Station, Meteorological, Climatological and Geophysical Agency, Indonesia * Corresponding author: iktri.madrinovella@universitaspertamina.ac.id Received: May 20, 2023. Revised : May 31, 2023. Accepted: June 10, 2023. Published: July 31, 2023 DOI: 10.25299/jgeet.2023.8.02-2.13879 Abstract This research was carried out to make a map of the dominant frequency (f0), amplification factor (A0), seismic susceptibility index (Kg), Vs30 , Sediment Layer Thickness (H) and Peak ground acceleration (PGA). Microtremor measurements were carried out with a three - component seismometer of the TDL-303S type as many as 27 measurement points. The data was analyzed by the Horizontal to Vertical Spectral Ratio (HVSR) method. The PGA calculation was carried out using the Kanai equation with a reference to the Palu -Donggala earthquake on September 28, 2018. The results showed that the distribution of the dominant frequency value (f0) ranged from 0.4149 Hz-0.8869 Hz, the soil amplification factor (A0) ranged from 2,199–4,884, the seismic vulnerability index (Kg) ranged from 8.79 s2/cm- 41.41 s2/cm, the shear wave velocity to a depth of 30 meters ( Vs30) ranged from Vs30 197.7 m/s-320.2 m/s , the thickness of the sedimentary layer ranges from 260.3 m-291.1 m and the peak ground acceleration (PGA) of Kanai ranges from 137.3 gal – 234.2 gal by using Mw 7.4 earthquakes with an intensity scale (MMI) VI to VII. The coastal area of West Palu bay has an intermediate seismic vulnerability II to a high seismic vulnerability IV so that it will be vulnerable in the event of an earthquake disaster. Are as that have a very high vulnerability index are in the upper western and easternmost regions while those with a lower level tend to have a lower vulnerability index value. Keywords: Passive seismic, Microtremor, HVSR, Seismic vulnerability index, Peak ground acceleration, Palu-Donggala 1. Introduction The motion of the Palu-Koro Fault activity is the main cause of the earthquake disaster that occurred in Palu City. The Palu–Koro fault is approximately 240 km, which leads from northern Palu to Southern Palu (Malili) until it reaches Bone Bay. The fault is cytostral and very active with a fault movement speed of about 25 to 30 mm per year, causing frequent earthquakes in the Palu city area. The Palu-Koro Fault passes through The West Palu District, where in West Palu District has many important buildings that become the economic center for the people of Palu City, such as universities, schools, hotels, malls, banks, health centers, village offices and sub-district offices. The geological map shows that the research area (West Palu Bay coastline) has the lithology of alluvium and coastal deposit, which means the area has less dense rock. The West Palu Bay coastline is likely vulnerable if an earthquake occurs near the area. Fig 1. Geological Map Review Palu Sheet, Central Sulawesi (BMGK Kota Palu, 2018) 24 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 The Horizontal to Vertical Spectral Ratio (HVSR) method from Nakamura is one of the methods that can be used to utilize microtremor data. It can generate the dominant frequency value and amplification of the soil in an area, to compute the seismic vulnerability index and the thickness of the sedimentary layer. The purpose of this study was to determine the level of seismic vulnerability in the study area by making a map of the distribution of the dominant frequency value, amplification factor, seismic vulnerability index, shear wave velocity and sediment layer thickness. 2. Material and Methods 2.1 Microtremor Microtremor is a continuous subsurface vibration that has a low amplitude ranging from 0.2-20 Hz, the source of vibration comes from various kinds of surface activities such as wind interactions with trees, ocean waves, vehicles, and human activities that can cause below-surface vibrations (Kanai, 1983). These harmonic vibrations are retained in the sedimentary layer on the surface and are constantly reflected by the boundary plane of the layer. 2.2 HVSR method Nakamura (2008) The spectral ratio of horizontal components to vertical component called the HVSR method with the assumption that shear waves predominate in microtremor and surface waves (Rayleigh and Love waves) are ignored. The transfer function between wave vibrations in sediments and bedrock is considered as the ratio of the horizontal and vertical spectra of surface vibrations. This means that the peak period of H/V and the peak value itself correspond to the dominant period and its amplification factor. 2.3 Amplification Factor (A0) Nakamura (2008) states that the value of the soil strengthening factor (amplification) is related to the ratio of the impedance contrast of the surface layer with the layer below. If the contrast ratio of the impedance of the two layers is high, the value of the strengthening factor is also high, and vice versa. 2.4 Dominant Frequency (f0) The dominant frequency or the natural frequency of the soil or rock is the most dominant or most frequently appearing frequency value. The parameter of the dominant frequency is that it can indicate or be a reference to the type and characteristics of a rock. 2.5 Shear wave velocity at the depth of 30 m (Vs30) Vs30 is the average speed of the surface shear wave is up to a depth of 30 meters. It is one of the local site effect variables that contribute to variations in the level of earthquake shocks on the surface. To estimate the speed of the shear wave (Vs) can be done by various methods, including drilling (borehole). The calculation of Vs30 is using topography slope, which is published by USGS (United States Geological Survey) as open-source data based on Allen and Wald (2007). Other reference, Herak (2008) introduced a technique for estimating the speed of shear waves (Vs) using HVSR microtremor inversions based on body waves by the equation: 𝑉𝑠30 = ℎ30 ∑ ℎ𝑖 𝑉𝑠𝑖 (1) Vs30 is the shear wave velocity at the depth of 30 m (m/s), h30 is the depth of 30 m, hi is the depth of the i-th layer (m) and Vsi is the shear wave velocity of the i-th layer (m/s). 2.6 Sediment thickness (h) The characteristics of sedimentary rocks that are loose, are not compact with relatively low mass density so that the propagation speed of shear waves is slower than that of the bedrock below. The value of f0 is related (inversely proportional) to the thickness of the sediment layer where the greater the dominant frequency value, the thinner a sedimentary layer, on the contrary, when the dominant frequency value is smaller, the thicker the sedimentary layer. The thicker the sediment layer results in an increased risk or impact of soil susceptibility to earthquakes. Here is the formula for the relationship of the thickness of the sedimentary layer to the dominant frequency. ℎ = 𝑉𝑠30 4𝑓0 (2) h is the sediment thickness (m), Vs30 is the shear wave velocity at the depth of 30 m (m/s), fo is the dominant frequency (Hz). 2.7 Seismic vulnerability index (Kg) Seismic vulnerability index is used to identify the vulnerability degree of the subsoil that deformations due to earthquakes occur. High seismic vulnerability index values will generally be found on soils with soft sedimentary rock lithologies. A high seismic vulnerability index value can identify that the area is prone to earthquakes and will experience strong shocks, and vice versa, small seismic vulnerability index values are generally found in soil layers with solid building rock lithologies so that when an earthquake occurs it does not experience much strong shaking. 𝐾𝑔 = 𝐴0 2 𝑓0(𝜋 2𝑉𝑏 ) (3) Ao is the amplification factor, fo is the dominant frequency (Hz) and Vb is the shear wave velocity of the bedrock (m/s). 2.8 Peak Ground Acceleration (PGA) Peak ground acceleration (PGA) is the largest soil vibration acceleration value that has ever occurred somewhere caused by an earthquake wave. Earthquake shocks felt on the surface can be described descriptively with an illustration scale based on the level of shock and impact. The formula of PGA uses Kanai method (1966) modified by Douglas (2004). 𝑃𝐺𝐴 = 5 √𝑇0 10 0.61𝑀−(1.66+ 3.6 𝑅 ) log 𝑅+0.67 1.83 𝑅 (4) PGA in gal, M is the magnitude of an earthquake near the site, R is the distance between the hypocenter and the site (km). 2.9 Data Acquisition The acquisition of microtremor data was carried out on the coastline of West Palu Bay, Palu City, Central Amirudin et al./ JGEET Vol 08 No 02-2 2023 25 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Sulawesi accompanied by BMKG Class I Palu on March 10-14, 2022. The location of each station was determined by using GPS (coordinates and elevations) by leveling or adjusting the position, connecting the seismometer, measuring in 60 minutes for each station, saving the data in trace (*trc) format, and reformatting to miniseed (MSD) format. Fig 2. Microtremor data acquisition equipment The instrument used for data acquisition is Seismometer DS-4A to measure the ground motion, digital portable seismograph TDL-300S with sampling frequency 100 Hz, geological compass to define the sensor direction, GPS Garmin for the measurement point location, GPS antenna to define the time and location, connector cable to connect the instruments, laptop to monitor and check the recorded signal in digital portable seismograph, and logbook to write the coordinates, time of measurements and field conditions. Fig 3. Microtremor data acquisition location map 2.10 HVSR Analysis with Geopsy Software The microtremor data is obtained from the 3 (three) component short-period digital seismometer, which consists of 2 (two) horizontal components (north-south direction and east-west direction) and 1 (one) vertical component. The data is stored in digital waveform with sampling rate 100 Hz, and it must be transformed into frequency domain so the spectral ratio can be analyzed. The process of transforming the data from time domain into frequency domain is called Fast Fourier Transform (FFT). 2.11 Fast-Fourier Transform: 𝐹(𝑡) = ∫ 𝑓(𝜔)𝑒𝑖𝜔𝑡𝑑𝑡 ∞ −∞ (5) 26 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 4. 3-Component Microtremor Signal The first step of HVSR analysis is to define the parameter of filtering, windowing, smoothing. The filtering used is low pass filter of 20 Hz, according to the characteristics of the low frequency microtremor signal. The windowing length is 60 s which means the frequency spectrum is computed for every 60 s of data. However, there are high frequency noise that will disturb the data. It must be removed by determining the data by using the ratio of short-term average (STA) to long-term average (LTA) of the data amplitude. The length for STA determination used is 1 s and LTA is 30 s. Fig 5. Windowing Every signal in every time window is transformed into frequency spectrum. The frequency spectrum from the horizontal component (North-South and East-West) is being computed using quadratic mean. The ratio of H/V component means the ratio of the frequency spectrum of the vertical component (Z) to the quadratic mean of the horizontal components. Amirudin et al./ JGEET Vol 08 No 02-2 2023 27 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 6. The calculation of H/V spectral ratio (https://www.geopsy.org/documentation/geopsy/hv.html) The maximum number of the H/V of the average spectrum means the value of the amplification factor (Ao), and the frequency of the Ao is the value of the dominant frequency. The average H/V curve is determined by the length of window (60 s), so the standard deviation value is used to validate the curve reliability and peak clarity using the SESAME criteria (2004). https://www.geopsy.org/documentation/geopsy/hv.html 28 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 7. HVSR curve Fig 8. H/V curve criteria for reliable curve and clear peak (SESAME criteria) The H/V curves derived by all the 27 measurement points generally meet the SESAME criteria (2004) almost 100% for the reliable curve and 80% for the clear peak. The seismic vulnerability index is computed using the value of the dominant frequency and the amplification factor of the H/V curve. And the each point data of the H/V curve is used to generate the Vs30 value using depth inversion method. 2.12 HVSR Inversion with Dinver Software The Dinver software is used to calculate the shear wave velocity (Vs) using the ellipticity curve method, which is one of the HVSR inversion methods that can be used to identify subsurface structures to obtain a grounf profiles curve that shows the Vs value of each depth. The H/V curve model is constructed using model parameters in the form of primary wave velocity (Vp), Amirudin et al./ JGEET Vol 08 No 02-2 2023 29 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 secondary wave velocity/shear wave velocity (Vs), and layer depth (H), as well as input condition parameters (poisson ratio) and ρ (rock density). The matching process is carried out repeatedly, each iteration (iteration) is corrected to each of the initiation parameters, until the final model curve has the smallest mismatch of the observed HVSR curve. The result of the smallest misfit is obtained in the form of Vp and Vs each depth corresponding to the HVSR curve from the observation data 𝜎. The parameters used in the inversion of the HVSR curve can be seen in table 1 where the number of layers used in this study is 6 layers, the model parameters are P-wave velocity (Vp), shear wave velocity (Vs), Poisson’s ratio (𝜎), and density (ρ). The maximum depth used is 300 m. Table 1. Initial Model Parameter Vs (m/s) Vp (m/s) Poisson’s Ratio Density (Kg/m3) Maximum Depth (m) Layer 1 Layer 2 Layer 3 Layer 4 100-400 300-460 460-580 480-760 160-1000 480-1150 576-1450 768-1900 0.2-0.5 2000 300 Layer 5 660-1000 1056-2500 Layer 6 900-1300 1440-3250 The results of the HVSR inversion are the ellipticity curve and the ground profile curve. The ellipticity curve serves to check the similarity between the initial parameter of the model and the observation. The ground profile curve generates the Vs value to the depth below the surface. Fig 9. The ellipticity curve and the ground profile curve. The Dinver software used the algorithm introduced by Wathelet (2008) based on Sambridge (1999) and Wathlet (2005). 30 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 10. The flow chart of HVSR analysis and inversion 3. Results and Discussion The distribution map shows the dominant frequency value in range 0.4149 – 0.8869 Hz. Based on the classification of the dominant frequency of Nakamura, the soil type is classified as the alluvial formation with frequency less than 2.5 Hz and the sediment thickness above 30 m. The geology study shows that the area is formed by the delta sedimentation, top soil, sand, mud, and gravel. The eastern part of the fault plane shows the significantly low dominant frequency, which indicates the high potential of damage if an earthquake occurs. Amirudin et al./ JGEET Vol 08 No 02-2 2023 31 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 11. The dominant frequency map The amplification factor shows the strength of the shaking after the seismic wave propagates in the less dense rock is 2.199-4.884 times stronger than in the bedrock. The amplification factor of the research area shows the medium level, which the western and the center area shows the biggest amplification factor number. The amplification number can be increased if the rock has undergone deformation (weathering, folding, roughing) that changes the physical properties of the rock. Fig 12. The amplification factor map The shear velocity at the depth of 30 m is called Vs30, and it ranges between 197.7 – 320.2 m/s. The Vs30 is generated by HVSR curve inversion, and it shows the lower value in the center part of the research area or the center of Palu City. This area has lower Vs30 value (<200 m/s) and fo value (0.4 - 0.5 Hz), also has the higher amplification factor value (>4), which needs to be concerned. Fig 13. The Vs30 map 32 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Based on the sediment layer thickness map in range 260.3 – 291.1 m, it shows the center of Palu City also has the thicker layer (> 280 m). In the north-western part of the fault plane, it also shows the thicker sediment layer (> 280 m), lower fo value (0.5 – 0.6) and higher amplification factor value (>4). It has crowded infrastructure in this area, such as malls, hotels, markets, hospital and schools. Fig 14. The sediment layer thickness map The result of the seismic vulnerability index by using the amplification factor and the dominant frequency data shows the range between 8.791 – 41.41 s2/cm which is classified as high index. The center area of Palu City and the north-western part of the fault plane has the highest index (more than 30 s2/cm). Table 2. The Classification of the seismic vulnerability index (BMKG Palu, 2018) Kg Index < 0.16 0.16 - 6 6 – 9 9 – 12 12 – 16 16 – 18 Low Middle I Middle II Middle III Middle IV High I 18 – 21 High II 21 – 24 High III > 24 High IV Based on the seismic vulnerability index result and Table 2, the West Palu Bay coastline has the vulnerability index middle to high. The geological map shows the entire area consists of alluvium and coastal deposit, and the highest index is in the western and the center part of the research area (High III – High IV). The highest index shows the most vulnerable area of the West Palu Bay coastline, which can be damaged severely if an earthquake occurs near the area. The western part (West Palu) and the center part (Central Palu) are the location of houses, malls, hotels, markets, hospitals. The review of ground motion after Palu-Donggala earthquake 2018 reported by BMKG shows the response of citizens that has significantly corrected the earthquake intensity in Palu City. It is consistent to the seismic vulnerability index, the area of Palu City has greater damage when the earthquake occurs. Fig 15. The seismic vulnerability index map Amirudin et al./ JGEET Vol 08 No 02-2 2023 33 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 Fig 16. BMKG Shakemap of the Palu earthquake, 2018 before (left) and after (right) correction The calculation of Peak Ground Acceleration (PGA) is also conducted to view the impact of the dominant frequency. Thera are many empirical ways to calculate PGA by using parameter magnitude and distance between the source and the site. By using the Kanai formula, the PGA calculation can consider the geological effect that is represented by the dominant frequency value. The calculation uses the Palu-Donggala earthquake (Mw 7.4, September 28, 2018) as the magnitude and the source location parameter. The PGA result ranges between 137.3 – 234.2 gal, which the highest PGA located in the south-western part of the fault plane. The PGA result is consistent to the dominant frequency, but inconsistent to the other parameter (Ao, Vs30, Kg, H). It can be concluded that this PGA formula does not associate with the result effectively, but it has strong relationship with the dominant frequency. Fig 17. The PGA map (left) and the BMKG shakemap of the Palu-Donggala earthquake on September 28, 2018 (BMKG Palu, 2018) The result of PGA using Kanai formula and BMKG shakemap shows the similarity in the south-western part of the fault plane, but inconsistent to the other part. 7. Conclusion 1. The south-western part of the fault plane and the center part of the Palu City are the most potentially vulnerable area with the lower dominant frequency (fo 0.4 – 0.6 Hz), lower shear velocity (Vs30 < 200 m/s) higher amplification factor (A0 > 4), higher seismic vulnerabity index (Kg 30 – 41 s2/cm), thicker sediment layer (H 280 – 290 m). It is pretty crowded in the south-western part area which has malls, hotels, markets, hospital and schools. The center of the Palu City also has many houses, university, and markets. It is also consistent to the shakemap of BMKG (The review of Palu-Donggala earthquake 2018) which shows the greater intensity after the correction (felt by people). 2. The result shows that the area of West Palu and Central Palu are the most vulnerable area, 34 Amirudin et al./ JGEET Vol 08 No 02-2 2023 Special Issue from The 1st International Conference on Upstream Energy Technology and Digitalization (ICUPERTAIN) 2022 they have the potential to be severely damaged if a great earthquake occurs. 3. The PGA value is affected majorly by the dominant frequency factor, so the result of the PGA value does not associate effectively with other parameters (Ao, Vs30, Kg, H). 4. The research area does not cover the map of the liquefaction or affected area of the 2018 earthquake. But the vulnerable area is still counted as the area that can be affected by the motion of the Palu-Koro fault. 5. These results and methods can be implemented in other study area, as long as the area has similar geological and seismological conditions. However, the validation used in this research is limited, so it can have potential error in analysis. To improve the confidence of this result, other similar method or computation should be conducted. Acknowledgment We acknowledge the Meteorological, Climatological and Geophysical Agency (BMKG), especially Palu Class I for providing the data and the opportunity to conduct the field acquisition. References Allen, T. I., Wald, D. J. 2007, Topographic slope as a proxy for global seismic site conditions (VS 30) and amplification around the globe: U.S. Geological Survey Open-File Report 2007-1357, 69. https://doi.org/10.1785/0120060267 BMKG Kota Palu., 2018. Laporan Akhir Analisis Indeks Kerentanan Seismik Kota Palu. BMKG Kelas I Palu. Kota Palu. Douglas., 2004, Ground Motion Estimation Equation 1964- 2003, South Kensington Campus Press, London. Geopsy Tutorial: H/V Measurement (Modified 2008-10- 21) https://www.geopsy.org/documentation/geopsy/hv.ht ml Herak, M., 2008, ModelHVSR: a Matlab tool to model horizontal-to-vertical spectral ratio of ambient noise. Computers and Geosciences 34, 1514-1526. https://doi.org/10.1016/j.cageo.2007.07.009 Kanai, K., 1983, Seismology in Engineering, Tokyo University, Japan. Kanai., 1966, Improved Empirical Formula for Characteristics of Stray Earthquake Motion. Proceedings of The Japanese Earthquake Symposium (1-4). Japan: Reported in Trifunac & Brandy (1975). Nakamura, Y., 2008. On The H/V Spectrum, The 14th World Conference on Earthquake Engineering: Beijing, China. Sambridge, M., 1999, Geophysical Inversion with a neighbourhood algorithm: I. Searching a parameter space. Geophysical Journal International 138, 479- 494. https://doi.org/10.1046/j.1365- 246X.1999.00876.x SESAME., 2004, Guidelines for the Implementation of the H/V Spectral Ratio Technique on Ambient Vibratio Measurement, Processing and Interpretation, European Commission - Reseacrh General Directorate. Wathelet, M., 2005, Array recordings of ambient vibrations: surface-wave inversion. PhD thesis Universite de Liege Belgium. Wathelet, M., 2008, An improved neighborhood algorithm: parameter conditions and dynamic scaling. Geophysical Research Letters 35, L09301. https://doi.org/10.1029.2008GL033356.pdf © 2023 Journal of Geoscience, Engineering, Environment and Technology. All rights reserved. This is an open access article distributed under the terms of the CC BY-SA License (http://creativecommons.org/licenses/by-sa/4.0/). https://doi.org/10.1785/0120060267 https://www.geopsy.org/documentation/geopsy/hv.html https://www.geopsy.org/documentation/geopsy/hv.html https://doi.org/10.1046/j.1365-246X.1999.00876.x https://doi.org/10.1046/j.1365-246X.1999.00876.x http://creativecommons.org/licenses/by-sa/4.0/ http://creativecommons.org/licenses/by-sa/4.0/