Crustal Deformation and Fault Strength of the Sulawesi Subduction Zone

This paper investigates the seismicity and rheology of the North‐Sulawesi subduction zone. Body‐wave modeling is used to estimate focal mechanisms and centroid depths of moderate magnitude (M5–M6.5) earthquakes on the North Sulawesi megathrust and surrounding region. The slip vectors of megathrust earthquakes radiate outward from Sulawesi, indicating motion that is incompatible with the relative motion of two rigid plates. Instead, the observed deformation implies lateral spreading of high topography, controlled by gravitational potential energy contrasts. This finding suggests that the observed deformation of Sulawesi results from stresses transmitted through the lithosphere, rather than basal tractions due to circulation in the mantle. Our modeling of the force balance on the megathrust shows that the subduction megathrust is weak, with an average shear stress of ∼13 MPa and an effective coefficient of friction of 0.03. Elsewhere in Sulawesi, slip vectors of other earthquakes suggest similar potential‐energy‐driven deformation is present, but at significantly slower rates. Our results show the importance of lateral rheology contrasts in determining deformation rate, and hence seismic hazard, in response to a given driving force.


Methods
We use teleseismically recorded P and SH body waveforms downloaded from the IRIS datacenter to evaluate the focal mechanisms and depths of earthquakes in North Sulawesi and the Celebes Sea. To reduce the effect of near source/receiver structure, and complexities arising from the arrival of other seismic phases, we restrict our analysis to stations at epicentral distances of 30-90°. The decision of which analysis technique to use on each earthquake is primarily determined by the magnitude of the earthquake. Earthquakes with a magnitude greater than ∼5.5 are typically large enough to produce sufficiently clear waveforms at teleseismic distances to use body-waveform modeling techniques to invert for the source depth, sourcetime-function and focal mechanism (Section 2.1). For earthquakes with magnitudes lower than ∼5.5, we initially attempt to model them using the same body-waveform modeling technique described in Section 2.1. Provided there are sufficient stations with high signal-to-noise ratios, we can sometimes model these earthquakes well. However, if the signal-to-noise ratio is not high enough, or stations do not sufficiently cover the focal sphere to result in a well-constrained inversion, we resort to a simpler technique to determine the earthquake depth only, using the near-source surface reflections pP and sP (Section 2.2). As there are no estimates of the crustal velocity beneath Sulawesi, we use 6.5 and 3.5 km s −1 for P-and S-wave speed, respectively, because they are a reasonable average for global basement velocities. These velocities are used to generate the Green's functions. Synthetic seismograms are then aligned with the observed arrivals, which means that the absolute timing of arrivals is removed from our analysis (so whole-Earth structure does not play a role), and only the relative timings of direct and depths phases are analyzed (which allows the earthquake depths to be estimated). This methodology gives us no constraints on the latitude and longitude of the events, only on the depth and mechanism, so we use the prior latitude and longitude location estimates from the ISC (ISC, 2020). Using a plausible range of velocities in the crust does not significantly alter the modeled focal mechanism (Craig et al., 2011) but does affect the estimated depth.

