ALBERT

Description: Summary On 24th August 2016 at 01:36 UTC a ML6.0 earthquake struck several villages in central Italy, among which Accumoli, Amatrice and Arquata del Tronto. The earthquake was recorded by about 350 seismic stations, causing 299 fatalities and damage with macroseismic intensities up to 11. The maximum acceleration was observed at Amatrice station (AMT) reaching 916 cm/s2 on E-W component, with epicentral distance of 15 km and Joyner and Boore distance to the fault surface (RJB) of less than a kilometre. Motivated by the high levels of observed ground motion and damage, we generate broadband seismograms for engineering purposes by adopting a hybrid method. To infer the low frequency seismograms, we considered the kinematic slip model by Tinti et al. (2016). The high frequency seismograms were produced using a stochastic finite-fault model approach based on dynamic corner-frequency. Broadband synthetic time series were therefore obtained by merging the low and high frequency seismograms. Simulated hybrid ground motions were compared both with the observed ground motions and the ground-motion prediction equations (GMPEs), to explore their performance and to retrieve the region-specific parameters endorsed for the simulations. In the near-fault area we observed that hybrid simulations have a higher capability to detect near source effects and to reproduce the source complexity than the use of GMPEs. Indeed, the general good consistency found between synthetic and observed ground motion (both in the time and frequency domain), suggests that the use of regional-specific source scaling and attenuation parameters together with the source complexity in hybrid simulations improves ground motion estimations. To include the site effect in stochastic simulations at selected stations, we tested the use of amplification curves derived from HVRSs (horizontal-to-vertical response spectra) and from HVSRs (horizontal-to-vertical spectral ratios) rather than the use of generic curves according to NTC-18 Italian seismic design code. We generally found a further reduction of residuals between observed and simulated both in terms of time histories and spectra.

Description: Summary We explore here the benefits of using constraints from seismic tomography in gravity data inversion and how inverted density distributions can be improved by doing so. The methodology is applied to a real field case in which we reconstruct the density structure of the Pyrenees along a southwest-northeast transect going from the Ebro basin in Spain to the Arzacq basin in France. We recover the distribution of densities by inverting gravity anomalies under constraints coming from seismic tomography. We initiate the inversion from a prior density model obtained by scaling a pre-existing compressional seismic velocity Vp model using a Nafe-Drake relationship : the Vp model resulting from a full-waveform inversion of teleseismic data. Gravity data inversions enforce structural similarities between Vp and density by minimizing the norm of the cross-gradient between the density and Vp models. We also compare models obtained from 2.5D and 3D inversions. Our results demonstrate that structural constraints allow us to better recover the density contrasts close to the surface and at depth, without degrading the gravity data misfit. The final density model provides valuable information on the geological structures and on the thermal state and composition of the western region of the Pyrenean lithosphere.

Description: Summary The relatively short duration of the early stages of subduction results in a poor geological record, limiting our understanding of this critical stage. Here, we utilize a 2D numerical model of incipient subduction, that is the stage after a plate margin has formed with a slab tip that extends to a shallow depth and address the conditions under which subduction continues or fails. We assess energy budgets during the evolution from incipient subduction to either a failed or successful state, showing how the growth of potential energy, and slab pull, is resisted by the viscous dissipation within the lithosphere and the mantle. The role of rheology is also investigated, as deformation mechanisms operating in the crust and mantle facilitate subduction. In all models, the onset of subduction is characterized by high lithospheric viscous dissipation and low convergence velocities, whilst successful subduction sees the mantle become the main area of viscous dissipation. In contrast, failed subduction is defined by the lithospheric viscous dissipation exceeding the lithospheric potential energy release rate and velocities tend towards zero. We show that development of a subduction zone depends on the convergence rate, required to overcome thermal diffusion and to localise deformation along the margin. The results propose a minimum convergence rate of ∼ 0.5 cm yr−1 is required to reach a successful state, with 100 km of convergence over 20 Myr, emphasizing the critical role of the incipient stage.

