A Secondary Zone of Uplift Due to Megathrust Earthquakes

The 1960 M9.5 Valdivia and 1964 M9.2 Alaska earthquakes caused a decimeters high secondary zone of uplift a few hundred kilometers landward of the trench. We analyze GPS data from the 2010 M8.8 Maule and 2011 M9.0 Tohoku-Oki earthquakes to reveal the persistent existence of a secondary zone of uplift due to great earthquakes at the megathrust interface. This uplift varies in magnitude and location, but consistently occurs at a few hundred kilometers landward from the trench and is likely mainly coseismic in nature. This secondary zone of uplift is systematically predicted by our 2D visco-elasto-plastic seismo-thermo-mechanical numerical simulations, which model both geodynamic and seismic cycle timescales. Through testing hypotheses in both simple and realistic setups, we propose that a superposition of two physical mechanisms could be responsible for this phenomenon. First, a wavelength is introduced through elastic buckling of a visco-elastically layered fore-arc that is horizontally compressed in the interseismic period. The consequent secondary zone of interseismic subsidence is elastically rebound during the earthquake into a secondary zone of relative uplift. Second, absolute and broader uplift is ensured through a mass conservation-driven return flow following accelerated slab penetration due to the megathrust earthquake. The dip and width of the seismogenic zone and resulting (deep) coseismic slip seem to have the largest affect on location and amplitude of the secondary zone of uplift. These results imply that subduction and mantle flow do not occur at constant rates, but are rather modulated by earthquakes. This suggests a link between deep mantle and shallow surface displacements even at time scales of minutes.

earthquakes at the megathrust interface. This uplift varies in magnitude and location, but consistently occurs at a few hundred kilometers landward from the trench and is likely predominantly coseismic in nature. This secondary zone of uplift is systematically predicted by our 2D continuum visco-elasto-plastic seismo-thermomechanical (STM) numerical simulations, which physically-consistently model the dynamics at both geodynamic and seismic cycle timescales. Through testing hypotheses in both simple and realistic setups, we propose that a superposition of two physical mechanisms could be responsible for this phenomenon. First, a wavelength is introduced through elastic buckling of a visco-elastically layered fore-arc that is horizontally compressed in the interseismic period. The consequent secondary zone of subsidence is elastically rebound during the earthquake into a secondary zone of relative uplift. Second, absolute and broader uplift is ensured through a mass conservation-driven return flow following accelerated slab penetration due to a megathrust earthquake. The dip and width of the seismogenic zone and resulting ers et al 2018). In the interseismic period the overriding plate is coupled to the subducting plate along the seismogenic zone. Subduction thus drags the overriding plate landward and down. This causes subsidence from the trench throughout a large part of the seismogenic zone and interseismic compression causes uplift beyond. This uplift slowly tappers to zero in the far field. In the coseismic period surface displacements typically show elastic rebound of these interseismically accumulated displacements (e.g., Reid 1910;Moreno et al 2010). This thus leads to strong uplift as the overriding plate slips seaward, while the coastal regions typically located above the hypocenter manifest subsidence. Again farther land inward, standard (visco-)elastic models models show zero vertical displacements (e.g., Wang 2007;Meltzner et al 2006). This classical conceptual model is contrasted by two great megathrust earthquakes (M w > 8.5) in the 1960s in Chile and Alaska, where a distinct secondary zone of uplift (SZU) was measured landward of the hypocenter (Plafker 1969;Plafker and Savage 1970). These static measurements were, however, made years after the earthquakes. They thus not allowed one to separate contributions from coseismic or early postseismic deformation. The classical interpretation of a very gradual tapering to zero uplift is also contrasted by more recent seismo-thermomechanical (STM) models, which predict the presence of a secondary zone of uplift (van Dinther et al 2013b). These models self-consistently simulate both subduction dynamics and seismogenesis in a setup where visco-elastic structure is governed by conservation laws and a visco-elasto-plastic rheology based on laboratory experiments.
Following the two great megathrust earthquakes in the 1960s several papers identified specific settings and physical mechanisms that would allow for a sec-ondary zone of uplift. For the 1960 M9.5 Valdivia earthquake, Plafker and Savage (1970) reproduced the secondary zone of uplift by introducing a downward curving fault that steepens suddenly. Linde and Silver (1989) reanalyzed the same dataset and suggested that slip must have also occurred until depths up to 65-80 km, while a strong kink in the interface was required below the peak of the secondary bulge to reproduce this feature. Vita-Finzi and Mann (1994) explained the deformation pattern in Valdivia by elastic flexure of a continuous elastic beam following displacements of massand resulting buoyancy effects. For the 1964 M9.2 Alaska earthquake, Plafker (1969) speculated that it could be caused by a sudden increase in horizontal compressional strain and termed it a 'Poisson bulge', while noting this feature as a major unresolved problem. Alternatively, Plafker (1972) shortly postulated a hypothesis of transverse crustal buckling resulting from horizontal compression of the continental plate.
The occurrence of great megathrust earthquakes in the last decade allowed for major advances in understanding horizontal displacements, particularly with respect to the contribution of afterslip, visco-elastic relaxation, and relocking to postseismic deformation (e.g., Wang et al 2012;Sun et al 2014;Klein et al 2016;Govers et al 2018). However, models following the 2010 and 2011 megathrust earthquakes do typically not reproduce a secondary zone of uplift (e.g., Govers et al 2018). Interestingly some models do reproduce (parts of) it, but do not describe it as such (e.g., Miyashita 1987;Sun and Wang 2015).
This literature overview shows that there is no consensus on wether a secondary zone of uplift is a universal physical phenomenon. Additionally, there is no consensus on the physical mechanisms responsible for such a secondary zone of uplift. Through re-analyzing high-quality data from the last decade and dedicated numerical models (Section 2) we aim to understand wether the classical conceptual model of surface displacements should be extended with a secondary zone of uplift.
Our analysis of published data for four great megathrust earthquakes confirms the existence of a secondary zone of uplift (Section 3.1). We then study STM models of different degrees of complexity to propose two physical mechanisms working together to form a secondary zone of uplift (Section 3.2-3.4). Finally, we discuss the limitations, implications and predictions of our findings (Section 4).

