Search

The estimation of the slope (b -value) of the frequency–magnitude distribution of earthquakes is based on a formula derived by Aki decades ago, assuming a continuous exponential distribution. However, as the magnitude is usually provided with a limited resolution, its distribution is not continuous but discrete. In the literature, this problem was initially solved by an empirical correction (due to Utsu) to the minimum magnitude, and later by providing an exact formula such as that by Tinti and Mulargia, based on the geometric distribution theory. A recent paper by van der Elst showed that the b -value can be estimated also by considering the magnitude differences (which are proven to follow an exponential discrete Laplace distribution) and that in this case the estimator is more resilient to the incompleteness of the magnitude data set. In this work, we provide the complete theoretical formulation including (i) the derivation of the means and standard deviations of the discrete exponential and Laplace distributions; (ii) the estimators of the decay parameter of the discrete exponential and trimmed Laplace distributions and (iii) the corresponding formulas for the parameter b. We deduce (iv) the standard 1 σ intervals for the estimated b. Moreover, we are able (v) to quantify the error associated with the Utsu minimum-magnitude correction. Furthermore, we have discussed the formulas to produce statistically independent magnitude differences. We tested extensively the b -value estimators on simulated synthetic data sets including complete catalogues as well as catalogues affected by a strong incompleteness degree such as aftershock sequences where the incompleteness is made to vary from one event to the next. We have also analysed the real aftershock sequence of the 30/10/2016 Norcia (central Italy) to integrate the finding of the simulations. To judge the performance of the various estimators we have introduced an index p that can be seen as a non-parametric extension of the Student's t index. The main outcomes of this paper are that (1) the b -value estimators devised for continuous magnitude data are not adequate for binned magnitudes, (2) for complete data sets, estimators based on magnitudes and on magnitude differences provide substantially equivalent results, (3) for incomplete magnitude data sets, estimators based on magnitude differences provide better results and (4) for incomplete aftershock sequences there is no evidence that methods based on positive magnitude differences are superior than other methods using differences. This conclusion is further confirmed by our analysis of the above-mentioned Norcia seismic sequence. This last finding contrasts with the van der Elst's claim that the so called |${{b}_ + }$| method is the most adequate to treat real aftershock sequences. [ABSTRACT FROM AUTHOR]

Full-waveform inversion (FWI) is an effective method for imaging subsurface properties using sparsely recorded data. It involves solving a wave propagation problem to estimate model parameters that accurately reproduce the data. Recent trends in FWI have seen a renewed interest in extended methodologies, among which source extension methods leveraging reconstructed wavefields to solve penalty or augmented Lagrangian (AL) formulations have emerged as robust algorithms, even for inaccurate initial models. Despite their demonstrated robustness on synthetic data, challenges remain, such as the lack of a clear physical interpretation and reliance on difficult-to-compute least-squares (LS) wavefields. Moreover, the literature lacks a general and through comparison of these source extension methods with each other and with the standard FWI. This paper is divided into three critical parts. In the first, a novel formulation of these methods is explored within a unified Lagrangian framework. This novel perspective permits the introduction of alternative algorithms that use LS multipliers instead of wavefields. These multiplier-oriented variants appear as regularizations of the standard FWI, are suitable to the time domain, offer tangible physical interpretations, and foster enhanced convergence efficiency. The second part of the paper delves into understanding the underlying mechanisms of these techniques. This is achieved by solving the associated non-linear equations using iterative linearization and inverse scattering methods. The paper provides insight into the role and significance of Lagrange multipliers in enhancing the linearization of the equations. It explains how different methods estimate multipliers or make approximations to increase computing efficiency. Additionally, it presents a new physical understanding of the Lagrange multiplier used in the AL method, highlighting how important it is for improving algorithm performance when compared to penalty methods. In the final section, the paper presents numerical examples that compare different methods within a unified iterative algorithm, utilizing benchmark Marmousi and 2004 BP salt models. [ABSTRACT FROM AUTHOR]

More accurate inversion of source fault geometry and slip parameters under the constraint of the Bayesian algorithm has become a research hotspot in the field of geodetic inversion in recent years. In nonlinear inversion, the determination of the weight ratio of the joint inversion of multisource data is more complicated. In this context, this paper proposes a simple and easily generalized weighting method for inversion of source fault parameters by joint geodetic multisource data under the Bayesian framework. This method determines the relative weight ratio of multisource data by root mean square error (RMSE) value and can be extended to other nonlinear search algorithms. To verify the validity of the method in this paper, this paper first sets up four sets of simulated seismic experiment schemes. The inversion results show that the joint inversion weighting method proposed in this paper has a significant decrease in the large residual value compared with the equal weight joint inversion and the single data source joint inversion method. The east–west deformation RMSE is 0.1458 mm, the north–south deformation RMSE is 0.2119 mm and the vertical deformation RMSE is 0.2756 mm. The RMSEs of the three directions are lower than those of other schemes, indicating that the proposed method is suitable for the joint inversion of source parameters under Bayesian algorithm. To further verify the applicability of the proposed method in complex earthquakes, the source parameters of the Maduo earthquake were inverted using the method of this paper. The focal depth of the inversion results in this paper is closer to the focal depth released by the GCMT agency. In terms of strike angle and dip angle, the joint inversion in this paper is also more inclined to the GCMT results. The joint inversion results generally conform to the characteristics of left-lateral strike-slip, which shows the adaptability of this method in complex earthquakes. [ABSTRACT FROM AUTHOR]

SUMMARY: Before inverting Moho topography, the traditional Parker–Oldenburg method requires the determination of two important hyperparameters, the average Moho depth and Moho density contrast. The selection of these two hyperparameters will directly affect the inversion results. In this paper, a new method for estimating hyperparameters is proposed which is used to improve the Parker–Oldenburg method. The new method is improved by using simulated annealing to accurately estimate the average Moho depth and Moho density contrast based on the relationship between Moho depths and corresponding gravity anomalies at seismic control points. Synthetic tests show that compared to the improved Bott's method and the trial and error method, our method reduces the error in Moho density contrast and average Moho depth by 0.83 and 1.81 per cent, respectively. In addition, compared with the trial and error method, our method greatly improves the computational efficiency. In a practical example, we apply this method to invert the Moho topography in the northern South China Sea. The inversion results show that the Moho topography in the northern South China Sea ranges from 8.2 to 33 km. The root mean squared error between our Moho topography and the seismic validation points is 0.94 km. Compared with the CRUST 1.0 model, our Moho topography is more accurate. [ABSTRACT FROM AUTHOR]

SUMMARY: For accurate modelling of groundwater flow and transport processes within an aquifer, precise knowledge about hydraulic conductivity K and its small-scale heterogeneities is fundamental. Methods based on pumping tests, such as hydraulic tomography (HT), allow for retrieving reliable K-estimates, but are limited in their ability to image structural features with high resolution, since the data from time-consuming hydraulic tests are commonly sparse. In contrast, geophysical methods like induced polarization (IP) can potentially yield structural images of much higher resolution, but depend on empirical petrophysical laws that may introduce significant uncertainties to the K-estimation. Therefore, this paper presents a joint inversion procedure for both HT and IP data, which allows for combining the complementary abilities of both methods. Within this approach, a traveltime inversion is applied to the HT data, while the IP inversion is based on a full-decay time-domain forward response, as well as a reparametrization of the Cole–Cole model to invert for K directly. The joint inversion is tested on a synthetic model mimicking horizontally layered sediments, and the results are compared with the individual HT and IP inversions. It is shown that jointly inverting both data sets consistently improves the results by combining the complementary sensitivities of the two methods, and that the inversion is more robust against changes in the experimental setups. Furthermore, we illustrate how a joint inversion approach can correct biases within the petrophysical laws by including reliable K-information from hydraulic tests and still preserving the high-resolution structural information from IP. The different inversion results are compared based on the structural similarity index (SSIM), which underlines the robustness of the joint inversion compared to using the data individually. Hence, the combined application of HT and IP within field surveys and a subsequent joint inversion of both data sets may improve our understanding of hydraulically relevant subsurface structures, and thus the reliability of groundwater modelling results. [ABSTRACT FROM AUTHOR]

