Geophysical Journal International

SummaryMuography is a promising and rapidly developing physical prospecting technique based on the attenuation of muon flux. The feasibility and potential of applying muography to mining were presented in a small number of previous case studies. In this work, we aimed to address three challenges that limit the applicability and efficiency of muography in mineral exploration: (1) application to low-density-contrast ore body exploration, (2) analysis of altitudinal impacts on measurements, and (3) precise and efficient reconstruction. We conducted the first case of applying muography to the exploration for low-density-contrast ore bodies. Six muon imaging systems were placed underground to collect surviving muons for roughly half a year. We analyzed the altitudinal impact on the data measurements and proposed a simplified method to eliminate it. We also developed the seed algorithm, a novel three-dimensional reconstruction algorithm, that can reconstruct anomalies located far away from the detectors and avoid their elongation along the observed directions. Benefiting from the seed algorithm, a low-density-contrast orebody and a limonitic siliceous slate structure were reconstructed, demonstrating the sensitivity of this technique in density distinction; discoveries of several mined-out areas are important for accident avoidance; and reconstruction of the stope and scarps served as strong circumstantial evidence of the reliability of the result. The success of this experiment shows the great value of muography in the economic, research and safety aspects of mineral exploration and inspection. Moreover, the overcoming of challenges is helpful for the development of muography, making it an effective and competitive technique in this field.

SummaryThree-dimensional finite element models, which can handle the stress perturbations caused by subsurface mechanical heterogeneities and fault interactions, have been combined with the finite source inversion to estimate the coseismic slip distribution over the fault plane. However, the mesh grid for discretizing the governing equations in the finite element model significantly affects the numerical accuracy. In this study, we performed kinematic finite source inversion with idealized (regular observation point array; M1A–M1D) and regional (GEONET stations in Japan; M2A–M2H) models with different mesh sizes to quantitatively analyse the effect of the mesh grid size around the fault plane on the inverted fault slip distribution. Synthetic observation data vectors obtained from the finest models (M1A and M2A) are compared with those from the coarser models (M1B–M1D and M2B–M2H), which were adopted to construct Green's function matrix. We found that the coarser mesh models derived a smaller surface displacement, leading to a decrease in the norm of Green's function matrix, which in turn influences the fault slip magnitude from the finite source inversion. Finally, we performed the source inversion for the fault slip distribution of the 2011 Mw 9.0 Tohoku-Oki earthquake using the coseismic surface displacements recorded at the GEONET and seafloor stations and finite element modeling. By reducing the mesh size on the fault, we confirmed that the estimated magnitude of fault slip converged to approximately 80 m, which is consistent with the range of fault slip amounts from previous studies based on the Okada model. At least 0.88 million total domain elements and a 6.7 km2 mesh size on the fault plane with an area of 240 × 720 km2 are required for the convergence of the fault slip. Furthermore, we found that the location of the maximum fault slip is less sensitive to the mesh size, implying that source inversion based on a coarse mesh model (i.e. less than 0.5 million elements and > ∼60 km2 mesh size) can quickly provide the rough fault slip distribution.

SummaryA similar waveform search plays a crucial role in seismology for detecting seismic events, such as small earthquakes and low-frequency events. However, the high computational costs associated with waveform cross-correlation calculations represent bottlenecks during the analysis of long, continuous records obtained from numerous stations. In this study, we developed a deep-learning network to obtain 64-bit hash codes containing information on seismic waveforms. Using this network, we performed a similar waveform search for ∼35 million moving windows developed for the 30 min waveforms recorded continuously at 10 MHz sampling rates using 16 acoustic emission transducers during a laboratory hydraulic fracturing experiment. The sampling points of each channel corresponded to those of the 5.8-year records obtained from typical seismic observations at 100 Hz sampling rates. Of the 35 million windows, we searched for windows with small average Hamming distances among the hash codes of 16 channel waveforms against template hash codes of 6057 events that were catalogued using conventional autoprocessing techniques. The calculation of average Hamming distances is 1430–1530 times faster than that of the corresponding network correlation. This hashing-based template matching enabled the detection of 23,462 additional events. We also demonstrated the feasibility of the hashing-based autocorrelation analysis, where similar event pairs were extracted without templates, by calculating the average Hamming distances for all possible pairs of the ∼35 million windows. This calculation required only 15.5 h under 120 thread parallelisation. This deep hashing approach significantly reduced the required memory compared with locality-sensitive hashing approaches based on random permutations, enabling similar waveform searching on a large-scale dataset.

