Interseismic deformation at subduction zones
investigated by 2D numerical modeling:
case study before the 2010 Maule earthquake

Marcelo Contreras1, 2, *Andrés Tassara1, 2, Muriel Gerbault4,
Rodolfo Araya5, Klaus Bataille3

1 Programa de Doctorado en Ciencias Geológicas, Facultad de Ciencias Químicas, Universidad de Concepción, Víctor Lamas 1290, Concepción, Chile.
marcont@udec.cl

2 Facultad de Medicina, Universidad San Sebastián, General Cruz 1257, Concepción, Chile.
mcontrerask@docenteuss.cl

3 Departamento de Ciencias de la Tierra, Facultad de Ciencias Químicas, Universidad de Concepción, Víctor Lamas 1290, Concepción, Chile.
andrestassara@udec.cl; bataille@udec.cl

4 Geosciences Environnement Toulouse, CNRS UMR5563, IRD, Université P. Sabatier, 14 av. Edouard Belin, 31400 Toulouse, France.
gerbault@get.obs-mip.fr

5 Departamento de Ingenería Matemática, Universidad de Concepción, Víctor Lamas 1290, Concepción, Chile.
rodolfo.araya@udec.cl

* Corresponding author: andrestassara@udec.cl

Fig. 1. GPS-derived surface velocity vectors across a trench-perpendicular transect at Arauco Peninsula (aprox 37.3º S) in south-central Chile (topography/bathymetry of the area is color-coded). Vectors where calculated by Ruegg et al. (2009, blue dots) and Moreno et al. (2010, yellow dots) from campaign measurements done before the Mw8.8 2010 Maule earthquake. Contours of 3, 6 and 9 m of co-seismic slip of the 2010 event are shown in red (Moreno et al., 2012). The convergence velocity vector is shown as a white arrow.

 

The Back Slip Model (Savage, 1983) assumes that the megathrust is a finite fault embedded into an elastic media and simulates the interseismic period by a downwards slip of the upper plate’s base along this fault, simulating the coupling of the forearc wedge with the descending oceanic slab. The kinematic downdip movement of this fault produces an elastic deformation that translates into a surface displacement, which compared to geodetic observations allows to constrain the backslip velocity and therefore the degree of fault locking. Although the BSM approach is an oversimplification of the geometrical and rheological configuration of a subduction zone, which cannot be used to gain direct understanding on the mechanical interaction between plates leading to earthquakes, its simplicity has allowed successful applications to geodetic data analysis and seismic hazard assessment in a number of subduction zones worldwide (e.g., Prawirodirdjo et al., 1997; Mazzotti et al., 2000; Khazaradze and Klotz, 2003; Chlieh et al., 2008; Moreno et al., 2010).

In the attempt to construct a conceptual model that could describe the subduction process in a more realistic way, and yet preserve the simplicity of the kinematic link between geodetic observations and plate locking, Kanda and Simons (2010) proposed the Elastic Subducting Plate Model (ESPM). This model is based on a previous conceptualization of Zhao and Takemoto (2000) and considers a subducting elastic slab of finite thickness coupled to an elastic upper plate along the seismogenic megathrust zone that moves with respect to the surrounding asthenospheric mantle. The ESPM also includes the kinematic effect of slab flexure near the trench axis by considering a decreasing velocity gradient along the slab’s radius of curvature. Plate flexure and its associated bending stresses have been shown to be a fundamental process that connects the dynamic force balance of tectonic plates with the seismic cycle (Conrad et al., 2004; Capitanio and Morra, 2012; Schellart, 2009; Buffett and Rowley, 2006). As demonstrated by Kanda and Simons (2010), the ESPM mimics the deformation field caused by the BSM when the thickness of the subducting plate tends to zero. However, if at least a portion of flexural stresses is not continuously released in the shallow portion of the subduction zone, then the predicted surface velocities of both models can differ significantly. This conclusion implies that the interpretation of geodetic observations within the framework of more (ESPM) or less (BSM) realistic approaches for assessing the degree of interseismic plate locking (and the maturity of a seismic cycle), remains to be developed together with an understanding of its mechanical implications.

With the main aim of analyzing the mechanical implications of the ESPM on the kinematics and mechanics of plate subduction and its influence on megathrust seismogenesis, we present here a series of 2D numerical models that range from a simple kinematically-constrained elastic rheology at short time scales (a year) to a complex elasto-visco-plastic rheology at large time-scales (millions of years). First we present our own version of an ESPM as implemented using the finite element method (FEM), with validation against analytical solutions of elastic dislocations. Our method further allows to account for realistic geometries of a subduction zone and the introduction of gradients in the velocity distribution at the base of the subducted plate. In section 3, we apply our modeling method to GPS-derived surface velocities obtained years to months before the Mw8.8 2010 Maule earthquake (Moreno et al., 2010; Ruegg et al., 2009) along a nearby trench-orthogonal transect in south-central Chile (Fig. 1). At this latitude, the proximity of the coastline to the trench (ca. 100 km) allows appreciating surface deformation by land-based GPS stations. The exercise shows that in order to fit the geodetic observations, the velocity of the subducting slab at the base of the plate near the trench axis must be significantly reduced with respect to the far-field plate velocity. In section 4, we use other numerical modeling approaches (ADELI: Hassani et al., 1997, and PAROVOZ: Gerbault et al., 2009) to investigate the meanings of this result in terms of the mechanical interaction between the subducting slab and the underlying astenospheric mantle. We compare the instantaneous stress field generated by kinematic and geometric configurations similar to those applied to our initial elastic model, with larger time-scale subduction models. Such long-term models show that when bending builds up naturally as a consequence of continuous subduction of a competent plate, the high deviatoric stresses that accumulate at the bent base of the slab are located in a position equivalent to the region where a reduction of slab velocity was required in the EPSM to fit the geodetic data. We finally discuss the implications of these results in terms of the relationship between the state of stress in the long-term in a subduction zone, and short-term seismogenic energy release.

