Kinematic Representations of Linear and Power‐Law Viscoelastic Deformation Through the Earthquake Cycle

Viscoelastic deformation below the Earth's elastic crust is modulated by stresses generated by both plate tectonic and earthquake cycle processes. Rapid near‐fault deformation following large earthquakes has been interpreted as the signature of viscoelastic stress diffusion in the upper mantle and earthquake cycle models have been developed that integrated this effect throughout the duration of the earthquake cycle. Here we develop a surrogate method to approximate the upper crustal kinematics associated with time‐dependent viscoelastic deformation that does not require knowledge of nor provide constraints on subcrustal rheology. To do this we show how the effects of deformation in the viscoelastic upper mantle can be emulated by introducing a set of effective dislocations at crust‐mantle interface. The approach is shown to approximate both linear and power‐law viscoelastic rheologies and offers a way to accurately represent viscoelastic kinematics even in the case where upper mantle rheology is poorly constrained.


Introduction
Earthquakes generate elastic stress changes over time scales associated with the propagation of seismic waves. This includes stresses that propagate from the largely elastic upper crust into the viscoelastic mantle below. In this lower layer, stresses may slowly diffuse and evolve due to viscoelastic processes. This process further exerts stresses on the base of the crust, which causes geodetically detectable motions both immediately following large earthquakes (e.g., Nur & Mavko, 1974;Pollitz et al., 2001) and throughout the earthquake cycle (e.g., Savage and Lisowski, 1998;Thatcher, 1983). The fundamental approach to modeling these motions is to develop mechanical models that explicitly represent the elastic crust and viscoelastic mantle ( Figure 1). For linear rheology cases (e.g., Maxwell, Burgers rheologies), solutions are generally based on exploiting the correspondence principle which allows for the solving an elastic problem in the Laplace domain that is equivalent to the viscoelastic problem in the time domain (e.g., Cohen & Kramer, 1984;Hetland & Hager, 2005, 2006Savage and Prescott, 1978). For the cases where non-linear rheologies are considered or where linear rheologies with lateral variations in material properties are considered, numerical methods are used (Freed, 2005;Freed & Bürgmann, 2004;Freed et al., 2010). Common to both the semi-analytic and numerical methods is the assumption that the physics of the problem can be reasonably specified. At a practical level, this requires being able to specify a rheology and its parameters (e.g., power-law exponent, effective viscosity). Applications have relied on a wide range of parameters and have been used to explain geodetic observations during the postseismic epoch in both two- (Ryder et al., 2007) and three-dimensions, (Wen et al., 2012) as well as to estimate material properties of the upper mantle (Pollitz, 2003).
Perturbation style approaches have also been advanced to model the viscoelastic effects throughout the earthquake cycle in three-dimensions (Devries et al., 2017;Hearn et al., 2013;Pollitz and Evans, 2017). These methods rely on the assumption of linear rheologies so that the viscoelastic effects can be superimposed independently. Here, we develop a kinematic approach to viscoelasticity that does not require knowledge of the of viscoelastic rheology. A significant goal is to enable the modeling of viscoelastic effects without having to know or be able to model the physics of viscoelastic process in detail so that other parameters such as fault slip rates and interseismic Abstract Viscoelastic deformation below the Earth's elastic crust is modulated by stresses generated by both plate tectonic and earthquake cycle processes. Rapid near-fault deformation following large earthquakes has been interpreted as the signature of viscoelastic stress diffusion in the upper mantle and earthquake cycle models have been developed that integrated this effect throughout the duration of the earthquake cycle. Here we develop a surrogate method to approximate the upper crustal kinematics associated with time-dependent viscoelastic deformation that does not require knowledge of nor provide constraints on subcrustal rheology. To do this we show how the effects of deformation in the viscoelastic upper mantle can be emulated by introducing a set of effective dislocations at crust-mantle interface. The approach is shown to approximate both linear and power-law viscoelastic rheologies and offers a way to accurately represent viscoelastic kinematics even in the case where upper mantle rheology is poorly constrained.

Plain Language Summary
The Earth's surface deforms at variable rates in the time between large earthquakes. A reason for this is that the motions within the earth, including its viscoelastic mantle, continue even when faults are locked. Here, we describe a method for representing the effects of viscoelastic motion that does not require detailed knowledge of viscoelastic properties. Additionally, we provide an application of this method to geodetic measurements of the interseismic deformation along the North Anatolian fault. coupling rates can be estimated. A corollary of this approach is that it provides no constraints on the rheology of the sub-crustal viscoelastic layer.

