Climatic windows for human migration out of Africa in the past 300,000 years

Late Quaternary climate reconstructions

Precipitation

Our reconstructions of Late Quaternary precipitation are based on outputs from a statistical emulator of the HadCM3 general circulation model⁶². The emulator was developed using 72 3.75° × 2.5° resolution snapshot climate simulations of HadCM3, covering the last 120k years and in 2k year time steps from 120k to 22k years ago and 1k-year time steps from 21k years ago to the present, where each time slice represents climatic conditions averaged across a 30-year post-spin-up period^28,63. The emulator is based on grid-cell-specific linear regressions between the local time series of HadCM3 climate data and four time-dependent forcings, given by the mean global atmospheric CO₂ concentration and three orbital parameters: eccentricity, obliquity, and precession. The values of these four predictors are known well beyond the last 120k years; thus, applying them to the calibrated grid-cell-specific linear regressions allows for the statistical extrapolation of global climate up to 800k years into the past⁶². The emulated climate data have been shown to correspond closely to the original HadCM3 simulations for the last 120k years, and to match long-term empirical climate reconstructions well⁶².

Here, we used precipitation data from the emulator, denoted ({bar{{{{bf{P}}}}}}_{{{{{rm{HadCM3}}}}}_{{{{rm{em}}}}}}(t)), of the last 300k years at 1k-year time steps, (tin {{{{bf{T}}}}}_{300k}). The data were spatially downscaled from their native 3.75° × 2.5° grid resolution, and subsequently bias-corrected, in two steps, similar to the approach described in ref. ⁶⁴, whose description we follow here. Both steps use variations of the delta method⁶⁵, under which a high-resolution, bias-corrected reconstruction of precipitation at sometime t is obtained by applying the difference between lower-resolution present-day simulated and high-resolution present-day observed climate—the correction term—to the simulated climate at time t. The delta method has been used to downscale and bias-correct palaeoclimate simulations before (e.g. for the WorldClim database⁶⁶), and, despite its conceptual simplicity, has been shown to outperform alternative methods commonly used for downscaling and bias-correction⁶⁷.

A key limitation of the delta method is that it assumes the present-day correction term to be representative of past correction terms. This assumption is substantially relaxed in the dynamic delta method used in the first step of our approach to downscale ({bar{{{{bf{P}}}}}}_{{{{rm{HadCM}}}}{3}_{{{{rm{em}}}}}}(t)) to a ~1° resolution. This method involves the use of a set of high-resolution climate simulations that were run for a smaller but climatically diverse subset of T_300k. Simulations at this resolution are computationally very expensive, and therefore running substantially larger sets of simulations is not feasible; however, these selected data can be very effectively used to generate a suitable time-dependent correction term for each (tin {{{{bf{T}}}}}_{300k}). In this way, we can increase the resolution of the original climate simulations by a factor of ∼9, while simultaneously allowing for the temporal variability of the correction term. In the following, we describe the approach in detail.

We used high-resolution precipitation simulations from the HadAM3H model⁶³, generated for the last 21,000 years in 9 snapshots (2k year time intervals from 12k to 6k years ago, and 3k year time intervals otherwise) at a 1.25° × .83° grid resolution, denoted ({bar{{{{bf{P}}}}}}_{{{{rm{HadAM}}}}3{{{rm{H}}}}}(t)), where (tin {{{{bf{T}}}}}_{21k},)represents the nine time slices for which simulations are available. These data were used to downscale ({bar{{{{bf{P}}}}}}_{{{{rm{HadCM}}}}{3}_{{{{rm{em}}}}}}(t)) to a 1.25° × 0.83° resolution by means of the multiplicative dynamic delta method, yielding

The ⊞-notation indicates that the coarser-resolution data were interpolated to the grid of the higher-resolution data, for which we used an Akima cubic Hermite interpolant⁶⁸, which, unlike the bilinear interpolation, is continuously differentiable but, unlike the bicubic interpolation, avoids overshoots. The time (hat{t}in {{{{bf{T}}}}}_{21k}) is chosen as the time at which climate was, in a sense specified below, close to that at time (tin {{{{bf{T}}}}}_{300k}). In contrast to the classical delta method (for which (hat{t}=0) for all (t)), this approach does not assume that the resolution correction term, ((frac{{bar{{{{bf{P}}}}}}_{{{{rm{HadAM}}}}3{{{rm{H}}}}}(hat{t})}{{bar{{{{bf{P}}}}}}_{{{{rm{HadCM}}}}{3}_{{{{rm{em}}}}}}^{ boxplus }(hat{t})})), is constant over time. Instead, the finescale heterogeneities that are applied to the coarser-resolution ({bar{{{{bf{P}}}}}}_{{{{rm{HadCM}}}}{3}_{{{{rm{em}}}}}}(t)) are chosen from the wide range of patterns simulated for the last 21k years. The strength of the approach lies in the fact that the last 21k years account for a substantial portion of the glacial-interglacial range of climatic conditions present during the whole Late Quaternary. Following ref. ⁶⁴, we used global CO₂, a key indicator of the global climatic state, as the metric according to which (hat{t}) is chosen; i.e. among the times for which HadAM3H simulations are available, (hat{t}) is the time at which global CO₂ was closest to the respective value at the time of interest, t.

