ALBERT

Description: Summary Landslides can cause devastating damage. In particular, heavy rainfall-triggered landslides pose a chain of natural hazards. However, such events are often difficult to detect, leaving the physical processes poorly understood. Here we apply a novel surface-wave detector to detect and locate landslides during the transit of Typhoon Talas 2011. We identify multiple landslides triggered by Typhoon Talas, including a landslide in the Tenryu Ward, Shizuoka prefecture, Japan, ∼400 km east from the typhoon track. The Tenryu landslide displaced a total volume of 1.2 − −1.5 × 106 m. The landslide is much smaller than those detected by using globally recorded surface waves, yet the event generated coherent seismic signals propagating up to 3000 km away. Our observations show that attributes of small and large landslides may follow the same empirical scaling relationships, indicating possible invariant failure mechanisms. Our results also suggest an alerting technology to detect and locate landslides with a sparse seismic network.

Description: Summary We studied the broad-band spectra of the 8 largest earthquakes that have occurred in Chile in the last 25 years using strong-motion records and 1-Hz high-rate GNSS (cGNSS) data. To avoid the numerical instability problem with the double integration of the accelerograms, we computed velocity spectra integrating the acceleration time series in the spectral domain and compared them to time-differentiated the cGNSS displacement records. To compute the velocity spectrum, we used a multitaper algorithm so as to provide stability over the entire spectral band. We found that the velocity spectra of records obtained close to the main rupture of the earthquakes are different from classical Aki and Brune spectra. The velocity spectrum of large events in Chile presents a flat trend at low frequencies produced by the near-field waves. This trend converges at low frequencies to the static displacement as determined from GNSS data. For different magnitude earthquakes, we observe a transition in the ground-velocity spectrum from a decay of ${f^{ - 1}}$ at high frequencies and a flat trend at low frequencies to a more classical model with a peak at the corner frequency. The source-station distance influences the shape of the velocity spectrum at low frequencies, but there is no simple rule for the records available at present. At intermediate frequencies, the spectra are controlled by surface waves and S waves. We found a transition in the velocity spectrum for the 2014 Iquique earthquake, which indicates a change in the decay of the spectrum for stations at distances greater than ∼200 km. Finally, we show that the flat low-frequency trend of the velocity spectra determined from accelerograms, and the peak ground-displacement (PGD) determined from GNSS data scales with the moment to the power 2/3.

Description: Summary Seismic imaging techniques such as elastic full waveform inversion (FWI) have their spatial resolution limited by the maximum frequency present in the observed waveforms. Scales smaller than a fraction of the minimum wavelength cannot be resolved, and only a smoothed, effective version of the true underlying medium can be recovered. These finite-frequency effects are revealed by the upscaling or homogenization theory of wave propagation. Homogenization aims at computing larger scale effective properties of a medium containing small-scale heterogeneities. We study how this theory can be used in the context of FWI. The seismic imaging problem is broken down in a two-stage multiscale approach. In the first step, called homogenized full waveform inversion (HFWI), observed waveforms are inverted for a smooth, fully anisotropic effective medium, that does not contain scales smaller than the shortest wavelength present in the wavefield. The solution being an effective medium, it is difficult to directly interpret it. It requires a second step, called downscaling or inverse homogenization, where the smooth image is used as data, and the goal is to recover small-scale parameters. All the information contained in the observed waveforms is extracted in the HFWI step. The solution of the downscaling step is highly non-unique as many small-scale models may share the same long wavelength effective properties. We therefore rely on the introduction of external a priori information, and cast the problem in a Bayesian formulation. The ensemble of potential fine-scale models sharing the same long wavelength effective properties is explored with a Markov chain Monte Carlo algorithm. We illustrate the method with a synthetic cavity detection problem: we search for the position, size and shape of void inclusions in a homogeneous elastic medium, where the size of cavities is smaller than the resolving length of the seismic data. We illustrate the advantages of introducing the homogenization theory at both stages. In HFWI, homogenization acts as a natural regularization helping convergence toward meaningful solution models. Working with fully anisotropic effective media prevents the leakage of anisotropy induced by the fine scales into isotropic macro-parameters estimates. In the downscaling step, the forward theory is the homogenization itself. It is computationally cheap, allowing us to consider geological models with more complexity (e.g. including discontinuities) and use stochastic inversion techniques.

Description: Summary Particle Image Velocimetry (PIV), a method based on image cross-correlation, is widely used for obtaining velocity fields from time series of images of deforming objects. Rather than instantaneous velocities, we are interested in reconstructing cumulative deformation, and use PIV-derived incremental displacements for this purpose. Our focus is on analogue models of tectonic processes, which can accumulate large deformation. Importantly, PIV provides incremental displacements during analogue model evolution in a spatial reference (Eulerian) frame, without the need for explicit markers in a model. We integrate the displacements in a material reference (Lagrangian) frame, such that displacements can be integrated to track the spatial accumulative deformation field as a function of time. To describe cumulative, finite deformation, various strain tensors have been developed, and we discuss what strain measure best describes large shape changes, as standard infinitesimal strain tensors no longer apply for large deformation. PIV or comparable techniques have become a common method to determine strain in analogue models. However, the qualitative interpretation of observed strain has remained problematic for complex settings. Hence, PIV-derived displacements have not been fully exploited before, as methods to qualitatively characterize cumulative, large strain have been lacking. Notably, in tectonic settings, different types of deformation - extension, shortening, strike-slip - can be superimposed. We demonstrate that when shape changes are described in terms of Hencky strains, a logarithmic strain measure, finite deformation can be qualitatively described based on the relative magnitude of the two principal Hencky strains. Thereby, our method introduces a physically meaningful classification of large 2D strains. We show that our strain type classification method allows for accurate mapping of tectonic structures in analogue models of lithospheric deformation, and complements visual inspection of fault geometries. Our method can easily discern complex strike-slip shear zones, thrust faults and extensional structures and its evolution in time. Our newly developed software to compute deformation is freely available and can be used to post-process incremental displacements from PIV or similar autocorrelation methods.

Description: Summary It is generally accepted that melt extraction from the mantle at mid-ocean ridges is concentrated in narrow regions of elevated melt fraction called channels. Two feedback mechanisms have been proposed to explain why these channels grow by linear instability: shear flow of partially molten mantle and reactive flow of the ascending magma. These two mechanisms have been studied extensively, in isolation from each other, through theory and laboratory experiments as well as field and geophysical observations. Here, we develop a consistent theory that accounts for both proposed mechanisms and allows us to weigh their relative contributions. We show that interaction of the two feedback mechanisms is insignificant and that the total linear growth rate of channels is well-approximated by summing their independent growth rates. Furthermore, we explain how their competition is governed by the orientation of channels with respect to gravity and mantle shear. By itself, analysis of the reaction-infiltration instability predicts the formation of tube-shaped channels. We show that with the addition of even a small amount of extension in the horizontal, the combined instability favours tabular channels, consistent with the observed morphology of dunite bodies in ophiolites. We apply the new theory to mid-ocean ridges by calculating the accumulated growth and rotation of channels along streamlines of the solid flow. We show that reactive flow is the dominant instability mechanism deep beneath the ridge axis, where the most unstable orientation of high-porosity channels is sub-vertical. Channels are then rotated by the solid flow away from the vertical. The contribution of the shear-driven instability is confined to the margins of the melting region. Within the limitations of our study, the shear-driven feedback does not appear to be responsible for significant melt focusing or for the shallowly dipping seismic anisotropy that has been obtained by seismic inversions.

Description: Summary To understand earth processes, geoscientists infer subsurface earth properties such as electromagnetic resistivity or seismic velocity from surface observations of electromagnetic or seismic data. These properties are used to populate an earth model vector, and the spatial variation of properties across this vector sheds light on the underlying earth structure or physical phenomenon of interest, from groundwater aquifers to plate tectonics. However, to infer these properties the spatial characteristics of these properties need to be known in advance. Typically, assumptions are made about the length scales of earth properties, which are encoded a priori in a Bayesian probabilistic setting. In an optimisation setting, appeals are made to promote model simplicity together with constraints which keep models close to a preferred model. All of these approaches are valid, though they can lead to unintended features in the resulting inferred geophysical models owing to inappropriate prior assumptions, constraints or even the nature of the solution basis functions. In this work it will be shown that in order to make accurate inferences about earth properties, inferences can first be made about the underlying length scales of these properties in a very general solution basis. From a mathematical point of view, these spatial characteristics of earth properties can be conveniently thought of as “properties” of the earth properties. Thus, the same machinery used to infer earth properties can be used to infer their length scales. This can be thought of as an “infer to infer” paradigm analogous to the “learning to learn” paradigm which is now commonplace in the machine learning literature. However, it must be noted that (geophysical) inference is not the same as (machine) learning, though there are many common elements which allow for cross-pollination of useful ideas from one field to the other, as is shown here. A non-stationary trans-dimensional Gaussian Process (TDGP) is used to parameterise earth properties, and a multi-channel stationary TDGP is used to parameterise the length scales associated with the earth property in question. Using non-stationary kernels, i.e., kernels with spatially variable length scales, models with sharp discontinuities can be represented within this framework. As GPs are multi-dimensional interpolators, the same theory and computer code can be used to solve geophysical problems in 1D, 2D and 3D. This is demonstrated through a combination of 1D and 2D non-linear regression examples and a controlled source electromagnetic (CSEM) field example. The key difference between this and previous work using TDGP is generalised nested inference and the marginalisation of prior length scales for better posterior subsurface property characterisation.

Description: Summary The Z-transform of a complex time signal (or the analytic signal of a real signal) is equal to the Z-transform of a prediction error divided by the Z-transform of the prediction error operator. This inverse is decomposed into a sum of partial fractions, which are used to obtain impulse response operators formed by non-causal filters which complex-conjugate symmetric coefficients. The time-components are obtained by convolving the filters with the original signal, and the peak frequencies, corresponding to the poles of the prediction error operator, are used for mapping the time-components into frequency components. For non-stationary signals, this decomposition is done in sliding time windows, and the signal component values, in the middle of each window, are attributed to the peak value of its frequency response which corresponds to the pole of this partial fraction component. The result is an exact, but non-unique, time-frequency representation of the input signal. A sparse signal decomposition can be obtained by summing along the frequency axis in patches with similar characteristics in the time-frequency domain. The peak amplitude frequency of each new time component is obtained by computing a scalar prediction error operator in sliding time windows, resulting in a sparse time-frequency representation. In both cases, the result is a time-frequency matrix where an estimate of the frequency content of the input signal can be obtained by summation over the time variable. The performance of the new method is demonstrated with excellent results on a synthetic time signal, the LIGO gravitational wave signal, and on seismic field data.

Description: Summary Physical properties of near-surface soil and rock layers play a fundamental role in the seismic site effects analysis, being an essential element of seismic hazard assessment. Site-specific mechanical properties (i.e. shear- and compressional-wave velocities and mass density) can be inferred from surface wave dispersion and horizontal-to-vertical or ellipticity data by non-linear inversion techniques. Nevertheless, results typically exhibit significant inherent non-uniqueness as different models may fit the data equally well. Standard optimization inversion techniques minimize data misfit, resulting in a single representative model, rejecting other models providing similar misfit values. An alternative inversion technique can be formulated in the Bayesian framework, where the posterior probability density on the model space is inferred. This paper introduces an inversion approach of surface wave dispersion and ellipticity data based on a novel multizonal transdimensional Bayesian formulation. In particular, we parameterize one-dimensional layered velocity models by the varying number of Voronoi nuclei, allowing us to treat the number of layers as an unknown parameter of the inverse problem. The chosen parameterization leads to the transdimensional formulation of the model space, sampled by a reversible jump Markov chain Monte Carlo algorithm to provide an ensemble of random samples following the posterior probability density of model parameters. The used type of the sampling algorithm controls a model complexity (i.e. the number of layers) self-adaptively based on the measured data's information content. The method novelty lies in the parsimonious selection of sampling models and in the multizonal formulation of prior assumptions on model parameters, the latter allows including additional site-specific constraints in the inversion. These assumptions may be based on, e.g. stratigraphic logs, standard penetration tests, known water table, and bedrock depth. The multizonal formulation fully preserves the validity of the transdimensional one, as demonstrated analytically. The resultant ensemble of model samples is a discrete approximation of the posterior probability density function of model parameters and associated properties (e.g. VS30, quarter-wavelength average velocity profile, and theoretical SH-wave amplification function). Although the ultimate result is the posterior probability density function, some representative models are selected according to data fit and maximum of the posterior probability density function. We first validate our inversion approach based on synthetic tests and then apply it to field data acquired from the active seismic survey and single-station measurements of ambient vibrations at the SENGL seismic station site in central Switzerland.

Description: Summary Until now, the polar motion resonance (PMR) complex frequency has been determined in the seasonal and retrograde diurnal band of the polar motion. In this study this resonance is studied in the prograde diurnal band, where polar motion is mainly composed of periodic terms caused by the diurnal oceanic tide. The resonance parameters (period and quality factor) are encompassed in the frequency transfer function between generating tidal potential and polar motion, and can be estimated accordingly. To this aim, we gather three published sets of prograde diurnal terms determined from GNSS and VLBI, to which we append our own estimates based upon a processing of the VLBI delays over the period 1990-2020. Then, by fitting the PMR parameters so that the prograde diurnal terms match the corresponding components of the tide generating potential, we obtained a resonance period of about 401 days and an equivalent quality factor of −22, differing from the ones reigning in the seasonal band (PPMR ≈ 431 days; QPMR ≈ 56 − 255) and in the retrograde diurnal band (PPMR ≈ 380 days; QPMR ≈ −10). Our estimates confirm strikingly the theoretical prediction derived from the tidal ocean angular momentum derived from the FES 2014 ocean tide model.

Description: Summary Bayesian inversion of electromagnetic data produces crucial uncertainty information on inferred subsurface resistivity. Due to their high computational cost, however, Bayesian inverse methods have largely been restricted to computationally expedient 1D resistivity models. In this study, we successfully demonstrate, for the first time, a fully 2D, trans-dimensional Bayesian inversion of magnetotelluric data. We render this problem tractable from a computational standpoint by using a stochastic interpolation algorithm known as a Gaussian process to achieve a parsimonious parametrization of the model vis-a-vis the dense parameter grids used in numerical forward modeling codes. The Gaussian process links a trans-dimensional, parallel tempered Markov chain Monte Carlo sampler, which explores the parsimonious model space, to MARE2DEM, an adaptive finite element forward solver. MARE2DEM computes the model response using a dense parameter mesh with resistivity assigned via the Gaussian process model. We demonstrate the new trans-dimensional Gaussian process sampler by inverting both synthetic and field magnetotelluric data for 2D models of electrical resistivity, with the field data example converging within 10 days on 148 cores, a non-negligible but tractable computational cost. For a field data inversion, our algorithm achieves a parameter reduction of over 32x compared to the fixed parameter grid used for the MARE2DEM regularized inversion. Resistivity probability distributions computed from the ensemble of models produced by the inversion yield credible intervals and interquartile plots that quantitatively show the non-linear 2D uncertainty in model structure. This uncertainty could then be propagated to other physical properties that impact resistivity including bulk composition, porosity and pore-fluid content.

