## 1 Introduction: The concept of an earthquake deformation cycle

The cyclic buildup and release of elastic strain around earthquake faults is accompanied by transient deformation at the Earth's surface, especially in the years to decades following large earthquakes. Measured by modern satellite geodetic techniques, these signals are both winnowing out competing mechanisms and refining our knowledge of the underlying processes. To characterize post-seismic transients in both space and time, we will interpret them within the context of the earthquake deformation cycle.

### 1.1 Temporal characteristics

Fig. 1 shows the general features of surface displacement as a function of time as observed near the fault. The sudden abrupt drops in displacement correspond to the co-seismic stress release that occurs in the few seconds to several minutes as the fault slips. After the rupture ceases, post-seismic transient movements can continue for years to decades. Steady inter-seismic deformation occurs at a fairly uniform rate throughout the cycle, even while transient motions persist. This inter-seismic strain accumulation reloads the fault, which will slip again when the stress on it exceeds the frictional strength of the fault, and the cycle recurs, typically on a time scale of several decades to thousands of years.

This concept has guided attempts to predict the date of an earthquake in the future. The time-predictable model states that the time to the next earthquake is the co-seismic slip in the most recent earthquake divided by the fault-slipping rate [107]. An equivalent expression also applies in stress since it is linearly proportional to displacement in an elastic medium such as the Earth's crust. A rigorous statistical application of the linear time predictable model, however, failed to predict the most recent earthquake at Parkfield, California [74]. Nonetheless, the concept of an earthquake deformation cycle provides a contextual framework for interpreting geodetic measurements [111–113].

According to this concept, the vector position **x** of a benchmark on the ground varies as a function of time t:

$$\mathbf{x}(t)=x({t}_{0})+H({t}_{\mathrm{q}})\mathbf{u}+(t-{t}_{0})\mathbf{v}+\mathbf{w}(t)$$ | (1) |

**u**is its co-seismic displacement, H is the Heaviside ‘step’ function that switches from 0 to 1 at time ${t}_{\mathrm{q}}$,

**v**is the inter-seismic velocity, and $\mathbf{w}(t)$ is a function of time describing the temporal evolution of the post-seismic transient. In the simplest case, if the post-seismic transient is negligible ($\mathbf{w}=0$) and the co-seismic displacement

**u**is the same in each ‘characteristic’ earthquake, then the cycle is periodic (and thus predictable) with a period (or recurrence time) $T=u/v$. The Earth, however, does not always follow such idealized behavior.

To facilitate the discussion, we will make two simplifying assumptions. First, we will assume that the temporal evolution and spatial distribution of the post-seismic transient displacement at time t and vector position coordinate ${\mathbf{x}}_{i}$ is separable (e.g., [36])

$$\mathbf{W}(t,\mathbf{x})=f(t)\times \mathbf{G}(\mathbf{x})$$ | (2) |

In other words, we can completely describe the spatio-temporal dependence of the post-seismic signal with one map and one time series. This assumption does not apply to models using similarity variables that explicitly couple **x** and t, such as the stress diffusivity parameter in the Elasser viscous model (e.g., [37]).

Second, we assume that the function $f(t)$ is zero just after the mainshock at time ${t}_{\mathrm{q}}^{+}$. It then rises to $f(\tau )=1-1/e=0.63$ after a relaxation time τ and finally to 1 after an infinitely long time interval. Accordingly we can write the total, or fully relaxed, post-seismic displacement field as:

$$\underset{t\to \infty}{\mathrm{lim}}\mathbf{W}(t,\mathbf{x})=\underset{t\to \infty}{\mathrm{lim}}[f(t)]\mathbf{G}(\mathbf{x})=\mathbf{G}(\mathbf{x})$$ | (3) |

This displacement field records the motion of a point on the ground along its trajectory from the moment ${t}_{\mathrm{q}}^{+}$ just after the earthquake to its final, fully relaxed position.

### 1.2 Spatial characteristics

When a fault slips seismically, a continuous marker curve, such as a road, becomes discontinuous. Two points that were adjacent at positions ${\mathbf{x}}_{1}$ and ${\mathbf{x}}_{2}$ before the earthquake (at time ${t}_{1}$) are permanently separated in space after the earthquake (at time ${t}_{2}$). The difference in their displacements ${\mathbf{u}}_{1}$ and ${\mathbf{u}}_{2}$ defines the slip vector:

$$\mathrm{\Delta}\mathbf{u}={\mathbf{u}}_{2}-{\mathbf{u}}_{1}=[{\mathbf{x}}_{2}({t}_{2})-{\mathbf{x}}_{2}({t}_{1})]-[{\mathbf{x}}_{1}({t}_{2})-{\mathbf{x}}_{1}({t}_{1})]$$ | (4) |

