Fault source investigation of the 6 December 2016Mw 6.5 Pidie Jaya, Indonesia, earthquake based on GPS and its implications of the geological survey result

On 6 December 2016 at 22:03 UTC, a devastating magnitude 6-class strike-slip earthquake occurred along an unidentified and unmapped fault in Pidie Jaya, northern Sumatra. We analysed the possible fault using continuous Global Positioning System (GPS) observation available in the region. In our investigation,we searched for the fault source parameters of the northand south-dipping left-lateral faults and the westand east-dipping rightlateral faults. We identified that the fault responsible for the earthquake was located offshore, with a southwestnortheast direction. We also computed the Coulomb failure stress and compared the result with the distribution of the aftershocks. In this study, we demonstrated that the result of the geological field survey conducted soon after the mainshock was attributed to the secondary effects of ground shaking and near-surface deformation, and not surface faulting. The newly identified offshore fault proposed by this study calls for further investigation of the corresponding submarinemorphological attributes in this particular region.


Introduction
Along Sumatra, Indonesia, lies the tectonically active Sumatran fault of 1900-km in length divided into multiple segments with slip rates ranging from 12 to 16 mm/yr [1,2,3]. Previous studies have reported additional fault segments in southern Sumatra based on geological investigations [4], and the Pantan Terong and Celala segments in northern Sumatra based on the analysis of an earthquake event along these segments [5,6].
On 6 December 2016 at 22:03 UTC, a devastating earthquake occurred in Pidie Jaya, northern Sumatra, causing the deaths of 102 people, more than 7,000 damaged houses, and an economic loss of USD 230 million, as reported by the Indonesian National Board for Disaster Management (BNPB). The seismic analyses of the German Research Centre for Geoscience (GFZ), the United States Geological Survey (USGS), the Global Centroid Moment Tensor (GCMT) [7,8], and the Indonesian Agency for Meteorology, Climatology, and Geophysics (BMKG) [9] have indicated that the earthquake moment magnitude was Mw 6.5-6.6, with a depth between 10.0 and 17.5 km. Interestingly, the focal mechanism of the 2016 Pidie Jaya earthquake (PJE) was a strike-slip mechanism, located ∼50 km east of the currently well-identified Sumatran fault ( Figure 1). As the 2016 PJE occurred at an unidentified and unmapped fault, the fault nodal plane ruptured by the mainshock is debatable.
Using information from the CMT catalogue, Muzli et al. [10] used one-month aftershock data starting from eight days after the mainshock, hereinafter referred to as the later aftershocks, to constrain the coseismic slip distribution of the 2016 PJE. The distributed one-month aftershocks indicate a left-lateral rupture, with the fault trace from the southwest-northeast direction, called the Panteraja fault. In their analysis, unfortunately, the right-lateral fault as another possible fault source was not investigated. The use of aftershock locations from eight days after the Figure 1: Tectonic setting of the study area. Aftershocks from 6 to 13 December 2016 denoted by the yellow circles are from Supendi et al. [9], while aftershocks from 14 December 2016 to 14 January 2017 indicated by the green circles are from Muzli et al. [10]. The segments of the Sumatran fault are from Sieh and Natawidjaja [1] and Gunawan et al. [5]. (a) The epicentre of the mainshock showed by a star symbol, with the dark green colour indicating the location from the GFZ catalogue; red colour from the USGS catalogue; dark brown from the GCMT catalogue; and blue from the BMKG catalogue. Cyan squares indicate the location of the GPS stations. (b) The epicentre difference between the two catalogues of Muzli et al. [10] and Supendi et al. [9] from [14][15][16][17][18][19] December 2016. The inset shows the global location of our study with the blue circle indicating the location of the IGS stations used in the analysis. mainshock to define the mainshock fault source location may be considered as misleading information because the static stress changes from one rupture could trigger subsequent ruptures [11,12,13,14,15]. Moreover, Supendi et al. [9] showed that one-week aftershocks just after the mainshock, hereinafter referred to as early aftershocks, are estimated based on the BMKG seismic network, have the tendency of the northwest-southeast direction ( Figure 1).
The motivation of this study was to identify the fault source of the 2016 PJE based on the available GPS (Global Positioning System) observation. To the best of our knowledge, thus far, this GPS observations have not been used for analysing the fault slip of the 2016 PJE. We will introduce the aftershock, coseismic displacement, and fault slip modelling in Section 2. In Section 3, we will show the results of the fault slip and our investigation using the Coulomb failure stress. Section 4 will discuss the preferred fault model and the implications of the geological survey result. Finally, Section 5 will summarize our findings.

