Sunda and Sumatra Block Motion in ITRF2008

The Southeast Asia region is mostly surrounded by active subduction zones in which the Australian plate, the Indian plate, and the Philippine Sea plate submerges beneath the continental plates and blocks. The Sunda block covers the large part of the Southeast Asia region, which comprises of Indochina, the South China Sea, the northeastern part of Sumatra, Borneo, the northern part of Java, and the shallow seas in between. We collect the GPS data in the whole Southeast Asia region for the period from 1994 to 2016, and process the original carrier phase data of GPS using GAMIT/GLOBK 10.6 to obtain the velocity field in the International Terrestrial Reference Frame, ITRF2008. The velocity field thus obtained is utilized to update the Euler rotation parameters of the Sunda block in ITRF2008, and model the long-term slip rates between the adjacent plate and blocks. In this study, we model the Sunda block and the Sumatra block together with the Australian plate by using TDEFNODE. The estimated Euler pole parameters of the Sumatra and Sunda blocks are estimated as their locations at (37.4°S, 106.8°E) and (46.2°N, 89.4°W), respectively, and their angular velocities of 0.371°/Myr clockwise, and 0.327°/Myr counter clockwise, respectively. These parameters result in the slip rate of the Sumatra fault with magnitude of ~9 mm/yr.


Introduction
The Southeast Asia region is mostly surrounded by active subduction zones in which the Australian plate, the Indian plate, and the Philippine Sea plate submerges beneath the continental plates and blocks. The largest part of the area is called as the Sunda block which covers Indo-china, some part of Sumatra, Borneo, some part of Java, and the shallow seas located in between (Simons et al., 2007). In the north, the Sunda block is bounded by the collision zone between the Eurasian plate. In the west, the Indian plate subducts beneath the Sunda block. In the south, the Sunda block is bounded by the Sumatran fault and the Baribis-Kendeng fault in Java Island. In the east, the West Makassar block and the Philippine Sea plate is the boundary of the Sunda block. Hall and Morley (2004) suggested that the Sunda block presently moves as a coherent lithospheric block although its geological origin must not be monolithic. The interior of the Sunda block has been characterized by a very low seismicity consistent with strain rate lower than ~7 × 10 -9 /yr estimated by Simons et al. (2007). By using high precision GPS measurements, the crustal motions can be accurately determined. Many GPS studies recently published have defined the Sunda block as an independent block (Bock et al., 2003;Simons et al., 2007). However, there are still significant discrepancies in the location of its boundaries and its rotation pole. The latest study of the Sunda block motion (Simons et al., 2007) yielded the Euler rotation parameters in the International Terrestrial Reference Frame, ITRF2000 (Altamimi et al., 2002). In this study, the GPS velocity field thus obtained is utilized to update the Euler rotation parameters of the Sunda block in ITRF2008.

Methodology
In this study, GNSS data from several sources for the period from 1994 through 2016 are used: 68 IPGSN stations for the period from January 2007 to the end of 2016, 28 International GNSS Service (IGS) stations for the period from January 2007 to December 2016, 7 SuGAr stations for the period from January 2002 through December 25, 2004 (before the Sumatra-Andaman Islands Earthquake, M9.1).

GNSS Data Processing Strategy
GNSS carrier phase observables from each station are used to estimate the 3-dimensional Cartesian coordinates and the integer phase ambiguities by applying loose constraints to the parameters: 6 orbital elements, semi-major axis, eccentricity, inclination, longitude of ascending node, argument of perigee, and mean anomaly are fixed and the Earth orientation parameters (pole position and rate and UT1 rate). The atmospheric zenith delay and its horizontal gradients at each station, the IGS elevation dependent antenna models, the ocean loading and solid Earth tides are given from the models. The GNSS stations were grouped into several networks with at least four overlapping tie stations to combine the individual network solutions. The maximum number of stations per network was set to 20. This step produced the RMS (root mean square of residuals per degree of freedom) of coordinate components for each station below 15 mm and the normalized RMS for all solutions are ~0.2.
The unconstrained parameters including the covariance matrices from GAMIT solutions were input into the GLOBK Kalman filter, which combines those loosely constrained parameters and applies uniform constraints to all parameters to derive final combined solutions in the reference frame ITRF2008 for the station coordinate estimates for each daily solution (Dong et al., 1998;Herring, 2015) and to estimate six components of a Helmert transformation (3 translation and 3 rotation). It has been done by using the stable stations in the list included in GAMIT software out of 28 IGS stations for each epoch. Basically, it was an iterative process, in which four iterations were used to eliminate bad stations and to calculate station weights for the reference frame stabilization (Nikolaidis, 2002). Each time series was cleaned using an outlier detection that excluded the outliers exceeding 95 % of confidence level (2 sigmas) and with uncertainties larger than 10 mm for the horizontal components and 20 mm for the vertical ones.

