We present the first example of a short-term uplift transient within the long-term stagnant forearc of the Arica Bend, northernmost Chile. We date a storm deposit embedded into fan delta sediments, which were deposited close to sea level and are now located ∼40 m higher. The radiocarbon age is ca. 10 ka, yielding an average uplift rate of ∼4 mm/yr. During the Holocene, the section has been dissected by the Lluta River, the long profile of which shows geomorphic response to coastal uplift. However, the coast of the Arica Bend is characterized by million-year-scale relative subsidence and current global positioning system surveyed quiescence. We therefore interpret the inferred uplift signal to represent a transient increase in uplift rates throughout the Holocene. Similar short-term transients have been globally documented in relation to various subduction zone settings, although their causes still remain a matter of debate and need further investigation. We propose that periods of increased surface uplift may also occur in long-term quiescent forearcs, possibly due to temporal variations of plate interface properties (i.e., strength and/or seismic locking) that induce thickening beneath the coast.

Raised terraces and paleoshorelines are widely used to identify uplifted coasts and estimate their uplift rates over a variety of time scales (e.g., Regard et al., 2010). In forearc settings, the comparison between uplift rates inferred from geological records and those directly measured through global positioning system (GPS) monitoring allowed us to unravel the complexity of subduction kinematics, where rapid uplift and subsidence due to stress accumulation and release throughout the seismic cycle are thought to result in slow forearc uplift over a million-year scale (Savage, 1983; Saillard et al., 2017). During the past two decades, however, periods of vertical displacement at rates mismatching both long-term and short-term observations have been directly measured through GPS monitoring of subduction zones worldwide (Liu and Rice, 2007, and references therein) and inferred from the stratigraphy of emerging forearcs (Natawidjaja et al., 2007; Sawai et al., 2004). These transient pulses of uplift (Liu and Rice, 2007), apparently misaligned from the megathrust earthquake cycle, disclosed a further complexity of the link between seismicity and surface uplift. Therefore, reporting undocumented uplift transients from diverse settings is crucial to better understand this process.

Here we present for the first time a well-dated case, where transient uplift pulses have produced several meters of surface displacement in a forearc setting undergoing long-term tectonic stagnation. We report evidence for an uplift pulse at the Arica Bend, northernmost Chile (Fig. 1), which is the only coastal area of the Central Andean forearc that has undergone million-year-scale quiescence (Madella et al., 2016; Mortimer and Saric, 1972), as evidenced by the occurrence of the ∼2.7-m.y.-old Lauca ignimbrite found at 60 m a.s.l. (above sea level) ∼3.2 km from the coast (Fig. 2) (Wörner et al., 2002). Long-term tectonic stagnation at Arica contrasts with long-term uplift farther north and south, where the Coastal Cordillera evidences forearc emergence at a million-year time scale (Fig. 1) (Melnick, 2016; Regard et al., 2010; Saillard et al., 2009). Next to Arica, we investigated the sedimentological architecture of a raised conglomerate suite, which we will interpret to have been deposited in a fan delta environment. We dated these deposits through 14C applied to wood fragments embedded in the sediments, and investigated the geomorphic response of the Lluta River (Figs. 1 and 2) to possible tectonic perturbations. Results suggest that, despite both million-year-long (Wörner et al., 2002) and GPS-inferred (Melnick, 2016) quiescence, the Arica Bend underwent transient uplift during the Holocene, uplift that we propose resulted from a period of intensified forearc thickening through upper plate fault activity.

Sedimentology and Chronology

