How Droughts Influence Earthquakes

Earthquakes result from strain build-up from without and weakening from within faults. A generic co-seismic condition is presented that includes just three angles representing, respectively, fault geometry, fault strength, and the ratio of fault coupling to lithostatic loading. Correspondingly, gravity fluctuations, bridging effects, and granular material production/distribution form an earthquake triad. As a dynamic constituent of the gravity field, groundwater fluctuation is the nexus between the triad components. It is pivotal in regulating major seismic irregularity, by reducing natural (dry, or purely tectonic, stationary seismicity) inter-seismic periods and by lowering magnitudes. Specifically, to exert stress on the fault, groundwater does not need to reside deep in proximity to the locked fault interface, as it can work remotely. It can act mechanically-direct (MD), by a differential de-loading and superimposing a seismogenetic lateral stress field, thereby aiding plate-coupling, from without, or mechanically-indirect (MI) by enhancing fault fatigue, and hence weakening the fault, from within. To verify this hypothesis, gravity measurements, and a numerical model, are used. The remote action hypothesis is globally applicable. Detailed results are presented for the Himalayan and New Zealand regions. The gravity recovery and Climate experiment (GRACE measurements) reveals that major earthquakes (Mw 5 and above) always occur in the dry stage, indicating drought and associated groundwater extraction is an important trigger for major earthquakes. By exploring 73 historical records successfully reproduced by the model, it is found that for collisional (e.g., the peri-Tibetan Plateau) and strike-slip (e.g., the San Andreas Fault) systems, the MD mechanism dominates, because the orographically induced spatially highly variable precipitation is channeled into greater depth by through-cut faults. Droughts elsewhere also are seismogenetic, but likely through MI effects. In a warming future climate, mechanisms identified here play a greater role in increasing the recurrence frequency of major earthquakes, but also in slightly reducing their severity.


Introduction
ing crust, bridging effects, and fatigue [4]. With the fault strength set as constant, the fault geometry and plate-coupling stress determine a natural, limiting, earthquake occurrence frequency [5]. However, the observed earthquake occurrence seldom obeys this limit and typically occurs well ahead of schedule, due to various fault-weakening factors. Among many other first order mechanisms explainable by Eq. (1), the present study singles out the weight loss (i.e. de-loading) of the upper crust as seismogenetic. As recently reported, for the central Californian valley [6] the weight loss of the upper crust can result from extended droughts that deplete the soil moisture content residing in the rhizosphere and/or by groundwater loss (extraction). By synthesizing geodetic dislocations, obtained from radar and GPS measurements of surface deformation, using an elastic dislocation model, Gonzalez, et al. [7] found that the areas of fault slip, along the Alhambra de Murcia Fault of southeast Spain, correlated well with positive Coulomb stress changes induced by groundwater extraction from nearby aquifers. Water table decreases of hundreds of meters were detected. Thus, the nucleation and main rupture of the Mw 5.1 Lorca earthquake (May 11, 2011) were assumed to be attributable to a loss of groundwater. In the following discussion, groundwater is used as a generic term for both soil moisture, and for drainage into a groundwater reservoir. For relatively weak faults (e.g. subduction megathrusts), the lateral stresses exerted by groundwater extraction are of the same order of magnitude as the fault strength. The normal stress component exerted by groundwater transfers readily into lateral stress earthquakes, especially the role of groundwater, and water residing in the rhizosphere, in affecting seismic activity. Tectonic earthquakes are frictional phenomena between plates in contact [1]. At the wedge junction, the upper plate's weight (or more generally, the compressive stress) aids friction in resisting the increasing plate-coupling stress τ resulting from differential plate motion. Figure 1 is an idealized configuration of uniform fault zone geometry, which assumes a constant slope angle, θ , unlimited width in the transverse direction, and uniform material with a constant maximal frictional coefficient In this cases, the co-seismic condition is: where the repose angle is = arctan( /G) φ τ , and G is the compressive stress on the shear (fault) zone. The subscript 'c' means the critical value, θ is the slope of the fault zone, determined by the geometry of the plate interface, and f θ is the maximum static friction angle corresponding to fault strength. Values of f θ depend on wall rock material properties and lithology [2,3]. Between seismic events, the repose angle increases steadily and approaches the sum of the fault slope angle and the maximum static friction angle. The driving stress ( ) τ increases and the overall resistive stress decreases (i.e., f θ decreases); both are seismogenetic [1]. In reality, erosion from within and the loading from without likely are simultaneous. Earthquake cycles are stress adjustment cycles. Closely related to the three angles in Eq. (1) are the earthquake triads: i.e., weight fluctuations of the overrid- Figure 2, the upper (usually continental) crust is represented as four "carriages" of a train, each carriage signifying different segments of the upper crust, which is in force balance between the frictional force at the bottom and the drag in the front and therefore it is stationary. Depending on how each carriage is footed (analogous to the friction between the upper and lower crust plates), there are several possible configurations/scenarios (the four rows in Figure 2). If the elastic plate can maintain rigidity over the entire shear zone, all segments of the shear zone must reach stress saturation for coherent sliding to occur. Prior to that critical moment, all segments 'saturate' at the same rate, as the compressive driving stress increases (the top train in Figure 2). In reality, earthquakes involving the 'saturation' of the entire shear zone seldom occur, because tectonic plates have difficulty in maintaining rigidity over extended distances, manifested as the existence of through-cut faults, heave bumps and many small magnitude earthquakes within close proximity of each other. Consequently, many locations do not contribute their shares of resistance and rely on neighbors to stay in place.
Depending on how τ is distributed on the fault zone, there are various possible combinations of saturated and partially saturated segments along the fault zone (trains 2-4 in Figure  2). In Figure 2 it is assumed the distance that plate maintains its rigidity also is the horizontal dimension of each segment. Thus, once saturated, a corresponding co-located earthquake can occur. For example, in the second train (second row of the scenarios) the first and the last segment are saturated. As they are separated, two small magnitude earthquakes (or one event with two small ruptures [8]) are generated. The third train has two adjacent segments saturated. In this case, in the fault environment. The following discussions use subduction faults and collisional systems with thrust faults as prototypes. By substituting compressive stress for lithostatic loading, the same reasoning also applies to strikeslip faults, as explained in the Supplementary Material of Ref. [6].
Displacements caused by a major earthquake (Mw of 5 or greater) usually range between 0.5-3 m, i.e., they are at least 3 orders of magnitude smaller than the dimension of the seismically locked fault zone, which can be several hundred kilometers. Consequently, material in the locked fault zone has experienced numerous such wearing and tearing events. When simulating near-future earthquakes, the geological backgrounds, including topography, material properties and fault geometry, can reasonably be assumed to be constant. Thus, the natural limiting frequency of seismic slip/rupture occurrence of a fault is highly stable. The irregularity in earthquake occurrence is primarily due to fault weakening mechanisms and to external forcing from the upper boundary. This study focuses on more dynamic fluctuations affecting the earthquake triad. In the Supplementary Material of this study (SM), a stability index is introduced, based on Eq. (1) to explain, simply, the interaction of the triads and the role played by groundwater fluctuations.
As an under-determined problem, multi-solutions (i.e., uncertainties in the causal mechanisms) in seismogenesis can result from insufficient information about the lateral interactions among different segments of the overriding crust ( Figure  2). Rather than the self-weight, the groundwater is influential as it affects the lateral interaction between crust segments. In  (-). Depending on the distribution of lateral pressure (i.e., how the 4.8 MPa is partitioned), the bottom frictional forces are of very different magnitudes and even in directions. Compressional pressure is shown as positive only because here most numbers are compressional. However, in the numerical model, compression stress is negative by convention. The degree of friction saturation is color shaded. Filled Heptagrams (convex seven pointed stars) indicate the epicenters.

