Locating Fully Locked Asperities Along the South America Subduction Megathrust: A New Physical Interseismic Inversion Approach in a Bayesian Framework

The largest earthquakes in subduction zones occur where significant interseismic slip deficit has accumulated on the plate interface. Slip deficit accumulates most quickly in mechanically locked regions, and these also cause the regions around them to accumulate slip deficit; therefore, large earthquakes are typically expected to rupture in and around locked areas. The locations and dimensions of these locked zones have been difficult to resolve using standard techniques and available data sets. We develop a new statistical interseismic inversion approach that incorporates the physical interactions between nearby fault areas to directly determine the distribution of locking on the subduction plate interface (simultaneously with rigid forearc motions) from interseismic surface velocities. Because we include physical prior information in the inversion procedure, this approach reduces uncertainties in the rate of slip deficit accumulation, even in locations (such as near the trench) where kinematic inversions of onshore data have relatively low resolution. Applying the inversion to the South America subduction zone, we find that the pattern of locking and corresponding slip deficit rates correlate well with recent and historical large earthquake ruptures. Locked patch dimensions are <40 km and account for no more than 30% of the area of the plate interface. The small size of the imaged locked zones is a natural outcome of our physical assumptions and implies that mechanical locking is caused by correspondingly small geological features. Despite their small dimensions, locked zones generate substantial slip deficit on the surrounding plate interface, consistent with the slip patterns of large megathrust earthquakes.


Introduction
Landward interseismic velocities on the upper plate in subduction zones, such as those recorded at Global Navigation Satellite System (GNSS) stations, are a consequence of the mechanical coupling between the subducting slab and upper plate Savage, 1983). These geodetic observations indicate that there are regions where the plate interface ("megathrust") either remains completely locked or slides at a rate substantially less than the convergence rate ("asperities"; e.g., Gombert et al., 2018;Li & Freymueller, 2018;Loveless & Meade, 2016;Metois et al., 2016;Nocquet et al., 2014). The same set of observations suggests that regions outside of asperities slip aseismically at a faster rate, up to the relative plate velocity (via a mechanism such as stable frictional sliding or ductile creep).
The location and behavior of asperities are important for understanding the occurrence of subduction zone earthquakes because asperities accumulate both slip deficit (equal to the total plate convergence minus the actual fault slip; Savage, 1983) and shear stresses at the highest rates of any region on the megathrust. The total accumulated slip deficit (given the past history of earthquakes along the megathrust) defines the theoretical maximum possible coseismic slip that can occur during an earthquake. At the same time, in order to nucleate an earthquake, shear stress on the plate interface must reach a brittle strength threshold, controlled by its frictional properties (e.g., Marone, 1998;Scholz, 1998). Most previous interseismic studies have focused on mapping the observed surface motions directly to the distribution of slip deficit rates and interpreting the areas accumulating slip deficit at high rates to represent the locations and dimensions of asperities. The asperities are subsequently associated with large earthquake ruptures (e.g., Chlieh et al., 2008;Metois et al., 2016;Moreno et al., 2010;Nocquet et al., 2017) and major geological features in the subduction zone (e.g., incoming sediment thickness: Heuret et al., 2012, Olsen et al., 2020; roughness of the incoming plate: Carena, 2011;van Rijsingen et al., 2018;Wang & Bilek, 2014; plate interface geometry: Bletery et al., 2016;or upper plate geology: Béjar-Pizarro et al., 2013;Liu & Zhao, 2018). However, there are still significant uncertainties in the distribution of slip deficit in subduction zones, due in large part to observation geometries. These uncertainties preclude the straightforward identification of mechanical asperities and their correlation with seismicity and geological features. In this study, we circumvent this problem by explicitly introducing the physical effects of asperities into the interseismic analysis, which allows us to directly estimate the locations and dimensions of these asperities. This approach has the additional benefit of improving constraints on the magnitudes of slip deficit rates by introducing physical a priori information on the modeled slip rates.
An implicit assumption in most previous kinematic slip deficit inversions is that the distribution of fault slip does not need to be consistent in the context of continuum mechanics (Figures 1a and 1b). In other words, slip in each area of the plate interface is determined independently from its neighbors, other than perhaps a smoothing operator which is not based on physical constraints. However, several studies have demonstrated that low-slip-rate asperities restrict the regions around them from slipping at the convergence rate because of the mechanical continuity of the plates, implying that both asperities and the plate interface around them should accumulate slip deficit (Almeida et al., 2018;Bürgmann et al., 2005;Gans et al., 2003;Herman et al., 2018;Johnson & Fukuda, 2010;Malservisi et al., 2003;Wang & Dixon, 2004;Figure 1c). We follow these previous studies and model the plate interface as being either completely locked or unlocked ("binary locking"); this is an end-member representation of the megathrust rheology. As Herman et al. (2018) explicitly tested, a mechanically locked boundary condition remains a good approximation for asperities as long as they slip at <5% of the convergence rate, irrespective of the specific rheology or frictional properties of the megathrust. Therefore, we will simply refer to these asperities as "mechanically locked" from this point onward. In this physical framework, shear stresses increase in mechanically locked regions such that they accumulate slip deficit at the plate convergence rate. The evolution of shear stresses outside of locked asperities remains less well constrained; however, these regions are generally thought to slide more stably, meaning that the resistive shear tractions remain relatively constant with continuing fault slip (Figures 1c and 1d; e.g., ). We will therefore simply refer to these areas as "unlocked," and model them as sliding stably. We later discuss how our analysis would change if the shear tractions outside of locked asperities do not remain constant. These unlocked regions accumulate interseismic slip deficit due to their proximity to locked patches ( Figure 1c); in the absence of nearby locked asperities, these unlocked regions would be able to slide at or near the convergence rate.
Although these effects have been described conceptually and previously modeled, they have never been systematically incorporated into an inversion to constrain the distribution of mechanical locking along the subduction megathrust. We therefore develop a new procedure to include these mechanical effects in the inversion, allowing us to statistically determine which areas of a plate boundary are locked or unlocked based on interseismic velocity observations (section 2). We then apply this new analysis to the well instrumented and seismically active South America subduction zone (sections 3 and 4).