SummaryIrazú and Turrialba are a twin volcanic complex that marks a distinct stop in volcanism along the Central America volcanic arc. We present a new travel-time velocity model of the crust beneath Irazú and Turrialba volcanoes, Costa Rica, and interpret it considering the results of previous ambient noise tomographic inversions. Data were acquired by a temporary seismic network during a period of low activity of the Irazú-Turrialba volcanic complex in 2018-2019. Beneath the Irazú volcano, we observe low P-wave velocities (VP = 5 km s-1) and low velocity ratios (VP/VS = 1.6). In contrast, below the Turrialba volcano, we observe a low S-wave velocities (VS = 3 km s-1) and a high VP/VS (= 1.85) anomaly. We found that locations of low VP and VS anomalies (- 15 per cent) correspond well with shear wave velocity anomalies retrieved from ambient noise tomography. At shallower depths, we observe high VP and VS anomalies (+ 15 per cent) located between the summits of the volcanoes. Sub-vertical velocity anomalies are also observed at greater depths, with high VP and VS anomalies appearing at the lower limits of our models. We propose a complex structure of an intermediate magmatic reservoir, presenting multi-phase fluid states of a liquid-to-gas transition beneath Irazú and a juvenile store of magmatic fluid beneath Turrialba, while shallow fluid transport provides evidence of magmatic-hydrothermal interactions.

SummaryIncreasing magnetisation within mature hydrocarbon reservoirs provides a new technique in identifying oil-water contacts (OWCs) in cored wells with the potential to assess yield thereby reducing the need for further exploration. Authigenic precipitation of magnetic minerals at OWCs may also help locate paleocontacts (PCs), where structural changes to the petroleum system have caused hydrocarbon remigration. This study determines the magnetic characteristics of magnetic enhancements at OWCs and possibly PCs in silliclastic and carbonate reservoirs at the Wytch Farm oil field, Wessex Basin, UK. Increases in saturation magnetisation and susceptibility are observed at the OWC in 11 of the 12 analysed cored reservoirs owing to the increased presence of magnetite and vivianite. Geochemical analysis and shallow reservoirs suggest biogenic and inorganic mineral precipitation is extensive at the OWC depending on iron, sulphur, and phosphorus availability. Similar magnetic characteristics have been observed in magnetic enhancements above the OWC in numerous wells which may represent OWCs before a basin-wide easterly tilt caused hydrocarbon remigration in the Cenozoic. Multiple magnetic enhancements above the OWC in westerly onshore wells, suggest this remigration may have occurred as numerous phases.

SummaryThe drift history of the Lhasa terrane plays an essential role in understanding the tectonic evolution of the Bangong-Nujiang Tethyan Ocean and the Neo-Tethyan Ocean, as well as the evolutionary history of the Tibetan Plateau. Here, a combined rock magnetic, petrographic, and palaeomagnetic study is performed on the Upper Permian–Lower Triassic limestones (∼259–251 Ma) in the western Lhasa terrane. The site-mean direction for the 28 sites is Dg = 32.1°, Ig = 50.3°, kg = 47.9, and α95 = 4.0° in situ and Ds = 342.9°, Is = 32.7°, ks = 43.2, and α95 = 4.2° after tilt-correction, yielding a palaeopole at 68.9° N, 314.4° E with A95 = 4.3°, corresponding to a palaeolatitude of 18.0° ± 4.3° N. The fold tests are not significant because the sampling section shows monoclinic features with minor variations in their bedding attitudes. The palaeopoles for the directions before and after tilt-correction are compared with reliable Late Permian–Palaeogene palaeopoles obtained from the Lhasa terrane. Based on these comparisons, the studied limestones were remagnetized prior to tilting and this remagnetization most likely occurred during the Early Cretaceous. The depositional environment of the limestones may have changed from anoxic to suboxic and oxic during the Early Cretaceous, leading to the oxidation of iron sulphide to authigenic magnetite. Meanwhile, the Late Jurassic–Early Cretaceous convergence between the western Lhasa and Qiangtang terranes may have resulted in tectonic fluid migration and the formation of calcite veins and stylolites in the limestones. This is supported by the presence of small calcite veins and stylolites in some samples, as well as the fact that the framboidal oxides were formerly sulphides (mostly pyrite), implying that the majority of the iron oxides observed in the limestones were authigenic. These processes indicate that chemical remanent magnetization caused by the growth of magnetic minerals related to tectonic fluid migration was most likely the mechanism for the limestone remagnetization.