We mapped a conglomerate suite located along the Lluta River 2.8 km away from its outlet (Figs. 1 and 2). The sedimentary fabric of the deposit was described and samples for 14C dating were collected. The average postdepositional uplift rate is calculated through the equation
where the uplift rate u equals the height difference between current elevation hc and paleoelevation hd, minus sea-level elevation sd at the time of deposition td, divided by td. The depositional environment of the deposit has been constrained through its sedimentological properties (backshore; see following), and hd is considered to equal that of today’s backshore environment, which we estimate based on a 4-km-wide swath profile across the Lluta River delta (Fig. 2), extracted from the ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model) with a 30-m-resolution grid. The vertical accuracy of the DEM for barren, poorly urbanized, open spaces is considered not to exceed 2 m of standard deviation (ASTER GDEM Validation Team, 2011). We therefore neglect this vertical error, because a much larger uncertainty for hd arises from the 4-km-wide swath’s standard deviation (Fig. 2). The site location and height were measured with a handheld GPS, for which we assume an uncertainty of ±5 m. Sea-level elevation at the time of deposition was estimated using the Holocene sea-level curve proposed by Dura et al. (2016, and references therein), which has been determined for the locality of Quintero, central Chile (33°S). Thus far, it is the best-constrained quantification of relative sea-level variations for the western central Andes at the Holocene time scale.

Radiocarbon dating was performed with a MICADAS (Mini CArbon DAting System) accelerator mass spectrometer (AMS) at the Laboratory for the Analysis of Radiocarbon with AMS (LARA) of the University of Bern, following the protocol described in Szidat et al. (2014). Analytical results are reported in uncalibrated 14C ages, from which the average calendar age was inferred using the SHCal13 calibration curve (Hogg et al., 2013).

Stream Profile Analysis

Once a stream has fully adjusted its long-stream profile to external perturbations, the river gradient S at any point of its course is best represented by a power law function (Flint, 1974; Whipple and Tucker, 1999):
where the coefficient ks and exponent θ denote the steepness and concavity index, respectively, and A is the drainage area. Accordingly, any deviations from this graded profile can be used to identify transient effects related to base-level fluctuations, tectonic perturbations, or climatic changes (Wobus et al., 2006). We applied these concepts to the lowermost 80 km of the Lluta River, following the calculation procedure reported by Valla et al. (2011), thereby predicting a steady-state equilibrium profile for the identical source elevation, river nodes spacing, and catchment area as the observed profile. These parameters were extracted from the ASTER GDEM. The best-fit outlet elevation and concavity were determined through a least-square minimization of the elevation difference between observed and predicted profiles (Δh). We set the source elevation at 1400 m a.s.l. (80 km from the outlet), thereby excluding the main knickzone of the Lluta River (e.g., Veit et al., 2016). This approach has been chosen in order to detect perturbations that would otherwise be obliterated by the prominent convexity of the knickzone. In addition, the Lluta River flows entirely on a gravel bed for its lowermost 80 km, within a deeply incised canyon. Along this reach, the stream is almost devoid of tributary channels, which possibly alter the stream erosive power due to along-strike changes of sediment supply (Sklar and Dietrich, 1998). The analyzed reach, moreover, is located in a desert environment with constant hyperarid to arid climate conditions throughout the Holocene (Grosjean et al., 2003). Because of these conditions, we consider the lowermost 80 km of the Lluta River as suitable for the outlined stream profile analysis, and we thus anticipate vertical deviation from a steady-state profile to represent the geomorphic response to tectonic forcing along the channel. We further corroborate the predicted deformation pattern through field observations of the channel morphology.

Sedimentology and Chronology

