Geophysical Journal International

SummarySeismic tomography is a principal method for studying mantle structure, but imaging of Earth’s wavespeed anomalies is conditioned by seismic wave sampling. Global models use misfit criteria that may strive for balance between portions of the data set but can leave important regional domains underserved. We evaluate two full-waveform global tomography wavespeed models, GLAD-M25 and SEMUCB-WM1, in the mantle below the Pacific Ocean. The region of the South Pacific Superswell contains multiple hotspots which may be fed by plumes anchored in the Large Low Shear-Velocity Province at the base of the mantle. The uneven distribution of seismic receivers worldwide leaves several candidate plumes beneath various hotspots poorly resolved. We assess the regional quality of GLAD-M25 relative to its global performance using a partition of the seismic waveform data used in its construction. We evaluate synthetic waveforms computed using the spectral-element method to determine how well they fit the data according to a variety of criteria measured across multiple seismic phases and frequency bands. The distributions of travel-time anomalies that remain in GLAD-M25 are wider for trans-Pacific paths than globally, suggesting comparatively insufficiently resolved seismic velocity structure in the region of interest. Hence, Pacific-centered regional inversions, based on (augmented) subsets of the global data set have the potential to enhance the resolution of velocity structure. We compare GLAD-M25 and SEMUCB-WM1 by cross-validation with a new, independent, data set. Our results reveal that short- and long-wavelength structure is captured differently by the two models. Our findings lead us to recommend focusing future model iteration on and around the Pacific Superswell and adding data that sample new corridors, especially using ocean sensors, to better constrain seismic velocity structure in this area of significant geodynamic complexity.

SummaryWe present a new 3D crustal P-wave velocity (VP) model for the greater Alpine region (GAR). We use and merge three different high-quality datasets for local earthquake tomography covering 24 years, starting from January 1st, 1996, up to December 31st, 2019. We processed and repicked the waveforms from the events reported by the European-Mediterranean Seismological Centre with M > 3.0 inside the greater Alpine region for the period between May 2007 and December 2015 using a recently developed automated arrival time-picking procedure (ADAPT framework). This allows bridging the data gap between previously published (pre-2007) datasets and the recently published AlpArray research seismicity catalogue and thus provides a high-quality, highly consistent set of P-wave arrival times covering 24 years. With this data set we derived a new minimum 1D VP model and associated station delays covering the entire GAR. Subsequently, we performed a series of local-earthquake-tomography (LET) inversions obtaining a 3D VP model with a horizontal node spacing of 20×20 km and between 7 to 15 km variable vertical spacing in the well-resolved area of investigation, thus improving the spatial and uniformly high-resolution coverage compared to previous LET studies in the area. For well-known major crustal structures, such as, e.g. the geophysical Ivrea body, deep foreland basins and main orogenic crustal roots, our tomographic results correlate well with features documented by various previous seismic studies in the region. This correlation increases our confidence in the model's accuracy throughout the well-resolved area. Additionally, our model reveals previously poorly known, or unknown crustal features and it documents details in the Moho topography throughout the region. Eventually, we present a LET-Moho map (VP isoline of 7.25 km/s) for the GAR with spatially nearly uniform resolution and document its comparison with previously published Moho maps. The new regional 3D VP crustal model also correlates well with a previously published VS crustal model obtained by ambient noise tomography. These comparisons document the new LET results of combined 3D VP crustal velocities and Moho topography being intrinsically consistent and reliable within the region of high resolution. Hence, in addition to further improving our understanding of crustal structure geometries in the GAR, our results also provide pivotal information for a future reference seismic 3D crustal model of the region.

