The cryptic seismic potential of blind faults revealed by off-fault geomorphology, Pichilemu, Chile.

In seismically-active regions, mapping capable faults and estimating their recurrence time is the first step to assess seismic hazards. Fault maps are commonly based on geologic and geomorphic features evident at the surface; however, mapping blind faults and estimating their seismic potential is challenging because on-fault diagnostic features are absent. Here, we study the Pichilemu Fault in coastal Chile, unknown until it generated a M7.0 earthquake in 2010. The lack of evident surface faulting suggests a partly-hidden blind fault. Using off-fault deformed marine terraces, we estimate a slip-rate of 0.42±0.04 m/ka, which when integrated with deformation estimated from satellite geodesy during the 2010 earthquake suggests a 2.5±0.25 ka recurrence time for M6.6-6.9 extensional earthquakes. We propose that extension is associated with stress changes during megathrust earthquakes and accommodated by sporadic slip during upper-plate earthquakes. Our results have implications for assessing the seismic potential of cryptic faults along seismically-active coasts.


Introduction
The unexpected rupture of previously-unmapped faults during recent large-magnitude earthquakes emphasizes a major lacuna in our knowledge of the location and seismic potential of tectonically-active structures 1  with past earthquakes [7][8][9] . Nevertheless, such evidence may be completely absent or difficult to identify in areas affected by active blind faults e.g. 10,11 . Blind faults are geological structures whose ruptures do not reach the surface 12 , and therefore may have a hidden seismic potential. Such structures are common in sedimentary basins and have been identified from geophysical images e.g. 13,14,15 and indirect geomorphic observations 16,17 . However, estimating the seismic potential of blind faults is difficult using typical geologic and geomorphic evidence of surface ruptures during past earthquakes. Here we show that quantifying deformation using off-fault geomorphic strain markers provides valuable insight into potentially-active faults in coastal areas.
In subduction zones, the state of stress of the upper continental plate changes throughout the seismic cycle, both in space and polarity 18 . As a response to the drastic polarity change associated with a megathrust earthquake, the upper plate is commonly affected by enhanced extension promoting an increase in seismicity and occasionally triggered crustal earthquakes e.g. 19,20 . Slip on upper-plate faults triggered by megathrust earthquakes has been reported from Japan, Alaska, and Chile [20][21][22] , and has been inferred along most subduction zones e.g. 20,23,24,25 . However, historical and paleoseismic observations suggest that slowly-slipping crustal faults are characterized by recurrence times involving thousands of years, and therefore may not be reactivated during every megathrust earthquake, which commonly recur over several decades to a few centuries e.g. 20,24,26 .
Because upper-plate faults are widespread along coasts bordering subduction zones, they may pose local hidden hazards. In addition, such crustal faults may produce higheramplitude seismic waves at local scales than megathrust earthquakes and may locally enhance the amplitude and shorten arrival times of tsunamis in the near-field e.g. 27 .
Therefore, mapping crustal structures along subduction zones and quantifying their sliprate and relation with the megathrust earthquake cycle is a fundamental task to adequately assess the spatiotemporal characteristic of earthquake and tsunami hazards.
In this study we focus on the Pichilemu Fault (PIF), a partly-blind, hidden fault, which was unknown until it generated two shallow Mw 7 and 6.9 earthquakes 11 days after the M8.8 Maule megathrust event that affected central Chile in 2010. We show that on-fault displaced geomorphic markers are absent along the fault trace, but that off-fault strain markers may be used to estimate a long-term slip-rate, which when integrated over the 2010 coseismic deformation pattern yields a recurrence rate of such earthquakes. Our results show the hidden potential of blind faults with implications for seismic hazards along coastlines bordering subduction zones.