Data collection from literature
A secondary zone of uplift in nature can be detected by surveying land elevations before and after an earthquake. Decades ago methods as described in e.g. Plafker (1965) and Plafker and Savage (1970) provided estimates with measurement uncertainties on the order of a few decimeters. Near the coast line relative sea level changes were mapped using local markers, such as high-tide lines or vertical growth limits up to which specific sessile marine organisms or plants can live.
Inland elevation changes were obtained by differencing results from leveling methods obtained in two subsequent surveys. Nowadays, land-based GPS data provide widespread information on vertical displacements with an accuracy on the order of centimeters.
We analyze the megathrust earthquakes for which a decent amount of measurements exists at a few hundred kilometers landward from the trench. This requires a coastline for relative sea level change measurements or land for GPS measurements. Accordingly, we identified four megathrust earthquakes: the 1960 M9.5 Valdivia (Plafker and Savage 1970), 1964M9.2 Alaska (Plafker 1969), 2010 M8.8 Maule (Vigny et al 2011M9.0 Tohoku-Oki (Ozawa et al 2011Sato et al 2011) earthquakes. For these earthquakes, we combined published data into trench-perpendicular profiles of vertical displacements. Relevant aspects regarding the origin of these datasets, specific values for the resulting secondary zone of uplift, and the corresponding tectonic parameters for each subduction zone are summarized in Table 1. The values and uncertainties in the data are adopted from the referenced studies (Plafker and Savage 1970;Plafker 1969;Vigny et al 2011;Ozawa et al 2011;Sato et al 2011). In general, uncertainties decrease as measurements are done latter in time with more accurate acquisition methods. Furthermore, the time interval between the earthquake and survey (∆t) gives a rough estimate of the amount of postseismic deformation that is potentially included in the displacement data (see discussion in section 4.3).

Numerical model
We use the seismo-thermo-mechanical (STM) numerical models developed and detailed in van Dinther et al (2013a), van Dinther et al (2013b) and . These models are based on the continuum-mechanics framework of I2ELVIS, which is a 2-D implicit, conservative finite difference thermo-mechanical code (Gerya and Yuen 2007). The fully staggered Eulerian grid is combined with a Lagrangian markerin-cell technique to allow for large deformation through advecting properties (e.g. lithology, stress) along with the particles they are attached to. The code solves for the pressure as well as horizontal and vertical velocity assuming conservation of mass in an incompressible medium (equation 5) and conservation of momentum Table 1 An overview of the differences between studied earthquakes and data acquisition methods. The time span between earthquake and survey (∆t) provides an indication of the potential amount of postseismic data included, while the uncertainties amongst others depend on the measuring techniques. We include estimates of peak slip, rupture width (defined as downdip width of the zone with slip >5 m), average interface dip (defined here as the average dip for 100 km depth), and downdip limit of the seismogenic zone from the literature. The lowermost block characterizes the surface displacements: HP 1 is the transition from primary uplift to subsidence (first hinge point), whereas HP 2 is the second hinge point (the transition from subsidence to secondary uplift). S 1,max is the maximum subsidence of the primary zone of subsidence and U 2,max denotes maximum uplift of the secondary zone of uplift. Data sources:  The constitutive equations relate strain ratesε ij to deviatoric stresses σ ij using a non-linear visco-elasto-plastic rheology according tȯ This represents a Maxwell visco-elastic body in series with a frictional plastic slider, where η is effective viscosity, G is shear modulus,

Dσ ij
Dt is the objective co-rotational time derivative solved using a time explicit scheme, g p lastic is the plastic flow potential, χ is the plastic multiplier connecting plastic strain rates and stresses, and σ yield is the plastic yield strength. The amount of elastic versus viscous deformation is determined by the viscoelasticity factor (G∆t)/(G∆t+η vp ) (e.g., Moresi et al 2003;Gerya 2010). The non-linear viscosity η is largely follows experimentally determined dislocation creep flow laws as where R is the gas constant (8.314 J/(mol· • C). Stress exponent n, pre-exponential factor A D , activation energy E a and activation volume V a are experimentally determined parameters set for each lithology.
Brittle deformation is modeled using non-associative plasticity with a pore fluid pressure-effective pressure-dependent yield strength Earthquake-like events result from a strongly slip rate-dependent frictional formulation with where µ s is the static friction coefficient, V c is the characteristic velocity, and γ In summary, the resulting code handles both long-term subduction dynamics and short-term seismogenesis in a physically consistent manner. In the large-scale setup this means that the slab and seismogenic zone geometries together with its thermal and viscosity structures evolve autonomously. They influence the corresponding stress and strength distributions, which ultimately lead to the generation of spontaneous earthquake-akin events. Hence we model the interaction between the lithosphere, slab and mantle through spontaneously developing faults and events. An important disadvantage of the current version is the fact that events durations are very long (on the order of a ∼100 years) due to the constant time step of 5 years (van Dinther et al 2013b). This makes inertial dynamics negligible and prevents us from uniquely distinguishing co-from postseismic processes.
However, inertial dynamics in terms of shear wave propagation is resolved in the simplified analogue model setup, although waves are somewhat slow due to low scaled shear wave speeds of gelatin.

Model setups
To better narrow down the physical mechanisms governing a secondary zone of uplift we use two model setups. These setups vary in degree of lithological, rheological and geometrical complexity.   The simplified setup is adapted from van Dinther et al (2013a) and is based on the upscaled analogue modeling setup of Corbi et al (2013). In this setup a rigid, straight slab subducts beneath a (visco-)elastic wedge, which is bounded by a rigid backstop ( Figure 1b). For ease of numerical computation we rotate the setup and gravity by a slab dip of 10 • to align the megathrust interface and slab with the lower boundary. The wedge-shaped fore-arc is confined by a backstop, which is moved further away from the trench to reduce its influence on simulated surface displacements to a minimum (see section 3.3.1). The fore-arc wedge deforms elastically (99.7%). In this study we add lower crustal and upper mantle layers that largely deform viscously (≈98%). The megathrust interface additionally features plastic deformation, as controlled by a seismogenic zone with slip rate weakening friction bounded by slip rate strengthening friction regions.
In each setup a sticky air that deforms viscously at all time steps approximates the free surface (Crameri et al 2011;van Dinther et al 2013b). This allows for unhampered evolution of both temporal and permanent topography.

Results and Analysis
We first compile published vertical displacement data to understand how universal a secondary zone of uplift is (section 3.1). Second we study the universality of a secondary zone of uplift in seismo-thermo-mechanical models with a realistic setup tailored to Southern Chile (section 3.1). This section also describes the evolution and characteristics of region of the SZU throughout the seismic cycle. Section 3.3 analyzes the physical mechanisms responsible for such a SZU through studying both a realistic Southern Chile setup and a simplified wedge model. Finally, we discuss some parameters influencing a SZU (section 3.4).  To understand wether a secondary zone of uplift exists in all great megathrust earthquakes we compile the available vertical surface displacements resulting from great megathrust earthquakes. The collected data for four out of four earthquakes show a secondary zone of uplift (Figures 2 and 3). These secondary zones of uplift are remarkably spatially coherent with 164 out of 167 measurements indicating uplift ( Figure 2). Two measurements that show subsidence are obtained near the second hinge point of the 1960 Valdivia earthquake, which is the location where subsidence changes to uplift (HP 2 ). The location and magnitude of this secondary zone of uplift, however, vary significantly from one tectonic region to another. This can be appreciated quantitatively by studying Table 1), which contains the available data on the secondary zone of uplift and relating earthquake characteristics for each event.

Natural Data
Within the comparable Chilean tectonic region, we observe a correlation in  The M9.2 Alaska earthquake caused a secondary uplift U 2,max of maximum 0.3 m with a contrast between the eastern profile (recorded along the Richardson Highway) and the western profile (Alaska railroad, Figure 3D). The second hinge point at the western transect was measured at ∼500 km, while it occurred at ∼350 km for the eastern transect ( Figure 2D). This difference could arise from the location with respect to the lateral limit of the rupture or from the sharp bend of the slab and trench in Alaska, which influences the slab dip. Such large differences in location of the SZU within one earthquake are not observed for the Chilean and Tohoku earthquakes, where the trenches and slabs are rather straight.  Table 2 in Hashimoto et al 2006). Station SAMP at the eastern side of Sumatra is located at roughly 400 km from the trench and recorded an uplift of 6.2 ± 8.5 mm. Additionally, an uplift of 12.5 ± 7.3 mm was recorded at 600-700 km from the trench in Phuket (Thailand). Interestingly, levelling data following the 1946 M8.2 Nankaido earthquake also shows three locations with uplift in a secondary zone beyond 250 km from the trench (see Figure 10 in Miyashita (1987)). These two and three uplift measurements suggest a secondary zone of uplift could also be present for the 2004 M9.2 Sumatra-Andaman and 1946 M8.2 Nankaido earthquakes. Nonetheless, due the limited statistical meaning of two or three data points in space, we exclude these two earthquakes from our analysis.

