Sediment routing systems (SRSs) are a critical element of the global response to ongoing climate change. However SRS response to climate forcing is complex, fragmentary, and obscured when viewed over short, human time scales (10−1–102 yr). Over long time scales (>102–103 yr), the aggregated, system-wide response of SRSs to climate forcing can be gleaned with more confidence from the sedimentary record, but the nature and time scales of this aggregated response to abrupt climate change are still poorly understood. Here, we investigate the aggregated temporal response of a SRS in northern Spain to abrupt climate warming at the Paleocene-Eocene thermal maximum (PETM). Our results show that terrestrial sites in northern Spain record a temporal lag of 16.5 ± 7.5 k.y. between the onset of the PETM, defined by an abrupt negative excursion in the δ13C profile, and the onset of coarse-grained deposition. Within the same SRS at a deep marine site 500 km to the west, we observe a temporal lag of 16.5 ± 1.5 k.y. using an age model that is independent of that used for the terrestrial sites. These results suggest that the aggregated, system-wide response of SRSs to present-day global warming—if we take the PETM as an appropriate modern-day analogue—may persist for many millennia into the future.


The Paleocene-Eocene thermal maximum (PETM) is the most informative geological analogue for understanding the impact of rapid (<5 k.y.) and large-magnitude (>4 °C) warming on global hydrology and sediment routing systems (SRSs) (Haywood et al., 2011; Foreman et al., 2012; Carmichael et al., 2017). The PETM is associated with a large negative carbon-isotope excursion (CIE) driven by the release of >4000 Pg of isotopically light carbon into the global carbon cycle (Gutjahr et al., 2017). This event is recorded in the carbon isotope (δ13C) values of calcium carbonate minerals and organic matter from both terrestrial and marine strata (McInerney and Wing, 2011). During the PETM, changes in atmospheric and ocean-surface temperature caused dramatic and globally nonuniform change of the hydrological cycle (Bowen and Bowen, 2008; Carmichael et al., 2017). Studies of mid-latitude areas in the central United States (Foreman et al., 2012; Kraus et al., 2013; Baczynski et al., 2013) and northern Spain (Schmitz and Pujalte, 2007; Manners et al., 2013) indicate increased seasonal precipitation during the PETM, within generally dry climate regimes. In northern Spain, there is a clear association between the CIE and an increase in the amount, and caliber, of detrital material transported to both terrestrial and marine environments (Schmitz and Pujalte, 2003, 2007; Dunkley Jones et al., 2018; Pujalte et al., 2015). A laterally extensive (>500 km2), 1–7-m-thick conglomerate unit (“Claret conglomerate”) marks this change in the terrestrial environment (Schmitz and Pujalte, 2003). Similarly, in Wyoming, USA, a 30-m-thick and laterally extensive channel-belt sandstone (“boundary sandstone”) is concomitant with the PETM (Foreman, 2014); and in Colorado, USA, a similar change in sedimentary architecture is associated with the PETM (Foreman et al., 2012). The short temporal duration of PETM onset (<5 k.y.) (Bowen et al., 2015; Dunkley Jones et al., 2018), and the potential for precise, high-resolution stratigraphic correlation between sections using δ13C data, enable PETM Earth system leads and lags to be resolved at the millennial scale (Knight and Harrison, 2013). The preservation and exposure of terrestrial and deep-marine PETM sequences in the Tremp-Graus Basin of northern Spain (Schmitz and Pujalte, 2003, 2007; Domingo et al., 2009; Manners et al., 2013; Dunkley Jones et al., 2018) allows the study of system-wide response of a SRS to rapid and abrupt climate warming. We use these sequences in northern Spain (Fig. 1) to explore the aggregated, system-wide response of a PETM SRS.


To quantify the temporal relationship between the onset of the PETM climate perturbation and the response of the Spanish PETM SRS, we first assume that a close coupling existed between the onset of the CIE and the PETM climate pertubation. This is supported by high-resolution studies of expanded PETM sections from the New Jersey margin (USA) that find no significant lag between δ13C and δ18O records (Zeebe et al., 2016). Here we focus on quantifying the lag time (tlag) between the onset of the CIE and the onset of coarse-grained detrital deposition (OCD) in a study section. To do this, we first calculate a mean sedimentation rate for each section using the latest estimates of PETM CIE duration (Westerhold et al., 2018) and the measured stratigraphic thickness of the CIE at each location. To account for uncertainty in the precise identification of the CIE “core” (∼100 k.y. duration) versus “core and recovery” (∼180 k.y.) in the study sections, we calculate lower and upper estimates for sedimentation rate using both of these durations (Westerhold et al., 2018).

