Earthquake ruptures and topography controlled by plate interface deformation

What controls the location and segmentation of mega-earthquakes in subduction zones is a long-standing problem in earth sciences. Prediction of earthquake ruptures mostly relies on interplate coupling models based on Global Navigation Satellite Systems providing patterns of slip deficit between tectonic plates. We here investigate if and how the seismic and aseismic patches revealed by these models relate to the distribution of deformation along the plate interface, i.e. basal erosion and/or underplating. From a mechanical analysis of the topography applied along the Chilean subduction zone, we show that 5 extensive plate interface deformation takes place along most of the margin. We show that basal erosion occurs preferentially at 15 km depth while underplating does at 35 ± 10 and 60 ± 5 km depth, in agreement with P-T conditions of recovered underplated material, expected pore pressures, and spatial distribution of marine terraces and uplift rates. Along southern Chile, large sediment input favors shallow accretion and underplating of subducted sediments, while along northern Chile, extensive basal erosion provides material for the underplating. We then show that all major earthquakes of southern Chile are 10 limited along their down-dip end by underplating while, along northern Chile, they are surrounded by both basal erosion and underplating. Segments with heterogeneously distributed deformation largely coincide with lateral earthquake terminations. We therefore propose that long-lived plate interface deformation promotes stress build-up and leads to earthquake nucleation. Earthquakes then propagate along fault planes shielded from this long-lived permanent deformation, and are finally stopped by segments of heterogeneously distributed deformation. Slip deficit patterns and earthquake segmentation therefore reflect 15 the along-dip and along-strike distribution of the plate interface deformation. Topography acts as a mirror of distributed plate interface deformation and should be studied systematically to improve the prediction of earthquake ruptures.


Methodology
To do so, we applied the critical taper theory (CTT, Davis et al., 1983, Fig. 2a). The theory predicts the topographic slope of a 60 wedge on the verge of failure as well as dips of its internal faults, from the frictional and pore fluid pressure properties of the wedge and megathrust (Dahlen, 1984). We used the solution of Dahlen (1984) for a non cohesive wedge, describing the critical taper as a function of the angle Ψ B formed by the maximum principal stress σ 1 and the base of the wedge and the angle Ψ 0 formed by σ 1 and the top of the wedge. The solution for the lower compressional branch is: The angles φ int and φ b are the internal and basal coefficients of friction defined as µ int = tanφ int and µ b = tanφ b , and The internal and basal Hubbert-Rubbey fluid pressure ratios λ and λ b are defined in (Davis et al., 1983) as: where ρ and ρ w are the wedge material and water densities and D is the water depth. The solution is exact if λ = λ b and the approximation is valid for small tapers as used in this study (Wang et al., 2006).
The relationship between the topographic slope α and the slab dip β forms an envelope separating different mechanical states 80 (Fig. 2a). The wedge is at critical state if the taper formed by α and β follows a critical envelope. In that case, activation of the megathrust requires internal faulting. In contrast, if the taper exceeds this critical limit, the wedge enters in a stable domain (in the sense of the CTT), where the only possible active fault is the megathrust.
Dynamic effective frictions of seismogenic zones have been shown to reach extremely low values (µ ∼ 0.01 − 0.03, Fulton et al., 2013;Gao and Wang, 2014). As a consequence, a forearc above a seismogenic megathrust, with standard α and β, will 85 systematically fall in the CTT stable domain (Fig. 2a). In contrast, aseismic megathrusts are characterized by larger effective friction (µ ∼ 0.1 − 0.15, Cubas et al., 2013;Gao and Wang, 2014), allowing the wedge to reach critical conditions (Cubas et al., 2013) (Fig. 2a). Furthermore, at critical state, if the basal effective friction of the wedge approaches the internal effective friction, the dip of internal faults within the wedge will decrease and become parallel to the plate interface. A small difference of friction will thus lead to deformation either by basal accretion (in the lower plate) or by basal erosion (in the upper plate), 90 instead of standard accretion characterized by thrust faults reaching the surface (Dahlen, 1984) (Fig. 2a). If the basal effective friction reaches the internal one, a highly fractured forearc is even expected.
Following the method of (Cubas et al., 2013), we used ETOPO 1 for the bathymetric and topographic slope, and slab 2.0 for the slab dip (Hayes et al., 2018). We built swath profiles perpendicular to the trench every 0.1 longitudinal degree (Suppl. Mat. Fig. 1), smoothed by a rectangular window, to get average topographic slopes and slab dips with their standard deviations 95 necessary to the inversion procedure (Fig. 2,Suppl. Mat. Fig. 2,5). Along these profiles, we selected segments parallel to critical envelops. For each of these segments, we retrieved, by inversion, probability density functions for the friction of the wedge, the pore pressure ratio of the wedge and the effective friction of the megathrust (Suppl. Mat. Fig. 3). The CTT inversion only provides the internal pore pressure, which can be considered as a lower bound for the megathrust pore pressure. Moreover, close effective internal and basal frictions implies close internal and megathrust pore pressure. The inversion procedure is 100 described in (Cubas et al., 2013). We only kept segments with extremely low misfits, and with values consistent with standard frictions (from 25 to 43 o for φ int , Byerlee (1978), from 1 to 42.9 o for φ ef f b , and from 0.35 to 0. 95 for λ). Sensitivity tests were run by Cubas et al. (2013), showing that an error of ±5 o on β implies a horizontal translation of the taper and a 3 o variation for the effective basal friction, without affecting the critical state of the forearc. Best misfits for each parameter are provided in Supplementary Material (Table 1 & Fig. 6).

