ALBERT

Description: Summary Studies of glacial isostatic adjustment (GIA) provide important constraints on the Earth's mantle viscosity. Most GIA models assume Newtonian viscosity through the mantle, but laboratory experimental studies of rock deformation, observational studies of seismic anisotropy, and modeling studies of mantle dynamics show that in the upper mantle non-Newtonian viscosity may be important. This study explores the non-Newtonian effects on the GIA induced variations in mantle stress and viscosity and on surface observables including vertical displacement, relative sea level (RSL) and gravity change. The recently updated and fully benchmarked software package CitcomSVE is used for GIA simulations. We adopt the ICE-6G ice deglaciation history, VM5a lower mantle and lithospheric viscosities, and a composite rheology that combines Newtonian and non-Newtonian viscosities for the upper mantle. Our results show that: 1) The mantle stress beneath glaciated regions increases significantly during deglaciation, leading to regionally reduced upper mantle viscosity by more than an order of magnitude. Such effects can be rather localized at the periphery of glaciated regions. However, non-Newtonian effects on far-field mantle viscosity are negligibly small. GIA induced stress is also significant in the lithosphere (∼30 MPa) and lower mantle (∼2 MPa). 2) The predicted RSL changes from non-Newtonian models display distinct features in comparison with the Newtonian model, including more rapid sea-level falls associated with the rapid deglaciation at ∼14,000 years ago followed by a more gradual sea-level variation for sites near the centers of formerly glaciated regions, and an additional phase of sea-level falls for the last ∼8000 years for sites at the ice margins. Similar time-dependence associated with the deglaciation is also seen for rate of vertical displacement, suggesting a relatively slow present-day rates of vertical displacement and gravity change. These features can be explained by the non-Newtonian effects associated with a loading event which manifest a fast relaxation stage followed by a relative slow relaxation stage. Our results may provide GIA diagnoses for distinguishing non-Newtonian and Newtonian rheology.

Description: Summary Northeastward subduction of the oceanic Rivera and Cocos plates in western Mexico poses a poorly understood seismic hazard to the overlying areas of the North America plate. We estimate the magnitude and distribution of interseismic locking along the northern ∼500 km of the Mexico subduction zone, with a series of elastic half-space inversions that optimize the fits to the velocities of 57 GPS stations in western Mexico. All velocities were corrected for the coseismic, afterslip and viscoelastic rebound effects of the 1995 Colima-Jalisco and 2003 Tecomán earthquakes. We explore the robustness of interseismic locking estimates to a variety of mantle Maxwell times that are required for the viscoelastic corrections, to the maximum permitted depth for locking of the subduction interface, and to the location assigned to the Rivera-Cocos-North America plate triple junction offshore from western Mexico. The best fitting locking solutions are associated with a maximum locking depth of 40 km, a triple junction location ∼50 km northwest of the Manzanillo Trough, and a mantle Maxwell time of 15 yr (viscosity of 2 × 1019 Pa·s). Checkerboard tests show that the locking distribution is best resolved at intermediate depths (10-40 km). All of our inversions define a gradual transition from strong locking (i.e. 70-100 percent) of most (70%) of the Rivera-North America subduction interface to strong but less uniform locking below the Manzanillo Trough, where oceanic lithosphere transitional between the Cocos and Rivera plate subducts, to weak to moderate locking (averaging 55 percent) of the Michoacán segment of the Cocos-North America interface. Strong locking of the ∼125-km-long trench segment offshore from Puerto Vallarta and other developed coastal areas, where our modelling indicates an average annual elastic slip-rate deficit of ∼20 mm yr−1, implies that ∼1.8 m of unrelieved plate slip has accrued since the segment last ruptured in 1932, sufficient for a M ∼ 8.0 earthquake.

Description: Summary We invert ∼25 years of campaign and continuous Global Positioning System daily positions at 62 sites in southwestern Mexico to estimate coseismic and postseismic afterslip solutions for the 1995 Mw = 8.0 Colima-Jalisco and the 2003 Mw = 7.5 Tecomán earthquakes, and the long-term velocity of each GPS site. Estimates of the viscoelastic effects of both earthquakes from a 3-D model with an elastic crust and subducting slab, and linear Maxwell viscoelastic mantle are used to correct the GPS position time series prior to our time-dependent inversions. The preferred model, which optimizes the fit to data from several years of rapid postseismic deformation after the larger 1995 earthquake, has a mantle Maxwell time of 15 years (viscosity of 2 × 1019 Pa·s), although upper mantle viscosities as low as 5 × 1018 Pa·s cannot be excluded. Our geodetic slip solutions for both earthquakes agree well with previous estimates derived from seismic data or via static coseismic offset modelling. The afterslip solutions for both earthquakes suggest that most afterslip coincided with the rupture areas or occurred farther downdip, and had cumulative moments similar to or larger than the coseismic moments. Afterslip thus appears to relieve significant stress along the Rivera plate subduction interface, including the area of the interface between a region of deep non-volcanic tremor and the shallower seismogenic zone. We compare the locations of the seismogenic zone, afterslip and tremor in our study area to those of the neighboring Guerrero and Oaxaca segments of the Mexico subduction zone. Our newly derived interseismic GPS site velocities, the first for western Mexico that are corrected for the coseismic and postseismic effects of the 1995 and 2003 earthquakes, are essential for future estimates of the interseismic subduction interface locking and hence the associated seismic hazard.

