Computational earthquake sequence models provide generative estimates of the time, location, and size of synthetic seismic events that can be compared with observed earthquake histories and assessed as rupture forecasts. Here we describe a three-dimensional probabilistic earthquake sequence model that produces slip event time series constrained across geometrically complex non-planar fault systems. This model is kinematic in nature, integrating the time evolution of geometric moment accumulation and release with empirical earthquake scaling laws. The temporal probability of event occurrence is determined from the time history of geometric moment integrated with short-term Omori-style rate decay following each earthquake achieving long-term time-averaged moment balance. Similarly, the net geometric moment monotonically controls the probability of event localization, and seismic events release geometric moment with spatially heterogeneous slip on three-dimensional non-planar fault surfaces. We use this model to generate a synthetic earthquake sequence on the Nankai subduction zone over a 1,250-year-long interval, including 700+ M_W = 5.5-8.5 coseismic events, with decadal-to-centennial scale quiescent intervals quasi-periodic great earthquake clusters followed by aftershock sequences.

Non-inertial afterslip has been inferred to occur following large earthquakes. An explanation for this slow slip phenomenon is that coseismically generated stresses induce sliding on parts of a fault surface with velocity-strengthening frictional properties. Here we develop an alternative mesoscale heuristic explanation for afterslip based on the idea that afterslip may occur on any portion of a fault that exhibits positive residual geometric moment following an earthquake, including sections that ruptured coseismically. Following a large earthquake, this model exhibits exponential time decay of afterslip and allows for variable sensitivity to coseismic event magnitude and residual geometric moment. This model provides a partial explanation for the spatial relationship of co- and post-seismic slip associated with the 2011 M_W = 9 Tohoku-oki earthquake.

Forecasting the timing of earthquakes is a long-standing challenge. Moreover, it is still debated how to formulate this problem in a useful manner, or to compare the predictive power of different models. Here, we develop a versatile neural encoder of earthquake catalogs, and apply it to the fundamental problem of earthquake rate prediction, in the spatio-temporal point process framework. The epidemic type aftershock sequence model (ETAS) effectively learns a small number of parameters to constrain the assumed functional forms for the space and time correlations of earthquake sequences (e.g., Omori-Utsu law). Here we introduce learned spatial and temporal embeddings for point process earthquake forecasting models that capture complex correlation structures. We demonstrate the generality of this neural representation as compared with ETAS model using train-test data splits and how it enables the incorporation additional geophysical information. In rate prediction tasks, the generalized model shows >4% improvement in information gain per earthquake and the simultaneous learning of anisotropic spatial structures analogous to fault traces. The trained network can also be used to perform short-term prediction tasks, showing similar improvement while providing a 1000-fold reduction in run-time.

2022

Bem2d.jl: A quadratic co-location two-dimensional boundary element approach to quasi-static faulting problems with gravity and non-planar topography

Linear elastic boundary element models are commonly used tools to understand the mechanics of earthquake cycle processes and their contribution to the growth of tectonic structures. Here we describe a two-dimensional plane strain linear elastic boundary element approach to earthquake cycle and tectonic processes based on the displacement discontinuity method. This approach integrates an analytic solution for coincident interactions using three node quadratic elements and the classic particular integral approach to uni-directional gravity. Three node quadratic elements are more accurate than classical constant displacement elements and enable the exact representation displacements, stresses and tractions on elements subject to slip gradients. We demonstrate the recovery of analytic solutions and illustrate the combined effects of faulting and gravitational body forces in the presence of topographic relief.

Interseismic Masking of Fault Slip Deficit Rates by Earthquake Cycle Processes and Local Block Rotations

Interseismic geodetic observations at active boundaries contain information about tectonic motions and earthquake cycle processes. The combined effects of these two processes may serve to mask tectonic signals near active faults. Motivated by the observation of coseismic rotational motion during the 2016 M_W=7.8 Kaikōura New Zealand earthquake we develop idealized models for the masking of tectonic motions near blocks with local rotation poles. We show that differential interseismic velocities may reach local, rather than far field, maxima contributing to suppression of apparent slip deficit rates. This effect is most pronounced for small ( < 100 km) tectonic blocks with deep ( > 10 km) effective locking depths.

Kinematic Representations of Linear and Power-Law Viscoelastic Deformation Through the Earthquake Cycle

Meade, Brendan J.,
Mallick, Rishav,
and Carrero-Mustelier, Emily

Abstract Viscoelastic deformation below the Earth’s elastic crust is modulated by stresses generated by both plate tectonic and earthquake cycle processes. Rapid near-fault deformation following large earthquakes has been interpreted as the signature of viscoelastic stress diffusion in the upper mantle and earthquake cycle models have been developed that integrated this effect throughout the duration of the earthquake cycle. Here we develop a surrogate method to approximate the upper crustal kinematics associated with time-dependent viscoelastic deformation that does not require knowledge of nor provide constraints on subcrustal rheology. To do this we show how the effects of deformation in the viscoelastic upper mantle can be emulated by introducing a set of effective dislocations at crust-mantle interface. The approach is shown to approximate both linear and power-law viscoelastic rheologies and offers a way to accurately represent viscoelastic kinematics even in the case where upper mantle rheology is poorly constrained.

On the Choice and Implications of Rheologies That Maintain Kinematic and Dynamic Consistency Over the Entire Earthquake Cycle

Mallick, Rishav,
Lambert, Valère,
and Meade, Brendan

Journal of Geophysical Research: Solid Earth
Sep
2022

Viscoelastic processes in the upper mantle redistribute seismically generated stresses and modulate crustal deformation throughout the earthquake cycle. Geodetic observations of these motions at the surface of the crust-mantle system offer the possibility of constraining the rheology of the upper mantle. Parsimonious representations of viscoelastically modulated deformation through the aseismic phase of the earthquake cycle should simultaneously explain geodetic observations of (a) rapid postseismic deformation, (b) late in the earthquake cycle near-fault strain localization. To understand how rheological formulations affect kinematics, we compare predictions from time-dependent forward models of deformation over the entire earthquake cycle for an idealized vertical strike-slip fault in a homogeneous elastic crust underlain by a homogeneous viscoelastic upper-mantle. We explore three different rheologies as inferred from laboratory experiments: (a) linear Maxwell, (b) linear Burgers, (c) power-law. The linear Burgers and power-law rheologies are consistent with fast and slow deformation phenomenology over the entire earthquake cycle, while the single-layer linear Maxwell model is not. The kinematic similarity of linear Burgers and power-law models suggests that geodetic observations alone may be insufficient to distinguish between them, but indicate that one may serve as an effective proxy for the other. However, the power-law rheology model displays a postseismic response that is non-linearly dependent on earthquake magnitude, which may offer a partial explanation for observations of limited postseismic deformation near some magnitude 6.5–7.0 earthquakes. We discuss the role of mechanical coupling between frictional slip and viscous creep in controlling the time-dependence of regional stress transfer following large earthquakes and how this may affect the seismic hazard and risk to communities living close to fault networks.

2021

A global set of subduction zone earthquake scenarios and recurrence intervals inferred from geodetically constrained block models of interseismic couplingc distributions

The past 100 years have seen the occurrence of five M > 9 earthquakes and 94 M > 8 earthquakes. Here we assess the potential for future great earthquakes using inferences of interseismic subduction zone coupling from a global block model incorporating both tectonic plate motions and earthquake cycle effects. Interseismic earthquake cycle effects are represented using a first-order quasistatic elastic approximation and include 1e7 square kilometers of interacting fault system area across the globe. We use estimated spatial variations in decadal-duration coupling at 15 subduction zones and the Himalayan range front to estimate the locations and magnitudes of potential seismic events using empirical scaling relationships relating coupled area to moment magnitude. As threshold coupling values increase, estimates of potential earthquake magnitudes decrease, but the total number of large earthquakes varies non-monotonically. These rupture scenarios include as many as 14 recent or potential M > 9 earthquakes globally and up to 18 distinct M > 7 events associated with a single subduction zone (South America). We also combine estimated slip deficit rates and potential event magnitudes to calculate recurrence intervals for large earthquake scenarios, finding that almost all potential earthquakes would have a recurrence time of less than 1,000 years.

