Can Plate Bending Explain the Observed Faster Landward Motion of Lateral Regions of the Subduction Zone After Major Megathrust Earthquakes?

Greater landward velocities were recorded after six megathrust earthquakes in subduction zone regions adjacent to the ruptured portion. Previous explanations invoked either increased slip deficit accumulation or plate bending during postseismic relaxation, with different implications for seismic hazard. We investigate whether bending can be expected to reproduce this observed enhanced landward motion (ELM). We use 3D quasi‐dynamic finite element models with periodic earthquakes. We find that afterslip downdip of the brittle megathrust exclusively produces enhanced trenchward surface motion in the overriding plate. Viscous relaxation produces ELM when a depth limit is imposed on afterslip. This landward motion results primarily from in‐plane elastic bending of the overriding plate due to trenchward viscous flow in the mantle wedge near the rupture. Modeled ELM is, however, incompatible with the observations, which are an order of magnitude greater and last longer after the earthquake. This conclusion does not significantly change when varying mantle viscosity, plate elasticity, maximum afterslip depth, earthquake size, megathrust locking outside of the rupture, or nature and location of relevant model boundaries. The observed ELM consequently appears to reflect faster slip deficit accumulation, implying a greater seismic hazard in lateral segments of the subduction zone.

the rupture, with no other nearby stations. Landward velocities up to ∼4 mm·yr −1 greater than before the event were observed in the 5 years after the Iquique earthquake, at stations ∼300-400 km along-trench on either side of the rupture centroid (Hoffmann et al., 2018;Yuzariyadi & Heki, 2021). Hoffmann et al. (2018) found landward increases, with respect to preseismic values, as high as 10 mm·yr −1 in the second year after the event.
The magnitudes of landward velocity changes after all six earthquakes initially vary over time in the period shortly after the earthquake (Yuzariyadi & Heki, 2021). This transient period largely coincides with the previously inferred duration of substantial postseismic transients (particularly afterslip) and lasts ∼5 years after the Tohoku earthquake and ∼2 years after the other, smaller events. The transient behavior of the trench-perpendicular velocity changes is not consistent across all observations. It includes increases, decreases, and even transitions from trenchward to landward changes within the first 2 years after the Oaxaca (Yuzariyadi & Heki, 2021) and Iquique (Hoffmann et al., 2018) earthquakes. However, in all cases, after this transient period, velocity changes stabilize and remain constant, except for a moderate increase in the following 3 years after the Iquique earthquake (Yuzariyadi & Heki, 2021).
An increase of the landward velocity can signify changes in the magnitude or timing of the next earthquake in the area, for instance indicating an increase in slip deficit accumulation and seismic hazard. Ascertaining the mechanism responsible for the landward velocity changes can thus clarify what changes to seismic hazard should be expected where the changes are observed. One interpretation of the observed increase in landward velocities is that it results from an increase in frictional interplate coupling on the megathrust in the vicinity of these observations (Loveless & Meade, 2016). The hypothesis stems from the kinematic inversion of observed velocities from different time spans into interplate coupling. This interpretation implies that the frictionally coupled area on the interface would increase due to a megathrust event hundreds of km away.
Another explanation for the increased landward velocities is that the subducting slab accelerates over a wide portion of the margin as a result of the megathrust unlocking in the rupture zone (Heki & Mitsui, 2013). The hypothesis is consistent with marine GPS-acoustic (GPS-A) observations showing increased Pacific plate velocities close to the rupture zone following the 2011 Tohoku-Oki earthquake (Tomita et al., 2015). However, slab acceleration due to an altered force balance resulting from the coseismic unlocking of asperities can only occur until the ruptured asperities are relocked. Relocking is sometimes inferred to have occurred within a few months to a year after the 2010 Maule, 2011 Tohoku, and other large megathrust earthquakes (Bedford et al., 2016;Govers et al., 2018;Remy et al., 2016;Sato et al., 2011). In that case, transient slab acceleration cannot explain postseismic ELM observed over several years. Nevertheless, no consensus exists on whether relocking is universally rapid, or even whether rapid relocking is needed to explain observations after the Tohoku earthquake (Peña et al., 2019;Watanabe et al., 2014). Both increased coupling and slab acceleration require additional postseismic changes to the subduction system other than well-established postseismic processes (e.g., asperity relocking, visco-elastic relaxation, afterslip, poroelastic rebound). Melnick et al. (2017) proposed another mechanism that would be intrinsic to deformation after large megathrust events. In their mechanical models of combined coseismic and postseismic deformation, over a 100-year timestep, resulting from a drop in megathrust friction, they saw a pattern of velocity changes in the far-field similar to what was observed at GNSS stations following the Maule earthquake. Melnick et al. (2017) and Loveless (2017) proposed that the elastic bending of the plates, in response to postseismic relaxation, explains the far-field landward increases in landward velocities associated with the Maule earthquake. They also argued that the bending producing the velocity changes could cause temporal clustering of megathrust earthquakes by triggering ruptures of asperities. The 2015 Illapel, and 2016 Chiloé earthquakes, which followed the 2010 Maule earthquake in Chile, were interpreted as an example of such clustering (Loveless, 2017;Melnick et al., 2017). This interpretation implies that landward velocity changes may also be responsible for increased shortening rates between clustered historical megathrust earthquakes (Melnick et al., 2017), evidenced for instance by increased subsidence rates recorded by Sumatran microatolls (Meltzner et al., 2015;Philibosian et al., 2014). However, Melnick et al. (2017) and Loveless (2017) compared the effect of the simulated bending with observations only qualitatively, without analyzing the amplitude or temporal evolution of the velocity changes resulting from relaxation with the observed ones. Furthermore, their investigation of the parameter sensitivities and driving mechanisms of the process is incomplete in that it does not include, for instance, bulk rheological parameters, model domain extent, and preseismic interplate locking pattern.
In this paper, we investigate how far-field enhanced landward motion (ELM) may be produced as part of the megathrust earthquake cycle, assuming no variations in the interplate locking pattern or slab acceleration. More specifically, we study whether the expected acceleration produced by plate bending falls within the observed range and under what conditions this bending driven by postseismic relaxation may occur. As part of this goal, we aim to establish the sensitivity of this bending mechanism to key aspects of the megathrust earthquake cycle, such as the earthquake magnitude, the downdip extent of afterslip, and the rheology and extent of the plates and mantle.
We use numerical models of the earthquake cycle, with physically consistent stresses, strains and megathrust slip, to quantitatively simulate the postseismic deformation field. As ELM is observed at several subduction margins, we build generic seismic cycle models, not tailored toward any specific geographical region or megathrust earthquake. We conclude that in-plane bending of the plate probably does produce some ELM, but by itself cannot explain the global observations of ELM.