Description: Summary Although many studies have revealed that the atmospheric effects of electromagnetic wave propagation (including ionospheric and tropospheric water vapor) have serious impacts on Interferometric Synthetic Aperture Radar (InSAR) measurement results, atmospheric corrections have not been thoroughly and comprehensively investigated in many well-known cases of InSAR focal mechanism solutions, which means there is no consensus on whether atmospheric effects will affect the InSAR focal mechanism solution. Moreover, there is a lack of quantitative assessment on how much the atmospheric effect affects the InSAR focal mechanism solution. In this paper, we emphasized that it was particularly important to assess the impact of InSAR ionospheric and tropospheric corrections on the underground nuclear explosion modeling quantitatively. Therefore, we investigated the 4th North Korea (NKT-4) underground nuclear test using ALOS-2 liters-band SAR images. Because the process of the underground nuclear explosion was similar to the volcanic magma source activity, we modeled the ground displacement using the Mogi model. Both the ionospheric and tropospheric phase delays in the interferograms were investigated. Furthermore, we studied how the ionosphere and troposphere phase delays could bias the estimation of Mogi source parameters. The following conclusions were drawn from our case study: The ionospheric delay correction effectively mitigated the long-scale phase ramp in the full-frame interferogram, the standard deviation decreased from 1.83 cm to 0.85 cm compared to the uncorrected interferogram. The uncorrected estimations of yield and depth were 8.44 kt and 370.33 m, respectively. Compared to the uncorrected estimations, the ionospheric correction increased the estimation of yield and depth to 9.43 kt and 385.48 m while the tropospheric correction slightly raised them to 8.78 kt and 377.24 m. There were no obvious differences in the location estimations among the four interferograms. When both corrections were applied, the overall standard deviation was 1.16 cm, which was even larger than the ionospheric corrected interferogram. We reported the source characteristics of NKT-4 based on the modeling results derived from the ionospheric corrected interferogram. The preferred estimation of NKT-4 was a Mogi source located at 129°04′22.35‘E, 41°17′54.57″N buried at 385.48 m depth. The cavity radius caused by the underground explosion was 22.02 m. We reported the yield estimation to be 9.43 kt. This study showed that for large-scale natural deformation sources such as volcanoes and earthquakes, atmospheric corrections would be more significant, but even if the atmospheric signal did not have much complexity, the corrections should not be ignored.

Description: SUMMARY Self-consistent modelling of magmatic systems is challenging as the melt continuously changes its chemical composition upon crystallization, which may affect the mechanical behaviour of the system. Melt extraction and subsequent crystallization create new rocks while depleting the source region. As the chemistry of the source rocks changes locally due to melt extraction, new calculations of the stable phase assemblages are required to track the rock evolution and the accompanied change in density. As a consequence, a large number of isochemical sections of stable phase assemblages are required to study the evolution of magmatic systems in detail. As the state-of-the-art melting diagrams may depend on nine oxides as well as pressure and temperature, this is a 10-D computational problem. Since computing a single isochemical section (as a function of pressure and temperature) may take several hours, computing new sections of stable phase assemblages during an ongoing geodynamic simulation is currently computationally intractable. One strategy to avoid this problem is to pre-compute these stable phase assemblages and to create a comprehensive database as a hyperdimensional phase diagram, which contains all bulk compositions that may emerge during petro-thermomechanical simulations. Establishing such a database would require repeating geodynamic simulations many times while collecting all requested compositions that may occur during a typical simulation and continuously updating the database until no additional compositions are required. Here, we describe an alternative method that is better suited for implementation on large-scale parallel computers. Our method uses the entries of an existing preliminary database to estimate future required chemical compositions. Bulk compositions are determined within boundaries that are defined manually or through principal component analysis in a parameter space consisting of clustered database entries. We have implemented both methods within a massively parallel computational framework while utilizing the Gibbs free energy minimization program Perple_X. Results show that our autonomous approach increases the resolution of the thermodynamic database in compositional regions that are most likely required for geodynamic models of magmatic systems.

Description: SUMMARY We present a statistical rock physics inversion of the elastic and electrical properties to estimate the petrophysical properties and quantify the associated uncertainty. The inversion method combines statistical rock physics modeling with Bayesian inverse theory. The model variables of interest are porosity and fluid saturations. The rock physics model includes the elastic and electrical components and can be applied to the results of seismic and electromagnetic inversion. To describe the non-Gaussian behaviour of the model properties, we adopt non-parametric probability density functions to sample multimodal and skewed distributions of the model variables. Different from machine learning approach, the proposed method is not completely data-driven but is based on a statistical rock physics model to link the model parameters to the data. The proposed method provides pointwise posterior distributions of the porosity and CO2 saturation along with the most-likely models and the associated uncertainty. The method is validated using synthetic and real data acquired for CO2 sequestration studies in different formations: the Rock Springs Uplift in Southwestern Wyoming and the Johansen formation in the North Sea, offshore Norway. The proposed approach is validated under different noise conditions and compared to traditional parametric approaches based on Gaussian assumptions. The results show that the proposed method provides an accurate inversion framework where instead of fitting the relationship between the model and the data, we account for the uncertainty in the rock physics model.

