Geometric Variation in the Surface Rupture of the 2018 Mw7.5 Palu Earthquake from Subpixel Optical Image Correlation

: We obtained high-resolution (10 m) horizontal displacement ﬁelds from pre- and post-seismic Sentinel-2 optical images of the 2018 Mw7.5 Palu earthquake using subpixel image correlation. From these, we calculated the curl, divergence, and shear strain ﬁelds from the north-south (NS) and east-west (EW) displacement ﬁelds. Our results show that the surface rupture produced by the event was distributed within the Sulawesi neck (0.0974–0.6632 ◦ S) and Palu basin (0.8835–1.4206 ◦ S), and had a variable strike of 313.0–355.2 ◦ and strike slip of 2.00–6.62 m. The NS and EW displacement ﬁelds within the Palu basin included ﬁne-scale displacements in both the near- and far-fault, the deformation patterns included a small restraining bend (localized shortening), a distributed rupture zone, and a major releasing bend (net extension) from the curl, divergence, and shear strain. Surface rupture was dominated by left-lateral strike-slip from initiation to termination, with a localized normal slip component peaking at ~3.75 m. The characteristics and geometric variation of the ruptured fault controlled both the formation of these surface deformation patterns and sustained supershear rupture.


Introduction
Over the past three decades, our ability to measure ground displacements produced by tectonic processes (e.g., earthquakes, landslides, volcanic activity, etc.) has improved dramatically due to improvements in space geodesy [1][2][3]. Using space geodesy measurements (i.e., global navigation satellite system (GNSS), interferometric synthetic aperture radar (InSAR), and satellite images)), we can restructure the ground displacement produced by earthquakes, and further reveal the geometrical parameters of the ruptured faults, coseismic deformation patterns, and even deep structural features [4][5][6]. However, GNSS can only provide sparse coseismic displacements and only in regions with a monitoring station [7]. In contrast, InSAR can measure all or most coseismic ground displacement; however, it cannot provide the near-fault displacement of ruptured faults due to incoherence of the interference phase caused by large amplitude displacement and/or intense ground shaking ( Figure 1b) [5,8]. Much significant information is contained within near-fault displacement, including ruptured fault traces, characteristics of coseismic deformation, and the shallow surface slip model [9][10][11]. Subpixel optical image correlation [7] provides an opportunity to reconstruct coseismic horizontal displacement fields from high-resolution pre-and post-seismic optical satellite images, which is required to measure horizontal coseismic displacement in both the near-fault (within hundreds of meters from the ruptured fault traces) and far-fault (> 1 km away from the ruptured faults) regions. In theory, the subpixel correlation method using a pair of optical satellite images (before and after the event) is able to identity 0.1 pixel of ground displacement with an accuracy of 0.05 pixel (subpixel detection capability); therefore, when using optical satellite images with a ground pixel resolution of 0.5-10 m, the smallest displacements that may be measured are 0.05-1.00 m [12,13]. Hence, sub-pixel optical image correlation can be used to reconstruct detailed near-and far-fault coseismic ground displacement, and allows us to quantify complex surface ruptures using a precise and convenient method. A number of previous studies have utilized subpixel optical image correlation technology to process multi-source optical image data, such as SPOT-5 imagery [14], SPOT-7 imagery [15], Pleiades imagery [11], Landsat-8 imagery, Worldview imagery [16], and Sentinel-2 imagery [10], in order to reconstruct coseismic horizontal displacement fields. Horizontal displacement fields provide valuable information on near-fault deformation patterns [10,17] and fault rupture characteristics, and provide robust constraints on the shallowest part of fault slip models [18]. Documented examples of complex coseismic surface ruptures using subpixel optical image correlation include the 2013 Mw7.7 Balochistan earthquake [11,15,16], the 2016 Mw7.8 New Zealand earthquake [18], the 2018 Mw7.5 Palu earthquake [10,19,20], and the 2019 Mw7.1 Ridgecrest earthquake [21].
At 10:02:43 (Coordinated Universal Time, UTC) on 28 September 2018, an Mw7.5 earthquake event caused by a strike slip fault struck the city of Palu in Sulawesi, Indonesia [22]; a few minutes later, a 4-7 m-high tsunami was produced. The earthquake was also responsible for the generation of multiple landslides around eastern Palu city, triggered by the strong ground motion in Palu basin (Figure 1b), where the main fault ruptured [23,24]. A number of studies have investigated the dynamic rupture process of this earthquake, including surface rupture kinematics, the complex multiple-fault model, and how a tsunami and multiple landslides could have been produced by strike-slip event. Bao et al. [25] and Fang et al. [26] presented strong evidence of early and persistent supershear rupture at sub-Eshelby speed (propagated at a sustained velocity of 4.1 km/s) using seismological far-field Rayleigh Mach waves. Socquet et al. [10] identified supershear rupture characteristics on the southern Palu coastlines using space geodetic imagery of the Advanced Land Observing Satellite (ALOS)-2 and Sentinel-2 satellite. Slip partitioning of the shallow surface was documented by Song et al. [27], who found that a left-lateral strike-slip rupture controlled slip on the onshore fault, while slip on the offshore fault was dominated by an east-dipping normal fault, these observations explain how this strike-slip earthquake was able to produce a tsunami.
The Palu fault is considered to make contributions to the region crustal deformation and mass lateral extension in the area of Sulawesi (Figure 1). Although the triple junction zone where the Philippine Sea, Australian and Sunda plates meet (Figure 1a) is highly seismically active, the Palu fault has a low level of seismicity. Obviously, the seismicity of these active faults developing eastern of the Palu basin is much higher than the Palu fault (Figure 1b). Previous studies have documented that the Palu fault, with a left-lateral strike slip rate of 39 mm/y and an extension rate of 11-14 mm/y, develop along the bedrock mountain. However, the rupture trace of the earthquake documented by Socquet et al., indicated that the initial rupture trace within the Palu basin did not reproduce the Palu fault. Hence, the 2018 Palu earthquake provides a good opportunity to investigate the structural features of the Palu fault, and factors controlling surface rupture [28][29][30][31].
In this study, we first reconstructed the coseismic horizontal displacement fields for both the north-south (NS) and east-west (EW) components using Sentinel-2 optical images, we then computed the curl, divergence, and shear strain of both components to quantitatively analyze the geometric complexity of the surface rupture produced by the 2018 Mw7.5 Palu earthquake. An additional aim of the study was to demonstrate the use of a simple method (i.e., horizontal displacement fields, curl, divergence, and shear strain) to investigate complex surface deformation patterns.

