Shallow Shear-Wave Velocity Beneath Jakarta, Indonesia Revealed by Body-Wave Polarization Analysis

Abstract

Noting the importance of evaluating near-surface geology in earthquake risk assessment, we explored the application to the Jakarta Basin of a relatively new and simple technique to map shallow seismic structure using body-wave polarization. The polarization directions of P-waves are sensitive to shear-wave velocities (V_s), while those of S-waves are sensitive to both body-wave velocities. Two dense, temporary broadband seismic networks covering Jakarta city and its vicinity were operated for several months, firstly, from October 2013 to February 2014 consisting of 96 stations, and secondly, between April and October 2018 consisting of 143 stations. By applying the polarization technique to earthquake signals recorded during these deployments, the apparent half-space shear-wave velocity (V_s^ahs) beneath each station is obtained, providing spatially dense coverage of the sedimentary deposits and the edge of the basin. The results showed that spatial variations in V_s^ahs obtained from polarization analysis are compatible with previous studies, and appear to reflect the average V_s of the top 150 m. The low V_s that characterizes sedimentary deposits dominates most of the area of Jakarta, and mainly reaches the outer part of its administrative margin to the southwest, more than 10 km away. Further study is required to obtain a complete geometry of the Jakarta Basin. In agreement with previous studies, we found that the polarization technique was indeed a simple and effective method for estimating near-surface V_s that can be implemented at very low-cost wherever three-component seismometers are operated, and it provides an alternative to the use of borehole and active source surveys for such measurements. However, we also found that for deep basins such as Jakarta, care must be taken in choosing window lengths to avoid contamination of basement converted phases.

1. Introduction

The growth of the global population over the past century, combined with the accelerating pace of urbanization, has resulted in the explosive growth in the number of megacities (population over 10 million). Especially in the developing world, many of the buildings in these cities have not adhered to earthquake resilient construction practice with the result that “an epicentral hit on a megacity has the potential to cause 1 million fatalities” [1]. Also, many of these cities are concentrated in sedimentary basins, because of their flat topography and access to fertile soils, sources of water and maritime commerce. Improving earthquake risk assessment is one of the most important ways of alerting policymakers to the danger earthquakes pose. This requires not only knowledge of the population and building exposure and “hard rock” seismic hazard, but also an understanding of the potential for sedimentary basins to cause amplification and resonance of seismic wave motion. The observation of such effects in recent destructive earthquakes—the 2015 Kathmandu [2] and 2017 Mexico City earthquakes [3]—imparts some urgency to the evaluation of basin effects in other megacities. Such an evaluation can be conducted by profiling the shear-wave velocity (V_s) of sedimentary deposits down to bedrock.

Many approaches have been developed to estimate V_s depth profiles based on geotechnical and geophysical methods. Direct approaches such as borehole drilling, vertical seismic profiling (VSP), and standard penetration testing (SPT) can give accurate information on the velocity profile. However, these are often too expensive to cover the entire area of a large city. Other approaches exist such as seismic refraction or reflection surveys, but these use active sources like explosives and require large, regularly-spaced sensor arrays that are impractical in built-up urban areas. These problems are avoided using passive approaches such as microtremor measurements [4,5,6,7] and interferometry studies [8,9,10,11], which have become very popular for use in densely-populated areas.

Recently, Park and Ishii [12] introduced an alternative approach for estimating “near-surface” V_s by using measurements of body-wave polarization. The term polarization is used here to address the wave’s particle motion (below we explain that “near-surface” V_s is actually apparent half-space V_s, which we denote as V_s^ahs). They utilized the well-known relationships between the V_s and compressional-wave velocity (V_p), and the polarization of body waves interacting with the free surface of a half-space. The incoming polarization directions of P-waves are sensitive to V_s, while those of S-waves are sensitive to both V_p and V_s. Therefore, it is possible to estimate near-surface V_s and V_p by observing the polarization directions of P-waves and S-waves generated by earthquakes at a seismic station. This approach requires no artificial source or other expensive equipment, and it can be applied with a minimum of computational effort.