Description: SUMMARY The horizontal-to-vertical spectral ratio (HVSR) of ambient noise is commonly used to infer a site's resonance frequency (${f_{0,site}}$). HVSR calculations are performed most commonly using the Fourier amplitude spectrum obtained from a single merged horizontal component (e.g. the geometric mean component) from a three-component sensor. However, the use of a single merged horizontal component implicitly relies on the assumptions of azimuthally isotropic seismic noise and 1-D surface and subsurface conditions. These assumptions may not be justified at many sites, leading to azimuthal variability in HVSR measurements that cannot be accounted for using a single merged component. This paper proposes a new statistical method to account for azimuthal variability in the peak frequency of HVSR curves (${f_{0,HVSR}}$). The method uses rotated horizontal components at evenly distributed azimuthal intervals to investigate and quantify azimuthal variability. To ensure unbiased statistics for ${f_{0,HVSR}}$ are obtained, a frequency-domain window-rejection algorithm is applied at each azimuth to automatically remove contaminated time windows in which the ${f_{0,HVSR}}$ values are statistical outliers relative to those obtained from the majority of windows at that azimuth. Then, a weighting scheme is used to account for different numbers of accepted time windows at each azimuth. The new method is applied to a data set of 114 HVSR measurements with significant azimuthal variability in ${f_{0,HVSR}}$, and is shown to reliably account for this variability. The methodology is also extended to the estimation of a complete lognormal-median HVSR curve that accounts for azimuthal variability. To encourage the adoption of this statistical approach to accounting for azimuthal variability in single-station HVSR measurements, the methods presented in this paper have been incorporated into hvsrpy, an open-source Python package for HVSR processing.

Description: SUMMARY This study focuses in the analysis of the internal structure of the upper 3 km of Los Humeros (LH) caldera and the relation of electrical and hydrothermal anomalies. For this purpose, we measured, processed and interpreted 78 broad-band magnetotelluric (MT) soundings. We performed a 3-D inversion of the data set (ModEM) using all MT soundings, although only half of the available frequencies per sounding due to limited computed power. We also carried out the 2-D inversions (NLCG) of the invariant determinant along two orthogonal profiles (EW and NS) crossing the caldera structure; their comparison yields similar resistivity and structural models results. The resistivity modelling is complemented with the results of a joint 3-D inversion of an accurate gravity database of 720 stations, and total field aeromagnetic data (SGM) from the caldera crater. The combined results provide novel details about the structure of the shallow geothermal reservoir of the resurgence caldera complex hosting the active hydrothermal system. Density and resistivity models show the existence of a composed crater basin structure separated by an EW high-density structure; the northern basin is associated to the LH crater, whereas the southern basin associates to the emergent Los Potreros (LP) caldera basin. The magnetization model indicates that there is a common source for the magnetic volcanic products observed at the caldera surface, and that the LP fault is the more magnetized fault of the geothermal system. The propylic zoning under the geothermal field, which according to the MT model results has resistivities above ∼100 Ω-m, was extrapolated using this and additional criteria to obtain the distribution of other hypothetical propylitic zones of hydrothermal potential.

Description: SUMMARY For the time stationary global geomagnetic field, a new modelling concept is presented. A Bayesian non-parametric approach provides realistic location dependent uncertainty estimates. Modelling related variabilities are dealt with systematically by making little subjective a priori assumptions. Rather than parametrizing the model by Gauss coefficients, a functional analytic approach is applied. The geomagnetic potential is assumed a Gaussian process to describe a distribution over functions. A priori correlations are given by an explicit kernel function with non-informative dipole contribution. A refined modelling strategy is proposed that accommodates non-linearities of archeomagnetic observables: First, a rough field estimate is obtained considering only sites that provide full field vector records. Subsequently, this estimate supports the linearization that incorporates the remaining incomplete records. The comparison of results for the archeomagnetic field over the past 1000 yr is in general agreement with previous models while improved model uncertainty estimates are provided.

