Surface Wave Tomography Research Articles

Ambient noise surface-wave imaging in a hardrock environment: implications for mineral exploration

Summary The advancement of seismic methods is vital for mineral exploration in the ongoing energy transition. In this study, we investigate the application of ambient noise seismic interferometry and surface-wave analysis to characterize the subsurface in a mineral exploration context. We then confirm the results of the passive seismic investigation through an active source experiment. We collected ambient noise data using a 2D seismic line initially deployed for an active source reflection seismic study. By cross-correlating the signals, we retrieved the surface waves and constructed a 2D shear-wave velocity profile using conventional surface-wave analysis. We utilized the active source data to establish initial assumptions about the surveyed medium and then validated the passive seismic experiment. The passive seismic results are concordant with the active source results and allow for the interpretation of geological contacts and fault zones. Our work demonstrates the potential of passive seismic methods for investigating local tectonic settings and their role in hardrock mineral exploration.

Using contrasts in horizontal P-wave reflectivity to map the base of the continental lithosphere: Results for the central and eastern U.S.

Vertical-incidence seismic reflection profiles generated from global phases for 16 earthquakes recorded by stations of the Transportable Array (TA) show distinctive patterns of P-wave reflectivity in the uppermost mantle beneath the central and eastern United States. The overall distribution of reflections identified objectively using the sign test statistic applied to bootstrapped stacks is consistent with a westward increase in depth of the lithosphere-asthenosphere boundary (LAB) from roughly 110 to 250 km that is marked within the lower lithosphere by piecewise continuous segments of elevated horizontal P-wave reflectivity. For some profiles, the onsets of zones of increased reflectivity closely match depths corresponding to the maximum negative S-wave velocity gradients found by surface-wave tomography, suggesting that P-wave reflectivity can be used to help characterize properties of the lower lithosphere. We suggest that the vertical change in horizontal reflectivity straddles the lithosphere-asthenosphere transition, encompassing a broad zone of layering caused by increased strain in the lower lithosphere as well as drag-induced flow in the asthenosphere. Some of the lines also show waveforms that fall within the depth range of arrivals identified as midlithospheric discontinuities (MLDs) in overlapping Sp receiver-function profiles. The reflection waveforms observed in the TA lines are mostly multicyclic with a mix of polarities indicating a layered transition, consistent with previous observations and model studies that show the breakup of single Sp waveforms into a series of less prominent, shorter-period P-wave reflections as the dominant frequency of incident energy is increased.

Surface Wave Imaging Method Based on Spatial Autocorrelation Principle for Mining Microseismic Monitoring

Surface Wave Imaging Method Based on Spatial Autocorrelation Principle for Mining Microseismic Monitoring

Frequency-Bessel multi-mode surface wave tomography utilizing microseisms excited by energetic typhoons

Frequency-Bessel multi-mode surface wave tomography utilizing microseisms excited by energetic typhoons

Combined application of FTAN and cross-spectra analysis to ambient noise recorded by a microseismic monitoring network

A case study of seismic interferometry applied to a small microseismic monitoring network is here presented. The main objectives of this study are (i) to quantify the lateral variability of shear-wave velocities in the studied area, and (ii) to investigate the bias produced by noise directionality and non-stationarity in the velocity estimate. Despite the limited number of stations and the short-period character of the seismic sensors, the empirical Green’s functions were retrieved for all station pairs using two years of passive data. Both group and phase velocities were derived, the former using the widespread frequency-time analysis, the latter through the analysis of the real part of the cross-spectra. The main advantage of combining these two methods is a more accurate identification of higher modes, resulting in a reduction of ambiguity during picking and data interpretation. Surface wave tomography was run to obtain the spatial distribution of group and phase velocities for the same wavelengths. The low standard deviation of the results suggests that the sparse character of the network does not limit the applicability of the method, for this specific case. The obtained maps highlight the presence of a lower velocity area that extends from the centre of the network towards southeast. Group and phase velocity dispersion curves have been jointly inverted to retrieve as many shear-wave velocity profiles as selected station pairs. While the average model can be used for a more accurate location of the local natural seismicity, the associated standard deviations give us an indication of the lateral heterogeneity of seismic velocities as a function of depth. Finally, the same velocity analysis was repeated for different time windows in order to quantify the error associated to variations in the noise field. Errors as large as 4% have been found, related to the unfavorable orientation of the receiver pairs with respect to strongly directional noise sources, and to the very short time widows. It was shown that using a one-year time window these errors are reduced to 0.3%.