Park and Ishii [12] concluded that their technique could be used very widely to study near-surface velocity structure wherever three-component seismometers are deployed. Since it is applicable to urban areas, we applied this technique based on P-waves polarization to map shallow V_s beneath Jakarta, the capital city of Indonesia. Because it lies on a thick young sedimentary basin [13,14,15,16], seismic risk in Jakarta is likely to be enhanced by amplification and resonance effects. In late 2013 and again in 2018, two temporary seismic networks consisting of three-component broadband seismometers were deployed for approximately 3 months and 5 months, respectively, covering most of the city and its vicinity. While Park and Ishii [12] presented their application using the Hi-net array in Japan, the combination of seismicity, dense seismic network, and sedimentary basin structure in Jakarta provides a good opportunity to evaluate this technique on a more local scale. The results are benchmarked against borehole data and shallow V_s estimated from other studies.

2. Geologic Setting and Seismometer Deployments

Located on the northern coast of the island of Java, the surface geology of Jakarta mainly consists of alluvial deposits. According to Turkandi et al. [17], the city’s surface geology from the coastline to around 6 km southward to the center of the city consists of Holocene sand dunes and alluvial fan, while deposits of Pleistocene alluvial fan cover the southern part of the city with little trace of Tertiary volcanic deposits. The sedimentary deposits are estimated to thicken northward, with alluvial fan sediments reaching thicknesses of 300 m or more in the city center. Unfortunately, although the geological setting of Indonesia is described generally in Van-Bemmelen [18], there are only a few recent studies of the detailed geological setting of the Jakarta Basin (e.g., Turkandi et al. [17]).

Nevertheless, the existence of a thick sedimentary basin overlying tertiary bedrock strongly suggests that amplification and resonance could enhance the damaging effects of earthquakes. After seismic waves propagate through the transition between stiff rock and soft soil, their energy is trapped within the soil layer by internal reflections [19]. The reverberating seismic waves add constructively at certain frequencies but interfere destructively at others, with the former being resonant frequencies that are determined by the thickness and V_s of the sedimentary layer. Because of resonance, the energy and amplitude of seismic waves, which normally decay rapidly in intensity with distance from the earthquake, can instead increase dramatically in both amplitude and duration. Not accounting for such basin resonances may lead to an underrated earthquake risk assessment.

To better understand such basin effects, three-component broadband seismometers (Trillium Compact) and digitizers built by the Australian National University were installed at various sites (Figure 1) in two separate deployments. Stations spaced at 3–5 km, were deployed temporarily on a concrete slab floor in schools throughout the city. The first deployment comprised 96 stations and operated from October 2013 to February 2014, covering most of the area of Jakarta. Twenty-six stations were maintained for 3 months of recording as semi-permanent stations and the other 26 stations were redeployed in three phases, each of one-month duration. The second deployment, comprising 143 stations was deployed in 2018 and aims to include coverage just outside Jakarta in order to reveal the extent of the basin edge. The 30 stations were maintained and redeployed in five phases with at least one-month duration for every site.

Using data from the first deployment, Saygin et al. [15] extracted Rayleigh wave Green’s functions from cross-correlograms of ambient noise at the different station pairs and imaged the basin structure using Ambient Noise Tomography. They found that the sedimentary basin covers most of the area of the city with a thickness up to 1500 m below central Jakarta [15]. Considering this evidence, they concluded that basin effects in Jakarta are likely to enhance the damaging effects of earthquake-generated seismic waves. This conclusion was also reached by Cipta et al. [16,20], who used the same data to invert Horizontal-to-Vertical Spectral Ratio (HVSR) curves to achieve better resolution of the basement architecture. However, neither of these models revealed the basin edges, which extend outside the city of Jakarta beyond the extent of the 2013–2014 seismometer deployment. For this reason, the 2018 seismometer deployment was undertaken to extend the coverage beyond the city of Jakarta itself, and hopefully resolve the basin edges. In addition, the extent to which either of these studies resolves the very shallow (<100 m depth) V_s structure is unclear.

