Remarkable ground uplift and reverse fault ruptures for the 2013 Bohol earthquake (Mw 7.1), Philippines, revealed by SAR pixel offset analysis

Abstract

By applying a pixel offset analysis using RADARSAT-2 SAR data to an inland crustal earthquake that occurred on Bohol Island, Philippines on 15 October 2013, we succeeded in mapping a ground displacement associated with the earthquake. The most concentrated crustal deformation with ground displacement exceeding 1 m is located in the northwest part of the island. The crustal deformation is zonally distributed and extends a length of approximately 50 km in the ENE–WSW direction. The ground in the mountainous area moved toward the satellite, while the ground in the northern coastal zone moved away from the satellite. A clear displacement discontinuity with a length of about 5 km, probably corresponding to earthquake surface faults, can be identified in the northeastern region. Our fault model consisting of two rectangular planes shows nearly pure reverse-fault motion on south-southeast-dipping planes with moderate dip angles. A local rupture occurs in the northeast at shallow depths and produces surface ruptures. By applying an additive color process using SAR amplitude images, significant changes in backscatter intensity were detected along the coast from Maribojoc to Loon; these changes suggest that the seafloor uplifted and the shoreline resultantly shifted seaward. The area showing the shoreline change is in good spatial agreement with the locally distributed large ground uplift predicted from our fault model. We identified a good correlation between the ground upheaval produced by the reverse-fault motion and elevation in the mountainous area, which is consistent with the idea that repeated historical reverse faulting developed the present-day topography.

Background

A devastating inland earthquake with moment magnitude (Mw) of 7.1 struck Bohol Island, Philippines on 15 October 2013 [1]. Within the Philippine islands, which straddle a region of complex tectonics, Bohol Island is located in the Sunda block, beneath which the Philippine Sea plate subducts from the east westward. GPS data analyses and focal mechanism solutions suggest that the island is regionally subjected to E–W to NW–SE compression [2, 3]. The Philippines have suffered from considerable historic seismic activity, mainly along the subduction zones and the Philippine fault zone (Figure 1), but the seismicity in and around Bohol Island has been relatively low [4]. The largest previous earthquake that had occurred during the past several decades was the M6.8 event that occurred east of the island in 1990 (Figure 1), and no seismic events with a magnitude exceeding 7 had occurred in the crust. The existence of the East Bohol fault in the southeast part of Bohol Island has been well known, however, the 2013 seismic event did not involve this fault but occurred along a previously undiscovered fault (hereafter called North Bohol fault).

Figure 1

Tectonic setting of Bohol Island, Philippines. Red stars represent the epicenters of the 2013 Bohol earthquake determined by the USGS, the PHIVOLCS, and the Global CMT Project, respectively. Gray and white circles are earthquakes with magnitudes larger than 4 that occurred in the crust after and before the 2013 Bohol event since 1980 (derived from the USGS earthquake archive). A beachball diagram is the Global CMT solution. The blue-colored dotted line is the East Bohol fault. The left panel represents the tectonic setting of the Philippine Islands.

One of the notable features of the seismic event was remarkable ground surface changes, which were observed in field surveys [5]. These changes included the appearance of earthquake surface faults with vertical offsets of several meters and shoreline changes caused by ground uplift. According to reports by the U. S. Geological Survey (USGS) and the Philippine Institute of Volcanology and Seismology (PHIVOLCS), among others, a reverse fault mechanism with a NW–SE compressive axis was inferred from seismic wave analyses, thus vertical ground movements associated with the reverse motion must have been involved in the ground surface changes. However, it remains unclear where and how the fault rupture contributed to these changes. Measurement of ground displacements around the epicentral area certainly plays a key role in answering these questions, but there is no geodetic data from which we could obtain the detailed crustal deformation.

Satellite synthetic aperture radar (SAR) data can provide detailed and spatially comprehensive ground information. Interferometric SAR (InSAR) analysis has the advantage of detecting ground deformation in a vast region with high precision (e.g., [6, 7]). However, for the Bohol event, the standard InSAR approach is not helpful for determining the details of the seismic rupture. Only C– or X-band SAR data are available for the event, thus it is not suitable to apply an InSAR method to measure ground displacement on Bohol Island, which is covered by forest [8–10]. We actually conducted an InSAR analysis using the C-band data we handle in this study, but a coherent loss area resultantly spread over the island. Thus, in order to reveal the unknown surface displacements, we conducted a pixel offset method that enabled us to robustly detect large ground deformation even in incoherent areas [11–13].

Methods

