Fast, accurate solutions for curvilinear earthquake faults and anelastic strain
Abstract
Imaging the anelastic deformation within the crust and lithosphere using surface geophysical data remains a significant challenge in part due to the wide range of physical processes operating at different depths and to various levels of localization that they embody. Models of Earth’s elastic properties from seismological imaging combined with geodetic modeling may form the basis of comprehensive rheological models of Earth’s interior. However, representing the structural complexity of faults and shear zones in numerical models of deformation still constitutes a major difficulty. Here, we present numerical techniques for highprecision models of deformation and stress around both curvilinear faults and volumes undergoing anelastic (irreversible) strain in a heterogenous elastic halfspace. To that end, we enhance the software Gamra [37] to model triangular and rectangular fault patches and tetrahedral and cuboidal strain volumes. This affords a means of rapid and accurate calculations of elastostatic Green’s functions for localized (e.g., faulting) and distributed (e.g., viscoelastic) deformation in Earth’s crust and lithosphere. We demonstrate the correctness of the method with analytic tests, and we illustrate its practical performance by solving for coseismic and postseismic deformation following the 2015 Mw 7.8 Gorkha, Nepal earthquake to extremely high precision.
1 Introduction
The deformation of the Earth’s lithosphere at time scales relevant to the earthquake cycle is accommodated by slip on faults and distributed strain in finite regions of anelastic deformation. Faults exhibit complex dynamics governed by nonlinear friction laws [23, 55] that involve a wide range of thermally activated chemical and physical processes [66, 46]. In addition, fault networks can exhibit great geometrical intricacies, and single fault segments are often characterized by morphological gradients [34, 33]. The combination of constitutive and geometrical complexities constitutes a fundamental challenge for earthquake prediction [48].
Crustal and lithosphere dynamics is also characterized by distributed anelastic deformation. For example, in the lower crust or the mantle asthenosphere, it is responsible for loading faults and accommodating transient strains [39]. The kinematics of crustal deformation can be inferred from geodetic or seismic data by discretization of fault surfaces and linear inversions of deformation data during the interseismic, coseismic, and postseismic phases of the earthquake cycle [41, 45, 9, 5]. The development of these techniques in the last few decades has led to an explosion of knowledge on fault behavior [54, 4, 14, 71, 1]. These methods have recently been extended to incorporate the distributed deformation of large domains of the lithosphere [36, 8, 44, 53], such that it is now possible to build models of Earth’s deformation that represent fault slip and distributed strain consistently using elastostatic Green’s functions.
Numerical simulations of fault dynamics using Green’s function and the integral method have produced predictive scenarios of seismicity. Forward models of crustal dynamics are important to reveal the frictional or rheological properties of the Earth, and to make useful predictions about the longterm behavior of the mechanical system. A large body of work focuses on the interaction between localized and distributed deformation [35, 10, 36, 22, 21]. These studies rely on analytic or numerical methods that treat fault slip and viscoelastic flow in fundamentally different ways. The introduction of Green’s functions for distributed deformation allow us to simplify the treatment of rheological heterogeneities, nonlinear anelastic behavior, and strong mechanical feedbacks (for example, between shear heating and rheology). Analytic solutions for rectangular and triangular dislocations [19, 20, 58, 57, 49, 50, 47, 27, 42] and cuboidal strain volumes [36, 8] allow us to build kinematic and dynamic models of lithosphere deformation on curved surfaces with distributed, offfault anelastic strain. Despite being of fundamental importance, all these solutions share an important caveat as they do not represent the lateral variations of elastic properties of the ambient rocks.
The goal of this paper is to describe a reliable numerical method for constructing displacement and stress Green’s functions for curvilinear faults and strain volumes. The method builds upon a threedimensional, adaptive, multigrid elasticity solver for embedded dislocations described in a recent publication [37]. This manuscript presents functional improvements of the method to model localized and distributed deformation in an arbitrary heterogeneous elastic half space. We adopt the approach used in the geodesy community and consider fault surfaces discretized in triangular or rectangular elements. Similarly, we consider volumes discretized in tetrahedral or cuboidal volumes. We provide a short treatment of the modeling approach in Section 2. We demonstrate the relevance of the technique by modeling the coseismic and postseismic deformation on the 2015 Mw 7.8 Gorkha, Nepal earthquake in Section 3. The algorithms described in this paper are implemented in Gamra, a freely available tool for realistic earthquake modeling (see Section 6).
2 Methods
2.1 Elastostatics with embedded dislocations
We assume that Earth’s deformation at static equilibrium is governed by the equations of linear elasticity
(1) 
where is the Cauchy stress and is a forcing term. This assumption is thought to be valid for seismotectonic activity with infinitesimal strain. The stress components are defined using Hooke’s law in terms of the displacement components , Lame’s first parameter , and the shear modulus as
(2) 
We use Einstein summation notation, where each index , , is understood to stand for , , and in turn, repeated indices are summed, and commas (,) denote derivatives. Eqs. 1 and 2 govern the static displacement field induced by fault slip [60, 61], or the quasistatic velocity field driven by viscoelastic flow or other anelastic deformation processes [7, 6]. In the case of quasistatic deformation, the forcing term is per unit time, and is the velocity field. Because a large breadth of physical processes can be captured within the same functional form, solving equations 1 and 2 with realistic material properties can find a broad range of applications.
The method we use to solve these equations is largely described in the previous paper [37]. To summarize, we use a parallel multigrid solver on a staggered, adapted finite difference grid as in Figure 1. For dislocations, we add carefully constructed correction terms that depend on the mesh size. For example, when Eq. 1 is expanded out, it includes the term . Expressing this derivative at point in Figure 1 using standard finite differences gives
(3) 
where is the width of each finite difference cell and the notation denotes the value of at the coordinate . We model dislocations as step functions in displacement . In Figure 1, a dislocation with a discontinuity in passes between the points and . As the resolution goes to zero, the dominant term in the finite difference derivative from Eq. 3 at point is
(4) 
which diverges. This divergence is a numerical artifact. The derivatives are well behaved on both sides of the dislocation, it is only the numerical approximation that is divergent.
To fix this issue, we consider dislocations as a boundary between different domains. To use finite difference expressions like Eq. 3, we have to transport all of the quantities to the same domain as where we want to evaluate the expression. For the example in Eq. 3, the expression is being evaluated at . So the terms and are already on the same side. The term has to be transported by subtracting the slip . This results in a correction term which is equal and opposite to the divergent term in Eq. 4. This is a constant term that does not depend on the solution, so we can fold this correction into a modified forcing term
(5) 
For interpolation stencils for multigrid, we can use the same logic to transport all of the quantities to the same domain. The end result is that the interpolation stencils acquire some extra constant terms.
2.2 Dislocations on triangular patches
Our previous work [37] was applied to dislocations defined on rectangular patches. To extend the method to curvilinear faults, we assume that a fault surface can be represented by a triangular mesh where slip is piecewise uniform. This representation is commonly utilized in geodetic imaging [41, 38, 9, 52] and earthquake cycle modeling [40, 59, 52]. Higherorder discretization methods, such as the nonuniform rational basis spline, are not yet widely adopted in the field of crustal dynamics.
This representation allows us to apply most of the same methods used for rectangular patches. The correction terms of the finite difference scheme for Eq. 5 are calculated using intersections between a triangle and stencil elements as in Figure 1. Figure 2 shows a cross section of a solution for a single triangular fault. Figure 3 shows that, even for a line crossing near a singular corner, the computed stress is well behaved and converges to the analytic solution for a homogeneous halfspace [47]. This result gives us confidence that our numerical approach can capture displacement and stress in the near field of triangular fault elements with high accuracy.
2.3 Strain volumes
We extend the method to strain volumes to represent distributed anelastic deformation, such as viscoelasticity, poroelasticity, and thermoelasticity. Under the infinitesimal strain assumption, the plastic deformation of the lower crust or the asthenosphere can be represented by a linear combination of finitevolume elements. This allows the simulation of complex forward models of deformation with nonlinear constitutive properties using the integral method [36] or directly imaging plastic flow in large domains of Earth’s interior using surface geodetic data [44, 53]. We are interested in modeling the displacement and stress in the halfspace surrounding these elements while incorporating a realistic distribution of elastic moduli. To allow realistic volume meshes, we consider tetrahedra as the basic element, but we also include cuboids for convenience. We assume that these finite volumes undergo anelastic deformation with a piecewise uniform applied strain . A realistic distributed strain can be modeled by adding up many of these individual elements.
The implementation of internal dislocation and anelastic strain differ because for the latter there are no longer step functions in displacement, but rather gradients. Consider Figure 1. The difference between at and at depends on how much of the strain volume is traversed by the stencil. In this case, the difference for is . Using this, we compute a correction for strain volumes in a manner similar to that used to compute Eq. 5, but modified as follows
(6) 
Comparing this with Eq. 5, the only difference is replacing with . This means that after computing for all of the strain volumes, we can reuse all of the machinery for adaptive multigrid solutions that we previously employed for dislocations.
To test our implementation, there is no analytic solution for the deformation of tetrahedral volumes in a half space. However, there is a solution for a cuboid aligned with the surface [36, 8]. So we construct a cuboid out of 6 tetrahedra. Figure 4 shows the computed solution for this arrangement, and Table 1 details the convergence of the error in . We also implement cuboids directly, and the solutions are identical to the solutions of cuboids made up of tetrahedra.
0.1  26.54 
0.05  3.191 
0.025  1.545 
0.0125  1.222 
0.00625  0.4033 
To test more complicated volumes with variable moduli, we create a model using one of the tetrahedra from Figure 4. Figure 5 shows the dimensions of the tetrahedra, and Figure 6 shows a solution for the stress. We do not have an analytic solution for this setup, but we can approximate a tetrahedral strain volume with a number of small triangular faults as in Figure 5. This approach works very well at coarse resolution. At higher resolution, the mesh can see the spaces between the faults, so the details start to differ more significantly. Table 2 demonstrates that the two methods converge as the number of faults increases. This gives us some confidence in the correctness of our implementation of strain volumes with heterogeneous elasticity.
128  0.4129 