SummaryNumerical simulations of infinite Prandtl number convection in Cartesian domains have shown that a combination of internal and basal heating allows for behaviour not observed in either end-member cases of pure basal or pure internal heating. In particular, these mixed heating systems exhibit a decrease in the upper boundary layer velocity as internal heating increases. This leads to an inverse relationship between surface heat flow and boundary layer velocity. The inverse relationship has been attributed to boundary layer interactions, leading to deviations from classic boundary layer theory. Herein, we extend that work by presenting results from numerical experiments for mixed-heated convection in an isoviscous fluid in a fully 3D spherical domain. We show that an increase in internal heating causes a decrease in surface velocity, consistent with previous Cartesian results. We confirm that boundary layer interactions decrease with increased internal heating, which correlates with decreasing surface velocities. A scaling theory, previously applied to Cartesian geometry, is modified for spherical geometries and tested against the results of the numerical solutions. The modified scalings lead to good fits for temperature and heat flux trends. The scalings predict that velocities can decrease with increased internal heating from low to moderate internal heating rates and become constant at higher heating rates, consistent with numerical results. The quantitative match between velocity scalings and numerical results is not as good as observed for heat flow and temperature trends. We attribute this to surface velocities being more strongly affected by observed changes in convective wavelengths and planform transitions from sheet-like to plume-like downwellings as the rate of internal heating and/or basal heating increases.

SummarySeismic tomography is used to image subsurface structures at various scales, accomplished by solving a nonlinear and nonunique inverse problem. It is therefore important to quantify velocity model uncertainties for accurate earthquake locations and geological interpretations. Monte Carlo sampling techniques are usually used for this purpose, but those methods are computationally intensive, especially for large datasets or high-dimensional parameter spaces. In comparison, Bayesian variational inference provides a more efficient alternative by delivering probabilistic solutions through optimization. The method has been proven to be efficient in 2D tomographic problems. In this study, we apply variational inference to solve 3D double-difference (DD) seismic tomographic system using both absolute and differential travel time data. Synthetic tests demonstrate that the new method can produce more accurate velocity models than the original DD tomography method by avoiding regularization constraints, and at the same time provides more reliable uncertainty estimates. Compared to traditional checkerboard resolution tests, the resulting uncertainty estimates provide a better measure for the reliability of the solution. We further apply the new method to data recorded by a local dense seismic array around the San Andreas Fault Observatory at Depth (SAFOD) site along the San Andreas Fault (SAF) at Parkfield. Similar to previous studies, the obtained velocity models show significant velocity contrasts across the fault. More importantly, the new method produces velocity uncertainties of less than 0.34 km/s for ${{\rm{V}}}_p$ and 0.23 km/s for ${{\rm{V}}}_s$. We therefore conclude that variational inference provides an effective tool for solving 3D seismic tomographic problems and quantifying model uncertainties.

SummaryWe present a new seismic shear wave velocity model of the upper mantle of the Antarctic Plate region, AP2024. It includes the lithosphere and underlying mantle down to 660 km depth beneath both the continental and oceanic portions of the plate. To augment the limited seismic station coverage of Antarctica, we assemble very large regional and global data sets, comprising all publicly available broadband seismic data. The model is built using 785 thousand seismograms from over 27 thousand events and 8.7 thousand stations. It is constrained by both body and Rayleigh surface waves, ensuring the dense data sampling of the entire upper mantle depth range. The tomographic inversion is global but focused on the Antarctic Plate, with the data sampling maximised in the Southern Hemisphere, with elaborate automated and manual outlier analysis and removal performed on the regional data, and with the regularisation tuned for the region. The upper mantle of the Antarctic continent exhibits a bimodal nature. The sharp boundary along the Transantarctic Mountains separates the cratonic eastern from tectonic western Antarctica and shows a shear-velocity contrast of up to 17% at ∼100 km depth. The bimodal pattern is also seen in the oceanic part of the plate, with the older oceanic lithosphere beneath the Indian sector of the Southern Ocean showing higher shear velocities. The continental lithosphere in East Antarctica shows high velocity anomalies similar to those beneath stable cratons elsewhere around the world. It is laterally heterogeneous and exhibits significant thinning in the near-coastal parts of Dronning Maud Land and Wilkes Land. A low velocity channel is observed along the southern front of the West Antarctic Rift System and is probably related to Cenozoic rifting. High seismic velocity anomalies are detected beneath the Antarctic Peninsula and are likely to indicate fragments of the recently subducted Phoenix Plate Slab. Low velocity anomalies beneath Marie Byrd Land extend into the deep upper mantle and are consistent with a deep mantle upwelling feeding West Antarctica intraplate magmatism.