2019

Machine learning in seismology: Turning data into insights

Kong, Q.,
Trugman, D. T.,
Ross, Z. E.,
Bianco, M. J.,
Meade, B. J.,
and Gerstoft, P.

This article provides an overview of current applications of machine learning (ML) in seismology. ML techniques are becoming increasingly widespread in seismology, with applications ranging from identifying unseen signals and patterns to extracting features that might improve our physical understanding. The survey of the applications in seismology presented here serves as a catalyst for further use of ML. Five research areas in seismology are surveyed in which ML classification, regression, clustering algorithms show promise: earthquake detection and phase picking, earthquake early warning (EEW), ground-motion prediction, seismic tomography, and earthquake geodesy. We conclude by discussing the need for a hybrid approach combining data-driven ML with traditional physical modeling.

2018

Deep learning of aftershock patterns following large earthquakes

DeVries, P. M. R.,
Viégas, F.,
Wattenberg, M.,
and Meade, B. J.

Aftershocks are a response to changes in stress generated by large earthquakes and represent the most common observations of the triggering of earthquakes. The maximum magnitude of aftershocks and their temporal decay are well described by empirical laws (such as Bath’s law and Omori’s law), but explaining and forecasting the spatial distribution of aftershocks is more difficult. Coulomb failure stress change3 is perhaps the most widely used criterion to explain the spatial distributions of aftershocks, but its applicability has been disputed9. Here we use a deep-learning approach to identify a static-stress-based criterion that forecasts aftershock locations without prior assumptions about fault orientation. We show that a neural network trained on more than 131,000 mainshock-aftershock pairs can predict the locations of aftershocks in an independent test dataset of more than 30,000 mainshock-aftershock pairs more accurately (area under curve of 0.849) than can classic Coulomb failure stress change (area under curve of 0.583). We find that the learned aftershock pattern is physically interpretable: the maximum change in shear stress, the von Mises yield criterion (a scaled version of the second invariant of the deviatoric stress-change tensor) and the sum of the absolute values of the independent components of the stress-change tensor each explain more than 98 per cent of the variance in the neural-network prediction. This machine-learning-driven insight provides improved forecasts of aftershock locations and identifies physical quantities that may control earthquake triggering during the most active part of the seismic cycle.

The rotations of tectonic plates provide a partial description of the total observed displacements at the Earth’s surface. The estimated number of kinematically distinct plates has increased from 12 in 1990 to 56 in 2010 as a result of the increase in the number of kinematic observables. At length scales 1,000 km, rotation-only plate models are inaccurate because geodetic signals of long-term plate motions are complicated by earthquake cycle effects. Here we present results from a global block model that unifies large-scale plate motions and local earthquake cycle effects at plate boundaries. Incorporating the rotations of 307 distinct plates, elastic strain accumulation from 16 subduction zones and 1.59 x 1e7 square kilometers of fault system area, this model explains 19,664 interseismic GPS velocities at a resolution of 2.2 mm/year. Geodetically constrained fault slip deficit rates yield a cumulative global moment accumulation rate of 1.09 × 1022 N·m/year, 12% larger than the average annual coseismic moment release rate from 1900 to 2013. The potential contribution to the total moment rate budget can be estimated from the frequency distribution of the modeled fault slip-deficit rates, which follow an exponential distribution. Integrating this frequency distribution over all possible slip rates indicates that the geologic structures included in this reference global block model account for 98% of the global moment budget. Comparing our results with population distribution, we find that 50% of the world’s population lives within 200 km of an active fault with a slip rate >2 mm/year.

2017

A comparison of geodetic and geologic rates prior to large strike slip earthquakes: A diversity of earthquake cycle behaviors?

Comparison of preevent geodetic and geologic rates in three large-magnitude (Mw 7.6–7.9) strike-slip earthquakes reveals a wide range of behaviors. Specifically, geodetic rates of 26–28 mm/yr for the North Anatolian fault along the 1999 Mw = 7.6 Izmit rupture are 40% faster than Holocene geologic rates. In contrast, geodetic rates of 6–8 mm/yr along the Denali fault prior to the 2002 Mw 5 7.9 Denali earthquake are only approximately half as fast as the latest Pleistocene-Holocene geologic rate of 12 mm/yr. In the third example where a sufficiently long pre-earthquake geodetic time series exists, the geodetic and geologic rates along the 2001 Mw = 7.8 Kokoxili rupture on the Kunlun fault are approximately equal at 11 mm/yr. These results are not readily explicable with extant earthquake-cycle modeling, suggesting that they may instead be due to some combination of regional kinematic fault interactions, temporal variations in the strength of lithospheric-scale shear zones, and/or variations in local relative plate motion rate. Whatever the exact causes of these variable behaviors, these observations indicate that either the ratio of geodetic to geologic rates before an earthquake may not be diagnostic of the time to the next earthquake, as predicted by many rheologically based geodynamic models of earthquake-cycle behavior, or different behaviors characterize different fault systems in a manner that is not yet understood or predictable.

Block motion changes in Japan triggered by the 2011 great Tohoku earthquake

Plate motions are governed by equilibrium between basal and edge forces. Great earthquakes may induce differential static stress changes across tectonic plates, enabling a new equilibrium state. Here we consider the torque balance for idealized circular plates and find a simple scalar relationship for changes in relative plate speed as a function of its size, upper mantle viscosity, and coseismic stress changes. Applied to Japan, the 2011 Mw 9.0 Tohoku earthquake generated coseismic stresses of 102–105 Pa that could have induced changes in motion of small (radius \100 km) crustal blocks within Honshu. Analysis of timedependent GPS velocities, with corrections for earthquake cycle effects, reveals that plate speeds may have changed by up to \3 mm/yr between $3.75 year epochs bracketing this earthquake, consistent with an upper mantle viscosity of 3e18 Pa s, suggesting that great earthquakes may modulate motions of proximal crustal blocks at frequencies as high as 1028 Hz.

Enabling large‐scale viscoelastic calculations via neural network acceleration

DeVries, P. M. R.,
Thompson, T. B.,
and Meade, B. J.

One of the most significant challenges involved in efforts to understand the effects of repeated earthquake cycle activity is the computational costs of large-scale viscoelastic earthquake cycle models. Computationally intensive viscoelastic codes must be evaluated at thousands of times and locations, and as a result, studies tend to adopt a few fixed rheological structures and model geometries and examine the predicted time-dependent deformation over short (<10 years) time periods at a given depth after a large earthquake. Training a deep neural network to learn a computationally efficient representation of viscoelastic solutions, at any time, location, and for a large range of rheological structures, allows these calculations to be done quickly and reliably, with high spatial and temporal resolutions. We demonstrate that this machine learning approach accelerates viscoelastic calculations by more than 50,000%. This magnitude of acceleration will enable the modeling of geometrically complex faults over thousands of earthquake cycles across wider ranges of model parameters and at larger spatial and temporal scales than have been previously possible.

Interactive visualization of spatially amplified GNSS time-series position fields

Meade, B. J.,
Freeman, W. T.,
Wilson, J.,
Viégas, F.,
and Wattenberg, M.

Global Navigation Satellite System (GNSS) position time series are used pervasively in earthquake science to measure the surface response to earthquake cycle deformation. Characteristic usage cases are focused on the temporal windowing of position data to isolate coseismic, postseismic, or interseismic deformation. Here, we present an interactive visualization approach for the temporal evolution of GNSS time session in 2D in which the position estimates are amplified relative to their true positions, or are amplified relative to a reference state. This approach enables a rapid visual assessment of deformation patterns across all phases of the earthquake cycle in relation to topographic structure and active faults including azimuth reversal of coseismic and interseismic deformation.

Intermittent granular dynamics at a seismogenic plate boundary