Body-Waveform Source Modeling
We first filter the seismograms using a long-period filter which emulates the response of a WWSSN 15-100 s instrument. By removing the shorter periods, we minimize the signals due to source complexity and nearsource velocity heterogeneity. This period range allows us to approximate the earthquake as a point-source, for events up to ∼M7. During filtering, the waveforms are detrended and resampled at one sample-per-second. The horizontal components are then rotated into a radial-tangential reference frame between source and receiver.
We inverted the recorded waveforms using the MT5 algorithm Zwick et al., 1994). This algorithm uses both the P and SH waveforms (recorded on the vertical and tangential components) to determine the best fitting double-couple focal mechanism, moment, source time function, and centroid depth. It uses an iterative least-squares solver to minimize the misfit between synthetic seismograms and the recorded data within user-determined windows ( Figure 2). A complete description of the technique can be found in Taymaz et al. (1991) and Nábělek (1984). Typical errors for this technique are ±4 km in centroid depth, ±10° in strike, ±5° in dip, and ±5° in rake (Molnar & Lyon-Caen, 1989;Taymaz et al., 1991). The error in slip vector is typically ±15° (McCaffrey, 1991 Sulawesi's northern arm, as discussed in this study, is highlighted by the black rectangle. Seismicity from the GCMT catalog, with a depth of less than 50 km, a magnitude greater than 6 and which are greater than 70% double-couple are plotted as focal mechanisms with gray compressive quadrants. Earthquakes reanalyzed using the techniques outlined in this paper are plotted with black compressive quadrants. Specific earthquakes, as discussed in the text are labeled with letters. (c) Normal-faulting earthquakes re-analyzed in this study (green compressive quadrants) and from the GCMT catalog (gray compressive quadrants), respectively. (d) Earthquake mechanisms in the outer-rise of the North-Sulawesi subduction zone. Surface velocities calculated using the GPS system and relative to stable Sundaland (Simons et al., 2007) are indicated by the black arrows. (Kouali et al., 2016;Mustafar et al., 2017;Socquet et al., 2006). The 95% confidence interval is indicated by the oval. (e) Events in the EHB (E. R. Engdahl et al., 2020) catalog colored by depth. Displayed focal mechanisms are from the CMT catalog with a depth greater than 50 km, a magnitude greater than 6 and which are greater than 70% double-couple are plotted as focal mechanisms with compressive quadrants colored by earthquake depth. (f) Events located along the Section A-B. The red and blue lines delineate the surface of the south-dipping Sulawesi slab and the west-dipping Halmahera (or Sangihe) slab, respectively.

Depth-Phase Modeling
Source modeling using the MT5 technique described above requires good azimuthal coverage of the focal sphere to produce high-quality solutions. However, when the azimuthal coverage is poor, or the number of stations which have an acceptable signal-to-noise ratio is too few, we opt to fix the focal mechanism and source time function to the global centroid moment tensor (CMT) solution and attempt to constrain the source depth only. We use broadband records and focus on the P-wave and its near-source surface reflections (i.e., pP and sP). Green's functions are generated at a suite of depths for each station that records the P-wave arrival, using the WKBJ method (Chapman & Orcutt, 1985). We convolve these Green's functions with the broad-band instrument response, an attenuation function (t* = 1.0), and the source time function from the GCMT solution. We then align the observed P-waveform and the synthetic on the first arrival, before normalizing each trace by dividing by its variance (Figure 2). A root-mean-square misfit in a window covering the depth phases is then used to compare the fit between the observed waveform and the synthetic at each depth ( Figure 2). The typical errors in depth using this technique are ±2 km for a given velocity model (Craig et al., 2011).
Our inversion approach initially focuses on identifying one or two stations with clear depth phases. We find the minimum misfit for these stations using a window which changes dynamically for each depth (window length = sP time -P time + 5 s). These results are used to define a static window to use over the entire data set of stations which record the earthquake (i.e., those with a P-wave arrival picked). We then compare the observed waveforms and synthetics for all of the stations simultaneously. We remove those stations which do not have identifiable depth phases. This may be the case if: (1) the depth phase has a take-off angle close to a nodal plane, (2) there are near source/receiver noise sources, (3) there is geometrical or temporal complexity in the source not described by the GCMT model, or (4) there is sufficient attenuation along the ray path.
For both of the techniques we use, a prominent feature of some of the data are waves which reverberate in the water layer, producing long-period signals after the depth phases (i.e., stations SYO and P124,Figure 2a). We attempted to model these reverberations to improve our model fits by including a water layer in our velocity model. However, despite the fits to some seismograms being improved, others were degraded due to the laterally varying water depth as a function of event-station azimuth, leading to differences in the period of the reverberations between seismograms of the same event. We instead opted to use variable data window lengths to remove the reverberations from the inversion window, or the removal of the affected stations for the final solutions where this was not possible. Authors who have used similar methods (e.g., Craig et al., 2011) found that the incorporation of a water layer did not change the final solution. This effect arises because the direct wave and the depth-phase from the sea-bed arrive earlier than the reverberations from the water layer, so those reverberations do not appreciably affect our estimate of centroid depth below the sea floor, or our estimated focal mechanism.

Earthquake Distribution and Mechanisms
The minimum-misfit solutions for the earthquakes we study are tabulated in the Supplementary Information along with the inversion results. They are plotted in Figures 1 and 3, along with the solutions from the CMT catalog for events that we were unable to study due to data availability. Low-angle thrust-faulting earthquakes (with nodal plane dips of 5-32° and depths of ∼10-30 km) occur along the megathrust offshore northern Sulawesi, along which the oceanic lithosphere of the Celebes Sea subducts southwards. They are located close to the coastline, and by analogy with other subduction zones probably represent earthquakes close to the down-dip limit of the locked parts of the megathrust (Hyndman, 2013;Ruff & Tichelaar, 1996;Simoes et al., 2004). These thrust events display a distinct spatial pattern in the azimuth of the slip vector.
GREENFIELD ET AL.
10.1029/2020TC006573 5 of 18 From east to west, slip vector azimuths rotate from −15° to 15° (Figure 4). We discuss the possible causes for this pattern in Section 4.
In the west of the subduction zone, where the arc terminates against the Palu-Koro strike-slip fault, we observe north-south-oriented, sinistral, strike-slip earthquakes ( Figure 1). To the east, the Sulawesi arc abuts against the easterly dipping Halmahera (or Sangihe) subduction zone. Here the observational earthquake record is poor but two northeast-southwest striking, high-angle (∼45°), thrust events are observed (labeled as E, Figure 3). The high-angle nature of these events suggests they are not likely to be due to failure of the megathrust. Normal-faulting earthquakes occur along the North-Sulawesi subduction zone with 14 observed in the GCMT catalog (Ekström et al., 2012) (Figure 1c). Most are in a single cluster midway along . Shallow megathrust structure constrained using teleseismically recorded earthquakes. (a) shows the seismicity which has been reanalyzed in this and previous studies (Gómez, et al., 2000;Vigny et al., 2002) plotted as focal mechanisms. Events with black compressive quadrants are interpreted to be on the megathrust, purple compressive quadrants on the Palu-Koro strike slip fault and red compressive quadrants in the down-going plate on the outer rise. The event indicated by the blue compressive quadrants is a normal event within the accretionary prism and is discussed in the text. Green compressive quadrants indicate normal events. The red lines contour the depth of the southward-dipping North Sulawesi Slab at 20 km intervals. Contours with depths greater than 50 km, which, due to over-smoothing and a paucity of earthquakes to map the slab, probably do not represent the true slab, are delineated as dashed lines (Hayes et al., 2018). Panel (b) and (c) show cross sections along A-B and C-D, respectively. The red line delineates the surface of the slab from the Slab2 model (Hayes et al., 2018). The black line is the surface of the slab delineated from this study. the arc from which we select a representative two for refinement. In our catalog of refined solutions three normal earthquakes are observed close to the northern coast of Sulawesi's northern arm (green mechanisms in Figure 3).
Earthquakes which have low-angle planes and are likely to have occurred on the megathrust have centroid depths between 15 and 30 km (with a single exception at 36 km). We identify no earthquakes in the far eastern part of the arc (longitude > 124°) occurring on the megathrust. The two high-angle thrust events in this region, which we associate with the accretionary prism, occur at depths of 10-17 km. A recent geometrical model of the Sulawesi subduction zone was compiled in the Slab2 model (red line, Figure 3) (Hayes et al., 2018). The earthquakes we identify as occurring along the megathrust have depths 7-20 km shallower than the model (Figures 3b and 3c) suggesting that, at least at shallow depths, the Slab2 model is a poor representation of the Sulawesi slab.
The Sulawesi outer rise is seismically quiet. The GCMT catalog records only nine earthquakes in this region ( Figure 1d). Bending stresses at subduction zone outer rises can cause normal and thrust earthquakes to be generated (Chapple & Forsyth, 1979;Craig et al., 2014). A single earthquake in the GCMT catalog is normal and its location is in the same place as a thrust event ( Figure 1d). The ISC catalog shows a large spread in GREENFIELD ET AL.
10.1029/2020TC006573 7 of 18  Mustafar et al., 2017;Socquet et al., 2006) relative to stable Sundaland (Simons et al., 2007) with black arrows. The 95% confidence interval is indicated by the oval. Slip vectors from earthquakes analyzed in this study and previous studies (Gómez et al., 2000;Vigny et al., 2002) are indicated by the red vectors. Gray arrows indicate the slip vector directions from earthquakes in the GCMT catalog. estimated locations from different agencies for these events, and it is unclear whether they occur in the subducting or over-riding plate. We therefore do not interpret these events.
The remaining events within the subducting plate have low magnitudes (M < 5.3) and we estimated their depths using a combination of the MT5 and the WKBJ methods. The three events for which we obtained well-constrained solutions (red compressive quadrants, Figure 3a) have depths of 20-33 km. Two of the events are thrust events, typical for the compressional stresses caused at depth by the bending of the oceanic plate. The remaining event is strike-slip, but we note that the mechanism is from the GCMT catalog and could have large errors due to the small size of the event.
In the following sections, we investigate the possible causes of the lateral variation in slip vector azimuth, and then combine our earthquake source models with free-air gravity anomalies and simple dynamic models to infer the forces driving the deformation and estimate the rheology of the region.