Computing Fault Slip
In this study, we assume that every location on the plate interface is in one of two possible states: fully locked with zero slip or fully unlocked with no resistance to slip (Figures 1c and 1d). Previous studies used finite element models (FEMs) to determine the distribution of fault slip rates under these binary locking conditions Herman et al., 2018;Malservisi et al., 2003), but FEMs are too computationally expensive to run many times in a statistical search. For example, the converged 3-D FEM of Herman et al. (2018) takes up to 60 s to run on a single processor. If instead the Green's functions relating fault slip to shear tractions are precomputed (using an FEM or analytical solutions), the fault slip can be calculated much more rapidly (~1 s for a comparable model). Therefore, in this study, similar to Almeida et al. (2018), the slip deficit rate in unlocked faults ( _ u j ) surrounding locked faults is computed by solving the following equation, which directly minimizes the shear stresses on those faults: Figure 1. Conceptual model of the fault slip rate and shear stressing rate along a subduction megathrust (location shown in inset) for a standard kinematic inversion approach and the mechanical framework incorporated in this study. (a) In standard kinematic inversions, a slip rate is prescribed in a fault section. Implicitly, the plate interface outside of these areas is assumed to slip at the convergence rate, unless observations or regularization require slip in those regions. (b) The shear stressing rate is nonzero both inside and outside of asperities in the standard approach. This stress distribution is inconsistent with respect to continuum mechanics. (c) In this study, we impose mechanical consistency on the fault slip rates. Therefore, locked asperities restrict the regions around them from sliding at the full convergence rate. (d) We assume that the plate interface outside of locked regions slides stably, so that it accumulates zero shear stresses. Shear stresses only accumulate within locked zones so that these areas do not slip.
where L is the set of locked faults, G ij are the Green's functions relating slip rate to shear tractions for all unlocked fault pairs, and t k i are the shear tractions produced by locked faults that accumulate on unlocked faults (full derivation in the supporting information Text S1). We verify that under similar boundary conditions, an FEM representation of fault slip around a locked zone converges to the same solution ( Figure S1).
Although binary locking may be a simplification of the plate interface rheology, it produces similar results as more complex representations of the megathrust (e.g., Corbi et al., 2017;Herman et al., 2018;Hetland & Simons, 2010;Kaneko et al., 2010). This is because the continuum mechanics underlying Equation 1 are valid independent of the specific megathrust rheology. One important scenario to consider is when the unlocked faults do not slide stably, so that they also accumulate shear tractions. Equation 1 can be modified to include this behavior by simply adding the shear traction accumulation rate in each unlocked fault to the right-hand-side (Text S1). Although it is beyond the scope of this study, these tractions could be investigated as a free parameter in the inversion.
A benefit of incorporating physics-based mechanical continuity into the forward model is that the distribution of fault slip rates is automatically smoothed and less than the convergence rate, so there is no need for additional regularization constraints. However, it also results in a nonlinear pattern of slip around the locked asperities Figure 1c), making it more difficult to use in standard least squares inversion approaches. By running forward models in a direct search for the distribution of locked patches, we avoid the complications associated with slip rate nonlinearity altogether.

