Geophysical Journal International

SummarySubsurface temperature measurements play a crucial role, for instance, in optimizing geothermal power plants and monitoring heat-storage systems. Previous studies have demonstrated that time-lapse variations in temperature can be correlated with variations in seismic wave speeds, offering the potential for temperature monitoring via seismic surveys. However, an apparent discrepancy has emerged between field and laboratory experiments. Field studies predominantly report positive correlations between temperature and seismic wave speeds, while laboratory experiments often show anti-correlations. This inconsistency underscores the need for a more comprehensive, physics-based understanding of temperature-induced wave speed changes. In this study, we strive to bridge the gap between field and laboratory findings by examining several mechanisms governing temperature-induced seismic wave speed changes, namely the intrinsic temperature dependency of elastic parameters and thermally-induced elasticity. We present a physics-based modelling approach to identify the primary mechanisms responsible for temperature-induced seismic wave speed changes. By considering several end-member models, we find that intrinsic temperature dependency of elastic parameters (negative correlation) compete with thermal pressure effects (positive correlation). The precise initial and boundary conditions and physical parameters of the system under consideration will determine the weight of both effects. Temperature related dilatation does not seem to play an important role. We apply our approach to loosely consolidated sediments in the shallow subsurface of the Groningen region, where subsurface temperature fluctuations are driven by seasonal atmospheric temperature fluctuations roughly between −5 and 30 ○C. For these models, we predict seasonal temperature-induced changes in body wave speeds of up to 8% in the first few meters of the subsurface, high-frequency (above 2 Hz) surface wave phase velocity variations in the range of 1-2%, and relative changes in site amplification on the order of 4%. These findings contribute to a more comprehensive understanding of the intricate relationship between temperature and near-surface seismic properties, offering insights for applications as subsurface temperature monitoring systems.

SummaryIn this work, we revisit the seminal concept of Johnson-Koplik-Schwartz (JKS) length Λ, i.e. a characteristic length representing an effective pore size which controls various transport-related properties of porous media, such as, the permeability and the electrical conductivity. We present a novel closed-form equation that predicts the behavior of Λ in partially saturated media, for different saturation states. Using previous models in the literature that predict the intrinsic and relative electrical conductivities under partially saturated conditions, we infer the JKS length Λ and the electrical formation factor F as functions of water saturation and properties associated with the pore-size distribution of the probed porous medium. The proposed method permits to estimate the effective permeability and the relative permeability directly from electrical conductivity measurements, thus opening new-avenues for the remote characterization of partially saturated media. We believe that this new model will prove useful for various characterization and modeling applications from reservoir (CO2 or hydrogen storage) to vadose zone studies.

SummaryThe possibility of a transient rheological response to ice age loading, first discussed in the literature of the 1980s, has received renewed attention. Transient behavior across centennial to millennial time scales has been invoked to reconcile apparently contradictory inferences of steady state (Maxwell) viscosity based on two distinct data sets from Greenland: Holocene sea-level curves and GNSS-derived modern crustal uplift data. To revisit this issue, we first compute depth dependent Fréchet kernels using 1-D Maxwell viscoelastic Earth models and demonstrate that the mantle resolving power of the two Greenland data sets is highly distinct, reflecting the differing spatial scale of the associated surface loading: the sea-level records are sensitive to viscosity structure across the entire upper mantle while uplift rates associated with post-1000 CE fluctuations of the Greenland Ice Sheet have a dominant sensitivity to shallow asthenosphere viscosity. Guided by these results, we present forward models which demonstrate that a moderate low viscosity zone beneath the lithosphere in Maxwell Earth models provides a simple route to simultaneously reconciling both data sets by significantly increasing predictions of present-day uplift rates in Greenland whilst having negligible impact on predictions of Holocene relative sea-level curves from the region. Our analysis does not rule out the possibility of transient deformation, but it suggests that it is not required to simultaneously explain these two data sets. A definitive demonstration of transient behavior requires that one account for the resolving power of the datasets in modeling the GIA process.