Earthquakes at seismogenic plate boundaries are a response to the differential motions of tectonic blocks embedded within a geometrically complex network of branching and coalescing faults. Elastic strain is accumulated at a slow strain rate on the order of 1e-15 s, and released intermittently at intervals >100 yr, in the form of rapid (seconds to minutes) coseismic ruptures. The development of macroscopic models of quasistatic planar tectonic dynamics at these plate boundaries has remained challenging due to uncertainty with regard to the spatial and kinematic complexity of fault system behaviors. The characteristic length scale of kinematically distinct tectonic structures is particularly poorly constrained. Here, we analyze fluctuations in Global Positioning System observations of interseismic motion from the southern California plate boundary, identifying heavy-tailed scaling behavior. Namely, we show that, consistent with findings for slowly sheared granular media, the distribution of velocity fluctuations deviates from a Gaussian, exhibiting broad tails, and the correlation function decays as a stretched exponential. This suggests that the plate boundary can be understood as a densely packed granular medium, predicting a characteristic tectonic length scale of 91 ± 20 km, here representing the characteristic size of tectonic blocks in the southern California fault network, and relating the characteristic duration and recurrence interval of earthquakes, with the observed sheared strain rate, and the nanosecond value for the crack tip evolution time scale. Within a granular description, fault and blocks systems may rapidly rearrange the distribution of forces within them, driving a mixture of transient and intermittent fault slip behaviors over tectonic time scales.

Viscoelastic block models of the North Anatolian fault: A unified earthquake cycle representation of pre- and post-seismic geodetic observations

DeVries, P. M. R.,
Krastev, P. G.,
Dolan, J. F.,
and Meade, B. J.

Bulletin of the Seismological Society of America
2017

Along the North Anatolian fault (NAF), the surface deformation associated with tectonic block motions, elastic strain accumulation, and the viscoelastic response to past earthquakes has been geodetically observed over the last two decades. These observations include campaign-mode Global Positioning System (GPS) velocities from the decade prior to the 1999 Mw 7.4 İzmit earthquake and seven years of continuously recorded postseismic deformation following the seismic event. Here, we develop a 3D viscoelastic block model of the greater NAF region, including the last 2000 yrs of earthquake history across Anatolia, to simultaneously explain geodetic observations from both before and after the İzmit earthquake. With a phenomenologically motivated simple two-layer structure (schizosphere and plastosphere) and a Burgers rheology (with Maxwell viscosity log10 ηM ≈ 18:6–19:0 Pa·s and Kelvin viscosity log10 ηK ≈ 18.0–19.0 Pa·s), a block model that incorporates tectonic plate motions, interseismic elastic strain accumulation, transient viscoelastic perturbations, and internal strain can explain both the pre- and post-İzmit earthquake observations with a single unified model. Viscoelastic corrections to the interseismic GPS velocity field with the unified model reach magnitudes of 2-9 mm/yr. Geodetically constrained slip-deficit rate estimates along the central NAF and northern strand of the NAF in the Sea of Marmara vary nonmonotonically with Maxwell viscosity and change by up to 23% ( 4 mm=yr) for viscosities ranging from 1018 to 1023 Pa·s. For the best-fit viscosity structures, central NAF slip-deficit rates reach 22 mm=yr, increasing to 28 mm=yr in the Sea of Marmara. Along the central NAF, these rates are similar to the fastest geologic slip-rate estimates. The fastest slip-deficit rate estimates along the entire fault system ( 27–28 mm=yr) occur less than 50 km from Istanbul, along the northern strand of the NAF in the Sea of Marmara.

What is better than Coulomb failure stress? A ranking of scalar static stress triggering mechanisms from 10e5 mainshock-aftershock pairs

Meade, B. J.,
DeVries, P. M. R.,
Faller, J.,
Viégas, F.,
and Wattenberg, M.

Aftershocks may be triggered by the stresses generated by preceding mainshocks. The temporal frequency and maximum size of aftershocks are well described by the empirical Omori and Bath laws, but spatial patterns are more difficult to forecast. Coulomb failure stress is perhaps the most common criterion invoked to explain spatial distributions of aftershocks. Here we consider the spatial relationship between patterns of aftershocks and a comprehensive list of 38 static elastic scalar metrics of stress (including stress tensor invariants, maximum shear stress, and Coulomb failure stress) from 213 coseismic slip distributions worldwide. The rates of true-positive and false-positive classification of regions with and without aftershocks are assessed with receiver operating characteristic analysis. We infer that the stress metrics that are most consistent with observed aftershock locations are maximum shear stress and the magnitude of the second and third invariants of the stress tensor. These metrics are significantly better than random assignment at a significance level of 0.005 in over 80% of the slip distributions. In contrast, the widely used Coulomb failure stress criterion is distinguishable from random assignment in only 51–64% of the slip distributions. These results suggest that a number of alternative scalar metrics are better predictors of aftershock locations than classic Coulomb failure stress change.

2016

Geodetically constrained models of viscoelastic stress transfer and earthquake triggering along the North Anatolian fault

DeVries, P. M. R.,
Krastev, P. G.,
and Meade, B. J.

Over the past 80 years, 8 Mw > 6.7 strike-slip earthquakes west of 408 longitude have ruptured the North Anatolian fault (NAF) from east to west. The series began with the 1939 Erzincan earthquake in eastern Turkey, and the most recent 1999 MW 5 7.4 Izmit earthquake extended the pattern of ruptures into the Sea of Marmara in western Turkey. The mean time between seismic events in this westward progression is 8.5 6 11 years (67% confidence interval), much greater than the timescale of seismic wave propagation (seconds to minutes). The delayed triggering of these earthquakes may be explained by the propagation of earthquake-generated diffusive viscoelastic fronts within the upper mantle that slowly increase the Coulomb failure stress change (DCFS) at adjacent hypocenters. Here we develop three-dimensional stress transfer models with an elastic upper crust coupled to a viscoelastic Burgers rheology mantle. Both the Maxwell (1e19 Pa s) and Kelvin (1e19 Pa s) viscosities are constrained by studies of geodetic observations before and after the 1999 Izmit earthquake. We combine this geodetically constrained rheological model with the observed sequence of large earthquakes since 1939 to calculate the time evolution of DCFS changes along the North Anatolian fault due to viscoelastic stress transfer. Apparent threshold values of mean DCFS at which the earthquakes in the eight decade sequence occur are between 0.02 to 3.15 MPa and may exceed the magnitude of static DCFS values by as much as 177%. By 2023, we infer that the mean time-dependent stress change along the northern NAF strand in the Marmara Sea near Istanbul, which may have previously ruptured in 1766, may reach the mean apparent time-dependent stress thresholds of the previous NAF earthquakes.

Kinematically consistent models of viscoelastic stress evolution

Following large earthquakes, coseismic stresses at the base of the seismogenic zone may induce rapid viscoelastic deformation in the lower crust and upper mantle. As stresses diffuse away from the primary slip surface in these lower layers, the magnitudes of stress at distant locations (>1 fault length away) may slowly increase. This stress relaxation process has been used to explain delayed earthquake triggering sequences like the 1992 Mw = 7.3 Landers and 1999 Mw = 7.1 Hector Mine earthquakes in California. However, a conceptual difficulty associated with these models is that the magnitudes of stresses asymptote to constant values over long time scales. This effect introduces persistent perturbations to the total stress field over many earthquake cycles. Here we present a kinematically consistent viscoelastic stress transfer model where the total perturbation to the stress field at the end of the earthquake cycle is zero everywhere. With kinematically consistent models, hypotheses about the potential likelihood of viscoelastically triggered earthquakes may be based on the timing of stress maxima, rather than on any arbitrary or empirically constrained stress thresholds. Based on these models, we infer that earthquakes triggered by viscoelastic earthquake cycle effects may be most likely to occur during the first 50% of the earthquake cycle regardless of the assumed long-term and transient viscosities.

Two decades of spatiotemporal variations in subduction zone coupling offshore Japan

