Multi vegetation model evaluation of the Green Sahara climate regime

During the Quaternary, the Sahara desert was periodically colonized by vegetation, likely because of orbitally induced rainfall increases. However, the estimated hydrological change is not reproduced in climate model simulations, undermining confidence in projections of future rainfall. We evaluated the relationship between the qualitative information on past vegetation coverage and climate for the mid‐Holocene using three different dynamic vegetation models. Compared with two available vegetation reconstructions, the models require 500–800 mm of rainfall over 20°–25°N, which is significantly larger than inferred from pollen but largely in agreement with more recent leaf wax biomarker reconstructions. The magnitude of the response also suggests that required rainfall regime of the early to middle Holocene is far from being correctly represented in general circulation models. However, intermodel differences related to moisture stress parameterizations, biases in simulated present‐day vegetation, and uncertainties about paleosoil distributions introduce uncertainties, and these are also relevant to Earth system model simulations of African humid periods.

In coupled atmosphere-ocean general circulation model simulations of the mid-Holocene (6 kyr B.P. during the last Green Sahara episode) the North African monsoon migrates a few degrees northward, but models submitted to the Paleoclimate Model Intercomparison Project (PMIP) phases 2 and 3 fail to reproduce a convincing Green Sahara as inferred from the paleoclimate reconstructions [Braconnot et al., 2007;Perez-Sanz et al., 2014]. Part of this is undoubtedly caused by the idealized and unrealistic nature of the model boundary conditions employed: modern observed soils and vegetation either are prescribed where they are not interactive or are initialized from modern-day conditions. When some models are forced with extensive vegetation [Claussen and Gayler, 1997;Renssen et al., 2006], or an idealized increase in dust loading [Pausata et al., 2016], there is a better agreement with proxy data. However, including either dynamic vegetation or an interactive dust cycle (e.g., MIROC-ESM and HadGEM2-ES in Coupled Model Intercomparison Project phase 5 (CMIP5)) or using the stronger early-Holocene orbital forcing [e.g., Marzin and Braconnot, 2009;Singarayer and Valdes, 2010] does not appear to lead to substantial positive feedbacks in the region.
More generally, it is regarded that the Green Sahara represents a long-standing problem for Earth system models [Harrison et al., 2015] which undermines confidence in projections of monsoon change in future states [Bony et al., 2015]. This calls for a greater understanding of the hydrological regime and the ability of models to reproduce the coupled ecophysical transition.
To our knowledge no study has systematically evaluated the climate change that is required to stimulate a transition to a Green Sahara using dynamic vegetation models (DGVMs). We present ensembles of simulations with three different DGVMs. This study therefore provides a missing link between pollen-based climate reconstructions and general circulation model (GCM) simulations of African humid periods. It also addresses whether model-based inferences are compatible with the range of rainfall derived from fossil pollen or sedimentary leaf wax. Our simulations contribute to understanding the representation of biophysical feedbacks in the Sahel region [Zeng et al., 1999] and will provide a basis for improved understanding of vegetation feedbacks in existing and upcoming (CMIP6/PMIP4) mid-Holocene and Last Interglacial Earth system model simulations.