Kinematic Representations of Viscoelastic Deformation
Viscoelastic earthquake cycle models are predicated on explicitly modeling the interaction between what is assumed to be a finite thickness linear elastic upper crust coupled to a viscoelastic layer beneath it. For the simple class of non-inertial models that assume the periodic occurrence of earthquakes (e.g., Cohen & Kramer, 1984;Hetland & Hager, 2005;Savage and Prescott, 1978) stresses are generated that are instantaneously propagated into the viscoelastic layer. These stresses are then allowed to diffuse in the viscoelastic layer, while maintaining continuity of displacement and tractions at the interface with the elastic crust. This continuity constraint gives rise to the prediction of continued surface deformation in the immediate postseismic era and throughout the earthquake cycle.
The central idea in the kinematic representation argument presented here is that we can solve for an effective set of boundary conditions at the interface between the elastic and viscoelastic layers that mimic the effect of the explicitly mechanically coupled solution. To do this, we exploit the fact that displacements in the elastic layer are entirely a function of the applied boundary conditions. Specifically for the case where there are no body forces, the stresses and tractions at the boundaries of a linear elastic body can always be written as the result of a set of effective displacements. This is the fundamental approach in displacement-based boundary element modeling (e.g., Banerjee & Butterfield, 1981;Brebbia & Dominguez, 1994;Crouch & Starfield, 1983;Gaul et al., 2013) and the relationship between the effective displacements and physical displacements and tractions can be written as, where T*, and H*, are kernel matrices that relate effective displacements, at source coordinates, q, to displacements and tractions at observation coordinates, p, along a boundary C.
To demonstrate how this general principle can be applied to the specific earthquake cycle problem, we consider the example of an infinitely long strike-slip fault that ruptures completely through an elastic crust with a homogeneous viscoelastic layer below it. We consider the different cases for the lower crustal rheology: (a) Linear Maxwell, (b) power-law, n = 3, and power-law, n = 6 ( Figure 3). For each of these cases, we run forward viscoelastic earthquake cycle models to generate predicted surface displacements through the earthquake cycle after it has been spun up to a cycle invariant state (Mallick et al., 2022). We perform a velocity decomposition (Figure 2) to isolate the viscoelastic contribution to the interseismic surface velocities, where v t is the total interseismic velocity (representing what is geodetically measured in between large earthquakes), v b , is the block velocity (here assumed to be purely translational), v sd is the slip deficit velocity (i.e., the elastic only component), and v p is the viscoelastic perturbation velocity.
We let each forward model velocity profile through the earthquake cycle be, v t , and find the best-fit block motion and slip deficit rates. This represents the expected interseismic deformation in the case where there is no time-dependent viscoelastic deformation. The remaining perturbation to the steady-state velocity expectation, v p , can be interpreted as a contribution from exclusively viscoelastic effects (Devries et al., 2017;Hearn et al., 2013;Pollitz and Evans, 2017). It is this component of the total velocity field that we represent with a set of effective displacements at the base of the elastic crust to represent the kinematic effects of the fully coupled mechanical model. For this idealized demonstration, we discretize the basal elastic interface as a set of 20 horizontal 22 km wide dislocation elements spanning the interval from −220 to 220 km ( Figure 4). Anti-plane slip on any one of these buried dislocation patches gives rise to fault parallel surface displacements that reach a localized fault parallel velocity maxima directly above the center of the dislocation patch and decay symmetrically with distance. Variable slip on these dislocation patches can produce a wide range of different surface velocity patterns, including those that very precisely mimic v p . Best-fit basal slip distributions can be determined with linear least squares, and here we consider examples for two different times (immediately after one large earthquake and just before a second large earthquake) through the earthquake cycle are shown in Figure 5. At both times, and for all three rheology cases, v p appears similar for both the forward model and basal dislocation approximation. The best-fit dislocation slip rates are smoothly varying in this noise-free example with the exception of the dislocation patches closest to the fault which are the locus of the sign change.
While these examples demonstrate a simple approach to developing kinematic representations of time-dependent, we can also observe only small distinctions in the distribution of basal slip between different n values (Figure 5) which may limit the extent to which this method might provide constraints on rheology. In the synthetic cases the basal slip distributions generally, but not always, decay with distance from the fault. This decay is more gradual for the n = 1 case than for the non-linear cases, but the n = 3 and n = 6 appear quite similar, consistent with the fact that their surface velocities differ little.