In the second step of our approach, we used the classical multiplicative delta method to bias-correct and further downscale ({{{{bf{P}}}}}_{ sim 1^circ }(t)) to a hexagonal grid⁶⁹ with an internode spacing of ~55 km (~0.5°),

where ({{{{bf{P}}}}}_{{{{rm{obs}}}}}(0)) denotes present-era (1960–1990) observed precipitation⁷⁰.

We reconstructed land configurations for the last 300k years using present-day elevation⁷¹ and a time series of Red Sea sea level⁷². For locations that are currently below sea level, the delta method does not work. For these locations, precipitation was extrapolated using a inverse distance weighting approach. With the exception of a brief window from 124–126k years ago, sea level in the past was lower than it is today; thus, present-day coastal patterns are spatially extended as coastlines move, but not removed. For all (tin {{{{bf{T}}}}}_{300k}), maps of annual precipitation ({bar{{{{bf{P}}}}}}_{ sim 0.5^circ }(t)) with the appropriate land configuration are available as Supplementary Movie 1.

Based on these data representing 30-year climatological normals at 1k-year time steps between 300k years ago and the present, we generated, for each millennium, 100 maps representing 10-year average climatologies as follows. We used 3.75° × 2.5° climate simulations from the HadCM3B-M2.1 model, providing a 1000-years-long annual time series of annual precipitation for each millennium between 21k years ago and the present⁷³. Millennia were simulated in parallel; thus, the 1000-years-long time series representing each millennium is in itself continuous, but the beginnings and ends of the time series of successive millennia generally do not coincide. For (tin {{{{bf{T}}}}}_{21k}), we denote the available 1000 successive maps of annual precipitation by ({{{{bf{P}}}}}_{{{{rm{HM}}}}}^{(1)}(t),ldots ,{{{{bf{P}}}}}_{{{{rm{HM}}}}}^{(1000)}(t)). We used these data to compute the relative deviation of the climatic average of each decade within a given millennium, and the climatic average of the 30-year period containing the specific decade as

Finally, we applied these ratios of 10-year to 30-year climatic averages to the previously derived 1k-year time step climatologies to obtain, for each (tin ,{{{{bf{T}}}}}_{300k}), 100 sets of 10-year average annual precipitation,

where, analogous to our approach in Eq. (1),(, boxplus ) denotes the interpolation to the ~55 km hexagonal grid, and where (hat{t}) is chosen as the time at which global CO₂ was closest to the respective value at time t.

Aridity

The Köppen aridity index used here is defined as the ratio of annual precipitation (in mm) to the sum of mean annual temperature (in °C) and a constant of 33 °C (cf. Eq. (8)). This measure of aridity was found to be the most reliable one of a set of alternative indices in palaeoclimate contexts³⁰.

Decadal-scale mean annual temperature data between 300k years ago and the present were created using analogous methods to those previously applied to reconstruct precipitation. 3.75° × 2.5° resolution emulator-derived simulations of mean annual temperature of the past 300k years at 1k time steps⁶², denoted ({bar{{{{bf{T}}}}}}_{{{{rm{HadCM}}}}{3}_{{{{rm{em}}}}}}(t)), were first downscaled by means of the additive dynamic delta method, using 1.25° × 0.83° HadAM3H simulations of mean annual temperature of the past 21k years, denoted ({bar{{{{bf{T}}}}}}_{{{{rm{HadAM}}}}3{{{rm{H}}}}}(t)), yielding, analogous to Eq. (1),

Analogous to Eq. (2), Next, present-day observed mean annual temperature, ({bar{{{{bf{T}}}}}}_{{{{rm{obs}}}}}(0)), was used to further downscale and bias-correct the data by means of the additive delta method to obtain

For all (tin {{{{bf{T}}}}}_{300k}), maps of mean annual temperature ({bar{{{{bf{T}}}}}}_{ sim 0.5^circ }(t)) with the appropriate land configuration are available as Supplementary Movie 1.

