A Secondary Zone of Uplift Measured After Megathrust Earthquakes: Caused by Early Downdip Afterslip?

A secondary zone of surface uplift (SZU), located from 200 to 400 km landward of the trench, has been measured after several megathrust earthquakes. The SZU reached a few centimeters hours to days after the 2011 Mw 9.1 Tohoku (Japan) and 2010 Mw 8.8 Maule (Chile) earthquakes. Published coseismic finite‐fault models for these events do not reproduce the measured SZU. One interpretation is that this SZU is universal, driven by volume deformation around the slab interface (van Dinther et al., 2019, https://doi.org/10.1007/s00024-019-02250-z). In contrast, we demonstrate the SZU may instead result from slip on the slab interface, and suggest it might be caused by rapid afterslip. We can reproduce the SZU with fault slip if elastic heterogeneities associated with the subducting slab are accounted for, as opposed to assuming homogeneous or layered elastic lithospheric structures.

. Synthetic and observed trench perpendicular profiles of vertical surface displacements. (a) Vertical surface displacement induced by a ∼40-km-deep primary slip patch, by a secondary downdip patch (∼90-km-depth), and the sum of the two. The zoomed inset (c) shows that the sum of these two patches induces a ∼10 cm secondary zone of uplift (SZU). (b) Cross section of the synthetic subduction zone, with the location of the primary and downdip slip patches. (d) Coseismic static vertical displacement measured by survey and continuous GNSS for the 2010 M w 8.8 Maule earthquake for profile A (three near-trench points are from profile B to mimic (a)), and predictions from several published models. Location of the profiles are in Figure S1 in Supporting Information S1 and Figure 5, data from Vigny et al. (2011). (e) The zoomed inset shows the inability of published finite fault slip models to explain the measured SZU. Predictions from Delouis et al. (2010), Luttrell et al. (2011), Pollitz et al. (2011), and Lin et al. (2013) have been produced using models from the SRCMOD database (Mai & Thingbaijam, 2014); others have been reproduced from published material (Langer et al., 2020;Lorito et al., 2011;Moreno et al., 2012;Yue et al., 2014). Note that these models were derived using different datasets (sometimes including only a subset of the data shown here). Vertical bars indicate measurement errors, typically ∼10-20 mm and therefore smaller than the size of the dot. expect that a finite-fault model could infer a downdip slip patch to explain any observed SZU. However, existing published slip models do not.
In the following, we investigate under which assumptions the SZU can, or cannot, be predicted with fault slip. We begin by assuming that the SZU is coseismic by default. We explore the effect of assuming homogeneous elastic properties versus a stiffer subducting slab and more compliant forearc, on predicted surface displacements. We first investigate the effect of 3D elastic heterogeneities for a synthetic subduction case. We then focus on the Maule event, for which the SZU reaches a few cm and could not be reproduced ( Figure 1c) even with added complexity in crustal properties: curved and deeper slab geometries, topography, heterogeneous crustal elastic properties, etc. (Langer et al., 2020;Lin et al., 2013;Moreno et al., 2012). While we do not discard the possibility that the SZU might be affected or result from deformation of the volume around the slab interface, we show it may be caused by slip on this interface. We end with a discussion of the geodetic datasets that have recorded the SZU for the 2010 M w 8.8 Maule and 2011 M w 9.0 Tohoku earthquakes, and discuss the timing of the SZU relative to the mainshocks.

A Synthetic Example: Secondary Zone of Uplift Caused by Downdip Slip
We begin by designing a synthetic subduction zone, where the lithosphere is divided in domains of different elastic properties, generic trench-perpendicular topographic variations and a curved slab interface whose architecture varies slightly along strike ( Figure 2f). This subduction zone is characterized by a stiff plunging slab overlain by a compliant oceanic crust; the continental domain consists of a 35-km-thick crust, more compliant than the underlying mantle whose density increases with depth (domain properties detailed in Text S2, Table S1, Figures S2 and S3 in Supporting Information S1). We apply slip on a limited region of the slab interface ( Figure 1b). Because of the inhomogeneous elastic structure, we rely on a finite element approach (Pylith, Aagaard et al., 2013) to calculate surface displacements.
We first compare the strain produced by a ∼40-km-deep slip patch on the assumed fault, embedded either in a model with 3D variations of elastic properties or in a layered crust ( Figure 2). The layered crust is derived from the continental domain of its 3D counterpart and does not incorporate variations in topography ( Figure 2g). Relative to the layered elastic models, the 3D-heterogeneous models produce a primary zone of subsidence, located 150-200 km from the trench, that is smaller in amplitude and tapers to zero closer to the trench. In the region of primary subsidence, the impact of elastic heterogeneity is ∼5 times larger for vertical displacements than for horizontal ones (Figure 2, Figure S4 in Supporting Information S1, 25% of peak amplitude vs. 5%, respectively).
We then assume two slip patches, the primary patch peaks at 17 m of slip while the secondary downdip patch has 3.5 m of slip ( Figure 1b). We here consider all slip to be coseismic. With the heterogeneous elastic model, we calculate the induced displacement offsets at 50 locations randomly distributed at the surface, to which we add two E-W profiles. The profiles mimic the spatial distribution of the GNSS data of the Maule event ( Figure  S1 in Supporting Information S1, Vigny et al., 2011). Induced displacements are similar to the ∼15-cm-uplift measured 250-300 km away from the trench after the Maule earthquake (Figures 1a and 1d, Figure S1 in Supporting Information S1). We add white and spatially correlated noise to these synthetic data, and try to recover the target slip patches assuming the correct fault geometry with subfaults larger than those used to produce synthetic data. We assume an elastic structure that is different from the one used to calculate synthetic surface displacements. The assumed structure is either layered (Figure 2g), or with 3D variations (Figure 3d). We use a Bayesian sampling approach to infer fault slip from the synthetic displacement (detailed in Section S1 in Supporting Information S1, Minson et al., 2013).
When the crust is assumed layered (or homogeneous), the secondary uplift cannot be fit (and is not within posterior uncertainty, Figures 3a and 3c, Figure S6 in Supporting Information S1, respectively). Relative to the model with heterogeneous elastic properties, a layered crust produces wider and larger primary zone of subsidence, while the horizontal displacements are only slightly impacted ( Figure 2). The amount of slip required to explain the horizontal displacements is incompatible with the slip required to explain the vertical ones. Most inversions will typically favor fitting the horizontal measurements, since they are larger and usually have smaller measurement uncertainty.
With a layered or homogeneous crust, limited downdip slip is imaged, as required by the horizontal displacements, if the fault is deep enough ( Figure S6 in Supporting Information S1). Assuming a fault model that is too shallow, and/or subject to unphysical spatial smoothing, can prevent resolution of the downdip patch ( Figure S7 in Supporting Information S1). The SZU can be produced with incorrect inferred slip, and to the detriment of the fit to the horizontal displacements, if assuming very low measurement errors for the vertical displacements only (1 mm, i.e., very strongly favoring their fit) and a fault geometry that extends to great depths ( Figure S8 in Supporting Information S1).
In contrast, adopting a relatively realistic lithospheric structure (e.g., with 3D heterogeneities in elastic properties for a typical subduction zone, even if the properties are imperfectly known, detailed in Table S2 in Supporting Information S1), allows one to reproduce the SZU, and to recover the downdip slip patch (Figures 3b and 3c). Accounting for uncertainties in elastic properties (Figures 3c and 3d) improves the fit to the data. With this simple synthetic example, we show that a SZU can be produced by downdip slip on the slab interface by accounting for 3D variations in elastic properties.

Recovering the Secondary Uplift of the 2010 M w 8.8 Maule Earthquake
We build a realistic elastic structure across the subduction zone for the calculation of the Green's functions (Figures S10 and S11 in Supporting Information S1, slab geometry from Slab2, elastic properties from LITHO1.0, topography from ETOPO1, Amante & Eakins, 2009;Hayes et al., 2018;Pasyanos et al., 2014). While more detailed velocity models and datasets are available, our goal is to explore the secondary uplift, not to image the slip in detail. We also account for potential uncertainties in the assumed fault geometry and elastic properties (following the methodology presented in Ragon and Simons, 2021). Uncertainties in fault geometry are calculated by varying the dip of the assumed slab geometry while keeping the location of the trench and elastic properties fixed. Note that changing the fault geometry to fit the SZU has already been attempted by several authors (Langer et al., 2020;Lin et al., 2013), without success, and therefore we expect the uncertainties in fault geometry will have little impact here. The inferred slip model reproduces the SZU (Figure 4). We image a primary zone of fault slip in most of the offshore region, with a large uncertainty of 2-to-4 m (and up to 10 m in the near-trench domain, Figure S15 in Supporting Information S1). Downdip of this primary region of slip, at ∼90-km-depth, we infer a well-constrained slip zone with an amplitude of 2.5-3 m, equivalent to M w = 7.2, which is responsible for the secondary uplift. Models assuming a layered or homogeneous crust do not image this downdip slip and do not reproduce the SZU (Figure 1c and enclosed references, Figures S12-S14 in Supporting Information S1). Models assuming an heterogeneous elastic structure, but neglecting related epistemic uncertainties, are able to reproduce the SZU albeit not as well as when epistemic uncertainties are accounted for (Figures S13 and S14 in Supporting Information S1).  Figure  S5 in Supporting Information S1). In panels (b) and (c), the assumed fault replicates the true geometry shown in panel (a), but extends to greater depths. In panel (c), uncertainties in elastic properties are accounted for: Note the difference in the spatial distribution of posterior uncertainties. (d) Assumed 3D elastic properties, μ 0 = 52 GPa, which differ from the properties used to calculate synthetic observations (displayed in Figure 2f). (e) Trench perpendicular profile of the target synthetic data and predicted vertical displacements (at 0-km-along-strike). Vertical error bars indicate the posterior uncertainty. Predictions in light red are for the model shown in Figure S9 in Supporting Information S1.
Our results suggest that previously published models for the Maule earthquake were not able to reproduce the SZU (Figure 1c) because most of them were inferred assuming a layered crust. While Moreno et al. (2012) assumed 3D heterogeneous elastic properties, the shallow fault geometry they used and the impact of spatial regularization likely prevented a downdip patch to be imaged. Note that some authors do infer downdip slip as required by horizontal displacements (e.g., Bedford et al., 2013;Delouis et al., 2010;Vigny et al., 2011;Yue et al., 2014), but that the inferred slip could not cause a SZU because of the incompatibility shown in our synthetic example (Figure 3a). The combined effect of strong assumptions on the crustal elastic structure and fault geom etry, and the common use of unphysical regularization (e.g., Ortega-Culaciati et al., 2021), probably prevented published models from producing the mild secondary uplift of the Tohoku earthquake. As with the Maule earthquake, some authors do infer downdip slip as required by horizontal displacement without being able to fit the vertical measurements (e.g., Periollat et al., 2022).

What Is the Secondary Zone of Uplift?
That we image downdip slip does not mean slip is uniquely the cause of the SZU. Challenges in modeling highly disparate time-scales (from seconds to years) prevented van Dinther et al. (2019) from confirming the universal process they invoke is coseismic, rather than lasting several weeks after the mainshock. For the 2010 Maule earthquake, we infer downdip slip at ∼90-km-depth, where only a few aftershocks occurred, none with  horizontal survey GNSS postseismic offsets measured 12 days after the mainshock ( Figure S16 in Supporting Information S1).
The SZU observed after megathrust earthquakes other than the Maule event is located 200 from the trench in Chile, 350 in Japan, and 400 km in Alaska (van Dinther et al., 2019). Assuming that the SZU finds its origin in slip downdip of the coseismic rupture, because of the various slab geometries, the downdip slip would have consistently occurred at ∼90-120-km-depth. For the 1960 M w 9.5 Valdivia and 1964 M w 9.2 Alaska earthquakes, leveling data measured a few months to years after the mainshocks will probably contain a large postseismic component (e.g., Plafker, 1965;Plafker & Savage, 1970;van Dinther et al., 2019). In contrast, that coseismic geodetic data for the Maule and Tohoku earthquakes recorded the SZU would suggest it is coseismic or has been produced by very rapid postseismic slip (hours to weeks after the mainshock), which is often included in coseismic geodetic offsets.
For the Maule earthquake, the SZU has been mostly recorded by survey GNSS which were acquired several days to weeks after the mainshock (Vigny et al., 2011), and therefore contain some postseismic signal. At two continuous stations located in the region of the SZU (MAUL and ANTC, Figure 5), coseismic vertical offset is of 5 ± 9 and −16 ± 11 mm, while offsets at collocated survey stations (CT70 and LLA0, Figure 5) reach 102 ± 14 and 120 ± 13 mm (Vigny et al., 2011). Such difference would indeed suggest that the ∼10 cm SZU has been caused by postseismic deformation in the weeks following the mainshock. However, daily time-series estimated at the same continuous GNSS station do need 2 years to reach 10 cm uplift (MAUL and ANTC, Figure 5, Klein et al., 2022). Coseismic survey offsets published by Vigny et al. (2011) were calculated by extrapolating interseismic velocities over 10 years. The few measurements available (only three data points in 1996, 1999, and 2002, e.g., Ruegg et al., 2009) and non-negligible seasonal variations (>20 mm in amplitude at MAUL and ANTC, Figure 5), have likely altered the sparse interseismic velocity measurements, and make the errors on survey vertical offsets larger than those on continuous data. 3D displacement fields extracted from InSAR data also contain some postseismic deformation and show, in the SZU location, from −50 to +20 cm of vertical offset (Xiong et al., 2022) and are therefore not reliable to investigate the SZU.
For the Tohoku earthquake, Ozawa et al. (2011) measured up to 4 cm of coseismic uplift, and up to 5 cm of uplift in the 14 days after the mainshock, but not necessarily at the same locations as the coseismic offsets. Daily time series of Periollat et al. (2022) show up to 45 ± 5 mm uplift in the 1-3 days following the mainshock ( Figure 5). While some daily positions could suggest a 3-day transient postseismic uplift, the vertical component of their 30-s time series has a poor signal-to-noise ratio and cannot be exploited. 3D displacement field derived from InSAR data do reproduce some far-field uplift, but is not independent from measured GNSS offsets (J. Hu et al., 2013).
We demonstrate that the SZU is likely caused by downdip slip happening in the hours to days following the mainshock. While, for the Maule earthquake, the SZU remains ambiguous and a coseismic cause cannot be discarded, it is clearly measured in the hours to days after the Tohoku earthquake. Any further conclusion cannot be made without a thorough examination of early postseismic GNSS time series with a rate higher than 1 day, which is beyond the scope of this study.

Discussion and Conclusion
A SZU has been observed after several megathrust earthquakes. In this study, we investigate if (and which) assumptions in the foward and/or inverse approach could prevent the SZU to be reproduced with slip on the slab interface. We show that neglecting variations in elastic properties due to the plunging slab induces an incompatibility in the amount of slip required to explain the measured horizontal, or vertical, displacements, preventing models from reproducing the SZU. In contrast, we demonstrate that assuming realistic heterogeneous elastic properties, a sufficiently deep fault geometry, and discarding any non-physical regularization of the inverse problem, we infer the SZU as caused by slip downdip of the main coseismic rupture.
Inconsistencies in the fit to vertical versus horizontal measurements have already been discussed for various subduction zones and processes. For instance, Klein et al. (2018) report an inconsistency in the amount of slow slip needed to fit horizontal versus vertical observations a few hundreds of km from the trench. Some postseismic models of the Maule event (e.g., Lin et al., 2013), or synthetic tests performed for an infinitely long megathrust (Hsu et al., 2006), report similar inconsistencies. It is common practice to discard or down-weight vertical data 9 of 11 because of such inconsistencies and larger measurement errors. We show that by accounting for heterogeneities in elastic structure, we can reconcile vertical and horizontal observations.
With synthetic tests and a study of the 2010 M w 8.8 Maule earthquake, we suggest that the SZU is likely caused by deep postseismic slip (or localized strain) happening within the first hours to days following the mainshocks. For both the Maule and 2011 M w 9.1 Tohoku earthquakes, the ambiguity of the SZU measurements highlights the difficulty to accurately evaluate the contribution of very early deformations occurring after large earthquakes (Twardzik et al., 2019). Our results motivate further study of the postseismic phase as early as possible after the mainshock, as already emphasized by several authors (e.g., Jiang et al., 2021;Ragon et al., 2019;Twardzik et al., 2019).
While the occurrence of very early deep afterslip (hours to days after the mainshock, ∼100-km-deep) remains to be further investigated, it is consistent with a rate strengthening frictional behavior of the megathrust. For instance, numerical simulations of Muto et al. (2019) and Barbot (2020) showed that stress-driven aseismic afterslip can occur at great depths (60-100-km-depth) when considering rate-and-state friction laws. However, the amount of afterslip simulated with velocity-strengthening friction by van Dinther et al. (2019) was too low to reproduce the measured uplift. Alternatively, viscous flow could also explain such early postseismic deformation (e.g., Montési & Hirth, 2003). Mallick et al. (2022) show that power-law viscous flows are of greater amplitude at shorter time-scales for large earthquakes, which might explain why the SZU has only been observed for megathrust earthquakes. The occurrence of early viscous flow (hours to days) is consistent with longer-term steady-state viscoelastic relaxation invoked for the Maule earthquake (e.g., Klein et al., 2016;Peña et al., 2020Peña et al., , 2021, However, the similarity of surface displacements produced either by afterslip or by viscous flows prevents discriminating potential processes driving the SZU (e.g., Mallick et al., 2022;Weiss et al., 2019). In contrast, that the early SZU measured after the Tohoku earthquake is followed by subsidence may suggest the physical mechanisms driving the secondary uplift are different than for Chile, although viscous flows have been invoked (e.g., Agata et al., 2019;Luo & Wang, 2021;Sun et al., 2014).

Data Availability Statement
Materials (without data sets) presented in this paper are archived and available on Zenodo (Ragon, 2022). Static GNSS offsets for the 2010 Maule earthquake have been published in Vigny et al. (2011) and Lin et al. (2013). GNSS time series have been process by Klein et al. (2022), with data provided by Centro Sismológico Nacional (CSN) of the Universidad de Chile (Báez et al., 2018), that can be retrieved from the GNSS database (http://GNSS. csn.uchile.cl). GNSS time series for the 2011 Tohoku earthquake have been processed by Periollat et al. (2022). Static GNSS time series are attributed to GNSS Products, CNRS, OSUG, ISTERRE (2022). The Bayesian simulations were performed with the AlTar2 package (https://github.com/lijun99/altar2-documentation). The Classic Slip Inversion (CSI) Python library (Jolivet et al., 2014) developed by Romain Jolivet was used to build inputs for the Bayesian algorithm. The mesh for the FEM simulations was built using Coreform Cubit (v2022.4, retrieved from http://coreform.com). We used the finite element code Pylith (Aagaard et al., 2013) to perform the simulations. Slab geometry, topography, and crustal elastic properties from Slab2, LITHO1.0, and ETOPO1 models are available in Hayes et al. (2018), Pasyanos et al. (2014), and Amante and Eakins (2009). 3D data were visualized using the open-source parallel visualization software ParaView/VTK (Ahrens et al., 2005). Figures were generated with the Matplotlib (Hunter, 2007) and Seaborn (Waskom, 2021) Python3 libraries.