Numerical models for slip on the subduction interface motivated by field observations

Subduction interfaces show variations in plate-boundary behavior as a function of thermal structure, with changes in slip characteristics downdip along a single convergent margin, as well as global differences in seismicity between warm and cold subduction zones. Along an individual subduction interface, earthquakes are typically generated within a depth range associated with temperatures of 150–350 °C, with creep occurring updip and downdip of the seismogenic zone (Moore and Saffer, 2001). The size distribution of earthquakes is impacted by the age of the subducting slab, with younger subducting plates associated with a greater proportion of large earthquakes (Nishikawa and Ide, 2014). Some of the warmer accretionary margins with thick incoming sediment piles such as Sumatra and Cascadia experience supercycles defined by hundreds of years of subsidence separating sequences of earthquake clusters (Sieh et al., 2008; Goldfinger et al., 2013). Each rupture sequence evolves uniquely in terms of the order and grouping of asperities (Philibosian et al., 2017). The causes of cyclic earthquake clustering and temperature control on subduction interface slip behavior are enigmatic, but they may be related to periodic variability in the frictional properties of the interface (Mouslopoulou et al., 2016) through temperature-controlled healing. Here, we argue that geochemical processes play a role in modifying the behavior of the subduction interface through thermally activated healing of fractures and microslip surfaces accomplished by redistribution of silica during the time period between seismic ruptures (Fisher and Brantley, 2014; Saishu et al., 2017).

Evidence indicating that silica redistribution is an important controller of plate-boundary behavior comes from laboratory experiments that demonstrate the impact of silica kinetics on the rate-state friction characteristics of fault gouge material (Nakatani and Scholz, 2004; den Hartog et al., 2012; den Hartog and Spiers, 2013) and from observations of microstructures in ancient subduction-related fault zones (e.g., Fisher and Byrne, 1987; Fagereng et al., 2011; Kimura et al., 2012; Yamaguchi et al., 2012; Raimbourg et al., 2015). In this paper, we describe features that can be used to place constraints on a numerical model of the subduction plate interface. We show that many characteristics of subduction seismicity, such as earthquake clustering, an updip limit to the seismogenic zone (Fig. 1; Moore and Saffer, 2001), and a relationship between subducting plate age and earthquake size distributions (Nishikawa and Ide, 2014), can result from the failure of portions of the interface that nucleate and strengthen during interseismic periods at rates determined by silica kinetics.

Ancient accretionary complexes in the forearcs of active subduction zones expose paleosubduction interfaces that provide important insights into the deformation mechanisms and geochemical processes occurring along active plate boundaries. These ancient plate boundaries are characterized by regionally extensive zones of tectonic mélange (hundreds to thousands of kilometers along strike) juxtaposing accreted packages that decrease in age seaward, with consistent kinematic indicators of top-to-the-trench simple shear that occurs in the footwall of the subduction interface during underthrusting (Byrne, 1984; Fisher and Byrne, 1987; Taira et al., 1988; Kitamura and Kimura, 2012; Fagereng et al., 2011; Kimura et al., 2012; Yamaguchi et al., 2012; Hashimoto et al., 2012; Raimbourg et al., 2015). We base this discussion on examples of ancient fault zones that were accreted within a temperature range that spans the seismogenic zone (150–350 °C), for example, the Uyak mélange (Byrne and Fisher, 1990), the Ghost Rocks mélange (Vrolijk et al., 1988), the Mugi mélange (Matsumura et al., 2003), the Yokonami mélange (Hashimoto et al., 2012), the Okitsu mélange (Ikesawa et al., 2003), the Kure mélange (Mukoyoshi et al., 2006), the Central belt of the Kodiak Formation (Brantley et al., 1997; Fisher and Brantley, 2014), the Makimine mélange (Kitamura and Kimura, 2012), and the Hyuga mélange (Kondo et al., 2005; Raimbourg et al., 2015).