2D Sedimentary structures at the southeast margin of the Tarim Basin, China, constrained by love wave ambient noise tomography

Summary Imaging the detailed structure of sedimentary basins is crucial for natural resource exploration and essential for better analysis and correction of sediment responses when studying deeper interior earth structures using seismic data. In the Tarim Basin, previous studies on the sedimentary structures are mostly obtained by active-source seismic surveys, which can provide high-resolution underground interface information but are highly costly or environmentally unfriendly. In this paper, ambient noise tomography, an efficient and economical method based on background vibration, is employed to construct the sedimentary velocity structure at the southeast margin of the Tarim Basin. Based on ambient noise data collected from a linear dense short-period seismic array, we extract Love wave signals from T-T component cross-correlation functions (CCFs) and measure Love wave dispersion curves at a period band of ∼0.3–11.5 s. Then, we utilize a one-step direct surface wave tomography method to image a fine 2D sedimentary shear wave velocity structure with a depth reaching 10 km. Our results reveal a clear layered sedimentary structure, the paleo Tadong uplift, and the thrust Cherchen fault. Our study provides reference sedimentary velocity models for the southeastern Tarim Basin, focusing on depths shallower than 10 km. This model is intended to serve as crucial input for studies requiring detailed shallow sedimentary velocity data. Moreover, our research demonstrates that the application of the ambient noise tomography method with dense arrays has great potential for enhancing resource exploration efforts in sedimentary basins.

Shallow crustal structure of the northern Longmen Shan fault zone revealed by a dense seismic array with ambient noise analysis

The Longmen Shan (LMS) fault zone located in the southwest of China is not only the boundary tectonic zone of the Tibetan Plateau and South China block, but also an important part of the north–south seismic belt of China. To investigate the detailed crustal structure of the LMS fault zone, we deployed a dense seismic array of 151 short-period temporary stations in its northern section from March 16 to April 6, 2023. Continuous vertical component ambient noise waveforms recorded by these stations were cross-correlated for all station pairs, enabling the extraction of Rayleigh wave phase velocity dispersion curves across periods ranging from 0.5 s to 6 s. We then inverted these data to derive the crustal shear wave velocity using direct surface wave tomography. Meanwhile, a linear sub-array consisting of 54 stations and crossing the surface rupture of the 2008 Wenchuan earthquake was used to image the fault structure by applying the Transmitted Surface Wave Reverse Time Migration method (TSW-RTM). The ambient noise tomography results show that strong lateral heterogeneity exists in the shallow crust in the study area, with a high shear velocity in the northwest and a low shear velocity in the southeast. The low shear velocity is generally distributed between the Yingxiu-Beichuan fault (YBF) and the Guanxian-Jiangyou fault (GJF). The results of TSW-RTM indicate that the Wenchuan-Maoxian fault (WMF) and the YBF are both characterized with a steep dip to the northwest. Our results could serve as the important basis for the study of segmentation characteristics of the LMS fault and seismic hazard preparation and mitigation.

Crustal structure of the Bushveld complex, South Africa from 1D shear wave velocity models: Evidence for complex-wide crustal modification

Thirty-nine 1D shear wave velocity profiles, obtained by jointly inverting receiver functions and Rayleigh wave group velocities, are used to investigate the crustal structure of the Bushveld Complex in northern South Africa. Data from teleseismic earthquakes recorded on broadband seismic stations between 1997 and 1999 and 2015–2020 were used to compute P-wave receiver functions. Rayleigh wave group velocities between 5 and 30 s period were obtained from an ambient noise tomography and combined with group velocities between 30 and 60 s period from a published continental-scale surface wave tomography model. Moho depths of 45–47.5 km are found under the center of the complex compared to 40 km thick crust, on average, surrounding the complex, indicating ∼5–7 of crustal thickening. The bottom ∼10 km or more of the lower crust across much of the Bushveld Complex has a Vs ≥ 4.0 km/s, consistent with a mafic composition, whereas in most areas around the margins of the complex the thickness of the mafic lower crust is much less than 10 km. In the upper crust higher velocity structure (Vs > 3.6 km/s) above 15 km depth underlain by lower velocity structure is seen in many locations, suggesting the presence of mafic/ultramafic layering. These results favor the continuous-sheet model for the structure of the Bushveld Complex because the ensemble of 1D models is characterized by three diagnostic features consistent with that model: (1) thicker crust under the center of the complex than away from the complex; (2) a greater thickness of high-velocity (i.e., mafic) layering in the lower crust under the complex compared to away from the complex; (3) high-velocity (i.e., mafic/ultramafic) layering in the upper crust beneath much of the complex. The lack of upper crustal mafic/ultramafic layering beneath some parts of the complex is consistent with the post-emplacement tectonic and magmatic history of the complex.