Deciphering tectonic control
In an attempt to decipher which tectonic features influence this secondary zone of uplift, we compare values for this admittedly too low number of four earthquakes (Table 1). The amount of secondary uplift U 2,max seems somewhat correlated to earthquake magnitude and thus the total amount of slip. Total slip can for this be approximated as slip times area, which is derived from moment magnitude scaled to seismic moment (Blaser et al 2010) and assumes shear moduli are roughly equal.
More total slip or larger magnitude leads to a higher secondary zone of uplift with the exception of the M9.0 Tohoku earthquake. However, the amount of slip on the shallow portion of the Tohoku megathrust interface was exceptional (e.g., Fujiwara et al 2011). The distance from the second hinge point HP 2 to the trench seems to increase with the downdip width of each rupture and decrease with the dip of the megathrust interface (Table 1). Alaska with the flattest subduction zone shows the most horizontally stretched pattern, whereas the Chilean slab dips most steeply and show a more compressed uplift pattern (Figure 4 and 3). This suggests that if earthquake slip penetrates farther away from the trench, the secondary zone of uplift is shifted accordingly.
In summary, all four megathrust earthquakes studied show a similar displacement pattern including a secondary zone of uplift. Differences in amplitude and hinge point position are considerable and are likely related to slab geometry and rupture size.

