The analysis of Coulomb stress changes has become an important tool for seismic hazard evaluation because such stress changes may trigger or delay subsequent earthquakes. Processes that can cause significant Coulomb stress changes include coseismic slip and transient postseismic processes such as poroelastic effects and viscoelastic relaxation. However, the combined influence of poroelastic effects and viscoelastic relaxation on co- and postseismic Coulomb stress changes has not been systematically studied so far. Here, we use three-dimensional finite-element models with arrays of normal and thrust faults to investigate how pore fluid pressure changes and viscoelastic relaxation overlap during the postseismic phase. In different experiments, we vary the permeability of the upper crust and the viscosity of the lower crust or lithospheric mantle while keeping the other parameters constant. In addition, we perform experiments in which we combine a high (low) permeability of the upper crust with a low (high) viscosity of the lower crust. Our results show that the coseismic (i.e., static) Coulomb stress changes are altered by the signal from poroelastic effects and viscoelastic relaxation during the first month after the earthquake. For sufficiently low viscosities, the Coulomb stress change patterns show a combined signal from poroelastic and viscoelastic effects already during the first postseismic year. For sufficiently low permeabilities, Coulomb stress changes induced by poroelastic effects overlap with the signals from viscoelastic relaxation and interseismic stress accumulation for decades. Our results imply that poroelastic and viscoelastic effects have a strong impact on postseismic Coulomb stress changes and should therefore be considered together when analyzing Coulomb stress transfer between faults.

Large earthquakes affect the stress state on faults in the vicinity of the source fault that ruptured during the seismic event in such a way that subsequent earthquakes on neighboring faults (called “receiver faults”) may be promoted or delayed. It is therefore crucial for seismic hazard evaluation to analyze the earthquake-induced stress changes, which are usually expressed in terms of Coulomb stress changes, ∆CFS (King et al., 1994; Stein, 2003; Freed, 2005; Scholz, 2010):

formula

where ∆τ is the change in shear stress (positive in direction of slip of the source fault), μ is the friction coefficient, ∆σn is the change in normal stress (positive if fault is clamped), and P is the pore pressure. An increase in the Coulomb stress implies that receiver faults are brought closer to failure, whereas a stress decrease indicates that the next earthquake may be delayed (King et al., 1994; Stein, 1999). Hence, the analysis of Coulomb stress changes can be used to identify possible locations of future earthquakes and has become an important tool for the evaluation of seismic hazard in a region (e.g., Stein, 2003; Toda et al., 2005; Freed et al., 2007; Parsons et al., 2008; Field et al., 2009; Luo and Liu, 2010; Wan and Shen, 2010; Serpelloni et al., 2012; Wang et al., 2014; Bagge et al., 2019).

Coulomb stress changes on receiver faults may be generated by a variety of earthquake-induced mechanisms, including static stress changes caused by the coseismic slip (King et al., 1994; Stein, 1999; Ryder et al., 2012) and dynamic stress changes caused by seismic waves (Belardinelli et al., 1999; Gomberg et al., 2001; Oglesby et al., 2003; Pollitz et al., 2012). Transient postseismic Coulomb stress changes may be caused by afterslip, pore fluid pressure changes, and postseismic viscoelastic relaxation. Afterslip may occur when unrelieved stress leads to aseismic slip on the fault plane (Freed, 2005). Pore fluid pressure changes are induced by the coseismic slip on the source fault and may trigger aftershocks by reducing the effective normal stress (Cocco et al., 2000; Piombo et al., 2005; Chiarabba et al., 2009; Chiaraluce, 2012; Malagnini et al., 2012). During the postseismic phase, coseismically induced pore fluid pressure gradients are relaxed by fluid flow if rocks are sufficiently permeable (e.g., Albano et al., 2017; Nespoli et al., 2018; Tung and Masterlark, 2018). Postseismic viscoelastic relaxation refers to the process by which the coseismically imposed stress is relaxed due to viscous flow in the lower crust and lithospheric mantle (Nur and Mavko, 1974). Depending on the viscosity of the lower crust and lithospheric mantle, postseismic viscoelastic relaxation may cause considerable stress changes on faults in the brittle upper crust (Pollitz, 1997; Freed and Lin, 1998; Kenner and Segall, 1999; Nostro et al., 2001; Masterlark and Wang, 2002; Gourmelen and Amelung, 2005; Bagge and Hampel, 2017; Verdecchia et al., 2018). In this study, we focus on the effects of pore fluid pressure changes and viscoelastic relaxation on postseismic Coulomb stress change evolution.

Since the recognition that the stress transfer between source and receiver faults may trigger earthquakes (King et al., 1994; Stein, 1999, 2003), static Coulomb stress changes have been routinely calculated after major earthquakes (e.g., Wang and Chen, 2001; Parsons et al., 2008; Wan and Shen, 2010; Serpelloni et al. 2012). In addition, static Coulomb stress changes were used to explain past earthquake sequences (e.g., Nostro et al., 1997; Stein et al., 1997; Hardebeck et al., 1998). Analyses solely based on static Coulomb stress changes neglect, however, that transient processes as well as interseismic strain accumulation may considerably affect the magnitude and distribution of Coulomb stress changes on receiver faults during the postseismic phase (e.g., Cocco et al., 2000; Freed and Lin, 2001; Masterlark and Wang, 2002; Pollitz and Sacks, 2002; Cianetti et al., 2005; Bagge and Hampel, 2017; Bagge et al., 2019; Verdecchia et al., 2018). So far, studies that have considered transient processes in Coulomb stress calculations typically have accounted for either poroelastic effects (e.g., Albano et al., 2017, 2019; Nespoli et al., 2018; Tung et al., 2018; Miao et al., 2021) or viscoelastic relaxation (e.g., Luo and Liu, 2010; Verdecchia et al., 2018; Bagge et al., 2019) based on the assumption that the two processes act on sufficiently different temporal and spatial scales. Only a limited number of studies have considered transient Coulomb stress changes arising from both poroelastic and viscoelastic effects (Freed and Lin, 2001; Masterlark and Wang, 2002; Barbot and Fialko, 2010; Tung and Masterlark, 2018).

The validity of the assumption that poroelastic effects and viscoelastic relaxation act on different time scales has been challenged by inversions of geodetic data from intra-continental earthquakes and numerical models, which have shown that viscoelastic relaxation may already occur in the early postseismic phase (e.g., Ryder et al., 2007, 2010; Barbot and Fialko, 2010; Hampel and Hetzel, 2015; Bagge and Hampel, 2017; Mandler et al., 2021) and hence on the time scale of poroelastic effects (Dempsey et al., 2013; Albano et al., 2017, 2019; Tung et al., 2018). This temporal overlap of poroelastic and viscoelastic effects is reflected in the postseismic velocity and stress fields (Peikert et al., 2022), but the consequences for Coulomb stress changes on receiver faults have not been systematically studied so far. A better understanding of how transient processes overlap and affect postseismic Coulomb stress changes is, however, crucial for evaluating the stress state on receiver faults because stress changes cannot be measured directly for natural faults. Instead, Coulomb stress changes are typically calculated from the geodetic records of surface displacements and therefore strongly depend on the resolution of the data, the characteristics of the specific earthquake, and the geometry of the fault that was ruptured during the earthquake. To achieve a better understanding of principal Coulomb stress changes and the role of non-transient and transient processes during the postseismic phase, conceptual models are needed that allow investigation of stress change patterns independent of a specific earthquake or fault. So far, however, only few studies have investigated such principal Coulomb stress changes (King et al., 1994; Nostro et al., 2001; Cocco and Rice, 2002; Bagge and Hampel, 2016, 2017), but none of them considered poroelastic effects and viscoelastic relaxation together.

