Coseismic Slip Distribution of the 2019 Mw 7.5 New Ireland Earthquake from the Integration of Multiple Remote Sensing Techniques

Abstract

The 2019 Mw 7.5 New Ireland earthquake occurred at an equatorial area where the dense vegetation prevents remote sensing techniques such as C- or X-band interferometric synthetic aperture radar (SAR) from acquiring coherent phase measurements. Therefore, in this paper, multiple remote sensing techniques including the L-band interferometric SAR, the range and azimuth offset tracking of SAR intensities, and the offset tracking of optical images were employed to map its co-seismic deformation field and to determine the slip distribution. The surface rupture was clearly and consistently captured by all offset observations, with the ground fault trace striking at an angel of 315° and extending over 10 km. An iterative weighting strategy based on the residual root mean square of inversions using individual datasets was developed to determine the relative weight of each dataset, allowing for the joint inversion of the fault geometry, the refinement of the dip angle, and the determination of the best fitting slip distribution. The resultant model indicates a nearly left-lateral strike-slip motion on the Weitin fault that ruptured to the surface with a maximum slip of 6.10 m, occurring at a depth of ~10 km, and a geodetic moment release of 1.03 × 1020 Nm, corresponding to a magnitude of Mw 7.31. The distribution of aftershocks shows about 70% of aftershocks were located in the area with increased Coulomb failure stress and few aftershocks in the subduction zone to the south of the Weitin fault were triggered by this event.

1. Introduction

On 14 May 2019 (UTC 12:58:25), an Mw 7.5 earthquake occurred in New Ireland, eastern Papua New Guinea (Figure 1). It ruptured the Weitin fault, a strike-slip fault across the south of New Ireland, along the boundary between the South Bismarck and North Bismarck microplates. The epicentre initially estimated by the U.S. Geological Survey (USGS) was 4.081° S, 152.569° E, with a focal depth of ~10.0 km, indicating a shallow strike-slip fault. The earthquake was followed by 43 aftershocks with a magnitude > Mb 4.0 within five days of the main shock. Two moderate aftershocks on 14 May, both with a magnitude of Mb 5.0, respectively occurred ~40 min and ~100 min after the main event and were located ~18 km southeast and ~70 km northwest of the main shock. However, the largest aftershock occurred on 17 May 2019 with a magnitude of Mb 5.9, located ~74 km southeast of the main event. The earthquake catalogue from USGS show that the focal depth of the three aftershocks was 10.0, 12.1, and 21.0 km, respectively.

Papua New Guinea, located in a complex tectonic setting between the Pacific plate and Australian plates (Figure 1a), is one of the most seismically active regions in the world. Several microplates in its centre and west including the Solomon Sea, South Bismarck, and North Bismarck microplates compose part of the edge of the Pacific and Australian plates [1,2]. As a transform boundary [3], more than 35 earthquakes with a magnitude greater than Mw 7 have occurred since 1970. These include a well-documented Mw 8.0 earthquake that occurred on 16 November 2000, 45 km to the northwest of the 14 May 2019 event, which was followed by two Mw 7.8 aftershocks on the subduction zone between the Solomon Sea and Pacific plates [4]. Tregoning et al. [5] used teleseismic data to relocate the aftershocks and concluded that the accumulated strain on the Weitin fault was fully released during the main shock because no aftershocks were located in the upper 15 km of the Weitin fault. Determining the coseismic slip distribution of the 16 November 2000 event constrained by the teleseismic wave data, Geist and Parsons [6] modelled the changes in the Coulomb failure stress, which they inferred had contributed to the first Mw 7.8 thrust aftershock. However, Park and Mori [7] carried out inversions using teleseismic P waveforms and showed that the static stress triggering mechanism hardly explains all of the triggered events in the earthquake sequence. To conclude, uncertainties exist due to limited observations for the investigation of the interaction of the faults in New Ireland and subduction trenches (e.g., New Britain trench), especially during large strike-slip earthquakes that occurred on the Weitin fault. Therefore, the 14 May 2019 event, covered for the first time by abundant geodetic observations with a high spatial resolution, provides a great opportunity to study the detailed fault behaviour and potential seismic hazard around New Ireland.

