Pattern scaling of simulated vegetation change in northern Africa during glacial cycles

Over the last hundreds of millennia, natural rhythms in Earth's astronomical motions triggered large-scale climate changes and led periodically to humid conditions in much of northern Africa. Known as African Humid Periods (AHPs), such times sustained vast river networks, vegetation, wildlife, and prehistoric settlements. The mechanisms, extent, and timing of the changes remain poorly constrained. Although AHPs along glacial cycles are recognizable in marine sediment records, the related land cover changes are difficult to reconstruct due to the scarcity of proxy data over the continent. Moreover, most available information covers only the latest AHP during the Holocene. Here we use a comprehensive Earth system model to look at additional, much earlier, possible cases of AHPs. We simulate the full last glacial cycle, aiming to reproduce the last four AHPs as seen in available proxies. The simulated AHPs seem in broad agreement with geological records, especially in terms of timing and relative strength. We focus on the simulated vegetation coverage in northern Africa, and we detect a dominant change pattern that seems to scale linearly with known climate forcing variables. We use such scaling to approximate northern African vegetation fractions over the last eight glacial cycles. Although the simple linear estimation is based on a single mode of vegetation variability (that explains about 70 % of the variance), it helps to discuss some broad-scale spatial features that had only been considered for the Holocene AHP. Extending the climate simulation several millennia into the future reveals that such (palaeo-based) pattern scaling breaks when greenhouse gases (GHGs) become a stronger climate change driver.

1 Introduction

Broad swings in environmental conditions over northern Africa have happened for at least the last 11 million years (Crocker et al., 2022). What is today the hyper-arid Sahara was formerly at times a much wetter region, adequate for sustaining perennial water bodies and vast expanses of vegetation cover that significantly reduced dust output to surrounding areas (Gasse, 2000; deMenocal et al., 2000). Such times, also known as African Humid Periods (AHPs), were recorded as quasi-periodic intervals of enhanced moist conditions in regional geological archives (Grant et al., 2017, 2022). Humid climates set in, along with changes in the land surface, which, via feedback loops, amplified the effects of changes in the Earth's orbit on the distribution of incoming solar radiation (Claussen et al., 2017; Pausata et al., 2020). The extent and timing of the changes remain uncertain due to proxy data limitations (Yacoub et al., 2023), while understanding of the causes and mechanisms is incomplete, as evinced in persistent discrepancies between existing geological records and climate simulations (Brierley et al., 2020).

The Holocene interglacial had the best-known case of an AHP, spanning roughly 11.7 to 4.2 ka (ka is thousands of years ago) (Claussen et al., 2017). Evidence collected for this period around northern Africa revealed the changes in past lake highstands (Lézine et al., 2011), river networks (Drake et al., 2011), dust fluxes (McGee et al., 2013), vegetation cover (Hély et al., 2014), and rainfall (Tierney et al., 2017b). In addition, numerical climate models were able to reproduce some of these changes during the Holocene AHP and hinted at basic mechanisms underpinning AHP dynamics (Claussen et al., 2017). Proxy data and models spurred an interesting debate that highlighted the spatial and temporal complexities of an AHP, since spatial heterogeneity resulted in marked differences in localized rates of change and asynchronous AHP termination across northern Africa (Brovkin and Claussen, 2008; Shanahan et al., 2015; Dallmeyer et al., 2020). Nevertheless, what happened before the Holocene and how known mechanisms played a role in even deeper times is relatively unknown. Repeated climate and land cover changes prevented the preservation of geological records older than the Holocene over land (Drake et al., 2018), while the costs of complexity and spatial and temporal resolutions limited in many cases the reach of simulations with comprehensive climate models (Tierney et al., 2020a).

What is known about pre-Holocene AHPs is mainly informed by signals recorded in marine sediments off the northern African coastline. Marine proxies, combined with the few surviving terrestrial signatures, and taking the Holocene AHP as a baseline, make it possible to examine the cases of AHPs along glacial cycles (Blome et al., 2012; Drake et al., 2013; Larrasoaña et al., 2013; Scerri et al., 2014). For instance, the multiple lines of geological evidence suggest that the previous interglacial during Marine Isotope Stage (MIS) 5e, around 127 ka, most likely had a much stronger AHP than the Holocene, in agreement with modelling results about the influence of changes in Earth's orbit. Likewise, marine proxies from around northern Africa show other potential AHPs during the last few glacial cycles (e.g. Rossignol-Strick, 1983; Ehrmann and Schmiedl, 2021a). However, limitations in projecting marine proxy signals far inland, compounded with scant terrestrial sources, restrict our understanding about the amplitude of northern African climate variability and AHP development. A clear outline of AHP changes is crucial to explain, for instance, regional geomorphology (Drake et al., 2022), ecological corridors during hominin evolution (Scerri et al., 2018), and long-term climate–land surface interactions (Pausata et al., 2020).

Climate models offer a theoretical window into possible spatial patterns of AHPs. For instance, early experiments with an atmosphere-only model already coarsely simulated four AHPs since 150 ka (Prell and Kutzbach, 1987), as many as indicated by marine proxies. However, modelling long glacial cycles to include multiple AHPs requires a compromise between model complexity, simulation length, and spatial resolution, which hinders the amount of detail that can be reproduced and evaluated. Studies that kept relatively complex models and higher spatial resolutions chose to either simulate discontinuous time slices (Singarayer and Burrough, 2015) or accelerate climate forcing rates (Kutzbach et al., 2008, 2020), in both cases reducing the effects of high-frequency variability. Other studies opted for intermediate-complexity models and large grid spacings to focus mostly on temporal patterns of average changes of AHPs (e.g. Tuenter et al., 2003; Tjallingii et al., 2008; Menviel et al., 2021). Although simulations broadly agreed with trends in the proxy records, underestimation of Saharan shrinking during AHPs and the associated rainfall increase is a consistent problem across models, leading us to question current understanding on the mechanisms and scale of the changes. Simulating a full glacial cycle with interactive climate and vegetation components at higher spatial resolution may reveal dominant AHP patterns underlying geological records or behind data–model mismatches.

We use the comprehensive Max Planck Institute for Meteorology Earth System Model (MPI-ESM; Giorgetta et al., 2013) to perform a transient global climate simulation since 134 ka at an affordable spatial resolution in order to see subcontinental changes in the vegetation cover of northern Africa. The main goal is to find general trends in simulated patterns during AHPs of the last glacial cycle. We evaluate the vegetation cover of different AHP responses and relate it explicitly with the large-scale climate drivers: (1) Earth's orbital parameters, (2) radiative effects of levels of atmospheric greenhouse gases (GHGs), and (3) high-latitude ice sheets. These drivers are known to set multiple controls on localized rain-bringing mechanisms (Dallmeyer et al., 2020; Blanchet et al., 2021), which could lead to diverse AHP biogeographies. Assuming local environments always reacted similarly to climate forcing, it is possible for AHP vegetation patterns to be approximated as a function of known forcing values. In addition, we also simulate and briefly assess the potential future evolution of northern African vegetation. In the following sections we describe the climate model, explain the experimental setup, present results focused on northern African vegetation cover during AHPs, and discuss the possibility of a forcing-based pattern scaling.