Here, we use three-dimensional (3-D) finite-element models with a generalized setup of intra-continental normal and thrust fault arrays to gain insights into principal Coulomb stress changes arising from coseismic slip, poroelastic effects, postseismic viscoelastic relaxation, and interseismic stress accumulation. Due to this implementation of transient and non-transient processes in the models, we are able to evaluate their relative importance for the magnitude and distribution of the Coulomb stress changes on the receiver faults in the model fault array. We investigate the spatial and temporal scales of poroelastic and viscoelastic effects and the resulting impact on the Coulomb stress changes by conducting experiments with variable permeability and viscosity structure of the crust and the lithospheric mantle. Therefore, our models provide comprehensive insights into the role of the different processes that may generate postseismic Coulomb stress changes. Due to their generalized setup, our conceptual models also help in differentiating between principal stress change patterns and stress changes arising from the specific characteristics of an earthquake or fault, which will help in interpreting stress change patterns calculated for earthquakes on natural faults.

Principal Model Setup

The 3-D finite-element models of this study are created with the commercial software Abaqus (version 2018). All models consist of a 200 × 200–km–wide and 100-km-thick lithosphere, which is subdivided into an elastic upper crust, a viscoelastic lower crust, and a viscoelastic lithospheric mantle. The principal model setup of the normal and thrust fault reference models including the elastic parameters (Poisson’s ratio ν, Young’s modulus E), density ρ, viscosity η, and permeability K of the layers is shown in Figure 1. Viscoelastic behavior is implemented as linear, temperature-independent Maxwell viscoelasticity. Although this rheology represents a simplification of the actually depth-dependent and possibly nonlinear viscoelastic behavior of the lower crust and lithospheric mantle (e.g., Freed and Bürgmann, 2004; Ellis et al., 2006), the implementation of viscoelastic layers itself is a considerable advantage because other models commonly used the homogeneous, purely elastic half-space models based on Okada (1992). Furthermore, values for linear viscosity are derived from the postseismic deformation after a number of earthquakes (e.g., Masterlark and Wang, 2002; Hetland and Hager, 2003; Nishimura and Thatcher 2003; Gourmelen and Amelung, 2005; Ryder et al., 2014). The viscosity values applied in our models (Fig. 1; Table 1) reflect the range of viscosity estimates of continental lithosphere (Kaufmann and Amelung, 2000; Masterlark and Wang, 2002; Hetland and Hager, 2003; Nishimura and Thatcher, 2003; Gourmelen and Amelung, 2005; Burov and Watts, 2006; Klemperer, 2006; England et al., 2013; Ryder et al., 2014; Shi et al., 2015; Henriquet et al., 2019).

In Abaqus, the permeability is defined as the hydraulic conductivity k that can be calculated from k = K × ρfluid × g / ηfluid (where ρfluid is fluid density, 1000 kg/m3; g is acceleration due to gravity, 9.81 m/s2; and ηfluid is fluid viscosity, 998 × 10−6 kg/m/s). To evaluate the effect of a variable permeability on the model results, we use permeabilities (K) between 10−10 m2 and 10−16 m2 for the upper crust in different experiments (Fig. 1; Table 1). These values reflect the range of permeabilities derived for continental upper crust (Manning and Ingebritsen, 1999; Ingebritsen and Manning, 2010; Stober and Bucher, 2015). For the lower crust, we apply a value of 10−18 m2 (cf. Manning and Ingebritsen, 1999). For the lithospheric mantle, we apply a permeability of 10−35 m2 to simulate an impermeable layer (Tung et al., 2018). The initial pore pressure distribution in the models is hydrostatic; as a boundary condition, a pore pressure (Ppore) of 9.8 × 108 Pa is applied at the model bottom. We simulate the coupling between solid and fluid phase by applying the coupled pore fluid diffusion and stress analysis in Abaqus (2018). To specify the proportion between volumes of voids and solids in the medium and the volume of fluid trapped in the medium, Abaqus uses the void ratio, which is a few percent for crystalline basement rocks (e.g., Masterlark and Wang, 2002; Masterlark, 2003). In our models, we use a void ratio of 0.06 and a saturation of 1. No pore fluid flow occurs across the model boundaries, and the faults are treated as impermeable, i.e., fluid cannot flow across them (cf. Rudnicki, 1986; Dempsey et al., 2013; Albano et al., 2017, 2019; Peikert et al., 2022). The implementation of impermeable faults is supported by observations from natural faults and experimental results, which show that faults act as a barrier to fluid flow once an impermeable fault gouge layer has developed (Scholz, 1987; Parsons et al., 1999; Piombo et al., 2005; Ingebritsen and Manning, 2010). Gravity is included in all models as a body force. Isostatic effects are implemented by adding a lithostatic pressure (Plitho) of 3 × 109 Pa and an elastic foundation to the model bottom (cf. Hampel et al., 2019). The property of the elastic foundation represents an asthenosphere with a density of 3200 kg/m3. The model bottom is free to move in vertical and horizontal directions; model sides in the xz-plane are fixed in the y-direction. By applying a velocity boundary condition to the model sides in the yz-plane, the model domain is either extended or shortened to simulate the interseismic strain accumulation.

In the model center, a source fault (SF), which experiences the coseismic slip during the analysis, and 10 surrounding receiver faults (RF1–RF5 and RF7–RF11) are embedded in the upper crust as frictional contact interfaces between footwall and hanging wall of each fault (Fig. 1). All faults extend from the top to the bottom of the upper crust and are 40 km long, which is a typical length for faults that may be ruptured by Mw ~7 earthquakes (Wells and Coppersmith, 1994; for details, see Model Phases section below). Distances between the faults are ≥15 km in the x-direction and ≥5 km in the y-direction, following natural spatial configurations of faults in, for example, the Basin and Range province (western United States; Wesnousky et al., 2005), the Aegean region (Roberts and Michetti, 2004), and the foreland of the Tibetan Plateau (Meyer et al., 1998; Hetzel et al., 2004). For normal and thrust faults, typical dips of 60° and 30° are used, respectively. In order to capture the Coulomb stress changes in the surroundings of the source fault, four receiver faults are located in the footwall and hanging wall of the source fault (RF4, RF5, RF7, and RF8), two faults are located along strike of the source fault’s tips (RF2 and RF10), and four other faults are located outside of the immediate hanging wall and footwall of the source fault (RF1, RF3, RF9, and RF11). With this distribution of the receiver faults, our model fault array also reflects that natural receiver faults are commonly located along strike of the source fault or in its footwall or hanging wall (e.g., Nostro et al., 2001; Hetland and Hager, 2003; Lin et al., 2011; Albano et al., 2021).

Slip on the model faults is initiated by extending or shortening the model domain at a total rate of 6 mm/a and controlled by a Mohr-Coulomb criterion |τmax| = c + μσn, where τmax is the shear stress, c is cohesion (zero in the models), µ is the friction coefficient (0.6 in the models), and σn is the normal stress. Note that although the fault planes are predefined with a rectangular shape, the spatial slip distribution on the fault plane is not prescribed and develops freely in response to the extension or shortening of the model.

Model Phases