Locking Search
We test many combinations of locked and unlocked fault areas ("sub-faults") to determine which distribution of locked zones produces good fits to observed interseismic velocities. This direct search approach ensures that every model has a physically consistent distribution of slip deficit rate and yields robust statistical estimates of the locking distribution. As the number of sub-faults increases, it becomes computationally intractable to test every locked/unlocked combination. Therefore, we use a Markov Chain Monte Carlo approach; specifically, we apply the Metropolis-Hastings algorithm (Hastings, 1970;Metropolis et al., 1953), which is becoming ubiquitous in the Earth Sciences for its relative ease of implementation and its utility in high-dimensional parameter spaces. This algorithm produces a set of models (each consisting of a set of locked and unlocked sub-faults) whose parameter values (e.g., locking and slip deficit) are distributed proportional to their probability (Sambridge & Mosegaard, 2002). The general form of the Metropolis-Hastings algorithm is as follows: In iteration zero, an initial model is selected from the parameter space. This model is set to be the "current" model. In each iteration of the search, a new ("proposed") model is generated from the current model by a transition operation (our discrete problem requires a different transition operator than the typical random sample from a multivariate Gaussian: details provided below). If this proposed model produces a better fit to the data than the current model (i.e., a higher objective function, φ), then it is accepted as the current model. If the proposed model produces a worse fit to the data (lower objective function), then it may still be accepted as the current model, with probability: Once the current model has been either updated or retained from the previous iteration, the next iteration begins.
These competing options (always taking better proposed models and allowing worse proposed models to be accepted) allows the algorithm to both find local optima and also potentially emerge from them to discover new local optima. After an initial "burn-in" set of iterations, in which the search can discover a region of parameter space that produces good fits to the data, the Metropolis-Hastings algorithm samples the parameter space proportionally to the probability of the parameters, as defined by the choice of φ, the objective function. In the context of Bayes' theorem, the algorithm samples from the posterior probability density function, p(M| D) ("the probability of the model parameters-in this case, the sub-faults being lockedgiven the observed data set"): where p(D| M) ("the probability of obtaining the data set, given a set of model parameters" or "likelihood") is proportional to the data fit and p(M) ("the probability of the model parameters") contains prior information about the model parameter values.
For our analysis, we initialize the search with a random distribution of locked patches. The advantage of randomly initializing the search (instead of starting with a model that produces better fits to the data) is that the algorithm can discover optimal locked/unlocked combinations without user bias. The cost of this choice is that more initial iterations may be needed to reach the locked distributions with relatively good fit. In each iteration, the proposed model is generated by taking the current model and flipping faults to/from locked and unlocked. The number of faults to flip is randomly selected from a uniform distribution between one and some user-set maximum value in every iteration. Although this maximum value can be set to any number and the algorithm will still (in theory) sample p(M| D) given enough iterations, a better proposal choice can make the algorithm more efficient (Bagnardi & Hooper, 2018;Minson et al., 2013;Sambridge & Mosegaard, 2002). Our preliminary tests of the search suggest that flipping a maximum of 2-3% of the total number of sub-faults provides a balance between large enough perturbations to explore the parameter space and small enough to map the probability distribution in sufficient detail, thus reducing the number of forward models that need to be run. We find that the algorithm becomes particularly inefficient (sometimes never discovering good-fitting locking distributions) when the number of faults to flip is too large, greater than~15% of the total.
In this study, we maximize the −χ 2 misfit between the model-predicted and observed GNSS velocities: where v⃑ o are observed velocities, v⃑ p are predicted velocities (computed using the direct shear stress minimization approach; section 2.1), and C −1 is the inverse of the data covariance matrix. This corresponds to Gaussian data uncertainties, so that the likelihood is where N is the number of velocity components and |C| is the determinant of the covariance matrix.
In synthetic tests of this algorithm, where the model parameterization is known perfectly, C contains only the covariances of the observed interseismic velocities (C obs ). In real applications, uncertainties due to the modeling choices can be as high as 15-20% of the signal magnitude (the "fractional error"; Minson et al., 2013). In these cases, the covariance matrix can incorporate both observation and model uncertainties, that is, C ¼ C obs +C model . For efficiency in the calculation, we set C model to be diagonal: Here, α is the fractional error and v o,i is the ith observed velocity component. We take the approach of setting the fractional error to a fixed value at the beginning of the search.

10.1029/2020GC009063
We also account for the possibility of forearc motions superimposed on the megathrust locking signal by adding Euler poles describing rigid plate motions to the search. In each iteration of the algorithm, a new proposed location and angular velocity of the pole are sampled from a user-specified Gaussian distribution around their values in the current model. Rigid velocities are calculated for each GNSS station in the corresponding forearc sliver and added to the predicted locking-generated velocity; this vector sum is now v⃑ p .