SummarySeismic-wave scattering is observed, to variable degrees, on Earth, its moon, and Mars. Current scattering models and data processing typically rely on two end-member phenomena: weak single- or multiple-scattering events (ballistic) on the one hand, or intense scattering such that the wavefield retains no path information or bearing on its origin (diffuse).This study explores the existence of scattering behaviour intermediate between these end-members, as well as the properties of heterogeneous media that facilitate a transition between them. We apply multi-scale entropy and frequency-correlation analysis to seismic coda, and observe a distinct transitional behaviour is present within a part of the investigated parameter space. Analysis of terrestrial, lunar and Martian seismograms further demonstrate the applicability of these new methods across a wide range of scattering behaviours, while also highlighting their shortcomings. Results from the planetary data indicate partially non-diffuse behaviour and low complexity within specific bandwidths of lunar wavefields, potentially contradicting the current paradigm that lunar wavefields are diffuse, and require continued study. Further, Martian seismograms are shown to possess greater statistical entropy than lunar seismograms and diffuse energy properties, yet still display distinct phase arrivals, suggesting substantial scattering and transitional scattering behaviour on Mars. The robust, comparative nature of multi-scale entropy and frequency-correlation analysis, applied to idealised simulation as well as three separate planetary bodies, therefore provides a promising framework for future exploration of scattered wavefields across ballistic, transitional, and diffuse regimes, that complements existing methods.

SummarySeismic interferometry (SI) retrieves new seismic responses, e.g., reflections, between either receivers or sources. When SI is applied to a reflection survey with active sources and receivers at the surface, non-physical (ghost) reflections are retrieved as well. Ghost reflections, retrieved from the correlation of two primary reflections or multiples from two different depth levels, are only sensitive to the properties in the layer that cause them to appear in the result of SI, such as velocity, density, and thickness. We aim to use these ghost reflections for monitoring subsurface changes, to address challenges associated with detecting and isolating changes within the target layer in monitoring. We focus on the feasibility of monitoring pore-pressure changes in the Groningen gas field in the Netherlands using ghost reflections. To achieve this, we utilise numerical modelling to simulate scalar reflection data, deploying sources and receivers at the surface. To build up subsurface models for monitoring purposes, we perform an ultrasonic transmission laboratory experiment to measure S-wave velocities at different pore pressures. Applying SI by auto-correlation to the modelled datasets, we retrieve zero-offset ghost reflections. Using a correlation operator, we determine time differences between a baseline survey and monitoring surveys. To enhance the ability to detect small changes, we propose subsampling the ghost reflections before the correlation operator and using only virtual sources with a complete illumination of receivers. We demonstrate that the retrieved time differences between the ghost reflections exhibit variations corresponding to velocity changes inside the reservoir. This highlights the potential of ghost reflections as valuable indicators for monitoring even small changes. We also investigate the effect of the sources and receivers’ geometry and spacing and the number of virtual sources and receivers in retrieving ghost reflections with high interpretability resolution.

SummaryMaterial density remains poorly constrained in seismic imaging problems, yet knowledge of density would provide important insight into physical material properties for the interpretation of subsurface structures. We test ambient noise wavefield sensitivities to subsurface density contrasts through spatial and temporal wavefield gradients via Wave Equation Inversion (WEI), a form of seismic gradiometry. Synthetic results for 3D acoustic media suggest that it is possible to estimate relative density structure with WEI by using a full acoustic formulation for wave propagation and gradiometry. We show that imposing a constant density assumption on the medium can be detrimental to subsurface velocity images, whereas the full acoustic formulation assuming variable density improves our knowledge of both material properties. It allows us to estimate density as an additional material parameter, as well as to improve phase velocity estimates by accounting for approximations to the density structure. In 3D elastic media, severe approximations in the governing wave physics are necessary in order to invert for density using only an array of receivers on the free surface. It is then not straightforward to isolate the comparatively weak density signal from the influence of phase velocity using gradiometric WEI. However, by using receivers both at the surface and in the shallow subsurface we show that it is possible to estimate density using fully elastic volumetric WEI.

