North African humid periods over the past 800,000 years

The Sahara region has experienced periodic wet periods over the Quaternary and beyond. These North African Humid Periods (NAHPs) are astronomically paced by precession which controls the intensity of the African monsoon system. However, most climate models cannot reconcile the magnitude of these events and so the driving mechanisms remain poorly constrained. Here, we utilise a recently developed version of the HadCM3B coupled climate model that simulates 20 NAHPs over the past 800 kyr which have good agreement with NAHPs identified in proxy data. Our results show that precession determines NAHP pacing, but we identify that their amplitude is strongly linked to eccentricity via its control over ice sheet extent. During glacial periods, enhanced ice-albedo driven cooling suppresses NAHP amplitude at precession minima, when humid conditions would otherwise be expected. This highlights the importance of both precession and eccentricity, and the role of high latitude processes in determining the timing and amplitude of the NAHPs. This may have implications for the out of Africa dispersal of plants and animals throughout the Quaternary.

Supplementary Figure 3. Composite seasonal precipitation during periods of precession maxima (PMax), the 20 North African Humid Periods (NAHPs), and the difference between the two (NAHP-PMax).a) Summer (June, July, August; JJA) precipitation (mm/JJA), b) winter (December, January, February; DJF) precipitation (mm/DJF).Only the differences that are considered 95% confident according to a student t-test are shown.showing only the most abundant plant function type (PFT) in each grid square, c) vegetation reconstruction utilising the North African vegetation classification of Larrasoaña, et al. 3 which demonstrates regional vegetation types.Only the differences that are considered 95% confident according to a student t-test are shown.Grassland (<10% woody cover) Supplementary Figure 7. Timeseries of modelled Saharan (15:30°N, 15°W:35°E) annual precipitation (mm/yr) and the Atlantic Meridional Overturning Circulation (AMOC) index.The AMOC index is defined as the mean annual strength of the AMOC between 40°N:50°N at 800m water depth.Supplementary Figure 8. Composite latitude/height plots for the Western Sahara (-15:30°N, 10°W:20°E) during the skipped beats, for the NAHP_All experiment, the Orb_GHG experiment, and the difference between the two (NAHP_All-Orb_GHG).The NAHP_All experiment includes all boundary conditions and the Orb_GHG experiment incorporates orbital and greenhouse gas boundary conditions.The difference panels therefore represent the impact of the ice sheet on the climatology.a) Summer (JJA) lagrangian tendency of air pressure (filled colour contours; pa/s) where negative (positive) values represent upward (downward) motion, and zonal wind (black line contours; m/s) where solid/positive (dashed/negative) contours represent westward (eastward) air flow.b) Summer (JJA) temperature (filled colour contours; °C), geopotential height (black line contours; m) and relative humidity (green line contours; %).Solid (dashed) contours represent positive (negative) differences.Only the differences that are considered 95% confident according to a student t-test are shown.