Seismotectonic and geologic setting of the Pichilemu Fault
The PIF is located at the central Chile margin, where the Nazca plate is subducted underneath South America at 66 mm/yr 28 . This region comprises the Coastal Range that reaches maximum elevations of ~600 m and formed almost exclusively by crystalline metamorphic rocks, related to a Paleozoic accretionary prism overprinted by brittle deformation during Mesozoic and Cenozoic exhumation [29][30][31] . The seawards slope of the range has a sequence of uplifted marine terraces, some of them overlain by shallow marine deposits 32 . The dense vegetation and thick soil cover have hampered mapping of geological and geomorphic features in this area, resulting in different interpretations regarding the presence of neotectonic structures e.g. 3,32 .
The Pichilemu area was affected by the 2010 Maule earthquake (M8.8), which ruptured a ~500-km-long portion of the megathrust with a peak slip of 17 m south of the PIF (Fig.   1A) e.g., 33 . The Maule earthquake triggered slip along the Santa Maria and Pichilemu coastal faults, located in the southern and northern parts of the rupture zone, respectively 3,21 (Fig. 1A). Before the 2010 earthquake, only a few faults affecting the crystalline basement and Cenozoic sedimentary cover had been previously mapped in this area 31,34 ( Fig. 1B); the PIF was unknown.
Eleven days after the Maule earthquake, the PIF slipped during two M6.9 and 7 normalfaulting earthquakes (Fig. 1B), followed by ~12,000 aftershocks located between the plate interface and ~4 km depth 35 (Fig. 1C). Aftershocks delineated a NW-striking and SWdipping 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 earthquakes were associated with a maximum slip of ~3 m along the main strand of the PIF, extending from 0 to 22 km depth 36 . However, despite the large magnitude of these earthquakes, no evidences of surface ruptures were neither found during field surveys nor detected in radar interferometry images 3,36 , which suggests the PIF is a blind, yet active structure.

On-fault tectonic geomorphology from bare-earth LiDAR topography
We used a Digital Terrain Model (DTM) derived from airborne Light Detection And Ranging (LiDAR) data at 1-m resolution to estimate fluvial metrics 37 and analyse the surface expressions of the PIF (See methods section). The coastal reaches of the PIF comprise two catchments of ~100 km 2 each (C1 and C2; Fig. 2A). A conspicuous set of NW-SE-oriented, subparallel lineaments exist in catchment C1 and extend for ~9 km; these are associated with aligned drainages and slope breaks ( Fig. 2A and Fig. S1). The catchments are developed entirely on the metamorphic bedrock, except near the coastline where the valley floors are filled by marine, fluvial and eolian deposits. The catchment asymmetry is manifested in the deviation between the position of the main river trunk stream and the catchment centreline ( Fig. 2A). The asymmetry of C1 and C2 is highlighted by trunk streams that converge along a section parallel to the edge of densely aligned aftershock seismicity and diverge from the catchment centreline in opposite directions ( Fig. 2A). Catchment C1 is characterized by higher relief, with drainages reaching ~500 m above mean sea level (amsl) and 16  Catchment and drainage metrics suggest variable surface deformation in this area. For instance, the marked catchment asymmetry suggests local block tilting of the PIF footwall and hanging-wall blocks. In addition, differences in Ksn, knickpoint locations, and drainage elevations may be indicative of differential vertical displacements ( Fig. 2A-D).
However, a detailed analysis of the PIF lineaments suggests that these features cannot be directly interpreted as fault scarps associated with surface breaks during recent earthquakes. Slope and elevation breaks are neither observed on maximum elevations derived from swath profiles nor on ridge-crest profiles (Figs. 2E and 2F). In addition, no clear relationship exists between the trace of these lineaments and the fringe of higher Ksn 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 38 (See Section 4.2 for further details). We found no relationship between chi-plot (c) values at the location of the knickpoints and their distance to the mapped lineaments (Figs. 2B and S2, Methods).
Furthermore, the distance between knickpoints and lineaments is rather uniform (Fig. 2B), suggesting that knickpoints are not related to active fault scarps, and rather reflect baselevel changes associated with relative sea-level variations. The results of morphometric analyses of LiDAR topography suggests that the area has been affected by surfacedeformation processes with a certain degree of spatial asymmetry, which may reflect tilting and differential uplift. However, these deformation patterns cannot be directly associated with any structure that has a marked surface expression. We therefore conclude that this area is not characterized by localized deformation at the surface, but rather by strain distributed over a 10-km-wide region.