2 Methods

2.1 Model description

We use version 1.2 of MPI-ESM (Mauritsen et al., 2019). The model couples dynamical components representing the general circulations of the atmosphere (ECHAM6.3; Stevens et al., 2013) and ocean (MPIOM1.6; Jungclaus et al., 2013), along with changes in the land surface and vegetation cover (JSBACH3.2; Reick et al., 2013). The standard release of this model version participated in the Coupled Model Inter-comparison Project Phase Six (CMIP6; Eyring et al., 2016). Implementation of the coarse-resolution setup, T31GR30 (about 3.75° horizontally), with 31 atmospheric and 40 oceanic levels, enables long climate simulations at a maximum of about 700 simulated years per day (Mikolajewicz et al., 2018). Simulations covering several glacial cycles require updates during the runtime of the land–sea mask, glacier mask, and river routing that align with reconstructions of ice-sheet and sea-level changes. This is done once every decade, as described in Meccia and Mikolajewicz (2018) and Riddick et al. (2018). The specific model setup is actively maintained and performs well within the range of current general circulation models (Kapsch et al., 2022).

How the model computes natural vegetation cover is particularly relevant for this study. A detailed description is given in Reick et al. (2013), and a performance assessment is given in Brovkin et al. (2013). Briefly, the model uses a “mosaic” approach, where each land surface grid cell is tiled to be occupied by fractions of predefined plant functional types (PFTs). PFTs include eight cover classes between woody (tropical and extra-tropical trees plus shrubs) and herbaceous (C₃ and C₄ grasses) vegetation, each with different turnover timescales. When one vegetation type dies via natural causes (age) or via disturbances (fire or windthrow), the “void” it leaves is filled (or “colonized”) by the growth of the most competitive PFT in terms of net primary productivity (NPP), which is a function of climate and atmospheric CO₂. Seeds for all PFTs are assumed to be universally present. Fractions of PFTs can only occupy the portion of a grid cell that is hospitable for growth, which is determined based on having sufficiently productive (NPP) vegetation over consecutive years. The distribution of inhospitable fractions across the globe corresponds to hot and cold deserts. Vegetation cover in the model has immediate effects on the exchange processes between land surface and the atmosphere, which can subsequently alter climate patterns and feedback on vegetation (Brovkin et al., 2009).

2.2 Experiments

We simulate a transient Earth system response to prescribed forcing for the last 134 kyr (kyr is thousand years as in duration), almost entirely as a single continuous run (120 to 0 ka), except for a 14 kyr interval spliced a posteriori to extend analysis until 134 ka (details about splicing are in Appendix A). Earth's orbital parameters are computed from Berger (1978). Atmospheric levels of GHGs (CO₂, CH₄, and N₂O) are taken from the ice core record compilation of Köhler et al. (2017) and are updated every 10 model years as decadal means. Ice sheets from the GLAC-1D dataset (Tarasov and Peltier, 2002; Tarasov et al., 2012, 2014; Briggs et al., 2014; Abe-Ouchi et al., 2013) are interpolated spatially and temporally to suit model resolution and are also prescribed at decennial time steps, with corresponding adjustments in topography, bathymetry, and river routing (Meccia and Mikolajewicz, 2018; Riddick et al., 2018). Details about the GLAC-1D dataset extension until before the last interglacial are given in Menviel et al. (2019).

Figure 1 Transient simulation of the last glacial cycle. Forcing variables (left) represent prescribed conditions driving the experiment: (a) orbital parameters following Berger (1978) are shown using the orbital monsoon index of Rossignol-Strick (1983), (b) ice sheets from the GLAC-1D dataset (references in main text) are shown with the total ice volume in the Northern Hemisphere (NH), and (c) greenhouse gases (GHGs) from Köhler et al. (2017) are shown as radiative forcing (RF) changes from the pre-industrial (PI) era. Simulation results compared with proxy records (right): (d) changes in global mean near-surface temperature relative to the PI, compared to isotopic data from Lisiecki and Raymo (2005a); (e) Atlantic meridional overturning circulation (AMOC) in the Atlantic Ocean (1000 m depth, basin mean at 26° N); and (f) changes in precipitation in the Tropical North Atlantic (TNA; defined in the main text) compared to marine sediment reflectance (L^∗; smoothed with a third-order low-pass Butterworth filter with a cutoff frequency of $\mathrm{1}/\mathrm{5}$ kyr⁻¹) data from Deplazes et al. (2013a), both related to shifts in the Intertropical Convergence Zone (ITCZ). We include labels for marine isotope stages (MISs) according to Railsback et al. (2015).

Figure 1 (left) compiles a set of variables directly linked to forcing factors: orbitally driven changes in low-latitude insolation in terms of the so-called orbital monsoon index (Fig. 1a), the Northern Hemisphere total volume of ice sheets (Fig. 1b), and the GHG radiative forcing change relative to the pre-industrial (PI) era (1850 CE; Fig. 1c). The monsoon index of Rossignol-Strick (1983) is an insolation difference between the Northern Tropic (23.4° N) and the Equator during the caloric summer season defined by Milankovitch (1941). Changes in the radiative forcing of GHGs are based on an “equivalent-CO₂ concentration” explained in Ganopolski et al. (2010). We ignore any land use, and atmospheric aerosols remain at pre-industrial values (Kinne et al., 2013). Figure 1 also includes MIS labels according to Railsback et al. (2015), which we use throughout the text simply to refer to approximate ages of the late Quaternary when an AHP may have occurred.

Additionally, we run two transient experiments 10 kyr into the future. In both, the ice sheets remain unchanged from present-day conditions, while the orbital forcing continues to change according to Berger (1978). Such an approximation for the ice sheets is based on the study of Ganopolski et al. (2016), where they suggest no substantial changes are expected within the next several millennia. For GHG concentrations, we choose two alternative Shared Socioeconomic Pathway (SSP) scenarios from the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC): SSP1-1.9 and SSP3-7.0 (IPCC, 2023). The two scenarios are selected simply as contrasting cases of potential future GHG dominance. The GHG data for these scenarios come from Meinshausen et al. (2020) until the year 2500 CE, and data from 2500 CE onwards come from the previous modelling output of Brovkin et al. (2012), similarly to Kleinen et al. (2021).

3 Results

