An Earthquake Model Based on Fatigue Mechanism—A Tale of Earthquake Triad

Diandong Ren^1,2

¹School of Electrical Engineering, Computing and Mathematical Sciences, Curtin University, Perth, Australia.
²Australian Sustainable Development Institute, Curtin University, Perth, Australia.

DOI: 10.4236/gep.2019.78021                

Abstract

Earthquakes are the result of strain build-up from without and erosion from within faults. A generic co-seismic condition includes merely three angles representing respectively fault geometry, fault strength and the ratio of fault coupling to lithostatic load. Correspondingly, gravity fluctuation, bridging effect, and granular material production/distribution form the earthquake triad. As a dynamic component of the gravity field, groundwater fluctuation is the nexus among the three intervened components and plays a pivotal role in regulating major earthquake irregularity: reducing natural (dry) inter-seismic periods and lowering magnitudes. It may act mechanical-directly (MD) through super-imposing a seismogenic lateral stress field thus aiding plate-coupling from without; or mechanical-indirectly (MI) by enhancing fault fatigue, hence weakening the fault from within. A minimum requirement for a working earthquake prediction system is stipulated and implemented into a well-vetted numerical model. This fatigue mechanism based modeling system is an important supplement to the canonical frictional theory of tectonic earthquakes. For collisional systems (e.g., peri-Tibetan Plateau regions), MD mechanism dominates, because the orographically-induced spatially highly biased precipitation is effectively channeled into deeper depth by the prevalence of through-cut faults. Droughts elsewhere also are seismogenic but likely through MI effects. For example, ENSO, as the dominant player for regional precipitation, has strong influence on the gravity field over Andes. Major earthquakes, although bearing the same 4 - 7 years occurrence frequency as ENSO, have a significant hiatus, tracing gravity fluctuations. That granular channels left behind by seamounts foster major earthquakes further aver the relevance of MI over Andes. Similarly, the stability of the Cascadia fault is found remotely affected by Californian droughts (2011-15), which created a 0.15 kPa/km stress gradient along the Pacific range, which also is the wave guide.

Keywords

Earthquakes, Erosion, Co-Seismic

1. Earthquakes as Frictional Phenomena and Earthquake Triad

Tectonic earthquakes are frictional phenomena between contacting plates (Scholz, 1998; Wang & Hu, 2006; Wang & He, 1999) . It is a game of gravity- (or more generally compressive stress-) aided friction resisting against the ever-increasing plate-coupling stress ( $\tau$ ) resulted from differential plate motions. For an ideal configuration (Figure 1(a)) of uniform fault geometry (i.e., constant slope angle $\theta$ and unlimited width in the transverse direction) and uniform material (so that maximum frictional coefficient ${\mu }^{\prime }=\tan \left({\theta }_f\right)$ being a constant), the co-seismic condition is

${\varphi }_c=\theta +{\theta }_f$ . (1)

where repose angle $\varphi =\arctan \left(\tau /G\right)$ , with G being compressive stress load on the shear (fault) zone. Subscript c means critical value. $\theta$ is the slope of the fault zone, determined by the geometry of the plate interface, and ${\theta }_f$ is the maximum static friction angle corresponding to fault strength. Values of ${\theta }_f$ depend on wall rock material properties and lithologies. During inter-seismic stage, the repose angle increases steadily and approaches the sum of fault slope angle and the maximum static friction angle. Driving stress increase and the overall resistive stress decrease (weakening of the fault, or ${\theta }_f$ decrease) both are seismogenesis (Scholz, 1998) . In reality, erosion from within and burden from without likely are simultaneous. Earthquake cycle is a stress adjustment cycle. Closely related to three angles in Equation (1) are the earthquake triads: Gravity fluctuation of the overriding plate, bridging effect and fatigue. With fault strength set as constant, fault geometry and plate-coupling stress determine a natural (limiting) earthquake occurrence frequency. However, actual occurrence seldom obeys this limit and usually occurs far ahead of schedule, due to various fault-weakening factors.

Since displacements caused by each earthquake (ranging between 0.5 - 3 m) are at least 3 orders of magnitude smaller than the dimension of the seismically locked fault zone, which could well be several hundred kilometers, material in the locked fault zone is well-seasoned: all having experienced numerous times of wearing and tearing. For predicting near future earthquakes, the geological backgrounds, including material properties and fault geometry, can be safely assumed to be constant. Thus, the natural limiting occurrence frequency over a fault is rather stable. The irregularity in earthquake occurrence is primarily due to fault weakening mechanism. This study focuses on more dynamic fluctuations affecting the triad of earthquake. Following discussions use subduction fault as prototype. By substituting compressive stress for lithostatic loading, same reasoning is viable for strike-slip faults. Strictly, simulating an earthquake event includes estimates of its timing (Section 2), epicenter location and magnitude (Section 2.1). Actual full 3D numerical model, with detailed parameterization outlined (Section 3) is used in the simulations discussed in Section 4. In addition to the parameterizations of material rheology (visco-elastic mantle, granular debris as produced by seismic events of all kinds, and the brittle elastic crust), Section 3 also is dedicated to the setting of the fault-following model grids, static geologic properties (e.g., plate interface geometry and physical property of each layer of medium).

2. Complexity of the Earthquake Problem: Under-Determination and Interactions of the Triad

Different segments of a realistic heterogeneous fault zone (Figure 1(b)) have different slope angles and maximum friction coefficients and also experience different loads. The driving stress distributed to each segment approaches its maximum affordable resistive stress gradually, like the saturation process of porous media. Fault “stress saturation” is thus analogously defined here as the ratio of driving stress to maximum affordable resistive stress.

$S_i={\tau }_i/\left[G_i\tan \left({\theta }_i+{\theta }_{f,i}\right)\right]=\tan \left({\varphi }_i\right)/\tan \left({\theta }_i+{\theta }_{f,i}\right)$ (2)

Consequently, an S value of 0 corresponds to full locking while an S of 1 corresponds to creeping. The overall stability of a realistic fault becomes ${\bar{S}}_i=\left[{\displaystyle \underset{i}{\sum }{\tau }_i}\right]/\left[{\displaystyle \underset{i}{\sum }{G}_i\left({{\mu }^{\prime }}_i+\tan {\theta }_i\right)/\left(1-{{\mu }^{\prime }}_i\tan {\theta }_i\right)}\right]\le 1$ . S is thus a useful diagnostic index for fault stability.

2.1. Multi-Solution in Seismogenesis

In reality, earthquakes involving the saturation of the entire shear zone seldom occur, because actual tectonic plates have difficulty in maintaining rigidity over extended distance, manifested as the existence of through-cut faults and many small magnitude earthquakes within close proximity to each other. Consequently, many places do not contribute their shares of resistance and depend on neighbors to stay in place. Depending on how $\tau$ is distributed on the fault zone, there can be various “legal” combinations of saturated and partially saturated segments along the fault zone. In reality, earthquake occurrence can be compounded by any factors affecting $\tau$ , G, and ${\theta }_f$ (or fault strength ${\mu }^{\prime }$ ) in Equation (2). Multi-solution in seismogenesis essentially is resulted from the lack of enough information as constraints on, or an under-determination issue of lateral interactions. Realizing that earthquakes, as a result of relative motion of plates involved, inevitably have footprints on the mass (Wallace et al., 2009; Ren et al., 2015) , energy (Gao & Wang, 2014) and momentum fields ( Wang et al., 2012 and references therein), efforts for earthquake prediction are increasingly based on observational facts of these physical parameters and their fields.

(a)(b)(c)(d)

Figure 1. Diagram of an ideal crust front/fault (a) and a plane view of a more realistic plate interface (b). Panel (a) is a highly idealized continental plate overlying oceanic plate configuration. Black bold arrows in (a) show the opposing oceanic plate and the overlain continental plate. Stress analysis is performed for an arbitrary point (P) at the shear zone. $\mu$ is the static friction; G is the gross gravity of the overlain plate section; $\tau$ is the driving stress resulted from tectonic motion. Static friction and gravity work together to contest against continental driving stress. Once the driving stress overcomes gravity-aided maximum static friction, a co-seismic sliding occurs. Red-toned color shades illustrate stress saturation. Granular material, often in the form of fault gouge, inter-laces the plate-interface (b). Either because geometrical irregularity of the contacting interface, or because the existence of density anomaly (e.g., cavity) above, there are places (“pillars” or strong contact regions) that take a disproportionally larger load than the neighboring regions (Bridging effect). Plates are “riveted” together by these “sticky” spots. These strong contact regions tend to wear down faster, through generating granular debris more actively. “Asperities” and “barriers” (Wallace et al., 2004) can both be the “pillars” defined here. Panel (c) shows the vertical resistive stress σ_zz fields around a region with significant bridging effects (i.e., with cavity caused gravity anomaly). Note that the under-belly of the overlain plate is in weak contact with lower plate and also is a region with active granular material production (that region experiences strong tensile stress). The fault creeps primarily by wearing geometrical irregularities or heterogeneous “rises/bumps/sticky spots”. The epicenter is the ensemble (weighted by frictional heat released at these regions) mean location of the sticky points being co-seismically unlocked. Real cavity usually occurs at much shallower depths. Fault interface underneath still feels the bridging effects by producing a corresponding uneven granular material generation pattern. Panel (d) is SEGMENT simulated global distribution of inter-pate granular material thickness (m).

The co-seismic release of potential energy is transformed mainly into heat through inter-plate frictional force. A large, deep-burying area approaching saturation simultaneously, un-locking systematically and releasing the stored elastic energy coherently generally produces large magnitude earthquakes (Wang & Bilek, 2011) . The contact surface of actual tectonic plates can have sophisticated topography and may involve material heterogeneity (Figure 1(b)). Therefore, it is difficult to predict the future rupture propagation without knowing the details at the interface. Recent studies (Gao & Wang, 2014; Wang et al., 2012; Wang & Bilek, 2011) indicate that widespread and smooth contact between two relatively weaker plates harbors large earthquakes—a situation usually satisfied at the megathrusts. As an opposite case of rough contact, the subducting seamounts generally discourage the generation and propagation of large ruptures. This is because the seamounts consume a very small fractional of the plate interface, acting more effectively as a grinder for the overlain plate rather than as a blocker (Figure 1(b); Wang & Bilek, 2011 ). Downstream of these seamounts, there are stripes of granular debris composed by the material scrapped off both the sea mounts and the overlain continental plates. These existing findings learn further support to the importance of lateral coherence among different segments of the fault in determining the magnitude of the earthquakes (Equation (2)).

2.2. Intertwining Earthquake Triad

One important factor causing non-uniform distribution of driving stress among segments is the existence of macroscopic crevasses of voids/cavity within the plates (Rowe et al., 2012) . In reality, bridging effects (Ren, 2014a) —load directly above can be partially carried by neighboring segments (or vise-versa)—add yet another layer of complexity. Bridging effect contributes to gross resistive stress only when the fault interface has non-uniform friction coefficients. Otherwise, if the interface has same roughness ${\mu }^{\prime }$ , the net contribution to resistive stress from bridging effects may be minuscule due to the canceling effects from neighboring regions that is relieved portion of their lithostatic stress.

From Equation (2), gravity of the overriding plate always serves as a stabilizing factor. The reduction of gravity (e.g., from extended droughts) contributes to instability. Fluctuations in gravity field not only have local stability consequence, they also affect neighboring regions through bridging effect and a unique fatigue process of the fault zone: granular material (GM, Jop et al., 2006; Ren et al., 2008 ) generation and transportation, a common fault-weakening mechanism. Production of granular debris is the primary form of fatigue of the fault interface. Regions of strong contact, likely being a result of bridging effect, encourage active granular material generation. Ironically, granular material has much smaller viscosity and acts as lubricant for the plates involved. Regions carrying more loading are now footed on slippery floor (i.e., more weight is loaded on smooth contacts and less is loaded on rough contacts along the shear zone). Thus, bridging effect, through affecting granular material generation and redistribution, influences the frictional properties of the fault interface and achieves a negative spatial covariance between loading and fault strengths, being seismogenic (i.e., contributing to fault instability) according to Equation (2) (i.e. reduction of the ${\displaystyle \underset{i}{\sum }{\tau }_i}$ ). Sources of bridging effects are plenty. For example, large scale structures in the overriding plate such as mountains and valleys signifies bridging effects at depth on the seismogenic fault zone (Bejar-Pizarro et al., 2013) . Even for regions without topographic features, existence of cavities/crevasses at depth (not necessarily all the way to the locked interface; the vertical change in loading transfers readily to horizontal stress in the fault circumstances) causes non-negligible bridging effects on faults (to be further detailed in Section 3.1.1). Bridging effects from static geological background, although may cause variations in the frictional properties of the seismogenic zone, vary slowly and are not responsible for short term variations in major earthquake occurrence. Large scale variations in groundwater, however is more dynamic and is a viable candidate for explaining earthquake occurrence irregularity on decadal to centennial time scales. Cyclic forcing is most effective in causing fatigue (Ren & Leslie, 2014) . Cyclic groundwater fluctuations, especially when acting in concord with the resonance frequency of the plates, assist granular material generation.