Each model run consists of three model phases. In the first model phase, viscoelastic behavior and pore fluid flow are activated and lithostatic and hydrostatic pressure distributions as well as isostatic equilibrium are established in the model. During the second model phase, slip on the faults is initiated by extension or shortening of the model in the x-direction (Fig. 1). The faults are unlocked and allowed to accumulate slip until they reach a constant slip rate, which depends primarily on the applied extension or shortening rate, the fault dip, and the viscosity structure of the lithosphere (cf. Hampel et al., 2010; Bagge and Hampel, 2017; Peikert et al., 2022). Once all faults achieve a constant slip rate, they are in equilibrium with the regional stress field that results from extension or shortening and the viscosity structure of the model. The earthquake cycle is simulated in the third model phase, which comprises the preseismic, coseismic, and postseismic phases. In the preseismic phase, all faults are locked and slip accumulation stops while extension or shortening of the model continues. At the beginning of the coseismic phase, the source fault in the center of model fault array is unlocked, which causes a sudden slip (i.e., model earthquake). All receiver faults remain locked during the coseismic phase. The coseismic slip on the source fault is controlled by the applied far-field extension or shortening rate, the rheological properties of the model, and the length of the preseismic phase (cf. Hampel and Hetzel, 2012). In this study, we have chosen a duration of the preseismic phase such that the source fault experiences a maximum coseismic slip of 2 m (Fig. 2A), which represents a typical value of a Mw ~7 earthquake on a 40-km-long fault (Wells and Coppersmith, 1994). Note that the source fault fully relaxes during the coseismic phase (cf. Ellis et al., 2006; Bagge et al., 2019), i.e., the stress drop is complete, and no afterslip occurs on the source fault during the subsequent postseismic phase. During the postseismic phase, all faults including the source fault are locked again while extension or shortening continues to simulate interseismic strain accumulation. Coseismic and postseismic Coulomb stress changes are calculated for all predefined faults from the Abaqus output of the normal and shear stress for frictional contact interfaces, which takes into account the regional stress field.

Coseismic Displacement and Coulomb Stress Changes

Figure 2 illustrates the model state at the end of the coseismic phase (cf. Hampel et al., 2013; Bagge and Hampel, 2016, 2017). Note that the coseismic slip, surface displacements, and Coulomb stress changes are similar in all models of our study (Table 1) because the permeability of the crust or the viscosity structure of the lithosphere have an only negligible influence on the coseismic deformation. In all models, the source fault experiences a total coseismic slip of ~2 m (Fig. 2A). The vertical surface displacements show coseismic footwall uplift and hanging-wall subsidence in the normal fault model and primarily hanging-wall uplift in the thrust fault model (Fig. 2B), in accordance with geological and geodetical observations from normal and thrust fault earthquakes in nature (e.g., Yu et al., 2001; Chen et al., 2006; Lin et al., 2009; Liu-Zeng et al., 2009; Cheloni et al., 2010; Yu et al., 2010; Serpelloni et al., 2012). The horizontal displacement fields in the x-direction indicate surface movements away from the normal fault and toward the thrust fault (Fig. 2B), which implies extension and shortening across the normal and thrust source faults, respectively. Within the footwall of the normal fault, shortening occurs because the magnitude of the horizontal displacements in the positive x-direction decreases with distance from the source fault (cf. Hampel and Hetzel, 2015; Peikert et al., 2022). In the hanging wall, only minor horizontal movements occur perpendicular to the strike of the fault because the coseismic subsidence of the hanging wall mainly induces horizontal movements and shortening parallel to the fault’s strike (cf. Hampel and Hetzel, 2015, their figure 2; Hampel et al., 2013). In contrast, the direction and magnitude of the horizontal movements in the thrust fault model are such that extension occurs within the hanging wall and footwall (cf. Hampel and Hetzel, 2015; Peikert et al., 2022). Both coseismic shortening in normal fault footwalls and coseismic extension near thrust faults are in accordance with geological and geodetic observations from natural earthquakes (e.g., Slemmons, 1957; Myers and Hamilton, 1964; King and Vita-Finzi, 1981; Philip and Meghraoui, 1983; Crone et al., 1987; Meghraoui et al., 1988; Yu et al., 2001; Chen et al., 2006; Lin et al., 2009; Liu-Zeng et al., 2009; Cheloni et al., 2010; Yu et al., 2010; Serpelloni et al., 2012). Figure 2C shows the distribution of the coseismic Coulomb stress changes on the source and receiver faults in our models (cf. Bagge and Hampel, 2016). For our model earthquake with a coseismic slip of ~2 m and a magnitude of Mw ~7, the coseismic stress drop is on the order of 20–30 MPa and hence within the range of coseismic stress drop derived for natural intra-plate earthquakes (Kanamori and Anderson, 1975; Hanks, 1977; Scholz, 2002). In both the normal and thrust fault models, the receiver faults located in the footwall and hanging wall of the source fault (RF4, RF5, RF7, and RF8) mostly show negative Coulomb stress changes. The largest positive Coulomb stress changes are obtained for receiver faults RF2 and RF10, which are located along strike of the source fault tips, as well as on RF5 located in the hanging wall of the source fault. Notably, the positive stress changes occur in the upper part of RF5 in the normal fault model but in the lower part of the fault plane in the thrust fault model. Mixed patterns with both positive and negative stress changes occur on receiver faults RF1, RF3, RF9, and RF11.

Experiments Conducted for Parameter Study

For our parameter study, we define three reference models for each fault array (Table 1). Reference models R1nf (normal fault array) and R1tf (thrust fault array) include both poroelastic effects and viscoelastic relaxation. The second type of reference model (R2nf and R2tf) considers poroelastic effects but no viscoelastic relaxation during the postseismic phase. Models R3nf and R3tf include postseismic viscoelastic relaxation but no poroelastic effects. For our sensitivity analyses, we first perform a series of experiments in which we vary one parameter a time, i.e., we vary either the permeability of the upper crust (P1nf–P3nf and P1tf–P3tf) or the viscosity of the lower crust (V1nf–V2nf and V1tf–V2tf) (Table 1). The second set of experiments consists of four end-member configurations for which we combined high and low permeability and viscosity values, respectively (Table 1; PV1nf–PV4nf and PV1tf–PV4tf).

In this section, the postseismic Coulomb stress changes for the different normal and thrust fault models are shown for different time points between the first months and 50th year after the earthquake. First, we present the results for the normal and thrust fault reference models (R1nf–R3nf and R1tf–R3tf) in Figure 3 and Figure 4, respectively. Afterwards, we show the results from models with variable permeability (Figs. 5 and 6), variable viscosity (Figs. 7 and 8), and end-member configurations (Figs. 9 and 10).

Reference Models

