Search

Detailed imaging of accretionary wedges reveals splay fault networks that could pose a significant tsunami hazard. However, the dynamics of multiple splay fault activation during megathrust earthquakes and the consequent effects on tsunami generation are not well understood. We use a 2‐D dynamic rupture model with complex topo‐bathymetry and six curved splay fault geometries constrained from realistic tectonic loading modeled by a geodynamic seismic cycle model with consistent initial stress and strength conditions. We find that all splay faults rupture coseismically. While the largest splay fault slips due to a complex rupture branching process from the megathrust, all other splay faults are activated either top down or bottom up by dynamic stress transfer induced by trapped seismic waves. We ascribe these differences to local non‐optimal fault orientations and variable along‐dip strength excess. Generally, rupture on splay faults is facilitated by their favorable stress orientations and low strength excess as a result of high pore‐fluid pressures. The ensuing tsunami modeled with non‐linear 1‐D shallow water equations consists of one high‐amplitude crest related to rupture on the longest splay fault and a second broader wave packet resulting from slip on the other faults. This results in two episodes of flooding and a larger run‐up distance than the single long‐wavelength (300 km) tsunami sourced by the megathrust‐only rupture. Since splay fault activation is determined by both variable stress and strength conditions and dynamic activation, considering both tectonic and earthquake processes is relevant for understanding tsunamigenesis., Plain Language Summary: In subduction zones, where one tectonic plate moves beneath another, earthquakes can occur on many different faults. Splay faults are relatively steep faults that branch off the largest fault (the megathrust) in a subduction zone. As they are steeper than the megathrust, the same amount of movement on them could result in more vertical displacement of the seafloor. Therefore, splay faults are thought to play an important role in the generation of tsunamis. Here, we use computer simulations to study if an earthquake can break multiple splay faults at once and how this affects the resulting tsunami. We find that multiple splay faults can indeed fail during a single earthquake due to the stress changes from trapped seismic waves, which promote rupture on splay faults. Rupture on splay faults results in larger seafloor displacements with smaller wavelengths, so the ensuing tsunami is bigger and results in two main flooding episodes at the coast. Our results show that it is important to consider rupture on splay faults when assessing tsunami hazard., Key Points: Multiple splay faults can be activated during a single earthquake by megathrust slip and dynamic stress transfer due to trapped waves. Splay fault activation is facilitated by their favorable orientation with respect to the local stress field and their closeness to failure. Long‐term geodynamic stresses and fault geometries affect dynamic splay fault rupture and the subsequent tsunami., Volkswagen Foundation (VolkswagenStiftung) http://dx.doi.org/10.13039/501100001663, Royal Society (The Royal Society) http://dx.doi.org/10.13039/501100000288, EC | H2020 | H2020 Priority Excellent Science | H2020 European Research Council (ERC) http://dx.doi.org/10.13039/100010663, Deutsche Forschungsgemeinschaft (DFG) http://dx.doi.org/10.13039/501100001659, National Science Foundation (NSF) http://dx.doi.org/10.13039/100000001, https://github.com/TUM-I5/SWE, https://doi.org/10.5281/zenodo.6969455

High-resolution computer simulations of earthquake sequences in three or even two dimensions pose great demands on time and energy, making lower-cost simplifications a competitive alternative. We systematically study the advantages and limitations of simplifications that eliminate spatial dimensions in quasi-dynamic earthquake sequence models, from 3D models with a 2D fault plane down to 0D or 1D models with a 0D fault point. We demonstrate that, when 2D or 3D models produce quasi-periodic characteristic earthquakes, their behavior is qualitatively similar to lower-dimension models. Certain coseismic characteristics like stress drop and fracture energy are largely controlled by frictional parameters and are thus largely comparable. However, other observations are quantitatively clearly affected by dimension reduction. We find corresponding increases in recurrence interval, coseismic slip, peak slip velocity, and rupture speed. These changes are to a large extent explained by the elimination of velocity-strengthening patches that transmit tectonic loading onto the velocity-weakening fault patch, thereby reducing the interseismic stress rate and enhancing the slip deficit. This explanation is supported by a concise theoretical framework, which explains some of these findings quantitatively and effectively estimates recurrence interval and slip. Through accounting for an equivalent stressing rate at the nucleation size h* into 2D and 3D models, 0D or 1D models can also effectively simulate these earthquake cycle parameters. Given the computational efficiency of lower-dimensional models that run more than a million times faster, this paper aims to provide qualitative and quantitative guidance on economical model design and interpretation of modeling studies., Journal of Geophysical Research: Solid Earth, 127 (8), ISSN:2169-9313, ISSN:0148-0227, ISSN:2169-9356