Thus, the earthquake triad is interlinked and hence reinforces one another through groundwater fluctuation as the nexus. Groundwater-aided granular material generation (G³) is the mechanism modulating earthquake cycles. Detailed discussion on fatigue in rock interface, bridging effects and their physical parameterizations in the Scalable Extensible Geofluid Modeling framework for ENvironmenT Intelligence System (SEGMENT, Ren, 2014a; Ren et al., 2008; Ren et al., 2012 ) can be found in Section 3, together with the sources of input datasets. Through solving conservation equations of mass, momentum and energy, the SEGMENT provides 3D distribution of strain, stress, and full three components of velocity in the simulation domain. How a major earthquake is identified using these prognostics is described in Section 3, irrespective of the physical domain they reside (on the main thrust interpolate or not) and the morphology (e.g., megathrust earthquakes or simultaneous rupture of several strike-slip faults as the 2012 Indian Ocean earthquake). It is noteworthy that variables in Equations (1) and (2), such as the degree of stress saturation, are merely diagnostics of SEGMENT instead of being prognostics. Through these diagnostics, one can vividly view the stress build-up during the inter-seismic stage. Section 3 also includes detailed discussions of recondite approaches in obtaining present global distribution of granular material thickness and the future evolution.

3. Earthquake Triad (3.1) and Their Representation in SEGMENT Modeling System (3.2)

Fault geometry, although being the far-largest factor in determining the timing of major earthquakes, are looked upon as geological background here. The super-imposed gravity fluctuations are primarily responsible for earthquake irregularity. Following discussion focuses on the bridging effect and the fault weakening consequences from granular material (GM) formation. Groundwater is not singled out as a subsection because it is the nexus among the components and it is natural to discuss it as an integrated component throughout the text. The following is a description of a minimum requirement for a working earthquake prediction system for realistic 3D fault configuration—the Groundwater aided Granular material Generation (G3) framework and its implemented in the well-vetted SEGMENT system.

3.1. Earthquake Triad

3.1.1. Bridging Effects at the Fault Zone Make Instability “Non-Local”

That different portions of the fault/shear zone contribute different resistive stresses primarily is due to uneven loading (the variations in G_i along the shear zone). There are many causes for loading deviation from the lithostatic values. For convenience, a generic term “bridging effects” (Figure 1(c); Ren, 2014a ) is assigned to this phenomenon. It is a disturbance to the vertical resistive stress field (R_Z). The “G” in Equations (1) and (2) is a special case of R_Z when bridging effect is insignificant. In principle, disturbances at any depth above the fault zone can be felt at the plate interface. It must be noted that, for a point source (of horizontal dimension much less than fault dimension), the influence range of bridging effect in each horizontal direction expands linearly. So the footprint increases at a square power of depth. Consequently, the amplitude of the disturbance decreases at a rate of depth-squared. Thus, the closer the disturbance source is to the shear zone, the more effective in superposing onto R_Z and contributing to motion tendency of this segment of the shear zone. According to this reasoning, localized load disturbances such as man-made infrastructure (e.g., dams, high raisers) are far less important than a cavity of similar size existing at 20 km depth in disturbing the load distribution of a shear zone at ~30 km depth. The limiting size (L) of cavity at depth H depends on rock tensile strength (f_c) and density ( $\rho$ ) according to a square root law: $L=b\sqrt{f_{c}H/\rho g}$ , with b a cavity geometry dependent constant. It is apparent that larger dimension cavities are allowed at deeper depths of the brittle crust. At deeper depth (>~30 km), the thermodynamic reduction in tensile strength of rock dominates and causes drastic reduction in L. As a result, the cavities on the locked interface are of sub-meter scale at most (Ujjie et al., 2007; Cowan et al., 2003) . Thus, bridging effects from inter-plate geometry irregularity may be only non-negligible for regions surrounding sea mounts for earthquake evolution and lifecycles. These features have clear remote-sensing signature, especially when experiencing transition in size as a result of cave-in (usually a result of continued granular material formation). Loading irregularity is thus the most prevalent form of bridging effects. Load fluctuations of shallower level, only when occurring at large enough spatial scale, affect fault stability. In this sense, regional scale groundwater fluctuations, even only resides in upper couple kilometers depth, may have significant impact on fault stability. The deeper they reach (such as around the through-cut faults in the Longmenshan Fault system (to be detailed in the discussion section) the more direct and efficient they are in affecting fault stability.

The bridging effect contributes to gross resistive stress only when the plate interface has non-uniform frictional properties. To have a net effect on resistive stress, the extra loading should be negatively coupled with the spatial variation of fault strength (i.e., their spatial co-variance being negative). Granular material (Jop et al., 2006; Ren et al., 2008) formation and transportation is such a nexus. Through bridging effects, neighbor’s weight is partially transferred over regions with strong contacts (like “pillars” of a bridge). There may not be any macroscopic voids; by “pillars” it is meant only of strong compressive stress region, or strong contact regions. These regions of strong contact encourage granular material formation (a consensus of yielding criteria of brittle rocks; Rowe et al., 2012 ). With the formation of granular material, which has viscosity (<10⁴ Pa S) many orders of magnitude smaller than polycrystalline rocks (>10²¹ Pa S), the shear zone with strong contacts weakens and this is seismogenic. For example, if only 5 mm granular material is put between rock plates under the confining pressure similar to the locked fault zone, the equivalent frictional coefficient, or fault strength drops to <10⁻⁴, orders of magnitude smaller than solid-against-solid rock friction.

Sea mounts are a type of inter-plate “pillars/asperities” that grind upper plate and also wear out themselves. “Bridging effect” actually is the controlling factor for granular material generation. Bridging effect caused co-seismic and interseismic stress irregularities has strain/deformation consequences at all depths up to the earth surface, where the near surface displacements (directly associated with strain buildup and release) can now be accurately measured by Global Positioning System (GPS; Chlieh et al., 2011; Sun et al., 2014; Loveless & Mead, 2010 ) and Synthetic Aperture Radar interferometry (InSAR; e.g., Weston et al., 2012; Cavalie & Jonsson, 2014 and references therein) instruments. From Equation (2), it contributes to stability if more weight is loaded on rough contacts and less is loaded on smooth contacts along the shear zone. Bridging effect, through affecting granular material production rates, influences the frictional properties of the fault interface and achieves the seismogenic effect of “more weight is loaded on smooth contacts and less is loaded on rough contacts along the shear zone”. Sources of bridging effects are plenty. For example, large scale structures in the overriding plate such as mountains and valleys signify bridging effects at depth on the seismogenic fault zone (Bejar-Pizarro et al., 2013) . However, topographic features caused bridging effects and hence variations in the frictional properties of the seismogenic zone cannot explain the irregularity of earthquakes that occur on decadal to centennial scales. Large scale variations in groundwater, however is a viable candidate. Thus, fluctuations in gravity field not only have local stability consequence (Equation (2)), they also affect neighboring regions (laterally) through bridging effect and the unique fatigue process of the fault zone.

3.1.2. Fatigue in Rock Interface—Generation of Granular Material

Two objects, when being pressed against each other, wear and tear at the interface, a manifestation of material fatigue. The interfaces of tectonic plates are of no exception. One proof is the existence of widespread fault gouge and pseudotachylyte melt in wall rocks (Rowe et al., 2012) . Seismic activities, irrespective magnitudes, produce granular material. The generation rate of GM however is a functional of at least three factors: 1) the degree of imbalance of the three principle stress components, in a form of yielding criteria for brittle material; 2) the existing amount of granular material or the current state of GM; and 3) forming granular material consumes elastic energy accumulated in the plates, through forming new surfaces under huge confining pressure of the fault environment. Point (1) is supported by recent findings of frictional melt along the fault surface. It happens at only a small fraction of the total rupture surface area but always begins at asperities, the effectively strongest points on the fault surface. Cyclic forcing is most effective in causing fatigue. Groundwater fluctuations, especially when acting in concord with the resonance frequency of the plate segment, enhance granular material generation. The second requirement is because, as granular material builds up at the interface, it hinders further production of granular material in exactly the same way accumulated saw dust hinders further cutting of a piece of wood. This is especially relevant for the inter-plate granular material because, compared with its parent rock it has larger porosity and easily fills up the limited inter-plate voids. Earthquakes are such an effective mechanism for removing and re-distributing existent granular material, specifically by thermal pressurization or acoustic fluidization. Granular material production and redistribution achieve a negative spatial covariance between loading and effective frictional coefficients, being seismogenic according to Equation (2).

In addition to the critical role played in granular material generation, groundwater, as an extra loading on the overriding plate, is by itself a stabilizing factor, according to Equation (2). Reduction of groundwater, usually as a result of extended wide-spread droughts, is however seismogenic. Although the weight of groundwater itself is negligible compared with the weight of the overriding plate, it is of the same order of magnitude as the residual stress (between plate-coupling compressive stress and gravity-aided friction). That is, the inter-seismic stage is actually in a very tricky balance (Wallace et al., 2017) . Fluctuation of groundwater, through superimposing a lateral/horizontal stress gradient, may play a critical role in the timing of major earthquakes. Fortuitously, fluctuations in groundwater can now be remotely sensed by gravity satellites such as the Gravity Recovery and Climate Experiment (GRACE; Famiglietti & Rodell, 2013; Ren, 2014b; Voss et al., 2013 ) and by InSAR measurements (Liu et al., 2017) .

Centered on the master co-seismic criterion (Equation (1)), the triad determining earthquake irregularity is discussed. It is now clear that, because of the heterogeneity in driving stress distribution, and the overlain loading condition, unlocking is often not completed at once. Rather, different geological locations may unlock in certain orders. The ensemble center of the ruptures is the epicenter and the magnitude is directly related to the rupture size. Geological backgrounds such as fault geometry and the plate motion speeds evolve rather slowly. The irregularity in earthquake occurrence mainly is resulted from the inter-plate granular material production and transportation. Large spatial scale fluctuations in groundwater (i.e., those arising from large scale droughts and floods), through bridging effects, are a suitable mechanism in granular material generation. In a certain sense, the climate fluctuations, through affecting groundwater loading pattern, affect the earthquake cycles. Bridging effects by themselves may not cause fault weakening. But with the associated granular material generation, they weaken faults through forming a negative spatial covariance pattern between loads and friction coefficients. The generation and transportation of granular material at the plate contact reconcile above-mentioned counter-intuitive characters of major earthquakes (e.g., as enumerated in Wang (2015) ) and also are critical in orchestra of small-scale ruptures to generate major earthquakes.

This study focuses on subduction zones (e.g., Andes and Cascadia) because they are responsible for the largest earthquakes and also because they have significant footprints on the gravity (mass) and thermal (Gao & Wang, 2014) fields. The continental collisions (continental against continental, e.g., Himalayas) here are not differentiated from subductions (oceanic against continental plates). Despite the improved measurements of the Earth’s structure and rupture geometry, the timing of major earthquakes remains a scientific puzzle (Pritchard & Simons, 2006; Sawai et al., 2004) . In this research, first order schemes representing earthquake triad are implemented in a sophisticated 3D modeling system, to estimate the spatio-temporal scales of major earthquakes.

3.2. Numerical Model Setup, Parameterization of the Earthquake Triads, and Data Acquisition

