Prolonged post‐seismic deformation of the 1960 great Chile earthquake and implications for mantle rheology

Contemporary crustal deformation of the southern Andean margin shows an interesting feature: While nearly all coastal GPS sites move landward, consistent with inter‐seismic deformation near a locked subduction fault, sites 300–400 km landward of the rupture region of the Mw9.5 1960 Chile earthquake are moving in the opposite direction. We attribute the seaward motion of these inland sites to a prolonged crustal deformation due to mantle stress relaxation following the 1960 great earthquake. In order to reproduce the observed seaward motion using a three‐dimensional finite element model we need to incorporate a mantle viscosity of about 3 × 1019 Pa s. The possibility that the seaward motion is caused by a silent slip event on the plate interface at large depths cannot be completely excluded, and our analysis provides a working model for future field tests.


Introduction
[2] Large plate boundary earthquakes are known to induce time-dependent post-seismic crustal deformation. For timescales of months to a few years, the deformation is commonly attributed to the continuing aseismic slip of the co-seismically ruptured fault segment [e.g., Heki et al., 1997] and/or its downdip extension [e.g., Barrientos and Ward, 1990]. The afterslip model is simple to use: usually the Earth is represented as an elastic half-space; but the viscoelastic rheology of the earth cannot always be ignored. For example, postglacial rebound analysis by James et al. [2000] indicates that the mantle viscosity in a forearc/backarc environment can be of the order of 10 19 Pa s. Regardless of the presence or absence of afterslip, a viscosity of this magnitude will give rise to mantle stress relaxation that may cause post-seismic deformation, particularly for time-scales of decades or longer.
[3] In the proximity of the rupture zone and very shortly after the earthquake, it is usually difficult to distinguish between the contributions to post-seismic deformation from fault afterslip and mantle relaxation [Pollitz et al., 1998]. However, very large earthquakes, such as great subduction earthquakes of M > 9, may induce large stress in the mantle over a very broad region, and the relaxation of this stress may cause prolonged crustal deformation far away from the coseismic rupture zone. There have been only two earthquakes of this type recorded in the past century, and in both cases, seaward crustal motion has been observed landward of the rupture region several decades after the occurrence of the events. Freymueller et al. [2000] reported similar motions for sites on the western Kenai Peninsula on Alaska, located 300 to 400 km from the trench. Freymueller et al. [2000] and Zweck et al. [2002] attributed these unusual velocities to the continued post-seismic response to the 1964 M w 9.2 Alaska earthquake and the absence of a shallow locked interface. Another example of trenchward oriented crustal motion was reported by Klotz et al. [2001] for the stations located landward of the 1960 M w 9.5 Chile earthquake rupture area. In this paper we interpret these observations in terms of mantle stress relaxation, and estimate the mantle viscosity of the Chilean forearc/backarc using a three-dimensional (3-D) viscoelastic finite element model.
[4] The May 22, 1960 M w 9.5 Chilean earthquake ruptured a 900 km long segment of the plate boundary fault (Figure 1, inset). Using co-seismic uplift observations Barrientos and Ward [1990] estimated the distribution of spatially variable co-seismic fault slip, ranging from 0 to 40 m. If the slip was assumed to be uniform along a rectangular fault of 850 Â 130 km, they obtained an average slip of 17 meters.

GPS Observations
[5] The GPS data used in this study have been obtained under the ongoing South American Geodynamic Activities (SAGA) project initiated in 1993 by the GeoForschungsZentrum (GFZ) Potsdam. Data from 29 stations between 35°S and 43°S latitudes were obtained during two GPS campaigns in 1994 and 1996. Every site was occupied for at least three consecutive days with daily observations exceeding 20 hours. All campaign data, together with data from selected IGS sites, were processed with GFZ software EPOS, using the combined IGS satellite orbits and the Earth orientation parameters. No other constraints were imposed. The fiducial-free solution was transformed to the ITRF97 reference frame using a Helmert transformation constrained by the IGS stations. Detailed descriptions of our data processing procedure and a table of GPS velocities can be found in Klotz et al. [2001].
[6] The velocity vectors presented in Figure 1 are given in a reference frame fixed to a nominal ''stable'' South America (SA) that is defined using IGS stations KOUR, FORT, and LPGS ( Figure 1, inset). The computed rms residuals for the horizontal motion between these three stations do not exceed 2 mm/yr, supporting our assumption that they represent a relatively stable and rigid SA plate. Errors in the campaign GPS horizontal velocities do not exceed 4 mm/yr. The only available velocity estimate in our study area not from the SAGA project is for site ANTC [Kendrick et al., 1999] located several km away from our marker ANTU. This independently derived velocity vector falls within the 95% confidence interval of our velocity estimate for ANTU ( Figure 1).
[7] As reported in [Klotz et al., 2001], most of our GPS observations along the south-central Andean fore-arc are characteristic of inter-seismic forearc deformation due to a locked subduction fault. The GPS sites are moving landward in the direction of plate convergence with their rates decreasing gradually away from the trench. However, the velocity field near the area affected by the 1960 Chile earthquake exhibits an unusual behavior. Inland sites 300 -400 km from the trench are moving seaward in a coherent fashion, opposite to the type of motion expected for the inter-seismic phase of the earthquake deformation cycle ( Figure 1). This motion can not be an artifact due to the choice of reference frame. GPS data from all sites along the margin were processed using the same procedure and referenced to the same nominal stable SA as explained above. The motion is not caused by recent large local earthquakes, as no events of M > 7 have been registered in this area since 1975 ( Figure 1). We are led to the hypothesis that this seaward motion represents post-seismic crustal deformation in response to the 1960 great earthquake by the following facts: 1) the area of the seaward motion spatially correlates with the region affected during the 1960 great earthquake 2) a contiguous set of GPS sites move in a coherent fashion, and 3) the sense of motion is exactly opposite of the direction of plate subduction, hence in the direction of the co-seismic slip during the 1960 earthquake.

