# The cryptic seismic potential of the Pichilemu blind fault in Chile revealed by off-fault geomorphology

J Jara-Muñoz, D Melnick, S Li, Anne Socquet, J Cortés-Aranda, D Brill, M R Strecker

► To cite this version:

J Jara-Muñoz, D Melnick, S Li, Anne Socquet, J Cortés-Aranda, et al.. The cryptic seismic potential of the Pichilemu blind fault in Chile revealed by off-fault geomorphology. Nature Communications, 2022, 13 (1), 10.1038/s41467-022-30754-1. http://dx.doi.org/10.1038/s41467-022-30754-1

HAL Id: hal-04979004

https://hal.science/hal-04979004v1

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-022-30754-1

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# The cryptic seismic potential of the Pichilemu blind fault in Chile revealed by off-fault geomorphology

J. Jara-Muñoz ID <fnref n="1" />,⁷✉, D. Melnick ID ², S. Li ID ³, A. Socquet ID ⁴, J. Cortés-Aranda ID ⁵, D. Brill ID ⁶ & M. R. Strecker<fnref n="1" />

The first step towards assessing hazards in seismically active regions involves mapping capable faults and estimating their recurrence times. While the mapping of active faults is commonly based on distinct geologic and geomorphic features evident at the surface, mapping blind seismogenic faults is complicated by the absence of on-fault diagnostic features. Here we investigated the Pichilemu Fault in coastal Chile, unknown until it generated a Mw 7.0 earthquake in 2010. The lack of evident surface faulting suggests activity along a partly-hidden blind fault. We used off-fault deformed marine terraces to estimate a fault-slip rate of 0.52 ± 0.04 m/ka, which, when integrated with satellite geodesy suggests a 2.12 ± 0.2 ka recurrence time for Mw~7.0 normal-faulting earthquakes. We propose that extension in the Pichilemu region is associated with stress changes during megathrust earthquakes and accommodated by sporadic slip during upper-plate earthquakes, which has implications for assessing the seismic potential of cryptic faults along convergent margins and elsewhere.

<footnote n="1">¹Department of Earth Sciences, University of Potsdam, Karl-Liebknecht-Str. 24-25, Potsdam, Germany. ²Austral University of Chile, Institute of Earth Sciences, Edificio Pugín, Campus Isla Teja, Valdivia, Chile. ³State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, No. 19, Beitucheng Western Road, Beijing, China. ⁴University of Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, UGE, ISTerre, 38000, Grenoble, France. ⁵Department of Earth Sciences, Universidad de Concepción, Víctor Lamas 1290, Concepcion, Chile. ⁶Institute of Geography, University of Cologne, Otto-Fischer-Straße 4, Cologne, Germany. ⁷Present address: Department of Civil Engineer, Hochschule Biberach, Karlstraße 9-11, Biberach, Germany. ✉email: jara@geo.uni-potsdam.de</footnote>
Unexpected ruptures that occur during large-magnitude earthquakes along previously unmapped faults emphasise a major lacuna in our knowledge concerning the location and seismic potential of tectonically active structures¹. Over the past decade, earthquake ruptures have occurred along previously unidentified or not fully mapped faults during at least five Mw > 6 earthquakes. These include the 2010 Pichilemu earthquakes in Chile, the 2016 Kaikoura earthquake that ruptured the unmapped Papatea Fault, the 2010 Darfield and 2011 Christchurch earthquake sequence in New Zealand, and the 2019 Ridgecrest events in California²⁻⁶. These earthquakes highlight an important gap in our understanding of the seismogenesis of hidden faults in a variety of geodynamic environments. Whereas mapping potentially active faults commonly relies on identifying geomorphic and geologic features indicative of surface rupture and deformation that could be associated with past earthquakes⁷,⁸, in areas with active blind faults such evidence may either be completely absent or difficult to identify⁹,¹⁰. Blind faults are geological structures whose ruptures do not reach the earth's surface¹¹, thereby hiding their seismogenic potential. Such structures are common in sedimentary basins and have often been identified on the basis of geophysical imagery¹²⁻¹⁴ and indirect geomorphic observations¹⁵. However, estimating the seismic potential of blind faults is difficult in the absence of any on-fault geologic or geomorphic evidence from past earthquakes. Here we demonstrate that quantifying deformation using off-fault geomorphic strain markers provides valuable insight into cryptic, potentially active faults in coastal areas. We use the terms on-fault and off-fault to distinguish between surface deformation that has occurred along the fault trace and within the zone surrounding the fault, respectively.

The polarity and spatial distribution of stresses in the upper plate of subduction zones changes throughout the seismic cycle¹⁶. For instance, in response to the nearly instantaneous polarity change associated with a megathrust earthquake, the upper plate is commonly affected by enhanced extension that results in increased seismicity and the triggering of occasional crustal earthquakes¹⁷,¹⁸. Slip on upper-plate faults triggered by megathrust earthquakes has been reported in Japan, Alaska, and Chile¹⁸⁻²⁰, and has been inferred to be a common feature along most subduction zones¹⁸,²¹⁻²³. However, historical and paleoseismic observations suggest that crustal faults with a low slip rate are characterised by recurrence times involving thousands of years, and may therefore not be reactivated during every megathrust earthquake, since these commonly recur over periods ranging from a number of decades to a few centuries¹⁸,²²,²⁴. Because both exposed and hidden upper-plate faults are widespread along coastlines bordering subduction zones, they pose significant local hazards. Such crustal faults can produce higher amplitude seismic waves at local scales than megathrust earthquakes and may locally increase the amplitude and shorten the arrival times of tsunamis in the near-field²⁵. Mapping crustal structures along subduction zones and quantifying their slip rates and relationships to megathrust earthquake cycles is therefore a fundamental requirement for obtaining an adequate assessment of the spatiotemporal characteristic of earthquake and tsunami hazards.

In this study we focus on the Pichilemu Fault (PIF), a hidden fault that was unknown until it generated two shallow Mw 7 and 6.9 earthquakes 11 days after the Mw 8.8 Maule megathrust event that affected central Chile in 2010. Our study combines geomorphic and morphometric analyses using high-resolution LiDAR topography, luminescence dating, radar interferometry, and numerical modelling. We show that while on-fault displaced geomorphic markers are absent along the surface fault traces, off-fault strain markers can be used to estimate a long-term slip rate,

which, when integrated over the 2010 coseismic deformation pattern and assuming a characteristic slip behaviour, allows a recurrence rate to be inferred for such earthquakes. Our results demonstrate the hidden seismogenic potential of blind faults, with implications for seismic hazard along coastlines bordering subduction zones.

## Results and discussion

**Seismotectonic and geologic setting.** The PIF is located at the coast of the central Chile margin, where the Nazca plate is subducting beneath South America at 66 mm/yr²⁶ (Fig. 1A). This region comprises the Coastal Range, which reaches maximum elevations of ~600 metres above sea level (MASL) and consists almost exclusively of crystalline metamorphic rocks and scattered intrusive bodies; the metamorphic rocks are related to a Paleozoic accretionary prism overprinted by brittle deformation during Mesozoic and Cenozoic exhumation²⁷⁻²⁹ (Fig. S1A). The seaward slope of the range is sculpted by a sequence of uplifted marine terraces, some of them overlain by shallow marine deposits³⁰. Dense vegetation and thick soil cover have hampered the mapping of geological and geomorphic features in the area, resulting in different interpretations regarding the presence of tectonically active structures e.g., refs. 4,30.

The Pichilemu area was affected by the 2010 Maule earthquake (Mw 8.8), which ruptured a ~500 km-long portion of the megathrust, with a peak slip of 17 m to the south of the PIF³¹ (Fig. 1A). The Maule earthquake triggered instantaneous slip along the Santa Maria fault and non-instantaneous slip along the PIF, which are located in the southern and northern parts of the rupture zone, respectively⁴,¹⁹ (Fig. 1A). Prior to the 2010 earthquake, only a few faults affecting the crystalline basement and Cenozoic sedimentary cover were known in this area²⁸,³² (Fig. 1B); the PIF was unknown.

Eleven days after the Maule earthquake, the PIF slipped during two Mw 6.9 and 7 normal-faulting earthquakes (Fig. 1B), followed by ~12,000 aftershocks located between the megathrust and a depth of ~4 km³³ (Fig. 1C). The aftershocks delineated a NW-striking, SW-dipping structure extending for ~80 km along strike. Surface displacements estimated from GPS and ALOS/PALSAR Interferometric Synthetic Aperture Radar (InSAR) collected 2 days before and 44 days after the Maule earthquake suggest that the 2010 PIF earthquakes were associated with a maximum slip of ~3 m along the main strand of the PIF, extending from 5 to 22 km depth³⁴. However, despite the large magnitude of the PIF earthquakes, no evidence of surface ruptures was found during field surveys nor were any detected in radar interferometry images⁴,³⁴, which suggests that the PIF is a blind, yet tectonically active structure.