As a base for comparing the stress evolution in the different experiments of our parameter study, we first present the results from the normal and thrust fault reference models with both poroelastic effects and viscoelastic relaxation (R1nf, Fig. 3A; R1tf, Fig. 4A). In the first month after the earthquake, most normal and thrust faults show Coulomb stress changes that are inverse to the coseismic stress changes, i.e., areas with positive coseismic stress changes show negative postseismic stress changes and vice versa. In the normal fault model (Fig. 3A), most faults exhibit positive Coulomb stress changes, with values of as much as 7 MPa on the source fault and 0.15 MPa (RF1 and RF9) to 2.5 MPa (RF7) on the receiver faults. On the receiver faults located close to the source fault (RF2, RF5, RF7, and RF10), both positive and negative stress changes occur. In the thrust fault model (Fig. 4A), all faults except RF3 and RF11 show positive and negative stress changes along their lower parts of as much as 9 MPa on the source fault and 0.25 MPa (RF3 and RF11) to 4 MPa (RF5) on the receiver faults. In the following months, the magnitude of Coulomb stress changes decreases on all faults, and the spatial distribution of positive and negative stress changes gradually changes in both the normal and thrust fault reference models (Figs. 3A and 4A). Stress changes in the first year are generally dominated by the stress changes of the first month in all reference models. The normal fault model shows similar stress change magnitudes on the source fault and the receiver faults in its footwall (RF3, RF7, RF8, and RF11) but higher magnitudes on the other receiver faults (Fig. 3A). In the thrust fault model, the stress change magnitudes are similar to those of the first month, but the spatial distribution of negative stress changes is different on some faults compared to the first month (Fig. 4A). Until the second year after the earthquake, the Coulomb stress changes decrease by one to two orders of magnitude. Most normal faults show only positive stress changes between 0.2 MPa (source fault) and 0.025 (RF1 and RF9); only RF5 still exhibits positive and negative stress changes. In the thrust fault model, all receiver faults in the hanging wall (RF1, RF4, RF5, and RF9) and RF7 and RF8 in the footwall show zones of negative stress changes along the upper part of the fault plane. All other faults experience positive stress changes of as much as 0.25 MPa on the source fault, 0.025 MPa on RF1, RF3, RF9, and RF11, and 0.2 MPa on RF5. From the fifth year onward, the magnitude of the positive stress changes on the source fault slowly decreases in all reference models. Most receiver faults show a uniform Coulomb stress distribution with an average stress increase of 0.02 MPa (normal fault model) and 0.015 MPa (thrust fault model) until the 50th year. Exceptions are normal receiver fault RF4, which shows a non-uniform distribution of positive stress changes, and RF5, which exhibits positive and negative stress changes until the 50th year in both models. In the thrust fault model, all faults located parallel to the source fault in its hanging wall and footwall continue to show positive and negative stress changes.

The second set of reference models involves poroelastic effects but no viscoelastic relaxation (R2nf, Fig. 3B; R2tf,Fig. 4B). Until the fifth year after the earthquake, these models show the same evolution as the reference models with both processes. Afterwards, only positive Coulomb stress changes are observed on the receiver faults in both the normal and thrust fault model, which indicates that poroelastic effects have decayed and interseismic stress accumulation prevails.

The third set of reference models is computed with viscoelastic relaxation but without pore fluid pressure changes (R3nf, Fig. 3C; R3tf, Fig. 4C). The postseismic Coulomb stress changes are two to three orders of magnitude lower than in the reference models with both processes from the first month onward. In the normal fault model, only positive Coulomb stress changes occur on all faults in the first three months except for receiver fault RF5, which shows negative stress changes in its upper part. From the first year until the 50th year, the positive stress changes on all other receiver faults remain constant at a value of 0.02 MPa. In the thrust fault model, receiver faults RF1, RF2, RF3, RF9, RF10, and RF11 show positive stress changes with a maximum value of 0.002 MPa in the first three months and 0.02 MPa from the first to 50th year after the earthquake. RF4 exhibits a zone of negative stress changes in the first months after the earthquake, which turns into a zone of positive stress changes during the first year. RF5 shows both positive and negative stress changes over the model time.

To illustrate the postseismic surface deformation in the reference models with and without poroelastic effects and viscoelastic relaxation, Figures 3D and 4D show the vertical and horizontal velocity fields for all three normal and thrust fault reference models (R1nf–R3nf and R1tf–R3tf). In the first month, the normal fault models with poroelastic effects (R1nf and R2nf) show horizontal velocities of as much as 500 mm/a and shortening across the source fault, whereas extension occurs within the hanging wall and footwall. Both hanging wall and footwall of the source fault subside at rates of as much as 1000 mm/a. In the following months, the zones of horizontal shortening and subsidence shift toward receiver fault RF5. Both models show the same velocity patterns until the second year after the earthquake. Compared to the two reference normal fault models with poroelastic effects (R1nf and R2nf), R3nf (with viscoelastic relaxation only) shows a different evolution during the early postseismic phase, which is characterized by extension across the source fault, uplift of the footwall and subsidence of the hanging wall, and velocities that are as much as two orders of magnitude lower. From the fifth year onward, velocity patterns in reference models R1nf and R3nf become gradually similar, indicating that the influence of poroelastic effects decreases while the signal from viscoelastic relaxation increases. Regarding the horizontal and vertical velocity fields, the thrust fault reference models (Fig. 4D; R1tf–R3tf) generally show a similar evolution compared to the normal fault reference models, except that all movements occur in the opposite direction, i.e., zones of shortening or uplift in the normal fault model are zones of extension and subsidence in the thrust fault models and vice versa.

Models with Variable Permeability

For the first model set in our parameter study, we varied the permeability of the upper crust between 10−10 m2 and 10−16 m2 (cf. Manning and Ingebritsen, 1999; Ingebritsen and Manning, 2010; Stober and Bucher, 2015) while keeping the viscosity of the lower crust constant (P1nf–P3nf, Fig. 5; P1tf–P3tf, Fig. 6). In models with a permeability of 10−10 m2 (P1nf, Fig. 5A; P1tf, Fig. 6A), the stress distribution in the first month and first year after the earthquake resemble the stress distribution of the second month in the reference models, i.e., it shows positive and negative stress changes on all faults except on normal faults RF1, RF4, and RF9 and thrust fault RF8. The stress changes reach values of 4 MPa on the source fault and as much as 2.5 MPa on the receiver faults in both models. Already in the second month, the stress changes in P1nf and P1tf decrease by two orders of magnitude and the Coulomb stress patterns are similar to those of reference models R1nf and R1tf in the second year, in which only normal fault RF5 and thrust faults RF1, RF4, RF5, RF7, RF8, and RF9 show negative stress changes. In the thrust fault model, these zones turn into zones of stress increase in the third month, and new zones of negative stress changes appear on the lower part of RF4 and RF5. Normal receiver faults show an average stress increase of 0.02 MPa over the years and a uniform Coulomb stress distribution, except RF5, which also exhibits negative stress changes. In the thrust fault model, RF1–RF3 and RF9–RF11 experience a uniform stress increase of as much as 0.02 MPa, which remains constant over the following years. RF4 and RF5 still show a negative stress zone in the lower part of the fault planes, which completely disappears on RF4 and becomes smaller on RF5 until the 50th year.

In the normal and thrust fault models with a permeability of 10−14 m2 (P2nf, Fig. 5 B; P2tf, Fig. 6B), almost all receiver faults experience both positive and negative stress changes, except normal faults RF1 and RF9 and thrust faults RF3 and RF11, which show only positive stress changes. The highest stress changes occur on RF7 (2 MPa) and RF5 (2 MPa) and on the source fault, where values of 12 MPa and 8 MPa are reached in the normal and thrust fault model, respectively. The negative zone on normal fault RF4 disappears in the second month and the magnitudes of the stress changes decrease, while the stress changes do not considerably alter in the following months in both models. In the second year, normal faults RF2, RF5, RF7, and RF10 show both positive and negative stress changes of as much as 0.5 MPa, similar to the first year in R1nf. In the thrust fault model, the stress changes of the second year are similar to those of the reference model R1tf, with magnitudes of as much as 1.5 MPa. In the fifth year, the zones of negative stress changes become smaller on most receiver faults in both the normal and thrust fault models, with the exception of RF7, where a second zone of negative stress changes develops. From the 10th year onward, the stress change distribution in the normal fault model shows similar patterns but higher magnitudes compared to R1nf. In the thrust fault model, the stress changes in the 10th year resemble the stress changes of the second year in R1tf, but afterwards they evolve in a different way. Whereas a uniform stress change distribution develops on RF1–RF3 and RF9–RF11, negative stress changes can be found on RF5 in the 50th year and RF7 still shows higher values of 0.04 MPa in the 50th year.