Both the 2013–2014 and 2018 seismometer deployments were intended to make use of noise interferometry to study the basin structure, especially the first seismic network [15]. Nevertheless, earthquake signals were recorded during both deployments. We have evaluated recorded signals from 56 earthquakes with a good signal-to-noise ratio (SNR), varying from local to regional and teleseismic earthquakes. Applying body-wave polarization analysis to these signals seems feasible and worthwhile, and this is what we report on in this study. However, not all stations record the same events due to the different phases of deployment. The numbers of observed data at every station are summarized by color representation in Figure 1b.

3. Method

3.1. Apparent Incident Angles of Body-Waves

When a body wave arrives at Earth’s free surface, the incident wave is both reflected and converted. In particular, the incident P-wave generates a reflected P-wave and a converted SV wave. This means the particle motion of P-waves recorded by a three-component seismometer on the free surface is determined by the combination of the incoming and the two outgoing waves, which is defined by the apparent incidence angle ($\overline{\theta }$).

Considering a P-wave incident on the free surface of a uniform half-space with P(S)-velocity V_p(V_s), the apparent incidence angle is different from the true incidence angle ($\theta$), with their relationship derived in Wiechert [21] as:

$$p=\frac{\sin \theta }{V_p}=\frac{\sin \left(\frac{1}{2}\overline{\theta }\right)}{V_s}$$

(1)

where p is the ray parameter (or horizontal slowness) of the P-wave.

The derivation of their relationship is also shown by [12] using the free surface boundary conditions. Equation (1) can be rewritten as:

$$V_s=\frac{\sin \left(\frac{1}{2}\overline{\theta }\right)}{p}$$

(2)

$$\overline{\theta }=2\arcsin \left(p{V}_s\right)$$

(3)

which defines the half-space V_s if the ray parameter of the P-wave and the apparent angle are known [22].

In this study, Equation (2) is used to estimate a “near-surface” V_s for the Jakarta Basin. However, since the true V_s profile in the basin is not that of a half-space but has V_s increasing with depth, our estimates of “near-surface” V_s are actually an estimate of apparent half-space V_s, which we denote as V_s^ahs in what follows. V_s^ahs should be representative of the actual V_s averaged over some depth range. This closely follows Svenningsen and Jacobsen [22], who use a different notation, V_s,app to denote the apparent half-space V_s.

3.2. Calculating Polarization

The apparent incidence angles of P-waves are measured from the particle motion of the observed body-waves, with particle motion in the direction of apparent incidence for P-waves. This particle motion can be observed using a “particle motion” plot of the curve connecting particle position in the vertical and radial plane at successive times (see Figures 3c and 4c). Therefore, it is natural to project the recorded three-component seismograms (vertical, north-south, and east-west signal amplitudes) onto the vertical-radial-transversal plane using a priori information of the earthquake source and seismic station location. Herein, the discretized vertical and radial component time series data are defined by the column vectors z and r, respectively.

Polarization of particle motion can be measured in the time domain using principal component analysis (PCA) [23]. For the selected signals, a data covariance matrix is arranged as follows:

$$z={\left[z_1,z_2,\dots ,z_M\right]}^T;r={\left[r_1,r_2,\dots ,r_M\right]}^T$$

$$C=\frac{1}{M}\left[\begin{array}{cc}{z}^{T}z& z^{T}r\\ r^{T}z& r^{T}r\end{array}\right]$$

(4)

where M is the number of data points for each component. The particle motion is then specified by the eigenvectors of this covariance matrix. The eigenvalues $\lambda$ and eigenvectors $v$ of the covariance matrix C are calculated by solving:

$$\left(C-\lambda I\right)v=0$$

(5)