Off-fault tectonic geomorphology: deformed marine terraces
Marine terraces -geomorphic markers of past relative sea-level positions 39 -are ubiquitous along the coast of central Chile 32 . We mapped terraces at Pichilemu using a LiDAR DTM and the TerraceM-2 software 40 32 ; thus, we correlate the surface morphology and geometry of these deposits with MIS 5e and MIS 9, at 125 and 320 ka respectively (Figs. 3A, 3B, S4, and S6). We tentatively correlate two additional terrace levels with MIS 7 and MIS 11, between 250 and 380 ka (See Fig. S6). South of the PIF, the lower marine terrace level is continuously exposed between Punta de Lobos and La Puntilla, with widths of 1 to 3 km.
These surfaces have been interpreted as a rasas 32 , i.e. terrace surfaces formed by reoccupation during successive highstands. Marine sediments that cover this rasa level are exposed at 16 and 32 m amsl and yielded post-IR IRSL ages of 328 ± 33 and 314 ± 30 ka, corresponding to MIS 9 (Figs. S3, and S6). The upper terrace levels are characterized by sharply defined paleo-cliffs and narrow paleo-platforms, with mean shoreline angle elevations at 42, 60, and 80 m and whose elevation decreases southward. We interpret the age of the three upper terrace levels to range between MIS 11 and 19 by correlating the sequence with a global sea-level curve of Rohling et al. 41 and Bintanja et al. 42 (See methods and Fig. S6).
The estimated uplift rates vary from 0.06-0.15 m/ka to 0.37-0.46 m/ka across the PIF ( Fig   3C), with associated 2s 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 of ~10 km wavelength, which is compatible with faster uplift and back tilting along the PIF foot wall block, and monoclinal roll-over folding in the hanging wall, which is consistent with a NW-striking, SW-dipping normal fault (Figs. 3B and 3C). 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 two faults (F1 and F2; Figs. 4C and 4D), which evidently offset the deformation pattern inferred from marine terraces. We generated two models, the first one was based on uplift rates estimated from measured shoreline angles, and the second ones was based on a smoothed surface derived from the interpolation of these measurements, which increases the area for comparing measured and modelled results. In the first model we obtained a slip-rate of 0.42 +0.05 -0.03 m/ka for F1 (90% confidence) ( Fig. 4C and 4D), for a fault extending from 26 to 0.4 km depth (Fig. S9B, S10D and S10E), and 0.06 +0.07 -0.03 m/ka for F2 with an updip depth of 1.0 km. The second model yielded a slip-rate of 0.42 ± 0.04 m/ka for F1 ( Fig. 4C and 4D) and 0.08 +0.07 -0.02 m/ka for F2. The best-fit updip slip depths are 0.8 and 1.2 km for F1 and F2, respectively (Figs. S9C, S10G and S10H). Importantly, both models predict similar deformation patterns and slip rates, with slip not reaching the surface. We select the second as the best-fit model that better reproduces the deformation pattern of marine terraces with lower uncertainty (Fig. 4D).

Coseismic slip and long-term slip-rate of the PIF
Our coseismic and long-term slip models suggest blind faulting with similar surface deformation patterns. The differences in updip depths between both models may be related to the relatively simple model setups, which assumes heterogeneous slip along a planar fault slip with homogeneous rheology. However, both models reproduce well the pattern and magnitude of surface displacements well, suggesting that deformation associated with the PIF occurs along a 10-km-wide area.