For the third set of experiments, we used a low permeability of 10−16 m2 (P3nf, Fig. 5 C; P3tf, Fig. 6C). In the first month, these models show positive and negative stress changes of as much as 2 MPa on the normal source fault, as much as 6 MPa on the thrust source fault, between 0.2 MPa (RF7) and 0.01 MPa (RF1 and RF9) on the normal receiver faults, and between 0.25 (RF5) and 0.02 (RF3 and RF11) on the thrust receiver faults. After the first month, the negative stress changes on the outer receiver faults (RF1, RF3, RF9, and RF11) disappear while they become smaller on the other faults. Stress change magnitudes slightly decrease. In the first year after the earthquake, all faults parallel to the source fault experience positive and negative stress changes in the normal fault model, with the highest stress changes being found on RF7 (0.7 MPa). The thrust fault model shows positive and negative stress changes on faults parallel and along strike to the source fault, with the highest value occurring on RF5 (0.9 MPa). Over the years, the zones with negative stress changes become smaller and turn into positive stress changes, except for normal fault RF5 and thrust faults RF4 and RF5, which continue to show negative stress changes in the 50th year.

Figures 5D and 6D illustrate that the permeability of the upper crust strongly influences the horizontal and vertical surface velocities, in particular during the early postseismic phase. Models with a permeability of 10−10 m2 (P1nf and P1tf) show high velocities in the first month, which then decrease by two orders of magnitude in the second month. Both the normal and thrust fault models exhibit pronounced zones of extension and shortening, which disappear in the third month. In contrast, vertical and horizontal velocities are initially lower for permeabilities of 10−14 m2 (P2nf and P2tf) and 10−16 m2 (P3nf and P3tf) and then show a more gradual decrease over time. Shortening across the normal source fault and extension across the thrust source fault prevail until approximately the fifth postseismic year. As in the reference models with poroelastic effects, both the hanging wall and footwall of the source fault initially experience subsidence in the normal fault model and uplift in the thrust fault model. This pattern gradually changes to the interseismic distribution of subsidence and uplift, with the time scale of the change depending on the permeability of the upper crust.

Models with Variable Viscosity

For the second model set in our parameter study, we changed the viscosity of the lower crust to 1018 Pa·s (V1nf and V1tf) and 1022 Pa·s (V2nf and V2tf), respectively, while keeping the permeability of the upper crust constant (Figs. 7 and 8). In the normal fault model with a viscosity of 1018 Pa·s (V1nf, Fig. 7A), the Coulomb stress changes reach values of 5 MPa on the source fault, 1.7 MPa on receiver fault RF7, and only 0.1 MPa on RF1 and RF9 in the first month, i.e., the values are slightly lower than in the reference normal fault model (cf. Fig. 3). In contrast, the thrust fault model V1tf (Fig. 8A) shows stress magnitudes similar to those of the thrust fault reference model (cf. Fig. 4). In both the normal and thrust fault models, the stress changes evolve largely similarly to those of the reference models in the first two months. In the third month, normal faults RF1, RF9, and the source fault still show only positive stress changes, whereas all other faults experience both positive and negative stress changes. Compared to the normal fault reference model, the stress magnitudes are higher and the negative stress changes occur at different positions. In the thrust fault model, one or two zones of negative stress changes with higher stress magnitudes than in the reference model are found on the receiver faults in the second and third month. In the first year after the earthquake, the stress change magnitudes are generally higher on the source fault and on most of the receiver faults than in the reference models. Receiver faults RF1, RF4, RF9, RF5, and RF7 show negative stress changes in their lower or upper parts, respectively. All other receiver faults are characterized by positive stress changes. Until the second year, the magnitudes of the stress change decrease, but they remain one order of magnitude higher than in the reference model. Between the first and 10th year, the receiver faults evolve differently from those of the reference models with respect to the spatial distribution of the negative stress changes. On some receiver faults, the zones with stress decrease become larger and occur in the lower part and/or the upper part of the fault plane (e.g., RF1, RF4, and RF9). On other faults, new zones of negative stress changes develop with time (e.g., RF3, RF8, and RF11). In the 10th year, both positive and negative stress changes can be observed normal receiver faults RF1, RF3–RF5, RF7–RF9, and RF11 and on all thrust faults, including the source fault. In the 50th year, a uniform stress change distribution with an average stress increase of as much as 0.02 MPa is found on all faults.

In the models with a viscosity of 1022 Pa·s (V2nf, Fig. 7B; V2tf, Fig. 8B), the stress field patterns and magnitudes evolve similarly to those of the stress field of the reference models in the first years. From the fifth year onwards, all normal and thrust faults show a constant stress increase of 0.02 MPa until the 50th year.

The surface deformation (Figs. 7C7D and 8C8D) reveals how the viscosity influences the horizontal and vertical velocities for a constant permeability. In the first month, the surface velocity fields in both normal and thrust fault models with high and low viscosities are dominated by the poroelastic effects and resemble the stress changes of the first month in the reference models. In the second month, the velocity patterns start to deviate from each other because in the model with a low viscosity (1018 Pa·s; V1nf and V1tf), the signals from poroelastic effects and viscoelastic flow begin to overlap. The surface velocity patterns in the models with low viscosity remain different from those of the respective reference models until the 50th year. In contrast, poroelastic effects continue to dominate the surface velocity fields until the second year in the model with a high viscosity (1022 Pa·s; V2nf and V2tf). Afterwards, the signal from slow viscoelastic flow prevails.

End-Member Models with Variable Permeability and Viscosity

In the final model set, we combine a high (low) permeability of the upper crust with a low (high) viscosity of the lower crust to maximize or minimize the effects from the overlap of poroelastic effects and viscoelastic relaxation (PV1nf–PV4nf, Fig. 9; PV1tf–PV4tf, Fig. 10).

The Coulomb stress changes in the model PV1nf (Fig. 9A) and PV1tf (Fig. 10A) with a high permeability (10−10 m2) and a low viscosity (1018 Pa·s) reveal a dominance of poroelastic effects in the first month and a strong overlap of poroelastic effects and viscoelastic relaxation from the second postseismic month onward. From the first to the second month, the stress changes on the receiver faults of both the normal and thrust models decrease by one order of magnitude, and one or two zones of negative stress changes develop on most faults. From the first to the 50th year after the earthquake, the poroelastic effects decay and the Coulomb stress changes evolve similar to those of models V1nf and V1tf with a low viscosity (cf. Figs. 7 and 8) but with different magnitudes.

In models PV2nf (Fig. 9B) and PV2tf (Fig. 10B), we combined a low permeability (10−16 m2) with a low viscosity (1018 Pa·s). In the first month, both positive and negative stress changes occur on most normal receiver faults (except RF3 and RF11) as well as on thrust receiver faults RF7 and RF8 that are located in the source fault hanging wall. The stress patterns in both models change only slightly over the next months. Stress change magnitudes are lower than in models V1nf and V1tf in the first month, but the subsequent decrease in magnitudes is slower. In the first year in the normal fault model, RF8 and the receiver faults on the hanging wall and along strike of the source fault show positive and negative stress changes, with the magnitudes being similar to those of model V1nf on most faults and lower on RF7. The evolution of the stress change pattern is largely similar to that of model V1nf, but the magnitudes always remain higher. In the thrust fault model, all receiver faults except RF2 and RF10 experience both positive and negative stress changes in the first year, with magnitudes being lower compared to model V1tf. From the second year onward, the spatio-temporal evolution resembles that of model V1tf, but with higher magnitudes. In the 50th year, the source fault still shows zones of negative stress changes.