Physics-based dynamic rupture models capture the variability of earthquake slip in space and time and can account for the structural complexity inherent to subduction zones. Here we link tsunami generation, propagation, and coastal inundation with 3D earthquake dynamic rupture (DR) models initialized using a 2D seismo-thermo-mechanical geodynamic (SC) model simulating both subduction dynamics and seismic cycles. We analyze a total of 15 subduction-initialized 3D dynamic rupture-tsunami scenarios in which the tsunami source arises from the time-dependent co-seismic seafloor displacements with flat bathymetry and inundation on a linearly sloping beach. We first vary the location of the hypocenter to generate 12 distinct unilateral and bilateral propagating earthquake scenarios. Large-scale fault topography leads to localized up- or downdip propagating supershear rupture depending on hypocentral depth. Albeit dynamic earthquakes differ (rupture speed, peak slip-rate, fault slip, bimaterial effects), the effects of hypocentral depth (25–40 km) on tsunami dynamics are negligible. Lateral hypocenter variations lead to small effects such as delayed wave arrival of up to 100 s and differences in tsunami amplitude of up to 0.4 m at the coast. We next analyse inundation on a coastline with complex topo-bathymetry which increases tsunami wave amplitudes up to ≈1.5 m compared to a linearly sloping beach. Motivated by structural heterogeneity in subduction zones, we analyse a scenario with increased Poisson's ratio of ν = 0.3 which results in close to double the amount of shallow fault slip, ≈1.5 m higher vertical seafloor displacement, and a difference of up to ≈1.5 m in coastal tsunami amplitudes. Lastly, we model a dynamic rupture “tsunami earthquake” with low rupture velocity and low peak slip rates but twice as high tsunami potential energy. We triple fracture energy which again doubles the amount of shallow fault slip, but also causes a 2 m higher vertical seafloor uplift and the highest coastal tsunami amplitude (≈7.5 m) and inundation area compared to all other scenarios. Our mechanically consistent analysis for a generic megathrust setting can provide building blocks toward using physics-based dynamic rupture modeling in Probabilistic Tsunami Hazard Analysis.

The 1960 M9.5 Valdivia and 1964 M9.2 Alaska earthquakes caused a decimeters high secondary zone of uplift a few hundred kilometers landward of the trench. We analyze GPS data from the 2010 M8.8 Maule and 2011 M9.0 Tohoku-Oki earthquakes to reveal the persistent existence of a secondary zone of uplift due to great earthquakes at the megathrust interface. This uplift varies in magnitude and location, but consistently occurs at a few hundred kilometers landward from the trench and is likely mainly coseismic in nature. This secondary zone of uplift is systematically predicted by our 2D visco-elasto-plastic seismo-thermo-mechanical numerical simulations, which model both geodynamic and seismic cycle timescales. Through testing hypotheses in both simple and realistic setups, we propose that a superposition of two physical mechanisms could be responsible for this phenomenon. First, a wavelength is introduced through elastic buckling of a visco-elastically layered fore-arc that is horizontally compressed in the interseismic period. The consequent secondary zone of interseismic subsidence is elastically rebound during the earthquake into a secondary zone of relative uplift. Second, absolute and broader uplift is ensured through a mass conservation-driven return flow following accelerated slab penetration due to the megathrust earthquake. The dip and width of the seismogenic zone and resulting (deep) coseismic slip seem to have the largest affect on location and amplitude of the secondary zone of uplift. These results imply that subduction and mantle flow do not occur at constant rates, but are rather modulated by earthquakes. This suggests a link between deep mantle and shallow surface displacements even at time scales of minutes.