Data Processing
We used Sentinel-2 optical images (near-infrared band 8, spatial resolution of 10 m) and the offset tracking method to reconstruct the horizontal displacement field of the earthquake. First, the pre-(17 September 2020) and post-earthquake (2 October 2018) optical images were processed by orthorectification and co-registration using the COSI-corr software (Co-registration of Optically Sensed Images and Correlation) [32,33]. We then calculated the EW and NS displacement components of the coseismic horizontal displacement field using the MicMac package based on subpixel image correlation technology. The MicMac package, which was developed by the French National Geographic Institute, can be used to compute the ground displacement field with an accuracy of 0.1 pixel [15]. When utilizing MicMac to compute ground displacement, a non-linear cost is used to limit the impact of noise on displacement measurements. When dealing with optical images with significant inhomogeneity, the processing step becomes particularly important [34]. The non-linear cost can be computed as: where Cor is the normalized cross-correlation coefficient, C min is the correlation threshold (below this value, the correlation has no influence), and γ determines the influence of the correlation scores. Detailed correlation parameter settings are shown in Table 1 [34][35][36]. We then identified the influence of clouds in the coseismic area, and improved the quantity of coseismic displacements using a normalized cross-correlation coefficient (NCCC; Figure 2c), only those NS ( Figure 2a) and EW ( Figure 2b) displacements with an NCCC greater than 0.7 were retained to ensure that our displacement fields had high correlation (i.e., to avoid the interference of image noise and clouds on the measurement results). The high-correlation NS and EW displacements indicate that the surface rupture trace (dotted line in Figure 2a) produced by the 2018 Palu Mw7.5 earthquake is consistent with the trace of Palu fault to the south of Palu city (0.8835 • S), the location of which is well-documented based on investigations using tectonic geomorphology and geodetic data. When propagating offshore of Palu Bay, the rupture trace disappeared (length of~28 km) from the displacement fields due to a loss of correlation at sea. The rupture, which had a length of 75 km, mean strike of N173 • , and surface slip of~3 m, reappeared further north within the Sulawesi neck (0.0974-0.6632 • S). Compared with the Sulawesi neck, coseismic displacement fields within the Palu basin (0.8835-1.4206 • S) showed finer-scale near-and far-fault displacements. The rupture in this region was almost linear with a strike near south-north, and formed a number of interesting rupture patterns, including a major releasing bend (1.1866-1.2303 • S).