Materials and Methods
The SM outlines the interactions of the earthquake triads in the simplest possible manner. From Eq. (A1), the weight of the overriding plate always serves as a stabilizing factor. A weight reduction (e.g., from extended droughts) contributes to instability as it is seismogenetic. Fluctuations in the weight of the overlaying plate not only have local stability consequences, they also affect neighboring regions through the bridging effect and a unique fatigue process of the fault zone [13,14], by GM generation and transportation, a common fault-weakening mechanism. Two objects, when pressed against each other, wear and tear at the interface, as a manifestation of material fatigue. The interfaces of tectonic plates are no exception. Production of granular debris is the primary form of fault interface fatigue. Regions of strong contact, likely as a result of the bridging effect, enhance active granular material generation. GM has a much smaller viscosity and acts as a lubricant for the plates involved. Regions carrying more loading now effectively rest on a 'slippery floor' (i.e., more weight is loaded on smooth contacts and less is loaded on rough contacts, along the shear zone). Thus, the bridging effect, by affecting granular material generation and redistribution, influences the frictional properties of the fault interface. It therefore achieves a negative spatial covariance between loading and fault strengths, thereby being seismogenetic according to Eq. (A1). Sources of bridging effects are numerous. For example, large scale structures in the overriding plate, such as mountains and valleys, signify bridging effects at depth on the seismogenetic fault zone [15]. Bridging effects from static geological backgrounds, although they may cause variations in the frictional properties of the seismogenetic zone, change only slowly and are not responsible for short-term variations in major earthquake occurrence. Large scale variations in groundwater, however, are more dynamic and thus are viable candidates for explaining earthquake occurrence irregularity on decadal to centennial time scales. Cyclic forcing is the most effective means of causing fatigue [16]. Cyclic groundwater fluctuations, especially when acting in concord with the resonance frequency of the plates, assist GM generation.