2) Model Validation
A broad validation of the HadCM3B model against a range of observational datasets is provided in section 5 of Valdes, et al. 5 .Their analysis shows that the model accurately reproduces key aspects of the climate system, including surface air temperatures (SATs), sea surface temperatures (SSTs), ocean circulation patters and precipitation.Valdes, et al. 5 also shows that HadCM3B outperforms many of the CMIP5 models particularly with respect to surface air temperatures (see their Figure 2).
The snapshot modelling approach has previously been shown to accurately simulate Quaternary palaeoclimate 6,7 , however there are a number of limitations.Firstly, the approach assumes that the climate system is in equilibrium at each time period and that these simulations are a reasonable approximation of the climate at each time period.This assumption has been justified in previous studies [7][8][9] .Alternative methodologies, such as transient simulations with accelerated forcing, face similar assumptions.The model biases that might affect our simulated climatology are outlined in Valdes, et al. 5 , and are comparable to other models.The methodology does not capture millennial scale climate variability, specifically the Dansgaard-Oeschger and Heinrich events which are evident in North Atlantic climate over the Quaternary 10 .In Armstrong, et al. 6 , millennial scale variability was added in post processing, by scaling the results of a North Atlantic freshwater hosing experiment with the Greenland Ice core temperature record in order to determine the amplitude and spatial pattern of the DO/HE events.However, the uncertainty associated with this methodology increases with locations further away from Greenland as it utilises the Greenland Ice core record, and so we do not include that methodology here.The model also does not incorporate atmospheric aerosol/dust feedbacks.Despite the uncertainties, our approach is highly efficient as the model simulations can be run simultaneously.A fully continuous simulation would have taken many years to complete on a high-performance computer.
In order to validate the HadCM3BB-v1.0model, we first compare modelled pre-industrial (0kyr BP, equivalent to the year 1950 without anthropogenic greenhouse gas emissions) precipitation and wind vectors against ERA-5 reanalysis data 11 (Supplementary Figure 10).Although the data is not directly comparable as it represents different time periods, past work has shown that these biases are small compared to the model biases 5 .Note that the model data has been downscaled to 1° resolution via bilinear interpolation.Broadly, the model accurately reproduces synoptic-scale wind patterns, the position of the winter and summer ITCZ, and the spatial pattern/intensity of the West African Monsoon (WAM) over land.In the summer, the model simulates too much precipitation across the Sahara on the order of 0.19mm/day, and very slightly underestimates winter precipitation (-0.04 mm/day).The average annual Saharan precipitation bias is 0.06mm/day equivalent to 21 mm/yr.This positive bias is because the model is prone to drizzle like many GCMs 5 , which has a relatively large impact in hyper-arid regions such as the Sahara.This bias will impact the amplitude of any simulated North African Humid periods (NAHPs) in this study, although it is small compared to the magnitude of those events.Outside of the Sahara, both summer and winter precipitation is of comparable magnitude across tropical and subtropical Africa.This is except for precipitation in the region of the Gulf of Guinea which is too high in the model.This may be a result of the stronger simulated regional summer westerlies, which may reflect inaccurate SAT/SSTs in a climatologically complex region which experiences strong upwelling.
However, over land the spatial imprint and intensity the WAM is comparable to observations.The spatial biases in the modelled precipitation, including across the Sahara, are comparable to other CMIP5 models as discussed in Section 5.1.2 of Valdes, et al. 5 (see their Figure 4).
Validating palaeo-model data is challenging due to the lack of reliable observational datasets which are themselves subject to extensive uncertainties and biases.Here we compare model timeseries for SATs, SSTs and precipitation against numerous global proxy datasets covering the past 0-800kyr (Supplementary Figures 11-13) in order to validate in the temporal domain.
The associated references are labelled in the Figure captions.
Supplementary Figure 11 shows 12 land-based temperature records including five Greenland ice-core records (b-f) [12][13][14][15] , two Antarctic ice core records (k, l) 16,17 one of which is a pan-Antarctic records comprised of 5 stacked cores (k) 17 , and five records from locations away from the ice sheets (g -j, m) [18][19][20][21][22] .Note that because our modelling approach does not incorporate millennial scale variability, the DO events that dominate the Greenland ice core records are not represented in our simulated dataset.This is particularly evident for the Bolling-Allerod/Younger Dryas interstadial/stadial transition in panels d-f.Despite this, there is generally good agreement in the broad scale patterns in modelled climate with the scale of temperature change in Greenland.However the modelled temperatures are generally too cold during the Holocene and warm beyond.The lack of a Holocene Climatic Optimum in the model has also been shown in a number of other climate models 23 , and may reflect biases across the current generation of models used and/or biases with the reconstructions themselves.In Antarctica, there is good agreement with the ice core datasets with regards to absolute temperatures, the timing of the glacial/interglacial cycles, and the scale of temperature change throughout the past 800kyrs.In Africa, the model does well at recreating the observed temperature in Lake Malawi and the Limpopo basin, particularly with regards to the timing of the glacial/interglacial changes.In Lake Tanganyika however the modelled temperatures are apparently too cold.Note however that lake derived temperature reconstructions are subject to a range of uncertainties.Similarly, a LOVECLIM model reconstruction of the Albanian Lake Ohrid indicates that the timing of large-scale temperature changes is well represented, however absolute temperatures are warmer in our model.
Supplementary Figure 12 shows a range of SST records with a wide spatial coverage, albeit generally focused at lower latitudes.In the North Atlantic (b-f) [24][25][26][27][28][29] , the model does generally well at simulating absolute temperatures and the broad scale timing of large-scale temperature changes.However, variability in the observational data is generally greater than that simulated due to the lack of millennial scale variability in the model dataset.There is good agreement with timing of the glacial/interglacial cycles and the extent of temperature change in panels d and f.In the South Atlantic and off Angola (g-i) [30][31][32] , modelled temperatures are too warm, but again the timing of temperature change is consistent.The west coast of Africa is a complex region climatologically due to the presence of strong regional upwelling which may not be properly resolved.In the Indian and Arabian ocean (j-l) [33][34][35] there is good agreement with modelled absolute temperatures and the patterns of temperature change, albeit with too little variability.There is very good agreement in the SST records in the South China Sea and off Australia (m-o) 33,36,37 .In the Pacific (p-u) 33,[38][39][40][41][42][43][44] modelled temperature are generally too warm, however the patterns of temperature change are broadly consistent with the observations.
Validating palaeo-precipitation is more challenging that SATs and SSTs due to the lack of comparable observational datasets.Palaeo-precipitation proxies such as lake levels, geochemical composition of sediments, and vegetation composition (pollen, δ 13 C) do not easily translate into quantitative precipitation estimates, but instead track relative changes in the precipitation or moisture availability.These proxies may also be influenced by additional factors other than precipitation which can hamper direct comparison.Furthermore, precipitation is problematic due to its significant spatial heterogeneity; therefore, observational records may demonstrate localised trends which may not be captured by a coarse scale model.For example, precipitation is sensitive to orography and the position of the coastlines, which is particularly subject to model resolution and may introduce biases.This may influence direct comparison in regions such as the Eastern African rift valley and around the complex coastline of the Mediterranean.
Supplementary Figure 13 shows a collection of Mediterranean and African proxies for moisture availability.In the Mediterranean and North Africa (b-c) 21,45 , the observation timeseries show Quercus (deciduous Oak) abundance (b) from Lake Ohrid in Albania and Ba/Al which tracks sapropel intervals from Eastern Mediterranean sediments (c).Both records indicate clear precession and eccentricity forcing, with periods of enhanced moisture availability at precession minima during interglacials and suppression during glacials.In both panels (b-c), the timing of enhanced moisture coincides with the NAHPs identified in our model timeseries (Figure 1), and both observations lend support to the conclusions of this study.There is a good correlation between the modelled precipitation in the Nile River Basin and high Mediterranean Ba/Al intervals.In the case of Lake Ohrid, comparison with model data is complicated because this region is not resolved in the model due to low resolution of the coastline, therefore the model timeseries shown in panel b is from a region to the North of Lake Ohrid.This is likely why there is disagreement with the model data, such as the lack of an eccentricity signal.However, precession timing is broadly consistent for both model and observations.
In the Western Sahara (d-e) 46,47 , modelled precipitation shows good agreement with patterns of dust deposition in marine sediments (d) and humidity change inferred from grain-size analyses of siliciclastic marine sediments (e).There is consistent precession timing in both the model and observations, with strong peaks in the observations during the interglacials and broadly weaker peaks during glacial precession minima.In tropical West Africa (f) 48 , pollen abundance from a sediment core off Liberia is of poor temporal resolution but demonstrates enhanced peaks at around 100kyr BP and 200kyr BP which coincides with interglacials, indicative of increased precipitation during these periods and increased aridity during glacials.
Lake Bosomtwe woody cover data (g-h) has recently been published with two different age models 49,50 yielding different temporal patterns.Both records show variability on a precession timescale, but the different age models appear to create significant differences in the larger scale patterning between the two versions of the proxy.This highlights the problems of model-proxy comparisons.In panel g, the interglacial timing corresponds more accurately with the model data, and there is also generally better correlation with the strong precession variability demonstrated by the model in this region.
In central Africa, terrigenous mineral composition of sediment cores from the Zaire fan (jk) 51,52 shows very good correlation with the model data, although precession is more dominant in the observation record in panel j.Both records and model demonstrate precession and obliquity orbital signals, and generally enhanced river run-off and moisture availability during precession minima and interglacials.
Eastern Africa, which includes the East African Rift valley (i, l-o) 19,[53][54][55][56] , is a complex climatological region influenced by highly variable topography and a number of atmospheric systems including the walker circulation, the East African monsoon and El-Nino.A record of Ca/Ti from Lake Tana in Ethiopia is indicative of lake level (i), and the authors state 53 that it shows evidence for a more stable wet environment during the Holocene and the last interglacial between ~132-70 ka BP, although this is not clear in the observations.The record shows some agreement with our model data, but it does not demonstrate the same dominant precession forcing.
A potassium-based aridity record from the Chew Bahir basin in southern Ethiopia (l) shows good correlation with the model data (note that this data has been obtained from the study of Kaboth-Bahr, et al. 57 as the original data is not available).The observations show clear dominant precession peaks modulated by eccentricity, the timing of which agree well with the model data.A 400kyr eccentricity signal is also apparent in both the model and observation data, and a trend towards increasing aridity since the last interglacial.A pollen record from Lake Magadi in Kenya (m) supposedly shows aridification of the region over the past 500kyr 55 , however the temporal resolution of the record is poor prior to 250kyr BP.This long-term aridification is not apparent in the model data, and the proxy does not show a strong trend either.Overall there is generally poor correlation between the model and Lake Magadi data.This is also apparent for the three precession minima peaks during the last interglacial between 70-120kyr, which are too closely grouped together in the observation dataset compared to the model and orbital data.
A leaf wax carbon isotopic record from the Kenyan Koora Basin (n) shows precession forcing strongly modulated by eccentricity, with particularly high amplitude variability during the penultimate interglacial between 180-270kyr BP.Although there is some disagreement between the model and the timing of some of the precession peaks, there is very good modeldata agreement with the eccentricity signal which indicates increased/decreased moisture availability during interglacials/glacials, and the enhanced impact of precession during the penultimate interglacial.A lake level record from Lake Malawi (o) shows, according to the authors 19 , evidence for 100kyr eccentricity with moist interglacial and dry glacials.However, the lake level data does not show precession forcing which is in contrast to the model data.However, the model does capture some of the megadroughts identified in the Lake Malawi data of Lyons, et al. 58 .
In Southern Africa (p-q) 59,60 , precession variability is antiphased with that in Northern and Central Africa, with enhanced precipitation during precession maxima (enhanced austral summer insolation).A leaf-wax δD record from a marine core off Namibia shows dominant precession forcing which correlates very well with model variability.Both the model and data indicate aridification in the region since the last glacial maximum.A Fe/Ca record from the Limpopo basin in South East Africa (q) show clear precession peaks that correlate well with our model data.There is also good agreement between the model and observations of stronger precession forcing during eccentricity maxima.However the model does not indicate increasing aridity between 800 to 600kyr as implied by the authors 59 , but this trend does not appear to be strong in the proxy data either.
In summary, the HadCM3BB-v1.0model accurately simulates present-day precipitation and wind patterns over Africa with biases comparable to other CMIP5 models, including the position of the ITCZ and WAM.Temporally, the model and the snapshot modelling approach that we utilise is shown to produce a climatology that is comparable to proxy datasets for SATs, SSTs and precipitation, in terms of absolute values and patterns of change over the past 800kyr.This model/data agreement lends strong support to the overall approach that we have employed.on the atmosphere has a dominant impact in reducing glacial precipitation, with a secondary impact from glacial sea-surface temperatures.The impact of glacial SSTs primarily reduces precipitation between 10-20°N whereas the impact of albedo change is primarily between 10-35°N.