SummaryThis study investigates the effectiveness of inversion methods using tsunami waveforms to analyse volcanic tsunami sources, which are a type of non-seismic tsunami source. We focused on the 2018 Anak Krakatau tsunami triggered by a volcanic eruption. This study developed a static initial sea surface displacement model based on tsunami waveform inversion with data recorded at tide gauge stations using a Gaussian-shaped sea surface displacement for the unit source. A key characteristic of our model is that all initial velocity components of the tsunami were zero. We tested 12 scenarios for accuracy to determine the most plausible sea surface displacement. The optimal displacement model reasonably reproduced the observed tsunami waveforms. The calculated water volume at the initial sea surface displacement was reasonably consistent with the total collapse volume of the Anak Krakatau eruption by magnitude. These findings suggest that our approach to developing a static source model can effectively apply to non-seismic tsunami events. Although this approach offers simplified tsunami source modeling for tsunami estimation during volcanic eruptions with complex source dynamics, further validation is required for its application to other non-seismic tsunami events.

SummaryDespite significant progress in paleomagnetic research over the last century, the origin, evolution, and long-term behavior of the geomagnetic field remains poorly understood. One significant open question is when and how the inner core nucleated. Since geomagnetic field behavior is intrinsically linked to the thermal evolution of the core, scientists have turned to the global paleointensity record to search for proxies for inner core nucleation. From this record, two signals have been identified as possible indicators of inner core nucleation: (1) A spike in magnetic field strength between 1.5–1.0 Ga, and (2) an initially strong, but gradually decreasing field strength that resulted in a weak dynamo in the Ediacaran. Although both these hypotheses are vastly different, they do have one common challenge hindering rigorous testing: A paucity of paleointensity data. This is especially true for the Precambrian time period for which well-preserved outcrops are scarce and weathering/alteration is nearly inescapable. Despite making up almost 90% of Earth's history, data from this super eon comprises <10% of the global paleointensity database. This lack of data for most of Earth's history represents a considerable gap in our knowledge and greatly impedes our ability to understand the origin and evolution of our planet and its magnetic field. In an effort to fill in this gap, we performed paleointensity experiments on Precambrian-aged mafic dikes from India (Malani Igneous Suite and Bastar, Dharwar, and Bundelkhand cratons) with ages ranging from ∼740 Ma to ∼2.36 Ga. To monitor thermal alteration and minimize the effects of non-ideal grain sizes, the Thellier method following the IZZI protocol was used. Successful results were obtained for samples from the Bundelkhand (∼740 Ma) and Bastar (∼1.89 Ga) cratons. The Bastar results fall in a ∼40 Myr gap in the database and corroborate field trends predicted by the Monte Carlo axial dipole moment model (MCADAM), which suggests that intensity values were moderately low (2–4 × 1022 Am2) in the middle Paleoproterozoic. The Bundelkhand result suggests that the field may have been rapidly decaying in the late Tonian to early Cryogenian.

