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

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) in-ter-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 orographical-ly-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


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 ( τ ) resulted from differential plate motions.
For an ideal configuration (Figure 1(a)) of uniform fault geometry (i.e., constant slope angle θ and unlimited width in the transverse direction) and uniform material (so that maximum frictional coefficient 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).

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.
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 ( ) ( ) . S is thus a useful diagnostic index for fault stability.

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 τ 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 τ , G, and f θ (or fault strength µ′ ) 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.  µ is the static friction; G is the gross gravity of the overlain plate section; τ 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 me- 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)).

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 µ′ , 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 i 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).

D. D. Ren Journal of Geoscience and Environment Protection
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 3 ) 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.

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.

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 ( ρ ) according to a square root law: 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 macros- Journal of Geoscience and Environment Protection copic 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 4 Pa S) many orders of magnitude smaller than polycrystalline rocks (>10 21 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 −4 , 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".  (2)), they also affect neighboring regions (laterally) through bridging effect and the unique fatigue process of the fault zone.

Fatigue in Rock Interface-Generation of Granular Material
Two objects, when being pressed against each other, wear and tear at the inter- 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; 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 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.

Numerical Model Setup, Parameterization of the Earthquake Triads, and Data Acquisition
Above discussed theoretical concepts, or the Groundwater aided Granular material Generation-G 3 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 4 km 2 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.

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  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 D. D. Ren 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 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.,  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.

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: where a is the thickness of the granularized rock layer, starting from 0, t is time, c 1 is the e-fold time scale, and c 2 is the Paris coefficient (Timoshenko & Gere, 1963). Apparently, c 1 signifies the negative feedback as granular material accumulates. The unit surficial energy density of rocks, J (<300 mJ/m 2 , Gilman  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 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.  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 1 ). 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 1 , we did not involve the healing mechanism in Equation (3) (Rowe et al., 2012). For the inter-seismic periods, Equation (3)  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 µ′ value of 10 −3 , 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.

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.

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 . 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 w 7.5 and >M w 5.5. Without explicit statements, >M w 5 over other world regions is regarded as major earthquakes.

Earthquake Activity in the TP Region: Correlations with Groundwater Fluctuations
Large scale (e.g., thousand kilometers) groundwater fluctuations exert extra  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 2 saturated area (as defined in Equation (2) 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 w 7.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 w 7. 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 D. D. Ren Journal of Geoscience and Environment Protection 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. Journal of Geoscience and Environment Protection gence 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)).

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). . 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). 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 w 8.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 sec-

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 w 7, 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 W 5.5 and M W 7, all major earthquakes greater than M w 7.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 w 7.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.

Cascadia Subduction Zone's Stability Affected by Californian Droughts
Detailed discussion of Cascadia subduction fault, such as effective strength ( µ′ ), geometrical characteristics and unique geological conditions can be found in Wang and He (1999 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 13 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 w 7.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 w 7.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).

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 sev-  (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 interpo- where in geoscience convention, θ is longitude, λ is latitude and r is the distance from the earth center ( Figure A1). Still in right-handed coordinate system, Notice also that transpose of I 1 also is inverse of I 1 (I 1 is thus an orthogonal matrix). In addition, I 1 also is unitary and normal matrix (det(I 1 ) = I).
From Equation (A9), 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 ( In Figure A2, the Burgers rheology is represented by a serial connection of a