2. Finite element modeling of an elastic subduction zone during the interseismic phase

In order to simulate the elastic deformation produced by interseismic plate locking at subduction zones we use the finite-element method (FEM) to develop a 2D model in which displacement along a fault plane is implemented based on a split node technique (Melosh and Raefsky, 1981). We start with a conceptual model of an idealized subduction zone whose general configuration setup is shown in figure 2. This setup can indistinctly reproduce the BSM or the ESPM. For both models the subduction geometry is characterized by two parameters: Θ is the dip angle of the slab below the locked seismogenic zone and DDL is the depth to the Downdip Limit at which both plates are locked. For sake of simplicity as a first step, the DDL also corresponds to the limit between the upper plate and the underlying asthenospheric mantle (Fig. 2). Different domains (upper plate, subducting plate, asthenospheric mantle) can be recognized with eventually different values of the Young Modulus. The upper surface of the model is assumed to be stress free (i.e., it freely deforms) whereas along the lateral and bottom boundaries (extending down to 300 km in the mantle) no displacement is allowed (anchored boundaries). For the ESPM setup, a subducting plate of thickness H is defined and its motion is imposed as a boundary condition at its left-hand border. On the opposite edge of the slab at 300 km depth, free motion is allowed. In order to simulate the kinematic conditions of the interseismic phase, we impose a uniform slip along the boundaries of the elastic plate with the asthenospheric mantle (red lines in figure 2), and assume a 100%-coupled seismogenic zone (continuous deformation and stress in between the upper and lower plate from the top surface to the DDL depth). In contrast, in the simpler geometric setup corresponding to the BSM, a homogeneous slip is imposed along the seismogenic zone (blue line in figure 2) and all other boundaries are fixed.

 

Fig. 2. Model setup for the FEM experiments of subduction models. The figure define dimensions of the modeled space, boundary conditions, geometric parameters and different material domains (W1-3). For the ESPM configuration a slip SESPM is imposed along the red slab boundaries and the interplate boundary (blue) is kept 100% locked, whereas for the BSM a backslip SBSM is imposed along the interplate boundary. See text for details.

 

3. Testing subduction models and FEM validation

3.1. BSM, an end-member case of ESPM

In order to first validate our model setup and to reproduce the analysis of Kanda and Simons (2010) with our finite element tool, we assign the values of 40º and 40 km for Θ and DDL respectively, and create two different models as shown in figure 3; ESPM (Fig. 3a) and BSM (Fig. 3b). For the ESPM we assign 1 m of slip along the moving boundaries of the slab (red lines in figure 2) whereas for the BSM the same displacement is applied along the plate interface. In the limit when the slab thickness (H) of a ESPM tends to zero, the resulting deformation field should approach the one generated by a BSM, since the imposed  displacements cancel each other except below the locked seismogenic zone, where slip at the base of the slab subsists and mimics backslip as for a BSM (Kanda an Simons, 2010). The displacement field resulting from the ESPM with a plate of thickness (H) of 5 km (approaching the case of a zero-thickness BSM) is shown in figure 3c and can be compared with the field generated by the BSM in figure 3d. We see that both models generate an almost identical deformation within the volume of the upper plate: the norm of the displacement field decreases from a maxima at the trench to about 50% at location x=50 km, and to about 25% at location x=80 km.

 

Fig. 3. Comparison of BSM and ESPM with thickness H=5 km. Panels A and B show the model setup for a ESPM and BSM respectively. Lower panels depicts the calculated displacement field after applying 1 m slip to the boundaries of the slab in a ESPM (C) and the interplate boundary for a BSM (D).

 

The similarity between both models can also be appreciated in figure 4, which displays horizontal and vertical components of the surface displacement obtained from a BSM and from a ESPM with different plate thicknesses H between 5 and 50 km. Landward from the trench the solutions of the BSM and the 5 km thick slab ESPM are almost identical, while they both differ significantly with increasing slab thickness.

 

Fig. 4. Comparison of horizontal (A) and vertical (B) surface displacement for BSM and ESPM with changing slab thickness H between 5 and 50 km.

 

Thus, the end-member ESPM with a nearly zero slab thickness can be efficiently modeled by a simple BSM, as shown by Kanda and Simons (2010) with their analytical formulation. This exercise serves as an independent validation of our numerical model.

3.2. Modeling interseismic deformation by BSM across the Arauco Peninsula before the 2010 Maule earthquake.