**On-fault tectonic geomorphology.** We used a Digital Terrain Model (DTM) derived from airborne Light Detection And Ranging (LiDAR) data at 1 m resolution to estimate fluvial metrics³⁵ and analyse the surface expression of the PIF (see Methods Section, Analysis of on-fault geomorphic features, and Figs. S1A–F and S2A–E). The coastal reaches of the PIF comprise two catchments of ~100 km² each (C1 and C2; Fig. 2A) developed almost exclusively on metamorphic bedrock, except for the western part of C2, which contains isolated intrusive outcrops, and near the coast, where the valley floors are filled by marine, fluvial, and aeolian deposits (Fig. S1A). The median annual rainfall is similar in both catchments, with values of 622 mm/yr in C1 and 618 mm/yr in C2, based on data from the Tropical Rainfall Measuring Mission (TRMM)³⁶ (Fig. S2E). The catchment asymmetry is evident in each of these catchments and manifested in the deviation between the main trunk stream and the
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Tectonic setting of the Pichilemu Fault. **A** Slip distribution for the 2010 Maule earthquake³¹ and active faults from the Chilean Active Fault database⁵⁸. The Santa María Fault (SMF) and Pichilemu Fault (PIF) were reactivated by the earthquake⁴,¹⁹. YF: El Yolki Fault; CF: Carranza Fault. **B** Topography (SRTM data from www2.jpl.nasa.gov), faults, and focal mechanisms of the PIF earthquakes⁴ in the Pichilemu region. Contours show aftershock density after the Maule earthquake. Thick black lines denote the PIF branches (F1 and F2). Also shown are the Topocalma Fault (TOF), a Quaternary normal fault, and Vichuquén Fault (VF), an inferred fault affecting the crystalline basement⁴. **C** Crustal profile showing aftershocks (Mw > 1) recorded between March 15 and September 30, 2010³³, used to infer the subsurface geometry of the PIF. Megathrust from Slab 2.0 model⁴⁷.</figcaption></figure>

catchment centreline, which is quantified using the symmetry factor (Ts)³⁷ (see Methods Section, Analysis of on-fault geomorphic features). The asymmetry of C1 and C2 is highlighted by trunk streams that converge along a section parallel to the trace of the PIF (Ts > 0.6); further east, the catchments become

progressively symmetric with Ts~0.1 (Fig. 2A and Fig. S2B). The area to the west of the centreline in C2 includes NE-SW-elongated parallel drainages that are probably associated with local surface tilting. Catchment C1 is characterised by higher local relief than C2 with median values of 111 and 135 m, respectively, and reaching 300 m in C1 and 250 m in C2 (Fig. S2C). The median slope of both catchments is similar (13° for C1 and 14° for C2), but its distribution is slightly biased towards lower values in catchment C2 (Fig. S2D). The drainage network of catchment C1 reaches ~500 MASL and includes 16 knickpoints distributed between 100 and 350 MASL (Fig. 2C). Steepness index ($K_{sn}$) values reach up to 100, with the higher values forming a fringe at elevations between 300 and 400 MASL. C2, in contrast, includes eight knickpoints distributed between 100 and 150 MASL and has $K_{sn}$ values of up to 70 (Fig. 2A). A conspicuous set of NW-SE and NE-SW-oriented lineaments can be identified from aligned drainages and small slope breaks. The westernmost lineaments (L1 and L2; Fig. 2A and Fig. S1), are ~3 km-long and have ~50 m-high scarps partly degraded by river incision and associated with contacts between metamorphic and intrusive rock units (Fig. S1A). Their traces are oblique to fault geometries inferred from crustal seismicity (F1 and F2). Farther east, lineaments L3, L4 and L5 are highlighted by aligned valleys and a trellis drainage pattern (Fig. 2A, E and Fig. S1). These lineaments extend for ~7 km along strike, are associated exclusively with metamorphic rocks, and are subparallel to local metamorphic foliations²⁷,³⁸ (Fig. S1A), but oblique to the trace of F1.

Catchment and drainage metrics indicate variable degrees of surface deformation in this area. For instance, the high symmetry factors of both catchments together with the subparallel drainages in the western part of C2 suggest local tilting of the PIF footwall and hanging-wall blocks in opposite directions (Fig. 2 and Fig. S2B). The difference in local relief and catchment slope between C1 and C2 reflect the varying degrees of river incision, which, together with the differences in $K_{sn}$ values, knickpoint locations, and drainage elevations, may indicate differential vertical displacements. The minimum differences between the catchments in terms of rainfall and metamorphic-bedrock lithology (100% of the area in C1 and ~75% in C2; Fig. S1A) suggest that these two factors exert a negligible control on catchment asymmetry and drainage metrics.

A detailed analysis of the PIF lineaments suggests that these features cannot be directly interpreted as fault scarps associated with surface ruptures during recent earthquakes. The drainages crossing lineaments L1 and L2 are associated with low $K_{sn}$ values, suggesting that if these lineaments were related to active faults, either the faults are slipping at low rates that have a little effect on river incision, or they are associated with a blind fault that does not reach the surface. Furthermore, the lineaments coincide with lithological contacts, indicating that their morphology might respond to differential erosion rather than the effects of surface faulting. Considering their orientations and morphologies, lineaments L3 to L5 are the best candidates for fault scarps. However, even though they are associated with subtle breaks in slope, a well-developed scarp cannot be distinguished on scarp-perpendicular swath or ridge-crest profiles (Fig. 2E, F). Furthermore, no clear spatial relationship exists between the trace of these lineaments and the fringe of higher $K_{sn}$ values of C1 regarding overlap or orientation (Fig. 2A).

Chi-plots are commonly used to identify transient signals propagating upstream along fluvial systems, such as tectonically generated knickpoints³⁹ (see Methods Section, Analysis of on-fault geomorphic features for further details). However, we found no relationship between chi-plot ($\chi$) values at the knickpoint locations in C1 and the distance from the knickpoints to the mapped lineaments (Fig. 2B and Fig. S2A). Furthermore, the
<figure><img src="image_5.png" /><figcaption>**Fig. 2** On-fault tectonic geomorphology from bare-earth LiDAR topography. **A** Shaded-relief map of the PIF coastal reach showing the two analysed catchments (C1 and C2), the thick dashed black lines are the faults (F1 and F2) inferred from spatial alignments of crustal seismicity. Note asymmetry of the main trunk streams with respect to the catchment centreline. The black lines labelled L1 to L5 denote lineaments that may represent surface-breaching faults (see text for discussion). Steepness index ($K_{sn}$) determined for a reference concavity of 0.45. **B** Scatter plots and histograms comparing fluvial metrics with the lineament metrics within C1. Scatter plot and histograms of chi-values ($\chi$) at knickpoints with respect to distance to lineament (See chi-map Fig. S2A). **C** Longitudinal river profiles of both catchments showing knickpoint locations. **D** Slope-area plots for both catchments. Note that C1 has a higher $K_{sn}$ value suggesting more rapid uplift. **E** Red-relief map of lineaments, ridge profiles (green lines) and swath boxes (black rectangles) shown in **F**. The arrows indicate the trace of the lineaments. **F** Swath and ridge-crest profiles. Profiles 1, 2, a, and b are centred on the trace of lineament L3, and profiles 3 and c are centred on the trace of L5. Note the lack of evident scarps suggesting an absence of any recent surface-breaching ruptures.</figcaption></figure>
linear distances between knickpoints and the L3 and L5 lineaments do not display a clear trend (Fig. 2B), suggesting that knickpoints are not related to the potential active fault scarps but rather reflect the effects of base-level changes associated with relative sea-level variations. The results of our morphometric analyses of LiDAR topography indicate that the area has been affected by surface deformation with a degree of spatial asymmetry, which may be a result of tilting and differential uplift. However, these deformation patterns cannot be directly related to any particular structure with a marked surface expression. We therefore conclude that this area is not characterised by localised deformation at the surface, but rather by strain distributed over a 10 km-wide region.

**Off-fault tectonic geomorphology.** Uplifted marine terraces—geomorphic markers of past relative sea-level positions⁴⁰—are ubiquitous along the coast of central Chile³⁰ and can be used as regionally correlatable strain markers⁴¹. We mapped terraces at Pichilemu using a LiDAR DTM and the TerraceM-2 software⁴² (Fig. 3A, see Methods Section, Analysis of off-fault geomorphic features). To the south of the PIF there are four levels of wave-cut terraces sculpted into the bedrock reaching up to 100 MASL (Fig. 3B and Fig. S3); these are occasionally covered by a thin veneer of marine and aeolian sediments. In contrast, to the north of the PIF six distinct sedimentary units corresponding to wave-built terraces are exposed; they consist of shallow marine sandstone bodies onlapping against the crystalline bedrock at elevations between 50 and 170 MASL, each comprising a single regressive cycle (Fig. 3B, Fig. S4A–C). From the lower and intermediate terrace levels at 51 and 115 MASL, respectively, we obtained post-IR IRSL (post-infrared infrared stimulated luminescence, see Methods Section, Post-IR IRSL dating) ages of 106 ± 9.3 and 297 ± 29 ka, which correspond to Marine Isotope Stages (MIS) 5d and 8 (Fig. 3A, Figs. S4A, S5, and S6, Table 1). These sediments were deposited above bedrock during shoreline progradation shortly after the corresponding MIS highstand e.g., refs. ³⁰,⁴³. We have correlated the surface morphology and geometry of these deposits with MIS 5e and MIS 9, at 125 and 320 ka respectively (Fig. 3A, B, Figs. S4, and S6). We have also tentatively correlated two additional terrace levels with MIS 7 and MIS 11 based on a composite sea-level curve (Fig. S6), these relative sea-level highstands correspond to ages between 250 and 380 ka.

South of the PIF, the lower marine terrace level is continuously exposed between Punta de Lobos and La Puntilla, with widths ranging between 1 and 3 km. This surface has been previously interpreted as a rasa³⁰, i.e., a terrace surface formed by marine reoccupation during successive highstands. Marine sediments that cover this rasa level between 16 and 32 MASL have yielded post-IR IRSL ages of 328 ± 33 and 314 ± 30 ka, corresponding to MIS 9 (Fig. 3A, Figs. S5 and S6). The upper terrace levels are characterised by well-defined paleo-cliffs and narrow paleo-platforms, with mean shoreline-angle elevations of 42, 60, and 80 MASL decreasing southward. By correlating the sequence with global sea-level curves⁴⁴,⁴⁵, we interpret the age of the three upper terrace levels to range between MIS 11 and MIS 19 (Fig. S6).