Synthetic Tests
Synthetic surface velocities are generated from a model with a plate interface dipping at 15°and extending 500 km downdip and 1,000 km along strike. This fault is divided into 50 square sub-faults, each 100 km to a side. Four sub-faults are locked by assigning them to have backslip of 70 mm/year (the synthetic plate velocity; Figure S2). The slip deficit rate is computed in the 46 unlocked sub-faults. The synthetic surface velocities are calculated for this distribution of backslip at 60 locations on the upper plate with a better offshore distribution than is typically available at subduction zones. Gaussian noise with standard deviation 0.2 mm/year (comparable with the lowest uncertainties reported for interseismic velocity observations) is added to all three components of these synthetic velocity observations. Both shear stressing rates and surface velocities are computed using the equations for rectangular dislocations in an elastic half-space (Okada, 1992).
The inverse models use the same plate interface geometry and all three velocity components. The Metropolis-Hastings algorithm is run for 100,000 iterations (with the first 1,000 iterations for burn-in) to build the locked/unlocked probability density function. We test the effect of varying the inverse model sub-fault dimensions (with 50-, 100-, and 200-km squares) and shapes (100-150 km wide triangles; Meade, 2007). In all of these tests, regions with highest locking probability overlap with the locked zones of the input model. The corresponding inverted surface velocities also fit the synthetic velocities well.
Larger sub-faults result in greater confidence in locking probability at the expense of spatial detail. Smaller sub-faults provide better spatial resolution, but the locking probability is spread over many adjacent sub-faults.
Using the same model geometry, we also investigate the size of a single locked patch that can be resolved ( Figure S3). The locked patch is placed on the megathrust centered at a depth of 50 km, backslip rates are computed on the surrounding interface, and synthetic surface velocities are computed (with noise) at the same grid of points as in the previous synthetic test. Using an inverse model with triangular sub-faults 100-150 km wide, we test its ability to constrain square locked patches with side lengths of 20, 40, 60, and 80 km. We find high locking probabilities that overlap with the locked zone down to a size of 60 km. For a 40-km locked zone, the search identifies its presence, but mislocates it. A 20-km locked zone is not seen in this test search.
Based on these synthetic tests, we use a minimum sub-fault size of 40 km in our analysis of the South America subduction zone. We did not consider additional uncertainties (related to either the measurements or the modeling choices) in these proof-of-concept synthetic tests. Therefore, the locked patches in real applications may be more difficult to resolve. Nevertheless, these tests provide a baseline spatial scale for the dimensions of locked patches that can be identified along a subduction megathrust.

South America Seismotectonics
In the South America subduction zone, the Nazca plate descends eastward under the South America plate along a trench extending from southern Chile to Colombia (Figure 2). The convergence velocity varies along strike, from 75 mm/year at the southern end of the subduction zone to 60 mm/year at the northern end, with a nearly constant azimuth (75°-78°) (DeMets et al., 2010). There have been 12 moment magnitude (Mw) 7.5+ megathrust earthquakes in this subduction zone since 1995 (Hayes, 2017). There are many published slip models for these events, and on the scale of the subduction zone all of these generally have similar epicenters and rupture areas as the Hayes (2017) models (e.g., the models in the SRCMOD database: Mai & Thingbaijam, 2014; or many other published models, e.g., Duputel et al., 2015;Herman et al., 2017;Klein et al., 2016;Lay et al., 2010;Nocquet et al., 2017). This subduction zone also has historical records of large magnitude megathrust earthquakes over the past 400-500 years (Kelleher, 1972;Lomnitz, 2004).

Interseismic Velocity Observations
The data set we use to constrain the locking distribution in the South America subduction zone consists of interseismic horizontal velocities relative to fixed South America at 632 GNSS sites on the upper plate measured over the past 20 years (Blewitt et al., 2016;Brooks et al., 2003Brooks et al., , 2011Chlieh et al., 2004;Kendrick et al., 2001;Khazaradze & Klotz, 2003;Klein et al., 2018;Klotz et al., 2001;Metois et al., 2012;Metois et al., 2013;Metois et al., 2014;Nocquet et al., 2014;Ruegg et al., 2009;Villegas-Lanza et al., 2016; Figure S4). We do not use vertical motions because of their sensitivity to model rheology. These velocities reflect different interseismic measurement periods uncontaminated by coseismic and postseismic motions. The stations extend along the entire length of the subduction zone from southern Chile (37°S) to Colombia (4°N), with few along-strike gaps greater than 100 km. The width of this combined data set measured perpendicular to the coast varies from <20 to >250 km.
Along much of the length of the subduction zone, the interseismic velocities are oriented ENE, with magnitudes of 20-30 mm/year near the coast decreasing to 0 mm/year in or east of the Andes. This is generally consistent with a locked plate interface and the relative motion of the Nazca plate driving South America eastward. There are also velocities oriented nearly parallel to the subduction margin (perpendicular to the plate motion) in northern Peru and southern Ecuador. These can be fit reasonably well with a single Euler pole and have been interpreted to reflect the motion of a rigid forearc sliver in this region (Nocquet et al., 2014). Similar forearc slivers have been interpreted in northern Ecuador and Colombia (Nocquet et al., 2014; directed to the northeast) and in Chile (Metois et al., 2016; directed to the east, parallel to the direction of shortening from locking).