Description: SUMMARY Rapid development of time-lapse seismic monitoring instrumentations has made it possible to collect dense time-lapse data for tomographically retrieving time-lapse (even continuous) images of subsurface changes. While traditional time-lapse full waveform inversion (TLFWI) algorithms are designed for sparse time-lapse surveys, they lack of effective temporal constraint on time-lapse data, and, more importantly, lack of the uncertainty estimation of the TLFWI results that is critical for further interpretation. Here, we propose a new data assimilation TLFWI method, using hierarchical matrix powered extended Kalman filter (HiEKF) to quantify the image uncertainty. Compared to existing Kalman filter algorithms, HiEKF allows to store and update a data-sparse representation of the cross-covariance matrices and propagate model errors without expensive operations involving covariance matrices. Hence, HiEKF is computationally efficient and applicable to 3-D TLFWI problems. Then, we reformulate TLFWI in the framework of HiEKF (termed hereafter as TLFWI-HiEKF) to predict time-lapse images of subsurface spatiotemporal velocity changes and simultaneously quantify the uncertainty of the inverted velocity changes over time. We demonstrate the validity and applicability of TLFWI–HiEKF with two realistic CO2 monitoring models derived from Frio-II and Cranfield CO2 injection sites, respectively. In both 2-D and 3-D examples, the inverted high-resolution time-lapse velocity results clearly reveal a continuous velocity reduction due to the injection of CO2. Moreover, the accuracy of the model is increasing over time by assimilating more time-lapse data while the standard deviation is decreasing over lapsed time. We expect TLFWI-HiEKF to be equipped with real-time seismic monitoring systems for continuously imaging the distribution of subsurface gas and fluids in the future large-scale CO2 sequestration experiments and reservoir management.

Description: SUMMARY Seven years after the beginning of a massive wastewater injection project in eastern Colombia, local earthquake activity increased significantly. The field operator and the Colombian Geological Survey immediately reinforced the monitoring of the area. Our analysis of the temporal evolution of the seismic and injection data together with our knowledge of the geological parameters of the region indicate that the surge of seismicity is being induced by the re-injection of produced water into the same three producing reservoirs. Earthquake activity began on known faults once disposal rates had reached a threshold of ∼2 × 106 m3 of water per month. The average reservoir pressure had remained constant at 7.6 MPa after several years of production, sustained by a large, active aquifer. Surface injection pressures in the seismically active areas remain below 8.3 MPa, a value large enough to activate some of the faults. Since faults are mapped throughout the region and many do not have seismicity on them, we conclude that the existence of known faults is not the only control on whether earthquakes are generated. Stress conditions of these faults are open to future studies. Earthquakes are primarily found in four clusters, located near faults mapped by the operator. The hypocentres reveal vertical planes with orientations consistent with focal mechanisms of these events. Stress inversion of the focal mechanisms gives a maximum compression in the direction ENE-WSW, which is in agreement with borehole breakout measurements. Since the focal mechanisms of the earthquakes are consistent with the tectonic stress regime, we can conclude that the seismicity is resulting from the activation of critically stressed faults. Slip was progressive and seismic activity reached a peak before declining to few events per month. The decline in seismicity suggests that most of the stress has been relieved on the main faults. The magnitude of a large majority of the recorded earthquakes was lower than 4, as the pore pressure disturbance did not reach the mapped large faults whose activation might have resulted in larger magnitude earthquakes. Our study shows that a good knowledge of the local fault network and conditions of stress is of paramount importance when planning a massive water disposal program. These earthquakes indicate that while faults provide an opportunity to dispose produced water at an economically attractive volume–pressure ratio, the possibility of induced seismicity must also be considered.

Description: Summary Virtual Deep Seismic Sounding (VDSS) uses the arrival time of post-critical SsPmp relative to the direct S wave to infer Moho depth at the Pmp reflection point. Due to the large offset between the virtual source and the receiver, SsPmp is more sensitive to lateral variations of structures than near-vertical phases such as Ps, which is used to construct conventional P receiver functions. However, the way post-critical SsPmp is affected by lateral variations in lithospheric structure is not well understood, and previous studies largely assumed a 1D structure when analyzing SsPmp waveforms. Here we present synthetic tests with various 2D models to show that lateral variations in lithospheric structures, from the lithosphere-asthenosphere boundary to sedimentary basins, profoundly affect travel time, phase, and amplitude of post-critical SsPmp, and that a 1D approximation is usually inappropriate when analyzing 2D data. Despite these strong effects we show, with synthetic examples and the ChinArray data from the Ordos Block in northern China, that a simple ray-theory-based back-projection method can retrieve the geometry of the crust-mantle boundary given array observations in cases with moderate lateral variations in the crust-mantle boundary and/or the lithosphere-asthenosphere boundary. The success of our back-projection method indicates that ray-theory approximations are sufficient in modeling SsPmp travel times in the presence of moderate lateral heterogeneity. In contrast, we show that the ray theory is generally insufficient in modeling SsPmp phase shifts in a strongly heterogeneous lithosphere due to non-planar down-going P waves incident at the crust-mantle boundary. Nonetheless, our results demonstrate the feasibility of direct imaging of the crust-mantle boundary with post-critical SsPmp even in the presence of 2D variations of lithospheric structure.