Earthquake recurrence of the Pichilemu Fault
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 triggered slip by megathrust earthquakes, such as during the events observed after the Maule earthquake. Historical and paleoseismic records suggest a recurrence time of ~0.1-0.2 ka for Maule-type events 44,45 , whereas paleo-tsunami records suggest recurrence times of 0.16 ka at the southern tip of the rupture 46 . If we consider a recurrence time of 0.2 ka, then a triggered slip of 0.084 m would account for the long-term PIF slip-rate (Fig. 5A); such slip would be equivalent to a seismic moment between 7.943 × 10 16 and 1.258 × 10 18 Nm or a Mw 5.2 to 6 earthquake, based on empirical relationships [47][48][49] (Fig. 5B). In turn, a recurrence time of 0.5 ka would imply 0.21 m of slip along the PIF (Fig. 5A) and an earthquake magnitude between Mw 5.7 and 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, in order to trigger a Mw 7 earthquake, 1.65 m of slip per event would be necessary. In this case, the recurrence time required to account for the long-term slip-rate would be 3.9 ka. If the probability distributions of the long-term PIF slip-rate and coseismic slip intersect, a recurrence time of 2.5 ka (2.1 to 3.3 ka at 95% confidence interval) may be inferred (Fig. 5A). Such a recurrence time would generate a Mw 6.6 to 6.9 earthquake, equivalent to a seismic moment release between 1.0 and 2.8 × 10 19 Nm (Fig. 5B), similar to 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 M>8 megathrust earthquakes in the Maule segment, implying sporadic slip triggering.

Pichilemu Fault triggering by megathrust earthquakes
Stress transferred to the upper plate after a megathrust earthquake may induce reactivation of crustal faults optimally oriented with respect to the new stress field e.g. 50  It has been suggested that the along-dip segmentation of megathrusts may control the magnitude and characteristics of subduction earthquakes 53 . We compared CFS at the PIF induced by earthquakes occurring at variable megathrust down-dip depths by simulating synthetic slip scenarios (Fig. S12). Previous studies have shown that stress-field perturbations of 0.1 to 0.3 MPa may be sufficient to trigger seismicity, while reductions of the same amount may prevent them [54][55][56] . Certainly, the threshold of stresses for fault reactivation will mostly depend on the amount of elastic strain accumulated by the fault; however, for the purpose of comparison we set a tentative CFS threshold at 0.2 MPa, which is below the minimum induced at the PIF during the Maule earthquake. We observe that when the locus of slip is located deeper down-dip along the megathrust (between 40 and 85 km depth), CFS values decrease (from 0.18 to 0.01 MPa) below the reactivation threshold (Fig. 5D). When the locus of slip is located at shallower depths (15 to 40 km depth), CFS values increase reaching 0.28 to 0.35 MPa, within the tensional stress scenario that may promote reactivation of the PIF (Fig. 5D). However, megathrust slip at depths shallower than ~15 km decreases CFS on the PIF, reaching values near the proposed CFS threshold for reactivation (Fig. 5D). The comparison of different scenarios suggests that triggered reactivation of the PIF will be favoured when megathrust slip is located in the vicinity of its down-dip extent.
The CFS modelling shows that different slip distributions along the megathrust may control the reactivation potential of the PIF, promoting or inhibiting its reactivation. These results may explain why the estimated 2.5 ka recurrence time of the PIF is longer than expected for megathrust earthquakes in the Maule segment. We consider two possible scenarios to explain this: 1) The Maule segment may be characterized by variable slip during megathrustearthquake ruptures. In this case megathrust earthquakes may trigger the reactivation of PIF when the PIF is favourably oriented to the megathrust slip pattern, which may occur occasionally. During the interseismic phase of the megathrust earthquake cycle, the PIF would accumulate elastic strain.
2) The Maule segment may be characterized by persistent earthquake ruptures. In this case, positive CFS values induced by successive earthquakes may build up tensional stresses on the PIF, which may be reactivated only when enough elastic strain has been accumulated and when the fault is favourably oriented to the stress field. A minimum amount of elastic strain accumulates during the interseismic phase, counteracting tensional stresses.
Regarding the first scenario, the relation between the coseismic slip distribution during the Maule earthquake and long-term deformation and the observed geophysical characteristics suggests three persistent segments 32,57 ; however, the magnitude of coseismic slip on these segments may vary between cycles 57 . Furthermore, the CFS values on the PIF are negative during the interseismic period that preceded the Maule earthquake 33 (Figs. 5A and B) suggesting no elastic strain is expected to be accumulated during the interseismic phase. Regarding the second scenario, a similar interpretation has been proposed to explain the reactivation of crustal faults in northern Chile 58 , where negative CFS during the interseismic stage may counterbalance the CFS induced by coseismic slip, extending the recurrence time of crustal fault earthquakes and providing an explanation for the long recurrence times that we estimated for PIF earthquakes.