Description: Summary Ambient wavefield data acquired on existing (so-called “dark fiber”) optical fiber networks using distributed acoustic sensing (DAS) interrogators allow users to conduct a wide range of subsurface imaging and inversion experiments. In particular, recorded low-frequency (

Description: Summary The stability of partial differential equations determines the properties of their solutions. This study focuses on the stability analysis of the equations describing wave propagation in fluids-saturated porous media. We briefly introduce the stability analysis method for the wave propagation equations and discuss the adverse effects on the solutions. In this way, the first part of this paper is mainly devoted to the analysis of the Tuncay and Corapcioglu's (TC) model, which describes the dynamic behaviour of porous media saturated with two immiscible fluids. It is pointed out that the TC model allows spatially bounded but time-exponentially exploding solutions and may yield unstable numerical results. Based on the deduced unstable factors, we construct a stable equivalent fluid (SEF) model. We rigorously analyze the stability of the SEF model using the energy method. For predicting the influence of saturation on wave velocity, the robustness of this model is preserved due to its consistency with the original TC model. Furthermore, the numerical simulations of the wave fields show that the results of the TC model exponentially increase with time after the initial effective wave signal, which does not occur in the SEF model curves. This indicates the necessity of considering the stability from a mathematical point of view during the construction of physical model. It could be useful to merge the mathematical stability theory with the geophysical wave propagation modelling theory.

Description: Summary Uncovering the distribution of magnitudes and arrival times of aftershocks is a key to comprehending the characteristics of earthquake sequences, which enables us to predict seismic activities and conduct hazard assessments. However, identifying the number of aftershocks immediately after the main shock is practically difficult due to contaminations of arriving seismic waves. To overcome this difficulty, we construct a likelihood based on the detected data, incorporating a detection function to which Gaussian process regression (GPR) is applied. The GPR is capable of estimating not only the parameters of the distribution of aftershocks together with the detection function, but also credible intervals for both the parameters and the detection function. The property that the distributions of both the Gaussian process and aftershocks are exponential functions leads to an efficient Bayesian computational algorithm to estimate hyperparameters. After its validation through numerical tests, the proposed method is retrospectively applied to the catalog data related to the 2004 Chuetsu earthquake for the early forecasting of the aftershocks. The results show that the proposed method stably and simultaneously estimates distribution parameters and credible intervals, even within t ≤ 3h after the main shock.

Description: Summary We reconstruct the post-52 Ma seafloor spreading history of the Southwest Indian Ridge at 44 distinct times from inversions of ≈20,000 magnetic reversal, fracture zone, and transform fault crossings, spanning major regional tectonic events such as the Arabia-Eurasia continental collision, the Arabia Peninsula’s detachment from Africa, the arrival of the Afar mantle plume below eastern Africa, and the initiation of rifting in eastern Africa. Best-fitting and noise-reduced rotation sequences for the Nubia-Antarctic, Lwandle-Antarctic, and Somalia-Antarctic plate pairs indicate that spreading rates everywhere along the ridge declined gradually by ≈50 percent from ≈31 Ma to 19-18 Ma. A concurrent similar-magnitude slowdown in the component of the Africa plate’s absolute motion parallel to Southwest Indian Ridge spreading suggests that both were caused by a 31-18 Ma change in the forces that drove and resisted Africa’s absolute motion. Possible causes for this change include the effects of the Afar mantle plume on eastern Africa or the Arabia Peninsula’s detachment from the Somalia plate, which culminated at 20-18 Ma with the onset of seafloor spreading in the Gulf of Aden. At earlier times, an apparently robust but previously unknown ≈6-Myr-long period of rapid kinematic change occurred from 43 Ma to 37 Ma, consisting of a ≈50 percent spreading rate slowdown from 43-40 Ma followed by a full spreading rate recovery and 30-40 ○ clockwise rotation of the plate slip direction from 40-37 Ma. Although these kinematic changes coincided with a reconfiguration of the paleoridge geometry, their underlying cause is unknown. Southwest Indian Ridge abyssal hill azimuths are consistent with the slip directions estimated with our newly derived Somalia-Antarctic and Lwandle-Antarctic angular velocities, adding confidence in their reliability. Lwandle-Antarctica plate motion has closely tracked Somalia-Antarctic plate motion since 50 Ma, consistent with slow-to-no motion between the Lwandle and Somalia plates for much of that time. In contrast, Nubia-Somalia rotations estimated from our new Southwest Indian Ridge rotations indicate that 189±34 km of WNW-ESE divergence between Nubia and Somalia has occurred in northern Africa since 40 Ma, including 70-80 km of WNW-ESE divergence since 17-16 Ma, slow to no motion from 26-17 Ma, and 109±38 km of WNW-ESE divergence from 40 Ma to ≈26 Ma absent any deformation within eastern Antarctica before 26 Ma.

Description: Summary Most of the existing three-dimensional (3-D) electromagnetic (EM) modeling solvers based on the integral equation (IE) method exploit fast Fourier transform (FFT) to accelerate the matrix-vector multiplications. This in turn requires a laterally-uniform discretization of the modeling domain. However, there is often a need for multi-scale modeling and inversion, for instance, to properly account for the effects of non-uniform distant structures, and at the same time, to accurately model the effects from local anomalies. In such scenarios, the usage of laterally-uniform grids leads to excessive computational loads, both in terms of memory and time. To alleviate this problem, we developed an efficient 3-D EM modeling tool based on a multi-nested IE approach. Within this approach, the IE modeling is first performed at a large domain and on a (laterally-uniform) coarse grid, and then the results are refined in the region of interest by performing modeling at a smaller domain and on a (laterally-uniform) denser grid. At the latter stage, the modeling results obtained at the previous stage are exploited. The lateral uniformity of the grids at each stage allows us to keep using the FFT for the matrix-vector multiplications. An important novelty of the paper is a development of a “rim domain” concept which further improves the performance of the multi-nested IE approach. We verify the developed tool on both idealized and realistic 3-D conductivity models, and demonstrate its efficiency and accuracy.

Description: Summary Denoising and onset time picking of signals are essential before extracting source information from collected seismic/microseismic data. We proposed an advanced deep dual-tasking network (DDTN) that integrates these two procedures sequentially to achieve the optimal performance. Two homo-structured encoder-decoder networks with specially designed structures and parameters are connected in series for handling the denoising and detection of microseismic signals. Based on the similarity of data types, the output of the denoising network will be imported into the detection network to obtain labels for the signal duration. The procedures of denoising and duration detection can be completed in an integrated way, where the denoised signals can improve the accuracy of onset time picking. Results show that the method has a good performance for the denoising of microseismic signals that contain various types and intensities of noise. Compared with existing methods, DDTN removes the noise with a minor waveform distortion. It is ideal for recovering the microseismic signal while maintaining a good capacity for onset time picking when the signal-to-noise ratio is low. Based on that, the second network can detect a more accurate duration of microseismic signals and thus obtain more accurate onset time. The method has great potential to be extended to the study of exploration seismology and earthquakes.

Description: Summary The ongoing InSight mission has recently deployed very broad band seismometers to record the Martian seismic activity. These recordings constitute the first seismic data set collected at the surface of Mars. This unique but sparse record compels for the development of new techniques tailored to make the best use of the specific context of single station-multiple events with several possible ranges of uncertainties on the event location. To this end, we conducted sets of Markov chain Monte-Carlo inversions for the 1-D seismic structure of Mars. We compared two inversion techniques that differ from the nature of the parameterization on which they rely. A first classical approach based on a parameterization of the 1-D seismic profile using Bézier curves. A second, less conventional approach that relies on a parameterization in terms of quantities that influence the thermo-chemical evolution of the planet (mantle rheology, initial thermal state, and composition), which accounts for 4.5 Gyr of planetary evolution. We considered several combinations of true model parameters to retrieve, and explored the influence of the type of seismic data (body waves with or without surface waves), the number of events and their associated epicentral distances and uncertainties, and the presence of potential constraints on Moho depth inferred from independent measurements/considerations (receiver functions and gravity data). We show that due to its inherent tighter constraints the coupled approach allows a considerably better retrieval of Moho depth and the seismic structure underneath it than the classical inversion, under the condition that the physical assumptions made in coupled approach are valid for Mars. In addition, our tests indicate that in order to constrain the seismic structure of Mars with InSight data, the following independent conditions must be met: (1) The presence of surface waves triggered by an internal source to constrain the epicentral distance. (2) The presence of just a few well-localized impact sources, with at least one located at close epicentral distance (

Description: Summary In the field of seismic interferometry, cross-correlations are used to extract Green’s function from ambient noise data. By applying a single station variation of the method, using auto-correlations, we are in principle able to retrieve zero-offset reflections in a stratified Earth. These reflections are valuable as they do not require an active seismic source and, being zero-offset, are better constrained in space than passive earthquake based measurements. However, studies that target Moho signals with ambient noise auto-correlations often give ambiguous results with unclear Moho reflections. Using a modified processing scheme and phase-weighted stacking, we determine the Moho P wave reflection time from vertical auto-correlation traces for a test station with a known simple crustal structure (HYB in Hyderabad, India). However, in spite of the simplicity of the structure, the auto-correlation traces show several phases not related to direct reflections. Although we are able to match some of these additional phases in a qualitative way with synthetic modelling, their presence makes it hard to identify the reflection phases without prior knowledge. This prior knowledge can be provided by receiver functions. Receiver functions (arising from mode conversions) are sensitive to the same boundaries as auto-correlations, so should have a high degree of comparability and opportunity for combined analysis but in themselves are not able to independently resolve VP, VS, and Moho depth. Using the timing suggested by the receiver functions as a guide, we observe the Moho S wave reflection on the horizontal auto-correlation of the north component but not on the east component. The timing of the S reflection is consistent with the timing of the PpSs-PsPs receiver function multiple, which also depends only on the S velocity and Moho depth. Finally, we combine P receiver functions and auto-correlations from HYB in a depth-velocity stacking scheme that gives us independent estimates for VP, VS, and Moho depth. These are found to be in good agreement with several studies that also supplement receiver functions to obtain unique crustal parameters. By applying the auto-correlation method to a portion of the EASI transect crossing the Bohemian Massif in central Europe, we find approximate consistency with Moho depths determined from receiver functions and spatial coherence between stations, thereby demonstrating that the method is also applicable for temporary deployments. Although application of the auto-correlation method requires great care in phase identification, it has the potential to resolve both average crustal P and S velocities alongside Moho depth in conjunction with receiver functions.

Description: Summary By combining scaled laboratory experiments and numerical simulations, this study presents a quantitative analysis of the bending radius (RB) of subducting slabs within the upper mantle, taking into account the effects of age (A). Based on a half-space cooling model, we constrain the density (ρ), viscosity (η) and thickness (h) of slabs as a function of A, and develop representative models to estimate RB for different A. Laboratory subduction models produce visually contrasting bending curvatures for young (A = 10 Ma), intermediate (A = 70 Ma) and old (A = 120 Ma) slabs. Young slabs undergo rollback, resulting in a small bending radius (scaled up RB ∼ 150 km), whereas old slabs subduct along a uniformly dipping trajectory with large bending radius (RB ∼ 500 km). Equivalent real scale computational fluid dynamic (CFD) simulations reproduce similar bending patterns of the subducting slabs, and yield RB versus A relations fairly in agreement with the laboratory results. We balance the buoyancy driven bending, flexural-resistive moments and viscous flow induced suction moment to theoretically evaluate the rate of slab bending. The analytical solution suggests an inverse relation of the bending rate with A, which supports our experimental findings. Finally, slab geometries of selected natural subduction zones, derived from high-resolution seismic tomographic images have been compiled to validate the experimental RB versus A regression. We also discuss the subduction settings in which this regression no longer remains valid.

Description: Summary As recovering the crust-mantle/Moho density contrast (MDC) significantly depends on the properties of the Earth's crust and upper mantle, varying from place to place, it is an oversimplification to define a constant standard value for it. It is especially challenging in Antarctica, where almost all the bedrock is covered with a thick layer of ice, and seismic data cannot provide a sufficient spatial resolution for geological and geophysical applications. As an alternative we determine the MDC in Antarctica and its surrounding seas with a resolution of 1° × 1° by the Vening Meinesz-Moritz gravimetric-isostatic technique using the XGM2019e Earth Gravitational Model and Earth2014 topographic/bathymetric information along with CRUST1.0 and CRUST19 seismic crustal models. The numerical results show that our model, named HVMDC20, varies from 81 kg/m3 in the Pacific Antarctic mid-oceanic ridge to 579 kg/m3 in the Gamburtsev Mountain Range in the central continent with a general average of 403 kg/m3. To assess our computations, we compare our estimates with those of some other gravimetric as well as seismic models (KTH11, GEMMA12C, KTH15C, and CRUST1.0), illustrating that our estimates agree fairly well with KTH15C and CRUST1.0 but rather poor with the other models. In addition, we compare the geological signatures with HVMDC20, showing how the main geological structures contribute to the MDC. Finally, we study the remaining glacial isostatic adjustment effect on gravity to figure out how much it affects the MDC recovery, yielding a correlation of the optimum spectral window (7≤ n ≤12) between XGM2019e and W12a GIA model of the order of ∼ 0.6 contributing within a negligible $pm 14$ kg/m3 to the MDC.

Description: Summary Spectral information obtained from induced polarization (IP) measurements can be used in a variety of applications and is often gathered in frequency domain (FD) at the laboratory scale. In contrast, field IP measurements are mostly done in time domain (TD). Theoretically, the spectral content from both domains should be similar. In practice, they are often different, mainly due to instrumental restrictions as well as the limited time and frequency range of measurements. Therefore, a possibility of transition between both domains, in particular for the comparison of laboratory FD IP data and field TD IP results, would be very favourable. To compare both domains, we conducted laboratory IP experiments in both TD and FD. We started with three numerical models and measurements at a test circuit, followed by several investigation at different wood and sandstone samples. Our results demonstrate that the differential polarizability (DP), which is calculated from the TD decay curves, can be compared very well with the phase of the complex electrical resistivity. Thus, DP can be used for a first visual comparison of FD and TD data, which also enables a fast discrimination between different samples. Furthermore, to compare both domains qualitatively, we calculated the relaxation time distribution (RTD) for all data. The results are mostly in agreement between both domains, however, depending on the TD data quality. It is striking that the DP and RTD results are in better agreement for higher data quality in TD. Nevertheless, we demonstrate that IP laboratory measurements can be carried out in both TD and FD with almost equivalent results. The RTD enables a good comparability of FD IP laboratory data with TD IP field data.

Description: Summary Isostasy explains why observed gravity anomalies are generally much weaker than what is expected from topography alone, and why planetary crusts can support high topography without breaking up. On Earth, it is used to subtract from gravity anomalies the contribution of nearly compensated surface topography. On icy moons and dwarf planets, it constrains the compensation depth which is identified with the thickness of the rigid layer above a soft layer or a global subsurface ocean. Classical isostasy, however, is not self-consistent, neglects internal stresses and geoid contributions to topographical support, and yields ambiguous predictions of geoid anomalies. Isostasy should instead be defined either by minimizing deviatoric elastic stresses within the silicate crust or icy shell, or by studying the dynamic response of the body in the long-time limit. In this paper, I implement the first option by formulating Airy isostatic equilibrium as the linear response of an elastic shell to a combination of surface and internal loads. Isostatic ratios are defined in terms of deviatoric Love numbers which quantify deviations with respect to a fluid state. The Love number approach separates the physics of isostasy from the technicalities of elastic-gravitational spherical deformations, and provides flexibility in the choice of the interior structure. Since elastic isostasy is invariant under a global rescaling of the shell shear modulus, it can be defined in the fluid shell limit, which is simpler and reveals the deep connection with the asymptotic state of dynamic isostasy. If the shell is homogeneous, minimum stress isostasy is dual to a variant of elastic isostasy called zero deflection isostasy, which is less physical but simpler to compute. Each isostatic model is combined with general boundary conditions applied at the surface and bottom of the shell, resulting in one-parameter isostatic families. At long wavelength, the thin shell limit is a good approximation, in which case the influence of boundary conditions disappears as all isostatic families members yield the same isostatic ratios. At short wavelength, topography is supported by shallow stresses so that Airy isostasy becomes similar to either pure top loading or pure bottom loading. The isostatic ratios of incompressible bodies with three homogeneous layers are given in analytical form in the text and in complementary software.

Description: Summary The Pollino range is a region of slow deformation where earthquakes generally nucleate on low-angle normal faults. Recent studies have mapped fault structures and identified fluid-related dynamics responsible for historical and recent seismicity in the area. Here, we apply the coda-normalization method at multiple frequencies and scales to image the 3D P-wave attenuation (QP) properties of its slowly-deforming fault network. The wide-scale average attenuation properties of the Pollino range are typical for a stable continental block, with a dependence of QP on frequency of $Q_P^{-1}=(0.0011pm 0.0008) f^{(0.36pm 0.32)}$. Using only waveforms comprised in the area of seismic swarms, the dependence of attenuation on frequency increases ($Q_P^{-1}=(0.0373pm 0.0011) f^{(-0.59pm 0.01)}$), as expected when targeting seismically-active faults. A shallow very-low-attenuation anomaly (max depth of 4-5 km) caps the seismicity recorded within the western cluster 1 of the Pollino seismic sequence (2012, maximum magnitude MW = 5.1). High-attenuation volumes below this anomaly are likely related to fluid storage and comprise the western and northern portions of cluster 1 and the Mercure basin. These anomalies are constrained to the NW by a sharp low-attenuation interface, corresponding to the transition towards the eastern unit of the Apennine Platform under the Lauria mountains. The low-seismicity volume between cluster 1 and cluster 2 (maximum magnitude MW = 4.3, east of the primary) shows diffuse low-to-average attenuation features. There is no clear indication of fluid-filled pathways between the two clusters resolvable at our resolution. In this volume, the attenuation values are anyway lower than in recognized low-attenuation blocks, like the Lauria Mountain and Pollino Range. As the volume develops in a region marked at surface by small-scale cross-faulting, it suggests no actual barrier between clusters, more likely a system of small locked fault patches that can break in the future. Our model loses resolution at depth, but it can still resolve a 5-to-15-km-deep high-attenuation anomaly that underlies the Castrovillari basin. This anomaly is an ideal deep source for the SE-to-NW migration of historical seismicity. Our novel deep structural maps support the hypothesis that the Pollino sequence has been caused by a mechanism of deep and lateral fluid-induced migration.

Description: Summary Perturbation method applied for the phase velocities in the monoclinic model with a horizontal symmetry plane gives insight into the seismic signatures of P, S1 and S2 wave modes associated with the monoclinic stiffness coefficients. The perturbation-based approximations are very accurate and can be used for modeling and inversion purpose.

Description: Summary We present a novel approach for imaging global mantle discontinuities based on full-waveform inversion (FWI). Over the past decades, extensive research has been done on imaging mantle discontinuities at approximately 400 km and 670 km depth. Accurate knowledge of their topography can put strong constraints on thermal and compositional variations and hence geodynamic modelling. So far, however, there is little consensus on their topography. We present an approach based on adjoint tomography, which has the advantage that Fréchet derivatives for discontinuities and measurements, to be inverted for, are fully consistent. Rather than working with real data, we focus on synthetic tests, where the answer is known in order to be able to evaluate the performance of the developed method. All calculations are based on the community code SPECFEM3D_GLOBE. We generate data in fixed 1-D or 3-D elastic background models of mantle velocity. Our ‘data’ to be inverted contain topography along the 400 km and 670 km mantle discontinuities. To investigate the approach, we perform several tests: (i) In a situation where we know the elastic background model 1-D or 3-D, we recover the target topography fast and accurately, (ii) The exact misfit is not of great importance here, except in terms of convergence speed, similar to a different inverse algorithm, (iii) In a situation where the background model is not known, the convergence is markedly slower, but there is reasonable convergence towards the correct target model of discontinuity topography. It has to be noted that our synthetic test is idealised and in a real data situation, the convergence to and uncertainty of the inferred model is bound to be larger. However, the use of data consistent with Fréchet kernels seems to pay off and might improve our consensus on the nature of mantle discontinuities. Our workflow could be incorporated in future FWI mantle models to adequately infer boundary interface topography.

Description: Summary This study investigates the seismic structure and anisotropy in the crust beneath Madagascar and south-eastern Africa, using receiver functions. The understanding of seismic anisotropy is essential for imaging past and present deformation in the lithosphere-asthenosphere system. In the upper mantle, seismic anisotropy mainly results from the orientation of olivine, which deforms under tectonic (fossil anisotropy) or flow processes (in the asthenosphere). In the crust, the crystallographic alignment of amphiboles, feldspars(plagioclase) or micas or the alignment of heterogeneities such as fractures, add to a complex geometry, which results in challenges to understanding the Earth's shallow structure. The decomposition of receiver functions into back-azimuth harmonics allows to characterize orientations of lithospheric structure responsible for azimuthally-varying seismic signals, such as a dipping isotropic velocity contrasts or layers of azimuthal seismic anisotropy. By analysing receiver function harmonics from records of 48 permanent or temporary stations this study reveals significant azimuthally-varying signals within the upper crust of Madagascar and south-eastern Africa. At 30 stations crustal anisotropy dominates the harmonics while the signature of a dipping isotropic contrast is dominant at the remaining 18 stations. However, all stations’ back-azimuth harmonics show complex signals involving both dipping isotropic and shallow anisotropic contrasts or more than one source of anisotropy at shallow depth. Our calculated orientations for the crust are therefore interpreted as reflecting either the average or the interplay of several sources of azimuthally-varying signals depending of their strength. However, comparing information between stations allows us to draw the same conclusions regionally: in both southern Africa and Madagascar our measurements reflect the interplay between local, inherited structural heterogeneities and crustal seismic anisotropy generated by the current extensional stress field imposed by the southward propagation of the East-African Rift System. A final comparison of our crustal orientations with SKS orientations attributed to mantle deformation further probes the interplay of crustal and mantle anisotropy on SKS measurements.

Description: Summary During the Cenozoic, the Earth experienced multiple first-order geologic events that are likely mantle flow related. These include the termination of large-scale marine inundation in North America in the Paleocene, the late Tertiary rise of Africa relative to other continents, and the long wavelength tilting of Australia since the late Cretaceous, which occurred when the continent approached the South East Asia subduction systems on its northward passage from Antartica. Here we explore a suite of eight high resolution, compressible, global mantle flow retrodictions going back 50 million years ago, using an adoint method with ≈670 million finite elements. These retrodictions show for the first time that these events emerge jointly as part of global Cenozoic mantle flow histories. Our retrodictions involve the dynamic effects from an upper mantle low viscosity zone (LVZ), assimilate a past plate motion model for the tangential surface velocity field, probe the influence of two different present-day mantle state estimates derived from seismic tomography, and acknowledge the rheological uncertainties of dynamic Earth models by taking in four different realizations for the radial mantle viscosity profile, two of which were published previously. We find the retrodicted mantle flow histories are sensitive to the present-day mantle state estimate and the rheological properties of the Earth model, meaning that this input information is testable with inferences gleaned from the geologic record. For a deep mantle viscosity of 1.7 × 1022 Pa.s and a purely thermal interpretation of seismic structure, lower mantle flow velocities exceed 7 cm/yr in some regions, meaning they are difficult to reconcile with the existence of a hotspot reference frame. Conversely, a deep mantle viscosity of 1023 Pa.s yields modest flow velocities (

Description: Summary The rate of earthquakes with magnitudes Mw ≤ 7.5 in the Ometepec segment of the Mexican subduction zone is relatively high as compared to the neighboring regions of Oaxaca and Guerrero. Although the reason is not well understood, it has been reported that these earthquakes give rise to a large number of aftershocks. Our study of the aftershock sequence of the 2012 Mw7.4 Ometepec thrust earthquake suggests that it is most likely due to two dominant factors: (1) The presence of an anomalously high quantity of over-pressured fluids near the plate interface, and (2) the roughness of the plate interface. More than 5,400 aftershocks were manually detected during the first ten days following the 2012 earthquake. Locations were obtained for 2,419 events (with duration magnitudes Md ≥ 1.5). This is clearly an unusually high number of aftershocks for an earthquake of this magnitude. Furthermore, we generated a more complete catalog, using an unsupervised fingerprint technique, to detect more smaller events (15,593 within one month following the mainshock). For this catalog, a high b-value of 1.50 ± 0.10 suggests the presence of fluid release during the aftershock sequence. A low p-value (0.37 ± 0.12) of the Omori law reveals a slow decaying aftershock sequence. The temporal-distribution of aftershocks shows peaks of activity with two dominant periods of 12h and 24h that correlate with the Earth tides. To explain these observations, we suggest that the 2012 aftershock sequence is associated with the presence of over-pressured fluids and/or a heterogeneous and irregular plate interface related to the subduction of the neighboring seamounts. High fluid content has independently been inferred by magneto-telluric surveys and deduced from heat flow measurements in the region. The presence of fluids in the region has also been proposed to explain the occurrence of slow slip events, low frequency earthquakes, and tectonic tremors.

Description: Summary The increasing density of geodetic measurements makes it possible to map surface strain rate in many zones of active tectonics with unprecedented spatial resolution. Here we show that the strain tensor rate calculated from GPS in the India-Asia collision zone represents well the strain released in earthquakes. This means that geodetic data in the India-Asia collision zone region can be extrapolated back in time to estimate strain buildup on active faults, or the kinematics of continental deformation. We infer that the geodetic strain rates can be assumed stationary through time on the timescale needed to build up the elastic strain released by larger earthquakes, and that they can be used to estimate the probability of triggering earthquakes. We show that the background seismicity rate correlates with the geodetic strain rate. A good fit is obtained assuming a linear relationship ($dot{N} = lambda cdot dot{epsilon }$ where $dot{N}$ is the density of the rate of Mw ≥ 4 earthquakes, $dot{epsilon }$ is strain rate and λ = 2.5 ± 0.1 × 10−3 m−2), as would be expected from a standard Coulomb failure model. However, the fit is significantly better for a non-linear relationship ($dot{N} = gamma _1 cdot dot{epsilon }^{gamma _2}$ with γ1 = 2.5 ± 0.6 m−2 and γ2 = 1.42 ± 0.15). The b-value of the Gutenberg-Richter law, which characterize the magnitude-frequency distribution, is found to be insensitive to the strain rate. In the case of a linear correlation between seismicity and strain rate, the maximum magnitude earthquake, derived from the moment conservation principle, is expected to be independent of the strain rate. By contrast, the non-linear case implies that the maximum magnitude earthquake would be larger in zones of low strain rate. We show that within areas of constant strain rate, earthquakes above Mw4 follow a Poisson distribution in time and and are uniformly distributed in space. These findings provide a framework to estimate the probability of occurrence and magnitude of earthquakes as a function of the geodetic strain rate. We describe how the seismicity models derived from this approach can be used as an input for probabilistic seismic hazard analysis. This method is easy to automatically update, and can be applied in a consistent manner to any continental zone of active tectonics with sufficient geodetic coverage.

Description: SUMMARY The Kohistan Island Arc (KIA) occupies the northwestern region of the Himalayan Mountains, sandwiched between Asia and India plates. Its formation, collision with plate boundaries, and evolution has been controversially discussed for a couple of decades. To better understand this, a palaeomagnetic study has been conducted on the Jutal dykes (ca. 75 Ma), intruded in the northeastern part of the KIA. Comprehensive rock magnetic investigations reveal that the magnetic carrier minerals are pyrrhotite and magnetite. An intermediate temperature component (ITC) predominates the natural remanent magnetization and shows good coincidence within-site; it is carried by pyrrhotite and is considered reliable, yielding a mean direction at Dg/Ig = 11.5°/39.9° (kg = 28.4, α95 = 3.5°) before and Ds/Is = 8.6°/12.1° (ks = 5.1, α95 = 9.1°) after tilt correction. A high-temperature component that is carried by magnetite exhibits random distribution within-site. The fold test for the ITC is negative, indicating a post-folding origin. Scanning electron microscopy combined with energy-dispersive X-ray spectroscopy indicates that the magnetic carrier minerals were influenced by metamorphism or thermochemical fluids. The comparison of mean palaeolatitude (22.6 ± 3.5°N) of the ITC with the collisional settings and thermal history of the study area implies that the remagnetization occurred at ∼50–35 Ma, consistent with the previous reported palaeomagnetic data of the KIA. We propose a tectonic model that shows the evolution of the Jutal dykes, supporting the concept that India collided with the KIA first, followed by a later collision with Asia.

Description: SUMMARY Interfaces are important part of Earth’s layering structure. Here, we developed a new model parametrization and iterative linearized inversion method that determines 1-D crustal velocity structure using surface wave dispersion, teleseismic P-wave receiver functions and Ps and PmP traveltimes. Unlike previous joint inversion methods, the new model parametrization includes interface depths and layer Vp/Vs ratios so that smoothness constraint can be conveniently applied to velocities of individual layers without affecting the velocity discontinuity across the interfaces. It also allows adding interface-related observation such as traveltimes of Ps and PmP in the joint inversion to eliminate the trade-off between interface depth and Vp/Vs ratio and therefore to reduce the uncertainties of results. Numerical tests show that the method is computationally efficient and the inversion results are robust and independent of the initial model. Application of the method to a dense linear array across the Wabash Valley Seismic Zone (WVSZ) produced a high-resolution crustal image in this seismically active region. The results show a 51–55-km-thick crust with a mid-crustal interface at 14–17 km. The crustal Vp/Vs ratio varies from 1.69 to 1.90. There are three pillow-like, ∼100 km apart high-velocity bodies sitting at the base of the crust and directly above each of them are a low-velocity anomaly in the middle crust and a high-velocity anomaly in the upper crust. They are interpreted to be produced by mantle magmatic intrusions and remelting during rifting events in the end of the Precambrian. The current diffuse seismicity in the WVSZ might be rooted in this ancient distributed rifting structure.

Description: Summary The recently constructed diffuse field theory from isotropic energy equipartition has been well developed in elasticity for full-wave interpretation of Horizontal-to-Vertical Ratio (HVSR), which links the signal autocorrelation with the imaginary part of Green's function. Here, the theory is extended to the saturated layered medium within the framework of Biot's theory to account for the offshore environment. The imaginary parts of Green's functions are obtained using direct stiffness method accompanied with Fourier-Hankel transform. In particular, the up-going wave amplitudes are modified to tackle the overflow during wavenumber integral and allow for fast calculations. After validating the method from the perspectives of Green's function calculation, emphasis is laid on evaluating the inaccuracies of HVSR calculation induced by model misuses in the lack of prior geological and geotechnical information. The numerical results considering the effects of layer sequence, impedance ratio, porosity and drainage condition show that the predominant frequency of the one-phase medium is slightly less than the two-phase medium with the maximum shift no more than 0.1Hz, while their amplitude differences can be prominent as impedance ratio and porosity increase, with the maximum difference up to 29 per cent. The shallowest soft layer has the dominant effects on HVSR amplitudes, whereas the buried low-velocity layer at depth over one-wavelength contributes little to the peak amplitude. Finally, the method is applied to a realistic case at Mirandola, Northorn Italy, which suffered extensive liquefaction-induced damages in 2012 Emilia earthquake. The well identified pattern of the experimental HVSR using the two-phase medium model illustrates the application potential of our method to further assist the subsurface geology retrieval.

Description: Summary Delineating spatial variations of seismic anisotropy in the crust is of great importance for the understanding of structural heterogeneities, regional stress regime and ongoing crustal dynamics. In this study we present a three-dimensional (3-D) anisotropic P-wave velocity model of the crust beneath northern California (2-22 km) by using the eikonal equation-based seismic azimuthal anisotropy tomography method. The velocity heterogeneities under different geological units are well resolved. The thickness of the low-velocity sediment at the Great Valley Sequence is estimated to be about 10 km. The high-velocity anomaly underlying Great Valley probably indicates the existence of ophiolite bodies. Strong velocity contrasts are revealed across the Hayward Fault (2-9 km) and San Andreas Fault (2-12 km). In the upper crust (2-9 km), the fast velocity directions (FVDs) at different depths are generally fault-parallel in the northern Coast Range, which may be caused by geological structure; while the FVDs are mainly NE-SW in Great Valley and the northern Sierra Nevada possibly due to the regional maximum horizontal compressive stress. In contrast, seismic anisotropy in the mid-lower crust (12-22 km) may be attributed to the alignment of mica schists. The anisotropy contrast across the San Andreas Fault may imply different mechanisms of crustal deformation on the two sides of the fault. Both the strong velocity contrasts and the high angle (∼45o or above) between the FVDs and the strikes of faults suggest that the faults are mechanically weak in the San Francisco bay area (2-6 km). This study suggests that the eikonal equation-based seismic azimuthal anisotropy tomography is a valuable tool to investigate crustal heterogeneities and tectonic deformation.

Description: Summary Surface wave retrieval from ambient noise records using seismic interferometry techniques have been widely used for multiscale shear wave velocity (Vs) imaging. One key step during Vs imaging is the generation of dispersion spectra and the extraction of a reliable dispersion curve from the retrieved surface waves. However, the sparse array geometry usually affects the ability for high-frequency (〉1 Hz) seismic signals acquisition. Dispersion measurements are degraded by array response due to sparse sampling and often present smeared dispersion spectra with sidelobe artifacts. Previous studies usually focus on interferograms domain (e.g., cross-correlation function) and attempt to enhance coherent signals before dispersion measurement. We propose an alternative technique to explicitly deblur dispersion spectra through use of a phase-weighted slant-stacking algorithm. Numerical examples demonstrate the strength of the proposed technique to attenuate array responses as well as incoherent noise. Three different field examples prove the flexibility and superiority of the proposed technique: the first dataset consists of ambient noise records acquired using a nodal seismometer array; the second dataset utilizes Distributed Acoustic Sensing (DAS) and a marine fiber-optic cable to acquire a similar ambient noise dataset; the last dataset is a vibrator-based active-source surface wave data. The enhanced dispersion measurements provide cleaner and higher-resolution spectra without distortions which will assist both human interpreters as well as ML algorithms in efficiently picking curves for subsequent Vs inversion.

Description: Summary Mt. Merapi, which lies just north of the city of Yogyakarta in Java, Indonesia, is one of the most active and dangerous volcanoes in the world. Thanks to its subduction zone setting, Mt Merapi is a stratovolcano, and rises to an elevation of 2968 m above sea level. It stands at the intersection of two volcanic lineaments, Ungaran–Telomoyo–Merbabu–Merapi (UTMM) and Lawu–Merapi–Sumbing–Sindoro–Slamet, which are oriented north-south and west-east, respectively. Although it has been the subject of many geophysical studies, Mt Merapi's underlying magmatic plumbing system is still not well understood. Here, we present the results of an ambient seismic noise tomography study, which comprise of a series of Rayleigh wave group velocity maps and a 3-D shear wave velocity model of the Merapi-Merbabu complex. A total of 10 months of continuous data (October 2013–July 2014) recorded by a network of 46 broadband seismometers were used. We computed and stacked daily cross-correlations from every pair of simultaneously recording stations to obtain the corresponding inter-station empirical Green's functions. Surface wave dispersion information was extracted from the cross-correlations using the multiple filtering technique, which provided us with an estimate of Rayleigh wave group velocity as a function of period. The group velocity maps for periods 3–12 s were then inverted to obtain shear wave velocity structure using the neighbourhood algorithm. From these results, we observe a dominant high velocity anomaly underlying Mt. Merapi and Mt. Merbabu with a strike of 152° N, which we suggest is evidence of old lava dating from the UTMM double-chain volcanic arc which formed Merbabu and Old Merapi. We also identify a low velocity anomaly on the southwest flank of Merapi which we interpret to be an active magmatic intrusion.

Description: Summary The crustal attenuation structure can effectively reveal the rheology and thermal properties of different geological blocks, and can provide seismological constraints on regional tectonic evolution. Based on 11,306 vertical-component Lg-wave seismograms recorded at 111 broadband stations from 891 crustal earthquakes that occurred between 1994 and 2020, a broadband Lg-wave attenuation model is obtained for Southeast Asia. This study demonstrates the capability of applying crustal Lg wave attenuation inversion in a complex region mixed with continents, islands and marginal seas. The resolution approaches 2 degrees in most parts in the study region. Lg blockages are observed at places with sharp Moho depth changes. The resultant Q models are consistent with regional geologic structures provided by previous studies. Prominent low attenuation anomalies are located in the Sundaland Core containing stable ancient crust, including Indochina, Malay Peninsula, East Sumatra, Sunda Shelf and Borneo Core. Regions with strong attenuation are associated with complex tectonic conditions, such as the Indo-Australian subduction zone, sutures in Sarawak and Sabah. The observed Lg-wave attenuation characteristics provide constraints on the tectonic affinities and evolutions of the geological blocks. The results show that the Borneo Core remained stable since its accretion with the Sundaland Core. Ancient blocks are characterized by weak Lg attenuation, whereas geologically younger terranes are often characterized by strong Lg attenuation, which can be exploited to better understand the separation and convergence of plates during the tectonic evolution.

Description: Summary Shallow seismic sources excite Rayleigh wave ground motion with azimuthally dependent radiation patterns. We place binary hypothesis tests on theoretical models of such radiation patterns to screen cylindrically symmetric sources (like explosions) from non-symmetric sources (like non-vertical dip-slip, or non-VDS faults). These models for data include sources with several unknown parameters, contaminated by Gaussian noise and embedded in a layered half-space. The generalized maximum likelihood ratio tests that we derive from these data models produce screening statistics and decision rules that depend on measured, noisy ground motion at discrete sensor locations. We explicitly quantify how the screening power of these statistics increase with the size of any dip-slip and strike-slip components of the source, relative to noise (faulting signal strength), and how they vary with network geometry. As applications of our theory, we apply these tests to (1) find optimal sensor locations that maximize the probability of screening non-circular radiation patterns, and (2) invert for the largest non-VDS faulting signal that could be mistakenly attributed to an explosion with damage, at a particular attribution probability. Lastly, we quantify how certain errors that are sourced by opening cracks increase screening rate errors. While such theoretical solutions are ideal and require future validation, they remain important in underground explosion monitoring scenarios because they provide fundamental physical limits on the discrimination power of tests that screen explosive from non-VDS faulting sources.

Description: Summary We present a theory of modern, thermally-induced deformation in a realistic Earth. The heat conduction equation is coupled with standard elastic deformation theory to construct a boundary-value problem comprised of eighth-order differential equations. The accurate and stable dual variable and position propagating matrix technique is introduced to solve the boundary-value problem. The thermal load Love numbers are defined to describe the displacements and potential changes driven by thermally-induced deformation. The proposed analytical method is validated by comparing the present results with exact solutions for a homogeneous sphere, which are also derived in this paper. The analytical method is then applied to a realistic Earth model to evaluate the effects of layering and self-gravitation of the Earth on displacement and changes of potential. Furthermore, the frequency-dependence in the thermal load is illustrated by invoking different thermal periodicities in the computation. Thermal anisotropy is also considered by comparing the results obtained using isotropic and transversely isotropic Earth models. Results show that, when simulating thermally-induced deformation, invoking a homogeneous spherical Earth leads to results that substantially differ from those obtained using a more realistic Earth model.

Description: Summary Mantle convection induces dynamic topography, the lithosphere's surface deflections driven by the vertical stresses from sub-lithospheric mantle convection. Dynamic topography has important influences on a range of geophysical and geological observations. Here, we studied controls on the Earth's dynamic topography through three-dimensional spherical models of mantle convection, which use reconstructed past 410 Myr global plate motion history as time-dependent surface mechanical boundary condition. The numerical model assumes the extended-Boussinesq approximation and includes strongly depth- and temperature-dependent viscosity and phase changes in the mantle. Our results show that removing the chemical layer above the CMB and including depth-dependent thermal expansivity have both a limited influence on the predicted present-day dynamic topography. Considering phase transitions in our models increases the predicted amplitude of dynamic topography, which is mainly influenced by the 410 km exothermic phase transition. The predicted dynamic topography is very sensitive to shallow temperature-induced lateral viscosity variations (LVVs) and Rayleigh number. The preservation of LVVs significantly increases the negative dynamic topography at subduction zones. A decrease (or increase) of Rayleigh number increases (or decreases) the predicted present-day dynamic topography considerably. The dynamic topography predicted from the model considering LVVs and with a Rayleigh number of 6 × 108 is most compatible with residual topography models. This Rayleigh number is consistent with the convective vigor of the Earth as supported by generating more realistic lower mantle structure, slab sinking rate, and surface and CMB heat fluxes. The evolution of the surface heat flux pattern is similar to the long-term eustatic sea-level change. Before the formation of Pangea, large negative dynamic topography formed between the plate convergence region of Gondwana and Laurussia. The predicted dynamic topography similar to that of present day has already emerged by about 262 Ma. Powers for degrees 1–3 dynamic topography at 337 Ma and 104 Ma which correspond to times of higher plate velocities and higher surface heat fluxes are larger.

Description: Summary We investigated the induced seismicity, source mechanisms and mechanical responses of a decameter-scale hydraulic stimulation of a pre-existing shear zone in crystalline rock, at the Grimsel Test Site, Switzerland. The analysis reveals the meter-scale complexity of hydraulic stimulation, which remains hidden at the reservoir-scale. High earthquake location accuracy allowed the separation of four distinct clusters, of which three were attributed to the stimulation of fractures in the damage zone of the shear zone. The source mechanism of the larger-magnitude seismicity varied by cluster, and suggests a heterogeneous stress field already prevailing before stimulation, which is further modified during stimulation. In the course of the experiment, stress redistribution led to the aseismic initiation of a tensile-dominated fracture, which induced seismicity in the fourth of the identified seismic clusters. The streaky pattern of seismicity separated by zones without seismicity suggests fluid flow in conduits along existing fracture planes. The observed sub-meter scale complexity questions the forecasting ability of induced seismic hazard at the reservoir scale from small-scale experiments.

Description: Summary Ground observatory and satellite-based determinations of temporal variations in the geomagnetic field probe a decadal to annual time scale range where Earth’s core slow, inertialess convective motions and rapidly propagating, inertia-bearing hydromagnetic waves are in interplay. Here we numerically model and jointly investigate these two important features with the help of a geodynamo simulation that (to date) is the closest to the dynamical regime of Earth’s core. This model also considerably enlarges the scope of a previous asymptotic scaling analysis, which in turn strengthens the relevance of the approach to describe Earth’s core dynamics. Three classes of hydrodynamic and hydromagnetic waves are identified in the model output, all with propagation velocity largely exceeding that of convective advection: axisymmetric, geostrophic Alfvén torsional waves, and non-axisymmetric, quasi-geostrophic Alfvén and Rossby waves. The contribution of these waves to the geomagnetic acceleration amounts to an enrichment and flattening of its energy density spectral profile at decadal time scales, thereby providing a constraint on the extent of the f−4 range observed in the geomagnetic frequency power spectrum. As the model approaches Earth’s core conditions, this spectral broadening arises because the decreasing inertia allows for waves at increasing frequencies. Through non-linear energy transfers with convection underlain by Lorentz stresses, these waves also extract an increasing amount of energy from the underlying convection as their key time scale decreases towards a realistic value. The flow and magnetic acceleration energies carried by waves both linearly increase with the ratio of the magnetic diffusion time scale to the Alfvén time scale, highlighting the dominance of Alfvén waves in the signal and the stabilising control of magnetic dissipation at non-axisymmetric scales. Extrapolation of the results to Earth’s core conditions supports the detectability of Alfvén waves in geomagnetic observations, either as axisymmetric torsional oscillations or through the geomagnetic jerks caused by non-axisymmetric waves. In contrast, Rossby waves appear to be too fast and carry too little magnetic energy to be detectable in geomagnetic acceleration signals of limited spatio-temporal resolution.

Description: Summary The squirt flow model, proposed by Mavko & Jizba, has been widely used in explaining the frequency-related modulus and velocity dispersion between ultrasonic and seismic measurements. In this model, the saturated bulk modulus at high frequency is obtained by taking the so-called unrelaxed frame bulk modulus into Biot's or Gassmann's formula. When using Gassmann's formula, the mineral bulk modulus is taken as matrix bulk modulus. However, the soft pores (cracks) in rocks have a weakening effect on the matrix bulk modulus. The saturated bulk modulus at high frequency calculated with mineral bulk modulus as matrix bulk modulus is higher than the real values. To overcome this shortcoming we propose a modified matrix bulk modulus based on the Betti-Rayleigh reciprocity theorem and non-interaction approximation. This modification takes the weakening effect of soft pores (cracks) into consideration and allows calculating the correct saturated bulk modulus at high frequency under different soft-pore fractions (the ratio of soft porosity to total porosity) or crack densities. We also propose an alternative expression of the modified matrix bulk modulus, which can be directly obtained from laboratory measurements. The numerical results show that the saturated bulk modulus at high frequency using the original matrix bulk modulus (i.e. mineral bulk modulus) is approximated to that using the modified one only for rocks containing a small amount of soft-pore fraction. However, as the soft-pore fraction becomes substantial, using the original bulk matrix modulus is not applicable, but the modified one is still applicable. Furthermore, the results of the modified squirt flow model show good consistency with published numerical and experimental data. The proposed modification extends the applicable range of soft-pore fraction (crack density) of the previous model, and has potential applications in media having a relatively substantial fraction of soft pores or almost only soft pores, such as granite, basalt, and thermally-cracked glasses.

Description: Summary We report 24 palaeomagnetic directions and 10 high-quality Thellier-derived palaeointensity (PI) values, obtained from 27 sites located in Baja California Peninsula, northwestern Mexico. Sampling was done in four rock units (magnesian andesites, calc-alkaline lavas, ignimbrites, adakites) belonging to San Borja and Jaraguay monogenetic volcanic fields. These units have erupted between ∼ 15 and 2.6 Ma (previous K-Ar and 40Ar/39Ar data); hence results are presented in two consecutive periods: middle-late Miocene and Pliocene. The identified main magnetic minerals in the sampled sites are titanomagnetite, magnetite, and minor hematite, of variable grain size, present as intergrowths or surrounding grains, which reflect varying oxidation/reduction conditions during emplacement of high-temperature magmas. Based on previous geological and geophysical records, the kinematic evolution was carefully considered in the region, allowing for the independent restoration of the palaeoposition of each sampled site. Previous palaeodirections were also evaluated and corrected for tectonic motion in order to combine them with present data. Accordingly, a number of 15 and 36 directional data are used to calculate palaeopole position for Pliocene and middle-late Miocene periods, respectively, selected from a total of 74 data points. Pliocene (Plat = 87.8°, Plong = 147.5°, K = 41.06, A95 = 6.0°) and middle-late Miocene (Plat = 86.0°, Plong = 172.7°, K = 41.08, A95 = 3.8) palaeopole positions, calculated after tectonic corrections, are not statistically different from expected North American reference pole. Tectonic correction for Middle-late Miocene virtual geomagnetic poles plays an important role in reducing the resultant tilting from 2.7° to -0.8°. PI mean were calculated for Pliocene and middle-late Miocene periods at 29.2 ± 9.1 μT and 23.2 ± 6.3 μT, respectively. Compiling global filtered PI data, together with our results, indicates that the strength of the geomagnetic field during middle-late Miocene was weak (virtual dipole moment = 5.0 ± 2.2 × 1022 Am2) compared to Pliocene (6.4 ± 2.8 × 1022 Am2), and also relative to the present-day value (7.6 × 1022 Am2). This indicates the global nature of the low dipole moment during the middle-late Miocene period. However, issues related to the spatio-temporal distribution of PI data still present an obstacle to validating these suggestions; therefore, more reliable data are still needed.

Description: Summary We recompute the 26-year weekly Geocenter Motion (GCM) time series from 1994 to 2020 through the network shift approach using Satellite Laser Ranging (SLR) observations to LAGEOS1/2. Then the Singular Spectrum Analysis (SSA) is applied for the first time to separate and investigate the geophysical signals from the GCM time series. The Principal Components (PCs) of the embedded covariance matrix of SSA from the GCM time series are determined based on the w-correlation criterion and two PCs with large w-correlation are regarded as one periodic signal pair. The results indicate that the annual signal in all three coordinate components and semi-annual signal in both X and Z components are detected. The annual signal from this study agrees well in both amplitude and phase with those derived by the Astronomical Institute of the University of Bern and the Center for Space Research, especially for the Y and Z components. Besides, the other periodic signals with the periods of (1043.6, 85, 28), (570, 280, 222.7), and (14.1, 15.3) days are also quantitatively explored for the first time from the GCM time series by using SSA, interpreting the corresponding geophysical and astrodynamical sources of aliasing effects of K1/O1, T2 and Mm tides, draconitic effects, and overlapping effects of the ground-track repeatability of LAGEOS1/2.

Description: Summary The inverse problem of finding the slowness vector from a known ray direction in general anisotropic elastic media is a challenging task, needed in many wave/ray-based methods, in particular, solving two-point ray bending problems. The conventional resolving equation set for general (triclinic) anisotropy consists of two fifth-degree polynomials and a sixth-degree polynomial, resulting in a single physical solution for quasi-compressional (qP) waves and up to 18 physical solutions for quasi-shear waves (qS). For polar anisotropy (transverse isotropy with a tilted symmetry axis), the resolving equations are formulated for the slowness vectors of the coupled qP and qSV waves (quasi-shear waves polarized in the axial symmetry plane), and independently for the decoupled pure shear waves polarized in the normal (to the axis) isotropic plane (SH). The novelty of our approach is the introduction of the geometric constraint that holds for any wave mode in polar anisotropic media: The three vectors—the slowness, ray velocity and medium symmetry axis—are coplanar. Thus, the slowness vector (to be found) can be presented as a linear combination of two unit-length vectors: the polar axis and the ray velocity directions, with two unknown scalar coefficients. The axial energy propagation is considered as a limit case. The problem is formulated as a set of two polynomial equations describing: a) the collinearity of the slowness-related Hamiltonian gradient and the ray velocity direction (third-order polynomial equation), and b) the vanishing Hamiltonian (fourth-order polynomial equation). Such a system has up to twelve real and complex-conjugate solutions, which appear in pairs of the opposite slowness directions. The common additional constraint, that the angle between the slowness and ray directions does not exceed ${90^{ m{o}}}$, cuts off one half of the solutions. We rearrange the two bivariate polynomial equations and the above-mentioned constraint as a single univariate polynomial equation of degree six for qP and qSV waves, where the unknown parameter is the phase angle between the slowness vector and the medium symmetry axis. The slowness magnitude is then computed from the quadratic Christoffel equation, with a clear separation of compressional and shear roots. The final set of slowness solutions consists of a unique real solution for qP wave and one or three real solutions for qSV (due to possible triplications). The indication for a qSV triplication is a negative discriminant of the sixth-order polynomial equation, and this discriminant is computed and analyzed directly in the ray-angle domain. The roots of the governing univariate sixth-order polynomial are computed as eigenvalues of its companion matrix. The slowness of the SH wave is obtained from a separate equation with a unique analytic solution. We first present the resolving equation using the stiffness components, and then show its equivalent forms with the well-known parameterizations: Thomsen, Alkhalifah and ‘weak-anisotropy’. For the Thomsen and Alkhalifah forms, we also consider the (essentially simplified) acoustic approximation for qP waves governed by the quartic polynomials. The proposed method is coordinate-free and can be applied directly in the global Cartesian frame. Numerical examples demonstrate the advantages of the method.

Description: Summary Seismic inversion is one of the most commonly used methods in the oil and gas industry for reservoir characterization from observed seismic data. Deep learning (DL) is emerging as a data-driven approach that can effectively solve the inverse problem. However, existing deep learning-based methods for seismic inversion utilize only seismic data as input, which often leads to poor stability of the inversion results. Besides, it has always been challenging to train a robust network since the real survey has limited labeled data pairs. To partially overcome these issues, we develop a neural network framework with a priori initial model constraint to perform seismic inversion. Our network uses two parts as one input for training. One is the seismic data, and the other is the subsurface background model. The labels for each input are the actual model. The proposed method is performed by log-to-log strategy. The training dataset is firstly generated based on forward modeling. The network is then pre-trained using the synthetic training dataset, which is further validated using synthetic data that has not been used in the training step. After obtaining the pre-trained network, we introduce the transfer learning strategy to fine-tune the pre-trained network using labeled data pairs from a real survey to acquire better inversion results in the real survey. The validity of the proposed framework is demonstrated using synthetic 2D data including both post-stack and pre-stack examples, as well as a real 3D post-stack seismic data set from the western Canadian sedimentary basin.

Description: Summary Probing seismic anisotropy of the lithosphere provides valuable clues on the fabric of rocks. We present a 3-D probabilistic model of shear wave velocity and radial anisotropy of the crust and uppermost mantle of Europe, focusing on the mountain belts of the Alps and Apennines. The model is built from Love and Rayleigh dispersion curves in the period range 5 to 149 s. Data are extracted from seismic ambient noise recorded at 1521 broadband stations, including the AlpArray network. The dispersion curves are first combined in a linearised least squares inversion to obtain 2-D maps of group velocity at each period. Love and Rayleigh maps are then jointly inverted at depth for shear wave velocity and radial anisotropy using a Bayesian Monte-Carlo scheme that accounts for the trade-off between radial anisotropy and horizontal layering. The isotropic part of our model is consistent with previous studies. However, our anisotropy maps differ from previous large scale studies that suggested the presence of significant radial anisotropy everywhere in the European crust and shallow upper mantle. We observe instead that radial anisotropy is mostly localized beneath the Apennines while most of the remaining European crust and shallow upper mantle is isotropic. We attribute this difference to trade-offs between radial anisotropy and thin (hectometric) layering in previous studies based on least-squares inversions and long period data (〉30 s). In contrast, our approach involves a massive dataset of short period measurements and a Bayesian inversion that accounts for thin layering. The positive radial anisotropy (VSH 〉 VSV) observed in the lower crust of the Apennines cannot result from thin layering. We rather attribute it to ductile horizontal flow in response to the recent and present-day extension in the region.

Description: Summary The analysis of surface wave dispersion curves is a way to infer the vertical distribution of shear-wave velocity. The range of applicability is extremely wide: going, for example, from seismological studies to geotechnical characterizations and exploration geophysics. However, the inversion of the dispersion curves is severely ill-posed and only limited efforts have been put in the development of effective regularization strategies. In particular, relatively simple smoothing regularization terms are commonly used, even when this is in contrast with the expected features of the investigated targets. To tackle this problem, stochastic approaches can be utilized, but they are too computationally expensive to be practical, at least, in case of large surveys. Instead, within a deterministic framework, we evaluate the applicability of a regularizer capable of providing reconstructions characterized by tunable levels of sparsity. This adjustable stabilizer is based on the minimum support regularization, applied before on other kinds of geophysical measurements, but never on surface wave data. We demonstrate the effectiveness of this stabilizer on: i) two benchmark—publicly available— datasets at crustal and near-surface scales; ii) an experimental dataset collected on a well-characterized site. In addition, we discuss a possible strategy for the estimation of the depth of investigation. This strategy relies on the integrated sensitivity kernel used for the inversion and calculated for each individual propagation mode. Moreover, we discuss the reliability, and possible caveats, of the direct interpretation of this particular estimation of the depth of investigation, especially in the presence of sharp boundary reconstructions.