Concept
We develop three-dimensional mechanical models of the full earthquake cycle. The model geometry involves a realistic slab profile and is uniform in the trench-parallel direction ( Figure 2). Deformation is driven by velocities imposed at the ends of the plates and by the locked portions of the megathrust (asperities), which we define and which accumulate slip deficit during the interseismic period.
It is customary to model the earthquake cycle using the backslip approach, in which the downgoing plate (slab) is assumed to deform only according to the imposed distribution of interseismic slip accumulation and coseismic slip, which must be determined independently of the model. Additionally, afterslip is either also imposed arbitrarily or assumed to be driven by coseismic stress changes and arbitrary background stresses only. In contrast, in our methodology the imposed plate velocity and the mechanical continuity of the two plates determine the pattern of interseismic slip deficit accumulation beyond the asperities, on the rest of the megathrust and on the shear zone downdip of it. The plate velocity, asperities, and plate continuity also determine the shortening of the slab updip of the asperity and its downdip lengthening. Deformation in the model during the coseismic and postseismic phases, after the asperity is ruptured, thus reflects the recovery of elastic deformation of the slab, including faster downdip motion of its middle portion as it catches up with the top and bottom. It also includes Figure 2. Model setup geometry, subdomains, boundary conditions and dimensions. The colors on the external surfaces indicate the displacement boundary conditions (BCs): light orange-free slip along x and y at the lateral sides; cyan-velocity BCs at the top and bottom of the downgoing plate; dark blue-free slip along z at the landward end. The colors on the top and bottom of the slab distinguish the asperity (red), rest of the brittle megathrust (dark fuchsia), shear zone (bright fuchsia), and interfaces where we impose relative motion at the interplate convergence rate (90 mm·yr −1 ). Model boundaries without colors have no displacement BCs, and therefore use zero-traction BCs.
coseismic and postseismic slip that is fully physically consistent with, and determined by, the background stresses as well coseismic and postseismic stress changes.
As the far-field overriding plate is fixed horizontally, all displacements and velocities, both imposed as boundary conditions and resulting from the models, are expressed with respect to the overriding plate. The megathrust is represented by a discrete fault. We consider the two principal mechanisms of postseismic relaxation, afterslip and viscous relaxation (Bürgmann & Dresen, 2008;Diao et al., 2014;Ozawa et al., 2004Ozawa et al., , 2011. We focus on the postseismic period of repeating earthquake cycles.
The models have a tetrahedral finite element mesh with a variable resolution, with nodes as little as 4 km apart in high-strain areas close to the edges of the megathrust and asperities. The reference model includes 533,755 nodes and 3,114,252 elements and contains 6,000 time steps with a size (Δt) of 1 year, corresponding to 20 earthquake cycles. A visualization of the mesh is shown in Figure S1 in Supporting Information S1. Posterior estimates of the model error (Verfürth, 1994) show that the selected mesh is dense enough to support our conclusion that our results are accurate within a few %.
Following each coseismic phase and each afterslip phase, the system is mechanically re-balanced via multiple model iterations. After model spin-up, earthquake cycles are near-identical. There is a difference in surface displacement of less than a few mm between equivalent stages of one cycle and the preceding or following one, while 27 m of interplate convergence occurs over a cycle. We show results from the 19th to 20th cycle.

Model Domain and Geometry
The model geometry extends for 2,000 km along-trench (in the y direction) and 2,200 km in the trench-perpendicular horizontal (x) direction ( Figure 2). The lateral extent of the model domain is chosen so that regions where ELM is expected are not affected by the model edges. We verified that extending the domain further along-trench changes surface motion only minorly and close to the lateral edges. The trench is located at x = 0 and the oceanward model boundary at x = 212 km. The positive x direction thus points oceanward. The domain has a vertical extent of 338 km, with z positive upward and z = 0 at the top of the overriding plate. The distance between the trench and the landward edge of the model is 1,988 km. We used pilot models to verify that enlarging the domain in the along-trench or landward directions does not alter the surface deformation of the overriding plate. We deal with enlargement in the other directions in Section 3.2.6.
The downgoing plate has a thickness of 80 km, consistent with the seismologically detected depth of the lithosphere-asthenosphere boundary of oceanic plates (e.g., Kawakatsu et al., 2009;Kumar & Kawakatsu, 2011), especially for older lithosphere such as on the margins of the Pacific plate (Liu et al., 2017). The top of the downgoing plate follows a trench-perpendicular cross-section through the Slab2 (Hayes et al., 2018) model geometry for the Japan subduction zone, taken to be representative of a typical subduction zone. The overriding plate is 40 km thick with a flat top surface, except for a taper to the trench (at z = −8 km) over a horizontal distance (along x) of 18 km.

Rheology
The model consists of two elastic plates and two asthenospheric domains with isotropic viscoelastic rheological properties ( Figure 2). The constitutive equations (Govers & Wortel, 2005) are based on compressible elastic deformation and incompressible viscous deformation. Here we use a linear viscosity so that the viscoelastic properties follow a Maxwell model with a characteristic stress relaxation time τ ("Maxwell time"; Appendix A1 in Govers et al. (2018)). Most models have a Young's modulus E = 100 GPa and a shear modulus G = 40 GPa (corresponding with bulk modulus K = 66.7 GPa, compressibility β = 1.5 ⋅ 10 −2 GPa −1 , and Poisson's ratio ν = 0.25). These elastic parameters are chosen to be consistent with seismological observations (Dziewonski, 1984) as well as spatially uniform, for the sake of simplicity in studying model sensitivity to their value. Below we discuss how a PREM elasticity profile (Dziewonski & Anderson, 1981) affects the results.
The mantle wedge and sub-slab asthenosphere in most of our models have a viscosity η = 10 19 Pa s. This value is roughly consistent with viscosities determined from observations of postseismic deformation after the 2011 Tohoku-Oki  and 2010 Maule (Klein et al., 2016) earthquakes. These viscosity and shear modulus values correspond to a Maxwell time τ = η/G of 7.92 yr (e.g., Melosh & Raefsky, 1983;Spence & Turcotte, 1979). In Section 3.2 we investigate the sensitivity of the results to material properties.

Boundary Conditions
We impose horizontal and trench perpendicular velocity boundary conditions on the oceanic side of the subducting plate ( Figure 2). The rest of the vertical oceanic side, bounding the sub-slab asthenosphere, is allowed to move only in the vertical direction, because we do not model long-term convective motions. For the same reason, we allow only vertical motion along the vertical continental backside of the model. Slab parallel velocity boundary conditions are imposed where the slab passes through the model bottom boundary. No displacement boundary conditions are applied along the rest of the basal model boundary, corresponding to a zero-traction boundary condition. We apply free-slip boundary conditions at the lateral sides of the model, that is, we allow no displacement perpendicular to these boundaries. We investigate the effect of altering the boundary condition and of moving the vertical location of this boundary in Section 3.2.6.
Isostatic restoring pressures counteract vertical motions of the free surface of both plates (Govers & Wortel, 1993;Wu, 2004). These pressures have a magnitude proportional to vertical displacement. The constant of proportionality is the gravitational acceleration (9.8 m s −2 ) times the density contrast (3,250 kg m −3 at the top of the overriding plate, 2,200 kg m −3 at the top of the oceanic plate).