SummaryThe current crustal stress field is of key importance to understand geodynamic processes and to assess stability aspects during subsurface usage. To provide a 3D continuous description of the stress state, linear elastic forward geomechanical-numerical models are used. These models solve the equilibrium of forces between gravitational volume forces and surfaces forces imposed mainly by plate tectonics. The latter are responsible for the horizontal stress anisotropy and impose the inverse problem to estimate horizontal displacement boundary conditions that provide a fit best to horizontal stress magnitude data within the model volume. However, horizontal stress magnitude data have high uncertainties and they are sparse, clustered, and not necessarily representative for a larger rock volume. Even when Bayesian statistics are incorporated and additional stress information such as borehole failure observations or formation integrity test are used to further constrain the solution space, this approach may result in a low accuracy of the model results, i.e. the result is not correct. Here, we present an alternative approach that removes the dependence of the solution space based on stress magnitude data to avoid potential low accuracy. Initially, a solution space that contains all stress states that are physically reasonable is defined. Stress magnitude data and the additional stress information are then used in a Bayesian framework to evaluate which solutions are more likely than others. We first show and validate our approach with a generic truth model and then apply it to a case study of the Molasse foreland basin of the Alps in Southern Germany. The results show that the model's ability to predict a reliable stress state is increasing while the number of likely solutions may also increase, and that outlier of stress magnitude data can be identified. This alternative approach results in a substantial increase in computational speed as we perform most of the calculations analytically.

SummaryWe present a kinematic slip model and a simulation of the ensuing tsunami for the 2020 Mw 7.0 Néon Karlovásion (Samos, Eastern Aegean Sea) earthquake, generated from a joint inversion of high-rate GNSS, strong ground motion and InSAR data. From the inversion, we find that the source time function has a total duration of ∼20 s with three peaks at ∼4, 7.5 and 15 s corresponding to the development of three asperities. Most of the slip occurs at the west of the hypocenter and ends at the northwest down-dip edge. The peak slip is ∼3.3 m, and the inverted rake angles indicate predominantly normal faulting motion. Compared with previous studies, these slip patterns have essentially similar asperity location, rupture dimension and anti-correlation with aftershocks. Consistent with our study, most published papers show the source duration of ∼20 s with three episodes of increased moment releases. For the ensuing tsunami, the eight available gauge records indicate that the tsunami waves last ∼18-30 hours depending on location, and the response period of tsunami is ∼10-35 min. The initial waves in the observed records and synthetic simulations show good agreement, which indirectly validates the performance of the inverted slip model. However, the synthetic waveforms struggle to generate long-duration tsunami behavior in simulations. Our tests suggest that the resolution of the bathymetry may be a potential factor affecting the simulated tsunami duration and amplitude. It should be noted that the maximum wave height in the records may occur after the decay of synthetic wave amplitudes. This implies that the inability to model long-duration tsunamis could result in underestimation in future tsunami hazard assessments.

SummaryAmbient noise tomography is a well-established tomographic imaging technique but the effect that spatially variable noise sources have on the measurements remains challenging to account for. Full waveform ambient noise inversion has emerged recently as a promising solution but is computationally challenging since even distant noise sources can have an influence on the inter-station correlation functions and therefore requires a prohibitively large numerical domain, beyond that of the tomographic region of interest. We investigate a new strategy that allows us to reduce the simulation domain while still being able to account for distant contributions. To allow nearby numerical sources to account for distant true sources, we introduce correlated sources and generate a time-dependent effective source distribution at the boundary of a small region of interest that excites the correlation wavefield of a larger domain. In a series of 2D numerical simulations, we demonstrate that the proposed methodology with correlated sources is able to successfully represent a far-field source that is simultaneously present with nearby sources and the methodology also successfully results in a robustly estimated noise source distribution. Furthermore, we show how beamforming results can be used as prior information regarding the azimuthal variation of the ambient noise sources in helping determine the far-field noise distribution. These experiments provide insight into how to reduce the computational cost needed to perform full waveform ambient noise inversion, which is key to turning it into a viable tomographic technique. In addition, the presented experiments may help reduce source-induced bias in time-dependent monitoring applications.

SummarySeismic and electromagnetic explorations are two of the most successful geophysical applications for understanding the subsurface earth, and the joint interpretation of seismic and electromagnetic survey data can help to better characterize the rocks because they contain independent and complementary information about the rocks. However, the successfulness of the joint interpretation depends on the understanding of the correlations between the elastic and electrical rock properties and their influencing factors. Confining pressure is an important geological parameter that has been found to give rise to linear elastic-electrical correlations in sandstones. However, it is still poorly known about what controls the slopes of the pressure dependent linear correlations, even though slope is one of the most important parameter determining the linear correlation. We make artificial sandstones with controlled porosity and permeability, respectively, and measure their pressure dependent elastic (electrical resistivity) and electrical (P-wave velocity) properties simultaneously as well as porosity. We show that the slopes of the measured electrical resistivity versus P-wave velocity as an implicit function of confining pressure correlate positively with the compliant porosity in all the samples. The results not only reveal the petrophysical parameter that controls the slopes of the pressure dependent linear elastic-electrical correlations in sandstones, but also provide a basis for the discrimination of the slope-controlling parameter from the simultaneously measured elastic and electrical properties.

