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

. We study interseismic deformation preceding the Mw8.8 2010 Maule earthquake by means of two-dimensional finite-element modeling. Our goal is to gain insight into the fundamental factors controlling elastic strain build-up and release in subduction zones, and to evaluate different modeling approaches of surface displacement as observed by GPS. We developed a linear elasticity solver that allows us to implement a realistic subducting plate geometry constrained by geophysical data. We test the influence of subducting plate thickness, variations in the updip and downdip limit of a 100% locked interplate zone, elastic parameters, and velocity reduction at the base of the subducted slab. We compared our modeled predictions with interseismic GPS observations along an EW profile crossing the Maule earthquake rupture area, in order to determine best fitting parameters. Our results indicate little influence of the subducting plate thickness at a given downdip limit, which itself has a strong influence on surface deformation. However, the fit to observations is achieved only after reducing the velocity at the base of the subducted slab below the trench region to 10% of the far-field convergence rate. We link this novel result to complementary numerical models that gradually evolve toward considering longer time-scales and complex rheology in order to evaluate the mechanical meaning of the above mentioned inferred kinematic conditions. This allowed us to link the velocity reduction at the base of subducting slabs with a long-term state of high flexural stress resulting from the mechanical interaction of the slab with the underlying mantle. Even a small amount of theses high deviatoric stresses may transfer towards the upper portion of the slab as strain energy that could participate into the mechanical loading of the megathrust and therefore in triggering large earthquakes there.


Introduction
The megathrust interplate shear zone along subduction zones is seen to behave frictionally locked during the time lapse between two large earthquakes. Interplate locking during this interseismic period can be modeled by an elastic deformation of the upper plate that translates into the displacement of the surface in the direction of convergence, a mechanism that has been widely used to explain observations of conventional geodetic measurements (see for example Thatcher and Rundle (1984) and Hashimoto and Jackson (1993) who modeled century-long geodetic data in Japan). Modern satellite geodesy, mostly the Global Positioning System (GPS), has improved the measurement of surface displacements at centimeterlevel accuracy, allowing a better understanding of the mechanisms of interplate coupling and the generation of great destructive earthquakes (e.g., Smith and Turcotte, 1993;Scholz, 2002;Avouac, 2015). In order to connect geodetic observables at the Earth's surface and plate coupling along the megathrust fault at depth, different approximations can be used. In this contribution we use 2D numerical models of subduction zones to analyze the influence of different parameters and several plate boundary approximations (including the Back Slip Model, BSM and, Elastic Subducting Plate Model, ESPM) on kinematic and mechanical conditions accompanying the seismic cycle. To do this we constrain the deformation caused by different models using interseismic surface velocity vectors measured by GPS previous to the Mw8.8 Maule 2010 earthquake along a trench-normal transect located near the Arauco Peninsula in south-central Chile (Fig. 1).
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 RESUMEN. Deformación intersísmica en zonas de subducción investigada por modelos numéricos 2D: estudio de caso antes del terremoto del Maule 2010. Estudiamos la deformación intersísmica que precedió al terremoto Mw8.8 del Maule 2010 por medio de modelos de elementos finitos en dos dimensiones. Nuestro objetivo es comprender los factores fundamentales que controlan la acumulación y liberación de strain elástico en zonas de subducción, y evaluar distintas aproximaciones en la modelación del desplazamiento superficial observado por GPS. Desarrollamos un código de elasticidad lineal que permite implementar geometrías realistas de la placa subductada según datos de control geofísicos. Mediante esta herramienta testeamos la influencia del espesor de la placa subductada, variaciones en el límite superior e inferior de una zona 100% acoplada del contacto interplacas, parámetros elásticos del material, y reducción de velocidad en la base de la placa subductada. La comparación de las predicciones del modelo contra observaciones de GPS intersísmicas en un perfil EW que cruza el área de ruptura del terremoto del Maule permite determinar los parámetros que mejor reproducen las observaciones. Nuestros resultados indican poca influencia del espesor de la placa subductada para un límite inferior de acople dado, el que a su vez tiene una gran inf luencia sobre la deformación superficial. Sin embargo, el ajuste a las observaciones se logra solo luego de reducir la velocidad en la base del slab bajo el área de la fosa a un 10% de la velocidad de convergencia. Nosotros ligamos este novedoso resultado con modelos numéricos complementarios que gradualmente incorporan escalas de tiempo largas y complejidades reológicas con el fin de evaluar el significado mecánico de este resultado cinemático. Esto nos permite relacionar la reducción de velocidad en la base de la placa subductada con un estado de alto estres f lexural en el largo plazo resultante de la interacción mecánica de la placa con el manto circundante. Incluso una pequeña proporción de estos altos estreses deviatóricos podrían ser transferidos a la parte superior de la placa subductada como energía elástica que pudiera a su vez participar en la carga mecánica de la falla de subducción y por tanto en la generación de grandes terremotos. 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 landbased 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, andPAROVOZ: 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.

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.

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.
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.
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.

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.

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.5x10 11 Pa, respectively.
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.   ) 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.
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).

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 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.
section, e.g., thought to be appropriate for slab-mantle interaction at the short-scale, we progressively evolve towards boundary conditions appropriate for longterm slab-mantle interaction. We thus discuss step by step the implications of some assumptions, and how our reasoning brings insight on interseismic seismogenic behavior.

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).
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 f low, 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., 2009or 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 slabmantle 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 f lexural 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 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.
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), f lexural 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.
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.

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 selfconsistently 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 lowviscosity 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 andMorra, 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.

Conclusions
The well known Back Slip Model (BSM; Savage, 1983) solution of an elastic dislocation in a semiinfinite 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 f lexure, 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.