Geophysical Journal International

SummarySeismic full-waveform inversion (FWI) or waveform inversion (WI) has gained extensive attention as a cutting-edge imaging method, which is expected to reveal the high-resolution images of complex geological structures. In this paper, we regard each 1-D signal in the inversion system as a 1-D probability distribution, then use the Jensen-Shannon divergence (JSD) from information theory to measure the discrepancy between the predicted and observed signals, and finally implement a novel 2-D multiparameter shallow-seismic waveform inversion (MSWI). Essentially, the novel approach achieves an implicit weighting along the time-axis for each 1-D adjoint source defined by the classical waveform inversion (CWI), thus enhancing the extra illumination for a deeper medium compared with the CWI. By evaluating the inversion results of the two-layer model and fault model, the reconstruction accuracy for S-wave velocity and density of the new method is increased by about 30% and 20% compared with that of the CWI under the same conditions, respectively. The reconstruction performance for P-wave velocity of these two methods is almost equal. In addition, the new 2-D MSWI is also resilient to white Gaussian noise in the data. Numerically, the inversion system has almost the strongest sensitivities to the S-wave velocity and density, performing the poorest sensitivity to the P-wave velocity. Finally, we test the novel method with a detection case for a power tunnel.

SummaryWe apply the adjoint method to efficiently calculate sensitivity kernels for long-period seismic spectra with respect to structural and source parameters. Our approach is built around the solution of the frequency-domain equations of motion using the Direct Solution Method (DSM). The DSM is currently applied within large-scale mode coupling calculations and is also likely to be useful within finite-element type methods for modelling seismic spectra that are being actively developed. Using mode coupling theory as a framework for solving both the forward and adjoint equations, we present numerical examples that focus on the spectrum close to four eigenfrequencies (the low-frequency mode, 0S2, and higher frequency modes, namely 2S2, 0S7, and 0S10 for comparison). For each chosen observable, we plot sensitivity kernels with respect to 3D perturbations in density and seismic wave-speeds. We also use the adjoint method to calculate derivatives of observables with respect to the matrices occurring within mode coupling calculations. This latter approach points towards a generalisation of the two-stage splitting function method for structural inversions that does not rely on inaccurate self-coupling or group-coupling approximations. Finally, we verify through direct calculation that our sensitivity kernels correctly predict the linear dependence of the chosen observables on model perturbations. In doing this, we highlight the importance of non-linearity within inversions of long-period spectra.

SummaryRapid and accurate characterization of earthquake sources is crucial for mitigating seismic hazards. In this study, based on 18000 scenario ruptures ranging from Mw 6.4 to Mw 8.3 and corresponding synthetic high-rate Global Navigation Satellite System (HR-GNSS) waveforms, we developed a multi-branch neural network framework, the Continental Large Earthquake Agile Respond (CLEAR), to simultaneously determine the magnitude and slip distributions. We apply CLEAR to recent large strike-slip events, including the 2021 Mw 7.4 Maduo earthquake and the 2023 Mw 7.8 and Mw 7.6 Turkey doublet. The model generally estimates the magnitudes successfully at 32 s with errors of less than 0.15, and predicts the slip distributions acceptably at 64 s, requiring only approximately 30 ms on a single CPU. With optimal azimuthal coverage of stations, the system is relatively robust to the number of stations and the time length of the received data.