The Megathrust
We use the slippery node technique (Melosh & Williams, 1989) to model slip along the megathrust in response to shear tractions that develop in the rest of model. The megathrust is infinitely thin in this formulation, and we impose resistive shear tractions to lock parts of the interface during periods between earthquakes. Herman and Govers (2020) demonstrated that interseismic GPS velocities along the South America subduction margin can be well reproduced using a physical model of fully locked asperities with dimensions of ≈50 km on a megathrust that can slip freely otherwise. Low shear tractions up-and downdip of seismogenic asperities is consistent with stable sliding at low friction (Hardebeck, 2015;Ikari et al., 2011;Lindsey et al., 2021;Scholz, 1998). Between earthquakes we therefore consider portions of the megathrust as either locked (asperities) or unlocked.
The megathrust outside the asperities is continuously unlocked and can thus slip freely between earthquakes. However, the mechanical continuity of the plates adjacent to the fault results in accumulation of slip deficit within a distance of ∼50 km from the asperity . To discourage slip on the uppermost portion of the megathrust (Fujiwara et al., 2011;Kanamori, 1972;Moore & Saffer, 2001;Sladen & Trevisan, 2018), we apply small shear tractions at depths shallower than 15 km. Their direction is opposite to coseismic slip and their amplitude is directly proportional to it, with a spring constant of 200 Pa m −1 .
We use asperities that are circular in map view and that have a diameter of 50 km. In all models, the center of one asperity is located 120 km landward from the trench (x = −120 km) in the middle of the model (y = 0). Some models have additional asperities where landward velocity accelerations may be expected (Section 3.2.5). A model "earthquake" occurs by slip on the megathrust when the central asperity is unlocked, which is imposed to happen every 300 years. Unlocking relaxes all shear tractions on the asperity, and the numerical model finds a solution to the new force balance and stresses using 10 iterations. We assume that nothing resists slip deficit release on the megathrust, thus maximizing earthquake magnitude for a given asperity and interplate convergence rate, consistently with the free slip allowed interseismically on the non-locked megathrust. The moment magnitude of the model earthquake is determined by the total accumulated slip deficit on the megathrust. The asperity relocks immediately at the end of the coseismic phase.

Shear Zone Downdip of the Megathrust
The contact between the mantle wedge and the slab, downdip of the brittle megathrust that releases slip deficit coseismically, hosts slow slip, tremors and low-frequency earthquakes (Behr & Bürgmann, 2021;Lay et al., 2012;Tichelaar & Ruff, 1993). Geodynamic models show that a viscoelastic shear zone develops on geological time scales that facilitates differential motion between the slab and the mantle wedge (van Keken et al., 2002). The maximum depth extent of rapid postseismic relative motion (afterslip) on the slab-wedge interface is incompletely constrained but is commonly taken to extend to ∼80-100 km (Diao et al., 2014;Freed et al., 2017;Hu, Bürgmann, Uchida, et al., 2016;Klein et al., 2016;Sun et al., 2014;Yamagiwa et al., 2015) based on postseismic relaxation observations. We simplify the rheological complexity of the contact zone (Perfettini & Avouac, 2004) by representing it by a thin viscoelastic shear zone with a very low viscosity and with the same elastic properties as the surrounding rocks Muto et al., 2019). During the (instantaneous) coseismic motion on the megathrust, there is no differential motion (slip) on the shear zone. Immediately after the coseismic phase,the asperity relocks and very rapid viscous shear stress relaxation occurs in the shear zone. We refer to such rapid postseismic shearing as afterslip. Afterslip is effectively instantaneous in our models. We compute it by rebalancing forces and stresses, using 10 iterations, immediately following the coseismic phase, during which no differential motion is allowed on the shear zone downdip of the megathrust. Model afterslip is consequently complete before the onset of bulk viscous relaxation in the wedge and sub-slab asthenosphere Muto et al., 2019). The shear zone is represented in the numerical model by an infinitesimally thin interface using slippery nodes . Additional relative motion occurs on the shear zone during postseismic and interseismic periods as a result of viscous relaxation and continued convergence. At depth, beyond the downdip end of the shear zone, the wedge and slab are modeled as fully viscously coupled, in agreement with the classical geodynamic view of the subduction system (e.g., Conder, 2005;Kneller et al., 2005;Leng & Mao, 2015;Long, 2013;van Keken et al., 2002). In the context of our earthquake cycle models we are not interested in the steady-state convective flow ("corner flow") in the wedge that is driven by slab motion. We therefore use a similar approach to Savage (1983) along the deeper slab-wedge interface, as follows. The total flow field is the response to both steady subduction and perturbations due to the earthquake cycle. By imposing a steady differential slip rate on the part of the interface where the slab and wedge are fully coupled, we isolate the viscoelastic response to the earthquake cycle only. Using the split node technique (Melosh & Raefsky, 1981), we impose a differential slip equal to the imposed subduction rate.

Slab-Asthenosphere Boundary
We are also uninterested in modeling the steady, long-term, Couette convective flow due to the fact that the slab and underlying asthenosphere are viscously coupled. We thus isolate the response of the sub-slab asthenosphere to the earthquake cycle. Faulted nodes impose the long term subduction velocity as a relative displacement rate along the base of the downgoing plate.

Surface Motion Due to Postseismic Relaxation
Postseismic relaxation in our models involves bulk viscous relaxation and afterslip. Since afterslip is effectively instantaneous in our models, only bulk viscous relaxation produces changes in surface velocities. We compute these velocity changes as Δ ⃗ −pre = ⃗ − ⃗pre , the difference between postseismic velocities at time t after the earthquake and the velocities pre at the last timestep before the earthquake. The latter velocities are taken to represent the near-steady-state contribution of continued convergence with stable coupling at the asperity. When considering cumulative displacement due to both relaxation mechanisms up to a certain time t after the earthquake (Section 3.1.3), we remove the contribution of continued convergence by subtracting ⋅ ⃗pre .
Before computing the velocity changes and displacement due to postseismic relaxation, we correct the velocities and displacement for the small effect of deformation due to long-term slab bending and unbending under the applied boundary conditions. The correction is computed by subtracting velocities from an identical model Since the model geometry has reflection symmetry about a trench-perpendicular plane through the middle of the model (y = 0), we only plot half of the model (y ≥ 0) when showing surface velocity or displacements.