The estimated uplift rates vary between 0.06–0.15 m/ka and 0.37–0.46 m/ka across the PIF (Fig. 3C), with associated 2 s errors between 0.01 and 0.08 m/ka, suggesting protracted emergence over the past ~620 ka. Overall, the marine terrace sequence displays a broad warping pattern with a wavelength of ~10 km, which is compatible with rapid uplift and back tilting along the PIF footwall block, and monoclinal rollover folding in the hanging wall, which is consistent with a NW-striking, SW-dipping normal fault at depth (Fig. 3B, C).

**Coseismic slip and long-term slip rate of the PIF.** Using a combination of coseismic displacements derived from GPS and INSAR together with fault geometries inferred from aftershock seismicity³³,⁴⁶, we estimated slip along the PIF during the 2010 earthquakes (Figs. 3D, 4A–B, Figs. S7A–E and S8A–F). The INSAR data comprised two Envisat scenes acquired 2 days before and 7 days after the PIF earthquakes, obtaining an ascending interferogram (See Methods Section, Estimating coseismic slip during the PIF earthquakes). The aftershocks clearly mark the down-dip termination of the PIF at ~26 km depth where it intersects the megathrust⁴⁷; in contrast, the up-dip limit is more diffuse and most likely located at a depth between 5 and 8 km (Fig. S7B–E). We carried out a set of forward models by iterating up-dip depths and slip magnitudes for F1 and F2 within pre-defined ranges (see Table S1 and Methods Section, Estimating coseismic slip during the PIF earthquakes). We varied the along-dip extent of the fault by moving the up-dip limit up and down along fixed dips (55° for F1 and 72° for F2), strikes (N38°W for F1 and N16°W for F2), and down-dip depths (26.2 km for F1 and 20 km for F2). The model allowed us to estimate a slip magnitude and up-dip depth for each fault, along with their corresponding uncertainties, using the normalised root mean squared error (NRMSE, Figs. S9A and S10A–D). The model results suggest that F1 slipped 1.1 m (with 90% confidence interval (CI) between 0.95 and 1.15 m) between the megathrust and a depth of 4.6 km (CI: 3.6–5 km, Figs. S9A and S10B); and F2 with 0.1 m of slip (CI: 0–0.4 m) and 5.2 km up-dip (CI: 2.2–5.4 km). The 1.1 m coseismic slip of F1 obtained by forward modelling is similar to the median slip of 1.2 m predicted by the inverse model (Fig. S11A–F).

In order to evaluate a possible dependency between the slip of the F1 and F2 faults and to rule out a trade-off, we tested their fault-slip variability as independent variables with respect to the NRMSE. We observed variabilities of ~400% and ~10% for F1 and F2 as independent variables, respectively, suggesting a strong dependency of F2 with respect to slip values of F1 lower than 1.1 m (Fig. S12). However, our best-fit model suggests that the slip magnitude of F2 is one order of magnitude smaller than that of F1, and therefore we conclude that the F2 trade-off has a minimal effect on our results.

The forward model residuals are <0.05 m in the near-field (within 10 km of the fault), increasing to ~0.1 m in the far-field. (Fig. S8F). The higher far-field residuals are probably associated with post-seismic deformation following the 2010 Maule earthquake, which was not considered in our model. However, because of its deep source (at the megathrust and continental mantle⁴⁸), post-seismic deformation is associated with wavelengths of ~100 km⁴⁸,⁴⁹, which is an order of magnitude greater than deformation related to the shallower PIF (Fig. 4A, B). We therefore consider that our model adequately reproduces the near-field surface deformation resulting from coseismic slip along the PIF. Our model results suggest that the PIF earthquakes occurred along blind faults, whose deformation pattern differ from those associated with surface-breaching faults (see comparison in Fig. S13A, B).

We estimated a long-term slip rate for the PIF (over the past ~500 ka) by forward modelling the spatial pattern of uplift rates estimated from marine terraces. We modelled surface deformation associated with the same two faults used in the coseismic model (F1 and F2; Fig. 4C, D), which clearly offset the elevations of marine terraces. We generated two models, one based on uplift rates estimated from measured shoreline angles, and a second one, based on an interpolated surface derived from these measurements to increase the area available for comparing measurements with model results (Fig. 4D, see Methods Section, Analysis of off-fault geomorphic features). In the first model, we
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Coseismic and long-term surface deformation associated with the PIF. A Map of the six levels of marine terraces in the PIF area and ages assigned to Marine Isotope Stages (MIS). Boxes show the location of swath profiles used to map shoreline angles (colour-coded points). Thick black lines denote the two PIF branches (F1 and F2). B Shoreline angles projected along profile P-Q; arrows indicate relative displacements of F1 and F2. C Uplift rates deduced from shoreline angles along profile P-Q, grey lines are standard error estimates of uplift rates (se calculation details in Methods Section Analysis of off-fault geomorphic features). D Coseismic Line-of-Sight (LOS) displacements in the Pichilemu area determined from Envisat® radar interferometry of the 2010 PIF earthquakes. Red and orange arrows represent horizontal and vertical displacements, respectively, estimated at the continuous GPS station PICH (red dot) as a result of the 2010 PIF earthquakes³⁴.</figcaption></figure>

obtained a slip rate of 0.48 m/ka for F1 (CI: 0.46–0.54 m/ka), for a fault extending from 26.2 to 1.2 km depth (CI: 0.8–2.6 km) and 0.2 m/ka for F2 (CI: 0.08–0.22 m/ka) with an up-dip depth of 3.6 km (CI: 1.2–3.8 km). The second model yielded a slip rate of 0.52 m/ka (CI: 0.48–0.56 m/ka) for F1 and 0.1 m/ka for F2 (CI: 0.08–0.2 m/ka; Figs. S9B and S10E–H). The best-fit up-dip slip depths were 1.8 km for F1 (CI: 1–2.6 km) and 1.2 km for F2 (CI: 0.8–3.4 km F1; Figs. S9C and S10I–L). Importantly, both models predict similar deformation patterns and slip rates, and these differed from those predicted for a surface-breaching fault (Figs. S13C, D). We selected the second model as the best-fit that optimally reproduces the deformation pattern of marine

terraces, with lower uncertainties than the first model and convergence towards low NRMSE values (Fig. 4D and Fig. S9C).

Marine terraces in the Maule region record two wavelengths of deformation; long wavelength deformation patterns (>100 km) associated with deformation along the coast controlled by the megathrust, and short-wavelength patterns associated with crustal faults³⁰. The pattern of deformed marine terraces in the PIF area exhibits a short deformation wavelength of ~10 km, which is one order of magnitude smaller than the deformation related to the megathrust. The effect of megathrust deformation is therefore negligible compared to the shorter wavelength of crustal faults.
<table><caption>Table 1 Post-IR IRSL samples.</caption><thead><tr><th>Sample</th><th>Long (deg)</th><th>Lat (deg)</th><th>Z (m)</th><th>Depth (m)</th><th>Paleodose (Gy)</th><th>No. aliquots</th><th>Over dispersion (%)</th><th>U (ppm)</th><th>Th (ppm)</th><th>K (%)</th><th>Water (%)</th><th>Dose rate (Gy/ka)</th><th>Age (ka)</th></tr></thead><tbody><tr><td>PICH-1</td><td>-71.97</td><td>-34.38</td><td>51</td><td>2.80</td><td>242 ± 13</td><td>5</td><td>2</td><td>1.44 ± 0.07</td><td>4.94 ± 0.28</td><td>0.86 ± 0.03</td><td>5 ± 5</td><td>2.3 ± 0.1</td><td>106 ± 9.3</td></tr><tr><td>LOBO-1</td><td>-72.04</td><td>-34.42</td><td>14</td><td>3.20</td><td>748 ± 58</td><td>14</td><td>22</td><td>0.96 ± 0.05</td><td>3.61 ± 0.22</td><td>1.09 ± 0.04</td><td>12 ± 5</td><td>2.3 ± 0.1</td><td>328 ± 33</td></tr><tr><td>PICH-2</td><td>-71.96</td><td>-34.33</td><td>115</td><td>2.20</td><td>706 ± 51</td><td>19</td><td>22</td><td>1.13 ± 0.07</td><td>5.32 ± 0.34</td><td>1.05 ± 0.01</td><td>6 ± 5</td><td>2.4 ± 0.1</td><td>297 ± 29</td></tr><tr><td>PICH-4</td><td>-72.03</td><td>-34.42</td><td>38</td><td>6.5</td><td>784 ± 54</td><td>17</td><td>18</td><td>1.18 ± 0.07</td><td>5.28 ± 0.29</td><td>1.20 ± 0.02</td><td>9 ± 5</td><td>2.5 ± 0.1</td><td>314 ± 30</td></tr><tr><td>SM mb</td><td>-73.51</td><td>-37.04</td><td>2</td><td>0.5</td><td>3.4 ± 0.1</td><td>15</td><td>9</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Radionuclide analyses include Uranium (U), Thorium (Th) and Potassium (K). Sample SM mb is a present-day beach berm sample from Santa Maria Island used to evaluate the completeness of signal resetting.

<figure><img src="image_7.png" /></figure>

**Fig. 4 Coseismic and long-term fault-slip models.** A Best-fit model forward of Line of Sight (LOS) displacements as a result of the 2010 PIF earthquakes, green and yellow arrows denote horizontal and vertical displacements, respectively, estimated at continuous GPS station PICH³⁴. B Swath profile showing Envisat* LOS displacements and best-fit forward model. Red dot shows PICH displacement projected to Envisat LOS. Insets in B are the probability-density distributions of coseismic slip of the forward (black line) and inverse (red line) models (See Figs. S10 and S11), blue lines denote best-fit slip and median slip of both models, respectively. C Best-fit dislocation model based on interpolated marine terrace uplift rates (yellow dots denote measurement locations). D Swath profile comparing fault-slip model based on shoreline angles (red line) with a model based on interpolated uplift rates from shoreline angles (black line). Insets show probability-density distributions of slip rate values used to estimate confidence interval (See Fig. S10). The green polygons in D are the interpolated uplift rate surfaces projected along the profile.