Description: Summary The spatial correlation of earthquake ground motion intensity can be measured from strong motion data; however, the data used in past studies is sparsely sampled in space, and only the inter-station distance was considered as a correlation variable. These limitations mean that we have only weak constraints on the true correlation structure of ground motion and that potentially important aspects of spatial correlation are unconstrained. In this study, we combine a large-N seismic array and graph analytics to explore this issue at a local scale using small local and regional earthquakes. Our result suggests site conditions, and how they interact with the incident seismic wavefield, strongly condition the spatial correlation of ground motion. Future progress in characterizing ground motion spatial variability will require dense wavefield measurements, either through nodal deployments, or perhaps distributed acoustic sensing (DAS) measurements, of seismic wavefields. Aftershock sequences of major earthquakes would provide particularly data-rich targets of opportunity.

Description: SUMMARY Mining-induced seismic events can be followed by aftershocks which increase the risk associated with the exploitation. The understanding of the aftershock generation process in induced seismicity may improve post-earthquake safety procedures applied in mines. Rudna copper ore mine in southwestern Poland commonly experiences intense and strong seismic activity accompanying the room-and-pillar exploitation of copper ore. Some strong (magnitude 〉 2) mining events are followed by numerous aftershocks and some are not followed by any. In this study we seek to find whether there is any geological, technological or seismological cause of this diversity. We study 46 strong mining events and focus on their aftershock productivity. We analyze the geological and mining setting of the studied events, their signal similarity, stress drops and the ground motion effect using data from three different seismic networks. Our results show that seismic events producing large aftershock sequences may share similar focal mechanisms and have larger ground effects than events with no aftershocks. The results also indicate the potential differences in stress drops. This interesting observation may help to better evaluate the aftershock hazard in mines. It also indicates the need for a more detailed analysis of the focal mechanisms of strong events and their relationship to the exploitation technique.