We average the simulation results as annual means at centennial scale. In the following sections, we present our results firstly from a simulated global context, and then we focus on the past and future of AHPs. We analyse the simulated AHPs of the last glacial cycle and develop a pattern-scaling approach that is also useful for exploring much older AHPs of the late Quaternary. Then we look at simulated future AHPs and how they compare to the pattern-scaling method. Whenever possible, we compare model output to proxy results and insights from previous modelling experiences.

3.1 Earth system response

Figure 1 (right) shows some globally relevant mean climate variables. Global climate is the warmest during the MIS 5e interglacial around 127 ka (Fig. 1d). About 0.7 K warmer than simulated pre-industrial temperature, the last interglacial warming matches proxies in the sign of change, but the amount could be underestimated by up to 1.3 K (Turney et al., 2020), similarly to other climate models (Lunt et al., 2013; Otto-Bliesner et al., 2021). Minimum average temperature occurs during the Last Glacial Maximum (LGM) around 24 ka, roughly 4.4 K colder than during the PI, which is not far from recent estimates ranging from 5 to 7 K (Tierney et al., 2020b; Osman et al., 2021). The amplitude and timing of temperature changes seem in agreement with the isotopic data from Lisiecki and Raymo (2005a).

Ocean dynamics maintain an Atlantic meridional overturning circulation (AMOC) around 19 Sv (Fig. 1e), including millennial-scale fluctuations related to freshwater input from glacier run-off as the ice sheets change. The strongest AMOC variability occurs during late deglaciation phases around 130–128 and 15–9 ka. The last deglaciation includes an especially strong meltwater pulse around 14.4 ka that nearly collapses the AMOC in the simulation (Kleinen et al., 2023). After this meltwater pulse, the remainder of the simulation includes some unexpected AMOC variability during the early Holocene. The amplitude of this early Holocene AMOC variability is not represented in the meltwater release from the ice-sheet reconstruction, and it was not seen in previous deglaciation experiments with the same model (e.g. Kapsch et al., 2022; Kleinen et al., 2023). In our simulation, it appears to be related to a slight difference in the initial conditions of ocean volume that results in some shallow regions in the Arctic Ocean. The AMOC changes are associated with drops in global temperature and in precipitation in the tropics (see Fig. 1d and f).

As a proxy for large-scale circulation changes in the atmosphere, we also show in Fig. 1f the simulated precipitation changes in the Tropical North Atlantic (TNA; 10 to 20° N, 80 to 20° W), where changes are coupled to movements in the average position of the Intertropical Convergence Zone (ITCZ). Unlike global temperature, these changes follow more closely the oscillations of the orbital monsoon index, and the timing of the minima of precipitation (when the ITCZ is southernmost) broadly agrees with the Cariaco Basin record of Deplazes et al. (2013a). The Cariaco Basin record shows total reflectance values of dark/light segments in a sediment core extracted off the coast of Venezuela which relate to rather enhanced (low L^*) or rather suppressed (high L^*) rainfall rates and hence tracked displacements of the ITCZ (Deplazes et al., 2013a). However, the simulation is slightly out of phase with the proxy record, and not all excursions from mean values match, probably because of regional- or local-scale signals in the data and the issues in smoothing both time series.

3.2 African Humid Periods

3.2.1 Timing and strength

Marine sediments off the coast of northern Africa reveal four AHPs in the last 134 kyr (Fig. 2a–d). AHPs appear as times of pronounced depletion of deuterium (D) ratios in leaf waxes from around eastern Africa (Tierney et al., 2017a) as a result of increased rainfall (Fig. 2a). AHPs can also be seen as times of reduced Saharan dust output to the Atlantic Ocean (Skonieczny et al., 2019) and Mediterranean Sea (Ehrmann and Schmiedl, 2021a) because increased soil-stabilizing vegetation cover restricted dust deflation (Fig. 2b, c). Likewise, during AHPs, increased river discharge into the Mediterranean stratified seawater and promoted the formation of sapropel bands (Fig. 2d). In spite of notable differences in location and methods amongst proxy records, AHPs stand out rather consistently.

Figure 2 Northern African climate in the last glacial cycle as seen in (a–d) marine proxies offshore and (e) average simulated vegetation coverage (woody, herbaceous, and total) in the Sahara (defined in main text). Proxy data (grey) are interpolated evenly at 1 kyr resolution and smoothed (black) with a third-order low-pass Butterworth filter with a cutoff frequency of $\mathrm{1}/\mathrm{10}$ kyr⁻¹: (a) deuterium ratio changes (δD) in terrestrial plant lipids (wax) in core RC09-166 near eastern Africa (Tierney et al., 2017a); (b) ²³⁰Th-normalized Saharan dust flux into the Atlantic Ocean according to core MD03-2705 near western Africa (Skonieczny et al., 2019); and (c) influx of fine-sized dust to the eastern Mediterranean Sea according to the ratio of kaolinite and chlorite minerals in core SL99, with (d) its sapropel sequence (Ehrmann and Schmiedl, 2021a). In panel (e), we also show a secondary y axis with the orbital monsoon index. Proxy sites are shown in an inset globe.

Focusing on the kaolinite/chlorite proxy (Fig. 2c), it reveals details about the intensity of the AHPs, according to the amplitude of dust pulses, which depends on the degree of (rain-fed) weathering that accumulates the kaolinite mineral in water bodies while an AHP is active (Ehrmann et al., 2017). Accordingly, in Fig. 2c, the AHP of MIS 5e stands out clearly as the strongest (largest dust pulse around 118 ka), while the estimated magnitude of other AHPs based on kaolinite/chlorite has been indicated to vary with the core-to-shore distance (Ehrmann and Schmiedl, 2021a).

Figure 2e shows simulation results of total vegetation fraction in the Sahara (20–30° N, 20° W–25° E). The simulation appears well in phase with trends in the geological records. We find four AHP-like changes in northern African climate, seen as large peaks in average Saharan vegetation. The vegetation peaks track the orbital monsoon index but with a small lag of about 1 kyr, which is longer than required for the vegetation to expand over the region. Most of this lag probably appears because the orbital monsoon index is tied to the summer season only and to the Northern Tropic (23.4° N), while the simulated vegetation response is also affected by what happens in the extra-tropics, not only in the summer but also during the winter season. Saharan vegetation peaks happen sometime between peak boreal summer temperatures and mild boreal winter temperatures.

Simulated AHPs also appear as times of herbaceous cover expansion and woody-type invasion of the Sahara. The AHP of MIS 5e is the strongest event with more than 15 % Saharan vegetation cover, followed by two weaker AHPs during MIS 5c and 5a (around 103 and 83 ka, respectively), both clearly over 10 % total vegetation fraction. In the simulation, the Holocene has the weakest AHP, barely approaching 10 % of vegetation in the Sahara. Such a weak Holocene AHP is unexpected according to previous applications of the same model to the last deglaciation (e.g. Dallmeyer et al., 2022).