SummarySite-specific Probabilistic Tsunami Hazard Assessment (PTHA) is a powerful tool for coastal planning against tsunami risk. However, its typically high computational demands led to the introduction of a Monte Carlo Stratified Importance Sampling (SIS) approach, which selects a representative subset of scenarios for numerical inundation simulations. We here empirically validate this sampling approach, for the first time to our knowledge, using an existing extensive dataset of numerical inundation simulations for two coastal sites in the Mediterranean Sea (Catania and Siracusa, both located in Sicily, Italy). Moreover, we propose a modified importance sampling function to prioritise seismic tsunami scenarios based on their arrival time at an offshore point near the target site, in addition to their wave amplitude and occurrence rate as leveraged in the previous work. This sampling function is applied separately in each earthquake magnitude bin, and allows denser sampling of near-field earthquakes to whose variations tsunamis are very sensitive. We compare the confidence intervals of the offshore PTHA estimates obtained with the new and the original importance sampling functions. Then, we benchmark our onshore PTHA estimates obtained with both functions against the inundation PTHA calculated using the full set of scenarios. We also test the assumption that onshore random errors follow a normal distribution, as found previously for the offshore case. As a result of the benchmarks, we find that the SIS approach works satisfactorily. Introducing the arrival time as an additional sampling factor enhances the precision of the estimates of both the mean and the percentiles for the two coastal sites considered. With this modification it is possible to deal efficiently with heterogeneous near-field earthquake sources involving coastal deformation at Catania and Siracusa, in addition to regional crustal and subduction sources. By comparing the sampling errors with the model (epistemic) uncertainty, an optimal trade-off between the number of simulations employed and the uncertainty of the PTHA model can be found, even for such a complex situation. A relatively small number of scenarios, on the order of a few thousand, is sufficient to perform site-specific PTHA for practical applications. These numbers correspond to 4–8% of the already reduced ensembles used in previous assessments at the same sites.

SummaryThe Pamir tectonic zone originates from the intense collision of the Indo-Eurasian plate. Identifying the faults in the Pamir region region is essential for elucidating the collision mechanism and seismic characteristics. This paper compares the effect of the two-dimensional discrete wavelet transform (DWT2D) and the non-subsampled shearlet transform (NSST) on gravity field separation through synthetic model gravity field experiments. The results show that NSST can avoid the Gibbs phenomenon of DWT2D and better maintain the gravity field distribution. The surface gravity disturbances data of the Global Gravity Model Plus (GGMplus) with a high-spatial resolution (7.2 arcsec or approximately 200 m) is employed to separate the region-residual gravity fields in the neighbouring domain of the Pamir region based on the NSST. Furthermore, the gravity gradient tensor (GGT) is computed, and the correspondence between the GGT and the location and strike of the surrounding faults is analyzed. The results show that the GGT component and its various combinations can effectively identify shallow and deep faults, the residual field GGT and its combinations can effectively identify the distribution and direction of shallow faults, and the regional field GGT and its combinations can effectively identify the distribution and direction of deep faults. The existence of north-south trending faults in the Pamir-Hindu Kush region is widely accepted. However, our study has revealed an east-west trending concealed fault in the deep areas of the Hindu Kush (Depth > 200 km). This finding provides significant insights for studying the bidirectional subduction of the Indian and Eurasian plates. This research not only helps us to analyze the tectonic characteristics of the shallow and deep parts of the region separately but also provides complementary information for investigating the distribution of deep underground faults, especially when fault inversion of intermediate to deep source earthquakes is limited by factors such as uncertainty in source depth and complexity of seismic wave velocities.

SummaryIn Pakistan, the relative displacement between the Indian and Eurasian plates is accommodated by a left lateral transpression zone comprising the Chaman and Ghazaband faults and the Sulaiman Range. The current tectonic deformation of the Sulaiman Range is known only from some focal mechanisms and a few neotectonic studies. In this study, we propose an InSAR quantification of current tectonic deformation using the Sentinel 1 satellite. Velocity maps for the ascending and descending tracks enabled us to locate active faults affected by creep: the Harnaï and Kingri strike-slip faults, and the Gwal-Bagh thrust. We propose a numerical simulation that considers these faults as well as the level of detachment fold-and-thrust belt. Our results suggest the existence of out-of-sequence deformation along the Gwal-Bagh thrust, creep along the Harnaï and Kingri strike-slip faults, and slip along the décollement of the Sulaiman Range. The eastern part of the Sulaiman Range is characterized by a partitioning of the deformation with a left lateral strike-slip along the N170° Kingri fault and an eastward thrust. In contrast, the western part is characterized by north-south compressive deformation associated with right lateral strike-slip on the Harnaï N120° fault. Modelling of the co-seismic deformation of the 21 October 2021 earthquake shows that this earthquake occurred on a fault with a ramp geometry but affected by a strike-slip motion.