Description: Summary The present study introduces a novel method for computing post-seismic crustal internal deformation in a layered earth model. The surface dislocation Love number (DLN) calculated by the reciprocity theorem was implemented as the initial value. Furthermore, numerical integration of the value from the Earth's surface to the interior was undertaken to obtain the internal DLN. This method does not require a combination of the general solution and particular solution for the calculation of internal deformation above the seismic source, thus avoiding the loss of precision. When the post-seismic deformation within a certain period is calculated, the particular solutions at the beginning and end of the considered period cancel each other. This simplifies the calculation of post-seismic internal deformation. The numerical results depict that as the degrees increase, the post-seismic DLN reaches stability in a shorter interval of time. Thus, for improved efficiency of the post-seismic internal deformation calculation, the post-seismic DLNs should be calculated within 2000 degree and integrated with the co-seismic results. As an application, the post-seismic Coulomb failure stress changes (ΔCFS) induced by the 2011 Tohoku-Oki earthquake in the near field around the Japanese archipelagos and two major faults in Northeast China were simulated. The results exhibit that the ΔCFS values in the near field agree well with those simulated by the method in a half-space layered earth model, thus verifying the present method. The co-seismic ΔCFS on the Mishan-Dunhua fault in Northeast China, as an example, is only 0.094–0.668 KPa. However, the ΔCFS caused by the viscoelastic relaxation of the mantle within 5 years following the 2011 Tohoku-Oki event on the same fault exceeds the co-seismic results. Therefore, the cumulative effect of the viscoelastic relaxation of the mantle is deserving of attention.