A raised conglomerate unit occurs 2.8 km from the coast of Arica along the Lluta River, in a gravel pit located within Holocene fluvial deposits, ∼1 km inland from the limit between the current littoral deposits and Pleistocene fluvial deposits (Fig. 2). The analyzed section (Fig. 3) comprises a 1-m-thick clast-supported conglomerate bed (unit 1) with a massive-bedded fabric and well-rounded clasts. The pebble b-axes are between centimeters and decimeters long; they are imbricated with an eastward dip. Interbedded coarse-grained black-brownish sand lenses in unit 1 are thin (decimeters) and parallel to wavy laminated at the millimeter scale. Following a sharp erosive contact with decimeter-deep scours, a 1-m-thick light grayish-yellow suite of fine- to coarse-grained sandstone beds (unit 2) is ∼50 m long and gently dips 2° toward the west, where it disappears at the base of the section. To the east, this unit is cut by the modern land surface (Fig. 3A). The sands of unit 2 display an alternation of tabular and trough cross-beds, wavy laminations with meter-wide wavelengths, hummocky cross-stratifications, asymmetric ripples and sandwaves with amplitudes of as much as 10 cm and wavelengths of as much as 50 cm, and water-escape structures (Fig. 3). In addition, unit 2 contains millimeter-scale shell fragments and interbedded decimeter-thick and meter-wide ash-rich sand lenses with a massive-bedded fabric. Embedded decimetric wood fragments are aligned parallel to the bedding surfaces (Fig. 3). The encountered sedimentary structures suggest rapid deposition at the transition between the lower and upper flow regimes, with inland current direction (Fig. 3). Unit 2 is cut and topped by unit 3, a 1-m-thick conglomerate displaying an analogous fabric to unit 1. Three wood fragments sampled from unit 2 (Fig. 3) yielded three consistent ages within error, which averaged and calibrated equal 10.440 ± 0.650 cal ka B.P. (Table 1).

The ensemble of sedimentary structures observed in units 1 and 3 (Fig. 3) suggests fluvial deposition by the Lluta River, according to seaward imbrications, scours, and laminated sand channel-fill structures. In contrast, the eastward paleocurrent and the occurrence of shell fragments suggest that sands of unit 2 have been deposited during a storm, possibly forming a washover fan within the backshore environment (e.g., Wang and Horwitz, 2007). In particular, the lack of any freshwater shells in all the mapped fluvial deposits of the study area (García et al., 2004) corroborates our interpretation of the marine origin of unit 2. Here we acknowledge that discerning between storm and tsunami deposits can be complicated (Kortekaas and Dawson, 2007). The preservation potential of tsunami deposits along the Peru-Chile margin, however, has been shown to be extremely low, due to the hyperarid climate that facilitates eolian erosion (Spiske et al., 2013a). The well-sorted sedimentary fabric, the thickness of the suite, the lack of pebbles and rip-up clasts, and especially the occurrence of hummocky cross-stratifications and current ripples in unit 2 suggest deposition during a storm rather than a tsunami. Among many possible features, tsunami deposits tend to display (1) bidirectional imbrication of pebbles and cross-laminae in response to run-up and backwash, (2) boulders and rip-up clasts evidencing very high energy, (3) normally graded sand resulting from suspension transport, and (4) poorly sorted fabric (Kortekaas and Dawson, 2007; Spiske et al., 2013b), all inconsistent with our observations. In addition, storm deposits exhibiting a fabric analogous to that of unit 2 were mapped nearby (García et al., 2004). For these reasons, although we cannot completely discard the possibility of tsunamigenic deposition, we interpret unit 2 to constitute a washover fan, the paleoelevation of which ca. 10 ka was 12.0 ± 7.0 m a.s.l. (Fig. 2). Accordingly, the ensemble of observations, including fluvial conglomerates with a backshore interbed, suggests that the analyzed sediments were part of a fan delta environment deposited by the Lluta River. Furthermore, considering 11.1 ± 6.1 m of sea-level rise since time of deposition (Dura et al., 2016), the fan delta currently located at 40 ± 5 m a.s.l. has undergone a cumulative uplift of 39.1 ± 18.1 m. Therefore, the inferred average coastal uplift rate of the Lluta River fan delta equals 3.9 ± 2.0 mm/yr for the past 10440 ± 650 yr, which is a conservative estimate. Here we note that reported wave heights for the largest historical tsunami in the region of Arica are as great as ∼20 m (Spiske et al., 2013b), within the uncertainty of our paleoelevation estimate. The resulting Holocene uplift rate of 3.9 ± 2.0 mm/yr would therefore remain unaltered even considering tsunamigenic deposition.

Stream Profile Analysis

