Human and climate global-scale imprint on sediment transfer during the Holocene.

Accelerated soil erosion has become a pervasive feature on landscapes around the world and is recognized to have substantial implications for land productivity, downstream water quality, and biogeochemical cycles. However, the scarcity of global syntheses that consider long-term processes has limited our understanding of the timing, the amplitude, and the extent of soil erosion over millennial time scales. As such, we lack the ability to make predictions about the responses of soil erosion to long-term climate and land cover changes. Here, we reconstruct sedimentation rates for 632 lakes based on chronologies constrained by 3,980 calibrated 14C ages to assess the relative changes in lake-watershed erosion rates over the last 12,000 y. Estimated soil erosion dynamics were then complemented with land cover reconstructions inferred from 43,669 pollen samples and with climate time series from the Max Planck Institute Earth System Model. Our results show that a significant portion of the Earth surface shifted to human-driven soil erosion rate already 4,000 y ago. In particular, inferred soil erosion rates increased in 35% of the watersheds, and most of these sites showed a decrease in the proportion of arboreal pollen, which would be expected with land clearance. Further analysis revealed that land cover change was the main driver of inferred soil erosion in 70% of all studied watersheds. This study suggests that soil erosion has been altering terrestrial and aquatic ecosystems for millennia, leading to carbon (C) losses that could have ultimately induced feedbacks on the climate system.

Changes in human land use have increased soil erosion rates globally (1, 2) and have led to substantial alterations in soil stability and productivity, functioning of river ecosystems, and lateral transfers of sediments and carbon, affecting ultimately long-term change of climate (3–5). The impacts of human land uses and climate on soil erosion amplitude and spatial extent are closely tied to the timescale considered. On geological timescales of millions of years, erosion is controlled by climate and tectonic cycles (6), which, in turn, regulates atmospheric CO2 and climate through chemical weathering (7) and organic C burial in sediments (8). During the course of the Holocene, anthropogenic modifications of catchments, including vegetation clearance and burning, as well as agricultural and urban expansion, have led to rapid fluctuations in local to regional soil erosion rates (9–11), and ultimately this has accelerated erosion 10- to 100-fold in some regions (1). However, it is still unclear to what extent climate and human factors have interacted and controlled long-term (i.e., >100 y) soil erosion rates, and when a significant portion of the Earth surface shifted to human-driven soil erosion rates (3).

Natural archives preserved as lake sediments provide a key source of evidence for assessing soil erosion that occurs in lake catchments and is integrative of all fluxes and processes that remove soil, rock, or dissolved material from the watershed, including gully, till, or rill erosion. For many lakes, sediment accumulation rates (SARs) from cores collected in profundal depositional zones can provide an estimation of gross changes in the rate of deposition of allochthonous sediment related to erosion (12), although SARs variability can also be related to other processes (). Due to their relatively high temporal resolution and often continuous history of deposition, lake cores are valuable archives of catchment-scale soil erosion dynamics, (i.e., temporal variability in catchment erosion intensity). Furthermore, the abundance of Holocene lake records and the diversity of catchment systems around the globe have established lake sediment cores as strategic data repositories for a better understanding of transfer of matter on the surface of Earth.

The objectives of the present study are to identify first increases in sediment supplies to lakes to infer soil erosion temporal dynamics across globally distributed sites to quantify spatiotemporal patterns and to investigate the effects of climate and land cover change (LCC) during the Holocene using SAR records from lakes. Syntheses such as this one have the potential to make key advances in our understanding of critical events for Earth system processes. A major limitation for identifying the transition from an overarchingly climatically controlled system to one in which human drivers explain changes in soil erosion rates is the lack of spatio-temporally resolved empirical data, which are needed to discriminate between local stressors mediated primarily by human activities and climatic factors (13, 14). Furthermore, our gap in knowledge about changes in global soil erosion rates limits our ability to clearly define of the onset of the Anthropocene (15). Finally, the lack of synthesis also limits our understanding of the timing and extent of increased soil erosion and hinders predictions about long-term total carbon losses from soils and land-atmosphere exchanges (3).