With the intention of offering an independent validation of our FEM approach, now we reproduce the BSM model proposed by Ruegg et al. (2009) to fit the GPS interseismic velocities measured before the Mw8.8 Maule 2010 earthquake. For this we use velocity vectors reported by Ruegg et al. (2009) and Moreno et al. (2010) that were acquired during the decade preceding the Maule earthquake. These measurements are roughly aligned in a profile across the center of the Arauco Peninsula from the coastline extending east to the Andean volcanic arc near the Chile-Argentina boundary (Fig. 1). Ruegg et al. (2009) defined their BSM with a planar fault geometry characterized by the following parameters: Strike=0°, rake=−90°, DIP=17°, Up-Dip Limit UDL=0 km and width of the fault W=180 km. We constructed a 2D finite element mesh with these geometric parameters and imposed a backslip velocity of 65.5 mm/yr equal to the trench-normal convergence velocity.

Figure 5 compares the surface displacement predicted by the BSM analytical solution, which assumes the classical formulation of Okada (1985), with that computed with our ESPM FEM approach. It can be observed that both solutions are almost identical, mostly for the horizontal component, and that they fit the GPS data equally well.

 

Fig. 5. Modeled horizontal (A) and vertical (B) velocity and GPS observations (points with error bars) with the analytical solution of Okada (1985) and our FEM solution reproducing the same BSM parameters used by Ruegg et al. (2009).

 

4. Exploring parameters of the ESPM for interseismic deformation before the 2010 Maule earthquake

As recognized by Kanda and Simons (2010) and many others, a BSM approximation with a planar fault is an incorrect physical model of the process of plate subduction, and therefore we attempted to model the interseismic GPS data across the Arauco Peninsula with our ESPM approach. We first constructed a finite element mesh that interpolates representative points of the local topography and bathymetry, as well as a realistic geometry of the slab upper surface (derived from the geophysically-constrained 3D model of Tassara and Echaurren, 2012). The upper plate has a constant thickness of 40 km.

With this FEM mesh (Fig. 6), we first study the influence of the slab plate thickness (which we assume constant across the model), the downdip limit (DDL) of the seismogenic zone, and variations of the elastic parameters between different domains. We assume the same boundary conditions as those defined for the conceptual ESPM (section 2) and apply a displacement of the oceanic plate from the left-hand side equal to 65 mm/yr, simulating its trench-normal motion during one year. The locked zone between the Nazca and South American plates is considered fully locked in-between the trench and its downdip limit (DDL), which can vary between 26 and 60 km depth in our simulations. Additionally, we considered variations of the elastic parameters of each domain Ω (see figure 2) in order to quantify their effect on surface displacements. Following the values given by Turcotte and Schubert (2002) we consider these domains to be dominantly composed by granite, granodiorite and/or peridotite with a Young modulus E 0.5, 1.0 and/or 1.5x1011 Pa, respectively.

 

Fig. 6. Model setup for a ESPM across the Arauco Peninsula near. The slab geometry was constructed according to the 3D geophysical model of Tassara and Echaurren (2012). The figure shows applied boundary conditions and defines main parameters of the model (see text for details).

 

With this setting, we try to reproduce the interseismic horizontal surface displacement observed with the GPS data. In a first attempt, we express a potential reduction of velocity at the bottom of the slab (v_bcs) with respect to the applied convergence velocity (v_c), which results from plate bending at the trench (Kanda and Simons, 2010). This velocity can be evaluated geometrically and it is non-uniform below the seismogenic zone as defined by v_bcs =γ( (1-H/R)/(1-H/2R)), with H the thickness of the slab and R its local radius of curvature. γ is a parameter between 0 and 1 that allows us to express this reduced velocity at the base of the slab.

Figure 7a (horizontal component) and 7b (vertical component) show that, for a granodioritic upper plate with peridotitic slab and mantle, with γ=1 and at a given DDL, varying the slab thickness H between 10 and 50 km has little effect on horizontal surface displacements, which appear significantly smaller than the uncertainty in GPS measurements. Variations in DDL in turn have a much larger effect. A DDL at 27 km depth produces deformation concentrated towards the trench, whereas a DDL near 60 km depth shifts the surface deformation towards the volcanic arc with an excess of deformation in the forearc near the coast. The best-fitting combination of DDL and H under the described conditions is accomplished for values of 42 km and 30 km respectively (Fig. 7a   and b), although we note that none of the models is really able to give an acceptable fit to the GPS data. This difficulty in reproducing the observed interseismic surface displacements with an ESPM approach remains independent from the choice of elastic parameters for each model domain.

 

Fig. 7. ESPM Model results compared with GPS velocities along the Arauco Peninsula. Left-hand and right-hand panels show horizontal and vertical velocities respectively. With respect to uncertainties of GPS observations, surface displacement produced by the ESPM with realistic subduction geometry is almost insensitive to slab thickness H (as shown in A and B). Young Modulus E of the upper plate has a minor effect (as shown in C and D). Displacement is mostly controlled by the Downdip Limit (DDL) of interplate locking. Even our best fitting model (H=30 km, Granodiorite, DDL=42 km) is unable to reproduce the observed GPS velocities.

 

We then show that this misfit can be corrected by introducing a reduction in the velocity at the bottom of the subducting plate below the coupling zone. Now we consider the same configuration as described in the previous section, but we modify the velocity distribution along the bottom boundary of the subducting plate, as illustrated with arrows of different colors in figure 8. One configuration considers that this basal velocity  is constantly reduced by a factor γ=0.1 with respect to. When we apply this velocity reduction to the entire curved bottom of the slab, the resulting surface deformation is largely reduced well below the level given by GPS measurements (black curve in figure 9). However, if we apply this reduced velocity along the bottom of the slab from the trench only down to the DDL depth, then the match to GPS observations is strongly improved (blue curve in figure 9).

 