SummaryLaunched on December 16, 2022, the Surface Water and Ocean Topography (SWOT) satellite, using synthetic aperture radar interferometric techniques, measures sea surface heights (SSHs) across two 50-km-wide swaths, offering high-resolution and accurate two-dimensional SSH observations. This study explores the efficiency of SWOT in seamount detection employing the vertical gravity gradient (VGG) derived from simulated SWOT SSH data. Simulated circular and elliptical seamounts (height: 900-1500 m) are integrated within the South China Sea's 4000 m background depths. Geoid perturbations induced by these seamounts are extracted through the residual depth model principle, subsequently merged with the DTU21MSS model for simulating SWOT SSH observations. For comparative assessment, SSH data from Jason-2 and Cryosat-2 are included. An automatic algorithm (AIFS) is presented to identify seamount centers and base polygons using VGG derived from simulated altimeter SSH data. The analysis reveals SWOT-derived VGGs precisely locate all seamount centers, base polygons, and elliptical seamount azimuths. The merged Jason-2 and Cryosat-2 data face challenges with identifying small circular and elliptical seamounts. Detecting long narrow elliptical seamounts remains arduous; however, SWOT-derived VGGs successfully elucidate the approximate shapes and major axis azimuths of the elliptical seamounts. Validated against ‘true values’ of VGG, the root-mean-squared deviation (RMSD) of SWOT-derived VGG stands at 1.33 Eötvös, whereas the merged Jason-2 and Cryosat-2 data exhibit an RMSD of 1.93 Eötvös. This study shows the enhanced capability of SWOT from its high-resolution two-dimensional SSH observations in advancing seamount detection via satellite-derived VGG. We identify challenges and recommend improved detections using data integration and machine learning.

SummaryBayesian methods provide a valuable framework for rigorously quantifying the model uncertainty arising from the inherent non-uniqueness in the magnetotelluric (MT) inversion. However, widely-used Markov chain Monte Carlo (MCMC) sampling approaches usually require a significant number of model samples for accurate uncertainty estimates, making their applications computationally challenging for 2D or 3D MT problems. In this study, we explore the applicability of the Hamiltonian Monte Carlo (HMC) method for 2D probabilistic MT inversion. The HMC provides a mechanism for efficient exploration in high-dimensional model space by making use of gradient information of the posterior probability distribution, resulting in a substantial reduction in the number of samples needed for reliable uncertainty quantification compared to the conventional MCMC methods. Numerical examples with synthetic data demonstrate that the HMC method achieves rapid convergence to the posterior probability distribution of model parameters with a limited number of model samples, indicating the computational advantages of the HMC in high-dimensional model space. Finally, we applied the developed approach to the COPROD2 field dataset. The statistical models derived from the HMC approach agree well with previous results obtained by 2D deterministic inversions. Most importantly, the probabilistic inversion provides valuable quantitative model uncertainty information associated with the resistivity structures derived from the observed data, which facilitates model interpretation.

SummaryThe Rogaland Igneous Complex (RIC) in southern Norway intruded into Sveconorwegian granulite crust beginning ∼930 Ma. Three massif anorthosite bodies, Egersund-Ogna, Helleren and Åna-Sira, were intruded some 10 myr later by the Bjerkreim-Sokndal layered intrusion. The Garsaknatt leuconorite and the ilmenite-rich Tellnes norite, one of the youngest rock in the complex at ∼920 Ma, intrude the anorthosite or nearby country rock. Magnetic mineralogy and paleomagnetic studies carried out on the Tellnes norite, the Garsaknatt leuconorite, and the surrounding Åna-Sira anorthosite, indicate the magnetization of all three bodies are dominated by hemo-ilmenite carrying the remanence as a thermo-chemical remanent magnetization (TCRM), although magnetite is present in some samples. The three bodies yield steep negative inclinations with northwesterly declinations (Tellnes, I = -71.9°, D = 305.0°, α95 = 10.6°; Garsaknatt, I = -73.1°, D = 312.7°, α95 = 4.7°; Åna-Sira, I = -81.2°, D = 326.3°, α95 = 6.7°). When combined with data from other bodies in the RIC, the older anorthosites have steeper inclinations, and higher paleolatitudes, while the younger units have less steep inclinations and shallower paleolatitudes by nearly 10°, indicating northward plate motion during cooling of the intrusions. Age of the remanence is difficult to determine precisely, however, best estimates are ∼910 Ma for the older anorthosites and ∼900 Ma for the younger intrusions. Although these differences are significant, a unified pole position (35.6°N, 215.1°E), combining all the 111 sites from the RIC, strongly supports the assumed position of southern Baltica in Rodinia at ∼900 MA.