where I is the 2 × 2 identity matrix. From the data comprising vectors z and r, solving Equation (5) yields two eigenvalues ${\lambda }_1$ and ${\lambda }_2$, which have respective eigenvectors $v_1$ and $v_2$. In contrast to Park and Ishii [12], we only utilized eigenvector $v_1$ related to ${\lambda }_1$ (${\lambda }_1$ > ${\lambda }_2$), which defines the maximum energy in the data. The principal polarization of the selected signals in the vertical and radial component is given by the eigenvector $v_1={[v_z,v_r]}^T$.

In estimating the polarization, first of all, the data are selected by windowing to isolate the direct P-wave signal, with window length chosen to include as much of the respective waveform as possible without including coda that is contaminated by arrivals of different wave type or incidence angle (see below). For every window, Equations (6) and (7) are applied to estimate the principal P-wave polarization ${v^p}_1$. Then, the apparent incidence angles are defined by:

$$\overline{\theta }=arctan\left(\frac{{v^p}_r}{{v^p}_z}\right)$$

(6)

3.3. Estimating Apparent Half-Space Velocities

The procedure to invert the observed apparent incidence angles $\overline{\theta }$ for estimates of V_s^ahs is described as follows. For the observed data, Equation (6) describes the observed apparent incidence angle of the ray from a particular earthquake to the seismic station. Equation (3) is used for forward modeling the apparent incidence angles for given V_s^ahs and p. Arranged as an objective function for a single station, we use the misfit modified from [12]:

$$f\left(V_s\right)=\frac{{\sum }_{i=1}^N\left[w_i^p{\left({\overline{\theta }}_i^{obs}-{\overline{\theta }}_i^{cal}\left(V_S^{ahs}\right)\right)}^2\right]}{{\sum }_i^N{w}_i^p}$$

(7)

where superscripts obs and cal denote observed and calculated apparent angles for the i-th earthquake, respectively, and the summation is over N earthquakes. The weighting values $w_i^p$ are given based on the quality of measured data. In this case, we use the total variance in the measured data of particle motion, which is:

$$w=\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}$$

(8)

According to Equation (3), we require the ray parameters of P-waves for a certain station and earthquake geometry. The program Travel Time Toolbox (TTBox) [24] was utilized to compute seismic ray paths and travel times using a 1-D spherical velocity model. Then, ${\overline{\theta }}_i^{cal}$ can be computed as a function of V_s^ahs. A grid-search over V_s^ahs was used to determine the values that minimize the misfit in Equation (7).

4. Application in Jakarta

4.1. Window Selection

After the 56 earthquakes have been identified during the recording, the time information for the P-waves is required to calculate the time windows used for the polarization calculations.

At each station and for each of the 56 earthquakes, the arrival times (or onsets) of P-waves are automatically picked using a kurtosis based algorithm [25]. Unfortunately, this automatic step does not always work well for low SNR signals. Therefore, we manually checked the low SNR data and the P-waves onsets were refined manually to get more reliable times, while no manual re-picking was needed for signals with high SNR.

In order to calculate the polarization using PCA as described in Equations (4) and (5), the signal is selected and windowed to isolate the direct P-wave, starting from its onset. Band-pass filtering can be applied to strengthen the signals, as long as it still preserves the source frequency content. However, a problem arises when choosing the ideal length of the time window. The aim is to include as much of the waveform of the direct P-wave as possible without including contamination from arrivals of different wave type or horizontal slowness, as might be expected, for example, from phases converted at the basement of the sedimentary basin. Such contamination may bias or increase the uncertainty in estimates of seismic velocities. In addition, teleseismic earthquakes will require longer time windows compared to regional and local earthquakes.

Figure 2 shows examples of a regional and a teleseismic earthquake recorded at one of the stations. The band-pass filters of 0.1–2 Hz and of 0.1–1 Hz are used for the regional earthquake and the teleseismic earthquake, respectively. Their P-waves can be clearly distinguished. Obviously, the duration of the teleseismic earthquake is much longer. Then, the calculation of their polarizations using PCA is illustrated by Figure 3 and Figure 4. The onsets are used as the starting point of the time window.