In this study, we bring together a global synthesis of pollen samples (n = 43,669) and 14C dates (n = 3,980) to estimate temporal trends in land cover and lake SARs in 632 lake-watersheds (Fig. 1 and Dataset S1). The ages of the first increase in soil erosion in lake-watersheds were then derived from the SAR trends (). The lake pollen records (retrieved from three pollen databases [European (EPD), Global (GPD), and African Pollen Database (APD)]) were grouped into 5 plant-functional types and 2 land-cover types (arboreal and nonarboreal). Our a priori prediction is that decreases in the proportion of arboreal pollen (AP) would reflect land cover changes (i.e., land clearance) that, in turn, would lead to soil degradation and erosion. Trends in long-term precipitation, air temperature, and wind were assessed using MPI-ESM-1.2LR transient Holocene simulations (Dataset S1). Generalized additive models (GAMs) were used to estimate the long-term contribution of land cover and climate on erosion. Our approach is summarized in .

Fig. 1.

(A) Location of the 632 study sites and changes in lake sediment rates over the last 10,000 y. Increases in lake sediment rates (red dots) were observed in 35.1% of the watersheds according to M-K τ coefficients; the rest of the sites show a decrease (green dots) or steady (white dots) rates. (B) Distribution of M-K trend τ coefficients. (C) Number of dated cores containing a pollen record and distribution of 14C dates (intervals = 12,000 y) used to calculate SARs.

Sampled lake-catchments capture a comparable range of morphometric properties (i.e., site elevation, watershed area, lake area, and lake depth) relative to the 14.3 million lakes in the global HYDROLAKE database (16), as well as a wide range of soil types, land cover, and climatic conditions ( and Datasets S2–S5). Based on our analyses and reconstruction of SARs ( and ), we found that ∼35% of the long-term records of SARs for 632 sites were essentially increasing (Kendall’s tau coefficient τ > 0.2, P value <0.05), while 26% were decreasing (τ < −0.2, P value <0.05) and 39% were stationary or showed no significant trends (Fig. 1 and Table 1). Piecewise regression modeling on each time series () showed that the global signal of SAR and mass accumulation rates (i.e., MAR, corrected from compaction, see ) increased since 4,000 calendar years before the present (cal. BP) (i.e., mean breakpoint year = 4,150 cal. BP [minimum-maximum range 3,250–4,547]). The fraction of lakes recording positive anomalies in SARs also increase circa 4,000 cal. BP, at rate of 12 ± 0.2% ka−1 (P < 0.01; Fig. 2 ).

Table 1.

Results of M-K tests and multiple regression models

	All sites (n = 632)	GAM results
		GAM sites		P value LCC < P value Air T°		P value LCC < P value Prec.		P value LCC < P value climate
τ (M-K test)	%	n	%	n	%	n	%	n	%
>0.2 (Ps < 0.05)	35.1	209	57.3	141	67	144	69	122	58
>0.5 (Ps < 0.05)	20	116	31.8	86	74	86	86	78	67
<0.2 (Ps < 0.05)	23.8	155	42.5	115	74	111	111	102	66
<0.5 (Ps < 0.05)	12.8	79	21.6	64	81	67	67	63	80
−1 to 1	100	365	100	256	70	255	70	224	61

The M-K trend coefficients (τ) show that 35.1% of all our study sites have recorded an increase in SARs over the last 10,000 y. GAM showed that the probability of temporal variation in SARs of lakes depends on changes in the proportion of AP over the past 12,000 y more than changes in climate variables in 58–80% cases. GAM sites correspond to the number of sites used for the GAM models. GAM formula to assess controls on SARs is a function of changes in land cover, air temperature, and precipitation: SAR ∼ s(AP) + s(Prec) + s(T°).

Fig. 2.

Trends in lake SARs (proxy of erosion), land cover as well as climate change during the Holocene. (A and B) Density plots of the anomalies in SARs and AP percentages (43,669 pollen samples), respectively, for 632 14C-dated sites (3,980 calibrated 14C ages). Global trends for SARs are shown in 632 lakes based on an GAM, with 95% confidence intervals on the predicted means (C), as SARs medians trends calculated every 50 y (n = 632 sites) and point breaks detected on the SAR medians curve global trend based on a bootstrap distribution (D). (E) Changes in land cover are shown as trends AP percentage for the subset of sites (n = 116) recording an increase in SARs over time (solid green line), and for a subset of sites (n = 287) recording no significant changes in SARs over time, denoted here as benchmark sites (dashed green line). (F) From ref. 15: trends in global temperature anomalies (relative to the average over the period 1961–1990 in blue) and atmospheric CO2 (red) from the GRIP ice core59.

