Main

The Nile River is an iconic waterway that connects the humid tropics of Africa to the semi-arid Mediterranean coast and has sustained the establishment of complex societies for many millennia1. Today, the Nile discharge is controlled mainly by seasonal migrations of the African summer monsoon, bringing high amounts of rainfall to the Ethiopian Highlands2. Recent episodes of droughts and floods bear serious consequences for communities living in densely populated areas along the Nile course3,4. Global climate models participating in the Coupled Model Intercomparison Project (CMIP) forecast wetter conditions in eastern Africa as well as strongly variable and intermittent Nile river run-off5 (Extended Data Fig. 1 and Extended Data Table 1). The combination of higher discharge and more-variable Nile floods poses serious challenges for water management (storage and distribution) and social stability6. However, the model simulations of the Nile response to global warming show a wide range of possible scenarios5, leading to highly uncertain forecasts and risks assessment for potential adaptation measures.

To test model simulations and evaluate the response of the Nile River to wetter conditions, we used an annually resolved sediment archive from offshore the Nile mouth deposited during the early Holocene ‘North African Humid Period’ (NAHP)7. This period provides useful benchmarks for future climates because it is characterized by a large increase in rainfall in northeastern (NE) Africa (Fig. 1a and Extended Data Figs. 2 and 3) due to stronger insolation and a larger ocean–continent temperature gradient. The estimated amount of rainfall brought by the monsoon during the NAHP was higher by ~75% (refs. 8,9) and is thought to have generated a three- to fivefold increase in Nile River run-off10. The morphology of the Nile River itself was largely modified, with a dense network of tributaries and a much larger drainage area11,12. Geomorphologic evidence also points to the occurrence of episodes of violent rainfall and torrential floods during the NAHP, which probably led prehistoric populations to abandon the Nile valley13,14,15. So far, the frequency and recurrence of these extreme run-off events have not been examined, which has limited our ability to identify the main driving factors.

Fig. 1: Precipitation in the Nile River watershed during the NAHP and varved record of core P362/2-33.
figure 1

a, Map of humidity anomalies at 9 kyr bp (9 K) as compared with pre-industrial conditions estimated by the Alfred Wegener Institute-Earth System Model (AWI-ESM) model (Methods and Extended Data Figs. 2 and 3), with the Nile River (pink) and the location of core P362/2-33 offshore the mouth of the Nile River (star). b, Thickness record of four different sublayers corresponding to seasonal deposits (see legend) along the lithology and photograph of core P362/2-337.

Full size image

Core P362/3-33 is ideally located to record annual changes in Nile flood dynamics during the NAHP (Fig. 1a). This 6-m-long sediment core was retrieved in 2008 on the western Nile deep-sea fan and contains a pristine 5-m-long laminated sequence7,16 (Fig. 1b). The chronology was constrained using a combination of radiocarbon dating and annual layer (varve) counting. The radiocarbon ages were obtained on fossil planktic foraminifera Globigerinoides ruber (white), which live in subsurface waters and provide a reliable estimation of depositional ages7,16. Varve counting and measurements of sublayer thickness were realized under the microscope. Microfacies analysis revealed a repeating sequence of sublayers that was interpreted as a seasonal depositional cycle7, which was used to determine the boundaries of annual layers. The annual deposition of the layers in core P362/3-33 was confirmed by fitting an autoregressive gamma-walk sedimentation model based on a million Monte Carlo simulations through eleven radiocarbon ages (Extended Data Fig. 4). The excellent match between layer thickness expressed as accumulation rates and the probability distribution of radiocarbon ages confirms the presence of varves in core P362/3-33. In addition, this Bayesian model provides us with a very precise and accurate age determination with errors ranging between 90 and 250 years.

A continuous record of summer floods