Dynamic Vegetation Model Experiments
Three different DGVMs were used: the Joint UK Land Environment Simulator (JULES) version 4.1 [Clark et al., 2011;Best et al., 2011], the Sheffield Dynamic Global Vegetation Model (SDGVM) [Woodward et al., 1995;Beerling and Woodward, 2001], and the Lund-Potsdam-Jena dynamic vegetation model (LPJ version 2.1) [Sitch et al., 2003;Gerten et al., 2004]; see Table S1 in the supporting information. These models calculate the fluxes of carbon and water at the land surface using a discrete set of plant functional types (PFTs) to represent the global distribution and structure of vegetation. Modeled vegetation distribution, structure, and productivity are dependent on the imposed soil texture, climate fields, and atmospheric CO 2 . The main differences in terms of process representations and a comparison of the modeled response to future climate scenarios have been documented before [Sitch et al., 2008]. JULES 4.1 employs nine rather than five plant functional types (PFTs) and updated relationships between leaf nitrogen and photosynthesis, as well as the introduction of a trade-off between leaf lifespan and leaf mass per unit area [Harper et al., 2016]. JULES has been extensively evaluated with observations [Blyth et al., 2010;Harper et al., 2016] and includes improvements derived from simulations of the Last Glacial Maximum [Hopcroft and Valdes, 2015a]. SDGVM incorporates dynamic vegetation and coupled carbon and nitrogen cycles, with six PFTs [Woodward et al., 1995;Beerling and Woodward, 2001], and has shown a reasonable representation of global vegetation dynamics [Beerling et al., 1997;Beerling and Woodward, 2001]. LPJ uses 10 PFTs and has been extensively evaluated against observations of global terrestrial productivity and vegetation and soil moisture [Sitch et al., 2003;Wagner et al., 2003;Gerten et al., 2004;Hickler et al., 2005]. Each dynamic vegetation model was run at a resolution of 1.875 ∘ × 1.25 ∘ in longitude-latitude, with an atmospheric CO 2 concentration of 280 ppmv for the preindustrial and with no anthropogenic land use. Soil types are prescribed based on the default data sets used within JULES and LPJ. In SDGVM desert soils had unrealistically low wilting points and so these values are read from a global soil database [Global Soil Data Task Group, 2000]. The soil data used in each model are described in more detail in the supporting information [Dharssi et al., 2009;Global Soil Data Task Group, 2000;Best et al., 2011;Gerten et al., 2004;Woodward et al., 1995;Zobler, 1986;Poulter et al., 2015].

10.1002/2017GL073740
The three models require different driving climate fields (see Table S2). For the control "preindustrial" simulations, JULES is forced with CRU TS 3.22 [Harris et al., 2014] and NCEP reanalysis [Kalnay et al., 1996] 6-hourly fields for . LPJ and SDGVM were forced with monthly fields from CRU TS 3.22 for 1961-1990 and spun up for 1000 and 500 years, respectively.

Mid-Holocene Climate Drivers
Monthly simulated climate fields for the mid-Holocene and preindustrial from PMIP3/CMIP5 coupled general circulation model simulations with CCSM4 [Gent et al., 2011], GISS-E2-R [Schmidt et al., 2014], IPSL-CM5A-LR [Kageyama et al., 2013], and MPI-ESM-P [Giorgetta et al., 2013], shown Figure S1, were used to construct multimodel mean anomalies for each required driving variable. See Table S3 for a comparison of the resolution and climate sensitivity values (as calculated by Hopcroft and Valdes [2015b]) of these models. The climate fields are shown in Figures S2-S4 for North Africa. With the exception of wet days (only required by LPJ), the anomalies were combined with the observational data sets from CRUTS3.22, to create mid-Holocene climate fields. Wetdays fields were subsequently calculated from the resultant monthly precipitation fields using a logarithmic relationship between the two variables calibrated against the observed relationship over the present-day north African monsoon region.
Mid-Holocene simulations are forced with an atmospheric CO 2 concentration of 264 ppm, and since LPJ and SDGVM are not forced with radiation fluxes, insolation is modified according to the orbital parameters at 6 kyr B.P. for these models. JULES is driven with surface downward short-wave radiation taken from the GCMs which already implicitly includes this astronomical forcing. For JULES the monthly fields from the mid-Holocene GCM simulations were interpolated to the required 6-hourly time step. Anomalies were combined in a relative sense for radiation fluxes, to avoid biasing the diurnal cycle.
Nine sensitivity experiments were conducted by perturbing the location or magnitude of the mid-Holocene climate anomalies imposed on the DGVMs. To alter the latitudinal extent of the monsoon anomaly, the mid-Holocene anomaly (GCM mid-Holocene minus GCM preindustrial) was artificially extended by +5 and +10 ∘ N. To perturb the magnitude, the anomalies in each variable (temperature, humidity, radiation fluxes, or cloud cover) are multiplied by an amplification factor and the multimodel regression between rainfall and that variable. These anomalies are then combined with the observed climate fields to produce the mid-Holocene driving variables. In the northern part of the domain, where monsoon changes are minimal, the temperature, relative humidity, and surface short-wave radiation responses are a direct result of the increased Northern Hemisphere summer insolation. These variables are therefore only perturbed by the amplification factor over grid cells with a concurrent significant precipitation change of greater than 0.5 mm d −1 , otherwise the climate changes in these regions are assumed to be caused directly by orbitally induced changes in solar insolation. The monthly mean precipitation, temperature, and cloud cover anomalies over North Africa in each simulation are shown in supporting information Figures S2-S4.