In our simulation, the Holocene climate is unexpectedly largely affected by pronounced AMOC variability (see Fig. 1e), whose amplitude is greater than prescribed with the ice-sheet reconstruction changes (unlike MIS 5e, where AMOC fluctuations are according to ice-sheet changes). The large AMOC fluctuations appear connected with a small ocean volume difference in the initial conditions of the experiment. To diagnose the weak Holocene AHP, we include Appendix B, where we show the global climate anomalies with respect to pre-industrial temperature and precipitation. The effects of the AMOC variability are clearly seen in the early Holocene having a much colder North Atlantic region (Fig. B1f). This ensures that the ITCZ and related rain bands move southward and explains the anomalously weak Holocene AHP (see Tjallingii et al., 2008).

Nevertheless, the strength of AHPs in Fig. 2e is meant only relative to other AHPs, since an overall Sahara greening below 20 % is still substantially less than previous estimates from intermediate-complexity models (e.g. Menviel et al., 2021; Duque-Villegas et al., 2022) and less than estimates from the same model but implemented at higher spatial resolution (Dallmeyer et al., 2020). In Appendix C, there is a brief comparison of the simulated climate around the present day (mean values for 1990–2020 CE) with current observational products. We find the model simulates less vegetation in the Sahel (despite neglecting land use in the simulation) and is overall much drier (about 1 mm d⁻¹) over northern Africa than inferred from the observational data. Explaining the mismatch between different-complexity models or why there is a dry bias in the model at this spatial resolution is beyond the scope of this work; thus we focus on relative changes rather than absolute values.

3.2.2 Vegetation patterns

The intensity of an AHP should also manifest in the spatial extent of changes in vegetation coverage. Figure 3 shows maximum vegetation cover fraction in northern Africa during AHPs and changes from the PI. In this case, AHP values are selected from time slices around the vegetation optimum in the Sahara (see guide in figure). It is important to note that it is not trivial to define AHPs for the whole northern Africa, since it is a large region that responds with spatiotemporal heterogeneity (see Dallmeyer et al., 2020). For this reason, we show maximum values that may happen at different times (a few centuries earlier or later) for different sites (grid points) within an AHP time slice. Spatial heterogeneity leads to an interesting distinction about the changes in Fig. 3 between areas that are already “green” for pre-industrial conditions (like below the Sahel), which become “greener” during AHPs, and areas that are barren at the PI which turn “green” during AHPs (like in the Sahara).

Figure 3 Northern African vegetation coverage. (a–e) Maximum vegetation coverage during the pre-industrial (PI) era and AHPs and (f–i) the relative changes from the PI. A time series inset (top right) has marine isotope stage (MIS) labels to highlight mapped AHP time slices spanning about 4 kyr. On the maps, a horizontal line shows the Northern Tropic (23.4° N). Only the PI map shows modern coastlines, while others show the coastlines according to reconstructed sea-level changes. Lake Chad is shown in all maps as a point of reference.

Pre-industrial conditions in the simulation (Fig. 3a) already seem to underestimate zonal patterns of vegetation (like present-day dry bias in Appendix C), resulting in little vegetation in the Sahel (around 15° N) and an exaggerated Dahomey Gap (dividing the coast of Guinea) that is much less conspicuous in higher-resolution versions of the same climate model (e.g. Bathiany et al., 2014). During AHPs, there is a continent-wide increase in vegetation, closely approaching the Northern Tropic (23.4° N) and effectively shrinking the Sahara (Fig. 3b–e). The largest changes occur around the Sahel, within 10 to 20° N, with increases of over 40 % (Fig. 3f–i). The smallest changes seem to be cornered in the northeastern Sahara. Near the Equator, in places where there is already close to 100 % vegetation cover under pre-industrial forcing, there are also small or no changes during AHPs.

AHP patterns in Fig. 3b–e seem to become weaker sequentially, having the most intense changes during the AHP of MIS 5e on the left, followed by two relatively weaker intermediate cases, and ending with the mildest AHP case during the Holocene (around 10 ka during MIS 1) on the right. Most apparent is the fading of intense changes around the Sahel and a less pronounced northward push along the west coast. Such step-wise (ignoring long time gaps between AHPs) weakening of AHPs is also evident in the average vegetation fraction in the Sahara (see Fig. 2e). Incremental change occurs in the same direction as the orbital forcing (see Fig. 1a), thereby suggesting the possibility of a linear scaling of a baseline pattern with climate forcing variables. In the following section, we explore such a linear relationship for vegetation cover fractions.

3.3 Spatiotemporal analysis

AHPs seem quasi-periodic in response to natural rhythms of the Earth system. Possibly, there are preferred modes (patterns) of AHP variability (i.e. climate and biogeographies) that regularly appeared in the past during AHP development. In such cases, it may be assumed that the broadest patterns would tend to scale proportionally to large-scale climate forcing conditions. We investigate this possibility of pattern scaling in the simulated vegetation fractions in northern Africa during the last glacial cycle, with the main focus on the potential to capture with this method the simulated AHP patterns of Fig. 3.

We find a baseline pattern of northern African vegetation cover variability with an analysis of empirical orthogonal functions (EOFs) and their related principal components (PCs) for the entire last glacial cycle. In the analysis, the AHPs are simply the cases of maximum possible vegetation coverage. For this, we construct a covariance matrix with changes from the PI of annual maximum vegetation cover fraction for grid cells inside northern Africa (m=126) at centennial temporal resolution (n=1340). The covariance matrix (126×1340) is centred about its mean to remove trends and is factorized with a singular value decomposition to obtain its eigenvectors ordered according to explained variance. After decomposition we use only the first eigenvector (EOF₁), which, together with its corresponding time series (PC₁), explains the greatest variability.

Figure 4 EOF analysis of northern African changes in vegetation fraction from the pre-industrial (PI) era: (a) leading variability mode (EOF₁) and (b) corresponding PC₁ time series scaled to unit variance. PC₁ is smoothed with a third-order low-pass Butterworth filter with a cutoff frequency of $\mathrm{1}/\mathrm{5}$ kyr⁻¹. Panel (b) also includes a forcing-based estimate for the smooth PC₁ (central trend and spread due to uncertainty of regression parameters in parentheses) whose linear regression model equation is shown, where O_lag1k is the orbital forcing monsoon index lagged 1 kyr (cf. lag in Fig. 2e), G stands for the GHG radiative forcing change from the PI, and I_NH is the ice-sheet volume in the Northern Hemisphere. Modern coastlines are shown in panel (a).