The lowermost 80 km of the Lluta River profile are characterized by an overall slight concavity and a prominent convexity between 3 and 15 km from the outlet (Fig. 4). The observed profile plots above the theoretical steady state at 0–13 km and 48–65 km from the outlet, suggesting relative uplift, whereas it plots below the theoretical steady-state profile along the rest of the considered reach, denoting relative subsidence (Fig. 4). This pattern of relative uplift can be further corroborated through field observations. If neither water discharge nor bedrock lithology vary along the considered reach, narrow valley flats have been shown to correlate with ongoing uplift and incision (Brocard and van der Beek, 2006), while broad floodplains are thought to develop along subsiding reaches. In our case, field observations show that near the coast (0–13 km), the Lluta River is currently incising a 10-m-deep gorge into the dated fan delta deposits (Fig. 2), where the stream profile is characterized by a prominent convexity (Fig. 4). In contrast, between 13 and 48 km the model predicts relative subsidence, where the floodplain is as much as ∼1.5 km wide (Fig. 4). Between 48 and 65 km from the outlet, relative uplift is inferred where the valley flat narrows to <250 m and where perched fluvial gravels and hanging tributary fans are found (Fig. 4). Nevertheless, the interpretation of relative uplift along this latter upstream reach is less well constrained. At ∼54 km, the only tributary catchment of the lowermost 80 km of the Lluta River joins the canyon from the southeast (Fig. 2). This confluence makes the drainage area increase by 13% and produces a concave kink in the predicted profile (Fig. 4). This tributary catchment, however, is at elevations lower than 3400 m a.s.l. (Fig. 2), where seasonal precipitation is strongly inhibited by the Andean rain shadow (Houston and Hartley, 2003). For these reasons, the inference of relative uplift in the reach at 48–65 km might not reflect a local tectonic forcing, but rather the effect of a dry tributary confluence on the model.

In summary, stream profile analysis, channel morphology, and field observations suggest that the Lluta River shows transient geomorphic response to uplift and incision near the coast, subsidence and aggradation between 12 and 40 km from the mouth, and possibly uplift between 48 and 65 km distance from the outlet. These inferences, moreover, are consistent with the 14C age and the sedimentary fabric of the fan delta deposits described here.

Possible Causes for Fluvial Incision Near the Coast

The Lluta River currently flows 10 m beneath the dated paleodelta (Fig. 2), implying an average net incision rate of 1.0 ± 0.1 mm/yr during the Holocene. We note that Dura et al. (2016) proposed a Holocene eustatic curve that indicates >5 mm/yr of sea-level rise from an elevation of ∼−11 m a.s.l. at 12 ka to a maximum of 3 m a.s.l. at 7 ka. According to Dura et al. (2016), sea level dropped steadily to the current base level after 5 ka. This mid-Holocene 3 m highstand likely coincided with marine abrasion of the coast and subsequent marine regression, thereby inducing a wave of incision upstream of the backshore environment. However, this mechanism alone cannot explain the observed 10 m incision or the profile convexity found at km 13 (Fig. 4), because the amplitudes of both morphological changes are much larger than that of the sea-level drop.

Rivers that drain the western Andean forearc are known to respond to El Niño–driven floods by strong sediment pulses. These processes are capable of depositing several meters of fluvial gravel (Bekaddour et al., 2014), temporally altering the channel morphology and obstructing the previously established river course. Although floods are likely to be recorded by conglomerate sheets, we discard such a mechanism alone to be recorded by the fan delta deposit, because the 3-m-thick analyzed section comprises records of both marine and fluvial processes. Based on a chronology of paleosols in the Lluta valley, Veit et al. (2016) showed that in the late Pleistocene, the Central Andean Pluvial Event is likely to have caused a wave of incision 30–50 km further upstream. The phase of valley-floor lowering, however, terminated at 15.4 ka at the latest (Veit et al., 2016), 5 k.y. before the deposition of the fan delta conglomerates at Arica. Accordingly, we cannot use this mechanism to explain the 10-m-deep incision observed ∼3 km onshore from the current coastline. We are thus left with a scenario where the phase of fluvial downcutting was most likely initiated by coastal uplift at the Arica Bend during the Holocene. These mechanisms and the possible implications for the seismic behavior of the region are elaborated in the following.