Description: Summary The piezomagnetic effect is defined as a change in magnetisation with applied stress. Changes in the geomagnetic field caused by the piezomagnetic effect, referred to as the piezomagnetic field, have been theoretically estimated and compared by previous studies to interpret observed variations in the geomagnetic field. However, the piezomagnetic field estimated in previous studies may not provide an accurate estimation because they ignored spatial variations in elasticity, leading to only a rough approximation of the properties of Earth's crust. In this paper, a semi-analytical procedure for calculating the piezomagnetic field arising from a point dislocation source embedded in a layered elastic medium is derived. Following a well-established method of the vector surface harmonic expansion, all of the governing equations written in partial differential equations in a real domain, together with the linear constitutive law of the piezomagnetic effect, are converted to a set of ordinary differential equations in a wavenumber domain. Equations in the wavenumber domain are solved analytically, and each component of the piezomagnetic field in the real domain is obtained after applying the Hankel transform. By using the derived procedure, the piezomagnetic and displacement fields due to a finite fault with strike-slip, dip-slip, and tensile-opening mechanisms are estimated for media with layered elasticity structures. Results for a finite fault are obtained by integrating the point source solution over the fault plane. The results of the numerical analysis allow the effect of heterogeneities in rigidity on the piezomagnetic effect to be examined and implications for the findings of previous investigations to be drawn. In cases where the moment-release at the dislocation source is fixed, the effect of the rigidity differences between upper and lower layers on the piezomagnetic field is minor even in the case where the Curie point depth is near the source of dislocation. This result is in contrast to a previous study that assumed the Mogi model and suggested that heterogeneities in the horizontal direction may be of importance when combined with layered rigidity structures. A contrast is seen between the piezomagnetic and displacement fields corresponding to models with layered rigidity structures: the piezomagnetic field is roughly proportional to the moment-release on a source fault, whereas the displacement field is proportional to slip or opening of the fault. Provided that the rigidity of the crust increases with increasing depth, the calculated piezomagnetic field is likely to have been underestimated in many earlier studies, which assumed uniform rigidity and a geodetically inverted size of slip.

Description: Summary Unravelling relations between lateral variations of mid-crustal seismicity and the geometry of the Main Himalayan Thrust system at depth is a key issue in seismotectonic studies of the Himalayan range. These relations can reveal along strike changes in the behavior of the fault at depth related to fluids or the local ramp-flat geometry and more generally of the stress build-up along the fault. Some of these variations may control the rupture extension of intermediate, large or great earthquakes, the last of which dates back from 1505 CE in far western Nepal. The region is also associated to lateral spatio-temporal variations of the mid-crustal seismicity monitored by the Regional Seismic Network of Surkhet-Birendranagar. This network was supplemented between 2014 and 2016 by 15 temporary stations deployed above the main seismic clusters giving new potential to regional studies. Both absolute and relative locations together with focal mechanisms are determined to gain insight on the fault behavior at depth. We find more than 4000 earthquakes within 5 and 20 km-depth clustered in three belts parallel to the front of the Himalayan range. Finest locations reveal close relationships between seismic clusters and fault segments at depth among which mid-crustal ramps and reactivated tectonic slivers. Our results support a geometry of the Main Himalayan Thrust involving several fault patches at depth separated by ramps and tear faults. This geometry most probably affects the pattern of the coseismic ruptures breaking partially or totally the locked fault zone as well as eventual along strike variations of seismic coupling during interseismic period.