Our coseismic and long-term slip models both suggest blind faulting with similar surface deformation patterns. The differences in up-dip depths between both models may be associated with a partial rupture during the PIF earthquakes, as has been observed elsewhere in other crustal earthquakes e.g., refs. ⁵⁰,⁵¹, or to the relatively simple model setups that assume heterogeneous slip along a planar fault with homogeneous rheology. The inverse model predicts a median slip of 1.2 m arranged into an irregular patch of slip centred at ~12.5 km depth that extends upwards to ~6 km depth (Fig. S11D). This is consistent with blind faulting, which is also indicated by both the coseismic and long-term forward models. Taken together, our models closely reproduce the pattern and magnitude of observed surface displacements, suggesting that deformation associated with the PIF occurs over a ~10 km-wide area.

**Earthquake recurrence of the Pichilemu fault.** Both the coseismic and long-term surface deformation patterns of the PIF are similar, and both are consistent with extensional fault

kinematics. We therefore propose that the PIF accrues permanent deformation only during slip triggered by megathrust earthquakes, such as during the events observed after the Maule earthquake. Historical, paleoseismic, and paleo-tsunami records suggest a recurrence time of ~0.1–0.2 ka for Maule-type events⁵²⁻⁵⁴. If we consider a recurrence time of 0.2 ka, then a triggered slip of 0.1 m would account for the long-term PIF slip rate (Fig. 5A). The inferred slip would be equivalent to offsets during a Mw 5.3–6.1 earthquake based on empirical relationships⁵⁵⁻⁵⁷ (i.e., a seismic moment between 1.1 × 10¹⁷ and 1.8 × 10¹⁸ Nm, Fig. 5B). In turn, a recurrence time of 0.5 ka would imply 0.26 m of PIF slip (Fig. 5A) and an earthquake magnitude between Mw 5.9 and Mw 6.3 (Fig. 5B). In both cases the estimated slip per event would be five to ten times less than during the 2010 PIF earthquakes. Instead, the amount of slip per event required to trigger a Mw 7 earthquake would be 1.65 m; in this case, the recurrence time required to account for the long-term slip rate would be 3.2 ka. By combining the probability distribution of the long-term PIF slip rate with that of the coseismic slip we infer a recurrence time of 2.12 ka (1.85–2.28 ka at 90%
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Recurrence time for earthquakes on the PIF. **A** Recurrence times estimated for different values of coseismic slip values and long-term slip rates, assuming that PIF slip occurs only during megathrust earthquakes. **B** Equivalent seismic moment for each recurrence time in **A**, indicated by colour-coded points. Lines denote different empirical relationships (see references in text). The black-to-grey areas in **A**, **B** represent the probability-density of the PIF slip. The dashed red lines in **A**, **B** show the 90% confidence interval. The black line in **A** and white line (B) show the best-fit slip per event and slip rate, respectively.</figcaption></figure>

confidence interval, Fig. 5A) for Mw 6.7 to Mw 6.8 earthquakes (seismic moment between $1.4 \times 10^{19}$ and $2.0 \times 10^{19}$ Nm; Fig. 5B), similar to that of the 2010 PIF earthquakes. Our results suggest that the recurrence time of the PIF may be over an order of magnitude greater than that of Mw > 8 megathrust earthquakes in the Maule segment, implying sporadic triggering of fault slip.

The forearc along the 2010 Maule rupture zone includes ten normal faults similar to the PIF⁵⁸ and associated with similar throws³⁰, slip rates²⁴,⁵⁹ and lengths of their surface traces. If all of these faults would be associated with a similar ~2 ka recurrence time to that of the PIF, then extensional slip on each fault could be sequentially triggered by the megathrust earthquakes that occur approximately every ~0.2 ka in this region. However, triggering would depend on the relationship between the normal fault geometry and the locus of the megathrust slip, as discussed below.

**Megathrust earthquakes and Pichilemu fault slip.** Megathrust earthquakes induce an instantaneous reversal in the polarity of the stress field in the upper plate. This change may induce the seismic rupture of crustal faults that are optimally oriented with respect to the new stress field e.g., refs. ⁶⁰,⁶¹. To explore possible rupture scenarios, we modelled the Coulomb failure stress (ΔCFS) induced by different slip distributions during megathrust earthquakes and considering the PIF as a receiver normal fault. The Maule earthquake involved slip within two subsegments at the northern and southern parts of the rupture³¹, and induced positive ΔCFS values of 1.3–6 MPa in the PIF area (Fig. 6A). Model results based on geophysical and paleo-seismological observations along convergent margins suggest that the heterogeneous frictional properties of the subduction megathrust may result in temporally variable slip patterns e.g., refs. ⁶². It is therefore possible that these two Maule subsegments may be associated with different slip behaviour and stress transfer to the PIF. For this reason, we generated synthetic megathrust-slip distributions at different locations along the Maule rupture zone (Fig. 6A–D and Fig. S14A–D). Slip along the southern Maule subsegment would result in a ΔCFS of −0.2 MPa at the PIF, inhibiting slip (Fig. 6D). Conversely, megathrust slip along the northern subsegment would result in a ΔCFS of 0.66 MPa in the

PIF area (Fig. 6B), promoting slip (Fig. 6B). We noticed that ΔCFS values on the PIF increase progressively as the locus of synthetic slip shifts northward, reaching a maximum of 0.23–0.34 MPa for friction coefficients of between 0.4 and 0.75 (Fig. S14A, B).

It has been suggested that the down-dip segmentation of megathrusts may control the magnitude and characteristics of subduction earthquakes and the activity of crustal faults⁶³,⁶⁴. We compared the ΔCFS induced by earthquakes at variable megathrust down-dip depths by simulating synthetic slip scenarios (Fig. S14C). We observed that when the locus of slip is in the deeper part of the megathrust (between 40 and 50 km depth), ΔCFS values in the PIF area are negative (−2 MPa), thus inhibiting slip (Fig. 6C). In contrast, when the locus of slip occurs at shallower depths or directly below the PIF (20–30 km depth), ΔCFS values on the PIF are above 0.2 MPa, and thus promote slip (Fig. S14D). However, megathrust slip at depths of less than ~15 km reduced the ΔCFS values below 0.2 MPa (Fig. S14D), indicating a ΔCFS that is still positive, but less favourable for triggering slip along PIF.

We explore the potential triggering of PIF during the interseismic phase using the degree of plate locking estimated from GPS measurements obtained over the decade preceding the Maule earthquake³¹. The interseismic phase induced negative ΔCFS values at the PIF (−1.44 MPa, Fig. 6E) inhibiting normal fault slip. This is consistent with the lack of extensional crustal earthquakes in the area during the 25 years prior to the Maule earthquake³⁴, suggesting that normal faulting along the PIF is unlikely during the interseismic phase. Farías et al.⁴ suggested that Mesozoic N-S-striking structures, such as the Vichuquén Fault (Fig. 1A), are favourably oriented for reverse reactivation during interseismic contraction and could contribute to the long-term build-up of topography. However, these N-S-striking faults are not included in the Chilean Database of Active Faults⁵⁸ and we found no field evidence supporting their potential Quaternary activity. Further studies are clearly needed to assess the influence of interseismic contraction in the reactivation of upper-plate faults.

The Coulomb stress models show that different slip distributions during megathrust earthquakes can either promote or inhibit slip along the PIF. These results may explain why the
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Coulomb stress changes and PIF slip. Mean ΔCFS maps calculated for a depth range between 5 and 25 km using the fault parameters of the PIF and a friction coefficient of 0.4. **A** Mean ΔCFS map derived from the coseismic slip of the Maule earthquake³¹ on the PIF. **B–D** Mean ΔCFS produced by synthetic slip scenarios and the same fault parameters as used in **A**. **E** Mean ΔCFS calculated using the degree of interseismic plate locking estimated from GPS measurements obtained over the decade preceding the Maule earthquake³¹. Note that positive ΔCFS that promote normal slip on the PIF occur only in scenarios depicted in **A**, **B**, **F** and **G** Profiles across the PIF showing ΔCFS during the 2010 Maule earthquake and during the pre-Maule interseismic period. White contours in **A–D** denote coseismic slip magnitude and the degree of plate locking in **E**.</figcaption></figure>

estimated 2.1 ka recurrence time for large PIF earthquakes exceeds the recurrence time of megathrust earthquakes in this region. We propose two possible scenarios: First, if the Maule segment was characterised by ruptures with different lengths, magnitudes, and slip distributions, then the PIF could have been triggered by only one in 10–12 megathrust earthquakes, and only when megathrust slip involved the northern Maule sub-segment at shallow depths (20–30 km depth). Second, Mw > 8 historical megathrust earthquakes such as the 1835 and 1751 events suggest that the Maule segment was characterised by similar earthquake ruptures⁶⁵. Previous studies of crustal faults in northern Chile have proposed that build-up of tensional stresses after successive seismic cycles may be counterbalanced by negative ΔCFS during the interseismic phase⁶¹, extending the recurrence time of crustal faults and providing an explanation for the long recurrence times that we estimated for the PIF earthquakes.

**Origin of the PIF and implications for blind faults.** Blind faults usually occur in sedimentary basins when fault-tip propagation fails to reach the surface across a thick sedimentary cover⁶⁶. This is, however, not the case for the PIF, which mostly affects crystalline basement rocks. Blind structures can also be controlled by mechanical and rheological heterogeneities in the upper crust⁶⁶. The PIF host rocks are part of a Palaeozoic accretionary wedge characterised by pervasively deformed high-pressure and high-temperature metasedimentary and metavolcanics rocks (Fig. S1A) that locally form melanges²⁷⁻²⁹,³⁸. These rocks have been overprinted by brittle deformation during Cenozoic exhumation, resulting in a rheologically heterogeneous and fragmented upper