256  0.1645 
512  0.0750 
1024  0.0520 
3 The 2015 Mw 7.8 Gorkha, Nepal earthquake
As an illustration of the relevance and performance of the proposed modeling approach, we construct realistic models for the coseismic deformation and initial postseismic relaxation of the 2015 Mw 7.8 Gorkha, Nepal earthquake. This earthquake [3, 26, 25] took place on the Main Himalayan Front, a megathrust that separates the Indian and Eurasian plates [62]. The continental collision is responsible for the highest mountain ranges on Earth and has the potential for very destructive earthquakes [12, 56, 15, 16]. The 2015 Gorkha earthquake only broke a fraction of the megathrust, prompting questions about the dominant role of fault morphology [52, 33] or friction [43] in partial ruptures. The main shock was followed by a detectable transient deformation [72, 28] that was compatible with accelerated fault slip in the downdip extension of the rupture and viscoelastic flow in the lower crust of Southern Tibet [74]. The absence of coseismic and postseismic slip updip of the rupture, i.e., south of Katmandu, has raised serious concerns about the remaining seismic hazard in the region [13]. We simulate the static deformation that accompanies coseismic fault slip on the megathrust and the velocity induced by viscoelastic flow in the lower crust of the downgoing plate. In both cases, the solution encompasses a box km. The boundary conditions on the side and bottom are free slip: zero shear stress and zero normal displacement. The boundary conditions on the top are free surface: zero shear and normal stress. The adapted mesh ranges from a resolution of 9,300 m down to 146 m. This is equivalent to a fixed resolution of . A uniform mesh would require about elements, but the adapted mesh only requires elements. The models took between 16 and 26 hours to solve on a Dell R720 with 16 physical cores (Intel Xeon CPU ES2670).
3.1 Coseismic slip on a curved fault
We use the slip model from [52], which consists of 2,841 triangular patches with individual slip amplitudes and directions. The slip model incorporates a realistic megathrust geometry constrained by structural and geological observations [33] and explains a comprehensive geodetic dataset made of spaceborne radar observations and GPS measurements [52]. The elastic moduli are set with PREM [24]. Figure 7 shows the slip model and computed displacements. Because of the lowangle décollement, the deformation is confined above the hanging wall and above the rupture. Despite the level of detail packed into the model, the simulation required no manual meshing. The regions of highly variable displacement strain are remeshed adaptively [37]. To push the limits of the method, we set the refinement critera to refine when the error in the displacement is greater than 0.1 cm. This is about 0.01% of the maximum displacement, and about 10 times smaller than the error from having the boundaries only 300 km away.
3.2 Distributed postseismic deformation
The stress change induced by the main shock can accelerate viscoelastic flow in ductile regions. Remarkably, the governing partial differential equations for the instantaneous velocity field due to viscoelastic flow are the same as for elastostatic deformation when equivalent body forces (per unit time) are used to represent the irreversible deformation [7, 6]. This equivalence can be exploited to map out the distribution of strain rate in Earth’s interior using the surface velocity field to reveal, for example, the spatial distribution of effective viscosity [44, 53].
We illustrate the capability of computing the instantaneous velocity field due to viscoelastic deformation in the lower crust of the downgoing plate. To generate the strain rate volumes, we assume that the lower crust dives 35 km down a ramp centered on the Main Himalayan Front [17]. We create a mesh of the lower crust using km cuboids. Using the slip model of [52], we use the analytic solution [47] to compute the change in stress at the center of each cuboid. We then set the strain rate of each cuboid assuming an effective viscosity of Pas.
Using these applied strain rate volumes as input, we use Gamra to compute the induced strain rate in the bulk (Figure 8). We set the refinement criteria to refine when the error in the induced velocity is greater than cm/s. This is about 20 times smaller than the error due to the boundaries being only 300 km away. Viscoelastic deformation is more distributed than fault slip owing to the greater depth of the source. The resulting strain rate is more diffuse, requiring more mesh refinement over a larger volume. This results in a model that takes longer to solve (26 hours rather than 16 hours for the coseismic slip model).
4 Conclusions and future work
We have demonstrated that we can robustly and accurately model arbitrary distributions of slip on a curvilinear fault and strain in tetrahedral and cuboidal volumes. A particular advantage of the adaptive mesh is that computing the velocity field due to any single strain volume is faster than for the overall lower crust model, because the mesh can remain coarse far away from the source. This makes the approach well suited for the calculation of elastostatic Green’s functions for localized (e.g., faulting) and distributed (e.g., viscoelastic) deformation.
While the approach is sufficient for a large class of problems, there are still some significant limitations. Topography can play a major influence on the surface displacements when there are large topographic gradients or if the deformation source is shallow (such as during volcanic unrest [18, 73] or underground explosions). In addition, the effect of the Earth’s curvature can play an important role for longwavelength deformation [51]. Also, coupling of deformation with the gravity field [51] is detectable for particularly large earthquakes [30, 29, 67]. These areas will be the focus of future work.
5 Acknowledgements
This research was supported by the National Research Foundation of Singapore under the NRF Fellowship scheme (National Research Fellow Award No. NRFNRFF201304) and by the Earth Observatory of Singapore and the National Research Foundation and the Singapore Ministry of Education under the Research Centres of Excellence initiative. This is EOS publication 184.
6 Computer Code Availability
The algorithms described in this paper are implemented in Gamra [37], which is available at https://bitbucket.org/wlandry/gamra. Gamra is written in C++, uses MPI for parallelism, and depends on a number of packages: HDF5 [64], SAMRAI [31, 32, 65], FTensor [68, 70], libokada [69], muparser [11], and Boost [2]. Gamra runs on on everything from laptops to supercomputers. While Gamra is still under active development, the version associated with this paper has the Mercurial [63] changeset ID