Above discussed theoretical concepts, or the Groundwater aided Granular material Generation—G³ hypothesis is implemented into a global spherical crown version of the SEGMENT, with the stress decomposition into driving stress and resistive stress follows the convention of Van der Veen and Whillans (1989) . The SEGMENT solves the equations of mass, momentum and energy conservation in a 3D spherical geometry. The model prognostics include three components of velocity, temperature, and six components of the deviotoric stresses. Their evolution is depicted by the time marching of the governing equations. Variables in Equations (1) and (2), such as the degree of stress saturation, are just diagnostics of SEGMENT instead of being prognostics. Through these diagnostics, one can vividly view the stress build-up during the inter-seismic stage. A major earthquake in model is identified as more than 10⁴ km² adjacent model grids in the horizontal dimension possess macroscopic motion velocities. Once an event is identified, the released energy is estimated within a time window fully covers the time period with macroscopic motion. An epicenter is diagnosed as the geological locations of the saturated grids, weighted by the respective energy release. Degree of stress saturation S is not directly used to identify earthquakes also because the grid elements are not isolated; there are interactions among adjacent grids. If the unlocking/rupture at a grid point releases a large amount of energy, a neighboring grid point, even still quite away from saturation, can be unlocked and set into motion. In fact, some very large earthquakes are formed because there are some sections of the fault (not necessarily adjacent to each other) get simultaneously ruptured and the released energy triggered neighboring sections, resulting in an upward spiral that causes a domino effect in propagating away of the rupture region.

3.2.1. Grid Stencil of SEGMENT-Fault Zone Oriented Coordinate System

For this study, the original 3D spherical model of the upper ~4 km medium is now applied into a 500 km thick lithospheric material (including crust and upper mantle). The horizontal resolution is generally 30 km at the earth surface and refined to 2 km at subduction zones (Figure 2). Vertical resolution varies as 101 stretched vertical layers are used to represent the 500 km depth. The stretching scheme is used so that the shear zone gets better represented. In setting up grids in the model simulation domain, special attention is paid so that the upper surfaces of the oceanic plates are locally parallel to the model grids spanned surface (there always is a vertical level that has surface parallel to the oceanic plate’s upper surface).

3.2.2. Static Geological Parameters and Viscosity Parameterizations