crust²⁸. We propose that this rheological heterogeneity may have preconditioned blind faulting by favouring strain diffusion across a broad zone of deformation within the uppermost crustal levels, thus preventing localised surface faulting. This assessment is further supported by the distribution of aftershock seismicity following the 2010 earthquakes across a ~10 km-wide zone (Figs. 1C and S7), and by low P-wave velocity inferred from seismic tomography, which defines an area characterised by intense fracturing⁴⁶. We suggest that the lineaments observed in the LiDAR topography (Fig. 2) represent secondary fault-tip bending structures inherited from old crustal fabrics that might control their oblique orientation with respect to the PIF, but that are now either inactive or characterised by very low displacement rates.

The possibility that heterogeneities in upper-crustal rheology either favour or hinder surface-breaching faults is further highlighted by the characteristics of the El Yolki Fault, located 125 km farther south, which shares similar kinematics and orientation with the PIF⁵⁹. In contrast to the PIF, the El Yolki Fault breaches the surface and is surrounded by metamorphic basement rocks consisting exclusively of moderately deformed metapelites displaying subhorizontal foliation²⁸. This more homogeneous rheological character may therefore have promoted localized strain along the fault zone and towards the surface.

Unlike the megathrust, crustal-scale faults in the upper plate of subduction zones can remain inactive for long periods and in most cases not detectable by geodetic or seismic monitoring⁶⁷. Furthermore, as some of these faults may be blind and lack clear geomorphic expression (such as the PIF), their unambiguous identification by field observations or remote-sensing methods is
challenging. Many of these cryptic crustal faults at subduction margins were only recognised following their seismic reactivation, which raises the question of how many other active faults with the potential to produce large earthquakes may exist in these environments. For example, a statistical analysis of crustal earthquakes in New Zealand has shown that more than half of the historical large earthquakes ruptured along previously unidentified faults¹. Further cryptic and partly hidden active faults have also been documented at the Cascadia⁶⁸,⁶⁹ margin as well as in the Kanto region of Japan⁷⁰, and similar conditions are likely to exist along other subduction zones. We conclude that regional-scale acquisition of LiDAR topography, combined with the quantitative analysis of off-fault and on-fault geomorphic markers and geodetic data is a promising way forward that will ultimately improve our knowledge of the location, deformation mechanisms, and recurrence time of hitherto unidentified seismogenic faults at convergent margins. The approach proposed in this study to detect blind faults and recognise their potential as cryptic seismogenic sources will contribute to the development of more accurate active-fault maps and hazard assessments, with far-reaching implications for seismic risk management in coastal areas.

## Methods

**Estimating coseismic slip during the PIF earthquakes.** We processed two Envisat® radar scenes acquired 2 days before and 7 days after the March 11 2010 earthquake doublet, obtaining line-of-sight (LOS) displacements (Fig. S8A–F, see details on the radar images and interferogram in Table S2). Because the interferometric coherence is relatively low (Fig. S8B), we applied filtering algorithms based on a moving window using an adaptive range filter⁷¹ to facilitate the unwrapping (Fig. S8C). The interferograms were processed using Roi-Pac software⁷² and unwrapped using the branch-cut algorithms. To corroborate our results, we compared the INSAR LOS displacement with those estimated using the permanent GPS station PICH projected to the LOS vector (location in Fig. 3D). To estimate coseismic slip along the PIF, we searched for the parameters that best reproduced the distribution of LOS displacements using forward elastic dislocation modelling⁷³. The elastic dislocation models were programmed in Matlab® using the function okada⁸⁵⁶⁷⁴. We inferred fault geometries from alignments in the cluster of crustal aftershocks associated with the PIF earthquakes (55° dip, N38°W strike and 26 km length for F1 and 72° dip, N16°W strike and 16 km length for F2; Fig. S7). The strike inferred from aftershocks is similar to those from focal mechanisms³⁴ and interpretations based on aftershock distributions³⁴,⁴⁶,⁷⁵. Aftershocks extended continuously down-dip from between ~2 and ~6 km depth to the megathrust at ~26 km (Fig. S7); however, the up-dip limit is rather diffuse. We used forward modelling to estimate the slip and up-dip depth of each fault for defined ranges and constant increments (Table S1). In addition, we included a hypothetical up-dip scenario of surface-breaching faults (Fig. S13). We set the range of slip values from 0.6 to 1.5 m for F1 and 0 to 0.65 m for F2, respectively, based on previous fault-slip inversions³⁴. Using these ranges, we generated 85,918 elastic models (Table S1 and Fig. S9A) and searched for the parameters that minimised the NRMSE.

The root mean squared error (RMSE) is defined as the difference between observation (yi) and model (y), with n being the number of observations (Eq. 1).

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n} (y_i - y)^2}{n}} \quad (1)$$

The NRMSE facilitates comparing models with different scales by normalising the RMSE (Eq. 2), where y max-y min is the range of observations. The uncertainties in model results were estimated using the lower 5% tail of the NRMSE distribution (Fig. S10). We estimated slip values using this tail distribution and defined the model uncertainty as the interval between the 5 and 95% percentiles (Fig. S10).

$$NRMSE = \frac{RMSE}{(ymax - ymin)} \quad (2)$$

To evaluate the consistency of our coseismic slip estimates we compared results from our forward model with those of an inverse model for a heterogeneous slip distribution. The inverse model was based on the automated fault model discretization method of Barnhart and Lohman⁷⁶, which can resolve scattered surface displacement observations by varying the model resolution with depth. We resampled the LOS displacements based on a Delaunay triangulation algorithm allowing to generate a dataset tractable for fault-slip inversions⁷⁶ (Fig. S11A–C). We inverted the resampled interferogram using the geometry of fault F1, a rake of ~90, a fault length of 57 km to avoid boundary artefacts, and a depth of 26 km to the intersection with the megathrust. The fault discretization method is based on

an iterative approach that resizes the triangular fault-slip patches according to the model resolution. This generates smaller patches near the surface and larger patches at depth and offshore. The resulting model geometry of F1 comprised 678 patches (Fig. S11D). The green function matrix calculated for the patches may lead to unstable solutions that can be resolved using regularisation. The regularisation parameter (lambda) is dependent of the data noise and the Laplacian smoothing. The most suitable regularisation represents the correct balance between smoothing and data noise to fit the underlying noise-free signal; while high or low regularisation values result in increasingly smooth or increasingly complex models that do not represent a good fit to the noise-free signal. We selected the best fit slip distribution using the L-curve method⁷⁷ that compares graphically the squared root of model and data misfits, and the optimum lambda is then determined at the point of maximum curvature (Fig. S11E). We additionally checked the jRi value as a function of the lambda and confirmed that the optimum lambda was selected in the jRi-low range⁷⁶ (Fig. S11F). The jRi value is a metric of the quality of the inversion and is related to the input data noise, the fault regularisation and fault parameters, therefore, a lower value of jRi represents a more suitable regularisation value and smoothing⁷⁶. We finally analysed the distributions of slip and model misfits statistically by re-discretizing the triangular meshes into a grid of equal area cells. The contribution of each triangle on each cell was calculated using a triangle area weighted average, this procedure allowed reducing the spatial bias produced by triangles of different sizes.

**Analysis of on-fault geomorphic features.** We performed a detailed geomorphic and morphometric analysis of off-fault and on-fault geomorphic features in the PIF area. On-fault features were identified by analysing topographic and fluvial metrics using the Topobox-2 software³⁵. We created red-relief maps (RRMs)⁷⁸ that combine the terrain openness with surface slope. RRM are useful for identifying lineaments and fault scarps together with changes in the fluvial network, because they lack the potential bias of light source direction common in shaded-relief maps (Fig. S1B). We also generated local relief maps (Fig. S2C) using a 500 m roving window, and slope maps (Fig. S2D) calculated using the 8-connected neighbourhood gradient algorithm of Topobox®⁸³.

Drainage and catchment morphology have the potential to be used as tectonic markers in the quantification of regional strain, uplift, and tilting. Ideally, catchments should be symmetric about the main trunk stream if they have incised a horizontal surface of uniform lithology under homogeneous climatic conditions⁷⁹. To estimate the catchment asymmetry, we extracted the main catchments, the drainage networks and the main trunk streams in the Pichilemu area using the flow-routing and accumulation algorithms of Topobox®⁸³. We subsequently estimated the catchment centreline by creating a distance buffer from the border of each catchment polygon. Then the distance matrix was skeletonised using plain curvature to obtain the centreline and was compared with the main trunk stream to estimate the catchment symmetry factor (Fig. S2B). The symmetry factor Ts³⁷ (Eq. 3) is used to highlight changes in catchment asymmetry and areas of lateral tilting and is defined as:

$$Ts = Da/Dd \quad (3)$$

where Da represents the distance from the catchment centreline to the main catchment trunk and Dd is the distance from the catchment centreline to the catchment boundary. The Ts values range between 0 and 1, where 0 represents a perfectly symmetric catchment and 1 represents a tilted catchment. Da orientations are presented as vectors that are colour-coded by Ts and incorporated into rose diagrams to interpret the dominant tilting direction (Fig. S2B).

Knickpoints were identified using the knickpointfinder algorithm in Topobox®⁸⁰, which iteratively maps the location of knickpoints by fitting a concave upward profile to the river-profile elevations, with the condition that the fitted curve must run below the profile elevations. The maximum vertical distance between the observed and modelled profile is defined by the tolerance parameter. We used a tolerance of 18.5 m, based on the maximum difference between the upstream maxima and downstream minima of the main trunk river profile⁸⁰.