SummaryIn the immediate vicinity of a source there are strong gradients in the seismic wavefield that are tamed and modified in DAS recording due to combined effects of gauge-length averaging and local stacking on the local strain field. Close to a source broadside propagation effects are significant, and produce a characteristic impact on the local DAS channels. In the presence of topography, of surface or cable, additional effects are introduced that modify the expected signal. All these influences mean that the results of tap tests used to calibrate the channel positions along a DAS cable may give a distorted view of the actual geometry. Such effects can be important for detailed mapping of faulting processes and comparable features.

SummaryAccurate characterization and monitoring strategies are essential for designing and implementing remedial programs for sites polluted with dense non-aqueous phase liquids (DNAPLs). Electrical resistivity tomography (ERT) is a widely used geophysical technique for mapping subsurface features and processes of interest and exhibits desirable characteristics for DNAPL sites due to its ability to gather large volumes of continuous subsurface information in a non-invasive, cost-effective, and time-efficient manner. However, ERT from the surface suffers from poor imaging quality with depth. Enhanced ERT imaging can be obtained via electrodes deployed on the surface and within horizontal boreholes, but so far it has only been investigated for two-dimensional (2D) imaging. This study evaluates the potential of three-dimensional (3D) surface-to-horizontal borehole (S2HB) ERT configurations for imaging 3D DNAPL source zones. Laboratory tank experiments were first conducted with a 3D S2HB ERT configuration, which consisted of a surface grid and a single borehole line of electrodes, being used to monitor DNAPL migration within porous media. Results demonstrate that 3D S2HB ERT with a single borehole provides improved sensitivity at depth, and therefore enhanced imaging compared to conventional 3D surface ERT. Further tank experiments were performed to assess the performance of single borehole S2HB ERT when (i) the distance between surface and borehole is increased, and (ii) additional horizontal boreholes are included. The S2HB ERT with a single borehole significantly outperforms surface ERT at larger depths and performs comparably to S2HB ERT with multiple boreholes. This study suggests that 3D S2HB ERT with a single borehole can provide the enhanced imaging ability needed to map DNAPLs, while also being relatively practical for implementation at field sites.

SummaryAccurate seismic phase detection and onset picking are fundamental to seismological studies. Supervised deep-learning phase pickers have shown promise with excellent performance on land seismic data. Although it may be acceptable to apply them to OBS (Ocean Bottom Seismometer) data that are indispensable for studying ocean regions, they suffer from a significant performance drop. In this study, we develop a generalised transfer-learned OBS phase picker – OBSTransformer, based on automated labelling and transfer learning. First, we compile a comprehensive dataset of catalogued earthquakes recorded by 423 OBSs from 11 temporary deployments worldwide. Through automated processes, we label the P and S phases of these earthquakes by analysing the consistency of at least three arrivals from four widely used machine learning pickers (EQTransformer, PhaseNet, Generalized Phase Detection, and PickNet), as well as the AIC (Akaike Information Criterion) picker. This results in an inclusive OBS dataset containing ∼36,000 earthquake samples. Subsequently, we employ this dataset for transfer learning and utilize a well-trained land machine learning model – EQTransformer as our base model. Moreover, we extract 25,000 OBS noise samples from the same OBS networks using the Kurtosis method, which are then used for model training alongside the labelled earthquake samples. Using three groups of test datasets at sub-global, regional, and local scales, we demonstrate that OBSTransformer outperforms EQTransformer. Particularly, the P and S recall rates at large distances (>200 km) are increased by 68% and 76%, respectively. Our extensive tests and comparisons demonstrate that OBSTransformer is less dependent on the detection/picking thresholds and is more robust to noise levels.