The Cause of Along-Strike Variations in Slip Vector Azimuth
One option to explain this variation in slip-vector-azimuth is to attribute it to consequences of the 3D geometry of the slab. Along-strike variability of megathrust slip vectors can occur if the slab curvature in the horizontal plane (referred to below as 'horizontal curvature') changes as it is subducted (Figure 5c) (i.e., for a planar slab there would be no change in slip vector orientation along the length of the trench, Figure 5a). If the slab does become more horizontally curved as it is subducted, we would expect to see a significant variation with longitude in the convergence direction between the Celebes Sea and Sulawesi. Based on the EHB earthquake catalog, there is no clear change in the horizontal curvature of the slab at shallow depths (<40 km, where the shallow deformation will control the slip vector azimuth on the megathrust) (Figures 1  and 3). Scatter in the earthquake hypocenters in this location means there could be some changes in horizontal curvature that we are not able to observe. However, any deformation associated with a progressive increase in horizontal curvature would be characterized by along-strike compression within the subducting slab ( Figure 5c). The earthquakes in the North Sulawesi slab indicate down-dip compression (Figure 1e), inconsistent with such a mechanism. These observations imply that along-strike shortening due to the de-GREENFIELD ET AL.
10.1029/2020TC006573 8 of 18 A second geometrical effect that could lead to lateral variations in megathrust slip vector azimuth would be the instantaneous bending of the slab into the curved subduction zone, even if its curvature then does not continue to develop at greater depths. In this case, the degree of along-strike variability of slip vector azimuth that can be developed can be obtained by simple trigonometry using the megathrust dip and the along-strike variability in trench azimuth. In this case, the possible along-strike variation in slip vector azimuth is <1 degree, which is smaller than that visible in the earthquake slip vector azimuths. This line of logic is equivalent to the geometrical analysis of Schettino and Tassi (2012), which predicts little deformation of the slab as it enters the subduction zone if the surface trace of the trench approximates a small circle.
Other possible causes of the variation in slip-vector azimuth are summarized in Figure 5. There could be along-strike compression in the down-going plate seawards of the trench (Figure 5b), or along-strike extension in the over-riding plate (Figure 5d), either of which would generate along-strike variation in the megathrust slip vector orientation. The lack of seismicity and bathymetric features in the Celebes Sea suggests that no such along-strike compression is present ( Figure 1). However, east-west extension on normal faults is observed to be active in northern Sulawesi, as shown on Figure 1c. This along-strike extension can accommodate the along-strike variation in megathrust slip vector azimuth. The cause of this extension is also suggested by the orientation of the megathrust slip vectors. These point radially to the local topography, down the local topographic slope. Such patterns have been observed elsewhere, and are thought to be due to the influence of gravitational potential energy contrasts leading to shortening along the topographic gradient (Copley & McKenzie, 2007;England & Houseman, 1986;Hatzfeld et al., 1997). Such a potential energy contrast exists in northern Sulawesi, between the continental crust of the island with its surface above sea level, and the old oceanic lithosphere of the Celebes Sea (much like the potential energy contrast between the Aegean Sea and the Eastern Mediterranean that is thought to play an important role in Aegean dynamics; England et al., 2016). We therefore infer that the along-strike variability of slip vector azimuth on the megathrust offshore northern Sulawesi is the result of gravitational potential energy contrasts driving shortening perpendicular to the local strike of the trench, accommodated by east-west extension in the over-riding continental crust.
In some subduction zones (e.g., The Sunda Arc, Figure 1a), along-strike variability in earthquake slip vectors is thought to be due to the partitioning of oblique convergence between the frontal thrust and strikeslip faulting within the overlying plate (McCaffrey, 1991). A continuous change in the slip vector can thus be caused by an increase in the proportion of slip accommodated on the arc parallel strike-slip system. However, no significant trench-parallel strike-slip systems are known in northern Sulawesi, implying such a mechanism does not play an important role in this region.
GREENFIELD ET AL.