Inversion Setup
The geometry of the subducting South America slab is defined by the Slab2 model (Hayes et al., 2018), divided into triangular sub-faults (Shewchuk, 1996; Figure S5). Sub-faults between the trench and the coast have side lengths of 40-50 km (areas of 700-1,200 km 2 ). The sub-fault dimensions increase with depth to minimize computation time, since these sub-faults will have low slip deficit rates ( Figure S2). This mesh size provides a compromise between capturing spatial variations in locking and analyzing the entire subduction zone. The convergence velocity at each sub-fault is calculated from the MORVEL plate motion model between the Nazca and South American plates (DeMets et al., 2010).
The velocity and shear stress Green's functions are computed using the analytical solution for uniform slip triangles in an elastic half-space (Meade, 2007) with a Poisson's ratio of 0.25 and a shear modulus of 40 GPa. We opt to use elastic half-space Green's functions because we can calculate them rapidly; however, because of this choice, we only use the horizontal velocity observations in order to avoid artifacts caused by mapping vertical motions into locking . The biases and uncertainties associated with the model rheology are incorporated into the data covariance matrix as a 20% fractional data covariance for every velocity.
We randomly initialize the locking distribution in the search by setting each sub-fault to be locked with a probability of 0.5. We only allow faults shallower than 80 km to become locked throughout the search.  (Hayes, 2017). The estimated lengths of historical ruptures are shown as red lines to the west of the trench, with distance from the trench indicating years before present; the largest events are labeled with their date (Kelleher, 1972;Lomnitz, 2004). The locations of proposed forearc slivers (Metois et al., 2016;Nocquet et al., 2014) are indicated by dashed lines.
Deeper faults are still included because they can contribute to the surface motion. In each iteration, up to 15 faults (~2% of the faults in the seismogenic zone) can flip between locked and unlocked.
We also search for the Euler poles describing the motions of forearc slivers in Chile, Peru, and Ecuador-Colombia (Figure 2; Figure S4). The initial locations and angular velocities of the Euler poles are taken to be near the published values for the Colombia-Ecuador sliver (83.4°W, 15.2°N, 0.287°/Ma) and the Peru sliver (63.8°W, 22.5°N, 0.092°/Ma) (Nocquet et al., 2014), and for the Chile sliver (138.7°E, 56.4°N, 0.120°/Ma) (Metois et al., 2016). The widths of the proposal distributions arẽ 25 km for the location and~0.005°/Ma for the angular velocity.
Each instance of the search runs for 50,000 iterations, with the first 5,000 iterations treated as burn-in. To better map out the locking and Euler pole probabilities, we run 20 separate searches (each starting with a different distribution of locked patches). These parallel chains may not necessarily be log-likelihood comparable, for example, if one dominantly samples an area of p(M| D) with low probabilities and another samples a region with high probabilities, these chains cannot be directly compared with each other. We therefore verify that our final locking distribution is robust by evaluating the results from different search configurations: (a) a single 50,000 iteration chain; (b) a simple combination of all 900,000 nonburn-in models; (c) a random resampling from the 900,000 non-burn-in models proportional to each model p(D| M) (e.g., Ching & Chen, 2007;Minson et al., 2013); and (d) a single, 500,000 iteration chain ( Figure S6). All of these results have similar locking probabilities, suggesting that the chains are converging to the same locking probability distribution. In the following sections, we show and discuss the results from the resampled ensemble of models.

Locking Probability
The fundamental result of the search is the probability that each sub-fault is locked (Figure 3). We calculate this probability by taking the number of iterations in which the sub-fault is locked and dividing by the total number of iterations in the search. We can also use the results of our search to determine how much of the plate interface area is locked. For five sections of the subduction megathrust that have hosted recent large earthquakes, as well as the region offshore of northern Peru, we compute the fraction of the megathrust area that is locked in each iteration of the search. Combining all iterations yields a probability distribution for that region showing the fraction of the interface that is locked (Figure 3, insets).
Most of the Chilean section of the subduction zone has a nonzero probability of being locked, but there is significant along-strike variability. Starting in the south, the sub-faults on interface from 38°S to 35°S all have a moderately high locking probability (0.4-0.6). From 35°S to 34°S, this high locking transitions to a section of the subduction zone with broadly low locking (34°S-29°S), except for a narrow section with high locking probability at~31°S. North of 29°S, the locking probabilities remain in the range 0.3-0.6 up to northern Chile at 18°S. Similar moderate to high locking probability continues through the bend in the subduction zone into the area of the plate interface offshore southern Peru (up to 12°S). At this latitude, there is another transition to low locking probability over a distance of~100 km. This section of low plate interface locking extends over 600 km along the northern Peru coast, from 10°S all the way  Figure 2) and low locking probabilities are associated with no recent earthquakes and the endpoints of historical ruptures. The histograms indicate the fraction of the outlined regions that are locked. No region has more than~30% of its megathrust area locked.
to 4°S. From this point northward, there is a higher probability of locking on the plate interface in the section of the plate interface offshore Ecuador up to 2°N.
The highest probability of any individual sub-fault being locked is 0.6 because the observations (plus uncertainties) cannot be used to distinguish between the signal of locking in that sub-fault or one of its neighbors  (Nocquet et al., 2014); (b) the Peru forearc sliver (Nocquet et al., 2014); and (c) the Chile forearc sliver (Metois et al., 2016). 10.1029/2020GC009063 with greater confidence. This is analogous to the resolution limitations inherent in kinematic slip deficit inversions at subduction zones. Based on our resolution tests, we estimate that we can detect the presence of an isolated locked zone down to~25 km wide. Although such a small locked zone does not generate enough of a signal by itself to be imaged, it produces slip deficit over a large enough area (a region~80 km wide) for an Mw~7.5 earthquake .
The median fraction of the plate interface area that is locked in the selected regions ranges from 0.01 to 0.24, with uncertainties ≤0.08. At most,~0.30 of the plate interface area is locked. The Maule (33°S-38°S) and Iquique/Pisagua (18°S-22°S) sections of the subduction zone offshore Chile have 0.20-0.25 of their area locked, whereas the Illapel/Coquimbo section offshore Chile (29°S-33°S), the Arequipa/Pisco regions offshore southern Peru (17°S-12°S), and the Pedernales region offshore Ecuador (4°S-2°N) have lower locked fractions of 0.08-0.16. The long, unlocked section of the interface offshore of northern Peru has 0.01 of its area locked, consistent with the near-zero locking probability in all sub-faults in this region.