The block-in-matrix fabric records stratal disruption across a shear zone tens to hundreds of meters thick (Fig. 2A; Fisher and Byrne, 1987), with pervasive scaly fabrics in mudstone lithologies and penetrative systems of veins in sandstone blocks (Figs. 2B and 2C). Microstructures record the development of an anastomosing network of scaly microfaults in mudstones, with shear failure and subsequent dissolution, while stronger sandstone blocks experience pervasive opening-mode cracking and precipitation of quartz (with minor albite, calcite, and chlorite). Crystallization in veins involves periodic cracking and sealing to form solid (antitaxial) and fluid (syntaxial) inclusion bands (Fig. 2C; e.g., Fisher and Brantley, 2014). The observation of pervasive cracking at high angle to blocks is interpreted as evidence for a system that is buffered at fluid pressures that exceed σ₃ (Fisher, 1996). This footwall system experiences cyclic fluctuations in crack porosity (Fisher and Brantley, 1992). Healing of cracks occurs by local diffusive transport from scaly microfaults with high mean stress to open fractures with lower mean stress, as well as by precipitation from cooling fluids derived externally (e.g., Raimbourg et al., 2015) through migration along a fracture compartment parallel to the interface (Sibson, 2017).

The association of earthquake epicenters and areas of large moment release on megathrusts suggests that earthquake recurrence is controlled by asperities (Thatcher, 1990). For the purposes of the modeling, we treated asperities as patches along the interface that strengthen at rates determined by silica kinetics. These patches exhibit behavior similar to geodetic asperities capable of developing slip deficits in the interseismic period (e.g., Bürgmann et al., 2005), and seismic asperities capable of large coseismic moment release (Lay and Kanamori, 1982). The potential for portions of the interface to strengthen through mineral redistribution during the interseismic period forms the basis for a two-dimensional (2-D) numerical model for the behavior of the subduction interface, where aging follows a temperature-dependent rate law.

The plate interface is represented in our analysis by a 2-D block-slider cellular automaton (CA) model (e.g., Brown et al., 1991; Huang et al., 1992; Christensen and Olami, 1992; Wang, 2012). The model depicts a plate boundary as a critically stressed interface held in frictional equilibrium. Similar to classic CA models, our model produces wide ranges of earthquake magnitudes and clustered temporal sequences like those observed in nature (Huang et al., 1992) when run without any consideration of asperities. Under these conditions, the model fails to account for two important patterns that are apparent in natural earthquake data. First, the model produces magnitude distributions better fit by exponential or lognormal equations than power laws, as are commonly seen in nature (e.g., Gutenberg and Richter, 1944; Nishikawa and Ide, 2014). Second, although the model produces high time-frequency cyclicity, ruptures are not transmitted long distances in space. That is, a given cell will fail, together in cascading failures with its neighbors, with a near-characteristic periodicity, but there is no coordination of failures across the model domain. In contrast, natural earthquake histories (e.g., Philibosian et al., 2017) depict plate-scale megacycles of strain energy accumulation and release, in which large, infrequent earthquakes—which make up the “long tail” of a power-law magnitude distribution—release strain energy over large portions of the plate interface.

More recent numerical modeling of earthquakes has enhanced our ability to investigate the effects of rate state and velocity-dependent friction as well as radiation damping, viscoelastic stress transfer, and layered elasticity (Tullis et al., 2012a, 2012b). However, none of these models appears to produce the Gutenberg-Richter power-law scaling of natural earthquake sizes.

The aim of the modeling in this study was to incorporate the effect of healing through geochemically driven mineral redistribution along the interface. We revisited the CA model of Huang et al. (1992), modifying it to produce areas of greater strength upon the interface during plate motion. Our approach was to use Arrhenius equations for silica kinetics to nucleate and strengthen patches on the interface due to a transition analogous to the inferred variations in fracture porosity inherent in crack-seal behavior. The model is based on the special case of the modern frictional model where aging occurs on a stationary interface, and static friction is greater than dynamic friction (Rabinowicz, 1951), but the model does not include the supercritical crack growth, rate-state friction, and viscoelastic stress transfer included in some earthquake simulators (e.g., Richards-Dinger and Dietrich, 2012; Tullis et al., 2012b).