SummarySeismic travel time anomalies of waves that traverse the uppermost 100-200 km of the outer core have been interpreted as evidence of reduced seismic velocities (relative to radial reference models) just below the core-mantle boundary. These studies typically investigate differential travel times of SmKS waves, which propagate as P-waves through the shallowest outer core and reflect from the underside of the core-mantle boundary m times. The use of SmKS and S(m-1)KS differential travel times for core imaging are often assumed to suppress contributions from earthquake location errors and unknown and unmodelled seismic velocity heterogeneity in the mantle. The goal of this study is to understand the extent to which differential SmKS travel times are, in fact, affected by anomalous mantle structure, potentially including both velocity heterogeneity and anisotropy. Velocity variations affect not only a wave's travel time, but also the path of a wave, which can be observed in deviations of the wave's incoming direction. Since radial velocity variations in the outer core will only minimally affect the wave path, in contrast to other potential effects, measuring the incoming direction of SmKS waves provides an additional diagnostic as to the origin of travel time anomalies. Here we use arrays of seismometers to measure travel time and direction anomalies of SmKS waves that sample the uppermost outer core. We form subarrays of EarthScope's regional Transportable Array stations, thus measuring local variations in travel time and direction. We observe systematic lateral variations in both travel time and incoming wave direction, which cannot be explained by changes to the radial seismic velocity profile of the outer core. Moreover, we find a correlation between incoming wave direction and travel time anomaly, suggesting that observed travel time anomalies may be caused, at least in part, by changes to the wave path and not solely by perturbations in outer core velocity. Modelling of 1-D ray and 3-D wave propagation in global 3-D tomographic models of mantle velocity anomalies match the trend of the observed travel time anomalies. Overall, we demonstrate that observed SmKS travel time anomalies may have a significant contribution from 3-D mantle structure, and not solely from outer core structure.

SummaryIncreasing evidence from seismic methods shows that anisotropy within subduction zones should consist of multiple layers. To test this, we calculate and model shear wave splitting across the Alaska-Aleutians Subduction Zone (AASZ), where previous studies have argued for separate layers of anisotropy in the subslab, slab, and mantle wedge. We present an updated teleseismic splitting catalog along the span of the AASZ, which has many broadband seismometers recently upgraded to three components. Splitting observations are sparse in the Western Aleutians, and fast directions are oriented generally trench parallel. There are significantly more splitting measurements further east along the AASZ. We identify six regions in the Central and Eastern Aleutians, Alaskan Peninsula, and Cook Inlet with a high density of splits suitable for multilayered anisotropy analyses. These regions were tested for multilayer anisotropy, and for five of the six regions we favor multiple layers over a single layer of anisotropy. We find that the optimal setup for our models is one with a dipping middle layer oriented parallel to paleospreading. A prominent feature of our modeling is that fast directions above and below the dipping layer are generally oriented parallel to the strike of the slab. Additionally, we lay out a framework for robust and statistically reliable multilayer shear wave splitting modeling.

SummaryAmbient noise tomography on the basis of Distributed Acoustic Sensing (DAS) deployed on existing telecommunication networks provides an opportunity to image the urban subsurface at regional scales and high-resolution. This capability has important implications in the assessment of the urban subsurface’s potential for sustainable and safe utilization, such as geothermal development. However, extracting coherent seismic signals from the DAS ambient wavefield in urban environments at low cost remains a challenge. One obstacle is the presence of complex sources of noise in urban environments, which may not be homogeneously distributed. Consequently, long recordings are required for the calculation of high-quality virtual shot gathers, which necessitates significant time and computational cost. In this paper, we present the analysis of 15 days of DAS data recorded on a pre-existing fiber optic cable (dark fibers), running along an 11 km long major road in urban Berlin (Germany), hosting heavy traffic including vehicles and trains. To retrieve virtual shot gathers, we apply interferometric analysis based on the cross-correlation approach where we exclude low-quality virtual shot gathers to increase the signal-to-noise ratio of the stacked gathers. Moreover, we modify the conventional ambient noise interferometry workflow by incorporating a coherence-based enhancement approach designed for wavefield data recorded with large-N arrays. We then conduct Multichannel Analysis of Surface Waves (MASW) to retrieve 1D velocity models for two exemplary fiber subsegments, and compare the results of the conventional and modified workflows. The resulting 1D velocity models correspond well with available lithology information. The modified workflow yields improved dispersion spectra, particularly in the low-frequency band (< 1 Hz) of the signal. This leads to an increased investigation depth along with lower uncertainties in the inversion result. Additionally, these improved results were achieved using significantly less data than required using conventional approaches, thus opening the opportunity for shortening required acquisition times and accordingly lowering costs.