At any single point ${\mathbf{x}}_{i}$ on the ground, a geodetic benchmark, for example, the cumulative displacement over times much longer than the duration of a single earthquake cycle is

$${\mathbf{u}}_{i}^{(\mathrm{plate})}={\mathbf{u}}_{i}^{(\mathrm{inter})}+{\mathbf{u}}_{i}^{(\mathrm{co})}+{\mathbf{u}}_{i}^{(\mathrm{post})}$$ | (5) |

Over long time intervals, it is convenient to think in terms of rate of displacement, i.e. velocity $\mathbf{v}=\partial \mathbf{u}/\partial t$. Accordingly, ${\mathbf{v}}^{(\mathrm{plate})}$ represents the long-term velocity averaged over the last several million years as in the conventional NUVEL-1A plate motion model [26,27]. Note, however, that the above sum applies in displacement **u**, but not velocity **v**, because the instantaneous co-seismic displacement ${\mathbf{u}}^{(\mathrm{co})}$ is not differentiable in time. For benchmarks located on the stable, interior parts of continents, where earthquakes are extremely infrequent and fairly small, we may neglect co- and post-seismic displacements. Under these conditions, we expect the long-term velocity ${\mathbf{v}}^{(\mathrm{plate})}$ to approximately equal the inter-seismic velocity ${\mathbf{v}}^{(\mathrm{inter})}$ that we can measure by repeated geodetic surveys. The validity of this approximation has been well confirmed by two global combinations of geodetic velocity estimates including hundreds of benchmarks: ITRF-2000 [3] and REVEL [105]. Except for a few plates (e.g., Nazca and India), differences between the long-term and short-term estimates are within measurement uncertainties. Improved estimates of the terrestrial reference frame are forthcoming and will provide more refined comparisons [28,97].

### 1.3 Capturing post-seismic transients with satellite geodesy

After having written an expression for the post-seismic deformation field, we now discuss how to measure it using satellite geodesy. Navigation systems such as the Global Positioning System (GPS) estimate the coordinates of a benchmark in some well-defined reference frame. For our purposes, we can consider these three coordinates as an absolute vector position $\mathbf{x}(t)$ valid at time t. With multiple surveys, we can determine a time series composed of the three components of the position vector at distinct times. To describe one component of a position time series including an earthquake, we must specify at least three scalar parameters: the initial position $x({t}_{0})$, the inter-seismic velocity v and the co-seismic displacement u. In addition, we would like to estimate at least one more parameter to describe the temporal evolution $f(t)$ part of the post-seismic transient w. For clarity, we will consider several real-world cases.

The first case involves a benchmark that was installed and measured only twice: once before the earthquake, and once again during the excitement following the earthquake. In this case, we have only two data points. Consequently, we can only estimate two parameters, typically $x({t}_{0})$ and u. We must assume or neglect the inter-seismic velocity v and the post-seismic transient w. Such a benchmark does not constrain models for post-seismic deformation.

In the second case, we have four data points: two before, and two after, the earthquake. Now we can estimate four parameters: $x({t}_{0})$, u, v, and w. To describe the post-seismic transient w with a single parameter requires assuming a simplest form for $f(t)$, such as a piece-wise linear ‘triangle’ function of time:

$${f}_{\mathrm{lin}}(t)=\frac{(t-{t}_{\mathrm{q}})}{{\tau}_{\mathrm{lin}}}\phantom{\rule{1em}{0ex}}\text{when}{t}_{\mathrm{q}}t\u2a7d{\tau}_{\mathrm{lin}}$$ | (6) |

In the third case, we assume a benchmark that we have been lucky enough to install and survey at least twice before an earthquake occurs. After the earthquake, an increase in funding allows at least three surveys in the post-seismic interval. In this case, we can estimate at least two more parameters to describe a curved shape for the temporal evolution $f(t)$ part of the post-seismic transient **w**. The shape of the curve depends on the material properties of the Earth's crust, as we will now describe.

## 2 Rheological analogies

To simplify the discussion, we draw analogies between the continental lithosphere and three mechanical analogs. First, we consider elastic behavior like that of a rubber eraser. Then we consider the effect of pore fluids in analogy with a sponge impregnated with water. Finally, we consider a flowing substrate analogous to honey.

### 2.1 Rubber eraser analog