SummaryEarthquake hypocenters are routinely obtained by a common inversion problem of P- and S-phase arrivals observed on a seismological network. Improving our understanding of the uncertainties associated with the hypocentral parameters is crucial for reliable seismological analysis, understanding of tectonic processes, and assessing seismic hazards. However, current methods often overlook uncertainties in velocity models and variable trade-offs during inversion. Here, we propose to unravel the effects of the main sources of uncertainty in the location process using techniques derived from the framework of global sensitivity analysis. These techniques provide a quantification of the effects of selected variables on the variance of the earthquake location using an iterative model that challenges the inversion scheme. Specifically, we consider the main and combined effects of (1) variable network geometry, (2) the presence of errors in analyst observations and (3) errors in velocity parameters from a 1D velocity model. These multiple sources of uncertainty are described by a dozen of random variables in our model. Using a Monte Carlo sampling approach, we explore the model configurations and analyze the differences between the initial reference location and 100,000 resulting hypocentral locations. The GSA approach using Sobol's variance decomposition allows us to quantify the relative importance of our choice of variables. It highlights the critical importance of the velocity model approximation and provides a new, objective and quantitative insight into understanding the sources of uncertainty in the inversion process.

SummaryThe Gulf of Guinea exhibits a continuous emission of narrow-band and long-period signals (16, 26, and 27 seconds) on teleseismic records, yet the underlying excitation mechanism remains unclear. This study establishes a connection between these tremors and the vibration of thin, decoupled crustal plates at unexplored volcanoes in the gulf. We first formulate the damped plate oscillation equation, by incorporating the vibration of the thin surface crustal plate and magma flow in the subsurface sill. The findings reveal that a fundamental-mode vibration with a period of several dozen seconds can be induced by a crustal plate that is less than 1.0 km thick but extends over tens of kilometers in both length and width, given a subsurface sill depth exceeding 10.0 cm. The thin plate hypothesis also allows for excitation of a few overtone modes, but such waves in higher frequencies diminish over long distances, leaving only the monotonous fundamental-mode vibration at teleseismic stations. The long duration of Guinea tremors at each recurrence is attributed to the presence of low viscosity basaltic magma, which influences the damping factor. Direct wave loads at the shallow gulf serve as the primary vibration source, accounting for seasonal variations and recurring patterns. Sporadic energy bursts may also occur due to large storms. Radiation patterns of Guinea tremors are linked to the geometric structure of the thin plate. Our theoretical estimates of tremor spectra closely align with observed data, confirming the model’s accuracy in capturing reported Guinea tremor characteristics. This study provides valuable insights into the origins of VLP tremors at continental volcanoes.

SummaryThe Global Navigation Satellite System—Acoustic ranging combination technique (GNSS-A) is a seafloor geodetic observation technique that achieves an accuracy of centimeters by combining high-rate GNSS data with acoustic ranging. The technique determines the seafloor position by acoustic ranging between the sea surface and multiple seafloor stations, using GNSS data from the sea surface station. Here, the gradient state of the underwater sound speed structure (SSS) is a significant source of error. The open-source software GARPOS can reduce the effect from underwater gradient structures by retrieving the underwater disturbance as a parameter projected onto the sea surface and seafloor. To evaluate the effects of underwater disturbances, a quantitative comparison of the model parameters is necessary. In this study, we developed a representation method to evaluate features of the ocean field. Here, the expression method was described in the order of a formulation and an interpretation in the case of a two-dimensional cross section and extension to the case of an actual three-dimensional field. This method makes it possible to evaluate SSS states in the GNSS-A observations. As an example, we showed the correlation between the anomaly of the expressed ocean state and the anomaly of the seafloor station position, showing that this expression method is an effective index for correcting bias errors. Additionally, we used the data from sites located in the Kuroshio area, a strong current near Japan, to show that the ocean state can be quantitatively interpreted using this expression method.