SUMMARY: The amplitude, frequency and polarization of ground motion at the surface can be affected by the local geology. While low-velocity sediments and fill can amplify ground motions in certain frequency ranges, the low velocities found in fault zones can also produce prominent wavelets. In this paper, we provide further evidence that polarization of ground motion can be affected by the geological fabric in fault zones that have sustained significant brittle deformation. Aside from the well-known effect of fault-trapped waves in the low-velocity zone with polarization azimuths parallel to the fault strike, the effect of stiffness anisotropy was recently recognized with polarization azimuths at high-angle to the fault strike and orthogonal to the locally predominant fracture field in the fault damage zone. To clarify further such features, we investigate directional amplification effects across the San Jacinto fault zone in Southern California using seismic data recorded by permanent seismic stations and dense across-fault arrays. We observe three main polarization trends. The first trend parallel to the fault strike is ascribed to fault-trapped waves along the low-velocity zone, in agreement with several studies in the last decade in the same region. The second and third trends are orthogonal to the orientation of R and T Riedel planes, respectively. They are related to the stiffness anisotropy in densely fractured rocks in the damage zone, which are more compliant orthogonal to their fractures. At some locations the two effects are superimposed, occurring in different and distinct frequency ranges. Directional amplification at rock sites can be important for expected ground motion and seismic hazard. However, in seismic engineering the current prescriptions of seismic codes do not account for amplification effects at rock sites at frequencies of engineering interest. [ABSTRACT FROM AUTHOR]

SUMMARY: Atmospheric pressure changes on Earth's surface can deform the solid Earth. Sorrells derived analytical formulae for displacement in a homogeneous, elastic half-space, generated by a moving surface pressure source with speed $c$. Ben-Menahem and Singh derived formulae when an atmospheric P wave impinges on Earth's surface. For a P wave with an incident angle close to the grazing angle, which essentially meant a slow apparent velocity $c_a$ in comparison to P- ($\alpha ^{\prime }$) and S-wave velocities ($\beta ^{\prime }$) in the Earth ($c_a \ll \beta ^{\prime } \lt \alpha ^{\prime }$), they showed that their formulae for solid-Earth deformations become identical with Sorrells' formulae if $c_a$ is replaced by $c$. But this agreement was only for the asymptotic cases ($c_a \ll \beta ^{\prime }$). The first point of this paper is that the agreement of the two solutions extends to non-asymptotic cases, or when $c_a /\beta ^{\prime }$ is not small. The second point is that the angle of incidence in Ben-Menahem and Singh's problem does not have to be the grazing angle. As long as the incident angle exceeds the critical angle of refraction from the P wave in the atmosphere to the S wave in the solid Earth, the formulae for Ben-Menahem and Singh's solution become identical to Sorrell's formulae. The third point is that this solution has two different domains depending on the speed $c$ (or $c_a$) on the surface. When $c/\beta ^{\prime }$ is small, deformations consist of the evanescent waves. When $c$ approaches Rayleigh-wave phase velocity, the driven oscillation in the solid Earth turns into a free oscillation due to resonance and dominates the wavefield. The non-asymptotic analytical solutions may be useful for the initial modelling of seismic deformations by fast-moving sources, such as those generated by shock waves from meteoroids and volcanic eruptions because the condition $c / \beta ^{\prime } \ll 1$ may be violated for such fast-moving sources. [ABSTRACT FROM AUTHOR]

SUMMARY: The geocentric axial dipole (GAD) hypothesis is key to many palaeomagnetic applications, for example plate-tectonic reconstructions; however, the validity of this hypothesis at high latitudes is not fully resolved. To address this, in this paper we determined the palaeomagnetic directional data of 156 lava units in Eyjafjarðardalur, Iceland, with the aim of determining the validity of the GAD hypothesis at high latitudes using time-averaged field (TAF) analysis. In addition to the palaeomagnetic directional data, we constructed an age model for the sequences using new 40Ar/39Ar dates, magnetostratigraphy and field data. The sequence age range is 2.6–8.0 Ma. We show that the mean virtual geomagnetic pole (VGP) for our data does not agree with the GAD theory at 95 per cent confidence, when only the standard tilt and tectonic corrections are made; however, when inclination-shallowing processes are accounted for, for example thermoremanence (TRM) anisotropy and refraction effects, the mean VGP can align with GAD at 95 per cent confidence. These inclination-shallowing processes are shown to reduce the inclination by up to 14° for some of the basaltic units. Applying the inclination-shallowing correction also reduces VGP dispersion to levels that agree with global model predictions. We propose that much of the scatter within the palaeomagnetic directional databases are due to inclination-shallowing process effects, which become more important as the natural remanent magnetization (NRM) intensity is high, for example >2 A m−1. We propose that inclination-shallowing processes can be identified and corrected for by examining the NRM intensity and dispersion. [ABSTRACT FROM AUTHOR]

This 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. [ABSTRACT FROM AUTHOR]

A 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 multimode 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 modelling 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. [ABSTRACT FROM AUTHOR]

In this paper, we use a new workflow to substantiate the characterization of a prominent, deep sediment conductor in the hyperextended Bjørnøya Basin (SW Barents Sea) previously identified in smooth resistivity models from 3-D deterministic inversion of magnetotelluric data. In low-dimensionality environments like layered sedimentary basin, 1-D Bayesian inversion can be advantageous for a thorough exploration of the solution space, but the violation of the 1-D assumption has to be efficiently handled. The primary geological objectives of this work is therefore preceded by a secondary task: the application of a new machine learning approach for handling the 1-D violation assumption for 21 MT field stations in the Barents Sea. We find that a decision tree can adequately learn the relationship between MT dimensionality parameters and the 1-D–3-D residual response for a training set of synthetic models, mimicking typical resistivity structures of the SW Barents Sea. The machine learning model is then used to predict the dimensionality compensation error for MT signal periods ranging of 1–3000 s for 21 receivers located over the Bjørnøya Basin and Veslemøy High. After running 1-D Bayesian inversion, we generated a posterior resistivity distribution for an ensemble of 6000 1-D models fitting the compensated MT data for each 21 field stations. The proportion of 1-D models showing ρ < 1 Ω·m is consistently beyond 80 per cent and systemically reaches a maximum of 100 per cent in the Early Aptian–Albian interval in the Bjørnøya Basin. In hyperextended basins of the SW Barents Sea, the dimensionality compensation workflow has permitted to refine the characterization of the deep basin conductor by leveraging the increased vertical resolution and optimal used of MT data. In comparison, the smooth 3-D deterministic models only poorly constrained depth and lateral extent of the basin anomaly. The highest probability of finding ρ < 1 Ω·m is robustly assigned to the syn-tectonic Early Aptian–Albian marine shales, now buried at 6–8 km depth. Based on a theoretical two phase fluid-rock model, we show that the pore fluid of these marine shales must have a higher salinity than seawater to explain the anomaly ρ < 1 Ω·m. Therefore, the primary pore fluid underwent mixing with a secondary brine during rifting. Using analogue rift systems in palaeomargins, we argue that two possible secondary brine reservoir may contribute to deep saline fluid circulation in the hyperextended basin: (1) Permian salt-derived fluid and, (2) mantle-reacted fluid from serpentinization. [ABSTRACT FROM AUTHOR]