Model Characteristics
We first present a "reference model," so called as its parameters and features will be the reference point for the sensitivity study of Section 3.2. The reference model (Ref) has uniform elastic moduli with realistic yet generic values, not aimed at approximating any specific locality: Young's modulus E = 100 GPa and shear modulus G = 40 GPa. We use a single, central asperity. This way, we prevent additional asperities and their interseismic, coseismic and postseismic signals from interfering with the postseismic relaxation that we study. In later models (Section 3.2.5) we discuss the effect of additional coupling in the form of other asperities, placed laterally along-trench. The asperity is located on the megathrust between the depths of 19.5 and 30.2 km from the surface (11.5 and 22.2 km from the trench). Its unlocking causes coseismic slip corresponding to a moment magnitude M W of 8.9. Afterslip occurs between 40 km (the lower limit of the megathrust interface) and 100 km depth along the slab-wedge interface. Figure 3a shows the cumulative surface displacement due to afterslip on the shear zone separating the slab from the asthenospheric wedge. The trench-perpendicular component of surface displacement of the overriding plate is entirely trenchward (positive). Its amplitude is highest (∼9 m) between the asperity and the trench and decreases with distance, in both the trench-perpendicular and the trench-parallel directions. . These velocity changes are landward as close as 700 km along-trench from the middle of the asperity. The maximum amplitude of the landward velocity change occurs around 110 km from the trench and 1,054 km from the middle of the asperity (Table 1). The trench-perpendicular gradient in landward velocity changes is small in the offshore, near-trench region ( Figure S2 in Supporting Information S1). The velocity changes are highest immediately after the earthquake and decay with time. For instance, the maximum landward velocity change ( −Δ −pre ) is 0.67 mm·yr −1 at t = 1 yr, 0.62 mm·yr −1 at t = 2 yr, and 0.58 mm·yr −1 at t = 3 yr. Figure 3c shows the temporal evolution of trench-perpendicular displacement of one point on the surface of the overriding plate. This point (x = −170 km, y = 1,060 km) is located at the lowest (most landward) Δ 1 −pre at the coastline, taken to have the same horizontal location as the downdip end of the megathrust. Displacement is measured as 0 at the end of coseismic slip. Afterslip, instantaneous in the model, produces the trenchward (i.e., positive) displacement at time 0. Landward (i.e., negative) displacement then occurs due to viscous relaxation. At this location, the trenchward displacement due to afterslip is greater than the cumulative ELM due to viscous relaxation at any time. In the 5 years after the earthquakes, the cumulative landward displacement due to viscous relaxation is everywhere smaller than the trenchward displacement due to afterslip. We expect the viscosity of the asthenosphere to control the rate at which viscous relaxation occurs and thus the temporal evolution of the resulting landward displacement. We later explore the effect of different viscosities (Section 3.2.3).

Maximum Depth Extent of Afterslip
We evaluate the sensitivity of our model results by varying the maximum depth at which the relative motion between the slab and mantle wedge can deviate from the interplate convergence rate. This restricts afterslip and associated slip deficit accumulation on the deep shear zone. This parameter is the major mechanical constraint on material deformation, for a given rheological structure and megathrust locking pattern. In the cutout, the color scale is clipped at 50 mm to show the displacement in the far-field along-trench region. The cyan contour marks 0 trench-perpendicular displacement, separating landward from oceanward motion. The black barbed line shows the location of the trench. The outline of the asperity is shown in red. The dashed orange lines are 2.5 m contours of slip on the shear zone and megathrust due to afterslip. The approximate location of the coastline, taken to be directly above the downdip limit of the locked asperity, is shown in green. The blue dot marks the point at which displacement is plotted in panel (c). Only half the model is shown because of symmetry about the middle (y = 0). (b) Velocity changes (postseismic minus pre-seismic), 1 year after the earthquake, due to viscous relaxation. The color field shows the amplitude of trench-perpendicular velocity, while the vectors show the direction and magnitude of horizontal velocity. The color scale is clipped at ±5 mm·yr −1 to show landward velocity changes. The cyan contour marks 0 trench-perpendicular velocity. The black barbed line shows the location of the trench. The outline of the asperity is shown in green. The dashed orange lines are 2.5 m contours of coseismic slip on the megathrust. The approximate location of the coastline is shown in green. The purple dot marks the point at which displacement is plotted in panel (c). Only half the model is shown. (c) Temporal evolution of total trench-perpendicular surface displacement (dots) at one point in the model (x = −170 km, y = 1,060 km), minus the contribution of the velocity at the end of the interseismic stage, beginning immediately after the coseismic stage.
First, we restrict afterslip to moderate depths, shallower than 75 km (model Aft75). The maximum landward velocity change 1 year after the earthquake is slightly lower than that produced in the reference model with a maximum afterslip depth of 100 km (Table 1). Landward velocity changes also occur ∼50 km along-trench closer to the middle of the asperity. We then restrict afterslip on the shear zone (downdip of the megathrust and thus deeper than 40 km) to very shallow depths, less than 45 km (model Aft45). The landward displacement due to afterslip is greatly reduced, but so is the maximum landward velocity change due to viscous relaxation (Table 1 and Figure 4 and Figure S3 in Supporting Information S1). Next, we allow afterslip to occur at greater depths, as much as 150 km (model Aft150). Compared to the reference model, the landward velocity changes at time t = 1 y after the earthquake have a near-identical maximum amplitude, occurring next to the trench and at a greater along-trench distance from the middle of the asperity (Table 1). Lastly, we completely remove any restriction on afterslip, allowing the relative velocity of the mantle wedge and slab to vary at any depth in response to postseismic deformation (model AllAft). Removing the restriction on afterslip completely eliminates any landward velocity changes due to viscous relaxation. In our models, not allowing time-variable slip rates in the deep shear zone is necessary for enhanced landward velocities to result from postseismic viscous relaxation. The spatial extent of this restriction determines the specific pattern of velocity changes produced.
To better understand the mechanism responsible for ELM generation in our models, we further investigate the relationship between the restriction of motion and the production of ELM by viscous relaxation. We take the model with no limits on afterslip (AllAft) and we introduce a backstop in the overriding plate. We do this by imposing no trench-perpendicular displacement,  at all depths within the plate, at a horizontal distance of 400 km from the trench. This model (AllAftB1) produces landward surface velocity changes due to postseismic viscous relaxation (Table 1). The far-field portion of the plate has an opposite pattern of trench-perpendicular motion, with landward velocity changes in the central part of the model and lower trenchward velocities farther along-trench. Increasing the horizontal distance from the trench to the free-slip boundary to 700 km (model AllAftB2) decreases the maximum landward velocity change 1 year after the earthquake and increases the minimum along-trench distance from the middle to landward velocity changes at that time.

Earthquake Magnitude
We examine the robustness of our results when the size of the earthquake changes. To this end, we reduce the interplate convergence rate, uniformly lowering the slip deficit accumulated and released over an earthquake cycle without varying its spatial pattern. Halving the convergence rate, and thus the seismic moment M 0 , reduces the moment magnitude M W from 8.92 to 8.71 and halves the displacement due to afterslip and the velocity changes due to viscous relaxation at any time. Similarly, reducing M 0 by an order of magnitude (and M W from to 8.25) also reduces the velocity changes and displacement to a tenth. Therefore, with a given interplate locking pattern, ELM produced by postseismic relaxation scales linearly with seismic moment M 0 . This is unsurprising, given the linear nature of the rheologies used in the model. Given the amplitude of the ELM in the reference model, even an earthquake larger than any ever recorded would produce smaller landward velocity changes than the largest values observed at GNSS stations.

