Conference Agenda

Many conjugate rifted margins portray a marked structural asymmetry, with the wide side characterized by the occurrence of a large offset low-angle detachment fault. One such structure, S, is imaged by 3D seismic data in West Iberia Margin, an archetypical magma-poor margin. Our understanding of S detachment formation and its role in generating large-scale asymmetry remains incomplete. Here, using a novel numerical modeling technique, we show that S consists of sequentially active high-angle faults propagate as narrow shear zones at depth. The shear zone active at a low angle of 20°- 25° in a restricted lower crustal channel overlying and merging with the Moho. Contrary to conventional ideas suggesting that S is the causing of asymmetry, we demonstrate that S detachment is the culmination of an asymmetry formation process that started earlier in the rifting history. Its sub-horizontal geometry is reproduced by the back-rotation of the shear-zones once inactive. We conclude that the distal configuration of magma-poor margins can be explained by the lateral coexistence of ductile low-angle shear zones leading to detachment formation, and brittle faulting providing water for mantle serpentinization.

The strength of faults is subject of an important debate throughout various Earth Scientific disciplines. Different scientific communities have different perspectives with respect to appropriate values for friction coefficients μ. Geodynamicists with a long-term perspective require very low effective strengths (μ<0.05), while at the same time realizing mountains need to be sustained as well. Geologists and seismologists typically start from Byerlee friction coefficients of 0.6<μ<0.85, whereas rock mechanics experiments at high seismic slip rates show short-term low dynamic friction values of 0.03<μ<0.3. Here I show that both long- and short-term approaches can be made more compatible through considering that a regional or global frictional strength should be approached as a strain-averaged quantity. Doing this accounts for large variations of strain in both time and space. What matters for large-scale models is that most deformation occurs over a very small space and time during which friction is exceptionally low, thus making the representative long-term strength low. This is supported by seismo-thermo-mechanical models that self-consistently simulate the dynamics of both long-term subduction and short-term seismogenesis. The latter sustain mountain building, while representative earthquake-like events occur on faults with pore fluid pressure-effective static friction coefficients between 0.125 and 0.005 (or 0.75<Pf/Ps<0.99). These low friction values suggest faults are weak and suggest the dominant role of fluid pressures in weakening faults in subduction zones. This is confirmed in analytical considerations based on mechanical energy dissipation, which provide an equation to calculate the long-term fault strength as a strain-average quantity. Constraining the four parameters in this equation by observations confirms that fluid weakening is more important for long-term weakening than dynamic frictional weakening and low static friction coefficients. From the short-term perspective of modeling earthquake rupture dynamics it is now also becoming evident that fluid overpressured faults are preferable. They namely facilitate the incorporation of laboratory-observed dynamic weakening (70-90%) by limiting the stress drop to reasonable values. In summary, this cross-scale perspective supports long-term effective friction values in the range of about 0.03 to 0.2.

Induced seismicity as a result of natural gas production is a major challenge from both an industrial and a societal perspective. The compaction caused by gas production leads to changes of the effective pressure fields in the reservoir and stress redistributions occur particularly in the surrounding faults. In addition, the strong coupling between fluid flow and solid rock deformations and the role of fluid flow regarding the frictional properties of the faults necessitate a coupled and comprehensive modeling framework. A general and fully coupled thermo-hydro-mechanical finite difference formulation is developed herein and the results are verified against numerical benchmarks. A visco-elasto-plastic rheological behavior is assumed for the bulk material and a return-mapping algorithm is implemented for accurate simulation of the stress evolution. The geometrical features of the faults are incorporated into a regularized continuum framework, while the response of the fault zone is governed by a rate-and-state-dependent friction model. Numerical simulations are provided for large-scale problems and their efficiency is assured through the evaluation of the consistently linearized systems of equations along with the use of advanced numerical solvers and parallel computing. Although the proposed framework is a step towards the modeling of earthquake sequences for induced seismicity applications, the features of the numerical model are highlighted for other applications, including seismic events in subduction settings where the role of fluid flow inside the faults is considerable. Another application of the present, fully coupled hydro-thermo-mechanical formulation is the prediction of the fluid pressurization phenomena, where the frictional heating increases the magnitude of the pore fluid pressure inside the faults, and the resultant degradation of dynamic frictional strength is naturally captured.