Groundwater time series within and outside the gravity satellite measurements period
Fluctuations in groundwater, which are remotely sensed by gravity satellites such as the Gravity recovery and Climate experiment (GRACE) [17][18][19] are used to highlight the strong correlations with the timing of major earthquakes. Launched in 2002, GRACE has revolutionized the detection of large-scale mass changes on Earth. Water changes, irrespective of whether they result from accelerated ice sheet melt, extreme precipitation events, or extended droughts, are the dominant drivers of regional mass fluctuations. In addition to measurement noise, the global monthly one-degree data used in this study have the following missing months: a larger magnitude earthquake is triggered. The fourth train is unique in that it involves interlaced compression and dilation among the segments, corresponding to a saddle-shaped strain distribution. Consequently, the three adjacent segments are simultaneously saturated and a super earthquake is formed.
Even in the case of this highly idealized fault, with its uniform maximum friction coefficients and uniform lithostatic loading, very different spatial patterns are possible, fostering earthquakes of drastically different magnitudes and at different spatial locations. These examples represent only a fraction of the myriad possibilities that can unfold in reality. In particular, earthquake occurrence can be compounded by any factors that affect τ , G, and f θ (or fault strength µ') in Eq. (A1). Multi-solutions in seismogenesis result largely from inadequate information as constraints on, or under-determination of lateral interactions.
For example, historical data should be treated with caution when predicting future earthquake scenarios because they usually do not provide information on the rupture length, width and slip; they provide essentially no indication of the horizontal coherence of the different plate segments. Earthquakes, because of the relative motion of plates involved, inevitably have footprints on the mass [9,10], energy [11] and momentum fields (Ref. [12] and references therein). Hence, earthquake predictions are increasingly based on observations of these physical parameters and their fields. The mechanism discussed here, i.e. the seismic consequences of groundwater fluctuations, fits this situation exactly. It is a hitherto largely neglected factor that plays a far larger role than anticipated in the existing approaches, which include only local weight effects as some extra weight (negligible compared with the ~20 km thick overriding crust plate). The effect of the groundwater fluctuations is magnified by the bridging effects and the related granular material (GM) generation, which is a fatigue for brittle material. Groundwater fluctuation caused fatigue accumulated temporally because there is no healing mechanism for the granular fatigue and each leaves an accumulative memory on the plate interface. Unlike the reported good correspondence between minor seismic activity on the San Andreas fault and the central California droughts and groundwater depletion [6], there seems to be no one-to-one correspondence between extreme precipitation and drought cycles and the large earthquakes. It is a multiple drought-wet cycles preceding one major earthquake phenomenon. More critically, the footprint in the lateral stress field accumulates spatially, especially along the frontier mountain chains that act as waveguides by maintaining plate segments' rigidity for extended distances. As it is a grand challenge, it is not surprising that there is no simple equation available that predicts earthquakes. All expressions and equations here are intended to convey the concepts of the earthquake triads. A sophisticated numerical modeling approach, as elaborated in the next section, is needed for reliable seismic hazards predictions.
nomena. However, the prediction of earthquakes is, however, very challenging and cannot be summarized in just one or a few simple equations. Formulae are used here to explain only the physics involved. At present, the only viable approach is a coordinated solving of a numerical modeling system that includes energy, momentum, and mass conservation principles, under prescribed rheological properties. SEGMENT is such a realization. The current rubric on earthquake prediction and forecasting (e.g., Ref. [24]) comprises four steps: The fault model provides physical geometry of the faults (especially large and active ones) considered; the deformation model provides slip-rate and creep estimates for each fault section, and derives the deformation rate of the modelled faults; the earthquake rate model provides long-term rate of all earthquakes throughout the region of interest; and, finally, an earthquake probability model provides the actual earthquake prediction. SEGMENT essentially covers the first three models (with the fault physical geometry and mechanical parameters as inputs to the modeling system). There is no application of probabilistic methods in earthquake forecasting and the epistemic uncertainties are assumed to be model uncertainty (e.g., model parameter setting and input data inaccuracies), which reflect the still incomplete understanding of earthquake mechanisms. In actual predictions, SEGMENT modeling is a set of parallel code and requires super computers for execution.
An understanding of the basic earthquake mechanics presented here is based on numerical simulations from SEG-MENT. SEGMENT is applied as a global spherical crown model, with the stress decomposed into driving stress and resistive stress components, following Van der Veen and Whillans [25]. By solving the equations for conservation of mass, momentum and energy, SEGMENT provides a 3D distribution of strain, stress, and the three components of velocity in the simulation domain. Identifying a major earthquake from these prognostics is described in Section 2.3. In addition to numerical details such as fault-following grid stencils, Section 2.3 also describes approaches in obtaining the current global distribution of GM thickness and its future evolution.

Numerical model setup and parameterization of the earthquake triads
The theoretical concepts/hypothesis discussed above and also in the SM, namely, GM generation facilitated by groundwater fluctuation, is implemented in a global spherical crown model SEGMENT. As a dynamic model unlike kinematic and empirical statistical models, SEGMENT solves the equations of mass, momentum (under assumed multiple rheology), and energy conservation in 3D spherical geometry. The model prognostic equations include the three components of velocity; temperature; and six components of the deviatoric stresses. Their evolution is depicted by forward time marching of the governing equations. Variables in Eqs. (1) and (A1), such as the degree of stress saturation, are diagnostic equations of SEGMENT rather than prognostic equations. Through these diagnostics, the stress buildup during the inter-seismic stage can be followed. A major earthquake in the model is identified from the predicted from neighboring months. In total, there are 15 years of GRACE measurements of time variations in gravity changes over the globe, before the quality of the data degrades. Before 2002, or for areas with land surface areas that are too small to be resolved by GRACE data, a land surface scheme within a model referred to as Scalable Extensible Geofluid Modeling framework for ENvironmenTal issues (SEGMENT [14,20,21], will be detailed in Subsection 2.2), driven by atmospheric parameters, is used to simulate groundwater fluctuations. The earthquake triad is interlinked with groundwater fluctuations. Groundwater-aided GM generation is the mechanism modulating earthquake cycles. As an initial attempt at supporting this hypothesis, global evidence is provided of the earthquake events and gravity fluctuations in the neighboring regions. The present study uses the global seismic records from the USGS Significant Earthquakes Archive (https://earthquakes.usgs.gov).
From Eqs. (1) and (A1), groundwater, by affecting the gravity field, contributes to fault instability. Its spatial fluctuation, through enhancing GM generation, also weakens the fault from within (Section 2.2 of the SM). To obtain the groundwater fluctuations over the past years from 2000-2015, a continuous time marching of the land surface scheme is performed with a 30-minute time step. Atmospheric parameters such as near surface air temperature, humidity, pressure, winds, and radiative components are from the NCEP/NCAR reanalysis. Precipitation data from the Global historical climatology network (GHCN, 1900(GHCN, -2015, GPCP, and The tropical rainfall measuring mission (TRMM, 1998(TRMM, -2015 are taken as inputs 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 (2002 onward), GRACE measurements provide a reliable means for the groundwater recharge/discharge. Model simulations are found to be satisfactory when verified against GRACE measurements. GRACE has a very coarse resolution and SEGMENT can provide the spatial distribution of groundwater at much finer resolution. Obtaining groundwater distribution in a warmer climate is challenging. A small ensemble of 27 climate models, listed in Table S1, and an innovative approach in estimating extreme precipitation [20] is used for the prediction of the future occurrence of earthquakes. The partitioning into surface runoff, evapotranspiration and ground water percolation also is performed within SEGMENT.
In addition to the geodetic measurements and atmospheric parameters, the earthquake model used here also requires some static data, such as fault geometry, material property and thermodynamic parameters of the entire simulation domain. In summary, CRUST1.0 [22] is used to set up the geometry and physical properties (e.g., density values) over the grids within the crust layer. It refers to the Advanced Solver for Problems in Earth's ConvecTion (ASPECT) [23] for setting up the thermal and initial flow fields. Model spin up runs are constrained by present GPS measurements of the surface elastic deformation.
setting up the grids above present sea level, the Advanced space-borne thermal emission and reflection radiometer (ASTER) global digital elevation map (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, 102-106°E; 30-33°N) of the eastern margin of the Tibetan Plateau, there was a detailed gravity survey after the Mw 7.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 [26] (and data pointed out by references therein) also provided information that is merged into the coarse resolution CRUST 1.0.
Based on the fact that the co-seismic ruptures consume only a portion of the plate interface [27] to the center and the down-dip section creeps, the upper and lower portions of the plate should be treated respectively as elastic 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 viscoelastic [12] but becoming more rigid upwards intending a smooth transition in the vertical domain. This is achieved by considering the thermal structure [11] 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 (~75-260 kJ/mol for natural rocks), specific volume V and universal gas constant R (8.314 J/mol/K) all are assumed to be 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 S2). In most cases, it is only a small fraction of the shear zone layer that is granular material (also recorded in arrays). The parameterization of GM viscosity follows Ref. [13] as implemented in SEGMENT [14]. The shear-thinning viscosity structure (Eq. 8 in Ref. [14]) implies a shear-localizing effect from GM accumulation. The present amount of GM ( Figure S2c) is deduced from surface velocity using the method in Ref. [21]. Future variations of GM are diagnosed (relevant equation is detailed in the next subsection on GM generation that immediately follows). The granular ensemble size is deduced from its parent rock type and the associated grain sizes. Re-mapping of the grid stencils is automatically performed only if the accumulated GM exceeds 20 m or the earthquake caused displacement is greater than 20 m.
The lower boundary stress conditions are critical for the simulation of stress build up 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 motion velocities. Once an event is identified, the released energy is estimated within a time window that covers fully the time period with macroscopic slip. An epicenter is diagnosed as the geological locations of the saturated grids, weighted by the respective energy release. The degree of stress saturation S is not directly used to identify earthquakes, because the grid elements are not isolated; there are interactions among adjacent grids. If the unlocking/ rupture at a grid point releases large amount of energy, a neighboring grid point, even if distant from (stress) 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) are 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.
Grid stencil of SEGMENT: For this study, the original 3D spherical model of the upper ~4 km medium (suffices for ice sheets and landslides applications) is now applied to the 500 km thick lithospheric material (including elastic crust and visco-elastic upper mantle). The horizontal resolution is generally 30 km at the Earth's surface and refined to 2 km at the subduction zones ( Figure S1). The vertical resolution varies using 101 stretched vertical layers to represent the 500 km depth. A stretching scheme is used so that the shear zone gets higher resolution and is better represented. In setting up the model grids in the simulation domain, it is ensured that there always is a vertical level that has a surface parallel to the oceanic plate's upper surface. That is, SEGMENT uses a fault zone-oriented coordinate system.

Static geological parameters and viscosity parameterizations:
The SEGMENT had been used satisfactorily for landslide research in the past. As frictional phenomena, landslides and tectonic earthquakes share high degree of similarity. Both involve a shear zone and weakening of the shear zone leads to instability. Reducing the weight of the overlying portion also is prone to instability for landslides. For example, around manmade reservoirs formed by occluding a river branch with a dam, landslides along the banks always occur during the discharge stage of the reservoir, same de-loading mechanism as presented in Eq.
(1) for earthquakes. Although the similarity between landslides and tectonic earthquakes makes the conversion of SEGMENT into an earthquake model feasible, the pressure and thermal environments of the fault interface and the shear zone of landslides are very different. Special attention to the viscosity parameterization is required. Also, direct verification of the schemes become difficult for earthquake applications as the plate interface is deep-buried that direct measurements of either stresses or dislocations are rarely available. Inverse methods are necessary for model verification, spin-up and initialization.
The CRUST1.0 module is used to set up the density values over the grids within the crust layer (both the oceanic and continental crusts). Underneath the bottom of the crust layer is assumed to be upper mantle material (silicates), and is given a uniform reference viscosity of Compared with GPS measurements on orogens and plates away from the faults, globally, errors in meridional velocity are < 0.2 mm/yr and in the longitudinal direction are less than 0.3 mm/yr. For the Andes region, the error is slightly larger than the global average (0.35 and 0.3 mm/yr, respectively in the two directions). Over the Cascadia region, however, the errors are smaller than global average (0.27 and 0.15 mm/yr respectively). Based on the optimized solution of neotectonics, 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 similar to those of the Advanced solver for problems in earth's convection (ASPECT) [23], it is noted that the ASPECT model, due to its intended research purpose, has coarser 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 between model top-level velocity and GPS measurements. By using careful numerical settings, SEGMENT avoided 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 3a).
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, it is the Yellowstone area), does the boundary layer flow structure indicate 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 is 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., . 0 n τ → = 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, (weather signal pre-filtered) gravity signal and thermal flux measurements available.
Accurate flow fields are critical for accurate estimations 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 averages. The geodetic data over land mass in panel (a) were not collected simultaneously (actually over 1991-2003) but spanned a period of many major earthquakes (e.g., the Mw 6.7 Northridge earthquake). Because major earthquakes during the time window of geodetic data collection may leave coseismic rates. The estimated inherent resonant frequencies are labeled in Figure S3.
The inter-plate GM chute system, as left behind by seamounts, removes granular deposits through an advection process (term ADV in Eq. 2) exactly the same way as surface runoff moves debris [32]. Transport of GM is left to the Navier-Stokes solver of SEGMENT. As in Figure S2b, the strong contact regions have the most active GM generation. Its existence in the locking zone assists the unlocking because the shear resistance GM 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. However, with the associated GM 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 [1]. GM formation is one feasible explanation for slide weakening. The generation and transport of GM at the plate contact can reconcile many counter-intuitive characteristics of major earthquakes [33] and also is critical in organizing small-scale ruptures to generate major earthquakes. Correct treatment of present GM ( Figure S2c) between plates also may provide a timeframe for the next major earthquake on the same fault. From Figure  S2c, it appears that, globally, the GM thickness distribution is positively correlated with major earthquakes occurrence frequencies. The following sections first discuss the correlation of earthquakes and fluctuations in gravity fields as a result of groundwater extraction, using two faults of different geometry and strength, and located in very different climate zones. Using another two subduction faults, the upcoming earthquakes as a result of hydrological cycle changes with The above discussion gives the global setting of SEG-MENT. For applications to specific regions, existing models over that orogen are also referred to in the parameter setting. For example, in the peri-Himalayan regions, the static model of Ref. [5] on Persia-Tibet-Burma orogen is referenced. The stationary creeping fields (as they use a 2D model, only one level of flow fields is available) are used to constrain SEGMENT's rheological parameters setting in the fault vicinity, especially near fault junctions. Similarly, Ref.
[29] is referenced in setting the tectonic environment of New Zealand, especially the numerous fault traces.

Granular material (GM) generation and transportation:
In analogous to the state equation in Ref. [1], the following force-restore relationship is proposed here for GM generation: where α is the thickness of the granularized rock layer, starting from 0, t is time, c 1 is the e-fold time scale, c 2 is the Paris coefficient [30], and c 1 signifies the negative feedback as GM accumulates. The unit surficial energy density of rocks, J (< 300 mJ/m 2 , Refs. [1,31]), is inversely proportional to c 2 . To form new surfaces at high confining pressure (~ 0.1 GPa), the energy requirement far surpasses the surficial energy. Thus c 2 also is inversely proportional to the confining pressure. Thermal effects also are included in c 2 . The coefficient δ is 0 when the stress intensity change falls within the elastic range of the parent rock and is 1 when exceeding the elasticity range. DE is the energy input rate when stress perturbation caused imbalance in principal stress exceeds the elastic range. External (e.g. groundwater fluctuations) tides with resonance frequencies are most effective in transferring distortion energy into the plates to cause fatigue and influence the GM production errors. As with all prediction systems, the uncertainty grows exponentially with forecasting lead time. In the following discussion, the uncertainty level will be presented as necessary.
Although global results have been obtained for earthquake correlations with fluctuations in gravity fields resulting from groundwater discharge and extraction, the following discussion focuses on two regions: The Hikurangi subduction fault-Alpine fault-puysegur mega thrusts of New Zealand (HAP), and the Tibetan Plateau and its surrounding regions (TP, a continental collision earthquake zone). Both regions are adequately covered by GPS measurements and have adequately sizable land areas that can be sensed by GRACE for the mass fluctuations around the faults. As a result of different fault geometries, plate opposing motion speeds, and the distribution patterns of inter-plate GM, earthquakes over the selected regions are of various morphologies. Major earthquakes over the TP and HAP regions refer respectively to earthquakes of magnitudes > Mw 5.5 and Mw > 5. In the Supplementary Material (SM), applications to subduction earthquake zones are further discussed for two other earthquake-prone regions: the Andes and the Cascadian subduction zones, along with global maps of time scales of plate resonance responses to external forcing by groundwater fluctuations.

Earthquake activity in the New zealand region
The weight of groundwater alone is unimportant as a normal pressure component in a fault. However, the lateral climate warming are discussed, as actual applications of the numerical earthquake model discussed in this section. Limited by size, these two subduction faults, namely the Andes and Cascadian faults, are addressed in detail only in the SM.

Results
There are three interlaced mechanisms working together to make the earthquake prediction difficult for present existing rubrics. The mechanisms identified here are tested positively against 73 historical cases, all within the NCEP-NCAR reanalyses period, with high quality precipitation data. The reproduced earthquakes are close to those observed, especially for cases after 2003 (the start time of GRACE measurements) and before 2015 (the aging of the current GRACE satellites). By simulating the time scales of the evolution of strain and stress in the fault zone, using SEGMENT, with appropriate parameter settings, assisted by remote sensing information, the timing and location of major earthquakes (Mw > 5 globally) can be estimated. The uncertainties in timing are reduced to a year, and the location errors to within 120 km (the root mean squared error for historical records during 1979-2015). The primary limitation on the prediction accuracy is the resolution of the CRUST-1.0 model, depending on which SEGMENT resolves fault geometry and composing material properties. Figure 4 is the hind-casts of historical earthquakes greater than Mw 8. It is apparent that simulation of cases within the remote sensing era has significantly smaller spatio-temporal errors. Over the same fault zone, those cases within the GRACE observational period have the smallest simulation terface. This is a mechanism for the temporal accumulation of groundwater effects. From the two major mechanisms described above, many seemingly distinct earthquake events become explicable. A further examination of the 73 historical records worldwide suggests the above groundwater hypothesis likely is universally applicable.

Earthquake activity in the Tibetan plateau region
Large scale (e.g., 1000 km) groundwater fluctuations exert additional stress on a fault interface. From Figure 6, almost all major (> Mw 5.5) earthquakes during 2003-2016 over the Tibetan Plateau also occurred during low groundwater phases (Figure 6b), similar to the New Zealand cases of Sec 3.1. Although the fluctuation of groundwater is secondary to plate motion in causing stress build-up, its addition can trigger an earthquake. Analyzing the stress field clearly indicates that once two plates are coupled together, the underlying plate drags the lighter overriding plate downward and they move in tandem deeper into the Earth, pushing aside the higher density upper mantle material. This causes mass loss around the fault region. This process, although only causing mass loss at a rate of two orders of magnitude smaller than the 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 stress it exerts on the surrounding region is critical for fault instability. The second critical point is that, through fatigue mechanisms, groundwater extraction leaves a memory on the fault interface and forms a mechanism for temporal accumulation. Figure 5 shows the occurrence of major earthquakes over New Zealand since 1995. All major earthquakes occurred during the relatively dry periods (low groundwater levels as indicated by the low spikes in the blue curve of gravity anomalies). This high correlation is not by chance. Using SEGMENT, the transient stress field exerted by the groundwater deficit (relative to a long-term mean) is simulated for the Kaikoura Mw7.8 earthquake on November 16, 2016 (Figure 5b). Kaikoura resides in a maximum region of the R xz stress field exerted by groundwater deficit, which reaches 5 kPa. This additional stress, even if was not the effective trigger for the major earthquake, is sufficient for triggering the following slow slip events (SSEs, personal communication with L Wallace [9]). For the events shown in Figure 5, it appears that most of the high correlations result from the groundwater deficit being an effective trigger rather than being a direct cause, even though the shear stress exerted by gravity fluctuation (groundwater extraction through net evaporation) can be orders of magnitude greater than the weight of the groundwater in the fault conditions (e.g., horizontal pushing as in the train carriage illustration in Figure 2). Also notable is that, through bridging effects and GM generation, each cycle of groundwater fluctuation leaves a 'mark' on the plate in- . On average, the groundwater reduction rate was ~31 Gt/yr. This super-imposed stress field produces a saturated area (as defined in Eq. (A1)) of ~85 km 2 annually. The matured (saturated) area is not spatially contiguous on the fault interface; it scattered across the inter-seismically locked interface in a complicated-connected manner. Connection of these sporadically scattered saturated patches is critical for generating earthquakes. For example, the occurrence of the first earthquake (April 25, 2015; Mw 7.8; epicenter at 28.147°N, 84.708°E) redistributed the granular material due underneath and established a stronger rock-rock contact, lessening the degree of saturation of the second major earthquake (May 12, 2015 of M7.3 magnitude and epicenter located 200 km to the southeast at 27.973°N, 85.963°E) region, so another month had passed before its occurrence. The earthquakes over Burma (Mw 6.9 on April 12, 2016 at 23.13°N, May, 2015, and another, third, earthquake in Burma a year later all reside within the mass loss region. The temporal sequence of these three major earthquakes and the ruptures' spatial pattern (starting from northwest and extending to southeast) can be explained by the GM production and its transportation, 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 follows the GM accumulation in the shear zone.
The three epicenter-consisting patches (comprising the three respective epicenters) all were approaching saturation by April 2015, after ~72 years of stress accumulation [34-36]. The western-most epicenter ruptured first because it was aided by groundwater reduction. Based on SEGMENT simulations, the groundwater recharge over the past 80 years has  33°N and 35°N) to the north east, ~ 120 ± 45 years later of magnitude ~ Mw7. For the far-larger TP region, the India Eurasian collisional system, while thickening the elastic lithosphere, has created numerous vertically extended crevasses. These crevasses are very effective in channeling groundwater to greater depth. Combined with the TP's orographic effects, precipitation variation transfers readily into extra lateral stress on the faults. Therefore, groundwater fluctuations become a dynamic factor affecting the epicenter location and the magnitude of major earthquakes. Because of 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 has formed and has persisted over the past half century. The natural, geological, recurrence frequency (background-deduced "dry" frequency) is greatly reduced over the region, especially over the Southeast Asia peninsula. The fault zone over Myanmar likely will become unstable in approximately two years. Of great interest is that the epicenters of the upcoming earthquakes, affected by the weakening of monsoonal precipitation, are located along the 96 ± 1°E parallel between 17-28°N (most frequently in the 17-23°N sector), roughly parallel to the Ayeyarwady River basin.
In addition to the pre-Himalayan and the New Zealand regions, the authors have shown that there are strong, statistically significant correlations between tectonic plate earthquake irregularities and groundwater fluctuations over many other global faults, agreeing with several recent studies (e.g., Refs. [6,7]). For example, the central Italy earthquake (August 24, 2016) and the Lorca earthquake (May 11, 2011) both occurred at low groundwater level stages. The main value of SEGMENT is its predictive ability. Previous models only involve the groundwater induced crustal unloading, or the direct mechanical effect of groundwater as in Eq. (1). Without considering the fault weakening effects from groundwater, these models explain reasonably well the observed fault slip patterns resulting from groundwater depletion. However, predicting earthquakes is beyond the capability of such models. The SEGMENT, by placing groundwater fluctuations at the nexus of all prongs of the earthquake triads, especially emphasizing the effective fault weakening from the cyclic variations of the induced crustal loading, possesses the ability to simulate the role of groundwater in future earthquakes. Using the SEGMENT modeling system, driven by precipitation scenarios from climate model projections, the future occurrence of earthquakes over the Cascadian and Andes subduction zones is investigated (in the SM). It was found that fatigue caused by groundwater fluctuations, associated with increased precipitation volatility in 21 st Century California, will likely be decisive in shortening the time scale of the natural earthquake cycle of the Cascadian subduction fault (SM). However, as a consequence, their magnitudes will be reduced. For the Andes region, in addition to the regular precipitation pattern from ENSO, the granular chute system left by seamounts on the Nazca and Jun Fernandez Ridges foster major earthquakes. 94.9°E and Mw 6.8 on August 25, 2016 at 20.9°N, 94.6°E), further to the southeast, were not directly related to the first two earthquakes. They were a result of the lightening of the Bay of Bengal region, to be further addressed, below. Unlike some recent studies 24 predicting that the future earthquake after the Nepal earthquake of 2015 would be in the sector to the northwest, the SEGMENT simulation indicates that the sector to the south-east also is becoming unstable. In actuality, the maturation pattern jumps between the eastern and western sectors, with cluster centers at (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 LFZ is a classic example of the 'non-unique solution' configuration presented in Figure 2. The LFZ defines the eastern margin of the Tibetan Plateau. A very recent synthetic approach, using a gravity survey, seismic reflection profiles and earthquake focal mechanism data, and sediment analysis indicated that the LFZ essentially is composed of two segments of different formations. The central-northern section, which was formed ~ 40 Myr ago [38,39], is a lithospheric throughcut fault zone that is of low elastic strength but has 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 favourable conditions for instability, also has a profound influence on the groundwater reservoir. The Mw 7.8 Wenchuan earthquake offers further evidence of the importance of gravity field changes in fostering major earthquakes.
The mass change fields averaged over the period 2003-2008 indicate clearly that Wenchuan resides at the center of a saddle region (Figure 3 in Ref.
[10]). As illustrated in the 4 th scenario train in Figure 2, the superimposed stresses form a strong compressive stress in the horizontal direction, aiding the existing geological compression. For earthquakes occurring along the same fault (but temporally consecutive), the compression caused by tectonic plate motion usually produces very similar saturation patterns. The primary reason for the randomness of epicenters is the increase in transient perturbations, such as groundwater fluctuations. The superimposed field causes a specific point to reach saturation first and then propagates away to cause a large rupture area. The rupture can grow because, to halt the motion of the unstable plates, neighboring unsaturated regions may be activated. If the propagation encounters too many unsaturated regions, the momentum is quickly lost, strongly limiting the magnitude of the earthquake. However, if the propagation starts from a region surrounded by a nearly saturated neighborhood, a positive feedback forms and is supportive of large magnitude earthquake formation. By initiating the rupture near Wenchuan, a more efficient energy release is achieved. From the SEGMENT simulations with future precipitation scenarios under the IPCC SRES A1B scenario provided by a small-volume climate model ensemble (as listed in Table 1), the next major earthquake is expected to occur in the Ganshu province (epicenter located along the fluctuation in the gravity field. Similarly, the stability of the Cascadia fault is found to be 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. Thus, climate fluctuations can and do determine the groundwater loading pattern through extended droughts and extreme precipitation storms, which in turn affect the earthquake cycles. The actual impacts of a warming climate are fault specific, but for the San-Andreas Fault, located in a Mediterranean climate zone, an enhanced hydrological cycle will further reduce the magnitude of earthquakes but enhance the occurrence frequency.
Existing geological data limited SEGMENT only represent major faults on orogens. Large number of minor faults [41,42] are not included in the input data, and their seismic activity in those locations are beyond the capability of our modeling system. New higher resolution geodetic information on plate interface and boundaries, sea floor spreading, fault systems, and geodesy all are necessary for further improvements of the predicting score, by improved description of continuum deformation and subsequent rupture/yielding.

Supplementary Materials
Due to size limitation, important details are provided in the Supplementary Material (SM). The author also will provide assistance if readers wish to re-produce the results.

Discussion
Attempts to understand earthquakes have a very long history and are globally widespread. The novel earthquake theory presented here follows from a confluence of several factors: Advances in remote sensing technology; recent, tectonically active stages motivated many researchers to re-evaluate previous assumptions [33]. Most critically, a well-documented landslide numerical model is used as a prototype of an advanced numerical modeling system to represent the mechanics of inter-plate earthquakes, and to verify the original hypotheses proposed here. This study shows, through data analysis and modeling, that variations in groundwater loading and occurrence of drought conditions together influence the occurrence and properties of earthquakes. Otherwise, the stationary seismicity of neotectonics is very stable for the majority of major faults mapped on existing neotectonic models. The temporal scale separation, i.e., the geological conditions vary at a much longer time scale than the occurrence of major earthquakes, which underscores why groundwater fluctuations constitute a viable mechanism for affecting earthquakes.
Earthquakes occur when the repose angle determined by plate-coupling stress reaches the sum of two angles representing, respectively, the fault geometry and fault strength. Accurate representation of the evolution of fault strength and repose angle therefore is critical for earthquake simulation and prediction. The earthquake triad, as defined in this study, is a generic framework that has universal applicability. Groundwater fluctuations, although relatively small, nevertheless play a pivotal role in linking the three triad components together and thereby determining earthquake irregularity. Groundwater influences fault instability directly by super-imposing a seismogenetic lateral stress field working in synergy with plate coupling stress, and also indirectly, by weakening the fault from within by enhancing fault fatigue from GM generation. The direct effect alone, the groundwater fluctuation superimposed stress field is of the same order as the residual of gravity-aided frictional stress and plate coupling stress; hence it affects the timing of major earthquakes. These findings agree with recent studies by Amos, et al. [6] and Gonzalez, et al. [7], and present the results in a coherent theoretical framework. The present study addresses the concerns of Avouac [40] and is a significant advance in understanding earthquake dynamics and thermodynamics for reliable earthquake prediction.
By exploring existing records of earthquakes using SEG-MENT, it is found that for collisional systems, direct effects from groundwater control the irregularity of major earthquakes. Essentially, the orographically highly variable precipitation is effectively channeled into deeper depth by the prevalence of through-cut faults. Droughts elsewhere also are seismogenetic but likely their effects are indirect. For example, the Elnino-southern oscillation (ENSO), as the dominant climate driver of regional precipitation, has a strong influence on the gravity field over the Andes region of South America. The occurrence of major earthquakes, although bearing the same 4-7 year occurrence frequency of ENSO, has a significant hiatus when compared with  Groundwater is not singled out as a subsection because it is the nexus between the components, and is discussed as an integrated component throughout the text.

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, the generic term "bridging effects" (Figure S2b; Ref. [20]) is assigned to this phenomenon. It is a disturbance of the vertical resistive stress field (R. Z ). The "G" in Eqs. (1) and (A1) is a special case of R. Z when the bridging effect is not salient. In principle, disturbances at any depth above the fault zone can be felt at the plate interface. It is noted that, for a point source (of horizontal dimension much less than the fault dimension), the influence range of the bridging effect in each horizontal direction expands linearly. So, the footprint increases as the square of the 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 it is in superposing onto R. Z and contributing to the motion tendency of this segment of the shear zone. According to this reasoning, localized load disturbances such as human-made infrastructure (e.g., dams) 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 a cavity at depth H depends on rock tensile strength (f c ) and density (ρ), following a square root law: where b is a cavity geometry dependent constant. It is apparent that larger dimension cavities are allowed at greater depths in the brittle crust. At greater depth (>~ 30 km), the thermodynamic reduction in the tensile strength of rock dominates and causes a drastic reduction in L. As a result, the cavities on the locked interface are sub-meter scale at most [44,45]. Bridging effects for regions surrounding sea mounts and these deep-burying cavities are non-negligible for an earthquake evolution and life cycle. These features have clear remote-sensing signatures, especially when experiencing transition in size as a result of collapses (usually a consequence of continued GM formation), or through an underground chute system channeling groundwater into deeper depths. Load fluctuations that are shallower, affect fault stability only when they occur on sufficiently large spatial scales. In this sense, regional scale groundwater fluctuations, even those residing in the upper few kilometers of depth, can have a significant impact on fault stability. The deeper they are, the more direct and efficient they are in contributing to fault stability.
The bridging effect contributes to gross resistive stress

Supplementary Material (SM) of "How Droughts Influence Earthquakes": Gravity Footprints of Major Earthquakes Are An Earthquake Triad Tale
The SM contains a definition of stability factor (Section 1), details of the earthquake triad (Section 2), and the application of SEGMENT to the Andean and Cascadian faults (Section 3). Five illustrations (Figures) for generic numerical model verifications and grid stencil representation also are placed in this SM.

Stability Index and the Intertwined Earthquake Triad
Different segments of a realistic heterogeneous fault zone ( Figure S2a) have different slope angles and maximum friction coefficients and also experience different loads. The driving stress distributed to each segment approaches its maximum supportable resistive stress gradually, similar to the saturation process of porous media. Fault 'stress saturation' is 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 value of 1 corresponds to creeping. Hence, the overall stability of a realistic fault becomes Whatever the actual form of an earthquake model, Eq. (A2) is an important co-seismic diagnosis. One important factor causing non-uniform distribution of driving stress among segments is the existence of macroscopic crevasses of voids/ cavity within the plates [43]. The illustrations in Figure 2 assume uniform lithostatic loading among all segments. In reality, the bridging effects [20], that is, the load directly above a plate segment can be partially carried by neighboring segments (or vice-versa), thereby adding yet another layer of complexity. Bridging contributes to gross resistive stress only when the fault interface has non-uniform friction coefficients (or, non-uniform fault strength). Otherwise, if the interface has the same roughness, the net contribution to resistive stress from bridging effects may be negligible due to the canceling effects from neighboring regions that are relieved of part of their lithostatic stress. The following subsection outlines the complex interaction among the earthquake triads, either occurring in sequence or concurrently. The seismic cycle is controlled largely by the compounding of these factors.

Gravity Fluctuation as the Nexus of the Earthquake Triad
Fault geometry, although by far the largest factor in determining the timing of major earthquakes, is treated here as the geological background. The super-imposed gravity fluctuations are primarily responsible for earthquake irregularity. The following discussion focuses on the bridging effect and Two objects, when pressed against each other, suffer wear and tear at the interface, a manifestation of material fatigue. The interfaces of tectonic plates are no exception. One validation/proof is the existence of widespread fault gouges and pseudotachylyte melt in wall rocks [33]. The generation rate of GM results from at least three factors: (1) The degree of imbalance in the three principle stress components, in a form of yielding criteria for brittle material; (2) The existing amount of GM or the current state of GM; and (3) The formation of GM consumes elastic energy accumulated in the plates, by producing new surfaces under the huge confining pressure of the fault environment. Point (1) Is supported by recent findings of frictional melt along the fault surface; it occurs only at a small fraction of the total rupture surface area, but always begins at asperities, the effectively strongest points on the fault surface. Cyclical forcing is most effective in causing fatigue. Groundwater fluctuations, especially when acting in concert with the resonance frequency of the plate segment, aid GM generation. The second requirement is that, as GM builds up at the interface, it hinders further production of GM in exactly the same way saw dust hinders further cutting of a piece of wood. This is especially relevant for the inter-plate GM because, compared with its parent rocks, it has larger porosity and easily fills the limited inter-plate space. Earthquakes are very effective mechanisms for removing and re-distributing existent granular material, specifically by thermal pressurization or acoustic fluidization. Granular material production and redistribution create a negative spatial covariance between loading and effective frictional coefficients, thus are seismogenic (i.e., contributing to instability) according to Eq. (A1). Figure S2c is the retrieved GM thickness, using the scheme in Ref. [23] and surface GPS displacements as constraints. The co-seismic release of potential (elastic as well as gravitational) energy is transformed mainly into heat through inter-plate frictional force (the destructive waves consume only a portion of the total energy). This possibly is a healing mechanism counteracting the granular fatigue. A large, deeply buried area approaching saturation simultaneously, un-locking systematically and releasing the stored elastic energy coherently, usually produces large magnitude earthquakes [51]. The contact surface of actual tectonic plates can have sophisticated topography and may involve material heterogeneity ( Figure S2a). Therefore, it is difficult to predict the future rupture propagation without knowing the interface details. Recent studies [11,12,51] indicate that widespread and smooth contact between two relatively weaker plates harbors large earthquakes; this is 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 fraction of the plate interface, acting more effectively as a grinder for the overlain plate rather than as a blocker/obstacle ( Figure S2a; Ref. [51]). Downstream of these seamounts, there are strips of granular debris composed of material scraped off both the sea mounts and the overlying continental plates. These findings further support the importance of lateral coherence between different segments of the fault in determining earthquake magnitudes (Eq. A1).
only when the plate interface has non-uniform friction coefficients. If the interface has the same roughness µ', the net contribution to resistive stress from bridging effects may be minuscule due to the canceling effects from neighboring regions that are relieved of part of their lithostatic stress. 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 is negative). Granular material [13,14] formation and transportation is such a nexus. Through bridging effects, a neighbor's weight is partially transferred over regions with strong contacts (similar to the "pillars" of a bridge). These regions of strong contact encourage GM formation (a consensus of yielding criteria of brittle rocks; Ref. [43]). With the formation of granular material with viscosity (< 10 4 Pa•S), many orders of magnitude smaller than polycrystalline rocks (> 10 21 Pa•S), the shear zone with strong contacts weakens. Such a situation is seismogenic. For example, if there is only 5 mm GM between rock plates under a confining pressure similar to the locked fault zone (~ 0.7 GPa), the equivalent frictional coefficient, or fault strength, drops to < 10 -4 , orders of magnitude smaller than solid-against-solid rock friction.
At the fault interface, cavities are rare. The 'pillars' are surrounded by rock materials. Sea mounts are a type of inter-plate "pillars" that grinds the upper plate and in turn suffers wear. The bridging effect therefore 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. The near surface displacements (directly associated with strain buildup and release) can now be accurately measured by Global Positioning System (GPS; Refs. [46-48]) and Synthetic Aperture Radar interferometry (InSAR; e.g., Refs. [49,50] and references therein) instruments. From Eq. (A1), it contributes to stability if more weight is loaded on rough contacts and less weight is loaded on smooth contacts along the shear zone. Bridging effects, by affecting GM production rates, influence the frictional properties of the fault interface and achieve a seismogenic effect caused by more weight being loaded on smooth contacts and less on rough contacts along the shear zone. Sources of bridging effects are numerous. For example, large scale structures on the overriding plate (e.g., mountains and valleys) signify bridging effects at depth on the seismogenic fault zone [15]. However, bridging effects caused by topographic features, and hence variations in the frictional properties of the seismogenic zone, cannot explain the irregularity of earthquakes that occur on decadal to centennial scales. However, large scale variation in groundwater is a viable candidate.
From Eq. (A1), the weight of the overriding plate always serves as a stabilizing factor, that is, its reduction is seismogenic. Fluctuations in the gravity field not only have local stability consequence, they also affect neighboring regions through bridging effects and the unique fatigue process of the fault zone.

Fatigue in rock interface-generation of granular material tion Faults
The following sections discuss earthquakes over two selected regions addressed in detail only in the SM, namely, the Andes and the Cascadian subduction faults, as specific applications of the numerical modelling systems.
Andes earthquakes: Controlled jointly by granular material pattern and groundwater fluctuation Figure S4 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 counter clockwise 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. Reference [54] suggests that the subduction of oceanic lithosphere coupled with underplating and a brief episode of gravity spreading contributed to crustal thickening in the back arc of the Central Andes. There remains no uniform theory that explains how very thick crusts develop and evolve (e.g., Ref. [55]).

Footprints of seamounts-role of granular material in
Andes earthquakes spatial pattern: Granular accumulation along the Peru coast (South of Lima) extends far inland. 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 are oriented at an angle to the Nazca Ridge's moving direction (as shown on Figure S5). This is similar to chiseling a rotating disk. The GM 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 GM to the west (A') ages as old as 14 Ma, whereas GM 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.2 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 the Nazca ridge at present but was very different 3Ma ago [56,57]. This chain of seamounts plowed beneath Puna at ~ 10-6 Ma ago like a hockey stick. At present, the 'head' of the stick had melted into the upper mantle and only the 'handle' is left. However, the GM 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 20ᵒS; 70ᵒW). 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 well with the granular stripes. The granular stripes left by the Juan Fernandez Ridge foster larger (but In addition to the critical role it plays in GM generation, groundwater, as an additional loading on the overriding plate is a stabilizing factor, according to Eq. (A1). 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 delicate balance. Fluctuation of groundwater, by superimposing a lateral/horizontal stress gradient, can play a critical role in the timing of major earthquakes. Fortuitously, fluctuations in groundwater have been remotely sensed in recent decades by gravity satellites such as the Gravity Recovery and Climate Experiment (GRACE) [17][18][19].
Centered on the co-seismic criterion (Eq. (1)), the triad determining earthquake irregularities are discussed. It is clear that, because of the heterogeneity in the 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 (energy released over all ruptures) 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 results from the inter-plate GM 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 fact, except for groundwater fluctuations, few other natural cyclic factors possess spatial wave lengths and temporal frequencies of the same orders as major tectonic earthquake irregularities. The climate fluctuations, by affecting groundwater loading patterns, impact the earthquake cycles. Bridging effects by themselves might not cause fault weakening. But, with the associated GM generation, they weaken faults through the formation of a negative spatial covariance pattern between loads and friction coefficients. The generation and transportation of GM at the plate contact reconcile the above-mentioned, counter-intuitive, characters of major earthquakes (e.g., as enumerated in Ref. [33]) and also are critical in orchestrating 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 [11] fields. The continental collisions (continental against continental; e.g., Himalayas) are not differentiated here 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 major, unsolved scientific problem [52,53]. In this study, first order schemes representing the earthquake triad are implemented in a sophisticated 3D modeling system, referred to as SEGMENT, to estimate the spatio-temporal scales of major earthquakes.

Applications of SEGMENT over Two Subduc-
droughts and the seismogenic effects (it is not a one to one correspondence: fatigue from multiple cycles of droughts leads to a reduced inter-seismic period). Rather than a direct temporal correspondence, there is a multiple year hiatus before the seismogenic effects from droughts to take effect. For example, the extended Californian drought may contribute to fault instability 2 (the local Riverside and Parkfield area to southern side Mendocino fault transition zones, or the Oregon-Washington coast) to some 10 years (northern Brish Columbia coast) later. In addition, the Cascadia fault's locking is distinctive in that the oceanic plate still is creeping, which is seen clearly by examining the velocity profile. It is locked in the sense that the neighboring regions have large gradients in motion speeds, signifying a stress build-up. For example, the Olympic National Park coast of the Cascadia subduction zone is now inter-seismically locked (Figure not shown for clarity). At 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 × 1013 J/yr. The estimated apparent frictional coefficient is 0.12. If there is no perturbation of the loading field, the natural frequency of seismic cycle is ~ 600 year. By releasing the accumulated elastic energy, an earthquake of Mw 7.8 would be produced (of energy of ~ 45.7 PJ). However, from SEGMENT simulations, it is apparent that the Californian drought exerts a cyclical stress on the interface and enhanced its fatigue. If the same hydrological cycle over the past 50 years was repeated for (hence is representative of) the past 800 years, the recurrence frequency would be reduced to ~ 446 ± 70 years. As a result, a > Mw 7.5 (~34 PJ) earthquake is expected within the upcoming decade. From the GM distribution pattern, the rupture likely would start from the southern sector, that is, between the Oregon coast and the coast of British Columbia.
The logic of relating California droughts to future Cascadia earthquakes might appear to be counter-intuitive as Cascadia has many of its own hydrological processes (e.g., it has 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 the Cascadia region of Canada does not have a large inter-annual to inter-decadal variation, unlike the California region which is influenced by both ENSO and the PDO climate drivers. More importantly, the Mediterranean climate is highly sensitive to climate warming [58]. Local Canadian precipitation fluctuation would be more direct but the magnitude of fluctuations truly is smaller than that of Californian precipitation. The North American plate cannot maintain full rigidity over such a distance. The waveguide for propagating the bridging effect zig-zags ~5.5 wave lengths along the Pacific coastal ranges before reaching the Cascadian fault. Thus, the fatigue caused by groundwater fluctuations associated with increased precipitation volatility in twenty-first century California likely will play a decisive role in shortening the timescale of the natural earthquake cycle and, at the same time, reducing the magnitude, for the entire coastal region (from Mexico to the south to Canada to the north). less frequent) earthquakes. For both seamount ridges, their left (northern) flanks foster larger earthquakes, because of the granular accumulation patterns.
Groundwater's indirect control (with apparent hiatus) of Andes earthquakes' spatio-temporal patterns: For the Peru and Chile coasts, ENSO ocean-atmospheric teleconnection patterns play a critical role in determining the phase of groundwater. Earthquakes, especially those greater than Mw 5.5, but less than Mw 7, are not necessarily occur in phase with the ENSO. This means that, unlike for the TP region, the direct mechanical effect, as that indicated by Eq. (1), may not be the dominant factor. However, the occurrence frequency of the major earthquakes is ~ 4-7 years, about the same as ENSO. A model simulation confirmed that it is the fluctuations, rather than the actual weight of groundwater that are critical for the GM generation and the ensuing earthquakes. Hence, for this region, groundwater-induced loading fluctuation affects earthquakes through granular material generation and evolution (subsections 2.1 and 2.2 of the SM and Section 2.3.3), rather than through the direct mechanical mechanism.
Despite the apparent temporal hiatus for most earthquakes between MW 5.5 and MW 7, all major earthquakes greater than Mw 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. ENSO seemingly lessened the severity of the earthquakes over the Andes. Without ENSO, all the scattered < Mw 7.5 earthquakes would 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 a strong correlation with seismogenic zone segmentation [15]. However, it simply is an illusion that surface topographic feature corresponds to ruptures at depth. Instead, the precipitation and the ensuing groundwater distribution that 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 groundwater fluctuations therefore is a practical means of understanding earthquake cycles. Thus, gravity mapping satellites such as GRACE provide essential assistance in predicting and verifying the occurrence of these earthquakes.

Cascadia subduction zone's stability affected by Californian droughts
Detailed discussion of the Cascadia subduction fault, such as its effective strength (µ'), geometrical characteristics and unique geological conditions are given in Ref. [3]. Here the additional stress fields super-imposed by precipitation fluctuations are discussed. Californian droughts (for example, that of 2012-2016), even though they were hundreds kilometers away, significantly enhanced the instability of the Cascadia subduction zone. Due to the lack of vertically extended crevasses, such as those resulting from the India-Eurasia collisional system, groundwater fluctuations take effect through aiding GM formation. Therefore there is a hiatus between Figure S3: SEGMENT diagnosed time scales of plate resonance frequency to groundwater fluctuation forcing. Background map was obtained from http://sideshow.jpl.nasa.gov/mbh/series.html and has been cleared for public release. Figure S4: The present plates interface geometry of the coastal Andes. a) Panel is crust thickness distribution; b) Panel is a vertical crosssection 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; c) Panel 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 in terms of the compressive stress between continental and oceanic plates. In particular, note that ocean-continental subduction here does not lead to thickening of the continental plate. Geometry data obtained from CRUST1.0 [21].