Interferometric synthetic aperture radar (InSAR) has been widely used for tectonic and volcanic studies since the early 1990s [8,9,10]. However, their phase measurements may significantly lose coherence in this heavy vegetation (as shown in Figure 1b,c) area, especially for sensors with short wavelengths such as X- and C-bands, limiting their usages in tropical regions. The decorrelation and unwrapping problems intensify for areas with large displacement gradient. It is therefore reasonable to employ multiple remote sensing observations that are feasible over vegetation-covered areas such as the L-band interferometric phase measurement (ALOS-2), the C-band SAR pixel offset (Sentinel-1), and the optical image pixel offset (Sentinel-2). SAR pixel offset tracking based on SAR amplitude images, also referred to as incoherent speckle tracking, can provide unambiguous surface displacement in both the line of sight (LOS) and azimuth directions by cross-correlating intensities [11,12,13]. The accuracy can reach 1/10 of one single-look pixel, which indicates its dependency on data pixel size [14,15]. Optical based offset tracking can measure ground deformation with a subpixel accuracy [16] in cloudless areas regardless of displacement gradient. Observations from these methods have been individually used for earthquake modelling [17,18,19], but few studies have considered combining them all for joint modelling due to their different satellite geometries, observation accuracies, and a proper weighting strategy.

In this study, the observations above-mentioned are combined using an iterative weighting strategy for a joint earthquake modelling. First, satellite data from Sentinel-1, Sentinel-2, and ALOS-2 PALSAR-2 in both ScanSAR and strip map (SM) mode were collected and processed to map the coseismic deformation field and determine the surface trace of the ruptured fault. With the assumption that the Earth crust model is elastic half-space homogeneous [20], finite-fault slips for the earthquake were then jointly inverted from these satellite images. Coulomb failure stress changes were finally calculated and discussed for the uncertainty of the interaction between the Weitin fault and subduction trenches.

2. Data and Processing Strategy

Sentinel is a series of Earth observation missions developed by the Copernicus initiative and operated by the European Space Agency (ESA). The Sentinel missions are operated day and night, providing complete, free, and open-access products through the Copernicus Open Access Hub. As the first of five Sentinel missions, Sentinel-1 performs C-band SAR imaging with a constellation of two satellites, Sentinel-1A and Sentinel-1B, running on the same orbital plane. Since the launch of the first satellite in April 2014, Sentinel-1 can offer repeated wide swath (~250 km) coverage and acquire imagery globally every six or 12 days regardless of the weather with the mode of Terrain Observation with Progressive Scans (TOPS) [21]. Another mission, the Sentinel-2, aims at providing multi-spectral and high-resolution optical imagery, carrying an optical instrument payload that samples 13 spectral bands: four bands at 10 m, six bands at 20 m, and three bands at 60 m spatial resolution [22]. With the launch of the first Sentinel-2 satellite in June 2015, the mission can offer systematic global coverage of land surfaces with an orbital swath width of 290 km and a high revisit frequency (e.g., five days at the Equator) [23].

As shown in

Table 1. Satellite data used in this study.

Satellite	Flight Direction	Reference Date	Repeat Date	Perpendicular Baseline (m)
Sentinel-1	Descending, right looking	13 May 2019	25 May 2019	1.1
Senitnel-2	Orthogonal looking	11 May 2019	26 May 2019	-
ALOS-2	Descending, right looking	12 May 2019	23 June 2019	114.6
ALOS-2	Ascending, right looking	09 March 2019	01 June 2019	598.9

, the Sentinel-1 image pair we used in this study was captured on 13 May 2019 and 25 May 2019 from descending track 16. The interferogram was mosaicked from two consecutive frames to cover a larger spatial extent. Considering the ruptured fault as a left-lateral northwest-striking fault with long-standing motion [24], we discarded the usage of the ascending track of Sentinel-1. GAMMA [25] was used to process SAR images in the single look complex format (level 1). To suppress speckle noise, a 20 × 4 multi-looking factor was applied in the range and azimuth. Compared with traditional SAR interferometry, processing Sentinel TOPS data requires a much more stringent image registration due to the rotation of antenna during the observation of each burst. The accuracy of coregistration in azimuth should be better than 0.001 of a pixel to avoid phase jumps at the interface between adjacent bursts [26]. To achieve such a high accuracy, an iterative amplitude matching procedure was conducted on the SLC after estimating terrain-induced pixel offsets with precise orbits from ESA and a 30 m digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) [27]. Once the azimuth offset correction is smaller than 0.02 pixels [28], a spectral diversity method considering interferometric phase in the burst overlap areas was used to further reduce coregistration errors. Upon the completion of high-accuracy coregistration, the conventional two-pass SAR interferometry method was followed. The first step was to simulate the topographic phase with the SRTM DEM, from which the differential interferogram was generated. The second step was to filter the interferogram using adaptive spectral filtering [29] and unwrap it using the minimum cost flow (MCF) method [30].