Postseismic vertical displacements
These data include different amounts of postseismic deformation mainly due to the different delay times when measuring the displacements (Table 1) Figure 3E). Events are detected based on the slip velocity, i.e. the coseismic phase starts when the markers located just above the interface start to move significantly seaward and end when they return to pre-event levels. The resulting displacements and hinge points in Figure 3E agree well with the reference event of van Dinther et al (2013b, Figure 5 therein) and all show a secondary zone of uplift. The consistency of the location of the SZU for different events is remarkable and seems fairly independent of event details and rupture size. This is even more remarkable, when you consider these do affect the primary zone of uplift.
The universal nature of the SZU is supported by extensive tests in which this realistic setup was more and less drastically changed. In these tests we have taken  care to ensure that the SZU is not influenced by numerical modeling parameters (e.g., domain size and boundary conditions). None of the tested numerical parameters influences the location and magnitude of the SZU in a noteworthy manner (e.g., compare red and blue lines with black lines in Figure 5). Additionally, all models with a wide range of tectonic and material parameters reveal a secondary zone of uplift. These models will be discussed in more detail in section 3.3.2 and 4.2) to better understand the physical mechanisms governing it.  To better understand what happens we analyze the spatiotemporal evolution of one quasi-characteristic event and one seismic cycle in detail ( Figure 6). Panel 6A portrays the vertical surface velocities as a function of distance to the trench (X) and time in time steps (Y). During the interseismic period the fore-arc within 100 km from the trench is dragged landward, since it is coupled to the landward subducting plate (also see spatial snapshots in panels 6B and D). This compression causes uplift of the overriding plate from about 100 to almost 300 km. Interestingly, a secondary zone of interseismic subsidence occurs at distances beyond 300 km from the trench. This model thus predicts very slow subsidence (i.e., less than mm/yr) at a few hundred kilometers landward of the trench.
During the coseismic period accumulated displacements within the overriding plate are largely elastically rebound (see Figure 5 in same models of van Dinther et al (2013b), where on average more than 90% of interseismic displacements is rebound ). This leads to a primary zone of uplift, which propagates seaward along with the rupture that nucleated just above the brittle-ductile transition.
Behind this propagating uplift, the extending overriding plate experiences primary subsidence. Beyond about 250 km the secondary zone of uplift occurs (also see spatial snapshots in panels 6C and E). This uplift occurs over the same time period during which primary coseismic uplift occurs, although it lasts slightly longer. This suggests that the secondary zone of coseismic uplift at least has a distinct coseismic component, although exact distinctions are not allowed due to the low temporal resolution of the numerical model.
Analyzing velocities at depth shows that the secondary zone of uplift is connected to a broad pattern of uplift occurring throughout the whole lithosphere and mantle landward of a line connecting the secondary hinge point at the sur-face and the start of ductile deformation at the megathrust interface (T>450 • C in Figure 6E). The second hinge point at the surface is located roughly above the interplate decoupling point, which is located at the megathrust interface at depths of about 90-95 km. The interplate decoupling point marks the depth below which motions of the hanging wall decouple from those in the footwall, which from then on induce corner flow in the mantle wedge (definition updated from Furukawa (1993) starting with van Dinther et al (2013b)). This is evident from seaward flow in the astenosphere, as opposed to interseismic landward displacement of the lithosphere ( Figure 6D). This is facilitated by significantly decreased overriding plate viscosities due to ambient temperatures approaching 1300 • C. This roughly corresponds to the thermal definition of the lithosphere-astenosphere boundary within the overriding plate. At depth, corner flow is driven by slab penetration continuously, although at variable speeds ( Figure 6D-F). Subduction and slab penetration are slow within the interseismic period. As the megathrust interface decouples during an earthquake, subduction and slab penetration are considerably accelerated. Displacements of the subducting plate are thus not rebound to their original position, but rather subduct faster to catch up with long-term subduction rates.
After the rupture arrests and the seismogenic zone relocks itself, postseismic velocities largely change back toward their interseismic pattern (also see spatial snapshot in panel 6F w.r.t. panel D). This reversal is delayed within 50 km from the trench due to afterslip in the velocity-strengthening domain and beyond about 400 km due to viscous relaxation of the lithosphere-astenosphere. Additionally, the secondary zone of interseismic subsidence starts at 250 km, while it moves progressively downdip until it starts at 300 km in the late interseismic period.

Physical mechanisms governing a secondary zone of uplift
Both observations from nature and results from our numerical models suggest that a secondary zone of uplift is a universal characteristic of great megathrust earthquakes. Here we attempt to identify possible physical mechanisms that are applicable to all subduction zones. To identify and directly test the proposed physical mechanisms we utilize two model setups (section 2.3, Figure 1). We aim to add a SZU to the simple elastic laboratory wedge setup, which simulates the classical model for surface displacements and thus does not include a SZU ( Figure 3F). Excluding an identified potential physical mechanism should also be able to remove the SZU from the realistic, albeit complex setup for Southern Chile ( Figure 3E).
Satisfying both these criteria provides a valid test for potential mechanisms, which discarded several potential mechanisms from acting at all or acting alone.