SummaryThis study attempts to interrogate the upper mantle deformation pattern beneath the Kumaon-Garhwal region, located in the western Himalaya, using shear wave splitting (SWS) analysis of core-refracted (XK(K)S) phases recorded at 53 broadband stations. The fast polarisation azimuths (FPAs) revealed by 338 well constrained measurements are dominantly clustered around ENE-WSW, with a few along the NE and E-W directions. The delay times vary from 0.2 to 1.4 s, with an average of 0.6 s that is smaller than that for the Indian shield (∼0.8 s), central and eastern Himalayas. The northern part of the lesser Himalaya shows a slightly smaller delay time compared to the southern part, which is attributed to the weakening of azimuthal anisotropy caused by the dipping of the Indian lithosphere. In order to understand the crustal contribution, its anisotropy is measured by analysing the splitting of Ps conversions from the Moho (Pms), akin to that of the XK(K)S phases. However, reliable results for crustal anisotropy could be obtained only at 10 stations. The average delay time due to crustal anisotropy is 0.47 s, with a variation from 0.2 to 0.9 s. Although the dominant period of Pms is smaller than that of SK(K)S, crustal anisotropy contributing to splitting of the latter phases cannot be ruled out. The orientation of FPAs obtained from Pms phases is found to be parallel or sub-parallel to those from XK(K)S phases, suggesting a similar deformation mechanism in the mid- to lower-crust and upper mantle. On the basis of FPAs derived from XK(K)S measurements, the Kumaon-Garhwal Himalaya (KGH) region can be divided into four sub-regions. In the western and eastern parts, the FPAs are mostly aligned along NE and ENE-WSW, and NE, respectively. In the central and south-eastern parts, their orientation is along ENE-WSW and NW, respectively. The strong ENE-WSW orientation in the central part could result from a slightly variable anisotropy in the crust to the upper part of the lithosphere or basal topography causing deflection of mantle flow. Also, the NW orientation in the south-eastern part of KGH is associated with a shallow source within the lithosphere. Application of the spatial coherency technique to single-layered anisotropic parameters results in a depth of 220-240 km, implying that the dominant source of anisotropy could lie in the upper mantle.

SummaryModeling seismic wavefields in complex 3D elastic media is the key in many fields of Earth Science: seismology, seismic imaging, seismic hazard assessment, earthquake source mechanism reconstruction. This modeling operation can incur significant computational cost, and its accuracy depends on the ability to take into account the scales of the subsurface heterogeneities varying. The theory of homogenization describes how the small-scale heterogeneities interact with the seismic waves and allows to upscale elastic media consistently with the wave equation. In this study, an efficient and scalable numerical homogenization tool is developed, relying on the similarity between the equations describing the propagation of elastic waves and the homogenization process. By exploiting the optimized implementation of an elastic modeling kernel based on a spectral element discretization and domain decomposition, a fully scalable homogenization process, working directly on the spectral-element mesh, is presented. Numerical experiments on the entire SEAM II foothill model and a 3D version of the Marmousi II model illustrate the efficiency and flexibility of this approach. A reduction of two orders of magnitude in terms of absolute computational cost is observed on the elastic wave modeling of the entire SEAM II model at a controlled accuracy.

SummaryThe increasing need for petrophysical parameter inversion has spurred the advancement of full-waveform inversion (FWI). This study introduces a new inversion method that can directly estimate petrophysical parameters from seismic waveform data. Specifically, we modified the regular elastic wave equation by including Biot poroelastic parameters and designed an inversion workflow based on the framework of elastic FWI. The newly introduced inversion algorithms have demonstrated their effectiveness and accuracy using the Marmousi model and a synthetic model related to the gas hydrate deposits based on the Shenhu area of the South China Sea. Our work demonstrates the feasibility of directly inverting petrophysical parameters and assessing if a deposit contains gas hydrate.