Elastic Thickness
The elastic thickness of the oceanic lithosphere can be estimated by modeling the flexure, and resulting gravity anomaly, in response to loading by seamounts or in subduction zones (Bry & White, 2007;Contera-Reyes & Osses, 2010;McNutt, 1984). The simple structure of the oceanic plates (i.e., a single load-bearing layer) makes such modeling relatively simple, in contrast to the difficulties and controversies involved in estimating the elastic thickness on the continents (e.g., McKenzie & Fairhead, 1997;Watts & Burov, 2003).
We use bathymetry from SRTM15+ bathymetry (Tozer et al., 2019) to estimate the elastic thickness by modeling the flexure of the oceanic plate using the equations outlined in Turcotte and Schubert (2002). Swaths of bathymetry were generated perpendicular to the trench and aligned at the location of the trench axis. The profiles are then stacked to reduce noise and calculate a standard deviation ( Figure 6). We invert the bathymetry between the trench and a distance to the north of 270 km ( Figure 6). Our inversion approach follows that of Bry and White (2007) and includes parameters to account for regional bathymetric anomalies and any long wavelength topographic features. We initially loop over a suite of trial elastic thicknesses and allow the other fitting parameters (bending moment, load magnitudes and location, and a linear slope and offset in the data) to vary to generate a misfit curve ( Figure 6). The parameters from the minimum misfit solution are then used as starting parameters when we invert for all free parameters at the same time. Our best-fitting elastic thickness is   10 6 24 km, where the errors refer to the range of values with a misfit less than 150% of the minimum misfit value ( Figure 6). This error estimate partly arises from nonflexural signals in the bathymetric data, but mostly due to the tradeoffs between the elastic thickness, bending moment, and long-wavelength tilt ( Figure 6).
The subducting oceanic plate has an age of 35-45 Ma and our estimates of the elastic thickness lie well within the typical distribution of elastic thicknesses estimated for mature oceanic lithosphere (age > 25 Ma) (see compilations in A. B. Watts, 2001, p. 458 andCopley, 2014). In the continents, there is controversy surrounding the estimates of the elastic thickness, especially in cratonic regions. Using spectral methods, the estimated elastic thickness in cratonic regions varies between 20 and 30 km and >80 km (e.g., Audet & Bürgmann, 2011;Kirby, 2014;Maggi et al., 2000;McKenzie et al., 2015;Watts & Burov, 2003). The highest values are calculated using the coherence between Bouguer gravity anomalies and topography (Forsyth, 1985). In many cratonic areas, there is little coherence between the free-air gravity anomaly and topography. In these regions, McKenzie (2016) argues that the Bouguer gravity anomaly is controlled by the ratio of the power spectra of the free air gravity anomalies and the uncompensated topography, and contains no information about the value of the elastic thickness. Our elastic thickness for the Celebes Sea matches the lower range of cratonic elastic thicknesses calculated from the free-air admittance in regions with sufficient coherence between topography and gravity (McKenzie et al., 2015), and is consistent with space-domain estimates based on the flexure in foreland basins (Jackson, 2008). These comparisons suggest that in relation to other oceanic and continental estimates globally, the Celebes Sea contains relatively strong lithosphere.