Elastic rebound following interseismic buckling of a visco-elastically layered lithosphere
We start from the most simple setup, where a rigid slab subducts beneath a wedgeshaped fore-arc (van Dinther et al 2013a) (Figures 1B and 7A). When the fore-arc is homogeneous and virtually elastic, spontaneous cycles of megathrust events confirm results of elastic models with near trench uplift followed by subsidence  Model set 2 introduces a lower crustal layer with reduced viscosity in the original analogue model setup, which introduces a gentle secondary zone of uplift around 250-300 km from the trench (cf. sets 1 and 2 in Figure 7B). This modulation is related to the presence of a thin elastic beam (i.e., the upper crust) separated by a viscously deforming layer (i.e., the lower crust). The thin elastic beam is then free to buckle during the interseismic period in response to horizontal compression due to an end load (e.g., Turcotte and Schubert 2002). Elastic buckling due to horizontal forces typically introduces a secondary zone of subsidence at large distances from the trench (e.g., Figure 6A). This zone of secondary subsidence is approximately rebound in the coseismic period due to reversal of elastic deformation (e.g., Reid 1910), thereby introducing a secondary zone of uplift. Buckling is not observed in the original model with the thick overriding plate (model 1 in Figure 7B), since the horizontal compressional forces are not large enough to buckle a thick beam with a very large elastic thickness.
Model set 3 shows a gentle secondary zone of uplift reduces to a secondary zone of slight relative uplift, albeit absolute subsidence, when the backstop is moved further inland (Figures 7A3 and B). This suggests that, for buckling to be effective in introducing secondary uplift, a means to generate more localized compression is needed. A backstop corresponds to a region that is significantly stronger than the region just trenchward of it. This can result from a transition from sedimentary to magmatic/metamorphic rocks (Byrne et al 1993a) or from thinner, warmer and weaker arc lithosphere to colder, thicker and stronger (e.g., cratonic) lithosphere This is added to reproduce the mechanical structure of the overriding plate in the large-scale model. A weak lower crust and weak mantle wedge create a doublebeam system made of the rigid upper crust (upper beam) and mantle lithosphere (lower beam). This facilitates additional buckling of the lithospheric mantle and thus leads to a higher-order wavelength of surface displacements. Note here that the peak of the secondary bulge is roughly located above the mantle wedge tip, which colocates with the interplate decoupling point. This peak, however, remains below the zero vertical displacement level and the surface at 200-300 km is thus still subsiding.

Need for a second mechanism
Numerical experiments in the analogue model setup show elastic rebound following interseismic horizontal compression of a visco-elastically layered lithosphere is able to introduce a high-order buckling wavelength. However, in several dozens of attempts widely varying mechanical structure through geometries and parameters, we are not able to generate a consistent secondary zone of uplift with amplitudes above the zero level. Our summarized results in Figure 7B2 show that a mechanism to generate broad uplift over larger wavelengths is missing. A rigid backstop in the near field (i.e., within 300 km) might be able to facilitate this for some events, although adding that to an analogue setup with complete visco-elastic layering still only partially ensures elevations above zero. However, the disputable presence of a backstop at those locations in both nature and our large-scale models suggests that a second mechanism is needed.
The need for a second mechanism is confirmed in dozens of experiments in a large-scale setup. A secondary zone of uplift can not be removed by eliminating the thin rigid beam structure that facilitates buckling of the overriding plate.
Increasing viscosities of the lower crust and/or upper mantle on the landward side to 10 25 Pa·s shows that surface displacement patterns are affected to a minor extent, so that a secondary zone of uplift remains present and largely unchanged (green line in Figure 5). This experiment also supports our observations that, at least in the the analogue setup, visco-elastic relaxation within the hanging wall is not a key mechanism. Besides this experiment we ran numerous experiments in the large-scale model aimed to remove the secondary zone of uplift. None of these experiments lead to removing of the secondary zone of uplift, while generally its amplitude and location remained largely unchanged (a few selected experiments are shown in Figure 5). This wide range of experiments rather shows that the additional contributing process an important mechanisms that is a very basic feature of our model, which can not be removed. One option that is difficult to quantify for us is the presence of a smoothly curved interface. However, since we do not observe a secondary zone of uplift in other smoothly curved models (e.g.,

Moreno et al 2009) we suspect curvature is not the key component. If this holds,
then the only other mechanism not present in our analogue model involves the slab and mantle kinematics and dynamics.

Mass conservation following slab penetration
To analyze the impact of subduction dynamics, we return to the spatial cross sections ( Figure 6D-F) and schematically represent our interpretation of the physical mechanism in Figure  supplementary movie S1), which is likely due to catching up of interseismic loading and slab pull. Coseismic slab displacements are thus considerable and on the order of a few tens of meters. For a 150 km wide downgoing region in two dimensions this amounts to a displacement of mantle material of a few square kilometers, which needs to be displaced somewhere to conserve mass (and momentum). These large amounts of mantle are preferentially forced upward beneath the overriding plate ( Figures 6E and 10B). This generates uplift on the landward side of the interplate decoupling point, which is mostly focused within 100-150 km landward of that (i.e., the area of the SZU). These motions also uplift the lithosphere, where at the same time space is being created by the seaward displacements of the lithosphere.
A simple calculation as suggested by Figure 8 indicates that conserving the mass displaced as calculated above can lead to uplift on the order of tens of centimeters.
In addition a smaller portion of mantle is displaced upward beneath the oceanic plate ( Figure 10B). These displacement patterns form two (tilted) convective cells similar to what is observed if mass is conserved due to a sinking object (Figure 1.2b in Gerya (2010)). These two convective cells become more narrow and more clear in the postseismic phase ( Figure 6F). Finally, this mechanism can be better understood through an analogy in which you put your finger into a pot with honey and as a results of this the surface surrounding your finger is uplifted.
In summary, both interseismic buckling of a double-beam overriding plate and upward mass-conservation driven flow in the mantle wedge are both necessary to produce a secondary zone of uplift due to megathrust earthquakes.

