Seismic Imaging of Lithospheric Structure Beneath Central-East Java Region, Indonesia: Relation to Recent Earthquakes

The Central and East Java region, which is part of the Sunda Arc, has an important role in producing destructive earthquakes and volcanic complexes in Indonesia as a result of the convergence between the Indo-Australian plate that subducts under the Eurasian plate. In this study, the 3-D seismic velocity structure (Vp, Vs, and Vp/Vs) of the crust and upper mantle was determined to reveal the presence of the subducted slab, volcanic sources, and seismogenic features in the Central and East Java region. We have manually re-picked P- and S-arrivals of 1,488 events from January 2009 to September 2017 recorded at 27 stations of the BMKG network. An iterative damped least-squares inversion method was applied to simultaneously calculate both hypocenter relocations and velocity structure beneath this high-risk region to a depth of 200 km. We then compare the tomographic results and seismicity to interpret structural features in the seismic zones. The subducted slab is dipping toward the north, imaged by the high-velocity regions with low Vp/Vs at depths of about 50–100 km. Low-velocity anomalies with high Vp/Vs above the slab at a depth of ∼100 km, imply the possible location of partial melting from slab dehydration. Fluids and melts are ascending to feed the volcanoes i. e., Merapi-Merbabu, Wilis, Pandan, Semeru, Bromo, and Ijen that also have similar low-velocity anomalies at 10–30 km depths, suggesting the presence of the sedimentary basin or magma reservoir. We also have redetermined the hypocenter location of the 2021 (Mw 6.1) Malang earthquake at 8.94oS, 112.45oE, with a depth of 59.7 km. The location error in the x, y, and z directions are 3.08, 6.39, and 11.91 km, respectively. This intraslab event with a thrusting mechanism is located in the high-velocity region and close to the intermediate-depth seismic clusters, which indicates the geometry of the oceanic slab. In the region of 1994 (Mw 7.8) Banyuwangi earthquake, we found a low-velocity anomaly at ∼ 50 km depth that might be associated with the presence of subducting seamount that is more hydrated than the surrounding slab. The slip over this subducting seamount caused the tsunamigenic earthquake.