SummaryA robust, in-situ estimate of shear-wave velocity VS and the small-strain damping ratio DS (or equivalently, the quality factor QS) is crucial for the design of buildings and geotechnical systems subjected to vibrations or earthquake ground shaking. A promising technique for simultaneously obtaining both VS and DS relies on the Multichannel Analysis of Surface Waves (MASW) method. MASW can be used to extract the Rayleigh wave phase velocity and phase attenuation data from active-source seismic traces recorded along linear arrays. Then, these data can be inverted to obtain VS and DS profiles. This paper introduces two novel methodologies for extracting the phase velocity and attenuation data. These new approaches are based on an extension of the beamforming technique which can be combined with a modal filter to isolate different Rayleigh propagation modes. Thus, the techniques return reliable phase velocity and attenuation estimates even in the presence of a multi-mode wavefield, which is typical of complex stratigraphic conditions. The reliability and effectiveness of the proposed approaches are assessed on a suite of synthetic wavefields and on experimental data collected at the Garner Valley Downhole Array and Mirandola sites. The results reveal that, under proper modeling of wavefield conditions, accurate estimates of Rayleigh wave phase velocity and attenuation can be extracted from active-source MASW wavefields over a broad frequency range. Eventually, the estimation of soil mechanical parameters also requires a robust inversion procedure to map the experimental Rayleigh wave parameters into soil models describing VS and DS with depth. The simultaneous inversion of phase velocity and attenuation data is discussed in detail in the companion paper.

SummaryThis paper deals with in-situ characterization of the small-strain shear-wave velocity VS and damping ratio DS from an advanced interpretation of Multi-channel Analysis of Surface Waves (MASW) surveys. A new approach based on extracting Rayleigh wave data using the CFDBFa method has been discussed in the companion paper. This paper focuses on mapping the experimental Rayleigh wave phase velocity and attenuation into profiles of VS and DS versus depth, which is achieved through a joint inversion procedure. The joint inversion of phase velocity and attenuation data utilizes a newly-developed Monte Carlo global search algorithm, which implements a smart sampling procedure. This scheme exploits the scaling properties of the solution of the Rayleigh eigenvalue problem to modify the trial earth models and improve the matching with the experimental data. Thus, a reliable result can be achieved with a limited number of trial ground models. The proposed algorithm is applied to the inversion of synthetic data and of experimental data collected at the Garner Valley Downhole Array site, as described in the companion paper. In general, inverted soil models exhibit well-defined VS profiles, whereas DS profiles are affected by larger uncertainties. Greater uncertainty in the inverted DS profiles is a direct result of higher variability in the experimental attenuation data, the limited wavelength range at which reliable values of attenuation parameters can be retrieved, and the sensitivity of attenuation data to both DS and VS. Nonetheless, the resulting inverted earth models agree well with alternative in-situ estimates and geological data. The results stress the feasibility of retrieving both stiffness and attenuation parameters from active-source MASW testing and the effectiveness of extracting in-situ damping ratio estimates from surface wave data.

SummarySeismic attributes are derived measures from seismic data that help characterize subsurface geological features and enhance the interpretation of subsurface structures: we propose to exploit the hidden layers of Long-Short Time Memory neural network predictions as possible new reflection seismic attributes. The idea is based on the inference process of a neural network, which in its hidden layers stores information related to different features embedded in the input data and which usually are not considered. Neural network applications typically ignore such intermediate steps because the main interest lies in the final output, which is considered as the exclusive exploitable feature of the process. On the contrary, here we analyse the possibility to exploit the intermediate prediction steps, hereafter referred as “Deep Attributes” because they are produced by a deep learning algorithm, to highlight features and emphasize characteristics embedded in the data but neither recognizable by traditional interpretation, nor evident with classical attributes or multi-attribute approaches. Nowadays, classical signal attributes are numerous and used for different purposes; we here propose an original strategy to calculate attributes previously never exploited, which are potentially complementary or a good alternative to the classical ones.We tested the proposed procedure on synthetic and field 2-D and 3-D reflection seismic data sets to test and demonstrate the stability, affordability and versatility of the entire approach. Furthermore, we evaluated the performance of deep attributes on a 4-D seismic dataset to assess the applicability and effectiveness for time-monitoring purposes and comparing them with the sweetness attribute.

SummaryThis article presents a new methodology to assist geophysicists in determining the first break event in a 3D seismic dataset using the well-known technique called Dynamic Time Warping algorithm (DTW), which is usually used to find the optimal alignment between two-time series. We used the optimal path from the cost matrix to identify the first-break in the seismogram using a few picks (seeds) made by an interpreter as a reference to perform this task. Furthermore, the data was pre-conditioned by the topographic and linear moveout to improve the method’s accuracy. To demonstrate the technique’s robustness, first, we applied the methodology in a synthetic seismic data. After demonstrating the efficiency of the algorithm, we applied the aforementioned methodology in the Polo-Miranga 3D seismic cube located in the Recóncavo sedimentary basin, Bahia-Brazil, and in the seismic data acquired from the Blackfoot field in Alberta, Canada. The high-quality results showed consistency in determining the first-break in all ranges of offsets, demonstrating an alternative way to accelerate this seismic processing step. Furthermore, we compared the results obtained by the proposed methodology with an algorithm based on comparing the short-time averages (STA) with long-time averages (LTA). Finally, we performed the static correction calculation to ensure that the time distortion resulting from the terrain and the low-velocity layer (LVL) was mitigated in shoot gathers and in the stacked section.