Using Seismicity to Constrain the Plate Strength
The elastic thickness estimated above provides a simple measure of the strength of the plate and suggests the subducting Celebes Sea plate is relatively strong. However, the oceanic lithosphere is not purely elastic where it bends into subduction zones, as shown by permanent strain in earthquakes (e.g., Chapple & Forsyth, 1979;Craig and Copley, 2014). In this section, we therefore use the seismicity itself to place constraints on the intraplate fault strength in the subducting Celebes Sea plate.
We follow the approach of Craig et al. (2014), who used the mathematical framework outlined in McAdoo et al. (1978) to estimate the strength of faults in a suite of subduction zones around the globe. In this framework, the curvature of the plate in a vertical plane ( 2 2 / d w dx , where w is plate deflection and x is the horizontal distance) can be related to a depth distribution of the yield stress,    z , through the lithosphere, where  is the dip of the failing faults,  is the effective coefficient of friction along these faults,  is the density of the lithosphere (3,200 kg m −3 ), g is the acceleration due to gravity (9.81 m/s 2 ), and the  2 sin term is positive for thrust faults and negative for normal faults.
We calculate the vertical curvature of the subducting plate directly from the SRTM15+ bathymetry (Tozer et al., 2019) by differentiating the bathymetry twice along an azimuth perpendicular to the local trench strike (Figure 7). We then extract profiles along the length of the trench and stack to reduce noise and GREENFIELD ET AL.  generate a single profile of the curvature of the plate in the vertical plane (Figure 7). Despite the appearance of a smooth Celebes Sea plate, filtering of the topography prior to differentiation is required to reduce the noise present in the data set. We tested Gaussian filter widths between 10 and 150 km, and settled on 50 km as the optimum value. The tested filter widths are somewhat higher than those used by Craig et al. (2014) (maximum of 40 km), and the shallowing bathymetry landward of the trench can alter the recovered curvature. However, it is clear from Figure 7 that the most negative curvature of the Sulawesi slab is between 70 and 150 km from the trench, so our chosen filter width has little effect on this result. Using our suite of curvature profiles, the most negative curvature ranges between −0.20 × 10 −6 and −0.12 × 10 −6 m −1 (Figure 7). These values are similar to the most negative curvature calculated from the flexural model in Section 4.1, suggesting our range of values is robust.
The approach of Craig et al. (2014) uses the difference in depth between the reverse and normal earthquakes found at the outer-rise of subduction zones to estimate the width of the elastic core. The Sulawesi outer-rise does not have any normal-faulting earthquakes recorded teleseismically, despite the observed strain at the surface due to curvature in the vertical plane (∼10 −3 ) being an order of magnitude greater than the typical elastic limit of rocks (10 −4 ). The lack of normal-faulting earthquakes is most likely due to the repeat-time of the outer-rise seismicity compared to our observation time (∼50 years). As an upper bound, if we assume the elastic core extends from the surface to the minimum depth of compressional earthquakes, the largest possible thickness of the elastic core is 22 km (Table 1). Using this value and the curvature of the subducting slab we can calculate the yield stresses using Equations 1 and 2. Assuming the yield stress is constant with depth, we calculate that  0 is 110-190 MPa. Using our second model, in which the yield stress varies linearly with depth and a fault dip of 45°, the effective coefficient of friction is estimated to be between 0.14 and 0.21, with the yield stress at the depth of the compressional earthquakes being 230-380 MPa (Table 1).
As discussed above, the lack of normal earthquakes in our catalog of outer-rise earthquakes is probably a result of the short observation time relative to seismic cycles rather than the oceanic plate being elastic all the way to the surface. It is therefore useful to consider what the strength of the plate would be if the elastic core was thinner than we can constrain using the earthquake catalog as it currently exists. In the global survey of outer-rise seismicity in Craig et al. (2014), typical elastic core widths are less than ∼10 km. The results of our calculations using elastic core widths of 5 and 10 km are presented in Table 1.
As expected, the reduced elastic core thickness used in these calculations lowers the estimated strength of the plate. The values we obtain are similar to the estimates Craig et al. (2014) made for the Middle America, Kuril, and Tonga trenches for the upper range of our considered elastic core widths, but lower for the smallest values considered here. An implication of this result is that there are either significant lateral variations in oceanic plate strength between this region and those studied by Craig et al. (2014) or that the elastic core here is likely to be wider, and the plate strength more similar (which arises from the lower curvature we observe here compared to the locations studied by Craig et al. (2014)).
Our calculated strengths are also lower than those calculated for the oceanic lithosphere using seismological techniques (e.g., >500 MPa; Choy & McGarr, 2002;Wei et al., 2013). Those studies use estimated earthquake depths and lithospheric strength profiles based on Byerlee's law for the brittle layer with a high coefficient of friction (0.6) and a dry crust and mantle (Kohlstedt et al., 1995). The difference between our estimates may arise from faulting in the outer rise acting to channel water into the lithospheric mantle beneath the oceanic crust and thus reducing the strength relative to the oceanic lithosphere close to mid-ocean ridges (Choy & McGarr, 2002). It is notable that both our estimated elastic thickness, and the strength of the faults, is similar to the those estimated for cratonic regions in the continents (i.e., McKenzie et al., 2015). This similarity suggests that although oceanic and continental lithosphere have different buoyancy relative to the underlying mantle, and so behave differently in convergent settings, differences in mechanical properties probably do not play a major role in these differences in behavior.
GREENFIELD ET AL. Note. z1 is the maximum depth of normal earthquakes, z2 the shallowest reverse earthquake and Δz is the difference in depth between the normal and reverse earthquakes and thus the width of the elastic core. σ 0 is the yield stress calculated assuming a constant yield stress with depth and σ(z1) and σ(z2) are the yield stresses in a depth dependent yield stress model evaluated at z1 and z2, respectively. μ is the effective coefficient of friction estimated for a depth-dependent yield stress model. a Depth calculated assuming a fixed elastic core width. b Range of values calculated from the range in most negative plate curvature in the outer rise and trench slope.