The steepness index (Ksn) can be useful for obtaining information on tectonic and/or climatic perturbations in a fluvial network. Ksn values allow the quantification of deviations from the steady-state concavity of a river profile and the detection of these deviations along the stream network. The estimation of Ksn values is based on the power law (Eq. 4) of detachment-limited incision into bedrock⁸¹.

$$\frac{dz}{dt} = U - kA^m \left(\frac{dz}{dx}\right)^n \quad (4)$$

Where U is the uplift rate, A is the upslope area, dz/dx is the channel slope and m, k and n are constants. In steady-state conditions dz/dt = 0; hence, we can rearrange the equation as [Eq. 5]:

$$\frac{dz}{dx} = \left(\frac{U}{k}\right)^{\frac{1}{n}} A^{-\frac{m}{n}} \quad (5)$$

were (U/k)¹/² represents the channel steepness and m/n is the channel concavity
(θ); thus, the equation can be written as [Eq. 6]:

$$K_{sn} = \frac{S}{A^{-\theta}} \qquad (6)$$

In contrast to the Ksn, chi-plots are based on the horizontal transform of the upslope area to linearise the concave upward profile for a well-chosen reference concavity. This spatial transform makes chi-plots useful for identifying transient erosional signals with a common origin propagating upstream along the drainage network, such as tectonically generated knickpoints or changes in base level. Furthermore, the chi dimension allows comparing these signals between different catchments, irrespective of their size or shape.

Estimating chi-values requires rearranging Eq. [4] to convert dx to a distance measured from the catchment outlet and assuming that U and k are spatially uniform (Eq. 7).

$$\int \frac{dz}{dx} dx = z(x0) + \left( \frac{U}{kA0^m} \right)^{\frac{1}{m}} \int \frac{A0}{A(x)^m} dx \qquad (7)$$

The integral of dz/dx is calculated using a reference area (A0), where x0 is the catchment outlet. The chi-plot values (χ) are estimated based on the right hand-integral (Eq. 8).

$$\chi = \int \left( \frac{A0}{A(x)} \right)^{\frac{1}{m}} dx \qquad (8)$$

To calculate the Ksn and chi-values in the catchments of the Pichilemu area we used Topobox®³⁵, which includes all the above equations. We used a reference concavity of 0.45.

**Analysis of off-fault geomorphic features.** To analyse off-fault geomorphic features we studied marine terraces using two methods, depending on the type of terraces and their origin, with the aim to identify and map marine shoreline angles. The shoreline-angle is a geomorphic marker located at the intersection between the paleo-platform and paleo-cliff that represent the maximum reach of the sea level during a highstand period that can be used to estimate vertical deformation and coastal uplift rates⁴⁰.

We analysed wave-built marine terraces following the principles outlined in Jara-Muñoz and Melnick⁴³ taking into account the morphology of the bedrock unconformity, the number of sedimentary cycles within the wave-built terrace, and the thickness of the sequence (Fig. S4). We mapped the surface morphology of the wave-built terraces using swath profiles to detect any breaks in slope. We also measured the depth to the crystalline bedrock in incised valleys and generated an isopach map of sedimentary sequence thickness (Fig. S4A, B). This allowed us to differentiate sedimentary-sequences and to improve the estimation of shoreline-angle elevations.

We studied the surface morphology of marine terraces using LiDAR topography and swath profiles in order to measure the locations and elevations of shoreline angles (see Supplementary Data 1). To estimate uplift rates (u), we correlated terrace levels with sea-level highstands using the IRSL ages and a composite sea-level curve for the southern hemisphere spanning the last 700 ka⁴⁴,⁴⁵ (Eq. 9, Fig. S6).

$$u = \frac{(E - e)}{T} \qquad (9)$$

where E is the elevation of shoreline angles, e is the elevation of the corresponding highstand and T the age of the terrace level. Uplift rate errors Se(u)² were estimated following Gallen et al.⁸² as [Eq. 10]:

$$Se(u)^2 = u^2 \left( \left( \frac{\sigma^2 H}{H^2} \right) + \left( \frac{\sigma^2 T}{T^2} \right) \right) \qquad (10)$$

where σH is the error in relative sea-level, defined as [Eq. 11]:

$$\sigma H = \sqrt{\sigma E^2 + \sigma e^2} \qquad (11)$$

where σT is the age uncertainty in the sea-level curve (7 ka), σE is the error of the shoreline-angle assessments and σe is 12 m uncertainty of the highstand elevation based on Rohling et al.⁴⁵.

**Post-IR IRSL dating.** Four sediment samples were analysed by thermoluminescence dating were analysed at the University of Cologne using the post-IR IRSL signal of K-feldspar obtained from marine terrace sediments. The sedimentary units sampled were (1) deposited as shallow marine sediments, ideally in berm or swash-zone environments; (2) comprised of sandy sediments with medium grain sizes and more than 20% of feldspars and (3) from the base of the sequence, as close as possible to the bedrock wave-cut platform (Fig. S3). We analysed medium-sized K-Feldspar sand grains (100–250 µm) following the post-IR IRSL₂⁹⁰ SAR protocol⁸³. The dated sediments were generally characterised by high feldspar signals with adequate reproducibility in dose-recovery tests performed after signal resetting in a solar simulator for 12 h (satisfactory ratios between measured and laboratory dose between 0.9 and 1.1). Burial doses were based on 5–19 aliquots of 8 mm diameter, using the central age model for calculation⁸⁴ (See Abanico-type plots in Fig. S5). In addition, we evaluated the completeness of signal resetting using a sample from a

modern beach berm (Sample SM mb, Table 1), which revealed an insignificant residual dose compared to the burial doses of the marine terrace samples (i.e., ~3.4 Gy). Radionuclide analysis (uranium, thorium and potassium) for dose-rate estimation was carried out using high-resolution gamma-spectrometry (see details in Table 1). We used a potassium content of 12.5 ± 0.5% to estimate internal dose rates, based on the method proposed by Huntley and Baril⁸⁵.

**Slip rate of the PIF from deformed marine terraces.** We estimated the PIF slip rate by searching for the best-fitting input parameters (i.e., the lowest NRMSE value) in a set of 116,964 model runs to find the minimum NRMSE value (Fig. S9B, C, Eq. 2). Uplift rates derived from shoreline-angle elevations of marine terraces were reproduced by forward elastic dislocation models by varying up-dip and down-dip slip depths, as well as the slip rate, of each fault (Table S1). The elastic model setups include the same geometry and down-dip depth used in the coseismic models. We performed preliminary modelling experiments using one and two faults, obtaining a better result using two faults (Fig. 4D). To perform the comparisons, we iterated the slip magnitude and up-dip depths of F1 and F2 using constant increments (Table S1). The best-fitting models are those with the lowest NRMSE value (Fig. S9B, C). To estimate the confidence interval of the best-fit models, we used the lower 5% tail of the NRMSE distribution (Fig. S10E, I). From these distributions, we defined the 5% and 95% percentiles as the confidence limits, equivalent to a 90% confidence interval, as the uncertainty of the best-fit model results (Fig. S10E–L).

Because marine terraces are exposed along the coast, they mostly reflect deformation along a 2D profile. To increase the area available for comparing dislocation-model results with observations we carried out a natural neighbour interpolation using a Delaunay triangulation of the scattered shoreline angles in 100 m bins (Fig. 4D). Because the spatial distribution of marine terraces is not suitable for estimating fault slip with a 3D inverse model, we used forward dislocation models in order to produce comparable results at both coseismic and long-term timescales.

**Modelling Coulomb stress failure.** We used the ΔCFS to evaluate the potential for a megathrust earthquake to induce slip along the PIF. We modelled the ΔCFS in the PIF area and directly on the PIF using the fault geometry inferred from aftershocks (Fig. S7) and the Coulomb 3.4 algorithms⁸⁶. The model setup and rheological parameters were based on previous studies of the South American margin and in the Maule region e.g., refs. 4,3,4,6¹, we used a Poisson’s ratio of 0.25 and Young modulus of 75 GPa both based on average values for upper-crustal materials⁸⁶. The friction coefficient of crustal faults commonly varies over a wide range⁶¹, and our analysis therefore considered a range of friction coefficients between 0.4 and 0.75.

Coulomb stresses imparted by slip on a source fault can either promote or inhibit slip along a receiver fault⁶⁰. For this analysis we modelled the slip imparted during the 2010 Maule earthquake and synthetic slip scenarios in elliptical patches of 100 km length and 50 km width (Fig. S14A, C). We set the area and amount of slip of these synthetic slip patches based on the average slip distribution of the 2010 Maule earthquake. In a first set of modelling experiments, the locus of 28 synthetic slip patches was shifted along the Maule rupture zone between 38° and 35°S and within its up-dip and down-dip depth limits³¹. In a second set of experiments we shifted the locus of slip of 15 synthetic patches to different down-dip depths values of between 10 and 50 km, beyond the down-dip depth limit of aftershock seismicity³³. We also modelled the potential effect of continental contraction induced by interseismic plate-locking on the PIF using a back-slip model based on the pre-2010 Maule earthquake locking distribution of Moreno et al.³¹. We considered a plate convergence rate of 66 mm/yr²⁶ and assumed complete elastic strain release during the 1835 earthquake, which was similar in rupture area and magnitude to the 2010 Maule earthquake³¹.

## Data availability

The data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Information and Supplementary Data 1. Additional data related to this paper such as the LiDAR topography may be requested from the authors.

## Code availability

The Matlab® codes used for marine terrace mapping and dislocation modelling are part of TerraceM available at www.terracem.com.

Received: 30 September 2020; Accepted: 6 May 2022;

Published online: 11 June 2022

## References