After the above processing, however, we found that the low coherence in this region made the unwrapping extremely difficult. The coherence loss was mainly caused by the heavy vegetation from tropical rainforests. Considering that the pixel spacing of Sentinel-1 in the range direction is about 2.3 m, it is possible to capture range offsets in the near field with SAR pixel offset tracking. The method will first search for the maximum cross-correlation between master and slave image windows after the high-accuracy coregistration, and then it will calculate the offsets between corresponding pixels. Following the offset tracking module in GAMMA software [31], we used a SLC offset search window of 300 × 60 pixels and a cross-correlation function window of 32 × 32 pixels. The cross-correlation coherence threshold for the acceptance of offsets was set as 0.1. After the offset tracking, we further used a median filter (9 × 9) to reduce the noise. Figure 2a shows the Sentinel-1 range offsets, where the red star indicates the epicentral location. Although the offsets were not clean enough and terrain-related residuals still seemed to exist, near-field coseismic deformation could clearly be seen. The maximum surface displacement in the range direction reached 2 m. Another finding was that the ground trace of the ruptured fault, extending over 10 km, could be easily recognized from the offset map.

The pre- and post-earthquake Sentinel-2 images (11 May 2019 and 26 May 2019) we used in the study were relatively cloud-free in the near field. We chose two images in band 8 with a resolution of 10 m and processed them using the COSI-Corr software package [16]. We used a sliding multi-scale window (initial 64 × 64 and final 32 × 32 pixels), a step size of four pixels, and four robustness iterations to optimize the masking of noise frequencies [17]. To further reduce noise, we discarded outliers greater than 5 m, detrended the displacement map with a linear ramp estimated from a spatial subset that excluded the area near the rupture, and denoised the results with a non-local means filter [32]. Figure 2b,c show the final displacement map in the north and east direction, respectively, where cloud-covered areas on land are masked. The surface trace of the fault rupture is visible on both maps and is consistent with that displayed in the Sentinel-1 range offset map (Figure 2a). To compare Sentinel-1 and Sentinel-2 deformation measurements in the same direction, we converted the Sentinel-1 range and azimuth offsets to horizontal displacements with Equation (1). The vertical component was not included in the equation since the fault rupture was mainly controlled by strike-slip motion. As shown in Figure 3a,b, the horizontal deformation map from the Sentinel-1 offsets was noisy due to the added azimuth offsets, but the deformation pattern approximately matched those in the Sentinel-2 deformation maps (Figure 2b,c). We further compared the deformation on the two profiles (denoted in Figure 3a), as shown in Figure 3c–f and found that the dispersion of the Sentinel-2 offsets was higher, indicating a lower precision compared to Sentinel-1.

Inherited from the Advanced Land Observing Satellite (ALOS) in 2014, ALOS-2 carries a phased array type L-band synthetic aperture radar-2 (PALSAR-2) sensor to acquire SAR images. PALSAR-2 has three observation modes: spotlight, strip map (SM), and ScanSAR, where the resolution successively decreases from 3 m to 100 m while the observation width increases from 25 km to 350 km. The revisit time also varies from several months to two weeks, depending on the observation mode and sensing areas. Compared with the C-band wave, the L-band wave can better penetrate vegetation to obtain ground information due to a longer radar wavelength [33], which facilitates the maintenance of coherence in rainforest areas.

According to the data availability for our study area, one descending ALOS-2 image pair from track 6 (12 May 2019 to 23 June 2019) in ScanSAR mode and one ascending image pair from track 108 (09 March 2019 to 01 June 2019) in SM mode were collected and processed. The processing for ALOS-2 ScanSAR images follows the conventional InSAR two-pass method after the mosaic of different swaths. We masked out those decorrelation areas with a coherence threshold of 0.4 to conduct phase unwrapping. The unwrapped displacement map, as shown in Figure 2d, is difficult to completely interpret the coseismic deformation in the near field because of coherence loss. As for those isolated areas such as the island to the southwest of the epicentre, phase unwrapping may not work well due to the lack of effective phase linking. Coherent pixels in the map can still be used in the slip modelling considering the high precision of SAR interferometry.