Megathrust Strength
As described in Sections 3 and 4, megathrust earthquakes on the North-Sulawesi subduction zone have slip vectors which are not parallel, as would be expected for the relative motion of rigid plates. Instead, they radiate outward, away from the high topography (Figures 3 and 4). Such a deformation pattern would be expected to cause arc-parallel extension behind the frontal thrusts (e.g., Armijo et al., 1986), which would manifest as normal faults striking perpendicular to the arc (Figure 5d). Such normal-faulting earthquakes are observed in northern Sulawesi (Figure 3). It is likely that we do not observe normal earthquakes further north, closer to the trench, because it is not mechanically feasible for along-strike extension to occur where the underlying megathrust is locked as the North-Sulawesi subduction zone is (Socquet et al., 2006). We suggest that it would be likely that such earthquake will be observed after megathrust events. This change in strain would be similar to that observed after the Tohoku earthquake (Wang et al., 2019), but with the difference that rather than the resulting extension direction being arc-perpendicular, we would expect it to be arc-parallel (or both, if the effects described for the Tohoku event were also to occur in this region).
As discussed above, slip vectors radiating outward from regions of high topography are indicative of deformation being driven by gradients in gravitational potential energy, that is, motion from high topography to low topography. This configuration allows us to model the forces driving the deformation using two end-member approaches (Figure 8). The first approach is to assume that the long-term strength of the crust lies in the brittle portion of the megathrust (Figure 8a). The second approach considers the brittle portion of the megathrust to be extremely weak, so that the long-term strength lies within the viscously deforming portion of the megathrust, downdip of the brittle-ductile transition (Figure 8b). In this approach, the movement along the megathrust is limited by the viscous strength of the shear zone. These two approaches represent end-member states, a combination of which is likely to exist on the Sulawesi subduction zone. However, by considering the end members we can put bounds on the strength of the megathrust regardless of what mechanism (or combination) is valid. As we shall see, the calculated shear stress is similar in either case.

If Strength Lies in the Brittle Part of the Crust
If we assume the strength of the megathrust lies within the brittle portion of the interface, we can model the magnitude of the shear stress by following the method outlined by Lamb (2006), which balances the forces on the wedge of material overlying the subduction megathrust. Three forces are considered as acting on the subduction wedge: the weight of the wedge, the frictional resistance to sliding along the megathrust, and a force acting horizontally on the back of the wedge (Figure 8a). We calculate the weight of the wedge using a constant density of 2,800 kg m -3 and consider a simple brittle rheology with a constant effective coefficient of friction to model the shear stress along the megathrust. The force acting on the rear of the wedge is calculated assuming no deviatoric stresses in the arc. In this case, the horizontal force can be equated with the summed pressure of the rock column beneath the arc (Lamb, 2006), which can be calculated by assuming a GREENFIELD ET AL.