Vegetation Regimes for Modern Conditions
We compared the simulated preindustrial fractional vegetation coverage with the European Space Agency Climate Change initiative (ESA-CCI) vegetation coverage [Poulter et al., 2015] in Figure 1. Over Western Africa the observed desert-steppe transition (taken to coincide with 20% total vegetation coverage) occurs at 281 mm yr −1 of rainfall. The simulated transition values are 378, 515, and 182 mm yr −1 in JULES, LPJ, and SDGVM, respectively. JULES and SDGVM are equidistant from the observations, while LPJ predicts much less vegetation coverage in low-rainfall areas. A similar picture emerges for total fractional vegetation coverage in the Sahara, which is 14%, 10% and 20% in the JULES, LPJ and SDGVM, respectively, compared to the observed coverage of 23% [Poulter et al., 2015]. Reconstructed potential natural vegetation coverage which is derived from a combination of historical cropland inventory data with observed and modeled vegetation coverage [Ramankutty and Foley, 1999] shows a desert-steppe boundary that is located approximately 2 ∘ farther south than the satellite observations but with a similar overall vegetation coverage. Thus, land use change (which is not included in these simulations) does not play a major role in the model-observation discrepancies.
The relationship between total (C3 + C4) grass coverage and precipitation (Figure 2) also shows that coverage is confined to high-precipitation regions in JULES and LPJ, while SDGVM has too much vegetation coverage in arid regions. For the grass fraction the 20% coverage transition rainfall values are 263 mm yr −1 in the observations and 442, 579, and 182 mm yr −1 in the three models. Thus, SDGVM is closer to the observations than JULES and LPJ which both predict too little vegetation coverage in low-rainfall areas. The relationship between vegetation coverage and precipitation ( Figure 2) shows considerable variation in the observations. The peak grass coverage is not well replicated in any model, though SDGVM looks qualitatively the most similar. Experiments with the latest version of SDGVM (which includes PFT age classes but not dynamic vegetation) [Woodward and Lomas, 2004] show similar productivity fields as the version of SDGVM used here, and so we expect that the results presented are robust to the inclusion of age class distributions. A trend toward tree-dominated landscapes at higher precipitation is seen in each model, with a peak tree coverage at 1600, 1100, and 1600 mm rainfall in JULES, LPJ, and SDGVM, respectively, compared with a peak of 1600 mm in the observations.