The ALOS-2 SM images covering the seismic region is in an ascending track, so the azimuth observations that are almost parallel with the ruptured fault should be visible. We measured the azimuth offsets with the offset tracking module in GAMMA, where the size of the SLC offset search window was set as 128 × 128 pixels. The fault trace can also be seen from the ALOS-2 azimuth offsets (Figure 2e), although within a smaller spatial extent. However, the offsets seem to be noisier than the Sentinel-1 range offsets, especially in the northeast edge, which may be caused by the larger azimuth pixel spacing (3.8 m) in ALOS-2 SM images.

$$\left(\begin{array}{l}East\\ North\end{array}\right)={\left(\begin{array}{cc}-\cos \alpha \sin \theta & \sin \alpha \sin \theta \\ \sin \alpha & \cos \alpha \end{array}\right)}^{-1}\left(\begin{array}{l}Range\\ Azimuth\end{array}\right)$$

(1)

3. Geophysical Modelling

The moment tensor solution from USGS [4] based on teleseismic waveforms suggests two nodal planes. Investigating the deformation patterns from the above images indicates a left-lateral slip component. Therefore, we used a fault plane dipping to the northeast for the inverse of the detailed fault geometry and its slip distribution. According to the Sentinel-1 (Figure 2a), Sentinel-2 (Figure 2b,c), and ALOS-2 (Figure 2e) offset maps, we estimated a surface trace with a strike angle of 315°, which best separates the directions of deformation.

Before inversion, each dataset was subsampled with different downsampling methods. We masked the island in the southwest of the epicentre to avoid possible unwrapping errors and used the quadtree method [34] to downsample the ALOS-2 ScanSAR measurements. For the remaining four offset sets, only pixels around the ruptured fault were extracted to retain a high signal-to-noise ratio. We then used the quadtree method to downsample Sentinel-1 and ALOS-2 SM offsets, and applied a subsampling scheme depending on the distance to the fault trace to Sentinel-2 measurements. The sampling interval for pixels within 5 km was only one, then increased to two for those pixels from 5 km to 20 km, and finally increased to four for those with a distance farther than 20 km.

3.1. Inversion Method

Assuming an elastic homogeneous half-space with a Poisson ratio of 0.25, a two-step inversion strategy can be used to invert for the fault geometry and the slip distribution. The first step is a nonlinear inversion which estimates the fault geometry by assuming a uniform slip on a rectangular fault plane [35]. The second step is a linear inversion to solve for the finite-fault slip distribution in a least-square sense. In this study, we fixed the strike angle of the fault based on its surface rupture, and then used multipeak particle swarm optimization (M-PSO) with a hybrid minimization algorithm [36] to search for the other fault geometry parameters including the dip angle, the width, length, and depth of the fault plane, and the upper boundary of the ruptured fault. After the determination of the optimal fault geometry, we discretized the fault plane into rectangular sub-patches (2 km × 2 km). Since the fault dip angle from the uniform model may not be the optimal for spatial-variable slip distribution [37,38], we refined the dip angle based on the relationship between weighted model misfits and multiple dip angles. To limit variations in the slip solution, the Laplacian smoothing was applied with a smoothing factor determined from a trade-off curve between the slip weighted misfit and roughness (as in Figure 4a).

3.2. Weighting Strategy

Following the above inversion strategy, we initially performed individual inversions using each of the five datasets (i.e., ALOS-2 ScanSAR interferometric phases, Sentinel-1 range offsets, ALOS-2 SM azimuth offsets, Sentinel-2 north–south offsets, and Sentinel-2 east–west offsets). The residuals of individual inversion were used to estimate the initial weights of these datasets in the joint inversion, which followed a ratio of 1:0.35:0.16:0.06:0.06. We then used M-PSO to obtain an optimal fault geometry with a length of 68 km, a width of 16 km, a dip angle of 89°, and a strike angle of 315°.