(a) (b)
constant crustal density of 2,800 kg m −3 . The normal-faulting seen across Sulawesi's northern arm indicates that there is a re-orientation of the principal stress axes in that region compared to on the megathrust, and we follow Lamb (2006) in suggesting that the region, therefore approximates the underlying assumption of negligible deviatoric stress in the arc. Any deviatoric stresses in the arc that are present will cause a change in the stress along the megathrust equal to   Δ 2 sin (where  Δ is the strength of faults in the arc and  is the dip of the slab). Assuming all the stress is supported in a seismogenic region 15 km thick in the region of the arc, a slab dip of 9º, and a maximum representative effective coefficient of friction equal to 0.2, we can expect a maximum error of <5 MPa in our final estimate of stress along the megathrust due to unmodelled deviatoric stresses in the arc. The other major source of error in these calculations is the dip of the subducting slab. We use the mean dip of the megathrust earthquakes (9°) as the dip of the slab. As the earthquakes are clustered beneath the coastline this dip is likely a maximum value because subduction slabs typically bend convex-landwards (Hayes et al., 2018). An error in slab dip of 2° leads to an error of 3 MPa in our estimate of the shear stress on the megathrust. The possibility of unmodeled deviatoric stresses in the arc are therefore the largest contribution to the possible error in our megathrust shear stress estimates.
Using the parameters and method described above, the calculated average shear stress on the Sulawesi megathrust to a depth of 30 km is 13 ± 5 MPa. To calculate , we assume that shear stress increases linearly with depth and therefore the peak shear stress (at the depth of the deepest interface earthquakes) is twice the mean shear stress. Following this method, we estimate  to be 0.03, similar to that of other megathrusts (Lamb, 2006), and to results from critical wedge theory applied to accretionary prisms (e.g., Lallemand et al., 1994). If we assume the cause of the low coefficient of friction is entirely due to high pore fluid pressure, we can expect values up to 95% of lithostatic pressure. However, it is likely that a combination of both lubricating rock types and pore fluid pressures act to lower the coefficient of friction (Copley, 2018).
In situ measurements of the coefficient of friction on megathrust fault rocks are rare. The static friction coefficient calculated for fault rocks sampled at a depth of ∼800 m below the sea floor at the site of the 2011 Tohoku-oki earthquake (Brodsky et al., 2020;Ikari et al., 2015) was 0.26. However, the relatively shallow depth of the sample location does not necessarily make this value representative of the entire megathrust. However, if we assume this value for the intrinsic coefficient of friction, then a pore fluid pressure equal to 88% of lithostatic pressure would be required to generate our estimated effective coefficient of friction in Sulawesi.

If Strength Lies in the Ductile Crust, Down-dip of the Brittle Megathrust
If we instead assume the brittle part of the megathrust is weak, the strength of the interface must lie in the viscously shearing part of the plate boundary, down-dip of the brittle megathrust, and be controlled by the viscosity in that region (Figure 8b). We estimate the viscosity using the approach outlined in Copley and McKenzie (2007) to model the observed surface velocities resulting from viscous flow over a rigid base ( Figure 8b) (Huppert, 1982;McKenzie et al., 2000). We assume the lower boundary is gently sloping, and represents the upper surface of the subducting plate. Using the observed surface velocities of 45 mm yr −1 (Socquet et al., 2006) the calculated viscosity is 5 × 10 19 -3 × 10 20 Pa s. The range in values is primarily caused by uncertainty in the dip of the subducting plate. Our quoted values are for dips in a range of 5-15°. This value is similar to the viscosities calculated for Tibet, and the Indo-Burman and Makran accretionary wedges (Ball et al., 2019;Copley & McKenzie, 2007) and is consistent with laboratory flow laws (Bürgmann & Dresen, 2008).
In this model, the viscous stresses on the shear zone at the base of the crust vary between 7 and 21 MPa. This range encompasses the strength estimates calculated above, assuming that the strength lies within the brittle part of the megathrust. Thus, our two approaches (viscous or brittle) have yielded approximately the same megathrust strength. Therefore, although we are not able to establish in which depth interval the dominant stresses are supported, we are able to estimate the stresses transmitted across the interface between the subducting lithosphere of the Celebes Sea and the over-riding crust of Sulawesi to be in the range ∼7-21 MPa.

Discussion
Radial motions in response to gravitational potential energy contrasts have been suggested for many mountain ranges to explain the observed pattern of deformation (e.g., England & Houseman, 1986;Hatzfeld et al., 1997;Copley & McKenzie, 2007). For example, the under-thrusting Indian shield acts as a rigid boundary over which the thickened (and therefore hot) Tibetan crust spreads out. The pattern of deformation in Tibet, as revealed using earthquake focal mechanisms and GPS derived velocities, is similar to that of the cartoon in Figure 5d. Shallow-dipping, frontal thrust earthquakes show radiating slip vectors away from the high topography. Behind them, a band of normal earthquakes accommodate arc-parallel extension (Armijo et al, 1986).
The pattern of seismicity we observe around Sulawesi's northern arm is similar to the pattern observed in the Himalaya, suggesting that the same processes cause the deformation in both regions. Our estimate of the yield strength of the subducting oceanic lithosphere (30-80 MPa) is higher than our estimated strength of the megathrust (<20 MPa and <0.03, respectively), confirming that the region acts as a region of weak crust spreading out over the strong subducting plate. In this situation, the deformation is governed by forces transmitted through the lithosphere, rather than stresses due to circulation in the underlying mantle. The force system in the underlying slab is implicit in our analysis because it plays a role in determining the geometry of the down-going plate (i.e., slab dip or slab curvature), but does not directly contribute to driving the deformation.