Description: Summary The Mw 6.0 earthquake that hit central Italy on 24 August 2016 caused an abnormally high level of destruction in the town of Amatrice. In order to clarify the role of site response in causing such a disaster, a series of ambient noise recordings acquired in the aftermath of the event are analysed here to identify site resonance properties from the ratios H/V between horizontal and vertical amplitudes of ground motion. Although the noise data acquisition was limited by the emergency management activities, the use of a new analysis technique, which isolates the contribution of Rayleigh waves to the noise wavefield and averages instantaneous estimates of H/V ratios, provided more stable results compared to the standard Nakamura's technique based on mean spectral ratios. The results demonstrated the occurrence of significant resonance phenomena, but without an obvious correlation with the spatial distribution of damage severity. It is apparent that the damage severity was also influenced by some additional local factors related to building vulnerability. Moreover, the time series analysis revealed seasonal variations in the Rayleigh wave ellipticity curves likely related to the water content changes in the surface deposits and their influence on the Poisson coefficient. Finally, the new method proved capable of recognizing time-varying directions of Rayleigh wave propagation. This capability could be exploited to support other passive seismic methods (e.g. ReMi), whose results reliability is limited by the lack of control on wave direction origin.

Description: SUMMARY This paper explores the effects on tsunami simulations of the level of detail of the bathymetric grid in use. For this purpose, we expand available bathymetric data sets of the Pacific Basin in spherical harmonics. For realistic scenarios of tsunamis generated by earthquake dislocations, we conclude that an expansion to a maximum degree lmax = 40, corresponding to wavelengths of 1000 km, is sufficient to reproduce the main features of the tsunami wavefield synthesized in deep water, that is, without considering final shoaling and interaction with coastal features.

Description: SUMMARY Full-waveform inversion (FWI) has demonstrated increasing success in estimating medium properties, but its computational cost still poses challenges in moving towards high-resolution imaging of targets at depth. Here, we propose a target-oriented FWI method that inverts for the medium parameters confined within an arbitrary region of interest. Our method is novel in terms of both local wavefield modelling and data redatuming, in order to build a target-oriented objective function which is sensitive to the target medium only without further assumptions about the medium outside. Based on the convolution-type representation theorem, our local forward modelling operator propagates wavefields within the target medium only while providing full acoustic coupling between the target medium and the surrounding geology. A key requirement of our local FWI method is that the subsurface wavefields surrounding and inside the target be as accurate as possible. As such, the subsurface wavefields are retrieved by the Marchenko method, which can redatum the single-sided surface reflection data to the target zone while preserving both primary and multiple reflections, with minimal a priori knowledge of the full-domain medium. Given a sufficiently accurate initial velocity macromodel, our numerical examples show that our local FWI method resolves the reservoir zone of a 2-D Barrett Unconventional P-wave velocity model much more efficiently than the conventional full-domain FWI without significantly sacrificing accuracy. Our method may further enable FWI approaches to high-resolution imaging of subsurface targets.

Description: SUMMARY Determining rock microstructure remains challenging, since a proper rock-physics model is needed to establish the relation between pore microstructure and elastic and transport properties. We present a model to estimate pore microstructure based on porosity, ultrasonic velocities and permeability, assuming that the microstructure consists on randomly oriented stiff equant pores and penny-shaped cracks. The stiff pore and crack porosity varying with differential pressure is estimated from the measured total porosity on the basis of a dual porosity model. The aspect ratio of pores and cracks and the crack density as a function of differential pressure are obtained from dry-rock P- and S-wave velocities, by using a differential effective medium model. These results are used to invert the pore radius from the matrix permeability by using a circular pore model. Above a crack density of 0.13, the crack radius can be estimated from permeability, and below that threshold, the radius is estimated from P-wave velocities, taking into account the wave dispersion induced by local fluid flow between pores and cracks. The approach is applied to experimental data for dry and saturated Fontainebleau sandstone and Chelmsford Granite.