Time Series Analysis Method
In general, a time series of GNSS site coordinates consists of three major signals, inserseismic site velocities, seasonal variations, and earthquake-related signals including coseismic and postseismic deformation (Feng et al., 2015). The interseismic velocities can be described as a linear term before or even after earthquakes. The seasonal variations can be expressed by superposition of sinusoidal terms with constant amplitudes and may be attributed to some different sources such as gravitational excitation, thermal and hydrodynamics, local site errors, and model errors (Dong et al., 2002). Since the seasonal signals are modeled in the GAMIT solutions as the solid Earth tides and/or the atmospheric zenith delay to some extent, they should be already small. The largest impact of the seasonal signals exists in vertical component.
The coseismic signals can be represented by offsets in the time series. The postseismic deformation can be attributed to after slip, viscoelastic response, or poroelastic response (e.g. Feigl and Thatcher, 2006). Detailed studies are necessary to distinguish which source is dominant and it is beyond the scope of this study. Instead an exponential formula is used to represent postseismic deformation in the time series in order to quantify the pattern of postseismic processes, even though this formula cannot examine the driving mechanisms. Some offsets are found in the time series obtained at the stations nearby significant earthquakes especially those of the IPGSN stations. At least two major earthquakes, the M7.0 Java Earthquake of 2 September 2009 and the M8.6 off the west coast of northern Sumatra Earthquake of 11 April 2012 affected the time series in the data period. Beside of these earthquakes, the possibility of ongoing postseismic deformation due to the M 7.7 Java Earthquake of 17 July 2006 was examined.

Time Series Analysis Method
The GNSS time series represent the motion of the Earth surface. This motion can be related to the kinematics of block rotation and the locking on the block-bounding faults. The motion of blocks can be described by Euler poles in spherical coordinates. An Euler vector in a spherical coordinate system Ωi = (λp, φp, ω) gives the displacement at a surface observation point (λ, φ) caused by the rotation of a block i in the reference frame ITRF2008, where λp and φp are the longitude and latitude of the Euler pole, respectively and ω is its angular velocity. The horizontal velocity at a point in the i-th block is V i = Ωi × X, where X is the vector pointing from the center of the Earth to the surface point (λ, φ). A geodetic inversion code, TDEFNODE is used to model simultaneously the block rotations, homogenous strain rates, and locking distributions on block-bounding faults (McCaffrey, 2009).

Results and Discussion
We carried out the kinematic modeling of the Sumatra block, the Sunda block, and the Australian plate using the GPS data from 1994 to 2016. These blocks are separated into several faults known as the Sumatran fault and the subduction fault along the Sumatra trench (the Sumatra subduction fault). The Sumatran fault is the boundary between the Sunda block and the Sumatra block, and running along the southwestern coast of Sumatra Island. The Sumatra subduction fault is the boundary between the Sumatra block and the Australian plate, where the activity of the interplate thrust earthquakes has been very high. I used the slab model from Hayes et al. (2012) to represent the Sumatra subduction fault. In order express the slab surface, several nodes were located every 200 km along the isodepth contours of the slab model with 20 km depth interval. Coupling factors on these nodes are set to be free parameters except for the nodes on eastern and western edges where the coupling factors are set to be fully creeping considering the depth of the seismogenic zone in this area. In the inversion, the elastic deformation is calculated by integrating over small patches assumed in the regions between the nodes. Each small patch has a length of 10 km along strike direction and a width of 5 km along dip direction in this study.
Smoothing of the locking is applied by imposing a damping factor of 5 × 10 10 derived from the trade-off curve result.
The Sumatran fault has set to be full coupling from 0 km to the depth of 15 km and nine nodes along each iso-depth contour with an interval of 5 km are located based on the bending of the fault trace. The rest of boundaries are assumed to have no relative motion from one another, and to be fully creeping. In this first model,  Figure 1 shows the fault slip rates decomposed into the fault-parallel and the fault-normal motions. The fault-normal convergence rates across the Sumatra trench increase from 46 mm/yr in the west to 54 mm/yr in the east. These estimates are generally similar with the previous study by McCaffrey (2009) who noted that the convergence rate at the Sumatra trench is likely to be 40-50 mm/yr. The fault-parallel velocities along the Sumatra trench show an increase from 18 mm/yr to 32 mm/yr from the east to the west. Right lateral slip rate is found along the Sumatran fault as ~9 mm/yr. These rates are slower than that predicted by Bradley et al (2017). Small compressions exist along the Sumatra fault at rates of 0.3 -0.5 mm/yr. The locations of the Euler pole of the Sunda block estimated in this study and the previous ones are shown in Figure 2. The smallest uncertainty of the Euler pole parameters for the Sunda block is delivered in Simons's study. However, the uncertainty of the Euler pole parameters in this study are acceptable considering the reweighting that has been applied in this estimation.

Conclusion
The Euler pole parameters of the Sumatra and Sunda blocks are estimated showing the locations at (37.4°S, 106.8°E) and (46.2°N, 89.4°W), respectively, and their angular velocities of 0.371°/Myr clockwise, and 0.327°/Myr counter clockwise, respectively. These parameters result in the slip rate of the Sumatra fault with magnitude of ~9 mm/yr. The distribution of slip deficit rate on the Sumatra subduction fault is characterized by fully locking in the area between 97°E and 102°E in the depth range from 0 km to 60 km, which overlaps with the rupture area of the 1797 Sumatra earthquake (M8.4) and the 2007 Mentawai earthquake (M 8.4).