We obtained a continuous record of ~1,500 varves between 9.47 and 7.94 kyr bp. Another 70 varves occur above a non-laminated (bioturbated) interval between 7.69 and 7.62 kyr bp (Fig. 2a). Microfacies analyses identified the flood layers as specific detrital sublayers related to terrestrial sediment-laden hyperpycnal deposits during the annual flood season of the Nile7,17. The flood-layer thickness varies between 0.3 and 10.0 mm, with pronounced changes in the average thickness at ~8.08, 8.40, 8.62, 8.80, 8.88, 9.17 and 9.34 kyr bp (±0.05 kyr) as indicated by change points in the sedimentation model (Fig. 2e and Methods). These changes in flood-layer thickness occur within a few decades (~30–40 years) and indicate large shifts in the volume of sediments deposited on the margin.

Fig. 2: Fluvial dynamics of the Nile River during the North African Humid Period (7500–9500 yr bp).
figure 2

a, Thickness of flood layers (grey curve), 30-year moving average (dark blue curve) and 2σ confidence (blue curves). b, Number of strong and weak floods (respectively, blue and orange curves) per 30 years (Methods and Extended Data Fig. 6). c, Median (thick dark blue line) and standard deviation (light blue line) of flood-layer thickness between change points. d, Median grain sizes (log(D50), thick gold line) between change points. e, Change points shown as probability density functions of sedimentation rate (N, number of iterations deviating from a linear accumulation rate), hashed line (median) and grey shading (2σ uncertainty). Insert: relationship between flood-layer thickness and median grain size (bi-plot), which is considered significant since the calculated P value of the null-hypothesis significance test is 0.005. The non-varved part of the record is indicated by a light orange shading and the mention ‘Not varved’.

Full size image

At our core site, the annual to centennial variability in flood-layer thicknesses tracks sediment discharge and fluvial activity in the Nile watershed rather than spatial migrations of the depocentre or sea-level changes. The western part of the Nile deep-sea fan is a stable depocentre during the Quaternary period18 and was particularly active between 15 and 6 kyr bp19. Local sea-level reconstructions obtained using an ensemble of three-dimensional mantle viscosity structures combined with the ICE-6G glaciation history model for the Nile deep-sea fan show a near-linear rise from about −18 to −2 m between 10 and 8 kyr bp20, with limited glacial isostatic adjustments (Extended Data Fig. 5a). No significant relationship could be observed between flood-layer thickness and rates of sea-level rise during the early Holocene (Extended Data Fig. 5b), which suggests that the post-glacial sea-level rise did not exert a substantial control on sediment dynamics over the western Nile deep-sea fan.

The correlation between flood-layer thickness and particle size provides a solid physical foundation for using our data to reconstruct variations in flood strength (Fig. 2d and Extended Data Fig. 5c,d). Indeed, the volume of sediment transported by floods might be influenced by complex, nonlinear interactions among rainfall, discharge and sediment availability and mobilization within large river systems. Such complex interactions might lead to heavy rainfall events causing minor erosional episodes, and vice versa21,22. By contrast, the size of particles carried by a stream is directly related to the discharge. This link is well established for the Nile River, where particulate matter during summer floods is predominantly silt-sized, while lower winter flow is characterized by finer, mud-sized particles2,23. The consistent relationship between layer thickness and particle size therefore provides robust evidence that our record of flood-layer thickness tracks past flood dynamics on interannual to decadal timescales.

Extreme floods and flickering fluvial regimes

The NAHP is characterized by the occurrence of extremely strong and variable Nile floods. In particular, the dominance of thick flood layers between 9.2 and 8.6 kyr bp (Fig. 2a) depicts a period of strong erosional activity in the Nile River basin and the deposition of large amounts of particulate matter offshore. The volume of sediments transported in this time interval was on average two to three times higher than after 8.6 kyr bp (Fig. 2c). This order of magnitude is similar to previous estimations of increases in heavy rainfall frequency and Nile run-off during the NAHP10,15,24, suggesting that rainfall intensity might exert a substantial control on the erosional activity in the Nile watershed.