A common observation across our sites recording accelerated SARs during the last 4,000 y was the trend of deforestation (10), which coincided with the increase of SARs (Fig. 2 ). Our statistical analyses support the conclusion that land cover change was the leading driver for the onset of accelerated lake sediment deposition. Assessment of the relative contribution of variables in a linear regression model show that variations in SARs were explained best by land cover, here expressed as arboreal pollen (AP) percentage (). Using general additive models on each lake-watershed time series, we found that the probability of positive anomalies in SARs increased as the anomalies in the proportion of AP decreased (P < 0.0001) in 70% of our studied watersheds (Fig. 3, , and Table 1). In contrast, SAR time series were less related to the changes in anomalies of precipitation and temperatures (P < 0.01) (Fig. 3 , Table 1, and ) (mean variances explained around 0.6). Sites displaying relatively constant sedimentation rates through time (i.e., benchmark sites, Mann–Kendall (M-K) test ε[−0.2;0.2]) had no clear signs of deforestation in our records (Fig. 2).

Fig. 3.

Holocene-scale trends in 12 annually resolved varved sites. (A) Density plot normalized SARs trends for 12 varved lakes. (B) Negative correlation between SARs and AP fraction of the varved sites. (C) Deviance explained in the SAR trends by changes in land cover, precipitation, and temperatures in the GAM models run for each of the 632 14C-dated sites. (D) GAMs were run 632 times, and P value results inform the contribution of LCC, precipitation, and air temperature to SARs.

Trend analyses show that despite the heterogeneity in the signal of SARs among sites, some regional patterns were also observed, highlighting different socio-economical histories. For instance, the rise in SARs is later in North America compared to Europe () and likely corresponds to the widespread establishment of European agricultural methods shortly after colonization (14). The decrease in sediment fluxes recorded in 23% of the lakes was not explained by land cover and climate variables, and instead may be a result of increased sediment retention due to water uses and river management. Human modification of rivers about 3,000 y ago increased with embankments, dams, watermills, especially in the Chinese and Roman Empires, which could have increased sediment retention (17).

The long-term reconstructions of sediment flux dynamics—i.e., relative changes in sedimentation rates—at global and regional levels presented herein involve some uncertainties due to the nature of 14C dating and the relatively poor dating resolution at some sites, but nonetheless provide a comprehensive global assessment of the timing of the “low frequency” variations (Fig. 2 ). To investigate how short-term variations in land cover could affect SARs, we conducted further analyses on varved records from 12 lakes around the world that had annually dated SARs and pollen records extending over the Holocene (Fig. 3 and ). Consistent with our larger but lower temporal resolution dataset, we found that the SAR trends were anticorrelated with the percentage of AP (Fig. 3). High rates of change in SAR and pollen records were generally evident in the Early Holocene, followed by an increase of AP percentage associated with decreasing SARs, and then finally during the Late Holocene, we observed increases in SARs and a decrease in AP percentage ().

Interestingly, both the varved and lower temporal resolution 14C-dated cores highlight the importance of climate driven change in SARs for the Early Holocene. This Early Holocene change (i.e., between 12,000 and 8,000 cal. BP) can be explained by the retreat of North Hemisphere ice sheet that progressively gave rise to the creation of most contemporaneous lakes and new lake SAR records. Long records of SARs in lakes found in previously glaciated parts of Europe and North America tended to show high pre-Holocene rates declining to minimum values during the Early to Middle-Holocene (14). For example, AP fraction shows an increase and SARs showed a decrease from 12,000 cal. BP to 10,000 cal. BP in Holzmaar (Germany), or between 10,000 cal. BP and 8,000 cal. BP in the Black Sea as a result of the transient response of vegetation to climate changes (). During the early to mid-Holocene, the lake-watershed systems seem to progressively reach a stabilization phase associated with the vegetation development, soil maturation, and lake ecosystem stabilization, as evidenced by relatively constant SARs between 8,000 and 4,000 cal. BP. Since the middle of the Bronze Age (3,500 cal. BP), the growing human influence upon erosion is the emergent global trend (9) (Fig. 3 and ).