If a high permeability (10−10 m2) and a high viscosity (1022 Pa·s) are combined (PV3nf, Fig. 9C; PV3tf, Fig. 10C), the Coulomb stress changes evolve similar to those of models P1nf and P1tf and show the same stress change magnitudes in the first two months and in the first year after the earthquake. In the following months and decades, only positive Coulomb stress changes are observed on both normal and thrust faults, which compares to models V2nf and V2tf from the fifth year onward. The source faults still show higher values of 0.1 MPa and 0.2 MPa in the normal and thrust fault models, respectively, which decrease over the years, whereas the receiver faults exhibit a uniform distribution of positive Coulomb stress changes between 0.015 MPa and 0.025 MPa.

In models PV4nf (Fig. 9D) and PV4tf (Fig. 10D), which combine a low permeability (10−16 m2) with a high viscosity (1022 Pa·s), the Coulomb stress changes show a slow, gradual decrease over the 50 years after the earthquake. In the early postseismic phase, the stress evolution is largely similar to that of models P3nf and P3tf. In the normal fault model, zones of negative stress changes on the receiver faults parallel to the source fault turn into positive stress zones between the second and tenth year. In the thrust fault model, negative Coulomb stress changes slowly disappear over time; in the 50th year, only RF5 and the source fault still show negative stress changes in their lower parts.

Our models with different permeabilities and viscosities provide insights into principal Coulomb stress changes arising from coseismic slip, poroelastic effects, postseismic viscoelastic relaxation, and interseismic stress accumulation. In this section, we first evaluate the main findings of our models and the relative importance of poroelastic effects and postseismic viscoelastic relaxation for the generation and evolution of postseismic Coulomb stress changes. We also link the postseismic Coulomb stress changes to the postseismic surface deformation. Afterwards, we compare our modeled Coulomb stress change patterns with stress change analyses for natural faults and earthquakes.

Relative Importance of Poroelastic Effects and Viscoelastic Relaxation

The model results show that the static Coulomb stress changes are significantly altered through space and time by both poroelastic effects and viscoelastic relaxation in the first months after the earthquake. As our parameter study reveals, the distribution, magnitude, and evolution of the postseismic Coulomb stress changes are controlled by the permeability of the upper crust, the viscosity of the lower crust, as well as the position of the receiver faults relative to the source fault. In the following, we discuss the transient postseismic stress changes separately from the coseismic stress changes to evaluate the relative importance of poroelastic effects versus viscoelastic relaxation. It should be noted, however, that for stress triggering, both coseismic and postseismic stress changes are relevant.

In the first months after the earthquake, models including poroelastic effects with (R1nf and R1tf) and without (R2nf and R2tf) viscoelastic relaxation reach stress change magnitudes of as much as two orders of magnitude higher compared to models that include only viscoelastic relaxation (R3nf and R3tf) (Figs. 3 and 4). The highest stress changes occur on the receiver faults closest to the source fault (as much as 4 MPa, R1tf), but the outermost receiver faults also show higher values of as much as 0.4 MPa (R1nf), which could bring the faults closer to failure (King et al., 1994). Both types of models (R1nf and R1tf, R2nf and R2tf) show the same stress field evolution in the first two years after the earthquake, which indicates that poroelastic effects have larger impact on the stress field in the early postseismic phase than viscoelastic effects in these models (Figs. 3 and 4). A strong impact of poroelastic effects during the early postseismic phase is observed in almost all models (Figs. 310). Note that the signal from poroelastic effects is recognizable as much as ~50–75 km away from the source fault and that the stress transfer to the receiver faults is not hindered by the impermeability of the model faults, as indicated by the spatially continuous surface velocity fields across the receiver faults (e.g., Figs. 3D and 4D). The reason is that the coseismically induced fluid flow occurs close to the source fault (cf. Peikert et al., 2022) but the resulting stresses can be transmitted farther away. The signal from poroelastic effects may overlap with the signals from viscoelastic relaxation and interseismic stress accumulation for months to several years. How long the poroelastic effects persist and cause Coulomb stress changes that may bring faults closer to failure depends on the permeability of the crust and the corresponding pore pressure dissipation time. While poroelastic effects decay within two months in models with high permeabilities of 10−10 m2 (P1nf and P1tf, PV1nf and PV1tf, PV3nf and PV3tf) due to the fast fluid flow, their signal can be observed for two and more than ten years in models with a permeability of 10−12 m2 (R1nf and R1tf, R2nf and R2tf, V2nf and V2tf) and 10−14 m2 (P2nf and P2tf), respectively. During these time periods, the Coulomb stress changes on receiver fault RF5 are still high enough (>0.1 MPa) to bring the faults closer to failure. The magnitudes of stress changes due to poroelastic effects decrease with decreasing permeabilities by one order of magnitude between permeabilities of 10−12 m2 and 10−16 m2. In models with a permeability of 10−16 m2 (P3nf and P3tf, PV2nf and PV2tf, PV4nf and PV4tf), poroelastic effects still overlap with viscoelastic relaxation and interseismic stress accumulation for as long as 50 years but are only high enough to bring RF5 or RF7 closer to failure to about the 10th year. Therefore, the poroelastic effects observed in our models act on time scales that overlap with that of the spatio-temporal evolution of the postseismic viscoelastic relaxation process. The Coulomb stress distribution of the models including only viscoelastic relaxation but no fluid flow (R3nf and R3tf) indicates that the influence of viscoelastic relaxation on the stress field is already recognizable in the first month, but the Coulomb stress changes due to viscoelastic relaxation are lower and hence overlapped by the stress changes caused by poroelastic effects. As soon as the poroelastic effects decay, the viscoelastic relaxation signal starts to dominate the stress field (Figs. 3 and 4). High viscosities in the lower crust of 1022 Pa·s (V2nf and V2tf, PV3nf and PV3tf, PV4nf and PV4tf) lead to increasing Coulomb stress changes, which remain constant over decades on all receiver faults with average values of 0.02 MPa and thus slightly higher than the values of the interseismic stress accumulation (0.01–0.02 MPa in our models). Models with a viscosity of 1020 Pa·s additionally show negative stress changes on RF5 for decades. Both viscosity values lead to Coulomb stress changes, which still outweigh the continuous interseismic stress increase in the 50th year after the earthquake. In models with a low viscosity of 1018 Pa·s (V1nf and V1tf, PV1nf and PV1tf, PV2nf and PV2tf), the viscoelastic relaxation causes positive and negative Coulomb stress changes on all faults with higher values and a significant change of distribution and magnitude over time. Until the fifth year (normal fault) on RF2, RF5, RF7, and RF10 and until 10th year (thrust fault) on RF5 and RF7, the Coulomb stress changes are high enough to bring the receiver faults nearest to the source fault closer to failure if the permeability is sufficiently high at the same time (V1nf and V1tf, PV1nf and PV1tf). In these models, the viscosity is low enough that the signal from viscoelastic relaxation already starts to dominate the stress field from the third month onwards and overlaps the signal from poroelastic effects. By the 50th year, the interseismic stress increase dominates the stress field.