Figure 4 shows the results of the EOF analysis. The dominant pattern of variability (EOF₁) explains about 73 % of the variance in vegetation fraction changes (Fig. 4a). It is a mostly zonally coherent pattern with an action centre largely inside 10–20° N. Near eastern Africa, the pattern adopts a more meridional structure between 30 and 45° E. According to its corresponding PC₁ (smoothed in Fig. 4b), this mode is related to uniform changes in vegetation fraction with respect to pre-industrial conditions, which are positive and of great amplitude usually during warm interstadials (shaded MIS boxes). PC₁ also seems in phase with the orbital forcing monsoon index (compare Fig. 1a). The phasing with the monsoon index is a good indicator of the possibility of estimating PC₁ as a function of external climate forcing.

Figure 4b also shows a linear model that combines the orbital monsoon index (lagged 1 kyr and symbolized as O_lag1k), the GHG radiative forcing (symbolized as G), and the Northern Hemisphere ice-sheet volume (symbolized as I_NH). The model is fit to the smooth PC₁ using an ordinary least-squares method. We choose to lag the orbital monsoon index as a predictor because of the lag in the vegetation peaks explained in Sect. 3.2.1. Fitting is done several times in different parts of the time series, using a sliding (every 10 kyr) window of 63 kyr (about three precession cycles), which explains the uncertainty in the regression parameters.

Overall, the linear model (_{forcing-based}PC₁ in Fig. 4b) estimates a similar trend to the smooth PC₁, with some uncertainty in the parameters that leads to some spread. It is possible that other more complicated fits could reproduce the smooth PC₁ pattern better, or with less spread, but we restrict our analysis to the simplest fit. In the linear equation, it is interesting to note that, for the last glacial cycle, the standard deviation of O_lag1k is close to 10.1 W m⁻², while that of G is around 0.8 W m⁻² and that of I_NH is around 11.0× 10⁶ km³. Multiplying these values by the absolute value of their respective fit parameters (0.07, 0.18, and 0.02), it shows that orbital, GHG, and ice-sheet variables influence about 66 %, 13 %, and 21 %, respectively, of the outcome of the linear model. The dominance of the orbital forcing is likewise clearly visible in the oscillating shape of the linear fit. The fitting seems to perform most poorly after 60 ka, when the orbital forcing becomes much weaker.

3.4 Pattern scaling

We can use the results of the previous section to obtain forcing-based estimates of the vegetation changes (relative to the PI) in northern Africa. Our method resembles a pattern scaling, where an invariant spatial pattern is assumed to scale linearly with a known forcing factor. The relationship to the climate forcing is based on the curve fit to the smooth PC₁ time series in Fig. 4b, using the linear sum of the orbital monsoon index (lagged 1 kyr), GHG radiative forcing, and Northern Hemisphere ice-sheet volume. Because of the use of the principal component time series, another name for this method is also principal component emulator (Wilkinson, 2010), which has also been used before in palaeoclimate modelling contexts (e.g. Holden et al., 2019). Using EOF₁ and _{forcing-based}PC₁, we can approximate mean vegetation changes having forcing variables as predictors. We are especially interested in estimating the mean changes during AHPs.

Figure 5 Assessment of the pattern-scaling method for estimation of change in vegetation fraction from the pre-industrial (PI) era: (a–d) mean spatial errors during the AHP time slices highlighted in panel (e); (e) domain-wide average performance. This assessment uses the MPI-ESM simulated output as a reference (i.e. errors are scaled minus simulated patterns). In panels (a)–(d), the root-mean-square errors (RMSEs) are shown. Maps include modern coastlines only as a spatial reference.

In Fig. 5, we evaluate the skill of the pattern scaling to predict the simulated vegetation changes (from the PI) during the same AHP time slices as in Fig. 3. This means the climate model output acts as reference and that errors come from scaled minus simulated fields. We find that, for the AHPs of MIS 5 (Fig. 5a, b, and c), the pattern scaling underestimates the vegetation changes inside the Sahara, which is most noticeable in western Africa. Only for the Holocene (Fig. 5d) is there a slight overestimation of vegetation growth during the AHP. That our approach projects only limited changes in the Sahara (e.g. north of 20° N) is expected because the EOF₁ pattern (Fig. 4b) has more pronounced loadings in the Sahel (around 15° N). We found other EOF patterns (not shown) that could have provided more changes in the Sahara; however, their associated temporal patterns were not simple to link to the external climate forcing. Nonetheless, considering the magnitude of some of the changes in northern Africa (above +40 % in Fig. 3), root-mean-square errors (RMSEs) below ±10 % vegetation cover change in Fig. 4a–d may be seen as acceptable.

Furthermore, Fig. 5e shows the average performance of the pattern scaling. We find a similar outcome: that the method predicts intermediate values between the rather strong AHPs during MIS 5 and the weaker Holocene AHP. We hypothesize that, if there had not been the unexpectedly strong AMOC variability during the early Holocene (also seen in Fig. 5e), which weakened the Holocene AHP, our method would find a curve fit able to predict slightly larger values closer to those during MIS 5. It is also interesting to see that the scaling spread, due to uncertain regression parameters, keeps the climate simulation largely within the reach of the pattern scaling.

3.5 Late Quaternary AHPs

The pattern-scaling method opens up the interesting possibility of quickly obtaining a rough estimate of the vegetation fraction distribution in much of northern Africa for other times for which reliable forcing data are available (and yet may lie outside of current modelling capabilities). It seems best suited for the Sahelian vegetation response (where the main loading of the EOF pattern lies); however, it still provides part of the vegetation change everywhere in northern Africa (e.g. the Sahara), such as a minimum probable estimate of the change. Two main assumptions underlying this approach are that (1) the EOF₁ pattern captures much of the large-scale pattern of regional variability and that (2) the forcing-based linear model effectively emulates the long-term time dependence of variation at that scale and the climate forcing. Although it may likely be biased towards a single glacial cycle, proxy data suggest the last glacial cycle could be representative of large-scale climate changes happening for at least the last eight glacial cycles (Lüthi et al., 2008a). Therefore, it is plausible that the pattern scaling generalizes to many other cycles.

Figure 6 Northern African vegetation fraction changes over the last 800 kyr. (a–c) Set of forcing factors used for scaling, including (a) in the Rossignol-Strick (1983) threshold for sapropel formation at 19.8 W m⁻². (d) Composite sapropel record (cores ODP160-966/7), with wide, dark bands for clear sapropels and narrow, light bands for oxidized or missing sapropels (Emeis et al., 2000a). (e) Inferred AHPs (colour bands) from bulk elements in sediment core ODP160-967 (Grant et al., 2022). (f) Pattern-scaling results for area average changes from the pre-industrial (PI) era. A horizontal line at 14 % is an empirical threshold for counting potential AHPs based on AHPs of the last glacial cycle. (g, h) Scaled patterns during the potential ancient AHP of MIS 15a showing modern coastlines, even though they may differ at that time.

