Early Holocene greening of the Sahara requires Mediterranean winter rainfall

Significance Explaining the greening of the Sahara during the Holocene has been a challenge for decades. A strengthening of the African monsoon caused by increased summer insolation is usually cited to explain why the Sahara was vegetated from 14,000 to 5,000 y ago. Here, we provide a unique climate record of quantified winter, spring, and summer precipitation in Morocco over the past 18,500 y, and numeric simulations, which show that moisture contributions from the Mediterranean Sea and the North Atlantic Ocean in winter, were as important as the expanded summer monsoon for the greening of the Sahara during the African humid period. The findings of this study will help to better understand and simulate climate variability over northern Africa.

: Mean monthly precipitation (mm) of five main Mediterranean plant taxa identified in the Tislit pollen record. The plot shows a marked seasonal distribution of the precipitation throughout the year with clear wet winters and dry summers. Figure S4: January and July reconstructed temperature changes from the Tislit pollen record. The red marks on the Y axes indicate the modern values in the Imilchil closest weather station. The shaded blue area corresponds to the African Humid period.
Leaf wax analyses δ 13 C analyses of n-alkanes were conducted on a ThermoFisher Scientific MAT 252 isotope ratio mass spectrometer that was coupled, via a GC-C combustion interface with a nickel catalyzer operating at 1000 °C, to a ThermoFisher Scientific Trace GC equipped with a HP-5ms column (30 m, 0.25 mm, 0.25 m).
μ Each sample was measured in duplicate if sufficient material was available. δ 13 C values were calibrated against a CO2 reference gas of known isotopic composition, and are given in ‰ VPDB. Accuracy and precision were determined by measuring n-alkane standards, calibrated against the A4-Mix isotope standard (A. Schimmelmann, University of Indiana) every six measurements. The difference between the long-term means and the measured standard values yielded a 1 error of σ <0.3 ‰. The accuracy and precision of the squalane internal standard were both 0.1 ‰. The precision of the replicate analyses of the n-C29 and n-C31 alkanes was 0.2 and 0.1 ‰ on average, respectively (<0.1 -0.7 ‰ for n-C29 and <0.1 -0.6 ‰ for n-C31 alkane). For samples which could only be analyzed once, the long-term precision of the standards (0.3 ‰) was assumed. As δ 13 C values of the n-C29 and n-C31 alkanes co-vary strongly (r=0.81), we calculated a weighted average δ 13 C value for δ 13 Cwax using their relative abundance. δ 13 Cwax varies from -33.3 to -24.5 ‰ VPDB, with an average propagated precision of 0.2 ‰ (<0.1 -0.7). Mostly, δ 13 Cwax showed depleted 13 C values, indicative of C3 plant vegetation, with the exception of the period corresponding to Heinrich event 1 (HE1, 16.1 -14.9 ka BP), during which there were enriched 13 C values ( figure S6). This is consistent with the natural vegetation in Morocco, which is entirely of the C3 type in modern conditions. The elevated δ 13 Cwax values during HE1 coincide with an abundance of Poaceae and could thus indicate an expansion of C4 grasses around Lake Tislit under the extreme climatic conditions of HE1.
δD analyses of n-alkanes were conducted on a Thermo Fisher Scientific MAT 253™ Isotope Ratio Mass Spectrometer that was coupled, via a GC IsoLink operating at 1420°C, to a Thermo Fisher Scientific TRACE GC equipped with a HP-5ms column (30 m, 0.25 mm, 1 m). Each sample was μ measured in duplicate if sufficient material was available. δD values were calibrated against an H2 reference gas of known isotopic composition and are given in ‰ VSMOW. Accuracy and precision were controlled by the laboratory's internal n-alkane standard, calibrated against the A4-Mix isotope standard every six measurements, and by daily determination of the H3 + factor. Measurement precision was determined by calculating the difference between the analyzed values of each standard measurement and the long-term mean of standard measurements, which yielded a 1 error of σ <3 ‰. H3 + factors varied between 5.0 and 5.2. The accuracy and precision of the squalane internal standard were 4 ‰ and 3 ‰, respectively. The precision of the replicate analyses of the n-C29 and n-C31 alkanes was 1 ‰ on average (<1 -5 ‰). For samples which could only be analyzed once, the long-term precision of the standards (3 ‰) was assumed. As D values of the δ n-C29 and n-C31 alkanes co-vary (r=0.78), we calculated a weighted average D value for D δ δ wax, using their relative abundance. D δ wax varies between -162 to -193 ‰ VSMOW, with an average propagated precision of 2 ‰ (<1 -5 ‰) (figure S6). As different vegetation types can cause offsets in D δ wax, due to different hydrogen isotope fractionation factors, we estimated D δ precip using δ 13 Cwax as an estimate for vegetation type, according to established procedures 5 . For comparability, we used δ 13 Cwax endmembers and fractionation factors reported by Tierney 6 and based on Garcin 7 and Sachse 8 . D δ precip estimates range from -54 to -87 ‰ VSMOW (figure 5 in SI), which are consistently more depleted than the D δ precip estimates for GC27 6 . This is due to the high altitude setting of Lake Tislit and reflects the input of plant waxes from local sources. Additionally, the D δ precip estimates were compensated for global ice volume changes affecting isotopes in the hydrological cycle. For this purpose, Bintanja's marine oxygen isotope compilation 9 was used, with the LGM-present δ 18 O range scaled to 1 ‰ VSMOW. This ice-volume correction has an insignificant effect for the Holocene and leads to slightly depleted D δ precip estimates for the period over 10,000 years ago (figure S6). In addition, we compensated the D δ precip estimates for temperature changes, as changes in condensation temperature lead to changes in D δ precip 76 . In the absence of a specific local temperature lapse rate for Morocco, we used the global average of 5.6 ‰/°C for D δ precip 10 . We used the averaged winter temperature (Dec -Feb) from the pollen reconstructions. The essential effect of the temperature compensation is more enriched D δ precip estimates when winter temperatures were lower, i.e. prior to the Holocene (figure S6). Figure S6: Stable isotope compositions of plant wax in Tilsit core. a) δ 13 C compositions of the n-C29 (blue) and n-C31 alkane (orange) and weighted average δ 13 Cwax (green), b) δD compositions of the n-C29 (blue) and n-C31 alkane (orange) and weighted average δDwax (green), c) δDprecip estimates (see main text) and effects of ice volume correction (δDprecip-ic) and temperature correction (δDprecip-ic-tc). Data in c) are shown without error estimates to compare effects of corrections.

