Ambient seismic noise has been widely used to study earth structure, using e.g., Ambient Noise Tomography (ANT) and Horizontal-to-Vertical Spectral Ratio (HVSR). Ambient seismic noise used in ANT studies is usually generated by microseisms, resulting from the interaction between ocean waves and the solid earth (Snieder and Wapenaar 2010). Microseisms generate ambient seismic noise at frequencies below 1 Hz, with the primary microseism energy, attributed to topographic coupling between ocean swell and seismic waves on the continental shelf, in the range 0.02–0.1 Hz, while secondary microseisms, attributed to non-linear forcing by ocean swell in both pelagic and coastal regions, is in the range 0.1-1 Hz (Nishida 2017). It has been shown that cross-correlating the ambient noise recorded simultaneously at a pair of seismographs can yield an estimate of the response of the medium (i.e., a point-source Green’s function) that would be measured if there were a source at one of the two seismographs and a receiver at the other (Lobkis and Weaver 2001; Shapiro and Campillo 2004; Snieder and Wapenaar 2010).
Previous studies have shown that ANT can be used to image upper crustal seismic velocity structure in Indonesia, e.g., in Lake Toba (Stankiewicz et al. 2010), Ijen Volcano (Spica et al. 2015), Central Java (Zulfakriza et al. 2014), East Java (Martha et al. 2017), Tarutung pull-apart basin, Sumatra Fault (Ryberg et al. 2016), Banda Arc (Porritt et al. 2016), Jakarta (Saygin et al. 2016), Bandung Basin (Pranata et al. 2020), and Agung-Batur Volcano, Bali (Zulfakriza et al. 2020).
Rosalia et al. (2020) utilized a data set collected from 85 temporary broadband seismographic stations deployed in western Java during 2016–2018 (Fig. 2), to estimate interstation Green’s functions using cross-correlations of ambient seismic noise. They applied a Trans-dimensional Bayesian approach to estimate fundamental-mode Rayleigh wave group velocity maps from these data in the period range 1–25 s and showed how these maps correlate with surface geology and centers of volcanic activity. Here we consider the inversion of these 2D group velocity maps to estimate shear velocity (Vs) structure as a function of depth in the upper crust. We consider a grid of points regularly distributed throughout the study area, at which the Rayleigh wave group velocity maps define dispersion curves over the frequency range 0.04-1.00 Hz. Each Rayleigh wave group velocity dispersion curves is inverted for a Vs profile over 0–20 km depth at the respective grid point, and these Vs depth profiles are interpolated between grid points to obtain a pseudo-3D model for Vs structure in the upper crust.
3.1 Depth Inversion using Neighbourhood Algorithm
We used the Neighbourhood Algorithm (NA) method to obtain Vs depth profiles at each point in a grid regularly distributed over the study area. The NA method is a Monte-Carlo global direct-search method which uses successive sampling of the model space to assemble an ensemble of models that preferentially samples the low-misfit regions of parameter space (Sambridge 1999a, b).
The NA utilizes a Voronoi cell tessellation if the model space in both a search and appraisal stage. The search stage consists of a direct search in multidimensional parameter space for an ensemble of models that preferentially sample its low-misfit region(s). In this stage, the NA first generates n0 models randomly distributed throughout the model-space. A mesh of Voronoi cells, each covering the nearest-neighbour portion of space associated with each model is created. For each model, a theoretical dispersion curve is computed and the misfit between the theoretical and observation curve is assigned to the corresponding cell. Then, ns new models are generated within each cell having lowest misfit. A new Voronoi cell tessellation is computed, new misfits are calculated, and the lowest-misfit cells are chosen to resample. This procedure is repeated, and the algorithm stops after a fixed number of iterations. The method results in fine sampling of portions of the model space associated with low misfit, while high-misfit regions are coarsely sampled. The appraisal stage consists of an algorithm for using the entire ensemble of models produced in the search stage to extract information about model uncertainty.
In this study, we produce a series of Vs depth profiles using the NA method to invert group velocity dispersion curves constructed from group velocity maps obtained using trans-dimensional tree inversion (Rosalia et al. 2020, Hawkins & Sambridge, 2015). We grid the study area with intervals of latitude and longitude of 0.1o, producing 557 grid points (Fig. 3a). In each point grid, a dispersion curve is extracted from the 2-D group velocity maps at different periods. Each dispersion curve is then inverted using the Dinver package of Geopsy (Wathelet 2008) to generate a 1-D shear-wave velocity profile at each corresponding grid point. Combining all 1-D Vs profiles from all grid points results in a 3-D model of the subsurface.
In the NA, layered models are selected by choosing parameters from a specific range bounded by a priori minimum and maximum values. The layer parameters used in the surface wave inversion can be described using: S-wave velocity (Vs), thickness (H), P-wave velocity (Vp), and density (Wathelet et al. 2004; Wathelet 2008). The surface waves are less sensitive to Vp and density than to Vs (Takeuchi and Saito 1972; Bache et al. 1978; Tanimoto 1991; Aki and Richards 2002). Therefore, in this study we will just focus on Vs structure, by constraining Vp and density as descried below.
We parameterize our model based on the CRUST1.0 model (Laske et al. 2013). The CRUST1.0 model is specified on a 1o x 1o grid and incorporates global estimates of sediment thickness, which is used as the guidance to determine the range of depth and velocity parameters. In this inversion, we used three horizontal layers and one half-space layer which is associated with one sediment, two crystalline crustal layers and a half-space, denoted Sediment and Crystalline Crust 1, 2, and 3, respectively, in Tables 1 and 2. Each layer is divided into sublayers that are evenly spaced, with velocities determined by linear interpolation of the velocities at the top and bottom of the layer. The Sediment layer is divided into three sub-layers and the Crystalline Crust 1 and 2 layers are each divided it into five sub-layers. The inversion allows the top and bottom velocities and thickness of each layer to vary.
Table 1
Model Crust 1.0 (Laske et al., 2013) in the western part of Java
Layer | P-wave (km/s) | S-wave (km/s) | Density (g/cm3) | Bottom Depth (km) |
Sediment | 2.0-3.6 | 0.6–1.9 | 1.9–2.3 | 1.39–4.30 |
Crystalline Crust 1 | 5.8-6.0 | 3.4–3.5 | 2.6–2.7 | 8.44–12.15 |
Crystalline Crust 2 | 6.3–6.6 | 3.6–3.8 | 2.7–2.9 | 20.56–21.92 |
Crystalline Crust 3 | 6.9–7.1 | 3.9–3.9 | 2.9–3.1 | 30–32 |
Table 2
The Range of Layer Parameters Used for Vs Inversion in Western Java.
Layer | P-wave (km/s) | S-wave (km/s) | Density (g/cm3) | Bottom Depth (km) |
Sediment | 1.0–4.0 | 0.5–2.5 | 2.0 | 0.5–4.5 |
Crystalline Crust 1 | 3.0–8.0 | 2.0-3.5 | 2.1 | 5.00-11.5 |
Crystalline Crust 2 | 4.0–8.0 | 2.8-4.0 | 2.3 | 12.5–22.00 |
Crystalline Crust 3 | 6.0–8.0 | 3.0–4.0 | 2.6 | Half-space |
We noted the range of values for the layer parameters Vp, Vs, and depth taken on by the Crust 1.0 model in western Java (Table 1), but found these ranges were too narrow to be used as the range of layer parameters in our inversion. We therefore chose the wider bounds on layer parameters indicated in Table 2, in order to allow for the inversion to explore models that were outside the range of Crust 1.0 model for western Java. In the inversion we assumed a fixed, uniform density in each layer of 2000 kg/m3, 2100 kg/m3, 2300 kg/m3, and 2700 kg/m3 for each layer, respectively. In this inversion, the relation of Vp and Vs is bounded by Poisson's ratio, ranging from 0.2 to 0.5, which, together with the prior information of geological structure, may help limit parameter space to be searched.
Other parameters that should be considered are the two tuning parameters that control the behavior of the NA’s exploitation (quick convergence of the misfit function) and exploration (investigation of nearly all local minima to find a global solution) of the model space. The parameters are the number of new models to sample at each iteration (ns) and the number of resampled Voronoi cells (nr). For the inversion parameters, we set the initial number of model samples (ns0) initially before starting the iteration to 200, the ns to 100, and the nr to 100. Based on those parameterizations, we generate 50200 models. In this inversion, we set the maximum number of iterations to 500.