Then, we fixed the fault plane to be 100 km × 20 km and performed the joint linear inversion. We determined the Laplacian smoothing factor as 2.0, where the fault model fit the observations well and exerted a relatively small roughness, regardless of the setting of the dip angles (Figure 4a). We further refined the dip angle to be 88.5°, which corresponded to a least model misfit (Figure 4b). To better balance the contributions of multiple datasets, we iteratively updated the weight ratio using the residual root mean square (RMS) of each dataset. The square of the residual RMS reciprocal was calculated to set the new weight ratio. As shown in Figure 5, after three iterations of joint inversion, the data weight proportion and residual RMS proportion of each dataset tended to be equal, and the final weight ratio was 1:0.14:0.07:0.04:0.05.

4. Results

Figure 6 shows the inverted slip distribution on the fault plane from the joint inversion and individual inversion. It can be found that compared to the slip distribution from joint inversion (Figure 6a), the slip distribution inverted from ALOS-2 ScanSAR data only (Figure 6b) can hardly illustrate the onshore surface rupture (12~24 km along strike), while that inverted from the Sentinel-1 (Figure 6c), ALOS-2 SM (Figure 6d), and Sentinel-2 (Figure 6e) offsets seem to overestimate the length of the offshore surface rupture. The checkerboard test was further used to examine the resolution of these slip solutions. We first generated synthetic ground observations based on checkerboard-like slip distributions (Figure 7a), and then inverted them for the slip solution from these synthetic observations using the same inversion method. Results showed that the joint inversion (Figure 7b) retrieved the fault slip slightly better than ALOS-2 ScanSAR only (Figure 7c) in the onshore part (0~31 km along strike), and performed much better than the other three datasets (Figure 7d–f), especially in the offshore part (>31 km along strike). However, the down-dip slip (>10 km along dip) of the offshore region was to some extent underestimated and smeared, since no deformation could be captured on the offshore by satellite sensing.

The results from joint inversion, which offers a comprehensive geophysical interpretation, suggests that the earthquake was mainly controlled by a left-lateral strike-slip component. The maximum strike-slip and dip-slip were 6.07 m and 0.49 m, respectively, with a rake angle of 4.6°. This preferred model indicates that the fault ruptured to the surface with a length beyond 50 km and a maximum surface slip of over 5 m. The main slip area extended 18 km along the down-dip and the peak sliding patches with a slip of 6.10 m were located at a depth of about 10 km, where the rupture propagated mostly along the strike direction. The estimated geodetic moment (M0) was 1.03 × 1020 N·m, corresponding to a magnitude of Mw 7.31.

The predicted displacements from the slip model fit the observations well. The RMS of the misfits were 8.7 cm, 23.8 cm, 33.8 cm, 45.6 cm and 41.2 cm for ALOS-2 ScanSAR interferometric phases, Sentinel-1 range offsets, ALOS-2 SM azimuth offsets, Sentinel-2 north-south (NS) offsets and Sentinel-2 east-west (EW) offsets, respectively. We further calculated the model predictions on the ground within each dataset and compared them with original observations. As shown in Figure 8, the deforming pattern of models basically followed that of the five datasets, and the residuals in the Sentinel-1 range offset map had the most uniform distribution, although not the smoothest. In the residual map of InSAR data from ALOS-2 ScanSAR (Figure 8c), there were two asperities on the west side of the fault trace, possibly caused by unwrapping errors in the original data. However, the causes of unwrapping errors were different. The island asperity was caused by pixel isolation while the other one is caused by low coherence. The relatively large residuals of ALOS-2 SM and Sentinel-2 (Figure 8i,l,o) were predictable because of their data quality, which is limited by resolution. Another reason for relatively large residuals in the edges of these images is that only the near-field observations were used in the joint inversion. Moreover, the atmospheric artefacts and early post-seismic deformation may also contribute to the residuals.

5. Discussion