Description: Summary Ferrimagnetic, monoclinic 4C pyrrhotite (Fe7S8) is the only iron sulfide with high relevance for paleomagnetism and rock magnetism that can be identified in rock materials by its characteristic low-temperature anomaly. Despite its relevance in natural magnetism and the many magnetic studies over the last decades, the physics and the crystallography behind this anomaly, also denoted Besnus transition, is a matter of debate. In this study we analyze the static and dynamic magnetization associated with the Besnus transition in conjunction with low-temperature structural data of 4C pyrrhotite reported in the literature. The correlation between the Fe–Fe bonds causing spin-orbit coupling and the dynamic magnetic properties show that the magnetic characteristics of the Besnus transition stem from the interaction of two magnetocrystalline anisotropy systems triggered by thermally induced structural changes on an atomic level in monoclinic 4C pyrrhotite. This refutes the widespread view that the Besnus transition is caused by a crystallographic change from monoclinic to triclinic.

Description: Summary Two-phase flow equations that couple solid deformation and fluid migration have opened new research trends in geodynamical simulations and modelling of subsurface engineering. Physical nonlinearity of fluid-rock systems and strong coupling between flow and deformation in such equations lead to interesting predictions such as spontaneous formation of focused fluid flow in ductile/plastic rocks. However, numerical implementation of two-phase flow equations and their application to realistic geological environments with complex geometries and multiple stratigraphic layers is challenging. This study documents an efficient pseudo-transient solver for two-phase flow equations and describes the numerical theory and physical rationale. We provide a simple explanation for all steps involved in the development of a pseudo-transient numerical scheme for various types of equations. Two different constitutive models are used in our formulations: a bilinear viscous model with decompaction weakening and a viscoplastic model that allows decompaction weakening at positive effective pressures. The resulting numerical models are used to study fluid leakage from high porosity reservoirs into less porous overlying rocks. The interplay between time-dependent rock deformation and the buoyancy of ascending fluids leads to the formation of localized channels. The role of material parameters, reservoir topology, geological heterogeneity and porosity is investigated. Our results show that material parameters control the propagation speed of channels while the geometry of the reservoir controls their locations. Geological layers present in the overburden do not stop the propagation of the localized channels but rather modify their width, permeability, and growth speed.

Description: Summary Using data from 3837 seismic stations deployed in or around China, we construct high-resolution models of crustal thickness (H) and seismic compressional and shear velocity ratio (Vp/Vs or κ) in continental China by analysis of 1,150,543 receiver functions. We group the receiver functions in cells with a spatial resolution of 0.25˚ × 0.25˚ in the North-South China Seismic Belt and parts of the North China Craton, and of 0.5˚ × 0.5˚ in other regions, classify the receiver functions based on their characteristics, and develop a modified H-κ stacking method to construct models in the regions where the receiver functions are significantly affected by sedimentary basins and by Moho architecture. The inferred crustal thickness model displays an eastward thinning trend from the thickest crust (〉 80 km) beneath the Qiangtang Block to the thinnest crust (〈 26 km) beneath the southern part of the Cathaysia Block. Crustal thickness is 26 – 50 km in several major basins and 26 – 55 km in the Precambrian cratonic blocks. The inferred Vp/Vs model in the crystalline crust displays moderate-to-high values (1.75 – 1.85) in the southeastern margin of the Tibetan Plateau, the Tengchong volcanic field, the Emeishan large igneous province, the north-central areas of the Bohaiwan and Songliao basins, the western margin of the Taikang Hefei Basin and the southeastern margin of the Cathaysia Block. Lower values (≤ 1.72) characterize the major regions of the Cathaysia Block and the Jiangnan Orogenic Belt, and the hinterlands of the Ordos Block and Sichuan Basin. We discuss possible tectonic processes, secular crustal evolution and crustal compositions that are consistent with our inferred crustal thickness and Vp/Vs structure in continental China. This study establishes a framework of seismic data sharing for future studies in the seismological community in one of the first steps of developing a China Seismological Reference Model.