SummaryThis study presents the findings of a splitting analysis conducted on core-refracted teleseismic shear waves (SKS, SKKS, and PKS, called together as XKS) and local shear waves, obtained from a dense seismological network spanning the Kamchatka Peninsula. The objective of the study is to examine the pattern of mantle flow beneath the study area through the investigation of seismic anisotropy. The peninsula is situated at the northeastern end of the Kuril-Kamchatka subduction zone, where the Kuril trench intersects with the western boundary of the Aleutian trench. The dataset utilized in this study comprises waveform data from a dense network of seismic stations (99 broadband and short period stations for the local shear wave splitting analysis and 69 broadband stations for the SKS splitting analysis). The seismograms were downloaded from publicly available data repositories including the IRIS Data Management Center and the GFZ Data Services (GEOFON program). The dense station coverage allows us to investigate the lateral variations in anisotropy, providing insights into the flow patterns within the mantle. The processing of the combined data sets of local shear wave and teleseismic XKS waves allowed us to partially decipher the source of anisotropy in the mantle. Small delay (splitting) times (∼0.35 sec) observed from the local-S data suggest that anisotropy in the mantle wedge is relatively weak with lateral variations. Larger splitting times (∼1.1 sec) observed for the XKS waves relative to local S suggest that the main part of splitting on the XKS waves occurs in the subslab mantle. On the other hand, the rotational pattern of seismic anisotropy observed by both the local S and XKS waves suggests the presence of a toroidal flow at the NE edge of the subducting slab, which affects both the mantle wedge and subslab mantle. For the regions away from the edge of the slab, the mantle flow seems to be governed mainly by the drag of the lithospheric plate over the underlying asthenosphere.

SummaryPlate-coupling estimates and previous seismicity indicate that portions of the Makran megathrust of southern Pakistan and Iran are partially coupled and have the potential to produce future magnitude 7 + earthquakes. However, the GPS observations needed to constrain coupling models are sparse and lead to an incomplete understanding of regional earthquake and tsunami hazard. In this study, we assess GPS velocities for plate-coupling of the Makran subduction zone with specific attention to model resolution and the accretionary prism rheology. We use finite element model-derived Green's functions to invert for the interseismic slip deficit under both elastic and viscoelastic Earth assumptions. We use the model resolution matrix to characterize plate-coupling scenarios that are consistent with the limited spatial resolution afforded by GPS observations. We then forward model the corresponding tsunami responses at major coastal cities within the western Indian Ocean basin. Our plate-coupling results show potential segmentation of the megathrust with varying coupling from west to east, but do not rule out a scenario where the entire length of the megathrust could rupture in a single earthquake. The full subduction zone rupture scenarios suggest the Makran may be able to produce earthquakes up to Mw 9.2. The corresponding tsunami model from the largest earthquake event (Mw 9.2) estimates maximum wave heights reaching 2 to 5 meters at major port cities in the Northern Arabian Sea region. Cities on the west coast of India are less affected (1-2 m). Coastlines bounding eastern Africa, and the Strait of Hormuz, are the least affected (< 1 m).

SummaryUsing a Bayesian approach we compare anecdotal tsunami runup observations from the 29 December 1820 Flores Sea earthquake with close to 200,000 tsunami simulations to determine the most probable earthquake parameters causing the tsunami. Using a dual hypothesis of the source earthquake either originating from the Flores Thrust or the Walanae/Selayar Fault, we found that neither source perfectly matches the observational data, particularly while satisfying seismic constraints of the region. Instead both posteriors have shifted to the edge of the prior indicating that the actual earthquake may have run along both faults.

AbstractWe introduce a novel probabilistic framework for the solution of non-linear geophysical inverse problems in complex variables. By using complex probability distributions, this approach can simultaneously account for individual errors of real and imaginary data parts, independently regularize real and imaginary parts of the complex model, and still take into account cross-sensitivities resulting from a complex forward calculation. The inverse problem is solved by means of optimization. An application of the framework to complex resistivity (CR) imaging demonstrates its advantages over the established inversion approach for CR measurements. We show that CR data, with real and imaginary parts being subject to different errors, can be fitted adequately, accounting for the individual errors and applying independent regularization to the real and imaginary part of the subsurface conductivity. The probabilistic framework itself serves as a basis for the future application of global sampling approaches, such as Markov chain Monte Carlo methods.