Spatial patterns of interplate coupling on global subduction zones can be used to guide seismic hazard assessment, but estimates of coupling are often constrained using a limited temporal range of geodetic data. Here we analyze 19 years of geodetic observations from the GEONET network to assess time-dependent variations in the spatial distribution of coupling on the subduction zones offshore Japan. We divide the position time series into five, 3.75-year epochs each decomposed into best-fit velocity, annual periodic signals, coseismic offsets, and postseismic effects following seven major earthquakes. Nominally interseismic velocities are interpreted in terms of a combination of tectonic block motions and earthquake cycle activity. The duration of the inferred postseismic activity covaries with the linear velocity. To address this trade-off, we assume that the nominally interseismic velocity at each station varies minimally from epoch to epoch. This approach is distinct from prior time-series analysis across the earthquake cycle in that position data are not detrended using preseismic velocity, which inherently assumes that interseismic processes are spatially stable through time, but rather the best-fit velocity at each station may vary between epochs. These velocities reveal significant consistency since 1996 in the spatial distribution of coupling on the Nankai subduction zone, with variation limited primarily to the Tokai and Bungo Channel regions, where long-term slow slip events have occurred, and persistently coupled regions coincident with areas that slipped during historic great earthquakes. On the Sagami subduction zone south of Tokyo, we also estimate relatively stable coupling through time. On the Japan–Kuril Trench, we image significant coupling variations owing to effects of the 1994 Mw = 7.7 Sanriku-oki, 2003 Mw = 8.2 Tokachi-oki, and 2011 Mw = 9.0 Tohoku-oki earthquakes. In particular, strong coupling becomes more spatially extensive following the 1994 event until 2011, coseismic-sense slip precedes the Tohoku-oki event, and coupling offshore northern Honshu is reduced after the 2011 earthquake. Despite the occurrence of the 2003 Tokachi-oki earthquake, persistent coupling offshore Hokkaido suggests ongoing seismic hazard, possibly similar to past Mw 9-class earthquakes interpreted from coastal paleoseismic records. This time-dependent analysis of interseismic deformation illuminates rich diversity in the distribution of subduction zone coupling, including spatiotemporal stability in coupling, effective reduction in strongly coupled regions due to aseismic thrust-sense slip events, and broad changes in the distribution of coupling following major earthquakes.

2015

Kinematic barrier constraints on the magnitudes of additional great earthquakes off the east coast of Japan

In the wake of the 2011 great Tohoku-Oki, Japan, earthquake, questions were asked about the extent to which this earthquake ruptured a part of the subduction zone that had been strongly coupled in the decades prior to the event (Avouac, 2011; Lay and Kanamori, 2011; Loveless and Meade, 2011). The primary means for investigating this question is the analysis of interseismic geodetic data. These data record both tectonic motions and the accumulation of elastic strain the upper crust. Combined with elastic dislocation and block models, Global Navigation Satellite System (GNSS) data from the GNSS Earth Observation Network System (GEONET) array had been used to estimate the spatial variation in interseismic coupling (defined as the fraction of relative motion between tectonic plates that contributes to elastic strain accumulation) along the Japan trench subduction interface prior to the occurrence of the 2011 event (Nishimura et al., 2004; Suwa et al., 2006; Hashimoto et al., 2009; Loveless and Meade, 2010). Each of these estimated interseismic coupling distributions identified a partially to strongly coupled region, extending for 400 km between 36° and 40° north latitude along the Japan trench, approximately bounding the lateral and down-dip extents of the coseismic rupture (Loveless and Meade, 2011). However, none of these estimates show a perfect spatial correspondence between the coseismic rupture extent and complete preseismic coupling. Rather, the coseismic rupture was laterally bounded by regions of the Japan subduction zone with a coupling fraction c of 0.3, in which 30% of the relative plate motion contributes to elastic strain accumulation during the interseismic regime (Loveless and Meade, 2011). This spatial relationship could result from insufficient elastic strain accumulation around portions of the interface characterized by lower, but nonzero, coupling values, which act as kinematic barriers to rupture propagation. In this article, we develop candidate rupture scenarios across the Japan, Sagami, and Nankai subduction zones, based on spatial variations in interseismic coupling and empirical earthquake parameter scaling laws. Specifically, we examine spatial clusters of the subduction interfaces within which the coupling exceeds a given fraction of plate convergence rate and estimate the moment magnitude of an earthquake that could occur within that region, based on a scaling relationship between cluster area and earthquake magnitude. We consider cases in which the shallowest portions of the plate interfaces, near the subduction trenches, are assumed to creep interseismically and in which they are allowed to be coupled. The assumed kinematic behavior affects potential rupture areas by up to 40%, though this difference may be reduced if coseismic dynamic overshoot (which permits ruptures to propagate into interseismically creeping regions) is prevalent near the trenches.

Rapid slip-deficit rates at the eastern margin of the Tibetan plateau prior to the 2008 Mw 7.9 Wenchuan earthquake

Thompson, T. B.,
Plesch, A.,
Shaw, J. H.,
and Meade, B. J.

The Longmen Shan is the steepest topographic front at the India-Asia collision zone and the site of the Mw 7.9 Wenchuan earthquake. Here to explain the interseismic GPS velocities across the greater Longmen Shan region, we develop a boundary element model including earthquake cycle effects, topography, the westward dipping Beichuan Fault and a 20 km deep, shallowly dipping, detachment, inferred from observed afterslip and from structural considerations. Previous analyses which neglected the detachment and earthquake cycle effects have found shortening rates near zero. In contrast, we find that interseismic GPS data are consistent with a shortening rate of 5.7 \pm 1.5 mm/yr and maximum surface slip-deficit rate of 9.5 \pm 2.5 mm/yr. This model unifies the interpretation of geodetic deformation throughout the earthquake cycle and suggests that the Longmen Shan is an active fold-and-thrust belt with of Wenchuan-like recurrence intervals as short as 600 years.

Total variation regularization of geodetically and geologically constrained block models for the Western United States

Geodetic observations of interseismic deformation in the Western United States provide constraints on microplate rotations, earthquake cycle processes, and slip partitioning across the Pacific–North America Plate boundary. These measurements may be interpreted using block models, in which the upper crust is divided into microplates bounded by faults that accumulate strain in a first-order approximation of earthquake cycle processes. The number and geometry of microplates are typically defined with boundaries representing a limited subset of the large number of potentially seismogenic faults. An alternative approach is to include a large number of potentially active faults bounding a dense array of microplates, and then algorithmically estimate the boundaries at which strain is localized. This approach is possible through the application of a total variation regularization (TVR) optimization algorithm, which simultaneously minimizes the L2 norm of data residuals and the L1 norm of the variation in the differential block motions. Applied to 3-D spherical block models, the TVR algorithm can be used to reduce the total variation between estimated rotation vectors, effectively grouping microplates that rotate together as larger blocks, and localizing fault slip on the boundaries of these larger block clusters. Here we develop a block model comprised of 137 microplates derived from published fault maps, and apply the TVR algorithm to identify the kinematically most important faults in the western United States. This approach reveals that of the 137 microplates considered, only 30 unique blocks are required to approximate deformation in the western United States at a residual level of 2 mm/yr.

2013

Earthquake cycle deformation in the Tibetan plateau with a weak mid-crustal layer

Geodetic observations of interseismic deformation across the Tibetan plateau contain information about both tectonic and earthquake cycle processes. Time-variations in surface velocities between large earthquakes are sensitive to the rheological structure of the subseismogenic crust, and, in particular, the viscosity of the middle and lower crust. Here we develop a semianalytic solution for time-dependent interseismic velocities resulting from viscoelastic stress relaxation in a localized midcrustal layer in response to forcing by a sequence of periodic earthquakes. Earthquake cycle models with a weak midcrustal layer exhibit substantially more near-fault preseismic strain localization than do classic two-layer models at short (<100 yr) Maxwell times. We apply both this three-layer model and the classic two-layer model to geodetic observations before and after the 1997 Mw = 7.6 Manyi and 2001 Mw = 7.8 Kokoxili strike-slip earthquakes in Tibet to estimate the viscosity of the crust below a 20 km thick seismogenic layer. For these events, interseismic stress relaxation in a weak (viscosity ≤ 1e18.5 Pa s) and thin (height ≤20 km) midcrustal layer explains observations of both preseismic near-fault strain localization and rapid (>50 mm/yr) postseismic velocities in the years following the coseismic ruptures. We suggest that earthquake cycle models with a localized midcrustal layer can simultaneously explain both preseismic and postseismic geodetic observations with a single Maxwell viscosity, while the classic two-layer model requires a rheology with multiple relaxation time scales.