Seismic cycles reflect the internal periodicity of the dynamic processes on faults and their surrounding materials. Understanding seismic cycles is fundamental for the research of earthquakes and may ultimately help to better assess long-term seismic hazard. Numerical models are well-suited to overcome limited spatiotemporal observations and improve our understanding on this topic. With the development of supercomputers and numerical techniques to reduce running time, we are able to model seismic cycles in larger models which contain more complexity. However, large models in 3D are still computational time and memory consuming. Moreover, this may not be optimal if the aspects of lateral or depth variations within the results are not needed to answer a particular objective. Usually it’s also hard to determine which dimension plays the most important role in a certain research question. This inspired us to investigate the advantages and limitations of various dimensional models by simulating seismic cycles in 0D, 1D, 2D and 3D.

We use a C++ numerical library GARNET which allows us to conduct numerical experiments of a strike-slip fault under rate-and-state friction, surrounded by an elastic medium with constant tectonic loading and, test them under different parameters and initial conditions. By adding dimensions, we figure a more detailed structure of the seismic cycle. The higher dimensional models show quantitative and qualitative similarity and difference to the lower dimensional ones thus present both the validity and the limitations of them. For example, the inertia waves are not possible to present in 0D while quasi-dynamic radiation damping term can be added here instead. Another example is that due to lack of grid extension along the fault, both 0D and 1D model fail to reveal earthquake nucleation phase. However, some observables such as the seismic cycle period, maximum/minimum stress and slip rates are calculated accurately in lower dimensional models which are much faster than higher dimensional models. To extensively understand the effect of the newly added dimension in higher dimensional models, we introduced non-uniform grid configuration and compared the models under variable resolution. Our results in 2D shows that only a thin layer near the fault need to be modeled in a fine resolution with coarser grid configured far away. We also implemented and compared quasi- and fully dynamic models in the same way. Our results indicate that both the size of simulated seismic events and their interval are reduced in quasi-dynamic models. This could provide us with guidance to identify the appropriate model complexity for various problems. Finally, these results are compared to other participating codes for the SCEC SEAS benchmarks.

Induced seismicity, as has occurred in Groningen over the past decades, is a complex phenomenon whose causes and processes are not well understood. Data assimilation provides a means to combine theory and observations to estimate key state variables, in this case velocity and shear stress. We investigate its relevance for tracking the evolution of a fault system with the ultimate goal of better understanding the processes involved in induced seismicity. This work is part of the InFocus project (An Integrated Approach for Estimating Fault Slip Occurrence) within the DeepNL research program aiming to better understand the consequences of human activities in the deep subsurface. Our framework consists of an Ensemble Kalman Filter (EnKF) implemented in a 1D model. We impose rate-and-state friction (RSF) behavior at one extreme of the model representing a stick-slip fault and assimilate observations of shear stress and slip rate measured in the homogenous elastic medium. We develop our framework using the Parallel Data Assimilation Framework (PDAF) of the Alfred Wegener Institute for a more efficient and faster implementation. The results suggest that ensemble data assimilation can estimate shear stresses and slip rates acting on the fault and can thus be used for forecasting events in this perfect model setup. Future experiments will investigate the application of data assimilation in a simulation of a laboratory setting that represents fault slip in Groningen. The ultimate goal is to extend this to a more realistic case that mimics the earthquakes of the Zeerijp area.

Splay faults in subduction zones pose a tsunami hazard, as they can accommodate more vertical displacements than the megathrust and are situated closer to the coast. Modelling dynamic splay fault rupture is complicated because of the lack of constraints on self-consistent initial conditions on both the splay and the megathrust in a dynamic rupture model.

