Omori‐like decay of postseismic velocities following continental earthquakes

Various mechanisms have been proposed to explain the transient, enhanced surface deformation rates following earthquakes. Unfortunately, these different mechanisms can produce very similar surface deformation patterns leading to difficulty in distinguishing between them. Here we return to the observations themselves and compile near‐field postseismic velocity measurements following moderate to large continental earthquakes. We find that these velocities have a remarkably consistent pattern, with velocity inversely proportional to time since the earthquake. This suggests that postseismic velocities show an Omori‐like decay and that postseismic displacements increase logarithmically over time. These observations are inconsistent with simple, linear Maxwell or Burgers body viscoelastic relaxation mechanisms but are consistent with rate‐and‐state frictional afterslip models and power law shear zone models. The results imply that near‐field postseismic surface deformation measurements are primarily the result of fault zone processes and, therefore, that the inference of lower crustal viscosities from near‐field postseismic deformation requires care.


Introduction
The rheology of the continental lithosphere remains poorly understood, with a number of different disciplines contributing sometimes conflicting observations. Glacial Isostatic Adjustment (GIA) studies in continental cratons suggest that the lithosphere has very high viscosities (>10 22 Pa s) or behaves elastically [e.g., Peltier and Drummond, 2008;Zhao et al., 2012] and generally resolve deep, mantle relaxation with viscosities on the order of 10 20 -10 21 Pa s. GIA studies examining plate boundary zones find much thinner elastic lithospheric thicknesses and mantle viscosities in the range 10 18 -10 19 Pa s [James et al., 2000;Ivins et al., 2011]. On shorter timescales, surface deformation following the draining of lakes [e.g., England et al., 2013] and strain concentration associated with interseismic deformation [e.g., Johnson et al., 2007] also suggest relatively high lithospheric viscosities in the range 10 19 -10 22 Pa s. Viscosity estimates obtained from early postseismic deformation are typically several orders of magnitude lower than those obtained from other disciplines (≲10 18 Pa s) [Watts et al., 2013]. Apparent viscosities are often seen to increase with time since the earthquake leading to inferences of various transient rheologies [e.g., Pollitz, 2003;Ryder et al., 2007;Freed et al., 2010]. Whether these inferred viscosities apply to the lower crust, upper mantle or simply the fault zone itself is not always clear [Bürgmann and Dresen, 2008;Wright et al., 2013]. Furthermore, the observed surface deformation can also be explained by ongoing fault slip (afterslip) and it is often challenging to distinguish between these competing hypotheses [e.g., Savage, 1990;Perfettini and Avouac, 2004;Ryder et al., 2007;Hao et al., 2012;Wright et al., 2013]. Despite these common challenges, it has been possible to determine the predominant postseismic deformation mechanism in a selection of cases [e.g., Jónsson et al., 2003;Freed, 2007;Freed et al., 2007;Jónsson, 2008].
Many studies do not test multiple models against their postseismic data, and so evaluating the relative importance of different deformation mechanisms is challenging. Different authors often reach very different conclusions about the lithosphere using similar observations from the same earthquake [Wright et al., 2013]. For example, postseismic deformation following the 1992 Landers earthquake has been attributed to deep afterslip [Savage and Svarc, 1997], poroelastic rebound and deep afterslip [Fialko, 2004], afterslip and fault zone contraction [Massonnet et al., 1996], lower crustal relaxation [Deng et al., 1998], and power law mantle flow [Freed and Bürgmann, 2004].
Here we compile postseismic observations from 34 moderate to large, continental, intracrustal earthquakes/earthquake sequences. The majority of these observations come from data sources within 25 km 10.1002/2017GL072865 of the coseismic fault. We show that the temporal behavior over timescales from hours to tens of years is remarkably consistent. The results support models in which rate-and-state frictional afterslip and/or viscoelastic relaxation of a power law shear zone are the primary causes of postseismic deformation. We discuss the implications for the rheology of the continental lithosphere.