b2d412042fb39cf780ccce431bfb4f476ac74bd7
References
 [1] E. Araki, D. M. Saffer, A. J. Kopf, L. M. Wallace, T. Kimura, Y. Machida, S. Ide, E. Davis, I. Expedition, et al. Recurring and triggered slowslip events near the trench at the nankai trough subduction megathrust. Science, 356(6343):1157–1160, 2017.
 [2] B. Authors. Libokada: A library for analytic solutions of fault patches and strain volumes. http://www.boost.org/, 20032018.
 [3] J.P. Avouac, L. Meng, S. Wei, T. Wang, and J.P. Ampuero. Lower edge of locked main himalayan thrust unzipped by the 2015 gorkha earthquake. Nature Geoscience, 8(9):708–711, 2015.
 [4] W. H. Bakun, B. Aagaard, B. Dost, W. L. Ellsworth, J. L. Hardebeck, R. A. Harris, C. Ji, M. J. S. Johnston, J. Langbein, J. J. Lienkaemper, A. J. Michael, J. R. Murray, R. M. Nadeau, P. A. Reasenberg, M. S. Reichle, E. A. Roeloffs, A. Shakal, R. W. Simpson, and F. Waldhauser. Implications for prediction and hazard assessment from the 2004 Parkfield earthquake. Nature, 437(7061):969–974, Oct. 2005.
 [5] S. Barbot, P. Agram, and M. De Michele. Change of Apparent Segmentation of the San Andreas Fault Around Parkfield from Space Geodetic Observations Across Multiple Periods. J. Geophys. Res., 118(12):6311–6327, 2013.
 [6] S. Barbot and Y. Fialko. A unified continuum representation of postseismic relaxation mechanisms: semianalytic models of afterslip, poroelastic rebound and viscoelastic flow. Geophys. J. Int., 182(3):1124–1140, 2010.
 [7] S. Barbot, Y. Fialko, and Y. Bock. Postseismic Deformation due to the Mw 6.0 2004 Parkfield Earthquake: StressDriven Creep on a Fault with Spatially Variable RateandState Friction Parameters. J. Geophys. Res., 114(B07405), 2009.
 [8] S. Barbot, J. D. Moore, and V. Lambert. Displacement and stress associated with distributed anelastic deformation in a halfspace. Bulletin of the Seismological Society of America, 107(2):821–855, 2017.
 [9] N. M. Bartlow, S. Miyazaki, A. M. Bradley, and P. Segall. Spacetime correlation of slip and tremor during the 2009 cascadia slow slip event. Geophysical Research Letters, 38(18), 2011.
 [10] Y. BenZion, J. R. Rice, and R. Dmowska. Interaction of the San Andreas Fault creeping segment with adjacent great rupture zones and earthquake recurrence at Parkfield. J. Geophys. Res., 98(B2):2135–2144, 1993.
 [11] I. Berg. muparser  fast math parser library. http://muparser.beltoforion.de/, 2014.
 [12] R. Bilham. Location and magnitude of the 1833 nepal earthquake and its relation to the rupture zones of contiguous great himalayan earthquakes. Current Science, 69(2):101–128, 1995.
 [13] R. Bilham. Seismology: raising kathmandu. Nature Geoscience, 8(8):582–584, 2015.
 [14] Q. Bletery, A. Sladen, B. Delouis, M. Vallée, J.M. Nocquet, L. Rolland, and J. Jiang. A detailed source model for the mw9. 0 tohokuoki earthquake reconciling geodesy, seismology, and tsunami records. Journal of Geophysical Research: Solid Earth, 119(10):7636–7653, 2014.
 [15] L. Bollinger, S. N. Sapkota, P. Tapponnier, Y. Klinger, M. Rizza, J. Van Der Woerd, D. Tiwari, R. Pandey, A. Bitri, and S. Bes de Berc. Estimating the return times of great himalayan earthquakes in eastern nepal: Evidence from the patu and bardibas strands of the main frontal thrust. Journal of Geophysical Research: Solid Earth, 119(9):7123–7163, 2014.
 [16] L. Bollinger, P. Tapponnier, S. Sapkota, and Y. Klinger. Slip deficit in central nepal: Omen for a repeat of the 1344 ad earthquake? Earth, Planets and Space, 68(1):12, 2016.
 [17] R. Cattin and J. Avouac. Modeling mountain building and the seismic cycle in the himalaya of nepal. Journal of Geophysical Research: Solid Earth, 105(B6):13389–13407, 2000.
 [18] V. Cayol and F. H. Cornet. Effects of topography on the interpretation of the deformation field of prominent volcanoes  application to etna. Geophysical Research Letters, 25(11):1979–1982, 1998.
 [19] M. Chinnery. The deformation of the ground around surface faults. Bulletin of the Seismological Society of America, 51(3):355–372, 1961.
 [20] M. Chinnery. The stress changes that accompany strikeslip faulting. Bull. Seism. Soc. Am., 53(5):921–932, 1963.
 [21] P. M. DeVries, P. G. Krastev, J. F. Dolan, and B. J. Meade. Viscoelastic block models of the north anatolian fault: A unified earthquake cycle representation of preand postseismic geodetic observations. Bulletin of the Seismological Society of America, 107(1):403–417, 2017.
 [22] P. M. DeVries and B. J. Meade. Kinematically consistent models of viscoelastic stress evolution. Geophysical Research Letters, 43(9):4205–4214, 2016.
 [23] J. H. Dieterich. Modeling of rock friction 1. experimental results and constitutive equations. J. Geophys. Res., 84(B5):2161–2168, 1979.
 [24] A. M. Dziewoński and D. L. Anderson. Preliminary reference earth model. Phys. Earth Plan. Int., 25(4):297–356, June 1981.
 [25] J. Elliott, R. Jolivet, P. González, J.P. Avouac, J. Hollingsworth, M. Searle, and V. Stevens. Himalayan megathrust geometry and relation to topography revealed by the gorkha earthquake. Nature Geoscience, 9(2):174–180, 2016.
 [26] J. Galetzka, D. Melgar, J. F. Genrich, J. Geng, S. Owen, E. O. Lindsey, X. Xu, Y. Bock, J.P. Avouac, L. B. Adhikari, et al. Slip pulse and resonance of the kathmandu basin during the 2015 gorkha earthquake, nepal. Science, 349(6252):1091–1095, 2015.
 [27] Z. Gimbutas, L. Greengard, M. Barall, and T. E. Tullis. On the calculation of displacement, stress, and strain induced by triangular dislocations. Bull. Seism. Soc. Am., 102(6):2776–2780, 2012.
 [28] A. Gualandi, J.P. Avouac, J. Galetzka, J. F. Genrich, G. Blewitt, L. B. Adhikari, B. P. Koirala, R. Gupta, B. N. Upreti, B. PrattSitaula, et al. Preand postseismic deformation related to the 2015, mw7. 8 gorkha earthquake, nepal. Tectonophysics, 2016.
 [29] S.C. Han, J. Sauber, and F. Pollitz. Broadscale postseismic gravity change following the 2011 tohokuoki earthquake and implication for deformation by viscoelastic relaxation and afterslip. Geophysical research letters, 41(16):5797–5805, 2014.
 [30] S.C. Han, C. K. Shum, M. Bevis, C. Ji, and C.Y. Kuo. Crustal Dilatation Observed by GRACE After the 2004 SumatraAndaman Earthquake. Science, 313:658–662, 2006.
 [31] R. D. Hornung and S. R. Kohn. Managing application complexity in the samrai objectoriented framework. Concurrency and computation: practice and experience, 14:347–368, 2002.
 [32] R. D. Hornung, A. M. Wissink, and S. R. Kohn. Managing complex data and geometry in parallel structured amr applications. Engineering with Computers, 2006.
 [33] J. Hubbard, R. Almeida, A. Foster, S. N. Sapkota, P. Bürgi, and P. Tapponnier. Structural segmentation controlled the 2015 mw 7.8 gorkha earthquake rupture in nepal. Geology, 44(8):639–642, 2016.
 [34] J. Hubbard, S. Barbot, E. M. Hill, and P. Tapponnier. Coseismic slip on shallow décollement megathrusts: implications for seismic and tsunami hazard. EarthScience Reviews, 141:45–55, 2015.
 [35] N. Kato. Seismic cycle on a strikeslip fault with rateand statedependent strength in an elastic layer overlying a viscoelastic halfspace. Earth, planets and space, 54(11):1077–1083, 2002.
 [36] V. Lambert and S. Barbot. Contribution of viscoelastic flow in earthquake cycles within the lithosphereasthenosphere system. Geophys. Res. Lett., 43(19):142–154, 2016.
 [37] W. Landry and S. Barbot. Gamra: Simple meshing for complex earthquakes. Computers and Geosciences, 90:49–63, 2016.
 [38] J. P. Loveless and B. J. Meade. Geodetic imaging of plate motions, slip rates, and partitioning of deformation in japan. J. Geophys. Res., 115(B2), 2010.
 [39] S. Masuti, S. D. Barbot, S.i. Karato, L. Feng, and P. Banerjee. Uppermantle water stratification inferred from observations of the 2012 indian ocean earthquake. Nature, 538(7625):373–377, 2016.
 [40] T. Matsuzawa, H. Hirose, B. Shibazaki, and K. Obara. Modeling short and longterm slow slip events in the seismic cycles of large subduction earthquakes. J. Geophys. Res., 115(B12301), 2010.
 [41] J. J. McGuire and P. Segall. Imaging of aseismic fault slip transients recorded by dense geodetic networks. Geophysical Journal International, 155(3):778–788, 2003.
 [42] B. J. Meade. Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space. Comp. Geosc., 33(8):1064–1075, Aug. 2007.
 [43] S. Michel, J.P. Avouac, N. Lapusta, and J. Jiang. Pulselike partial ruptures and highfrequency radiation at creepinglocked transition during megathrust earthquakes. Geophysical Research Letters, 44(16):8345–8351, 2017.
 [44] J. D. Moore, H. Yu, C.H. Tang, T. Wang, S. Barbot, D. Peng, S. Masuti, J. Dauwels, Y.J. Hsu, V. Lambert, et al. Imaging the distribution of transient viscosity after the 2016 mw 7.1 kumamoto earthquake. Science, 356(6334):163–167, 2017.
 [45] J. Murray and P. Segall. Spatiotemporal evolution of a transient slip event on the san andreas fault near parkfield, california. Journal of Geophysical Research: Solid Earth, 110(B9), 2005.
 [46] M. Nakatani. Conceptual and physical clarification of rate and state friction: Frictional sliding as a thermally activated rheology. J. Geophys. Res., 106(B7):13347–13380, 2001.
 [47] M. Nikkhoo and T. R. Walter. Triangular dislocation: an analytical, artefactfree solution. Geophysical Journal International, 201(2):1117–1139, 2015.
 [48] D. Oglesby. Rupture termination and jump on parallel offset faults. Bulletin of the Seismological Society of America, 98(1):440–447, 2008.
 [49] Y. Okada. Surface deformation due to shear and tensile faults in a halfspace. Bull. Seism. Soc. Am., 75(4):1135–1154, Aug. 1985.
 [50] Y. Okada. Internal deformation due to shear and tensile faults in a halfspace. Bull. Seism. Soc. Am., 82:1018–1040, April 1992.
 [51] F. F. Pollitz. Coseismic deformation from earthquake faulting on a layered spherical earth. Geophysical Journal International, 125(1):1–14, 1996.
 [52] Q. Qiu, E. M. Hill, S. Barbot, J. Hubbard, W. Feng, E. O. Lindsey, L. Feng, K. Dai, S. V. Samsonov, and P. Tapponnier. The mechanism of partial rupture of a locked megathrust: The role of fault morphology. Geology, 44(10):875–878, 2016.
 [53] Q. Qiu, J. D. P. Moore, S. Barbot, L. Feng, and E. M. Hill. Transient rheology of the sumatran mantle wedge revealed by a decade of great earthquakes. Nature Communications, 2018, in press.
 [54] G. Rogers and H. Dragert. Episodic tremor and slip on the cascadia subduction zone: The chatter of silent slip. Science, 300(5627):1942–1943, 2003.
 [55] A. Ruina. Slip instability and state variable friction laws. J. Geophys. Res., 88:10,359–10,370, 1983.
 [56] S. Sapkota, L. Bollinger, Y. Klinger, P. Tapponnier, Y. Gaudemer, and D. Tiwari. Primary surface ruptures of the great himalayan earthquakes in 1934 and 1255. Nature Geoscience, 6(1):71–76, 2013.
 [57] R. Sato and M. Matsu’ura. Strains and tilts on the surface of a semiinfinite medium. J. Phys. Earth, 22(2):213–221, 1974.
 [58] J. C. Savage and L. M. Hastie. Surface deformation associated with dipslip faulting. J. Geophys. Res., 71(20):4897–4904, 1966.
 [59] B. Shibazaki, T. Matsuzawa, and A. Tsutsumi. 3D modeling of the cycle of a great Tohokuoki earthquake, considering frictional behavior at low to high slip velocities. Geophys. Res. Lett., 38(L21305), 2011.
 [60] J. A. Steketee. On Volterra’s Dislocations in a SemiInfinite Elastic Medium. Can. J. Phys., 36:192–205, 1958.
 [61] J. A. Steketee. Some geophysical applications of the elasticity theory of dislocations. Can. J. Phys., 36:1168–1198, 1958.
 [62] P. Tapponnier, G. Peltzer, and R. Armijo. On the mechanics of the collision between india and asia. Geological Society, London, Special Publications, 19(1):113–157, 1986.
 [63] T. M. Team. Mercurial source control management. https://www.mercurialscm.org/, 20052018.
 [64] The HDF Group. Hierarchical data format version 5. http://www.hdfgroup.org/HDF5, 20002010.
 [65] The SAMRAI Team. SAMRAI: Structured Adaptive Mesh Refinement Application Infrastructure. https://computation.llnl.gov/projects/samrai/applications, 19972017.
 [66] G. D. Toro, D. L. Goldsby, and T. E. Tullis. Friction falls towards zero in quartz rock as slip velocity approaches seismic rates. Nature, 427:436–439, 2004.
 [67] M. Vallée, J. P. Ampuero, K. Juhel, P. Bernard, J.P. Montagner, and M. Barsuglia. Observations and modeling of the elastogravity signals preceding direct seismic waves. Science, 358(6367):1164–1168, 2017.
 [68] W. Landry. Implementing a high performance tensor library. Scientific Programming, 11:273–290, 2003.
 [69] W. Landry. Libokada: A library for analytic solutions of fault patches and strain volumes. https://bitbucket.org/wlandry/libokada, 20132017.
 [70] W. Landry and M. Sacristan. The FTensor Library. https://bitbucket.org/wlandry/ftensor, 20012017.
 [71] L. M. Wallace, Y. Kaneko, S. Hreinsdóttir, I. Hamling, Z. Peng, N. Bartlow, E. D’Anastasio, and B. Fry. Largescale dynamic triggering of shallow slow slip enhanced by overlying sedimentary wedge. Nature Geoscience, 10(10):765, 2017.
 [72] K. Wang and Y. Fialko. Observations and modeling of coand postseismic deformation due to the 2015 mw 7.8 gorkha (nepal) earthquake. Journal of Geophysical Research: Solid Earth, 2016.
 [73] C. A. Williams and G. Wadge. An accurate and efficient method for including the effects of topography in threedimensional elastic models of ground deformation with applications to radar interferometry. Journal of Geophysical Research: Solid Earth, 105(B4):8103–8120, 2000.
 [74] B. Zhao, R. Bürgmann, D. Wang, K. Tan, R. Du, and R. Zhang. Dominant controls of downdip afterslip and viscous relaxation on the postseismic displacements following the mw7. 9 gorkha, nepal, earthquake. Journal of Geophysical Research: Solid Earth, 122(10):8376–8401, 2017.