Mantle Viscosity
Mantle viscosity controls the rate of viscous relaxation, which produces enhanced landward velocity changes in our reference model. We alter the viscosity η, and thus the Maxwell relaxation time τ, to investigate its effect on our findings. First, in model LoEta1 we decrease η and τ in both the asthenospheric wedge and sub-slab asthenosphere by a factor of 5 compared to reference values, to 2 ⋅ 10 18 Pa s and ∼1.59 years, respectively. We decrease the timestep size by the same factor of 5 to accurately resolve the displacement. The earthquake size (M W = 8.91) and recurrence interval (T = 300 years) are unaltered. The resulting landward velocity changes are dramatically higher than in the model with reference rheology and earthquake size and a single asperity (Table 1 and Figure 5, Figure S4a in Supporting Information S1). However, the maximum amplitudes of the landward velocity changes are still smaller than observed (Section 1 Yuzariyadi & Heki, 2021). The velocity changes decay faster than with the reference viscosity, with the peak amplitude going from 2.5 mm·yr −1 at t = 1 year to 1.6 mm·yr −1 at t = 2 years. In a related experiment (LoEta2), we decrease the viscosity compared to the reference model to 2 ⋅ 10 18 Pa s in the mantle wedge only, keeping it at 10 19 Pa·s in the sub-slab mantle. The maximum landward velocity change after 1 year is more than 50% higher than in LoEta1 (Table 1 and Figure 5, Figure S4b in Supporting Information S1). However, these velocity changes are still lower than observed after the Tohoku-Oki, Tokachi-Oki and Maule earthquakes (Yuzariyadi & Heki, 2021). Also, the model velocities decay rapidly, having a maximum amplitude of 3.8 mm·yr −1 at t = 1 year and 2.0 mm·yr −1 at t = 2 years. The greater landward velocity changes due to viscous relaxation when the viscosity is lower in the mantle wedge only indicate that they are driven by viscous flow in the wedge itself, while flow in the sub-slab mantle opposes them.
Since the earthquake size and elastic properties have not changed, afterslip and the surface motion it causes, via elastic deformation, are the same as in the reference model. The displacement due to the instantaneous afterslip in the model is entirely trenchward. In reality, afterslip has a finite, relatively short duration (a few years following the Tohoku earthquake, for instance, per Muto et al., 2019;Yamagiwa et al., 2015). We compare the cumulative surface displacement due to bulk viscous relaxation in the 2 years after the earthquake (and thus after the instantaneous afterslip) with that due to the afterslip. The landward motion due to viscous relaxation does exceed the trenchward motion due to afterslip, in the along-trench far-field portions of the overriding plate, but by a very limited amount, only as high as ∼1.0 mm.
Increasing the viscosity of both asthenospheric domains by a factor of 5 to 5 ⋅ 10 19 Pa s (model HiEta1), decreases the maximum landward amplitude of velocity changes 1 year after the earthquake (Table 1 and Figure 5, Figure  S5a in Supporting Information S1). It also decreases the rate of decay with time of the velocity changes. For instance, the maximum landward amplitude after 10 years (0.12 mm·yr −1 ) is only 11.5% lower than after 1 year. Increasing the viscosity only in the mantle wedge has a small effect on the maximum landward velocity change at any time (Table 1 and Figure 5, Figure S5b in Supporting Information S1). However, it varies the spatial pattern of the velocity changes significantly, pushing the peak landward value far from the trench and at the lateral edge of the model (y = 1,500 km). This occurs because the relatively small contribution of sub-slab viscous relaxation to surface velocities on the overriding plate is increased.
We have shown how the viscosity of the mantle wedge controls the amplitude and temporal decay of the landward velocity changes. A low viscosity produces large velocity changes, which can even compensate for the trenchward motion due to afterslip and produce net ELM. However, the velocity changes decay rapidly with time as viscous relaxation proceeds and are much smaller already a few years after the earthquake. Higher viscosities produce long-lasting velocity changes due to viscous relaxation, but their amplitudes are very small. Furthermore, the occurrence of afterslip should lead to consistently landward average velocity changes in the months and years after the earthquake during which deep afterslip is occurring. In contrast, velocity changes have been observed to transition from trenchward to landward only after two earthquakes (Iquique and Oaxaca) and within the first year after the event (Hoffmann et al., 2018;Yuzariyadi & Heki, 2021).

Elastic Moduli and Compliance Contrast
We test the sensitivity of our reference model results to changing the elastic parameters of the overriding plate, where the enhanced landward velocities are observed. The effect on modeled ELM of varying the parameters within the realistic range for Earth materials is limited. Furthermore, tailoring the values and spatial distribution of model parameters realistically for specific settings and scenarios is outside the scope of this study. We thus vary the parameters uniformly, choosing extreme values to highlight their effect on ELM and help us investigate the mechanism that produces it. In model LoErefK, we reduce Young's modulus E by a factor of 5, from 100 to 20 GPa, and the shear modulus G from 40 to 6.9 GPa, without changing the bulk modulus K (66.7 GPa) and thus the compressibility = 1 (1.5 ⋅ 10 −11 Pa −1 ). This increases Poisson's ratio from 0.25 to 0.45, close to its uppermost possible value of 0.5. The resulting landward velocity changes are considerably greater and closer to the asperity than in the reference model (Table 1, Figure S6a in Supporting Information S1).
In a related but different experiment (RefEloK), we keep the reference E, bring ν to 0 (as low as possible while not negative) and halve K from 66.7 to 33.3 GPa. β is then twice as large (3.0 ⋅ 10 −11 Pa −1 instead of 1.5 ⋅ 10 −11 ) and G is 50 GPa. The resulting velocity changes 1 year after the earthquake have a very similar maximum amplitude as the reference model, although with a different pattern (Table 1, Figure S6b in Supporting Information S1). In particular, the maximum landward velocity change is closer to the trench but farther from the asperity. The minimum along-trench distance from the middle to the landward velocity changes is greater than in the reference model. The ELM produced by viscous relaxation, when trench-perpendicular displacement is restricted at a certain distance from the trench, is primarily due to the elastic stiffness G of the overriding plate.
We then introduce a contrast in elastic stiffness between the overriding plate within a few hundred km of the trench and the plate farther inland. This represents the contrast between the hot, intensely deformed, tectonically young arc and backarc region, trenchward of the contrast, and the more stable interior of the overriding plate, landward of the contrast. This contrast produces a steep decrease in trench-perpendicular interseismic velocities with distance from the trench in the first few hundred km adjacent to the coast, at the location of the locked asperity, compatibly with observations (e.g., Chlieh et al., 2008;Loveless & Meade, 2010;Métois et al., 2012;Ruegg et al., 2009;Weiss et al., 2016). We use values of Young's modulus E (150 GPa) and shear modulus G (60 GPa) five times greater at horizontal distances from the trench beyond 700 km than closer to the trench (where they are 30 and 12 GPa, respectively). This is roughly the minimum ratio of the contrast that produces a noticeable break in the trench-perpendicular gradient of interseismic velocities and allows for the use of elastic moduli near the bottom and top of the range of realistic values for consolidated rock materials. The surface velocity changes 1 year after the earthquake, have a maximum amplitude of ∼2.2 mm·yr −1 (Table 1 and Figure 6, Figure S7 in Supporting Information S1). This is considerably more than in the reference model, but still less than the observed landward velocity changes (Yuzariyadi & Heki, 2021, see Section 1), despite the model earthquake having a greater magnitude than all observed events but Tohoku-Oki. The peak landward velocity change at that time is located ∼520 km along-trench from the middle of the asperity, while the shortest distance from the middle to landward velocity changes then is ∼400 km. Afterslip still produces substantial displacement there (several tens of mm), causing the average cumulative velocity changes from both afterslip and viscous relaxation to be entirely landward over any length of time after the earthquake.