Fig. 8. Model setup for an ESPM with reduced velocity at the base of the slab due to bending. This is similar to figure 7 but here we introduce the possibility of reducing the displacement velocity at the base of the slab compared with the rest of plate as shown by arrows of different colors.

 

 

Fig. 9. Horizontal surface velocity resulting from a ESPM model that incorporates velocity reduction at the base of the slab. Our best ESPM notably improves the fitting to observed GPS rates if the velocity at the base of the slab is strongly reduced to 10% of the convergence velocity.

 

5. Mechanical interaction between subducted slab and mantle

In searching for a dynamic explanation to the necessity of a localized velocity reduction at the base of the subducted slab, we explore now with complementary numerical models the mechanical interaction between a subducting plate and its underlying mantle. Beginning with boundary conditions similar to those imposed in the previous section, e.g., thought to be appropriate for slab-mantle interaction at the short-scale, we progressively evolve towards boundary conditions appropriate for long-term slab-mantle interaction. We thus discuss step by step the implications of some assumptions, and how our reasoning brings insight on interseismic seismogenic behavior.

5.1. Stress field resulting from kinematic conditions imposed at the slab-mantle interface

To explore how the kinematic conditions applied in the previous ESPM models of sections 3 and 4 actually translate into mechanical interactions between the slab and the underlying asthenospheric mantle, we first use the FEM-based code ADELI (Hassani et al., 1997; Cattin et al., 1997; Provost and Chery, 2006), which has already been used to study various mechanical aspects of the seismic cycle as well as long-term subduction, and presents the advantage of handling pre-defined contact zones enabling relative motion between different space domains. In the models below we consider synthetic geometries, since the aim here is to identify the theoretical influence of kinematic conditions along the slab-mantle interface on the stress field.