Edge-driven mechanical microplate models of strike-slip faulting in the Tibetan plateau

The India-Asia collision zone accommodates the relative motion between India and Eurasia through both shortening and pervasive strike-slip faulting. To gain a mechanical understanding of how fault slip rates are driven across the Tibetan plateau, we develop a two-dimensional, linear elastic, two-stage, deformable microplate model for the upper crust based on the behavior of an idealized earthquake cycle. We use this approach to develop a suite of simple India-Asia collision zone models, differing only in boundary conditions, to determine which combination of edge forces and displacements are consistent with both the slip rate measurements along major Tibetan faults as well as the geodetically observed extrusion of crustal material toward Southeast Asia. Model predictions for the Altyn Tagh (1-14 mm/yr), Kunlun (3-10 mm/yr), Karakorum (5-12 mm/yr), and Haiyuan (3-5 mm/yr) faults are in agreement with geologically and geodetically inferred slip rates. Further, models that accurately reproduce observed slip rate gradients along the Altyn Tagh and Kunlun faults feature two critical boundary conditions: (1) oblique compressive displacement along the Himalayan range front west of the Shillong plateau, and (2) forcing in Southeast Asia. Additionally, the ratio of internal-block potency rate to the total potency rate for each microplate ranges from 28% to 79%, suggesting a hybrid view of deformation in Tibet as simultaneously localized on major faults and distributed at length scales <500 km.

Inference of multiple earthquake cycle relaxation timescales from irregular geodetic sampling of interseismic deformation

Meade, B. J.,
Klinger, Y.,
and Hetland, E. A.

Bulletin of the Seismological Society of America
2013

Characterizing surface deformation throughout a full earthquake cycle is a challenge due to the lack of high-resolution geodetic observations of duration comparable to that of characteristic earthquake recurrence intervals (250–10,000 years). Here we approach this problem by comparing long-term geologic slip rates with geodetically derived fault slip rates by sampling only a short fraction (0.001%–0.1%) of a complete earthquake cycle along 15 continental strike-slip faults. Geodetic observations provide snapshots of surface deformation from different times through the earthquake cycle. The timing of the last earthquake on many of these faults is poorly known, and may vary greatly from fault to fault. Assuming that the underlying mechanics of the seismic cycle are similar for all faults, geodetic observations from different faults may be interpreted as samples over a significantly larger fraction of the earthquake cycle than could be obtained from the geodetic record along any one fault alone. As an ensemble, we find that geologically and geodetically inferred slip rates agree well with a linear relation of 0.94 +/- 0.09. To simultaneously explain both the ensemble agreement between geologic and geodetic slip-rate estimates with observations of rapid postseismic deformation, we consider the predictions from simple two-layer earthquake-cycle models with both Maxwell and Burgers viscoelastic rheologies. We find that a two-layer Burgers model, with two relaxation timescales, is consistent with observations of deformation throughout the earthquake cycle, whereas the widely used two-layer Maxwell model with a single relaxation timescale, is not, suggesting that the earthquake cycle is effectively characterized by a largely stressrecoverable rapid postseismic stage and a much more slowly varying interseismic stage.

Revisiting the orogenic energy balance in the western Taiwan orogen with weak faults

Orogens at convergent margins must meet the energetic requirements necessary to lift rocks against gravity, allow for frictional sliding along basal detachments and accommodate internal deformation processes. The combination of critical taper and kinematic wedge theories predicts the partitioning between these energy sinks as a function of both fault and crustal strengths. Integrating contemporary estimates of both crustal pore fluid pressures and the coefficient of friction on major faults, we find that work associated with internal deformation processes is the dominant energy sink in the western Taiwan orogenic wedge. These processes consume 54% of the total work budget, while the dissipation rates associated with frictional sliding on the basal detachment and lifting rocks against gravity are lower, requiring only 11% and 35% respectively. These results suggest a mechanical dichotomy in orogenic wedges where the faulting dominates the kinematic deformation budget, but internal distributed deformation processes dominate the energy budget.

2012

Geodetic constraints on San Francisco Bay Area fault slip rates and potential seismogenic asperities on the partially creeping Hayward fault

The Hayward fault in the San Francisco Bay Area (SFBA) is sometimes considered unusual among continental faults for exhibiting significant aseismic creep during the interseismic phase of the seismic cycle while also generating sufficient elastic strain to produce major earthquakes. Imaging the spatial variation in interseismic fault creep on the Hayward fault is complicated because of the interseismic strain accumulation associated with nearby faults in the SFBA, where the relative motion between the Pacific plate and the Sierra block is partitioned across closely spaced subparallel faults. To estimate spatially variable creep on the Hayward fault, we interpret geodetic observations with a three-dimensional kinematically consistent block model of the SFBA fault system. Resolution tests reveal that creep rate variations with a length scale of <15 km are poorly resolved below 7 km depth. In addition, creep at depth may be sensitive to assumptions about the kinematic consistency of fault slip rate models. Differential microplate motions result in a slip rate of 6.7 ± 0.8 mm/yr on the Hayward fault, and we image along-strike variations in slip deficit rate at 15 km length scales shallower than 7 km depth. Similar to previous studies, we identify a strongly coupled asperity with a slip deficit rate of up to 4 mm/yr on the central Hayward fault that is spatially correlated with the mapped surface trace of the 1868 Mw = 6.9-7.0 Hayward earthquake and adjacent to gabbroic fault surfaces.

Geodetic imaging of coseismic slip and postseismic afterslip: Sparsity promoting methods applied to the great Tohoku earthquake

Geodetic observations of surface displacements during and following earthquakes such as the March 11, 2011 great Tohoku earthquake can be used to constrain the spatial extent of coseismic slip and postseismic afterslip, and characterize the spectrum of earthquake cycle behaviors. Slip models are often regularized by assuming that slip on the fault varies smoothly in space, which may result in the artificial smearing of fault slip beyond physical boundaries. Alternatively, it may be desirable to estimate a slip distribution that is spatially compact and varies sharply. Here we show that sparsity promoting state vector regularization methods can be used to recover slip distributions with sharp boundaries, representing an alternative end-member result to very smooth slip distributions. Using onshore GPS observations at 298 stations during and in the 2 weeks following the Tohoku earthquake, we estimate a band of coseismic slip between 30 and 50 km depth extending 500 km along strike with a maximum slip of 64 m, corresponding to a minimum magnitude estimate of Mw = 8.8. Our estimate of afterslip is located almost exclusively down-dip of the coseismic rupture, with a transition between 40 and 50 km depth and an equivalent moment magnitude Mw = 8.2. This depth may be interpreted as coincident with the transition from velocity strengthening to velocity weakening frictional behavior, consistent with the upper limit of cold subduction estimates of the thermal structure of the Japan trench.

2011

Partitioning of localized and diffuse deformation in the Tibetan plateau from joint inversions of geologic and geodetic observations