Vegetation Coverage in the Mid-Holocene Simulations
The dominant PFT in each grid cell was calculated and converted to a common set of five PFTs for intermodel comparison (broadleaf tree, needleleaf tree, C3 and C4 grasses, and shrubs; see Table S4). The dominant  PFT where the bare soil fractional coverage is less than 80% is shown in Figure 3 for the control preindustrial simulations and selected mid-Holocene sensitivity tests. The full set of nine simulations are shown in Figures S5-S7.
For the standard mid-Holocene simulations no model simulates a significant northward expansion of vegetation. This confirms previous work which concluded that the rainfall anomalies simulated by the GCMs are far too small to account for the vegetation change reconstructed for the early and mid-Holocene [Perez-Sanz et al., 2014;Harrison et al., 2015]. In fact, only when the applied rainfall is significantly amplified compared to the GCM predictions or shifted considerably farther northward do the models produce tangible land cover change across the domain.
The reconstructed vegetation types for the mid-Holocene Sahara are associated with some uncertainty. Hoelzmann et al. [1998] used pollen samples and charcoal and inferred a sharp transition from near 100% savanna to steppe at around 15 ∘ N. Hély et al. [2014] used a more recent and extensive set of pollen samples covering the Holocene and inferred humid forest (Guineo-Congolian), wooded savanna (Sudanian), and grassland (Sahelian) types as far as 20 ∘ N, 23 ∘ N, and 26 ∘ N, respectively, for 6 kyr B.P. We converted these results to the grid used in our simulations and assumed dominant PFTs based on observations of current biome properties. For Hély et al. [2014] we assumed latitudinal homogeneity, which is qualitatively consistent with independently derived hydrological reconstructions [Larrasoana et al., 2013].
The adapted reconstructions are compared with each simulation, and the best fitting of the nine mid-Holocene sensitivity simulations in each model is shown with the reconstructions in Figure 3. The three models turn out to be relatively consistent, so that the most extreme scenario considered provides the best model-data agreement for all but two cases, those being for JULES and SDGVM in comparison with the Hoelzmann et al. [1998] reconstruction. In this latter case a smaller rainfall anomaly (2xMH+shift10) provides a superior model-data comparison. The simulated vegetation cover is very different between the models. For example, Figure S5 shows grass expansion across the Sahara in SDGVM that is not evident in the other models, and Figure S6 shows more forest expansion in LPJ than in JULES or SDGVM. The different precipitation value at which peak forest coverage occurs (1100 mm yr −1 in LPJ versus 1600 mm yr −1 in the other models) as well as differences in key processes, which we evaluate further in the following section, contribute to this.
The results in Figures 3f-3h and 3j-3l show that the zonally averaged rainfall in the northern part of the domain (from 20 ∘ to 25 ∘ N) may have exceeded 800 mm yr −1 (see Figure S2), compared with just 20 mm yr −1 today and only 48 mm yr −1 in the PMIP3 GCM simulations analyzed here. This is consistent with a compilation of rainfall reconstructions from leaf wax hydrogen isotopes [Tierney et al., 2017]. However, this magnitude of change represents a more serious challenge for GCM simulations of this time period, compared to the pollen-based reconstructions [Bartlein et al., 2011], which suggest a zonally homogeneous rainfall anomaly of only 350 mm yr −1 .

Underlying Causes of the Intermodel Differences
The reasons for intermodel differences in vegetation distributions and carbon and water fluxes are complex. DGVMS are known to exhibit a range of sensitivities [Cramer et al., 2001] which are caused by varying importance of several interconnected processes. While a full analysis is beyond the scope of this study, moisture limitation is likely to be a dominant factor in this region. The moisture stress formulations are described in the supporting information. JULES and SDGVM both parameterize the moisture stress as a function of soil moisture, while LPJ uses the ratio of potential evapotranspiration and atmospheric water demand.
A moisture availability factor ( ) was calculated for each model as detailed in the supporting information. The annual mean values for each grid cell over north Africa (0-35 ∘ N, 20 ∘ W-40 ∘ E) are shown in Figure 4 for the preindustrial simulations. JULES has a near-linear relationship between soil moisture and , whereas in LPJ, shows little correlation with soil moisture, consistent with radiation being a more important controlling variable [Medlyn et al., 2016]. In SDGVM the factor is quasi-linear with soil moisture. In the relatively vegetated regions north of the equator (0-15 ∘ N) the models show reasonable agreement. North of this, SDGVM has almost universally higher values of compared with the other two models, while LPJ has extremely low values, mostly less than 0.2. This divergent behavior of moisture availability helps explain the differences in the control and the changes in the mid-Holocene simulations. The diversity in moisture stress parameterizations appears to introduce much of the intermodel uncertainty in the simulations. This uncertainty in modeled moisture stress could be addressed with field manipulation experiments [e.g., Medlyn et al., 2016], but this is beyond the scope of this work.