Extreme variations in flood strength during the 9.2–8.6 kyr bp interval is indicated by the much wider variance of flood-layer thickness compared with intervals with lower thickness (Fig. 2c). The frequency of stronger-than-normal and weaker-than-normal floods increases in parallel, although the number of strong floods is generally higher than the number of weak floods (Fig. 2b and Extended Data Fig. 6b,c). Our results suggest that highly erosional periods were associated with increased interannual variability in fluvial dynamics and lend credence to the forecasted stronger Nile flow variance under future warmer and wetter conditions5.

Another striking observation is the occurrence of rapid changes in fluvial regime as seen by the existence of change points in the accumulation rates (Fig. 2c and Extended Data Fig. 6e). These changes typically occur within 30–70 years and are characterized by substantial modifications of the erosion regimes in the Nile watershed, leading to the deposition of flood layers with thickness of different median and variance (Fig. 2c). Such rapid switches cannot be readily attributed to a single driver but might result from threshold-like responses of the Nile River to climatic or environmental changes. For example, factors such as changes in moisture source25, vegetation cover26 or fluvial morphology (for example, the activation and abandonments of tributaries)11,15 might lead to nonlinear erosional responses. In addition, intervals of stronger volcanic forcing might have triggered Nile flow failures, especially around 8.9–8.8 and after 8.6 kyr bp27,28 (Extended Data Fig. 5e).

Drivers of flood dynamics under intensive rainfall

The occurrence of seasonal floods during the early Holocene has also been reported at onshore locations in the Nile valley, for example, at Lake Faiyum29,30. However, our continuous record of annual Nile floods during the NAHP is unique in the region and allows us to investigate the drivers of interannual hydroclimates under wetter-than-present conditions. In addition, our new NAHP palaeo-flood record is of similar length and resolution as the early Islamic ‘nilometer’ records of past Nile levels at Roda Island (near Cairo) during the Common Era (CE) (622–1922 ce)28,31,32. Therefore, the comparison of flood frequency between these two records helps determine the drivers of flood variability under contrasting hydroclimatic conditions. Once detrended and normalized, the NAHP and CE records were analysed using wavelet and multi-taper methods to detect significant oscillatory regimes (Methods and Extended Data Figs. 7–9)31,33. In both records, a strong signal is found in the interannual range with persistent periodicities between 2 and 7 years (Fig. 3). This oscillatory mode is a clear fingerprint of El Niño/Southern Oscillation (ENSO), which has been shown to drive NE African hydroclimates both at present and during the CE31,34,35. This teleconnection operates through complex ocean–atmosphere interactions that control the zonal moisture (Walker) circulation35. La Niña phases are associated with enhanced rainfall and widespread flooding in the NE African monsoonal realm, while El Niño phases are associated with droughts and negative anomalies of precipitation35. The identification of the ENSO oscillatory mode in our record of the NAHP demonstrates both the existence of ENSO during the early Holocene and the stationarity of the Pacific–NE Africa teleconnection36. This result provides an important fundament for projections of ENSO-driven climatic and environmental change in a warmer world37.

Fig. 3: Oscillatory regimes during the CE and the NAHP.
figure 3

a, Multi-taper analysis of the log-transformed nilometer data with the periodogram (black), harmonics (blue) and significance levels against a red-noise background (90%: red; 95%: light blue; 99%: green). Significant periodicities above the 99% confidence level are indicated in blue and above 95% in green. b, Wavelet analyses of the detrended and log-transformed record of high Nile flow from the nilometers31 covering the past 1,300 years, with P = 0.05 significance levels against a red noise in orange (lag = 0.9), cone of influence overlain. The ENSO band is indicated in the wavelet plots by red rectangles. c, Multi-taper analyses of the log-transformed flood thickness record of core P362/2-33 covering the interval 7.5–9.5 kyr bp. d, Wavelet analyses of the detrended and log-transformed summer flood thickness record of core P362/3-33. See Methods and Extended Data Figs. 8 and 9 for more details on time-series detrending and normalization.