The spatial complexity of continental deformation in the greater Tibetan plateau region can be defined as the extent to which relative motion of the Indian and Asian plates is partitioned between localized slip on major faults and distributed deformation. Potency rates provide a quantitative metric for determining the magnitudes of diffuse and on-fault crustal deformation, which are proportional to strain rates within crustal micro-plates and fault slip rates, respectively. We simultaneously estimate micro-plate rotation rates, interseismic elastic strain accumulation, fault slip rates on major structures, and strain rates within 24 tectonic micro-plates inferred from active fault maps in the greater Tibetan plateau region using quasi-static block models constrained by interseismic surface velocities at 731 GPS sites and 9 Holocene–Late Quaternary geologic fault slip rates. The joint geodetic-geologic inversion indicates that geologic slip rates are kinematically consistent with differential micro-plate motions. Estimated left-lateral slip rates on the Altyn Tagh, west-central Kunlun, and Xianshuihe faults are relatively homogeneous along strike ( 11.5, 10.5, and 11 mm/yr, respectively) while segmentation of the eastern Kunlun fault by the intersecting Elashan and Riyueshan faults results in a decreased slip rate, consistent with geologic observations. The fraction, ϕ, of total potency rate associated with intrablock strain, uncorrected for observational noise, ranges from 0.27 in the Himalayan Range block to 0.91 in the Aksai Chin block. Monte Carlo simulations are used to quantify the likelihood that internal deformation is statistically distinguishable from the uncertainties in geodetic velocities. These simulations indicate that internal block deformation is statistically significant only within the Himalayan Range Front (where internal deformation accounts for ϕ = 0.11 of block potency rate budget), west-central plateau (ϕ = 0.68), Ganzi-Yushu/Xianshuihe (0.58), Burma (0.09), Jiali (0.38), and Aksai Chin (0.68) blocks. In the other 18 tectonic micro-plates within the plateau region, estimated internal block potency is not currently distinguishable from the expected contribution of observational noise to residual velocities. Of the total potency budget within the Tibetan plateau, 85% is taken up by slip on major faults, with the remaining 15% accommodated by internal processes at sub-block scale distinguishable from observational noise. The localization of the majority of plate boundary activity is also supported by the spatial distribution of modern and historical crustal earthquakes. Sixty-seven percent of the total moment released by earthquakes in the CMT catalog and 90% of historical moment since 1900 have been released within 25 km of the major faults included in the block model, representing only 12% of the characteristic half-block length scale of 215 km. The localization of deformation inferred from geologic, geodetic, and seismic observations suggests that forces applied to tectonic micro-plates drive fault system activity at the India–Asia collision zone over decadal to Quaternary time scales.

Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 Mw = 9.0 Tohoku-oki earthquake

Imaging the extent to which the rupture areas of great earthquakes coincide with regions of pre-seismic interplate coupling is central to understanding patterns of strain accumulation and release through the earthquake cycle. Both geodetic and seismic estimates of the coseismic rupture extent for the March 11, 2011 Mw = 8.9-9.0 earthquake Tohoku-oki earthquake may be spatially correlated (0.26 \pm 0.05 to 0.82 \pm 0.05) with a region estimated to be partially to fully coupled in the interseismic period preceding the earthquake, though there is substantial variation in the estimated distribution and magnitude of coseismic slip. The 400 km-long region estimated to have slipped ≥4 m corresponds to an area of the subduction zone interface that was coupled at ≥30% of long-term plate convergence rate, with peak slip near a region coupled ≥80%. The northern termination of rupture is collocated with a region of relatively low (<20%) interseismic coupling near the epicenter of the 1994 Mw = 7.6 Sanriku-oki earthquake, and near a region of potential long-term low coupling or ongoing slow slip. Slip on the subduction interface beneath the coastline (40–50 km depth) is best constrained by the land-based GPS data and least constrained on the shallowest portion of the plate interface due to the 230 km distance between geodetic observations and the Japan trench.

Stress modulation on the San Andreas fault by interseismic fault system interactions

During the interseismic phase of the earthquake cycle, between large earthquakes, stress on faults evolves in response to elastic strain accumulation driven by tectonic plate motions. Because earthquake cycle processes induce non-local stress changes, the interseismic stress accumulation rate on one fault is influenced by the behavior of all nearby faults. Using a geodetically constrained block model, we show that the total interseismic elastic strain field generated by fault interactions within Southern California may increase stressing rates on the Mojave and San Bernardino sections of the San Andreas fault within the Big Bend region by as much as 38% relative to estimates from isolated San Andreas models. Assuming steady fault system behavior since the C.E. 1857 Fort Tejon earthquake, shear stress accumulated on these sections due only to interaction with faults other than the San Andreas reaches 1 MPa, 3 times larger than the coseismic and postseismic stress changes induced by recent Southern California earthquakes. Stress increases along Big Bend sections coincide with the greatest earthquake frequency inferred from a 1500-yr-long paleoseismic record and may affect earthquake recurrence intervals within geometrically complex fault systems, including the sections of the San Andreas fault closest to metropolitan Los Angeles.

2010

Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan

Interseismic deformation in Japan results from the combined effects of tectonic processes including rotation of crustal blocks and the earthquake cycle process of elastic strain accumulation about upper plate faults and subduction zone interfaces. We use spherical linear block theory constrained by geodetic observations from densely spaced Global Positioning System (GPS) stations to estimate plate motions, fault slip rates, and spatially variable interplate coupling on the Japan-Kuril, Sagami, and Nankai subduction zones. The reference model developed in this paper consists of 20 blocks, produces a mean residual velocity magnitude of 1.84 mm/yr at 950 stations, and accounts for 96% of the observed interseismic deformation signal. We estimate fault slip rates in excess of 15 mm/yr along the Niigata-Kobe Tectonic Zone and Itoigawa-Shizuoka Tectonic Line through central Japan, confirming their hypothesized roles as major tectonic boundaries. Oblique convergence across the Nankai Trough is partitioned, with 3/4 of the 30 mm/yr of trench-parallel motion accommodated by strike-slip motion on the subduction interface and the remaining 1/4 accommodated by right-lateral slip on the Median Tectonic Line. In contrast, our models suggest negligible slip partitioning in eastern Hokkaido, where oblique slip on the Japan-Kuril subduction interface accommodates all of the trench-parallel component of relative plate motion. Inferred spatial variations in the rake and magnitude of slip deficit on subduction zone interfaces reflect the influences of megathrust geometry and earthquake cycle processes such as enhanced elastic strain accumulation about seismic asperities and coseismic sense fault motion indicative of silent slip events or afterslip following large earthquakes.

The signature of an unbalanced earthquake cycle in Himalayan topography?

Fifty percent of the relative motion between the Indian and Asian plates is accommodated by active convergence at the Himalayan Range Front (HRF). Earthquake cycle processes on shallowly dipping HRF thrust faults generate large earthquakes (Mw ≥ 7) and contribute to the growth of HRF topography. Interseismic rock uplift rates reach a maximum north of the active Main Frontal Thrust and have been suggested to significantly influence the collocated convex bulge in HRF topography. Using geodetically constrained models of interseismic rock uplift rates and simple channel erosion rate laws, we show that convex channel profiles are predicted when interseismic deformation outpaces coseismic deformation. Applying this model to the observed elevation profiles of 20 HRF-spanning channels in Nepal yields a minimum mean residual elevation (72 m) if interseismic deformation has outpaced coseismic deformation by a factor of four. The long-term earthquake deficit required for the application of this model is consistent with some estimates of historical moment imbalance but requires temporally variable fault system activity. The spatial correlation between nominally interseismic rock uplift and the HRF topographic bulge may be explained by (1) a noncausal geometric coincidence, (2) geodetic observations of significant deformation not directly related to earthquake cycle processes, or (3) an unbalanced earthquake cycle at the HRF.

2009

Block modeling with connected fault network geometries and a linear elastic coupling estimator in spherical coordinates

Meade, B. J.,
and Loveless, J. P.

Bulletin of the Seismological Society of America
2009

Geodetic observations of interseismic deformation provide constraints on the partitioning of fault slip across plate boundary zones, the spatial distribution of both elastic and inelastic strain accumulation, and the nature of the fault system evolution. Here we describe linear block theory, which decomposes surface velocity fields into four components: (1) plate rotations, (2) elastic deformation from faults with kinematically consistent slip rates, (3) elastic deformation from faults with spatially variable coupling, and (4) homogeneous intrablock strain. Elastic deformation rates are computed for each fault segment in a homogeneous elastic half-space using multiple optimal planar Cartesian coordinate systems to minimize areal distortion and triangular dislocation elements to accurately represent complex fault system geometry. Block motions, fault-slip rates, elastic coupling, and internal block strain rates are determined simultaneously using a linear estimator with constraints from both geodetically determined velocity fields and geologic fault-slip rate estimates. We also introduce algorithms for efficiently implementing alternative fault-network geometries to quantify parameter sensitivity to nonlinear perturbations in model geometry.

