Lithospheric Structure of the Central Andes Forearc from Gravity Data Modeling: Implication for Plate Coupling

Geodetic and seismological data indicates that the Central Andes subduction zone is highly coupled. To understand the plate locking mechanism within the Central Andes, we developed 2.5-D gravity models of the lithosphere and assessed the region’s isostatic state. The densities within the gravity models are based on satellite and surface gravity data and constrained by previous tomographic studies. The gravity models indicate a high-density (~2940 kgm) forearc structure in the overriding South American continental lithosphere, which is higher than the average density of the continental crust. This structure produces an anomalous pressure (20-40MPa) on the subducting Nazca plate, contributing to intraplate coupling within the Central Andes. The anomalous lithostatic pressure and buoyancy force may be controlling plate coupling and asperity generation in the Central Andes. The high-density forearc structure could be a batholith or ophiolite emplaced onto the continental crust. The isostatic state of the Central Andes and Nazca plate is assessed based on residual topography (difference between observed and isostatic topography). The West-Central Andes and Nazca ridge have ~0.78 km of residual topography, indicating undercompensation. The crustal thickness beneath the West-Central Andes may not be sufficient to isostatically support the observed topography. This residual topography may be partially supported by small-scale convective cells in the mantle wedge. The residual topography in the Nazca ridge may be attributed to density differences between the subducting Nazca slab and the Nazca ridge. The high density of the subducted Nazca slab has a downward buoyancy force, while the less dense Nazca ridge provides an upward buoyancy force. These two forces may effectively raise the Nazca ridge to its current-day elevation.


Introduction
Large earthquakes commonly occur in Pacific subduction zones where the greatest tectonic strain accumulates [1]. The Cascadia subduction zone, in the northwest of North America, has produced earthquakes greater than M w 8 every 500 years [2]. The Kamchatka subduction zone in the Northwest Pacific had a M w 9 earthquake in 1952 [3], and, southwest of this zone, 10 major earthquakes of M w ≥ 7:5 have occurred in the Japan-Kuril trench in the last century [4]. The area of this study, the Peru-Chile subduction zone (Figure 1), frequently experiences major earthquakes that are among some of the largest recorded (M w > 8:5; [5][6][7][8]). The largest earthquake ever recorded by instruments, M w 9.5, occurred near Valdivia, Chile, in 1960 [1]. In addition to the Valdivia earthquake, the South American subduction zone has had numerous large earthquakes in the past few centuries resulting in multiple full-length ruptures of the entire margin [7]. The subduction zone is well studied in the regions of asperity and high seismic moment release. Areas of low asperity, seismic gaps, are seen throughout the subduction zone [9]. Two noteworthy gaps are present at the Peru-Chile border (18°S, 71°W to 21°S, 71°W) and northern Peru (1°S, 81°W to 11°S, 79°W). These seismic gaps have experienced large earthquakes with magnitudes M w > 8:5 in 1746, 1868, and 1877 [10]. The 1877 Iquique megathrust, M w 8.5-9.0, ruptured ca. 450 km of the thrust fault [11]. Several earthquakes with magnitude M w < 8: 5  far less than the 1877 megathrust, ca. 450 km [12][13][14]. Even with the energy released from these recent earthquakes, it is expected that the zone still possesses enough elastic energy to experience similar events compared to the previous megaquakes in the region [15]. Earthquakes in this subduction zone are primarily caused by the Nazca plate subducting under the South American plate (Figure 1).
Geodetic and seismological measurements in the area indicate regions of high plate coupling, centered around the Central Andes seismic gap zone up to a depth of 50 km, assuming elastic behavior [15]. The correlation between highly coupled zones and megathrust earthquakes has been shown by previous studies [16][17][18]. Additionally, previous studies have shown mass redistribution or fluid pressure variations, seen in vertical gravity gradient changes over time, to potentially indicate areas of asperity generation [19].
What is not clearly understood is the mechanism that controls asperity generation and hence the seismic gaps. As shown by Gutscher et al. [20], geometry of the subducting Nazca plate plays a role in plate coupling. A shallow subduction leads to an increase of contact area resulting in greater compressional strain [21]. Previous research shows a shallow slab in the north that transitions to a steep slab in the south [20][21][22][23].
Two hypotheses have been suggested to explain the locking mechanism of the Peru-Chile subduction zone. The first hypothesis states that a high-density structure, within the forearc, produces a downward buoyancy force, locking the overriding and subducting plates [8,24]. The second hypothesis is that low-density oceanic features on the subducting slab pull the slab upwards during the subduction process and thereby lock the plate interfaces [9,25,26]. In this study, we examine the locking mechanism of the plate interface in the Central Andes based on gravity data modeling. Three 2.5-D gravity models, representative of the Central Andes subduction zone, are presented. The models show the crust and upper mantle structure within the region and help assess the plate locking mechanism in the Central Andes. The density models are based on terrestrial and satellite gravity data and constrained by velocity models from seismic tomography and receiver function.

