Earthquake-cycle deformation, which includes earthquake ruptures, interseismic strain, and transient slow slip events, spans spatial scales ranging from fractions of a meter to thousands of kilometers. Similarly, temporal scales range from seconds during an earthquake rupture to thousands of years of strain accumulation between earthquakes. We discuss results regarding the vertical crustal deformation associated with both slow and rapid crustal deformation across a transect of the central Mexican subduction forearc in the Guerrero seismic gap, where the Cocos plate underthrusts the North America plate. This sector of the subduction zone is characterized by a flat-slab geometry with zones of sharp bending-unbending of the slab, irregularly distributed seismicity, and exceptionally large slow slip events. We used the river network, topography, geomorphic features, and morphometry on a transect across the forearc to assess Quaternary crustal deformation. The Papagayo drainage network shows that the forearc has been uplifted since the late Cenozoic (∼25 Ma), and that rates of uplift increased since the beginning of the Holocene. Uplift is not homogeneous but shows a trend of increase away from the coast. This vertical deformation is strongly influenced by subduction processes. Thus, the Papagayo River network is strongly controlled by Holocene earthquake cycle processes. This is particularly true for the southern section of the drainage basin, where E-W–striking left-lateral strike-slip faults with a vertical component offset the course of the main river. These faults are accommodating part of the oblique plate convergence at the Mexican subduction zone. We measured the height of a series of terraces and dated quartz extracts by optically stimulated luminescence, and we calculated long-term rates of uplift ranging from 0.5 to 4.9 mm/yr. We discuss associations of forearc topography, faults, and long-term crustal deformation with the Cocos slab geometry, distribution of slow slip events, and earthquake-cycle deformation.
Spatial scales of earthquake cycles involving rapid earthquake rupture, interseismic deformation, and transient slow slip behavior span from fractions of a meter to thousands of kilometers at plate boundaries. Similarly, temporal scales range from seconds during an earthquake rupture to thousands of years of strain accumulation between earthquakes. The complexity of tectonic processes operating over a range of scales and the limited duration of observations often focus inquiry onto a limited space-time context that may not fully characterize or quantify tectonic variability (e.g., Kostoglodov et al., 2003, 2010; Larson et al., 2004; Franco et al., 2005; Marushima and Ishiyama, 2009; Pacheco and Singh, 2010; Walpersdorf et al., 2011). Thus, some of the most important questions concerning subduction zones remain partially addressed, such as: How is tectonic deformation distributed among various modes of fault slip, particularly for the unique conditions of subduction zones? Also, how spatially and temporally variable are megathrust slip patterns, and what processes and properties control variability? How do processes at depth and in the near surface vary? Lastly, what is the long-term expression of earthquake-cycle deformation processes observed in the geology and geomorphology of the subduction forearc?
Vertical deformation of the forearc during coseismic, postseismic, and interseismic stages of the earthquake cycle has been measured in several subduction zones along the Pacific rim, such as the Japan subduction zone (Aoki and Scholz, 2003; Ozawa et al., 2011; Wang et al., 2012), New Zealand (Beavan et al., 2007, 2010), Indonesia (Sieh et al., 1999; Sieh and Natawidjaja, 2000; Meltzner et al., 2006; Chlieh et al., 2007; Prawirodirdjo et al., 2010), Chile (Melnick et al., 2006, 2012b; Vargas et al., 2011; Wesson et al., 2015), and the Mexican subduction zone (Hutton et al., 2001; Azúa et al., 2002; Kostoglodov et al., 2003; Larson et al., 2004; Wei et al., 2011). Sense of vertical motion reverses between the stages of deformation (i.e., coseismic, postseismic, and interseismic). However, some areas are more complex, such as the coast of northeastern Japan, which experienced late-interseismic and coseismic subsidence related to the A.D. 2011 Tohoku earthquake (Japan) although associated landforms indicated long-term (103 to 106 yr) uplift (Sagiya, 2015; Hashima and Sato, 2017). Vertical deformation also varies spatially across and along the forearc. The A.D. 2004 Sumatra-Andaman earthquake (Indonesia) produced regions both of uplift and of subsidence depending on distance from the trench (Meltzner et al., 2006; Prawirodirdjo et al., 2010). Similarly, the 2011 Tohoku earthquake produced vertical coseismic and postseismic displacements across the forearc varying in magnitude and with either subsidence or uplift (Ozawa et al., 2011). These examples show how the vertical motions of the coasts can differ from the classical elastic model for stress accumulation and release in the upper plate, and they raise questions regarding (1) the pattern of net deformation that is accumulated through time over several seismic cycles, and (2) the driving mechanisms. We attempt to bridge the time scales of the earthquake cycle and vertical deformation by examining the geodetic record (101 yr) and the geologic record (≥103 yr) of vertical deformation of the overriding plate, and related earthquake events. Our goal is to understand the fundamental mechanism for uplift of the forearc inland through different time periods and to analyze individual geodetic results.
Here, to address how surficial and deep processes shape the forearc, we analyze topography and fluvial terraces on the overriding plate. An important source of information to determine vertical deformation of the forearc on time scales of 102 to 106 yr are uplifted fluvial terraces. The deformation of the surface of fluvial terraces, and in particular strath terraces, has been used as an indicator of displacement and long-term uplift of the upper plate in several subduction zones including the Cascadia (Oregon and Washington; Personius, 1995; Pazzaglia and Brandon, 2001; Wegmann and Pazzaglia, 2002), Peruvian (Pedoja et al., 2006; Hall et al., 2012), Chilean (Farías et al., 2005, and references therein), New Zealand (Berryman et al., 2010), and Japan (Yamamoto, 2005; Matsu’ura et al., 2008; Marushima and Ishiyama, 2009), and at the Mendocino triple junction (northern California; Merritts et al., 1994). Strath terraces can also be used to determine spatial variations in displacement and long-term uplift, because incision rates within a drainage basin can change both spatially and temporally (Burbank et al., 1996; Pazzaglia and Brandon, 2001; Righter et al., 2010).
However, few studies have focused on the long-term vertical deformation of the overriding plate in the Mexican subduction zone (Ramírez-Herrera and Urrutia-Fucugauchi, 1999; Ramírez-Herrera et al., 2007, 2009, 2011). The Mexican subduction zone is unusual, with shallow earthquakes close to the coast (Pacheco and Singh, 2010), an inland volcanic arc running oblique to the trench at ∼15° (García-Palomo et al., 2004), long-lasting slow slip events (SSEs) (Kostoglodov et al., 2003, 2010; Larson et al., 2004, 2007; Franco et al., 2005; Vergnolle et al., 2010; Walpersdorf et al., 2011; Cavalié et al., 2013), and erosional subduction (Clift and Vannucchi, 2004). Most studies on long-term vertical deformation have focused on accretion-dominated subduction zones, with few exceptions (e.g., Marshall and Anderson, 1995; Matsu’ura et al., 2009). Therefore, we focus on quantifying the vertical deformation rate for the past ∼5 k.y. for the Mexican subduction forearc within the Guerrero seismic gap to understand the long-term geologic uplift in contrast to the short-term, interseismic geodetically observed subsidence. Late Quaternary crustal deformation was assessed through analysis of the river network, topography, and geomorphic features across a transect of the Mexican subduction zone forearc, in the Guerrero gap, and along a profile of the Papagayo River. In turn, the height was measured of a vertical sequence of terraces, and optically stimulated luminescence (OSL) ages were obtained from quartz extracted from fluvial sediments, to estimate long-term rates of vertical deformation. The association amongst topography, faults, fluvial terraces, and river morphology is discussed and related to characteristics of the Cocos slab geometry, earthquake cycle deformation, and SSEs.
CENTRAL MEXICAN SUBDUCTION ZONE AND GUERRERO SEISMIC GAP
Geodynamic Setting and Seismic Activity
The study area is in the forearc of the Cocos–North America subduction zone, within the Guerrero seismic gap, a ∼200 km segment between 99.2° and 101.2° W (Fig. 1). No large thrust earthquakes have occurred in the northwestern part of the gap since the Ms 7.6 earthquake in 1911 (Anderson et al., 1989; Kostoglodov and Ponce, 1994; A.R. Lowry, personal commun., 2016) (Fig. 1). The region southeast of Acapulco has experienced mostly modest (Mw <7.1) earthquakes since 1957. The total seismic moment of these moderate earthquakes is insufficient to accommodate the expected strain accumulation in the Guerrero gap, suggesting substantial contributions of aseismic slip (Suárez et al., 1990; Anderson et al., 1994). Large SSEs, Mw ∼7.5, occurred on the subduction interface in 1998, 2001–2002, 2006, 2009–2010, and 2014 (Kostoglodov et al., 2003 and 2015; Walpersdorf et al., 2011). The width of the coupled, seismogenic plate interface in the Mexican subduction zone is narrow, ∼60 km, to a depth of 25 km based on the hypocenters of shallow thrust earthquakes (Pardo and Suárez, 1995). However, recent studies using continuous GPS observations on land and dislocation modeling suggest that the interplate coupling is broader and reaches 180–220 km inland from the trench in the northwestern section of the Guerrero seismic gap. The long-term coupling in the gap probably remains low (∼25%) on the shallow and narrow seismogenic zone while the downdip plate interface is strongly coupled (∼80–90%) only during short, 3–4 year inter-SSE epochs.
The subduction stress for the oblique Cocos–North America plate initiated ca. 22 Ma with the breakup of the Farallon plate into the Cocos and Nazca plates (Wortel and Cloetingh, 1981; Barckhausen et al., 2001). The oceanic lithosphere of the Cocos plate subducts beneath the North America continental plate with convergence rates from 6.4 to 6.7 cm/yr along the Guerrero coast (Fig. 1) (DeMets et al., 2010). Relative motion is slightly oblique (∼12°) to the trench, yielding a left-lateral component of ∼8 mm/yr (A.R. Lowry, personal commun., 2016). Convergence rates between the Cocos and North America plates along the Middle America Trench and the geometry of the subducting slab in the Guerrero gap appear to control the crustal deformation in this forearc region and the earthquake thrust process. A variety of models have been proposed for the geometry of the slab interface (Kostoglodov et al., 1996, 2003; Manea et al., 2004; Pérez-Campos et al., 2008; Kim et al., 2010) (Fig. 2). The best-fit model used in this study is based on hypocentral locations and gravity anomalies (Kostoglodov et al., 1996, 2003; Manea et al., 2004) (Fig. 2). In this model, the slab interface dips gently (∼5°–10°) down to ∼12 km depth and ∼77 km from the trench, steepens to ∼35° at ∼40 km depth and ∼115 km eastward from the trench, and becomes almost horizontal beneath central Mexico (Fig. 2).
GPS observations reveal uplift of the forearc on the Guerrero sector of the Mexican subduction zone during SSEs, while subsidence predominates between SSEs (Kostoglodov et al., 2003; Vergnolle et al., 2010; Walpersdorf et al., 2011; Radiguet et al., 2012). SSEs occur on a cycle of 3–4 yr and are recorded by the Mexican GPS network. The estimated SSEs duration (6–11 months), magnitude of slip (Mw ∼7.0–7.5), and slip rate (∼0.1–0.5 m yr−1) for the Guerrero sector of the Mexican subduction zone fault slip are amongst the largest reported (e.g., Ide et al., 2007).
The study area is located across the forearc in central Mexico within the Guerrero morphotectonic zone (Ramírez-Herrera and Urrutia-Fucugauchi, 1999). This coast is characterized primarily by deposition, with the landscape dominated by low deposition deltaic and alluvial plains. Coastal lagoons are commonly separated from the ocean by barrier beaches and rocky promontories. Most of the largest river basins drain into the Pacific Ocean, eroding granitoid bodies or metamorphic rocks (mainly gneisses), and upstream forming deep mountain valleys with steep slopes (Fig. S11). The elevations in the studied area range from sea level to ∼3500 m upstream.
We selected the Papagayo River basin to investigate tectonic processes across the forearc, from the coast and inland, because it is the largest catchment (>7500 km2) in this region. This catchment exhibits steep average slopes of ∼20° and a series of river terraces along the main course of the river, knickpoints, and apparent lateral offsets of the main course of the river. This area is also directly inland of the locked section of the slab of interest where the SSEs occur (Figs. 1, 2).
The coastal area in Guerrero has no marine terraces indicative of coastal tectonic uplift in the late Quaternary (Ramírez-Herrera and Urrutia-Fucugauchi, 1999). Recent studies of coastal sedimentary records indicate coastal subsidence during the Holocene of the Guerrero sector at rates of ∼1 mm/yr (Ramírez-Herrera et al., 2007, 2009, 2011). Although there appears to be coastal subsidence, further inland beyond the coastal-alluvial plains there are uplifted river terraces above bedrock channels of the Papagayo River.
A number of E-W–trending crustal faults in the study area are expressed clearly in the landscape and river morphology; here, we focus on the topographic expression of faults and potential correlation with processes at depth. Particularly well expressed in the landscape, downstream of the Papagayo basin, are two E-W structures that we named the Cacahuatepec fault and the Dos Arroyos fault (Gaidzik et al., 2016) (Figs. 3 and S1 [footnote 1]). Deflections of the Papagayo River suggest strike-slip movement of these structures. However, an abrupt change in elevation north and south of these faults (Fig. 3B) suggests a vertical slip component. Slickenlines and kinematic indicators recorded in the field along the trace of these crustal structures indicated a left-lateral strike-slip motion. In some locations are E-W left-lateral strike-slip meso-faults with a normal component. The slip mode for these faults remains unsolved despite the faults’ clear geomorphic expression (Gaidzik et al., 2016). It is unknown whether these faults remained inactive for the 20th and early 21st centuries. The focal mechanisms for this region encompass a relatively short period of time from 1995 to 2007 (Pacheco and Singh, 2010). The Global Centroid-Moment-Tensor (Dziewonski, et al., 1981; Ekström et al., 2012) catalog records focal mechanisms with reverse kinematics dipping gently toward the NNE nodal planes and steeply dipping toward the SSE at depths from 25 to 40 km for all events in the study area.
On the basis of structural field data and morphotectonic interpretation of a digital elevation model (DEM) and satellite images, we also identified a remarkable ∼200 km left-lateral strike-slip crustal fault, including a clear offset of the main channel of the Papagayo River; this we named the La Venta fault (Fig. 3; Fig. S1 [footnote 1]). Meso-faults with kinematic indicators recorded in the field (slickenlines, subordinate R, R′, P, and T fractures [Fossen, 2010], and displacement of veins) close to the dam in La Venta confirmed the left-lateral activity of one of the segments of the La Venta fault that we named the Agua del Perro segment. Even though the strike-slip motion predominates, measurements taken in the field of the spatial orientation of meso-fault surfaces with slickenlines indicated that the vertical component can be up to ∼35% of the total displacement along this fault (Gaidzik et al., 2016). The activity of the La Venta fault and Agua del Perro segment was confirmed by GPS measurements that show a 4–5 mm/yr left-lateral slip of the fault and SSEs that occur along this fault concomitantly with the subduction thrust SSEs (Kostoglodov et al., 2015).
Geodetic Measurements of Short-Term Crustal Deformation
We used published continuous GPS data from the Servicio Sismológico Nacional-Sismología-Universidad Nacional Autónoma de México permanent GPS network located on land to compare the rates of short-term vertical deformation with our long-term (Holocene) rates. GPS data come from time series starting in 1997 and ending in 2014 from the DOAR, CPDP, ACAP, ACYA, COYU, and CAYA permanent stations in the Guerrero zone. We focused on vertical displacements for both observed slow aseismic slip events and interseismic deformation (e.g., Kostoglodov et al., 2003, 2010; Larson et al., 2004, 2007; Franco et al., 2005; Vergnolle et al., 2010; Walpersdorf et al., 2011; Cavalié et al., 2013). According to these data, SSEs occurred on the subduction zone on the northwestern Guerrero seismic gap, and were observed every 3–4 yr in 1998, 2001–2002, 2006, 2009–2010, and 2014 (Kostoglodov et al., 2003, 2010; Franco et al., 2005; Vergnolle et al., 2010; Walpersdorf et al., 2011; Cavalié et al., 2013). The down-dip slip limits for the 2001–2002 and for the 2006 and 2009–2010 SSEs were 200 and 150 km from the coast, respectively. SSEs are reported at 15–20 km depth for an up-dip slip limit in the gap area, and deeper outside the gap (Radiguet et al., 2012). GPS measurements on land and modeling suggest high coupling during the interseismic period in the northwestern Guerrero seismic gap sector of the Mexican subduction zone, with a coupling ratio (i.e., the ratio of the slip deficit rate over the convergence rate) of >0.7 on the plate interface at a distance of 10–90 km from the coast, at approximate depths of 10–40 km (Radiguet et al., 2012). Further inland, the coupling decreases as the slab angle approaches horizontal (Radiguet et al., 2012). During the two largest SSEs (2001–2002 and 2006), uplift of as much as 70 mm near the coast and subsidence of 2–10 mm inland toward La Venta was reported by Vergnolle et al. (2010) and Walpersdorf et al. (2011). In the same area, range-change data from radar interferometry (InSAR) also indicated uplift of 40 mm during the 2006 SSE (Cavalié et al., 2013). In contrast, geodetically measured inter-SSE vertical surface displacements along the Pacific coast of Mexico range from about 20 mm/yr subsidence on the coast, to a few millimeters of uplift >150 km eastward from the coast (Kostoglodov et al., 2003). The horizontal deformation ranges from 10 to ∼25 mm/yr (Kostoglodov et al., 2003).
Published GPS displacement data indicate that secular GPS velocity vectors are oblique to the trench and decrease by 4–5 mm across the La Venta fault (i.e., north and south of the fault) (Gaidzik et al., 2016; this study), ∼100 km inland from the trench in the Guerrero sector. The decrease in velocity components across this fault suggests left-lateral motion of this fault. Kostoglodov et al. (2015) inferred that strike-slip SSEs occur simultaneously along this fault with the subduction thrust SSEs in the Guerrero sector of the Mexican subduction zone.
Geomorphic Indices and Landscape Analysis
Geomorphic indices were used to determine the relative tectonic activity and to differentiate between areas of varying uplift rates (e.g., Bull and McFadden, 1977; Bull, 1978; Rockwell et al., 1985; Keller, 1986; Burbank and Anderson, 2001; Keller and Pinter, 2002). Here, we calculate the following geomorphic indices: (1) stream length-gradient index (SL) (Hack, 1973; Keller and Pinter, 2002), (2) ratio of valley floor width to valley height (Vf) (Bull and McFadden, 1977; Keller and Pinter, 2002), and (3) minimum bulk erosion (Ebulk) (Giaconia et al., 2012; Bellin et al., 2014). We use a 15-m-resolution DEM and digitized topographic maps (scale 1:50,000; Instituto Nacional de Estadística y Geografía, 2014). We also constructed topographic swath profiles using maps of minimum, maximum, and average elevation from the DEM at 1 × 1 km resolution. The valley shape and the channel profile along the north-south swath topographic profile (50 × 15 km) were quantified to study changes in relief and elevation. Additionally, 70 cross-sections were used through the Papagayo River valley at 1 km intervals (within four defined sections of the river from La Venta to the coast) (Fig. 3) to better assess fault offset and deformation.
Stream Length-Gradient Index (SL)
Ratio Of Valley Floor Width To Valley Height (Vf)
Low values of Vf indicate deep V-shaped valleys in highly active areas, whereas high Vf values (>1.0) suggest broad flat-floored valleys, typical for tectonically quiescent areas, or regions with low uplift rates (Bull and McFadden, 1977; Bull, 1978; Rockwell et al., 1985; Keller, 1986; Keller and Pinter, 2002). We calculate for the Papagayo River a Vf ratio at 1 km intervals along the river using cross-sections from the DEM (Fig. 4D). The average of three measurements was used because of limited resolution of the DEM for valley width.
Minimum Bulk Erosion (Ebulk)
Minimum bulk erosion (Ebulk) (Giaconia et al., 2012; Bellin et al., 2014), or “incision” (Burbank and Anderson, 2001), is defined as the difference between a theoretical pre-incision surface and current topography. In that respect, it differs from the topographic residual, which is calculated by subtracting a smooth surface interpolated from the stream bottoms from a pre-incision surface (Hilley and Arrowsmith, 2000). Shown as a map, the minimum bulk erosion expresses the spatial distribution of the amount of material that eroded within the drainage basin (Fig. 5). Therefore, it allows highly eroded areas to be distinguished from weakly eroded surfaces (Giaconia et al., 2012). We calculated this index using ArcGIS software (http://www.esri.com/arcgis/about-arcgis), following the “sub-ridgeline relief” procedure proposed by Brocklehurst and Whipple (2002). For the difference between the DEM and the rasterized pre-incision surface, we used the raster calculator from ArcGIS software.
To identify, classify, and map river terraces, we used aerial photograph interpretation (scale 1:30,000), a DEM obtained from aerial photographs using a structure-from-motion procedure (e.g., Snavely et al., 2008; James and Robson, 2012; Johnson et al., 2014), satellite image interpretation, and field surveys. We also constructed the long profile for the Papagayo River from digitized topographic maps (available at 1:50,000 scale, 10–20 m contour lines [Instituto Nacional de Estadística y Geografía, 2014]) and long profiles for each terrace surface from field surveying data. We determined the elevation of each terrace above the stream channel by field survey using a laser rangefinder (accuracy ±1 m) and distinguished, where possible, the height of the bedrock on strath terraces and the surface elevation for both fill and strath terraces (Fig. 6). The position of each terrace was recorded in the field with a handheld GPS.
Optically Stimulated Luminescence (OSL) Dating
Optical dating is a suitable method for deciphering the late Quaternary record of fluvial deposition in the Papagayo River drainage. The luminescence emission is a measure of the burial period during which, following its exposure during grain transport and deposition, the material has been shielded from any further exposure to light (Aitken, 1998). The development of single-aliquot regeneration (SAR) protocols provides improved accuracy and precision for dating the burial time of quartz grains (e.g., Murray and Wintle, 2003; Wintle and Murray, 2006), particularly in fluvial environments, although there are concerns about full solar resetting in higher-energy depositional settings (cf. Rittenour, 2008; Meier et al., 2013; Forman, 2015). SAR protocols (Wintle and Murray, 2006) were used in this study to estimate the apparent equivalent dose of the 63–100 µm, 100–150 µm, 250–355 µm, or 355–425 µm quartz fraction for up to 60 separate aliquots. Each aliquot contained ∼50 to hundreds of quartz grains on a 1–2 mm-diameter plate area on a 1–164 cm-diameter circular aluminum disc. The quartz fractions were separated by density with the heavy liquid Na-polytungstate, and a 40 min immersion in HF (40%) etched the outer ∼10 µm of grains, which was affected by alpha radiation (Mejdahl and Christiansen, 1994). Quartz grains were rinsed in HCl (10%) to remove insoluble fluorides. The purity of the quartz separate was evaluated by petrographic inspection, point counting of a representative aliquot, and pre-screening of aliquots through infrared excitation.
An automated Risø TL/OSL-DA-15 system was used for SAR analyses. Laboratory irradiations used a calibrated 90Sr/90Y beta source coupled with the Risø reader. Optical stimulation for all samples was completed at an elevated temperature (125 °C) achieved at a rate of 5 °C/s. All SAR emissions were integrated over the first 0.8 s of stimulation out of 40 s of measurement, with the background based on emissions for the last 30–40 s interval. A series of experiments evaluated the effect of preheating at 160, 180, 200, 220, and 240 °C for isolating the most robust time-sensitive emissions and minimizing thermal transfer of the regenerative signal prior to the application of SAR dating protocols (Murray and Wintle, 2003). These experiments entailed giving a known dose of 22 Gy and evaluating which preheat resulted in recovery of this dose. There was concordance with the known dose (22 Gy) for preheat temperatures between 160 and 240 °C with an initial preheat temperature of 180 °C for 10 s, a second preheat heat of 160 °C for the test dose, and a final annealing temperature of 260 °C in the SAR protocols.
An assumption a priori is that quartz grains in this fluvial depositional system were not uniformly reset by sunlight because of a possible short distance of transport, turbid depositional conditions, and reworking of older deposits (cf. Rittenour, 2008; Forman, 2015). Equivalent dose (De) values for quartz separates were determined on 25–42 aliquots, and the age model used was dictated by overdispersion values. An overdispersion percentage of a De distribution is an estimate of the relative standard deviation from a central De value in the context of a statistical estimate of errors (Galbraith et al., 1999; Galbraith and Roberts, 2012). Overdispersion values ≤20% are assessed routinely for quartz grains that are well solar reset, such as aeolian sands (e.g., Olley et al., 2004; Wright et al., 2011; Forman et al., 2014), and this value is considered a threshold metric for calculation of a De value using the central age model (Galbraith et al. 1999; Galbraith and Roberts, 2012). In this study, one sample (site 7) yielded overdispersion values at one-sigma errors ≤20%, and thus the De was calculated by the central age model. Overdispersion values >20% may indicate mixing of grains of various ages or partial solar resetting of grains; the finite mixing model (Galbraith and Green, 1990) may be an appropriate statistical treatment for such data. The statistical analysis by the finite mixing model (Galbraith and Green, 1990) showed that De distributions for three of the four quartz extracts were skewed negatively. For these quartz separates, there was a clear lower-age population defined by at least six aliquots (Fig. S1 [footnote 1]) defining a De population that yielded an estimate of the depositional age.
River Terrace Correlation, River Incision, and Derived Uplift Rates
We used OSL ages and stratigraphy to correlate the terraces and other mapped surfaces. We calculated incision rates and derived uplift rates from differences in elevation between dated river terrace surfaces and current stream level (Fig. 7A).
Geomorphic Indices and Landscape Analysis
The topographic swath profile A-A′, 350 km long and 25 km wide from the coast inland to Mexico City, can be divided into three sectors (Fig. 2):
A section with the lowest elevation (0–1000 m above mean sea level [a.m.s.l.]) stretches from the coast to La Venta (∼50 km inland). Here, the relief gradually increases from the coast inland to Cacahuatepec town (∼20 km from the coast) at 20 m. Then from Cacahuatepec to La Venta, at ∼50 km from the coast, there is a significant increase in elevation up to La Venta, reaching a peak of 1000 m a.m.s.l. (Figs. 2, 3B).
A second sector in the profile shows a significant drop in elevation north of La Venta (Fig. 2). Then the elevation increases rapidly up to 2500 m a.m.s.l. and then stays at approximately the same elevation. A second significant decrease in elevation (to 500 m a.m.s.l.) coincides with the Mezcala River valley. North of this point the elevation again increases. There is another decrease in elevation (to ∼800 m a.m.s.l.) ∼200 km from the coast.
A third sector of the profile shows a significant increase in elevation up to 3500 a.m.s.l., and further north the elevation stays steady at ∼2500 m a.m.s.l. close to Mexico City (Fig. 2).
The topographic profile overall shows significant changes from the coast inland, from flat to high elevations with sudden decreases in between.
River Longitudinal Profile and Stream-Length Gradient Index (SL)
The river longitudinal profile was constructed to show changes in stream gradient along the river, therefore it shows the actual river length from the coast to inland for a total of 80 km. We used this profile to project the terrace data. Hereinafter, we refer to distances based on the longitudinal profile, i.e., along the river. The channel long profile has in general a concave shape. However, there are also a few breaks or steps from the general concave geometry (Fig. 4A). The first break is ∼28 km from the coast, then there is another break and increase in steepness of the profile at ∼33 km. The channel profile is less steep (close to concave) from 37 km upstream to La Venta.
Stream length-gradient index (SL) values calculated at 1 km intervals along the Papagayo River also show an abrupt transition close to Cacahuatepec (Fig. 4B). Upstream, SL values are generally high, reaching 800, and commonly >400, whereas downstream from Cacahuatepec the SL values are low, ∼100, reaching 200 in only one section. We also estimated average values for two larger sections of the Papagayo River. From a point upstream of La Venta (∼76 km from the coast) to the most prominent bend of the river close to Cacahuatepec (32 km from the coast) the SL is 388, whereas downstream from Cacahuatepec to the mouth of the river the SL is only 110 (Fig. 4B). For the same river sections, the channel gradient is 1.55 m/km and 0.55 m/km, respectively. The channel long profile also shows upstream knickpoints at ∼28, 37, and 76 km from the coast (Fig. 4A).
Valley Shape and Ratio of Valley Floor Width to Valley Height (Vf)
Sectors downstream from Cacahuatepec (sectors A and B; Figs. 3A, 3C) show a wide flat-floored valley in cross-sections constructed every 1 km along the Papagayo River. Upstream from the Cacahuatepec fault almost to the Dos Arroyos fault (sector C; Figs. 3A, 3C), the valley is V-shaped. The valley continues to be narrow upstream to La Venta (sector D; Figs. 3A, 3C).
The width of the river valley and the width/height Vf index along 80 km of the Papagayo River was analyzed from the dam in La Venta to the river mouth (Figs. 4C, 4D). High Vf values (>5.0 and up to 20) were observed from the river mouth to ∼33 km upstream, with a minor decrease (1 < Vf < 5) between 12 and 20 km, suggesting a broad flat-floored valley, typical of regions with low uplift rates. In contrast, upstream from ∼33 km, the valley narrows and the Vf ratio values are low, indicating deep V-shaped valleys common to areas with high uplift rates (Bull and McFadden, 1977; Bull, 1978; Rockwell et al., 1985; Keller, 1986; Keller and Pinter, 2002). Higher values of valley width upstream at ∼76 to ∼80 km from the coast are an artifact related to the effect of a reservoir behind a dam in La Venta (Fig. 4C).
Minimum Bulk Erosion (Ebulk)
Downstream from Cacahuatepec, 32 km inland, values for bulk erosion or minimum amount of material eroded range from 0 to 250 m (eroded-rock column), whereas upstream to La Venta (76 km inland), values range from 250 to 1000 m (eroded-rock column), particularly upstream from the Dos Arroyos fault (Fig. 5). This indicates an apparent fourfold increase of erosion upstream from the Dos Arroyos fault. The La Venta fault has a broad impact on bulk erosion values, reflected in higher values on the western part of the Papagayo drainage basin and north of the La Venta fault than elsewhere (Fig. 5; Fig. S1 [footnote 1]).
We identified multiple levels of fluvial terraces in the Papagayo River basin (Fig. 6). Along the Papagayo River, the spatial distribution of terraces appears to be related to valley floor width and distance upstream. The valley floor width is ∼400–500 m with wide fill terraces that dominate the valley from the coast to ∼9 km upstream. Upstream, ∼11–22 km from the coast, the valley narrows and strath terraces dominate the landscape, showing only narrow strips of filled terraces. From 22 to ∼33 km upstream, the valley widens with wider strath terraces preserved, which are more common in this segment of the river offset by the Cacahuatepec fault. Upstream of ∼35 km, the river valley narrows and has a bedrock channel floor, and strath terraces are absent or rare (Fig. 6A).
We classified river terraces into three groups according to their stratigraphy and age: Q1 is a fill terrace and the youngest; Q2 is older, a strath terrace; and Q3 is also a strath terrace, the oldest (Fig. 6).
Terrace Q1 is a cut-and-fill terrace that is spatially discontinuous in the study area, wide at the mouth of the river, and narrowed upstream with valley constriction (Fig. 6A). It is well preserved and overlain by alluvial deposits, including the most recent flood deposits produced by Hurricane Manuel on 15 September 2013. The thickness of the flood deposits varies according to the location, reaching >1–2 m in some locations downstream of latitude 16°55′N. Flood deposits are mostly composed of sand and in some cases gravel, cobbles, and other debris. The Q1 terraces were also subject to scour and erosion by the hurricane-induced flood event, producing cut terraces in some locations. In most places along the stream, the Q1 terrace is unpaired, reflecting continuous lateral erosion and vertical incision (Merritts et al., 1994).
The Q2 terrace is the most continuous surface in the study area, becoming intermittent upstream from latitude 16°55′N (Fig. 6A). This terrace surface is barren, and bedrock is exposed in some locations after the most recent flood by Hurricane Manuel (Fig. 6C). Similar to Q1, Q2 is covered by recent flood deposits in some locations (Fig. 6B). Q2 terraces are paired in most places along the stream profile.
The Q3 terrace was identified in only one location, ∼12–13 km upstream from the coast, on the eastern river bank (Figs. 6A, 6D). Q3 is a strath terrace with thick alluvial cover (massive sand) of at least 3 m thickness.
Terrace Correlation and Longitudinal Profiles
Terrace elevations and composition were determined along an 80 km stretch from the coast upstream and to close to La Venta (Fig. 6). Terraces were correlated on the basis of stratigraphy, type (fill, cut, and strath terraces), and age. Strath bedrock elevations were used to trace the long terrace profiles. The Q1 fill terrace and Q2 strath terrace long profiles parallel the current channel topography (Fig. 6E).
The channel and terrace cross-sections show an up-catchment increase in slope. The channel gradient was relatively low at 0.55 m/km from the river mouth to ∼30 km, and increased to 1.55 m/km from 30 to 40 km. The steepest gradient, 2.5 m/km, was confined to the stretch upstream from the La Venta fault (Fig. 6E).
Surface Terrace Ages and Correlation to Sea Level
We collected samples from four locations for OSL dating from terraces Q1 (site 9), Q2 (sites 2 and 7), and Q3 (site 8) (Fig. 6A), where the strath terraces were covered by a veneer of fine-grained alluvial sediments. We present equivalent dose (De) values and other chronologic details for those samples derived by different age models in Table 1. At sites 7 and 9, OSL dates provided historic ages (470 ± 30 yr and 65 ± 5 yr, respectively). The sample at site 7, which appears as a thinner, better-sorted deposit, may be local colluvium or eolian material (Fig. S2 [footnote 1]). The sample at site 9 was collected just below an overbank deposit, and it appears that we sampled part of the recent storm deposit from Hurricane Manuel (deposited in 2013). At site 2, the sample was collected at 1.20 m depth from the surface, just above the bedrock and regolith, and yielded an OSL age of 4670 ± 30 yr. Another sample was collected 1.8 m below the surface at site 8 on a strath terrace that was covered by overbank and recent sediments deposited during the 2013 flood. Fine sand to massive silty sediments covered the unexposed bedrock at site 8. The OSL analysis indicated an age of 18,635 ± 1325 yr. We used the dates and elevations of terraces Q2 and Q3 and correlated them with the Holocene and Pleistocene sea-level curves, respectively, to account for eustatic sea-level changes (Lambeck et al., 2014).
On the basis of the ages obtained, we calculated incision rates, ranging from 0.5 mm/yr downstream (Figs. 6E, 7A) to 2.4 mm/yr (Figs. 6E, 7A) close to Cacahuatepec (∼32 km from the coast). The stratigraphic correlation of strath terrace Q2 and the pattern of increase of incision rates upstream suggested an incision rate of ∼4.9 mm/yr for the Q2 terrace close to La Venta, ∼76 km from the coast (Figs. 6E, 7A). The upstream increase in incision rates estimated from OSL ages correlates well with the pattern of an upstream increase in exhumation rates (Fig. 7B), which are calculated from previously published (U-Th)/He thermochronology ages in the study area (Ducea et al., 2004). As further detail about the profile discontinuity, the valley is wide downstream and very narrow upstream. A discontinuity in the longitudinal river profile ∼30 km from the coast, near Cacahuatepec, also coincides with an increase in incision rates.
Incision Rates, Regional Faults, and Geometry of Subduction
The apparent increase upstream in incision rates, from 0.5 mm/yr to ∼5.0 mm/yr, coincides with increased bulk erosion upstream. The channel profile also indicates upstream knickpoints that suggest the adjustment of the river to either climatic change or tectonic uplift (Whittaker, 2012). We evaluate below potential mechanisms for the adjustment of the river to isostatic rebound, climate and/or sea-level changes, lithological changes, and/or tectonic uplift.
A number of factors could influence the spatially variable incision rates and knickpoints along the Papagayo River. First, on this continental coastline far from ice sheet loads, the mode of glacial isostatic deformation of the Fennoscandia shield clearly indicates that deformation is fully compensated on the regional scale, although not globally (Mörner, 2015). Therefore, global glacial adjustments of coasts and seafloors 1000 km beyond the ice sheet margin, as in western Mexico, are small if present; thus, the model proposing global adjustments (Peltier et al., 2015) may be unsustainable (Mörner, 2015). The isostatic loading models also predict high mid-Holocene sea levels in the Pacific Ocean (Peltier, 2004). This also seems not to be in accordance with observations (Mörner, 2015). Thus, we are careful not to apply model reconstruction and prediction based on proposed global loading mechanisms. This is particularly important in areas dominated by tectonics (Mörner, 2015). Second, some studies report that local climate changes in this region were not significant in the late Holocene (Metcalfe et al., 2000; Vázquez-Selem and Heine, 2004). However, others suggest that in the Central Balsas basin, north of our study area, during the late glacial period (14,000–10,000 yr B.P.), lake beds were dry and the climate was cooler and drier (Piperno et al., 2007). Also, climatic drying episodes are reported between 1800 yr B.P. and 900 yr B.P. (Piperno et al., 2007), suggesting that increased incision rates are not climate driven. Third, the local sea level has remained constant since the mid-Holocene, ca. 6000 yr B.P. (Curray et al., 1969; Sirkin, 1985), and the projected global model roughly coincides with the local data (Peltier et al., 2015). Fourth, no remarkable lithological changes occur along the Papagayo River, where the lithology is homogenous and mainly granites and gneisses. Thus, we can exclude the above-mentioned processes as mechanisms for the observed spatial changes in incision rates and presence of knickpoints on the Papagayo River channel profile.
A possible correlation between knickpoints in the river profile and the presence of potentially active faults is apparent. From the topography and river profile, it appears that the La Venta and Cacahuatepec faults are the most important faults and likely the most active. Marked changes in the elevation along the Papagayo River that coincide with the Cacahuatepec fault and Agua del Perro fault, and to a lesser extent the Dos Arroyos fault, also spatially coincide with first-order geometric changes in the plate-interface dip, a bend that marks the change from a steep to a flat plate interface. It also appears that the La Venta fault coincides with a difference in forearc seismicity to the north and south of this structure (Fig. 8), i.e., a decrease in seismicity from this fault to the north (Kostoglodov et al., 2015; Pérez-Campos et al., 2015). It seems that this fault is a barrier in terms of the overriding plate seismicity down to the megathrust. A consequence is that slip partitioning is occurring along the Guerrero sector of the forearc. Thus, for this fault, changes may be related to seismicity and slab geometry (Figs. 2, 8). The importance of the La Venta, Cacahuatepec, and Agua del Perro crustal faults is emphasized by plate interface activity (either stick slip or SSEs) that can trigger earthquakes and/or deformation on these faults (Wang et al., 2013; Melnick et al., 2012a).
The observed spatial increase of incision rates upstream may be related to (1) the presence of inferred left-lateral strike-slip faults with a vertical component (Cacahuatepec fault, Dos Arroyos fault, and Agua del Perro fault segment of the La Venta fault), and/or (2) the geometry of the subducting slab.
Late Quaternary Incision Rates versus Late Cenozoic Exhumation Rates
Incision rates calculated here for the Papagayo River basin transect are comparable to the incision rates observed at other subduction zones (Personius, 1995; Burbank et al., 1996; Pazzaglia and Brandon, 2001; Wegmann and Pazzaglia, 2002; Farías et al., 2005; Yamamoto, 2005; Matsu’ura et al., 2008; Marushima and Ishiyama, 2009; Berryman et al., 2010; Hall et al., 2012) (Table 2). The exhumation rates show a general tendency to increase inland (Fig. 7A). The exhumation rate (0.1 to ∼0.3 mm/yr) is ∼10× lower than our calculated late Quaternary incision rates. Wegmann and Pazzaglia (2002), studying the Clearwater Basin in the forearc of the Cascadia subduction zone, obtained similar inconsistency between low values of longer-term exhumation rates and much higher late Quaternary incision rates, but find similar increases inland. The work of Burbank et al. (1996) on the northwestern Himalayas also shows much lower values for long-term exhumation erosion rates compared to Quaternary incision rates. Thus, overall rates referring to different ages have a similar spatial pattern; however, the values are lower for the late Cenozoic erosion rates at the 106 yr time scale, whereas the younger, late Quaternary values, at a scale of a few thousand years, yield higher incision rates. This disparity in apparent rates measured on different time scales has been addressed previously regarding paleontology and episodic sedimentation (Sadler and Strauss, 1990). Most apparent geologic rates are lower at longer time scales because rates are not constant or linear but episodic, as with a fault slip with an earthquake. As the time scales shortened, rates commonly increase. This appears to be a common effect for sedimentation rates, speciation rates, and incision and erosion rates. We discuss this question further in the following section.
Short-Term (GPS Instrumental) Rates of Vertical Deformation
What we know about the earthquake cycle of deformation, including associated SSEs and megathrust in this sector, is limited to instrumental data. A large earthquake on the megathrust in the northwestern Guerrero seismic gap has not occurred since the Ms 7.8 event on 16 December 1911. However, sectors to the west-northwest and east-southeast of the Guerrero seismic gap have experienced large events, such as the earthquakes on 19 September 1985 in Michoacán (Mw 8.1), and to the east-southeast of the northwestern Guerrero seismic gap on 28 July 1957 (Ms 7.8) and 15 April 1907 (Ms 7.9). Coseismic uplift of the coast up to 1 m was recorded by observing the mortality of intertidal organisms after the 1985 earthquake (Bodin and Klinger, 1986). When the uplifted coastal sites were revisited after a survey in the 1990s, little (<10 cm) to no postseismic changes were found (R. Corona, personal commun. reported in Suárez and Sánchez, 1996). A maximum of 8.5 cm change and subsidence of the coast relative to preseismic time was recorded for the 1985 earthquake. These observations suggest that the re-leveled sites experienced little deformation, <10 cm subsidence up to 80 km inland, and are inland of the down-dip termination of the rupture zone. The dislocation models predict vertical deformation of the upper plate from coastal coseismic uplift to minor subsidence inland (Suárez and Sánchez, 1996).
We thus hypothesize that large megathrust earthquakes might produce coseismic uplift of the coast that decreases inland. To test this, we modeled coseismic vertical displacements, calculated using STATIC1D software (Pollitz, 1996) for three possible rupture scenarios in the Guerrero section of the Middle America Trench. The model assumes a 175 × 80 km rupture plane that is striking 290° and dipping 14°. Displacements are estimated along a 300 km transect from the coast that is perpendicular to the trench for a M 8.0 (2.5 m slip), M 8.3 (7.3 m slip), and M 8.5 (14.4 m slip) earthquake (Fig. S3 [footnote 1]). It is likely that during the postseismic and interseismic phases, only part of the observed uplift along the coast and inland is recovered. A part of this vertical elastic deformation may not be recovered in the postseismic phase, potentially leaving a signature in the topography, reflecting warping of the crust and contributing to the observed uplifted terraces inland. Thus, the development of the current topography and forearc uplift has necessitated many earthquake cycles over thousands and millions of years.
Alternatively, short-term published data derived from GPS and interferometry show uplift of ∼70 mm near the coast, which decreases inland to 2–10 mm near La Venta during the two largest SSEs. Walpersdorf et al. (2011) reported 20 mm of uplift at the Guerrero coast, 46 mm at 20 km inland, and zero at 90 km inland. Deformation between SSEs is characterized by subsidence at the coast (from 10 to 20 mm/yr) and inland (6 mm/yr) (Vergnolle et al., 2010; Walpersdorf et al., 2011; Cavalié et al., 2013). This spatial deformation has similarities between the late Holocene and the late Cenozoic. However, assuming that uplift of ∼40 mm occurs every 4 yr, and during the inter-SSE period there is a 10 mm/yr subsidence, the 40 mm uplift produced by one SSE is compensated in one SSE cycle. Thus, vertical deformation produced by SSEs has no apparent significance in the long-term (Holocene and late Cenozoic) crustal deformation.
Spatial and Temporal Variations of Vertical Deformation
The pattern of long-term (late Quaternary and late Cenozoic) tectonic uplift based on our results along an 80-km-long profile across the forearc is consistent with increasing rates from the coast to inland. However, late Quaternary rates (0.5, 2.4, 4.9 mm/yr) are 10× higher than those for the late Cenozoic (0.1–0.3 mm/yr) (Figs. 7A, 7B). A plausible explanation for this difference might be that for episodic events on centennial to millennial time scales, the signatures of individual events are captured, thus showing higher rates of deformation, whereas on 105 to 106 yr time scales signatures of several events are recorded and the time between events is longer (i.e., deformation during the interseismic time, SSEs, and several earthquake cycles) (Sadler and Strauss, 1990). Rates that are invariant on different time scales reflect less dominance of individual events and continuous processes such as SSEs. However, we currently still don’t know the rates of SSEs on a scale of tens of thousands of years.
Vertical crustal deformation produced during SSEs apparently does not leave a permanent signature on the topography, because most of the instrumentally recorded uplift produced by SSEs seems to be canceled during the interseismic period between SSEs. Therefore, when analyzing long-term effects of SSEs and inter-SSE periods on crustal deformation, their components over several events and periods should be integrated into a secular uplift rate to have significance (Fig. 9B).
What is then the main process inducing long-term (late Quaternary and late Cenozoic) crustal deformation? We hypothesize that not one but several multiple potential mechanisms contribute in different degrees and at different stages of the subduction earthquake cycle to the long-term deformation:
Deformation associated with moment release by large magnitude earthquakes in a subduction zone contributes to warping of the overriding plate that results in additions to (non-volcanic) topography of the forearc over several earthquake cycles (Dragert et al., 1994; Klotz et al., 2001; Melnick et al., 2006; Simpson, 2015; Wesson et al., 2015). A fraction of deformation produced during megathrust earthquakes, i.e., uplift of the forearc, as demonstrated by our models of vertical coseismic deformation, is not recovered and can be accumulated through time (Fig. 9A; Fig. S3 [footnote 1]).
The steady state of rock creep better explains the effects of mantle viscoelastic relaxation and tectonic erosion during the interseismic period. The viscoelastic response to megathrust earthquakes compensates for the short-term interseismic and coseismic deformation in the early interseismic stage (transient phase) (Wang et al., 2012; Hashima and Sato, 2017). Later, it is the locking effect that governs the short-term interseismic stage and causes the short-term deformation on land assisted by the steady slip, below the rupture area, and rock creep during the long term causing deformation during the interseismic period in a broader zone (steady state) (Fig. 9A). Tectonic erosion partly explains the long-term uplift by landward movement of arc topography (Hashima and Sato, 2017).
Although we know less about deep earthquake (intraslab) contribution to the long-term deformation, recent observations of the 19 September 2017 (Mw 7.1) Puebla deep earthquake and interferogram data suggest coseismic subsidence of 6 cm in a 50 km radius around the epicenter of the deep earthquake (Centre for Observation and Modelling of Earthquakes, Volcanoes and Tectonics–Natural Environment Research Council, Twitter post, 25 September 2017, https://twitter.com/NERC_COMET/status/912330067556134912/photo/1) (Fig. 9B). However, we don’t know yet the effect of these kind of events on the long-term expression of topography.
The available observations suggest that long-term late Quaternary and late Cenozoic rates of deformation reflect partially accumulated (inelastic) deformation produced by earthquakes over multiple earthquake cycles, incomplete coseismic release of strain that was accumulated during the interseismic stage, and (inelastic) crustal deformation produced by long-term subduction (Fig. 9C).
Although there is progress in deciphering the surface and deep process relations, such as the main processes inducing long-term (late Quaternary and late Cenozoic) deformation, and short-term versus long-term rate variations, our knowledge is still at its infancy and there are several questions that need to be addressed in future studies. Geodetic and seismic observations along the Mexican subduction zone are still insufficient to draw conclusions on short-term deformation, and the record of seismic data is insufficient to allow a better understanding of how deformation is distributed among various modes of fault slip. SSEs have been only recently discovered along the Mexican subduction zone to allow a better understanding of long-term subduction process, and even at the short time scale, our knowledge is not yet sufficient to explain the role of SSEs in the long-term subduction earthquake cycle.
We have compared our results on the long-term (late Quaternary) rates of deformation with the late Cenozoic and the short-term rates of deformation and their respective contribution to permanent deformation of the forearc in the Guerrero sector of the Mexican subduction zone. Among the different analyzed parameters and processes discussed, we concluded that multiple mechanisms contributed in different degrees and at different stages of the subduction earthquake cycle to the long-term deformation of the forearc (coseismic uplift by megathrust earthquakes over multiple earthquake cycles, and warping by rock creeping in the interseismic period). However, we did not find the significance of SSEs in the long-term permanent deformation and vertical deformation of the forearc. This offers an opportunity to better understand processes at depth and their contribution to the forearc permanent deformation among different modes of fault slip of subduction zones.
Ramírez-Herrera acknowledges funding by Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (PAPIIT) IN109117, Programa de Apoyos para la Superación del Personal Académico de la UNAM (PASPA)-2015, and Secretaria de Educación Pública (SEP)-Consejo Nacional de Ciencia y Teconología (CONACYT) grant 129456. Thank you to the Berkeley Seismological Laboratory, Department of Earth and Planetary Science, University of California Berkeley, for use of facilities during a stay as a Visiting Scientist. Ann Grant helped editing the final text manuscript. Karen González helped with Figure 9 drafting.