If the low viscosity is combined with a low permeability (PV2nf and PV2tf), viscoelastic relaxation dominates the stress field from the first month onward and the Coulomb stress changes are still high enough in the 10th year to bring RF7 (normal fault) and RF5 (thrust fault) closer to failure. This combination of permeability and viscosity in the crust highlights the possibility that viscoelastic relaxation may dominate already in the early postseismic phase. The signal caused by viscoelastic flow may be further enhanced if the lower crust has a power-law viscosity because such a rheology would lead to high strain rates in the early postseismic phase. Later, the strain rate would decrease more rapidly than for a linear Maxwell rheology. As shown by Freed and Bürgmann (2004), the power-law flow over a period of three years (2000–2003) after the Landers and Hector Mine earthquakes (southern California, USA) indicated that the change in the effective viscosity of the mantle was equivalent to an order of magnitude increase in linear viscosity.

In models that combine a low permeability with a high viscosity, the signals from both effects are weak but long lasting, with the result that the Coulomb stress changes remain to be a superposition of stress changes caused by poroelastic effects, viscoelastic effects, and interseismic stress accumulation in the 50th year. Based on our model results, poroelastic effects may affect the stress field for decades if the permeability is sufficiently low and viscoelastic relaxation is already recognizable in the first month after the earthquake. Hence, poroelastic effects and viscoelastic relaxation overlap over longer time scales than expected (e.g., Freed and Lin, 2001; Luo and Liu, 2010; Albano et al., 2017, 2019, 2021; Nespoli et al., 2018), which should be considered for the calculation of postseismic Coulomb stress changes.

In addition, poroelastic and viscoelastic effects may interfere with afterslip. Afterslip may occur during the early postseismic phase and typically affects the shallow and/or deep parts of the fault plane (e.g., Freed, 2005; Wang and Fialko, 2018; Liu-Zeng et al., 2020; Jin et al., 2022). Its signal is therefore rather limited in space. In some cases, afterslip has been inferred to be the only process that could be detected in the postseismic deformation field (Freed, 2007).

Coulomb Stress Changes and Surface Deformation

As shown by our model results (Figs. 38), both poroelastic effects and viscoelastic relaxation lead to considerable crustal movements around the source fault, which leave a signal in the vertical and horizontal surface displacement fields and ultimately cause the observed Coulomb stress changes. The magnitudes of the stress changes strongly correlate with the surface velocities. If a model shows high surface velocities, it also experiences high Coulomb stress changes in the same area. A velocity decrease by one order of magnitude generally leads to a decrease of the stress change magnitude by one order of magnitude. The movements in horizontal directions control the distribution of the Coulomb stress changes on the faults. Extension leads to positive stress changes on normal faults and negative stress changes on thrust faults, whereas areas of shortening experience negative stress changes on normal faults and positive stress changes on thrust faults (Bagge and Hampel, 2017). The velocities and spatio-temporal evolution of the surface deformations strongly depend on the permeability and viscosities of the crust (Peikert et al., 2022).

Models with high permeabilities (R1nf and R1tf, R2nf and R2tf, P1nf and P1tf, V2nf and V2tf) show high surface velocities in the early postseismic phase that are two orders of magnitude higher than surface velocities caused by only viscoelastic relaxation (R3nf and R3tf). The highest movements can be found on receiver faults RF5 and RF7, which experience the highest stress changes. In the early postseismic phase in the normal fault models, shortening between RF5 and the source fault as well as between RF7 and the source fault leads to negative stress changes on the upper part of RF5 and the lower part of RF7. In contrast, the other parts of these faults show positive stress changes due to the extension within the footwall and hanging wall. In most models, high velocity perturbations occur between the source fault and the receiver faults along strike to the source fault (RF2 and RF10) in the early postseismic phase, caused by poroelastic effects. This leads to higher Coulomb stress changes at the parts of RF2 and RF10 near the source fault. Due to the fast pore pressure diffusion, the velocities strongly change and the velocity field is dominated by the signal from viscoelastic relaxation and the regional extension and shortening from the third month (P1nf and P1tf) and the fifth year (R1nf and R1tf, R2nf and R2tf, V2nf and V2tf) onwards. The models considering only viscoelastic relaxation but no pore fluid pressure (R3nf and R3tf) shows velocity perturbations with low velocities from the first month onward, which negligible changes over the next 50 years. This underlines that a signal from viscoelastic relaxation is already recognizable in the first month.

A lower permeability leads to prolonged pore pressure diffusion and hence to a slower change in the velocity and stress field magnitudes and pattern. The velocity field is affected by poroelastic effects until at least the 10th year for permeabilities of 10−14 and 10−16 m2, respectively. As a result, the stress field shows a similar evolution. The areas without significant velocity perturbations generally do not experience high Coulomb stress changes. For example, in the normal fault model with a permeability of 10−14 m2, the velocity perturbations on the hanging wall of the source fault are limited to the area between the source fault and receiver fault RF5, resulting in low Coulomb stress changes on the outer receiver faults of the hanging wall (RF1 and RF9). Caused by the strong viscoelastic relaxation, the surface velocity field is highly disturbed in models with a low viscosity (V1nf and V1tf), showing higher velocities until the 10th year and a pattern of several areas of extension and shortening, resulting in a mixed pattern of positive and negative stress changes on all faults parallel to the source fault. The surface velocity field is dominated by the signal from viscoelastic relaxation already in the third month and strongly changes within 20 years, similar to the stress field. In normal fault models with a viscosity of 1020 Pa·s in the upper crust, an area of shortening occurs within the hanging wall around RF5 after the poroelastic effects decay, resulting in negative Coulomb stress changes on the upper part of RF5 until the 50th year. Both zones of shortening and extension only change negligibly over time and overlap the regional extension field over 50 years. The rest of the model domain indicates enhanced extension similar to the regional extension field and hence solely positive stress changes with magnitudes, which are typical for the interseismic stress accumulation controlled by the regional deformation (Bagge and Hampel, 2017). A high viscosity of 1022 Pa·s (V2nf and V2tf) leads to weak but long-lasting viscoelastic relaxation, which results in weak velocity perturbation and slow surface velocities over decades. The surface velocity field changes only negligibly over time. The Coulomb stress change magnitudes and distribution also do not change significantly over the decades. The velocity field pattern only slightly differs from the velocity field pattern caused by the regional extension and shortening and hence only experiences positive stress changes.

In most models, the surface velocity and stress field are dominated by poroelastic effects in the early postseismic phase (except V1nf and V1tf, and PV2nf and PV2tf). The highest movements and stress changes can be found in the first month after the earthquake, followed by a strong decrease in the velocities and stress change magnitudes in the following months if the permeability is sufficiently high. As a result, the postseismic velocity and stress field integrated over one year show elevated values and resemble the velocity and stress field of the first month if the permeability is higher than 10−14 m2 and hence include a strong signal from poroelastic effects, which overlaps and is much stronger than the signal from incipient viscoelastic relaxation. Therefore, it is important to choose the right time interval for the calculation of Coulomb stress changes and for the analysis of postseismic velocity fields derived from geodetic data. In some cases, the poroelastic effects have not disappeared in the early postseismic phase (Peikert et al., 2022) and the postseismic velocity field cannot be interpreted to reflect solely the signal from incipient viscoelastic relaxation (e.g., Aoudia et al., 2003; Liu-Zeng et al., 2020; Mandler et al., 2021).

Comparison with Coulomb Stress Changes for Natural Intra-Continental Earthquakes