Geologic Setting
The Andes mountains are a result of the convergence between the Farallon and South American plates, following the breakup of the supercontinent, Gondwana, in the Jurassic [27][28][29][30]. Due to the convergence of the two plates, four    Figure 2). Each of these arcs has a separation from the previous up to 100 km [31][32][33][34][35]. In the Late Oligocene, the Farallon plate split into the Cocos plate, in the north, and the Nazca plate, in the south, due to a change in the convergence angle and an increase in the speed of the southern part of the Farallon plate [8].
After the splitting of the Farallon plate into the Cocos and Nazca plates, the convergence rate of the Nazca plate with South America peaked in the Early-Mid Miocene, 15 cm/yr, [8,23] and has been on the downward trend since, with the current rate at 6.5 cm/yr [28,29].
It is estimated that nearly 200 km of the South American plate has been eroded around the Peru-Chile subduction zone [36]. This erosion is suspected because, around the Antofagasta area, the ancient Jurassic volcanic arc is currently 75 km east of the trench and 35 km above the subducting Nazca plate [5,37]. In addition to the location of the Jurassic volcanic arc, previous studies have shown that the area is lacking an accretionary prism and normal faults are observed from the Coastal Cordillera up to the trench [36,38].
The current formation of the Andes mountains encompasses an area of 800,000 km 2 and has elevations above 6 km, in the central region, and an average of 3.8 km elsewhere [39]. The subduction zone is segmented into four different latitudinal pieces: northern (10°N-5°S), central (5°S-33.5°S), southern (33.5°-46.5°S), and austral (46.5°-56°S) Andes, each varying in age, topography, and volcanism [39]. The oceanic crust is the youngest in the central region, less than 1 Myr at 20°S, increasing in age to the north, 28 Myr at 5°S, and the south, 48 Myr at 46.5°S [8,39]. The area of this study, the Central Andes, can be classified in a longitudinal sense as well ( Figure 2). The westernmost structure is the trench that runs approximately 5900 km down the coast and on average 64 km off the coast of South America. The deepest part of the trench reaches down to 8 km [25]. East of the trench is the Coastal Cordillera, an area of hills, less than 1.2 km in altitude, and valleys, composed of Mesozoic rock [40]. Eastward of the Coastal Cordillera is the Coastal Depression, an area composed of Lower Cretaceous volcanic sediments and plutons [23,41]. The Forearc Precordillera is east of the Coastal Depression and is composed of Paleozoic basement, Eocene magmatic rocks, and Mesozoic and Tertiary sedimentary and volcanic rocks [23,35]. The greatest elevation, 6 km, in the Andes is seen in the Western Cordillera, an area east of the Forearc Precordillera [39]. Flanking the Western Cordillera is the Altiplano-Puna plateau, this plateau is composed of Mesozoic-Cenozoic sedimentary infill that can be 10 km thick [41]. The Eastern Cordillera is the easternmost part of the Andes and is composed of Precambrian to Paleozoic sedimentary rocks within a system of thrust belts [40,41].
Within the region, the Nazca ridge subducts beneath the South American plate. This oceanic feature is less dense and thicker than the surrounded oceanic crust and experiences an upward buoyancy force [20,42,43]. This upward buoyancy force is strong enough to decrease the subduction angle of the plate [20,42,43]. In addition to the buoyancy force, the motion of the overriding plate, South America, has the potential to cause regional changes in the subduction angle [20,42,43]. The variation in the subduction angle can result in plate coupling within the region.