SummaryTime-lapse seismic full-waveform inversion (FWI) provides estimates of dynamic changes in the Earth’s subsurface by performing multiple seismic surveys at different times. Since FWI problems are highly non-linear and non-unique, it is important to quantify uncertainties in such estimates to allow robust decision making based on the results. Markov chain Monte Carlo (McMC) methods have been used for this purpose, but due to their high computational cost, those studies often require a pre-existing accurate baseline model and estimates of the locations of potential velocity changes, and neglect uncertainty in the baseline velocity model. Such detailed and accurate prior information is not always available in practice. In this study we use an efficient optimization method called stochastic Stein variational gradient descent (sSVGD) to solve time-lapse FWI problems without assuming such prior knowledge, and to estimate uncertainty both in the baseline velocity model and the velocity change over time. We test two Bayesian strategies: separate Bayesian inversions for each seismic survey, and a single join inversion for baseline and repeat surveys, and compare the methods with standard linearised double difference inversion. The results demonstrate that all three methods can produce accurate velocity change estimates in the case of having fixed (exactly repeatable) acquisition geometries. However, the two Bayesian methods generate significantly more accurate results when acquisition geometries changes between surveys. Furthermore, joint inversion provides the most accurate velocity change and uncertainty estimates in all cases tested. We therefore conclude that Bayesian time-lapse inversion using a joint inversion strategy may be useful to image and monitor subsurface changes, in particular where variations in the results would lead to different consequent decisions.

SummaryDensity is an important parameter for both geological research and geophysical exploration. However, for model-driven seismic inversion methods, high-fidelity density inversion is challenging due to seismic wave travel-time insensitivity to density, and crosstalk that density has with velocity. To circumvent the challenge of density inversion, some inversion methods treat density as a constant value or derive density from velocity through empirical equation. On the other hand, deep learning approaches are completely driven by data and have strong target-oriented characteristics, offering a new way to solve multi-parameter coupling problems. Nevertheless, the accuracy of the inversion results of data-driven algorithms is directly related to the amount and diversity of the training data, and thus, they lack the universality of model-driven algorithms. To achieve accurate density inversion, we propose a simultaneous inversion algorithm for velocity and density that combines the advantages of data- and model- driven approaches: A neural network model (U-T), combining the U-net and Transformer architectures, is proposed to construct nonlinear mappings between seismic data as inputs and the velocity and density as predictions. Next, the model-driven inversion algorithm uses the U-T prediction as the initial model to obtain the final accurate solution. In the model-driven module, envelope-based sparse constrained deconvolution is used to obtain full-band seismic data, while a variable dominant frequency full waveform inversion algorithm is employed to perform multi-scale inversion, ultimately leading to accurate inversion results of velocity and density. The performance of the algorithm on the Sigsbee2A and Marmousi models demonstrates its effectiveness.

SummaryThe Bjerkreim-Sokndal intrusion in southern Norway has been studied for decades due to the presence of magnetic remanence creating anomalies 12000 nT below background as measured by airborne magnetic surveys. The strong magnetic remanence also makes the Bjerkreim-Sokndal intrusion a good Earth analogue for remote studies of planets that have prominent magnetic signatures, such as Martian geologic environments. Although numerous geophysical surveys and samples have been collected in the area, there are limited 3D geological interpretations of the subsurface. Here, we used existing geophysical data to conduct forward and inversion modeling of the Bjerkreim lobe to investigate the subsurface geometry of the Bjerkreim-Sokndal intrusion. An extensive petrophysical property compilation was used as input data for the models, in combination with airborne magnetics and digital elevation models. This petrophysical compilation was initially analyzed using Principal Component Analysis to understand which variables would have the greatest impact on the models. Forward and inversion modeling show that crosscutting jotunite bodies, and small anorthosite blocks within the Bjerkreim lobe have a limited depth extent of 1 km. Massive and foliated anorthosites to the west of the Bjerkreim lobe extend to depths greater than 4 km indicating that the Bjerkreim-Sokndal intruded into these anorthosites. Complications in magnetic field fitting during the forward modeling of megacyclic units with strong magnetic remanence and the results from a new ground magnetic survey support the need to revisit mapped contacts of the cyclical units.