The length of the time window that can be exploited in selecting the signals may vary from very short to significantly long, and they affect the estimation. Depending on the time window length, the result of polarization direction can vary significantly, as shown in Figure 3b and Figure 4b. These differences may result in bias and/or high uncertainty in the estimation of seismic velocities. Choosing a stable and consistent polarization becomes crucial at this point.

A shorter window achieves a purely direct P-wave signal that is uncontaminated by phases of different wave type or incidence angle, but is not necessarily stable owing to its small number of data points. A longer window may contain phases of different wave type or incidence angle, which could bias the result. The synthetic tests (Appendix A) shows that window lengths greater than 1–2 seconds are unreliable due to the contamination of converted waves from the basin basement. Therefore, we analyzed each record to judge whether a mean estimate of P-wave incidence angle could be made, while removing outliers at the beginning and end of the window that might represent instability due to limited data or contamination by converted phases, respectively.

4.2. Apparent Half-Space Velocity Estimates

A station may observe several earthquakes during its recording duration, each of which contributes measurements of P-wave polarizations. These body-wave polarizations are then used to estimate apparent half-space velocities (V_s^ahs) beneath the station. A grid-search optimization was undertaken using Equation (7) for every station, while ray parameters were calculated from the AK135 1-D earth model [26]. Figure 5 shows plots of the apparent incidence angles as a function of ray parameter for station JKA12. Stations that record only a few earthquakes will have poorly constrained apparent half-space velocities, so we only utilized stations that recorded at least four clear earthquakes. Aiming to get the best fit model, V_s^ahs were searched from 50 m/s to 2000 and 4000 m/s using increments of 10 m/s.

The maps of apparent half-space velocities V_s^ahs within Jakarta and its vicinity are shown in Figure 6. Note that stations that lack data were not included in the results. We focus on the V_s^ahs results derived mainly from the relatively high-quality P-wave polarization measurements.

The estimated values of V_s^ahs range from 200 to 800 m/s. Within the boundary of Jakarta, most of the area is dominated by low V_s^ahs between 200 and 400 m/s. The area characterized by low V_s^ahs extends outward beyond the city, only to the southwest. On the other hand, the eastern edge of the city is characterized by higher V_s^ahs, especially outside the city’s eastern boundary. Unexpectedly, given the results of previous studies, higher V_s^ahs can also be observed in the northern part of the city very near to the coastline.

4.3. Comparison and Depth Estimation

The polarization technique appears to provide considerable information on the variation in the shallow structure within the Jakarta basin, as reflected in the measured V_s^ahs. Unfortunately, it is difficult to know what depth range these apparent half-space values correspond to. In their use of a similar technique for measuring apparent half-space V_s from receiver functions, it is noted by Svenningsen and Jacobsen [22] that the results depend on the frequency content of the signals used. They advocate narrow-band filtering centered on inverse period f = 1/T to estimate a curve V_s^ahs(T) that can be inverted for the shear-wave velocity profile. Given the complications that the basin structure may pose for receiver function computation, as well as the use of local and teleseismic events with varying frequency content, we adopted the simpler approach of Park and Ishii [12]. Comparison of their estimates for V_p^ahs and V_s^ahs at stations of Japan’s broadband network Hi-net, allowed Park and Ishii [12] to conclude that their results were representative of the top ~100 m of the V_p and V_s profiles, respectively.

Therefore, in order to assess whether our half-space velocity measurements are representative of shallow V_s structure and if so, to suggest a depth range to which they correspond, we compared our results to previous studies. Ridwan et al. [14] used the spatial autocorrelation (SPAC) method to estimate depth profiles of V_s for depths less than 1 km at 55 sites throughout Jakarta. They provide a good set of benchmark V_s profiles at their site KMAL, located in northwest Jakarta, where they compared their SPAC result with cone penetrometer (SPT) and downhole seismic measurements (Figure 7, modified from Figure 10c in Ridwan et al. [14]). From Ridwan et al. [14], the shear-wave velocity at a depth of 100 m beneath the site KMAL is 350 m/s. Our closest station to this site, which is 700 m away, is JK046, where we used the polarization technique to obtain a very similar estimate for V_s^ahs of 390 m/s. The small difference between these values could reflect either measurement error or be a genuine difference in V_s due to the slight difference in the locations of measurement.