Gravity Database and Methodology
To construct the 2.5-D gravity models of the Central Andes subduction zone, we use the free-air gravity anomaly data from the International Centre for Global Earth Models (ICGEM). These data are combined surface and satellite gravity field data from the EIGEN-6C4 (European Improved Gravity model of the Earth by New techniques) geopotential model [44]. The EIGEN-6C4 model is a combination of land gravity data, altimetry over the oceans, EGM2008 (Earth Gravitational Model; [45]), DTU 2′ × 2′ global gravity anomaly grid (Danish Technical University; [46]), and gravity data from three satellite missions (GOCE (Gravity Field and Steady-State Ocean Circulation Explorer), GRACE (Gravity Recovery and Climate Experiment), and LAGEOS (Laser Geometric Environmental Observation Survey)). These data have a spatial resolution of ca. 10 km and are based on the WGS84 (World Geodetic System) reference ellipsoid, thus allowing the development of a detailed representation of the subduction zone.
We applied standard gravity corrections to the free-air gravity anomaly, and the reductions include Bouguer, terrain, and Bullard-B corrections. The Bullard-B correction accounts for the curvature of the Earth. The Bouguer correction corrects for mass below the measurement point relative to the WGS84 reference ellipsoid. The terrain correction accounts for topographical effects around the measurement point within a 168 km radius. For the terrain correction, we used elevation data from the Global Relief Model of the Earth's surface (ETOPO-1; [47]). The standard reduction density for this study is 2670 kg m -3 . The 2.5-D gravity model was developed using GM-SYS Gravity and Magnetic Modelling software (Geosoft Oasis [48]) which makes use of Green's theorem to calculate the gravity anomalies of irregular structures.

Initial Model and Constraints
In this discussion, we present three gravity models, representative of the Central Andes subduction zone. The density models show the crust and upper mantle structure of the convergent zone with varying tectonic settings, as seen in Figures 1 and 2. The models run along latitudinal lines, span an approximate distance of 9 km, and reach a maximum depth of 250 km. Two southern models, at 20°S and 19°S, represent the southern Peru and northern Chile subduction zones. The northern model, 16°S, is located in the flat slab region where the Nazca ridge subducts beneath the South American plate.
Within these models, major structures are constrained by previous studies using varying seismic velocity models and data. The Moho depth and Lithosphere-Asthenosphere Boundary (LAB), ca. 40 km and 100 km, respectively, have been defined by velocity models developed for the region [9,20,[49][50][51][52][53][54][55]. Table 1 shows the P-wave velocities and densities of major structures in the Central Andes.
The densities of major tectonic structures above the Moho are derived from P-wave velocities at relevant pressure and temperature conditions using the Sobolev and Babeyko    4 Lithosphere [57] method. Sub-Moho structure density is calculated with the Nafe and Drake [58] method. The temperature used for density calculation is based on a subduction zone thermal model from Mahatsente and Ranalli [59]. Additionally, a pressure gradient of 30 MPa/km is used for density determination. These two assumptions, temperature and pressure, can introduce uncertainties into the density determination. We see a negligible amount of uncertainty in Table 2, thus confirming that our assumptions are within the acceptable range. The geometry and density of the subducting Nazca slab and the overriding South American plate are well defined by previous tomography and receiver function studies ( [20,52,54,55]., [50,53]). The oceanic and continental mantle lithospheres have velocity ranges of 7.8-8.15 km s -1 and 7.6-8.1 km s -1 , respectively [52][53][54]56]. The density of the asthenosphere is based on the P-wave velocity model [56].