The motion of a given block alters the stress upon its neighbors, such that they might fail as a consequence of the neighboring failure. Upon each failure event, the model iteratively searches for consequent failures until none occurs, and an individual rupture is complete. After each rupture finishes, the downgoing plate continues to move steadily to the location at which the next rupture begins—that is, where the next block fails—and the rupture process begins anew. Ruptures are understood to propagate at seismic velocity, such that the passage of time in the model corresponds to the motion of the downgoing plate, and the ruptures are treated as instantaneous.

To model the effect temperature on nucleation and growth of asperities along subduction zones, we used T_min = 25 °C at the updip end and a default T_max = 250 °C at the downdip end; temperature does not vary along strike and is constant for each block throughout each simulation.

The model of Huang et al. (1992) is explicitly nondimensional. We implicitly introduced a size to the model space by assigning a temperature field. Our default scenario had a difference of 225 °C between the updip and downdip limits, which could represent a total distance of 130 km, using typical interface-parallel temperature gradients (Syracuse et al., 2010). Our default simulation used 100 × 100 cells, for a cell separation of 1.30 km, and a cell area of 1.69 km². However, the slip length in the model was determined only by an arbitrary stress value (Eq. 3) and so was independent of the block dimensions. Despite the fact that we established areal dimensions for the interface based on natural temperature gradients, there is no necessary scale to the slip unit, d, because S (stress) in Equation 3 has arbitrary units, and φ is a dimensionless ratio. However, in assigning a scale to d, we also constrained (1) the model time scale, by setting a constant rigid-plate velocity of 5 cm/yr, and (2) the magnitude of the ruptures, by assuming a shear modulus of 32 GPa and using the relation of Kanamori (1977). The seismic moment was calculated based on the area and average slip of a rupture and converted into moment magnitude, M_w. Scaling the model’s arbitrary slip unit (d = 1) to 4 m produced an earthquake magnitude distribution in which the largest sizes were similar to those of active subduction zones over times scales of hundreds of years or longer. By increasing the resolution of our model to 300 × 300 cells, we were able to produce earthquakes down to M_w ∼4; our default grid size, used to explore behavior in the model, precluded reproducing earthquakes smaller than M_w ∼5. To evaluate the earthquake size distributions, we ran simulations until the cell grid experienced 400,000 earthquakes.

For case A, without asperities (Figs. 4A and 4B), the resulting model output was generally better fit by an exponential distribution than a power law (Fig. 4B). In such a case, the b* value has limited meaning; nevertheless, the best-fit b* values were relatively high (Fig. 4A), reflecting the paucity of large events within an exponential distribution relative to a power law distribution. Varying the spring constant and friction ratios had no systematic effect on size distribution parameters.

For case B, in which asperities form randomly, regardless of temperature, we still found poorer fits to power-law distributions (Fig. 4D). Changing the maximum strength of asperities did not change this result. Lowering the static to dynamic friction ratio for asperities did reduce the b* value (Fig. 4C), presumably by increasing the displacement of slipped cells (Eq. 3), but it did not improve the quality of the power-law fit (Fig. 4D).