SummaryWe expand the scope of HeFESTo by encompassing the rich physics of iron in the mantle, including the existence of multiple valence and spin states. In our previous papers, we considered iron only in its most common state in the mantle: the high-spin divalent (ferrous) cation. We now add ferric iron end-members to six phases, as well as the three phases of native iron. We also add low-spin states of ferrous and ferric iron and capture the behavior of the high-spin to low-spin transition. Consideration of the multi-state nature of iron, unique among the major elements, leads to developments of our theory, including generalization of the chemical potential to account for the possibility of multiple distinguishable states of iron co-existing on a single crystallographic site, the effect of the high-spin to low-spin transition on seismic wave velocities in multi-phase systems, and computation of oxygen fugacity. Consideration of ferric iron also motivates the addition of the chromia component to several phases, so that we now consider the set of components: Ca, Na, Fe, Mg, Al, Si, O, and Cr (CNFMASO+Cr). We present the results of a new global inversion of mineral properties and compare our results to experimental observations over the entire pressure-temperature range of the mantle and over a wide range of oxygen fugacity. Applications of our method illustrate how it might be used to better understand the seismic structure, dynamics, and oxygen fugacity of the mantle.

SummaryThe seismic b-value in Gutenberg–Richter law is an important parameter in earthquake science research and earthquake risk assessment. People have tried to use b-values for short-term earthquake forecasts, and this requires the premise of estimating reliable b-values for real-time seismic catalogs. However, estimating b-values for real-time catalogs, which are usually of poor qualities, is usually faced with many difficulties and problems. In this study, through a series of numerical tests, we investigate the performance of three methods, including the commonly used maximum likelihood estimation method and two relatively new b-value estimation methods, namely the b-positive and K − M slope methods, on calculating b-values for real-time seismic catalogs. We also apply these three methods to both observed seismic catalogs (the seismic sequence in the Costa Marchigiana, Italy, and a high-resolution early aftershock sequence of the 2023 two Mw ∼7.8 earthquakes in Türkiye) and synthetic real-time seismic catalogs. The results in this study show that it seems difficult to obtain accurate b-values for real-time earthquake catalogs by each of these three methods, but the combination of these methods may give a better judgment—if all three methods suggest that the change in b-value is significant, the probability of making a correct decision is very high. Facing the uncertainty of b-value estimation that always exists, we advocate exploring the effectiveness of standard b-value estimation strategies in practical applications.

SummaryThe Lacq area in southwest France has been associated with continuous moderate induced seismic activity since 1969. However, the mechanisms driving this induced seismicity are not fully understood: reservoir depletion has been proposed as the main factor, and more recently wastewater injection has been suggested to play a more important role (Grasso et al., 2021). The interpretation of these mechanisms relies heavily on the quality of earthquake locations, which we prove to be weak due to a lack of local instrumentation for several years. In order to provide the most complete and reliable induced event catalog for the studies of the Lacq induced seismicity mechanisms & seismic hazard, we made an exhaustive compilation, analysis and improvement of all available catalogs. We also provided new earthquake detections & relocations in a 3D velocity model from past and present temporary deployments never used for studying the Lacq area. Important remaining location uncertainties lead us to also carefully sort the events according to their location confidence, defining 3 classes of events (unconstrained location, location constrained within 2-3 km and 1-2 km respectively). This new harmonized catalog and the identification of well-constrained events, covering 50 years of induced seismicity, allow us to propose that wastewater injection is almost certainly the main mechanism driving the seismicity, with (i) most of the constrained events located within the reservoir boundaries and (ii) the released seismic energy variations following variations in injection operations at different scales. In particular, we have also highlighted a change in the injection-seismicity relationship around 2010–2013. From 2013, despite lower injection volumes, seismicity remained persistent and some clusters of earthquakes were detected predominantly in spring, summer, and early autumn, except in winter periods. From 2016, we observed a strong temporal relationship between days with higher rate/volume injections (approximately above 400m3/day) and both clustered events and higher magnitude earthquakes (greater than 2.4).