SummaryThe central Pamir plateau moves northward and collides into Eurasia at a rate that varies significantly over its 600 km-wide extension. However, the active structures accounting for such internal shear strain remain enigmatic. In this study, we use Interferometric Synthetic Aperture Radar (InSAR) data to investigate the coseismic and postseismic deformation of the ${M}_w6.9$ Sarez earthquake on 23 February 2023. Using a Bayesian framework, we find the most likely seismogenic fault geometry and explore the full solution space of slip distributions. Our results highlight the mainshock ruptures a nearly NNE fault dipping to the southeast. The finite-fault model exhibits a purely left-lateral strike-slip mechanism with little to no slip reaching the surface. Most of the coseismic slip remains confined to a depth of ∼5 to 20 km, consistent with a large shallow slip deficit. Postseismic afterslip, which decays rapidly within the month following the mainshock, cannot compensate for such coseismic shallow slip deficit. Integrating the analysis of coseismic slip, postseismic deformation, and regional seismic activity, we argue that in the central Pamir, significant north-south shear strain is accommodated along multiple parallel faults, often unmapped, hence posing a significant seismic hazard.

SummaryDetailed investigations of aseismic slow slip events (SSEs) are crucial for estimating the strain budget and SSE mechanisms within subduction zones. The Suruga Trough, which includes the Tokai seismic gap, is an important area in Japan from a hazardous perspective. However, the aseismic slip history of this trough following the 2011 Tohoku earthquake is difficult to determine as a result of post-seismic deformation caused by the earthquake. In this study, we provided detailed imaging of the interplate aseismic slip in the Suruga Trough after the 2011 Tohoku Earthquake by applying a network inversion filter to global navigation satellite system data and considering viscoelastic deformation and afterslip caused by the earthquake. The analysis revealed the 2012 Shima long-term SSE (l-SSE), 2013-2016 Tokai l-SSE, 2017-2020 Shima l-SSE, and 2023-2024 Atsumi+Tokai l-SSE, with the slip area expanding to the area adjacent to the Tokai seismic gap from July 2023, consequently changing the stress state to promote the anticipated Tokai earthquake. The findings of this study suggest that the recurrence interval of the Tokai slow slip ranges from 10 to 13 years, with a duration of approximately 4-5 years and a total magnitude ranging from 6.5 to 7.1. The l-SSE zone shows that the upper-limit temperature threshold, which is the temperature at the upper bound of the l-SSE zone aligning the 350°C isothermal line in the Tokai segment, does not hold in the Suruga Trough. The change in strike direction of the l-SSE zone suggests that a discontinuous factor controls the l-SSE occurrence, such as high pore pressure caused by fluid infiltration to the plate interface. Furthermore, we explored a gap between the short-term SSE (s-SSE) and l-SSE zones, and the findings indicated a non-continuous transition from l-SSE to s-SSE, thus providing insights into the discontinuous factors that regulate l-SSE and s-SSE generation. The recurrence interval (10–13 years) and duration (4-5 years) of the Tokai SSEs are long, and their moment rates (1015.8Nm/day) are low compared to those of the l-SSEs in other regions. The SSE parameters suggest that the scaling law may not apply to SSEs in the Suruga-Nankai Trough with the prolonged duration.

SummaryOn 9th January 2020, a Mw 6.4 earthquake struck the central Koryak Highlands of eastern Siberia, northeast of the diffuse triple junction between the North American, Pacific and Eurasian plates. The largest earthquake recorded in the central Koryak Highlands to date, it provides an excellent opportunity to study the little-known active tectonics of this remote, sparsely instrumented region. We mapped coherent, coseismic surface deformation with Sentinel 1 Interferometric Synthetic Aperture Radar (InSAR), making this one of the highest-latitude earthquakes to be captured successfully with satellite radar, in spite of the rugged, snow-covered terrain. Elastic dislocation modelling, teleseismic back-projections, calibrated hypocentral relocations, and teleseismic moment tensor solutions are used to resolve a left-lateral fault trending northwestwards, proximal but perpendicular to a regional geological suture zone, the Khatyrka-Vyvenka Thrust. The earthquake probably ruptured unilaterally northwestwards along a 20 km long segment that appears indistinct in the local topography, and likely generated no surface rupture. We interpret that these observations are indicative of a structurally immature fault zone and estimate a seismogenic zone thickness of 10–15 km. The Koryak Highlands earthquake illustrates how terrane boundaries within cordilleran belts may continue to accommodate tectonic strain long after accretion, resulting in significant earthquakes even along hidden faults.