From decades to epochs: Spanning the gap between geodesy and structural geology of active mountain belts

Allmendinger, R. W.,
Loveless, J. P.,
Pritchard, M. E.,
and Meade, B. J.

Geodetic data from the Global Navigation Satellite System (GNSS), and from satellite interferometric radar (InSAR) are revolutionizing how we look at instantaneous tectonic deformation, but the significance for long-term finite strain in orogenic belts is less clear. We review two different ways of analyzing geodetic data: velocity gradient fields from which one can extract strain, dilatation, and rotation rate, and elastic block modeling, which assumes that deformation is not continuous but occurs primarily on networks of interconnected faults separating quasi-rigid blocks. These methods are complementary: velocity gradients are purely kinematic and yield information about regional deformation; the calculation does not take into account either faults or rigid blocks but, where GNSS data are dense enough, active fault zones and stable blocks emerge naturally in the solution. Block modeling integrates known structural geometry with idealized earthquake cycle models to predict slip rates on active faults. Future technological advances should overcome many of today’s uncertainties and provide rich new data to mine by providing denser, more uniform, and temporally continuous observations.

Predicting the geodetic signature of Mw > 8 slow slip events

Elastic dislocation models of geodetic measurements above subduction zones have led to the identification of Mw ≈ 6.0-7.2 slow slip events (SSEs) that release elastic strain over periods of days to months, but great (Mw ≥ 8) SSEs have remained unidentified. We extrapolate observations of SSE duration and slip magnitude to show that slip velocity decreases with event magnitude and predict that the slip velocity of Mw ≥ 8 SSEs is ≤50 mm/yr. The slip velocity for great SSEs may never exceed the plate convergence rate and thus never produce a reversal in trench perpendicular motion. Instead, geodetically constrained estimates of apparent partial elastic coupling on subduction zone interfaces worldwide may be direct observations of ongoing Mw ≥ 8 silent earthquakes with durations of decades to centuries.

2008

Andean growth and the deceleration of South American subduction: Time evolution of a coupled orogen-subduction system

Present-day orography at the Andean margin is a result of isostasy, tectonic accretion, and erosional processes. The resulting excess mass of the Andes gives rise to frictional stresses on the seismogenic plate interface that resist the sinking of the subducting slab into the upper mantle. Thus, subduction rates should be sensitive to the time-dependent dynamics of a back-arc orogen, as well as erosional or accretional processes that affect orogen growth. Here we develop a two-dimensional coupled orogen–slab model that allows for the prediction of orogen size and plate motion in response to both tectonic and erosional forcing. We find that the frictional force exerted by the orogen on the subducting slab grows quadratically with orogen width and that the frictional resistance typically balances 10–50% of the slab pull force. The time evolution of the coupled orogen-subduction zone system is largely controlled by the rate of orogen growth, which is controlled by the rate of convergence and the erosivity of the climate state. In the case of the Andean margin, our models show that Miocene aridification leads to reduced erosion, increased orogen growth, greater frictional resistance to subduction, and, ultimately, to a 50% reduction in the convergence rate between the Nazca and South American plates.

Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments

To constrain the spatial distribution of erosion rates in the Waimea river watershed, on the western side of the island of Kauai, Hawaii, we calculate the frequency distribution of cosmogenic 3He concentrations ([3He]c) from helium isotopic measurements in olivine grains from a single sample of river sediment. Helium measurements were made in 26 aliquots of 30 olivine grains each. The average [3He]c from the 26 aliquots was used to estimate a basin-wide average erosion rate of 0.056 mm/yr, a value that is similar to erosion rates obtained from geochemical analyses of river sediments from tectonically stable landforms. However, forward models of cosmogenic nuclide production and sediment generation rates are inconsistent with the hypothesis that the observed [3He]c frequency distribution is the result of a homogeneous, basin wide, erosion rate. Instead, a distribution of erosion rates, from 0 to 4 mm/yr, may account for the observed frequency distribution. The distribution of erosion rates can be modeled by both non-linear slope- and curvature-dependent erosion rates with power law exponents ranging from 2.0 to 2.5. However, the spatial distribution of cosmogenic nuclides for slope- and curvature-dependent erosion rates are distinct, and we propose strategies to test further the extent to which erosion rates are controlled by the macroscale topographic features. These results demonstrate that the observed frequency distribution of cosmogenic nuclide concentrations in river sediments, combined with numerical modeling of erosion rates, can provide constraints on both the spatial variability of erosion rates in a drainage basin and the form of parameterized erosion laws.

Temporal variation in climate and tectonic coupling in the central Andes

McQuarrie, N.,
Ehlers, T. A.,
Barnes, J. B.,
and Meade, B. J.

Analog and numerical models predict a coupling between climate and tectonics whereby erosion influences the deformation of orogens. A testable prediction from modeling studies is the decrease in width of mountain ranges as a result of increased precipitation. Here we evaluate the effect of climate on a critically tapered orogen, the central Andes, using sequentially restored, balanced cross sections through wet (15\,^∘–16\,^∘S) and dry (21\,^∘S) regions of the orogen. In these regions, tectonics, basin geometry, and style of deformation are similar, allowing us to use variations in propagation (or changes in percent shortening) to evaluate whether alongstrike changes in width and morphology are climate driven in the north. Results indicate similar total percent shortening along the northern (40%) and southern (37%) sections, suggesting that a wetter climate has not limited the width (propagation) in the north. However, comparison of early (45–25 Ma) and recent (ca. 20–0 Ma) shortening indicates that early deformation produced 45% \pm 2% shortening of both sections, while recent deformation produced 41% \pm 2% (north) versus 32% \pm 2% (south) in the actively deforming Subandes. The latter suggests a coupling between climate and tectonics that began between ca. 19 and 8 Ma, and continues to 0 Ma, potentially limiting the width of the northern Subandes by 40 km.

2007

Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space

We present algorithms for analytically calculating the displacements, strains, and stresses associated with slip on a triangular dislocation element (TDE) in a homogeneous elastic half space. Following previous efforts, the solution is constructed as a dislocation loop where the deformation fields for each of the three triangle legs are calculated by the superposition of two angular dislocations. In addition to the displacements at the surface we derive the displacements and strains at arbitrary depth. We give explicit formulas for the strains due to slip on an angular dislocation, the calculation of angular dislocation slip components, a method for identifying observation coordinates affected by a solid body translation, and rules for internally consistent vertex ordering allowing for the superposition of multiple TDEs. Examples of surface displacements and internal stresses are given and compared with rectangular representations of geometrically complex fault surfaces.

Power-law distribution of fault slip-rates in southern California

The spatial partitioning of deformation in the continental crust and, in particular, at plate boundary zones is determined by the distribution of fault slip-rates. Analytic and numerical models of strain accumulation in the elastic upper crust have been divided into those that parameterize faulting as localized on a finite length fault system comprised of relatively few fast slip-rate faults, or as distributed throughout a continuum of relatively slow slip-rate faults. We demonstrate that in the southern California fault system, between the Pacific and North American plates, both geologically and geodetically constrained fault slip-rate catalogs obey a power-law frequency distribution. Using this empirically constrained scaling relationship we derive an analytic expression for the partitioning of potency accumulation rate, which determines the distribution and magnitude of slip localization. This model describes the kinematics of both micro-plate and continuum deformation models, and predicts that 97% of the deformation in southern California is accommodated on faults slipping at >1 mm/yr which is consistent with models of continental deformation which explicitly represent a large though finite number of deforming structures.

Present-day kinematics at the India-Asia collision zone