Parameter study
Lithosphere buckling and upward flow due to slab penetration suggest that the most influential parameters are related to the geometry of the fore-arc as well as slab and seismogenic zone geometry. This agrees with our preliminary findings based on the limited observations from nature, where the geometries of the slab interface and rupture seemed most important (Section 3.1).
To analyze the role of these geometries in a physically-consistent manner, we varied the slab age. This primarily affects the thermal structure within the slab and thus the megathrust geometry and seismogenic zone width. An older slab has more time to cool and thus leads to lower temperatures at the megathrust interface and a deeper and slightly wider seismogenic zone (especially for the young slabs; see legend Figure 9A). Ruptures in a wider seismogenic zone penetrate more landward and thus cause subsidence and extension further land inward, thereby leading to a landward shift of the spatial surface pattern and the second hinge point ( Figure 9A,B). Figure 9C shows that peak amplitude of the secondary zone of uplift hardly changes with slab age or downdip rupture limit. This is likely due to the fact that the seismogenic zone width only weakly increases with slab age (see legend Figure 9A), such that maximum earthquake size increases only weakly.
This minor increase in seismogenic zone width could potentially also explain why a low correlation between slab age and maximum earthquake magnitude has been observed recently (R=0.05 in Heuret et al 2011), as opposed to original suggestions by Ruff and Kanamori (1980, R=0.627).
We further analyze the role of event magnitude by also including other presented models to allow for a larger variation in earthquake size . We find that the location of the second hinge point is not correlated to event magnitude (Figure 9D). The amplitude of the secondary zone of uplift does appear to increase with increasing earthquake size and slip ( Figure 9E).
Summarizing observations from nature and our models, the geometry of the secondary zone of uplift is mostly influenced by the geometry of the seismogenic zone (i.e., its dip and resulting downdip rupture limit).

Discussion
We first confirm the existence of a secondary zone of uplift both from an observational and numerical modeling perspective. We then hypothesize on the physical mechanisms governing this SZU and propose that two complementary mechanisms