SummaryMagneto-Coriolis (MC) modes in Earth’s fluid core involve oscillations sustained by the combined effect of the Lorentz and Coriolis forces. Here, we investigate the properties of MC modes that involve purely axisymmetric flow, which we term axiMC modes. We provide a basic description of the wave dynamics of these modes, and simple predictions for the expected scalings of their frequency ω, decay rate λ, and quality factor Q based on a uniform ambient magnetic field. In particular, Q scales with the Elsasser number Λ, which depends on the square of the r.m.s. strength of the azimuthally averaged meridional field. When Λ > 1, Q > 1 and axiMC modes may be excited; when Λ ≪ 1, Q ≪ 1 and axiMC modes revert to quasi-free magnetic decay modes. We present computations of axiMC modes in an inviscid, electrically conducting sphere for two idealized ambient magnetic field configurations, a uniform axial field and an axial poloidal field. We show that a flow gradient in the axial direction is a key property of axiMC modes. For the uniform axial field, ω, λ and Q follow the scalings expected for a uniform field. For the axial poloidal field, the structure of the modes changes substantially when Λ ≳ 1, becoming more concentrated in regions of lower field strength. The combination of this structural change and advection of field lines by flow significantly increases λ, resulting in a Q that remains close to 1 even at high Λ. For a magnetic field strength inside the Earth’s core of a few mT, the gravest axiMC modes are expected to have periods in the range of one thousand to a few thousand years and a Q not substantially above 1. AxiMC modes may be connected to a part of the observed millennial changes in Earth’s magnetic field, may exchange axial angular momentum with the mantle, and hence may also explain a part of the observed millennial changes in length of day.

SummaryMulti-channel analysis of surface wave (MASW) is a non-destructive technique to characterize the sub-surface using the dispersive nature of Rayleigh waves. Field dispersion curves are inverted to predict the shear wave velocity structure of the ground and pavement profile. Adjusting the dynamic properties of the initially assumed soil profile necessitates information regarding the dominant sensitive layers. Therefore, a swift and accurate computation of the Jacobian of phase velocity is essential to generate an appropriate shear wave velocity profile and accelerate the inversion process. This is especially crucial for the 2D MASW survey, which requires hundreds of 1D inversions to create a high resolution 2D profile. Available numerical methods are computationally expensive and often suffer from instabilities for highly sensitive layers. The existing analytical methods involve mathematical complexities and require rigorous treatment. Furthermore, they are time-consuming and often found to be marginally faster than the numerical methods. Based on the fast delta matrix algorithm, the paper presents a new efficient analytical formulation of the Jacobian matrix of modal phase velocities concerning the layer parameters. The proposed algorithm leverages the simpler and fewer matrix elements of the fast delta matrix, thus significantly reducing the number of mathematical operations required. Additionally, it reduces the algorithm's cost by factorizing non-zero elements, thereby markedly reducing the computational time. Five different types of synthetic earth models are adopted from the published literature to validate the accuracy and efficacy of the newly developed algorithm. The presented work will significantly benefit the practicing engineers and geophysicists in processing field MASW test data.