Adjacent Megathrust Locking
Our previously presented models have a single locked asperity on the megathrust. The observed lateral velocity changes, however, occur in areas with non-zero preseismic landward velocities and thus inferred interplate locking (Loveless & Meade, 2016;Yuzariyadi & Heki, 2021). Therefore, in the LatAsp model we test the effect of locking the megathrust along most of its along-trench extent. Starting with the reference model, we add two intermediate lateral asperities extending from 150 to 650 km along-trench from the middle and two external lateral asperities extending from 800 to 1,300 km along-trench. All lateral asperities are identical to each other and ellipsoidal in map view. Their trench-perpendicular horizontal width (50 km) and distance from the trench (centered 120 km away) are the same as for the central asperity. All asperities are locked interseismically and need to periodically undergo coseismic phases for the model to have multiple earthquake cycle and thus develop background stresses. We use the same recurrence interval of 300 years for each asperity, and thus for the resulting earthquake supercycle. The first set of additional asperities has a coseismic phase 20 years after the central asperity and the second set after 20 more years. During each coseismic phase, the relevant asperities are unlocked (locking tractions are removed) and slip by as much as needed to release all the slip deficit they accumulated interseismically. At the end of each instantaneous coseismic phase, the relevant asperities are immediately relocked and undergo an instantaneous afterslip phase, during which the deep shear zone slips in response to coseismic clip on the megathrust. Afterward, convergence resumes and postseismic viscous relaxation occurs in the region affected by coseismic stress changes, while the other asperities continue their interseismic period. We look at the landward velocity changes due to viscous relaxation 1 year after the earthquake on the central asperity. The amplitude of velocity changes directly above the most external asperities and trenchward of them is decreased, compared to the reference model, to less than 0.5 mm·yr −1 ( Figure S8 in Supporting Information S1). The maximum landward amplitude is decreased and shifted farther from the middle (Table 1 and Figure 7). The overall area occupied by landward velocity changes is very similar, although it locally stretches closer to the middle of the central asperity. Overall, adding additional locked asperities on the lateral portions of the megathrust modifies the specifics of the ELM produced by postseismic viscous relaxation, without fundamentally altering it.

Boundary Conditions and Location of Boundaries
In our reference model (Ref), we imposed zero tractions on the bottom of the mantle, leaving it free to move. We now test the opposite end-member boundary condition (BC): a fixed (no displacement) BC at the base of the mantle (FixBot), which restricts motion more than is realistic. The resulting  landward velocity changes 1 year after the earthquake ( Δ 1 −pre ) are up to 2.5 times larger than in Ref and ∼400 km along-trench closer to the asperity (Table 1 and Figure 8). However, the peak amplitude (1.5 mm·yr −1 ) is still far smaller (by up to two orders of magnitude) than observed (Section 1; Yuzariyadi & Heki, 2021) and decays with time similarly to model Ref. In a separate test (FBThick), we expand the model domain to a depth of 660 km, covering the whole vertical extent of the upper mantle, while keeping the fixed BC at the base of the mantle. This makes the landward velocity changes Δ 1 −pre once again smaller (with a maximum value of 1.1 mm·yr −1 ) and farther along-trench from the middle of the asperity, although less so than in Ref (Table 1 and Figure 8).
We also test the effect of extending the horizontal, trench-perpendicular extent of the oceanic plate and sub-slab mantle. Extending it by 800 km, to a distance of 1,012 km from the trench (model LongOcean), results in landward velocity changes that are closer along-trench to the asperity and larger in amplitude than in the reference model (Table 1 and Figure 8). However, again, the peak amplitude (1.3 mm·yr −1 ) is still far smaller (by up to two order of magnitude) than observed (Section 1; Yuzariyadi & Heki, 2021). The cumulative moment release by coseismic slip and by afterslip before viscous relaxation are both larger than in the reference model, by 0.8% and 4.8%, respectively.

The Mechanism Producing ELM in Our Models
Our model results show that restricting the maximum depth of afterslip, by applying full viscous coupling between the slab and wedge farther downdip, is needed for ELM to be produced during viscous relaxation. Changing this depth affects the resulting ELM pattern, as does introducing a trench-parallel contrast in overriding plate compliance. These sensitivities suggest that the mechanism producing the ELM relies on restricting trench-perpendicular motion.
To better understand why ELM results from viscous relaxation when afterslip is restricted, we simulate the mechanical response of an elastic plate to trenchward tractions, such as those that occur at the base of the overriding plate as the mantle wedge relaxes viscously. Analytical models show that in-plane bending occurs in a semi-infinite elastic plate in response to an outward horizontal force on the free side of the plate (Landau et al., 1986, chapter 13). In the context of the overriding plate in a subduction zone, the free side would correspond to the trench and the force would result from a traction applied at its base. An important result of this conceptually simple analytical solution is that it produces landward displacement of the trench further in the far-field, but only if displacements are imposed to be zero at some distance from the rupture. Although this result is very interesting, it is of limited direct use to ELM because of simplifications in the model setup. Therefore, to identify the nature of the tractions that drive ELM, we explore a two-dimensional (2D) numerical model with a distributions of tractions and boundary conditions closer to the overriding plate in our 3D seismic cycle models.
The 2D model includes only a plate with a uniform thickness of 40 km and the same rheological parameters as in our reference earthquake cycle model. We ignore vertical motion and variation of horizontal motion with depth by using a plane-stress approximation (Govers & Meijer, 2001). We apply a free-slip boundary condition to the lateral and landward edges, while the trenchward edge is left free. A trenchward traction applied on a square patch at the bottom of the plate represents the trenchward tractions due to viscous relaxation in the mantle wedge in the vicinity of the rupture. In response to the traction and boundary conditions, the plate moves trenchward in the middle, but landward laterally. The trench-perpendicular width of the plate determines the location of the trenchward displacement. This conceptually simple model suggests that the ELM produced by viscous relaxation in the 3D earthquake cycle models is due to the fundamental in-plane elastic response to the trenchward viscous flow that occurs in the mantle wedge below. Figure 9 summarizes our understanding of the deformation mechanism that results in ELM due to postseismic relaxation. Viscous relaxation in the mantle wedge produces trenchward motion and applies a trenchward traction to the base of the overriding plate in the vicinity of the asperity, which is absent farther along trench. The viscous coupling of the slab and mantle wedge resists trench-perpendicular motion at a given distance from the trench (∼335 km) in the wedge and in the overlying overriding plate mechanically attached to it. The combination of trenchward forcing in the middle (along-trench) of the plate and of resistance to trenchward and landward motion all along the trench causes the elastic response of the plate to its in-plane (horizontal) forcings to consist of bending, that is, to include lateral motion in the opposite (landward) direction to the trenchward tractions it is subject to.
The location of the viscous coupling between slab and mantle wedge, and thus of the downdip limit of afterslip, determines the resulting pattern of motion, given a certain rheology, asperity size, slip deficit, and boundary conditions. The stress changes associated with coseismic slip and afterslip, and thus the force balance during viscous relaxation, are affected by applying different boundary conditions at the base of the mantle, representing the opposite end-member rheology of the deeper asthenosphere, and by extending the model domain in the depth or oceanward directions. However, the response of the plate to viscous relaxation is still of in-plane bending and produces fundamentally similar motion, albeit with different amplitudes and wavelengths. However, we find that, without any viscous coupling at depth, that is, without any downdip limit to afterslip, nothing prevents the entire overriding plate from moving trenchward, and no bending nor ELM occurs.
Our sensitivity study shows that the landward velocity changes depend much less on its compressibility (while the shear modulus is kept constant) than on the elastic stiffness of the plate (when the compressibility is kept constant). This suggests that bending of the plate is not controlled by the finite compressibility. In other words, the enhanced landward motion on either side of the rupture area is not a consequence of a tendency to conserve volume in response to the trenchward displacement of the plate. Rather, the landward motion is part of the elastic in-plane bending of the overriding plate.

