*2.2. Methods*

In NLLoc version 7.00 [33], the location algorithm adopts the inversion approach from Tarantola and Valette [34], and the event location procedure from Tarantola and Valette [34], Moser et al. [35], and Wittlinger et al. [36]. This algorithm generated the misfit function for each sample cell, estimated the probability density function (PDF), and determined the optimal hypocenter location. The location search process began by defining the search area in 3D by setting the boundaries of the x, y, z regions and then applying the 1-D AK135 velocity model [37,38]. The optimal hypocenter location for each event was determined by estimating the probability density function (PDF) using an oct-tree importance sampling.

The result of the hypocenter location from the NonLinLoc program with an equal amount of P- and S-waves used to process the travel time tomography inversion was performed using SIMULPS12 [39] to obtain the 3D structure of both Vp and Vp/Vs ratio and simultaneously hypocenters relocation by using pseudo-bending ray tracing [40], as shown in Figure S3. This program performs an inversion of Vp and Vp/Vs models from P- arrival times and S-P times directly, not Vp and Vs separately, due to differences in data quality between P- and S- arrival times [41,42]. The initial velocity model used in the inversion process is a 3D regional-global seismic velocity model (Vp) beneath Indonesia obtained from high-resolution tomographic imaging [23]. Meanwhile, the value of the Vp/Vs ratio was obtained from the Wadati Diagram (Figure S4), which is 1.76 and is consistently used for each layer.

We performed the inversion through three stages: Firstly, a damping analysis was performed to obtain the optimum damping parameter input by making variations in the data distribution, number, size, and grid nodes distance [41]. The damping value was obtained by comparing the variance model (km2/s2) and the variance data (s2) for both Vp and Vp/Vs values. The optimum value for Vp and Vp/Vs damping obtained was 150 and 130, respectively (Figure S5). This value was used to obtain a stable solution when performing the tomographic inversion and a checkerboard resolution test (CRT). Secondly, we ran a checkerboard resolution test using several grid spacings to define the best model parameterization used in this study. We decided to use a 60 × 60 km grid and a vertical grid size of 30 km after considering the results and dataset distribution. Thirdly, we performed the tomographic inversion using the damped least square where the norm of the perturbations model is weighted and combined with the misfit data squared. We set high residual data < 10 s at the early inversion. It is automatically discarded if a value greater than 10 s is found. We set the maximum number of iterations to 20 times and stopped at 15 iterations. The weighted Root Mean Square (RMS) residual value decreases each iteration until convergence. The stopping condition provided by an F-test indicates that variance reduction is insignificant after more iterations [39].