Our approach should give minimum tlag values for two reasons. First, sedimentation rates in most sections increase substantially during the core of the PETM (Dunkley Jones et al., 2018; John et al., 2007), and applying these higher sedimentation rate values to the onset phase will minimize calculated tlag values. Second, the basal surface of the Claret conglomerate in the Spanish sections truncates underlying strata that accumulated during the onset interval of the PETM. Therefore, the preservation of a stratigraphic offset between the onset of the PETM CIE and the base of the Claret conglomerate is strong evidence for a temporal lag of SRS response. On this basis, the tlag values calculated here are almost certainly minimum estimates. We note that, for our purposes, the calculation of tlag is not compromised by the Sadler effect (Sadler, 1981).

To further explore the stratigraphic response of SRSs to abrupt hydrological change, we utilize a one-dimensional (1-D) sediment transport model (see the GSA Data Repository1). The model solves for the change in topography due to the transport of sediment downslope, which is a function of both local slope and surface water discharge (Armitage et al., 2016). Precipitation rate is increased from a baseline of 0.5 m yr−1 following a box profile with a time to peak precipitation of 5 k.y. reflecting the time scale of PETM onset (Westerhold et al., 2018).


The analysis shows that the Tendrui, Claret, and Campo sections exhibit a stratigraphic offset, but that the Esplugafreda section does not (Fig. 2). The calculated durations of the stratigraphic offset (tlag) are: Tendrui, tlag = 13–24 k.y.; Claret, tlag = 9–16 k.y.; Esplugafreda, tlag = 0 k.y. The lower rates of sediment accumulation at the Esplugafreda section (0.1 m yr−1) has led to a reduced level of stratigraphic completeness (Sadler, 1981; Straub and Esposito, 2013), inhibiting proxy record preservation and precluding the identification of a resolvable tlag (Foreman and Straub, 2017). We tentatively estimate the Campo tlag to be 20–36 k.y.

A new age model for the deep marine segment of this SRS at Zumaia (Fig. 1), alongside a record of detrital mass accumulation rate (MARD) (Dunkley Jones et al., 2018), allows for a unique comparison of the terrestrial and marine response to the same event. The MARD data from the Zumaia section show an immediate response to the CIE, increasing from 1 g cm−2 k.y.−1 to 3 g cm−2 k.y.−1 over the first 5 k.y. (Fig. 2E). This elevated value is maintained for 10 k.y. before increasing abruptly to 7 g cm−2 k.y.−1 over a period of <4 k.y., an increase that lags behind the CIE by tlag = 15–18 k.y., while maximum MARD values are attained 25–30 k.y. after the onset of the CIE (Fig. 2E). The immediate response of MARD to the CIE is due to the greater advective length scales or transportability of finer-grained material (Ganti et al., 2014). However, the synchroneity of the arrival of coarse material at the terrestrial sites and the increase in MARD at the deep marine site suggests that a single causative mechanism is responsible for the observed time lag. The data presented above strongly suggest that SRSs may take ∼104 yr to respond to abrupt, large-magnitude climate change.

The results of the 1-D sediment transport model demonstrate that a tlag can be reproduced (Fig. 3). The greater the precipitation increase over the 5 k.y. onset duration, the shorter the value of tlag (Fig. 3). The truncation of individual grain size profiles (e.g., Fig. 3B) signifies sediment bypass and non-deposition (i.e. a time gap) of the coarsest grain-size populations. We suggest that this phase of bypass is recorded in northern Spain as a 20–30-m-thick succession of fine-grained floodplain sediment that directly overlies the Claret conglomerate. The model predicts tlag values of 0 k.y. at 50 km distance from the origin, 10–35 k.y. at 100 km, and 45–85 k.y. at 150 km (Fig. 3D), which is not dissimilar to the field estimates of tlag. We acknowledge that these model results are approximations.