Locating plate interface deformation
Since we seek to capture distribution of deformation along the plate interface, we compared the internal and megathrust effective frictions and mapped the difference (Fig. 1b). A one degree difference implies a dip difference < 10 o between the megathrust and the forward-verging thrusts. This representation allows inferring areas prone to either standard accretion, or to 110 basal erosion/accretion, the latter being herein defined as distributed deformation along the plate interface.
The transition from erosion to accretion is particularly well captured by the method despite its relative simplicity ( Fig.   1a-b). North of 25 o S (i.e., at the transition between erosive and accretionary subduction), the difference of effective frictions is systematically lower than 2 o , and mostly lower than 1 o , which means that deformation is essentially located along the plate interface. The difference increases notably south of the Juan Fernandez ridge, and progresses inland towards the south, 115 consistent with active faults mapped near the Arauco peninsula (Fig. 1d). The method also captures basal erosion linked to the few subducted seamounts imaged around the Iquique rupture (Geersen et al., 2015) (Fig. 1c).
To interpret the origin of plate interface deformation, we set the limit between accretion and interface deformation for a difference of friction angle of one degree and then searched for the depth distribution of each process ( Fig. 3a-b). Results for different thresholds are also presented.

120
Accretion mostly occurs at depths shallower than 20 km whereas distributed deformation along the plate interface is found at every depth with three favorable peaks: a first one at 15 km depth, a second at 35±10 km depth and a third at 60±5 km depth.
Basal erosion due to bathymetric features is known to be efficient at shallow depth (von Huene et al., 2004), but might also occur slightly deeper as testified by subsidence of forearc basins (Clift and Hartley, 2007). Underplating has been observed at shallow depth (Kimura et al., 2010;Tréhu et al., 2019;Bangs et al., 2020), even for seamount moats (Clarke et al., 2018), 125 but mostly occurs below 30-40 km, i.e. the long-term coupling-decoupling transition (Agard et al., 2018). The compilation of P-T conditions for recovered underplated material shows a similar distribution as CTT, with peaks at depths of 30±5 km and 50±5 km and a gap in between (Fig. 3b). Numerical modelling also suggests that erosion predominates over underplating at depths shallower than 20 km (Menant et al., 2020). As a consequence, we assume that the interface deformation documented here mostly relates to basal erosion for depths < 20 km, whereas underplating dominates at greater depths, in particular below Most of the critical segments attesting to distributed interface deformation have lengths smaller than 7 km, and do not exceed 20 km (Fig. 3c). This is in good agreement with the length of graben and horst and seamounts observed on the subducting plate (von Huene and Ranero, 2003;Geersen et al., 2015). This is also consistent with seismic observations (

Relating plate interface deformation with long-term coastal uplift 145
In order to approach internal effective friction, the basal friction of the wedge needs to be relatively high (Fig. 2a). Deformation captured at depth by CTT therefore corresponds to underplating in the making, when deformation is still distributed: once tectonic slicing is achieved, basal friction will reach low values characteristic of mature faults, bringing the taper in the stable domain and impeding detection by CTT.
We here show that active, incipient underplating takes place along most of the Chilean margin (Fig. 1, 4a), whether or not 150 strain accumulation ultimately turns these critical zones into tectonic slices (Agard et al., 2018;Bonnet et al., 2019;Kimura et al., 2007). Along the Tocopilla and Maule segments, incipient underplating takes place below the slab -continental Moho intercept, whereas it is located above and near the intercept between 25-33 o S. Along the northern erosive margin (Fig. 4a-b), the thin oceanic sediment cover together with eroded material and possibly pieces stripped from the subducting plate are being underplated. Southward, larger sediment input favors shallow accretion and underplating of subducted sediments ( Fig. 4a-b), To evaluate the possible link between coastal uplift and regions of distributed deformation along the plate interface, we compare our spatial distribution of underplating with documented marine terraces and uplift rates (Saillard et al., 2017) (Fig. 4a-c). While the irregular distribution of marine terraces may relate to fragmentary preservation rather than discontinuous underplating, most of them coincide with segments of underplating (or with standard accretion near the Arauco peninsula).

160
The typical length of the segments showing distributed deformation along the interface is modest (∼10-20 km). This advocates for propagation of deformation only into the topmost part of the subducting plate, e.g. into the stack of oceanic/eroded material or along a hydrated basalt layer within the crust. Although the comparison with uplift rates inferred from morphometric analysis and modelling of landscape evolution (Melnick, 2016) shows no particular pattern, underplating of relatively thin slices would be consistent with both the low uplift rates and the rock record (Agard et al., 2018). 165 Shallow normal faults identified along northern Chile have been related to potential underplating (Adam and Reuther, 2000;Clift and Hartley, 2007). We here show that they are located above underplating areas, and that normal faults transition southward to thrusts where accretion prevails (Fig. 4a-c).

Large megathrust earthquakes are surrounded by plate interface deformation
We now compare zones of distributed deformation with areas struck by recent large earthquakes (Fig. 1c) shows the diversity of published co-seismic models, emphasizing the uncertainties related to the proposed slip areas, that need to be taken into account for the comparison.
To summarize, along the accretionary part of the margin, mega-earthquakes are delimited down-dip by underplating whereas 180 along the erosive domain, the extent is controlled by both basal erosion, at rather shallow depth, and underplating.
We now investigate the relationships between the critical segments characterized by relatively high effective basal friction, hence expected to behave aseismically, both with the along-strike segmentation inferred for large historical earthquakes (Mw respectively: Fig. 4a-d), particularly in the southern region. The distribution of underplating also shows a long wavelength (4 to 8 o degrees latitude) corresponding to the three main seismic domains visible on Fig. 4d (respectively located north of Mejillones, between Mejillones and Punta Choros, and to the south of Punta Choros). No correlation exists with sediment input (Fig. 4b). Marine terraces are mostly found where underplating is maximum, i.e. in the central domain (Fig. 4c). The proportion of a swath profile at critical state is then compared to the average estimate of coupling along this profile (Fig. 4e). A  fig.7), an anti-correlation is observed with coupling: zones of high coupling correspond to domains where underplating is somewhat more limited (e.g., between 31 and 38 o S). As shown by Fig. 1d, areas marked by intense plate deformation are mostly located down-dip of the highly locked patches, in particular in the south. Continuous and 195 voluminous sediment influx provides spatially and temporally stable mechanical conditions (Olsen et al., 2020). This leads to stationary seismic asperities followed, once the Moho is crossed, by the regular underplating of very thin and relatively short slices (typically 300-500 m thick, 5-10 km long). This spatial and temporal stability promotes topographic build-up, explaining the spatial coincidence between strongly coupled patches and efficient underplating. In contrast, segments with