Results and Discussion
5.1. Gravity Anomaly Analysis. The Bouguer anomaly, shown in Figure 3, indicates significant differences between the ocean (ca. 440 mGal) and the continent (ca. -350 mGal). In addition to gravity differences, there is also an inverse relationship between topography and anomaly. A positive elevation correlates with a negative Bouguer anomaly, whereas a negative elevation correlates with a positive anomaly. This correlation is expected because the Bouguer anomaly is directly related to the density distribution in the crust and mantle. The Bouguer anomaly over the Western and Eastern Cordillera, where the maximum elevation and average crustal root depth are 6.3 km and 50 km, respectively, is on the order of -475 mGal (1 mGal = 10 −5 ms −2 ). The Bouguer anomaly decreases away from the cordillera, a result of the absence of a crustal root. Towards the coast, the elevation approaches sea level and the anomaly reaches a value of 0 mGal. Within the ocean, the elevation decreases to a maximum of 7.6 km depth in the trench. The dense Nazca plate and shallower mantle result in a maximum Bouguer anomaly of 541 mGal. Within this region, the Nazca ridge, ca. 250 mGal, is a prominent feature seen in the Bouguer anomaly. This is a result of the Nazca ridge being less dense than the surrounding oceanic crust. In addition to the change in density, the Nazca ridge also reaches a maximum thickness of 18 km, far larger than the average 6 km thickness of the surrounding oceanic crust [60]. The change in density and  Figure 4: (a) Regional Bouguer anomaly shows the large-scale structures. The crustal root is seen as the long negative anomaly trending NW-SE. (b) Residual Bouguer anomaly shows shallow structures. These anomalies can be attributed to volcanic centers and cordilleras. 5 Lithosphere thickness is a product of the Nazca ridge being an extinct spreading center and having an age of 31 ± 1 Ma at the trench [61].
To further analyze the signal content of the Bouguer anomaly, we applied upward continuation and wavelength filtering. The upward continuation was applied at heights ranging from 10 to 100 km. We chose a final height of 50 km because it effectively showed regional structures. We then applied a regional filter with a cutoff wavelength of 290 km. This wavelength effectively removed the small structure signatures in the data, thus leaving the regional features. The results of the upward continuation and low-pass filter were compared to select the cutoff wavelength. Each showed similar amplitude, thus indicating a regional structure map.
The regional gravity map (Figure 4(a)) shows a large continuous negative anomaly, -350 mGal, over the cordilleras. This negative anomaly represents the low-density crustal root of the cordilleras. The other key structure that can be seen in the regional gravity map is the Nazca ridge. The ridge has an average anomaly of 250 mGal and is surrounded by an oceanic crust that is on average 300 mGal. The gravity signal of the shallow structures in the region, which is the result of the high-pass (short-wavelength) filter (Figure 4(b)), shows positive anomalies in a parallel trend with the cordilleras. Some of these positive anomalies correlate with volcanic centers in the region.