SummaryMonitoring volcanic deformation is crucial for understanding volcanic behavior, but challenges like limited GNSS coverage, infrequent SAR data acquisitions, and coherence loss during eruptions complicate this task. Our study on the 2018 eruption of Sierra Negra utilizes Sentinel-1A/B images to track surface deformation patterns that revert to their initial state across three phases: before, during, and after the eruption. We implemented an adaptive workflow using the shortest temporal baseline of consecutive SAR image pairs, including InSAR, optical flow, and pixel offset tracking methods, to accurately capture surface displacement linked to the dynamics of the magma reservoir in the caldera and a nearby (sub-) horizontal dike. Results show that while the caldera subsided gradually over two months during lava flow (initially at a rate of several meters) until it began to uplift again, the northern region alternated between uplift and subsidence twice in the line-of-sight (LOS) direction. This pattern suggests repeated magma injections into the sub-horizontal dike sustained the lava flows from the northeastern fissure. The one-day difference between SAR images from ascending and descending tracks enabled us to estimate the underground magma transfer rate at approximately 60 m/h, which aligns with the magma migration trajectory indicated by seismic data. By integrating InSAR and offset tracking methods, we provide a comprehensive view of surface displacement throughout the volcanic eruption cycle.

AbstractCrystallographic preferred orientation (CPO) of peridotite minerals is frequently invoked to explain the widespread dependence of seismic wave speed on propagation direction in Earth’s mantle — a property known as seismic anisotropy. As established by rock mechanics experiments, CPO constitutes a direct signature of past and ongoing strain regimes experienced by rocks during mantle flow. Therefore, an improved understanding of CPO generation promises to yield valuable information on the rheology and corresponding deformation mechanisms activated through mantle dynamics. Simulating CPO in geodynamical models is computationally challenging and has often been restricted to steady-state mantle flows. However, within Earth’s vigorously convecting mantle the steady-state assumption is questionable, thus motivating the need to couple CPO simulations with time-evolving mantle flow models. Here, we present a new Python implementation of the D-Rex CPO model, called PyDRex, which predicts salient features of mineral grain size and orientation evolution whilst providing a well-documented, user-friendly interface that supports flexible coupling to geodynamical modelling frameworks. PyDRex also packages numerous post-processing routines for strain analysis and visualisation of grain orientation distributions. We provide a set of benchmark simulations based on previous D-Rex implementations that validate PyDRex and demonstrate sensitivities to model parameters for both steady-state and time-dependent flows. Analysis of benchmark results highlights the role of dynamic recrystallisation in controlling competing grain growth in both the softest and hardest crystallographic orientations. When employing a commonly used value for the grain boundary mobility parameter (M* = 125), we also find that transient CPO textures are generally not well resolved if crystals are represented by fewer than 5000 ‘grains’ (weighted orientation samples) — a configuration rarely employed in most previously published studies. Furthermore, kinematic corner-flow models suggest that CPO produced at mid-ocean ridges has a non-linear dependence on depth, which implies that even ostensibly simple mantle flows can result in complex distributions of seismic anisotropy. Our analyses motivate further experimental calibration of parameters controlling dynamic recrystallisation and potential improvements to the numerical treatment of subgrain nucleation.

SummaryForward modeling is crucial for seismic data processing, which is the core of reverse time migration and full-waveform inversion. Numerical simulation based on conventional elastic wave equations in stationary solids neglects the fluidity of fluids (e.g., seawater), making it difficult to simulate the propagation of seismic waves in moving fluids accurately. To solve the problem, we start with classical equations of fluid mechanics and derive a new set of elastic wave equations that can be used to simultaneously model wave propagation both in moving fluids and stationary solids. For high-precision numerical simulations, a staggered-grid finite-difference scheme is used to solve the proposed equations. Numerical tests on a homogeneous uniformly moving model demonstrate that the dynamic and kinematic characteristics (e.g., wavelength, amplitude) of elastic waves in moving fluids are quite different from those in stationary medium. Forward modeling for a two-layer model that has a flowing water layer and a stationary rock layer is used to study the reflection and transmission patterns of elastic waves in the solid-fluid interface. With the help of the superposition principle of vectors and Snell's law, the transmission angles can be easily calculated. A further test for a more complex stratified model indicates that the energy and travel time differences of reflected waves are expected to be evidence for the identification of moving fluids. Numerical experiments on the Marmousi II model demonstrate that the relative wavefield error is positively correlated with the maximum moving velocity and the wavelet dominant frequency.