As discussed above, Cipta et al. [16] applied HVSR analysis to the 2013–2104 deployment data to obtain the V_s structure beneath Jakarta. We compared our results with average values of Cipta et al.’s [16] V_s profiles taken over different depths from the surface. We found that our results best match Cipta et al.’s [16] V_s profiles when the latter are averaged over the top 150 m, as shown in Figure 8. Note that the different samples shown are due to the lack of observed earthquakes at some stations in our study.

By comparison to the benchmark site KMAL [14], our V_s^ahs matches the average of the V_s profile over the top 100 m reasonably well. Unfortunately, no data deeper than that is available and one site is insufficient to indicate an overall trend. Meanwhile, our V_s^ahs estimates agree with many of the 150 m depth values of the much wider dataset of average V_s profiles obtained in [16]. A few inconsistencies do exist, but these may be due to anomalies in the measurements such as misalignment or miscalibration of sensors. In any case, our comparisons with previous studies suggests that the V_s^ahs estimates obtained in this study approximately represent average V_s at around the top 150 m depth.

4.4. Correlation with Surface Geology

The distribution of low V_s^ahs agrees with the geologic mapping that shows sedimentary deposits comprise almost all areas of Jakarta’s surface geology. Although the contrast between Holocene and Pleistocene deposits is small and their boundary is not well resolved in our V_s^ahs map, their deposits obviously fill the Jakarta basin to a depth of at least 150 m. These deposits reach the outer part of the Jakarta administrative margin mostly to the southwest, more than 10 km away. The result is in accordance with Cipta et al. [20] who reported the estimation of the basin extension based on extrapolation of the basement depth.

The higher V_s^ahs found in the southern part of the basin may relate to Tertiary volcanic deposits, or even the edges of the basin, which is Tertiary rock. This feature gets more distinct in the eastern part of Jakarta. The gradation from low V_s^ahs to higher V_s^ahs at the eastern administrative boundary underlines the contrast between sedimentary deposits and denser rock. The sedimentary deposits diminish eastward, while the thinnest layer lies at the southeastern corner of our study area. We suggest that the basin edge emerges near the surface in this area.

5. Conclusions

We applied the simple body-wave polarization technique of Park and Ishii [12] to obtain the variation in apparent half-space S-wave velocity (V_s^ahs) over the Jakarta Basin. By measuring the apparent incidence angles of earthquake-generated P-waves using principal component analysis, we obtained estimates of V_s^ahs. Although care had to be taken in choosing window lengths to avoid contamination of the direct P-wave by basement converted phases, we found we were able to obtain stable estimates of polarization by choosing windows of 1–2 sec duration.

The spatial variations we observed in V_s^ahs estimated using the body-wave polarization technique seem sensible when compared to other studies. In particular, when comparing our V_s^ahs with the V_s profiles in [14] and [16], the mapping of V_s^ahs appears to be correlated to the average of V_s profiles over the top 150 meters. In further studies, it might be useful to investigate the frequency-dependence of body-wave polarization in an attempt to reveal further details of the shallow V_s profiles.

Our estimates of V_s^ahs reflect the shallow Vs structure obtained within the Jakarta city limits in earlier studies [15,16], but extended this information beyond the city limits of Jakarta to what is thought to be the basin edge [20]. Although the surface geology of the entire study area is composed of quaternary sediments (Figure 1b), we found that on average, V_s^ahs increases towards the outer edge of the study area (Figure 6). This may indicate that V_s^ahs is sensitive to a reduction in basement depth that indicates the effective edge of the basement. In future studies we hope to obtain a more complete geometry of the Jakarta sedimentary basin that will enable a more accurate ground-motion simulation for hypothetical earthquake scenarios that can characterize seismic risk in Jakarta.