For a wide range of model parameters in case C, we found a better fit to a power-law distribution than to an exponential distribution (Figs. 4F, 4H, and 4J). Modifying asperity friction and strength produced qualitatively similar results to case B (Fig. 4E), only with a significantly improved power-law fit (Fig. 4F). We also tested the parameters that controlled the rate of asperity nucleation and strengthening, i.e., C and A_str (Figs. 4G and 4H) and E and G (Figs. 4I and 4J). In both comparisons, we observed that if either nucleation or strengthening was severely inhibited, the effect of asperities disappeared, and the model behavior was similar to that of case A. At intermediate nucleation and strengthening rates, a power-law size distribution emerged, with local minima in b* values, which approached the values of natural earthquake catalogues (∼1). Finally, when both nucleation and strengthening were rapid, the power-law quality of fit declined again, and b* increased. As we discuss below, it appears that the wide range of sizes inherent in the power-law distributions here emerged from having an updip cool region that experienced small earthquakes and a downdip warm region that experienced larger earthquakes produced by the rupturing of stronger asperities. The behavior of each of the three cases is illustrated by animations in the GSA Data Repository Item.¹

Figure 5 shows variations through time in the slip deficit along the interface. In all cases, the model exhibited minor edge effects within the three rows of blocks at the updip and downdip ends (Fig. 5). The presence of larger ruptures along the boundary did not significantly change the earthquake statistics, and we restrict our attention to the temporal and spatial distribution of earthquakes in the interior of the model space. In the absence of asperities (case A), the model produced a statistically significant cyclicity of strain energy buildup and release in time, which was local and not synchronized through time updip and downdip (Figs. 5A and 6). Addition of the nucleation and growth of asperities independent of temperature (case B) resulted in larger-magnitude events, with greater variations in slip deficit during cycles than the asperity-free model (Figs. 4B, 5B, and 6), but no clustering through time.

In the case of nucleation and strengthening of asperities as a function of the thermal structure of the interface (case C), the cooler updip portion of the interface (i.e., top of the model) behaved much like the asperity-free case, with high-frequency, low-magnitude events that are not apparent given the scale bar used in Figure 5C. However, there was a transition downdip from low-magnitude, high-frequency background failures updip to development of larger-magnitude failures downdip due to asperity nucleation (Fig. 5C). Three maps of slip deficit upon the plate interface are shown for different times in Figure 7; times are delineated in Figures 5C and 6. This transition occurred at ∼150–160 °C in the models and is analogous to the updip limit of the seismogenic zone. Increasing temperature within the lower part of the interface led to a greater density of asperity cells, and cycles with hundreds of years of elastic strain buildup punctuated by clusters of large failures. Updip propagation of large events that initiated in the areas of asperity formation caused little change in the slip deficit within the shallow region.

The size distribution of earthquakes that results from thermally activated asperity development follows a power-law relationship up to large magnitudes, unlike the no-asperity case (Fig. 4; e.g., Brown et al., 1991; Huang et al., 1992). Moreover, there is an increase in b* value on cumulative frequency plots of earthquake magnitude in the models for interfaces with colder thermal structure (Fig. 8A). The same relationship is observed in real-world subduction zones. To illustrate, we consider the Tonga, northern Chile, and Mexico subduction zones. These are, respectively, the coldest, the closest to the average, and the third hottest of the subduction zone segments modeled by Syracuse et al. (2010), as defined by the slab surface temperature at 30 km depth. The two hottest sections, both in Cascadia, have too few events for a meaningful cumulative frequency plot. The cumulative frequency plots for these subduction zones (Fig. 8B) show a trend of decreasing slope with increasing temperature, much like that seen in Figure 8A. Figure 8B was constructed using data from the Global Centroid Moment Tensor (GCMT) catalog, taking events within 250 km to either side of the cross sections of Syracuse et al. (2010) and within 25 km of the slab surface when projected onto the plane of the cross section. This filtering eliminated earthquakes within each transect that do not fall within the dipping band of hypocenters that define the plate boundary.