Vegetation model simulations
CARAIB is a mechanistic model that describes stomatal regulation and photosynthesis, growth and respiration, competition for resources, and mortality of a set of plant objects, which can be plant functional types (PFT), bioclimatic affinity groups or plant species. The daily model inputs are mean air temperature T, diurnal thermal amplitude Tmax -Tmin, precipitation P, percentage of sunshine hours SH, air relative humidity RH, and wind speed W. The main model outputs are soil water amount and water fluxes, as well as the abundance, gross and net primary productivity (GPP and NPP), biomass and leaf area index of all PFTs. From these outputs, a biome map can be constructed in a post-processing subroutine [11][12][13] . PFTs have been used in accordance with Henrot's classification 14 , which includes 26 PFTs designed for large-scale or global applications. These PFTs are distributed between two vegetation storeys: grasses and shrubs (PFTs 1-11) in the understorey, and trees (PFTs 12-26) in the overstorey. The model also calculates a budget of water fluxes to evaluate the amount of soil water in the root zone. Soil water deficit can induce: (1) stomatal closure and reduction of photosynthetic assimilation; (2) reduction of the leaf area index (leaf desiccation); and, under most severe water stress, (3) mortality of the PFTs. Natural fires can also induce mortality during dry periods. The model time step is 1 day for updating soil water or carbon reservoirs, and 2 hours for the calculation of photosynthesis and respiration fluxes.
In simulation 1 (figure 3.1), we selected the HadCM3 general circulation model 15 , as it performed an almost transient simulation over the last 21,000 years (with 1,000-year snapshots), with atmosphere and ocean grids of 3.75° x 2.5° and 1.25° x 1.25°, respectively, over 20 levels each. For the 9 ka experiments of HadCM3, only monthly mean data were available. The monthly anomalies were then calculated as the absolute anomaly (i.e., the difference between the 9 ka simulation and the pre-industrial control simulation) for temperature, and as a weighted mean of the absolute and relative (i.e., the same difference between the 9 ka simulation and the control, multiplied by the ratio between the observed and the control simulation values) anomaly for precipitation. The weights were chosen in such a way that they favored the absolute anomaly whenever possible, except when negative precipitation tended to be produced for 9 ka. In the latter situation, the weight of the relative anomaly was progressively increased 16 . The monthly anomalies were then simply added to the 1901-1930 monthly mean values, yielding a 30-year sequence of monthly temperature and precipitation for 9 ka. CARAIB's weather generator was then used to produce the daily values necessary to force the simulations for each year, drawing on the monthly sequence. This weather generator is based on the currently observed day-to-day variability of temperature (minimum and maximum) and precipitation per (geo-)climatic zone. A renormalization process is performed to ensure that the monthly values are not altered 17 . Note that the day-to-day variability reconstructed by the generator differs from that of the GSWP3 1901-1930 signal. However, this should not significantly impact the results, since both variabilities are based on meteorological data collected in a similar climatic zone (hence with similar distribution functions), and a relatively long record (30 years) is generated.
For consistency, in simulations 2 and 3 ( figure 3.2 and 3.3), the air relative humidity r was increased and the percentage of sunshine hours s was decreased at the same time as the precipitation was increased. These changes in r and s are critical for maintaining a Green Sahara, since they significantly reduce evapotranspiration and thus lead to wetter soils. They were estimated on the basis of statistical relationships established in the simulated areas for the 1901-1930 period in our climate dataset.
Using the monthly climatological means for 1901-1930, a significant linear correlation was indeed obtained in a log-log plot of relative humidity r (or sunshine hours s) versus precipitation that was increased by 1 mm/month (P+1). The Pearson correlation coefficient was 0.772 for r and -0.759 for s, when all months and all pixels of the simulated domain were included. This yielded the following relationships, allowing us to calculate r and s as a function of precipitation P: where r and s are expressed in % and limited to the indicated ranges, precipitation P is in mm/month, and αr = 23,0233, βr = 0,2153, αs = 88,8467, βs = -0,0936193. These relationships allowed us to calculate the anomalies: to be added to each value of r and s in the 1901-1930 time series, to obtain a 30-year time series of r and s for simulations 2 and 3.   Tables  Table S1. Identified pollen taxa of trees, shrubs and herbs in the Tislit record (68 taxa) with those used for the pollen-based climate reconstruction (45 in bold character).