Curl, Divergence, and Shear Strain
The high-correlation displacement fields (NS and EW component) within the Palu basin provide an excellent opportunity to investigate the geometric variation controlling the surface rupture and deformation patterns; however, it is challenging to identify the characteristics and geometric variation of surface rupture when only utilizing horizontal displacement fields. Hence, the fault parallel ( Figure 3c) and perpendicular ( Figure 3d) displacement fields were calculated by combining the NS and EW displacements (Figure 3a,b) as follows: where D para and D perp represent the fault parallel and perpendicular displacement fields, respectively, U N and U E are the NS and EW displacements, respectively, and θ is the fault strike. We further extracted the curl (ω), divergence (δ), and shear strain (γ) from the NS and EW displacement components to precisely extract the surface rupture trace, and to analyze the geometric variation controlling the surface rupture and surface deformation patterns. Curl (ω), as a means of visualizing and analyzing surface ruptures, delineates fault traces clearly and conforms to the right-hand rule (left-lateral strike-slip corresponds to anticlockwise rotation and hence positive ω at the fault) [11]. Divergence (δ) can be considered as an operator to qualitatively analyze the uplift areas generated by vertical motion. Positive δ at the fault is consistent with the opening of cracks in the uplifted region (net extension) [37]. A shear zone can be identified from the shear strain (γ), which measures the variation of the angle between NS and EW displacements, where a rupture generated by pure shear stresses corresponds to an angle of 0 • (left-lateral and hence positive γ at the fault) or 180 • (right-lateral and hence negative γ at the fault); otherwise, the rupture is characterized by both shear stresses and tensional or compressional components [38,39]. Curl (ω), divergence (δ), and shear strain (γ) are computed as follows: where → F is the displacement field, u and v are the EW and NS displacement components from Sentinel-2 optical images, respectively, x is the east-west direction (E), and y is the north-south direction (N).
The rupture within the Palu basin extended southwards from the Palu coastline for~65 km (Table 2) Figure 4b) is consistent with the negative divergence (δ < 0) at the fault (unlike the mean divergence of segment E, which is close to zero). Minor normal slip may reflect mudslides and/or liquefaction occurring west of fault segments A and B [23]. In summary, segments A, B, and E experienced pure left-lateral strike-slip rupture.  From the divergence at the fault, we also identified geometric variation that controlled the surface rupture. The rupture moved westward for 0.5 km at the end of segment B, and then propagated south along the base of some bedrock mountains (segments C and D). Negative divergence at the termination of segment B transformed to positive divergence (i.e., overall extension) for segments C and D, as expected from the normal slip of the two segments. The majority of the distributed rupture zone (~1.085 • S, length of 3.6 km, mean width of 0.5 km) formed within segment C due to geometric variations of this segment. Interestingly, the rupture changed abruptly at the start of segment D, bending sharply southwest for 7 km with strike of~313.0 • , where it formed a releasing bend. The mean shear strain value (γ mean = 0.0058) at the releasing bend is much lower than the mean divergence (δ mean = 0.0372), which is intuitive given that the stresses dominating the segment rupture were most likely tensional [39,40].

Localized Shortening
Our observations within the Palu basin show geometric variation controlling the surface rupture. An example of localized shortening was mapped using the high-resolution displacement, curl, divergence, and shear strain ( Figure 5). We obtained a near-fault strike-slip of~4.5 m (Figure 5a,b) and localized shortening of~0.8 m (Figure 5c,d) from across-fault profile AB, the peak curl value (large positive value and hence left-lateral strike-slip) coincides with the main rupture (blue stripe in Figure 5f). In contrast, the small negative divergence of the main rupture indicates localized shortening (black polygon in Figure 5a) due to thrust faulting (upwards movement of the east-side of the main fault), we further explain the localized shortening as a restraining bend, where the rupture bent slightly west forming a barrier to the northward rupture and further producing the localized thrust component (at least at the surface). The rupture on the main fault was dominated by left-lateral strike-slip, as demonstrated in the shear strain profile of the main fault (maximum value of~0.02). In summary, as demonstrated in Figure 5, rupture on this segment of the fault was characterized by left-lateral strike-slip and localized thrust faulting (Figure 5k).