3-D FEM Visco-Elastic Model
[8] We test the above hypothesis of prolonged post-seismic deformation using a 3-D viscoelastic finite element model. The model consists of an elastic continental plate of 40 km thickness, an elastic oceanic plate (and slab) of 30 km thickness, and a viscoelastic mantle (Figure 2a). The Young's moduli of the elastic plates and mantle are assumed to be 120 GPa and 160 GPa, respectively, and the Poisson's ratio is assumed to be 0.25 for the entire system. The finite element mesh extends from 1500 km seaward to 1500 km landward of the trench and is 3000 km long in the along-strike dimension ( Figure 2b). The bottom of the model mesh is set at 500 km, in the middle of the mantle transition zone. There are a total of 31,520 tri-linear eight-node elements, with grid spacing ranging from 0.6 km around the subduction fault to over 150 km at the model boundaries (Figures 2b and 2c).
[9] According to our model, the present-day surface deformation is the combined effect of the current interseismic strain accumulation due to the locking of the subduction fault and the delayed response to the 1960 great earthquake due to mantle stress relaxation. The co-seismic rupture of 1960 is modeled by imposing a uniform slip of 20 m in the plate convergence direction [Angermann et al., 1999] along a fault of 900 km long and 120 km wide, plus Figure 1. GPS velocities (with 1s error ellipses) relative to the nominal stable South America defined by the IGS stations shown as diamonds in the inset. Velocity at station ANTC (slightly displaced to the North to avoid overlapping with ANTU) was reported by Kendrick et al. [1999]. All other velocities were reported by Klotz et al. [2001]. Thick dashed line outlines the area of geodetically measured coseismic crustal uplift and subsidence associated with the 1960 Chile earthquake [Plafker, 1972]. Plate convergence vector is according to Angermann et al. [1999]. The epicenter of the main event determined by Cifuentes [1989] is shown as a star. Thin black contour lines represent the depth of the subducting slab. Filled circles show the only two earthquakes with M > 7 that occurred in the past 30 years in this region. an 80 km wide downdip transition from full rupture to zero slip. Similar transition zones are added to the northern and southern edges of the fault to avoid abrupt termination of the co-seismic rupture in these places. The rupture model is consistent with the rupture size and average slip estimated by Barrientos and Ward [1990]. Details of the co-seismic slip distribution do not affect the far-field, long-term postseismic deformation that we consider. The full-rupture area is locked after the earthquake, with locking simulated using the backslip approach [Savage, 1983]. The northern and southern extensions of the subduction fault that did not rupture during the 1960 earthquake remain in a locked state throughout the co-seismic and post-seismic phases of the 1960 earthquake.
[10] For the purpose of studying decadal scale postseismic deformation, most of the short-term fault slip is lumped into our kinematic co-seismic rupture model. To allow a smooth connection between the fault slip and mantle relaxation, we use a thin (1 km) viscoelastic layer along the plate interface as shown in Figure 2a. This layer approximately accounts for the velocity strengthening behavior of the fault downdip from the seismogenic zone. Ideally, one could model the aseismic fault slip using a rate and state dependent friction law, but a 3-D finite element model incorporating both the friction law and viscoelastic rheology is yet to be developed. Further details of the modeling method can be found in Wang et al. [2001].
[11] The model depicts the following deformation process. At the time of the earthquake, the forearc region is instantaneously stretched in the roughly margin-normal direction, consistent with geodetic observations of Plafker [1972]. Note that this ''stretch'' is relative to the initial stress/strain field and represents a reduction of compressive margin-normal stress in the forearc. The instantaneous response of the entire system is elastic. The co-seismic stretch induces shear stress in the mantle, which resists seaward motion of the overlying plate. Therefore deformation occurs only in the proximity of the rupture zone. As the mantle stress relaxes over time, the resistance to motion decreases, and the inland region begins to move seaward. Locking of the fault after the earthquake causes the coastal sites to move landward, while the inland sites are still in the process of catching up with the earlier co-seismic motion. Figure 1 shows the surface velocities 35 years after the earthquake as predicted by the model using a uniform viscosity of 3 Â 10 19 Pa s. Given the simplicity of the model and uncertainties in the GPS data, the model reproduces the overall pattern of the observed velocities quite well. The model would better match the velocities of the coastal GPS sites if a wider locked zone is used north of the 1960 earthquake zone [Klotz et al., 2001], but such details are neglected in the present model.
[12] In the following, we test the sensitivity of the model results to a few important parameters by examining model predicted margin-normal surface velocities along a profile running in the same direction through the center of the rupture zone. Figure 3a shows the effects of changing the uniform mantle viscosity value. Positive velocities represent landward motion. The observed seaward motion of the sites 300 -400 km landward of the center of the rupture zone is about 5 -10 mm/yr (Figure 1), and the viscosity of 3 Â 10 19 Pa s provides a reasonable fit to these velocities. Increasing or decreasing the viscosity by 2 Â 10 19 Pa s results in much poorer fits.
[13] Results shown in Figure 3b are designed to test the importance of the time-dependent fault motion downdip from the rupture zone. We have used a thin viscoelastic layer to approximately account for the possible post-seismic fault slip. Varying the viscosity of this layer by 2 Â 10 19 Pa s while keeping the mantle viscosity at 3 Â 10 19 Pa s only slightly affects the velocities of the inland region 35 years after the earthquake. This supports the argument that details of the time-dependent frictional behavior of the fault do not have a large impact on far-field post-seismic deformation. In contrast, varying the mantle viscosity by 2 Â 10 19 Pa s while keeping the thin-layer viscosity at 3 Â 10 19 Pa s (results not shown) produces results very similar to what is shown in Figure 3a. The decadal scale, large wavelength deformation depends mainly on the mantle viscosity.
[14] Results in Figure 3c show that the prolonged postseismic deformation is a unique behavior of long-rupture earthquakes. If the along-strike length of the rupture is 200 km, there is no seaward motion of the inland sites 35 years after the earthquake. This is because a short rupture produces too small an initial perturbation (the seaward ''stretch'') to the mantle.