Description: Summary Crustal deformation caused by hydrological processes has long been detected using space geodetic techniques, yet questions remain about the relative contributions of surface water and groundwater to the geodetic signals in different regions. Here, we investigate forward models of elastic loading deformation caused by a variety of water-storage changes within the Great Lakes region, including fluctuations in lake-water volume, soil moisture, and snow load. We use lake-level data from the Great Lakes Environmental Research Laboratory, soil-moisture content from the North American Land Data Assimilation System (NLDAS), snow load from the Snow Data Assimilation System (SNODAS), and background hydrological load at the global scale from Gravity Recovery and Climate Experiment (GRACE). We compare the modeled surface deformation with estimates of hydrological loading deformation inferred from Global Navigation Satellite System (GNSS) measurements. We find that seasonal deformation measured by GNSS is dominated by regional-scale hydrological loading based on strong correlations with the modeled loading displacements. The mean correlation coefficient for the study network is 0.56. The correlation coefficients vary spatially within the study region and exceed 0.9 at some stations near to the Great Lakes. We assess the relative contribution of each individual hydrological component to the total integrated hydrological load. We find that soil moisture consistently explains the largest percentage (27 per cent-69 per cent) of the total vertical loading deformation for 87 per cent of GNSS stations in the Great Lakes region. Snow loading and soil moisture contribute relatively equally in the northern reaches of the study area (e.g. Canadian shield, northern Superior basin). Lake loading accounts for about 10–25 per cent of the total loading signal in the immediate vicinity of the lakes. We also investigate the sensitivities of the surface loading displacements to three different Earth models, including two with lateral variations in structure. The structural variations considered here have limited impact (

Description: Summary Convection in the liquid outer core of the Earth is driven by thermal and chemical perturbations. The main purpose of this study is to examine the impact of double-diffusive convection on magnetic field generation by means of three-dimensional global geodynamo models, in the so-called “top-heavy” regime of double-diffusive convection, when both thermal and compositional background gradients are destabilizing. Using a linear eigensolver, we begin by confirming that, compared to the standard single-diffusive configuration, the onset of convection is facilitated by the addition of a second buoyancy source. We next carry out a systematic parameter survey by performing 79 numerical dynamo simulations. We show that a good agreement between simulated magnetic fields and the geomagnetic field can be attained for any partitioning of the convective input power between its thermal and chemical components. On the contrary, the transition between dipole-dominated and multipolar dynamos is found to strongly depend on the nature of the buoyancy forcing. Classical parameters expected to govern this transition, such as the local Rossby number -a proxy of the ratio of inertial to Coriolis forces- or the degree of equatorial symmetry of the flow, fail to capture the dipole breakdown. A scale-dependent analysis of the force balance instead reveals that the transition occurs when the ratio of inertial to Lorentz forces at the dominant length scale reaches 0.5, regardless of the partitioning of the buoyancy power. The ratio of integrated kinetic to magnetic energy Ek/Em provides a reasonable proxy of this force ratio. Given that Ek/Em ≈ 10−4 − 10−3 in the Earth’s core, the geodynamo is expected to operate far from the dipole-multipole transition. It hence appears that the occurrence of geomagnetic reversals is unlikely related to dramatic and punctual changes of the amplitude of inertial forces in the Earth’s core, and that another mechanism must be sought.

Description: Summary The Ordos Block is located at the intersection of the Tibetan Plateau, the South China Plate, and the North China Craton. The deep deformation of its surrounding areas is geologically complicated. Seismic anisotropy provide information about past and present deformation of the upper mantle and help to better understand deep deformation processes. We obtained the anisotropic pattern at high resolution within and surrounding the Ordos Block by analyzing teleseismic data from ∼710 newly deployed temporary seismic stations based on SKS phases from ∼86 earthquake. The central and eastern Ordos Block with a relatively thick and rigid lithosphere is characterized by a weak anisotropy. In the western part of the Ordos Block, the fast-wave polarization direction is dominantly NW–SE. We believe the lithosphere in the western part of the Ordos Block may have undergone significant deformation caused by expansion and compression of the Tibetan Plateau. Near the Datong Basin, the delay time is 0.92 s and the fast-wave polarization direction is mainly oriented NW–SE, perpendicular to the NE-directed compression of the Tibetan Plateau and parallel to the basin extension direction. We speculate that this anisotropy is related to the long-range effect of NE compression of the northeastern Tibetan Plateau on the low-velocity upper mantle in this area. The weak anisotropy in the central part of the Shanxi Rift indicates that the thickness and mechanical strength of the high-velocity lithosphere are higher than in the southern and northern regions.

Description: Summary Correlation of the coda of Empirical Green's functions from ambient noise can be used to reconstruct Empirical Green's function between two seismic stations deployed different periods of time. However, such method requires a number of source stations deployed in the area surrounding a pair of asynchronous stations, which limit its applicability in cases where there are not so many available source stations. Here, we propose an alternative method, called two-station C2 method, which uses one single station as a virtual source to retrieve surface wave phase velocities between a pair of asynchronous stations. Using ambient noise data from USArray as an example, we obtain the interstation C2 functions using our C2 method and the traditional cross-correlation functions (C1 functions). We compare the differences between the C1 and C2 functions in waveforms, dispersion measurements, and phase velocity maps. Our results show that our C2 method can obtain reliable interstation phase velocity measurements, which can be used in tomography to obtain reliable phase velocity maps. Our method can significantly improve ray path coverage from asynchronous seismic arrays and enhance the resolution in ambient noise tomography for areas between asynchronous seismic arrays.

Description: Summary We calculate and analyse the coordinate time series of 282 permanent GPS stations located in Greece and 47 in surrounding countries. The studied period is 2000–2020. The average GPS time series length is 6.5 years. The formal velocity uncertainties are rescaled to be consistent with the velocity scatters measured at 110 pairs of stations separated by less 15 km. We remove the effect of the crustal earthquakes of Mw ≥ 5.3. We quantify and model the postseismic deformations. Two relaxation times are usually needed, one short of some weeks, one long of one year or more. For the large Mw = 6.9 events of Samothraki 2014 and Methoni 2008, the postseismic deformation equals or exceeds the coseismic one. We detect at three stations a deformation transient in May 2018 that may correspond to a slow earthquake beneath Zakynthos and north-west Peloponnese, with equivalent magnitude 5.8. The density and accuracy of the velocities make it possible to better quantify several characteristics of the deformation in the Aegean, in particular: (a) the transition from the Anatolian domain, located in the south-east, to the European domain through the western end of the North Anatolian fault; (b) the north-south extension in the western Aegean; (c) the east-west extension of the western Peloponnese; (d) the clockwise rotation of the Pindos; (e) the north-south extension in central Macedonia. Large parts of the central Aegean, eastern Peloponnese and western Crete form a wide stable domain with internal deformation below 2 nstrain yr−1. We build a kinematic model comprising ten crustal blocks corresponding to areas where the velocities present homogeneous gradients. The block boundaries are set to fit with known localized deformation zones, e.g. the rift of Corinth, the North Anatolian fault, the Katouna fault. When the velocities steps are clear but not localized, e.g. through the Peloponnese, the boundary line is arbitrary and represents the transition zone. The model fits the velocities with a root mean square deviation of ± 0.9 mm yr−1. At the boundaries between blocks we compare the predicted and observed deformations. We find shear rates of 7.4 and 9.0 mm yr−1 along the Movri and Katouna faults, 14.9 and 8.7 mm yr−1 along the North Anatolian fault near Lemnos and near Skopelos respectively, extension of 7.6, 1.5 and 12.6 mm yr−1 across the Gulf of Patras, the Trichonis Lake and the Ambracian Gulf. The compression across western Epirus is 3.7 mm yr−1. There is a dextral transtensional movement of 4.5 mm yr−1 between the Amorgos and Astypalea islands. Only the Ionian Islands region shows evidence of coupling along the subduction interface.

Description: SUMMARY Ekman layers develop at the boundaries of the Earth’s fluid core in response to precession. Instabilities in these layers lead to turbulence when a local Reynolds number, Re, based on the thickness of the Ekman layer, exceeds a critical value. The transition to turbulence is often assessed using experiments for steady Ekman layers, where the interior geostrophic flow is independent of time. Precessionally driven flow varies on diurnal timescales, so the transition to turbulence may occur at a different value of Re. We use 3-D numerical calculations in a local Cartesian geometry to assess the transition to turbulence in precessional flow. Calculations retain the horizontal component of the rotation vector and account for the influence of fluid stratification. The transition to turbulence in a neutrally stratified fluid occurs near Re = 500, which is higher than the value Re = 150 usually cited for steady Ekman layers. However, it is comparable to the nominal value for precessional flow in the Earth. Complications due to fluid stratification or a magnetic field can suppress the transition to turbulence, reducing the likelihood of turbulent Ekman layers in the Earth’s core.

Description: SUMMARY An underground nuclear explosion (UNE) couples mechanical energy into crustal rock, which propagates as seismic and acoustic waves. These different physical phenomena transport, by different pathways, to standoff detectors at varying distances. The transport pathways attenuate the original signal but in different ways. Enabled by correct statistical weighting, signal attenuation models can be used to combine these disparate sensor data to estimate the yield of an UNE. Contemporaneous statistical models, used in yield estimation, can be improved with an advanced partition of error for these physical signal propagation models. We present an advanced multivariate approach to error modelling of multiphenomenology physical signatures. In addition to measurement error, our error model represents physical model biases as random with a physics-based covariance structure. To illustrate this proposed framework, we demonstrate the estimation of explosion yield using openly available seismic and acoustic data from chemical single-point explosions.

Description: SUMMARY Following the installation of a temporary seismological network in western Greece north of the Gulf of Patras, we determined the quality of the sites of each of the 10 stations in the network. For this, we used the horizontal-to-vertical spectral ratio (HVSR) method and calculated an average curve over randomly selected days between 0 and 10 Hz. The daily HVSR curve is determined by the HVSR 12-hr calculation (1 hr every two) without distinction between seismic ambient noise and earthquake signal. The HVSR curves obtained can be classified in three categories: flat curves without amplification, curves with a amplification peaks covering a large frequency range, and curves with one or more narrow peaks. In this third category C3, one station has one peak, two have two and one has three. On the contrary of what it is commonly assumed, the amplitudes and the resonance frequencies of these narrow peaks are not stable over time in C3. We determined the maximum of the amplitude of each peak with the corresponding central frequency for each day during 2.5 yr. Except for the station with three peaks, which finally appears stable within the uncertainties, the principal peak exhibits a seasonal variation, with a maximum in winter and a minimum in summer, the observations being more dispersed during winter. The second peak, when it exists, varies in the same way except at one station where it varies oppositely. These variations are clearly correlated with the loading and unloading cycle of the underlying aquifers as shown by the comparison with water level and yield measurements from wells located close to the stations. Moreover, they are also correlated with the vertical surface displacements observed at continuously recording GPS stations. The dispersion of the observed maximum amplitude in winter is probably related to the rainfall and the soil moisture modifying the S-wave velocity as revealed by other studies. From this study, we would like to emphasize that the use the HVSR method to constrain the S-wave velocity and the thickness of the sediment layer over the bedrock in the basin, has to be done with caution. Upon further confirmation of its robustness, the HVSR methodology presented here could be a good and easy-to-use tool for a qualitative survey of the aquifer backdrop and its seasonal behaviour, and of the soil moisture conditions.

Description: Summary Squirt flow is an essential cause of wave dispersion and attenuation in saturated rocks. The squirt flow model, proposed by Gurevich et al. (2010), has been widely applied to explain the wave dispersion and associated attenuation for saturated rocks at sonic and seismic frequency bands. In this model, the saturated bulk modulus is obtained by taking the partially relaxed frame bulk modulus as the dry frame modulus into Gassmann's formula with the mineral bulk modulus as the matrix bulk modulus. However, because of the weakening effect of soft pores on rock matrix bulk modulus, the model cannot accurately predict the saturated bulk modulus when the soft-pore fraction (the ratio of the soft porosity to total porosity) becomes large. We modified this model following Gurevich et al. (2010) by setting a different boundary condition. The modified squirt flow model can obtain correct saturated bulk modulus for large soft-pore fractions in the full range of frequencies, showing excellent consistency with the predictions of Gassmann and Mavko & Jizba (modified) at both low- and high-frequency limits, respectively. Modeling results show that the saturated bulk moduli and their dispersions calculated by the original and modified models exhibit little difference when the soft-pore fraction is small. Under this condition, the original model is as effective and accurate as the modified one. When the soft-pore fraction becomes larger, the differences in the bulk moduli and their dispersions become substantial, suggesting the original model is not applicable any longer. Furthermore, the differences calculated for the intermediate frequency range is even more obvious than other ranges, suggesting that the modified model should be used to calculate the bulk modulus and the dispersion in this frequency range. In summary, the modified squirt flow model can extend the original model's applicable range in terms of soft-pore fraction and has a potential application in rocks having a relatively large amount of soft-pore fraction such as basalts.

Description: Summary The noise attenuation is a fundamental step in seismic data processing, especially when ground-roll suppression remains a challenge. Rank-reduction methods have become quite popular in recent decades, as they promote significant improvements in the quality of data, highlighting reflections in seismograms. We present a methodology for ground-roll filtering, which combines the application of a recursive-iterative Singular Spectrum Analysis method, in the time-domain, as a particular way to decompose seismic data, with the computation of the average instantaneous frequency of the signal components. This combination allows for a precise estimation and filtering of the ground-roll noise. The frequency values are used for determining, in each component, the low-frequency parts associated with the ground-roll. For every single component, the ground-roll is attenuated by zeroing, and stacking the data-components, where the average instantaneous frequency values match the ground-roll bandwidth of frequency. Also, in order to enhance the lateral coherence of the reflectors, we present an extension of the recursive-iterative algorithm for a multichannel case. The multichannel algorithm is applicable on a shot, or common mid-point family of seismic traces, after the normal move-out correction. The numerical results using real data show the effectiveness of the proposed methodology for ground-roll attenuation and for improving the velocity analysis.

Description: SUMMARY The costly power requirements of delivering seismic data back to Earth from planetary missions requires the development of algorithms for lander-side signal analysis for telemetry prioritization. This is difficult to explicitly program, especially if no prior seismic data are available from the planetary body. Deep learning computer vision has been used to generalize seismic signals on Earth for earthquake early warning problems but such techniques have not yet been expanded to planetary science. We demonstrate that Convolutional Neural Networks can be used to accurately catalogue planetary seismicity without local training data by building binary noise/signal classifiers from a single Earth seismic station and applying the models to moonquakes from the Apollo Passive Seismic Experiment (PSE) and the Lunar Seismic Profiling Experiment (LSPE). In order to promote generality and reduce the amount of training data, the algorithms use spectral images instead of time-series. Two- to five-layer convolution models are tested against a subset of 200 Grade-A events from the PSE and obtained station accuracy averages of 89–96 per cent. As the model was applied to an hour trace of data (30 min before and after the Grade-A event), additional detections besides the Grade-A event are unavoidable. In order to comprehensively address algorithm accuracy, additional seismic detections corresponding to valid signals such as other moonquakes or multiples within a particularly long event needed to be compared with those caused by algorithm error or instrument glitches. We developed an ‘extra-arrival accuracy’ metric to quantify how many of the additional detections were due to valid seismic events and used it to select the three-layer model as the best fit. The three-layer model was applied to the entire LSPE record and matched the lunar day–night cycle driving thermal moonquake generation with fewer false detections than a recent study using Hidden Markov Models. We anticipate that these methods for lander-side signal detection can be easily expanded to non-seismological data and may provide even stronger results when supplemented with synthetic training data.

Description: SUMMARY The challenge of ruling out potential rupture nodal planes with opposite dip orientations during interferometric synthetic aperture radar (InSAR)-based kinematic inversions has been widely reported. Typically, slip on two or more different fault planes can match the surface deformation measurements equally well. The ambiguous choice of the nodal plane for the InSAR-based models was thought to be caused by InSAR's 1-D measurement and polar orbiting direction, leading to its poor sensitivity to north–south crustal motion. Through synthetic experiments and simulations, this paper quantitatively demonstrates the main reason of the ambiguous InSAR-based models, which confuse researchers in the small-to-moderate thrust earthquake cases investigation. We propose the inherent 1-D measurement is not the principle cause of the fault plane ambiguity, since models derived from the same InSAR data predict similar, but not identical, 3-D deformation patterns. They key to differentiating between these different models is to be able to resolve the small asymmetry in the surface deformation pattern, which may be smaller in amplitude than the typical noise levels in InSAR measurements. We investigate the fault geometry resolvability when using InSAR data with different noise levels through ‘R’ value. We find that the resolvability does not only rely on the InSAR noise, but also on the fault geometry itself (i.e. depth, dips angle and strike). Our result shows that it is impossible to uniquely determine the dip orientation of thrust earthquakes with Mw  5.0 km with InSAR data at a noise level that is typical for mountain belts. This inference is independent from the specific data set (i.e. interferogram or time-series) and allows one to assess if one can expect to be able to resolve the correct fault plane at all.

Description: Summary Secondary microseisms are ubiquitous ambient noise vibrations due to ocean activity, dominating worldwide seismographic records at seismic periods between 3 and 10 s. Their origin is a heterogeneous distribution of pressure fluctuations along the ocean surface. In spherically symmetric Earth models, no Love surface waves are generated by such a distributed surface source. We present global-scale modeling of three-component secondary microseisms using a spectral-element method, which naturally accounts for a realistic distribution of surface sources, topography and bathymetry, and three-dimensional (3D) heterogeneity in Earth’s crust and mantle. Seismic Love waves emerge naturally once the system reaches steady state. The ergodic origin of Love waves allows us to model the horizontal components of secondary microseisms for the first time. Love waves mostly originate from the interaction of the seismic wavefield with heterogeneous earth structure, in which the mantle plays an important role despite the short periods involved. Bathymetry beneath the source region produces weak horizontal forces that are responsible for a weak and diffuse Love wavefield. The effect of the bathymetric force splitting in radial and horizontal components is overall negligible when compared to the effect of 3D heterogeneity. However, we observe small and well-focused Love-wave arrivals at seismographic stations in Europe due to force splitting at the steepest portion of the North Atlantic Ridge and the ocean-continent boundary. The location of the sources of Love waves is seasonal at periods shorter than about 7 s, while seasonality is lost at the longer periods. Sources of Rayleigh and Love waves due to the same storm may be located very far away, indicating that energy equipartitioning might not hold in the secondary microseism period band.

Description: Summary The seismic velocity of the oceanic crust is a function of its physical properties that include its lithology, degree of alteration, and porosity. Variations in these properties are particularly significant in young crust, but also occur with age as it evolves through hydrothermal circulation and is progressively covered with sediment. While such variation may be investigated through P-wave velocity alone, joint analysis with S-wave velocity allows the determination of Poisson's ratio, which provides a more robust insight into the nature of change in these properties. Here we describe the independent modelling of P- and S-wave seismic datasets, acquired along an ∼330 km-long profile traversing new to ∼8 Myr-old oceanic crust formed at the intermediate-spreading Costa Rica Rift (CRR). Despite S-wave data coverage being almost four-times lower than that of the P-wave dataset, both velocity models demonstrate correlations in local variability and a long-wavelength increase in velocity with distance, and thus age, from the ridge axis of up to 0.8 and 0.6 km s−1, respectively. Using the Vp and Vs models to calculate Poisson's ratio (σ), it reveals a typical structure for young oceanic crust, with generally high values in the uppermost crust that decrease to a minimum of 0.24 by 1.0–1.5 km sub-basement, before increasing again throughout the lower crust. The observed upper crustal decrease in σ most likely results from sealing of fractures, which is supported by observations of a significant decrease in porosity with depth (from ∼15 to 

Description: Summary Previous studies have demonstrated that finite-fault simulations of actual or hypothetical earthquakes using deterministic, physics-based simulation techniques constitute an effective tool for characterizing near-fault ground strains and rotations in the low-frequency range. The characteristics of these motions are further investigated in this study by performing forward ground-motion simulations of three well-documented strike-slip earthquakes (i.e. 2004 Mw 6.0 Parkfield, 1979 Mw 6.5 Imperial Valley, 1999 Mw 7.5 Izmit) using models of the seismic source and crustal structure available in the literature. Time histories of ground strains and rotations are numerically generated at near-fault stations and at a dense grid of observation points extending over the causative fault. This is achieved by finite differencing translational motions simulated at very closely spaced stations using a kinematic modeling approach. The simulation results show that the three strike-slip earthquakes produce large-amplitude pulse-like shear strain and torsion in the forward direction of rupture propagation. The time histories of specific components of displacement gradient, strain, and rotation at near-fault stations can be estimated from those of ground velocities using a phase velocity, whereas peak ground torsions in the near-fault region can be reasonably estimated from peak horizontal ground velocities using a scaling factor. However, both the phase velocity and the scaling factor exhibit significant variability in the near-fault region of the considered earthquakes. The concept of isochrones is also utilized to associate fault rupture characteristics with near-fault ground strains and rotations. The results indicate that the seismic energy radiated from the high-isochrone-velocity region of the fault—which encompasses areas of large slip locally driven by high stress drop—arrives at a near-fault station in a short time interval that coincides with the time window of the large-amplitude pulse-like shear strain and torsion.

Description: Summary Elongate inclusions immersed in a viscous fluid generally rotate at a rate that is different from the local angular velocity of the flow. Often, a net alignment of the inclusions develops, and the resulting shape preferred orientation (SPO) of the particle ensemble can then be used as a strain marker that allows reconstruction of the fluid’s velocity field. Much of the previous work on the dynamics of flow-induced particle rotations has focused on spatially homogeneous flows with large-scale tectonic deformations as the main application. Recently, the theory has been extended to spatially varying flows, such as magma with embedded crystals moving through a volcanic plumbing system. Additionally, an evolution equation has been introduced for the probability density function (PDF) of crystal orientations. Here, we apply this new theory to a number of simple, two-dimensional flow geometries commonly encountered in magmatic intrusions, such as flow from a dyke into a reservoir or from a reservoir into a dyke, flow inside an inflating or deflating reservoir, flow in a dyke with a sharp bend, and thermal convection in a magma chamber. The main purpose is to provide a guide for interpreting field observations and for setting up more complex flow models with embedded crystals. As a general rule, we find that a larger aspect ratio of the embedded crystals causes a more coherent alignment of the crystals, while it has only a minor effect on the geometry of the alignment pattern. Due to various perturbations in the crystal rotation equations that are expected in natural systems, we show that the time-periodic behavior found in idealized systems is probably short-lived in nature, and the crystal alignment is well described by the time-averaged solution. We also confirm some earlier findings. For example, near channel walls, fluid flow often follows the bounding surface and the resulting simple shear flow causes preferred crystal orientations that are approximately parallel to the boundary. Where pure shear deformation dominates, there is a tendency for crystals to orient themselves in the direction of the greatest tensile strain rate. Where flow impinges on a boundary, for example in an inflating magma chamber or as part of a thermal convection pattern, the stretching component of pure shear aligns with the boundary, and the crystals orient themselves in that direction. In the field, this local pattern may be difficult to distinguish from a boundary-parallel simple shear flow. Pure shear also dominates along the walls of a deflating magma chamber and in places where the flow turns away from the reservoir walls, but in these locations, the preferred crystal orientation is perpendicular to the wall. Overall, we find that our calculated patterns of crystal orientations agree well with results from analogue experiments where similar geometries are available.

Description: Summary A new global model of the present-day thermochemical state of the lithosphere and upper mantle based on global waveform inversion, satellite gravity and gradiometry measurements, surface elevation and heat flow data has been recently presented: WINTERC-G (Fullea et al., 2021). WINTERC-G is built within an integrated geophysical-petrological framework where the mantle seismic velocity and density fields are computed in a thermodynamically self-consistent framework, allowing for a direct parameterisation in terms of the temperature, pressure and composition of the subsurface rocks. In this paper, we combine WINTERC-G thermal and compositional fields along with laboratory experiments constraining the electrical conductivity of mantle minerals, melt and water, and derive a set of new global three dimensional electrical conductivity models of the upper mantle. The new conductivity models, WINTERC-e, consist of two end-members corresponding to minimum and maximum conductivity of the in-situ rock aggregate accounting for mantle melting, mineral water content and the individual conductivities of the main stable mantle mineral phases. The end-member models are validated over oceans by simulating the magnetic field induced by the ocean M2 tidal currents and comparing the predicted fields with the M2 tidal magnetic field estimated from six-year Swarm satellite observations. Our new conductivity model, derived independently from any surface or satellite magnetic data sets, is however able to predict tidal magnetic fields that are in good agreement with the Swarm M2 tidal magnetic field models estimated by Sabaka et al. (2018, 2020) and Grayver & Olsen (2019). Our predicted M2 tidal magnetic fields differ in amplitudes by about 5-20% from the Swarm M2 tidal magnetic field, with the high conductivity WINTERC-e end-member model accounting for mantle melt and water content capturing the structure of Swarm data better than the low conductivity end-member model. Spherically symmetric conductivity models derived by averaging new WINTERC-e conductivities over oceanic areas are slightly more conductive than the recent global conductivity models AA17 by Grayver et al. (2017) derived from Swarm and CHAMP satellite data in the 60-140 km depth range, while they are less conductive deeper in the mantle. The conductivities in WINTERC-e are about 3-4 times smaller than the AA17 conductivities at a depth of 400 km. Despite the differences in electrical conductivity, our spherically symmetric high conductivity end-member model WINTERC-e captures the structure of Swarm M2 tidal magnetic field almost the same as a state of the art 1D conductivity models derived entirely from magnetic data (AA17, (Grayver et al., 2017). Moreover, we show that realistic lateral electrical conductivity inhomogeneities of the oceanic upper mantle derived from the temperature, melt and water distributions in WINTERC-e contribute to the M2 tidal magnetic field up to ±0.3 nT at 430 km altitude.

Description: Summary We examine the potential of frequency-dependent Rayleigh wave ellipticity, derived from microtremors, for the investigation of heterogeneous subsurface structure. Based on numerical simulation, we analyze the effects of interference waves in microtremors, primarily the various propagation directions of the Rayleigh waves, linear polarization waves, and white noise, on the ellipticity frequency-dependent estimation of the Rayleigh waves. A data processing scheme to separate the Rayleigh waves from the interference waves is proposed and verified by synthetic data. We performed a field experiment in the mountainous areas of Southwest China to show that the ellipticity frequency-dependency of Rayleigh waves in the period range of 0.05 to 5 s can be estimated from the microtremor records with the proposed data processing scheme. In addition, the method is feasible for investigating lateral heterogeneity within the top several hundred meters in the mountain regions. The study also reveals that the features of the ellipticity anomaly of a local heterogeneity are related to the propagation directions of the Rayleigh waves, and to reduce the ambiguity of the anomaly, the propagation direction of the waves picked for the ellipticity estimation should be consistent with (along or opposite to) that of the survey line. Then, to eliminate the effects of the phase differences due to the propagation direction, or time, the ellipticity for each location should be estimated by a single event rather than multiple events from the derived Rayleigh wave arrivals.

Description: Summary The differential effective medium (DEM) theory studied in this article describes elastic moduli of a fractured medium with help of differential equations, where crack density is the independent variable and fluid saturation is a parameter. The effective medium is isotropic for randomly oriented flat ellipsoidal cracks and thus fully characterized by two elastic constants. In this article we derive an analytical solution of the equation for Poisson’s ratio and we transform the differential equation for Young’s modulus into a non-linear algebraic equation. Fluid saturation and crack density can then be determined from measured wave propagation velocities by a simple algorithm. We also derive approximate solutions for elastic moduli as a function of crack density and saturation, which allows to quantify the uncertainty of the result due to measurement errors. The DEM theory leads to higher crack densities than the self-consistent (SC) method and to lower crack densities than the non-interacting (NI) theory for measured elastic moduli, while all three methods give similar fluid saturation fractions. As an example of application of our theoretical results, we study weathered granite in the Strengbach water catchment in the Vosges mountains in France. We have performed full waveform sonic logging measurements in an 86 m deep borehole located at an altitude of 1130 m above sea level, which is used for hydro-geophysical and geochemical studies of a granitic aquifer. The logging data allows us to investigate P and S waves in the depth range between 40 and 80 m. The P and S wave propagation velocities take average values of 5.0 km/s and 2.7 km/s, respectively, with the highest values of 5.8 km/s and 3.2 km/s at 75–80 m depth. From these velocities we obtain a water saturation of 75 ± 25 per cent. The crack density describes the degree of weathering of the granite, which generally decreases with depth, but takes high values near layers of strongly weathered granite. Crack density is on average 0.5, with the highest value of 1.0 at 65 m and the lowest value of 0.2 at 75 m depth. The analysis of the full waveform logging data by the DEM method supports results from previous geochemical and hydrological studies in the Strengbach catchment which concluded that water is stored in deeper layers of the granitic aquifer.

Description: Summary The ∼M6 1799 Bouin earthquake is considered as one of the largest earthquakes to have struck Western France. However, the seismogenic source potentially responsible for this event remain marginally documented. We present results from a focused offshore-onshore multidisciplinary survey in its meizoseismal area in order to identify the fault segments that potentially ruptured during this earthquake. Based on macroseismic data and the geology, we focused our study on the so-called Machecoul Fault as a potential source of the 1799 Bouin event. Our survey includes extensive high-resolution seismic reflection, high resolution bathymetry and a one-year seismological survey. These data were combined with existing topography, onshore gravity data and drill data to document the geometry of the Marais Breton / Baie de Bourgneuf basin, the past tectonic activity and the current local microearthquakes at depth along its bounding faults. Offshore and onshore observations suggest a recent activity of the segmented Machecoul Fault bounding the basin to the North. Offshore, the planar contact between the Plio-Quaternary sediments and the basement along the fault trace as well as the thickening of these sedimentary units near this contact suggests tectonic control rather than erosion. Onshore, the recent incision of the footwall of the fault suggests a recent tectonic activity. The temporary local seismological experiment deployed between 2016 and 2017 recorded a diffuse micro-seismicity down to the depth of 22 +/-5 km along the southward dipping Machecoul Fault, associated with predominantly normal fault mechanisms. Altogether, these results suggest that the Machecoul Fault is a serious candidate for being the source of the historical Bouin 1799 earthquake.

Description: SUMMARY The timing of inner core nucleation is a hugely significant event in Earth's evolution and has been the subject of intense debate. Some of the most recent theoretical estimates for the age of nucleation fall throughout the Neoproterozoic era; much younger than previously thought. A young inner core requires faster recent core cooling rates and a likely hotter early core; knowledge of its age would be invaluable in understanding Earth's thermal history and total energy budget. Predictions generated by numerical dynamo models need to be tested against such data, but records are currently much too sparse to constrain the event to a precise period of time. Here, we present results from 720 Ma dolerite dykes (and one sill) from the Franklin Large Igneous Province, which fall within a crucial 300 Myr gap in palaeointensity records. This study uses three independent techniques on whole rocks from 11 sites spread across High Arctic Canada and Greenland to produce virtual dipole moments ranging from 5 to 20 ZAm2 (mean 11 ZAm2); almost one order of magnitude lower than the present-day field. These weak-field results agree with recent ultralow palaeointensity data obtained from Ediacaran rocks formed ∼150 Myr later and may support that the dynamo was on the brink of collapse in the Neoproterozoic prior to a young inner core formation date.

Description: Summary The classical Backus-Gilbert method seeks localized Earth-structure averages at the shortest length scales possible, given a dataset, data errors, and a threshold for acceptable model errors. The resolving length at a point is the width of the local averaging kernel, and the optimal averaging kernel is the narrowest one such that the model error is below a specified level. This approach is well suited for seismic tomography, which maps three-dimensional Earth structure using large sets of seismic measurements. The continual measurement-error decreases and data-redundancy increases have reduced the impact of random errors on tomographic models. Systematic errors, however, are resistant to data redundancy and their effect on the model is difficult to predict. Here, we develop a method for finding the optimal resolving length at every point, implementing it for surface-wave tomography. As in the Backus-Gilbert method, every solution at a point results from an entire-system inversion, and the model error is reduced by increasing the model-parameter averaging. The key advantage of our method stems from its direct, empirical evaluation of the posterior model error at a point. We first measure inter-station phase velocities at simultaneously recording station pairs and compute phase-velocity maps at densely, logarithmically spaced periods. Numerous versions of the maps with varying smoothness are then computed, ranging from very rough to very smooth. Phase-velocity curves extracted from the maps at every point can be inverted for shear-velocity (VS) profiles. As we show, errors in these phase-velocity curves increase nearly monotonically with the map roughness. We evaluate the error by isolating the roughness of the phase-velocity curve that cannot be explained by any Earth structure and determine the optimal resolving length at a point such that the error of the local phase-velocity curve is below a threshold. A 3D VS model is then computed by the inversion of the composite phase-velocity maps with an optimal resolution at every point. The estimated optimal resolution shows smooth lateral variations, confirming the robustness of the procedure. Importantly, the optimal resolving length does not scale with the density of the data coverage: some of the best-sampled locations display relatively low lateral resolution, probably due to systematic errors in the data. We apply the method to image the lithosphere and underlying mantle beneath Ireland and Britain. Our very large dataset was created using new data from Ireland Array, the Irish National Seismic Network, the UK Seismograph Network, and other deployments. A total of 11238 inter-station dispersion curves, spanning a very broad total period range (4–500 s), yield unprecedented data coverage of the area and provide fine regional resolution from the crust to the deep asthenosphere. The lateral resolution of the 3D model is computed explicitly and varies from 39 km in central Ireland to over 800 km at the edges of the area, where the data coverage declines. Our tomography reveals pronounced, previously unknown variations in the lithospheric thickness beneath Ireland and Britain, with implications for their Caledonian assembly and for the mechanisms of the British Tertiary Igneous Province magmatism.

Description: Summary Lamellar magnetism is a source of remanent magnetization in natural rocks different from common bulk magnetic moments in ferrimagnetic minerals. It has been found to be a source for a wide class of magnetic anomalies with extremely high Koenigsberger ratio. Its physical origin are uncompensated interface moments in contact layers of nanoscale ilmenite lamellae inside an hematite host, which also generate unusual low-temperature (low-T) magnetic properties, such as shifted low-T hysteresis loops due to exchange bias. The atomic-magnetic basis for the exchange bias discovered in the hematite-ilmenite system is explored in a series of articles. In this third article of the series, simple models are developed for lamellae interactions of different structures when samples are either cooled in zero-field, or field-cooled in 5 T to temperatures below the ordering temperature of ilmenite. These models are built on the low-temperature measurements described earlier in Paper II. The important observations include: a) the effects of lamellar shapes on magnetic coupling, b) the high-T acquisition of lamellar magnetism and low-T acquisition of magnetization of ilmenite lamellae, c) the intensity of lamellar magnetism and the consequent ilmenite magnetism in populations of randomly oriented crystals, d) lattice-preferred orientation of the titanohematite host crystal populations, and e) the effects of magnetic domain walls in the host on hysteresis properties. Based on exemplary growth models of lamellae with different geometries and surface couplings we here provide simple models to assess and explain the different observations listed above. Already the simplified models show that the shapes of the edges of ilmenite lamellae against their hematite hosts can control the degree of low-T coupling between ilmenite, and the lamellar magnetic moments. The models also explain certain features of the low-T exchange bias in the natural samples and emphasize the role of lattice-preferred orientation upon the intensity of remanent magnetization. The inverse link between ilmenite remanence and exchange-bias shift in bimodal low-T ilmenite lamellae is related to different densities of hematite domain walls induced by the clusters of ilmenite lamellae.