Natural fault networks are geometrically complex systems that evolve through time. The evolution of faults and their off-fault damage patterns are influenced by both dynamic earthquake ruptures and aseismic deformation in the interseismic period. To better understand each of their contributions to faulting we simulate both earthquake rupture dynamics and long-term deformation in a visco-elasto-plastic crust subjected to rate- and state-dependent friction. The continuum mechanics-based numerical model presented here includes three new features. First, a 2.5-D approximation is created to incorporate the effects of a viscoelastic lower crustal substrate below a finite depth. Second, we introduce a dynamically adaptive (slip-velocity-dependent) measure of fault width to ensure grid size convergence of fault angles for evolving faults. Third, fault localization is facilitated by plastic strain weakening of bulk rate and state friction parameters as inspired by laboratory experiments. This allows us to simulate sequences of episodic fault growth due to earthquakes and aseismic creep for the first time. Localized fault growth is simulated for four bulk rheologies ranging from persistent velocity weakening to velocity strengthening. Interestingly, in each of these bulk rheologies, faults predominantly localize and grow due to aseismic deformation. Yet, cyclic fault growth at more realistic growth rates is obtained for a bulk rheology that transitions from velocity-strengthening friction to velocity-weakening friction. Fault growth occurs under Riedel and conjugate angles and transitions towards wing cracks. Off-fault deformation, both distributed and localized, is typically formed during dynamic earthquake ruptures. Simulated off-fault deformation structures range from fan-shaped distributed deformation to localized splay faults. We observe that the fault-normal width of the outer damage zone saturates with increasing fault length due to the finite depth of the seismogenic zone. We also observe that dynamically and statically evolving stress fields from neighboring fault strands affect primary and secondary fault growth and thus that normal stress variations affect earthquake sequences. Finally, we find that the amount of off-fault deformation distinctly depends on the degree of optimality of a fault with respect to the prevailing but dynamically changing stress field. Typically, we simulate off-fault deformation on faults parallel to the loading direction. This produces a 6.5-fold higher off-fault energy dissipation than on an optimally oriented fault, which in turn has a 1.5-fold larger stress drop. The misalignment of the fault with respect to the static stress field thus facilitates off-fault deformation. These results imply that fault geometries bend, individual fault strands interact, and optimal orientations and off-fault deformation vary through space and time. With our work we establish the basis for simulations and analyses of complex evolving fault networks subject to both long-term and short-term dynamics.

Taking the full complexity of subduction zones into account is important for realistic modelling and hazard assessment of subduction zone seismicity and associated tsunamis. Studying seismicity requires numerical methods that span a large range of spatial and temporal scales. We present the first coupled framework that resolves subduction dynamics over millions of years and earthquake dynamics down to fractions of a second. Using a two-dimensional geodynamic seismic cycle (SC) method, we model 4 million years of subduction followed by cycles of spontaneous megathrust events. At the initiation of one such SC event, we export the self-consistent fault and surface geometry, fault stress and strength, and heterogeneous material properties to a dynamic rupture (DR) model. Coupling leads to spontaneous dynamic rupture nucleation, propagation and arrest with the same spatial characteristics as in the SC model. It also results in a similar material-dependent stress drop, although dynamic slip is significantly larger. The DR event shows a high degree of complexity, featuring various rupture styles and speeds, precursory phases, and fault reactivation. Compared to a coupled model with homogeneous material properties, accounting for realistic lithological contrasts doubles the amount of maximum slip, introduces local pulse-like rupture episodes, and relocates the peak slip from near the downdip limit of the seismogenic zone to the updip limit. When an SC splay fault is included in the DR model, the rupture prefers the splay over the shallow megathrust, although wave reflections do activate the megathrust afterwards.

