- Split View
-
Views
-
Cite
Cite
G. Hillers, S. A. Miller, Dilatancy controlled spatiotemporal slip evolution of a sealed fault with spatial variations of the pore pressure, Geophysical Journal International, Volume 168, Issue 1, January 2007, Pages 431–445, https://doi.org/10.1111/j.1365-246X.2006.03167.x
- Share Icon Share
Summary
A range of observations suggest the formation and maintenance of sealed and hence overpressured compartments in fluid-infiltrated fault zones. It is assumed that hydromechanical properties of regions with variable pore pressure states control the fault's stability and thus its characteristic response, that is, seismic or aseismic slip accumulation. We investigate in a systematic parameter space study the effects of spatial variations in pore pressure on spatiotemporal slip evolution along a hydraulically isolated fault plane. The 3-D continuum model is governed by rate-and-state friction and constitutive laws for porosity reduction. We show that the model response is sensitive to the degree of overpressurization and the efficiency of dilatant hardening mechanisms. Low pore pressures and small dilatancy effects result in unstable response types, whereas high pore pressures and large dilatant effects lead to stable and aseismic creep. Regions with an unstable response are shown to support most of the stresses accumulated during interseismic periods. Accelerated slip nucleates preferably in regions of low pore pressure. Statistical properties of model seismicity produce a wide range of event sizes for moderate and large earthquakes, in the case where dilatant mechanisms are inefficient. In case of efficient slip rate controlled porosity increase, less instabilities grow into large earthquakes. Final slip maps demonstrate the applicability of the chosen method to model seismicity controlled by frictional and hydraulic processes on a planar fault plane. The evolution of governing variables that depend on the pore pressure environment provide a conceptual basis for the interpretation of observed response characteristics.
1 Introduction
The presence of fluids in the Earth's crust and its mechanical effects on the seismogenic process is a well-documented but still not well-understood phenomenon because of feedbacks among several possible hydromechanical mechanisms that are difficult to separate. Numerous studies identified the potential importance of fluid migration following a main shock producing aftershock sequences (e.g. Nur & Booker 1972; Bosl & Nur 2002; Miller 2004; Piombo 2005, and references therein). Free fluids are expected to be involved in seismic swarm activity (Waite & Smith 2002), and they remotely triggered earthquakes in geothermal and volcanic areas (Husen 2004). Stabilizing as well as weakening fluid-related mechanisms during rapid seismic slip, such as dilatant hardening and pore fluid expansion due to shear heating, respectively, are assumed to alter rupture propagation in large earthquakes (Garagash & Rudnicki 2003a, b). These competing mechanisms act on different time and velocity scales, that is, hardening effects influence nucleation phases associated with small slip rates, whereas thermal pressurization has a negligible effect during initial stages but becomes important when the sliding velocity is high (Andrews 2002). The dominating effect of thermohydraulical processes during seismic slip was shown to control the evolution of the friction coefficient and thus the rupture propagation in recent dynamic numerical experiments (Bizzarri & Cocco 2006a, b). The observed rapid slip-weakening behaviour substantially controls the imposed friction law which is necessary to model observed seismograms (Andrews 2006).
A number of laboratory and field observations confirm the formation and maintenance of low-permeability seals between the fault gouge and the host rock (Sibson 1994) that effectively trap pore fluids and lead to the formation of a conduit in the fault's core (Moore 1994; Zhang & Tullis 1998; Zhang 1999, 2001). Porosity reducing mechanisms that can lead to overpressured fluid states are plastic pore closure, stress-induced dissolution and crack healing and sealing (Walder & Nur 1984; Nur & Walder 1992; Sleep & Blanpied 1992; Blanpied 1998; Lockner & Byerlee 1994). Alternatively, near-lithostatic fluid regimes can result from continued upwelling of overpressured fluids from the bottom of the seismogenic part of the crust (e.g. Rice 1992), or dehydration of hydrous minerals from frictional heating. Hydraulically isolated overpressured fault zones remain sealed indefinitely even though occasional fault valve behaviour (Sibson 1992; Cox 1995)—breaking the seals during large slip events—may temporarily allow the equilibration of pore pressure states in-and outside the fault zone (Blanpied 1992).
The San Andreas fault's strength is found to be anomalously low, based on the lack of a significant aberrant heat flow signal accompanied by borehole stress measurements indicating a steep angle between the fault's strike and the maximum compressive stress (Lachenbruch & Sass 1980; Hickman 1991). An explanation that satisfies both constraints are elevated pore pressures reducing sufficiently the frictional strength (Rice 1992). Heterogeneous fault zone materials with different physical properties exposed to pore compacting mechanisms is a likely cause of different degrees of overpressurization in a fault zone. Focusing on this mechanism Lockner & Byerlee (1995) conducted a numerical study of a 1-D Burridge–Knopoff type model (Burridge & Knopoff 1967), where most of the shear stress is supported by a small number of compartments with relatively low pore pressure. Miller (1996) extended this concept to a 2-D fault plane with a simple model for fluid pressure variation within the fault zone. In this study we adopt these concepts (Fig. 1) and investigate, in a systematic parameter space study of a 3-D continuous fault model, the effect of variable fluid pressures trapped in a sealed fault zone on spatiotemporal slip evolution.
Previous 3-D hydromechanical strike-slip fault models investigated the properties of spatiotemporal seismicity controlled by fluid-related mechanisms, such as the development of sealed pockets of high pore pressure and permeability changes during the seismic cycle (Miller 1996, 1999; Fitzenz & Miller 2001). These models demonstrated the simplicity of generating and maintaining complexity along a fluid-infiltrated fault zone by including dominant hydromechanical effects. Their inherent discreteness, however, led to generic results such as slip events with physically plausible properties that make it difficult to estimate the relative importance of fluid-related processes, since purely elastic models of the discrete class have been shown to generate complex slip evolution and seismicity over a wide range of sizes without fluid interaction (e.g. Burridge & Knopoff 1967; Langer 1996; Carlson & Langer 1989a, b; Carlson 1991; Bak 1987, 1988; Bak & Tang 1989; Ito & Matsuzaki 1990; Lomnitz-Adler 1993; Ben-Zion 1996; Zöller 2005a, b; Dahmen 1998). The quasi-dynamic continuous approach presented in this study aims to investigate qualitatively the nucleation, propagation and arrest of slip instabilities in a sealed, and thus overpressured, fluid-infiltrated fault in a consistent framework. The modelling emphasis is the long-term slip evolution and different response characteristics of different regions on the fault, stable or unstable, which are controlled by pore pressure dependent hydromechanical interactions. Although the quasi-dynamic approach produces a solution of the coseismic evolution of the governing variables, we recognize that neglecting inertia, the normal stress dependence of the state variable, and temperature dependence may limit the accuracy of the details of the rupture process. Nevertheless, much can be learned about the interaction of hydromechanical effects and additional mechanism can be investigated in future studies.
We use the 3-D elastohydraulic model of Hillers & Miller (2006) who investigated the origin of slip complexity due to dilatancy as a function of permeability structures in the off-fault dimension using homogeneous properties along strike. The use of rate-and-state friction allows multicycle simulations, and results are compared to purely elastic models to extract the contribution of hydraulic parameters. The approach is an extension of a 1-D spring-block model discussed by Segall & Rice (1995) to study the efficiency of dilatant hardening mechanisms to suppress unstable sliding. Hillers & Miller (2006) showed that homogeneous frictional and hydraulic properties can lead to non-uniform, dilatancy controlled slip histories for a realistic range of parameters. In this paper, we expand this work by using heterogeneous distributions of variable pore pressure regimes along the fault, ranging from hydrostatic to near lithostatic, which have been shown to describe realistic properties of weak, mature fault zones. We show that it is important to model the physically probable spatial variations in pore pressure along strike with an appropriate 3-D model to study the interaction of different pore pressure regimes on a fault. Because of the dominant role of the spatial heterogeneity of overpressure on the fault's stability, the results are discussed in the context of complexity generating processes to identify the controlling mechanisms.
The remainder of the paper is organized as follows. In Section 2, we discuss the conceptual model and its implementation in the numerical model, followed by an introduction of basic parameter setting (Section 3). Section 4 describes spatiotemporal slip evolution as a function of the degree of overpressurization and dilatant mechanisms. We introduce a procedure to extract physical quantities from our simulations in Section 5, before we discuss statistical properties and hypocentre locations of model seismicity in Sections 6 and 7, respectively. Finally, we focus on the temporal evolution of variables depending on the degree of overpressurization (Section 8).
2 Model
2.1 Conceptual model
Values of the permeability κ range from 10−19 m2 to 10−14 m2, even across a single fault line (Wibberly & Shimamoto 2003). Low permeabilities that lead to essentially undrained faults as considered here are on the order of 10−19 m2, prohibiting effectively any fluid flow between the fault's core and the crustal bulk. The parameter β= 5 × 10−4 MPa−1 is equal to the porosity, φ, times the sum of fluid- and elastic pore compressibility, respectively, and is assumed to be constant (Segall & Rice 1995), although the porosity is time dependent. The stability of the system is thus controlled by a constant effective hydraulic diffusivity, and by a slip rate induced pore compaction term discussed below. Likely interactions and feedback mechanisms of temperature changes, fault valve related fluid flow, crack formation in the off-fault dimension during rapid slip, and chemical cementation processes, which are usually difficult to separate, are not parametrized in this study. Assuming a viscosity of water at 60°C, ν≈ 5 × 10−4 Pa s, and a thickness of the cemented border, LD, of the order of 10 m (50 m) (Wibberly & Shimamoto 2003), c* is of the order of 10−2 yr−1 (5 × 10−4 yr−1).
To isolate the effects of variable degrees of overpressurization we assume the hydraulic diffusivity, c*, to be constant and homogeneous across the fault plane. The equilibration of p and p∞ is restricted using undrained conditions (c*→ 0) to account for sealed faults that develop spatially variable pore pressure regimes in their core, which can vary from hydrostatic (λ≈ 0.35) to near-lithostatic (λ≈ 0.9) conditions. We approximate heterogeneous pore pressure states by using 2-D λ distributions, and the pore pressure changes with respect to these initial conditions by slip rate controlled pore compaction. Hence, the degree of heterogeneity is quenched and its evolution and maintenance is not modelled explicitly by physical processes. Because the quasi-dynamic approach in tandem with heterogeneous hydraulic properties along the fault is shown to be limited in producing realistic coseismic slip rates, we follow Hillers & Miller (2006) interpreting intermittent accelerations as unstable, seismic slip events, even if the associated slip rates are much lower than natural velocities.
2.2 Modelling procedure
We solve the resulting set of six first-order ordinary differential equations in the variables using an implicit Runge–Kutta method for stiff systems with adaptive step-size control, RADAU5 (Hairer & Wanner 1996). We use the fast Fourier transform (FFT) to compute the convolution of the Green's functions matrix K with the velocity field, v∞−v (Rice 1993; Stuart & Tullis 1995; Rice & Ben-Zion 1996). The model faults are 100 km long and 24 km deep, and discretized into 512 × 128 numerical cells along strike and depth, respectively. Details of the numerical procedure and specific solver parameters can be found in the appendix of Hillers (2006).
3 Parameter Setting
3.1 Critical nucleation size
3.2 Stability regimes
Fig. 4(b) shows an alternative horizon of the multidimensional phase space, λ=λ(ɛ) (eq. 16 in Hillers & Miller 2006). A relatively large value of ɛ > 10−4 stabilize the system's slip evolution, whereas smaller ɛ values representing less-effective dilatant mechanisms leave the unstable stick-slip response unaffected. Since we use heterogeneous 2-D λ distributions with λ∈[0.35 0.9] in tandem with variable values of ɛ across the undrained fault, we anticipate different response characteristics of patches with different degrees of overpressurization (Hillers & Miller 2006). This shows that results depend not only on λ but also on ɛμ0/β, because this ratio of competing mechanisms of dilatancy and compressibility of the pore space controls the stability of overpressured compartments (Taylor & Rice 1998; Hillers & Miller 2006). That is, Segall & Rice (1995) showed that ɛ/β(b−a), weighted by the steady state friction coefficient, is the smallest value of σe for instabilities to develop. Unstable slip occurs only if a sufficient increase of the available pore space, φ, outperforms incompressibility, leading to a drop in pore pressure. Since μ0= 0.7 and β= 5 × 10−4 MPa−1 are constant throughout this study, we discuss slip evolutions in terms of the parameter ɛ.
4 Pore Pressure and Dilatancy Controlled Slip Evolution
In terms of the Lockner & Byerlee (1995) concept (Fig. 1), regions with low pore pressures, supporting most of the fault's stress in interevent times, fail developing accelerated slip rates. The slip event is terminated in neighboured regions of high-p values, because there unstable slip cannot occur because of effective dilatant hardening mechanisms. We observe a time delay between slip events at lower-p regions (Fig. 5a, x=[12, 63] km), with interceding creep in adjacent high-p regions. However, this pattern is not regular because sometimes the stable patch at x= 12 km slips prior to the one located at x= 63 km (e.g. u= 5–7 m), whereas in later sequences this patch becomes unstable first (e.g. u= 14–16 m). It is evident that heterogeneous λ distributions produce a persistent pattern of earthquake nucleation and variability in slip evolution that corresponds to the associated degree of overpressurization. Similar to Hillers (2006a), where earthquakes nucleate in regions of low L values, the creeping and locked parts in the present simulations depend linearly on the dominating degree of overpressurization. This is in contrast to Liu & Rice (2005), where small variations in the frictional parameters a and a−b along strike lead to non-linear and hence unpredictable slip pattern, and to Hillers & Miller (2006), where dilatant effects on a homogeneously parametrized fault produce non-stationary slip patterns.
Implicitly assumed seals bounding neighbouring patches against each other might break during seismic slip episodes and lead to equilibrated pore pressure regimes within the fault zone (Miller 1996). These boundaries possibly reseal in interseismic periods and different degrees of compacting mechanisms or direct fluid sources at depth are assumed to lead to different degrees of overpressurization prior to the next slip event. Hence, a more realistic treatment of these mechanism known as ‘fault valve behaviour’ (Sibson 1992; Cox 1995) might add further degrees of complexity to seismicity evolution, as well as the explicit modelling of the influence of off-fault damage, shear heating, and chemical alteration of the material's properties.
5 Earthquake Parameters
We generate seven different 2-D λ distributions with different degrees of heterogeneity, using 4 × 1, 8 × 2 and 16 × 4 patches along strike and depth, respectively (Fig. 6). To compare seismicity generated by this set of simulations we determine quantities that are listed in typical earthquake catalogues. We extract a seismic catalogue from the continuously simulated slip velocities generated by our numerical experiments using criteria for a seismic event introduced by Hillers (2006a). The algorithm measures the cumulative slip of a compact zone of cells as long as its slip rates are larger than a certain velocity threshold, separating seismic from aseismic slip. For details on the definition and interpretation of stable, aseismic response in contrast to unstable, accelerated slip in the present quasi-dynamic framework with limited spatial and temporal resolution see Hillers & Miller (2006). The event size is measured by the scalar potency P (sum of seismic slip times rupture area in [km2 cm]) (Ben-Zion 2003), and the corresponding event magnitude is obtained by an empirical scaling relation of Ben-Zion & Zhu (2002) for events larger than ML= 3.5, log10(P) = 1.34 ML− 5.22. Hypocentres are identified as the cell whose sliding velocity first exceeds the threshold slip rate at the onset of slip instability.
6 Statistical Properties
6.1 Frequency—magnitude distribution (FMD)
We perform two simulations for each of the seven λ distributions with ɛ= 5 × 10−5 and ɛ= 10−4, respectively. We stack the resulting catalogues extracted with the procedure explained in Section 5 and plot the number of hypocentres as a function of λ at the hypocentre location and the FMD (Fig. 7). The histograms show that earthquakes nucleate preferably in regions where λ is low. Differences between Figs 7(a) and (c) reveal that the localization of hypocentres at low λ is more pronounced for larger ɛ, which has already been implied by Fig. 5. The FMDs confirm the trend that different values of ɛ produce in general the same behaviour, but results differ in detail. If dilatant effects are relatively small we observe a clear break in scaling at ML= 5.5. For larger magnitudes, the slope is comparable to the reference slope of a frequency size statistics characterizing global strike-slip events shallower than 50 km (Frohlich & Davis 1993). Earthquakes with magnitudes smaller than ML= 5.5 occur far less often than an extrapolation of the slope to small magnitudes would imply. Thus, heterogeneous p distributions are sufficient to produce realistic FMD statistics only for moderate and large earthquakes.
Small slip events with magnitudes ML < 5 are produced less frequently compared to continuous elastic models with heterogeneous distributions of the critical slip distance L (Hillers 2006a, b), and to discrete models parametrizing structural heterogeneity (Ben-Zion 1996; Ben-Zion 2003; Zöller 2005a, b). This is partly due to the relatively large value of L used in the present study, which is an order of magnitude larger than in Hillers (2006a) and Lapusta (2000), resulting in an increased nucleation zone. Future studies using smaller L values are expected to generate seismicity with more realistic statistical properties over the entire magnitude range. Moreover, the linear range of λ is physically limited, compared to the logarithmic range of L in Hillers (2006a) and Hillers (2006b), indicating that the range of size scales is an effective tuning parameter (Ben-Zion 1996; Zöller 2005b). For ɛ= 10−4, the break in scaling in the FMD is less pronounced but shows a more gentle curvature of the FMD. It appears that only one ML≥ 6.5 event in the whole set of simulations occurred, revealing the ability of increased dilatancy mechanisms to stabilize fault slip more effectively (Fig. 5c; Segall & Rice 1995; Taylor & Rice 1998; Hillers & Miller 2006).
Note that the overall productivity of the present implementation cannot be compared to the generation of model seismicity from previous 3-D discrete elastic (e.g. Ben-Zion 1996; Zöller 2005b) and elastohydraulic models (e.g. Miller 1999; Fitzenz & Miller 2001), and continuous elastic models with heterogeneous frictional properties (Hillers 2006a, b). This is partly due to the numerically challenging implicit integration scheme, where the chosen spatial discretization of 512 × 128 cells along-strike and downdip allows the calculation of the solution only for a limited time interval. Furthermore, the present tuning parameter λ is relatively inefficient to produce seismicity with statistical properties similar to natural seismicity, compared to the effect of strong geometrical heterogeneities studied by Ben-Zion (2003), Zöller (2005b) and Hillers (2006a).
6.2 Synthetic slip maps
We show a set of final slip maps (Fig. 8) generated by our quasi-dynamic simulations using pore pressure regimes with different degrees of heterogeneity across a fault. The slip distributions of specific events show the same characteristics as those generated by the rate-and-state friction model with spatial variations of the critical slip distance, although they lack some of the persistent second-order features such as multiple high slip areas (Hillers 2006a). This proves that non-uniform pore pressure distribution on a fault are capable of producing heterogeneous slip maps of dynamic instabilities with properties comparable to examples of real (ML > 6) earthquakes (Mai 2004). However, modelled quasi-dynamic coseismic slip (umax≈ 2 m for ML7) generally tends to underrepresent slip values observed in large events (umax > 5 m, MW 7.2 Landers earthquake). Note that synthetic coseismic slip is larger in case the entire fault slips unstable (Hillers & Miller 2006). The visual similarity suggests the applicability of the chosen approach to study hydraulic mechanisms responsible for observed features of natural seismicity. Slip maps generated in a previous hydromechanical study of a 2-D fluid-infiltrated fault (Miller 2002) are strongly influenced by the governing discrete character of the numerical implementation and the binary static/dynamic treatment of the frictional resistance, whereas in the present model the numerical grid does not control the solution. Hence, slip maps presented here are direct consequences of the interaction between regions of different degrees of overpressurization.
7 Hypocentre Locations
7.1 Testing procedure
Following Hillers (2006b) we test hypocentre locations (HL) generated by various simulations against 2-D distributions of physical parameters prevailing at the HL to investigate the dependence of earthquake nucleation on spatially varying hydromechanical properties of the fault zone. First we compare the HL to the original 2-D λ distribution, normalized by the underlying distribution of λ values. The second parameter is the effective nucleation size, hsTR, based on the estimate of Taylor & Rice (1998), hnTR, which is derived using eq. (9) with σe(x, z), a(z) and b(z). See Hillers (2006b) for the derivation of hs. Parameters are normalized to the range [0 1] to be comparable among each other.
7.2 Hypocentre locations
Fig. 9 displays the normalized histograms of hypocentre locations as a function of normalized parameter range of λ and hsTR. Clearly, most slip instabilities nucleate in regions where λ and hence the pore pressure is low, as outlined in Fig. 7. For both model classes with ɛ= 5 × 10−5 and ɛ= 10−4, respectively, the maximum occurrence of HL can be found at smallest values, followed by an approximate 1/λ decay toward larger values. The same behaviour is observed by Hillers (2006a) and Hillers (2006b) in a number of models, where most earthquakes tend to nucleate in regions where the heterogeneity parameter L is small, followed by a similar decay. Hence, both models produce most earthquakes in regions that are most likely to become unstable because of smallest nucleation sizes. This is confirmed in the correlation of HL with hsTR, showing a pronounced localization at small values. Within the framework of the present study, and assuming the present ɛ values to be realistic, these observations imply that stress sustaining portions during interseismic periods can be interpreted with regions having a low pore pressure state that have the tendency to develop intermittent instabilities. Those low-p regions may be labelled strong, opposite to weak high-p regions. Note that this correspondence is justified in the framework of the present model implementation. In terms of the present nomenclature Rice (1993), Rice & Ben-Zion (1996), Lapusta (2000) and Hillers (2006a) used an effective normal stress at seismogenic depth of σe= 50 MPa, leading to a rather high λ and thus weak fault zone. Ben-Zion & Rice (1995) used heterogeneous but static pore pressure profiles to explore its effects on slip evolution. However, because these elastic models do not explicitly simulate the evolution of the pore space and associated pore pressure changes, and its effect on fault zone stability, an implication about the role of fluids cannot be deduced from their results.
A persistent observation in the present study is the negative correlation between earthquake nucleation and high pore pressures. However, the identification of high-p regions as weak zones seems to imply just the opposite. It is important to note that the weak zones tend to release tectonic load by aseismic, steady stable sliding (Fig. 5c), whereas strong low-p zones fail seismically, developing intermittent instabilities. This response pattern is an original outcome of the present parametrization, leading to stable response modes due to effective dilatant mechanisms in the undrained limit. Possible rapid in-plane fluid flow that might trigger earthquakes due to a sudden drop of σe is not considered here.
8 Evolution of Variables
We observe pore pressure controlled response functions along the fault. Fig. 10 shows stress and normalized velocity evolutions, τ=τ(u) and v=v(u), respectively, along a fault with a hydrostatic pore pressure regime. The model contains a high-p compartment deliberately placed slightly offset from the centre of the fault to investigate its effect on rupture spreading and the associated evolution of governing variables τ and v. Hydraulic barriers, that is, stably sliding overpressured regions, are imposed around x= 0 km (recall the use of periodic boundary conditions), such that the arrangement is intended to trigger instabilities around x≈ 75 km. However, the non-linearity of the system's response and its resulting slip evolution makes it difficult to trigger an event at a certain a priori specified place. In fact, very few instabilities emerge that develop a rupture front that passes the high-p patch at the fault's centre. Rupture speed and slip rates do not reach realistic values, hence the result is intended to illustrate qualitatively the effect of the overpressured patch. For example, slip velocity ranges only four orders of magnitude, whereas in fully dynamic treatments the difference in slip rates between locked and slipping episodes is as large as 10 orders of magnitude (Lapusta 2000). The origin of the abscissa, u= 0 m, corresponds to the slip horizon at the beginning of the instability, which is not necessarily balanced along the fault (Fig. 5). During the course of the system's spatiotemporal evolution, an instability nucleates around cells 280–310 where p is hydrostatic. The shape of the resulting velocity function (Fig. 10, cells 286, 306) show a gradual rise, and a distinct maximum before it drops to its interseismic regime. The corresponding stress curves reach their maximum prior to the maximum slip rates, and drop to the frictional sliding level, in agreement with kinematic and fully dynamic yet dry investigations (Tinti 2004; Piatanesi 2004).
At x > 125 km the instability soon dies out because of the stabilizing effect of the high-p patch centred at the model's origin. This behaviour might be a particular feature of the quasi-dynamic treatment where realistic slip rates of the order of a meter per second are not reached. A fully dynamic solution that considers heating effects is likely to produce different stopping mechanisms. We can isolate a relatively mature rupture front that passes the high-p patch left of the hypocentre. At the passage, the stress pulse remains significantly smaller compared to the surrounding low-p regions, since the compartment's response is stabilized by dilatant mechanisms and hence cannot support large shear stresses. The velocity, however, is of the same magnitude as the slip rates developed left of the patch. This indicates that the velocity evolution of the overpressured region does not differ qualitatively from the hydrostatic regions, but associated quasi-dynamic stress levels are significantly lower. At cells 186–246 the event produces prominent stress pulses. Here, too, the peak stresses are reached prior to the maximum slip rates, suggesting an initial slip hardening phase and a delayed initiation of dynamic fault weakening. Hence, acceleration is not accompanied instantaneously by weakening mechanisms, as suggested by Ohnaka (1996), Ohnaka & Shen (1999), and Ohnaka (2003). A similar observation has been made by Bizzarri (2001) in a purely elastic drained model, suggesting this effect is not caused by the explicit consideration of dilatant mechanisms, but by the use of rate-and-state friction. Except in the nucleation zone and the high-p patch, the maximum slip rates—and corresponding stress pulses—are reached towards the end of the slip event, and are caused by reflected slip waves at the instability's boundaries. To the left of cell 186, the slip event ceases approaching the stable barrier. The overall complexity and heterogeneity of the curves among the observation points is a consequence of the evolutionary character of the current approach, different to the aseptic and stylized approach of simulations of a single dynamic rupture with homogeneous properties and forced nucleation. Therefore, details of the response functions like the relative position of curve peaks are due to the precise state of the system prior to the instability. Nevertheless, the clear difference between the responses of the high- and low-p regions is of general validity.
We extend the specific analysis of this tuned example and present response functions obtained at computational points at seismogenic depth for a typical example of a heterogeneous 2-D λ distribution that consists of 8 × 2 patches (Fig. 11). We plot responses from two locations in each λ patch, indicating that they essentially coincide but are nevertheless influenced by the response of adjacent patches (cells 80/112 and 336/368). The variable's histories display the evolution over several non-uniform cycles (Fig. 5), confirming that modelled coseismic slip (Fig. 8) is smaller than what is observed in real seismicity. Although limited in its resolution, this example serves as a conceptual basis for interpreting different response characteristics along a fault zone based on non-uniform hydraulic properties. The essence of this and other simulations (not shown) is that regions of low-p exhibit larger stress drops than regions of high-p, where exact results depend on specific choices of the dilatancy coefficient (Fig. 7). Small patches of high-p surrounded by larger areas of low-p show a significantly lower stress pulse when a rupture is passing (e.g. cell 400, see Fig. 10). From this bimodal (low/high-p) point of view, τ(u) functions shown in Fig. 11 give a broader overview over possible response types. The stress functions confirm the result by Lockner & Byerlee (1995) that areas of low-p have a higher stress level than neighbouring regions of high-p (e.g. cells 368, 400).
Corresponding slip rates reveal a dependence of maximum slip velocities on the degree of overpressurization. We find that patches with low λ develop slip rates that are up to two orders of magnitudes larger than regions with a large λ (e.g. cells 144, 400). However, the spatial λ heterogeneity connected to a relatively large cell size leads to slip rates of these high-p regions that are significantly larger than the background driving rate, v∞. The particular response of overpressured regions is governed by the efficiency of dilatant mechanisms, as demonstrated in Figs 5(a) and c, for example, at x= 70 km, where small values of ɛ allow v≫v∞, but v≈v∞ for ɛ being relatively large. That is, the continuum limit approach leads to accelerated slip rates of high-p regions if ɛ is sufficiently small, as these regions are passed by growing ruptures, which nucleate in low-p areas (Fig. 10). Hence, the modelling provides a physical explanation for slow slip events preceded and/or followed by seismic slip events (Segall 1996; Reilinger 2000; Hearn & Reilinger 2002).
Porosity evolution illustrates the velocity dependence of φ (eq. 8), since v(u) and φ(u) show corresponding pattern. Moreover, the inverse relation between porosity and pore pressure changes reveals the effect of pore compacting mechanisms. Differences in p(u) depending on λ are less significant compared to stress and velocity changes, demonstrating the efficiency of porosity changes in controlling the fluid pressure in sealed fault zones. During interseismic compacting periods p is above the background level, whereas coseismic increase in φ causes p to drop below the corresponding value of p∞ (Hillers & Miller 2006). Note that present changes in p are small compared to those generated by Lockner & Byerlee (1995), Miller (1996, 1999), and Fitzenz & Miller (2001) which is due to different pore compacting constitutive laws (Sleep & Blanpied 1992; Sleep 1995). Eq. (8) was found to produce no significant changes in p (Segall & Rice 1995), but within the present framework it has stability controlling properties parametrized by the dilatancy coefficient ɛ (Fig. 4). Future models that produce a more fertile fluid environment could include fluid sources at depth, or dehydration mechanisms.
9 Discussion and Conclusions
Extending the analysis of Segall & Rice (1995), our simulations of a 3-D, quasi-dynamic implementation of the conceptual model by Lockner & Byerlee (1995) and Miller (1996) demonstrate that instabilities nucleate preferably in low pore pressure regions, and that their extension is controlled by high pressure areas. In particular, effective dilatant mechanisms lead to aseismic, stable response functions of high-p patches, whereas less-effective dilatant hardening allows accelerated slip in these regions but no substantial stress changes. The degree of seismic coupling is controlled by the relative importance of dilatant processes parametrized by the dilatancy coefficient, which is not very well resolved based on laboratory estimates (Segall & Rice 1995; Sleep 2000). However, key features of our results remain unaffected by particular choices of ɛ, as long as large values (ɛ > 3 × 10−4) do not stabilize the fault independent of the degree of overpressurization. Hence, the response types depend not only on λ, but on the mechanical interaction of fluid pressure, dilatant effects, the counteracting degree of compressibility, and the (nominal) friction coefficient. Furthermore, the continuous character of the solution—indicated by a sufficient h/h* ratio—ensures that specific features of the seismic response are not controlled by effects of the spatial discretization, which can not be ruled out in previous discrete studies of hydromechanical interactions (Miller 1996, 1999; Fitzenz & Miller 2001). In addition, the rate-and-state formulation does not require an a priori defined failure threshold, which has shown to be a powerful method to consistently simulate earthquake cycles with realistic accelerating phases (Dieterich 1992; Rice 1993; Lapusta 2000; Rubin & Ampuero 2005).
We observe that instabilities clearly tend to nucleate in regions with low pore pressures, which seems to contradict observations of elevated fluid pressures at hypocentre locations (Zhao 1996). This discrepancy between modelling and observational evidence can be assigned to our time-independent treatment of the overpressurization pattern. This computationally necessary simplification does not rule out the possibility that a short-term fluid intrusion triggers seismicity (Miller 2004).
Although our model makes major but substantially supported simplifications it allows the efficient investigation of non-uniform slip responses based on likely hydromechanical properties assumed to control the seismic response of sealed and hence weak strike-slip fault zones. However, important processes have not been addressed in the present study that might alter a fluid-infiltrated fault's dynamic slip evolution significantly. Amongst these are fault valve behaviour allowing diffusion processes in the fault and leakage into the host rock (Blanpied 1992; Sibson 1994; Sleep 1995; Miller 1996), frictional heating (Lachenbruch & Sass 1980; Mase & Smith 1985; Spray 1995; Rice 2006; Bizzarri & Cocco 2006a, b) and slip-induced formation of new cracks (Scholz 1990).
The three main controlling mechanisms and properties that govern certain aspects of the observed response characteristics are summarized in Table 1. It becomes obvious that the effects of each of the single ingredients, rate-and-state friction, hydromechanics, and the role of heterogeneity, are studied in more detail elsewhere. However, the present approach combines these governing principles in a single, consistent 3-D model. Because of previously isolated investigations and applications of the three main mechanisms, we explored their contribution to the response characteristics presented here. In particular, rate-and-state friction allows the efficient investigation of multicycle simulations, hydromechanical processes control the stability of otherwise unstable regions, and imposed heterogeneity leads to seismicity over several orders of magnitude.
It has been shown that regions of near-hydrostatic pore pressure show the largest coseismic slip during unstable slip episodes, compared to regions of near-lithostatic fluid pressure. Corresponding stress evolution reveal that these low-p regions are expected to excite more seismic energy into the crust. Based on simulations presented here we suggest interpreting high-slip regions during earthquakes with relatively dry and hence strong fault sections. Due to the time-independent character of in-fault seals prohibiting pressure equilibration in the fault's core, these strong areas remain spatially immobile in our model. Heterogeneous but static pattern of pore pressure distributions produce non-uniform slip pattern, but the non-uniformity can be related to the governing λ distribution in that low-p regions always respond unstable and high-p regions stable. Thus, the chosen parametrization is sufficient to carefully investigate overpressurization controlled slip pattern with appropriate temporal and spatial resolution, but it lacks the ability to generate highly non-linear slip pattern produced in previous continuous studies (Liu & Rice 2005; Hillers 2006a, b).
As has been demonstrated, a relatively large L value and the linear range of physically plausible λ values and the associated variability in the effective nucleation size explain the lack of frequent occurrence of small events. Ben-Zion (1996), Zöller (2005b) and Hillers (2006b) showed that the range of size scales of the imposed heterogeneity parameter controls the statistics of the associated seismicity pattern. Large ranges of size scales are shown to result in FMDs that show a Gutenberg-Richter type power-law behaviour and are associated with geometric complexity of fault zones, whereas small ranges of sizes scales lead to a characteristic event type response pattern that correspond to homogeneous fault structures. These studies, and to a lesser degree the present investigation, demonstrate that a linear or log-linear distribution of controlling parameters on a planar model fault is sufficient to generate seismicity over a broad range of scales, in contrast to its often advocated correspondence to fractal parameter distributions (e.g. Scholz & Madelbrot 1989; Turcotte 1997; Dimri 2000; Hergarten 2002; Sornette 2004). FMDs presented here do not reproduce one of the two end-member cases related to fault zone maturity (Wesnousky 1994). However, future use of smaller L values with heterogeneous pore pressure distributions is expected to generate more smaller events due to smaller nucleation sizes. Together with a sufficient range of ɛ values it is anticipated to lead to a more realistic FMD, indicating that heterogeneous pore pressure states along a fault are related to geometrical complexity.
Acknowledgments
The manuscript benefited from comments by Massimo Cocco and an anonymous reviewer. The work was sponsored by EC-Project RELIEF (EVG1-CT-2002-00069). This is contribution number 1441 of the Institute of Geophysics, ETH Zurich.
References
Author notes
Now at: Institute for Crustal Studies, University of California, Santa Barbara, USA. E-mail: gregor@crustal.ucsb.edu