Supplementary
Latitude/height plots (Supplementary Figure 8, Supplementary Figure 15a-b) indicate that during glacial precession minima periods, the ice sheets cool the atmosphere, weaken the West African Monsoon (WAM), strengthen (weaken) the African Easterly Jet (Tropical Easterly Jet), and cause the sub-tropical westerly jet (STWJ) to descend which inhibits precipitation.This is similar to the mechanism of precipitation suppression during precession maxima (Figure 2b).
The contribution of glacial SSTs (Supplementary Figure 15c-d) primarily weakens the WAM and alters the position of the African/Tropical Easterly jets, which decreases uplift and precipitation between 10-20°N (Supplementary Figure 15c).There is little impact on the dynamics of the sub-tropical westerly jet and only a small cooling impact on atmospheric temperatures (Supplementary Figure 15d).The broad impact of the ice sheets is to cool SSTs, specifically in the North Atlantic, sub-tropical gyre and the Mediterranean.This may influence the land-sea thermal contrast around West Africa and the Gulf of Guinea, which consequently influences the strength of the WAM.However the direct mechanism by which this occurs, and which oceanic regions are responsible for this, is beyond the scope of this study.Glacial SSTs also cause a more negative Atlantic meridional temperature gradient, although there is little impact on the global meridional gradient (Supplementary Figure 9).This may contribute in shifting the ITCZ to a more southerly position, as has been shown in other studies 61,62 .
The atmospheric anomalies due to the impact of glacial albedo (Supplementary Figure 15e-f) is very similar to that shown by the overall addition of the ice sheets (Supplementary Figure 15a-b).Albedo change weakens the WAM, alters the position and strength of the African/Tropical Easterly jets, and causes a northward shift, strengthening, and descending sub-tropical westerly jet.This inhibits uplift and convective rainfall between 10-35°N.This is due to the widespread cooling of the troposphere north of 10°N.This also results in a more negative global and Atlantic meridional temperature gradient (Supplementary Figure 9), which may act to shift the ITCZ to a more southerly position during glacials 61,62 .These results indicate that it is albedo change and the resultant tropospheric cooling which is responsible for much of the suppression of Saharan precipitation during glacials (Supplementary Figure 15f).
In summary, ice sheets suppress the NAHPs during period of precession minima primarily due to the thermodynamic impact of albedo change.This cools the atmosphere which reduces the meridional temperature gradient, weakens the WAM and alters the dynamics of the subtropical westerly jet which inhibits uplift.In addition, glacial SSTs suppress precession-induced strengthening of the WAM. the anomaly between the two (sensitivity experiment -Control).The region is identified in pink in Figure 2, and the specific experiments are Orb_GHG (panels a-b), glacial_noSST (panels c-d), glacial_noAlb (panels e-f), and glacial_noTopo (panels g-h).These are defined in each panel label.Panels a, c, e, g) the lagrangian tendency of air pressure for (filled colour contours; pa/s) where negative (positive) values represent upward (downward) motion, and zonal wind (black line contours; m/s) where solid/positive (dashed/negative) contours represent westward (eastward) air flow.The major components of the West African Monsoon (WAM) are identified in Figure 2b.Panels b, d, f, h) temperature (filled colour contours; °C), geopotential height (black line contours; m) and relative humidity (green line contours; %).Solid (dashed) contours represent positive (negative) anomalies.Anomalies are 95% confident according to a student t-test.