Our physical understanding of earthquakes, along with our ability to forecast them, is hampered by limited indications on the current and future state of stress on faults. Integrating indirect observations, laboratory experiments and physics-based numerical modelling to quantitatively estimate this evolution is crucial. However, quantitative integrations are tenuous in light of the scarcity and uncertainty of observations and the difficulty of modelling the physics governing earthquakes. We show that observations and prior physical knowledge, along with their errors, can be efficiently integrated through the statistical framework of ensemble data assimilation (EDA), which is adopted from weather forecasting. To evaluate whether fault stress estimation and forecasting is possible, we perform a perfect model test in a subduction zone setup that mimicks a scaled laboratory experiment. Synthetic noised data on velocities and stresses from one point near the surface are assimilated using an Ensemble Kalman Filter. These data update the velocity and stress states throughout 150 ensemble members, whose dynamics is governed by a seismic cycle model. This visco-elasto-plastic forward model forecasts the system’s evolution through solving Navier–Stokes equations with a strongly rate-dependent friction coefficient. The ensemble assimilation of data from a single location provides probabilistic estimates of fault stress and dynamic strength evolution, which capture the true solution exceptionally well. This is possible, because the sampled error covariance matrix contains prior information from the physics that relates velocities, stresses and pressure at the surface to those at the fault. In the analysis step, this covariance allows stress and strength distributions to be reconstructed. In the subsequent forecast step the physical equations are solved to propagate the updated states forward in time. This provides probabilistic information on the likelihood of occurrence of the next earthquake in this synthetic laboratory setting. Throughout the ensemble simulations the forecasting ability for large, quasi-periodic events turns out to be significantly better than that of a periodic recurrence model. For example, it only requires an alarm to sound for 17 per cent instead of 68 per cent of the time to forecast 70 per cent of 21 events. We show that combining our prior knowledge of physical laws with observations through a Bayesian framework provides distinct added value with respect to using observations or numerical models independently. This educational test thus shows vast potential for including physics-based information into probabilistic seismic hazard assessment using EDA. To analyze its real world potential assumptions on an exact representation of the physics in a 2-D simplified system remain to be explored., Geophysical Journal International, 217 (3), ISSN:0956-540X, ISSN:1365-246X

Orientations of natural fault systems are subject to large variations. They often contradict classical Andersonian faulting theory as they are misoriented relative to the prevailing regional stress field. This is ascribed to local effects of structural or stress heterogeneities and reorientations of structures or stresses on the long-term. To better understand the relation between fault orientation and regional stresses, we simulate spontaneous fault growth and its effect on the stress field. Our approach incorporates earthquake rupture dynamics, visco-elasto-plastic brittle deformation and a rate-and state- dependent friction formulation in a continuum mechanics framework. We investigate how strike slip faults orient according to local and far-field stresses during their growth. We identify two modes of fault growth, seismic and aseismic, distinguished by different fault angles and slip velocities. Seismic fault growth causes a significant elevation of dynamic stresses and friction values ahead of the propagating fault tip. These elevated quantities result in a greater strike angle relative to the maximum principal regional stress than that of a fault segment formed aseismically. When compared to the near-tip time-dependent stress field the fault orientations produced by both growth modes follow Anderson’s classical faulting theory. We demonstrate how the two types of fault growth may be distinguished in natural faults by comparing their angles relative to the original regional max- imum principal stress. A stress field analysis of the Landers-Mojave fault suggests that an angle greater than approximately 25° between two faults indicates seismic fault growth.

There is increasing evidence that the Himalayan seismicity can be bimodal: blind earthquakes (up to Mw ~ 7.8) tend to cluster in the downdip part of the seismogenic zone, whereas infrequent great earthquakes (Mw 8+) propagate up to the Himalayan frontal thrust. To explore the causes of this bimodal seismicity, we developed a two-dimensional, seismic cycle model of the Nepal Himalaya. Our visco-elasto-plastic simulations reproduce important features of the earthquake cycle, including interseismic strain and a bimodal seismicity pattern. Bimodal seismicity emerges as a result of relatively higher friction and a non-planar geometry of the Main Himalayan Thrust fault. This introduces a region of large strength excess that can only be activated once enough stress is transferred upwards by blind earthquakes. This supports the view that most segments of the Himalaya might produce complete ruptures significantly larger than the 2015 Mw 7.8 Gorkha earthquake, which should be accounted for in future seismic hazard assessments., Nature Communications, 10, ISSN:2041-1723

We present a 2-D numerical modeling approach for simulating a wide slip spectrum in a viscoelastoplastic continuum. The key new model component is an invariant reformulation of the classical rate- and state-dependent friction equations, which is designed for earthquake simulations along spontaneously evolving faults. Here we describe the methodology and demonstrate that it is accurate and stable in a setup consisting of a mature strike-slip fault zone. We show that the nucleation and propagation of an earthquake are well resolved, as supported by a good agreement with various analytical approximations, including those of the nucleation and cohesive zone lengths. Results generally converge with respect to grid size, time step, and other numerical parameters. The convergence rate with respect to grid size depends on the internodal averaging scheme, is influenced by wave reflections, and deteriorates for inclined faults. The simulated slip spectrum, ranging from stable sliding at the loading rate to periodic aseismic slip to periodic seismic slip as a function of nucleation size, is in general agreement with the literature. In this simple setup, dynamic pressure does not play a significant role. By analyzing the role of viscous deformation, we identify and confirm by our simulations a theoretical viscosity threshold below which earthquakes cannot nucleate. This threshold is shown to depend on the reference strength of rate- and state-dependent friction and the loading strain rate, which is in agreement with previous work on the brittle-ductile transition.

