Geodetic observation of strain accumulation in the Banda Arc region

Abstract The Banda Arc region has produced several large destructive earthquakes, some of which have been followed by tsunamis. To better understand the earthquake potential in this area, we performed a comparison between geodetic and seismic moment rates. Data were collected from 110 continuous and campaign GPS stations observed for approximately ten years. The results show that the derived velocity field indicates that the Banda Arc deformation is characterized mainly by crustal shortening caused by the interaction of the Australian, Pacific, and Philippine Sea plates. Meanwhile, the contraction strain pattern dominates the Banda Arc area except around Papuan Bird’s Head. Areas with high strain rates have a history of significant seismicity, such as the Flores-Wetar Back Arc, the area around Ambon, and the Papuan Bird’s Head. The ratio of the geodetic moment rate to the seismic moment rate in the Banda, Bird’s Head, South Sulawesi, and Sumba zones are ∼1.5–7.0, indicating a moment deficit rate. The moment deficit rate provides an equivalent earthquake potential of Mw 7.7–8.1. This potential may be related to an aseismic deformation or stress accumulation, the under-sampling of long-term earthquake rates within the seismic catalogs, or a composite of these factors.


Introduction
Eastern Indonesia is the region with the most complex tectonic configuration because of the convergence of four active plates, namely, the Sunda, Australian, Pacific, and Philippine Sea plates (Hamilton 1979;Nishimura and Suparka 1990;Milsom 2001). The convergence of the Australian plate and Sunda Block is partitioned between a megathrust and a continuous zone of back arc thrusting extending over 2000 km studied. The continuous geodetic strain rates is suitable for identifying areas with high deformation rates and can distinguish the type of tectonic structure (Hackl et al. 2009). Strain rates can also be translated into a geodetic moment rate by comparing the seismic energy released within the tectonic area. Several studies have conducted such a comparison to explain the potential seismic hazard of certain regions around the world, such as in some Mediterranean areas (D'Agostino et al. 2009;Chousianitis et al. 2015;Sparacino et al. 2020), the Himalayan Arc (Bungum et al. 2017;Sharma et al. 2020), the Zagros Fold Thrust Belt region (Palano et al. 2018), Northeast Tibet (Pan et al. 2020), and Egypt (Sawires et al. 2021). A similar study was also carried out in the Nusa Tenggara and Banda regions by Robiana et al. (2018). However, it was poorly explained due to the ongoing deformation process in the Banda region.
In this study, we combine geodetic and seismic data analyses to understand the earthquake potential in the Banda Arc region. First, we processed continuous and periodic GPS data to obtain velocity fields from 110 stations across the study area. Then, we modeled the present-day deformation of the study area by analyzing the accumulated strain rate derived from interpolated GPS velocity, which excluded the coseismic offset of the earthquakes. Next, we converted the strain rate into the geodetic moment rate. We also used several seismic catalogs which provides earthquake locations, occurance times, magnitudes, and depths of the seismicity in our study area. Finally, the seismic moment rate of the earthquake catalog was calculated and compared with the geodetic moment rate. Comparing the geodetic and seismic moment rates can explain current tectonic activity in the Banda Arc region and its implications for potential seismic hazards, which is the main contribution of this study.

GPS data and processing
We used GPS data collected from 49 permanent and 61 periodic stations adequately distributed across Eastern Indonesia, which were provided by the Geospatial Information Agency of Indonesia (BIG) (https://srgi.big.go.id). The GPS dataset covered different periods for each station, ranging from 2.5 to 9 years with an average duration of 5.4 years. GPS data were processed using GAMIT/GLOBK software version 10.7 (Herring et al. 2018). GAMIT's main inputs are GPS station observation data in RINEX format, orbital or g files (sp3/g-files) and Global Navigation files (SOPAC, sopac.ucsd.edu; CORS, www.ngs.noaa.gov). Furthermore, the differential phase observation method was used to estimate the station coordinates, tropospheric delays, satellite orbit parameters, and phase ambiguity estimated for each observation day. To estimate the daily position time series in the 2008 International Terrestrial Reference Frame (ITRF2008) (Altamimi et al. 2011), the loosely constrained GPS solutions from GAMIT were combined in GLOBK with the daily solutions of 16 International GNSS Service (IGS) stations.
Subsequently, the GPS time series obtained at each station was edited and checked for quality visually and computationally. Outlier data were eliminated by setting a two sigma (2r) limit from the average of all the data (m) at the 95% confidence level to obtain the interseismic velocity of the daily coordinate linear trend. The annual and semi-annual cycles and offsets caused by tool changes or earthquakes (Nikolaidis 2002;Raharja et al. 2016) were identified and eliminated by calculating the interseismic velocity values at each GPS station. Figure 2 shows four examples of GPS stations' daily horizontal and vertical position time series free from seasonal variations, offsets, and outliers. These stations, namely BANI, CSAU, CFAK, and CREO, are in Banda Neira (Maluku), Saumlaki (Maluku), Fak-Fak (West Papua), and Reok (East Nusa Tenggara), respectively. The GPS data used in this analysis were obtained from 2010 to 2018, with some covering a shorter period due to minor data gaps during the study period. All the stations in Figure 2 move in the N-E direction, except CFAK, which moves toward the N-W due to the very small influence of the Australian plate. In the next processing, the horizontal and vertical components were used to determine the strain tensors and time series analysis regarding offsets caused by equipment changes.

Strain rate analysis
Strain rate analysis is often used in deformation studies to identify the active areas with significant strain accumulation (Sagiya et al. 2000). This process was conducted by implementing the spline interpolation method of geodetic data (Hackl et al. 2009). This method can convert a GPS velocity field into a continuous strain rate covering the entire study area. Although the spline interpolation method cannot determine the fixed value of the strain, it is still beneficial for identifying tectonic areas with significant strain rates. This process was also applied to the latitude and longitude components of the geodetic velocity data with a grid size of 30 Â 30 km ($0.3 degrees), which covers approximately half of the mean station distance, specifically in inland areas. In addition, the tension (T) algorithm was applied when interpolating the velocity data, with T ¼ 0.3 (Wessel and Bercovici 1998). The spline interpolation method was applied to the easting and northing velocity components, obtained from the derivation of each component. This produced four continuous planes that represent the strain rate tensor as follows (Turcotte and Schubert 2014): where i and j denote the easting and northing components, with v at point x: Similarly, the rotation rate asymmetry tensor can be calculated using the following equation (Turcotte and Schubert 2014): After determining all the tensor components, the strain derivative products were estimated, including the principal strain rate, maximum shear strain rate, dilatation rate, and second invariant rate (Turcotte and Schubert 2014): The eigenvalues (k 1 and k 2 ) represent the magnitude of the strain, with positive and negative values denoting expansion and shortening, respectively. The maximum shear strain was computed as follows: Meanwhile, the direction for a maximum shear strain rate of 45 from the eigenvector direction, which corresponds to the largest eigenvalue, was determined using the following equation: The second-order strain rate tensors (i.e., dilatation rate (e dilatation ) and second invariant (I)) were estimated as follows: The principal strain rate was used to evaluate the dominant strain mechanism that occurs in an area, while the maximum shear strain rate was utilized to identify the area of deformation due to a strike-slip fault. Furthermore, the dilatation rate can be used to analyze the horizontal deformation caused by thrust or normal faults, while the second invariant estimated the magnitude of the strain rate (Sagiya et al. 2000;Hackl et al. 2009;Ashurkov et al. 2016;Susilo 2017). Figure 3 shows that the velocity values of the GPS stations in Eastern Indonesia vary widely in the ranges À87.32 to 18.05 mm/year and 0.13 to 81.92 mm/year for the easting and northing components, respectively. In general, we find that the GPS velocity vectors rotate counterclockwise, increasing in value toward the east. This is in line with the findings of the previous studies by Koulali et al. (2016) and Susilo et al. (2016).

GPS velocity field
The velocity values of the forearc of the Banda region increase towards the north, which accommodates the dominant convergence between the Australian plate and Sunda Block. By contrast, the velocity vectors in West Papua move to the northwest due to the force from the Pacific plate to Eastern Indonesia in the westward direction.
The velocity patterns in Seram Island are divided into two directions, namely N-E and N-W. These patterns suggest that this island has a strike-slip faulting mechanism (Meilano et al. 2021).

Strain rates
The results of geodetic strain analysis are demonstrated as maps of the principal strain, maximum shear strain, dilatation, and second invariant rates. However, irrespective of the inability of these results to estimate the absolute strain values, they are useful to picture the area accumulating a contemporary crustal deformation process quickly. These results are then utilized as a starting and supportive point for further analyses and other geological and seismic data, such as quantifying the geodetic earthquake moment rates and comparing them with the seismic energy release (Ward 1994;D'Agostino et al. 2009;Mazzotti et al. 2011;Chousianitis et al. 2015;Walpersdorf et al. 2015;Palano et al. 2020;Sharma et al. 2020). The investigation of this comparison investigation is also presented in Subsection 3.3 of this paper. Figure 4a shows the pattern estimate and value of the principal strain rate. The expansion and contraction values in Eastern Indonesia vary from 65 to 330 nanostrain/year and from À299 to À3 nanostrain/year, respectively. The strain value in Papuan Bird's Head is higher than those in the other areas and dominated by the expansion patterns around Nabire. In contrast, the shortening pattern dominates in the Nusa Tenggara and Seram Island areas. The strain rates are low in the area between Southern Kalimantan and Eastern Sulawesi and south of Nusa Tenggara, indicating that these regions are part of the Sunda Block and Australian plate, respectively (Kreemer et al. 2000). Figure 4b shows the maximum shear strain rate pattern in Eastern Indonesia, ranging from 0.4 to 200 nanostrain/year. This pattern indicates that the area experiencing the highest shear strain rate is around Cendrawasih Bay in Papua ($200 nanostrain/year). This fact is in line with the mapping of earthquake sources by the National Earthquake Study Center (PuSGeN 2017), which estimates several strike-slip faults in the Cendrawasih Fault, Randhawa Fault, and Wapoga Fault with a slip rate of 17-35 mm/year. In addition, the Nabire, Ambon, and Maumere areas experience shear at a value of $150 nanostrain/year. The significant shear rate in the Ambon region is suspected to be one cause of the M6.5 strike-slip fault earthquake in September 2019 (Meilano et al. 2021). Despite the Flores Back Arc Thrust upswing in the northern area of Flores Island, the shear strain rate is still enormous, at around 120 nanostrain/year. We suggested that this phenomenon is related to the possible presence of an unmapped fault system in the area.  Figure 4c shows that the dilatation rate ranges from À200 to 250 nanostrain/year. The red and blue indicate areas experiencing expansion and shortening. The resulting dilatation pattern indicates that Eastern Indonesia, specifically in Banda Sea and Nusa Tenggara, experiences significant and varied deformations due to thrust faults. A large negative dilatation ($150 nanostrain/year) pattern is also seen in the north of Flores Island, confirming that the area is experiencing a significant accumulation of energy that could produce earthquakes in the future. In contrast, high values of positive dilatation are observed in Cendrawasih Bay, indicating that this area is undergoing an extended deformation pattern. Relatively small rates of expansion (<50 nanostrain/year) are also found in the south and southeast of the left-hand and righthand of Sulawesi, respectively, as well as in the south of Flores Island, which can be attributed to local expansional deformation. Figure 4d presents the analysis of the second invariant rate for Banda Arc, which expresses the strain in the region. We can see that the closer the second variant is to the east, the larger is the strain magnitude, at a range of 0-200 nanostrain/year. Similar strain patterns to the maximum shear strain map are observed in the northern areas of Maumere-Flores Island, Ambon and Seram Island, Nabire, and Cendrawasih Bay. These areas manifest the most significant strain magnitude in Eastern Indonesia of >100 nanostrain/year.

Computing the geodetic and seismic moment rates
To compare the spatial distribution of the geodetic and seismic moment rates consistently, we must divide the study area into several seismogenic zones. Previous studies have used a uniform grid approach to define their seismogenic zones (Palano et al. 2018;Pan et al. 2020). Here, we divide Banda Arc into five zones, considering the limitation of the GPS distribution, as well as summary from previous studies (Argus et al. 2011;Koulali et al. 2016;Susilo 2017), which have defined these zones based on the main tectonic features, seismic data, and velocity patterns with respect to the Australian plate. The segmentation is coded BAND, BRHD, SSUL, SUMB, and TIMO, which represent the Banda Sea, Papuan Bird's Head, South Sulawesi, Sumba, and Timor areas, respectively ( Figure 5).
To compare the accumulated moment budgets with the released moment rates, the geodetic moment rate is calculated using Equation (9) following Ward (1994), which sum the largest eigenvalue of the strain rate tensors from each zone. Some zones that have a dense network of observations have better strain quality than areas that have a less dense network (offshore areas). Hence, to avoid the interpolation uncertainty in biasing the moment rate analysis, we exclude the interpolated strain values located at a distance greater than 2 (4 grid cells) from the nearest observation site (Ashurkov et al. 2016).
where M 0g is the geodetic moment rate, l is the rigidity, H is the depth of the seismogenic layer at which most of the earthquakes occur, A is the surface area, and _ e 1 and _ e 2 are the absolute eigenvalues of the strain rate tensor of the grid cells within each seismogenic zone. A l value of 3.5Â10 10 N/m 2 is used in this study (Kreemer et al. 2000). Meanwhile, the seismogenic depth is adopted from the global crustal model CRUST1.0 (Laske et al. 2013) by averaging the gridded depths that belong to each zone. Table 1 lists the geodetic moment rates for all the seismogenic zones of the study area, along with the variables used.
In contrast, the seismic moment rate M 0s is estimated by accumulating the moments of the earthquakes within each seismogenic zone. We use the seismic catalog that cover the time span between 1960 and 2018, with 13,703 events (including small earthquakes) in the area extending from 13 E to 136 E longitudes and from À14 S to 2 N latitudes. The seismic catalog is compiled from the ISC-EHB bulletin (http://www.isc.ac.uk/isc-ehb/; last accessed on 12 May 2022) within the period from 1960 to 2008 and the catalog developed by Supendi et al. (2020) within the period from 2009 until 2018. To unify the different magnitude scale of the ISC-EHB bulletin (mb, Ms, MD, and ML) into the moment magnitude (M w ), the following magnitudescaling relationships (Sawires et al. 2021) are used to convert them: Furthermore, the earthquake moment magnitudes (M w ) of the catalog are converted into scalar seismic moment rate (M 0 ) using the conversion formula of Hanks and Kanamori (1979). Then, we implement the formula of Molnar (1979) to estimate the seismic moment rate as follows: where M 0max is the moment of the largest earthquake that may occur within each seismogenic zone, defined from the calculated M max of PuSGeN (2017) or the historical M max of Utsu's catalog (Utsu 2002). Furthermore, a and b are the seismicity level and slope of the Gutenberg-Richter's recurrence relation, respectively: logNðMÞ ¼ a À bM), where N ðMÞ is the cumulative number of earthquakes with a magnitude above M; c ¼ 1.5 and d ¼ 9.1 are used from the Kanamori (1977) relation. The corrected values of the seismic moment rate for each zone in the study area are obtained from this calculation. The resulting values are arranged in Table 2, along with the a À b values, magnitude of completeness (Mc) of the catalog, and options of M max : To see whether the comparison between the geodetic and seismic moment rates dependens on the conversion method used, we calculate the moment rates using other methods ( Savage and Simpson 1997;Hyndman and Weichert 1983). The results show that the moment rate ratios are relatively consistent for all seismogenic zones. Therefore, we suggest that the comparison between the geodetic and seismic moment rates (moment rate ratio) is independent of the method used. The comparison results can be seen in the Appendix.

Comparison of the geodetic and seismic earthquake moment rates
The comparison between the geodetic and seismic quantification is presented as a ratio of the geodetic to the seismic moment rates in Table 2 and Figure 6. Based on this analysis, the BAND zone has high values for both the geodetic and the seismic moment rates, which are around 2 to 10 times higher than the others. This means that this region experiences high-frequency seismic activities as well as a large amount of deformation activity. BAND along with SSUL show significant levels of agreement, with a moment ratio of about two. Furthermore, satisfactory agreement is observed in SUMB, with a moment ratio relative to unity, suggesting that seismic activities almost completely release the geodetic budget.
Significant discrepancies between the moment rates are observed in TIMO, with a ratio below unity. The seismic moment rate is about twice as high as the moment accumulation in this zone. Some relative explanations to this phenomenon are because of the lack of geodetic data over the region, and/or the measured crustal deformation released by several moderate to large shallow events at depths < 40 km.  One of the events was the devastating M w 7.8 Flores earthquake on December 12, 1992, followed by a tsunami of 2.5-3.2 meters above mean sea level, causing 87 fatalities and 1400 deaths (Tsuji et al. 1995). Another possible reason why TIMO has the lowest moment rate ratio is that this region's relative motion gradually decreases along the Timor and Tanimbar Troughs, from 32 ± 2.0 mm/year to 1.0 ± 1.7 mm/year (Koulali et al. 2016). This decrease in rates is also related to the continental-arc collision beneath the Timor and Tanimbar Troughs, which is portrayed by the low strain accumulation in these zones in Figure 4. In contrast, BRHD has the highest of moment rate ratio, with the GPS-based moment rate higher by a factor of seven than the earthquake-based moment rate. This comparison generally indicates that the area is still accumulating the energy budget needed to be released by seismic activities, has not witnessed large earthquakes over the last few centuries, the long-term earthquake rates are under-sampled within the seismic catalog used, long-term aseismic crustal deformation has occurred, or a combination of these factors.

Earthquake potential
This study also evaluates the moment deficit by subtracting geodetic moment rate from the seismic moment rate in each seismogenic zone (Ward 1998;Li et al. 2018;Sharma et al. 2020). Considering the hypothesis that a total moment deficit rate is an earthquake potential that will be released as one or more earthquakes in the future, we estimated the maximum magnitude needed for the seismic activity to completely fulfill the discrepancy between the geodetic and seismic moment rates (Chousianitis et al. 2015;Bungum et al. 2017;Sharma et al. 2020). Our results obtain the potential of future significant earthquakes with a maximum magnitude equal to M w 8.1 and M w 7.8 over BAND and SSUL, respectively and M w 7.7 in BRHD and SUMB.
We find evidence of earthquake potential in the event of the 2019 M w 6.5 Ambon Earthquake, which ruptured within the area of BAND. This event was followed by hundreds of aftershocks with magnitudes M L < 5 at depths of < 30 km (Meilano et al. 2021). The earthquake potential of SSUL was also recently proven by the occurrence of two significant earthquakes that ruptured in Mamuju-Majene (M w 6.2) and East Nusa Tenggara (M w 7.4) in 2021. The epicenter of both earthquakes was located right over the high strain rates area, as imagined by the map of second invariant rates in Figure 4d. Furthermore, the focal mechanism of the 2021 East Nusa Tenggara Earthquake was leftlateral strike-slip (USGS); however, none of its sources have been identified in previous studies. The high maximum shear strain and second invariant rates clearly show the possibility that a strike-slip fault exists in the NW-SE direction (see Figure 4b and d).
In terms of the high ratios for BRHD, BAND, and SSUL, these regions have several known active faults (Irsyam et al. 2020), as shown in Figure 6, leading to the implication of seismic hazard. There are at least three major earthquake sources in BRHD, namely the Sorong Fault

Conclusions
This study provides an interseismic strain analysis derived from the velocity field of both continuous and campaign 110 GPS sites observed from 2010 to 2018. The results show that the Banda Arc deformation is characterized mainly by crustal shortening caused by the Australian plate, Pacific plate, and the Philippine Sea plate in the north, east, and north, respectively. The contraction pattern strain rates dominate the Banda Arc area except around Papuan Bird's Head, which shows the highest expansion pattern. Meanwhile, the highest maximum shear strain rates are found around the Ambon-Seram and Flores-Wetar Back Arc regions that comprise major thrust structures (i.e., the Flores Back Arc Thrust and Seram Trough). These patterns are assumed to be associated with either the existence of unknown strike-slip structures or the post-seismic effect of past earthquakes. Lately, these strange patterns have been confirmed by the occurrence of two significant strike-slip events, namely the M w 6.5 Ambon Earthquake in 2019 and the M w 7.2 Flores-Selayar Earthquake in late 2021.
Furthermore, we compare the geodetic result with the seismicity catalog from 1900 to 2018 to understand the earthquake potential in Banda Arc by estimating the moment rate ratio of the five seismogenic zones around it. The strain rate analysis is converted into the geodetic moment rates of each seismogenic zone using Ward's formula, while the earthquake catalog is translated into seismic moment rates using the method proposed by Molnar. The comparison of the geodetic to the seismic moment rates shows a consistent value. For the BAND, BRHD, SSUL, and SUMB areas, the geodetic moment is higher than the seismic moment rate (with ratios within range 1-7), which indicates deficit moment rates equal to the M w 8.1, M w 7.7, M w 7.9, and M w 7.7 earthquakes for each zone, respectively. This deficit in seismic moment rates is primarily due to the possibility of no major earthquakes recorded in this region. By contrast, the seismic moment rate ratios are below unity (<1) in TIMO. In this zone, the lower moment rate ratio can be primarily attributed to the release of crustal deformation seismically. Table A1. The comparison of moment ratio using some methods: #1 and #2 are geodetic moment formulas from Ward (1994); and Savage and Simpson (1997), respectively. #A and #B are seismic moment equations from Molnar (1979); and Hyndman and Weichert (1983) Figure A1. The interpolated velocity for each horizontal component: (a) northing, (b) easting. Red is positive and blue is negative.