Lithospheric Structure and Locking Mechanism.
Geodetic data indicates a highly coupled plate interface between the subducting Nazca and overriding South American plates beneath the Central Andes subduction zone (Chlieh et al., 2013). We have developed three models in the Central Andes to explore the effects of the forearc structures on plate coupling. We use these models to evaluate our hypothesis of high-density structures causing a downward buoyancy force that effectively locks the Nazca and South American plates.
The first two models shown in Figures 5(a) and 5(b) are located in the southern seismic gap. In this region of north Chile, the slab is subducting at a relatively steep angle, ca. 30°. The location of the slab is constrained by available earthquake data within 0.1 degrees north and south of the profile.
Both models show a high-density structure (2940 kg m -3 ) in the forearc, which is higher than the average density of the continental crust. The high-density forearc structure could be batholith or ophiolite. Within the study area, there is a collection of Middle to Late Cretaceous granodiorite and tonalite rock units, known as the Andean batholith [62]. Previous studies have found indication of batholitic and   [64,65]. This high-density structure could provide enough of a downward force to effectively lock the plate interfaces.
To determine the lithostatic load variations in the forearc, we calculated vertical stress anomalies on top of the subducting Nazca slab along an east-west transect based on the density model ( Figure 6). The vertical stress anomalies are determined relative to a reference lithospheric column, representing the average continental crust and upper mantle. The reference lithospheric column consists of the 15 km thick upper crust, 20 km thick lower crust, and lithospheric mantle. The densities of the reference upper crust, lower crust, and lithospheric mantle are 2670, 2900, and 3350 kg m -3 , respectively. The vertical stress anomaly, resulting from the highdensity forearc structure, on top of the subducting Nazca plate (Figures 5(a) and 5(b)), ranges from 20 to 40 MPa ( Figure 6). This is significantly higher than the pressure exerted on the average continental lithosphere at the same depths. Thus, the high-density forearc structure may be exerting extra pressure on the subducting Nazca slab and thereby locking the plate interfaces.
Additionally, the mantle wedge, resulting from subduction, has the potential to experience small-scale convection [66]. The small-scale convection in the mantle wedge can provide strong enough forces to help inhibit the motion of the subducting Nazca plate. The forces provided by these two features, small-scale convection in the mantle wedge and high-density structures, could effectively lock the subducting Nazca plate.
As we move north, out of the seismic gap, our next model (Figure 7) is located where the Nazca ridge subducts beneath South America. This location has similar densities for the crust and mantle as the two southern models. One substantial difference within the crust is the absence of a high-density structure. Contrary to the south, previous studies indicate a shallow, less than 25 km, slow-velocity (low-density) structure [53]. Even with the absence of a high-density structure, GPS data still indicates that the plates in this region are coupled (Figure 8, [15]).
The plate coupling in this region can be attributed to the subducting Nazca ridge. The low-density material of the ridge causes the subducting Nazca slab to experience an upward buoyancy force, counteracting the downward thermal buoyancy force and pulling the slab upward into the South American plate (Figure 7). The vertical stress anomaly on the subducting Nazca plate that is caused by the lowdensity ridge can be seen in Figure 9. A vertical stress anomaly of ca. 15 MPa is the least amount of stress among the three models, agreeing with the lesser coupling coefficient value seen for this region in Figure 8.
The coupling in both seismic gaps can be attributed to density anomalies in the crustal structure of the forearc. In the south, a high-density structure in the overlying forearc produces a downwards buoyancy force that effectively locks the plates together. Opposite of this, a low-density structure in the north, Nazca ridge, in the oceanic crust produces an upward buoyancy force. Thus, the trench parallel segmentation of the overriding South American and Nazca plates (density and crustal thickness variations) may control plate coupling and asperity generation in the Central Andes.

Model Assumptions and Variations.
For the models presented, various assumptions were made when constructing the initial model. First, we assumed that the subducting Nazca slab is a continuous structure with no slab tearing. This assumption, which is based on available earthquake data and seismic tomography models, ensures a smooth transition between the surface oceanic Nazca plate and the subducting Nazca slab, thus simplifying the model. However, model assumptions and simplifications can change with availability of new additional data. If slab tear is present within the region, then the lithospheric structure could drastically change, in terms of density and shape, when compared to what is presented here. This would be due to the differing geometry of the underlying slab.
Another assumption we made is that the geometry of the overriding continental crust is concave and not convex, as if it were being pulled down by the subducting Nazca plate. This assumption, which is based on seismic data, ensures simplistic structure in the overriding continental crust. If the overriding crustal geometry was convex, then the overriding plate structure would have different geometry and potentially result in differing densities.