Description: Summary Some geological configurations, like sedimentary basins, are prone to site effects. Basins are often composed of different geological layers whose properties are generally considered as spatially homogeneous or smoothly varying. In this study, we address the influence of small-scale velocity fluctuations on seismic response. For this purpose, we use the Spectral Element method to model the 2D SH wave propagation on a basin of 1.1 km long and ≈ 60 m deep, representing a 2D profile in the city of Nice, France. The velocity fluctuations are modeled statistically as a random process characterized by a Von Karman autocorrelation function and are superimposed to the deterministic model. We assess the influence of the amplitude and correlation length of the random velocities on the surface ground motion. We vary the autocorrelation function’s parameters and compute seismic wavefields in 10 random realizations of the stochastic models. The analyses of our results focus on the Envelope and Phase differences between the waveforms computed in the random and deterministic models; on the variability of ground motion intensity measures, such as the peak ground velocity (PGV), the pseudo-spectral acceleration response (PSA); and the 2D basin response (transfer function). We find that the amplitude of fluctuations has a greater effect on the ground motion variability than the correlation length. Depending on the random medium realization, the ground motion in one stochastic model can be locally amplified or de-amplified with respect to the reference model due to the presence of high or low velocity contrasts, respectively. When computing the mean amplification of different random realizations the results may be smaller than those of the reference media due to the smoothing effect of the average. This study highlights the importance of knowing the site properties at different scales, particularly at small scales, for proper seismic hazard assessment.

Description: Summary Joint inversion of multiphysics data is a practical approach to the integration of geophysical data, which produces models of reduced uncertainty and improved resolution. The development of effective methods of joint inversion requires considering different resolutions of different geophysical methods. This paper presents a new framework of joint inversion of multiphysics data, which is based on a novel formulation of Gramian constraints and mitigates the difference in resolution capabilities of different geophysical methods. Our approach enforces structural similarity between different model parameters through minimizing a structural Gramian term, and it also balances the different resolutions of geophysical methods using a multiscale resampling strategy. The effectiveness of the proposed method is demonstrated by synthetic model study of joint inversion of the P-wave traveltime and gravity data. We apply a novel method based on Gramian constraints and multiscale resampling to jointly invert the gravity and seismic data collected in Yellowstone national Park to image the crustal magmatic system of the Yellowstone. Our results helped to produce a consistent image of the crustal magmatic system of the Yellowstone expressed both in low-density and low-velocity anomaly just beneath the Yellowstone caldera.

Description: Summary Detection of seismic events at or below the noise level is enabled by the use of dense arrays of receivers and corresponding advances in data analysis methods. It is not only important to detect tectonic events, but also events from man-made, non-earthquake sources and events that originate from coupling between the solid Earth and the atmosphere. In urban environments with high ambient noise levels the effectiveness of event detection methods is unclear, particularly when deployment restrictions result in an irregular receiver array geometry. Here we deploy a dense nodal array for 1 month in the highly populated city state of Singapore. We develop a new detection method based on image processing that we call spectrogram stacking, which detects anomalous, coherent spectral energy across the array. It simultaneously detects multiple classes of signal with differing spectral content and aids event classification, so it is particularly useful for signal exploration when signal characteristics are unknown. Our approach detects more local events compared to the traditional STA/LTA and waveform similarity methods, while all methods detect similar numbers of teleseismic and regional earthquakes. Local events are principally man-made non-earthquake sources, with several events from the same location exhibiting repeating waveforms. The closest earthquake occurs in peninsular Malaysia, in an area where no earthquakes have previously been detected. We also detect ground motion over a wide frequency range from discrete thunder events which show complex coupling between acoustic and elastic wavefield propagation. We suggest that care should be taken deciphering local high-frequency tectonic events in areas prone to thunder storms.

Description: Summary In the context of seismic imaging, full waveform inversion (FWI) is increasingly popular. Because of its lower numerical cost, the acoustic approximation is often used, especially at the exploration geophysics scale, both for tests and for real data. Moreover, some research domains such as helioseismology face true acoustic media for which FWI can be useful. In this work, an argument that combines particle relabelling and homogenization is used to show that the general acoustic inverse problem based on band-limited data is intrinsically non-unique. It follows that the results of such inversions should be interpreted with caution. To illustrate these ideas, we consider 2-D numerical FWI examples based on a Gauss-Newton iterative inversion scheme, and demonstrate effects of this non-uniqueness in the local optimization context.

Description: Summary The GRACE and GRACE-FO missions have provided an unprecedented quantification of large-scale changes in the water cycle. However, it is still an open problem of how these missions’ data can be referenced to a ground truth. Meanwhile, stationary optical clocks show fractional instabilities below 10−18 when averaged over an hour, and continue to be improved in terms of stability and accuracy, uptime, and transportability. The frequency of a clock is affected by the gravitational redshift, and thus depends on the local geopotential; a relative frequency change of 10−18 corresponds to a geoid height change of about 1 cm. Here we suggest that this effect could be exploited for sensing large-scale temporal geopotential changes via a network of clocks distributed at the Earth’s surface. In fact, several projects have already proposed to create an ensemble of optical clocks connected across Europe via optical fibre links. Our hypothesis is that a clock network with collocated GNSS receivers spread over Europe – for which the physical infrastructure is already partly in place – would enable us to determine temporal variations of the Earth’s gravity field at time scales of days and beyond, and thus provide a new means for validating satellite missions such as GRACE-FO or a future gravity mission. Here, we show through simulations how glacial, hydrological and atmospheric variations over Europe could be observed with clock comparisons in a future network that follows current design concepts in the metrology community. We assume different scenarios for clock and GNSS uncertainties and find that even under conservative assumptions – a clock error of 10−18 and vertical height control error of 1.4 mm for daily measurements – hydrological signals at the annual time scale and atmospheric signals down to the weekly time scale could be observed.

Description: Summary Reflection seismic imaging usually suffers from a loss of resolution and contrast because of the fluctuations of the wave velocities in the Earth’s crust. In the literature, phase distortion issues are generally circumvented by means of a background wave velocity model. However, it requires a prior tomography of the wave velocity distribution in the medium, which is often not possible, especially in depth. In this paper, a matrix approach of seismic imaging is developed to retrieve a three-dimensional image of the subsoil, despite a rough knowledge of the background wave velocity. To do so, passive noise cross-correlations between geophones of a seismic array are investigated under a matrix formalism. They form a reflection matrix that contains all the information available on the medium. A set of matrix operations can then be applied in order to extract the relevant information as a function of the problem considered. On the one hand, the background seismic wave velocity can be estimated and its fluctuations quantified by projecting the reflection matrix in a focused basis. It consists in investigating the response between virtual sources and detectors synthesized at any point in the medium. The minimization of their cross-talk can then be used as a guide star for approaching the actual wave velocity distribution. On the other hand, the detrimental effect of wave velocity fluctuations on imaging is overcome by introducing a novel mathematical object: The distortion matrix. This operator essentially connects any virtual source inside the medium with the distortion that a wavefront, emitted from that point, experiences due to heterogeneities. A time reversal analysis of the distortion matrix enables the estimation of the transmission matrix that links each real geophone at the surface and each virtual geophone in depth. Phase distortions can then be compensated for any point of the underground. Applied to passive seismic data recorded along the Clark branch of the San Jacinto fault zone, the present method is shown to provide an image of the fault until a depth of 4 km over the frequency range 10-20 Hz with an horizontal resolution of 80 m. Strikingly, this resolution is almost one eighth below the diffraction limit imposed by the geophone array aperture. The heterogeneities of the subsoil play the role of a scattering lens and of a transverse wave guide which increase drastically the array aperture. The contrast is also optimized since most of the incoherent noise is eliminated by the iterative time reversal process. Beyond the specific case of the San Jacinto Fault Zone, the reported approach can be applied to any scales and areas for which a reflection matrix is available at a spatial sampling satisfying the Nyquist criterion.