We used RADARSAT-2 data from the ascending orbit acquired on 12 January 2013 and 27 October 2013 for data analysis. This data pair provided the shortest temporal baseline among the SAR images covering the source region. The data obtained were strip-map imagery (Wide Multi-Look Fine mode) with an incidence angle of 34.1° at the scene center. The area analyzed is indicated by the frame in Figure 1. We processed the SAR data from SLC products using a software package Gamma [14]. After conducting coregistration between two images acquired before and after the mainshock, we divided the single-look SAR amplitude images into patches and calculated the offset between corresponding patches by using an intensity tracking method. This method was performed by cross-correlating samples of backscatter intensity of a master image with those of a slave image [15]. A pixel offset is feasible even for such incoherent areas provided that a large correlation window is used to track similar speckle patterns [15]. We employed a nearly square search patch of 256 × 256 pixels (range × azimuth: ~680 m × ~640 m). We then reduced the artifact by applying an elevation-dependent correction incorporating NASA Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data with a 3-arcsec resolution [16] in the same manner as used by Kobayashi et al. [13]. The measured offset consisted of two components: (1) displacement along the line of sight (range offset) and (2) horizontal movement along the ground parallel to the satellite track (azimuth offset).

Results and discussion

Ground displacement field for the 2013 Bohol earthquake

Figure 2 shows the estimated range and azimuth offset fields over the entire analyzed area. The range offset analysis succeeded in mapping the ground displacement in the proximity of the source region (Figure 2a). Intensive deformation was revealed in the northwest part of the island with ground displacement exceeding 1 m, and the deformation extends with a length of approximately 50 km in the ENE–WSW direction. The ground on the south side of the crustal deformation area (warm-colored area) moved toward the satellite, while the ground on the north side (cold-colored area) moved away from the satellite.

Figure 2

Pixel offset fields. (a) Displacement field in range component. Warm and cold colors represent displacements toward and away from the satellite, respectively. Gray-colored frames indicate surface projections of our fault model; thick lines stand for the upper edges. (b) Same as Figure 2a but for azimuth component.

It is notable that a clear displacement discontinuity, across which the ground movement was in opposite direction, can be identified in the northeastern area (Figures 2a and 3a). Field surveys have discovered surface ruptures near the displacement boundary [5]. The length of the displacement boundary was estimated to be approximately 5 km, which is consistent with the field surveys. Figure 4 shows displacement profiles along the cross sections of lines 1–4 (shown in Figure 2). A sharp, large displacement offset suggesting a surface rupture can be seen in the line 3 profile. The observed offset is approximately 2 m, equivalent to a vertical movement of 2.4 m. The clearly identified displacement boundary terminates around the position of 124.13°E, 10.01°N, and large ground movement cannot be clearly seen farther northeast (line 4 in Figure 4a). On the other hand, except in the northeastern part, the displacement changes are relatively gradual (lines 1 and 2 in Figure 4a).

Figure 3

Enlargement views of the offset fields. An enlargement view surrounded by the white frame shown in Figure 2 for range offset field (a) and azimuth offset field (b). A clear displacement discontinuity, across which the ground movement component is opposite, can be clearly recognized.

Figure 4

Displacement profiles. (a) Displacement changes of range offset component along cross section lines 1 to 4 shown in Figure 2 (red lines). A displacement offset can be identified in line 3. Blue solid lines represent displacements calculated from our fault model. Shaded profiles indicate the corresponding topography. (b) Azimuth offset component of line 3.

The azimuth offset field shown in Figure 2b is rather noisy. Crustal deformation produced by reverse-fault motion, in which northward/southward ground movement should be dominant, cannot be identified clearly. However, significant ground movement can be observed near the displacement discontinuity area with a relatively high signal-to-noise ratio. In this area, ground movement opposite to the satellite flight direction (cyan), i.e., nearly southward horizontal displacement, of approximately 0.5 m can be recognized. The southward movement terminates at the displacement discontinuity observed in the range offset (Figure 3b).

Fault modeling

Based on the obtained displacement data, we tried to construct a fault model under the assumption of a rectangular fault with uniform slip in an elastic half-space [17]. A rectangular fault model has the advantage that it can represent a macroscopic feature of the source property with simple notation. A slip distribution model would have provided us with a more detailed picture regarding the fault rupture, but the result of the pixel offset analysis would not necessarily have enough measurement accuracy to estimate a more complex slip distribution. Thus, we here provide only a simple-shape fault model.