We use the scaling approach to estimate changes in northern Africa since 800 ka using available forcing data. The orbital monsoon index is computed for this interval based on the parameters of Berger (1978). The GHG radiative forcing and ice-sheet data come from Ganopolski and Calov (2011), who compiled the trace gas data from Antarctic ice core records (Petit et al., 1999; EPICA Community Members, 2004) and simulated the ice sheets with a coupled intermediate-complexity model and ice-sheet model. Figure 6a–c shows this longer set of forcing variables, where the new GHG and ice-sheet data agree with the forcing used in our simulation and to fit the linear model for the overlapping period of the last 134 kyr. In Fig. 6a, we also show, as a guide, the orbital monsoon index threshold at 19.8 W m⁻², which relates to the occurrence of AHPs as detected by Rossignol-Strick (1983) and discussed more recently in Duque-Villegas et al. (2022).

Evidence of the AHPs of the last 800 kyr is available in marine proxies. Figure 6d and e depict Mediterranean proxy records linked to past AHPs: the sapropel record of Emeis et al. (2000a) and the bulk composition elements in the sediment core of Grant et al. (2022), respectively. Both long records agree on many events, although the geochemical analysis of Grant et al. (2022) detects several more AHPs. We compare such records with our simulation and pattern-scaling approach.

Figure 6f shows the last eight glacial cycles of average vegetation changes in northern Africa according to the pattern scaling and according to the forcing variables shown in Fig. 6a–c. Figure 6f also includes a reference line set at 14 % of average vegetation change, based on the AHPs of the last glacial cycle, which is useful to guess potential more ancient AHPs. According to this reference guide, two other potential AHPs stand out, possibly of a similar strength to the AHP during MIS 5e: one happening around 580 ka, within MIS 15a, and another one around 330 ka, within MIS 9e. The Mediterranean records also clearly contain these two events (Fig. 6d and e).

The average vegetation changes in Fig. 6f also compare well with the BIOME4 dataset generated with the statistical method in Krapp et al. (2021), although the reconstructed spatial patterns have differences in the zonal extent of the changes (not shown). Moreover, Fig. 6g and h show the estimated changes from the PI and vegetation fraction pattern, respectively, during one of the potentially strongest AHPs (around 580 ka). It is only one example of the pattern scaling for beyond-simulated AHPs. When considering such an ancient AHP, caution is important, and it is useful to know, for instance, that the land–sea mask and topography may differ from the ones used to calculate the EOF₁ pattern.

Based on the empirical threshold line in Fig. 6f, the scaling approach can identify most clearly some 19 (potential and known) AHPs. This number agrees with the 21 possible AHPs based only on the orbital monsoon index threshold (see Fig. 6a) detected by Rossignol-Strick (1983). Considering the spread in the parameters of the linear fit, the number could be as low as 8 (including the Holocene AHP, despite its minimum spread value being under the empirical threshold) or as high as 21 (like the monsoon index threshold would predict). All 21 of the candidate AHPs would be in agreement with signals in Mediterranean sediment cores (combining Fig. 6d and e), neglecting small lags between monsoon index peaks and sedimentary signals. A similar number of 20 AHPs was reported from many independent equilibrium simulations along this time range (Armstrong et al., 2023). Likewise, a data synthesis by Larrasoaña et al. (2013) indicated 22 potential AHPs for this interval. The refined geochemical analyses by Grant et al. (2022), which read subtle humidity variations, estimated about 30 episodes of Saharan greening since 800 ka (see Fig. 6e). However, some of the events could perhaps be grouped into single cases of AHPs.

3.6 Future AHPs