Data Collection
We compiled 151 postseismic velocity data points from 45 publications. These studies cover 34 continental earthquakes/earthquake sequences. The use of the full 3-D deformation field through time provides the most information about postseismic processes and rheology. However, trying to combine this information from multiple earthquakes is challenging. Instead, we combine information about the temporal variation of velocities from multiple earthquakes. In order to collect a consistent data set across different earthquakes and time intervals, we extracted the maximum surface velocities observed in the studied time interval (see Table S1 in the supporting information). If displacements are given at different time intervals for the same earthquake, we use the same observation point (where possible), and this criterion supersedes the criterion of picking the maximum amplitude signal. In other words, if the location of maximal displacement moves with time, we do not follow this and instead extract measurements from the same location through time.
The distance of the maximum deforming surface location from the coseismic fault is a reasonable proxy for the depth range of maximum deformation. Most of the maximum surface velocities come from points in the near-field (Figure 2), which will primarily be affected by near-field, shallow deformation. However, some points are at slightly larger distances; this suggests deeper deformation (e.g., in the lower crust) is dominant.
A variety of measurement techniques were used in the original studies, primarily Global Navigation Satellite Systems (GNSS) (e.g., GPS), interferometric synthetic aperture radar (InSAR), and leveling. GNSS is routinely used to study postseismic deformation; GNSS measurements account for 86 of our 151 data points. In general, horizontal displacements/velocities are reported and we extracted the maximum values, excluding any clear outliers. Interferometric synthetic aperture radar (InSAR) has provided 55 of our 151 data points. InSAR measures the change in distance (range) between the satellite and the ground, equivalent to the projection of the 3-D displacement vector in the satellite line-of-sight (LOS) direction. We chose to extract the maximum postseismic LOS displacement change for each earthquake. This approach may underestimate postseismic deformation if the predominant deformation direction is at a large angle to the satellite look direction. In practice, InSAR is most sensitive to vertical and east-west displacements. Leveling data records changes in height at particular points along transects and is the source of 10 data points in our compilation. The method allows dense spatial coverage along transects but generally has sparse temporal coverage. Nevertheless, it allows us to study postseismic deformation from old earthquakes, for example, from the 1940 Imperial Valley earthquake [Reilinger, 1984]. This method is insensitive to horizontal motions, which is problematic for predominantly strike-slip earthquakes.
Isolating a postseismic signal in actively deforming regions of the earth requires the removal of other signals (e.g., interseismic deformation). Most of the studies in our compilation achieve this using data prior to the earthquake or well-established models of other deformation sources. We note those studies where this signal separation is not possible in Table S1 [e.g., Calais et al., 2002;Vergnolle, 2003]. The earliest velocity observations included in our compilation are 3.5 h after an earthquake and the latest are 91.5 years after. Postseismic velocities range from over 10,000 mm/yr immediately following an earthquake to less than 1 mm/yr decades later. Using a log-log plot allows us to show the data clearly from all timescales (Figure 1). We connect data points from earthquakes where multiple postseismic deformation rates from different times after the earthquake were available.
In order to examine postseismic velocities at very early times we use data from GPS stations POMM, LAND, and HUNT following the 2004 Parkfield earthquake [Langbein, 2006]. We use the 30 min solutions provided by Langbein [2006] to examine postseismic velocities very shortly after an earthquake. There is significant scatter in the individual 30 min solutions, but a clear trend is visible in the displacement time series. We perform a piecewise linear fit to the data and extract the postseismic surface velocity over different intervals following the earthquake.
Approximately 100 years after an earthquake, the postseismic signal is only just detectable with current measurement techniques. The latest postseismic velocity estimates we have are those from earthquakes in  Figure 1a. (c) Normalized results. Light grey circles represent the original data. Blue circles represent the data normalized using intercepts of linear trends through each earthquake time series (see text). Red line shows the best fitting Omori model and black solid lines show the region where 95% of Omori models plot. Data sources: Barbot et al. [2008], Atzori et al. [2008], Reddy et al. [2012], Calais et al. [2002], Hsu et al. [2002], Perfettini and Avouac [2004] Vergnolle, 2003]. These estimates are dependent on models to deconvolve secular, long-lived deformation and transient, postseismic velocities and as such do not represent "raw" observations. Despite this difficulty, the postseismic deformation estimates from these points agree with data from earlier in the postseismic phase.
The data show a remarkably simple trend in log-log space, with a gradient of approximately −1. Postseismic surface velocity (V) is inversely proportional to the time (t) since the earthquake. This trend is evident in individual earthquakes as well as for the whole data set. With these results, we can estimate the range of possible postseismic surface velocities at a given time after an earthquake, assuming a 1∕t relationship. Ninety-five percent of our postseismic velocities lie between the blue dashed lines in Figure 1, with a factor of approximately 55 separating the lower velocities and higher velocities at any time. This information may be particularly useful for targeting future studies of long-lived postseismic deformation or deploying GNSS stations.