Earthquake size distributions, recurrence intervals of earthquake supercycles, and the slip distributions for these models of the subduction interface are a response to a population balance algorithm for asperities, with expressions for nucleation, strengthening, and failure, and with a temperature-dependent rate law for nucleation and strengthening. In the absence of asperities, the model interface experiences high-frequency fluctuations in slip deficit, with small-amplitude variations and no spatial pattern in the location or size of rupture events. We produced a significant change in behavior when we increased the proportion of the interface that becomes stronger during the interseismic period due to temperature-controlled nucleation and strengthening of asperities, which we attribute to healing of the fracture porosity. Buildup of slip deficit between earthquake clusters in the models is heterogeneous, much like estimates of plate coupling based on global positioning system (GPS) velocity fields in forearcs (Bürgmann et al., 2005).

We find that supercycles of earthquake clustering and development of a seismogenic zone with an updip limit are best reproduced when the nucleation and growth of asperities occur at a rate that increases exponentially with temperature, as expected based on Arrhenius-law silica kinetics. Under these circumstances, the recurrence interval and amplitude of fluctuations in the slip deficit depend on two factors: (1) the thermal structure of the interface, and (2) the difference in failure criteria between the asperity cells and other “unhealed” cells on the interface. The nonasperity condition is intended as an analogue for areas of the plate interface with relatively high footwall fracture porosity, whereas asperity cells are areas with diminished porosity through mineral precipitation.

Hotter subduction zones experience faster nucleation and strengthening of asperities that lead to a longer period for supercycles and a larger slip deficit and asperity area at the time of failure, in part because closely spaced asperities interact strongly and can behave as a single asperity (Johnson, 2010). The model in this study produces a creeping section updip of the hotter region of the interface characterized by cycles of buildup and release of elastic strain, similar to characterizations of the seismogenic zone. The sharp break in downdip behavior is also consistent with the slope break in forearc systems, and the inferred connection to the change from velocity strengthening to velocity weakening in models of wedge mechanics (Wang and Hu, 2006).

In our model, earthquake size distributions reflect a power law similar to the Gutenberg-Richter relationship on cumulative frequency plots. In the case of a fault surface with asperities, power-law distributions of earthquakes result from fractal asperity size distributions (Johnson and Nadeau, 2002). CA models without heterogeneity in cell failure criteria do not scale properly for the largest earthquake sizes, in part because in nature, supercritical crack growth allows rupture propagation above a given size at constant stress (Ben-Zion and Rice, 1993). In our models, the b* value, or slope of the power-law relationship, is smaller for warmer interfaces because of the impact of temperature on asperity nucleation. This model output is consistent with observations of b* values from natural subduction systems (Fig. 8B), where plate boundaries with younger and therefore hotter subducting plates have smaller b* values than older plates—a characteristic that has been explained as the effect of slab buoyancy on the normal stress along the interface (Nishikawa and Ide, 2014).

We propose an alternative explanation for variations in b* values: Warmer interfaces have accelerated rates of geochemical redistribution that result in faster healing. Asperities in our models form stochastically due to geochemical processes rather than representing geometric roughness elements on a plate boundary that lead to repeating events in the same location. A geochemical mechanism of healing is consistent with the observation that the largest earthquakes occur along sedimented margins with relatively smooth subduction interfaces (Scholl et al., 2015).

The power-law distributions produced by our CA model emerge when there is a wide range of asperity sizes and strengths. In such a scenario, a positive feedback emerges between cell asperitization and stress release upon rupture. If a cell is strengthened by a rise in effective normal stress and decrease in crack porosity along the interface, it stores more elastic strain energy with progressive plate motion. When such a cell finally does slip during an earthquake, it will undergo greater displacement (Eq. 3) and thus reach a more relaxed position relative to that of its unasperitized neighbors. This more relaxed position means it will take a longer time to reach its critical failure stress, and so it will have a greater likelihood to nucleate and be strengthened by mineral redistribution and undergo a reduction in fracture porosity. In this way, earthquake scaling for the plate interface is nonlinear, a characteristic shared among many natural systems that exhibit power-law behavior (Barabási and Albert, 1999).