Discussions
[15] The mantle viscosity of 3 Â 10 19 Pa s deduced from our FEM modeling corresponds well with the findings of Piersanti [1999]  using a viscoelastic spherical Earth model and suggested asthenospheric viscosities between 3 -4 Â 10 19 to 1 -2 Â 10 20 Pa s. Studies of the fluctuation of the southern Patagonian ice-fields during the last 5000 years [Ivins and James, 1999] suggest a mantle viscosity for the southern tip of South America of 5 Â 10 18 to 5 Â 10 19 Pa s, consistent with our findings. Our attempt to reproduce the above mentioned uplift rates using the FEM model have not been successful; the simple treatment of the fault does not allow any afterslip that may be required to explain observed costal vertical motion after the earthquake, but the crustal motion 400 km inland of the trench is less sensitive to details of the fault behavior.
[16] Another possible explanation for the seaward motion of inland sites is a transient silent slip on the plate interface. Dragert et al. [2001] reported a silent slip event on the plate interface of the Cascadia subduction zone downdip from its currently locked seismogenic portion. A cluster of continuously monitoring GPS stations about 100 to 250 km from the trench were observed to briefly reverse their direction of motion. Similar transient motions have been reported for a few other subduction zones [Lowry et al., 2001;Ozawa et al., 2001]. Our area of seaward motion is farther away from the trench (200 to 400 km) than that reported by  for Cascadia. If the motion was due to a silent slip on the plate interface, the slip had to take place in the 50 to 150 km depth range. It had to be quite a coincidence that a transient motion occurred at the time of our 1994 and 1996 campaign surveys and only in the region of the 1960 earthquake. Nevertheless, we cannot completely exclude the possibility of a silent slip. Our model predicts that the seaward motion of those inland sites should continue for the next few decades with decreasing rates. The motion direction of these sites in the next few years will support or invalidate our hypothesis.

Conclusions
[17] GPS campaign surveys at the Chile subduction margin indicate landward motion of coastal sites but seaward motion of inland sites. Over the entire length of the SAGA GPS network of more than 2500 km along the Andean margin, this phenomenon was observed only in the region landward of the great 1960 Chile earthquake rupture area. We explain the seaward motion in terms of post-seismic crustal deformation due to mantle stress relaxation in response to this earthquake. If the hypothesis of post-seismic deformation is valid, our 3-D viscoelastic finite element modeling leads to the following conclusions: 1) The average mantle viscosity for this active continental margin is about 3 Â 10 19 Pa s. 2) Details of the short-time fault motion immediately after the earthquake do not significantly affect the decadal scale deformation due to relaxation. 3) Along-strike rupture length has a first order control on the long-term post-seismic behavior. Prolonged post-seismic deformation is a unique feature of long-rupture great earthquakes.