Inference of Viscoelastic Effects Along the North Anatolian Fault
To demonstrate the applicability of this approach to non-synthetic examples, we consider a dense set of interseismic velocity estimates near the NAF (Weiss et al., 2020). The east-west components of these velocities are approximately parallel to the central NAF at a longitude of 34.5° (Figure 6). We extract a set of 9,041 fault parallel velocity components and treat them as observations of total interseismic velocities, v t . While in the synthetic examples above we isolated v p here we simultaneously estimate v b , v sd , and v p using a classical kinematically consistent block model formulation (Matsu'ura et al., 1986) augmented with a set of three basal dislocations. This is a lower resolution parameterization than the 20 elements used for the synthetic examples, but serves effectively as an approximation to a smoothed solution with a larger number of elements. Further, in terms of the behavior of the coupled linear block model system of equations, there is covariance between model parameters. To focus our estimates on the basal dislocation slip rates p * that generate v p we apply geologically estimated, a strongly weighted prior slip rate constraint of 18.5 mm/yr (Kurt et al., 2013). Taken together this linear system to solve here is very small with only five parameters (velocities of two blocks and slip on three buried dislocation patches).
The observed velocities and the velocity field decomposition from the estimated block motion, slip deficit, and basal dislocation contributions are shown in Figure 6. We find that the total velocity reproduces the observed  Numerically calculated total velocities (v t , upper row) and viscoelastic velocity perturbations (v p , bottom row) for three different power-law exponent cases (n = 1 (left-column), 3 (center-column), and 6 (right-column)) with observation times at t = 0, 0.001, 1, 5, 10, 25, 50, 100, 150, and 200 years. Pink to magenta profiles are from early in the earthquake cycle, and blue to cyan profiles are from late in the earthquake cycle. Note that the v p profiles in the lower row use cube-root y-axis scaling to compress the range and maintain sign. The late in the earthquake cycle v p profiles are closer to the steady state expectation, v b + v sd , for the non-linear cases (n = 3 and 6) than they are for the linear case (n = 1).

Figure 4.
Example surface velocity contributions from basal dislocations. In each of the 20 subpanels we show the velocity profiles (solid blue lines) predicted by anti-plane slip on a basal dislocation patches located at a depth of 20 km. Each of these is 20 km wide (gray regions) and they are arranged end-to-end. The velocity profiles extend beyond the edges of each patch due to the effects of the elastic deformation in the finite thickness elastic crust. local maxima in interseismic velocity approximately 80 km south of the NAF fault trace due to the contributions from the basal slip patches where are representing the kinematics. Estimates of fault effective basal slip rates are nearly symmetric (−5.9, 0.7, and 5.2 mm/yr on the basal dislocation patches from south to north) but it's not possible to determine the existence of a potentially similar local maxima on the northern side of the NAF from these data as the region in question is beneath the Black Sea. This example demonstrates the kinematic approach to representing and approximating viscoelastic contributions to total interseismic velocity fields, and may also contribute to and may contribute to explaining the reported discrepancy between geologic slip rates and geodetically inferred slip deficit rates along the NAF (e.g., .  (Arpat & Şaroğlu, 1972(Arpat & Şaroğlu, , 1975 and the rectangle near 35° indicates the region where a subset of fault parallel velocity data are selected from (ranging from X to X′). Lower left panel: Observed velocity data (light gray dots, extracted from the from X to X′ profile), estimated block velocities, v b (solid red line), estimated slip deficit velocities, v sd (dotted red line), estimated viscoelastic perturbation velocities, v p (dashed red line) and v t estimated total model velocities (solid black line). The mean average error for this model is 1.70 mm/yr. Lower right panel: inferred slip rates on three basal dislocation patches spanning the interval from −120 to 120 km.

Conclusions
The kinematic approach to representing viscoelastic contributions to observations of interseismic deformation is designed as a generalizable approach that lies between analytic models with linear rheologies (e.g., Savage and Prescott, 1978) and numerical models with non-linear rheologies (Mallick et al., 2022). Challenges with this approach include: (a) it is not a mechanical approach that allows for the prediction of displacements and stresses in time and (b) while stresses in the upper elastic crust can be calculated, this approach offers no constraints on the stress evolution in the viscoelastic upper mantle. The things that are useful about this approach are: (a) it can effectively capture the effects associated with any viscoelastic rheology even those that are unknown and unconstrained, (b) it is straightforward to implement in two-and three-dimensions, (c) it is computationally lightweight compared to direct numerical simulations. While displacements, strains, and stresses in a linear elastic body are determined by the applied boundary conditions for any shaped body in we have only provided a demonstration of how this might be applied to a two-dimensional anti-plane strain faulting environment. This approach could be used for time-dependent kinematic imaging approaches (Fukuda et al., 2004;Segall & Matthews, 1997) to simultaneously capture the temporal evolution of viscoelastic deformation and fault coupling.