Full size image

Recent studies suggested that the spatial pattern of ENSO might have been different during the early Holocene but led to persistent ENSO variability38,39. This finding supports our observations and helps to reconcile contradictory records of past ENSO variability for this time interval40,41,42,43. In addition, the spatial characteristics of the so-called ‘coastal’ ENSO dominating the early Holocene is captured in our simulations (Extended Data Fig. 2c)44. However, it should be noted that due to higher levels of atmospheric moisture during the early Holocene45, even a weak ENSO signal would be amplified and lead to more-variable floods as observed in our record and projected for future climates5,46 (Fig. 2c).

Longer, decadal to centennial-scale oscillatory modes are also identified in both CE and NAHP records (Fig. 3). Beyond the periodicities related to the solar cycles (11–22 year sunspot and Hale cycles), two robust periodicities are detected at ~54–55 and 93–100 years, which might be related to the Pacific decadal oscillation47 (50–70 years) or the Atlantic multi-decadal oscillation48 (40–60 years) (Fig. 3a,c). Although the observed periodicities cannot be attributed to a single climatic mode, their identification provides insights on low-frequency modulation of flood dynamics in the Nile basin. We show that these modes have been operating under contrasting climatic conditions of the CE and NAHP (Fig. 3a,c) and probably modulated the occurrence of flood extremes (Fig. 4). Indeed, bandpass filters applied to the log-transformed flood records show that the number of extreme flood events is driven largely by multi-decadal oscillations. Despite being widely identified in climate records, the drivers of this multi-decadal variability remain elusive and difficult to capture in global climate models49. We show here that in the Nile fluvial system, this variability might be an essential driver of extreme flood events.

Fig. 4: Multi-decadal variability drives extreme flood events during the CE and the NAHP.
figure 4

a, Number of extreme flood events per 30-year bins (blue: strongest floods; orange: weakest floods) in the high Nile level records (nilometer) during the CE, compared with the filtered log-transformed record of high Nile flood levels (pink). b, Similar comparison for the flood thickness record (Nile deep-sea fan) of the NAHP. See Methods for more information on bandpass filtering and data processing; number of extreme floods as in Fig. 2b and Extended Data Figs. 6 and 7.

Full size image

In the context of a rapidly changing climate, annually resolved records from past wetter intervals help capture the full range of natural climatic variability in flood records, which is crucial to build reliable forecasting tools50. Regular flooding in the Nile River is essential for the subsistence of dense populations in NE Africa but is also the source of important environmental and political pressure (flooding, trans-border conflicts)6. We provide here sedimentary evidence of a highly variable river with extremes in flood intensity and strong erosive events during the NAHP. The stationarity of high- and low-frequency signals under contrasting climates of the CE and NAHP provides a strong basis for testing hydrological models incorporating this climatic variability. However, the existence of nonlinearities in the river response to multi-centennial oscillations leading to rapid switches in erosion regime calls for further research. Finally, our data show that the volume of sediments transported and deposited offshore the Nile River varies by a factor 2 to 3 throughout the record. Given that sediment transport is, on interannual to decadal timescales, directly related to river flow, our observations provide relevant benchmarks for scaling future infrastructure.

Methods

Microfacies description in core P362/2-33

Core P362/2-33 contains ~5 m of finely laminated sediments that consist of alternating dark- and light-coloured millimetre-thick layers. The microfacies and their chemical composition have been described in detail7, and we present here the whole-core thickness measurements for each sublayer type (Fig. 1b). A regular sequence of four types of sublayers occurs throughout the core, with two of them being always present. These sublayers have been associated with seasonal depositional regimes: (1) summer floods, (2) autumn blooms (plankton), (3) winter run-off and (4) authigenic carbonates. Summer flood sublayers are characterized by coarser detrital grains and high Ti/K ratios while winter run-off consists of clay-sized, low Ti/K sediments7. Autumn blooms occur after the summer flood and are characterized in core P362/2-33 by the presence of planktic foraminifera shells. Such blooms were observed during historical Nile floods51. Authigenic carbonates are fine-grained calcite deposits that are interpreted to form in the bottom waters due to anoxic to sulfidic conditions at the end of the spring7. Several event layers were also identified in the core as matrix-supported layers and associated with large-scale remobilization episodes or large flood events17. Their facies is clearly distinct from that of the summer flood sublayers.