Due to heavy vegetation in the study area, the C-band SAR interferometry hardly contributes to the mapping of the coseismic deformation field. The L-band ALOS-2 SAR images generate better interferograms, but still suffer coherence loss in areas along the ruptured fault. Such decorrelation may be caused by heavy vegetation or a large deformation gradient. SAR and optical offset fields that can be calculated from Sentinel-1 and Sentinel-2 data are free from coherence loss. Despite the lower precision, offset maps can provide enough near-field data constraints for slip inversion. Moreover, the surface trace of the ruptured fault is consistently displayed in these offset maps, which fixes the strike direction of the fault model. Our study has demonstrated the possibility of combining them for a joint inversion, which is beneficial in areas with limited observations. We also intended to jointly use satellite radar interferograms both in descending and ascending, but the deformation component in the LOS of the ascending ALOS-2/Sentinel-1 interferograms is limited since the strike direction of the ruptured fault conflicts with the satellite azimuth. In addition, there are no ascending ALOS-2 ScanSAR images available spanning the coseismic period.

For joint inversion with multiple data, the determination of the weight of each dataset imposes a large impact on the slip solution. If we ignore the data quality and weight of each dataset equally, a slip distribution with different motion patterns would be generated even with the same fault model and smoothing factor. The relative weight ratio estimated from our iterative weighting method show that the contributions of SAR and optical offsets were lower than 25%. Although our study combined multiple satellite data, most data points were still distributed on one side of the fault, which makes the model constraint unbalanced. The checkerboard test also verified that these satellite data hardly provided good constraints for the offshore slip. This is probably the reason for the lower geodetic moment magnitude.

The slip distribution results indicate a lengthy rupture on the surface. To evaluate the effect of such fault ruptures on the surrounding seismogenic environment, we calculated the Coulomb failure stress (CFS) changes resulting from this event. The CFS changes representing the transfer of stress have been widely used to characterize the evolution of seismicity and quantify the triggering effect of medium/large earthquakes [39,40,41,42]. We estimated the CFS on an optimally oriented fault at a depth of 10 km with a friction coefficient of 0.4. The shear modulus was assumed as 33 GPa for a Poisson’s ratio of 0.25. The aftershocks collected from the USGS were superposed on the map of CFS changes (Figure 9). Although the number of recorded aftershocks were few, most aftershocks (~70%) were located in the area with increased CFS. The map also shows that the 2019 New Ireland earthquake increased CFS by more than 5.0 bars along the strike of the ruptured fault, except for the northwest and southeast edges of the faults.

The 16 November 2000 Mw 8.0 event shares the same fault structure and initiates from a similar location with this event. Additionally, their aftershocks do not delimit a linear zone of seismicity as expected for near-vertical strike-slip events [5]. However, the aftershocks of the 2000 Mw 8.0 event were clustered to the east of the Weitin fault, while those of this event were more scattered around the fault. The shallow rupture of the 2000 Mw 8.0 event may have reactivated many older subduction fractures in the upper plate [5], which can explain the static stress changes triggering the two 2000 Mw 7.8 thrust earthquakes [7]. However, in this event, no strong subduction-related activities to the south of the Weitin fault were triggered. We could not deny that the seismicity history repeated on the Weitin fault with a similar rupture mechanism, but the triggering effect of this event was more limited along the plate boundaries and associated structures when compared to the 2000 earthquake sequence, which reflects the uncertainty of the interaction between the Weitin fault and the subduction trenches.

6. Conclusions

In this study, we used multiple remote sensing observations including ALOS-2 ScanSAR interferometric phases, ALOS-2 SM azimuth offsets, Sentinel-1 range offsets, and Sentinel-2 offsets to map the coseismic deformation field and to invert for the slip distribution. First, the surface trace of the ruptured fault was obtained from SAR (Sentinel-1 and ALOS-2 SM) and optical (Sentinel-2) offset fields, which agreed well with each other. Then, an iterative weighting method based on the residual RMS of each dataset was developed to better balance the contributions of multiple datasets, with the Laplacian smoothing factor as 2.0 and the fault dip angle as 88.5° being successively determined along with the best fitting slip model. The joint inversion technique is beneficial to areas where limited observations are available due to heavy vegetation (e.g., tropical regions) or ice covers (e.g., polar or plateau regions).

The preferred slip distribution suggests a nearly pure left-lateral strike-slip motion on the Weitin fault and a surface rupture length of ~50 km. The maximum slip was ~6.10 m with a geodetic moment of 1.03 × 1020 N·m, corresponding to a magnitude of Mw 7.31. The distribution of aftershocks and Coulomb failure stress changes by this event show that ~70% of aftershocks were located in the area with increased stress and almost no strong aftershocks were triggered in the subduction zone to the south of the Weitin fault.