Like other 3D tectonic models (e.g. subduction zone models), an accurate representation of the fault geometry (e.g., interface geometry such as subduction thrust geometry). This requires high resolution geodic and seismic data. The CRUST1.0 (Laske et al., 2013) is used to set up the density values over the grids within the crust layer (oceanic and continental crusts). Underneath bottom of the crust layer (shown in CRUST 1.0) is assumed to be upper mantle material (silicates) and given a uniform viscosity of ${\eta }_0=\nu ={10}^{19}\text{Pa}\cdot \text{S}$ and a density of 3.3 × 10³ kg/m³ as base values and a pressure- and temperature-dependent modifier (to be detailed very soon). In setting up the grids above present sea level, the Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation map (SATER DEM, asterweb.jpl.nasa.gov/) and ETOPO1 (downloadable from https://www.ngdc.noaa.gov/mgg/global/global.html) bathymetry are referenced to. For faults with intensive gravity measurements, for example the Longmenshan Fault Zone (LFZ) of the eastern margin of the Tibetan Plateau, there was a detailed gravity survey after the M_w7.9 Wenchuan earthquake on 12 May 2008. Finer resolution data on elastic plate thickness, density and rigidity profiles are available and used in replacement of the CRUST1.0-derived parameters. Similarly, the Andes survey (Andes Geological Survey, 2006) and data pointed out by references therein provided information also is merged into the coarse resolution CRUST 1.0. For the Hikurangi subduction zone (New Zealand), the revised interface geometry (Williams et al., 2013) are used instead of the CRUST 1.0. For the Puysegur trench, we adapted the geometry data from Liu and Bird (2002) . For the Cascadia, Mexico, and Japan Trench, the higher resolution interface geometry are obtained from Gao and Wang (2017) . Personal communication with Z. Liu also provided higher resolution geometry data along the west frontal range of N. America. Other regions we have to tolerate the lower resolution interface geometry provided by SLAB 1.0 (Hayes et al., 2012) .

Based on the fact that the co-seismic ruptures consume only a portion of the plate interface (Lay, 2015) to the center and the down-dip section creeps, the upper and lower portions of the plate should be treated respectively as elastic

Figure 2. Horizontally and vertically stretched (so that the fault region and especially the plate interface zone is better represented) grids in SEGMENT, as delineated by black vertical lines and horizontal curves. The “locked” region gets the finest spatial resolution. In the corresponding calculational space, two shape factors (i.e., Jacobian operators) are involved in the governing equations. Background map showing subduction fault geometry is adapted from Wang (2015) .

and visco-elastic material. Grids in the upper crust are considered elastic and elastic moduli are specified using canonical values. Grids in the lower crust and upper mantle away from the shear zone are considered to be visco-elastic (Wang et al., 2012) but getting more rigid upwards intending a smooth transition in the vertical domain. This is achieved by considering the thermal structure (Gao & Wang, 2014) in the parameterization of material viscosity

$\eta \left(T,P\right)={\eta }_{0}D\exp \left[\frac{Q+PV}{RT}\right]$ . Here P is confining pressure, T is absolute

temperature, and D depends on the rock composition, grain size and fluid content. Here Arrhenius’ activation energy Q, specific volume V and universal gas constant R were assumed constants in the model. In response to different stress conditions, wall rocks can be compressed or stretched within a narrow range before onset of brittle fatigue. Arrays are allocated to record the amount of granular material formed at each grid cube. The shear zone is a granular zone of variable thickness (Figure 1(d)). In most cases, it is only a small friction of the shear zone layer that is granular material (also recorded in arrays). The parameterization of granular material viscosity follows Jop et al. (2006) , as implemented in SEGMENT (Ren et al., 2008) . The shear-thinning viscosity structure (Equation 8 in Ren et al., 2008 ) implies a shear-localizing effect from granular material accumulation. The present amount of granular material (Figure 1(d)) is deduced from surface velocity using the method in Ren et al. (2012) . Future variation of granular material is diagnosed from Equation (3) (described in next subsection on granular material generation; immediately follows). The granular ensemble size is deduced from its parent rock type and the associated grain sizes. In addition to the brittle interface of oceanic and continental crusts, GM also exist in the underbelly of the overlying continental crust just updip of the mantle wedge corner, a byproduct of (forearc) mantle wedge serpentinization (e.g., Hyndman et al., 2015; Audet & Burgmann, 2014 ). The GM material amount in this zone can be estimated from the clustering of episodic tremor and accompanying slow slip events and their magnitudes. Re-mapping of the grid stencils are automatically performed only if granular material accumulated exceeds 20 m or earthquake caused displacement over 20 m.

Setup of the lower boundary stress conditions is critical for the simulation of stress buildup inside the simulation domain. The GPS measurements of plate motion speeds $V$ are used to fine tune the parameters of upper mantle layers so that lateral stress (basal drag) exerted on the lowest model layer at the active zones (mantle driven plate motion) and passive zones (plate driven mantle motion) nearly balance. From the vertical profile of the horizontal velocity components, it is clear that for most regions, the mantle motion is passively forced by the overlain plates and only in the active thermal hot spots (for Cascadia, the Yellowstone area), the boundary layer flow structure indicates that it is the mantle flow that exerts a strong drag to the overlain plate. Reflected in the flow structure, mantle flow speeds increase away from the lithosphere-asthenosphere boundary (LAB) rather than being a maximum at the LAB. Therefore, the driving stress comes primarily from lateral stress from neighboring grids within the lithosphere, rather than from the local basal drag from mantle (ultimately it is from the convective mantle but local to subduction zones it is not). The upper boundary condition is assumed to be stress-free (i.e., $\tau \cdot n=0$ at upper surface). The model spin-up from ad hoc initial conditions is aided by advanced data assimilation schemes using the remotely sensed plate motion, gravity (weather signal pre-filtered) signal and thermal flux measurements available. Accurate flow field is critical for an accurate estimation of the stress build-up between plates. Figure 3 shows present plate motion velocities at two vertical levels: 0 - 50 km and 100 - 250 km depth averaged. Compared with GPS measurements, globally, errors in meridional velocity is <0.02 mm/yr and in longitudinal direction is less than 0.03 mm/yr. For the interested Andes region, the error is slightly larger than global average (0.035 and 0.03 mm/yr). Over the Cascadia region, however, the errors are smaller than global average (0.027 and 0.015 mm/yr respectively). Based on the optimized solution, further experiments are performed on the sensitivity of groundwater perturbation and the dynamic evolution of the fault strength and the associated earthquake cycles. Although the governing thermos-dynamic equations of SEGMENT are very similar to those of the Advanced Solver for Problems in Earth’s ConvecTion (ASPECT, Kronbichler et al. (2012) ), it should be noted that the ASPECT model, due to its intended research purpose, has coarse resolution and cannot sufficiently represent surface topography. In fact, surface topography-induced creeping is a significant component for GPS sites in mountainous regions and on islands. This results in discrepancies in model top-level velocity and GPS measurements. Through skilled numerical setting, SEGMENT sidestepped this issue and significantly improved the velocity simulation of the near surface (<100 km depth) layers, where the stress build-up is critical for earthquakes (Figure 3(a)).

Figure 3. Present day near surface (a) and deep (b) level plate motion velocities (blue arrows, in m/yr). GPS measurements (black and red (for regions of interest only) arrows) are from ftp://data-out.unavco.org/pub/products/velocity/pbo.final_nam08.vel (for North American Plate) and http://sideshow.jpl.nasa.gov/post/series.html (for global coverage). In Panel (b), resonance frequencies of global tectonic plates at the interfaces are labeled (years in red font). Hollow red arrows indicate the rupture tendency for the upcoming major earthquakes that are temporally in sequence. This also is the relocating direction of their epicenters.

3.2.3. Groundwater Time Series

From Equations (1) and (2), groundwater, through affecting the gravity field, contribute to fault instability. Its spatial fluctuation, through enhancing granular material generation, also weakens the fault from within (Section 3.1.2). To obtain the groundwater fluctuations over the past 15 years (2000-2015), a continuous time marching of the land surface scheme is performed with 30-min time step. Atmospheric parameters such as near surface air temperature, humidity, pressure, winds, and radiative components are from NCEP/NCAR reanalyses, and precipitation information from the Global Historical Climatology Network (GHCN, 1900-2015), GPCP, and The Tropical Rainfall Measuring Mission (TRMM, 1998-2015) are taken as input during the observational periods to drive the land surface scheme in SEGMENT to diagnose how much precipitation is percolated into ground water reservoir. For the recent decade (2003 onward), GRACE measurements provide a reliable means for the groundwater recharge/discharge. Model simulations are satisfactory when verified against GRACE measurements. GRACE is of very coarse resolution and model can provide spatial distribution of groundwater at much finer resolution. Obtaining groundwater distribution in a warmer climate is challenging. A small ensemble of climate models (27 climate models listed in Table 1) and an innovative approach in estimating extreme precipitation (Ren, 2014a) is used for the prediction of future occurrence of earthquakes. The partition into surface runoff, evapotranspiration and ground water percolation also is performed within SEGMENT.

Table 1. Twenty seven GCM models used in this study to provide atmospheric forcing parameters to SEGMENT.

3.2.4. Granular Material Generation and Transportation

In analogous to the state equation in Scholz (1998) , the following force-restore relationship is proposed here for granular material generation:

$\frac{\partial a}{\partial t}=-c_{1}a+c_2\delta \Delta E+ADV$ (3)

where a is the thickness of the granularized rock layer, starting from 0, t is time, c₁ is the e-fold time scale, and c₂ is the Paris coefficient (Timoshenko & Gere, 1963) . Apparently, c₁ signifies the negative feedback as granular material accumulates. The unit surficial energy density of rocks, J (<300 mJ/m², Gilman (1960) ), is inverse proportional to c₂. To form new surfaces at high confining pressure (~0.1 GPa), the energy need far surpasses the surficial energy. The c₂ thus also is inversely proportional to confining pressure. Thermal effects also are included in c₂. The coefficient $\delta$ is 0 when the stress intensity change falls within the elastic range of the parent rock and is 1 when exceeding the elasticity range. $\Delta E$ is the energy input rate when stress perturbation caused imbalance in principal stress exceeds the elastic range. External (e.g. groundwater fluctuation) tides with resonance frequencies are most effective in transferring distortion energy into the plates to cause fatigue and influence the granular material production rates.

$\Delta E=b_1{\text{e}}^{b_0\left|\nu -{\nu }_p\right|}{\left(\Delta K\right)}^4$ (4)

where the Paris coefficient $b_1$ depends on material property such as unit surficial energy and porosity of fault rocks. The coefficient $b_0$ is a negative number indicating the exponential damping of the groundwater tidal wearing when it is not synchronized with the natural frequency of the crust plate, $\nu$ is the groundwater variation frequency, ${\nu }_p$ is the resonance frequency of the crust plate, and $\Delta K$ is the range of stress intensity change (proportional to the groundwater variations’ amplitude squared). From 2003 onward, GRACE measurements are available for improved SEGMENT estimates of the groundwater fluctuation amplitude. Plates at their interface assume complex configurations and determining the resonance frequency ${\nu }_p$ is achieved by solving the deflecting equation, with realistic plate configurations and distributed density and temperature from local geology maps and recent observational studies (e.g., the Andes Survey 2006). The estimated inherit resonant frequencies are labeled in Figure 4. The framework of Equations (3) and (4) also is applicable to the silica granules “floating” above the mantle wedge corner. In that circumstances, the $\Delta E$ is parameterized as thermal control: a thermal switch is designed as a function of the temperature gradients around the mantle wedge and the plate dipped in. In that manner, GM generation is shut-off for cold-slab subduction zones such as the Japan Trench and Hikurangi; whereas being active around warm-slab subduction zones such as the Mexico and Nankai subduction zones.

Figure 4. SEGMENT diagnosed time scales of plate resonance frequency to groundwater fluctuation forcing (numbers in red). Background map was obtained from sideshow.jpl.nasa.gov/mbh/series.html and has been cleared for public release.

The inter-plate granular material chute system, as left behind by seamounts, removes granular deposits through advection process (ADV term in Equation (3)) exactly the same way as the surface runoff moves debris (Dietrich & Perron 2006) . Transportation of granular material is left to the Navier-Stokes solver of SEGMENT. Although granular materials are ubiquitous between contacting plates, under the huge confining pressure, the mean free path, in analogous to Prandtl mixing length in fluid dynamics, of the granular particles are very limited for most places lacking the parallel granular stripes left behind by seamounts or strong-contact “asperities” (Figure 1(b)). In the absence of macroscopic tunnels between plates, granular material redistribution resembles the diffusion of heat: moving along the potential gradient (included in coefficient c₁). Earthquakes, primarily through p-wave shaking, increase the mean free path of the granular particle and make it more like fluids, and facilitate its downstream dispatch. Advection dominates the co-seismic down-dip transport and diffusion processes works steadily during the inter-seismic period. Earthquakes also are the mechanism shedding off the “saw dust” and facilitating generation of new granular material. Granular material’s formation is a fatigue/damage process (Bercovici & Ricard, 2014) . As cyclic variations in forcing favors fatigue (first two terms on the RHS of Equation (3)), especially when the forcing frequency is in resonance with the segment of the plate. In this sense, groundwater fluctuation, which possesses frequency on annual to multi-decadal scales, likely in concord with the fatigue of the inter-plate strong contacts. Thus, large spatial scale droughts and floods contribute to the irregularity of earthquake cycle.

Although in principle healing effects can be included into c₁, we did not involve the healing mechanism in Equation (3) because it is only applied in the brittle crust zone near subduction zones, where temperature are lower due to the double insulation from oceanic and continental plates. The best way to heal a fatigue, the phase changes (e.g., melt) are not relevant here. The released elastic energy during an earthquake is primarily dissipated as heat at the interface. The remainder propagates away as P and S-wave energy. Earthquake involves long distance displacements may cause melt and partially heal the fatigues interface (Rowe et al., 2012) . For the inter-seismic periods, Equation (3) is safely applicable. Clarifying these points about granular material also helps as numerical modeling of faults are now feasible with super-computers and more detailed observations (e.g., INSAR and GPS measurements of displacements; GRACE measurements of mass changes and numerous ways of thermal heat fluxes).

Transportation of granular material is left to the Navier-Stokes solver of SEGMENT. The motion of the granular material inside the chute system at the plate interface is primarily low Reynolds number (Re) laminar flow. The gravity potential associated with the self-weight of granular material is negligible compared with its dominance at the earth surface flow systems. Conversely, confining pressure, the granular viscosity, and the interaction of granular with the channel boundary become important. These forces cause a passive pumping of granular material in the inter-plate chute systems. The thermal process also is simplified by the laminar flow because convective mixing in regular free fluids are non-existent and only diffusion kinetics are dominant. In many aspects, it is similar to micro-fluids in the capillary pipes. Confining pressure becomes the driver in guiding the flow direction.

As in Figure 1(c), the strong contact regions have the most active granular material generation. Its existence in the locking zone assists the unlocking because the shear resistance granular material can provide at the strain rate of the locking stage is equivalent to a ${\mu }^{\prime }$ value of 10⁻³, four orders of magnitudes smaller than solid-solid contact, representing a lubricant between the two elastic plates. The bridging effects by themselves may not cause fault weakening. But with the associated granular material generation, they definitely weaken faults through forming a negative spatial covariance between loading and friction coefficients. During the past two decades, it is being gradually realized that the maximum apparent frictional coefficient decreases with time (Scholz, 1998) . Granular material formation is one feasible explanation for slide weakening. The generation and transport of granular material at the plate contact can reconcile many counter-intuitive characteristics of major earthquakes ( Wang, 2015 and references therein) and also is critical in organizing small-scale ruptures to generate major earthquakes. Proper treatment of present granular material (Figure 1(d)) between plates also may provide a time frame for the next major earthquake on the same fault. From Figure 1(d), it seems, globally, the granular material thickness distribution is positively correlated with major earthquakes’ occurrence frequency. Following sections are discussion of earthquakes over two selected regions, as actual applications of the numerical models.

4. Results

There are three interlaced mechanisms working together to make the earthquake prediction tricky for existing rubrics. The mechanisms we identified after sedulous essays using the singular SEGMENT model and remote sensing data are whetted against 73 historical cases, all within the reanalyses period with quality precipitation data. The reproduced earthquakes are viable, especially for cases after 2003, the starting point of the GRACE measurements, and before 2015, the aging of the current GRACE satellites. Basic model background verifications are available in Section 3, including simulated surface velocities against GPS measurements. Following sub-sections summarize model verification and simulations over three representative regions where groundwater fluctuations play significant role in affecting occurrence of major earthquakes. Although this study attempts a projection of future major earthquakes, the results cannot be regarded too literally. This is because major earthquakes are rare events and the existing information suitable for training the dynamical model still is scarce. Consequently, many ad-hoc assumptions are still involved in the parameterizations and interpretation of the input data, making the predictions of spatial and temporal of the future occurrences only of reference values.

4.1. Model Verification

By simulating the time scale of the evolution of strain and stress in the fault zone, using SEGMENT, with appropriate parameter setting, assisted by remote sensing information, the timing and location of major earthquakes (M_w > 5 globally) can be estimated, with uncertainties in timing reduced down to a year and location error within 120 km (root mean squared error for historical records during 1979-2015). The primary limitation on the prediction accuracy resides with the resolution of CRUST-1.0 model, depending on which SEGMENT resolves fault geometry and composing material properties. As with all prediction system, the uncertainty grows exponentially with forecasting time. In the following discussion, the uncertainty level will be presented as necessary.

Although global results are obtained, the following discussion focuses on three regions: The Tibetan Plateau and surrounding regions (TP), the Cascadia, and the Andes subduction faults. All three regions are adequately covered by GPS measurements and have large enough land areas that can be sensed by GRACE for the mass fluctuations around the faults. Primarily from different fault geometries, plate opposing motion speeds, and the distribution patterns of inter-plate granular material, earthquakes over the selected three regions are of various morphologies. Major earthquakes over the Andes and the TP refer respectively to those with magnitude > M_w7.5 and >M_w5.5. Without explicit statements, >M_w5 over other world regions is regarded as major earthquakes.

4.2. Earthquake Activity in the TP Region: Correlations with Groundwater Fluctuations

Large scale (e.g., thousand kilometers) groundwater fluctuations exert extra stress on fault interface. Almost all major (>M_w5.5) earthquakes during 2003-2016 over the Tibetan Plateau occurred during low groundwater phases (figure not shown). Although the fluctuation of groundwater is secondary to plate motion in causing stress build up, it is an “extra”, and usually “the last straw that breaks the camel”. Analyzing the stress field clearly indicates that once two plates are coupled together, the underlain plate drags the lighter overriding plate downward and they go in tandem deeper into the Earth, pushing away the higher density upper mantle material. This causes mass loss around the fault region. This process, although only causing mass loss at a rate two orders of magnitude smaller than groundwater discharge, is however aided when the region is experiencing mass loss due to groundwater reduction (e.g., when evaporation steadily exceeds precipitation during extended droughts). The recent Nepal earthquakes of April and early May, 2015, and a third earthquake in Burma a year later all resides inside the mass loss region. The temporal sequence of the three major earthquakes and the ruptures’ spatial pattern (starting from northwest and extending to southeast) can be explained by the granular material production and transportation condition as proposed here. They represent key stages of strain energy release in the lateral direction, along the transverse direction of the fault by unlocking the granular sticky spots scattered over the plate interface. The order of unlocking strictly followed the granular material accumulation on the shear zone.

The three epicenter-consisting patches (consisting the three respective epicenters) were all about to approach saturation by April 2015, after ~72-year’s stress accumulation (Bollinger et al., 2014; Men & Zhao, 2016) . The reason the west-most one ruptured first is because the aiding of groundwater reduction. Based on SEGMENT simulation, the groundwater recharge during the past 80 years has a decreasing trend over the Bay of Bengal region (74-104˚E; 15-35˚N as the region of interest), agree with the monsoonal precipitation trends as identified in Wang and Ding (2006) . On average, the groundwater reduction rate is ~31 Gt/yr. The super-imposed stress field produces ~85 km² saturated area (as defined in Equation (2)) annually. The matured area is not spatially contiguous on the fault interface. Rather they scatter sporadically over the inter-seismically locked interface. Connection of these sporadically scattered saturated patches is critical for generating earthquakes. For example, the occurrence of the first earthquake (April 25, 2015; M7.8; epicenter at 28.147˚N, 84.708˚E) redistributed the granular material due underneath and re-assumes a stronger rock-rock contact. This lessens the degree of saturation of the second major earthquake (occurred on May 12, 2015 of M7.3 magnitude and epicenter located 200 km to its southeast at 27.973˚N, 85.963˚E) region and it took another month for its occurrence. The earthquakes over Burma (M_w6.9 on April 12, 2016 at 23.13˚N, 94.9˚E and M_w6.8 on August 25, 2016 at 20.9˚N, 94.6˚E), further to the southeast, were not directly related to these two events. They were a result of the lightening of the Bay of Bengal region, to be further addressed later. Unlike some recent studies (e.g., Men & Zhao, 2016 ) that predicted the future earthquake (after the Nepal earthquake 2015) would be to the north-west sector, SEGMENT simulation indicates that the sector to the south-east also is becoming unstable. Actually the maturation pattern is a jumping between the eastern sector and the western sector, with cluster centers (30˚N, 84.4˚E) and (27.5˚N; 86.2˚E). Consequently, major earthquakes also occur alternatively (between these two regions in a bi-polar mode).

The Longmenshan fault zone (LFZ, 102-106˚E; 30-33˚N) is a vivid incarnation of the “non-unique solution” configuration as presented in Section 2. The LFZ defines the eastern margin of the Tibetan Plateau. Very recent synthetic approach (gravity survey, seismic reflection profile and earthquake focal mechanism data, and sediment analysis; personal communication with Z-X. Li, 2016) indicated that the LFZ is roughly composed of two segments of variant formations. The central-northern section, formed ~40 Myr ago (Jiang & Li, 2014; Richardson et al., 2008) , is a lithospheric through-cut fault zone that possesses low elastic strength but with strong dextral strike-slip motions. In contrast, the southern sector is a crustal thrust zone dominated by shallow-angle thrust motion of the TP over the moderately stiff South China Plate. This Mobius-twisted fault plane and contrasting behaviors along the LFZ, in addition to providing kinematically favorable conditions for instability, have also profound influence on groundwater reservoir. The M_w7.8 Wenchuan earthquake further lent evidence of the importance of gravity fields in fostering major earthquakes.

Examining the mass change fields averaged between 2003 and 2008, it is clear that Wenchuan resides at the center of a saddle region (Figure 3 in Ren et al. (2015) ). For earthquakes fostered along the same fault (but temporally consecutive), the tectonic plate motion caused compressing usually invokes very similar saturation patterns. The primary reason for the randomness of epicenters is due to more transient perturbations such as groundwater fluctuation. The super imposed field cause a specific point to reach saturation first and propagate away to cause a large rupture area. The rupture can grow because, to seize the motion of the unstable plates, neighboring unsaturated regions may be activated. If the propagation encountered too many unsaturated regions, the momentum is quickly lost and this strongly limits the magnitude of the earthquake. On the other hand, if the propagation starts from a region surrounded by a ripen neighborhood, a positive feedback forms and resulted into a large magnitude earthquake. By initiating the rupture near Wenchuan, a more efficient energy release is achieved. From SEGMENT simulation with future precipitation scenarios under A1B scenario provided by a small-volume climate model ensemble (as listed in Table 1), the next major earthquake will be within Ganshu province (epicenter located along 104.8˚E longitudinal arc between 33 and 35˚N) to the north east, ~120 ± 45 years later of magnitude ~M_w7. For the far-larger TP region, the India-Eurasin collisional system, while thickening the elastic lithosphere, created numerous vertically extended crevasses. These crevasses are very effective in channeling groundwater into depth. Combined with TP’s orographic effects, precipitation variation transfers readily into extra lateral stress on the faults. Therefore, groundwater fluctuation becomes a dynamic factor affecting the epicenter location and the magnitude of major earthquakes. Because the persistent gravity reduction in the Bay of Bengal area and Southwest China, coupled with insignificant changes over the plateau and the northern Tarim and very arid Qaidam basins, a saddle-shape super-imposed stress field is formed and maintained over the past half century. The natural, geological background-deduced recurrence frequency (“dry” frequency) is greatly reduced over the region, especially over the Southeast Asia peninsula. The fault zone over Burma likely becomes unstable in about two years. Very interestingly, the epicenters of the upcoming earthquakes, affected by the weakening of monsoonal precipitation, are located along 96 ± 1˚E parallel between 17-28˚N (most frequently 17-23˚N sector), roughly parallel the Ayeyarwady river basin.

Factors controlling earthquake irregularity over Cascadia and Andean share similarity with the TP region but have their respective specialties. Here provides a brief summary before detailed discussion. The granular material filled tunnels left behind by seamounts located on top of the Nazca and Juan Fernandez ridges foster many major Andean earthquakes. Because of the proximity of ENSO, the groundwater fluctuation over the Peru and Chilean coasts bear apparent 4 - 7 year fluctuations. However, unlike the collisional system of TP, there lack through-cut faults to channel precipitation into deeper depth (closer to the locked fault zone). From Section 3.1.1, the deeper the groundwater reservoir resides, the more efficient its mechanical effects. As a result, groundwater fluctuations have to exert effects mechanical-indirectly through aiding granular material generation and distribution. Consequently, earthquakes occur after a significant hiatus. Andes earthquakes are co-controlled by groundwater fluctuation and existing granular material distribution patterns.

Cascadia is similar to Andean subduction faults such as Nazca and Juan Fernandez in that they are weak inter-plate coupling and the natural occurrence frequency of earthquakes are rather low. Also, there lack secondary through-cut faults on the continental plates. Cascadia fault is unique in that local precipitation, although being of highest rate in North America, lacks annual to decadal variabilities. That explains the ~600 year natural major seismic cycle. To its south, along the pacific coast cordillera (ranges), lies the Mediterranean climate California, which precipitation is affected by teleconnection patterns such as ENSO and PDO and turned out to be sensitive to climate warming (Ren et al., 2011) . The Cascadia fault is thus more sensitive to Californian precipitation than local precipitation.

4.3. Andes Earthquakes: Co-Controlled by Granular Material Pattern and Groundwater Fluctuation

Figure 5 illustrates the present plates interface geometry at coastal Andes, as CRUST1.0 provided to SEGMENT. Panel (b) shows the ongoing base state creeping structure. The black arrows are horizontal flow speeds. The convergence in this direction is apparent but the total convergence is not apparent because the horizontal velocities rotate counterclockwise in the horizontal plane and the magnitudes are not changing so drastically. The coupling is weak as to the compressive stress between continental and oceanic plates. Especially, the ocean-continental subduction here does not lead to thickening of the continental plate. Baby et al. (1997) suggested that the subduction of oceanic lithosphere coupled with under-plating and a brief episode of gravity spreading contributed to crustal thickening in the back arc of the Central Andes. There still lacks a uniform theory that explains how very thick crusts develop (e.g., anomalies in the map of Christensen and Mooney (1995) ).

4.3.1. Footprints of Seamounts-Role of Granular Material in Andes Earthquakes Spatial Pattern

Granular accumulation along the Peru coast (South of Lima) extends deep inland (Figure 1(d)). The pattern reflects the way the Nazca Ridge passed through the continental plate. Because the thrust is at an acute angle to the continental plate motion, the channels on the overriding plate left by the seamounts is oriented at an angle to the Nazca Ridge’s moving direction (as shown on Figure 6).

Figure 5. Present plates interface geometry at coastal Andes (represented in the SEGMENT model using CRUST 1.0 data). Panel (a) is crust thickness distribution. Panel (b) is a vertical cross-section along the red line in panel (a). The black arrows are horizontal flow speeds. The convergence in this direction is apparent but the total convergence is not apparent because the horizontal velocities rotate counterclockwise and the magnitudes are not changing so drastically. Panel (c) is the density structure in the same vertical cross-section (as in (b)). Over this fault, the motion is driven by mantle flow from underneath. The coupling is weak as to the compressive stress between continental and oceanic plates. Especially, we would like to emphasize that ocean-continental subduction here does not lead to thickening of the continental plate. Geometry data obtained from CRUST1.0 (Laske et al., 2013) .

Figure 6. Granular material (shades are thickness in meter) around the subducting Nazca Ridge under the advancing South American plate. Currently locked regions (A, B to the left) with thick granular material accumulation are connected through macroscopic tunnels with seamounts (A, B and C on the right, over the Nazca Ridge). These seamounts grind the overriding continental plate and created the slanted tunnels. These tunnels collapsed and had earth surface footprints. Bold blue arrow is the moving direction of the flat oceanic slab with sea mounts A-C sitting on. The Hollow red arrow indicates the direction of the “tunnels” on the overriding continental plate plowed by seamounts. Note that the M_w8.1 Chilean 2014 event occurred exactly on the granular belt left by the “handle” of the Juan Fernandez Ridge.

This is just like chiseling a rotating disk. The granular material leftover also is distributed in stripes oriented the same way. In fact, all stripes originated from the sea mounts extends to the left flanks, due to the relative motion of the underlying and overriding plates. The granular material to the west (A’) ages as old as 14 Ma, whereas granular material at A (on the current flat lab) is newly generated. As the channels get older, they collapse and get wider but shallower, following the mechanisms as described in subsections 1.1 and 2.4 of the SM. These granular stripes become the soft belly of the fault and are the locales of many major earthquakes. Juan Fernandez to the south is very similar to Nazca ridge at present but are very different 3 Ma ago (Yanez et al., 2001; Yuan et al., 2002) . This chain of seamounts plowed beneath Puna at ~10-6 Ma ago like a hockey stick. At present, the “head” of the stick had melt into the upper mantle and only the “handle” left. However, the granular material grinded off by the “head and shaft” of the stick still acts as locales of earthquakes, especially at the seismic section of the subducting fault at the joint of Peru and Chile (at 20S; 70W). Also, the channels and the granular stripes left behind by the seamounts on the Juan Fernandez Ridge are more south-north oriented than those left behind by the Nazca Ridge. The epicenters of major earthquakes after 1800 correspond very well with the granular stripes. The granular stripes left by Juan Fernandez Ridge foster larger (but less frequent) earthquakes. For both seamount ridges, their left (northern) flanks foster larger earthquakes, because the granular accumulation patterns.

4.3.2. Groundwater’s Indirect Control (With Apparent Hiatus) of Andes Earthquakes’ Spatio-Temporal Patterns

For the Peru and Chilean coast, El Nino and the opposite phase La Nina ocean-atmospheric teleconnection patterns (ENSO) play a critical role in determining the phase of groundwater. Earthquakes, especially those greater than M_w 5.5 but less than M_w7, do not necessarily in phase with the ENSO. This means that, unlike for TP region, the direct mechanical effect, such as indicated by Equation (1), may not be the dominant factor. However, the occurrence frequency of the major earthquakes is at ~4 - 7 years, same as ENSO. Model simulation confirmed that it is the fluctuation, rather than the actual weight of groundwater that is critical for the granular material generation and ensuing earthquakes. In other words, for this region, groundwater caused loading fluctuation affects earthquakes through granular material generation and evolution, rather than through the direct mechanical mechanism (Equation (2)).

Even there are apparent temporal hiatus for most earthquakes between M_W5.5 and M_W7, all major earthquakes greater than M_w7.5 occurred at dry stages of the groundwater fluctuation cycles. For example, the 1998, 2003 and 2010 earthquakes all occurred at the nadir of the regional gravity field. In a certain sense, ENSO lessened the severity of the earthquakes over the Andes. Without ENSO, all the scattered < M_w7.5 earthquakes will occur as less frequent but more severe 1960-like super earthquakes. Through affecting precipitation partition into evaporation, runoff and infiltration into groundwater reservoir, topographic features have strong correlation of seismogenic zone segmentation (Bejar-Pizarro et al., 2013) . However, it simply is an illusion that surface topographic feature corresponds to ruptures at depth. Instead, the precipitation and the ensuing groundwater distribution contribute to earthquake irregularity. This is primarily because the input of strain energy from the periodic recharge and discharge of groundwater in a large regional area increases the prospect of locking formation, through granular fatigue. Monitoring the fluctuation tides of groundwater thus is a practical way for understanding earthquake cycle. Gravity mapping satellites such as GRACE provide essential assistance.

4.4. Cascadia Subduction Zone’s Stability Affected by Californian Droughts

Detailed discussion of Cascadia subduction fault, such as effective strength ( ${\mu }^{\prime }$ ), geometrical characteristics and unique geological conditions can be found in Wang and He (1999) . Here we would like to discuss the extra stress field super-imposed by precipitation fluctuations. Californian droughts (2011-2015), even though hundreds kilometers apart, significantly enhanced the instability of the Cascadia subduction zone. Due to the lack of vertically extended crevasses such as resulted from the India-Eurasia collisional system, groundwater fluctuations take effect through aiding granular material formation. Therefore there is a hiatus between droughts and the seismogenic effects. Rather than a direct temporal correspondence, there is a multiple year hiatus before the seismogenic effects from droughts to set on. For example, the extended Californian drought may contribute to fault instability 2 (southern side Oregon-Washington coast) to 10 years (northern BC coast) later. In addition, the Cascadia fault’s locking is special in that the oceanic plate still is creeping, which is clear from examining the velocity profile. It is locked in a sense that the neighboring regions have large gradients in motion speeds, signifying stress build-up. For example, the Olympic National Park coast of the Cascadia subduction zone is now inter-seismically locked (Figure 5). At the contact, the speed of the continental plate is close to zero. The oceanic plate, however, still is in motion. The continental plate is accumulating elastic potential energy at an annual rate of 7.6 × 10¹³ J/yr. The estimated apparent frictional coefficient is 0.12. If there is no perturbation on the loading field, the natural frequency of seismic cycle is ~600 year. By releasing elastic energy accumulated, an earthquake of M_w7.8 will be produced (of energy of ~45.7 PJ). However, from SEGMENT simulations, it is apparent that even the Californian drought exerted a cyclic stress on the interface and enhances the fatigue of the interface. If the same hydrological cycle over the past 50 years are representative for the past 800 years, the recurrence frequency is reduced to ~446 ± 70 years. As a result, a >M_w7.5 (~34 PJ) earthquake is expected within the upcoming decade. From the granular material distribution pattern, the rupture likely starts from the southern sector, or between the Oregon coast and the BC, Canada coast.

The logic of relating California droughts to future Cascadia earthquakes seems counter-intuitive since Cascadia has a lot of its own hydrological processes (it enjoys some of the highest precipitation rates in North America). The reason the fault is more sensitive to Californian precipitation is because of the fluctuation in precipitation, not the absolute amount. The high latitude precipitation over Cascadia region of Canada does not have large inter-annual to inter-decadal variation as in the Californian region, which is dually influenced by ENSO and PDO. More importantly, the Mediterranean climate is sensitive to climate warming (Ren et al., 2011) . Local Canadian precipitation fluctuation would be more direct but the magnitude of fluctuations truly is smaller than that from Californian precipitation. The North American plate cannot keep full rigidity over such a distance. The waveguide for propagating the bridging effect zig-zags ~5.5 wave lengths roughly along the pacific coast cordillera (ranges).

5. Discussion

The quest to understand Earthquakes is an ancient and universal question. The new earthquake theory and its precepts presented here are a convergence of several factors: recent advances in remote sensing technology; recent, tectonically active stage motivated many researchers to re-evaluate previous assumptions (Wang, 2015) , sometimes to a radical degree, and most critically, a well-vetted landslide model as prototype of an advanced numerical modeling system to handle the mechanics of inter-plate earthquakes and to verify the original hypotheses proposed here.

Earthquake occurs when the plate-coupling stress determined repose angle reaches the sum of two angles representing respectively fault geometry and fault strength. Accurate representation of the evolution of fault strength and repose angle thus is critical for earthquake simulation and prediction. The earthquake triad is a generic framework that applies unanimously. Groundwater fluctuations, although being small, play a pivotal role in connecting the three components together in regulating earthquake irregularity. Groundwater may influence fault instability mechanical-directly by super-imposing a seismogenic lateral stress field that works in synergy with the plate coupling stress; and mechanical-indirectly, or weakens the fault from within, achieved by enhancing fault fatigue through granular material generation. Mechanical-directly alone, groundwater fluctuation superimposed stress field is on the same order as the residual of gravity-aided frictional stress and plate coupling stress thus is non-negligible in affecting the timing of major earthquakes.

Closely centered on the earthquake triad, a minimum requirement for a working earthquake prediction system for realistic 3D fault configuration, the Groundwater aided Granular material Generation framework is proposed here and implemented in the well-vetted SEGMENT numerical modeling system. By exploring existing records of earthquakes using SEGMENT, it is found that for collisional systems, mechanical-direct effects from groundwater control the irregularity of major earthquakes. This is because the orographically spatially highly biased precipitation is effectively channeled into deeper depth by the prevalence of through-cut faults. Droughts elsewhere also are seismogenic but likely take effects mechanical-indirectly. For example, ENSO, as the dominant player for regional precipitation, has strong influence on the gravity field over the Andes region. The occurrence of major earthquakes, although bearing the same 4 - 7 years occurrence frequency of ENSO, has a significant hiatus when compared with the fluctuation in the gravity field. Similarly, the stability of the Cascadia fault is remotely affected by Californian droughts. The input of strain energy from the periodic recharge and discharge of groundwater in a large regional area increases the prospect of un-locking of seismically coupled plates. In a certain sense, climate fluctuations, through extended droughts determining groundwater loading pattern, affect the earthquake cycles.

This study is a complement to the prevailing Coulomb fracture criterion and Anderson’s theory of faulting (Jaeger, 1963; Fossen, 2010) . The earthquake triad concept reconciles and rectifies some existing conceptual errors in the field. All previous earthquake models considered one or at most two aspects of the triad hence lacked predictive capability. Solution non-uniqueness essentially roots in information paucity. Remotely sensed data found their use in deducing interpolate granular material distribution as well as verifying model simulated precipitation drainage into groundwater reservoir, greatly improved model performance. Although we still cannot make punctilious predictions in timing and location of epicenter, the performance of the model is trenchant enough and deserves further improvements, especially as we are entering a big data era.

Acknowledgements

We are thankful to Drs Z-X Li and K. Wang for discussion on earthquake mechanisms and various modern approaches for monitoring stress build up. Dr Z-X Li also is appreciated for providing the Andes Geological Survey details. Drs Lance M Leslie, M. Lynch and J Bornman proof-read the first draft of this manuscript. The data used are all publicly available.

Additional Information

The author declares no competing financial interests.

Appendix 1: Velocity Transformation between SEGMENT’s Sphere Crown Model and Cartesian System of ASPECT

For generic curvilinear system, we have the following identities.

$\text{d}\Re =H_i{e}_i\text{d}{q}_i$ (A1)

$V=\text{d}\Re /\text{d}t=H_i{e}_i\text{d}{q}_i/\text{d}t=V_{qi}{e}_i$ (A2)

$V=\text{d}\Re /\text{d}t=u{i}_0+v{j}_0+w{k}_0$ (A3)

Using (A2),

$\{\begin{array}{l}u=\text{d}x/\text{d}t=\frac{\partial x}{\partial q_1}\text{d}{q}_1/\text{d}t+\frac{\partial x}{\partial q_2}\text{d}{q}_2/\text{d}t+\frac{\partial x}{\partial q_3}\text{d}{q}_3/\text{d}t=\frac{\partial x}{\partial q_i}\frac{V_{qi}}{H_i}\\ v=\text{d}y/\text{d}t=\frac{\partial y}{\partial q_1}\text{d}{q}_1/\text{d}t+\frac{\partial y}{\partial q_2}\text{d}{q}_2/\text{d}t+\frac{\partial y}{\partial q_3}\text{d}{q}_3/\text{d}t=\frac{\partial y}{\partial q_i}\frac{V_{qi}}{H_i}\\ w=\text{d}z/\text{d}t=\frac{\partial z}{\partial q_1}\text{d}{q}_1/\text{d}t+\frac{\partial z}{\partial q_2}\text{d}{q}_2/\text{d}t+\frac{\partial z}{\partial q_3}\text{d}{q}_3/\text{d}t=\frac{\partial z}{\partial q_i}\frac{V_{qi}}{H_i}\end{array}$ (A4)

Applying to a spherical coordinate system (Figure A1) with the forward transformation equation of

Figure A1. Geometry representation of the spherical coordinate system as in SEGMENT. Panels (a)-(c) show respectively the projections of unit directional vector $i_0,j_0$ and $k_0$ into the spherical coordinates ${\theta }_0,{\lambda }_{\text{0}}$ and $r_0$ . The geometrical relationships in (d) is repeatedly used in deriving the cosine relationships of angles composed by the intersection line and two vectors residing respectively in the two mutually perpendicular planes. The projection relationship between the unit cosine directions are summarized in (a)-(c).

$\{\begin{array}{l}x=r\cos \lambda \cos \theta \\ y=r\cos \lambda \sin \theta \\ z=r\sin \lambda \end{array}$ (A5)

where in geoscience convention, $\theta$ is longitude, $\lambda$ is latitude and r is the distance from the earth center (Figure A1). Still in right-handed coordinate system, lets define $q_1=\theta ,q_2=\lambda ,q_3=r$ . The corresponding Lame operators:

$\{\begin{array}{l}{H}_{\theta }=\sqrt{{\left({x^{\prime }}_{\theta }\right)}^2+{\left({y^{\prime }}_{\theta }\right)}^2+{\left({z^{\prime }}_{\theta }\right)}^2}=r\cos \lambda \\ H_{\lambda }=\sqrt{{\left({x^{\prime }}_{\lambda }\right)}^2+{\left({y^{\prime }}_{\lambda }\right)}^2+{\left({z^{\prime }}_{\lambda }\right)}^2}=r\\ H_r=\sqrt{{\left({x^{\prime }}_r\right)}^2+{\left({y^{\prime }}_r\right)}^2+{\left({z^{\prime }}_r\right)}^2}=1\end{array}$ (A6)

From Equation (A2),

$\{\begin{array}{l}{V}_{\theta }=r\cos \lambda \frac{\text{d}\theta }{\text{d}t}\\ V_{\lambda }=r\frac{\text{d}\lambda }{\text{d}t}\\ V_r=\frac{\text{d}r}{\text{d}t}\end{array}$ (A7)

From Equation (A4), we get

$\{\begin{array}{l}u=\frac{\partial x}{\partial q_i}\frac{V_{qi}}{H_i}=\frac{\partial x}{\partial \theta }\frac{V_{\theta }}{H_{\theta }}+\frac{\partial x}{\partial \lambda }\frac{V_{\lambda }}{H_{\lambda }}+\frac{\partial x}{\partial r}\frac{V_r}{H_r}\\ =\left(-r\cos \lambda \sin \theta \right)\frac{V_{\theta }}{r\cos \lambda }+\left(-r\sin \lambda \cos \theta \right)\frac{V_{\lambda }}{r}+\left(\cos \lambda \cos \theta \right)\frac{V_r}{1}\\ =-\sin \theta V_{\theta }-\sin \lambda \cos \theta V_{\lambda }+\cos \lambda \cos \theta V_r\\ v=\frac{\partial y}{\partial q_i}\frac{V_{qi}}{H_i}=\frac{\partial y}{\partial \theta }\frac{V_{\theta }}{H_{\theta }}+\frac{\partial y}{\partial \lambda }\frac{V_{\lambda }}{H_{\lambda }}+\frac{\partial y}{\partial r}\frac{V_r}{H_r}\\ =\left(r\cos \lambda \cos \theta \right)\frac{V_{\theta }}{r\cos \lambda }+\left(-r\sin \lambda \sin \theta \right)\frac{V_{\lambda }}{r}+\left(\cos \lambda \sin \theta \right)\frac{V_r}{1}\\ =\cos \theta V_{\theta }-\sin \lambda \sin \theta V_{\lambda }+\cos \lambda \sin \theta V_r\\ w=\frac{\partial z}{\partial q_i}\frac{V_{qi}}{H_i}=\frac{\partial z}{\partial \theta }\frac{V_{\theta }}{H_{\theta }}+\frac{\partial z}{\partial \lambda }\frac{V_{\lambda }}{H_{\lambda }}+\frac{\partial z}{\partial r}\frac{V_r}{H_r}\\ =0\frac{V_{\theta }}{r\cos \lambda }+\left(r\cos \lambda \right)\frac{V_{\lambda }}{r}+\left(\sin \lambda \right)\frac{V_r}{1}\\ =\cos \lambda V_{\lambda }+\sin \lambda V_r\end{array}$ (A8)

Equation (A8) is forward transformation, or using velocity in spherical coordinates (SEGMENT) to express velocity components in Cartesian system.

Define matrix $\begin{array}{c}{I}_1=\left(\begin{array}{c}{i}_0\\ j_0\\ k_0\end{array}\right)\left({\theta }_0,{\lambda }_{\text{0}},r_0\right)=\left(\begin{array}{ccc}\cos \left(i_0,{\theta }_0\right)& \cos \left(i_0,{\lambda }_{\text{0}}\right)& \cos \left(i_0,r_0\right)\\ \cos \left(j_0,{\theta }_0\right)& \cos \left(j_0,{\lambda }_{\text{0}}\right)& \cos \left(j_0,r_0\right)\\ \cos \left(k_0,{\theta }_0\right)& \cos \left(k_0,{\lambda }_{\text{0}}\right)& \cos \left(k_0,r_0\right)\end{array}\right)\\ =\left(\begin{array}{ccc}-\sin \theta & -\cos \theta \sin \lambda & \cos \theta \cos \lambda \\ \cos \theta & -\sin \theta \sin \lambda & \sin \theta \cos \lambda \\ 0& \cos \lambda & \sin \lambda \end{array}\right),\end{array}$ then Equation (A8) can be expressed as

$\left(\begin{array}{c}u\\ v\\ w\end{array}\right)=I_1\left(\begin{array}{c}{V}_{\theta }\\ V_{\lambda }\\ V_r\end{array}\right)$ (A9).

Inverse matrix of I₁ is

$\begin{array}{c}{\left(I_1\right)}^{-1}=\left(\begin{array}{c}{\theta }_0\\ {\lambda }_{\text{0}}\\ r_0\end{array}\right)\left(i_0,j_0,k_0\right)=\left(\begin{array}{ccc}\cos \left({\theta }_0,i_0\right)& \cos \left({\theta }_0,j_0\right)& \cos \left({\theta }_0,k_0\right)\\ \cos \left({\lambda }_{\text{0}},i_0\right)& \cos \left({\lambda }_{\text{0}},j_0\right)& \cos \left({\lambda }_{\text{0}},k_0\right)\\ \cos \left(r_0,i_0\right)& \cos \left(r_0,j_0\right)& \cos \left(r_0,k_0\right)\end{array}\right)\\ =\left(\begin{array}{ccc}-\sin \theta & \cos \theta & 0\\ -\cos \theta \sin \lambda & -\sin \theta \sin \lambda & \cos \lambda \\ \cos \theta \cos \lambda & \sin \theta \cos \lambda & \sin \lambda \end{array}\right)\end{array}$ (A10).

Notice also that transpose of I₁ also is inverse of I₁ (I₁ is thus an orthogonal matrix). In addition, I₁ also is unitary and normal matrix (det(I₁) = I).

From Equation (A9), $\left(\begin{array}{c}{V}_{\theta }\\ V_{\lambda }\\ V_r\end{array}\right)={\left(I_1\right)}^{-1}\left(\begin{array}{c}u\\ v\\ w\end{array}\right)$ , that is

$\{\begin{array}{l}{V}_{\theta }=\left(-\sin \theta \right)u+\left(\cos \theta \right)v\\ V_{\lambda }=\left(-\cos \theta \sin \lambda \right)u-\left(\sin \theta \sin \lambda \right)v+\left(\cos \lambda \right)w\\ V_r=\left(\cos \theta \cos \lambda \right)u+\left(\sin \theta \cos \lambda \right)v+\left(\sin \lambda \right)w\end{array}$ (A11)

Note also that I1 and its inverse matrix also denote the projection of the unit vectors between the curvilinear and Cartesian coordinate system. For example, simply look at the geometry in Figure (A1), it is apparent that the projection of ${\lambda }_{\text{0}}$ in the $\left(i_0,j_0,k_0\right)$ directions is exactly the second row in the inverse matrix (Equation (A10)). Similarly, projection of $j_0$ in the curvilinear coordinate system is the second row in I₁, or, $j_0=\cos \theta {\theta }_0-\sin \lambda \sin \theta {\lambda }_{\text{0}}+\cos \lambda \sin \theta r_0$ . Notice the similarity with Equation (A8). Actually, plugging these unit vectors’ projection expressions into Equations (A2) and (A3), the velocity transformation equations can be similarly obtained. For example,

$\begin{array}{c}v=V\cdot j_0=\left(V_{\theta }{\theta }_0+V_{\lambda }{\lambda }_{\text{0}}+V_r{r}_0\right)\left(\cos \theta {\theta }_0-\sin \lambda \sin \theta {\lambda }_{\text{0}}+\cos \lambda \sin \theta r_0\right)\\ =V_{\theta }\cos \theta -\sin \lambda \sin \theta V_{\lambda }+\cos \lambda \sin \theta V_r.\end{array}$ This is

exactly the second equation of Equation (A8). Similarly, using the direction-cosine relationships, for example, using the projections of ${\theta }_0$ in the original Cartesian coordinate system:

$\cos \left({\theta }_0,i_0\right)=-\sin \theta ;\cos \left({\theta }_0,j_0\right)=\cos \theta ;\cos \left({\theta }_0,k_0\right)=0$ . It is straightforward to obtain that

$V_{\theta }=V\cdot {\theta }_0=\left(u{i}_0+v{j}_0+w{k}_0\right)\left(-\sin \theta i_0+\cos \theta j_0+0k_0\right)=-\sin \theta u+\cos \theta v$ , exactly the first equation of Equation (A11).

Appendix 2: Viscoelastic Relaxation after Major Earthquakes

Rheology of mantle material is generally assumed to be viscoelastic. A combination of Maxwell fluid and Kelvin solid in the form of Burgers rheology is commonly accepted in the geophysics community.

In Figure A2, the Burgers rheology is represented by a serial connection of a Maxwell fluid of viscosity ${\eta }_M$ and rigidity ${\mu }_M$ and a Kelvin solid of viscosity ${\eta }_K$ and rigidity ${\mu }_K$ . ${\tau }_M$ and ${\tau }_K$ are Maxwell and Kelvin relaxation times.

1) A maxwell fluid is a viscoelastic material having the properties of both elasticity and viscosity. It can be represented by a purely viscous damper and a purely elastic spring connected in series. In this configuration, under applied axial stress, the total stress, ${\sigma }_{Total}$ and the total strain ${\epsilon }_{Total}$ can be defined as follows. ${\sigma }_{Total}={\sigma }_{D\left(amper\right)}={\sigma }_{S\left(pring\right)}$ ; ${\epsilon }_{Total}={\epsilon }_{D\left(amper\right)}+{\epsilon }_{S\left(pring\right)}$ . Where the subscript D indicates the strain/stress in the damper and the substract S indicates the stress/strain in the spring. For Newtonian fluids, we have

$\frac{\text{d}{\epsilon }_{Dl}}{\text{d}t}={\sigma }_D{\eta }_M={\sigma }_T{\eta }_M$ ; and from Hook’s law for elastic modulus form $\frac{\text{d}{\epsilon }_S}{\text{d}t}=\frac{1}{E}\frac{\text{d}{\sigma }_S}{\text{d}t}=\frac{1}{E}\frac{\text{d}{\sigma }_T}{\text{d}t}$ . The total strain rate thus follows

$\frac{\text{d}{\epsilon }_T}{\text{d}t}=\frac{\text{d}{\epsilon }_D}{\text{d}t}+\frac{\text{d}{\epsilon }_S}{\text{d}t}=\frac{{\sigma }_T}{\eta }+\frac{1}{E}\frac{\text{d}{\sigma }_T}{\text{d}t}.$ (A12)

In a Maxwell material, stress $\sigma$ , strain $\epsilon$ , and their rates of change w.r.t. time are governed by equations of the form $\stackrel{\dot{}}{\epsilon }=\frac{\sigma }{\eta }+\frac{\stackrel{\dot{}}{\sigma }}{E}.$

Equation (A12) can be applied either to the shear stress or to the uniform tension in a material. In the former case, the viscous component sorresponds to that for a Newtonian fluid. In the latter case, it has a slightly different meaning in relating stress and rate of strain. The Maxwell model usually is applied to the case of small deformations. For large deformations one should include some geometrical non-linearity. The following we discuss two special scenarios: a) effect of a sudden deformation, or the stress evolution after a sudden deformation; and b) Effect of a sudden stress, or how the strain evolves upon a sudden applied stress in the axial direction.

a) Effect of a sudden deformation

If an M material is suddenly deformed and held to a strain of ${\epsilon }_0$ , then the stress decays with a characteristic time of $\frac{\eta }{E}$ .

A brief proof of this starts from the governing equation (Equation A12) with initial condition of $\epsilon \left(t=0\right)={\epsilon }_0$ and strain rate term $\stackrel{\dot{}}{\epsilon }=0$ . That is

$\{\begin{array}{l}0=\frac{\sigma }{\eta }+\frac{\stackrel{\dot{}}{\sigma }}{E}\\ {\epsilon |}_{t=0}={\epsilon }_0,\sigma =\sigma \left({\epsilon }_{\text{0}}\right)={\sigma }_0;{\epsilon }_0={\sigma }_0/E\end{array}$ (A13).

Using the generic solution for the first order ordinary differential equation (Appendix 3), it is straightforward to obtain the solution for Equation (A13).

$\sigma \left(t\right)={\sigma }_0{\text{e}}^{-\frac{E}{\eta }t}$ (A14)

Figure A2. Burgers fluid as a serial connection of a Maxwell fluid (subscripts M) and a Kelvin solid (subscripts K). The Maxwell fluid and Kelvin solid are respectively serial and parallel connection of pure spring and pure viscous damper filled with Newtonian fluids.

If we then free the material at time t₁, then the elastic element will sprong back by the value of

${\epsilon }_{back}=-\frac{\sigma \left(t_1\right)}{E}={\epsilon }_0{\text{e}}^{-\frac{E}{\eta }{t}_1}$ (A15)

This is because at time t₁, from Equation (A14), keeping the existing deformation need a confining pressure of $\sigma \left(t_1\right)={\sigma }_0{\text{e}}^{-\frac{E}{\eta }{t}_1}={\sigma }_0\exp \left(-\frac{E}{\eta }{t}_1\right)$ . If remove the confining pressure at this moment, the elastic recovery would be $\frac{\sigma \left(t_1\right)}{E}=\frac{1}{E}{\sigma }_0\exp \left(-\frac{E}{\eta }{t}_1\right)\stackrel{\text{usingEq}\text{.}\left(\text{A13}\right)’\text{sinitialconditionexpression}}{=}{\epsilon }_0\exp \left(-\frac{E}{\eta }{t}_1\right)$ . As the viscous element would not return to its original length, the irreversible component of deformation can be simplified as

${\epsilon }_{irreversible}={\epsilon }_0-{\epsilon }_{back}={\epsilon }_0\left(1-\exp \left(-\frac{E}{\eta }{t}_1\right)\right)$ (A16)

b) Effect of a sudden stress

If an M material is suddenly subject to a stress ${\sigma }_0$ , then the elastic element would suddenly deform and the viscous element would deform with a constant rate

$\epsilon \left(t\right)=\underset{reversible}{\frac{{\sigma }_0}{E}}+\underset{irreversible}{\frac{{\sigma }_0}{\eta }t}$ . (A17)

(Recall the definition of Newtonian fluid as $\stackrel{\dot{}}{\epsilon }=\frac{\sigma }{\eta }$ , part of Equation (A12)).

If at some time t₁ we would release the material, then the deformation of the elastic element would be the spring back deformation and the deformation to the viscous part would not change back. That is at any arbitrary time t₁,

${\epsilon }_{reversible}=\frac{{\sigma }_0}{E};{\epsilon }_{irreversible}=\frac{{\sigma }_0}{\eta }{t}_1$ .

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Andes Geological Survey (2006).
[2]	Audet, P., & Burgmann, R. (2014). Possible Control of Subduction Zone Slow Earthquake Periodicity by Silica Enrichment. Nature, 510, 389-392. https://doi.org/10.1038/nature13391
[3]	Baby, P., Rochat, P., Mascle, G., & Hérail, G. (1997). Neogene Shortening Contribution to Crustal Thickening in the Back Arc of the Central Andes. Geology, 25, 883-886. https://doi.org/10.1130/0091-7613(1997)025<0883:NSCTCT>2.3.CO;2
[4]	Bejar-Pizarro, M., Socquet, A., Armijo, R., Carrizo, D., & Genrich, J. (2013). Andean Structural Control on Interseismic Coupling in the North Chile Subduction Zone. Nature Geoscience, 6, 462-467. https://doi.org/10.1038/ngeo1802
[5]	Bercovici, D., & Ricard, Y. (2014). Plate Tectonics, Damage and Inheritance. Nature, 508, 513-516. https://doi.org/10.1038/nature13072
[6]	Bollinger, L., Sapkota, S., Tapponnier, P., Klinger, Y., Rizza, M., Van der Woerd, J., Tiwari, D., Pandey, R., Bitri, A., & Bes de Berc, S. (2014). Estimating the Return Times of Great Himalayan Earthquakes in Eastern Nepal: Evidence from the Patu and Bardibas Strands of the Main Frontal Thrust. Journal of Geophysical Research: Solid Earth, 119, 7123-7163. https://doi.org/10.1002/2014JB010970
[7]	Cavalie, O., & Jónsson, S. (2014). Block-Like Plate Movements in Eastern Anatolia Observed by InSAR. Geophysical Research Letters, 41, 26-31. https://doi.org/10.1002/2013GL058170
[8]	Chlieh, M., Perfettini, H., Tavera, H., Avouac, J. P., Remy, D., Nocquet, J. M., Rolandone, F., Bondoux, F., Gabalda, G., & Bonvalot, S. (2011). Interseismic Coupling and Seismic Potential along the Central Andes Subduction Zone. Journal of Geophysical Research: Solid Earth, 116, B12405. https://doi.org/10.1029/2010JB008166
[9]	Christensen, N., & Mooney, W. (1995). Seismic Velocity Structure and Composition of the Continental Crust: A Global View. Journal of Geophysical Research: Solid Earth, 100, 9761-9788. https://doi.org/10.1029/95JB00259
[10]	Cowan, D., Cladouhos, T., & Morgan, J. (2003). Structural Geology and Kinematic History of Rocks Formed along Low-Angle Faults, Death Valley, California. Geological Society of America Bulletin, 115, 1230-1248. https://doi.org/10.1130/B25245.1
[11]	Dietrich, W., & Perron, J. (2006). The Search for a Topographic Signature of Life. Nature, 439, 411-418. https://doi.org/10.1038/nature04452
[12]	Famiglietti, J., & Rodell, M. (2013). Water in the Balance. Science, 340, 1300-1301. https://doi.org/10.1126/science.1236460
[13]	Fossen, H. (2010). Structural Geology (463 p). Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511777806
[14]	Gao, X., & Wang, K. (2014). Strength of Stick-Slip and Creeping Subduction Megathrusts from Heat Flow Observations. Science, 345, 1038-1041. https://doi.org/10.1126/science.1255487
[15]	Gao, X., & Wang, K. (2017). Rheological Separation of the Megathrust Seismogenic Zone and Episodic Tremor and Slip. Nature, 543, 416-419. https://doi.org/10.1038/nature21389
[16]	Gilman, J. (1960). Direct Measurements of the Surface Energy of Crystals. Journal of Applied Physics, 31, 2208-2218. https://doi.org/10.1063/1.1735524
[17]	Hayes, G., Wald, D., & Johnson, R. (2012). Slab 1.0: A Three-Dimensional Model of Global Subduction Zone Geometries. Journal of Geophysical Research: Solid Earth, 117, B01302. https://doi.org/10.1029/2011JB008524
[18]	Hyndman, R., McCrory, P., Wech, A., Kao, H., & Ague, J. (2015). Cascadia Subducting Plate Fluids Channeled to Fore-Arc Mantle Corner: ETS and Silica Deposition. Journal of Geophysical Research: Solid Earth, 120, 4344-4358. https://doi.org/10.1002/2015JB011920
[19]	Jaeger, J. (1963). Extension Failures in Rocks Subject to Fluid Pressure. Journal of Geophysical Research, 68, 6066-6067. https://doi.org/10.1029/JZ068i021p06066
[20]	Jiang, X., & Li, Z. (2014). Seismic Reflection Data Support Episodic and Simultaneous Growth of the Tibetan Plateau Since 25 Myr. Nature Communications, 5, 5453. https://doi.org/10.1038/ncomms6453
[21]	Jop, P., Forterre, Y., & Pouliquen, O. (2006). A Constitutive Law for Dense Granular Flows. Nature, 441, 727-730. https://doi.org/10.1038/nature04801
[22]	Kronbichler, M., Heister, T., & Bangerth, W. (2012). High Accuracy Mantle Convection Simulation through Modern Numerical Methods. Geophysics Journal International, 191, 12-29. https://doi.org/10.1111/j.1365-246X.2012.05609.x
[23]	Laske, G., Masters, G., Ma, Z., & Pasyanos, M. (2013). Update on CRUST1.0—A 1-Degree Global Model of Earth’s Crust.
[24]	Lay, T. (2015). The Surge of Great Earthquakes from 2004 to 2014. Earth and Planetary Science Letters, 409, 133-146. https://doi.org/10.1016/j.epsl.2014.10.047
[25]	Liu, X., Wang, Y., & Yan, S. (2017). Ground Deformation Associated with Exploitation of Deep Groundwater in Cangzhou City Measured by Multi-Sensor Synthetic Aperture Radar Images. Environmental Earth Sciences, 76, 6. https://doi.org/10.1007/s12665-016-6311-0
[26]	Liu, Z., & Bird, P. (2002). Finite Element Modeling of Neotectonics in New Zealand. Journal of Geophysical Research: Solid Earth, 107, ETG 1-1-ETG 1-18. https://doi.org/10.1029/2001JB001075
[27]	Loveless, J., & Meade, B. (2010). Geodetic Imaging of Plate Motions, Slip Rates, and Partitioning of Deformation in Japan. Journal of Geophysical Research: Solid Earth, 115, B09301. https://doi.org/10.1029/2008JB006248
[28]	Men, K., & Zhao, K. (2016). The 2015 Nepal M8.1 Earthquake and the Prediction for M ≥ 8 Earthquakes in West China. Natural Hazards, 82, 1767-1777. https://doi.org/10.1007/s11069-016-2268-2
[29]	Pritchard, M., & Simons, M. (2006). An Aseismic Slip Pulse in Northern Chile and along-Strike Variations in Seismogenic Behavior. Journal of Geophysical Research: Solid Earth, 111, B08405. https://doi.org/10.1029/2006JB004258
[30]	Ren, D. (2014a). Storm-Triggered Landslides in Warmer Climates. New York: Springer.
[31]	Ren, D. (2014b). The Devastating Zhouqu Storm-Triggered Debris Flow of August 2010: Likely Causes and Possible Trends in a Future Warming Climate. Journal of Geophysical Research: Atmospheres, 119, 3643-3662. https://doi.org/10.1002/2013JD020881
[32]	Ren, D., & Leslie, L. (2014). Effects of Waves on Tabular Ice-Shelf Calving. Earth Interactions, 18, 1-28. https://doi.org/10.1175/EI-D-14-0005.1
[33]	Ren, D., Fu, R., Leslie, L. M., & Dickinson, R. (2011). Predicting Storm-Triggered Landslides. Bulletin of the American Meteorological Society, 92, 129-139. https://doi.org/10.1175/2010BAMS3017.1
[34]	Ren, D., Leslie, L. M., & Karoly, D. (2008). Landslide Risk Analysis Using a New Constitutive Relationship for Granular Flow. Earth Interactions, 12, 1-16. https://doi.org/10.1175/2007EI237.1
[35]	Ren, D., Leslie, L. M., & Lynch, M. (2012). Verification of Model Simulated Mass Balance, Flow Fields and Tabular Calving Events of the Antarctic Ice Sheet against Remotely Sensed Observations. Climate Dynamics, 40, 2617-2636. https://doi.org/10.1007/s00382-012-1464-3
[36]	Ren, D., Leslie, L., Shen, X., Hong, Y., Duan, Q., Mahmood, R., Li, Y., Huang, G., Guo, W., & Lynch, M. (2015). The Gravity Environment of Zhouqu Debris Flow of August 2010 and Its Implication for Future Recurrence. International Journal of Geosciences, 6, 317-325. https://doi.org/10.4236/ijg.2015.64025
[37]	Richardson, N., Densmore, A., Seward, D., Fowler, A., Wipf, M., Ellis, M., Yong, L., & Zhang, Y. (2008). Extraordinary Denudation in the Sichuan Basin: Insights from Low-Temperature Thermochronology Adjacent to the Eastern Margin of the Tibetan Plateau. Journal of Geophysical Research: Solid Earth, 113, B04409. https://doi.org/10.1029/2006JB004739
[38]	Rowe, C., Kirkpatrick, J., & Brodsky, E. (2012). Fault Rock Injections Record Paleo-Earthquakes. Earth and Planetary Science Letters, 335-336, 154-166.
[39]	Sawai, Y., Satake, K., Kamataki, T., Nasu, H., Shishikura, M., Atwater, B., Horton, B., Kelsey, H., Nagumo, T., & Yamaguchi, M. (2004). Transient Uplift after a 17th-Century Earthquake along the Kuril Subduction Zone. Science, 306, 1918-1920. https://doi.org/10.1126/science.1104895
[40]	Scholz, C. (1998). Earthquakes and Friction Laws. Nature, 391, 37-42. https://doi.org/10.1038/34097
[41]	Sun, T., Wang, K., Iinuma, T., Hino, R., He, J., Fujimoto, H., Kido, M., Osada, Y., Miura, S., Ohta, Y., & Hu, Y. (2014). Prevalence of Viscoelastic Relaxation after the 2011 Tohoku-Oki Earthquake. Nature, 514, 84-87. https://doi.org/10.1038/nature13778
[42]	Timoshenko, S., & Gere, J. (1963). Theory of Elastic Stability. New York: McGraw-Hill.
[43]	Ujjie, K., Yamaguchi, A., Kimura, G., & Toh, S. (2007). Fluidization of Granular Material in a Subduction Thrust at Seismogenic Depth. Earth and Planetary Science Letters, 259, 307-318. https://doi.org/10.1016/j.epsl.2007.04.049
[44]	Van der Veen, C., & Whillans, I. (1989). Force Budget: I. Theory and Numerical Methods. Journal of Glaciology, 35, 53-60. https://doi.org/10.3189/002214389793701581
[45]	Voss, K., Famiglietti, J., Lo, M., de Linage, C., Rodell, M., & Swenson, S. (2013). Groundwater Depletion in the Middle East from GRACE with Implications for Transboundary Water Management in the Tigris-Euphrates-Western Iran Region. Water Resources Research, 49, 904-914. https://doi.org/10.1002/wrcr.20078
[46]	Wallace, L., Beavan, J., & McCaffrey, R. (2004). Subduction Zone Coupling and Tectonic Block Rotation in the North Island, New Zealand. Journal of Geophysical Research: Solid Earth, 109. https://doi.org/10.1029/2004JB003241
[47]	Wallace, L., Ellis, S., & Mann, P. (2009). Collisional Model for Rapid Fore-Arc Block Rotations, Arc Curvature, and Episodic Back-Arc Rifting in Subduction Settings. Geochemistry, Geophysics, Geosystems, 10, Q05001. https://doi.org/10.1029/2008GC002220
[48]	Wallace, L., Kaneko, Y., Hreinsdottir, S., Hamling, I., Peng, Z., Bartlow, N., D’Anastasio, E., & Fry, B. (2017). Large-Scale Dynamic Triggering of Shallow Slow Slip Enhanced by Overlying Sedimentary Wedge. Nature Geoscience, 10, 765-770. https://doi.org/10.1038/ngeo3021
[49]	Wang, B., & Ding, Q. (2006). Changes in Global Monsoon Precipitation over the Past 56 Years. Geophysical Research Letters, 33, L06711. https://doi.org/10.1029/2005GL025347
[50]	Wang, K. (2015). Subduction Faults as We See Them in the 21st Century. AGU Birch Lecture.
[51]	Wang, K., & Bilek, S. (2011). Do Subducting Seamounts Generate or Stop Large Earthquakes. Geology, 39, 819-822. https://doi.org/10.1130/G31856.1
[52]	Wang, K., & He, J. (1999). Mechanics of Low-Stress Forearcs: Nankai and Cascadia. Journal of Geophysical Research: Solid Earth, 104, 15191-15205. https://doi.org/10.1029/1999JB900103
[53]	Wang, K., & Hu, Y. (2006). Accretionary Prisms in Subduction Earthquake Cycles: The Theory of Dynamic Coulomb Wedge. Journal of Geophysical Research: Solid Earth, 111, B06410. https://doi.org/10.1029/2005JB004094
[54]	Wang, K., Hu, Y., & He, J. (2012). Deformation Cycles of Subduction Earthquakes in a Viscoelastic Earth. Nature, 484, 327-332. https://doi.org/10.1038/nature11032
[55]	Weston, J., Ferreira, A., & Funning, G. (2012). Systematic Comparisons of Earthquake Source Models Determined Using In SAR and Seismic Data. Tectonophysics, 532-535, 61-81. https://doi.org/10.1016/j.tecto.2012.02.001
[56]	Williams, C., Eberhart-Phillips, D., Bannister, S., Barker, D., Henrys, S., Reyners, M., & Sutherland, R. (2013). Revised Interface Geometry for the Hikurangi Subduction Zone, New Zealand. Seismological Research Letters, 84, 1066-1073. https://doi.org/10.1785/0220130035
[57]	Yanez, G., Ranero, C., Huene, R., & Diaz, J. (2001). Magnetic Anomaly Interpretation across the Southern Central Andes (32-34 S): The Role of the Juan Fernandez Ridge in the Late Tertiary Evolution of the Margin. Journal of Geophysical Research, 106, 6325-6345. https://doi.org/10.1029/2000JB900337
[58]	Yuan, X., Sobolev, S., & Kind, R. (2002). Moho Topography in the Central Andes and Its Geodynamic Implications. Earth and Planetary Science Letters, 199, 389-402. https://doi.org/10.1016/S0012-821X(02)00589-7