Previous studies pointed out the differences between orbitally dominated AHPs in the past and predicted, GHG-modulated AHPs in the future (Claussen et al., 2003; D'Agostino et al., 2019; Duque-Villegas et al., 2022). For the next 100 kyr, Earth keeps a low eccentricity orbit (a cycle of about 400 kyr) that hinders the possibility of a new AHP inception. At most, there could be some mild greening in northern Africa around 66 kyr into the future (Duque-Villegas et al., 2022). The low eccentricity reflects in the small variation in the orbital monsoon index in Fig. 7a, which stays far from the threshold of about 19.8 W m⁻² for the onset of an AHP detected by Rossignol-Strick (1983) and Duque-Villegas et al. (2022) (compare Fig. 7a with Fig. 6a). However, when compounded with increasing emission rates of GHGs, the strengthening of the planetary radiative forcing (as shown for the scenarios in Fig. 7b) enhances humidity levels globally, possibly enabling an earlier greening in northern Africa (Duque-Villegas et al., 2022).

Figure 7 Future vegetation coverage change in northern Africa: (a) orbital forcing monsoon index and (b) GHG radiative forcing scenarios; (c) simulated and scaling-based average vegetation fraction change from the pre-industrial (PI) era, where grey bands show the scaling spread; and (d, e) mean simulated and (f, g) scaling-based vegetation fraction under forcing scenarios around 9 kyr after the PI.

Scenarios in Fig. 7b have peak emissions in the near future (a few centuries) and later maintain a constant positive forcing, leading to an average global warming of about 7 K (SSP3) and 1 K (SSP1). Figure 7c shows simulated northern African vegetation response to those scenarios. SSP3 seems to arrive at AHP conditions with a similar level to that during MIS 5e (about 30 % average change from the PI), while SSP1 at about 10 % average vegetation change (which is smaller than during the Holocene) does not resemble an AHP (note the threshold at 14 % in the previous section).

Figure 7c also shows the average performance of the pattern scaling for future times. It performs relatively well for the SSP1 scenario, but it greatly underestimates the SSP3 changes (even considering the parameter spread). This means the method is rather conservative in terms of the GHG forcing, which could be related to the fitting being done within the narrow range of GHG variability of the last 134 kyr or could be because the model also includes the glacial changes in the ice sheets, which, for these future scenarios, are assumed to remain invariant at present-day conditions. Moreover, our estimate method is bound to the EOF₁ pattern, which means it is not suited to account for greening related to possible future shifts in vegetation boundaries.

Figure 7d–g show the patterns of vegetation change (from the PI) that were dynamically simulated (with MPI-ESM), and those estimated with the pattern scaling, near the orbital monsoon index local maximum at 9 kyr after present. This local maximum in tropical insolation could see substantial greening in northern Africa, potentially developing the next AHP, depending on the GHG forcing. Under the relatively low GHG forcing scenario, SSP1, the pattern scaling (Fig. 7f) performs well in comparison to MPI-ESM (Fig. 7d). However, Fig. 7e shows that, under the high GHG forcing scenario, SSP3, the simulated vegetation change pattern differs from its AHP predecessors, since it is missing the northward vegetation push along the west coast. Such a difference agrees with first-order effects seen in previous comparisons of climate change projections and the Holocene AHP (e.g. Gaetani et al., 2017; D'Agostino et al., 2019), although such studies did not extend until 9 kyr after present. Evidently, our single-EOF pattern-scaling approach by design lacks in skill to predict vegetation patterns simulated under high GHG forcing, since it is based only on late Quaternary vegetation changes.

4 Discussion

We present results of a coupled atmosphere–ocean–vegetation simulation for the last glacial cycle. To the best of our knowledge, the simulation is unique regarding model complexity and length and the absence of any acceleration technique. Other studies with similar times of interest opted for atmosphere-only models (Prell and Kutzbach, 1987), lesser-complexity models (Blanchet et al., 2021; Menviel et al., 2021; Duque-Villegas et al., 2022), ensembles of discontinuous equilibrium simulations (Singarayer and Burrough, 2015; Armstrong et al., 2023), or acceleration methods (Kutzbach et al., 2008, 2020). In this case, the simulation includes a realistic representation of interactions across a broad spectrum of variability frequencies, which means disturbances can propagate across the globe to superimpose climate forcing and local effects. However, for the same reason, it is not possible to implement a highly resolved grid in this long palaeo-experiment, which would enable a more refined comparison with the geological record and would more clearly capture shifts in vegetation boundaries (e.g. Dallmeyer et al., 2020).

We focus our analysis on vegetation cover in northern Africa, where continental-scale climate changes produce four AHPs and a northward greening reaching into the Sahara. The simulated vegetation changes generally agree with the timing and direction of change seen in available terrestrial and marine proxies (Blome et al., 2012; Drake et al., 2013; Tierney et al., 2017a; Skonieczny et al., 2019; Ehrmann and Schmiedl, 2021a). However, we detect a dry bias in the model that means Saharan greening is likely underestimated in our simulations. Most of the simulated greening happens within the Sahel region. Especially for the Holocene AHP, vegetation cover appears largely underestimated in our simulation because of unexpectedly large AMOC fluctuations, which resulted in a much colder North Atlantic and Northern Hemisphere while this AHP developed. Nonetheless, we keep this weakened Holocene AHP in our analyses because we still see considerable northern African greening, even if only for a few centuries.

Compared with previous modelling and reconstruction studies, the pattern of vegetation changes in the four simulated AHPs also seems consistent with recognized aspects of the Holocene AHP (Brierley et al., 2020; Dallmeyer et al., 2020). There is a northward progression of mostly zonal vegetation bands that extend east until about 40° E, with a slight northwest–southeast tilt related to further inland penetration of the West African Monsoon system. A “tilted” AHP development pattern is also detected in the simulations and reconstructions of Armstrong et al. (2023), and our results agree on the possibility that desert cover remains in a northeastern corner of northern Africa during AHPs, thereby breaking the zonal pattern of northward vegetation growth. Moreover, it was shown that palaeo-conditioning, or re-tuning a climate model to improve palaeo-modelling results, also led to a similar tilted vegetation pattern during the Holocene (Hopcroft and Valdes, 2021). In our case, however, such tuning is not applied.

Nevertheless, we also see differences with some of the few available sources that include not only the AHP of the Holocene, but also that of MIS 5e. The reconstructions of Larrasoaña (2021) show much more extensive greening of northern Africa and a more zonal pattern than our simulation. Also, for the AHP of MIS 5e, we are unable to simulate enough greening inside the Sahara, which would better fit the precipitation changes in the model–data comparison of Scussolini et al. (2019). We believe that part of the reason why these differences stand out is the dry bias in the model at the implemented resolution, but it is also possible that uncertainties and the sparse coverage in the geological record skew the reconstructions (see Otto-Bliesner et al., 2021).

Besides producing and evaluating vegetation patterns for the AHPs of the last glacial cycle, a most interesting outcome is the linear relationship we establish for simulated vegetation in northern Africa and the climate drivers forcing the simulations. We had found such a strong correlation previously but with a different climate model of intermediate complexity and only for Saharan greening (Duque-Villegas et al., 2022). Similarly to our previous work, the dominance of the orbital forcing is clear in our results, yet now, with the more comprehensive model, we are also able to extend this dependence to the spatial pattern predicted over northern Africa by the dynamical model via a simple linear scaling. The result is a cost-efficient estimation of northern African vegetation patterns for the late Quaternary, since simulating glacial cycles even at coarse resolution can take excessively long computing times.

We can compare the results of our pattern scaling with the global biome dataset generated by Krapp et al. (2021) for the last 800 kyr. In Krapp et al. (2021), authors also use model output for the last glacial cycle to create a linear regression model that can extend simulated variables much further back in time until 800 ka. Similarly, the predictor variables are forcing terms like orbital parameters and atmospheric CO₂. However, in that case, the approach was implemented globally and the estimated changes in northern Africa seem more zonal in comparison to our results. Nonetheless, the average evolution of total vegetation changes (one minus the fraction of desert plus barren biomes) agrees well with the domain-average performance of our method, although the estimated number of AHPs seems to be lower in their case.

An important caveat must be considered in our pattern scaling of vegetation fraction. Because this quantity is bounded (0 % to 100 %), it can happen that, in some cases, the scaling “overshoots” the estimation of the changes, since there is some maximum possible amount of change (from the PI) leading to completely bare (0 %) or vegetated (100 %) grid cells. Overshooting could happen, for instance, under more intense forcing (than in the last glacial cycle) or in cases where there is a greater sensitivity to the forcing. This implies the linear scaling may work well for vegetation fraction in this study only because of a dry bias in the model that prevents simulation of stronger (non-linear) greening.

Another factor we neglect in the scaling approach is the role of non-linear effects. Critical thresholds related to the climate forcing may prevent the pattern-scaling approach from correctly representing parts of the changes in the domain (Brovkin et al., 1998; Claussen et al., 2003; Duque-Villegas et al., 2022), most likely near transition zones, such as around the southern fringe of the Sahara, or where there is steep topography, such as towards eastern Africa. Those are some of the places that stand out when we assess the errors in the pattern scaling (see Fig. 5). Moreover, the dynamical systems that govern northern African moisture changes are connected via Walker and Hadley circulations with remote locations, whose responses to changes in the background climate also could induce more or less greening during simulated AHPs (e.g. Kaboth-Bahr et al., 2021). In this case, we have not looked into such teleconnection mechanisms. Although it is possible that some of the other coherent modes of variance improve our estimates (i.e. reduce emulation minus simulation differences), we opt for the simplicity of capturing only the variance at the largest scales, at the cost of a mean error of about ±10 % vegetation fraction change. The fact that the dominant EOF explains about 73 % of the simulated variance justifies our unimodal approach.

In our simulations of millennial-scale future climate change, we see the contrast between GHG-enabled and orbital-driven vegetation patterns of AHPs. The vegetation pattern under high GHG forcing (Fig. 7e) is similar to 21st-century projections of Sahel rainfall changes (e.g. Bony et al., 2013; Gaetani et al., 2017). In both cases (high GHG and 21st-century projections), the largest positive anomalies concentrate near central and eastern Sahel. However, in our case (high GHG), a key distinction is that the orbital configuration is changed considerably from the present day; thus we also see some vegetation growth in western Africa and approaching the Northern Tropic. Still, the western and northward development of vegetation under high GHG forcing is not as extensive as in the orbital-driven AHPs of the late Quaternary. This is related to the differences in dynamical (circulation) and thermodynamical (water vapour) responses to each of these forcing factors (Claussen et al., 2003; D'Agostino et al., 2019). Such competing effects complicate AHP pattern predictions in future climates.

5 Conclusions

After simulating a realistic climate evolution for the last glacial cycle, we produce new estimates of the vegetation cover patterns during the last four AHPs. A comparison with available data sources from the geological record and relevant modelling studies indicates that the estimated patterns capture some known continental-scale features of AHPs reasonably well. Nevertheless, a dry bias in the model version at the resolution we use frames our results rather as probable minimum estimates. We reduce the high-dimensionality in the variability of the vegetation cover changes in northern Africa to a single dominant pattern that accounts for 73 % of the variations. We find this single dominant pattern correlates well with known forcing variables, and we use this linear relationship to approximate to an acceptable extent the simulated vegetation changes of glacial cycles. Considering the high computational costs of running glacial cycle simulations, we think our approach could be used as a quick diagnostic tool for palaeovegetation cover in northern Africa, which we demonstrate by estimating vegetation cover changes for the last eight glacial cycles with considerable success. Additionally, the extension of our analysis several millennia into the future suggests that patterns of northern African greening under projected strong GHG forcing should differ enough from the vegetation patterns of late Quaternary AHPs that it would not be useful to assume a simple linear scaling with climate forcing for predictions.

Appendix A: Backwards extension of transient simulation

The complete 134 kyr simulation is mostly a single experiment that was run continuously from 125 to 0 ka. This interval is planned for studies about global climate of the last glacial cycle. However, because for this study we are particularly interested in the regional climate of northern Africa, we choose to complement this long experiment with another shorter, independent, and overlapping experiment that has an initial date slightly earlier, spanning 134 to 120 ka. We could use this to extend the longer experiment backwards in time. The main reason is to include in our analysis the tropical insolation maximum around 127 ka, which, according to available evidence, is an important time of extensive humidity over northern Africa. We opt for such an extension due to the computational expense and long waiting time (around 2 years) that it would take to re-run an entirely new simulation. The overlap in the experiments is to measure their level of agreement near a merging point. They should neither drift away from each other nor be at completely different states.

We find the splicing together of both experiments to have only a small consequence, producing a low-amplitude change that can be ignored for the study of the long-term climate signals we focus on for the study of AHPs. Figure A1 shows an example of the state evolution of both experiments for one grid cell in a highly variable region in the North Atlantic (Fig. A1a) and at the centre of our region of interest (Fig. A1b). The experiments are spliced together at 120 ka to produce the long simulation starting at 134 ka.

Figure A1 Extending the simulation deeper in time by merging overlapping experiments. Comparison of equilibrium and transient experiments for single grid cells (a) over the North Atlantic (35.26° N, 33.75° W) and (b) over land in northern Africa (20.41° N, 7.5° E). Experiment “mdv0003” spans 134 to 120 ka, while experiment “pmt0590” spans 125 to 0 ka. We zoom time on the interval close to the merging point. Circle markers show the end state of the equilibrium experiments, which is also the starting point of the transient experiments.

Appendix B: Background climate near simulated AHPs

To be able to compare more objectively the differences between simulated AHPs, we include Fig. B1, where we show the global conditions of near-surface temperature and precipitation near the times of the AHPs. The pre-industrial (PI) conditions (Fig. B1a), around the year 1850 CE, show the typical poleward temperature gradient profile and the precipitation maximum around the Equator. In Fig. B1b, we see that the climate during MIS 5e is overall much warmer than during the PI, except around the Sahel in northern Africa, where there is more transpiration because of the extended vegetation cover during the AHP. MIS 5c and MIS 5a show very similar conditions, also with slightly warmer conditions in the Arctic Ocean and Siberia. In great contrast, the Holocene (MIS 1) shows much colder temperatures over Greenland and the North Atlantic Ocean due to an excessive amount of sea ice that remains in this region after a strong meltwater pulse around 14.4 ka. Because of this, the Holocene AHP we present is not as green as seen in previous applications of the same model at the same spatial resolution. The precipitation changes (Fig. B1c) clearly show the shift in moisture and ITCZ towards the Northern Hemisphere during the AHPs. The largest increments in precipitation happen not only in northern Africa, but also over the Indian Ocean and over India. Precipitation changes during the Holocene AHP are much weaker than during the other AHPs.

Figure B1 Background climates during AHPs: the pre-industrial (PI) era means of (a) near-surface temperature and (b) daily precipitation rate fields, (c–f) temperature changes from the PI, and (g–j) precipitation rate changes from the PI. Modern coastlines are shown for the PI, while AHP time slices show the coastlines according to sea-level changes. An inset axis with global mean surface temperature highlights the AHP time slices.

Appendix C: Comparison to present-day observations

For an assessment of the dry bias in this specific implementation of the MPI-ESM model, we also compared the simulated present-day conditions of vegetation cover fractions and precipitation with some observational products in Fig. C1. In this case, “present day” is defined as the climatological mean of 31 years of model output, corresponding to the interval from 1990 to 2020 CE. The observational products were remapped to the model grid using a bicubic interpolation and averaged in time to obtain climatologic values. It is important to note that the simulations did not include any anthropogenic land-use changes; therefore, in many places (like near the Equator), it is expected that the simulated vegetation cover is an overestimation (potential vegetation cover). Figure C1a–c show that the patterns of vegetation are similar but that the model underestimated the vegetation cover around the Sahel (in spite of not including land use changes). A dry bias is clearly seen in the precipitation differences in Fig. C1d–f, especially around the Sahel. The large positive anomalies at the Equator (almost the only positive anomalies) indicate that the maximum of precipitation in the simulation may be displaced southward, relative to the observational product data.

Figure C1 Comparison in northern Africa of simulated present-day (1990–2020 CE) conditions to observational products for (a–c) vegetation cover fractions and (d–f) daily precipitation rates. (b) Vegetation cover data are from the SPOT/PROBA-V and PROBA-V/S3-OLCI products of the land service of Copernicus (Baret et al., 2013; Camacho et al., 2013). (e) Precipitation data are from the Global Precipitation Climatology Project (GPCP; Adler et al., 2016).