Plate Bending Due To Postseismic Relaxation
Our model results indicate that viscous relaxation following a megathrust earthquake can, by itself, produce ELM as part of a rotational pattern of velocity changes. This is generally consistent with the model results and interpretations of Melnick et al. (2017). In fact, they obtain a pattern of opposing rotation about a vertical axis, including landward motion in the lateral portions of the subduction zone, during the combined coseismic and postseismic deformation caused by a drop in megathrust friction. They propose elastic bending of both plates as the responsible mechanism. Corbi et al. (2022) also obtain a similar pattern of opposing rotation about a vertical axes in an analog model with an elastic wedge and a frictional megathrust separating it from a rigid subducting plate with a constant velocity. In their model, the rotational deformation includes landward motion above one asperity that occurs as a result of frictional failure, specifically the "slip" phase of stick-slip behavior, which is analog of the combined coseismic and postseismic phase of the Earth's megathrust systems. We focus on postseismic relaxation with the same megathrust locking pattern as preseismically. We do not use variable friction that determines megathrust tractions, but use shear tractions large enough to completely lock the asperities, except during the coseismic phase, when they are completely unlocked. We find that lateral ELM is produced by viscous relaxation because the trenchward flow in the mantle wedge due to the latter generates an elastic response in the overriding plate that produces the former. We characterize this elastic response as consisting primarily of in-plane bending, in agreement with the inferences of Melnick et al. (2017) and Loveless (2017).
The ELM in our models relies on the motion of the overriding plate being restricted at a certain distance from the trench. The distance between the trench and this restriction determines the spatial pattern and amplitude of landward velocity changes in response to a given earthquake. Melnick et al. (2017) and Corbi et al. (2022) applied this restriction in their models in the form of a backstop: they allowed no trench-perpendicular horizontal motion on a vertical model boundary parallel to the trench and located close to it on the landward side (700 km from the trench in the model of Melnick et al. (2017); roughly 3 asperity lengths in the analog scale model of Corbi et al.). In contrast, our models extend for nearly 2000 km landward of the trench and instead rely on the conditions imposed on the slab-wedge interface to restrict the motion of the overriding plate. Specifically, the restriction is defined by the depth at which the shear zone that hosts afterslip ends and full viscous coupling begins (100 km in the reference model). There is no direct evidence of the depth at which this transition occurs, or even if there is such a depth. Afterslip has been inferred to occur deeper than 40 km, but there is no evidence of it taking place beyond 100 km depth at most (Diao et al., 2014;Freed et al., 2017;Hu, Bürgmann, Uchida, et al., 2016;Klein et al., 2016;Sun et al., 2014;Yamagiwa et al., 2015). It is plausible, although not certain, that substantially deeper afterslip is not only undetectable at the surface, but truly absent because of viscous coupling between the mantle wedge and slab, in the absence of a localized shear zone. In this case, postseismic viscous relaxation is expected to produce no ELM.

Incompatibility With Observations
The rate of ELM, in our models that produce it, is much smaller than in observations. The observed ELM generally increases with the magnitude of the associated earthquake, as does the ELM in our model. However, the largest observed landward velocity change, following the Tohoku earthquake (M W 9.1), is more than an order of magnitude greater than in our reference model. This is the case even accounting for the smaller magnitude of the model earthquake (M W 8.9) and for the linear scaling of modeled ELM with seismic moment M 0 . For the Maule earthquake, smaller in magnitude (M W ∼ 8.8) than our reference model (M W ∼ 8.9), the maximum observed landward velocity change (∼9 mm·yr −1 ) is an order of magnitude greater than in our reference model (∼0.7 mm·yr −1 ). For the smaller earthquakes, the scaling indicates that ELM should be as much as two orders of magnitude smaller (for the Oaxaca earthquake, M W 7.4). Instead, the observed ELM following those earthquakes is only an order of magnitude smaller than the maximum observed value for the much larger Tohoku-Oki event (Yuzariyadi & Heki, 2021). Furthermore, the observed along-trench location of the ELM is also closer to the middle of the rupture than in the reference model, especially after the Iquique, Bengkulu and Oaxaca earthquakes. Although we can modify the boundary conditions and location of model boundaries to increase the amplitude of landward velocity changes and decrease its along-trench distance from the ruptured asperity, these changes are far from bridging the gap with observed amplitudes.
Our sensitivity tests indicate that overriding plate rheology and restrictions on afterslip affect the amplitudes and spatial pattern of the ELM occurring during viscous relaxation. In particular, introducing a lateral contrast between a more compliant overriding plate lithosphere (in the arc and backarc) and a less compliant plate interior increases the landward velocity changes. Such a contrast has also been inferred to determine the localization of high gradients in horizontal interseismic velocities in the arc and backarc, observed in multiple subduction zones. It is thus likely that the same compliance contrast responsible for the distribution of interseismic velocities amplifies the ELM produced by viscous relaxation, making it at least partly responsible for the fluctuations in the landward velocity changes observed in the early postseismic transient period.
Decreasing the viscosity in the mantle wedge can also produce large landward velocity changes, even exceeding the trenchward motion due to afterslip early after the earthquake. However, if the velocity changes are large shortly after the earthquake, they also decay rapidly with time. Conversely, increasing the viscosity produces a slower rate of decay of the velocity changes, but also lower amplitudes. Either way, the results are not in agreement with the observations, which show consistently long-lasting landward velocity changes, starting right after the earthquake and stabilizing to values of several mm·yr −1 after the transient period of a few years, while afterslip occurs (Yuzariyadi & Heki, 2021). Different rheologies not used in our models, such as Burgers viscoelasticity, could modulate the decay of velocity changes in different ways. For instance, large landward amplitudes could be achieved in the short term while exhibiting long-term viscosities compatible with the geodynamics of subduction zones. However, such rheologies cannot provide both large amplitudes and slow decay to the velocity changes due to relaxation of the same stress changes. Furthermore, the along-trench vicinity to the rupture of the landward velocity changes observed after the Bengkulu, Tokachi-Oki and Oaxaca earthquakes cannot be reproduced by any of the models in our sensitivity testing.
We find that afterslip produces entirely trenchward motion of the overriding plate in all our models. This is in contrast with the hypothesis that the bending producing landward velocity changes is driven by afterslip, proposed by Loveless (2017). In our models, afterslip is modeled as instantaneous and viscous relaxation happens after it has finished. Our implementation of the two postseismic relaxation processes in our models captures the main features of interseismic and coseismic behavior and allows to easily distinguish the contribution of afterslip and viscous relaxation. At the same time, it avoids the computational demands and expanded parameter space caused by simulating viscous flow in a narrow channel. However, in reality, afterslip has a finite duration and interacts with bulk viscous flow (Agata et al., 2019;Masuti et al., 2016;Muto et al., 2019;Yamagiwa et al., 2015). The degree to which afterslip affects the observed velocity changes depends on its distribution through time, as well as on the observation period and method of computation of the velocity changes from the displacement time series. The lack of a realistic temporal distribution of afterslip and the resulting surface displacement is a limitation of our implementation and precludes a direct comparison with observed displacement time series. Nevertheless, the entirely trenchward motion due to afterslip implies that the observed trench-perpendicular velocity changes, with amplitudes of several mm·yr −1 , cannot be explained by afterslip supplementing the motion due to viscous relaxation. This conclusion should not be affected by the lack of two-way feedback between afterslip and viscous relaxation, as the mechanical interaction between the two postseismic relaxation mechanisms has a small effect on the cumulative amplitude of horizontal displacement and on its spatio-temporal evolution, compared to the two processes not interacting (Agata et al., 2019;Muto et al., 2019).
We find that the modeled velocity changes due to viscous relaxation decay with time as the stresses are relaxed (Figure 3c). The contribution of afterslip, when distributed in time, produces a trenchward signal in trench-perpendicular velocity changes, regardless of the locking pattern on the megathrust. The resulting total velocity change due to both relaxation mechanisms should exhibit highly transient behavior, becoming more landward with time as afterslip decays. It should only reach small values (less than a mm·yr −1 in the reference model) and then decay in time as viscous relaxation continues. A transition from trenchward velocity changes in the first year to landward velocity changes in the second year after the Iquique earthquake is indeed observed by Hoffmann et al. (2018). Yuzariyadi and Heki (2021) observe generally less drastic temporal evolution of the velocity changes for all the six earthquakes they consider, including Iquique. However, they only analyze the temporal evolution of velocity changes at one station per earthquake. They do observe a transition from trenchward to landward velocity change in the first and second years, respectively, after the Oaxaca earthquake, at the Puerto Escondido station (OXPE). These transitions likely reflect substantial deep afterslip occurring only shortly after the earthquake, ceasing after about 1 year. Both Hoffmann et al. (2018) and Yuzariyadi and Heki (2021) agree that the velocity changes remain landward after afterslip is inferred to have ceased. No decay in the amplitudes of the trench-perpendicular velocity changes is observed by Yuzariyadi and Heki (2021) after the transient period. Amplitudes are constant after 2 years, except for a slight decay up to 5 years after the Tohoku earthquake and for a moderate increase up to 5 years after the Iquique earthquake. The two longest sets of time series, after the Tohoku and Tokachi earthquakes, show constant velocity changes in the last 4 years. This lack of decay cannot be explained by postseismic relaxation in our models, regardless of the peak amplitudes of velocity changes produced, and constitutes the greatest obstacle in explaining the observed ELM as caused by plate bending in response to postseismic relaxation.
Overall, we find that the elastic response of the plate to viscous relaxation, proposed by Melnick et al. (2017) and Loveless (2017), can plausibly occur, although only if full viscous coupling between the slab and mantle wedge is assumed to occur at a certain depth. We confirm that this response consists primarily of in-plane bending caused by the trenchward flow in the mantle wedge during viscous relaxation. However, according to our simulations, it is extremely unlikely that the temporal and spatial pattern of observed landward velocity changes later described by Yuzariyadi and Heki (2021) is primarily produced by bending in response to postseismic relaxation.