Possible Drivers for a Holocene Uplift Pulse

In Madella et al. (2016), it was proposed that the growing seaward concavity of the Bolivian orocline may have fostered the initiation of a positive feedback mechanism, where focused trench sedimentation and intensified pore-fluid pressures in the subduction channel would have weakened the plate interface and inhibited coastal uplift of the Arica Bend during Miocene–Pliocene time. For the Holocene, however, our geochronological and geomorphological data suggest the occurrence of uplift with a minimum average rate of ∼1 mm/yr, as suggested by fluvial incision near the coast. The relative elevation of the storm deposit yields an even higher average uplift rate of 3.9 ± 2.0 mm/yr. The Holocene rates, in either case, exceed the 0–0.1 mm/yr estimated for the million-year scale (e.g., Melnick, 2016) by at least one order of magnitude. The decadal GPS signal recorded by the IPOC (Integrated Plate boundary Observatory Chile; https://kg3-dmz.gfz-potsdam.de/ipoc) shows no or very minor vertical displacement at rates between 0 and 0.1 mm/yr, with slight subsidence following the 2014 Iquique earthquake. The bulk signal over an entire seismic cycle, however, displays little or no uplift. This shows that the current rates of surface displacement are strongly related to the seismic cycle of the neighboring regions, and it is not comparable to that estimated for the Holocene. Therefore, our results for the Arica Bend require a mechanism where both million-year-scale and recent quiescence coexist with increased uplift rates in the Holocene. As this uplift pulse clearly deviates from the background tectonic quiescence, we consider it as a transient signal.

Several processes have the potential to cause such a transient uplift of the coast, including subduction erosion with related isostatic unloading, subduction of a seamount, and sediment underplating. These mechanisms, however, are expected to cause uplift signals at time scales longer than 100 k.y. (Mouslopoulou et al., 2016); this is incompatible with the observed change from rapid to no uplift within a few thousand years.

The degree of locking (or seismic coupling ratio) is derived from geodetic measurements, and is expressed as a number between 0 and 1, indicating the fraction of slip that occurs through seismic events (Wang and Dixon, 2004). At the Arica Bend, the plate interface is considered to be less locked than most of the remaining Peru-Chile margin (Li et al., 2015) and to have undergone no historical megathrust earthquake (Dorbath et al., 1990). Accordingly, if convergence does not occur through subduction earthquakes, it has to result in a combination of (1) shortening in the eastern Andes hinterland (which does not affect the Arica Bend), (2) forearc thickening, and (3) aseismic slip. Sawai et al. (2004), for example, proposed that large shallow subduction earthquakes (<40 km depth) might trigger post-seismic creep along the deeper segment of the plate interface (>60 km) situated below the uplifted coast. At the Arica Bend, however, the anomalous subduction seismicity (Dorbath et al., 1990) and the low seismic coupling ratio (Li et al., 2015) constitute the major arguments in favor of a mechanism that does not involve large ruptures.

Based on modeling and compilation of short- and long-term uplift data from various emerging forearcs, Mouslopoulou et al. (2016) showed that rapid uplift transients over a 10–20 k.y. scale are most likely accomplished through temporally clustered upper plate fault activity. Moreover, Mouslopoulou et al. (2016) showed that strain may be localized along major reverse faults, as well as diffused along minor structures. In either case, the absence of active surface faults at the Arica Bend (Madella et al., 2016) suggests activity along blind structures. Instrumental evidence exhibiting the relationship between activation of such faults and underlying megathrust events is poor toward the portion of plate interface at or below the downdip end of seismic coupling, in contrast to splay faulting near the updip end (McCaffrey and Goldfinger, 1995). Nevertheless, recent modeling studies (e.g., DeDontney et al., 2012; Li et al., 2014) suggest that triggering of upper plate faulting from subduction earthquakes may occur under appropriate conditions. Depending on megathrust depth and slip distribution, as well as on the relative location of the upper plate splay faults with respect to the site of megathrust rupture, both normal and reverse faulting have been predicted numerically. Accordingly, reverse faulting would be triggered by megathrust events at or deeper than the location of the branch line, a relationship supported by the occurrence of deep megathrust events in northern Chile (Schurr et al., 2012) and their correlation to topography (Melnick, 2016). In addition, for activation to occur, the frictional strength of such faults has to be very low, and possible transient changes in this friction may temporally decouple activation of megathrust events and upper plate deformation (Li et al., 2014). It is likely, therefore, that all (1) stress rearrangement during megathrust events or (2) downdip drop of frictional resistance of the interface, (3) a temporal variation of the latter, or even (4) a temporal variation of upper plate fault strength, for example, from hydraulic loading (Simpson, 2006), may play a role in the transient uplift pulse observed here. Today’s pattern of interseismic microearthquakes (Comte et al., 1999), although likely not representative of the bulk strain or of a 10 ka transient, suggests that diffused fault activity within the overriding plate beneath the Arica Bend is not uncommon (Fig. 5), in contrast to the lack of such seismicity farther south. Moreover, the observed predominance of contractional focal mechanisms for the deeper upper plate events (Fig. 5) evidences the presence of upper plate faults that are potentially able to produce surface uplift if activated at higher rates for several thousand years, which we propose to have occurred during Holocene time. The inferred contractional seismicity below the coastline of Arica (Fig. 5; Comte et al., 1999) coincides with the location where seismic coupling ends in the downdip direction. This kinematic observation, however, does not correspond to a change in strength (i.e., frictional properties) of the plate interface (Wang and Dixon, 2004). Accordingly, a decrease in locking may be sufficient to cause a stress concentration in the overriding plate activating weak faults, but may be superseded by an unresolved downdip change of frictional resistance. Simulation of forearc deformation using dynamically scaled analogue experiments (Rosenau and Oncken, 2009) has shown that significant, if not most, reverse faulting above the downdip end of seismic coupling, although at constant frictional resistance, occurs in the interseismic period from loading of the upper plate faults. This relationship bears a striking similarity to the current spatial pattern of locking, seismicity, and uplift settings around Arica, i.e., a seismically quiescent period with little or no uplift.

Although the extent of our inferences is limited to the age of the dated marker, we cannot exclude that the uplift transient at the Arica Bend might have initiated earlier than 10 ka, or that similarly short-term transients may have occurred before the Holocene. The geological archives of such older events, however, would have likely been lost, due to low preservation potential and the long-term relative subsidence characterizing this portion of forearc. In addition, the calculated uplift rate of 3.9 ± 2.0 mm/yr can only be considered a minimum estimate, because although we integrate the measured displacement over 10 k.y., we do not know when the uplift pulse might have ended before today’s stagnation.

The sedimentological properties and the geochronological record of the raised fan delta deposits imply rapid Holocene uplift. Such deformation has not been recorded by geological archives in the study area at a million-year scale or by the current GPS data, both of which suggest the occurrence of relative subsidence. We conclude that coastal uplift transients over an ∼10 k.y. scale may occur also in forearcs generally characterized by a plate interface exhibiting low locking and the absence of long-term emergence. Moreover, we propose that these transients may be accommodated by upper plate faulting, thereby indicating the likelihood of temporal variations of either the frictional properties of both the plate interface and the upper plate faults, or of the degree of seismic locking of the plate interface.

Research was financed by Swiss National Science Foundation grant 200020–121680 to Schlunegger. We thank Heinrich Bahlburg for the helpful discussion. We also thank Vincent Regard, Christoph Glotzbach, and the anonymous reviewers who greatly improved this manuscript with their comments.