1. Nicol, A., Van Dissen, R., Stirling, M. & Gerstenberger, M. Completeness of the paleoseismic active-fault record in New Zealand. *Seismological Res. Lett.* 87, 1299–1310 (2016).
2. Barnhart, W. D., Hayes, G. P. & Gold, R. D. The July 2019 Ridgecrest, California, earthquake sequence: Kinematics of slip and stressing in cross-fault ruptures. *Geophys. Res. Lett.* **46**, 11859–11867 (2019).

3. Beavan, J., Fielding, E., Motagh, M., Samsonov, S. & Donnelly, N. Fault location and slip distribution of the 22 February 2011 Mw 6.2 Christchurch, New Zealand, earthquake from geodetic data. *Seismological Res. Lett.* **82**, 789–799 (2011).

4. Farias, M., Comte, D., Roecker, S., Carrizo, D. & Pardo, M. Crustal extensional faulting triggered by the 2010 Chilean earthquake: The Pichilemu Seismic Sequence. *Tectonics* **30**, TC6010 (2011).

5. Hollingsworth, J., Ye, L. & Avouac, J. P. Dynamically triggered slip on a splay fault in the Mw 7.8, 2016 Kaikoura (New Zealand) earthquake. *Geophys. Res. Lett.* **44**, 3517–3525 (2017).

6. Quigley, M. et al. Surface rupture during the 2010 Mw 7.1 Darfield (Canterbury) earthquake: Implications for fault rupture dynamics and seismic-hazard analysis. *Geology* **40**, 55–58 (2012).

7. Kumar, S. et al. Paleoseismic evidence of great surface rupture earthquakes along the Indian Himalaya. *J. Geophys. Res.: Solid Earth* **111**, B03304 (2006).

8. Sherrod, B. L. et al. Holocene fault scarps near Tacoma, Washington, USA. *Geology* **32**, 9–12 (2004).

9. Burrato, P., Vannoli, P., Fracassi, U., Basili, R. & Valensise, G. Is blind faulting truly invisible? Tectonic-controlled drainage evolution in the epicentral area of the May 2012, Emilia-Romagna earthquake sequence (northern Italy). *Ann. Geophys.* **55**, 525–531 (2012).

10. Talebian, M. et al. The 2003 Bam (Iran) earthquake: Rupture of a blind strike-slip fault. *Geophys. Res. Lett.* **31**, L11611 (2004).

11. Watterson, J. Fault dimensions, displacements and growth. *Pure Appl. Geophys.* **124**, 365–373 (1986).

12. Baudon, C. & Cartwright, J. A. 3D seismic characterisation of an array of blind normal faults in the Levant Basin, Eastern Mediterranean. *J. Struct. Geol.* **30**, 746–760 (2008).

13. Melnick, D., Bookhagen, B., Strecker, M. R. & Echtler, H. P. Segmentation of megathrust rupture zones from fore-arc deformation patterns over hundreds to millions of years, Arauco peninsula, Chile. *J. Geophys. Res.: Solid Earth* **114**, B01407 (2009).

14. Shaw, J. H. & Shearer, P. M. J. S. An elusive blind-thrust fault beneath metropolitan Los Angeles. *Science* **283**, 1516–1518 (1999).

15. Yeats, R. S. Large-scale Quaternary detachments in Ventura basin, southern California. *J. Geophys. Res.: Solid Earth* **88**, 569–583 (1983).

16. Savage, J. C. A dislocation model of strain accumulation and release at a subduction zone. *J. Geophys. Res.: Solid Earth* **88**, 4984–4996 (1983).

17. Becker, T. W., Hashima, A., Freed, A. M. & Sato, H. Stress change before and after the 2011 M9 Tohoku-oki earthquake. *Earth Planet. Sci. Lett.* **504**, 174–184 (2018).

18. Toda, S. & Sutsumi, H. Simultaneous reactivation of two, subparallel, inland normal faults during the M w 6.6 11 April 2011 Iwaki earthquake triggered by the M w 9.0 Tohoku-oki, Japan, earthquake. *Bull. Seismological Soc. Am.* **103**, 1584–1602 (2011).

19. Melnick, D., Moreno, M., Motagh, M., Cisternas, M. & Wesson, R. L. Splay fault slip during the Mw 8.8 2010 Maule Chile earthquake. *Geology* **40**, 251–254 (2012).

20. Plafker, G. Surface faults on Montague Island associated with the 1964 Alaska earthquake. U.S. Geological Survey Professional Paper. **543-G**, 41 (1967).

21. Audin, L., Lacan, P., Tavera, H. & Bondoux, F. Upper plate deformation and seismic barrier in front of Nazca subduction zone: The Chololo Fault System and active tectonics along the Coastal Cordillera, southern Peru. *Tectonophysics* **459**, 174–185 (2008).

22. Barnes, P. M., Nicol, A. & Harrison, T. Late Cenozoic evolution and earthquake potential of an active listric thrust complex above the Hikurangi subduction zone, New Zealand. *Geol. Soc. Am. Bull.* **114**, 1379–1405 (2002).

23. Singh, S., Hananto, N. & Chauhan, A. Enhanced reflectivity of aftershocks in the recent great Sumatran earthquake rupture zones. *Geophys. Res. Lett.* **38**, L04302 (2011).

24. Jara-Muñoz, J. et al. Quantifying offshore fore-arc deformation and splay-fault slip using drowned Pleistocene shorelines, Arauco Bay, Chile. *J. Geophys. Res.: Solid Earth* **122**, 4529–4558 (2017).

25. Wendt, J., Oglesby, D. D. & Geist, E. L. Tsunamis and splay fault dynamics. *Geophys. Res. Lett.* **36**, L15303 (2009).

26. Angermann, D., Klotz, J. & Reigber, C. Space-geodetic estimation of the Nazca-South America Euler vector. *Earth Planet. Sci. Lett.* **171**, 329–334 (1999).

27. Hyppolito, T. et al. LP/HT metamorphism as a temporal marker of change of deformation style within the Late Palaeozoic accretionary wedge of central Chile. *J. Metamorphic Geol.* **33**, 1003–1024 (2015).

28. Willner, A. P., Richter, P. P. & Ring, U. Structural overprint of a late Paleozoic accretionary system in north-central Chile (34-35 S) during post-accretional deformation. *Andean Geol.* **36**, 17–36 (2009).