Isostatic State of Central Andes.
The continental margin is characterized by recent active volcanism related to volcanic arcs (e.g., Western Cordillera). There is growing evidence of low-velocity zones below the western margin of the Peru-Chile subduction zone caused by deep crustal magmatic reservoirs. This is attributed to partial melting occurring at a depth of 20-60 km [67][68][69]. The effects of low-velocity zones associated with magma chambers and small-scale convection cells in the mantle wedge on the overall isostatic compensation in the forearc of the Central Andes are not well understood. To assess the isostatic state of the Central Andes, we determined the residual topography (difference between observed and isostatic topography) using elevation and crustal thickness models. The elevation dataset is the ETOPO-1 Global Relief Model [47], and the crustal thickness model is the LITHO1.0 model [70]. The computation of isostatic topography is based on the Vening Meinesz isostatic model and includes crustal and mantle support ( Figure 10, [71]).
The Nazca plate and the western part of the Central Andes are characterized by positive residual topography (ca. 0.78 km, Figure 11), indicating that the two regions may not be isostatically supported. Our analysis indicates that the crustal thickness beneath the western part of the Central Andes may not be sufficient to isostatically support the observed topography. The western part of the Central Andes may be partly supported by small-scale convection cells in the mantle wedge.
The extra topography in the Nazca ridge may be attributed to extra forces present in the subducting oceanic lithosphere. One such force could be the result of drastic density changes in the subducting Nazca slab. The highly dense subducted Nazca slab, ca. 3365 kg m -3 at ca. 210 km depth, has a significant downward pull. Opposite of this, the now subducting Nazca ridge is less dense than the surrounding oceanic crust, due to the younger age, and has an upward buoyancy force. These two forces, acting opposite of each other, effectively work in unison to raise the ridge to current-day elevation. The Nazca ridge is an ancient spreading center that is no longer active. According to basalt ages calculated by Ray et al. [61], the ridge is roughly 31 ± 1 Ma at the trench.
The present estimate of residual topography for the Nazca ridge (ca. 0.78 km) is higher than the previously predicted value for the ocean floor based on a global model (0.3 km; [72]). A 0.428 km difference between our study and the global model can be explained by the datasets used. The global model uses data from the CRUST1.0 model [73] whereas our model utilizes the updated LITHO1.0 model [70]. The updated LITHO1.0 model incorporates the CRUST1.0 model parameters and furthers the accuracy by fitting high-frequency surface wave dispersion maps [70].

Conclusions
The most recent major earthquakes in the Central Andes (Iquique 2014 (M w = 8:2) and Illapel 2015 (M w = 8:3)) did not break the entire seismic gap as previously predicted. The plate interface in the seismic gap zone is still highly coupled. The seismic gaps are thought to be the location of the next megathrust earthquake. To understand what is causing these seismic gaps and assess the locking mechanism of plate interfaces, we developed 2.5-D gravity models of the Central Andes subduction zone.
The gravity models in the seismic gap zone (at 20°S and 19°S) indicate the presence of a high-density (2940 kg m -3 ) structure, which is higher than the average density of the continental crust. This high-density forearc structure may be providing a downward force and locking the plate interface within a major seismic gap zone (18°S, 71°W to 21°S, 71°W). Outside of this major seismic gap zone, the gravity model (at 16°S) is lacking in a high-density structure but incorporates the subducting Nazca ridge. The effect of the lowdensity Nazca ridge is opposite in the sense that the structure reduces the negative thermal buoyancy force of the slab and thereby locks the plate interface. Thus, the trench parallel In addition to exploring the lithospheric structure of the Central Andes, we also evaluated the isostatic state of the region. The Nazca plate and western part of the Central Andes are characterized by positive residual topography (ca. -0.20-0.78 km), indicating that the two regions are isostatically undercompensated. The extra topography on the Nazca ridge may be attributed to density changes in the subducting Nazca plate. The high density of the subducted Nazca plate causes a downward force and works in unison with the low density of the Nazca ridge to raise the ridge to modern-day elevation.
The western part of the Central Andes may be partly supported by dynamic processes in the mantle wedge. The smallscale convective cells in the mantle wedge may be providing partial support to the observed topography. The effects of small-scale convective cells on plate coupling have not been investigated in this study. We leave a more detailed evaluation of the effects of small-scale convective cells on plate coupling for future studies.

Data Availability
The underlying dataset used in my research comes from the International Centre for Global Earth Models (ICGEM): http://icgem.gfz-potsdam.de/home. The model I used is the EIGEN-6C4 model and can be found here: http:// icgem.gfz-potsdam.de/calcgrid?modeltype=longtime~ã mp;modelid=7fd8f e44aa1518cd79ca84300aef4b41ddb2364aef9e82b7cdaab db60a9053f1. Other data (not the main data) that have been used are previous research papers, which are referenced specifically in the paper.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

Acknowledgments
This work has been supported by the Department of Geological Sciences, the University of Alabama, and Geological Sciences Advisory Board (GSAB), the University of Alabama.    10 Lithosphere 12 Lithosphere