The summer flood and winter run-off deposits are the main contributors to the total layer thickness (Fig. 1b). Large changes in the layer thickness of all sublayers are observed throughout the record, with thicker sublayers between ~9.2 and 8.7 kyr bp.

Age modelling

To estimate the age–depth relationship in core P362/2-33, we used 11 radiocarbon ages measured on planktic foraminifera Globigerinoides ruber that were published in ref. 16. The radiocarbon ages were calibrated using the Marine20 calibration curve, which already contains corrections for reservoir ages52.

We constructed the age model using a Bayesian age–depth modelling approach, which leverages relative age information from varve counting as prior information to model absolute age information from radiocarbon dating. The basic aim of this approach is to integrate the two sources of information and obtain a more accurate and precise age–depth relationship. Specifically, we used the varve counts to inform the prior distribution of ages at each depth interval and then updated this distribution using the radiocarbon dating data. By combining these two sources of information, we were able to produce a robust and reliable age–depth model for core P362/2-33.

To incorporate varve counts into a parameterized model for accumulation rates, we employed autoregressive gamma walks, as proposed by ref. 53. The general form of autoregressive gamma walks is given by

$${x}_{j}=\omega {x}_{j+1}+\left(1-\omega \right)\gamma$$
(1)

where xj is the accumulation rate (in yr cm−1) in depth interval j (with constant spacing Δc), ω is the autocorrelation at lag 1 and γ is a gamma distributed random variable with mean a and shape b (\(\gamma \sim {\mathrm{Gamma}}\left({a}_{\gamma },\,{b}_{\gamma }\right){{\it{iid}}}\)). To take the pronounced changes in accumulation rates of core P362/2-33 into account, we subdivided the sequence of varve counts into sections with similar mean accumulation rates using the Pruned Exact Linear Time change-point detection algorithm by ref. 54. Subsequently, the three parameters ωd, ad and bd were estimated for each section d. The resulting model was realized 106 times and each realization randomly anchored to the probability distribution of the youngest radiocarbon age (Extended Data Fig. 3, light blue area). Subsequently, we updated the ensemble of model realization by selecting age–depth relations maximizing the correspondence to radiocarbon ages (dark blue areas in Extended Data Fig. 4). This optimization process allows us to identify the most probable varve stratigraphy given the absolute age information of radiocarbon dates and quantify the uncertainty in the age estimates.

The modelled age–depth relationship falls within the uncertainty of all the 11 radiocarbon dates (Extended Data Fig. 4). This allows us to build an accurate age model with a high precision (uncertainties are below 100 years) and to maintain the high precision level throughout the profile. The excellent match between accumulation rates derived from varve counting and radiocarbon measurements also confirms the hypothesis that the sequence of four sublayers represents an annual cycle and therefore real varves.

Data analysis

To analyse the data, we used well-established statistical parameters. Following ref. 43, we have calculated the amount of extreme events (high or low flood) per 30 years on both the flood-layer thickness record of core P362/2-33 (Extended Data Fig. 6c) and the nilometer (high Nile) levels (Extended Data Fig. 7c)31. However, because our data of flood thickness are constrained (cannot be negative), we have used the log-transformed records and applied a 9 year high-pass filter using the FIR filters in the Past4 software55 (Extended Data Figs. 6b and 7b). Extreme events are defined as values higher (lower) than the 95 (5) percentile (Extended Data Figs. 2b and 3b). We then calculated the moving number of extreme high and low floods per 30 year window (Extended Data Figs. 6c and 7c).