Spurious Rayleigh-wave apparent anisotropy near major structural boundaries: a numerical and theoretical investigation

SUMMARY The recent developments in array-based surface-wave tomography have made it possible to directly measure apparent phase velocities through wave front tracking. While directionally dependent measurements have been used to infer intrinsic $2\psi $ azimuthal anisotropy (with a 180° periodicity), a few studies have also demonstrated strong but spurious $1\psi $ azimuthal anisotropy (360° periodicity) near major structure boundaries particularly for long period surface waves. In such observations, Rayleigh waves propagating in the direction perpendicular to the boundary from the slow to the fast side persistently show a higher apparent velocity compared to waves propagating in the opposite direction. In this study, we conduct numerical and theoretical investigations to explore the effect of scattering on the apparent Rayleigh-wave phase velocity measurement. Using 2-D spectral-element numerical wavefield simulations, we first reproduce the observation that waves propagating in opposite directions show different apparent phase velocities when passing through a major velocity contrast. Based on mode coupling theory and the locked mode approximation, we then investigate the effect of the scattered fundamental-mode Rayleigh wave and body waves interfering with the incident Rayleigh wave separately. We show that scattered fundamental-mode Rayleigh waves, while dominating the scattered wavefield, mostly cause short wavelength apparent phase velocity variations that could only be studied if the station spacing is less than about one tenth of the surface wave wavelength. Scattered body waves, on the other hand, cause longer wavelength velocity variations that correspond to the existing real data observations. Because of the sensitivity of the $1\psi $ apparent anisotropy to velocity contrasts, incorporating such measurements in surface wave tomography could improve the resolution and sharpen the structural boundaries of the inverted model.

Frequency‐Bessel Transform Method for Multimodal Dispersion Measurement of Surface Waves From Distributed Acoustic Sensing Data

AbstractThe array‐based frequency‐Bessel transform method has been demonstrated to effectively extract dispersion curves of higher‐mode surface waves from the empirical Green's functions (EGFs) of displacement fields reconstructed by ambient noise interferometry. Distributed acoustic sensing (DAS), a novel dense array observation technique, has been widely implemented in surface wave imaging to estimate subsurface velocity structure in practice. However, there is still no clear understanding in theory about how to accurately extract surface‐wave dispersion curves directly from DAS strain (or strain rate) data. To address this, we extend the frequency‐Bessel transform method by deriving Green’s functions (GFs) for horizontal strain fields, making it applicable to DAS data. First, we test its performance using synthetic GFs and verify the correctness of extracted dispersion spectrograms with theoretical results. Then, we apply it to three field DAS ambient‐noise data sets, two recorded on land and one in the seabed. The reliability and advantages of the method are confirmed by comparing results with the widely used phase shift method. The results demonstrate that our extended frequency‐Bessel transform method is reliable and can provide more abundant and higher‐quality dispersion information of surface waves. Moreover, our method is also adaptable for active‐source DAS data with simple modifications to the derived transform formulas. We also find that the gauge length in the DAS system significantly impacts the polarity and value of extracted dispersion energy. Overall, our study provides a theoretical framework and practical tool for multimodal surface wave dispersion measurement using DAS data.

Lithospheric structure of the Paraná, Chaco-Paraná, and Pantanal basins: Insights from ambient noise and earthquake-based surface wave tomography