SummaryThis study acquires the co-seismic deformation field of the MW7.4 earthquake that occurred in Maduo, China, on May 22, 2021, based on the BeiDou Navigation Satellite System (BDS), and the comparison with the results obtained by the Global Positioning System (GPS) reveals that the two systems are certain differences in their ability to acquire the horizontal co-seismic deformation field. The maximum difference in the horizontal co-seismic deformation is < 5 mm, and the maximum difference in the vertical co-seismic deformation is 8.7 mm. The dynamic displacement waveforms of the 2021 MW 7.4 Maduo earthquake acquired by BDS and GPS are very similar, which confirms that BDS can acquire ground-shaking images with an accuracy comparable to that of GPS. Based on the empirical relationship equation of the peak ground displacement (PGD) and moment magnitude (MW), this study verifies and calculates both the MW of the 2021 MW 7.4 Maduo earthquake and the error and finds that the MW can be quickly and accurately obtained by using the empirical PGD and MW equations, and this MW value can be used as a supplementary means of calibrating the MW of the large earthquake early warning (EEW) systems, which can be quickly determined by seismic wave data. Finally, by comparing the slip distributions inverted from the BDS and GPS co-seismic deformation fields, this study finds that BDS is equally effective as GPS.

SummaryHyuga-nada, off the Pacific coast of Kyushu along the Nankai Trough in southwest Japan, is one of the most active slow earthquake regions around Japan. We estimated the energies of shallow tremors and moments of shallow very low frequency earthquakes (VLFEs) in Hyuga-nada using data from a permanent onshore broadband network and temporary ocean bottom seismometer observations. The energies and moments of these slow earthquakes have a similar along-strike variation and are generally higher south of the subducted Kyushu-Palau Ridge than near the top of the ridge. This spatial variation is also related to the characteristics of slow earthquake migration. The along-strike migration speed was faster at initiation in the south, where the moments of slow earthquakes are higher. After migration entered the subducted Kyushu-Palau Ridge, its speed was decelerated with a parabolic pattern and their moments became smaller. Assuming a constant patch size of slow earthquakes, we estimated that the stress drop of VLFEs in the south of the subducted ridge was approximately three times higher than that near the top of the subducted ridge. According to our observations and a physical model, this stress drop difference between adjacent regions may cause parabolic migration. We also estimated the scaled energy of slow earthquakes from the ratio of the seismic energy rates of tremors to the seismic moment rates of accompanying VLFEs. The spatial variation in scaled energy is not identified inside the Hyuga-nada. Since the range of scaled energy is similar between the south and near the top of the subducted ridge, the apparent stress may be similar if the rigidity is the same. The dominant range of scaled energy of slow earthquakes in Hyuga-nada is 10−11.5–10−8.5. In addition to having similar or one order smaller values compared to other slow earthquake regions, the range of scaled energy in Hyuga-nada is broader. This broader range suggests wide range of characteristic time and various spectral features of slow earthquakes in Hyuga-nada. Based on a Brownian slow earthquake model, the wide range of characteristic time in this area suggests width variations of slow earthquake source area.

SummaryIn seismic tomography practices the Earth's surface is sometimes assumed to be either spherical or flat for convenience in forward modeling calculations. The effect of irregular surface topography on seismic wave propagation is thus ignored, resulting in biases in the phases and amplitudes of synthetic seismograms, which contribute to the residuals that are mapped into velocity structures in tomography inversions. In this study, we conduct a series of inversion experiments based on the adjoint waveform tomography method to quantitatively assess the topography effect on waveform-based inversion results. We first employ models with simplified topography to better highlight and quantify the topography effect. Results show that when topography effect is ignored in the forward modeling, it is mapped into velocity perturbations, leading to spurious velocity anomalies in tomography models. The strength of the spurious velocity anomaly is quasi-linearly related to locally averaged topography gradient. Our inversion experiments demonstrate that in places of strong topography variation, such as the Longmenshan Fault Zone region where the 2008 Mw7.9 Wenchuan earthquake occurred, topography effect can lead to spurious relative velocity anomalies of up to 10 per cent, which cannot be ignored in waveform-based tomography inversions.