Extended Data Fig. 6e shows the distribution of probable ages for the estimated change points (‘Age modelling’) in the final age model. The histograms represent the number of change points detected by a million Monte Carlo models, which allows us to provide an estimation of uncertainty around each change point and a determination of the skewness of the change-point distribution. The raw flood-layer thicknesses were then split in seven parts between the change points, and the changes in median and standard deviation were calculated (Extended Data Fig. 6d). These analyses show a doubling of the median thickness between different parts of the record, with thicker layers associated with a higher standard deviation than thinner layers.

Glacial isostatic adjustment modelling

We used outputs from the three-dimensional glacial isostatic adjustment model VILMA (VIscoelastic Lithosphere and Mantle Model) to evaluate the changes in relative sea level (RSL) at the core site (29° 50′ E, 31° 36′ N)20. This model ensemble allows us to test the effect of upper-mantle and asthenosphere viscosity structure on post-glacial sea-level rise according to the ICE-6G glaciation history56. The computed RSL variations show a near-linear rise from −18 to −2 m with very little effect of mantle viscosity structures (Extended Data Fig. 5a). After averaging the 18 model outputs (ensemble mean), we calculated the first derivative, which tracks the changes in rates of RSL changes. Finally, we calculated the median values of log-transformed RSL rates between change points to compare them with mean changes in flood-layer thickness (Extended Data Fig. 5b).

Grain-size measurements

The distribution of siliciclastic grains was measured at the University of Innsbruck on a set of 80 samples. After being rinsed three times with ultrapure water (MilliPore system), sediments were decarbonated using buffered acetic acid and the organic matter was removed using concentrated hydrogen peroxide. All solutions were prepared fresh in the clean lab at GFZ Potsdam using ultrapure acids. Samples containing the siliciclastic fraction were then shipped to Innsbruck, where they were mixed with sodium pyrophosphate to avoid the formation of clay aggregates and measured in triplicate on a Malvern Mastersizer 3000.

The results show grain sizes ranging between 0.1 and 100.0 µm with main modes at ~0.5–0.6 and 3–4 µm (not shown), very similar to results obtained on this sediment core previously at low resolution16. For the present study, we report the median grain size (D50), which is modulated by hydrodynamic sorting and flow strength23. The data are presented as individual samples (Extended Data Fig. 5d) and as log-transformed averages (median) computed between change points (Extended Data Fig. 5c). Comparing data located in different data spaces (only positive, unconstrained) in a log–log space allows us to compute linear regressions and estimate their statistical significance.

Time-series analysis

To investigate the frequency domain of our annual record of past Nile flow, we applied time-series analyses on the nilometer records of the CE, which are also annually resolved31. We present here the multi-taper method results performed on the nilometer and Nile deep-sea fan records using the online SSA-MTM toolkit version 4.4 (ref. 57). The MTM analysis was performed on normalized, log-transformed (log; Extended Data Fig. 8) and differentiated log-transformed (log[xt/xt−1]) (log(diff); Extended Data Fig. 9) flood-layer thickness and high Nile level data with a resolution of 3 and 5 tapers. The SSA-MTM toolkit allows us to detect harmonics (in blue) and frequencies significant at the 90%, 95% and 99% confidence levels against a red-noise background. These analyses show that both records contain significant frequencies in the range 2–7 years, which is typical of ENSO. The results for the high frequencies (interannual range) are similar for the log and log(diff) data, but the log(diff) spectra are distinctly ‘blue’ (Extended Data Fig. 9a,c) compared with the ‘red’ spectra for log data (Extended Data Fig. 8a,c). This is because the log(diff) procedure filters lower frequencies out of the signal. We therefore favour the log-transformed data for interpreting the MTM spectra (Extended Data Fig. 8a,c).