To investigate the lithospheric seismic structure of southern South America beneath the Paraná, Chaco-Paraná, and Pantanal basins, we measure two distinct sets of dispersion curves of Rayleigh waves: the first corresponds to 25,986 phase velocity dispersion curves in the period range 5–50 s, extracted from cross-correlation of the ambient noise between pairs of stations and thus sensitive mainly to crustal structure; the second set is derived from group velocities extracted from earthquake recordings, resulting in 33,888 dispersion curves with periods ranging from 8 to 200 s, allowing us to probe deeper Earth structure. From these data sets, we derive two independent S-wave velocity models following a two-step inversion procedure, resulting in a crustal and a lithospheric mantle model with depths extending up to 50 and 200 km, respectively. Features imaged in our crustal model include low velocity anomalies in agreement with the surface position of the major sedimentary basins and high velocity anomalies within fold-thrust provinces and the basement of the Pantanal basin. In the uppermost mantle, we identify high velocities in the region of the Paraná basin, São Francisco and Amazonian cratons, and Luiz Alves block, and a low velocity feature beneath the Pantanal basin, which suggests a potentially thinner lithosphere under the basin. We also illuminate a strong velocity gradient boundary between the eastern border of the Pantanal basin with the Paraná basin, from crustal to lithospheric mantle depths, which seems to delineate the western and southern borders of the Paraná basin, coinciding with the Western Paraná Suture. Additionally, we construct a Moho depth proxy map for the study area, which has an overall good agreement with previous results, except for a thicker crust beneath the southern Chaco-Paraná basin. This area, together with a thicker crust counterpart in the Paraná basin, has a remarkable correlation with the position of a volcanic layer related to the Paraná Magmatic Province, which leads us to suggest that magmatic processes played an important role in the south Chaco-Paraná basin as well.

Potential and limitations of noise-based surface-wave tomography for numerical site effect estimation: a case study in the French Rhône valley

Potential and limitations of noise-based surface-wave tomography for numerical site effect estimation: a case study in the French Rhône valley

Uncertainty Quantification in Radial Anisotropy Models Based on Transdimensional Bayesian Inversion of Receiver Functions and Surface-Wave Dispersion: Case Study Sri Lanka

ABSTRACT In geophysical inference problems, quantification of data uncertainties is required to balance the data-fitting ability of the model and its complexity. The transdimensional hierarchical Bayesian approach is a powerful tool to evaluate the level of uncertainty and determine the complexity of the model by treating data errors and model dimensions as unknown. In this article, we take account of the uncertainty through the whole procedure, thus developing a two-step fully Bayesian approach with coupled uncertainty propagation to estimate the crustal isotropic and radial anisotropy (RA) model based on Rayleigh and Love dispersion as well as receiver functions (RFs). First, 2D surface-wave tomography is applied to determine period-wise ambient noise phase velocity maps and their uncertainty for Rayleigh and Love waves. Probabilistic profiles of the isotropic average VS and RA as a function of depth are then derived at station sites by inverting the local surface-wave dispersion and model errors and RFs jointly. The workflow is applied to a temporary seismic broadband array covering all of Sri Lanka. The probabilistic results enable us to effectively quantify the uncertainty of the final RA model and provide robust inferences. The shear-wave velocity results show that the range of Moho depths is between 30 and 40 km, with the thickest crust (38–40 km) beneath the central Highland Complex. Positive RA (VSH>VSV) observed in the upper crust is attributed to subhorizontal alignment of metamorphic foliation and stretched layers resulting from deformation. Negative RA (VSV>VSH) in the midcrust of central Sri Lanka may indicate the existence of melt inclusions and could result from the uplift and folding process. The positive RA in the lower crust could be caused by crustal channel flow in a collision orogeny.

Computation and analysis of surface wave dispersion and attenuation in layered viscoelastic-vertical transversely isotropic media by the generalized R/T coefficient method

SUMMARY In this study, we propose a systematic and effective method, that is, an extended version of the generalized reflection/transmission (R/T) coefficient method, for computing the phase-velocity (${c}_r$) dispersion curves, attenuation coefficient ($\alpha $) curves, and eigenfunctions of both Rayleigh and Love waves as well as the ellipticity of Rayleigh waves in layered viscoelastic-vertical transversely isotropic (VTI) media. The numerical scheme of combining the root-searching method with the local optimization method is designed for determining the complex-valued modal solutions (i.e. complex wavenumber $k = {\omega {/ {\vphantom {\omega {{c}_r - i\alpha }}}} {{c}_r - i\alpha }}$) of surface waves. The near-surface sedimentary geological environment is taken as the model example because it is typical viscoelastic-VTI media. Besides the anisotropic-viscoelastic (AV) media, our algorithm can also compute surface waves in isotropic-elastic (IE), isotropic-viscoelastic (IV) and anisotropic-elastic (AE) media by resetting the corresponding parameters. Using the six-layer half-space models and in these four media, we verify the correctness of our algorithm by benchmarking the modal solutions against those from other methods. In the four-layer half-space model, by comparing the results of IE, IV, AE and AV media, we analyse the effects of velocity anisotropy, viscoelasticity and attenuation anisotropy on the dispersion and attenuation characteristics of both Rayleigh and Love waves in detail. Our study can provide a theoretical basis and useful tool for surface wave imaging considering the anisotropy and/or viscoelasticity of the medium, which has the potential to better investigate the solid Earth's internal structure.