Seismic 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 from information theory to measure the discrepancy between the predicted and observed signals, and finally implement a novel 2-D multiparameter shallow-seismic WI (MSWI). Essentially, the novel approach achieves an implicit weighting along the time-axis for each 1-D adjoint source defined by the classical WI (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 per cent 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. [ABSTRACT FROM AUTHOR]

This study estimated the gravity change expected on the sea surface due to magma accumulation beneath submarine volcanoes. For calculation, a semi-analytical model describing the deformation of an infinite elastic cone was considered, in which a point spherical source was embedded on the axis. The expected gravity change exceeds 0.2 mGal in certain scenarios when |$\delta V = 10^8\, \mathrm{m^3}$| of source volume change occurs 2 km below the summit, and the summit is shallower than 300 m below the sea surface. The steeper the slope of a submarine volcano, the greater the expected gravity change. 0.2 mGal of gravity variation can be detected with recent marine gravimeters such as cold atom interferometers. The computation method for the gravity change based on the integral transform (Mellin type) is elaborated in this paper. [ABSTRACT FROM AUTHOR]

Shear wave velocity (Vs) is a fundamental property of elastic media whose estimation from PS converted waves is challenging and requires modelling the boundary where P to S conversion occurs. This paper presents a PS tomography where seismic wave conversion/reflection points correspond to reflectors modelled with the level-set function set to zero [ ϕ (x, z) = 0]. The proposed method aims for stable Vs inversion in a seismic acquisition setting using multicomponent receivers. Synthetic models simulating true Vs , Vp and the location of the geological reflector are used in the study. The inversion starts by locating a flat reflector, ϕ (x, z) = 0, which defines the zone Ω1 between the surface and the reflector, where the initial Vs and Vp fields are also set. To calculate the traveltimes of incident PT (P wave that propagates in Ω1 from source to the reflector), converted PS and reflected PP waves, for both observed and modelled data (forward problem), the methodology proposed by Rawlinson and Sambridge is adopted. This method uses the arrival times of the P waves, Tpt , from the seismic source at each reflector point as secondary sources generating the times Tps and Tpp. These times are calculated as a solution to the eikonal equation by using the Fast Marching method. The PS and PP residual times are minimized by updating Vs , Vp and ϕ (x, z) = 0 through adjoint variables designed from a formulation using Lagrange Multipliers in a variational context. The performance of the algorithm is evaluated for models with synclinal, sinusoidal and monoclinal reflector geometries using numerical tests considering the inversion of: (1) ϕ , given the true values of Vs and Vp ; (2) ϕ and Vs , given the true value of Vp ; (3) ϕ and Vp , given the true value of Vs and (4) the three parameters ϕ, Vs and Vp , simultaneously. Good results are obtained by inverting Vs and ϕ , given the true value of Vp. The simultaneous inversion of the three parameters exhibits promising results, despite the illumination problems caused by the different distribution of the PS, PP and PT time gradients due to the geometry of the reflectors and the acquisition setting (sources–receivers in the same plane). The proposed tomography estimates Vs and reflector positions which could help in statics corrections and improve the lithological characterization of near surface. [ABSTRACT FROM AUTHOR]

Constitutive laws are a necessary ingredient in calculations of glacial isostatic adjustment (GIA) or other surface loading problems (e.g. loading by ocean tides). An idealized constitutive law governed by the Maxwell viscoelastic model is widely used but increasing attention is being directed towards more intricate constitutive laws that, in particular, include transient rheology. In this context, transient rheology collectively refers to dissipative mechanisms activated in addition to creep modelled by the Maxwell viscoelastic model. Consideration of such viscoelastic models in GIA is in its infancy and to encourage their wider use, we present constitutive laws for several experimentally derived transient rheologies and outline a flexible method in which to incorporate them into geophysical problems, such as the viscoelastic deformation of the Earth induced by surface loading. To further motivate this need, we demonstrate, via the Love number collocation method, how predictions of crustal displacement depart significantly between Earth models that adopt only Maxwell viscoelasticity and those with transient rheology. Throughout this paper, we highlight the differences in terminology and emphases between the rock mechanics, seismology and GIA communities, which have perhaps contributed towards the relative scarcity in integrating this broader and more realistic class of constitutive laws within GIA. We focus on transient rheology since the associated deformation has been demonstrated to operate on timescales that range from hours to decades. With ice mass loss enhanced at similar timescales as a consequence of anthropogenically caused climate change, the ability to model GIA with more accurate constitutive laws is an important tool to investigate such problems. [ABSTRACT FROM AUTHOR]

Ambient 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 d of DAS data recorded on a pre-existing fibre optic cable (dark fibres), 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 to retrieve 1-D velocity models for two exemplary fibre subsegments, and compare the results of the conventional and modified workflows. The resulting 1-D 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. [ABSTRACT FROM AUTHOR]

Despite the increasing accuracies of GRACE (Gravity Recovery and Climate Experiment)/GRACE-FO (GRACE Follow-On) gravity field models through worldwide endeavours, 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 modelling 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 (Laser Ranging Interferometers) range-rate residuals as the accuracy of the LRI range-rate data, despite its magnitude of being less than 0.1 mm 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 to 1 hr 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 data sets. 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 equivalent water height derived by using P4M6 filtering show that the maximum root mean square value of spatial differences caused by improving the temporal resolution of the atmospheric de-aliasing models can reach ∼13.4 mm in the subregion of the Congo River Basin. However, using different atmospheric data sets can lead to a maximum difference of ∼16.5 mm. For the Amazon River Basin, the corresponding maximum discrepancy is ∼18.1 mm, and that caused by improving temporal resolution is ∼9.4 mm. We further divide the Congo River Basin into several subregions using a lat-lon regular grid with a spatial resolution of 3°. The subsequent time-series results of mass changes reveal that the maximum contribution of temporal resolution and changes in the atmospheric data sets can reach 11.09 and 21.24 per cent, respectively. This suggests that it is necessary to consider the temporal resolution of de-aliasing products when studying mass changes at a regional scale. [ABSTRACT FROM AUTHOR]

We 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 behaviour 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 multiphase 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. [ABSTRACT FROM AUTHOR]

The Grüneisen parameter is an important parameter for the thermal state and evolution of the core, but its uncertainties and their implications are sometimes overlooked. Several formalisms using different parameters values have been used in different studies, making comparison between studies difficult. In this paper, we use previously published data sets to test the sensitivity of modelling the thermal state of the early core to the different formalisms and parameter values used to describe the evolution of the Grüneisen parameter with density. The temperature of the core obtained in our models is less sensitive to the uncertainties of the parameters used in Al'Tshuler et al. formalism than the uncertainties of the parameters used in Anderson formalism. [ABSTRACT FROM AUTHOR]

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

The effects of pore fluids on the elastic properties of sedimentary rocks can be broadly categorized into two mechanisms: variations in pore compressibility due to the physical properties of the fluids and alterations in the stiffness of the rock frame resulting from rock–fluid interactions. Particularly, as rock–fluid interactions alter the stiffness of contacts between mineral grains, changes in the fluid properties around grain contacts can induce volumetric deformation of the rock as well as in variations of the associated elastic coefficients. Even though many previous studies have explored the influence of swelling of clay minerals, the relationships between changes in elastic moduli and deformation in rocks, which hardly contain swellable minerals, remain to date enigmatic. In this paper, to evaluate quantitatively these effects, drying rates, strains and ultrasonic velocities of small cylindrical Berea sandstone samples were measured as their water saturation was decreased by evaporative drying. The measurements clearly showed drying shrinkage and drastic increases in shear and compressional moduli for all the samples under an almost fully dried condition. Previous studies considered that the alteration in surface energy of grain minerals between wet and dry states affected the contact stiffness between them. Some of them used micromechanical model, uniting the Digby grain contact model with the effective medium theory, to interpret the changes in elastic moduli observed in the non-swellable sandstones during water adsorption on the grain surface. The conventional micromechanical model assumes that a mineral grain is a pure sphere and that the number of contacts between the grains is one. However, grains in a sedimentary rock are generally not purely spherical, and the contact surface is composed of several adhesive asperities. We therefore modified the conventional model by introducing the curvature radius and the number of asperities per contact surface. The modified model well reproduced the shear moduli under wet conditions using the strain and moduli measured under dry conditions. On the other hand, the predictions of the compressional moduli using the model were partially in agreement with the experimental results. Therefore, we attempted to qualitatively interpret the relationship by combining the model with the viscoelastic effect associated with wave-induced fluid flow. The deformation and changes in the elastic moduli of rocks resulting from multiple pore fluids within them can be better understood by the present combined model. [ABSTRACT FROM AUTHOR]

The M w 6.8 Murghob earthquake is the third earthquake in an M w 6.4+ sequence occurring in the Pamir initiated by the 2015 Sarez M w 7.2 earthquake. It is of great significance to investigate their interactions and to assess future seismic hazards in the region. In this paper, we use Sentinel-1 radar interferometric data to retrieve coseismic deformation, invert for the slip distributions of the four events, and then investigate their interactions. The cumulative Coulomb failure stress changes (ΔCFS) suggest that the 2023 Murghob earthquake was promoted by the three prior earthquakes in the sequence. Pre-stress from historical earthquakes is a key factor in explaining the triggering mechanism of the two 2016 M w 6.4+ earthquakes. Stress loading and unloading effects on major faults in the region indicate that future attention should be paid in (1) the segment of the Sarez-Karakul fault north of the Kokuibel Valley, (2) the segment of the Sarez-Murghab thrust fault west of the Sarez-Karakul fault and (3) the east segments of the Pamir thrust fault system, all with a large positive ΔCFS. [ABSTRACT FROM AUTHOR]

It is a well-accepted hypothesis that deep-mantle primary plumes originate from a buoyant source layer at the core–mantle boundary (CMB), where Rayleigh–Taylor instabilities (RTIs) play a key role in the plume initiation process. Previous studies have characterized their growth rates mainly in terms of the density, viscosity and layer-thickness ratios between the denser overburden and the source layer. The RTIs, however, develop in the presence of global flows in the overlying mantle, which can act as an additional factor in the plume mechanics. Combining 2-D computational fluid dynamic (CFD) model simulations and a linear stability analysis, this paper explores the influence of a horizontal global mantle flow in the instability dynamics. Both the CFD simulation results and analytical solutions reveal that the global flow is a dampening factor in reducing the instability growth rate. At a threshold value of the normalized global flow velocity, short- as well as long-wavelength instabilities are completely suppressed, allowing the entire system to advect in the horizontal direction. Using a series of real-scale numerical simulations, this paper also investigates the growth rate  as a function of the density contrast, expressed in Atwood number |${A}_T = ({{{\rho }_1 - {\rho }_2}})/({{{\rho }_1 + {\rho }_2}})$|⁠ , and the viscosity ratio |$\ {\mu }^* = \ {\mu }_1/{\mu }_2$|⁠ , where |${\rho }_1,\ {\mu }_{1\ }$| and |${\rho }_{2,}\ {\mu }_{2\ }$| are densities and viscosities of the overburden mantle and source layer, respectively. It is found that increase in either |${A}_T$| or |${\mu }^*$| promotes the growth rate of a plume. In addition, the stability analysis predicts a nonlinearly increasing RTI wavelength with increasing global flow velocity, implying that the resulting plumes widen their spacing preferentially in the flow direction of kinematically active mantle regions. The theory accounts for additional physical parameters: source-layer viscosity and thickness in the analysis of the dominant wavelengths and their corresponding growth rates. The paper finally discusses the problem of unusually large inter-hotspot spacing, providing a new conceptual framework for the origin of sporadically distributed hotspots of deep-mantle sources. [ABSTRACT FROM AUTHOR]

This paper proposes a procedure of forward modelling 3-D frequency-domain wire-source electromagnetic data using the meshless generalized finite-difference (MGFD) method. This method is based on Taylor series expansions and the weighted least-squares method, and its basic principle is to express the partial derivatives of the unknown function on a particular central point by a linear combination of function values on the adjacent points. The advantages of the method over mainstream forward-modelling methods, for example, the regular finite-difference (FD) method, or the finite-element (FE) method, is that mesh generation is not needed: a discretization in the form of just points is applied instead. This allows the points to be distributed freely to fit the arbitrary shape of the structures in the model, which is helpful in the modelling of complex earth structures. It makes the MGFD method more suitable to deal with complex model than FD method. Also, unlike that in the FE method, interpolation functions are not required and no integral needs to be calculated in MGFD method. This results in high computational efficiency and a concise forward-modelling process. In this paper, the particulars of the MGFD method are introduced, the discretized MGFD system of equations (for an |${\boldsymbol{A}} - {\rm{\ }}\varphi $| potential decomposition of the fields, with the Coulomb gauge condition enforced and a primary–secondary separation approach to deal with the singularity of the source) are solved using a direct solver, and the forward-modelling code are programmed. To test the method and code, we compare the MGFD solutions for three 3-D earth models with the equivalent solutions calculated by other methods, and verify the correctness of the MGFD solution by the good agreement between the corresponding results (with relative error of the electric field |${{\boldsymbol{E}}}_{\boldsymbol{x}}$| smaller than 4.89 per cent). We also investigate the performance of this method when applying different discretizations of points, and when using different weighting functions, to assess the influence of these two factors on the forward-modelling accuracy and efficiency. Results indicate that denser point distributions and straightforward weighting functions result in better accuracy and efficiency. [ABSTRACT FROM AUTHOR]

Plenty of studies have suggested that pore fluid may play an important role in earthquake rupture processes. Establishing numerical models can provide great insight into how pore fluid may affect earthquake rupture processes. However, numerical simulation of 3-D spontaneous ruptures in poroelastic mediums is still a challenging task. In this paper, it is found that a closed-form time-domain Green's function of Biot's poroelastodynamic model can be constructed when the source frequency and source-field distance are within a certain range. The time-domain Green's function is validated by being transformed into the frequency domain and comparing with the frequency-domain Green's functions obtained by former papers. Poroelastic wave propagation phase diagrams for various two-phase poroelastic mediums are then plotted to show the applicable range of frequency and source-field distance for the new time-domain Green's function. It is shown that the applicable range not only include the frequency and spatial range of concern in seismology but also overlap that in acoustics. Based on the time-domain Green's function, the boundary integral equations (BIEs) for modelling dynamic ruptures in elastic mediums are extended to fluid-saturated mediums. In the meantime, a functional relationship between the effective stress tensor and the total stress tensor in fluid-saturated mediums is also obtained, which allows us to directly obtain the effective stress by BIEs. The spontaneous rupture processes controlled by the slip-weakening friction law on faults in elastic mediums and in fluid-saturated mediums are compared. It is found that under the same conditions, fluid-saturated rocks are more prone to supershear rupture than dry rocks. This result suggests that pore fluid may promote the excitation of supershear rupture. The poroelastic wave propagation phase diagrams also suggest that simulating a coseismic phase in the real scale requires a certain sample length in laboratories. They also suggest that an undrained governing equation is suitable for seismic wave propagation simulation in poroelastic media. [ABSTRACT FROM AUTHOR]

The aim of this paper is gravitational inversion, that is, the determination of the mass density distribution of a central body when given the observed external gravitational field that is produced by the mass itself. This inversion for a 3-D mass distribution is known to be grossly non-unique; whereas that for a 2-D mass distribution on a spherical surface is known to be unique. The latter justifies the surface 'mascon', or equivalent-water-thickness, solutions when not considering any interior mass transports for the Earth's time-variable gravity as observed by the GRACE satellite. In this paper, using the gravitational multipole formalism cast in the framework of linear Hilbert space with the notion of inner product, we do two things further: (i) we prove mathematically that the 2-D gravitational inversion on an arbitrary surface is unique; and (ii) we derive the algorithm that leads to the unique exact 2-D mass distribution solution, up to a truncated finite spherical harmonic degree. This pertains directly to reaching refined GRACE mascon solutions that account for the actual Earth surface shape including the ellipsoidal figure and the topography. In the process, we also (re)formulate the exact equations, for spherical or arbitrary shape, unifying the formulae that have appeared in the GRACE-relevant literature. [ABSTRACT FROM AUTHOR]

High-resolution 3-D seismic simulation imposes severe demands for computational memory, making low-storage seismic simulation particularly important. Due to its high-efficiency and low-storage, the half-precision floating-point 16-bit format (FP16) is widely used in heterogeneous computing platforms, such as Sunway series supercomputers and graphics processing unit (GPU) computing platforms. Furthermore, the low-storage Runge–Kutta (LSRK) technique requires lower memory resources compared with the classical Runge–Kutta. Therefore, FP16 and LSRK provide the possibility for low-storage seismic simulation. However, the orders of magnitude of the physical quantities (velocity, stress and Lamé constants) in the elastic wave equations are influenced by the P -wave and S -wave velocities and the densities of the elastic media. This results in a huge order of magnitude difference between the stored velocity and stress values, which exceed the range of the stored values of FP16. In this paper, we introduce three dimensionless constants, Cv , Cs and Cp , into elastic wave equations, and new elastic wave equations are derived. The three constants, Cv , Cs and Cp , keep the orders of magnitude of the velocity and stress at a similar level in the new elastic wave equations. Thus, the stored values of these variables in new equations remain within the range of the stored values of FP16. In addition, we introduce the use of the LSRK due to its low-storage characteristic. In this paper, based on the FP16 and LSRK low-storage techniques, we develop 3 optimized multi-GPU solvers for seismic simulation using the curved grid finite-difference method (CGFDM). Moreover, we perform a series of seismic simulations to verify the accuracy, stability, and validity of the optimized solver coupled with the two techniques. The verifications indicate that through maintaining the calculation accuracy, the computational efficiency of the solver is significantly optimized, and the memory usage is remarkably reduced. In particular, under the best conditions, the memory usage can be reduced to nearly 1/3 that of the original CGFDM solver. [ABSTRACT FROM AUTHOR]

SUMMARY: In this paper, the 3-D theoretical seismogram at the free surface of an elastic half-space with 'vertical transverse isotropy (VTI)' and an 'elliptic anisotropy' is synthesized via semi-analytical formulation. To this end, the time-domain solution is obtained via the Hankel and Laplace integral transforms and the Cagniard contour integration scheme. The formulations include complete Green's functions expressed in compact and elegant formulations in terms of some elementary line integrals over the finite interval (0, π/2), due to the general point force and point double-couple source of arbitrary orientation, varying with time as Heaviside step function. This solution is further implemented to model seismic waves generated from moment tensor source with a ramp-type slip function. The solution clearly delineates compressional P waves, shear SV waves and SH waves, diffracted SPwaves and surface Rayleigh waves, with two critical distances in which the mode conversion happens. The first marks the conversion of the total reflection of SV wave into the diffracted SPwave travelling on the free surface with the velocity of compressional wave. The second, however, marks a location where the order of arrival times of SP and SH waves is reversed. The interesting phenomenon of shear wave splitting by the apparent time lag between the arrivals of SV and SHwaves is further demonstrated. Studied also is the effect of non-double-couple (non-DC) components of the moment tensor of the shear fault. Particularly, it is shown that non-DC components may lead to large amplitude of the P wave in the presence of anisotropy, resulting in the changed polarity of P wave at the onset of the motion. It is found that for the elliptic anisotropy, and for oblique-thrust faults, the isotropic moment (ISO) is always negative if the elastic constants satisfy C11 > C33 and C44 > C66, whereas it is always positive if C11 < C33 and C44 < C66 (for oblique-reverse faults, the situation is reversed). It is shown that for some sedimentary rocks, the decomposition of moment tensor in the general VTI and in its elliptic approximation are close to each other, and that the moment tensor decomposition in VTI media is very sensitive to the elastic coefficients, such that a small change in the elastic constants may lead to a remarkable change in the moment tensor decomposition. Particularly, the effect of elastic constant C13 is shown to be significant. While the time-domain solution obtained in this paper can be also degenerated to the solution in the corresponding isotropic half-space, the Mathematica code of the solution provided can be served as benchmark to other numerical solutions, applicable for the computation of the theoretical seismogram in the involved media. [ABSTRACT FROM AUTHOR]

This paper presents a new methodology to assist geophysicists in determining the first-break event in a 3-D seismic data set 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 were 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 3-D 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 with long-time averages. Finally, we performed the static correction calculation to ensure that the time distortion resulting from the terrain and the low-velocity layer was mitigated in shoot gathers and in the stacked section. [ABSTRACT FROM AUTHOR]

Correcting atmospheric effects in measurements of vertical acceleration is challenging, particularly at frequencies above 0.3 mHz. Corrections based on 4-D atmosphere models should be the most effective, but are limited in temporal and spatial resolutions. So these models are commonly combined with a constant admittance between air pressure and gravity, though the correction can be further improved by allowing a frequency-dependent admittance. This paper studies the effectiveness of such procedures given that air pressure admittance varies with time because of variability in the underlying local atmospheric process, and finds that admittances estimated by cross-spectral analysis gives better results than a constant admittance does. The best results are achieved by combining a 4-D model with admittances estimated from time-series that cover the event of interest. [ABSTRACT FROM AUTHOR]

Dense array observation and seismic interferometry have revolutionized the imaging schemes of the earth structure. It is becoming possible to directly obtain the lateral variation of the earth's structure by applying array-based methods such as the cross-correlation beamforming (CBF) of the ambient noise to the subsets of the dense array, without tomography. CBF has been proven to extract the azimuth-averaged multimode surface wave dispersion curves. However, the resolution of the dispersion image generated by conventional CBF is low at high frequencies in the frequency–velocity (f - v) domain. Moreover, the irregular array geometry and uneven source distribution would bias the result of CBF, especially for the estimation of azimuth-dependent velocity. In this paper, two beamforming (BF) es are suggested to improve the resolution of multimode dispersion images in the f - v domain. First, the geometrical spreading of the wavefield is corrected to enhance the amplitude at high frequency (or large distance) and thereby improve the resolution of the dispersion image at high frequency. We call this scheme weighted correlation beamforming (WCBF). The azimuth-averaged velocity can be estimated with sufficient resolution using WCBF by stacking the BF output at each azimuth. We show that WCBF is the 2-D Fourier transform of the spatial wavefield from the viewpoint of the wavefield transform. Secondly, a modified beamforming scheme (MCBF) is suggested to reduce the effect of uneven source and/or irregular array geometry. The delay and summation in MCBF are performed only for plane waves incident from the stationary phase region. The azimuth-dependent velocity can therefore be estimated by MCBF with less dependence on the array geometry, as well as on the uneven source distribution. In terms of the estimation of azimuth-averaged phase velocity, we show the F-J method, another array-based method for extracting multimode surface waves from ambient noise using the Fourier–Bessel transform, is the azimuth-averaged version of WCBF. The reliability of WCBF and MCBF is verified based on the synthetic and field data using the array with different geometry. The dispersion image of multimode Rayleigh wave phase velocity at local and regional scales can be generated by WCBF or MCBF with high resolution. In particular, multimode dispersion curves at the local scale can be measured by MCBF with sufficient accuracy using quite short recordings from hours to days. This offers the possibility of a rapid assessment of the media properties. [ABSTRACT FROM AUTHOR]

This paper explores the applicability of ensemble Kalman inversion (EKI) with level-set parametrization for solving geophysical inverse problems. In particular, we focus on its extension to induced polarization (IP) data with uncertainty quantification. IP data may provide rich information on characteristics of geological materials due to its sensitivity to characteristics of the pore–grain interface. In many IP studies, different geological units are juxtaposed and the goal is to delineate these units and obtain estimates of unit properties with uncertainty bounds. Conventional inversion of IP data does not resolve well sharp interfaces and tends to reduce and smooth resistivity variations, while not readily providing uncertainty estimates. Recently, it has been shown for DC resistivity that EKI is an efficient solver for inverse problems which provides uncertainty quantification, and its combination with level set parametrization can delineate arbitrary interfaces well. In this contribution, we demonstrate the extension of EKI to IP data using a sequential approach, where the mean field obtained from DC resistivity inversion is used as input for a separate phase angle inversion. We illustrate our workflow using a series of synthetic and field examples. Variations with uncertainty bounds in both DC resistivity and phase angles are recovered by EKI, which provides useful information for hydrogeological site characterization. Although phase angles are less well-resolved than DC resistivity, partly due to their smaller range and higher percentage data errors, it complements DC resistivity for site characterization. Overall, EKI with level set parametrization provides a practical approach forward for efficient hydrogeophysical imaging under uncertainty. [ABSTRACT FROM AUTHOR]

Clairaut's theory that relates the Earth's oblate figure and internal ellipticity to its gravity under rotational-hydrostatic equilibrium has reigned classical geodesy over the centuries. In this paper, we (i) derive from first principles the classical Clairaut's theory for the polar oblateness of a rotating planet under axi-symmetric rotational-hydrostatic equilibrium and (ii) extend the development to the triaxial case for the equatorial ellipticity of a tidally locked synchronous-rotating moon under rotational-tidal-hydrostatic equilibrium. Typical derivations of the classical Clairaut's theory presented in the literature being rather laborious even to first order, we instead exploit two concise forms of methodology: the gravitational multipole formalism on the physics side, and the Jacobian determinant for the Clairaut coordinate transformation on the mathematics side. The outcome is a logical and straightforward derivation of Clairaut's theory to first order in its entirety, encompassing all the equations and related formulas in geodesy bearing Clairaut's name. That further allows a natural extension to a tidally locked moon. In particular it is demonstrated that the same Clairaut's differential equation applies to both cases governing both the polar oblateness and the equatorial ellipticity. [ABSTRACT FROM AUTHOR]

In this paper, we propose a fully coupled two-phase poroelastic deformation theory for a spherically layered and self-gravitating Earth. The earth model consists of a solid inner core, a fluid outer core and poroelastic mantle and crust. The boundary-value problems are posed in the Laplace-transformed domain using the spherical system of vector functions, and analytical solutions are obtained in each layer using the dual-variable and position matrix method. As an application, the surface loading problem is considered. The undrained and drained limits are discussed with the corresponding governing equations, which are different from the conventional ones where body force is ignored. Numerical examples show that the poroelastic effect can be significant since the difference between the results for undrained and drained limits are large with obvious time-dependent deformation. This newly derived theory for the coupled boundary-value problem will have broad applications, such as displacement and gravity change due to groundwater depletion, poroelastic rebound after an earthquake and risk evaluation of earthquakes induced by water injection. [ABSTRACT FROM AUTHOR]

This paper proposes the use of the parametric bootstrap as a simple way of obtaining reliable confidence intervals for the parameters in the space–time epidemic-type aftershock sequence (ETAS) model. Using an earthquake data catalogue spanning almost 20 yr of seismic activity in the Pacific Ocean off the shores of Vancouver Island (British Columbia, Canada), the authors show through simulation that the confidence intervals based on asymptotic maximum-likelihood theory can sometimes be misleading. In contrast, confidence intervals based on empirical quantiles from bootstrap samples have reliable nominal coverage, provided that edge effects are properly taken into account, both in the bootstrap procedure and in simulation studies. An R package called ETASbootstrap , developed by the authors, facilitates the use of this resampling procedure. [ABSTRACT FROM AUTHOR]

In geophysical surveys or experiments, recorded data are used to constrain properties of the planetary subsurface, oceans, atmosphere or cryosphere. How the experimental data are collected significantly influences which parameters can be resolved and how much confidence can be placed in the results. Bayesian experimental design methods characterize, quantify and maximize expected information post-experiment—an optimization problem. Typical design parameters that can be optimized are source and/or sensor types and locations, and the choice of modelling or data processing methods to be applied to the data. These may all be optimized subject to various physical and cost constraints. This paper introduces variational design methods, and discusses their benefits and limitations in the context of geophysical applications. Variational methods have recently come to prominence due to their importance in machine-learning applications. They can be used to design experiments that best resolve either all model parameters, or the answer to specific questions about the system to be interrogated. The methods are tested in three schematic geophysical applications: (i) estimating a source location given arrival times of radiating energy at sensor locations, (ii) estimating the contrast in seismic velocity across a stratal interface given measurements of the amplitudes of seismic wavefield reflections from that interface, and (iii) designing a survey to best constrain CO2 saturation in a subsurface storage scenario. Variational methods allow the value of expected information from an experiment to be calculated and optimized simultaneously, which results in substantial savings in computational cost. In the context of designing a survey to best constrain CO2 saturation in a subsurface storage scenario, we show that optimal designs may change substantially depending on the particular questions of interest. We also show that one method, so-called DN design, can be effective at substantially lower computational cost than other methods. Overall, this work demonstrates that optimal design methods could be used more widely in Geophysics, as they are in other scientifically advanced fields. [ABSTRACT FROM AUTHOR]

In view of the problem of identifying palaeochannels with high concealment and complex geological structures, this paper proposes an intelligent recognition method for palaeochannels based on Frangi filtering and deep learning. The methodology makes use of Maximum Entropy Wigner–Ville Distribution (MEWVD) method to process the original instantaneous amplitude data, which enhance the distinct features of micro palaeochannels in different sensitive frequency bands. The sample 2-D stratigraphic images generated from these data is labelled for the pre-training process of Attention R2U-Net neural network model. Subsequently, Frangi filter is used to identify and enhance the linear structures of river channels in target stratigraphic images, improving the identification effect of the neural network. Finally, RGB data fusion and 3-D visualization carving are performed on the identification data. This method not only eliminates redundant information using the Frangi filter but also proves that the Attention R2U-Net network model structure with attention mechanism and residual convolution structure can effectively improve the segmentation effect for river channels at different scales. Experimental examples show that this method achieves pixel-level feature segmentation of 3-D seismic data for palaeochannels, accurately depicting their shape, width, thickness, flow direction and other features, thus providing support for subsequent well deployment and horizontal well fracturing selection. [ABSTRACT FROM AUTHOR]

Physics-based and statistic-based models for describing seismic occurrence are two sides of the same coin. In this paper, we compare the temporal organization of events obtained in a spring-block model for the seismic fault with the one predicted by probabilistic models for seismic occurrence. Thanks to the optimization of the parameters, by means of a Maximum Likelihood Estimation, it is possible to identify the statistical model which fits better the physical one. The results show that the best statistical model must take into account the non-trivial interplay between temporal clustering, related to aftershock occurrence, and the stress discharge following the occurrence of high magnitude main shocks. The two mechanisms contribute in different ways according to the minimum magnitude considered in the data fitting catalogue. [ABSTRACT FROM AUTHOR]

The joint ESA/NASA Mass-change And Geosciences International Constellation (MAGIC) has the objective to extend time-series from previous gravity missions, including an improvement of accuracy and spatio-temporal resolution. The long-term monitoring of Earth's gravity field carries information on mass change induced by water cycle, climate change and mass transport processes between atmosphere, cryosphere, oceans and solid Earth. MAGIC will be composed of two satellite pairs flying in different orbit planes. The NASA/DLR-led first pair (P1) is expected to be in a near-polar orbit around 500 km of altitude; while the second ESA-led pair (P2) is expected to be in an inclined orbit of 65°–70° at approximately 400 km altitude. The ESA-led pair P2 Next Generation Gravity Mission shall be launched after P1 in a staggered manner to form the MAGIC constellation. The addition of an inclined pair shall lead to reduction of temporal aliasing effects and consequently of reliance on de-aliasing models and post-processing. The main novelty of the MAGIC constellation is the delivery of mass-change products at higher spatial resolution, temporal (i.e. subweekly) resolution, shorter latency and higher accuracy than the Gravity Recovery and Climate Experiment (GRACE) and Gravity Recovery and Climate Experiment Follow-On (GRACE-FO). This will pave the way to new science applications and operational services. In this paper, an overview of various fields of science and service applications for hydrology, cryosphere, oceanography, solid Earth, climate change and geodesy is provided. These thematic fields and newly enabled applications and services were analysed in the frame of the initial ESA Science Support activities for MAGIC. The analyses of MAGIC scenarios for different application areas in the field of geosciences confirmed that the double-pair configuration will significantly enlarge the number of observable mass-change phenomena by resolving smaller spatial scales with an uncertainty that satisfies evolved user requirements expressed by international bodies such as IUGG. The required uncertainty levels of dedicated thematic fields met by MAGIC unfiltered Level-2 products will benefit hydrological applications by recovering more than 90 per cent of the major river basins worldwide at 260 km spatial resolution, cryosphere applications by enabling mass change signal separation in the interior of Greenland from those in the coastal zones and by resolving small-scale mass variability in challenging regions such as the Antarctic Peninsula, oceanography applications by monitoring meridional overturning circulation changes on timescales of years and decades, climate applications by detecting amplitude and phase changes of Terrestrial Water Storage after 30 yr in 64 and 56 per cent of the global land areas and solid Earth applications by lowering the Earthquake detection threshold from magnitude 8.8 to magnitude 7.4 with spatial resolution increased to 333 km. [ABSTRACT FROM AUTHOR]

Satellite altimetry provides major data sources for marine gravity recovery, and typical altimetry derived models, for example, DTU21 and SS V32.1, were usually released with 1 arcmin × 1 arcmin gridding interval. Their true resolution is much lower than the nominal ∼2 km level. By contrast, the in situ ship-borne measurements are considered to have better short-wavelength resolution. In this paper, we aim to propose a new method to fully utilize satellite altimetry data and ship-borne measurements, namely the frequency-domain fusion method, and give certain analysis of new method along with two spatial-domain fusion methods. Comprehensive analysis is focus on four aspects: gravity signals in fusion images, numerical verifications, power spectra, as well as coherence analysis. Initial evaluation indicates that, first, the frequency-domain fusion method has advantage in flexibility, since it can autonomously select dominant bands to fuse different data sets. Secondly, the new method retains medium-long wavelength signals from altimetry-derived model and effectively incorporate dominant short-wavelength signals of in-situ measurements, while the spatial-domain methods are essentially full-wavelength fusion and inevitably diminish the role of satellite altimetry. To some extent, the new method maximize the positive contribution of satellite altimetry measurements and efficiently exploit the benefits of ship-borne data. Finally, verification experiments were similarly designed in three regions with different amount and ratio of ship-borne data to thoroughly evaluate various methods. In two regions with gridded and dense along-cruise ship-borne data, the average accuracy of this frequency-domain fusion results is improved by 0.346 and 0.613 mGal, respectively. In a region with sparse ship-borne data, we still recommend using spatial-domain fusion methods since the new method is unable to align ship-borne data with model grid. It is concluded from the above analysis that the new method effectively incorporates the short-wavelength signals from ship-borne data into the altimetry-derived gravity field model, and it is significant that the new method maximizes the application of advantageous bands from different data sources. [ABSTRACT FROM AUTHOR]

More than 20 yr of measurement data of the gravity missions GRACE (Gravity Recovery And Climate Experiment) and GRACE-FO (GRACE-Follow-On) allow detailed investigations of long-term trends in continental terrestrial water storage (TWS). However, the spatial resolution of conventional GRACE/GRACE-FO data products is limited to a few hundred kilometres which restrains from investigating hydrological trends at smaller spatial scales. In this study GRACE and GRACE-FO data have been used to calculate TWS trends with maximized spatial resolution. Conventionally, GRACE/GRACE-FO is presented as a series of either unconstrained gravity fields post-processed with spatial low pass filters or constrained inversions commonly known as Mascon products. This paper demonstrates that both approaches to suppress spatially correlated noise are mathematically equivalent. Moreover, we demonstrate that readily inverting all available sensor data from GRACE/GRACE-FO for a single TWS trend map, together with annual variations and a mean gravity field, provides additional spatial detail not accessible from the standard products. The variable trade-off between spatial and temporal resolution as a unique feature of satellite gravimetry allows for gravity products that are tailored towards specific geophysical applications. We show additional signal content in terms of long-term water storage trends for four dedicated examples (Lake Victoria, Northwest India, Bugachany Reservoir and High Plains Aquifer) for which external information from other remote sensing instruments corroborates the enhanced spatial resolution of the new mean-field trend product. [ABSTRACT FROM AUTHOR]

This paper presents a weakly supervised machine learning method, which we call DAS-N2N, for suppressing strong random noise in distributed acoustic sensing (DAS) recordings. DAS-N2N requires no manually produced labels (i.e. pre-determined examples of clean event signals or sections of noise) for training and aims to map random noise processes to a chosen summary statistic, such as the distribution mean, median or mode, whilst retaining the true underlying signal. This is achieved by splicing (joining together) two fibres hosted within a single optical cable, recording two noisy copies of the same underlying signal corrupted by different independent realizations of random observational noise. A deep learning model can then be trained using only these two noisy copies of the data to produce a near fully denoised copy. Once the model is trained, only noisy data from a single fibre is required. Using a data set from a DAS array deployed on the surface of the Rutford Ice Stream in Antarctica, we demonstrate that DAS-N2N greatly suppresses incoherent noise and enhances the signal-to-noise ratios (SNR) of natural microseismic icequake events. We further show that this approach is inherently more efficient and effective than standard stop/pass band and white noise (e.g. Wiener) filtering routines, as well as a comparable self-supervised learning method based on masking individual DAS channels. Our preferred model for this task is lightweight, processing 30 s of data recorded at a sampling frequency of 1000 Hz over 985 channels (approximately 1 km of fibre) in <1 s. Due to the high noise levels in DAS recordings, efficient data-driven denoising methods, such as DAS-N2N, will prove essential to time-critical DAS earthquake detection, particularly in the case of microseismic monitoring. [ABSTRACT FROM AUTHOR]

The deployment of electromagnetic (EM) induction tools while drilling is one of the standard routines for assisting the geosteering decision-making process. The conductivity distribution obtained through the inversion of the EM induction log can provide important information about the geological structure around the borehole. To image the 3-D geological structure in the subsurface, 3-D inversion of the EM induction log is required. Because the inversion process is mainly dependent on forward modelling, the use of a fast and accurate forward modelling tool is essential. In this paper, we present an improved version of the integral equation (IE) based modelling technique for general anisotropic media with domain decomposition pre-conditioning. The discretized IE after domain decomposition equals a fixed-point equation that is solved iteratively with either the block Gauss–Seidel or Jacobi pre-conditioning. Within each iteration, the inverse of the block matrix is computed using a Krylov subspace method instead of a direct solver. An additional reduction in computational time is obtained by using an adaptive relative residual stopping criterion in the iterative solver. Using this domain decomposition scheme, numerical experiments show computation time reductions by factors of 1.97–2.84 compared to solving the full-domain IE with a GMRES solver and a contraction IE pre-conditioner. Additionally, the reduction of memory requirement for covering a large area of the induction tool sensitivity enables acceleration with limited GPU memory. Hence, we conclude that the domain decomposition method is improving the efficiency of the IE method by reducing the computation time and memory requirement. [ABSTRACT FROM AUTHOR]

During the accretion stage, large impacts provided sufficient energy to melt the entire mantle into a terrestrial magma ocean. Processes occurring in the magma ocean may have led to the formation of heterogeneities still found in modern ocean island basalts. So far, no definitive mechanism exists to explain the survival of early heterogeneities for approximately 4.5 Ga. Addressing this question requires understanding the efficiency of convective mixing during both the early molten and the solid-state stages experienced by the Earth's mantle. While mixing in the solid mantle and in an essentially crystallized magma ocean has been relatively well documented, the efficiency of convective mixing in a liquid magma ocean has received less attention. In this paper we characterized the mixing efficiency of a convecting fluid in a rotating spherical shell, accounting for inertial effects, by computing finite-time Lyapunov exponents (i.e. the Lagrangian strain rate). We conducted a series of numerical experiments for a regime where the influence of the buoyancy force dominates that of rotation and we derived scaling laws to predict the mixing efficiency. We found that for a terrestrial magma ocean, in its fully liquid state, mixing time is of the order of a few minutes or less, even for initially large (∼1000 km) heterogeneities. Therefore, passive early mantle heterogeneities cannot survive a fully molten magma ocean stage. This suggests that short-lived heterogeneities (e.g. 182Hf−182W) were either created at the end of the accretional stage, or were stored in deeper regions of the Earth. [ABSTRACT FROM AUTHOR]

Polar regions and Greenland in particular are highly sensitive to global warming. Impacts on Greenland's glaciers may be observed through the increasing number of calving events. However, a direct assessment of the calving activity is limited due to the remoteness of polar regions and the cloudy weather which makes impossible a recurrent observation through satellite imagery. To tackle this issue, we exploit the seismological network deployed in Greenland which actively records seismic signals associated with calving events, hereinafter referred to as glacial earthquakes. These seismic signals present a broad frequency range and a wide diversity of waveform which make them difficult to discriminate from tectonic events as well as anthropogenic and natural noises. In this study, we start from two catalogues of known events, one for glacial earthquake events which occurred between 1993 and 2013 and one for earthquakes which occurred in the same time period, and we implement a detection algorithm based on the STA/LTA method to extract signals' events from continuous data. Then, we train and test a machine learning processing chain based on the Random Forest algorithm which allows us to automatically associate the events respectively with calving and tectonic activity, with a certain probability. Finally, we investigate 844 selected days spanning time of continuous data from the Greenland regional seismic network which results in a new, more exhaustive, catalogue of glacial earthquakes expanded of 1633 newly detected glacial events. Moreover, we extensively discuss the choice of the features used to describe glacial earthquakes, in particular the 39 new features created in this study which have drastically improved our results with 7 of the 10 best features being in the added set. The perspective of further expansion of the glacial earthquake catalogue applying the processing chain discussed in this paper on different time spans highlights how combining seismology and machine learning can increase our understanding of the spatio-temporal evolution of calving activity in remote regions. [ABSTRACT FROM AUTHOR]

Concerns raised by Okazaki & Ueda (2022) on the paper by Sambridge et al. (2022) are addressed. Two issues are discussed and some new numerical results presented. The first concerns whether the properties of the Wasserstein time-series misfit introduced in our earlier paper will translate to model space non-uniqueness in a seismic waveform inversion setting. It is argued that this is unlikely, given the special conditions, which must exist between all observed/predicted seismic waveform pairs for non-uniqueness to result. The second issue discussed is the efficacy of using the Sliced Wasserstein algorithm of Bonneel et al. (2015) as an alternate to the marginal Wasserstein algorithm, as proposed by Okazaki & Ueda (2022). It is argued that for optimization-based waveform fitting, the Sliced Wasserstein algorithm is a viable alternate provided care is taken to ensure that conditions arise which do invalidate analytical derivative expressions of the resulting Wasserstein misfit. In practice, this would likely mean recasting the 2D Optimal Transport problem posed in our earlier paper onto unstructured grids. [ABSTRACT FROM AUTHOR]

The analysis of earthquake recordings from three-component instruments can be challenging due to overlapping events. Time–frequency (TF) polarization methods are efficient tools for this purpose, which can discriminate these events. Previous polarization methods did not consider all three components simultaneously while transferring data to TF domain, which can cause inaccuracies in the reconstruction of wave amplitudes. Therefore, the three-component sparse adaptive S transform (3C-SAST) algorithm is preferred to other TF decompositions since it is mainly developed for polarization analysis purposes, and outperforms other TF methods. In this paper, we developed the 3C-SAST by adding a parameter to adjust the sparseness of the solution and make the resolution flexible. The developed TF decomposition is then used to extend the Morozov & Smithson method to TF domain, and devise a new TF polarization filter whose invertibility and resolution flexibility make it a promising tool for wavefield separation. This filter can eliminate the out-of-plane arrival energies and extract the Rayleigh waves for multicomponent data, which has application in Rayleigh wave tomography and seismological studies. We demonstrated the efficiency of the proposed method for seismic surface waves separation using synthetic signals and three-component teleseismic earthquake recording. [ABSTRACT FROM AUTHOR]

Compared to existing schemes, the decomposition of moment tensors for mining-induced seismic events into closing-crack and double-couple components has the advantage that each can be interpreted in terms of a physical source process (excavation convergence and slip/shear, respectively). Obviously, not every moment tensor permits such a decomposition, and we translate existing bounds on those that do to the Hudson source-type plot. For moment tensors falling within these bounds, it has previously been noted that there will be a infinite set of possible decompositions in general. We derive an implicit equation defining this set and suggest physically motivated criteria that can be used to select a single decomposition from it. Furthermore, for moment tensors falling outside the source-type bounds, we present a simple geometric scheme for determining the closest moment tensor within them (which can then be decomposed). To demonstrate the methods developed in this paper, we apply them to two catalogues of mining-induced seismicity. [ABSTRACT FROM AUTHOR]

In this paper, we address the study of the seismoelectric response of an elastic medium in contact with a poroelastic half-space. In particular, we advance in the understanding of the generation mechanism of the interface response (IR) and the evanescent electromagnetic (EM) fields occurring at the contact between both media, by proposing a seismoelectric phenomenological model (SPM). Essentially, the model consists of a sequence of electric dipoles that are activated successively, simulating the seismic-to-EM energy conversion taking place with the arrival of a seismic wave at the interface separating the media. We obtained SPM responses for different scenarios and acquisition configurations and compared them with responses computed using a code based on the finite-elements method, which solves the seismoelectric equations in the compressional P and vertical shear SV waves coupled with the transverse-magnetic (TM) fields (PSVTM) mode. The SPM successfully represents not only the evanescent wave but also the IR within the elastic medium. In particular, we show that the SPM is able to faithfully reproduce the relative amplitudes of both events and their radiation patterns with a minimum computational cost. In this way, it provides a novel insight in the study of the physical phenomenon behind the seismoelectric conversions. [ABSTRACT FROM AUTHOR]