Normalized Results
There is some scatter in the data which can be attributed to a number of factors. These earthquakes have different magnitudes, depths, and mechanisms, have occurred in different parts of the world, and have been observed using different techniques at different distances from the coseismic fault. However, given the large number of variables which could control postseismic surface velocities, it is surprising how little scatter there is in the data.
A linear regression on each of the individual earthquakes where three or more postseismic velocity measurements were available gave gradients with a mean of −0.96 and a standard deviation of 0.24. Coefficient of determination (R 2 ) values for these fits were generally high, with most (16 out of 22) higher than 0.9. The slope of the relationship is remarkably similar for all earthquakes, but there is a scale factor that affects the magnitude of postseismic velocities. To remove this effect we normalized each individual earthquake time series. The normalization was based on the linear regression such that the linear fit for each earthquake had a value of 12.16 mm/yr (log V = 1.09) at a time of 1 year (log t = 0); this normalized the results to the mean linear fit to all the data. The normalized results show significant reduction in scatter and have a gradient of −0.92 ± 0.13 (95% confidence interval from Bayesian inversion) close to the average gradient of the individual earthquakes ( Figure 1).

Scaling Factors
We investigate whether these postseismic velocities are dependent on a number of physical properties of the earthquake (moment, depth, and focal mechanism) or are measuring artifacts (distance of the observation from the fault, measuring technique) ( Figure 2). Earthquake stress drop is expected to play a significant role in determining the magnitude of any postseismic response. Larger stress drops result in greater stress transfer to neighboring regions of the fault/lithosphere. However, the surface deformation resulting from regions releasing stress depends on how large these regions are (determined by magnitude), and where they are in the lithosphere (determined by depth and focal mechanism). Stress drops are not easily accessible or generally well resolved for all earthquakes in our compilation, and so we do not include this parameter in our analysis. We performed a Bayesian inversion to determine the linear regression coefficients and their errors for magnitude, distance from fault, centroid moment tensor (CMT) depth, and time since the earthquake. Time exerted the strongest control, with magnitude and distance from fault playing secondary roles. We found that larger earthquakes generally produced greater postseismic velocities and that velocities increased with distance from the fault (for distances up to 35 km). CMT depth appeared to play no significant role, but this may be due to large errors in the reported CMT depths.
On a larger scale, there may be scatter due to lateral rheological heterogeneity over the Earth's surface. To test this, we plot each earthquake's power law decay coefficient at the location of the earthquake (see Figures S1a and S1b). We find no clear patterns, although higher decay coefficients are seen in the Tibetan plateau, suggesting that postseismic relaxation may occur more rapidly there. We also examined whether the best fit velocity gradient or intercept varies with time after an earthquake by using the fits to individual earthquakes (see Figures S1c and S1d). There was no clear pattern, with both velocity gradients and intercepts seemingly independent of time since the earthquake.

Consistency With Omori's Law
Our data compilation suggests that postseismic velocities can be fit using a version of Omori's law. Utsu [1957] developed the modified Omori's law which describes the number of aftershocks, n, expected at any time, t, after an earthquake: Our postseismic velocity compilation can be fit using an equation of the same form, where n(t) is replaced with V(t). We find that p values range between 0.80 and 1.04 (at 95% confidence).

Modeling
We use the temporal constraints provided by these observations to test postseismic relaxation models for the crust/upper mantle and fault zones. These models are all based on viscoelastic relaxation or afterslip. We do not consider poroelastic rebound here since this effect is only expected to last a short time following an earthquake [Jónsson et al., 2003;Fialko, 2004;Wright et al., 2013].