Distributed Rupture Zone
Another example of geometric variation controlling surface rupture is shown in Figure 6. The original structural characteristics (straight without major structure complexity) of fault segments A and B became curved to maintain the shear rupture (strike-slip of~4.7 m, Figure 6b), and this change in fault geometry may have caused a distributed rupture zone with a length of 3.5 km and width of 500-700 m in the shallow surface (blue stripe in Figure 6). The main and secondary ruptures can be intuitively identified from the curl (Figure 6f) and shear strain profiles (Figure 6j), where they reach peaks on the fault. Positive divergence (Figure 6h) on fault-perpendicular displacement profile CD (Figure 6d) for the main rupture corresponds to~0.6 m of net extension. The left-lateral strike-slip and localized normal faulting in this area (Figure 6k) may reflect the structural characteristics of the fault segment under the continuous concentration of near north-south strain.

Releasing Bend
As discussed, the southwest bend in the rupture (to a strike of~313.0 • ) formed a releasing bend of 7 km in length at the junction of the southern Palu basin and the bedrock mountain (1.1866-1.2303 • S). From across-fault profile EF, strike-slip was reduced to~3 m (Figure 7a,b), while normal slip increased to~1.5 m (Figure 7c,d), in accordance with the peaking curl ( Figure 7d) and increasing divergence value (Figure 7h). The maximum shear strain was less than~0.0175 (Figure 7i,j), indicating that the shear rupture was weakened by the change in the fault geometry. Under the combined action of left-lateral strike-slip and normal slip, a main subsidence area appeared to the north of the rupture within the southern Palu basin.