29. Willner, A. P. et al. Time markers for the evolution and exhumation history of a Late Palaeozoic paired metamorphic belt in North-Central Chile (34°-35° 30'S). *J. Petrol.* **46**, 1835–1858 (2005).

30. Jara-Muñoz, J., Melnick, D., Brill, D. & Strecker, M. R. Segmentation of the 2010 Maule Chile earthquake rupture from a joint analysis of uplifted marine terraces and seismic-cycle deformation patterns. *Quat. Sci. Rev.* **113**, 171–192 (2015).

31. Moreno, M. et al. Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake. *Earth Planet. Sci. Lett.* **321–322**, 152–165 (2012).

32. Lavenu, A. & Encinas, A. Deformación frágil de los depósitos noógenos de la cuenca de Navidad (Cordillera de la Costa, 34 S, Chile central). *Rev. geológica de Chile* **32**, 229–248 (2005).

33. Lange, D. et al. Aftershock seismicity of the 27 February 2010 Mw 8.8 Maule earthquake rupture zone. *Earth Planet. Sci. Lett.* **317–318**, 413–425 (2012).

34. Ryder, I. et al. Large extensional aftershocks in the continental forearc triggered by the 2010 Maule earthquake, Chile. *Geophys. J. Int.* **188**, 879–890 (2012).

35. Schwanghart, W. & Scherler, D. TopoToolbox 2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences. *Earth Surf. Dyn.* **2**, 1–7 (2014).

36. Bookhagen, B. High resolution spatiotemporal distribution of rainfallfall seasonality and extreme events based on a 12-year TRMM time series. (UC Santa Barbara Geography, 2013).

37. Cox, R. T. Analysis of drainage-basin symmetry as a rapid technique to identify areas of possible Quaternary tilt-block tectonics: An example from the Mississippi Embayment. *Geol. Soc. Am. Bull.* **106**, 571–581 (1994).

38. Kato, T. T. & Godoy, E. Middle to late Triassic mélange exhumation along a pre-Andean transpressional fault system: coastal Chile (26–42 S). *Int. Geol. Rev.* **57**, 606–628 (2015).

39. Perron, J. T. & Royden, L. An integral approach to bedrock river profile analysis. *Earth Surf. Process. Landf.* **38**, 570–576 (2013).

40. Lajoie, K. R. in *Active Tectonics* (ed. Wallace, R.) 95–124 (National Academy Press, 1986).

41. Freisleben, R., Jara-Muñoz, J., Melnick, D., Martínez, J. M. & Strecker, M. R. Marine terraces of the last interglacial period along the Pacific coast of South America (1°N–40°S). *Earth Syst. Sci. Data* **13**, 2487–2513 (2021).

42. Jara-Muñoz, J., Melnick, D., Pedoja, K. & Strecker, M. R. TerraceM-2: A Matlab® Interface for Mapping and Modeling Marine and Lacustrine Terraces. *Front. Earth Sci.* **7**, 1–18 (2019).

43. Jara-Muñoz, J. & Melnick, D. Unraveling sea-level variations and tectonic uplift in wave-built marine terraces, Santa María Island, Chile. *Quat. Res.* **83**, 216–228 (2015).

44. Bintanja, R., van de Wal, R. S. & Oerlemans, J. Modelled atmospheric temperatures and global sea levels over the past million years. *Nature* **437**, 125 (2005).

45. Rohling, E. et al. Antarctic temperature and global sea level closely coupled over the past five glacial cycles. *Nat. Geosci.* **2**, 500–504 (2009).

46. Calle-Gardella, D., Comte, D., Farias, M., Roecker, S. & Rietbrock, A. Three-dimensional local earthquake tomography of pre-Cenozoic structures in the coastal margin of central Chile: Pichilemu fault system. *J. Seismol.* **25**, 521–533 (2021).

47. Hayes, G. P. et al. Slab2, a comprehensive subduction zone geometry model. *Science* **362**, 58–61 (2018).

48. Weiss, J. R. et al. Illuminating subduction zone rheological properties in the wake of a giant earthquake. *Sci. Adv.* **5**, eaax6720 (2019).

49. Bedford, J. et al. A high-resolution, time-variable afterslip model for the 2010 Maule Mw = 8.8, Chile megathrust earthquake. *Earth Planet. Sci. Lett.* **383**, 26–36 (2013).

50. Elliott, J. et al. Depth segmentation of the seismogenic continental crust: The 2008 and 2009 Qaidam earthquakes. *Geophys. Res. Lett.* **38**, L06305 (2011).

51. Iezzi, F., Roberts, G., Walker, J. F. & Papanikolaou, I. Occurrence of partial and total coseismic ruptures of segmented normal fault systems: Insights from the Central Apennines, Italy. *J. Struct. Geol.* **126**, 83–99 (2019).

52. Bookhagen, B., Echtler, H. P., Melnick, D., Strecker, M. R. & Spencer, J. Q. G. Using uplifted Holocene beach berms for paleoseismic analysis on the Santa Maria Island, south-central Chile. *Geophys. Res. Lett.* **33**, L15302 (2006).

53. Lomnitz, C. Major Earthquakes of Chile: A Historical Survey, 1935–1960. *Seismological Res. Lett.* **75**, 368–378 (2004).

54. Dura, T. et al. Subduction zone slip variability during the last millennium, south-central Chile. *Quat. Sci. Rev.* **175**, 112–137 (2017).

55. Papazachos, B., Scordilis, E., Panagiotopoulos, D., Papazachos, C. & Karakaisis, G. Global relations between seismic fault parameters and moment magnitude of earthquakes. *Bull. Geol. Soc. Greece* **36**, 1482–1489 (2004).

56. Richards, P. G. & Aki, K. *Quantitative seismology: theory and methods*. Vol. I and II 700 (Freeman and Co., 1980).

57. Wells, D. L. & Coppersmith, K. J. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bull. Seismological Soc. Am.* **84**, 974–1002 (1994).
58. Maldonado, V., Contreras, M. & Melnick, D. A comprehensive database of active and potentially-active continental faults in Chile at 1: 25,000 scale. *Sci. Data* **8**, 1–13 (2021).

59. Melnick, D. et al. Hidden Holocene slip along the coastal El Yolki fault in Central Chile and its possible link with megathrust earthquakes. *J. Geophys. Res.: Solid Earth* **124**, 7280–7302 (2019).

60. King, G. C., Stein, R. S. & Lin, J. Static stress changes and the triggering of earthquakes. *Bull. Seismological Soc. Am.* **84**, 935–953 (1994).

61. Cortés-Aranda, J., González, G., Rémy, D. & Martinod, J. Normal upper plate fault reactivation in northern Chile and the subduction earthquake cycle: From geological observations and static Coulomb Failure Stress (CFS) change. *Tectonophysics* **639**, 118–131 (2015).

62. Bilek, S. L. & Lay, T. Subduction zone megathrust earthquakes. *Geosphere* **14**, 1468–1500 (2018).

63. Lay, T. et al. Depth-varying rupture properties of subduction zone megathrust faults. *J. Geophys. Res.: Solid Earth* **117**, B04331 (2012).

64. Li, S., Moreno, M., Rosenau, M., Melnick, D. & Oncken, O. Splay fault triggering by great subduction earthquakes inferred from finite element models. *Geophys. Res. Lett.* **41**, 385–391 (2014).

65. Udias, A., Madariaga, R., Buforn, E., Muñoz, D. & Ros, M. The Large Chilean Historical Earthquakes of 1647, 1657, 1730, and 1751 from Contemporary Documents. *Bull. Seismological Soc. Am.* **102**, 1639–1653 (2012).

66. Coleman, A. J., Duffy, O. B. & Jackson, C. A.-L. Growth folds above propagating normal faults. *Earth-Sci. Rev.* **196**, 102885 (2019).

67. Morell, K. D., Regalla, C., Leonard, L. J., Amos, C. & Levson, V. Quaternary rupture of a crustal fault beneath Victoria, British Columbia, Canada. *GSA Today* **27**, 4–10 (2017).

68. Kelsey, H. M., Sherrod, B. L., Blakely, R. J. & Haugerud, R. A. Holocene faulting in the Bellingham forearc basin: Upper-plate deformation at the northern end of the Cascadia subduction zone. *J. Geophys. Res.: Solid Earth* **117**, 1–26 (2012).

69. Liberty, L. M., Hemphill-Haley, M. A. & Madin, I. P. The Portland Hills Fault: uncovering a hidden fault in Portland, Oregon using high-resolution geophysical methods. *Tectonophysics* **368**, 89–103 (2003).

70. Ishiyama, T., Sato, H., Kato, N., Nakayama, T. & Abe, S. Active blind thrusts beneath the Tokyo metropolitan area: Seismic hazards and inversion tectonics. *Geophys. Res. Lett.* **40**, 2608–2612 (2013).

71. López-Quiroz, P., Doin, M.-P., Tupin, F., Briole, P. & Nicolas, J.-M. Time series analysis of Mexico City subsidence constrained by radar interferometry. *J. Appl. Geophys.* **69**, 1–15 (2009).

72. Rosen, P. A., Hensley, S., Peltzer, G. & Simons, M. Updated repeat orbit interferometry package released. *Eos Trans. Am. Geophys. Union* **85**, 47–47 (2004).

73. Okada, Y. Surface deformation due to shear and tensile faults in a half-space. *Bull. Seismological Soc. Am.* **75**, 1135–1154 (1985).

74. Beauducel, F. Matlab/Octave tools for geophysical studies (IPGP, http://www.ipgp.fr/~beaudu/matlab.html, 2014).

75. Ruiz, J. A. et al. Seismological analyses of the 2010 March 11, Pichilemu, Chile M w 7.0 and M w 6.9 coastal intraplate earthquakes. *Geophys. J. Int.* **197**, 414–434 (2014).

76. Barnhart, W. & Lohman, R. Automated fault model discretization for inversions for coseismic slip distributions. *J. Geophys. Res.: Solid Earth* **115**, B10419 (2010).

77. Segall, P. & Harris, R. Earthquake deformation cycle on the San Andreas fault near Parkfield, California. *J. Geophys. Res.: Solid Earth* **92**, 10511–10525 (1987).

78. Chiba, T., Kaneta, S.-I. & Suzuki, Y. Red relief image map: new visualization method for three dimensional data. *Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.* **37**, 1071–1076 (2008).

79. Ramsey, L. A., Walker, R. T. & Jackson, J. Geomorphic constraints on the active tectonics of southern Taiwan. *Geophys. J. Int.* **170**, 1357–1372 (2007).

80. Schwanghart, W. & Scherler, D. Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques. *Earth Surf. Dyn.* **5**, 821–839 (2017).

81. Whipple, K. X. & Tucker, G. E. Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. *J. Geophys. Res.: Solid Earth* **104**, 17661–17674 (1999).

82. Gallen, S. et al. Active simultaneous uplift and margin-normal extension in a forearc high, Crete, Greece. *Earth Planet. Sci. Lett.* **398**, 11–24 (2014).

83. Thiel, C. et al. Luminescence dating of the Stratzing loess profile (Austria)—Testing the potential of an elevated temperature post-IR IRSL protocol. *Quat. Int.* **234**, 23–31 (2011).

84. Galbraith, R. F., Roberts, R. G., Laslett, G., Yoshida, H. & Olley, J. M. Optical dating of single and multiple grains of quartz from jinmium rock shelter, northern australia: part i, experimental design and statistical models. *Archaeometry* **41**, 339–364 (1999).

85. Huntley, D. & Baril, M. The K content of the K-feldspars being measured in optical dating or in thermoluminescence dating. *Anc. TL* **15**, 11–13 (1997).

86. Toda, S., Stein, R. S., Sevilgen, V. & Lin, J. *Coulomb 3.3 graphic-rich deformation and stress-change software for earthquake, tectonic, and volcano research and teaching-user guide*. Report No. 2331-1258, 63 (US Geological Survey, 2011).

## Acknowledgements

This study was supported by the Millennium Scientific Initiative (ICM) of the Chilean government through grant C160025 "Millennium Nucleus CYCLO The Seismic Cycle Along Subduction Zones", Chilean National Fund for Development of Science and Technology (FONDECYT) grants 1181479, 11180509 and 1190258, the ANID PIA Anillo ACT192169, and the MARISCOS (MAule eaRquake: Integration of Seismic Cycle Observations and Structural investigations) project financed by the German Science Foundation (DFG), grant STR373/30-1. J.J.-M. was supported by the DFG "LIFE" project, grant JA2860/1. M.S. was supported by the DFG "ARTE" project, grant STR373/41-1. LiDAR data was provided by Forestal Arauco to the CYCLO project. We would like to thank Kenneth Fisk and Christian Hillemann for help in the field and also Richard Arthur, Alejandro Araneda, Sonia Muñoz Belmar, Andres Tassara and Javier Morin for support before and during the field campaign.

## Author contributions

J.J.-M., D.M., and M.S. designed the study, participated in the fieldwork, analysed the field observations and models. S.L. carried out the inverse modelling, J.C.-A. carried out the Coulomb modelling, A.S. processed and analysed the radar interferometry, D.B. analysed the IRSL samples. J.J.-M. and D.M. wrote the paper. All authors commented on the paper.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-30754-1.

**Correspondence** and requests for materials should be addressed to J. Jara-Muñoz.

**Peer review information** *Nature Communications* thanks Laura Gregory and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022