Linear Maxwell Models
The simplest rheological model usually invoked to explain postseismic deformation is the linear Maxwell model. A Maxwell material can be conceptualized as an elastic spring and viscous dash pot in series. In an earthquake there is an instantaneous elastic response followed by a decaying response through time as the material relaxes. Postseismic velocity in a Maxwell material will decay exponentially as where V 0 and are constants [e.g., Montėsi, 2004]. The Maxwell relaxation time of the material ( ) is equal to the ratio of viscosity ( ) to shear modulus ( ), and the instantaneous velocity (V 0 ) is inversely proportional to the relaxation time. Although instantaneous postseismic velocities at any time can be matched by this linear Maxwell model, it cannot explain the temporal decay. Low viscosities are required to explain rapid early motions while higher viscosities would be required to explain sustained slow motion (Figure 3).
A slightly more complicated linear rheology is a Burgers body. The Burgers body has a Maxwell material in series with a Kelvin (Voigt) material and shows two relaxation times ( ) [e.g., Pollitz, 2003Pollitz, , 2005. A transient phase of deformation is observed as the Kelvin material relaxes which is superimposed on the longer timescale Maxwell relaxation. Postseismic motion in a Burgers rheology can be described using where i = i i , V 0 is a constant initial velocity and K is a scalar [Malkin and Isayev, 2012]. The curves produced by these analytical expressions have the same form as those produced using VISCO1D [Pollitz, 1992], suggesting that the use of these expressions is valid. We perform a Bayesian inversion to obtain the maximum likelihood solution for a Burgers body (Figure 3). We find that a 1-D Burgers body model can produce results within the scatter of the data. The acceptable ratios of relaxation times are shown in a histogram in Figure 3. If both relaxing elements have a similar shear modulus, then the ratio in relaxation times can be seen as a ratio in viscosities. We find viscosity ratios of approximately 100 best fit the data. This is larger than the typical ratios found in previous studies [Ryder et al., 2011, and references therein] and is likely due to the need to explain a greater time range of observations simultaneously. Overall, though, the model prediction is inconsistent with the linear log-log trends observed in both the overall data compilation and in individual earthquakes.
Linear Maxwell models and Burgers body models cannot adequately reproduce the temporal evolution observed for postseismic velocities over extended time periods. The addition of more viscoelastic relaxation elements (e.g., a continuum of viscosity values) is able to better reproduce the observed trends [Hetland and Hager, 2006]. This viscosity variation may also be in space, for example, with depth [Riva and Govers, 2009;Yamasaki and Houseman, 2012] or distance from the fault zone [e.g., Yamasaki et al., 2014].

Afterslip Models
The rate-and-state friction law [Dieterich, 1979] has been widely used to explain a number of fault-related phenomena including earthquakes, slow slip events, steady creep, and afterslip. A number of authors have applied the rate-and-state friction law to postseismic afterslip [Marone, 1998;Hearn, 2002;Barbot et al., 2009]. Rate-and-state friction is often simplified by assuming that the fault is close to the steady state regime. Close to the steady state regime, the state variable is constant in time (̇= 0) and there is no loading stress rate. Under these conditions Marone [1998] showed that postseismic velocities are where is a characteristic decay time. In this regime velocity is seen to decay with time following a 1∕t relationship from an initial value V 0 (Figure 3). This is identical to the modified Omori's law formulation (equation (1)) if c = , K = V 0 and p = 1. The model reproduces the overall temporal decay well. Montėsi [2004] investigated postseismic deformation in a power law shear zone and derived a general relaxation law describing postseismic relaxation:

Shear Zone Models
where n is the power law exponent. Nonlinear rheologies are those where strain rate is proportional to stress raised to the value of the power law exponent.
This relaxation law allows us to test various shear zone rheologies since the surface velocity is simply expected to be directly proportional to the shear zone velocity. The law simplifies to two end members. When 1∕n → 1, the general law tends toward the following: This equation defines Newtonian viscous flow in a shear zone and is identical to the linear Maxwell result (equation (2)). This is shown as the blue line in Figure 3d. When n ≫ 1 and therefore 1∕n → 0, we obtain the equation for rate-and-state afterslip (equation (4)) [Montėsi, 2004]. From postseismic deformation alone, it is impossible to distinguish between postseismic relaxation of a shear zone with high n and frictional afterslip (see red line in Figure 3d). These models are also capable of explaining the overall pattern seen in our postseismic deformation data, provided the power law exponent is greater than the usual range of experimentally determined values.

Agreement With Common Postseismic Observations
Our compilation of observations is in agreement with a number of common postseismic observations which we outline below. First, postseismic surface displacement time series are commonly fit using logarithmic equations [e.g., Donnellan and Lyzenga, 1998;Freed et al., 2006a;Mahsas et al., 2008;Dogan et al., 2014]. These equations take the general form where C is a constant and is a time constant. If this equation is differentiated then we obtain the equation predicted by rate-and-state afterslip and/or relaxation of a high-n shear zone (i.e., V ∝ 1∕t for high t). Our data also show that postseismic deformation takes an identical functional form to Omori's law describing the decay of aftershocks (equation (1)). Typical p values lie in the range 0.6-2.5 with a median of about 1.1 [Utsu et al., 1995]. Our value for the decay of postseismic surface deformation through time has a p value of about 1 which shows that both postearthquake processes decay at similar rates. Aftershock sequences for continental earthquakes can last for decades to centuries [e.g., Ryall, 1977;Utsu et al., 1995;Ebel et al., 2000;Stein and Liu, 2009], which is in agreement with the long postseismic deformation timescales included here. Many have noted the similar decay rates of aftershocks and afterslip velocity, leading to the proposition that aftershocks are primarily controlled by frictional afterslip [e.g., Perfettini and Avouac, 2004;Savage et al., 2007;Helmstetter and Shaw, 2009].
When postseismic relaxation is explained using a viscoelastic material, many authors have noted an apparent increase in crustal/mantle viscosity through time following earthquakes [e.g., Pollitz, 2003;Freed et al., 2006a;Ryder et al., 2007]. Our results support this observation and show that this trend continues over a period of 100 years (Figure 3). Maxwell viscoelastic models require an increase in effective viscosity through time in order to match the decreasing postseismic velocities. This changing effective viscosity is sometimes modeled using multiple layers with different viscosity [e.g., Jónsson, 2008;Hetland and Hager, 2006], a transient rheology like a Burgers body or standard linear solid [e.g., Pollitz, 2003Pollitz, , 2005Ryder et al., 2007Ryder et al., , 2011, power law rheologies [e.g., Freed and Bürgmann, 2004;Freed et al., 2006b], or a combination of the above. Our compilation shows that any successful model must have a continuously changing viscosity, e.g., power law models with high n or linear models with a large number of relaxation times.

Temporal Variation and Model Characteristics
We have sampled postseismic deformation over a long time period: from a number of hours after an earthquake up to almost 100 years later. All these data define a linear trend in log-log space with a gradient close to −1. This observation suggests that all continental earthquakes exhibit similar temporal postseismic deformation patterns and allows us to give bounds for expected postseismic velocities at any given time after an earthquake.
Univiscous linear Maxwell materials and Burgers bodies with just two relaxation times cannot reproduce the linear trends seen in our compilation. Viscoelastic models need to be more complex with variations of viscosity in time, space, or both. In sufficiently complex rheologies, the surface displacements can be described by a logarithmic function, which when differentiated will produce a 1∕t relationship [Hetland and Hager, 2006].
Power law creep can reproduce linear trends in log-log space as required by our data. This model is also supported by rock experimental results, but the exponents found in these studies are usually between 2 and 5 [Carter and Tsenn, 1987;Hirth and Kohlstedt, 2003;Freed and Bürgmann, 2004]. Our compilation requires higher power law exponents and suggest that either the experimental results cannot be scaled up to tectonic conditions or more likely, that afterslip is the dominant mechanism recorded by our data.

Localized Deformation
We find that the temporal variation of our near-field compilation is best explained by afterslip or power law creep at high n. Afterslip is a very localized deformation mechanism, involving continued slip on the fault plane. Moore and Parsons [2015] found that a power law rheology contributed to the narrowing of viscous shear zones. They found that narrow shear zones would develop in viscous materials where viscosity varies with depth, with shear stress heating further narrowing the shear zone. Our data are consistent with a localization of postseismic deformation as aseismic afterslip in the upper crust and in deep shear zones in the middle to lower crust.
Fault zone-related processes best explaining our compilation is perhaps not surprising, since our surface velocity measurements are mostly from within a few kilometers of the fault zone ( Figure 2). These length scales are certainly small enough to be affected by processes centered on the fault zone [e.g., Freed, 2007] and may be dominated by fault zone processes rather than lithospheric relaxation. Furthermore, other strands of evidence support localized deformation in fault zones. Geodesy reveals that the majority of continental fault zones show significant strain localization between earthquakes [Wright et al., 2013;Vernant, 2015]. Geological evidence from exhumed roots of faults shows that motion at depth is likely localized into shear zones up to a few kilometers wide [Hanmer, 1988;Norris and Cooper, 2003;Vauchez and Tommasi, 2003;Frost et al., 2011]. Seismic experiments have also shown deep narrow structures along the North Anatolian Fault in Turkey [Fichtner et al., 2013;Kahraman et al., 2015;Taylor et al., 2016], the San Andreas [Zhu, 2000], and dip-slip faults in Tibet .
Both afterslip models and power law creep models find support from other observations. Afterslip models can explain the temporal decay in aftershock frequency [Helmstetter and Shaw, 2009]. The rate-and-state friction laws which form the basis for the afterslip models used here can also explain a large number of other aspects of the seismic cycle [Dieterich, 1994;Marone, 1998;Liu and Rice, 2005;Helmstetter and Shaw, 2009]. At high temperatures and stresses, rocks deform by power law creep in laboratory experiments [Wilks and Carter, 1990;Kohlstedt et al., 1995;Montėsi and Hirth, 2003;Bürgmann and Dresen, 2008]. These conditions may be prevalent in the deeper portions of fault zones.
Fault zone relaxation processes are usually considered to be relatively short-lived, but our data spans decades of the postseismic period. Afterslip has been observed in a few examples decades after earthquakes [Reilinger, 1984;Kaneko et al., 2013;Copley and Reynolds, 2014;Copley, 2014] and often is not tested for on these long timescales. Studies examining postseismic deformation decades after an earthquake should consider the role of continued afterslip, especially for explaining near-fault observations.
While afterslip/power law creep can explain the temporal variation seen in our data, it is unclear whether it can explain the time varying spatial patterns of postseismic deformation. For example, immediately after an earthquake, poroelastic rebound may play a significant role in determining the spatial pattern of postseismic deformation [Joünsson et al., 2003;Fialko, 2004] in the near field. Other studies have suggested that broad viscoelastic relaxation of the mantle is required to match far field observations, for example, after the 2002 Denali (Alaska) earthquake [Freed et al., 2006a] or the Landers and Hector Mine earthquakes in California . While we argue that the temporal decay of postseismic deformation is a powerful discriminant between competing mechanisms, the spatial patterns of postseismic deformation have been enough to constrain the most important deformation mechanism in a selection of cases [e.g., Joünsson et al., 2003;Freed et al., 2006a;Freed, 2007;Copley et al., 2012].

Implications
Despite these caveats, our findings have some important implications. Our compilation suggests that fault zone processes (afterslip or high n shear zones) generate the largest near-field postseismic signals. These signals may dominate postseismic deformation fields for decades, particularly at near-field sites. As such, caution should be exercised when interpreting lower-crustal viscosities derived from postseismic studies using predominantly near-field data.  [Wessel and Smith, 1991]. We thank Andy Hooper and Gregory Houseman for useful discussions. We thank the Editor, Andrew Newman, and Roland Bürgmann and Sylvain Barbot for thoughtful reviews which significantly improved the manuscript. A table with the compilation of postseismic velocities can be found in the supporting information.