Description: Summary Traditional earthquake location relying on first arrival picking is challenging for microseismic events with low signal-to-noise ratio. Over the past years, alternative procedures have been explored based on the idea of migrating the energy of an earthquake back into its source position by stacking along theoretical traveltime curves. To avoid destructive interference of signals with opposite polarity, it is common to transform the input signals into positive timeseries. Stacking-based source location has been successfully applied at various scales, but existing studies differ considerably in the choice of characteristic function, the amount of pre-processing and the phases used in the analysis. We use a dataset of 62 natural microearthquakes recorded on a 2D seismic array of 145 vertical geophones across the glacially-triggered Burträsk fault to compare the performance of five commonly used characteristic functions: the noise filtered seismograms and the semblance, the envelope, the short term average/long term average ratio and the kurtosis gradient of the seismograms. We obtain the best results for a combined P- and S-wave location using a polarity-sensitive characteristic function, i.e., the filtered seismograms or the semblance. In contrast, the absolute functions often fail to align the signals properly, yielding biased location estimates. Moreover, we observe that the success of the procedure is very sensitive to noise suppression and signal shaping prior to stacking. Our study demonstrates the usefulness of including lower quality S-wave data to improve the location estimates. Furthermore, our results illustrate the benefits of retaining the phase information for location accuracy and noise suppression. To ensure optimal location results, we recommend carefully pre-processing the data and test different characteristic functions for each new dataset. In spite of the sub-optimal array geometry, we obtain good locations for most events within ∼30-40 km of the survey and the locations are consistent with an image of the fault trace from an earlier reflection seismic survey.

Description: Summary Solving the wave equation to obtain wavefield solutions is an essential step in illuminating the subsurface using seismic imaging and waveform inversion methods. Here, we utilize a recently introduced machine-learning based framework called physics-informed neural networks (PINNs) to solve the frequency-domain wave equation, which is also referred to as the Helmholtz equation, for isotropic and anisotropic media. Like functions, PINNs are formed by using a fully-connected neural network (NN) to provide the wavefield solution at spatial points in the domain of interest, in which the coordinates of the point form the input to the network. We train such a network by back propagating the misfit in the wave equation for the output wavefield values and their derivatives for many points in the model space. Generally, a hyperbolic tangent activation is used with PINNs, however, we use an adaptive sinusoidal activation function to optimize the training process. Numerical results show that PINNs with adaptive sinusoidal activation functions are able to generate frequency-domain wavefield solutions that satisfy wave equations. We also show the flexibility and versatility of the proposed method for various media, including anisotropy, and for models with strong irregular topography.

Description: Summary The elastoplastic Iwan model has been used since the end of the 1970s to simulate nonlinear soil behavior in seismic wave propagation. In this work, we present an automatic algorithm to efficiently sample the shear-modulus reduction curve in function of shear deformation, which constitutes the exclusive ingredient of the elastoplastic model. This model requires the data from the shear- modulus reduction as a function of shear deformation, which are readily available in the literature and from specific laboratory tests. The method involves a discretization and interpolation of these data to be used. The quality of the solution depends on the number of interpolated points. However, a larger number of them produces an increase of the computational time. To overcome this, we present an automatic algorithm to efficiently sample the shear-modulus reduction curve. We numerically prove that the chosen discretization of the curve has a strong impact on the calculation load, in addition to the well-known dependance on the input motion amplitude level. Two tests of nonlinear wave propagation in 1D and 3D media show the clear gain in computation time when using the proposed automatic sampling algorithm.

Description: Summary The geophysical detection of magma bodies and the estimation of the dimensions, physical properties, and the volume fraction of each phase composing the magma is required to improve the forecasting of volcanic hazards and to understand transcrustal magmatism. We develop an analytical model to calculate P waves velocity in a three-phase magma consisting of crystals and gas bubbles suspended in a viscous melt. We apply our model to calculate the speed of sound as a function of the temperature in three magmas with different chemical compositions, representative of the diversity that is encountered in arc magmatism. The model employs the coupled phase theory that explicitly accounts for the exchanges of momentum and heat between the phases. We show that the speed of sound varies non-linearly with the frequency of an acoustic perturbation between two theoretical bounds. The dispersion of the sound in a magma results from the exchange of heat between the melt and the dispersed phases that affects the magnitude of their thermal expansions. The lower bound of the sound speed occurs at low frequencies for which all the constituents can be considered in thermal equilibrium, whereas the upper bound occurs at high frequencies for which the exchange of heat between the phases may be neglected. The presence of gas in a magma produces a sharp decrease in the velocity of compressional waves and generates conditions in which the dispersion of the sound is significant at the frequencies usually considered in geophysics. Finally, we compare the estimates of our model with the ones from published relationships. Differences are largest at higher frequencies and are