Stratigraphic sections from the terrestrial and deep marine segments of the northern Spain PETM SRS show strong evidence for a time lag between the onset of the PETM CIE and the onset of coarse-grained deposition. The results of a 1-D sediment transport model support this. The observed tlag was generated by the internal response of the SRS; so what was this mechanism? It is possible that the natural avulsion of river channels could generate a stratigraphic offset, and so tlag, but the tlag value would be spatially variable. The maximum value of an avulsion-related time lag can be approximated by a compensation time scale, Tc = h/r, (where h is channel depth and r is aggradation rate). This defines the time window over which sediment can be delivered to the majority of a river valley width (Straub et al., 2009; Wang et al., 2011), which is calculated to be Tc = 6.5 ± 3.7 k.y. for the sections in northern Spain. This value of Tc represents a maximum value, given that the time taken for a channel to visit 95% of the river valley width could be of the order of ∼0.25Tc (Straub and Esposito, 2013), and channel mobility increases with increasing sediment flux (Wickert et al., 2013), and with a reduction in vegetation cover; both of which are experienced by PETM landscapes in the subtropics (Schmitz and Pujalte, 2007; Foreman, 2014). We note that Tc = 6.5 ± 3.7 k.y. also represents a minimum resolvable lag time using geochemical proxies (Foreman and Straub, 2017). Our preferred mechanism for the generation of tlag is a delay in sediment transport from source area(s) to down-system locations as transport slopes adjust to the increased transport capacity of the fluvial system. This mechanism can explain a SRS-wide response. However, given the inherent complexity of landscapes and their nonlinear response to external forcing, it is possible that mountainous landscapes can act as a buffer to hydrological change (e.g., Knight and Harrison, 2013). Once this “buffering capacity” threshold is exceeded, then whole-scale landscape response, and associated increased sediment flux to areas down-system, ensues. Given this scenario, we cannot completely rule out the possibility that the time it takes to reach this threshold is manifested as a tlag or is an important contributor to it. This requires further work.

As a comparison, tlag values were calculated for two well-studied sections in the Piceance Creek basin of western Colorado (Foreman et al., 2012) and the Bighorn Basin of Wyoming (Foreman, 2014), giving values of tlag = 22–40 k.y. and 14–25 k.y., respectively (see the Data Repository). The presence of a tlag in PETM sections in the United States and in northern Spain might suggest a common mechanism of generation. However, we advise slight caution on this because the value of tlag will depend on the input grain size of the sediment supply, location within the basin, and local basin and earth-surface conditions that dictate the degree of stratigraphic completeness (Foreman and Straub, 2017).

Given the lag times present in many of the studied PETM sections, our results suggest that it may take up to 15 k.y. for SRS to achieve a total and aggregated, system-wide response to modern abrupt climate forcing. This time scale probably represents an upper estimate given that anthropogenic carbon release rates are an order of magnitude greater than that associated with the PETM CIE (Zeebe et al., 2016). The aggregated, system-wide response time scale described here is akin to the time scale necessary for “system clearing” (Jerolmack and Paola, 2010; Foreman et al., 2012) to occur in an individual SRS following abrupt climate forcing. This temporal and spatial character of SRS response is a function of the nature of climate forcing and the size and sensitivity of the SRS (Jerolmack and Paola, 2010; Knight and Harrison, 2013). From a stratigraphic perspective, accurately decoding past climate from sedimentary caliber or type alone is problematic given the observed complex response of SRSs, where sedimentary layers may well be recording an event that took place 104 yr earlier. From the perspective of modern global change, it is clear that the response of SRSs to anthropogenically induced climate change has only just started. Based on the available PETM data, dramatic changes in sediment erosion, sediment flux, and sediment accumulation rates across both terrestrial and marine segments of SRSs are likely to evolve and persist for millennia to come.


A linked terrestrial to deep marine sediment routing system in northern Spain records a stratigraphic offset between the abrupt onset of warming and hydrological change at the PETM, and the onset of deposition of greater amounts of detrital grains at both the terrestrial and deep marine sections. The calculated duration of this stratigraphic offset (tlag) is on the order of 16.5 ± 7.5 k.y. for terrestrial sites and 16.5 ± 1.5 k.y. for the deep marine site, each using independent age models. This SRS-wide response and tlag value range are reproduced using a 1-D sediment transport model, which supports the mechanism of delayed sediment transport from the sediment source area to locations down-system. Our data provide new field and modeling constraints on the response of SRSs to abrupt climate change and highlight the protracted response of our landscape to current global warming, which may take millennia.


We thank two anonymous reviewers and the editor. Armitage acknowledges funding through the French Agence National de la Recherche, Accueil de Chercheurs de Haut Niveau call, grant “InterRift”. Dunkley Jones acknowledges support from Natural Environment Research Council grant NE/P013112/1.

1GSA Data Repository item 2019059, details of the 1-D sediment transport model and values used to calculate the duration of stratigraphic offset, is available online at www.geosociety.org/pubs/ft2019.htm, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY license.