Collectively, these findings suggest that the abundance of trees in watersheds was the leading factor explaining temporal variations in soil erosion, with anthropogenic deforestation explaining accelerated soil erosion over recent millennia. Our results highlight the importance of large-scale (i.e., in terms of distribution of the lake records globally, but not in terms of the total contributing area) and long-term processes on soil erosion inferred from SARs, and how human activities began to act on these processes much earlier than other signatures of the Anthropocene globally, e.g., human appropriation of the N cycle since 1860 (18). This global-scale analysis of paleolimnological records adds to the growing evidence that humans are simultaneously increasing the river transport of sediment through soil erosion and decreasing this flux to the coastal zone through sediment retention in reservoirs (19). Altogether, our study shows strong evidence that sedimentation in lakes recorded a widespread signal of early increases that are likely the results of land use clearance and increased soil erosion.

Methods

Reconstructing Soil Erosion Dynamic.

Paleolimnological records were used to reconstruct a signal of soil erosion at both the annual and centennial scale for the past 12,000 y. It is now well established that long-term erosion can be reconstructed from geochemical signatures and SARs in well-dated lake sediments archives (10, 13). For many lakes, SARs from central profundal position in a lake can provide an estimation of gross changes in the rate of allochthonous inputs related to erosion, especially when SAR estimates are converted to influx rates of allochthonous inorganic material (12). As lacustrine productivity also contributes to lake sediments, SAR data were only used to infer, rather than to measure, the rate and processes of erosion (20). Our study of SARs was furthermore constrained to before the widespread cultural eutrophication has occurred since the Industrial Revolution. Under such baseline conditions, the SARs can conceivably serve to inform the first increase in soil erosion or sediment retention on millennial timescales (21).

Radiocarbon (14C) dating predominates for long Holocene timescales and provides the temporal frameworks used in this paper. Radiocarbon (14C) dating of lake sediments can document the rates at which sediment accumulated over the past 50,000 y. With a radioactive half-life of 5,730 y, the radioactive decay of 14C is minimal for samples less than 300 y of age and needs to be completed with other dating techniques. The preservation of varves (i.e., annually laminated sediment) in some hypoxic lakes offer the considerable advantage to precisely date sediments and reconstruct annual SAR variations by measuring varves from the sediment/water interface down-core, potentially without interruption (22, 23). Sediment archives also preserve indications of past land cover transformations via the pollen records, allowing the study of effects of LCCs on soil erosion when pollen records are available.

Palynological and radiometric data of lake sediment records were collected in May 2017 from three open source databases: the GPD, the EPD, and the APD. Sites included in our analyses had to satisfy several conditions. From the initial 4,711 sites inventoried, 632 sites were accepted because the study site: (i) was a lake sensu stricto (documented by authors, and confirmed by satellite observations) and a not peatbog; (ii) had to contain sediment core profile(s) that were dated with radiometric methods (i.e., minimum of three dated intervals; median here = 6.8 ± 4.4 SD) and/or by varves counting; (iii) had to have temporally resolved pollen records (i.e., minimum of 10 intervals over the sediment core; median = 40 ± 45 SD).

14C Age/Depth Models.

A total of 4,581 radiocarbon ages were calibrated to calendar (±2 σ) years before present using the CLAM 2.2 package in R software (24) and the northern and southern hemispheric calibration curves IntCal13 and SHCal13, respectively. Negative radiocarbon ages were calibrated with the postbomb curves from Hua et al. (25). For each core, 14C dates were processed semiautomatically to obtain 632 age-depth models after 1,000 iterations for each profile. Hiatuses and events of instantaneous deposition were removed from the analysis when the original authors recorded them. Linear depth-age models were selected to derive SARs in order to reduce temporal bias in estimating SAR changes by smoothing methods. Accumulation rates were derived from depth-age models at 10- and 50-y resolution, SAR = (di − di+1)/(ti+1 − ti), where t are calibrated ages and p sediment depths (). The SAR signals were normalized for each profile: SARnorm = (SAR − SARmean)/(SARmean), where SARmean corresponds to the mean of SARs (related from 8,000 BP to top of the lake sediments) for each related site location. To run the GAM analyses, only intervals including pollen records have been kept (median of 28 per sites) in order to not overestimate the confidence in the temporal resolution.