Origin of the PIF and implications for blind faults
Blind faults mostly occur in sedimentary basins, when fault tip propagation fails to reach the surface across a thick sedimentary cover 59 . This is, however, not the case of the PIF, which mostly affects crystalline basement rocks. On the other hand, blind structures may be also controlled by mechanical and rheological heterogeneities in the upper crust e.g. 59 .
The PIF host rocks are part of a Palaeozoic accretionary wedge characterized with pervasively-deformed high-pressure metasedimentary and metavolcanics rocks locally forming melanges [29][30][31] . These metamorphic rocks are 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 wide deformation zone in the uppermost crustal levels, preventing localized surface faulting. This assessment is further supported by the distribution of aftershock seismicity following the 2010 earthquakes across a ~10-km-wide fringe ( Fig. 1C and S8), and by a sharp increase in p-wave velocity at ~1-km depth inferred from seismic tomography 60

Estimating coseismic slip
We processed two Envisat® scenes obtained 2 days before and 7 days after the March 11 2010 earthquake doublet, obtaining line-of-sight (LOS) displacements ( Fig. S7 and Table   S1). To corroborate our results, we compared the LOS surface displacements with displacements estimated using positions of the permanent GPS station PICH (location in Fig. 3D) projected to the LOS vector. To estimate coseismic slip along the PIF, we searched for the parameters that better reproduced the distribution of LOS displacements using forward elastic dislocation modelling 66 . The elastic dislocation models were programmed in Matlab® using the function okada85© 67 . We inferred fault geometries from alignments in the cluster of crustal aftershocks associated with the PIF earthquakes (76º dip, N41ºW strike and 80 km length) (Fig. S8). The strike inferred from aftershocks is similar than previous interpretations and fault planes from focal mechanisms 36 and to the strike of lineaments mapped using the LiDAR bare-earth topography. To compare between model and observations we iterate the slip, updip, and down-dip for defined ranges and constant increments, generating 68.921 models (Figs. S9A and S10), and searching from the parameters that minimized the Normalized Mean Root Square Error (NRMSE). Eq.1: The root mean squared error (RMSE) is defined by the difference between observation (yi) and model (y), n is the number of observations (Eq.1).

Eq.2:
The NRMSE facilitates comparing models with different scales by normalizing 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 define the model uncertainty as the interval between the 5% and 95% percentiles (Fig. S10).
To evaluate the consistency of fault slip estimates, we compared our result with a 3D slip inversion of the coseismic displacement field estimated from InSAR data (Fig. S11). For the inversion we used the elastic dislocation model and the surface deformation constrained by InSAR, the fault geometry from focal mechanisms of the Pichilemu earthquakes, and a Poisson ratio of 0.25 following previous studies in the region 57,58 . To solve for the slip distribution, we used a least squares minimization determining the bestfit solution that reproduce the InSAR deformation patterns.

Analysis of on-fault geomorphic features in the Pichilemu area
We performed a detailed geomorphic and morphometric analysis of off-and on-fault geomorphic features in the PIF area. On-fault features where searched by analyzing topographic and fluvial metrics using the Topotoolbox-2 software 37 . We created red-relief maps (RRM) based on Chiba et al. 68 , which uses the terrain openness and surface slope.
RRM are useful to search for lineaments and to identify fault scarps as well as changes in the fluvial network because these maps lack the potential bias of light source direction in common shaded relief images (Fig. S1).
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 in homogeneous climatic conditions 69