SummaryThe purpose of this study is to investigate the seismic ambient noise spectral anomalies that occur near gas reservoirs. These anomalies involve a significant spectral amplification of the vertical component for frequencies generally between 1.5 and 4 Hz and have been reported at various hydrocarbon sites worldwide. There are differing views on the mechanisms responsible for these anomalies. The guideline for this study is that many hydrocarbon reservoirs share a common geological feature: an anticline structure. It appears to cause site effects that influence the amplitude of the ambient noise wavefield. This research examines a dense real dataset of ambient noise recorded at the Chémery underground gas storage site in France. The analysis identifies stable spectral anomalies between 1.2 and 2.4 Hz that are correlated to the position of the anticline structure, which also corresponds to the position of the gas bubble. We use a beamforming technique to study the composition and the origin of the ambient noise, and show that the variations of the spectral anomalies over time are correlated to changes in the source wavefield. Finally, we perform numerical simulations of Rayleigh wave propagation within a realistic 3D velocity model of the Chémery site, while using source distributions directly extracted from real data analysis. The comparison of the simulated anomalies with real data yields a satisfactory qualitative fit. We conclude that the fundamental mode Rayleigh wave site effect on the anticline is the main mechanism of the spectral anomaly.

SummaryA new ray-tracing algorithm based on the shortest path method was designed and optimized to model seismic rays. The workflow is based on Dijkstra’s algorithm to find the shortest path iteratively on Self-Adapting Random Grids. These random grids adapt from one iteration to the next, storing in memory those nodes that minimize the travel time in each iteration, consequently automatically densifying the grid in areas with significant velocity variations. Different data structures and grid geometries were studied, and it was confirmed that using a combination of a binary heap and random grids prevents systematic errors generated by using regular grids while keeping the computational times low. Since the nodes that minimize the travel time are stored in each iteration, the memory requirements increase logarithmically, with the initial iterations having the most significant impact on the error minimization but plateauing afterwards. It was found that, albeit this increase in memory requirements, by focusing the densification of nodes in areas close to the actual trajectory of the rays, the error minimization requires fewer iterations than other implementations based on multiple realizations of random grids (requiring less than a tenth of the number of iteration of other methods). A function to handle extreme topographic variations was also designed, which uses a boolean test to reject rays crossing the topography, making the algorithm suitable for first arrival modelling in complex topography areas, such as the required for tomographic inversion of first arrivals in thrust-belt land seismic.

SummaryPore pressure oscillations induced by stress variations including propagating seismic waves from remote earthquakes have been widely observed in various groundwater systems. The monitored pressure change in wells shows significant water-level oscillations to volumetric strain as well as to S and Love waves. Recent observations demonstrated azimuthal dependence of the pore pressure oscillations with respect to stress indicators and fault zone orientation. Within the fault zone, damaged induced anisotropy is the result of the alignment and orientation of cracks and other internal flaws within the rock. In this work, we provide a complete quantitative description of the pore pressure changes induced by passing seismic waves associated with different orientations and values of principal stress and damage tensor components. The model quantifies the azimuthal dependency of the pore pressure response by a non-dimensional ratio defined as the amplitude of the pressure oscillations induced by a shear strain normalized to the volumetric strain. Three angles and two values are needed to calculate the azimuth dependency of the pore pressure response: the angle between the directions of the maximum horizontal stress and the seismic event; fault zone orientation; micro-crack orientation within the fault zone; damage and stress values. The model predicts that maximum pore pressure response occurs when micro-cracks and maximum horizontal stress are in the same orientation, high damage, and high stress anisotropy. By adjusting these quantities, we recalculate results of recent seismological studies in Arbuckle disposal well, Osage County, Oklahoma. The presented model successfully predicts the observed azimuthal dependence in wave-induced fluid pressure response, and relates the anisotropic response to tectonic indicators such as the orientations of the maximum horizontal stress, fault zone, and micro-fractures.

SummaryIn near-surface surveys, shallow-seismic and ground-penetrating radar (GPR) full-waveform inversions (FWIs) have received increasing attention because of their ability to reconstruct high-resolution subsurface models. However, they have different sensitivities to the same targets and thus may yield conflicting geophysical parameter models. To solve this issue, we have developed an indirect joint petrophysical inversion (JPI) integrating shallow-seismic and multi-offset GPR data. These data are used to reconstruct porosity and saturation whereby we use only strong sensitivities between petrophysical and geophysical parameters. To promote its field application, we proposed an input strategy to avoid measuring rock matrix parameters and make indirect JPI more robust. We apply indirect JPI to the field data acquired in Rheinstetten, Germany and find that it reveals the mechanical, electrical and petrophysical properties more reliably than individual inversions. The reconstructed models are assessed by direct push technology, borehole sample measurements and migrated GPR image. Indirect JPI can fit seismic and GPR observed data simultaneously and provide consistent multi-parameter models, which are hard to achieve by FWIs and individual petrophysical inversions. We also find that the method is robust when there are uncertainties in petrophysical a priori information. Overall, the field example proves the great potential of using indirect JPI to solve real-world problems.