Link between earthquake ruptures and plate interface deformation
We herein show two major results: (1) recent earthquakes are bounded by extensive plate interface deformation characterized, along the southern accretionary part of the margin by underplating at their down-dip edge and, along the northern erosive part of the margin by both basal erosion and underplating (Fig. 1); (2) along strike, segments characterized by minimal interplate 210 deformation largely coincide with lateral earthquake terminations (circles; Fig. 4d). These observations demonstrate a close link between distributed interplate deformation and earthquake segmentation.
We therefore propose that plastic deformation and stress build-up associated with interplate deformation along distributed fault planes of limited extent (Fig. 6 -step 1), eventually leads to mega-earthquake nucleation (Fig. 6 -step 2). This process may be accompanied by transient seismic/slip events hardly resolvable by geodetic means due to their size, depth or distance to the 215 coast, as expected in zones of distributed deformation and of high pore fluid pressure (Liu and Rice, 2007;Kimura et al., 2010;Collot et al., 2017). This is suggested by the location of hypocenters at boundaries with distributed deformation (Fig. 1) and concords with observations for the Iquique earthquake (Ruiz et al., 2014a;Meng et al., 2015;Socquet et al., 2017) and others elsewhere (Tohoku, Kato et al., 2012), (Guerrero, Radiguet et al., 2016). Once nucleated, large earthquakes propagate along well localized and smoothed rate-weakening fault planes (Bletery et al., 2016) bounded by elongate zones of underplating (or 220 basal erosion when present; Fig. 1, 6 -step 3). These earthquakes abut on regions of heterogeneously distributed deformation and stress concentrations (Fig. 6 -step 4), which inhibit the development of well localized slip zones (Wang and Bilek, 2014) and impede rupture propagation (Wesnousky, 2006). Some of these fractures might however slip seismically during large events and produce the observed high frequency radiations (Meng et al., 2018(Meng et al., , 2015.

Down-dip segmentation and plate interface deformation 225
Within this framework, domain-C (Lay et al., 2012) would correspond to the region of incipient underplating, and domain-B to the smoothed rate-weakening fault zone. This is consistent with partial coupling of domain-C, where deformation is dominated by creep mechanisms, and elastic strain is accumulated along undeformed segments of the megathrust, producing moderate slip earthquakes with a higher recurrence time and significant coherent short-period seismic radiations. Domain-C might more specifically be delimited by the two underplating peaks (Fig. 3), the upper one corresponding to the slicing of 230 oceanic sediments and/or eroded continental material, the second to the propagation of deformation into the hydrated basalt layer once the continental Moho is crossed (Agard et al., 2018), a process ultimately leading to the underplating of very thin and relatively short slices (typically 300-500 m thick, 5-10 km long, Agard et al., 2018).