Description: Summary Thermoremanent magnetization (TRM), the primary magnetic memory of igneous rocks, depends for its stability through geologic time on mineral carriers with high coercivities and high unblocking temperatures. The paleomagnetic record of past magnetic field directions and intensities is the key to unraveling Earth's tectonic history. Yet we still do not fully understand how the familiar mineral magnetite, in the micrometer grain size range typically responsible for stable TRM, acquires and holds its signal. Direct indicators of magnetite remanence-carrying capacity and coercivity at high temperature T are saturation remanence relative to saturation magnetization Mrs/Ms and coercive force Hc. This study is the first to measure the variation of these hysteresis properties for magnetite, from room temperature to the Curie point, across the entire size range from 25 nm to 135 µm, covering superparamagnetic, single-domain, vortex, pseudo-single-domain and multidomain magnetic behaviour. The paper focuses on: (1) Hc(T) and Mrs(T) observations and their reproducibility; (2) mathematical relationships of Hc(T) and Mrs(T) to Ms(T), used in modeling TRM and for unbiased comparisons of thermal variations; (3) the shapes of magnetite grains and the number of domains they contain, revealed by demagnetizing factors N = Hc/Mrs; and (4) the grain-size dependences of Hc and Mrs at ordinary and elevated T, delineating domain structure changes and mechanisms of coercivity.

Description: Summary Metropolitan France is a region of slow tectonic deformation with sparse seismicity. On 11 November 2019, the ML 5.4 Le Teil earthquake became the largest seismic event recorded in the last 16 years. This event was recorded by the national seismic networks and also by a wide variety of other geophysical techniques including infrasound and InSAR measurements. These complementary technologies offer the opportunity to investigate in detail the earthquake source characteristics and the associated ground motion attenuation. Both seismic waveform inversions and InSAR interferogram reveal a shallow rupture on a reverse fault with an associated moment magnitude of 4.8-4.9. Infrasound signals also provide fast evidences pointing towards the area of ground surface displacements, which coincides with La Rouvière fault, in the Cévennes fault system, known as a formerly active normal fault during the Oligocene. The very significant amount of seismic records also helps toward validating the GMPE laws available for the region. This multi-technology characterisation documents the kinematics of this rare example of shallow intraplate fault reactivation.

Description: Summary Precise real time estimates of earthquake magnitude and location are essential for early warning and rapid response. While recently multiple deep learning approaches for fast assessment of earthquakes have been proposed, they usually rely on either seismic records from a single station or from a fixed set of seismic stations. Here we introduce a new model for real-time magnitude and location estimation using the attention based transformer networks. Our approach incorporates waveforms from a dynamically varying set of stations and outperforms deep learning baselines in both magnitude and location estimation performance. Furthermore, it outperforms a classical magnitude estimation algorithm considerably and shows promising performance in comparison to a classical localization algorithm. Our model is applicable to real-time prediction and provides realistic uncertainty estimates based on probabilistic inference. In this work, we furthermore conduct a comprehensive study of the requirements on training data, the training procedures and the typical failure modes. Using three diverse and large scale data sets, we conduct targeted experiments and a qualitative error analysis. Our analysis gives several key insights. First, we can precisely pinpoint the effect of large training data; for example, a four times larger training set reduces average errors for both magnitude and location prediction by more than half, and reduces the required time for real time assessment by a factor of four. Second, the basic model systematically underestimates large magnitude events. This issue can be mitigated, and in some cases completely resolved, by incorporating events from other regions into the training through transfer learning. Third, location estimation is highly precise in areas with sufficient training data, but is strongly degraded for events outside the training distribution, sometimes producing massive outliers. Our analysis suggests that these characteristics are not only present for our model, but for most deep learning models for fast assessment published so far. They result from the black box modeling and their mitigation will likely require imposing physics derived constraints on the neural network. These characteristics need to be taken into consideration for practical applications.

Description: Summary We use data from the Cascadia Initiative (CI) amphibious array and the USArray Transportable Array to construct and compare Rayleigh wave isotropic and azimuthally anisotropic phase speed maps across the Juan de Fuca and Gorda Plates extending onto the continental northwestern U.S. Results from both earthquakes (28–80 s) as well as ambient noise two- and three-station interferometry (10–40 s) are produced. Compared with two-station interferometry, three-station direct wave interferometry provides $〉 50\%$ improvement in the signal-to-noise ratio (SNR) and the number of dispersion measurements obtained, particularly in the noisier oceanic environment. Earthquake and ambient noise results are complementary in bandwidth and azimuthal coverage, and agree within about twice the estimated uncertainties of each method. We, therefore, combine measurements from the different methods to produce composite results that provide an improved data set in accuracy, resolution, and spatial and azimuthal coverage over each individual method. A great variety of both isotropic and azimuthally anisotropic structures are resolved. Across the oceanic plate, fast directions of anisotropy with 180○ periodicity (2ψ) generally align with paleo-spreading directions while 2ψ amplitudes mostly increase with lithospheric age, both displaying substantial variations with depth and age. Strong ($〉 3\%$) apparent anisotropy with 360○ periodicity (1ψ) is observed at long periods (〉50 s) surrounding the Cascade Range, probably caused by backscattering from heterogeneous isotropic structures.

Description: Summary Temporal changes in subsurface properties, such as seismic wavespeeds, can be monitored by measuring phase shifts in the coda of two seismic waveforms that share a similar source-receiver path but that are recorded at different times. These nearly identical seismic waveforms are usually obtained either from repeated earthquake waveforms or from repeated ambient noise cross-correlations. The five algorithms that are the most popular to measure phase shifts in the coda waves are the Windowed Cross Correlation (WCC), Trace Stretching (TS), Dynamic Time Warping (DTW), Moving Window Cross Spectrum (MWCS), and Wavelet Cross Spectrum (WCS). The seismic wavespeed perturbation is then obtained from the linear regression of phase shifts with their respective lag times under the assumption that the velocity perturbation is homogeneous between (virtual or active) source and receiver. We categorize these methods into the time domain (WCC, TS, DTW), frequency domain (MWCS), and wavelet domain (WCS). This study complements this suite of algorithms with two additional wavelet-domain methods, which we call Wavelet Transform Stretching (WTS) and Wavelet Transform Dynamic Time Warping (WTDTW), wherein we apply traditional stretching and dynamic time warping techniques to the wavelet transform. This work aims to verify, validate, and test the accuracy and performance of all methods by performing numerical experiments, in which the elastic wavefields are solved for in various 2D heterogeneous halfspace geometries. Through this work, we validate the assumption of a linear increase in phase shifts with respect to phase lags as a valid argument for fully homogeneous and laterally homogeneous velocity changes. Additionally, we investigate the sensitivity of coda waves at various seismic frequencies to the depth of the velocity perturbation. Overall, we conclude that seismic wavefields generated and recorded at the surface lose sensitivity rapidly with increasing depth of the velocity change for all source-receiver offsets. However, measurements made over a spectrum of seismic frequencies exhibit a pattern such that wavelet methods, and especially WTS, provide useful information to infer the depth of the velocity changes.

Description: Summary As the largest and most active intracontinental orogenic belt on Earth, the Tien Shan (TS) is a natural laboratory for understanding the Cenozoic orogenic processes driven by the India-Asia collision. On 19 January 2020, a Mw 6.1 event stuck the Kalpin region, where the southern frontal TS interacts with the Tarim basin. To probe the local ongoing orogenic processes and potential seismic hazard in the Kalpin region, both interseismic and instantaneous deformation derived from geodetic observations are employed in this study. With the constraint of interseismic global navigation satellite system (GNSS) velocities, we estimate the décollement plane parameters of the western Kalpin nappe based on a two-dimensional dislocation model, and the results suggest that the décollement plane is nearly subhorizontal with a dip of ∼3° at a depth of 24 km. Then, we collect both Sentinel-1 and ALOS-2 satellite images to capture the coseismic displacements caused by the 2020 Kalpin event, and the interferometric synthetic aperture radar (InSAR) images show a maximum displacement of 7 cm in the line of sight near the epicentral region. With these coseismic displacement measurements, we invert the source parameters of this event using a finite-fault model. We determine the optimal source mechanism in which the fault geometry is dominated by thrust faulting with an E–W strike of 275° and a northward dip of 11.2°, and the main rupture slip is concentrated within an area 28.0 km in length and${ m{,,}}$10.3 km in width, with a maximum slip of 0.3 m at a depth of 6–8 km. The total released moment of our preferred distributed slip model yields a geodetic moment of 1.59 × 1018 N$cdot $m, equivalent to Mw 6.1. The contrast of the décollement plane depth from interseismic GNSS and the rupture depth from coseismic InSAR suggests that a compression still exists in the Kalpin nappe forefront, which is prone to frequent moderate events and may be at risk of a much more dangerous earthquake.

Description: Summary The Ice, Cloud, and land Elevation Satellite 2 (ICESat-2) laser altimetry mission, launched in September 2018, uses 6 parallel lidar tracks with very fine along-track resolution (15 m) to measure the topography of ice, land, and ocean surfaces. Here we assess the ability of ICESat-2 ocean data to recover oceanographic signals ranging from surface gravity waves to the marine geoid. We focus on a region in the tropical Pacific and study photon height data in both the wavenumber and space domain. Results show that an ICESat-2 single track can recover the marine geoid at wavelengths 〉 20 km which is similar to the best radar altimeter data. The wavelength and propagation direction of surface gravity waves are sometimes well resolved by using a combination of the strong and weak beams, which are separated by 90 m. We find higher than expected power in the 3 km to 20 km wavelength band where geoid and ocean signals should be small. This artificial power is caused by the projection of 2-D surface waves with ∼300 m wavelengths into longer wavelengths (5-10 km) because of the 1-D sampling along the narrow ICESat-2 profile. Thus ICESat-2 will not provide major improvements to the geoid recovery in most of the ocean.

Description: Summary We expand the application of spatial autocorrelation (SPAC) from typical 1D Vs profiles to quasi-3D imaging via Bayesian Monte-Carlo inversion (BMCI) using a dense nodal array (49 nodes) located at the Utah Frontier Observatory for Research in Geothermal Energy (FORGE) site. Combinations of 4 and 9 geophones in subarrays provide for 36 and 25 1D Vs profiles, respectively. Profiles with error bars are determined by calculating coherency functions that fit observations in a frequency range of 0.2–5 Hz. Thus, a high-resolution quasi-3D Vs model from the surface to 2.0 km depth is derived and shows that surface-parallel sedimentary strata deepen to the west, consistent with a 3D seismic reflection survey. Moreover, the resulting Vs profile is consistent with a Vs profile derived from distributed acoustic sensing (DAS) data located in a borehole at the FORGE site. The quasi-3D velocity model shows that the base of the basin dips ∼22o to the west and topography on the basement interface coincident with the Mag Lee Wash suggests that the bedrock interface is an unconformity.

Description: Summary Scattered seismic coda waves are frequently used to characterize small scale medium heterogeneities, intrinsic attenuation or temporal changes of wave velocity. Spatial variability of these properties raises questions about the spatial sensitivity of seismic coda waves. Especially the continuous monitoring of medium perturbations using ambient seismic noise led to a demand for approaches to image perturbations observed with coda waves. An efficient approach to localize spatial and temporal variations of medium properties is to invert the observations from different source-receiver combinations and different lapse times in the coda for the location of the perturbations. For such an inversion, it is key to calculate the coda-wave sensitivity kernels which describe the connection between observations and the perturbation. Most discussions of sensitivity kernels use the acoustic approximation in a spatially uniform medium and often assume wave propagation in the diffusion regime. We model 2-D multiple non-isotropic scattering in a random elastic medium with spatially variable heterogeneity and attenuation using the radiative transfer equations which we solve with the Monte-Carlo method. Recording of the specific energy density of the wavefield that contains the complete information about the energy density at a given position, time and propagation direction allows us to calculate sensitivity kernels according to rigorous theoretical derivations. The practical calculation of the kernels involves the solution of the adjoint radiative transport equations. We investigate sensitivity kernels that describe the relationships between changes of the model in P- and S-wave velocity, P- and S-wave attenuation, and the strength of fluctuation on the one hand and seismogram envelope, travel time changes and waveform decorrelation as observables on the other hand. These sensitivity kernels reflect the effect of the spatial variations of medium properties on the wavefield and constitute the first step in the development of a tomographic inversion approach for the distribution of small-scale heterogeneity based on scattered waves.

Description: Summary The Pacific Coast of Central North America is a geodynamically complex region which has been subject to various geophysical processes operating on different time scales. Glacial isostatic adjustment (GIA), the ongoing deformational response of the solid Earth to past deglaciation, is an important geodynamic process in this region. In this study we apply Earth models with 3D structure to determine if the inclusion of lateral structure can explain the poor performance of 1D models in this region. Three different approaches are used to construct 3D models of the Earth structure. For the first approach, we adopt an optimal 1D viscosity structure from previous work and add lateral variations based on four global seismic shear wave velocity anomalies and two global lithosphere thickness models. The results based on these models indicate that the addition of lateral structure significantly impacts modelled RSL changes, but the data-model fits are not improved. The global seismic models are limited in spatial resolution and so two other approaches were considered to produce higher resolution models of 3D structure: inserting a regional seismic model into two of the global seismic models and, explicitly incorporating regional structure of the Cascadia subduction zone and vicinity, i.e. the subducting slab, the overlying mantle wedge, and the plate boundary interface. The results associated with these higher resolution models do not reveal any clear improvement in satisfying the RSL observations, suggesting that our estimates of lateral structure are inaccurate and/or the data-model misfits are primarily due to limitations in the adopted ice-loading histories. The different realisations of 3D Earth structure gives useful insight to uncertainty associated with this aspect of the GIA model. Our results indicate that improving constraints on the deglacial history of the southwest sector of the Cordilleran ice sheet is an important step towards developing more accurate of GIA models for this region.

Description: Summary We develop and apply a method to constrain the space- and frequency-dependent location of ambient noise sources. This is based on ambient noise cross-correlation inversion using numerical wavefield simulations, which honour 3-D crustal and mantle structure, ocean loading, and finite-frequency effects. In the frequency range from 3 - 20 mHz, our results constrain the global source distribution of the Earth’s hum, averaged over the southern hemisphere winter season of 9 years. During southern-hemisphere winter, the dominant sources are largely confined to the southern hemisphere, the most prominent exception being the Izu-Bonin-Mariana arc, which is the most active source region between 12 - 20 mHz. Generally, strong hum sources seem to be associated with either coastlines or bathymetric highs. In contrast, deep ocean basins are devoid of hum sources. While being based on the relatively small number of STS-1 broadband stations that have been recording continuously from 2004 - 2013, our results demonstrate the practical feasibility of a frequency-dependent noise source inversion that accounts for the complexities of 3-D wave propagation. It may thereby improve full-waveform ambient noise inversions and our understanding of the physics of noise generation.

Description: SUMMARY In this study, we aim to solve the seismic inversion in the Bayesian framework by generating samples from the posterior distribution. This distribution incorporates the uncertainties in the seismic data, forward model, and prior information about the subsurface model parameters; thus, we obtain more information through sampling than through a point estimate (e.g. maximum a posteriori method). Based on the numerical cost of solving the forward problem and the dimensions of the subsurface model parameters and observed data, sampling with Markov chain Monte Carlo (MCMC) algorithms can be prohibitively expensive. Herein, we consider the promising Langevin dynamics MCMC algorithm. However, this algorithm has two central challenges: (1) the step size requires prior tuning to achieve optimal performance and (2) the Metropolis–Hastings acceptance step is computationally demanding. We approach these challenges by proposing an adaptive step-size rule and considering the suppression of the Metropolis–Hastings acceptance step. We highlight the proposed method’s potential through several numerical examples and rigorously validate it via qualitative and quantitative evaluation of the sample quality based on the kernelized Stein discrepancy (KSD) and other MCMC diagnostics such as trace and autocorrelation function plots. We conclude that, by suppressing the Metropolis–Hastings step, the proposed method provides fast sampling at efficient computational costs for large-scale seismic Bayesian inference; however, this inflates the second statistical moment (variance) due to asymptotic bias. Nevertheless, the proposed method reliably recovers important aspects of the posterior, including means, variances, skewness and 1-D and 2-D marginals. With larger computational budget, exact MCMC methods (i.e. with a Metropolis–Hastings step) should be favoured. The results thus obtained can be considered a feasibility study for promoting the approximate Langevin dynamics MCMC method for Bayesian seismic inversion on limited computational resources.