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.

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², respectively. Maps of seismicity distribution before and after relocation are shown in Supplementary Figure S9 and Supplementary Figure S10 in the Supplementary Material, showing that the relocated seismicity distribution is more clustered and statistically improved compared to the initial locations. The RMS residuals are mostly below 0.5 s with an average value of 0.216 s (Supplementary Figure S11A). In addition, the travel time residuals associated with the initial and final model are represented as histograms in Supplementary Figure S11B, in which most of the value of travel time residuals are in the range of ±0.75 s.

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 (Koulakov et al., 2007). 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 ratio, a physical property of elastic materials. Distributions of Vp/Vs values indirectly correspond to the variations of pressure, presence of fluid, strength, and temperature (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 Figure 3, Supplementary Figure S4 and Supplementary Figure S5), the regions where the results are considered well constrained are within the dashed polygons.

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 low-velocity 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 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., 2005, 2008; Waltham 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 (Koulakov et al., 2007; 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., 2005, 2008; Waltham 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 (Koulakov et al., 2007; 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 low-density 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 clastic-dominated 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 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 cross-sections 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).

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^oS, 110.36^oE, 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) and Diambama 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 in the region is the 1994 Banyuwangi earthquake (Mw 7.8), with the latitude, longitude, and depth of 10.48^oS, 112.91^oE, 32.1 km, 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 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^oS, 112.47^oE 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 non-linear 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^oS, 112.45^oE 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.