Choosing Appropriate Regularization Parameters by Splitting Data into Training and Validation Sets—Application in Global Surface-Wave Tomography

Abstract Many linear(ized) geophysical inverse problems cannot be solved without regularization. Finding the regularization parameter that best balances the model complexity and data misfit is often a key step in the inversion problem. Traditionally, this is done by first plotting the measure of model complexity versus data misfit for different values of regularization parameter, which manifests as an L-shaped curve, and then choosing the regularization parameter corresponding to the corner point on the L-curve. For this approach, the difference in units between model complexity and data misfit must be considered, otherwise the result will be strongly affected by the scaling between these two quantities. Inspired by the machine learning literature, we here propose an extension to the traditional L-curve method. We first split the raw dataset into training and validation sets, obtain a solution by performing inversion on the training set only, and calculate data misfits on the validation set. We demonstrate the efficacy of this approach with a toy example and with two synthetic datasets. In realistic global surface-wave tomography studies where sampling of the Earth is nonuniform, we devise a procedure to generate a validation dataset with sampling as uniform as possible. We then show that the regularization parameter can be determined using this validation set, and this determination is apparently robust to the ratio of data split between training and validation sets. For both synthetic tests and realistic inversions, we find that our procedure can produce a minimal point that can be easily identified on the misfit curves calculated on the validation sets, and avoids the nuances encountered in the traditional L-curve analysis.

Illuminating Urban Near-Surface with Distributed Acoustic Sensing Multimodal Noise Surface-Wave Imaging

Abstract Urban subsurface exploration requires high spatial and temporal resolution, cost-effective operation, and minimal interference with urban activities. Distributed acoustic sensing (DAS)—an innovative seismic observation tool—emerges as a promising solution for urban surveys. In this study, we repurposed a 7.9 km telecommunication cable traversing Hefei into a seismic observation array with 3850 channels spaced at 2 m intervals. Noise cross-correlation functions (NCFs) were constructed from recordings by iDAS2 and ZD-DAS interrogators along the entire cable. Spatial variation in the NCFs was observed and attributed to different traffic conditions. Employing the recently developed modified frequency–Bessel transform method to NCFs from the 2 km southern subsection of the optic cable, we extracted broadband, high-resolution multimodal dispersion curves. The inverted near-surface structure beneath the cable unveiled a sediment thinning trend from the center to the periphery of the Hefei basin, consistent with borehole inspections. The three-station interferometry (C3) method and beamforming with the Bessel kernel function are applied to mitigate challenges arising from the weak coupling between the cable and the Earth, as well as persistent localized noise sources. These techniques facilitated the acquisition of broadband surface waves. Distinct secondary scatters are observed in NCFs near channels 2090 and 2287, accompanied by a substantial velocity contrast of 30%–40%, suggesting the existence of a blind fault. The study reaffirms the significant potential of DAS arrays for high-resolution imaging of subsurface structures in challenging urban environments, emphasizing the importance of advanced processing techniques to enhance imaging accuracy and robustness.

Long Period Rayleigh Wave Focal Spot Imaging Applied to USArray Data

AbstractWe demonstrate the effectiveness of seismic dense array surface wave focal spot imaging using USArray data from the western‐central United States. We study dispersion in the 60–310 s period range and assess the image quality of fundamental mode Rayleigh wave phase velocity maps. We apply isotropic spatial autocorrelation models to the time domain zero lag noise correlation wavefield data at distances of about one wavelength. Local estimates of the phase velocity, its uncertainty, and the regression quality imply overall better ZZ relative to ZR or RZ results. The extension of the depth resolution compared to passive surface wave tomography is demonstrated by the inversion of three clustered dispersion curves from different tectonic units. We observe anisotropic surface wave energy flux and the influence of body wave energy, but sensitivity tests at 60 s targeting the data range, correlation component, and processing choices show that the ZZ focal spots yield consistent high‐quality images compared to regional tomography results in the 60–150 s period range. In contrast, at 200–300 s the comparatively small scales of the imaged structures and the imperfect agreement with low‐resolution global tomography results highlight the persistent challenge to reconcile imaging results based on different data sources, theories, and techniques. Our study shows that surface wave focal spot imaging is an accurate, robust, local imaging approach. Better control over clean autocorrelation fields can further improve applications of this seismic imaging tool for increased resolution of the elastic structure below dense seismic arrays.