Petrographic evidence from opening-mode fractures suggests that mineral precipitation introduces positive feedback between fracture size and susceptibility for future growth by preferentially sealing tiny fractures, preventing them from being reopened (Clark et al., 1995; Hooker et al., 2012, 2014). Hooker et al. (2014) also noted that if the cementation rate is faster than some threshold, the feedback loop is disrupted because all fractures become sealed, regardless of their size; consequently, large fractures are not able to outpace small fractures, and a characteristic fracture size distribution—as opposed to a power-law distribution—emerges. This case, where cementation keeps pace with crack opening, is analogous to the situation in the models when excessively fast nucleation and strengthening of asperities disrupt the development of a power law (Figs. 4H and 4J). When all cells become equally strong during an interseismic interval, no cells systematically outpace their neighbors in terms of coseismic displacement, and the system becomes linear.

The cell size of 1.69 km² in the default model dictates the minimum size of these modeled fault patches of higher effective normal stress, sets the minimum size earthquake that can be produced, and precludes us from considering variability in interface behavior at a smaller wavelength. In nature, this roughness in interface properties, tied to spatial variations in crack porosity and fluid pressure in the footwall of the plate interface, is likely due to heterogeneity in the sedimentary properties (i.e., grain size, sorting, cohesion, and sand-shale ratio) of the underthrusting material.

There are numerous factors important for plate-boundary behavior that are not considered in these models, such as variations in plate roughness (Wang and Bilek, 2011) and variations from margin to margin in the inputs that characterize the footwall side of the interface. Moreover, the geochemical reactions involved in mineral distribution are more complex than congruent dissolution and reprecipitation of silica, with incongruent reactions that occur on microfaults, where quartz in veins is one product, along with albite, calcite, and chlorite. We used a default value of 54 kJ/mol as the activation energy for silica redistribution, similar to estimates from experiments on H₂O-quartz systems (Rimstidt and Barnes, 1980; Lander et al., 2008), and in slide-hold-slide experiments for the healing of fault gouge (Nakatani and Scholz, 2004). Given the variability globally in subduction inputs, the potential differences in geochemical reactions, and the lack of knowledge about whether the nucleation and growth processes are rate limited by interfacial kinetics or diffusion, the activation energy may vary from margin to margin or even downdip along an individual margin, although the sensitivity analyses suggest a range of activation energies over which the temperature dependence of the model leads to a power-law size distribution. Our results indicate that a rate law with exponential temperature dependence for the development and hardening of asperities produces some of the first-order characteristics of shallow seismogenic behavior along the subduction interface. As temperature increases in the model, a greater proportion of the interface is stronger. In natural systems, there is a temperature further downdip along the subduction interface at which flow laws for ductile deformation allow steady deformation at a lower shear stress than the frictional strength.

The field observations used in this study are from sand-shale strata typical of thickly sedimented margins, and other lithologies could behave differently. Accretionary complexes are also incomplete recorders of subduction processes over tens of millions of years, with most exposed fault zones related to times of thick incoming sediment piles. Nevertheless, the models in this study demonstrate that many of the characteristics of subduction zone earthquakes, particularly along sedimented margins such as Sumatra and Cascadia, can be explained in the context of nucleation and growth of asperities as a thermally activated geochemical process. Over a temperature range typical of natural seismogenic zones, variations in effective mean stress and chemical potential along the plate interface lead to heterogenous frictional characteristics in the interseismic period that could control the location, recurrence time, and magnitude of earthquakes.

Fisher is supported by National Science Foundation grant EAR-1524530 from the Tectonics Program of the National Science Foundation. Hooker is supported as a fellow of the GDL Foundation. The manuscript benefited from discussions with Chris Marone and from the comments by Daniel Pastor-Galán and two anonymous reviewers. Fisher benefited from field trips to the Shimanto belt with Gaku Kimura, Asuka Yamaguchi, Rina Fukuchi, and Yoshi Hashimoto. Fisher contributed the field observations and conceptual model, Hooker contributed the cellular automaton numerical model, and Oakley contributed the analysis of earthquake size distributions.