INTRODUCTION
Plate subduction along the Java Trench plays a major role in producing destructive earthquakes and volcanic complexes for the region. Many earthquakes are situated in or around the dipping Wadati-Benioff zone and are extensively caused by crustal faulting. Moreover, the joint effects between heat convection in the mantle wedge and fluids from slab dehydration lead to partial melting beneath the volcanic front, producing arc magmatism and volcanism (Zhao et al., 2012). To locate the subduction slab and estimate the magma plumbing system beneath the region, reliable subsurface information is necessary. Seismic tomography is one geophysical tool that can provide information about the three-dimensional distribution of seismic velocities of the subduction zone. For example, the presence of cold subducted slab and hot magma inclusions in the crust and upper mantle is reflected in the strong contrast of seismic velocities, which have been observed in many subduction zones (e.g., Zhao et al., 1995;Nakajima et al., 2001;Hall and Spakman, 2002;Husen et al., 2003;Reyners et al., 2006;Tsuji et al., 2008;Zhao et al., 2012).
The study area covers the Central and East Java region, which is part of the Sunda Arc, along the subduction zone in southwestern Indonesia, where there is orthogonal subduction of the Indo-Australian plate beneath the Eurasian Plate. The convergence rate was estimated at ∼5.6 cm/yr in the western part of Java and increased to ∼6.5 cm/yr in the east (Koulali et al., 2017). In addition to the subduction zone, several active faults and volcanoes control the seismicity around the study area (Marliyani, 2016;Pusat Studi Gempa Nasional, 2017) ( Figure 1). The region has a seismic gap that has not shown significant seismicity between the southern coast of Java and the trench. This region holds the potential for future megathrust events, as Widiyantoro et al. (2020) suggested that the high slip deficit in Central and East Java would produce a Mw 8.8 earthquake with a return period of 400 years. Damaging events in the region include a tsunami earthquake that occurred near Banyuwangi (Mw 7.8) in 1994, caused by a slip over a subducting seamount (Abercrombie et al., 2001). In 2006, a moderate shallow earthquake occurred around Yogyakarta (Mw 6.3) due to the Opak fault activity. Recently, a Mw 6.1 intra-slab earthquake occurred in the south of East Java near Malang on April 10, 2021 with at least 334 houses damaged and seven fatalities, according to the report from National Disaster Management Authority (BNPB), Indonesia (https://www.bnpb.go.id/).
The geological setting includes a basement structure beneath Java dominated by ophiolitic and volcanic arc material that was stranded at the margin of the ancient subduction zone (Sundaland) during the Late Cretaceous (Simandjuntak and Barber, 1996;Smyth et al., 2007). From the middle Eocene to the middle Miocene, the Southern Mountain zone and a large basin called the Kendeng zone were formed along the present subduction zone. In the Eocene-Pliocene, the Rembang zone was formed in the northern part of the island. Subduction ceased in the Cretaceous and started again in the Late Miocene, forming the modern volcanic arc, with mainly strato-volcanoes in central Java (Smyth et al., 2005(Smyth et al., , 2007. Monogenetic volcanoes, where  Figure 3 and Figure 5. (B) Grid node parameterization used for the seismic tomographic inversion (plus signs). The colors of the circles represent earthquake focus depth, active fault lineaments (blue lines) and volcanoes (black triangles) (Marliyani, 2016;Pusat Studi Gempa Nasional (PuSGeN), 2017). the volcanism produces the formation of tens to thousands of edifices can be observed in a few places in Java with maar-type edifices (e.g., Ijen volcanic complex). However, lava dome-type polygenetic volcanoes are more dominated in Java (Marliyani et al., 2020).
Previous seismic tomography studies have been carried out in the Indonesia region, including the seismic velocity distribution beneath Central and East Java on both global (e.g., Hafkenscheid et al., 2001;Hall and Spakman, 2015;Zenonos et al., 2019) and regional scales using teleseismic events (e.g., Fukao et al., 1992;Puspito et al., 1993;Widiyantoro and Van Der Hilst, 1996). They suggested that the Java slab dips steeply at depths of 250-300 km. Both P-and S-wave tomographic models show a gap in the slab between 350-500 km and the absence of intense seismicity. The gap is correlated to a temperature increase, producing the lowvelocity anomaly at this depth (Hall and Spakman, 2015;Zenonos et al., 2019). On local scales, previous studies have been used to study small-scale structures of crust and upper mantle i.e., fault zones and volcanism. They observed strong low-velocity anomalies beneath volcanoes in Central Java (Koulakov et al., , 2016Wagner et al., 2007;Bohm et al., 2013;Luehr et al., 2013;Rohadi et al., 2013;Haberland et al., 2014;Zulfakriza et al., 2014;Ramdhan et al., 2019). They have found an extreme lowvelocity anomaly in the crust beneath Merapi and Lawu volcanoes (Merapi-Lawu anomaly (MLA)) with a high Vp/Vs ratio of 1.9. Based on these results, a model of magma supply to volcanoes in Central Java was proposed. The MLA is supposed to be an anomaly related to an area with high content of fluids and melts in the crust .
There are still few tomographic studies in East Java, in particular, models that can resolve the crust and upper mantle structure. One study by Martha et al. (2017) has successfully imaged the shallow structure beneath East Java, using ambient noise tomography (ANT) with both BMKG and portable seismographic networks. The study, which has a resolution depth of 12 km, found low shear wave velocities associated with basins and sedimentary layers, along with high velocities associated with igneous rocks and carbonates. The low-velocity anomalies correlate with the presence of volcanoes and oil/gas fields in East Java. The study results also concluded that the Kendeng basin has a thickness of 8-10 km. Another prominent feature in East Java is a sediment-hosted geothermal system, nicknamed Lusi (Lumpur Sidoarjo), which initially erupted in May 2006. The magma plumbing system feeding the Lusi eruption has been studied, using local earthquake tomography and ANT that resolve the structure to 25 km depth (e.g., Fallahi et al., 2017;Karyono et al., 2020;Lupi et al., 2022). They suggest that magma and fluids are coming from the volcanic complex around the region and migrating through the fault system to the sedimentary basin, promoting Lusi. Extensively, imaging the larger-scale structure may help the understanding of the relationship between subduction and volcanism in the region.
This study aims to determine the 3-D seismic velocity structure (Vp, Vs, and Vp/Vs ratio) of the crust and upper mantle and its relation to geological features around the seismogenic zone and volcanic arc of Central and East Java, Indonesia, which represents a high seismic hazard to the region and has a population of more than 73 million people.

DATA AND METHOD
We used P-and S-arrivals determined by Muttaqy et al. (2021-in review) that were manually re-picked from waveforms recorded at 27 stations of the Meteorological, Climatological, Geophysical Agency of Indonesia (BMKG) for earthquakes within the time period of January 2009 to September 2017 using Seisgram2K (Lomax and Michelini, 2009). The earthquake locations were then determined using the NLLoc program (Lomax et al., 2000) with the global velocity model of AK135 (Kennett et al., 1995) for the longitude and latitude ranges of 108 o -118 o E and 5 o -12 o S, respectively, with at least six high-quality phase picks and magnitude (Mw) >3. Most of the earthquakes are distributed in the south of the island, which results in a maximum azimuthal gap of greater than 180 o with the average hypocenter location errors are 5.6, 15.5 and 20.5 km in the x, y, and z directions, respectively. In total, there are 1,488 initial events with 11,192 Pand 11,192 S-phases which were used for the tomographic inversion ( Figure 1 and Supplementary Figure S1). The estimated picking errors for the observed P-and S-arrivals are 0.1886 s and 0.297 s, respectively. This results in better locations compared to the BMKG catalog, which uses Geiger method (Geiger, 1912) with global velocity model IASP91. The BMKG location, which use a fewer number of S-arrivals for quick determination have a fixed-depths of 10 km for the shallow earthquakes and larger residual times. A detailed description of the seismicity analysis for the events has been provided in Muttaqy et al. (2021-in review).
The choice of a particular model parameterization is always essential for assessing the quality of the 3-D tomographic images, considering the recovery of checkerboard resolution tests (CRT) over various grid size models. Several factors need to be considered, such as the size of features that can be resolved, distribution of available events and stations, and the computing limitations (Kissling et al., 2001). In this study, we used a 50 × 50 km 2 grid spacing horizontally and various grid spacings of 10-100 km in the vertical direction ( Figure 1B). The initial velocity model used in this study was the 1-D seismic velocity model of Central Java  with a constant Vp/ Vs ratio of 1.75, based on the plotted Wadati diagram of this study (Supplementary Figure S1).
For the local tomographic inversion, we applied the SIMULPS12 code (Evans et al., 1994), which performs an iterative damped least-squares inversion method to calculate both hypocenter relocation and velocity structure, simultaneously. P-arrival times and S-P times are inverted directly for Vp and Vp/Vs ratio. It is more reliable to invert for Vp and Vp/Vs than Vp and Vs separately when there are fewer high-quality S-arrivals (Eberhart-Phillips, 1986). The procedure starts from determining the calculated arrival times using trial hypocenters, origin times, an initial model of the seismic velocity structure, and the pseudo-bending technique of Um and Thurber (1987) for the ray tracing (Supplementary Figure S2). The residuals between observed and calculated arrival times can be related to perturbations to the velocity structure parameters that are iteratively performed in order to minimize the weighted root-mean-square (RMS) misfit.
During the inversion, there are the grid nodes that were both wellsampled and poorly sampled by the rays. Hence, damping regularization is required to stabilize the solution which is determined by comparing data variance with model variance, as shown in the trade-off curves of Supplementary Figure S3 for a series of one-step inversions within a range of damping values (Eberhart-Phillips, 1986). An appropriate damping value of 50 that balances the data and model variance was selected for both Vp and Vp/Vs The earthquakes are relocated at each iteration step using the current determined 3-D velocity model. We used a weighting scheme as a function of residual with the full weight (1.0) for residuals less than 0.15 s, zero-weight (0.0) for residuals greater than 10 s, and an intermediate value of 0.02 for residuals near 1.5 s. The linear tapers are from 0.15 to 1.5 s and 1.5-10 s, so that the high residual data above 10 s will be automatically discarded in early inversion steps. We set the maximum number of iterations up to 20. However, the inversion automatically stopped at the 15th iteration when the convergence and least weighted RMS were achieved. An F-test in SIMULPS is used to provide a stopping criterion, indicating that variance reduction has become insignificant with further iterations.

Resolution Tests
The capability of the parameterized model and the distribution of data to resolve the structure beneath the study area can be tested with synthetic data constructed using the same distribution of earthquakes and stations as in our data set. We used several resolution tests, i.e., ray hit count (RHC), derivative weighted sum (DWS) (Toomey and Foulger, 1989), diagonal resolution element matrix (DRE) (Menke, 1989), and checkerboard resolution test (CRT), which is the most common way to examine the resolution of seismic tomography (Supplementary Figure S4 and Supplementary Figure S5). For adjacent grid nodes of the model, we perturbed the velocity by ± 10% relative to the initial 1-D velocity model and then calculated the travel times to the stations for each earthquake location. In order to evaluate the resolution test more appropriately, we added Gaussian noise with a mean of −0.0401 and a standard deviation of 0.6149 s to synthetic data (Supplementary Figure S6 and Supplementary Figure S7). This synthetic data set of travel times was then inverted for the velocity structure using the same number of iterations as the real inversion to examine how well the checkerboard pattern was recovered and is usually consistent with the higher values of DWS, RHC, and DRE.
The CRT results highlight the resolved area within the dashed polygons in Figure 2 and Figure 3 for horizontal and vertical sections, respectively, which show a reasonably good resolution for the study area at depths of 10-90 km. Following most studies that use the checkerboard test, if the inversion results show that one grid/block has the same anomaly (positive or negative) as the input model, even with a reduction in amplitude due to the implementation of damping, we consider that block to be resolved by the data while considering the limitations of CRT   (Rawlinson and Spakman, 2016). The similar areas are consistent with relatively high values for DRE, DWS, and RHC. Fairly high resolution is obtained around the slab model 2.0 (Hayes et al., 2018) and the mantle wedge, since many earthquakes are situated in the plate boundary and many seismic rays travel to the stations onshore, providing a relatively high density of seismic rays around given nodes beneath the region. Below 150 km, the resolution becomes smeared due to sparse many deep earthquakes.

Vp, Vs and Vp/Vs Models
After the 15th iteration, the tomographic inversion was automatically stopped to update the velocity model and relocate the seismicity when the convergence and least weighted RMS had been reached (Supplementary Figure S8). It leads to a total weighted RMS reduction from 0.73022 to 0.32071 s with the P-and S-P data variances of 0.1483 and 0.0032 s 2 , respectively. Maps of seismicity distribution before and after relocation are shown in Supplementary Figure S9 and Supplementary Figure S10  The final tomographic images of Vp and Vs (Figure 4 and Figure 5) are shown as percent perturbations relative to the initial 1-D seismic velocity model of Central Java . The blue and red colors represent high and low-velocity changes, respectively. The Vp/Vs values are shown in absolute values (Figure 4 and Figure 5). The Vp/Vs ratio is related to Poisson's  (Christensen, 1996;Jo and Hong, 2013;Lin et al., 2016). For this reason, the Vp/Vs ratio has been widely used to infer structural properties of subducting slabs such as the presence of partial melts or fluid-saturated medium (e.g., Nakajima et al., 2001;Zhang et al., 2004;Eberhart-Phillips et al., 2005;Nugraha and Mori 2006;Reyners et al., 2006;Koulakov et al., 2007;Tsuji et al., 2008;Rohadi et al., 2013;Rosalia et al., 2019). Based on the results of the resolution tests ( Figure 2 and The horizontal sections of the tomographic results are shown from the depth of 10-90 km (Figure 4). The map of Bouguer gravity anomaly in Figure 6A shows correlated features, which is taken from Center for Geological Survey (PSG) (https://geology.esdm.go.id/). In general, we observed that the low velocities (Vp and Vs) dominate around the volcanic complex of the study area which may reflect some features such as sediments or magmatism beneath the volcanoes (Nakajima et al., 2001;Zhao, 2015). These strong lowvelocity anomalies can be observed from the depth of 10-30 km beneath Mt. Lawu, Mt. Wilis, and Malang volcanic complex; and from 20 to 30 km beneath Mt. Merapi-Merbabu complex and Mt. Ijen-Baluran. The estimated Vp/Vs ratio for these regions is also relatively high (>1.75), suggesting the presence of fluids or partial FIGURE 6 | (A) Map of Bouguer gravity anomaly Center for Geological Survey (PSG) overlayed by Vp perturbation at 10 km depth (black contour) and (B) marine free-air gravity anomaly (Sandwell et al., 2014) in Central and East Java region. Black lines indicate faults and trench delineation. Black triangles are the volcanoes. Inset shows the location of subducted seamount.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 756806 8 melts in the sedimentary layer of the Kendeng basin. It correlates with the location of the strong negative Bouguer gravity anomaly of more than −75 mGal, which is at least 400 km long, parallel to the mainland ( Figure 6A) (Smyth et al., 2005Waltham et al., 2008) and previous tomographic studies which also confirmed the low-velocity anomaly area beneath Mt. Merapi-Lawu in Central Java with a high Vp/Vs ratio of 1.9 Wagner et al., 2007;Rohadi et al., 2013;Zulfakriza et al., 2014).
Meanwhile, the high-velocity anomalies (Vp and Vs) are considerably found at a depth of 10 km in the southern and northern part of the island, which are associated with the positive Bouguer gravity anomaly of the Southern Mountains zone and Rembang zone, respectively (Smyth et al., 2005Waltham et al., 2008). It is also consistent with low Vs anomaly in the northern and high-Vs anomaly in southern Java, as indicated by ambient noise tomography (Martha et al., 2017). From 20 to 30 km, we observe the high-velocity anomaly in southern Java that probably corresponds to the rigid forearc, and below 50 km, is associated with the oceanic lithosphere that is subducted toward the north.
We present vertical cross-sections through several areas across the Central and East Java region ( Figure 5) to look at the vertical distribution of velocities with the overlayed seismicity. The locations of the vertical sections are shown in Figure 1. The presence of the subducted slab is recognizable, dipping toward the north at depths of about 50-100 km, by the high-velocity areas with a low Vp/Vs ratio (e.g., Zhang et al., 2004;Reyners et al., 2006) and is consistent with the delineation of slab model 2.0 (Hayes et al., 2018). The earthquakes are also located in this area with high Vp and Vs to a depth of 150 km along with the slab. The region below 150 km depth is not well-resolved in our model. According to the tomograms, seismicity, and slab model, the dip of the subducted slab is increasing with depth toward the north, as suggested by previous studies (Fukao et al., 1992;Koulakov et al., 2007;Zenonos et al., 2019). However, the thickness and geometry of the slab are hard to estimate precisely from tomographic images. Fine structures of the subduction slab can be further studied by a dense seismic network, in particular, the ocean-bottom-seismometer (OBS) network (Zhao, 2015). Meanwhile, low-velocity zones above the top limit of the subducted slab are likely to be found at ∼100 km depth as the possible location of partial melting zones. Thermal models of the Java subduction zone also estimate the transition depths (from the decoupling to coupling between the slab and the overriding mantle) of about 80-155 km with the maximum temperature at the transition zone and the mantle wedge beneath the volcanic arc up to 693 and 1,458°C, respectively (Syracuse et al., 2010). In such areas, the high-fluid and thermal conditions are probably responsible for this anomaly-the water released from marine sediments and the oceanic crust of the subducted slab before the downgoing slab reaches beneath the volcanic arc. Fluids from this slab dehydration and temperature increase in the mantle wedge due to the heat convection from the deeper part of structure, lead to partial melting. Thus, seismic velocity is significantly lower relative to surrounding dry areas (Syracuse et al., 2010;Zhao, 2015;Tonegawa et al., 2017).
As previously mentioned, low-velocity anomalies with high Vp/ Vs ratio are observed beneath the volcanic complexes, especially seen in the vertical sections. In Figure 5A, where the vertical section A-A' slices the Mt. Merapi-Merbabu complex, we see low-velocities (Vp and Vs) at a depth of 20 km. It is probably associated with the magmatism beneath Mt. Merapi or Kendeng sedimentary basin, as also represented by previous local tomographic studies in the same area Wagner et al., 2007;Luehr et al., 2013;Ramdhan et al., 2019). Mt. Merapi is one of the most active volcanoes in the world, with an eruption cycle of two to 6 years and large eruptions occur every 50-100 years (Luehr et al., 2013;Ramdhan et al., 2019). Ramdhan et al. (2019) has imaged the internal structures of Mt. Merapi and determined two anomalies that are interpreted as shallow and intermediate magma bodies at depths of 5 and 15 km, respectively. The presence of these two magma bodies cannot be distinguished by the resolution of our tomography. However, our velocity models may depict the partial melting zone at a depth of ∼100 km, indicated by low-Vp and Vs The ascending fluids and melts migrate upward from the slab and concentrate beneath the volcanic front, supplying magma to the active volcanoes. Koulakov et al. (2007), Wagner et al. (2007), and Luehr et al. (2013) also highlight the prominent feature of Merapi-Lawu anomaly (MLA), as the strong low-velocity anomaly situated close to the Mt. Merapi and extendings to Mt. Lawu in the east. They discuss the origin of the MLA feature that is partly caused by thick lava and sedimentary deposits in the Kendeng basin which correlates with a negative gravity anomaly.
In the central part of the study area, Figure 5B shows vertical section B-B' through Mt. Wilis-Pandan. A low-Vp and Vs anomaly with a high Vp/Vs ratio down to a depth 30 km beneath Mt. Pandan is probably associated with a location of magma, and its movement toward the surface can be related to the earthquake clusters. Mt. Pandan is defined as a dormant volcano north of Mt. Wilis and close to the Kendeng thrust fault. Mt. Pandan is not listed as an active volcano, and it is unclear whether the local earthquakes are associated with the Kendeng thrust activity or magmatic activity. A gravity survey has been conducted to delineate the internal structure of Mt. Pandan and link it to the earthquakes sequence around the region (Santoso et al., 2018). The map of Bouguer anomalies shows lowdensity value beneath Mt. Pandan that may be related to a weak zone, hot materials, or magma source. Santoso et al. (2018) suggested that the subduction process triggers the fault movement and, at the same time, triggers the ascending magma to the surface.
In the eastern part of the study area, the volcanoes are distributed relatively close to each other. Most of the volcanic areas show features of low velocities that are interpreted to be regions of fluids and magma. In Figure 5C, section C-C′ slices through Mt. Semeru-Bromo. Low-velocity anomalies with high Vp/Vs are observed at ∼ 50 km depth beneath the volcanoes. Section D-D' in Figure 5D, which passes through Mt. Ijen-Baluran in the easternmost Java, shows a low-velocity anomaly with a high Vp/Vs ratio from ∼50 km depth beneath the volcano, and the wedge also appears to be associated with a low Vp/Vs ratio. The interaction between the volcanic and sedimentary domains is concentrated at the surface in East Java, producing sediment-hosted geothermal system and clasticdominated geyser-like eruption (nicknamed Lusi-Lumpur Sidoarjo) (Fallahi et al., 2017;Lupi et al., 2022). The plumbing system that feeds Lusi has been proposed by several tomographic studies in the region (e.g., Fallahi et al., 2017;Karyono et al., 2020;Lupi et al., 2022), suggesting that the magmatic system of Lusi connects to the Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 756806 surrounding volcanoes (Mt. Penanggungan and Mt. Arjuno-Welirang complex. The plate subduction promotes the formation of thrust faults below the volcanic arc. Magma and fluids originating from volcanic complex migrate through the weakened rocks along the fault system to the sedimentary basin, manifesting in the sediment-hosted geothermal system of Lusi. Furthermore, our results depict the partial melts feature at ∼ 100 km depth above the slab that is indicated by a low-velocity anomaly with a high Vp/Vs ratio. The partial melting is promoted by joint effects between heat convection in the mantle wedge and fluids from slab dehydration beneath the volcanic arc, producing magmatism in the region (Zhao et al., 2012). Based on the Vp models, we present the interpretative crosssections in Figure 7 for the Central and East Java region. As the oceanic Indo-Australia plate moves northward and subducts beneath the Eurasian plate, partial melting occurs at a depth of about 100 km. Fluids and melts are ascending to the surface and feeding the volcanoes located above the low-velocity anomalies. We realize that our tomographic results might have limitations in producing smooth images. Some features can be observed in certain sections while not seen in other sections. We expect that the main factor is the limitation of the ray coverage. The station distribution used in this study is relatively sparse, and path distribution can also be highly heterogeneous. Model parameterization with an irregular grid may be adapted to improve the stability of inversion and extractions of structural information beneath Central and East Java in further studies, as described by Rawlinson et al. (2014).

The Relationship to Recent Large Earthquakes
Inland faults contribute to the potential seismic hazard in the region since they have been recognized as active structures based on geological and geodetic studies (Marliyani, 2016;Pusat Studi Gempa Nasional, 2017), but some faults do not show a significant rate of recent seismicity. On the other hand, we observed shallow clustered earthquakes near the Opak fault, Central Java region. The activity of Opak fault is considered to be the cause of the 2006 Yogyakarta earthquake (Mw 6.3), with latitude, longitude, and depth of 8.04 o S, 110.36 o E, 15 km, respectively (International Seismological Centre, 2021) (http://www.isc.ac.uk/). The focal mechanism from the Global Centroid Moment Tensor (GCMT) (https://www.globalcmt.org/) (Dziewonski et al., 1981;Ekström et al., 2012) shows a strike-slip faulting mechanism with strike in the NE-SW direction. Figure 8B represents a vertical section of the tomographic model across the fault, along with the clustered earthquakes in the period of our observation, with depths less than 20 km. These earthquakes are located in the contrast of seismic velocities, reflecting the existence of different types of lithology. It is consistent with previous tomographic studies (e.g., Rohadi et al., 2013;Ramdhan et al., 2017;Diambama et al., 2019), and they suggested that low-velocity anomaly is associated with carbonates or other sedimentary rocks that overly the Oligo-Miocene volcanic deposits characterized by the high-velocity. However, we need to conduct a higher resolution of tomography in order to better image small-scale anomalies around the fault. Thi can be achieved by increasing the density of ray path coverage between the local earthquake and seismic stations and grid nodes refinement. At least, we can infer from our results that the geometry of Opak fault is dipping to the east, which supports the previous result of Tsuji et al. (2009) andDiambama et al. (2019).
The recent study of Widiyantoro et al. (2020) has revealed the seismic gaps south of Java, one of which is in the Central and East Java region. The lack of recent megathrust events may indicate great tsunamigenic earthquakes (Mw8.8) can be possible with return periods of about 400 years. The latest tsunami earthquake  respectively (International Seismological Centre, 2021). A thrusting fault mechanism is indicated for the mainshock (Dziewonski et al., 1981;Ekström et al., 2012) (https://www. globalcmt.org/), but the aftershocks are dominated by normal faulting (Abercrombie et al., 2001). They interpreted the observations as a mainshock that slipped over a subducted seamount, which produced normal faulting aftershocks at the outer rise of the Indo-Australian plate. On April 10, 2021, a Mw 6.1 earthquake occurred near Malang, East Java, without triggering a tsunami for the coastal FIGURE 8 | (A) Map of the large earthquake epicenters (yellow stars) and the seismicity (gray circles) used in this study. Focal mechanism solutions are taken from GCMT ((https://www.globalcmt.org/) The rectangles E and F are the area used to plot the vertical cross-section in Figure 8B, Figure 8C. Two bottom figures are vertical sections for Vp, Vs, and Vp/Vs ratio with relocated seismicity across (B) the 2006 Yogyakarta earthquake (Mw 6.3) (E-E9), (C) the 1994 Banyuwangi earthquake (Mw 7.8) and the 2021 Malang earthquake (F-F9). Green lines indicate the slab 2.0 model (Hayes et al., 2018). Black triangles show locations of volcanoes.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 756806 11 area of southern Java. However, the event produced strong shaking with MMI V in the East Java area, according to the shakemap reported by BMKG, causing building damage and fatalities (http://shakemap.bmkg.go.id/). The preliminary report of the hypocenter location is at 8.84 o S, 112.47 o E with a depth of 86 km (http://repogempa.bmkg.go.id/) with a thrust fault mechanism. To better constrain the event location, we have manually re-picked the mainshock arrivals from 19 stations (Supplementary Figure S7 and Supplementary Figure S8) and determined the hypocenter by applying a nonlinear method (Lomax et al., 2000). We plotted the wadati diagram to see the linearity between the phases and obtained a Vp/Vs ratio of 1.75 (Supplementary Figure S9). Finally, we have updated the hypocenter location of the 2021 Malang earthquake (Mw 6.1) at 8.94 o S, 112.45 o E with a depth of 59.7 km and included the result on our tomographic models. The hypocentral errors for the event are 3.08, 6.39, and 11.91 km in x, y and z directions, respectively. Figure 8C shows the vertical section of the velocity model across the 1994 and 2021 events. Both events are located on the high-velocity region associated with the geometry of the oceanic slab that steepens towards the north. Seismic clusters of our data period are also located relatively close to the 2021 event and provide the location of the rupture area in this region. Any ground shaking caused by large magnitude earthquakes will be a threat to the population of Java, although without triggering a tsunami. Meanwhile, the area near the 1994 event lacks thrust events since it is the location of the seismic gap. At ∼ 50 km depth, a low-velocity anomaly exists at the boundary of the subducted slab with a ∼50 km width and extends to the depth of ∼25 km. This feature may exhibit the presence of subducted seamount that were previously acknowledged by Abercrombie et al. (2001). Seamounts off Java have been initially recognized in early side-scan sonar along the eastern Java Trench (Masson et al., 1990) and revisited by multichannel seismic reflection profiles and multibeam bathymetric data (Kopp et al., 2006;Kopp, 2011;Xia et al., 2021). Some of these seamounts are currently being subducted and causing uplift of the forearc to the shallower water depth (Supplementary Figure S16). In the area of 1994 earthquake, a large subducting seamount modulates the seafloor bathymetry (Xia et al., 2021) and is manifested in the isolated free-air gravity anomaly close to the trench (Sandwell et al., 2014) (black square in Figure 6B). Based on these observations, Xia et al. (2021) suggest that the seamount has a height of ∼2 km and a possible width of ∼40-60 km. However, our resolution with grid spacings of 50 and 10 km in lateral and vertical variations, respectively, does not allow us to estimate the true dimension of this subducted seamount. Similar features of subducting seamounts indicated by low-V zones are also observed in the Pacific slab beneath Kanto, Japan with a ∼30 km width (Nakajima and Hasegawa, 2006) and Cocos slab beneath Costa Rica at 30 km depth (Husen et al., 2003;Dinc et al., 2010) and are associated with more hydrated regions compared to the surrounding slab. A possible location of partial melting is observed at a depth of 100 km, indicated by a low-velocity anomaly with a high Vp/Vs ratio. Fluids from slab dehydration and the temperature increase in the mantle wedge promote the partial melting in the area. The fluids/melts then ascend from the mantle wedge to feed volcanoes at the surface.

CONCLUSION
We have produced tomographic images of Vp, Vs, and Vp/Vs ratio beneath the Central and East Java region. Along with the relocated seismicity, we discuss the relationship between the earthquakes, subduction process and volcanism. In general, low-velocity features with relatively high Vp/Vs ratio dominate around the volcanic complexes, which probably correspond to the presence of fluid material and possible partial melting in the upper mantle which ascends from a depth of about 100 km. Shallow low-velocity features with low Vp/Vs ratio are probably associated with sedimentary deposits, such as the Kendeng basin. The subducted slab is revealed by high-velocity features with low Vp/Vs ratio at 50-100 km depths. The dip of the subducted slab is increasing with depth toward the north.
The 2006 (Mw 6.3) Yogyakarta earthquake is located on the structural discontinuity as indicated by the contrast between highlow velocity around the source region. We also have redetermined the hypocenter location of the 2021 (Mw 6.1) Malang earthquake at 8.94 o S, 112.45 o E, with a depth of 59.7 km. The 2021 event, along with the 1994 (Mw 7.8) Banyuwangi earthquake, are located in the high-velocity region and close to the seismic clusters, which reflects the geometry of the oceanic slab. A low-velocity anomaly at ∼ 50 km depth is associated with the presence of a subducting seamount that is likely more hydrated than the surrounding slab.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
FM, AN, JM, and NP conceived the study; FM, AN, JM, and PS contributed to the writing of the manuscript. All authors contributed to the preparation of the manuscript. All authors have read and approved the final manuscript.