Finally, we incorporated HadCM3B-M2 simulations of mean annual temperature of the past 21k years, ({{{{bf{T}}}}}_{{{{rm{HM}}}}}^{(1)}(t),ldots ,{{{{bf{T}}}}}_{{{{rm{HM}}}}}^{(1000)}(t)) for (tin {{{{bf{T}}}}}_{21k}), to obtain 10-year average mean annual temperature,

Based on these data, the Köppen aridity index at the same spatial and temporal resolution is calculated as

Comparison with empirical proxies

Long-term proxy records

Long-term proxy records allow us to assess whether simulations capture key qualitative dynamics observed in the empirical data. The lack of direct long-term time series reconstructions of annual precipitation and mean annual temperature makes it necessary to use proxies related to these two climate variables. Proxies providing temporal coverage beyond the last glacial maximum are not only extremely sparse in North Africa and Southwest Asia, but even the few records that exist are affected by environmental factors other than the specific climate variables considered here. For example, reconstructions of past wetness and aridity use proxies that reflect not only rainfall conditions but also the interaction of precipitation with other local and non-local hydro-climatic variables, e.g. river discharge or hydrological catchment across a larger area. Here, we have not attempted to correct for such processes, but assumed that the simulated climate at the site where the empirical record was taken provide a suitable approximation of the potentially broader climatic conditions relevant for the proxy data. Realistic climate simulations would therefore be expected to match major qualitative trends of the empirical records, rather than exhibit a perfect correlation with the data. We compared our precipitation simulations against three long-term humidity-related empirical proxies (Fig. 4a). Proxy 1⁷⁴ provides a time series of Dead Sea lake levels, for which wet and dry periods are associated with high-stand and low-stand conditions, respectively. Proxy 2¹⁹ from the southern tip of the Arabian Peninsula was obtained from a marine sediment core that allows for reconstructing past changes in aridity over land from the stable hydrogen isotopic composition of leaf waxes (δD_wax). Proxy 3¹⁸ is an XRF-derived humidity index from a core near the Northwest African coast. Temperature simulations were compared against two long-term records of δ¹⁸O, which varies over time as a result of temperature fluctuations (in addition to other factors), from the Peqiin and Soreq caves in Northern Israel⁷⁵ (Fig. 4e). Overall, the simulated data capture key phases observed in the empirical records well for both precipitation (Fig. 4b–d) and temperature proxies (Fig. 4f–h).

Pollen-based reconstructions

Pollen records used to empirically reconstruct past climate do not reach as far back in time as the above-described proxy records and are not available at the same temporal resolution; however, in contrast to those proxies, they can be used to quantitatively estimate local annual precipitation and mean annual temperature directly. Here, we used the dataset of pollen-based reconstructions of precipitation and temperature for the mid-Holocene (6k years ago) and the last glacial maximum (21k years ago)⁷⁶ (Fig. 5a). Our precipitation and temperature data are overall in good agreement with the empirical reconstructions (Fig. 5b–e). During the mid-Holocene, our simulations suggest slightly less precipitation at low levels than most of the empirical records (Fig. 5d), while our data match the empirical reconstruction available from a very arid location during the last glacial maximum very well (Fig. 5e).

Interglacial palaeolakes on the Arabian Peninsula

Finally, we plotted time series of our precipitation simulations in three locations in which palaeolakes have been dated to the last interglacial period, following the approach in ref. ²⁴, in which the authors tested whether their climate simulations predicted higher rainfall during the last interglacial period than at present at palaeolake sites on the Arabian Peninsula. Figure 6 shows the locations of three palaeolakes in the northeast (western Nefud near Taymal; proxy 1), the centre (at Khujaymah; proxy 2), and the southwest (at Saiwan; proxy 3) of the peninsula²⁴ (described in detail in refs. ^23,77), and our precipitation data in these locations. In two out of the three locations, our data predict that more rainfall occurred at the estimated timings of the palaeolakes than at any point in time since; in the third location, slightly more rainfall than during the dated time interval is simulated only for a period around 8k years ago.

Determining the minimum precipitation and aridity tolerance required for out-of-Africa exits

We denote by ({{{bf{X}}}}={({lambda }_{1},{phi }_{1}),({lambda }_{2},{phi }_{2}),ldots },)the set of longitude and latitude coordinates of the hexagonal grid with an internode spacing of ~55 km (~0.5°)⁶⁹ that are contained in the longitude window [15°E, 70°E] and the latitude window [5°N, 43°N] (shown in Fig. 3). We denote by ({{{bf{E}}}}) the set of the present-day elevation values of the coordinates in ({{{bf{X}}}}) (in meters)⁷⁸, i.e. ({{{bf{E}}}}({x}_{i})) is a positive number in a point ({x}_{i}=({lambda }_{i},{phi }_{i})) if ({x}_{i}) is currently above sea level, and negative if ({x}_{i}) is currently below sea level. We denote by (s(t)) the sea level (in meters) at the time (tin {{{{bf{T}}}}}_{300k}) (where ({{{{bf{T}}}}}_{300k}) represents the last 300k years in 1k time steps), for which we used a long-term reconstruction of Red Sea sea level⁷². In particular, we have (s(0)=0) at present day. For each millennium (tin {{{{bf{T}}}}}_{300k}), we denote by (bar{{{{bf{X}}}}}(t)) the subset of points in (X) that are above sea level:

Based on the precipitation map ({{{{bf{P}}}}}_{ sim 0.5^circ }^{(d)}(t)) for a decade (d=1,ldots ,100) in millennium (t) (Eq. (4)), and a given precipitation threshold value (p) (in mm year⁻¹), we denote by ({mathop{{{{bf{X}}}}}limits^{=}}_{p}^{(d)}(t)) the subset of (bar{{{{bf{X}}}}}(t)) that would be suitable grid cells for humans assuming that they cannot survive in areas where precipitation levels are below (p):

We then determined whether there was a connected path in ({mathop{{{{bf{X}}}}}limits^{=}}_{p}^{(d)}(t)) between an initial point, for which we used ({x}_{{{{rm{start}}}}}=(32.6^circ {{{rm{E}}}},10.2^circ {{{rm{N}}}})), and any point in a set of coordinates outside of Africa, defined as ({{{{bf{X}}}}}_{{{{rm{end}}}}}mathop{=}limits^{{{{rm{def}}}}}{(lambda ,phi )in {{{bf{X}}}}:lambda, > , 65^circ {{{rm{E}}}},{{{rm{or}}}},phi , > , 37^circ N}). This was defined to be the case if there was a finite sequence

of points ({x}_{i}in {mathop{{{{bf{X}}}}}limits^{=}}_{p}^{(d)}(t)) such that the distance between any two successive points ({x}_{i}) and ({x}_{i+1}) was less or equal to the maximum internode spacing of the grid (X). Based on this approach, the critical precipitation threshold below which no connected path exists for the precipitation map ({{{{bf{P}}}}}_{ sim 0.5^circ }^{(d)}(t)) was determined using the following bisection method. Beginning with ({hat{p}}_{0}=1000) mm y⁻¹ and ({check{p}}_{0}=0) mm y⁻¹, for which a connected path between ({x}_{{{{rm{start}}}}}) and ({{{{bf{X}}}}}_{{{{rm{end}}}}}) exists, respectively, for all and for no (t) and (d), the values ({hat{p}}_{k}) and ({check{p}}_{k}) were iteratively defined as

For all (k), the sought critical precipitation threshold, denoted ({p}_{{{{rm{crit}}}}}^{(d)}(t)), is bounded above by ({hat{p}}_{k}) and bounded below by ({check{p}}_{k}). For (kto infty ), both values converge to ({p}_{{{{rm{crit}}}}}^{(d)}(t)). Here, we defined

which lies within 1 mm y⁻¹ of the true limit value.

To specifically determine the precipitation tolerance required for a northern (Fig. 1a) or southern (Fig. 1b) exit, we rendered the passage of the respective other route impassable by removing appropriate cells from the grid. When investigating the southern route, we additionally assumed that no sea level and precipitation constraints applied within a ~40 km radius around the centre of the Bab al-Mandab strait.

For aridity, the procedure is identical, with the exception that ({mathop{{{{bf{X}}}}}limits^{=}}_{p}^{(d)}(t)) is defined based on the relevant aridity map, ({{{{bf{A}}}}}_{ sim 0.5^circ }^{(d)}(t)), and the value 4.0 is used for the initial upper threshold (denoted ({hat{p}}_{0}) above).

Width of the Strait of Bab al-Mandab

Similar to ref. ⁵², we reconstructed the minimum distance required to cover on water in order to reach the Arabian peninsula (present-day west coast of Yemen) from Africa (present-day Djibouti and southeast Eritrea). We used a 0.0083° (~1 km at the equator) map of elevation and bathymetry⁷⁸ and a time series of Red Sea sea level⁷² to reconstruct very-high-resolution land masks for the last 300k years. For each point in time, we determined the set of connected land masses, and the distances between the closest points of any two land masses. The result can be graph-theoretically represented by a complete graph whose nodes represent connected land masses and whose edge weights correspond to the minimum distances between land masses. The path involving the minimum continuous distance on water was then determined by solving the minmax path problem whose solution is the path between the two nodes representing Africa and the Arabian Peninsula that minimises the maximum weight of any of its edges (Fig. 1b grey shades).

Analyses were conducted using Matlab R2019a⁷⁹.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com