The offset field has ground surface changes over a vast range, and so produces too many values to be easily assimilated in a modeling scheme. Thus, to reduce the number of data for the modeling analysis, we downsampled the data beforehand, using a quadtree decomposition method. Essentially, we followed an algorithm presented by Jónsson et al. [18]. For a given quadrant, if, after removing the mean, the residue is greater than a prescribed threshold (25 cm in our case), the quadrant is further divided into four new quadrants. This process is iterated until either each block meets the specified criterion or the quadrant reaches a minimum block size (8 × 8 pixels in our case). Upon application of the above-mentioned procedure, the size of the data set was reduced from 471070 to 405 for the range offset data. We also applied the method in the same manner to the azimuth offset data. The azimuth offset data, however, were too noisy to incorporate into the modeling. Therefore, we used only the data in the northeastern area, ranging from 124.00°E to 124.15°E and from 9.95°N to 10.10°N, where the signal-to-noise ratio was relatively high.

For the modeling, we applied a simulated annealing method for searching the optimal fault parameters [19, 20]. We assumed a two-segment model that consists of a main fault producing the majority of the crustal deformation and a local but essential fault located in the northeast producing clear surface offsets. For the main fault, we randomly assigned parameters within the search range of 123.75°–124.00° in longitude, 9.8°–9.9° in latitude, 0–20 km in depth, 0–70 km in length, 0–30 km in width, 0°–90° (180°–270°) in strike, 0–90° in dip, 0°–180° in rake, and 0.0–10.0 m in slip amount. For the northeast fault, the clear displacement offsets reflecting surface ruptures strongly suggested that the fault rupture is rather shallow, thus we here fixed the fault top to be near the ground surface. With knowledge of the displacement discontinuity line, the search range of strike could be strongly limited to within 55°–65° (235°–245°) so as to fit the boundary line. We randomly assigned the other parameters within the search range of 124.10°–124.15° in longitude, 9.85°–10.01° in latitude, 0–10 km in length, 0–30 km in width, 0–90° in dip, 0°–180° in rake, and 0–10.0 m in slip amount. To estimate the individual confidence of each inferred parameter, we employed a bootstrap method [21].

Figure 5 shows the range and azimuth offsets calculated using our preferred model and the residuals between the observations and calculations. The estimated fault parameters are listed in Table 1. Our model was able to reproduce well the observations (Figures 4 and 5). The root mean square of the residuals in range offset was estimated to be 19.1 cm, which is comparable to or less than the typical measurement accuracy of pixel-offset analysis [13]. Thus, our uniform slip model sufficiently accounts for the nature of the source. Figure 6 shows the three-component displacement calculated from the fault model. From the calculations, we can recognize the remarkable characteristic that vertical displacement, particularly upheaval, is predicted to be dominant for the seismic event. Our fault model shows (1) a south-southeast-dipping fault plane with a dip angle of ~60° for the main fault and ~40° for the northeastern fault, (2) ENE–WSW oriented strikes, and (3) nearly pure reverse-fault motions. The total seismic moment is estimated to be 5.7 × 10¹⁹ N m (Mw 7.1) with a rigidity of 34 GPa. According to the results of the USGS, the PHIVOLCS, and the Global Centroid Moment Tensor Project (GCMT project), the moment magnitude was estimated to be Mw 7.1 (5.6 × 10¹⁹ N m), Mw 7.2, and Mw 7.1 (5.6 × 10¹⁹ N m), respectively [1, 22, 23]. Our result is in good agreement with these.

Figure 5

Range and azimuth offset fields calculated from the constructed fault model. (a) Range offset field calculated from our fault model. Gray-colored frames indicate surface projections of our fault model; thick lines stand for the upper edges. The color scale is the same as Figure 2. (b) Residual between the observations (Figure 2a) and the calculations (Figure 5a). (c) and (d) Same as Figure 5a and 5b but for azimuth offset.

Table 1 Fault parameters of our preferred model for the 2013 Bohol event

Figure 6

Three-component displacement calculated from our fault model. (a) NS component displacement. (b) EW component displacement. (c) UD component displacement. The contour interval is 0.2 m. (d) Enlargement view of UD component for the area corresponding to Figure 7b. A dotted line indicates the coastal zone along which the backscatter intensity increased after the earthquake (Figure 7). The color scale is the same as Figure 6c.

Shoreline changes revealed by SAR analysis

Significant ground upheaval was reported to have occurred in the coastal zone from the Maribojoc to Loon municipalities, and the shoreline shifted seaward by tens of meters ([5]; Toto B. personal communication). Examining the range offset field, relatively large displacements of approximately 1.5 m can be observed in and around the area (Figure 2). The narrow shoreline changes, however, could hardly be extracted by the pixel offset analysis because of the low spatial resolution, thus detailed information about these could not be acquired. To overcome this difficulty, we attempted to obtain the ground surface changes by using the microwave backscatter intensity of SAR amplitude images. For the analysis, we basically followed an additive color process presented by Tobita et al. [24]. After coregisteration between the two SAR images, we first assigned intensity variations in the pre-seismic amplitude image to variations in cyan, (R, G, B) = (0%, 100%, 100%), and then variations in the post-seismic image to red, (R, G, B) = (100%, 0%, 0%). Combining these two images, areas where backscatter increased / decreased / remained unchanged turned out to be red, cyan, and gray, respectively. Figure 7 shows the result of the additive color process. A clear red-colored zone can be observed extending in the coastal area from Maribojoc to Loon, where the backscatter intensity increased significantly. The total distance of the red zone is estimated to be about 13 km. We interpret it to mean that the seafloor uplifted and the shoreline shifted seaward resultantly.