SummaryEquations of motion are derived for (visco)elastic, self-gravitating, variably-rotating planets. The equations are written using a decomposition of the elastic motion that separates the body’s elastic deformation from its net translational and rotational motion as far as possible. This separation is achieved by introducing degrees of freedom that represent the body’s rigid motions; it is made precise by imposing constraints that are physically motivated and that should be practically useful. In essence, a Tisserand frame is introduced exactly into the equations of solid mechanics. The necessary concepts are first introduced in the context of a solid body, motivated by symmetries and conservation laws, and the corresponding equations of motion are derived. Next, it is shown how those ideas and equations of motion can readily be extended to describe a layered fluid–solid body. A possibly new conservation law concerning inviscid fluids is then stated. The equilibria and linearisation of the fluid–solid equations of motion are discussed thereafter, along with new equations for use within normal-mode coupling calculations and other Galerkin methods. Finally, the extension of these ideas to the description of multiple, interacting fluid–solid planets is qualitatively discussed.

SummaryFaults in general, and in clay materials in particular, have complex structures that can be linked to both a polyphased tectonic history and to the anisotropic nature of the intact rock. Drilling through faults in shaly materials allows measuring properties such as the structure, mineralogical composition, stress orientation and physical properties. We combine different petrophysical measurements on core samples retrieved from a borehole drilled perpendicularly to a fault zone affecting Toarcian shales from the Tournemire underground research laboratory (France). The borehole is cross-cutting the entire fault thickness which is of the order of ten meters. We perform several types of measurements: density, porosity, saturation directly in the field, and P-wave velocities together with P-waves anisotropy on core samples taken at regular intervals. Special protocols were developed to preserve as much as possible the saturation state of the samples. From our measurements, we were able to track the increase of damage, characterized by a smooth decrease in elastic moduli from the intact zone to the fault core. We then calculated Thomsen's parameters to quantify the elastic anisotropy evolution across the fault. Our results show strong variations of the elastic anisotropy with the distance to the fault core as well as the occurrence of anisotropy reversal.

SummaryBorneo and Sulawesi are two large islands separated by the Makassar Strait that lie within the complex tectonic setting of central Indonesia. The seismic structure beneath this region is poorly understood due to the limited data availability. In this study, we present Rayleigh wave tomography results that illuminate the underlying crustal structure. Group velocity is retrieved from dispersion analysis of Rayleigh waves extracted from the ambient noise field by cross-correlating long-term recordings from 108 seismic stations over a period of 8 months. We then produce a 3-D shear wave velocity model via a two-stage process in which group velocity maps are computed across a range of periods and then sampled over a dense grid of points to produce pseudo-dispersion curves; these dispersion curves are then separately inverted for 1-D shear wave velocity (Vs), with the resultant models combined and interpolated to form a 3-D model. In this model, we observed up to ± 1.2 km/s lateral Vs heterogeneities as a function of depth. Our models illuminate a strong low shear wave velocity (Vs) anomaly at shallow depth (≤ 14 km) and a strong high Vs anomaly at depths of 20 – 30 km beneath the North Makassar Strait. We inferred the sediment basement and Moho depth from our 3-D Vs model based on iso-velocity constrained by the positive vertical gradient of the Vs models. The broad and deep sedimentary basement at ∼14 ± 2 km depth beneath the North Makassar Strait is floored by a shallow Moho at ∼22 ± 2 km depth, which is the thinnest crust in the study area. To the east of this region, our model reveals a Moho depth of ∼45 ± 2 km beneath Central Sulawesi, the thickest crust in our study area, which suggests crustal thickening since the late Oligocene. Moreover, the presence of high near surface Vs anomalies with only slight changes of velocity with increasing depth in southwest Borneo close to Schwaner Mountain (SM) confirm the existence of a crustal root beneath this region.