Deformation Elsewhere in Sulawesi
Having established that the northern arm of Sulawesi is deforming in response to gravitational potential energy contrasts, we now discuss whether similar processes may act on the other margins of Sulawesi. South of the Palu-Koro fault system, GPS observations (Socquet et al., 2006) reveal that the relative motion between Sundaland and Sulawesi is small (<2 mm yr −1 ). Nevertheless, in some places, significantly sized earthquakes have occurred. In particular, two large (M7.1 and M6.7) shallow-dipping thrust events which occurred in the Makassar Strait in 1969 and 1984 (earthquakes labeled A and B in Figure 1) (Fitch, 1972;ISC, 2020) and two normal-faulting events, which occurred in the Tolo Trough (M5.6-1998 -M5.9-2009) labeled earthquakes C in Figure 1 are probably indicative of deformation driven by gravitational potential energy contrasts. In the case of the thrust earthquakes, shortening is occurring in the orientation perpendicular to the topographic gradient. In C's case the normal-faulting earthquakes may be indicative of arc-parallel extension behind the main frontal thrust of the Tolo Trough, in a similar manner to that described above for North Sulawesi.
These observations highlight the importance of lateral rheology contrasts as a first-order control on both the style and magnitude of deformation. Topographic contrasts between the island and the bounding seafloor exist around the entire perimeter of Sulawesi; therefore, the magnitude of the resultant force must be similar. However, significant deformation is only observed in the north where there is a subduction zone and therefore a zone of weakness. In contrast, such a plane is not available in the continental collision between Borneo and Sulawesi and thus they are relatively strong, and the resulting deformation is slow compared to that on the megathrust. The slow over-thrusting on the margins of Sulawesi in regions where subduction is not currently occurring indicates the role of pre-existing potential energy contrasts in providing a driving force for the over-thrusting, which could eventually result in subduction initiation.
Our analysis has possible implications for estimates of seismic hazard in this and similar regions. In Sulawesi, the hot and weak continental fragment is spreading out under its own weight. The style and rate of deformation is determined by the rheology on the margins of the island. In the north, where our calculated values for the megathrust indicate that it is weak, deformation rates are high, and we observe relatively frequent (100s of year repeat interval) potentially tsunamigenic earthquakes. In contrast, where there is no weak fault (i.e., in the Makassar Strait) deformation rates could be slow, and earthquake repeat intervals could be much larger. However, there is potential for infrequent, large magnitude, and hazardous earthquakes in these regions because of the presence of the driving potential energy contrasts across the margins of the island.

Conclusions
We have analyzed earthquakes around the island of Sulawesi in Indonesia, to study the active deformation in this poorly studied region. The slip vectors of the megathrust earthquakes along Sulawesi's northern arm show a radiating pattern incompatible with the relative motion of two rigid plates. We instead suggest that the accretionary wedge above the subducting Celebes Sea deforms in response to gravitational potential energy contrasts, spreading from a region of high gravitational potential energy to a lower one.
We use the curvature of the subducting Celebes Sea plate in the vertical plane, and the depth of seismicity in the outer-rise, to show that the subducting oceanic lithosphere is strong. It has a high elastic thickness (24 km), more similar to continental cratonic lithosphere than actively deforming regions. The shear stress at failure on intraplate faults in the outer-rise is calculated to be between 20 and 110 MPa, and the effective coefficient of friction is between 0.02 and 0.08. In contrast, the interface between the subducting and over-riding plates is calculated to be very weak. Regardless of whether the stresses are supported in the brittle or ductile part of the fault, we estimate the shear stresses to be in the range 7-21 MPa, equivalent to an effective coefficient of friction of 0.03. Our results imply that the deformation is driven by stresses transmitted though the lithosphere, rather than tractions on the base of the lithosphere caused by circulation in the underlying mantle.

Data Availability Statement
Seismic data used in this study are available to download from the IRIS datacenter (https://www.iris.edu/ wilber3). Focal mechanisms from this study are tabulated in the Supplementary Information and are available from Zenodo with the https://doi.org/10.5281/zenodo.4421239. We have made extensive use of the ISC earthquake bulletin and EHB catalogs (E. R. Engdahl et al., 1998Engdahl et al., , 2020ISC, 2020;Weston et al., 2018) and the GCMT catalog of source mechanisms (Dziewonski et al., 1981;Ekström et al., 2012).