The relative contribution of tectonic and kinematic processes to seismic behaviour of mountain belts is still controversial. To understand the partitioning between these processes we developed a model that simulates both tectonic and seismic processes in a continental collision setting. These 2D seismo-thermo-mechanical (STM) models obtain a Gutenberg–Richter frequency–magnitude distribution due to spontaneous events occurring throughout the orogen. Our simulations suggest that both the corresponding slope (b value) and maximum earthquake magnitude () correlate linearly with plate convergence rate. By analyzing 1D rheological profiles and isotherm depths we demonstrate that plate convergence rate controls the brittle strength through a rheological feedback with temperature and strain rate. Faster convergence leads to cooler temperatures and also results in more larger seismogenic domains, thereby increasing both and the relative number of large earthquakes (decreasing b value). This mechanism also predicts a more seismogenic lower crust, which is confirmed by a transition from uni- to bi-modal hypocentre depth distributions in our models. This transition and a linear relation between convergence rate and b value and is supported by our comparison of earthquakes recorded across the Alps, Apennines, Zagros and Himalaya. These results imply that deformation in the Alps occurs in a more ductile manner compared to the Himalayas, thereby reducing its seismic hazard. Furthermore, a second set of experiments with higher temperature and different orogenic architecture shows the same linear relation with convergence rate, suggesting that large-scale tectonic structure plays a subordinate role. We thus propose that plate convergence rate, which also controls the average differential stress of the orogen and its linear relation to the b value, is the first-order parameter controlling seismic hazard of mountain belts.

Correlations between geodynamic parameters and interplate seismicity characteristics in subduction zones are generally weak due to the short instrumental record and multiparameter influences. To investigate the role of subduction velocity Vs and the width of the seismogenic zone W on maximum magnitude Mmax, seismic rate τ, characteristic recurrence rate τc, and moment release rate MRR, we use synthetic data sets from simplified analog and numerical models to gain insight into natural subduction zones seismicity. Our models suggest that Mmax increases with W and is unaffected by Vs, τ increases with Vs, τc increases with Vs/W, and MRR increases with Vs × W. In nature, only the positive correlation between Vs and τ is significant. Random sampling of our time series suggest that the positive correlation between Vs and τ can be observed with short observation time windows. Other correlations, including Mmax versus W, become clear only for time window lengths longer than 1/τc.

""At oceanic margins, syn-convergent exhumation, subduction erosion, and inter-plate coupling are intimately related, but ample questions remain concerning their interaction and individual mechanisms. To analyze these interactions for a thick-skinned, visco-elastic wedge, we focus on properly modeling stresses, energies, and topographies at the inter-plate and wedge bounding interfaces using a Coulomb frictional contact algorithm. In this innovative plane-strain, free surface, Lagrangian finite element model, fault dynamics is modulated by retreating subduction. Subduction is dynamically driven by slab-pull due to a slab sinking in a semi-analytic, computationally favorable approximation of three-dimensional induced mantle flow. Nodal trajectories show that continuous underthrusting of a slab induces a steady state corner flow through forced underplating and subsequent trenchward extrusion due to gravitational spreading. This flow pattern confirms early-proposed models of syn-orogenic deep-seated rock exhumation propelled by coexisting extension and continuous shortening at depth. A distinct reduction in upward flowing material and accompanying decrease of exhumation velocities, to millimeters per year as observed in nature, is induced by a diversion of orogenic wedge material toward the mantle once a subduction channel is formed. The key parameter affecting model evolution and spontaneous formation of a subduction channel is basal friction, which modulates the amount of erosion. However, formation of a subduction channel entrance needs to be ensured through the deformability of the overriding plate, which is influenced by applied pressure at the overriding plate tip and material properties. The down dragging of the overriding plate is sufficient above a threshold inter-plate shear stress of about 2–7 MPa. ""