Another useful method to detect ENSO fingerprints in time series is to apply wavelet analyses33. We used the software PAST to perform a continuous wavelet transform55, using a Morlet wavelet of wavenumber 6, on the log (Extended Data Fig. 8b,d) and log(diff) data (Extended Data Fig. 9b,d). Strong signals are detected in the frequency bands between 2 and 7 years in all wavelet analyses. The signal power corresponding to the correlation strength of the mother wavelet is indicated by the colour coding (yellow: high; blue: low)58. Signal powers below the P = 0.05 significance level above the null hypothesis of a red noise (lag = 0.9) are highlighted in red, and cones of influence (delimiting the regions with boundary effects) are superimposed in the lower corners of the plots. For further interpretations, we favour the log(diff) records (Extended Data Fig. 9b,d), which allow us to explore variability in the interannual range on unconstrained records where the autocorrelation effects were removed.

Numerical modelling

We used numerical simulations from several sources to estimate the changes in rainfall over NE Africa for the next century and during the NAHP. Simulations for end-of-the-century precipitation distribution (modelled ensemble mean precipitation) based on the sixth phase of the Coupled Model Intercomparison Model (CMPI6)59 have been analysed to elucidate the spatial distribution of precipitation for different warming scenarios. We compared simulations for the shared socioeconomic pathways (SSP) 1.0–2.6 W m−2 (CMIP6-SSP126), 2.0–4.5 W m−2 (CMIP6-SSP245), 3.0–7.0 W m−2 (CMIP6-SSP370) and 5.0–8.5 W m−2 (CMIP6-SSP585) with the 1850–2010 baseline (CMIP6-Historical) (Extended Data Fig. 1 and Extended Data Table 1). These four scenarios vary between the low-end emission (‘taking the green road’, SSP126) and the high-end emission (‘taking the highway’, SSP585) scenarios60,61. The scenario emissions are anticipated to produce a radiative forcing in 2100 of approximately 2.6 W m−2 in SSP126 due to an increasing shift towards sustainable practices, and 8.5 W m−2 in SSP585 due to a fossil-fuel-driven development. The mean seasonal precipitation (August–September–October) has been computed for the ensemble mean following four warming scenarios according to the SSP (Extended Data Fig. 1).

To determine whether the rainfall patterns during the NAHP can be used as analogues of predicted monsoon changes in warmer and wetter climates, we also used simulations performed using the AWI-ESM for climatologies at 9 kyr bp (9 K) to evaluate the rainfall dynamics over eastern Africa during the NAHP (Extended Data Fig. 2). These simulations were performed using the boundary conditions as defined in ref. 62, with orbital parameters, greenhouse gases and ice sheets set at 9 K and an integration time of 800 years. The simulations depict an increase in summer rainfall in the monsoon region over North Africa (Extended Data Fig. 2a,b) and a Pacific Ocean configuration similar to a coastal El Niño39,44. This untypical El Niño is characterized by a coastal warming in the eastern Pacific in otherwise basin-wide neutral or La Niña-type conditions and leads to sizeable rainfall in the Maritime Continent44 (Extended Data Fig. 2c).

To test the reliability of the modelling experiment for estimating the rainfall dynamics at the Nile sources, we computed monthly means of daily precipitation in control runs for both the AWI-ESM and the CMIP6 ensemble and compared these with the observed annual rainfall cycle obtained from the Climatic Research Unit (CRU TS4) (Extended Data Fig. 3)63. Both the AWI-ESM and the CMIP6 ensemble means capture well the seasonal rainfall dynamics but tend to underestimate rainfall in spring and autumn in the AWI-ESM and in summer in the CMIP6 ensemble mean (Extended Data Fig. 3a,b). The summer precipitation is largely enhanced in the 9 K simulation, whereas the CMIP6 ensemble mean shows a stronger rainfall anomaly at the end of the wet season (August–September) (Extended Data Fig. 3c–e). It is remarkable that the summer rainfall anomaly is very similar for the 9 K and for the warmest SSP585 scenario (Extended Data Fig. 3e). The cumulative annual rainfall falls in similar ranges for the simulations at 9 K and in all SSP scenarios (Extended Data Fig. 3d).