Comparison With Past Studies
Our systematic study of a Green Sahara provides the first model-based assessment of the required climatic conditions. These simulations also give a context for understanding dynamic vegetation responses when coupled within climate-vegetation model simulations [e.g., Rachmayani et al., 2015] and provides impetus to better understand differences between dynamic vegetation schemes. Pollen-based reconstructions for this 10.1002/2017GL073740 time period [Peyron et al., 2006;Bartlein et al., 2011] show a substantially smaller precipitation change than inferred in two of the models. The pollen reconstructions also show less rainfall than recent implied by the pollen taxa [Hély et al., 2014] or leaf wax biomarkers [Tierney et al., 2017], which both imply larger rainfall increases of at least 1000 mm yr −1 as far as 20 ∘ N.
These discrepancies could imply that some of the reconstruction methods are biased, perhaps because there is no modern-day analogue for a Green Sahara and because changes in humidity and incident radiation are not explicitly accounted for in the reconstruction methodology. For example, while a multiple regression of temperature and precipitation with sunlight was employed in the inverse vegetation modeling study of Wu et al. [2007], the use of globally averaged regression coefficients [Guiot et al., 2000] may not be suitable for the climate regime of mid-Holocene north Africa.

Potential Atmosphere-Vegetation Feedbacks
The Sahel region is thought to be a hot spot for atmosphere-land surface coupling [The GLACE Team, 2004]. This could be argued to also apply to the whole Sahel/Sahara region during the mid-Holocene. A "greening" of the Sahara will reduce the land surface albedo, thereby enhancing the land-ocean thermal contrast contributing a positive rainfall feedback [Charney, 1975;Kutzbach et al., 1996;Braconnot et al., 1999]. The DGVMs simulate an increase in leaf area index and vegetation coverage in response to increasing rainfall. Using standard formulae for land surface albedo as a function of leaf area index [Best et al., 2011], we calculate the surface albedo change for the near-vegetated state. Averaged over 10-25 ∘ N, the models predict albedo reductions of 0.05 for JULES and LPJ and 0.08 in SDGVM.
Increased vegetation coverage will also increase potential canopy water storage leading to an increase in local moisture recycling from plant transpiration [Braconnot et al., 1999], again contributing to enhanced rainfall. The proportion of transpiration relative to the total evapotranspiration is 59% and 46% for LPJ and SDGVM (a similar diagnostic is not available from this version of JULES), whereas isotopic evidence [Jasechko et al., 2013] and satellite-based calculations [Miralles et al., 2011] suggest higher values of 76% and 80-90%, respectively. Hence, the importance of vegetation feedbacks for moisture recycling are likely underestimated in current DGVMs and Earth system models.
Vegetation will also contribute to organic matter in the litter fall causing changes in soil properties with potential for further positive feedbacks from the resultant surface albedo [Kutzbach et al., 1996;Levis et al., 2004;Vamborg et al., 2011]. Other secondary effects include a reduction in dust cycle [deMenocal et al., 2000;Williams et al., 2016]. However, Pausata et al. [2016] showed that with an idealized increase in dust loading over the Sahara region, simulated rainfall in the West African monsoon is amplified. This only occurs when accompanied by significant surface albedo decrease from vegetation expansion. Pausata et al. [2016] showed that the dust change has a smaller impact on rainfall than the vegetation-induced surface albedo decrease alone.

