Geophysical Journal International

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

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

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

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

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

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

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

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

SummaryOn March 18, 2021 an earthquake of magnitude Mw=6.0 occurred offshore the Algerian coasts and generated a tsunami with offshore amplitudes smaller than a few millimetres crossing the western Mediterranean Sea. The objective of this study is threefold: first, to determine whether seismic sources calculated in the context of tsunami early warning are relevant; second, to determine whether tsunami simulations are able to reproduce tide-gauge observations; and third, to define the sensitivity of simulations to the grid resolutions and tsunami parameters.In the Mediterranean Sea, a very small number of available coastal tide gauges recorded the tsunami. Among them, a few French tide gauge stations recorded water waves with amplitudes smaller than a few centimetres and with periods ranging from five to twenty minutes associated to harbour or bay resonances.Numerical simulations of the tsunami are performed by the operational code Taitoko for seven different source fault models. Three of them allow for a rapid source detection and characterization in the framework of tsunami warning at CENALT (Centre National d'Alerte aux Tsunamis, France). The integrated code Taitoko uses a system of multiple nested grids. Standard Boussinesq equations are solved in the Mediterranean grid, whereas nonlinear shallow water equations are solved in coastal and harbour grids with 25 and 5 m resolutions, respectively. Whatever the fault model, the observed time series of water heights are reproduced satisfactorily both in phase and amplitude by the model at Nice and Monaco but poorly at Port Mahon (Minorca) and Toulon.

SummaryDistributed acoustic sensing (DAS) is a relatively new technology for recording the propagation of seismic waves, with promising applications in both engineering and geophysics. DAS's ability to simultaneously collect high spatial resolution waveforms over long arrays suggests that it is well-suited for near-surface imaging applications such as 2D MASW (multichannel analysis of surface waves), which require, at a minimum, long, linear arrays of single-component receivers. The 2D MASW method uses a large number of sensor sub-arrays deployed along a linear alignment to produce 1D shear-wave velocity (VS) profiles beneath each sub-array. The 1D VS profiles are then combined to form a pseudo-2D VS image beneath the entire linear alignment that can be used for the purpose of identifying and characterizing lateral variations in subsurface layering. Traditionally, 2D MASW is conducted using arrays consisting of either 24 or 48 geophones. While additional receivers could easily be incorporated into the testing configuration, it is rare for researchers and practitioners to have access to greater numbers of seismographs and geophones. When a limited number of geophones are available for deployment, there is a need to pre-determine the geophone spacing and sub-array length prior to field data acquisition. Studies examining how the choice of sub-array geometry impacts the resulting pseudo-2D VS cross-sections have been largely limited to synthetic data. In response, this study utilizes DAS data to examine the effects of using various sub-array lengths by comparing pseudo-2D VS cross-sections derived from active-source waveforms collected at a well-characterized field site. DAS is particularly useful for 2D MASW applications because the sub-array geometry does not need to be determined prior to field data acquisition. We organize the DAS waveforms into multiple sets of overlapping MASW sub-arrays of differing lengths, ranging from 11 m to 47 m, along the same alignment, allowing for direct comparison of the derived pseudo-2D VS results at the site. We show that the length of the individual MASW sub-arrays has a significant effect on the resulting VS cross-sections, including the resolved location of a strong impedance contrasts at our study site, and evaluate the results relative to ground truth from invasive testing. Our results suggest that the choice of sub-array length is important and should be carefully chosen to meet project-specific goals. Furthermore, analysts may consider using multiple sub-array geometries during the data processing stage, as is made possible by DAS, to properly evaluate the uncertainty of 2D MASW results. This study demonstrates the potential of using DAS to collect data for 2D MASW in a manner that is efficient and flexible, and can be easily scaled up for use with very long arrays.

SummaryNatural fractures play a significant role in oil and gas reservoirs. Accurate predictions of fracture parameters are vital in reservoir prediction and oil and gas development. The birefringent phenomenon of shear waves in fractured media makes shear wave splitting (SWS) analysis an important tool in formulating fracture predictions. The traditional SWS analysis method is based on an orthogonal assumption of fast and slow shear waves. However, in an orthotropic medium composed of a background vertical transversely isotropic medium and a set of vertical fractures, fast and slow shear waves are not necessarily orthogonal. This causes the traditional SWS analysis method to fail. To solve this problem, we proposed an SWS analysis algorithm with a non-orthogonal assumption of fast and slow shear waves in this study. First, we introduced a parameter (difference angle) to characterise the angle between slow shear waves and the normal polarisation directions of the fast shear waves. Subsequently, based on the traditional two-parameter scanning algorithm, a parameter was added to facilitate three-parameter scanning. In addition, we derived an expression for the two-parameter scanning objective function using the non-orthogonal assumption. Two-parameter scanning can accurately extract fast and slow wave time delay data, but it cannot determine an accurate fast shear wave polarisation direction. Therefore, we optimised the three-parameter scanning algorithm as follows: first, we used two-parameter scanning to obtain the fast and slow wave time delays and then performed further scanning to determine the polarisation direction of the fast shear wave and difference angle. The optimisation algorithm significantly improved the computational efficiency. Subsequently, we tested the accuracy of this method using synthetic single-trace and three-component vertical seismic profile data. We demonstrated the implementation process of the three-parameter scanning method using actual data, separated fast and slow shear waves, and predicted fracture parameters. The final fracture parameters were verified.

SummaryAlthough ambient noise interferometry has been extensively utilized for seismic velocity tomography, its application in retrieving attenuation remains limited. This study presents a comprehensive workflow for extracting Rayleigh wave amplitude and attenuation from ambient noise, which consists of three phases: (1) retrieval of empirical Green's Functions (EGFs), (2) selection and correction of amplitude measurements, and (3) inversion of attenuation, site amplification, and noise intensity terms. Throughout these processes, an ‘asynchronous’ temporal flattening method is employed to generate high-quality EGFs while preserving relative amplitudes between stations. Additionally, a novel “t-symmetry” criterion is proposed for data selection along with the signal-to-noise ratio. Furthermore, 2-D sensitivity kernels are utilized to estimate the focusing/defocusing effect, which is then corrected in amplitude measurements. These procedures are designed to deliver reliable attenuation measurements while maintaining flexibility and automation. To validate the effectiveness of the proposed noise-based attenuation tomography approach, we apply it to a linear array, NCISP-6, located in NE China. The obtained results correlate reasonably well with known geological structures. Specifically, at short periods, high attenuation anomalies delineate the location of major sedimentary basins and faults; while at longer periods, a notable rapid increase of attenuation is observed beneath the Moho discontinuity. Given that attenuation measurements are more sensitive to porosity, defect concentration, temperature, melt, and volatile ratio than seismic velocities, noise-based attenuation tomography provides important additional constraints for exploring the crustal and upper mantle structures.

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

SummaryThe northern Chile Triple Junction (CTJ) is characterized by the ongoing subduction of the Nazca plate beneath the South American plate. The geological structures within the subduction zone undergo complex changes, resulting in significant tectonic activities and intense seismicity along the western margin of South America. Based on the Gravity Recovery and Climate Experiment (GRACE) data and earthquake catalogues, this study selects the northern CTJ area (25°S-40°S, 75°W-65°W) as the research object, adopts the mathematical methods of Independent Component Analysis (ICA) and Principal Component Analysis (PCA) to separate the earthquake-related signals within the GRACE data, and fits the changes of seismic b-values through the frequency-magnitude relationship. The characteristics of gravity changes before and after seismic events, the seismic activity parameter b-values, and the relationship between the gravity signals and b-values are discussed. The results show that mathematical methods can effectively extract seismic-related gravity components from the GRACE data. ICA, compared to PCA, provides better results in capturing the temporal variations associated with b-value time series, which exhibit good consistency in long-term trend changes. The average change of b-values in the study area is 0.66 ± 0.003, fluctuating over time. Generally, prior to larger seismic events, b-values tend to decrease. Along the western margin of South America, b-values are low; this aligns with the active tectonic activities between subducting plates. Additionally, a certain correlation between b-values and gravity changes is observed, but due to the influence of tectonic activities, the correspondence between b-values and gravity anomalies may not be consistent across different areas. The b-value is highly consistent with the strain rate model. Low b-values correspond to high strain rates along the western edge of South America, which is in line with the tectonic characteristics of frequent seismic activity in this area. A gradual concentration of gravity anomalies before major earthquakes is observed, accompanied by the gradual accumulation of smaller seismic events. Meanwhile, several months before the two major earthquakes, the spatial distribution of gravity appears to be similar to the co-seismic signals, but the nature of its generation remains to be explored. These methods and results not only add to the applications of GRACE in seismic studies but also raise questions for further exploration.

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

SummaryA Bayesian linear regression to determine the bias in the Nafe-Drake relationship between compressional velocity and density provides an improved model for the density structure of Kīlauea volcano, Hawaiʻi. In previous work, we combined the results of seismic tomography with the Nafe-Drake relationship between compressional velocity and density to explain the large values of gravity disturbances overlying the summits and rift zones of the island's volcanoes. These results were used to determine mechanisms for gravitational instability of the island flanks. Here we use laboratory measurements of the relationship of velocity and density for a wide range of Hawaiʻi island rocks as a prior in a Bayesian regression, with seismic tomography, to refine the 3D density structure for Kīlauea volcano. This refined structure shows dense bodies (3220 kg/m3) between 5 and 8 km below sea level that underly regions of magma storage, found from geodetic and geophysical studies, beneath the summit and East Rift Zone of Kīlauea volcano. Above these bodies, density iso-surfaces surround and cradle sources of pressure change determined from geodetic models, both at the summit and along the East Rift Zone. Continued subsidence of the summit following the 2018 eruption is aligned with a bowl-shaped density structure, formed primarily by density isosurfaces between 2800 and 2900 kg/m3 at 4 to 6 km depth. These surfaces underly the ∼3 km depth at which dike injection initiates, are largely aseismic, and from their density values are inferred to contain high concentrations of olivine. Taken together, these density structures are consistent with an olivine-rich mush with variable porosity that increases in density with depth and provides a mechanism to form olivine cumulates both at the summit and along the rift zones. This structural framework for Kīlauea volcano is consistent with melt and mush transport occurring over a large range of depths to accommodate the growth and spreading of the volcano.

SummaryFull waveform inversion has emerged as the state-of-the art high resolution seismic imaging technique, both in seismology for global and regional scale imaging and in the industry for exploration purposes. While gaining in popularity, full waveform inversion, at an operational level, remains a heavy computational process involving the repeated solution of large-scale 3D wave propagation problems. For this reason it is a common practice to focus the interpretation of the results on the final estimated model. This is forgetting full waveform inversion is an ill-posed inverse problem in a high dimensional space for which the solution is intrinsically non-unique. This is the reason why being able to qualify and quantify the uncertainty attached to a model estimated by full waveform inversion is key. To this end, we propose to extend at an operational level the concepts introduced in a previous study related to the coupling between ensemble Kalman filters and full waveform inversion. These concepts had been developed for 2D frequency-domain full waveform inversion. We extend it here to the case of 3D time-domain full waveform inversion, relying on a source subsampling strategy to assimilate progressively the data within the Kalman filter. We apply our strategy to an ocean bottom cable field dataset from the North Sea to illustrate its feasibility. We explore the convergence of the filter in terms of number of elements, and extract variance and covariance information showing which part of the model are well constrained and which are not. Analyzing the variance helps to gain insight on how well the final estimated model is constrained by the whole full waveform inversion workflow. The variance maps appears as the superposition of a smooth trend related to the geometrical spreading and a high resolution trend related to reflectors. Mapping lines of the covariance (or correlation matrix) to the model space helps to gain insight on the local resolution. Through a wave propagation analysis, we are also able to relate variance peaks in the model space to variance peaks in the data space. Compared to other posterior-covariance approximation scheme, our combination between Ensemble Kalman filter and full waveform inversion is intrinsically scalable, making it a good candidate for exploiting the recent exascale high performance computing machines.

SummaryIn petrophysics, physical rock properties are typically established through laboratory measurements of individual samples. These measurements predominantly relate to the specific sample and can be challenging to associate with the rock as a whole since the physical attributes are heavily reliant on the microstructure, which can vary significantly in different areas. Thus, the obtained values have limited applicability to the entirety of the original rock mass. To examine the dependence of petrophysical measurements based on the variable microstructure, we generate sets of random 2D microstructure representations for a sample, taking into account macroscopic parameters such as porosity and mean grain size. For each microstructure produced, we assess the electrical conductivity and evaluate how it is dependent on the microstructure’s variability. The developed workflow including microstructure modelling, finite element simulation of electrical conductivity as well as statistical and petrophysical evaluation of the results is presented. We show that the methodology can adequately mimic the physical behaviour of real rocks, showing consistent emulation of the dependence of electrical conductivity on connected porosity according to Archie’s law across different types of pore space (micro-fracture, inter-granular, and vuggy, oomoldic pore space). Furthermore, properties such as the internal surface area and its fractal dimension as well as the electrical tortuosity are accessible for the random microstructures and show reasonable behaviour. Finally, the possibilities, challenges and meshing strategies for extending the methodology to 3D microstructures are discussed.

SummaryThe reduction of effective normal stress during earthquake slip due to thermal pressurization of fault zone pore fluids is a significant fault weakening mechanism. Explicit incorporation of this process into frictional fault models involves solving the diffusion equations for fluid pressure and temperature outside the fault at each time step, which significantly increases the computational complexity. Here, we propose a proxy for thermal pressurization implemented through a modification of the rate-and-state friction law. This approach is designed to emulate the fault weakening and the relationship between breakdown energy and slip resulting from thermal pressurization and is appropriate for fully-dynamic simulations of multiple earthquake cycles. It preserves the computational efficiency of conventional rate-and-state friction models, which in turn can enable systematic studies to advance our understanding of the effects of fault weakening on earthquake mechanics. In 2.5D simulations of pulse-like ruptures on faults with finite seismogenic width, based on our thermal pressurization proxy, we find that the spatial distribution of slip velocity near the rupture front is consistent with the conventional square-root singularity, despite continued slip-weakening within the pulse, once the rupture has propagated a distance larger than the rupture width. An unconventional singularity appears only at shorter rupture distances. We further derive and verify numerically a theoretical estimate of the breakdown energy dissipated by our implementation of thermal pressurization. These results support the use of fracture mechanics theory to understand the propagation and arrest of very large earthquakes.

SummaryThe Alaska-Aleutian subduction zone represents an ideal location to study dynamics within a mantle wedge. The subduction system spans several thousand kilometers, is characterized by a slab edge, and has ample seismicity. Additionally, the majority of islands along the arc house broadband seismic instruments. We examine shear wave splitting of local-S phases originating along the length of the subduction zone. We have dense measurement spacing in two regions, the central Aleutians and beneath Alaska. Beneath Alaska, we observe a rotation in fast splitting directions near the edge of the subducting slab. Fast directions change from roughly trench perpendicular away from the slab edge to trench parallel near the boundary. This is indicative of toroidal flow around the edge of the subducting Alaska slab. In the central Aleutians, local-S splitting is primarily oriented parallel to, or oblique to, the strike of the trench. The local-S measurements, however, exhibit a depth dependence where deeper events show more consistently trench parallel directions indicating prevalent trench parallel mantle flow. Our local-S shear wave splitting results suggest trench parallel orientation are likely present along much of the subduction zone excited by the slab edge, but that additional complexities exist along strike.