Aftershock
Supendi et al. [9] reported aftershocks of the 2016 PJE from 6-19 December 2016. Their relocated aftershocks were derived after relocating the hypocentre with the doubledifference technique based on the routine national seismograph network data. On the other hand, Muzli et al. [10] deployed a temporary seismic network for a one-month duration starting from eight days after the mainshock. Consequently, those two catalogues should have different accuracies and different biases. To evaluate the difference between those two catalogues, we compared the location difference of the recorded hypocentre. Between the two catalogues, the overlap period was from 14-19 December 2016. Figure 1 shows the epicentre difference between the two catalogues.
Using a denser seismic network located near the mainshock, the catalogue from Muzli et al. [10] recorded more aftershocks than those recorded by Supendi et al. [9]. While Supendi et al. [9] recorded seven aftershocks between the 14 th and the 19 th of December 2016, Muzli et al. [10] recorded 76 aftershocks. From those seven aftershocks recorded by the two catalogues, in general, the distance differences were ∼4 km, except that two aftershocks had a distance of ∼10 km and ∼20 km.

GPS
We estimated the coseismic displacement of the 2016 PJE by using three available continuous GPS observations in northern Sumatra ( Figure 1). These data were obtained from the GPS stations maintained by the Geospatial Information Agency of Indonesia (BIG) and the Ministry of Agriculture and Spatial Planning/Land Administration Agency (ATR/BPN). In our GPS data processing, we used GAMIT/GLOBK to derive the displacement time series of each of the GPS stations [16]. During the process, we used 20 IGS (International GNSS Service) stations, i. e. ALIC, BAN2, COCO, CUSV, CNMR, DARW, DGAR, GUAM, HYDE, IISC, KARR, KOUC, LAE1, NTUS, PIMO, TOW2, TNML, WUHN, XMIS, YARR, and tied the local network into ITRF2008 reference frame [17].
A three-step approach was used during data processing. First, daily positions with atmospherically used, loose-constraint, prior GPS phase observations with fixed orbit and earth orientation parameters were estimated. Second, we combined the positions and the covariance with GPS solutions computed as part of MIT's processing for the IGS. Then, we examined antenna changes. Third, we estimated daily solutions at each GPS station and subtracted the velocity of three days after to five days before the mainshock. We used the result as the coseismic displacements associated with the 2016 PJE at each GPS station. Detail of GPS data processing can also be found in Gunawan et al. [18]. We notice that the largest observed coseismic displacement was ∼2.5 cm for the GPS station located ∼32 km northwest from the mainshock centroid, called PDIE. In contrast, the far-field GPS stations ∼100 km away from the epicentre showed considerably smaller displacements of 0.1 cm at ACEH and 0.3 cm at LHMI.
Using the horizontal and vertical displacements of GPS data, we conducted a fault slip inversion to determine the fault responsible for the 2016 PJE. As only one GPS station detected significant coseismic displacements, i. e. PDIE, we performed a linear least-square inversion on a single fault to avoid large uncertainties in the slip result. The fault slip was estimated using Gs = d, where G is Green's function, s is the fault slip, and d is the coseismic displacement from the GPS observation [19]. Green's function was calculated based on an elastic half-space model following Okada [20]. Then, the fault slip was estimated as follows:

Fault source modelling
The focal mechanism of the 2016 PJE showed two possible fault orientations, i. e. the right-lateral fault with the northwest-southeast direction, and the left-lateral fault with the southwest-northeast direction ( Figure 1). In this study, we searched for the fault responsible for the 2016 PJE by estimating the possibilities of four fault sources. The first and the second fault models (model 1 and model 2) were constructed following the north-and southdipping left-lateral faults, whereas the third and the fourth fault models (model 3 and model 4) were constructed following the west-and east-dipping right-lateral faults.
In our estimation, we fixed the fault strike, dip angle, length, and width, to reduce unknown parameters. The fault strikes of model 1 and model 2 were fixed at 230°and 50°, respectively, following the aftershock distribution of Muzli et al. [10]. In contrast, model 3 and model 4 incorporated faults with a strike of 140°and 320°based on the relocated aftershock of Supendi et al. [9] and the identified cracks from the geological survey (see Section 4.2). For both model 1 and model 2, we set the fault top depth of 1 km and a dip angle of 80°following the seismic analyses (GFZ, USGS, GCMT, and BMKG catalogues), which indicated that the fault was not a purely vertical strike-slip fault.
To estimate the length of the fault, we used an empirical relationship between the magnitude and the fault length and width [21]. For a strike-slip fault categorised as a magnitude 6.5 earthquake, the relationship was as follows: M W = a + b log(A), with a = 3.99, b = 1.000, and A is the fault area. Considering the earthquake depth that varied between 10.0 and 17.5 km (GFZ, USGS, GCMT, and BMKG), we set the fault width to 15 km, and the resulting fault length was 20 km. In our fault slip estimation, we searched for two unknown parameters, i. e. rake and slip.
In addition to these two parameters, we searched for the longitude and the latitude of the fault location. Because these two parameters were non-linear, we used a grid search algorithm to search for the optimum fault location for models 1-2. Meanwhile, the fault location of fault models 3-4 were defined following identified surface offsets during the geological field survey (see Section 4.2). For each model, the misfit was calculated based on the displacement data and model as follows: where n is the number of data items. The optimum reference point of the fault location for fault models 1-2, i. e. the upper-left corner of the faults, is shown in Figure 2.
The minimum misfits for fault models 1-4 were obtained when the fault source was located offshore (Figure 3). Our analysis showed that fault model 1 incorporated a slip of 38.3 cm with a rake of 39.1°; the slip and the rake of fault model 2 were 58.1 cm and −73.3°; the slip and the rake of fault model 3 were 80.5 cm and 94.8°, and the slip and the rake of fault model 4 were 47.7 cm and 151.1°, respectively.

Coulomb failure stress
The available aftershocks data as reported by Supendi et al. [9] and Muzli et al. [10] provided important information for the investigation of the fault source. Previous studies have suggested that aftershocks can be explained by the Coulomb failure stress (ΔCFF) analysis of the mainshock [22,23,24]. The Coulomb failure stress, which quantitatively describes the static stress changes, is defined as follows: ΔCFF = Δτ + μ ὔ Δσ, where Δτ is the shear stress change on a given fault plane (positive in the direction of the receiver fault slip), Δσ is the normal stress change (positive for fault unclamping), and μ ὔ is the effective friction coefficient for the receiver fault [25]. A positive ΔCFF suggests that the stress in the region is acting to promote the slip, while a negative ΔCFF suggests preventing the slip.
Previous laboratory studies have suggested an effective friction coefficient, μ ὔ , of between 0.60 and 0.85 [26]. Cattin et al. [27] suggested that μ ὔ > 0.6 might explain the observations following the stress change analysis after the 2004 Sumatra Andaman earthquake and the 2005 Nias earthquake in the region. Furthermore, Gunawan et al. [5] showed that the use of μ ὔ = 0.8 could explain the aftershocks following the 2 July 2013 Mw 6.1 Aceh earthquake in northern Sumatra, Indonesia. Taking into account the fact that the 2016 PJE was located in the same tectonic region as the regions considered in the above mentioned studies, we used μ ὔ = 0.8 for further analysis.
While the early aftershocks showed a tendency of the northwest-southeast direction [9], the later aftershocks indicated a southwest-northeast direction [10]. This might suggest that two fault planes were acting as the failure. In this study, we used a combination of these two possible fault mechanisms to construct the receiver faults. For the first receiver fault, the strike is 230°, the dip is 80°, and the rake is 0°. Meanwhile, for the second receiver fault, the strike is 140°, the dip is 80°, and the rake is 180°. The ΔCFF result of the mainshock for fault models 1-4 at a depth of 5 km is shown in Figure 3. The ΔCFF result shows that, in general, most of the aftershocks are evident in a region of positive values.