In the evaluation of the seismic hazard of a region, Coulomb stress changes are routinely calculated after major earthquakes, mostly focusing on static Coulomb stress changes. Aftershocks, however, may occur a few days to months to years after the main earthquake, which requires consideration of postseismic transient processes, as demonstrated by our study. Most previous studies, however, focused on a specific time interval and therefore restricted their analysis to a combination of static stress changes with either pore fluid pressure changes (Antonioli et al., 2005; Albano et al., 2017, 2019; Tung et al., 2018) or postseismic viscoelastic relaxation (e.g., Luo and Liu, 2010; Wang et al., 2014). Only a limited number of studies have considered static and both types of transient Coulomb stress changes in their analyses (Freed and Lin, 2001; Masterlark and Wang, 2002; Barbot and Fialko, 2010; Zhu and Miao, 2015; Tung and Masterlark, 2018). In the following, we qualitatively compare our findings with Coulomb stress calculations from natural dip-slip earthquakes with respect to the relative importance and spatio-temporal evolution of poroelastic effects and viscoelastic relaxation for aftershock triggering.

For example, Zhu and Miao (2015) analyzed the six-month aftershock sequence of the 2013 Lushan earthquake (eastern margin of the Tibetan Plateau), concentrated in a 40 × 15–km–wide, NE-SW-oriented area around the mainshock at depths between 10 km and 22 km. They showed that most of the aftershocks occurred in areas with negative static Coulomb stress changes and hence were not directly triggered by coseismic static Coulomb failure stress changes due to the Lushan mainshock. The authors assumed that aftershocks may have been triggered by postseismic stress transfer produced by changes in pore pressure and fluid flow but excluded viscoelastic relaxation as a mechanism. Our models indicate that poroelastic effects strongly alter the stress field shortly after the earthquake and positive stress changes are found in areas on the receiver faults near the source fault, which experience negative static stress changes before (cf. Piombo et al., 2005). The area investigated by Zhu and Miao (2015) is comparable to the position of receiver faults RF2, RF5, RF7, and RF10 in our models. In most of our models, these faults experience a considerable stress increase caused by poroelastic effects in the first six months.

Other examples of dip-slip earthquakes, for which aftershock triggering by static versus transient stress changes was evaluated, include the 1997 Umbria-Marche and the 2016 Amatrice-Visso-Norcia earthquake sequences in the Central Apennines, Italy. As shown by Cocco et al. (2000), the three largest earthquakes of the 1997 Umbria-Marche sequence may have been triggered solely by static Coulomb stress changes. However, the agreement between modeled and observed spatial patterns of the entire earthquake series was further increased by considering enhanced fluid flow (Cocco et al., 2000). Later, Antonioli et al. (2005) confirmed that postseismic pore pressure relaxation played a crucial role in the evolution of aftershocks during the 1997 Umbria-Marche earthquake sequence, given that most of the main shocks and aftershocks occurred in areas where the pore pressure increased around the fault. The authors used a permeability of 7.4 × 10−12 m2 for this area. As shown by our models, permeabilities of this order of magnitude may lead to strong poroelastic effects, which dominate the stress field from the first month to until at least the second year after the earthquake (Fig. 3). Our finding agrees with the model results from Albano et al. (2017, 2019), who showed that considerable poroelastic effects with a strong impact on the stress field may occur in the early postseismic phase because earthquake-induced pore pressure gradients fully dissipate after a few days to a few months for sufficiently high permeability (Albano et al., 2017, 2019). In the Amatrice-Visso-Norcia earthquake sequence, the 26 October 2016 event occurred even before the fluid overpressure induced by the 24 August Amatrice event had fully dissipated (Albano et al., 2019).

In contrast, a low permeability of 10−16 m2 for the crust beneath the Central Apennines was derived by Tung and Masterlark (2018), who calculated Coulomb stress changes for the 2016 Amatrice-Visso-Norcia earthquakes and found that poroelastic effects dominate the postseismic stress field and are responsible for the aftershock triggering. For such a permeability value, our model results indicate a considerably slower dissipation of the pore pressure than for higher permeabilities, with the consequence. However, Tung and Masterlark (2018) only considered viscoelastic behavior in the mantle and assumed that the contribution from viscoelastic relaxation is negligible. Given that Riva et al. (2007) and Aoudia et al. (2003) derived a viscosity of 1018 Pa·s for the lower crust beneath the Central Apennines from the postseismic deformation after the 1997 Umbria-Marche earthquake sequence, we argue that viscoelastic relaxation has also contributed to aftershock triggering. This is illustrated by our model PV2nf, which has a low permeability of 10−16 m2 combined with a low viscosity of 1018 Pa·s. For such a combination, poroelastic effects influence the stress field for several years due to slow pore pressure dissipation but are overprinted by the signal from viscoelastic relaxation already in the first month to decades later. This leads to higher Coulomb stress changes over larger distances, especially on the receiver faults along strike and parallel to the source fault.

In summary, our findings imply that the analysis of static Coulomb stress changes may not be a reliable tool for predicting stress transfer after major earthquakes because poroelastic effects and viscoelastic relaxation may alter both the magnitude and spatial distribution of the coseismically induced Coulomb stress changes already during the first month after the earthquake. Both poroelastic effects and viscoelastic relaxation should be considered when calculating Coulomb stress changes and analyzing geodetic measurements of earthquake-induced surface deformation. Poroelastic and viscoelastic effects may be isolated from each other by analysis of the postseismic surface deformation with respect to the magnitudes and patterns of horizontal and vertical surface velocities. In particular, surface velocities tend to be significantly higher in the presence of poroelastic effects compared to viscoelastic relaxation alone (Figs. 3 and 4). Also, the direction of the vertical movement of the hanging wall and footwall of the source fault differs depending on whether poroelastic effects are present or not. To tease poroelastic and viscoelastic effects apart for natural earthquakes, high-resolution geodetic data on vertical movements in the entire area around the earthquake source fault would therefore be highly desirable.

Transient processes should also be considered in the analysis of paleo-earthquake sequences (e.g., Verdecchia et al., 2018; Bagge et al., 2019) when the subsequent events occur on time scales of years to decades. As shown by our model results, such models should also account for interseismic stress accumulation because this process dominates the Coulomb stress change patterns after the signal from the transient process has disappeared.

We used 3-D finite-element models of intra-continental normal and thrust faults including coseismic slip, poroelastic effects, postseismic viscoelastic relaxation, and interseismic stress accumulation to investigate the relative importance of these processes for the spatio-temporal evolution of postseismic Coulomb stress changes. The models show that coseismic stress changes do not persist through the early postseismic phase but are considerably altered by poroelastic effects and viscoelastic relaxation within the first month after the earthquake. Poroelastic effects cause high Coulomb stress changes and dominate in the early postseismic phase but may still influence the stress and surface velocity field several years after the earthquake for sufficiently low permeability. Postseismic viscoelastic relaxation can affect the surface deformation patterns and the stress field already in the first few months after the earthquake if the viscosity of the lower crust is low enough. Depending on the combination of upper-crustal permeability and lower-crustal viscosity, our results indicate that the signals from poroelastic effects and postseismic viscoelastic relaxation may overlap in the early postseismic phase for up to several years. Poroelastic effects and viscoelastic relaxation have a strong influence on the magnitudes and distribution of postseismic Coulomb stress changes and should be considered together with interseismic stress accumulation when analyzing Coulomb stress transfer between faults and analyzing geodetic data on postseismic surface deformation.

Science Editor: David E. Fastovsky
Associate Editor: Jason W. Ricketts

We thank Editor D. Fastovsky, Associate Editor J. Ricketts, and two anonymous reviewers for the thoughtful and constructive comments, which helped to improve the manuscript. Funding by the German Research Foundation (DFG grant HA 3473/11-1) is gratefully acknowledged.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.