Seismic Hazard Implications
If the observed velocity changes are not attributable to bending caused by viscous relaxation, they must be caused by other mechanisms. One explanation is that large megathrust earthquakes result in changes in the interplate coupling on the megathrust, specifically an increase in the area of strong frictional coupling on the megathrust in the region where ELM occurs (Loveless & Meade, 2016). An increased area of coupling is a straightforward possible interpretation for any landward change in velocity at subduction zones. However, no explanation has been proposed for a megathrust earthquake rupture causing friction increases hundreds of km away. Another explanation is that the velocity of the slab transiently increases postseismically due to the altered force balance caused by unlocking the megathrust during the earthquake (Heki & Mitsui, 2013). Yuzariyadi and Heki (2021) test the correlations between the ELM they describe and the earthquake features predicted by the transient slab acceleration hypothesis, and find the evidence inconclusive but compatible with the hypothesis.
Both increased coupling and slab acceleration invoke an increased slip deficit, compared to the far-field, steady-state plate convergence rate, under the lateral far-field areas where the ELM is detected. Therefore, regardless of which of the two explanations is correct, it is likely that the seismic hazard increases at the locations and time at which enhanced landward velocities are observed. Yuzariyadi and Heki (2021) observe that the enhanced landward velocities correlate with increases in the seismic rate in the relevant lateral portions of the megathrust, which they interpret as evidence of increased stressing rates. The in-plane plate bending in our models occurs regardless of the presence of any coupling on the lateral portions of the megathrust, and thus also when the rest of the megathrust cannot have shear traction changes. The observed increase in seismicity is probably caused by a the mechanism other than plate bending that is needed to explain the observed landward motion requires. This mechanism is likely related to increased stressing and slip deficit accumulation on the lateral portions of the megathrust.
Further research is needed to investigate frictional behavior of the megathrust interface possibly responsible for increased coupling. Discriminating between the two mechanisms proposed to produce this increased slip deficit accumulation is necessary to distinguish whether the increased hazard consists of a greater likelihood of rupture (implied by greater stressing rate due to slab acceleration) or greater peak slip during the future ruptures (a possible consequence of greater frictional coupling on the megathrust). Future studies should also look for further geodetic evidence of transient slab acceleration, including elsewhere in the megathrust subduction system.

Conclusions
Postseismic viscous relaxation can indeed produce enhanced landward motion (ELM). The mechanism producing ELM in our models is the elastic, in-plane response of the overriding plate to the trenchward viscous flow due to relaxation in the mantle wedge. This elastic response consists largely of in-plane elastic bending of the plate. This mechanism relies on the restriction of afterslip provided by the viscous mechanical coupling of the mantle wedge and slab beyond the maximum depth of afterslip. Megathrust coupling in the lateral portions of the interface, above which ELM is observed, is not needed nor interferes significantly with the production of ELM in the models by postseismic viscous relaxation.
Enhanced landward velocity changes as part of the models plate bending due to viscous relaxation are small compared to observations. They also decay noticeably with time. This behavior is inconsistent with the observations of large ELM (several mm·yr −1 to a couple cm·yr −1 ) persisting over multiple years and only exhibiting transient behavior shortly after the earthquake. Furthermore, the geodetically observed ELM also occurs at smaller along-trench distances from the rupture than produced by plate bending in our models. We conclude that the observed ELM requires mechanisms other than postseismic plate bending. The most plausible explanation is thus that slip deficit accumulates at greater rates at the locations and times at which lateral landward velocity changes are observed, increasing seismic hazard there and then.

Data Availability Statement
The model output files that were used for the figures in the main text are available in the FAIR-compliant Yoda repository of Utrecht University at https://doi.org/10.24416/UU01-D7MWAP. Finite element meshes for the models in this paper are generated using Gmsh (Geuzaine & Remacle, 2009)