The time-dependent compaction of the sediment was taken into consideration using the equation of Emerson et al. (26): MAR = SAR × d × (1 − porosity). Where MAR = mass accumulation rate (g·cm−2·y−1); SAR = Sediment accumulation rate (cm·y−1); d = sediment density of 2.5 g.cm−3 (fixed; we assumed that changes in sediment density were relatively low compared to changes in SARs during the Holocene); porosity = water content (%). The porosity factor has been adjusted against 15 water content profiles collected from literature search and GAM regression ().

Holocene Land Cover and Climate Data.

The 632 sediments records from the three pollen databases contained data on pollen and other microfossils and macrofossils, which enables the assessment of land cover change over 3–12 millennia (27). Altogether, 5 plant structural functional types (PFTs) associated to terrestrial groups were used to analyze the temporal changes in land cover during the last 12,000 y: Herbs (HERB), dwarf shrubs (DWAR), bryophytes (BRYO), vascular cryptogams (Pteridophytes, VACR) and trees and shrubs (TRSH). AP fractions were then calculated and normalized to the total terrestrial PFT according to standard procedure (27) in order to assess changes in soil cover protection: Arboreal Pollen (AP %) = [TRSH/(TRSH + HERBS + VACR + DWAR + BRYO)] × 100.

Climate can be reconstructed from pollen records, but those data were already used in this study for LCC reconstructions. To avoid circularity and standardize our approach, climate time series for the past 8,000 y (i.e., average, minimum and maximum of precipitation, air temperature, and wind) were collected for all sites from the transient Holocene simulation with MPI-ESM-1.2 LR, with land and atmosphere spatial resolution of 1.9° × 1.9°, developed at the Max Planck Institute for Meteorology references (28). Time series were resampled at 50 and 10 y using piecewise cubic Hermitean interpolation polynomials (Pchip-R) for GAM and M-K analysis, respectively (). Finally, modern site characteristics data were extracted from the HydroLAKES database and from the modeled area of the hydrological basins calculated using the flow accumulation and flow direction rasters made available from HydroSHEDS (16).

Sensitivity Analysis on Annual Records.

SARs have been described as a reliable proxy of soil erosion, but the contribution of lake productivity to temporal variations of SARs warrants closer investigation. We first addressed this gap by considering a subsample of our varved lakes where multiple proxies were measured and can be used to estimate the contribution of autochthonous and allochthonous deposition to the SAR variability. Specifically, we quantified the proportion of variation in the annual SAR measurements for each varved record that could be explained by the terrigenous vs. lacustrine productivity signals from the same record (). The titanium (Ti) and the magnetic susceptibility were used as proxies for erosion, whereas biogenic silica flux, silicium-titanium (Si:Ti) ratio, and calcium-magnesium (Ca:Mg) ratio were used as proxies for productivity (12, 21). Second, we analyzed the response of our 632 SAR records to fluctuations in air temperature, precipitation, cropland cover, and human population density; here, we used the latter two metrics as an indirect proxy of nutrient supplies to the lake, assuming that these measures were the main driver of lake productivity during the Holocene. Nutrient supplies were estimated indirectly using the HYDE data on cropland and human population density (29).

Numerical Analysis.

To disentangle the impacts of climate and land cover drivers on erosion dynamics and to apportion changes to specific drivers, we use the GAMs techniques that are also able to take account of uncertainties associated with the available data. Variance Inflation Factors (VIFs) were used to detect collinearity among predictors (cutoff VIF > 0.4). We ran 632 models using the lake SAR time series as the response variable, the normalized AP percentage, precipitation, and air temperatures as fixed effect explanatory variables, and computed P values (reported on Fig. 1) and the deviance explained () related to contributions of LCC and climate variables. Confidence intervals were derived using the SEs produced by the predict.gam function in R, with type = ‘response’ specified in the model, mgcv library (30). Nonparametric M-K tests for monotonic trends were used to quantify trends of land cover for each of the 1,607 watershed time series within the past 12,000 y. This analysis was based on the Kendall rank correlation coefficient and was conducted using the Kendall library (31). A positive score shows a monotonically increasing trend, whereas a negative value shows a monotonically decreasing trend.