Rigid Forearc Motion
In the Ecuador-Colombia forearc, our search shows a mean Euler pole location at (286°E, 21°N) with a mean angular velocity of 0.05°/Ma, corresponding to east-southeastward velocities of 1-2 mm/year (Figure 4a). The Euler pole describing the motion of the Peru forearc sliver is located at (304°E, 25°N), with an angular velocity of 0.03°/Ma, corresponding to southeast-directed sliver velocities of 2-3 mm/year (Figure 4b). In the Chile sliver, we find a mean Euler pole location of (132°E, 70°N), with an angular velocity of 0.05°/Ma (Figure 4c). This produces eastward velocities of 3-4 mm/year throughout the forearc.

Slip Deficit Rate
We compute the mean slip deficit rate in each sub-fault (averaged over all search iterations) as an estimate of the potential for that region to slip in an earthquake ( Figure 5). Although the mean model does not satisfy the physics imposed in each iteration of the search, it informs us about the magnitude of slip deficit accumulation. Generally, the areas of highest mean slip deficit rate overlap with the areas that have a high probability of being locked. The mean slip deficit rate is 60 mm/year in the Maule region, 55 mm/year in the Iquique/Pisagua section, and 50 mm/year in the Pisco/Arequipa region. Outside of the high-probability locked areas, the mean slip deficit rate is greater than zero because of proximity to the locked regions. This effect has a geographic footprint of several hundreds of km (Almeida et al., 2018;Herman et al., 2018); even in center of the 600 km long unlocked section offshore of northern Peru, the mean slip deficit rate is still 2 mm/year. We also show the probability distributions of slip deficit rates for individual sub-faults ( Figure 5, insets). Some of these distributions have two peaks: a peak at the relative plate velocity representing the probability that the sub-fault is locked, and another peak representing the slip deficit rate if the fault is unlocked. This is not simply a sampling artifact; the bimodal slip deficit rate distribution is a fundamental consequence of the pattern of slip deficit rates produced by the mechanical locking forward model. We expect that if we more finely mesh the subduction megathrust (e.g., for local studies), we could better resolve the probabilities of slip deficit rates between Figure 5. Mean megathrust slip deficit accumulation rate. The color of each sub-fault indicates its mean slip deficit rate, with recent large megathrust earthquake rupture zones plotted on top in white. The earthquake slip always occurs in regions accumulating slip deficit at~50% or more of the convergence rate. Black arrows show observed GNSS velocities with their 95% confidence interval. Blue vectors show the predicted velocities from the mean locking model plus the mean forearc sliver model. These generally fit the observed velocities to within 95% confidence. For selected sub-faults (i-viii), we show the probability distribution of the slip deficit accumulation rate, with the solid dark gray line indicating the convergence velocity at that location, the dashed red line indicating the mean slip deficit rate, and the dotted orange line indicating the median slip deficit rate. these two peaks, but our synthetic tests suggest the two peaks remain distinct. Sub-faults with low locking probabilities have single peaks below the relative plate velocity, indicating slip deficit accumulation primarily by proximity to a locked zone. Because of the smoothing intrinsic to the physics of the forward model, the peaks of the distributions all have standard deviations less than 10 mm/year (although there is no single, straightforward metric to characterize the bimodal sub-fault slip distributions). Even in unlocked sub-faults in areas of low resolution (e.g., near the trench), the slip deficit accumulation rate is 10-40% of the relative plate velocity.
The modeled velocities from the mean slip deficit rate distribution plus mean sliver motions fit the observed interseismic velocities very well in most locations ( Figure 5). Even better fitting individual models can be found in the search ensemble. The fit to observations is best in Chile, southern Peru, and northern Ecuador, where our predicted interseismic velocities are almost all within the 95% confidence interval of the observed vectors. The quality of fit (relative to the magnitude of the observations) is slightly worse in central Peru, adjacent to the section where locking is low, but generally the predicted velocities still fall within the 95% confidence level of the observations. Large, systematic misfits are only found in central Ecuador.
Here, the velocities change abruptly from <5 mm/year to the southeast to~20 mm/year to the east over a distance of~50 km across the Andes. We cannot explain this with a combination of plate interface backslip and rigid body rotations; a fault zone oblique to the subduction trench is likely necessary to explain these surface velocities (Alvarado et al., 2016).