In an elastic solid, stress is linearly proportional to strain (Hooke's law). Two coefficients suffice to describe the rheology. Although several different parameterizations appear in the literature, including the bulk modulus k, Young's modulus E, Lamé first coefficient λ, the rigidity μ, and Poisson's ratio ν, any two completely specify the elastic medium. The relationships between these quantities are simple algebraic derivations that appear in many exams and textbooks (e.g., [58]). For simplicity, most studies assume that $\lambda =\mu $, so that these parameters drop out of the expressions for surface displacement. Such a medium, called a Poisson solid, has a Poisson's ratio ν of 1/4. In this case, the mathematical expressions for displacements, strains and stresses due to fault slip [80] no longer depend on λ, μ, or ν. It is a reasonable approximation because estimates in the upper crust vary from $\nu =0.23$ to $\nu =0.28$, as inferred from P- and S-wave velocities [30,86]. Poisson's ratio also depends on the amount of fluid included in the rock [117].

The Lamé coefficient μ is called the rigidity in mechanics and the shear modulus in seismology. Typical values (assumed) for rigidity μ in the Earth's crust range from 30 to 36 GPa, but values as low as 10 GPa [25] and as high as 50 GPa [10] have been used. For a fluid, the rigidity is zero because the medium cannot transmit shear stress.

The simplest recipe for describing deformation due to fault slip approximates the Earth as an elastic half space. This formulation has become the conventional standard for modeling co-seismic deformation. Okada [80] derives the expressions for the co-seismic (permanent) displacement at the Earth's surface caused by a slip across a rectangular fault surface at depth in closed analytic form. Accordingly, the surface displacement field $\mathbf{u}(x,y)$ due to a dislocation Δ**U** across a rectangular fault surface with position coordinates (${x}_{\mathrm{f}},{y}_{\mathrm{f}}$), width W, length L, strike α, dip δ, buried at depth d in an isotropic elastic medium with Poisson's ratio ν and shear modulus (rigidity) μ below an observing station located at ($x,y$) is:

$$\mathbf{u}(x,y)=\mathit{Okada}(\mathrm{\Delta}\mathbf{U},\nu ,\mu ,L,W,{x}_{\mathrm{f}},{y}_{\mathrm{f}},x,y,d,\delta ,\alpha )$$ | (7) |

For the complete set of equations, see Okada's 1985 paper that also corrects previous derivations. In the usual case of a Poisson solid with $\nu =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$, the expression no longer depends on ν or μ. For points at depth in the half-space, Okada [81] generalizes the formulation to an arbitrary point located anywhere within in the half-space. Okada's formulation has been applied successfully to over a hundred earthquakes, as reviewed by Feigl [32].

Why is the elastic dislocation formulation so important?

First, it captures the essential physics of deformation due to fault slip in a simple form. During the co-seismic interval, that lasts several hundred seconds at most, the lithosphere behaves elastically. Accordingly, the expressions for displacement **u** at the surface are linear functions of the fault slip Δ**U** at depth. This linearity allows us to compute the displacement field for a geometrically complicated fault system as the sum of the displacement fields produced by many smaller, simpler fault elements, or ‘patches’. It also simplifies the inverse problem of estimating the distribution of slip Δ**U** at depth from geodetic measurements of **u** at the surface. The half-space approximation appears to be valid for all but the largest earthquakes because the co-seismic displacement exceeds the millimeter level only within a few hundred kilometers of the fault. In other words, the length of the displacement vector $|\mathbf{u}|$ falls off steeply with the distance from the station location to the dislocation source. Several recent studies compare the half-space approximation to a radially symmetric spherical Earth model [9,25,75,110].

Second, the parameters in the elastic dislocation model are quantities that we can estimate independently. For example, we can infer the slip ΔU, the length L, and width W of the fault rupture from the seismic moment ${\mathrm{M}}_{0}=\mu ULW$ estimated from seismograms [1,43]. This approach involves fewer assumptions than, say, imposing unknown or poorly known forces at far-field boundaries.

Third, the elastic formulation is straightforward to write as computer code. A public-domain computer program written in FORTRAN performs these calculations for an arbitrary observation point on the Earth's surface [33]. Other formulations produce equivalent results [2,40,88,114].

Finally, the elastic dislocation formulation can be extended from the co-seismic case, where it provides an accurate physical description of how faults slip (a ‘model’ in the sense of a prototype or proxy), to other cases, where it can provide an approximate description (a ‘model’ in the sense of a heuristic conceptual analogy). It has become the universally accepted method of choice for calculating co-seismic displacements and slip distributions and has also been widely used to model post- and inter-seismic slip, as we discuss below.

### 2.2 Eraser analogy for post-seismic deformation

In this mechanism for post-seismic deformation, slip continues in the same plane as the mainshock rupture, as sketched in Fig. 3. An early, well-recorded example appeared in California following the Imperial Valley earthquake of 15 October 1979. In the days following the ${M}_{w}=6.5$ mainshock, the fault continued to slip at rates of 5 to 20 mm d^{−1}, diminishing to about 1 mm d^{−1} ten weeks later [44]. The largest total post-seismic displacements and afterslip rates occurred near the midpoint of the co-seismic surface rupture [44]. Near the northern termination of the co-seismic rupture, the after slip rate averaged $1.08\pm 0.8\text{}\text{mm}\phantom{\rule{0.2em}{0ex}}{\text{d}}^{\mathrm{-1}}$ between 27 October and 13 December 1979 [59]. Two leveling stations, one 5 m west of the fault trace, the other 30 m east, changed their relative vertical elevation by 14 cm (about as much as the co-seismic slip) in the first hundred days after the mainshock [106]. Six weeks after the event, the total post-seismic slip was 8 and 11 cm at two places, representing 13 to 29% of the co-seismic slip, respectively [23]. The overall picture is that afterslip represents continued slip in the same direction as the co-seismic rupture. It seems to be a rapid process with an exponential time scale of 17 days [23], but a logarithmic temporal function also fits the data [59,106]. In many cases, the afterslip is confined to very shallow depths ($<1\text{\u2013}2\text{}\text{km}$), as the semi-consolidated material near the surface ‘catches up’ to the co-seismic slip at greater depths. The geodetic signature seems to resemble that of creep [24,119].

The 1999 earthquake in Izmit, Turkey provides a complementary example of the post-seismic afterslip, in this case occurring at depths greater than that of the co-seismic slip [15,17,18,31,46,92,120]. Maximum afterslip rates decayed from greater than 5 mm d^{−1} just after the earthquake to 3 mm d^{−1} some three months later, with an exponential relaxation time of 57 days [17,31]. In the case of Izmit, the afterslip seems to be considerably deeper than the co-seismic rupture zone. In neither case can the afterslip inferred from the geodetic observations be explained by the recorded aftershocks.

A simple, empirical description of the temporal evolution of post-seismic position is a logarithmic function

$${f}_{\mathrm{log}}(t)={c}_{\mathrm{log}}\mathrm{log}[1+\frac{t-{t}_{\mathrm{q}}}{{\tau}_{\mathrm{log}}}]\phantom{\rule{1em}{0ex}}\text{for}t{t}_{\mathrm{q}}$$ | (8) |

$${f}_{\mathrm{log}}(t)={c}_{\mathrm{log}}\mathrm{log}({\tau}_{\mathrm{log}}+t-{t}_{\mathrm{q}})-{c}_{\mathrm{log}}\mathrm{log}({\tau}_{\mathrm{log}})$$ | (9) |

In other words, the final term is already included as the co-seismic offset **u** in the complete expression (1) for the position coordinate vector $\mathbf{x}(t)$. In practice, we can estimate only one such additive parameter for a given co-seismic epoch in a single coordinate time series.

At longer times, the logarithmic function rises slowly without bound. This is a serious drawback because post-seismic deformation cannot go on forever. One way to handle this shortcoming is to subtract a linear trend from the long-term behavior. The difference between the logarithmic form ${f}_{\mathrm{log}}$ and the linear form ${f}_{\mathrm{lin}}$ can resemble the exponential form ${f}_{\mathrm{exp}}$ discussed below. Although the logarithmic functional form seems to fit several post-seismic time series (e.g., station ROCH following Landers [36]), it does not correspond to a simple rheology. It can, however, be related to afterslip in a fault zone with rate- and state-dependent friction [65].

Extending this idea, Perfettini and his colleagues derive a model of brittle creep in the middle layer of the crust above a subduction zone using a rate- and state-friction law [84,85]. Similarly, Kato and others simulate post-seismic deformation on a vertical, strike-slip fault with rate- and state-dependent friction overlying a Maxwell viscoelastic half space [54]. Other studies formulate models along similar lines [60–62]. In its simplest form, a rate- and state-dependent friction law leads to a temporal evolution that shares attributes of the logarithmic and exponential forms. During the early post-seismic interval, it rises quickly, like the exponential form. At later times, however, it continues to rise, like the logarithmic form [36,84].

In space, afterslip produces a deformation field with characteristic distance scale on the order of magnitude of the depth at which the post-seismic slip occurred, less than 5 km for Imperial Valley, but 20 to 30 km for Landers and Izmit. Overall, the deformation field produced by post-seismic afterslip resembles the co-seismic field in shape and direction, but with smaller amplitude, with displacement vectors typically $\sim 1/10$ as long as their co-seismic values (Fig. 4).

### 2.3 Sponge analog

Earthquakes change the state of stress in the material around the generative fault, and if the medium is fluid saturated the seismically-induced static stress gradients generate fluid flow that couples into the deformation field. This stress perturbation can be as large as 1–10 MPa (10 to 100 bar) on the fault and decays with increasing distance from each dislocation source. It can affect hydrological conditions, as observed in streams and water wells [71–73]. Hydrothermal activity seems to have changed prior to the 1989 Loma Prieta earthquake in California [108]. The level of water in wells changed during or following the 1999 earthquakes in Turkey [121]. Similar effects were observed in an abandoned mine over several months after a ${M}_{w}\text{}7.2$ earthquake in Japan [79]. Fountaining of water and sediments was observed during the 2001 earthquake in India [87].

Similarly, several water wells in Iceland began to flow freely (becoming artesian) in the days following the ${M}_{w}\text{}6.6$ earthquake of 17 June 2000 [12]. Indeed, by mapping the polarity of the change (increase or decrease) in water level over various wells, Icelandic geophysicists were able to recognize the four-quadrant pattern characteristic of a strike-slip earthquake [12,52]. Perturbing the pressure of fluids trapped in the pores, fissures and interstices of rocks alters both their effective elastic moduli and the mechanical behavior of the fault zone itself.

### 2.4 Poisson recipe for modeling poro-elastic deformation analogous to a sponge

Increasing the fractional fluid volume in a rock changes its rheological behavior. As a result, the ratio ${V}_{\mathrm{p}}/{V}_{\mathrm{s}}$ of seismic wave velocities also changes, in a phenomenon known as dilatancy [39,42,77,101,103,117]. Consequently, the Poisson's ratio and changes with permeability. Assuming that the available pore space in crustal rocks is completely filled with water, Rice and coworkers [21,22,93] derive expressions for Poisson's ratio ν under ‘undrained’ conditions. Such a perturbation in the Poisson's ratio can explain post-seismic deformation over time scales of months and distance scales on the order of 10 km. According to this model, the co-seismic perturbation in the stress field alters the hydrological conditions by changing the fluid pressure instantaneously at time ${t}_{0}$, as sketched in Fig. 5. Then, the hydrological conditions readjust as the fluids drain from the interstices and fissures. Eventually, the rock returns to its initial, ‘drained’ condition. The simplest description of this process allows the Poisson's ratio to vary as a function of time. Just after the co-seismic rupture, at time ${t}_{\mathrm{q}}^{+}$, ‘undrained’ conditions apply:

$$\nu ({t}_{\mathrm{q}}^{+})={\nu}_{\mathrm{undrained}}=0.31$$ | (10) |

$$\nu ({t}_{1})={\nu}_{\mathrm{drained}}=0.27$$ | (11) |

$${u}_{i}({t}_{1})-{u}_{i}({t}_{\mathrm{q}}^{+})=\mathit{okada}({\nu}_{\mathrm{drained}})-\mathit{okada}({\nu}_{\mathrm{undrained}})$$ | (12) |

**U**are identical in both terms. In other words, this model approximates the post-seismic displacement generated by poro-elastic rebound as the difference between two co-seismic models with different values of Poisson's ratio. Since the undrained Poisson's ratio is generally higher than the drained value, the post-seismic displacement vectors will point in the opposite direction from the co-seismic ones. In map view, Fig. 6 shows that the horizontal displacements produced by the poro-elastic model are similar in magnitude to those produced by the after-slip model. The vertical displacements, however, are opposite in sign. For displacement fields recorded by satellite images, the horizontal and vertical components combine to form an asymmetric range change field. In the Landers case, for example, the projections of the horizontal and vertical components onto the radar line of sight have opposite sign on the east side of the fault, but the same sign on the west side [36].

To approximate the displacement field produced by post-seismic deformation of a poroelastic medium, the Poisson recipe takes the difference of two elastic co-seismic models with different values of Poisson's ratio, as described above. This recipe approximates a poro-elastic medium (solid rock and pore fluids) described by four rheological parameters as an elastic medium (analogous to a rubber eraser) described by two elastic parameters. This approximation makes several strong assumptions: (1) that the distribution of co-seismic slip is correct, (2) that the change in Poisson's ratio is spatially uniform, and (3) that the Poisson's ratio changes as a known function of time. At the time of this writing, many modeling studies are focusing on the validity of these assumptions [6,7,14,21,22,36,42,49,52,53,56,67–69,83,93,95,96,102,109,115].

The consequences of this process are challenging to describe, because they depend on local permeability and porosity, which can vary greatly according with rock type, fissure density, porosity and stress. With many simplifying assumptions, one can relate the relaxation time observed in geodetic time series to the hydrologic permeability k. For a fault with width $\sim 1\text{}\text{km}$ and a relaxation time of ∼1 month, the permeability should be on the order of $k\sim {10}^{\mathrm{-18}}\text{}{\text{m}}^{2}$ [21]. A more complete treatment, however, requires describing the material parameters with four parameters: two for the solid matrix plus two for the pore fluid [53,116].

The co-seismic perturbation to the hydrologic pressure field also alters the stress conditions on faults near the mainshock, triggering aftershocks where increasing pressure unclamps faults near failure [78]. This mechanism has apparently triggered aftershocks in Iceland [5] and the Afar [76]. By extending the theory to incorporate time-dependent rheology, we should be able to calculate the location of aftershocks triggered by the poro-elastic mechanism of post-seismic relaxation. This approach might also explain changes in the ${V}_{\mathrm{p}}/{V}_{\mathrm{s}}$ ratio observed in the 50 days following the $Mw\text{}8.0$ Antofogasta earthquake in Chili in 1995 [56].

### 2.5 Fault-zone collapse recipe for the sponge analog

Another recipe for describing sponge-like deformation involves a change in volume in the fault zone. For example, Savage et al. [100] proposed a model of ‘fault collapse’ for the 1989 Loma Prieta earthquake. In this case, the source of the post-seismic deformation is modeled as a dislocation with a (negative) tensile component that reduces the volume of a rectangular prism or dike (Fig. 7). One possible physical explanation for this process is the change in widths of cracks and fissures. After opening during the co-seismic rupture, they may slowly heal over the post-seismic interval, with an exponential relaxation time $\tau \sim 5\text{}\text{yr}$, as suggested by observations of seismic scattering [8]. The spatial pattern of deformation produced by such a process is shown in Fig. 8. Like the viscoelastic models with a weak substrate (see below), it predicts horizontal displacements perpendicular to the fault plane. Like the Poisson recipe, it predicts short-wavelength features near the fault.

Later, Arnadottir et al. [4], concluded that the “geodetic data do not place useful constraints on the amount of dilatancy” for the Loma Prieta event, a conclusion reaffirmed by [16]. However, the fault zone collapse mechanism has not been definitively rejected and has also been proposed for the Landers earthquake [66].

### 2.6 Honey analog

The third analog considers the lower crust and/or upper mantle as a weak substrate that can flow smoothly. Its rheological behavior is analogous to that of a spoon in honey – movement of the spoon results in flow and permanent, non-recoverable deformation of the honey.

The simplest recipe for describing such flow is a viscous rheology. The stress depends on the strain rate, and if the dependence is linear, the coefficient of proportionality is called the viscosity. Viscosity varies greatly on Earth: $\sim {10}^{\mathrm{-5}}\text{}\text{Pa}\phantom{\rule{0.2em}{0ex}}\text{s}$ in air, $\sim {10}^{\mathrm{-3}}\text{}\text{Pa}\phantom{\rule{0.2em}{0ex}}\text{s}$ in liquid water, $\sim {10}^{9}\text{}\text{Pa}\phantom{\rule{0.2em}{0ex}}\text{s}$ in glass softened for blowing, and 10^{18} to ${10}^{20}\text{}\text{Pa}\phantom{\rule{0.2em}{0ex}}\text{s}$ (1 to 100 EPa), in the asthenosphere.

A viscoelastic medium combines viscous and elastic rheologies. There are several possible mechanical combinations that have been applied to study the rheology of Earth's lithosphere. The behavior of a Maxwell solid is analogous to that of a spring and a dashpot connected in series [91]. The spring mimics (recoverable) elastic deformation in the upper crust with an elastic rigidity μ while the dashpot mimics (irrecoverable) viscous flow in the lower crust with a (Newtonian) linear viscosity η. To first order, the displacement signal relaxes exponentially in time over a characteristic ‘Maxwell’ time ${\tau}_{\mathrm{m}}$:

$${\tau}_{\mathrm{m}}=\frac{\eta}{\mu}$$ | (13) |

Estimating the viscosity for rocks at the temperature and pressure conditions prevailing in the Earth is quite challenging. The earthquake deformation cycle is a good rheological experiment because we can measure both the impulse (the co-seismic rupture) and the response (the post-seismic deformation). Other ways of estimating the viscosity in rocks include: laboratory experiments [55], post-glacial rebound [57], gravity studies of lithospheric flexure caused by topographic loading [118] or lacustrine loading [11]. Of course, these estimates of viscosity depend on the assumptions in the models.

### 2.7 Viscoelastic recipe for honey analog (Fig. 9)

Viscoelastic models for the earthquake cycle generally consist of an upper elastic layer with finite thickness where seismic faulting occurs and one (or more) viscoelastic layer(s) that underlie it. To first order, the viscoelastic model describes the temporal evolution of post-seismic displacement as an exponential function:

$${f}_{\mathrm{exp}}(t)=1-\mathrm{exp}(-\frac{t-{t}_{\mathrm{q}}}{{\tau}_{\mathrm{m}}})$$ | (14) |

^{17}to ${10}^{20}\text{}\text{Pa}\phantom{\rule{0.2em}{0ex}}\text{s}$, we expect post-seismic displacements driven by viscoelastic processes to decay over a Maxwell time on the order of 1 to 1000 yr. To implement such a viscoelastic rheology for a spherical Earth model, one can sum an infinite series of decaying exponentials (e.g., [89]). If the Maxwell time for the different modes depends on location, however, this formulation will contradict our previous assumption that the time dependence and spatial distribution functions are separable.

The viscoelastic post-seismic displacement field has several distinctive features. For the horizontal components, the co-seismic and post-seismic displacements have the same sign if the viscoelastic source is fairly deep. For the vertical displacements, the co-seismic and post-seismic displacements have opposite signs. This model also tends to produce large displacements perpendicular to the along-strike extension of the fault [102]. The displacement field produced by the viscoelastic relaxation mechanism decays with distance from the fault over length scales on the order of the elastic layer thickness, often taken to be the thickness of part or all of the crust, roughly 10 to 30 km (Fig. 10).

Estimating viscosity at depth from geodetic surveys at the surface is difficult because the latter are not very sensitive to the former. Moreover, the trade-offs between depth and viscosity make comparisons between different studies and different areas difficult. Nonetheless, most studies suggest viscosities on the order of 1 to 10 EPa s in the lower crust and/or upper mantle beneath the continents. It seems likely, however, that the weak, flowing substrate obeys a non-linear constitutive relation. For example, aseismic creep on particular sections of some strike-slip faults seems to deform according to a power-law rheology [103,104]. More recently, several studies have suggested power-law rheologies to explain post-seismic deformation [38,90] (Fig. 11).

## 3 Tricks of the modeling trade

### 3.1 Backslip

Savage and Prescott [99] introduced an ingenious shortcut into modeling the surface deformation effects of steady (i.e. constant rate) aseismic fault slip on the downdip extension of the co-seismic fault in elastic and viscoelastic Earth models. They reasoned that the superposition of ‘backwards’ slip (‘backslip’) at constant rate on the shallow co-seismic fault in a sense opposite to co-seismic slip with constant-rate ‘forward’ block motion across the entire fault will produce the same deformation at the Earth's surface as that of a fault locked on its co-seismic segment and freely sliding below. Since the block motion produces no internal strain, the model deformation is due only to the ‘backslip’ element of the model. One calculation then serves to obtain the elastic or viscoelastic effects of any desired distribution of buried steady fault slip, and ‘backslip’ models have become a very popular way of modeling inter-seismic deformation measured by GPS methods, particularly at subducting plate boundaries (e.g., [63]). Johnson and Segall [51] introduced a slightly different formulation using incremental slip that they compare to the standard backslip approach in a helpful sketch (their Fig. 6).

### 3.2 Distinguishing buried slip from viscoelastic relaxation

Some form of viscoelastic coupling of the largely elastic seismogenic crust and the underlying ductile lower crust or upper mantle must occur and is potentially crucial to understand the earthquake deformation cycle. However, distinguishing the transient buried slip from viscoelastic relaxation from geodetic observations has always been difficult because the observable surface deformation fields produced by the two mechanism are often very similar [111,113]. Indeed, Savage [98] has shown that for a two-dimensional vertical strike-slip fault, the two models are formally equivalent. In this case, an appropriate distribution of elastic dislocations on the fault can superposed to exactly mimic the surface deformation due to relaxation of a viscoelastic half-space lying beneath a faulted elastic plate. For a three-dimensional strike-slip fault this equivalence breaks down: the end effects are of opposite sign for the two models. This diagnostic feature can distinguish the afterslip mechanism from viscoelastic relaxation using INSAR data gathered after the 1999 Hector Mines earthquake [90]. For dip-slip faulting, the two mechanisms are not formally equivalent, even for two-dimensional faults.

## 4 Synthesis and conclusions

We have considered different mechanisms for explaining post-seismic deformation and arranged them into three categories analogous to a rubber eraser, a sponge, and a spoon in honey. Distinguishing one class of mechanism from another requires considering the temporal evolution and the spatial distribution of the displacement fields they produce. If we assume that the time- and space-dependence of post-seismic deformation are separable, we can describe the observed time series and maps with fairly simple mathematical expressions. The parameters in these expressions vary greatly, however, making it challenging to reject one model in favor of another. Furthermore, more than one process seems to be responsible for producing the post-seismic deformation observed in the seven years following the 1992 Landers earthquake in California. Other earthquakes also suggest different mechanisms acting on different, possibly overlapping, scales in time and space.

Although challenging, the post-seismic problem is well worth pursuing because it represents one of the few rheological experiments that is possible to perform both in situ and in a human lifetime. In the future, we expect progress to come in the following areas:

Robust estimation procedures. We have seen that the parameters in the various models of post-seismic deformation trade off with each other. For example, co-seismic offset will trade off with the characteristic time in a logarithmic function used to describe a time-series of station position. Similarly, estimating viscosity at depth from crustal deformation data measured at the Earth's surface is challenging because the latter is not very sensitive to the former. To obtain robust estimates of these correlated parameters requires estimating them simultaneously in a single optimization procedure. The usual multi-step process of ‘correcting’ the time series for co-seismic offsets and then the inter-seismic rate before curve-fitting the post-seismic part of the time series will bias the estimates of all parameters. Furthermore, if we can assume time–space separability, then we should constrain the temporal evolution function $f(t)$ to be the same at all measurement points. Toward this end, we prefer estimation procedures that “smooth” both the temporal evolution function $f(t)$ and the spatial distribution map $\mathbf{G}(\mathbf{x})$ such as GLOBK [47], Network Inversion Filter [104], or Permanent Scatterers [34,35] approaches.

Models that incorporate the geometric complexity of real faults. A normal fault that shallows with depth, as Meyer et al. [69,70] suggest for Grevena, is difficult to approximate with the simple model of a single, planar fault patch. Yet the latter approximation is the most feasible geometric parameterization in the non-linear inverse problem of estimating the focal mechanism [19,20]. Rigo et al. [94] resolved most of this discrepancy by applying a model formulation that admits a non-planar fault surface [64,114]. An open-source computer code for this type of model would foster progress.

More continuous GPS stations. We have seen that time series of geodetic measurements can help to distinguish between various mechanisms for post-seismic deformation. As continuously operating GPS (CGPS) stations become more ubiquitous, they will capture non-linear time series after more earthquakes. The density of CGPS networks in California and Japan (1 station per 1000 km^{2}) should be extended to other fault zones such as the South American cordillera, the boundaries of the Tibetan plateau, and the southwest Pacific. When a large (${M}_{w}>7$) earthquake occurs in an area without dense CGPS coverage, the priority should be installing CGPS instruments in a profile perpendicular to the strike of the co-seismic rupture. For example, if ten instruments are available, five should be installed on each side of the fault with distance of 10 km between them [13,45,50]. If possible, benchmarks near the fault surveyed by campaigns before the earthquake should be equipped with continuously operating instruments. Placing GPS receivers and seismometers in the same locations would save considerable effort.

Routine application of INSAR. To date, several factors make INSAR measurements of transient deformation a hit-or-miss, opportunistic affair. First, the lack of closely spaced orbital trajectories can force a compromise between short time spans (good correlation) and small orbital separation. Second, capturing earthquakes with INSAR is a major challenge because we do not know where they will occur. This implies that each INSAR-capable satellite must acquire a catalog of pre-quake images over all the land areas likely to produce a measurable earthquake. We estimate this area to be approximately 66 million km^{2}. A SAR mission dedicated to interferometric measurements of earthquake faulting is needed.

Related phenomena. INSAR and CGPS open two new windows in the spatio-temporal spectrum of seismological metrology: INSAR at distance scales between ∼1 and $\sim 10\text{}\text{km}$; CGPS at time scales of days to years. Prior to the introduction of these two techniques, measurements at these scales were prohibitively expensive or prone to drift. Now that both techniques have entered the realm of operational, routine observations, we should expect to see interesting observations of other seismological phenomena, such as slow earthquakes, inter-seismic strain accumulation, and perhaps even an earthquake precursor.

## Acknowledgements

We thank Jim Savage and Fred Pollitz for critical reviews that significantly improved the manuscript. We also thank Thora Arnadottir, Freysteinn Sigmundsson, Fred Pollitz, Loic Dubois, Dimitri Komatitsch, Hugo Perfettini, Michel Rabinowicz, Brad Hager and Semih Ergintav for helpful discussions. Yuri Fialko, Fred Pollitz, Hugo Perfettini, and Eric Hetland kindly provided results in advance of publication. This work was partially financed by GDR STRAINSAR and ANR.