Limitations
The 2D STM models are able to model both long-term subduction dynamics and short-term seismogenesis. This innovation allows them to predict new features and processes, which previously went either not observed or unrecognized. One such example is this consistent presence of a secondary zone of uplift due to megathrust earthquakes. The challenge of bridging all relevant time scales of subduction zone processes has, however, not been fully completed and thus involves some limitations for the large-scale model (section 4.6, van Dinther et al (2013b)). The main limitation is the exceptionally long coseismic duration of our events (years instead of seconds), which does not allow us to seperate out all postseismic and steadystate subduction contributions from the coseismic contribution. The same problem holds for natural data from the two earthquakes in the 1960s, since measurements were done years after the earthquakes. A second, relevant modeling limitation lies within the applied rheologies. These models use the assumption of incompressibility, which is a typical assumption made in geodynamic models. This could decrease the surface response due to accelerated slab penetration, although elimination of this characteristic response is not expected (as supported by similar motions in Sun and Wang (2015) for a Poisson's ratio of 0.25). These limitations have been overcome using adaptive time stepping, rate-and-state friction and compressibility in large-scale models of strike-slip settings . However, bridging from millions of years of subduction to earthquake dynamics remains a challenge (Herrendoerfer 2018). This is thus far only partially accomplished by coupling two different models (van Zelst et al 2019) or for resolving the postseismic phase (Sobolev and Muldashev 2017).

Embedding of proposed mechanisms
The proposed physical mechanisms form a new universal explanation for a secondary zone of uplift and for related modeled displacement patterns. In this section we compare our mechanisms and displacement patterns to other studies and start discussing open questions. Other studies with (visco-)elastic do typically not discuss a secondary zone of uplift (e.g., Wang 2007). Interestingly, also large-scale models of Miyashita (1987) reveal a secondary zone of coseismic uplift, which is also followed by subsidence that propagates landward in the postseismic phase.
In terms of physical mechanisms, the fact that they are universal and apply to all subduction zones is important, because we demonstrated its universal existence for various types of subduction zones. That suggests that a mechanism should not be valid only for one or two subduction zones, but rather for all. Hence the specific geometrical causes to locally uplift material due to abrupt curved steepening of this interface (Plafker and Savage 1970) or distinct slab interface kinks (Linde and Silver 1989) are not deemed relevant. Moreover, a compilation of slab data rather shows smooth interfaces without such rapidly changing slab shapes (Figure The first mechanism of interseismic buckling due to horizontal compression falls in the same category as the transverse crustal buckling of a horizontally compressed continental plate mentioned by Plafker (1972). We, however, emphasize that compression occurred during the interseismic period (not coseismic) and resulted in a secondary zone of subsidence. This was subsequently rebound during the earthquake and hence showed secondary uplift. Second, we added that a realistic visco-elastic layering of the lithosphere is needed to have a thin enough elastic beam that can buckle (Turcotte and Schubert 2002). The need of a realistically layered visco-elastic lithosphere to show higher-order wavelengths in vertical displacements has also been shown by Pollitz (1997). Besides inducing lithospheric buckling through horizontal forces, Vita-Finzi and Mann (1994) modeled buckling and a primary and secondary zone of uplift of an elastic beam atop a viscous astenosphere due to vertical forces (i.e., flexure). These oscillations resulted largely from buoyancy effects following mass displacements within which accelerated subduction lead to overlap and a positive buoyancy force that uplifts the primary zone. In our model density displacements are different and it is rather elastic rebound from an overriding plate that is dragged down with the slab that causes primary uplift. Instead secondary uplift seems to result from a superposition of horizontal buckling and slab-induced return flow.
The second mechanism of accelerated slab penetration and resulting accelerated convection to conserve mass is new. However, the resulting convective displacement patterns throughout the lithosphere-mantle system are similar to those modeled in Sun et al (2014); Sun and Wang (2015) (see Figure 3A in Sun and Wang (2015)). They, however, explain all displacements using the general term of visco-elastic relaxation and appoint an asymmetric rupture inducing greater ten-sion in the upper plate as driving mechanism for these convective displacements with two cells. In our model the rupture is much more symmetric ( Figure 6E) as the slab moves more and thus induces stress changes that are comparable (see Fig-ure3D in van Dinther et al (2013b)). Larger slab displacements are likely caused by the (cyclic) loading of our self-consistent system (and presence of spontaneous low viscosity zones around the slab as in Figure ??A and ?? and off-fault plasticity (van Dinther et al 2013b)). This ensures the relatively compressed slab catches up with slab pull in the coseismic period as subduction was partially inhibited in the interseismic period. Additionally, long-term subduction causes a pre-stress state with large extensional stresses in the slab and near neutral (locally compressional) in the wedge. Nonetheless, we still interpreted that rapid seaward motion of the wedge contributed, since these convective displacements fill up the space created by the seaward displaced wedge. Instead, as a primary control, slab penetration creates accelerated uplift as mass (and momentum) are conserved. Secondary, this uplift is tunneled towards the displaced wedge and creates a secondary zone of uplift. This slab penetration mechanism might also be relevant for visco-elastic relaxation and should be considered and explored in physical explanations of it.
Our choice to not refer to this mechanism as visco-elastic relaxation per se is supported by experiments showing that a secondary zone of uplift when the overriding mantle deforms elastically (green line in Figure 5). Therefore the ability of the visco-elastic mantle to delay and relax overriding plate displacements does not seem critical for a secondary zone of uplift. Instead letting the astenospheric mantle largely deform viscously does not remove its presence either, although it distinctly shifts the secondary zone of uplift in the direction of the land (purple line in Figure 5). This could reflect a more distinct contribution from elastic rebound following interseismic buckling.
This second mechanism could be crudely related to a transient and coseismic version of the long-term thin viscous sheet model combined with corner flow that uplifts the volcanic arc (Wdowinski et al 1989), as remarked upon by Vita-Finzi and Mann (1994). It namely relates to the corner flow of the mantle as induced by subduction to uplift a deformable lithosphere. We observe corner flow during the interseismic period ( Figure 10A). However, the flow pattern of its accelerated version during coseismic slab penetration is distinctly skewed towards the seaside as the overriding plate is displaced in that direction ( Figure 10B). In addition our uplift is rapid and thus involves a larger elastic component.
It is possible that internal slab deformation and slab unbending due to the earthquake (around X=300 km in Figure 6C It is open for discussion wether the second mechanism of mass conservation driven upward flow (Figure 8) occurs at coseismic and/or (early) postseismic time scales. Traditionally one might think that the response from the mantle will be too slow. However, we know that the mantle behaves elastically during and just after the earthquake, because of the prolonged propagation of seismic waves. Moreover, we know that mass must always be conserved, also when the large and heavy slab inevitably penetrates rapidly into the mantle. Together with the apparent occurrence of a mainly coseismic secondary zone of uplift, it seems this could occur at least during -or within the few days following -the earthquake. This might be observed on the Korean peninsula, where uplift is only observed for the five days following the 2011 M9.0 Tohoku earthquake (Kim and Bae 2012).
Nonetheless, the penetration of the slab into the mantle will be at least partially delayed, since on both sides of the slab the mantle resists its penetration as it does for the overriding plate. The resulting shear stresses will need to be relaxed on the time scales mostly dictated by mantle viscosity. In our model these viscosities are low and on the order of 10 18 Pa·s (Figure ??A). This reduction with respect to the surrounding mantle with viscosities of about 10 20 Pa·s occurs due to increased strain rates around the slab, which feedback non-linearly via the stress dependence of dislocation creep viscosity (eq. ??). Consequently, Maxwell relaxation times would be around one or a few years, unless accelerated slab penetration on time scales of minutes can increase them even further in the vicinity of the slab (e.g., to around 10 15−16 Pa·s as in Sobolev and Muldashev 2017). In summary, we estimate mass conservation following accelerated slab penetration operates on both coseismic and (early) postseismic time scales, where it also affects visco-elastic relaxation. How much is coseismic and how much is postseismic might be affected by the tectonic setting.

Implications
These results imply that subduction is not a gradual process of continuous subduction, as typically envisioned within the long-term communities. Subduction rather proceeds in shocks following the brittle stick-slip behaviour of the shallow seismogenic zone. During the interseismic period some subduction can occur. However, locking across a 100 km or 200 km portion of the megathrust interface can partially stall the penetration of the slab. When the whole megathrust unlocks, during a great megathrust earthquake, subduction catches up and the earthquake-induced displacement of the slab induces a significant amount of mantle flow ( Figure 6E, Supplementary movie S1). This makes megathrust earthquakes an integral driver of mantle flow. Similar ideas exploring the interaction between earthquakes and mantle flow are explored in other numerical models, which feature modulation of astenospheric flow (Barbot 2018) and modulation of residual polar wander (Cambiotti et al 2016). This thus suggests a link between deep mantle and shallow surface displacements on timescales of minutes to decades, which is shorter than previously considered.
These numerical results also demonstrate implications for geodetic-based source inversions. The inclusion of visco-elastic layering and their detailed geometrical implementation significantly impacts the resulting coseismic surface displacements (e.g., Figure 7). Conversely, when using these surface displacements to estimate slip at a fault within a homogeneous elastic medium typically used for source inversions, one would artificially adapt fault slip to compensate for the missed partially viscous features. This is similarly observed for fitting interseismic velocity data, where elastic models require the presence of a rigid micro plate (e.g., Simons et al 2007) that is not required using viscoelastic models (Trubienko et al 2013). Additionally, slip artifacts might be introduced by the secondary zone of uplift present in data, but not in an elastic forward model. This might make it difficult to fit model results to the data (e.g., Lin et al 2013, for the Maule earthquake). This supports emerging results that it is important to include a realistic visco-elastic structure in inversion for interseismic crustal deformation and earthquake slip in-

Predicting future observations
Based on our observational and numerical findings, we make several predictions that can be tested as more accurate data becomes available during and prior to future, large megathrust earthquakes.
We predict more secondary zones of uplift will be observed in future great megathrust earthquakes (and maybe also for M>8 or smaller earthquakes). The location of the secondary hinge point, its wavelength, its amplitude, and decay with time will vary with tectonic setting (Figures 3 and 9). The location could move further inland for subduction zones with more shallowly dipping slabs, whose seismogenic zones are wider and earthquakes can thus penetrate more inland.
More deep slip likely also translates into higher amplitudes of secondary uplift.
The contribution from interseismic buckling would be enhanced by the presence of more effective backstops (Figure 7) and by compression of a thinner upper crust and/or mantle lithosphere (as eq. 3-124 in Turcotte and Schubert 2002, predicts lower amplitudes at larger distances ). The contribution of mass-conserving return flow due to slab penetration would be enhanced by increased uplift amplitudes when more slab material displaces more mantle and effectively tunnels it to just landward of the interplate decoupling point (Figures 6E and 8). This is anticipated for subduction zones with slabs (and events) that have a larger lateral and/or depth extent of the slab (as to a small extent occurs for an older and cooler slab), while an event with more slip should be effective as well.
Finally, we predict a secondary zone of interseismic subsidence to occur at similar distances of between 200 and 500 km from the trench ( Figures 6A). This interseismic subsidence will be very slow ( Figure 6D) and is the counterpart or cause that through elastic rebound leads to a secondary zone uplift due to megathrust earthquakes. It likely occurs just on the landward side of the interplate decoupling point, where mantle displacements beneath the overriding plate become dominant. This transition facilitates both interseismic buckling and rapid upward displacements following slab penetration.

Conclusions
We propose to extend the classical earthquake vertical displacement pattern for great megathrust earthquakes from one to two zones of uplifts that flank a zone of primary subsidence. A second, minor zone of uplift was first predicted by physically consistent models starting to bridge long-and short-term dynamics. Subsequently we observed it for all four great megathrust earthquakes studied. This secondary zone of uplift starts at distances between 200 km and 350 km (or 500 km) from the trench and varies in magnitude from 0.4 to 11 decimeter.
Extensive numerical experiments in both realistic and simple setups could not identify a single physical mechanism that is able to respectively remove and add a secondary zone of uplift to the two setups. Instead we hypothesize that a superposition of at least two mechanisms is needed to generate a secondary zone of uplift.
We need a visco-elastically layered fore-arc to form two thin rigid beams that can buckle elastically in response to horizontal compression due to end loading in the interseismic period. This introduces a higher-order wavelength with a secondary zone of very minor subsidence. Elastic rebound due to an earthquake then causes a secondary zone of relative uplift. This is uplifted above zero by displacements that conserve mass (and momentum) following the earthquake-triggered penetration of the slab into the mantle. These upward displacements particularly localize in the about 150 km's landward of the interplate decoupling point, which typically corresponds to the area of the secondary zone of uplift.
We estimate that the most important parameters affecting the secondary zone of uplift are the seismogenic zone dip and (deep) coseismic slip magnitude and limit (or earthquake size). Recent postseismic data and coincident, albeit unresolved, timing in our numerical model point to a for the largest part coseismic nature, although this remains to be confirmed. Predictions from our models in terms of verifiable observations include more secondary zones of uplift (potentially also for smalller earthquakes) and a secondary zone of very minor subsidence. Additionally, we propose a suite of tectonic influences that could start to explain variations in its size and location. In any case a more accurate representation of the viscoelastic structure of the fore-arc helps to understand and invert for inter-, co-and postseismic displacements. Finally, our results imply that subduction is not a gradual processes, but that is rather accelerated and decelerated through seismic cycles following the slab penetration during great megathrust earthquakes. This suggests a link between deep mantle and shallow surface displacements time scales as short as minutes to decades.
Acknowledgements All natural data used for this paper are properly cited and referred to in the reference list. The data from the numerical experiments are available from the authors upon request. We are grateful to Andreas Fichtner and Paul Tackley for additional funding to continue this MSc thesis project. This work was also supported by a computing resources grant from the Swiss National Supercomputing Centre (CSCS; s741). We thank Arnauld Heuret, Emilie Klein, Saulė Simutė and Christophe Vigny for providing us with data as well as Rob Govers for discussions. We also thank three anonymous reviewers, who helped to largely improve the readability of the manuscript, and guest editor Sylvain Barbot for his patience and understanding in accommodating our time restrictions. We acknowledge the use of the geolib package for Matlab (Karney 2013). Finally, we acknowledge our contributions for transparency. L.P. performed and analyzed the computations and wrote the first version of the manuscript with support from Y.D. during and after his MSc thesis. Y.D. designed the study, closely supervised the work, and rewrote the manuscript during major revisions. T.G.
supervised the project. All authors discussed the results.

A Conservation equations
To obtain horizontal velocity vx, vertical velocity vz, and pressure P we solve the conservation of mass and momentum as ∂vx ∂x + ∂vz ∂z = 0 (5) Here σ ij represents the 2D deviatoric stress tensor. The conservation of momentum includes gravitational acceleration g and the inertial term, represented by density ρ times the Lagrangian time derivative of the respective velocity components Dv Dt . The momentum equations include the inertial term to stabilise high coseismic slip rates at low time steps (van Dinther et al 2013a). A time step of five years, however, reduces our formulation to a virtually quasi-static one.
In the large-scale model we also solve the heat equation , where Cp is isobaric heat capacity, DT /Dt is the Lagrangian time derivative of temperature, and qx and qz are the horizontal and vertical heat flux, respectively. The equation includes contributions from conductive heat transport and volumetric internal heat generation H due to adiabatic (de-)compression Ha, shear heating during non-elastic deformation Hs and lithology-specific radioactive heat production Hr (e.g., Gerya andYuen 2003, 2007).