Figure 7

Shoreline changes derived by the additive color process. (a) Difference image of backscatter intensity. Red- and cyan-colored areas represent increased and decreased backscatter intensity after the earthquake, respectively. (b) Enlargement view of Figure 7a. A red-colored zone can be clearly identified along the coast, suggesting that the shoreline shifted seaward.

The major source of error of this method is a tide-level difference between the two SAR data acquisition times. Thus, to confirm the validity of our analysis, we calculated it using Some Programs for Ocean-Tide Loading (SPOTL) software [25]. As a result, the difference in tide level is approximately 9 cm in and around the sea area, thus there would be no serious affect in the analysis result.

Figure 6d shows the model-calculated vertical displacement in and around Maribojoc and Loon, where our fault model predicts a large ground uplift of approximately 1.2 m (at maximum). Compared with other coastal areas, this area is locally subjected to large upheaval. A dotted line indicates the coastal zone showing the increase of backscatter intensity (Figure 7b). The locally distributed large uplift can account for the fact that shoreline changes occurred only in this area.

Relationship between ground movement and topography

There is a good spatial correlation between the range shortening area (warm-colored area) and the mountainous area (Figure 2). In particular, the relationship is obvious in the area from the coast to the north-facing slope of the mountain (lines 1 and 2 in Figure 4). The correlation between the two caused us to think about the landform evolution and whether the spatial correlation could be true. As mentioned in Section Methods, we applied an elevation-dependent correction to reduce the artifact displacement correlating with elevation, but the correction was often insufficient because of the large amount of topographic relief and the accuracy of the DEM data [13]. Thus, we had to consider carefully whether or not the spatial correlation is true. We here confirmed the potential error caused by the topography. For the SAR data pair that we analyzed, the perpendicular baseline of satellite orbit was estimated to be approximately −47 m at the scene center. In this case, the artifact, calculated using a formulation presented by Kobayashi et al. [26], should be at most 8 cm at the elevation of 500 m. Thus, we concluded that the error produced by the topographic effect can be ignored in this analysis.

Figure 8b shows the range offsets as a function of elevation for the analyzed domain indicated by the frame in Figure 8a. There is a good correlation between the range offsets that consist mainly of vertical movement and the elevation, with a correlation coefficient of 0.71. We may thus state that the reverse fault motion is in quite good harmony with the orientation of the long-term cumulative displacement in that the mountainous area is on the side of the hanging wall producing ground uplift. This may support the idea that reverse-fault motions have occurred repeatedly on the North Bohol fault and may have contributed to the development of the present-day topography, although there is no clear evidence on the ground surface indicating it to be an active fault.

Figure 8

Relationship between range offsets and topography. (a) A frame stands for the analyzed area. (b) Range offsets as a function of elevation. The horizontal and vertical axes indicate the elevation and the range offset, respectively.

Conclusions

We applied a pixel offset method using RADARSAT-2 SAR data to the 2013 Bohol earthquake and succeeded in mapping the crustal deformation. The following conclusions were derived from the analyses.

• (1) Intensive deformation with ground displacement exceeding 1 m extends in the northwest part of the island.

• (2) The crustal deformation zone has a length of approximately 50 km in the ENE–WSW direction.

• (3) The ground on the southern side of the crustal deformation area moved toward the satellite, while the ground on the northern side moved away from the satellite.

• (4) A clear displacement discontinuity with a length of about 5 km, probably corresponding to earthquake surface faults observed in field surveys, can be identified in the northeastern part of the source region.

• (5) Our fault model shows nearly pure reverse-fault motions on south-southeast-dipping planes with moderate dip angle. In the northeast, a local fault rupture occurs at shallow depths, causing the appearance of surface ruptures.

• (6) Remarkable changes in backscatter intensity were detected along the coast from Maribojoc to Loon by applying an additive color process, and these changes suggest that the seafloor uplifted and the shoreline resultantly shifted seaward.

• (7) The ground displacements produced by the reverse-fault motion are correlated with elevation, possibly suggesting that reverse-fault motions on the North Bohol fault have repeated historically and have contributed to the development of the present-day topography.