Potential Avenues for Model Improvement
Improving the simulation of vegetation dynamics in the semiarid regions is important for better understanding the natural carbon cycle [Poulter et al., 2014;Ahlström et al., 2015], for predicting the coupled land surface atmosphere response to ongoing climate change in hot spots like the Sahel, and, as shown in this work, for understanding the large environmental changes in Sahel/Sahara region during the late Quaternary era.
Bare soil evaporation is likely overestimated in current land surface models, because the formulations used are only strictly valid over a few centimeters depth [Mahfouf and Noilhan, 1991;van den Hoof et al., 2013] and because the proportion of transpiration is generally underestimated [Lawrence et al., 2007]. The modeled response of moisture-stressed plants may be improved by moving from empirically based dependencies of stomatal conductance to a photosynthesis optimization approach [e.g., Bonan et al., 2014;Clark et al., 2015].

Conclusions
During the early to middle Holocene era (11-4 kyr B.P.) the region of the present-day Sahara desert transitioned to a vegetated state, with lakes and river systems that are not evident today [Hoelzmann et al., 1998;Skonieczny et al., 2015]. This transition was ultimately driven by changes in incident solar insolation caused by variations in Earth's orbit through time.
Although some GCMs configured with a vegetated Sahara show rainfall anomalies close to the pollen reconstructions [Claussen and Gayler, 1997;Renssen et al., 2006;Pausata et al., 2016], three generations of climate models Braconnot et al., 2007;Perez-Sanz et al., 2014] fail to reproduce this transition [Harrison et al., 2015], although the idealized nature of these simulations discussed earlier and shortcomings in snapshot versus fully transient integrations should also be considered. For these and other reasons, the climatic regime that caused the Green Sahara is not well constrained.
We performed an ensemble of simulations with three different DGVMs to address the relationship between the observed vegetation change and the associated climate regime. We find that none of these models satisfactorily reproduces the current dependence of vegetation in semiarid regions on climate; two models simulate grass coverage that is too low, while the third shows the opposite. In spite of this, all three DGVMs require significantly more rainfall to satisfy the constraints from two contrasting reconstructions of the mid-Holocene vegetation distribution, as shown in Figures 3f-3h for the reconstruction of Hoelzmann et al. [1998] and Figures 3j-3l for the reconstruction of Hély et al. [2014]. The inferred precipitation anomalies are closer to recent reconstructions from leaf wax biomarkers [Tierney et al., 2017] and significantly exceed inferences from fossil pollen [Peyron et al., 2006;Bartlein et al., 2011]. Thus, the model results tentatively support the conclusion that the early to middle Holocene was significantly wetter than previously assumed.
Much of the intermodel variability found in this study is traced to moisture stress in the semiarid/arid regions. The lack of a consensus on the best formulation for moisture limitation of photosynthesis and transpiration therefore needs to be tackled to reduce intermodel uncertainty in vegetation dynamics [Medlyn et al., 2016]. This also implies that state-of-the-art Earth system models incorporating coupled atmosphere-vegetation components are unlikely to produce realistic simulations of past African humid periods including the mid-Holocene and the Last Interglacial (e.g., as planned for CMIP6). Given the important role of semiarid regions in the variability of atmospheric CO 2 [Poulter et al., 2014;Ahlström et al., 2015], and in the case of the Sahel, for human impacts of climate change, this raises questions for other aspects of the Earth system modeling and associated uncertainties [Zeng et al., 1999;Wang et al., 2014].
The dynamic vegetation models examined here are similar to the land surface components included in Earth system models. The recent PMIP3 ensemble shows an increase in precipitation of around 30 mm yr −1 over the critical region of 20-30 ∘ N for the mid-Holocene [Perez-Sanz et al., 2014]. This is at least a factor of 20 to 25 smaller than required by all three dynamic vegetation models examined. This and the discrepancies between rainfall and vegetation reconstructions all warrant further study.