Comparison to Great Earthquakes
The sections of the plate interface with high locking probabilities correlate well with the locations of the largest historical earthquakes (Hayes, 2017;Kelleher, 1972;Lomnitz, 2004; Figure 3). The main exception is offshore Chile from 23°S to 25°S, where we resolve high locking, but there is no record of great earthquakes. This section did host an Mw 8.0 earthquake in 1995 and an Mw 7.7 event in 2007 (Hayes, 2017). Similarly, most regions with lower locking probability are correlated with few to no great earthquakes in the historical record. These areas of low locking generally appear to act as persistent barriers to the large earthquakes that rupture the more highly locked regions. The only zone of low locking that historical ruptures have propagated across is from 32°S to 34°S. This correlation between interseismic locking and great earthquakes suggests that the distribution of these asperities is stable over multiple earthquake cycles and therefore that accurate interseismic locking distributions can be used to better anticipate likely great earthquake locations.
Most epicenters of Mw 8.0+ megathrust earthquakes, and many epicenters of Mw 7.5-8.0 events, are located in or near sub-faults with high locking probabilities. We interpret this to indicate that these earthquakes nucleate in the regions where stresses accumulate most quickly (Lay et al., 2012;Marone, 1998). However, several epicenters of smaller, Mw 7.5-8.0 events occur up to hundreds of km from zones of high locking probability. These smaller events are at the threshold of our resolution in this inversion, which suggests that they nucleate in asperities that are below the resolution limit of our plate-boundary-scale analysis.
Coseismic slip in great earthquakes does not always correlate with high locking, but it always occurs in regions accumulating slip deficit at a high rate ( Figure 5). The 2010 Maule and 2015 Illapel earthquakes ruptured areas with relatively high locking probability (and by definition high slip deficit). Slip in the 2000 Arequipa, 2007 Pisco, and Iquique earthquakes occurred in areas with a lower probability of being locked but accumulating slip deficit at~50% of the convergence rate. We interpret this to reflect that these great earthquakes ruptured into regions accumulating slip deficit adjacent to the locked asperity, which are made accessible to slip when the asperity unlocks. However, unless the rupture is over~200 km long, significant slip deficit will remain on the plate interface (including in the ruptured asperity; Herman et al., 2018), with the potential to produce single-asperity-rupturing Mw~8 events or multiple-asperityrupturing Mw~9 events.
These great megathrust earthquakes also generate tsunamis when large magnitude coseismic slip propagates to shallow depths and displaces the seafloor. It is important to constrain the slip deficit accumulated near the trench to evaluate these potential tsunami hazards; however, the resolution of kinematic slip 10.1029/2020GC009063 Geochemistry, Geophysics, Geosystems deficit inversions using onshore data is lowest near the trench (e.g., Gombert et al., 2018). The smoothing inherent to our mechanical forward model reduces the uncertainty in the magnitude of the slip deficit accumulation rate ( Figure 5). In South America, the slip deficit accumulation rate in sub-faults adjacent to the trench is at least 20% of the local convergence rate almost everywhere along the length of the subduction zone (Almeida et al., 2018;Herman et al., 2018). The only trench-adjacent sub-faults with lower slip deficit rates are located in the long stretches with no locked sections, such as northern Peru and central Chile. Therefore, any section of the subduction zone where some part of the megathrust is mechanically locked or slipping at a low rate must also accumulate slip deficit at ≥20% of the plate convergence rate near the trench, even if the near-trench area is not also mechanically locked.

Asperity Dimensions
Our inversion provides new constraints on the dimensions of the geological features responsible for locking the plate interface, given the assumption we have made that the plate interface is either completely locked (in asperities) or completely unlocked (outside of asperities). Under these binary locking conditions, no more than~30% of the megathrust surface area needs to be locked to produce the observed interseismic velocities, but at least 10% must be locked to produce great earthquakes. These results suggest that the locked zones are relatively small (length scales <50 km) and that the physical features on the slab or within the upper plate that are responsible for producing locking also have similarly small spatial scales. Therefore, plate-scale characteristics like slab curvature (Bletery et al., 2016) or regional characteristics like mean sediment thickness (Heuret et al., 2012) have wavelengths that are likely too long to account for this distribution of mechanical locking. In contrast, smaller geological heterogeneities, such as those associated with incoming seafloor roughness (Carena, 2011;van Rijsingen et al., 2018;Wang & Bilek, 2014) or geological heterogeneity within the upper plate (Béjar-Pizarro et al., 2013;Liu & Zhao, 2018) can be closer to dimensions of locking that we observe. Despite their small size, these locked zones produce sufficient slip deficit to account for the regular occurrence of great earthquakes in the subduction zone.
Although the interpretation of these small asperities is based on our analysis assuming a binary locking framework, we expect that we would find similarly small asperities if we had used a different plate interface rheology. This is because of the mechanical continuity of the system, which is valid irrespective of the specific megathrust rheology. The consequence of this continuity is that asperities-even if they are not completely locked-must produce a larger zone of slip deficit rate around them-even if the surrounding megathrust also has some shear resistance. This implies that if the basic mechanics of the system are included in the interseismic analysis, the inferred asperities will be smaller than their slip deficit rate footprint, as we have shown.
This modeling approach using a binary locking scheme provides an upper limit for the number and size of locked patches. If there were additional shear resistance in the regions we have defined as unlocked, this would add to the slip deficit rate that our model predicts, which would add to the landward surface velocities. As a result, fewer, smaller fully locked patches (specifically, having less of the plate interface area being locked or with near-zero slip) would be required to account for the magnitude of the interseismic velocities. The low resolution of the subduction plate interface provided by onshore data currently makes distinguishing a binary locked/unlocked interface from a more complex distribution of shear resistance difficult. However, in some locations where the upper plate extends farther over the trench (e.g., in Costa Rica; Feng et al., 2012) or with the addition of offshore data, it may become possible to further test these hypotheses.