Comparison with Published Results
With the help of the Sentinel-2 horizontal displacements and curl, divergence and shear strain maps, we found that the rupture produced by the 2018 Mw7.5 Palu earthquake was dominated by left-lateral strike-slip with thrust and extension components, this is highly consistent with the characteristics of surface slip reported by Socquet et al. [10] and Wang et al. [41]. More importantly, the divergence values extracted at the main rupture (Figure 4e) change from negative (shortening on segments A and B) to positive (extension on segments C and D) near 1.08 • S, which corresponds to fault-normal slip profiles along strike (by latitude) measured by Socquet et al. [10], which indicates the divergence derived from the NS and EW displacements can be considered as a precise and simple method to identify the extensional or compressional components. These pure strike-slip earthquakes can produce the extensional or compressional components (at least at the shallow surface) due to the geometric variations of the surface rupture. In addition, the rupture trace within the southern Palu basin was precisely identified from curl and shear strain maps and was in accordance with the results of Socquet et al. [10] and Song et al. [27].
Feng et al. [20] reconstructed the 3D (three dimensional) deformation field of the 2018 Mw7.5 Palu earthquake using NS and EW displacements from Planet images, with the azimuth and range of deformation from ALOS/Phased Array type L-band Synthetic Aperture Radar (PALSAR) synthetic aperture radar (SAR) data. Their vertical deformation field indicated two obvious subsidence areas to the east of the main fault caused by east-dipping normal faults, this is consistent with the geometric variation of segments C and D (a main releasing bend) and with the overall positive divergence values obtained in this study. In addition, from the curl, divergence, and shear strain maps, we identified a distributed rupture zone in segment C (near 1.1 • S). This has not been described in previously published results possibly because the deformation patterns (i.e., cracks, uplift and subsidence, rupture zone) are too small to be visually observed in NS and EW displacement maps. However, curl, divergence, and shear strain obtained from NS and EW displacement fields highlight these features in detail [11,[37][38][39][40].

Supershear Rupture during the 2018 Mw7.5 Palu Earthquake
Bao et al. [25] used a slowness-enhanced back-projection (SEBP) of data from the Australian seismic network to image spatiotemporal characteristics of kinematic rupture processes, taking the 2018 Mw7.5 Palu earthquake as a supershear event. The rupture propagated at a sustained velocity of 4.1 km/s from its initiation to its end, Fang et al. [26] and Ulrich et al. [42] obtained similar results. Previous investigations of supershear ruptures, including the 1999 Izimit [43] and 2001 Kunlun [44] earthquakes, indicate that they rupture at the bedrock fault or within a few kilometers, have straight traces without structural complexities, and have smooth coseismic slip in pure mode II [45]. Previous investigations also indicated the supershear rupture usually occurs in strike-slip earthquakes with unilateral rupture (coseismic slip of hanging wall is much larger than that of footwall). The characteristic appears in the 2018 Palu earthquake again. The coseismic slip of hanging wall (east side of surface rupture trace, Figures 3a and 4a) is produced by the overall left-lateral strike-slip event peak at~5 m and the west side (footwall) is only 2 m. Hence, the Palu earthquake is considered as a strong piece of evidence to demonstrate the characteristics (strike-slip event and unilateral rupture) of the supershear rupture. These characteristics were also observed in the 1999 Izimit [39] and 2001 Kunlun [40] earthquakes. The 1999 Izimit earthquake ruptured a length of 150 km along the northern Anatolian fault, and what is more interesting is that the rupture in the west and east fault segments propagated by 3.0 km/s, but the middle fault segment (length of 50 km) ruptured by 4.8 km/s (supershear velocity) [46,47]. More importantly, the supershear rupture in the middle fault segment was dominated by strike-slip motion and the surface trace is described in the field or seen on satellite images as remarkably linear, continuous and narrow (clearly associated with mode II rupture) [44,45]. In fact, the surface rupture trace (segments A, B, C and E, Figure 3) of the Palu earthquake is remarkably linear and continuous. Compared to this event, the rupture of the 2001 Kunlun earthquake was remarkably linear and continuous on the whole supershear rupture, and did not produce complex surface patterns and major bends (segment D produced by this event). We infer that the supershear rupture did not necessarily produce linear and straight traces, and it will produce some geometric variations in the process of maintaining the supershear rupture.
In addition, our results show that segments A and B of the Palu earthquake fault ruptured within the Palu basin (2-5 km from the western bedrock mountains), and did not reproduce the long-term fault structures of the Palu fault. The location of the Palu fault is well-documented based on investigations using tectonic geomorphology and geodetic data. According to these investigations, the Palu fault developed in bedrock along the junction of basin (the Palu basin) and mountains, and ruptured southward from the Palu coastline. In theory, the active faults (e.g., with bedrock properties, small friction force, without big barriers, etc.) are more suitable for supershear rupture [48]. What fault segments A and B, ruptured within the Palu basin by supershear velocity, may indicate that the basement of the Palu basin at segments A and B have the conditions to produce and maintain the supershear rupture. From here, despite the fact that we are unable to obtain the coseismic rupture trace using satellite images within the Palu bay, we further infer that the rupture propagating southward from the Sulawesi neck is closer to segment A. In other words, the material medium suitable for maintaining the supershear rupture has been found at segment A, and hence it is a good explanation for why the supershear rupture did not rupture along the bedrock mountain.
In addition, segments A and B correspond to a cross-basin fault system within the Palu basin, for which the fault system was identified by tracing Palu River channels in six separate images from 2003 to 2015 [49]. Moreover, the rupture did not propagate southward along the cross-basin fault system, but instead moved westward for 0.5 km at the end of segment B, propagating along the bottom of the bedrock mountain in segments C and D (following the Palu fault). Since the rupture was a persistent supershear event [25,26,42], the movement of the rupture at the end of segment B is reasonable. The material properties of segment C is more conducive to sustain supershear velocity. Another possibility may also explain why the rupture moved westward at the end of segment B. We think a large obstacle here formed a strong barrier for supershear rupture. Meanwhile, it remains unclear why the rupture propagated along segments A and B within the Palu basin instead of moving west along the bedrock mountain. The geometric variation controlling the surface rupture of segments C and D may have sustained the supershear rupture by continually adjusting the strike and location [49,50]. Hence, we infer that the characteristics and geometric variation of the ruptured fault played a significant role in the surface deformation patterns produced, and were needed to sustain the supershear rupture.

Conclusions
The coseismic horizontal displacement fields of the 2018 Mw7.5 Palu earthquake were obtained from Sentienel-2 optical images using subpixel image correlation. The curl, divergence, and shear strain fields derived from the NS and EW displacement fields provide a precise and simple means to analyze the characteristics and geometric variation that controlled the surface rupture. The coseismic surface rupture had a length of 65 km and mean strike of 343.0 • , it was predominately a left-lateral strike-slip event (maximum strike-slip of~6.62 m), with a minor normal component on an east-dipping fault peaking at~3.75 m. Interesting deformation patterns were precisely identified from curl, divergence, and shear strain maps, including a small restraining bend (area of localized shortening), a distributed rupture zone, and a major releasing bend (area of net extension). Our results show that high-resolution horizontal displacement fields combined with curl, divergence, and shear strain allow for quantitative and automatic mapping of the characteristics and geometric variation of surface ruptures.