Seismic structure of the 2015 Mw7.8 Gorkha earthquake revealed by ambient seismic noise and teleseismic surface wave tomography

The destructive 2015 Mw7.8 Gorkha earthquake occurred in the Main Himalayan Thrust due to the collision of the Indian and Asian plates, which provides a unique opportunity to understand the deep dynamic processes and seismogenic mechanisms of strong earthquakes. We construct a regional-scale shear-wave velocity model of the crust and uppermost mantle using ambient seismic noise and teleseismic surface wave at periods of 5–100 s around the Gorkha earthquake region. The new shear-wave velocity model exhibits prominently lateral heterogeneities in the Gorkha earthquake areas. We observe a high-velocity (high-V) zone around the Gorkha main shock in the Main Himalayan Thrust, indicating the existence of a high-strength asperity that sustains the stress accumulating. The aftershocks are primarily located in the low-velocity (low-V) anomalies and enclosed by two high-V anomalies, which appear to act as structural barriers that influence the spread of the aftershocks. Prominent low-Vanomalies from the lower crust to the mantle lithosphere are observed along the north–south trending rifts, suggesting the hot materials upwelling due to the tearing of the northward subducting Indian lithosphere. These observations may indicate that seismic velocity heterogeneity could play an essential role in earthquake initiation and the rupture process.

Comparison of surface-wave techniques to estimate S- and P-wave velocity models from active seismic data

Abstract. The acquisition of seismic exploration data in remote locations presents several logistical and economic criticalities. The irregular distribution of sources and/or receivers facilitates seismic acquisition operations in these areas. A convenient approach is to deploy nodal receivers on a regular grid and to use sources only in accessible locations, creating an irregular source–receiver layout. It is essential to evaluate, adapt, and verify processing workflows, specifically for near-surface velocity model estimation using surface-wave analysis, when working with these types of datasets. In this study, we applied three surface-wave techniques (i.e., wavelength–depth (W/D) method, laterally constrained inversion (LCI), and surface-wave tomography (SWT)) to a large-scale 3D dataset obtained from a hard-rock site using the irregular source–receiver acquisition method. The methods were fine-tuned for the data obtained from hard-rock sites, which typically exhibit a low signal-to-noise ratio. The wavelength–depth method is a data transformation method that is based on a relationship between skin depth and surface-wave wavelength and provides both S- and P-wave velocity (Vs and Vp) models. We used Poisson's ratios estimated through the wavelength–depth method to constrain the laterally constrained inversion and surface-wave tomography and to retrieve both Vs and Vp also from these methods. The pseudo-3D Vs and Vp models were obtained down to 140 m depth over an area of approximately 900 × 1500 m2. The estimated models from the methods matched the geological information available for the site. A difference of less than 6 % was observed between the estimated Vs models from the three methods, whereas this value was 7.1 % for the retrieved Vp models. The methods were critically compared in terms of resolution and efficiency, which provides valuable insights into the potential of surface-wave analysis for estimating near-surface models at hard-rock sites.

Incorporating topographic effects in surface wave tomography based on shortest-path ray tracing

SUMMARY We propose a new method to obtain 3-D shear wave velocity structures for regions with undulating topography based on surface wave dispersion data. In our method, we assume that surface waves propagate along a curved free surface, and the dispersion effects of these waves can be modelled as a frequency-dependent ray tracing problem. We use the shortest-path method to address off-great-circle propagation arising from inhomogeneous elastic parameters and topographic variations. We then apply our method to both synthetic and real-data inversions and demonstrate that ignoring topographic effects may significantly distort the inverted images. Finally, we analysed the accuracy of our method and provided a rule-of-thumb principle to quantitatively assess the need to account for topographic effects in surface wave tomography for a region.