Preference of the fault model
In our fault models 1-4 estimation, we showed that the minimum misfit between displacement data and the model was obtained when the fault source was offshore. Because of limited GPS data and the uncertainties of aftershock distribution, it is difficult to conclude which model is the best. Looking in detail at the slip inversion result, the displacement model for fault model 3 failed to fit the data at the LHMI station. At this station, though the GPS data showed an eastward direction, the model showed a westward direction.
Among the three remaining preferred models, the displacement model at the closest and most impacted station, i. e. the PDIE station, exceeds the displacement error and worsens the fit for fault model 4 than fault models 1 and 2. Thus, our preferred candidates are fault models 1 and 2, where the fault has a southwest-northeast direction, and misfits are 1 mm. For these two models, the GPS data inversion cannot constrain the dipping direction of the fault. Our preferred fault model with a southwestnortheast direction is similar to the result of the seismic analysis [10,28], and InSAR analysis [29]. Our result emphasizes that the fault source was located offshore and the mainshock transferred stress to the surrounding fault, which led to noticeable inland aftershocks (Figure 3).
Assuming a rigidity value of 40 GPa, the fault models 1 and 2 yielded a geodetic moment of 4.60 × 10 18 N⋅m and 6.97 × 10 18 N⋅m, which was equivalent to Mw 6.4 and Mw 6.5. This suggested that the geodetic moment of fault model 1 was slightly lower than the seismic moment of Mw 6.5, as shown by the USGS catalogue. On the other hand, a similar result is found between fault model 2 and the USGS catalogue. If we are to assume that the earthquake moment from seismic and geodetic analyses are identical [19], then our preference would be the south-dipping left-lateral of fault model 2. In addition, the south-dipping fault was also proposed by previous studies [10,29].
A comparison between fault models 1 and 2 revealed that the ratios of aftershocks occurring at the positive ΔCFF values for fault model 1 and model 2 were ∼68 % and ∼56 %, respectively. These average percentages are comparable to those of other earthquake cases, such as the 2008 Wenchuan earthquake [30], the 2013 Lushan earthquake [31], and the 2018 Palu-Donggala earthquake [32].

Implications of the geological survey result
Soon after the 2016 PJE, a team from the Geological Agency of Indonesia conducted a geological field survey from 7 to 21 December 2016. The purpose of this survey was to identify the surface offset of the 2016 PJE. During the survey, the team did not find the surface offset, which could be associated with the southwest-northeast direction of the Panteraja fault. Instead, they identified cracks that are possibly associated with the right-lateral offset of the mainshock. The measured cracks were varied between 2 to 3 cm ( Figure 4). Reports of the observed geological field survey were then published by the National Center for Earthquake Studies of Indonesia [33]. The three observations of cracked roads as shown in Figure 4 are of cm-offset cracks with an opening component, which cut orthogonally across east-west oriented asphalt roads. In any case, these features cannot be explained as fault trace ruptures of a northwest-southeastoriented fault. The coastal areas of Pidie Jaya were affected by enhanced ground shaking due to shallow seismic amplification along the coastal plain setting with low Vs30. Thus, these cracks are attributed to secondary effects of ground shaking and near-surface deformation, and not surface faulting. It is very common for rigid asphalt roads to crack locally in response to more general ground deformation of the substrate. Previous studies have similarly suggested that the damage caused by the earthquake was due to the effect of ground-motion amplification caused by greater sediment thickness along the coastline [10,34].
The offshore fault proposed by this study was not documented in the national document of the Indonesia Earthquake Source and Hazard Map 2017 [35]. This new information, thus, would be very useful for updating the next national seismic hazard maps of Indonesia. Furthermore, it also amends an updated published national report based on the geological field survey. Finally, the newly identified offshore fault proposed by this study calls for further investigation of the corresponding submarine morphological attributes in this particular region. This is crucial, especially with the subsequent disasters that occur after earthquakes, such as landslides and tsunamis, which may have a broader impact on the regions near the west coast of Peninsular Malaysia.

Conclusion
In this study, we investigated the possible fault source of the 6 December 2016 Mw 6.5 Pidie Jaya, Indonesia earthquake using the available GPS observations. Our GPS analysis showed that the minimum misfit was obtained when the fault source was located offshore, with a southwestnortheast direction. We note, however, that by assuming the earthquake moment derived from seismic and geodetic analyses are identical, our preference would be the southdipping left-lateral fault. Our investigation on Coulomb failure stress suggested that aftershocks occurring at the positive ΔCFF values for the interpreted fault models are comparable to other earthquake cases. Finally, we showed that our findings amend the updated published national report based on the geological field survey and we call for further investigation of the corresponding submarine morphological attributes in this particular region.