9 .
Annual mean Saharan (15:30°N, 15°W:35°E) and West Saharan (15:30°N, 15°W:15°E) precipitation, and annual mean global (poleward of 24°N and 24°S) and North Atlantic (poleward of 24°N and 24°S, 70°W:20°E) extratropic meridional temperature gradient (°C; North-South) for the 176kyr BP HadAM3B sensitivity experiments.The associated sensitivity experiments (see Methods) are identified on the Figure and clarify the thermodynamic and topographic impact of the ice sheet on annual average Saharan precipitation and the meridional temperature gradient.

Figure 10 .
A comparison of pre-industrial HadCM3BB-v1.0model output with ERA-5 reanalysis data.a) Winter and summer pre-industrial (0 kyr BP) HadCM3BB-v1.0output downscaled to 1° resolution for precipitation, wind and the position of the Intertropical Convergence Zone (ITCZ).b) Same as above but for ERA-5 reanalysis data 11 averaged over the period 1951-2020.The ITCZ is defined as where the northern and southern trade winds converge.

Figure 14 .Supplementary Figure 15 .
Thermodynamic and topographic Impacts of ice sheet addition on summer (JJA) precipitation at 176kyr BP. a) Summer (JJA) precipitation (mm/day) for HadCM3BB-v1.0NAHP_All, HadCM3BB-v1.0Orb_GHG, and the anomaly which represents the change in precipitation due to the addition of the ice sheet.b) same as above but utilising the Control HadAM3BB simulation to confirm the precipitation anomaly for the atmosphere only model, c) summer precipitation for Control, glacial_noSST and the anomaly representing change in precipitation due to addition of glacial SSTs, d) summer precipitation for control, glacial_noAlb and the anomaly representing change in precipitation due to addition of glacial albedo and its effect on the atmosphere, and e) summer precipitation for control, glacial_noTopo and the anomaly representing change in precipitation due to addition of glacial topography.West Saharan summer (JJA) latitude/height climatologies at 176kyr BP for the Control simulation, the thermodynamic/topographic sensitivity experiments, and

Supplementary Figure 13.
Validation timeseries plots for precipitation (mm/day).The same as Supplementary Figures 11 and 12 but for precipitation, with observation locations shown in panel a.All left y-axis labels are for model data in mm/day and right axis units are outlined in the panel label.The associated references are : b 21 , c 45 , d 46 , e 47 , f 48 , g 50 , h 49 , I 53 , j 52 , k 51 , l 56 ,m 55 , n 54 , o 19 , p 60 , q 59 .