Here, we build on a two-dimensional modelling framework first presented in Van Zelst et al. (2019) that considers the different temporal and spatial scales of the geodynamics of subduction, seismic cycles, dynamic ruptures, and tsunamis. A geodynamic seismic cycle model provides the dynamic rupture model with self-consistent initial stress and strength conditions, material properties, and fault geometries for the megathrust and six blind splay faults. The resulting surface displacements of the dynamic rupture model then serve as input for a tsunami propagation and inundation model.

We find that all six splay faults are activated when the megathrust ruptures by the main rupture front, dynamic stress transfer or stress changes induced by reflected waves from the surface. Splay fault rupture results in distinct peaks in the vertical surface displacements with a smaller wavelength and larger amplitudes. This translates into the resulting tsunami which has multiple, high-amplitude crests.

Our results suggest that larger-than-expected tsunamis could be attributed to rupture on large splay faults. Considering the significant effects splay fault rupture can have on a tsunami, it is important to understand splay fault activation and to consider them in hazard assessment.

Active faults are usually surrounded by narrow regions of localized deformation extending several hundred meters in width across the fault. This region of deformation consisting of a dense fracture network is macroscopically viewed as an elastic layer with low seismic wave velocities and referred to as a fault damage zone. The strength of the fault damage zone evolves throughout the seismic cycle, but the details of the mechanism and the nature of this evolution remain elusive. Understanding the structural evolution of this fault damage zone is key to unraveling the location, recurrence, the stressing history, and the probability of subsequent earthquakes in an active fault zone.

We use numerical simulations to understand the effects of the damage accumulation through multiple earthquakes and the rate at which the fault damage zone regains strength during the interseismic periods. We focus on the changes in the crack density of the fault damage zone over multiple earthquake sequences. We model this crack density change of the damage accumulation and the interseismic healing as changes in the shear wave velocity of an elastic layer surrounding a strike-slip fault. Using observations from Wenchuan, Landers, and Nojima, we constrain the changes in shear wave velocity and the rate of interseismic healing. We use two-dimensional earthquake cycle models of strike-slip faults with mode III rupture where the displacement is out of the plane of interest and stresses and friction vary along the depth. The fault damage zone is modeled as an elastic layer with a lower shear modulus compared to the surrounding host rock. We parameterize the shear modulus change using three variables: (1) the coseismic damage accumulation, which shows the amount of damage increase after an earthquake, (2) the healing time, which shows the interseismic duration it takes the fault zone to heal to its maximum level, and (3) the permanent damage, which shows the amount of damage that the fault zone never recovers. We will show the effect of this damage evolution on the recurrence and the surface slip expression of earthquake sequences. These simulations will provide a better insight into the partitioning of damage and healing during seismic cycles and the saturation of damage in mature fault zones.

Understanding spontaneous porous flow localisation in Earth's subsurface is of major importance for geological and industrial applications. These conductive pathways are of major importance to explain fast vertical fluid transfers and match observations provided by high-resolution geophysical monitoring of the subsurface. If the porous matrix undergoes time-dependent deformation, solitary waves of porosity can occur. Porosity then increases at the tip of the upwards travelling fluid pocket and decreases behind it. Resolving the underlying coupled and nonlinear hydro-mechanical processes is vital to account for accurate predictions of the time-dependent porosity modifications. Here, we apply an inversion approach that allows us to invert for the porosity distribution based on fluid velocity as observable. We use the adjoint method to compute the pointwise gradient which is used in a gradient descent framework to invert for the nodal porosity. Technically, the forward and the adjoint hydro-mechanical equations are solved using an implicit iterative pseudo-transient scheme on a Cartesian finite-difference grid. The solver executes in parallel on multiple graphical processing units at optimal and native performance relying on the Julia programming language.