SummaryThe Arunachal and Bhutan Himalaya, which are tectonically distinct from other regions of the Himalaya, have a structure that is quite intricate. The eastern Himalayan segment is a component of the region where the Indian and Eurasian plates collided 50 million years ago. The Indian plate goes beneath the Eurasian plate in the north, and in the eastern part of the region, the Indian plate subducts under the Burmese plate. Here, we studied the seismic attenuation of the uppermost mantle by measuring the quality factor of the Sn wave (SnQ) to understand the dynamics of the lithospheric mantle and the cause of the seismic anomalies found in this area. The upper mantle Q structure has significant lateral differences in Arunachal and the Bhutan Himalaya. Arunachal Himalaya’s central region is characterised by a very low Q ( ≤ 150). The successive low-high-low SnQ values in eastern Arunachal Himalaya near Siang region have been observed. The western Arunachal region, close to the Bhutan border, exhibits a contrast in Q values. We notice that low Q values ( ≤ 200) predominate in the central to eastern Bhutan Himalaya. The western part of Bhutan Himalaya exhibits relatively high Q ( ≥ 200) values, mostly near Paro and Thimpu. Interestingly, a clear boundary between low and high Q has been observed near Kakthang thrust (KT) in the Bhutan Himalaya. We found significant lateral variation of frequency dependent parameter (η) across the study region. They range from 0.25 to 0.75, with low values ( ≤ 0.5) found mostly in the central Bhutan Himalaya and in a few isolated areas of the Arunachal Himalaya. Low Q and a relatively higher η ( ≥ 0.5) might suggest that the scattering attenuation is the controlling mechanism for Sn wave attenuation in the upper mantle beneath Arunachal Himalaya. On the contrary, dominant low Q values across the central segment of the Bhutan Himalaya, along with a low to moderate body wave velocity and dominating low η values, subsequently corroborate that intrinsic attenuation is the dominant factor in the upper mantle of the central Bhutan Himalaya.

SummaryDespite the increasing accuracies of GRACE/GRACE-FO gravity field models through worldwide endeavors, the temporal aliasing effect caused by the imperfect background models used in gravity field modelling is still a crucial factor that degrades the quality of gravity field solutions. Since the important role of temporal resolution of atmospheric de-aliasing models, this paper specifically investigates the influence of temporal resolution on gravity field modeling from the perspectives of frequency, spectral, and spatial domains. To this end, we introduced the gravitational acceleration and geoid height derived from the static gravity field GOCO06s in the inner integral. The introduction of the static gravity field has a comparable impact on LRI range-rate residuals as the accuracy of the LRI range-rate data, despite its magnitude of being less than 0.1mm in the spatial domain. This finding also highlights the significance of error level in existing de-aliasing products as a crucial factor that restricts the current accuracy of gravity field solutions. Further analyses show that increasing the temporal resolution from 3 hours to 1 hour has an insignificant impact on the gravity solutions in both the frequency and spectral domains, which is also smaller than that caused by using different atmospheric datasets. However, in the spatial domain, LRI range-rate residuals can be effectively mitigated in certain regions of the Southern Hemisphere at mid- and high-latitudes by increasing the temporal resolution. Particularly, the discrepancies of mass change estimates brought about by enhancing temporal resolution have distinct characteristics, especially in the Congo River and the Amazon River Basins. The mass changes in terms of EWH derived by using P4M6 filtering show that the maximum RMS value of spatial differences caused by improving the temporal resolution of the atmospheric de-aliasing models can reach ∼13.4mm in the sub-region of the Congo River Basin. However, using different atmospheric datasets can lead to a maximum difference of ∼16.5mm. For the Amazon River Basin, the corresponding maximum discrepancy is ∼18.1mm, and that caused by improving temporal resolution is ∼9.4mm. We further divide the Congo River Basin into several sub-regions using a lat-lon regular grid with a spatial resolution of 3 degrees. The subsequent time series results of mass changes reveal that the maximum contribution of temporal resolution and changes in the atmospheric datasets can reach 11.09% and 21.24%, respectively. This suggests that it is necessary to consider the temporal resolution of de-aliasing products when studying mass changes at a regional scale.