Eq.3:
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.4: Were (U/K) 1/n represents the channel steepness and m/n is the channel concavity ( , thus the equation can be written as: Eq. 5: In contrast to Ksn, chi-plots are based on the horizontal transform of the upslope area to linearize the concave upward profile for a well-chosen reference concavity. This spatial transform makes chi-plots useful to identify transient erosional signals with a common origin propagating upstream along fluvial systems, such as tectonically-generated knickpoints or changes in base level. Furthermore, the chi dimension allows comparing these signals between different catchments independently of their size or shape.
Determining 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. 6: The integral of dz/dx is calculated using a reference area (A0), where x0 is the catchment outlet. The chi-plot values (c) are estimated based on the right hand-integral: Eq. 7: To calculate the Ksn and chi-values in the catchments of the Pichilemu area we used Topotoolbox® 37 , which includes all the indicated equations. We used a reference concavity of 0.45.

Analysis of off-fault geomorphic features in the Pichilemu area
To analyse off-fault geomorphic features, we studied marine terraces using two methods depending on their type and characteristics:

Analysis of wave-built terraces
We analysed wave-built marine terraces following 71 considering the morphology of the bedrock unconformity, number of sedimentary cycles within the wave-built terrace, and thickness of the sequence. We mapped the surface morphology of the wave-built terraces using swath profiles to detect slope breaks. In addition, we measured the depth to the crystalline bedrock in incised valleys and generated an isopach map of sedimentary sequence thickness (Fig. S4). This allows differentiating sedimentary sequences and improving the estimation of shoreline angle elevations.

Analysis of wave-cut terraces
We studied the surface morphology of marine terraces using LiDAR topography and swath profiles, in order to measure the location and elevation of the shoreline angle. The shoreline angle is a geomorphic marker located at the intersection between the paleoplatform 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 39 . To estimate uplift rates (u), we correlated terrace levels and sea-level highstands using the IRSL ages and a sea-level curve for the southern hemisphere spanning the last 700 ka 41,42 following Lajoie 39 as: where E is the elevation of shoreline angles, (e) the position of the respective highstand and T the age of the terrace level. Uplift rate errors are estimated following Gallen et al.,72 as: where σH is the error in relative sea-level defined as: Eq. 10: σH = where σT is the age uncertainty in the sea-level curve, σE is the error of the shoreline angle assessments and σe is 12 m uncertainty of the highstand elevation based on Rohling et al. 41 . We used an arbitrary uncertainty of 7 ka for the duration of the highstands (σT) following a previous study e.g. 32 .

Post-IR IRSL dating
Four samples for thermoluminescence dating were analyzed at the University of Cologne using the post-IR IRSL signal of K-feldspar obtained from marine terrace sediments. We selected sedimentary units for dating using three criteria: 1) sedimentary units formed by shallow marine sediments, preferable from berm or swash zone environments. 2) Sandy sediments of medium grain sizes with more than 20% of feldspars. 3) Samples from the base of the sequence as near as possible to the bedrock wave-cut platform (Fig. S3). We analyzed medium sand K-Feldspar grains (100-250 µm) following the post-IR IRSL290 SAR protocol 73 . The dated sediments were in general characterized 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 -1.1). Burial doses were based on 5 -19 aliquots of 8 mm diameter and using the central age model (CAM) of Galbraith et al. 74 for calculation (See radial plots inf Fig. S5). In addition, we evaluated the completeness of signal resetting using a sample from a modern beach berm (Sample SM mb, Uplift rates derived from shoreline angle elevations of marine terraces were reproduced by forward elastic dislocation models by varying up-and down-dip slip depths and the sliprate of each fault. The elastic model setups include the same geometry and down-dip depth used in 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 updip depths of F1 and F2 using constant increments. The best-fitting models were selected from the minimum NRMSE value. To estimate the confidence interval of the best-fit models, we used the lower 5% tail of the NRMSE distribution (Fig. S10) to calculate slip and slip rates. From these distributions, we define 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. S10).
Because marine terraces are exposed along the coast, they mostly reflect deformation along a 2D profile, and are therefore not suitable for estimating fault-slip with a 3D inverse model. In order to produce comparable results at the coseismic and long-term time-scales, and given that forward and inverse models of InSAR displacements produce similar results, we compare results using forward dislocation models.

Modelling Coulomb stress failure
We evaluate the reactivation potential of the PIF by a megathrust earthquake by analyzing

Acknowledgments
We would like to thank Kenneth Fisk and Christian Hilleman for help in the field and Juan