Locking and Forearc Motion
Our results also suggest a reexamination of the importance of forearc sliver motions to explain the interseismic velocity field in the South America subduction zone. We still find that forearc motions added to the locking signal improve the fit to data, but these rigid velocities are~50% smaller than previously estimated (Figure 4). The place where forearc motions (independent of the megathrust locking) appear to be most important is northern Peru and southern Ecuador. Here, we find systematically worse fits to the interseismic velocities (although generally still within the 95% confidence interval).
In central and southern Peru, southeastward velocities have previously been interpreted to reflect trench-parallel forearc motion. However, in our inversion incorporating mechanical continuity into the 10.1029/2020GC009063 megathrust locking pattern, we find that locking in southern Peru drives the unlocked part of central Peru to the southeast across the transition from a locked to unlocked interface. It is beyond the scope of this study to determine whether this trench-parallel interseismic motion is recovered or continues to accumulate over many earthquake cycles (leading to long-term forearc translation). A similar effect may contribute to other forearc sliver motions, for example, the westward velocities in the Central America subduction zone (LaFemina et al., 2009).

Model Uncertainties
For this study, where we are taking a step toward a more physical inversion framework, we deliberately used a simple elastic half-space and imposed binary locking on the megathrust. Although we included a 20% model uncertainty in the inversion, it is likely that changing this forward model would have an effect on the resulting locking and slip deficit rate estimates. One way to address this issue is by varying model parameters during the search phase and subsequently evaluating them as part of the model vector M in p(M| D). However, this added complexity to the forward model comes at the cost of longer model run-times and (potentially many) more samples needed to map the parameter space. An alternative (also used in sensitivity tests for standard kinematic inversion analyses) is to use a different model parameterization for each complete search and compare the resulting distributions. The disadvantage of this approach is that these results may not be quantitatively (log-likelihood) comparable. Here, we discuss aspects of the locking model that could benefit from further investigation.

Plate Interface Rheology
Our binary locking treatment of the plate interface excludes models in which slip in asperities is nonzero and slip outside of asperities accumulates shear tractions. We have discussed the major associated issues in more detail throughout this study. Here, we simply reiterate that our results are an improvement on standard kinematic models and that we expect more complex interface rheologies will primarily change the details of asperity locations, not our fundamental conclusions.

Model Rheology
We anticipate that changing the model from an elastic half-space to a more complex rheology will have a comparable effect as in standard kinematic slip deficit rate inversions (e.g., spatial variations in elastic moduli: Masterlark, 2003; viscous deformation at depth: Li et al., 2015). Based on these previous analyses, varying the model rheology is expected to have a larger effect on the downdip extent of the locking and the slip deficit rate distributions than on the along-strike patterns. However, the megathrust slip smoothing effects caused by mechanical continuity will still occur in these more complex systems.

Plate Motion Model
We used the MORVEL plate motion model (DeMets et al., 2010) to calculate the relative velocity between the Nazca and South America plates. Other plate motion models (e.g., GEODVEL: Argus et al., 2010) change the input convergence rate in South America by up to 10 mm/year. A larger convergence rate would imply less of the interface needs to be locked to produce the same surface velocities.

Conclusions
Incorporating the physical interactions between nearby sections of the megathrust into a statistical interseismic inversion allows us to directly image locked areas and robustly estimate (reduced) uncertainties in slip deficit accumulation rates. Applying the inversion to the South America subduction zone, the imaged locked zones correlate with recent and historical earthquakes. Less than 30% of the plate interface is locked, and this small total locked area produces enough slip deficit to account for the slip during great earthquakes. Slip deficit accumulates at 20% or more of the relative plate motion everywhere on the subduction megathrust at <50 km depth, except offshore of northern Peru, which is >250 km from the nearest locked patch. Including rigid forearc sliver motions in the inversion improves the quality of the fit to data. However, the magnitude of the forearc motion is smaller than previously estimated.

Data Availability Statement
The software for computing fault slip and surface velocities and performing the search is available online (https://doi.org/10.5281/zenodo.3894137). The interseismic velocity data are available through the references listed in section 3.1, and the velocities are compiled in Data Set S1. The triangulated subduction interface is available in Data Set S2. Inversion outputs are available in Data Sets S3 and S4.