SummaryDispersion curves of surface waves are widely used for the inversion of subsurface structures. To extract dispersion curves, many methods have been developed. Among them, multichannel analysis of surface waves such as slant stack and frequency-Bessel transform can extract not only the fundamental mode but also overtones. Inversion with overtones is proven to be more stable and has better resolution at greater depths. However, with a limited number of array receivers, artifacts and misfits due to array-layout effects arise in the dispersion spectra and impede the identification of dispersion curves. We evaluate the array-layout effects in the frequency-Bessel transform and calculate the array response functions which can help to mitigate artifacts and calibrate dispersion curves. We apply this technique to synthetic simulated, active source and ambient noise data. The artifacts caused by array-layout effects can be mitigated, which helps the identification of dispersion curves. We further calculate the Pearson correlation coefficient between the array response function and the dispersion spectrum section. It is used to calibrate the biases produced by the array-layout effects if we select dispersion curves by maximum values. The confidence intervals of the dispersion curves are then determined based on the correlation coefficients. It is helpful for the design of array layouts according to the investigation depths of interest.

SummaryPowerful infrasound (acoustic waves <20 Hz) can be produced by explosive volcanic eruptions. The long-range propagation capability, over hundreds to thousands of kilometers, of atmospheric infrasound motivates the development of regional or even global scale volcano-infrasound monitoring systems. Infrasound propagation paths are subject to spatiotemporal atmospheric dynamics, which lead to deviations in the direction-of-arrival (back-azimuth) observed at sensor arrays and contribute to source location uncertainty. Here we further investigate the utility of empirical climatologies combined with 3-dimensional ray-tracing for providing first-order estimates of infrasound propagation paths and back-azimuth deviation corrections. The intended application is in scenarios requiring rapid or precomputed infrasound propagation calculations, such as for a volcano-infrasound monitoring system. Empirical climatologies are global observationally based function fitting models of the atmosphere, representing robust predictors of the bulk diurnal to seasonal atmospheric variability. Infrasound propagation characteristics have previously been shown to have strong seasonal and diurnal components. At the International Monitoring System (IMS) infrasound station IS22, New Caledonia, quasi-continuous multi-year infrasound array detections show oscillating azimuthal variations for arrivals from volcanoes in Vanuatu, including Yasur (∼400 km range), Ambrym (∼670 km range), and Lopevi (∼650 km range). We perform 3-dimensional ray-tracing to model infrasound propagation from the Ambrym and Yasur volcano locations to IS22 every six hours (00:00, 06:00, 12:00, and 18:00 UTC) for every day of 2004 and 2019 for Ambrym and Yasur, respectively and evaluate the results as compared to the multi-year observations. We assess a variety of models and parameterizations, including both empirical climatologies and hybrid descriptions; range-independent and range dependent atmospheric discretizations; and unperturbed and perturbed range-independent empirical climatologies. The hybrid atmospheric descriptions are composed of ERA 5 reanalysis descriptions from the European Centre for Medium-Range Weather Forecasts (ECMWF) below ∼80 km altitude combined with empirical climatologies above. We propose and employ simple parametric perturbations to the empirical climatologies, which are designed to enhance the stratospheric duct and compensate for missing gravity wave perturbations not included in the climatologies, and thereby better match observations. We build year-long back-azimuth deviation interpolations from the simulations and compare them with three different multi-year array detection datasets from IS22 covering from 2003 up to 2022. Through a systematic comparison, we find that the range-independent empirical climatologies can capture bulk azimuth deviation variability and could thus be useful for rapid infrasound propagation calculation scenarios, particularly during favorable sustained propagation ducting conditions. We show that the hybrid models better describe infrasound propagation during periods of weak stratospheric ducting and during transient atmospheric changes such as stratospheric wind reversals. Overall, our results support the notion that climatologies, if perturbed to compensate for missing gravity wave structure, can improve rapid low-latency and precomputed infrasound source discrimination and location procedures.