The collision of the Indian subcontinent with Asia drives the growth and evolution of the greater Tibetan Plateau region. Fault slip rates resulting from the relative motion between crustal blocks can provide a kinematic description of the distribution of present-day deformation. I construct a three-dimensional, regional-scale elastic block model of the India-Asia collision zone that is consistent with geodetic observations of interseismic deformation, mapped fault system geometry, historical seismicity, and the mechanics of the earthquake cycle. This mechanical model of the elastic upper crust yields a set of kinematically consistent fault slip rates and block motions that may serve to constrain dynamic models of continental crustal dynamics.

2006

Orogen response to changes in climatic and tectonic forcing

Despite much progress, many questions remain regarding the potential dynamic coupling between atmospheric and lithospheric processes in the long-term evolution of mountain belts. As a complement to recent efforts to discover the interrelationships among climate, topography, erosion, and rock deformation under conditions of mass-flux steady state, we explore orogen response to changes in climate and tectonic influx. We derive an analytical model that predicts a powerful climatic control on orogen evolution and distinct, potentially diagnostic, responses to climatic and tectonic perturbations. Due to isostatic compensation, the near-surface rock uplift rate during transients is tightly coupled to climate-modulated erosional efficiency. System response is approximately exponential, with a characteristic response timescale that is inversely proportional to the climate- and lithology-modulated erosional efficiency, and is largely insensitive to initial conditions, tectonic influx, and both the sign and magnitude of perturbations.

2005

Block models of crustal motion in southern California constrained by GPS measurements

We estimate slip rates on major active faults in southern California using a block model constrained by Global Positioning System measurements of interseismic deformation. The block model includes the effects of block rotation and elastic strain accumulation consistent with a simple model of the earthquake cycle. Our estimates of the right-lateral strike-slip rate on the San Andreas fault vary by at least a factor of 5, from a high of 35.9 ± 0.5 mm/yr in the Carrizo Plain to a low of 5.1 ± 1.5 mm/yr through the San Bernadino segment. Shortening across the Puente Hills Thrust and left-lateral slip on the Raymond Hill fault are consistent with both thickening and escape tectonics in the Los Angeles Basin. Discrepancies between geodetic and geologic slip rate estimates along the San Andreas and San Jacinto faults, as well as in the Eastern California Shear Zone, may be explained by a temporal change in fault system behavior. We find no substantial evidence for long-term postseismic relaxation and infer that the viscosity of the lower crust/upper mantle may be relatively high (η > 1e19 Pa s)

Spatial localization of moment deficits in southern California

The balance between interseismic elastic strain accumulation and coseismic release defines the extent to which a fault system exhibits a surplus or deficit of large earthquakes. We calculate the regional moment accumulation rate in southern California based on a fault slip rate catalog estimated from a block model of interseismic deformation constrained by GPS measurements. The scalar moment accumulation rate, 17.8 ± 1.1 e18 N m/yr, is approximately 50% larger than the average moment release rate over the last 200 years. Differences between the accumulated and released elastic displacement fields are consistent with moment deficits that are localized in three regions: the southern San Andreas and San Jacinto faults, offshore faults and the Los Angeles and Ventura basins, and the Eastern California Shear Zone. The moment budget could be balanced by coseismic events with a composite magnitude of Mw ≈ 8.

2004

Controls on the strength of coupling among climate, erosion, and deformation in two-sided, frictional orogenic wedges at steady state

Many important insights regarding the coupling among climate, erosion, and tectonics have come from numerical simulations using coupled tectonic and surface process models. However, analyses to date have left the strength of the coupling between climate and tectonics uncertain and many questions unanswered. We present an approximate analytical solution for two-sided orogenic wedges obeying a frictional rheology, and in a condition of flux steady state, that makes explicit the nature and sensitivity of the coupling between climate and deformation. We make the simplifying assumption that the wedge grows in a self-similar fashion consistent with Airy isostasy such that topographic taper is invariant with orogen width, tectonic influx rate, and climate. We illustrate first how and why the form of the erosion rule matters to orogen evolution and then derive a physically based orogen-scale erosion rule. We show that steady state orogen width, crest elevation, and crustal thickness are controlled by the ratio of accretionary flux to erosional efficiency to a power dictated by the erosion process. Remarkably, we show that for most combinations of parameters in the erosion law, rock uplift rate is more strongly controlled by erosional efficiency than it is by the accretionary flux. Further, assuming frontal accretion with no underplating, the spatial distribution of erosional efficiency dictates the relative rock uplift rates on the pro-wedge and retro-wedge and the time-averaged trajectories of rocks through the orogen. The restriction to invariant frictional properties is conservative in these respects; systems subject to positive feedback between erosion and rheology will exhibit even stronger coupling among climate, erosion, and deformation than shown here.

Viscoelastic deformation for a clustered earthquake cycle

The clustering of earthquakes in time on the same fault affects the rate and pattern of interseismic deformation. We develop a simple analytic viscoelastic model of the surface velocity field through a clustered earthquake cycle by superposing the velocities of individual earthquake cycles of constant period but varying phase. Velocity profiles prior to and after an earthquake show a wider range of behavior than they do for earthquake cycles with constant repeat times. These new types of behavior provide possible explanations for discrepancies between geologic estimates of long-term fault slip rates and slip rate estimates from geodetic data.

2002

Estimates of seismic potential in the Marmara Sea region from block models of secular deformation constrained by global positioning system measurements

Meade, B. J.,
Hager, B. H.,
McClusky, S. C.,
Reilinger, R. E.,
Ergintav, S.,
Lenk, O.,
Barka, A.,
and Ozener, H.

Bulletin of the Seismological Society of America
2002

We model the geodetically observed secular velocity field in northwestern Turkey with a block model that accounts for recoverable elastic-strain accumulation. The block model allows us to estimate internally consistent fault slip rates and locking depths. The northern strand of the North Anatolian fault zone (NAFZ) carries approximately four times as much right-lateral motion ( 24 mm/yr) as does the southern strand. In the Marmara Sea region, the data show strain accumulation to be highly localized. We find that a straight fault geometry with a shallow locking depth of 6-7 km fits the observed Global Positioning System velocities better than does a stepped fault geometry that follows the northern and eastern edges of the sea. This shallow locking depth suggests that the moment release associated with an earthquake on these faults should be smaller, by a factor of 2.3, than previously inferred assuming a locking depth of 15 km.

2001

Present day kinematics of the Eastern California Shear Zone from a geodetically constrained block model

McClusky, S. C.,
Bjornstad, S. C.,
Hager, B. H.,
King, R. W.,
Meade, B. J.,
Miller, M. M.,
Monastero, F. C.,
and Souter, B. J.

We use Global Positioning System (GPS) data from 1993–2000 to determine horizontal velocities of 65 stations in eastern California and western Nevada between 35° and 37° N. We relate the geodetic velocities to fault slip rates using a block model that enforces path integral constraints over geologic and geodetic time scales and that includes the effects of elastic strain accumulation on faults locked to a depth of 15 km. The velocity of the Sierra Nevada block with respect to Nevada is 11.1±0.3 mm/yr, with slip partitioned across the Death Valley, (2.8±0.5 mm/yr), Panamint Valley (2.5±0.8 mm/yr), and Airport Lake/Owens Valley (5.3±0.7/4.6±0.5 mm/yr) faults. The western Mojave block rotates at 2.1±0.8°/My clockwise, with 3.7±0.7 mm/yr of left lateral motion across the western Garlock Fault. We infer 11±2 mm/yr of right lateral motion across the Mojave region of the Eastern California Shear Zone.

The current distribution of deformation in the western Tien Shan from block models constrained by geodetic data

We interpret Global Positioning System measurements of interseismic deformation throughout the western Tien Shan in the context of a block model which accounts for important geologic features (faults) and physical processes (elastic strain accumulation.) Through this analysis we are able to quantify the amount of deformation localized on active structures. In the central part of the belt the Dzhuanaryk fault zone appears to be the most important thrust fault, accommodating nearly five millimeters per year of north-south shortening across it. Conversely, the most widely recognized strike-slip fault in the region, the Talas-Fergana, is found to have very little of the previously estimated right lateral motion.