We define a first numerical setup (Adeli Model AM-O, Fig. 10a) in which the underlying mantle is considered as an inviscid fluid. Therefore it does not need to be included in the physical mantle domain (as in e.g., Hassani et al., 1997). This setup considers identical kinematic conditions along the plates to those assumed in the previous sections, with applied tangential velocity at the base of the subducting plate, and a fixed base of the overriding plate. The interaction with the overriding plate is simulated as a locked contact zone of high friction, simulating interseismic coupling (e.g., Cattin et al., 1997). The gravity force is included (to the difference of ESPM models of previous sections). Assumed to behave as an inviscid fluid. Note that the deviatoric shear stresses that develops inside both plates are of the order of 100 and 200 MPa (turquoise to green, identical in all three cases. Figure 10a displays the resulting internal stress field inside the slab (of the order of 100 MPa), resulting from the imposed velocity discontinuities at the base of the two plates. Two zones of discontinuous stresses develop in the upper plate in the forearc area (magnitude of about 200 MPa).

 

Fig. 10. Three models set with code ADELI, in order to illustrate the mechanical interaction between  the slab and the mantle, Each three cases impose a tangential velocity  along the base of the subducting plate, and no motion at the base of the overriding plate (as in the previous section). The deviatoric shear stresses that develop inside both plates reach ~100-200 MPa (turquoise to green), and are identical in all three cases. However, all 3 cases assume different conditions for the underlying mantle. A. Model AM-0 assumes that the mantle behaves as an inviscid fluid; B. AM-1 includes a mantle with temperature dependent viscosity (on average, of the order of 10² Pa.s): its lateral and bottom borders are set free-slip (similar to models in section 3). C. AM-2 includes a mantle of similar rheology but its boundaries react to Archimedes’ restoring force (as if embedded in an infinite half space). The large stress  field generated in the mantle results from the kinematic conditions applied on to the base of the plates, thus showing that they are out of mechanical balance.

 

If now we incorporate the mantle in our models, a conceptual issue appears: when one accounts for the flexed geometry of a slab inside the mantle, one must also think of the mechanical effects of this flexure. This mechanical effect is usually thought of as tension on the top of the bent slab and compression in its inner corner. However, the mantle underneath is also loaded, and deformed in association to this slab bending. Considering a plate that subducts from rigth down to the left, the surrounding mantle thus receives a clockwise torque from above, and flows accordingly to balance this torque (for maintenance of the mechanical equilibrium of the entire medium). Even though the mantle domain is visco-elastic and likely dissipates high stresses via low viscosity flow, mechanical balance must still be assessed. Consequently, how do the kinematical conditions imposed on the base of a subducting plate as defined in sections 3 and 4, fall within this frame-view of slab-mantle mechanical interaction?

When incorporating the mantle in numerical models of subduction, one must define the boundary conditions acting on its boundaries, either with kinematic or stress conditions. We test both options in the subsequent two models AM-1 and AM-2 (Fig. 10b and c). First, in AM-1 we prescribe free-slip boundary conditions on the lateral borders and bottom base of the mantle (as in models of previous sections). Upon applying a displacement of 1 meter of the subducting plate in this model (during 100 yrs), the resulting state of stress inside both subducting and overriding lithospheres is very similar to AM-0. However, to the difference of AM-0, now large deviatoric stresses also appear in the mantle underlying the overriding plate (much greater than 200 MPa, Fig 10b).

In the second model AM-2 (Fig. 10c), Archimede’s lithostatic stress reaction force is applied along the boundaries of  the modeled mantle domain (e.g., simulating that the modeled domain is embedded in an infinite half space, as in e.g., Gerbault et al., 2009 or Burov et al., 2005). The resulting state of stress displays even larger deviatoric stresses both in the eastern and western parts of the modeled mantle. This stress field results from the kinematic constraints that we imposed onto the base of the slab and the overriding plate, and which exert a flexural torque on the mantle. The resulting stress fields in both models AM-1 and AM-2 show that the flexure of natural subduction zones are actually not balanced in a simple way (e.g., not only geometrical) with the imposition of kinematic conditions on the slab-mantle interface. In other words, the application of the specific kinematic conditions proper to the ESPM produces large stresses in the underlying mantle. We show here with these three models (AM-0,1,2) that mechanical equilibrium is a significant issue, that is not resolved along the plate interface alone, and despite the fact that the stress distribution inside the slab remain similar in all three cases.

Let us now show another model AM-3 (Fig. 11a), in which we allow the slab to slide freely above the underlying mantle (thus no imposed velocity and no resisting friction). Now gravity is turned off again, and only boundary velocities are applied at the far-field right-end and at the bottom intersection of the slab with the modeled domain. A typical flexural pattern develops within 100 yr with bending stresses of the order of 200 MPa inside the slab (Fig. 11a). In the surrounding mantle in turn, no high stresses develop, in contrast to models AM-0,1,2.  One cannot avoid recognizing the existence of slab bending stresses within the frame of plate tectonics at the long-term scale of millions years. So now we are faced with the problem of deciding if the boundary conditions at the slab-mantle interface associated to slab bending at the long-term, also correctly simulate the same slab-mantle interface at the scale of the seismic cycle. In the quest for assessing the appropriate kinematic conditions that control interseismic deformation, one may (or may not) assume that slab bending is negligible when considering relative shear motion (from one earthquake to the other). In order to further clarify the importance of long term deformation and stress pattern, we present in the next section a last model that simulates the long-term flexural behavior of a subducting plate.

 

 

Fig. 11. Stress distribution inside an elastic plate when considering long-term subduction. A. Model AM-3, produced with code Adeli, with no gravity field, and with freely sliding plate-mantle boundaries (along which neither tangential nor normal motion are imposed). Within ~300 kyr of applied convergence from the right-hand side of the model (as in AM-2), flexural deviatoric stresses up to 260 MPa (in red) develop inside the elastic plate; B. Model LM built with code “Parovoz” (modified from Gerbault et al., 2009). Here gravity is included and the mantle and overriding plate both behave according to thermally dependent elasto-visco-plastic laws. Therefore their boundary is self-consistently defined and evolves throughout the time duration of the model. Within 4 Myr of applied far-field convergence, flexural stresses of the order of 1,000 MPa develop inside the elastic portion of the slab. Black and red lines display the orientation of principal tensile and compressional stress, respectively, typical of a flexural pattern; C. Profiles along the base of the slab in LM,  determined by the isotherms around 600 ºC: depth of profile (upper profile),  accumulated shear strain (middle), and shear stress (second invariant, bottom). At the slab’s corner around X=-50 km stresses exceed 1,000 MPa, whereas the shear strain remains low (red rectangle domain). Such a stress concentration is equivalent to a frictional resistance that explains the need to reduce artificially the eastward velocity in this bend in the FEM-based kinematic ESPM.

 

5.2. Long-term flexural stresses

Here we illustrate how flexural stresses develop in the long-term in a subduction zone. Therefore we display model LM, which is built with the finite-differences code Parovoz (Poliakov and Podladchikov, 1992), and assumes conditions similar to those defined by Gerbault et al. (2009), Dorbath et al. (2008) and Sobolev and Babeyko (2005) when studying the influence of rheology on deformation along the Chilean margin over several millions years. This model considers subduction over a time-scale of 4 Myrs and assumes self-consistently defined temperature-dependent elasto-visco-plasic rheology. Therefore the slab-mantle boundary is also self-consistently defined and consists in a continuous transition from competent (effectively elastic) plates to low viscosity asthenosphere. The setup is defined with the synthetic geometries of a bent slab which is pushed from its right-end side, eastward at the velocity of 4 cm/yr, while the overriding plate is remained fixed (free-slip far-field left-end of the model). At the base of the model domain the slab is also pulled downwards in the direction parallel to its inclination, at 4 cm/yr. As a result, flexural stresses build up and equilibrate within 4 Myr (Fig. 11b), using the flexure of the resistant slab to equilibrate the applied boundary conditions. Figure 11b displays the stress field via its 2^nd invariant. Highest values occur where slab curvature changes: these stresses approach 1 GPa at the bottom corner of the flexed lithosphere (sampling of the displayed values was made at depths corresponding to the 590-620 ºC isotherm interval, Fig. 11c). Thus we illustrate here that although the bent corner of the slab does not undergo much internal deformation (red rectangle domain in figure 11c), it actually stores large bending stresses capable of exceeding the brittle-yield strength.

Two side comments rise from these models. First, the large stresses that develop in this model result from our assumption of purely visco-elastic behaviour in this case, in order to remain consistent with the other models of the present study. If rheologies were assumed fully elasto-visco-plastic as in e.g., Gerbault et al. (2009), then the stresses would remain bounded by the yield strength and would not exceed ~400 MPa. However the strain and stress field still both develop a similar pattern to that displayed here. Second, if we compare more precisely the location of these large bending stresses with respect to the position of the trench, one notes that in model LM (Fig. 11b) they are offset by about 50 km further to the left of the trench and are generally more inclined with respect to the horizontal, in comparison with the AM models. This difference is due to the different designed initial geometries (plate length and curvature between Adeli and Parovoz), but does not affect the mechanical implications of our study.

In comparing the AM0,1,2 models with the AM3 and LM models we see that the imposition of specific velocities along the base of the slab, despite accounting for its realistic motion, actually impedes the development of flexural stresses inside the slab. That is the reason why high stresses develop in the mantle in models AM0,1,2 (Fig. 10), that counter balance the imposed velocities. In contrast in the long-term model LM, although the base of the bent slab is inherently kinematically locked with the low-viscosity underlying mantle, large stresses rather develop inside the slab itself since it is the solicited competent body under flexure.

Slab bending is actually involved in the energy balance that makes the dynamics of a subduction zone (Buffett, 2006; Ribe, 2010 Capitanio and Morra, 2012; Fourel et al., 2014. Bending stresses are intricated with the other driving and resisting forces and the energy balance does not provide similar quantities depending on the assumptions of the different authors. For instance downwards in the asthenosphere, bending energy is involved in the low viscosity asthenospheric flow that accompanies subduction a the scale of the entire upper mantle (e.g., viscous drag, Hager and O’Connell, 1978; Piromallo et al., 2006; Li et al., 2014). Downwards along the axis and inside the slab, bending is at least partly restored by slab unbending, as seen in figure 11b with a reversed orientation of the tensile/compressional stresses at ~100 km depth (e.g., Buffett, 2006; Hassani et al., 1997; Čížková et al., 2007; Faccenda and Mancktelow, 2010). Upwards above the slab, this energy must also, very likely, be partly released in some manner. Since pressure decreases, large stresses are no longer sustained (brittle failure occurs). However deformation is eased, and we can imagine that “stress” is transformed into “deformation”. Let us thus propose a simple back of the envelop exercise. Let us assume energy conservation between the bottom and top of a   lithospheric column set right in the bent corner of the slab: if the shear stresses are about 10 times greater in the deep domains with respect to the top domains at the slab’s bend (100 to 400 MPa as opposed to 10 to 40 MPa along the seismogenic zone, e.g., Lamb, 2006; Tassara, 2010), then energy conservation (expressed as the product of stress and strain rate) implies that the deformation that is to be released at the top of the slab should be 10 times greater than that along the base of the slab. Coincidentally, in order to fit the surface (top) interseismic deformation in the ESPM model of section 4, we had to reduce the velocity applied along the base of the slab by a similar factor of 10 (in the bent corner). This reasoning is consistent with the concept that stored (shear) stress at the bent base of the slab transfers towards its upper portion into (shear) strain. Unfortunately, deciphering the amplitude of this energy transfer with respect to other processes requires assessing non-linear temporally and spatially scaled parameters (such as speed of stress and strain guides in the deforming pores of the bent lithosphere, with respect to the above mentioned other modes of energy transfer), which falls out of the ambitions of the present paper.

6. Conclusions

The well known Back Slip Model (BSM; Savage, 1983) solution of an elastic dislocation in a semi-infinite half space (Okada, 1985), is widely used to predict interseismic deformation along a given fault zone and has been efficiently applied to model GPS observations before the 2010 Maule earthquake (Ruegg et al., 2009; Moreno et al., 2010). However, the BSM can not give insights into the physical and mechanical processes properly taking place at subduction zones because it exclude the presence of the subducting plate. Consequently we have developed a two-dimensional finite element approach simulating interseismic elastic deformation at subduction zones, which accounts for realistic interplate geometries and allows to apply kinematical constraints along the base of the subducting lithosphere. Similarly to Kanda and Simons (2010), our Elastic Plate Subducting Model (EPSM) confirms that the BSM solution is identical to the solution of an EPSM when plate thickness is negligible. Our EPSM geometrical setup was then constrained by independent geophysical data, and fitted with the observed GPS surface rates across the Arauco Peninsula. A satisfactory fit was obtained only when we reduced the basal velocity at the corner of the slab below the trench to 10% of the total convergence rate after correction of the geometrical effect of slab bending.

This velocity reduction likely represents the kinematic consequence of stress dissipation caused by plate bending near the trench axis, as suggested by complementary mechanical models that explore the slab-mantle mechanical interaction. These complementary models reveal that purely kinematic conditions applied to the slab impede the build-up of flexural stresses inside the slab and thus generate out of balance stresses in the adjacent mantle. The inherent flexural stresses associated with a naturally bending slab are shown to concentrate indeed at the basal corner of the subducting plate. Such high bending stresses thus explain the artificial reduction in basal velocity that we had to apply to our EPSM in order to fit GPS interseismic observations, and allows us to state that the slab-mantle boundary acts as an efficient frictional interface in the bent corner of the slab. This statement is consistent with the argumentation of Buffett (2006) that bending of the subducting plate exerts a basal friction that opposes to the motion of the downgoing plate, and that the driving forces of subduction include slab flexure (Forsyth and Uyeda, 1975), which is partially balanced by local resistance (Conrad et al., 2004).

Our study shows that in order to retrieve appropriately the mechanical processes in play during interseismic deformation, numerical models should account for the effects of plate flexure, thus pointing to the perspective of integrating both short-term and long-term approaches (e.g., Van Dinther et al., 2013). Meanwhile, the careful application of locally reduced basal velocities in an ESPM approach offers a more realistic treatment of GPS interseismic surface displacements than the conventional BSM approach.

Acknowledgements
We thank fruitful discussions with R. Hassani and M. Moreno. This research was funded by Fondecyt Projects 1101034, 1151175, ECOS/Conicyt C13U03, DFG/Conicyt PCCI130039, DAAD/Conicyt 2058-2010. The detailed review of F. Aron allows us to greatly improve the clarity and quality of this contribution.

References

Avouac, J.-P. 2015. From Geodetic Imaging of Seismic and Aseismic Fault Slip to Dynamic Modelling of the Seismic Cycle. Annual Review of Earth and Planetary Sciences 43: 233-271.

Buffett, B.A. 2006. Plate force due to bending at subduction zones. Journal of Geophysical Research: Solid Earth 111 (B9): B09405.

Buffett, B.A.; Rowley, D.B. 2006. Plate bending at subduction zones: Consequences for the direction of plate motions. Earth and Planetary Science Letters 245: 359-364.

Buffett, B.A.; Becker, T.W. 2012. Bending stress and dissipation in subducted lithosphere. Journal of Geophysical Research: Solid Earth 117: B05413. doi: 10.1029/2012JB009205.

Burov, E.; Guillou-Frottier, L. 2005. The plume head-continental lithosphere interaction using a tectonically realistic formulation for the lithosphere. Geophysical Journal International 161 (2): 469-490.

Capitanio, F. A.; Morra, G. 2012. The bending mechanics in a dynamic subduction system: Constraints from numerical modelling and global compilation analysis. Tectonophysics 522: 224-234.

Cattin, R.; Lyon-Caen, H.; Chéry, J. 1997. Quantification of interplate coupling in subduction zones and forearc topography. Geophysical Research Letters 24 (13): 1563-1566.

Chlieh, M.; Avouac, J.P.; Sieh, K.; Natawidjaja, D.H.; Galetzka, J. 2008. Heterogeneous coupling of the Sumatran megathrust constrained by geodetic and paleogeodetic measurements. Journal of Geophysical Research: Solid Earth 113 (B5): 2156-2202.

Čížková, H.; van Hunen, J.; van den Berg, A. 2007. Stress distribution within subducting slabs and their deformation in the transition zone. Physics of Earth and Planetary Interiors 161 (3-4): 202-214. doi: 10.1016/j.pepi.2007.02.002.

Conrad, C.P.; Bilek, S.; Lithgow-Bertelloni, C. 2004. Great earthquakes and slab pull: Interaction between seismic coupling and plate-slab coupling. Earth and Planetary Science Letters 218: 109-122.

Dorbath, C.; Gerbault, M.; Carlier, G.; Guiraud, M. 2008. Double seismic zone of the Nazca plate in northern Chile: High-resolution velocity structure, petrological implications, and thermomechanical modeling. Geochemistry, Geophysics, Geosystems 9 (7): Q07006. doi: 10.1029/2008GC002020.

Faccenda, M.; Mancktelow, N.S. 2010. Fluid flow during unbending: Implications for slab hydration, intermediate-depth earthquakes and deep fluid subduction. Tectonophysics 149 (1-2): 149-154. doi: 10.1016/j.tecto.2010.08.002.

Forsyth, D.W.; Uyeda, S. 1975. On the relative importance of the driving forces of plate motion. Geophysical Journal International 43 (1): 163-200.

Fourel, L.; Goes, S.; Morra, G. 2014. The role of elasticity in slab bending. Geochemistry, Geophysics, Geosystems 15 (11): 4507-4525.

Gerbault, M.; Cembrano, J.; Mpodozis, C.; Farías, M.; Pardo, M. 2009. Continental margin deformation along the Andean subduction zone: Thermo-mechanical models. Physics of the Earth and Planetary Interiors 177: 180-205.

Hager, B.H.; O’connell, R.J. 1978. Subduction zone dip angles and flow driven by plate motion. Tectonophysics 50 (2): 111-133.

Hashimoto, M.; Jackson, D. 1993. Plate tectonics and crustal deformation around the Japanese Islands.Journal of Geophysical Research: Solid Earth 98 (B9): 16149-16166.

Hassani, R.; Jongmans, D.; Chéry, J. 1997. Study of plate deformation and stress in subduction processes using two-dimensional numerical models. Journal of Geophysical Research: Solid Earth 102 (B8): 17951-17965.

Kanda, R.V.S.; Simons, M. 2010. An elastic plate model for interseismic deformation in subduction zones. Journal of Geophysical Research: Solid Earth 115 (B3): B03405.

Khazaradze, G.; Klotz, J. 2003. Short-and long-term effects of GPS measured crustal deformation rates along the south central Andes. Journal of Geophysical Research: Solid Earth 108 (B6): 2289. doi: 10.1029/2002JB001879.

Lamb, S. 2006. Shear stresses on megathrusts: Implications for mountain building behind subduction zones. Journal of Geophysical Research: Solid Earth 111 (B7): B07401. doi: 10.1029/2005JB003916.

Li, Z.H.; Di Leo, J.F.; Ribe, N.M. 2014. Subduction‐induced mantle flow, finite strain, and seismic anisotropy: Numerical modeling. Journal of Geophysical Research: Solid Earth 119 (6): 5052-5076.

Mazzotti, S.; Le Pichon, X.; Henry, P.; Miyazaki, S.-I. 2000. Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: An analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS. Journal of Geophysical Research: Solid Earth 105 (B6): 13159-13177.

Melosh, H.J.; Raefsky, A. 1981. A simple and efficient method for introducing faults intofinite element computations. Bulletin of the Seismological Society of America 71 (5): 1391-1400.

Moreno, M.; Rosenau, M.; Oncken, O. 2010. 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. Nature 467 (7312): 198-202.

Moreno, M.; Melnick, D.; Rosenau, M.; Báez, J.; Klotz, J.; Oncken, O.; Tassara, A.; Chen, J.; Bataille, K.; Bevis, M.; Socquet, A.; Bolte, J.; Vigny, C.; Brooks, B.; Ryder, I.; Grund, V.; Smalley, B.; Carrizo, D.; Bartsch, M.; Hase, H. 2012. Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake. Earth and Planetary Science Letters 321: 152-165.

Okada, Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America 75 (4): 1135-1154.

Piromallo, C.; Becker, T.W.; Funiciello, F.; Faccenna, C. 2006. Three dimensional instantaneous mantle flow induced by subduction. Geophysical Research Letters 33 (8): L08304. doi: 10.1029/2005GL025390.

Poliakov, A.; Podladchikov, Y. 1992. Diapirism and topography. Geophysical Journal International 109 (3): 553-564.

Prawirodirdjo, L.; Bocl, Y.; McCaffrey, R.; Genrich, J.; Calais, E.; Stevens, C.; Puntodewo, S.; Subarya, C.; Rais, J.; Zwick, P.; Fauzi, R.M. 1997. Geodetic observations of interseismic strain segmentation at the Sumatra subduction zone. Geophysical research letters 24 (21): 2601-2604. doi: 10.1029/97GL52691.

Provost, A.-S.; Chery, J. 2006. Relation between effective friction and fault slip rate across the Northern San Andreas fault system. Geological Society, London, Special Publications 253: 429-436.

Ribe, N.M. 2010. Bending mechanics and mode selection in free subduction: a thin-sheet analysis 180 (2): 559-576. doi: 10.1111/j.1365-246X.2009. 04460.x.

Ruegg, J.C.; Rudloff, A.; Vigny, C.; Madariaga, R.; de Chabalier, J.B.; Campos, J.; Kausel, E.; Barrientos, S.; Dimitrov, D. 2009. Interseismic strain accumulation measured by GPS in the seismic gap between Constitución and Concepción in Chile. Physics of the Earth and Planetary Interiors 175 (1-2): 78-85.

Savage, J.C. 1983. A Dislocation Model of Starins Accumulation and Release at Subduction Zone. Journal of Geophysical Research: Solid Earth 88 (B6): 4984-4996.

Schellart, W.P. 2009. Evolution of the slab bending radius and the bending dissipation in three-dimensional subduction models with a variable slab to upper mantle viscosity ratio. Earth and Planetary Science Letters 288 (1): 309-319.

Scholz, C.H. 2002. The mechanics of earthquakes and faulting. Cambridge University Press: 471 p.

Smith, D.; Turcotte, D. 1993. Contributions of space geodesy to geodynamics: earth dynamics. American Geophysical Union Geodynamics, Series 24: 213 p.

Sobolev, S.V.; Babeyko, A.Y. 2005. What drives orogeny in the Andes? Geology 33 (8): 617-620.

Tassara, A. 2010. Control of forearc density structure on megathrust shear strength along the Chilean subduction zone. Tectonophysics 495 (1-2): 34-47.

Tassara, A.; Echaurren, A. 2012. Anatomy of the Andean subduction zone: Three-dimensional density model upgraded and compared against global-scale models. Geophysical Journal International 189: 161-168.

Thatcher, W.; Rundle, J.B. 1984. A viscoelastic coupling model for the cyclic deformation due to periodically repeated Earthquakes at subduction zones. Journal of Geophysical Research 89: 7631.

Turcotte, D. L.; Schubert, G. 2002. Geodynamics. Cambridge University Press: 456 p. New York.

Van Dinther, Y.; Gerya, T.V.; Dalguer, L.A.; Mai, P.M.; Morra, G.; Giardini, D. 2013. The seismic cycle at subduction thrusts: Insights from seismo-thermo- mechanical models. Journal of Geophysical Research: Solid Earth 118: 6183-6202.

Watts, A.B.; Lamb, S.H.; Fairhead, J.D.; Dewey, J.F. 1995. Lithospheric flexure and bending of the Central Andes. Earth and Planetary Science Letters 134 (1): 9-21.

Zhao, S.; Takemoto, S. 2000. Deformation and stress change associated with plate interaction at subduction zones: a kinematic modelling. Geophysical Journal International 142 (2): 300-318.