Multi vegetation model evaluation of the Green Sahara climate regime
<![CDATA[By Peter O. Hopcroft, Paul J. Valdes, Anna B. Harper and David J. Beerling
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.
During at least the late Quaternary (over the past 1 Myr) the area of the present-day Sahara desert has periodically transformed to an area of lakes, rivers, and steppe and savannah vegetation. This is shown by fossil pollen [Prentice et al., 2000], paleolake levels [Kohfeld and Harrison, 2000], paleohydrological reconstructions [Pachur and Kropelin, 1987; Pachur and Holzmann, 1991; Hoelzmann et al., 1998, 2001; Skonieczny et al., 2015], and archeological evidence [e.g., Dunne et al., 2012; Manning and Timpson, 2014].
Generally, the increased rainfall required to stimulate a Green Sahara is taken to equal the transition threshold for desert to steppe biome types or around 200–300 mm yr−1 [Joussaume et al., 1999; Braconnot et al., 2007]. This is based on a modern climate and atmospheric CO2 (350–400 ppm). Given that CO2 was lower at 264 ppm [Elsig et al., 2009], this would have slightly decreased the water use efficiency of vegetation and so this threshold value would have been larger during the early to middle Holocene. Pollen-based precipitation reconstructions suggest that the mid-Holocene Sahara was wetter, with 354 ± 36 mm yr−1 (mean ± standard error) additional rainfall [Peyron et al., 2006; Bartlein et al., 2011]. Rainfall reconstructions derived from leaf wax deuterium (δD) recovered from marine sediment cores [Tierney et al., 2017] show a much greater rainfall change compared with the pollen data, with absolute values of up to 1280 mm yr−1 at 31°N. The rainfall increase therefore extends farther northward in this reconstruction. Wu et al.  used an inverse modeling approach with a single equilibrium vegetation model (BIOME4) to estimate Sahara precipitation to be 100–400 mm yr−1 during the mid-Holocene. Thus, existing estimates span a considerable uncertainty range.
Although the underlying mechanism behind the “Green Sahara” must derive from orbitally induced amplification of monsoon rainfall, multiple amplifying feedbacks are likely required to produce the regime as implied by reconstructions. These may include changes in surface vegetation and soil properties [Kutzbach et al., 1996; Claussen and Gayler, 1997; Bonfils et al., 2001; Levis et al., 2004; Vamborg et al., 2011], which can act to reinforce the monsoon anomaly following the mechanism outlined by Charney , enhanced local moisture recycling because of expanded early to middle Holocene coverage of wetlands and lakes [Coe and Bonan, 1997; Krinner et al., 2012], and changes in the sea surface temperatures [Braconnot et al., 1999; Zhao et al., 2005]. Expanded forest cover in Eurasia relative to today may also have contributed to a northward migration of the intertropical convergence zone during the mid-Holocene [Swann et al., 2014].
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.
2 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 CO2. 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 CO2 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].
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 1961–1990. 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.
2.1 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 CO2 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.
3.1 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.
3.2 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.
3.3 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.
4.1 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 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. , the use of globally averaged regression coefficients [Guiot et al., 2000] may not be suitable for the climate regime of mid-Holocene north Africa.
4.2 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.  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.  showed that the dust change has a smaller impact on rainfall than the vegetation-induced surface albedo decrease alone.
4.3 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].
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 [Joussaume et al., 1999; 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.  and Figures 3j–3l for the reconstruction of Hély et al. . 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 CO2 [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.
We thank the CRU at the University of East Anglia for the observational climate data and the Met Office for providing the JULES code on the Met Office Science Repository. We thank D. Gerten for providing the LPJ v2.1 model code and Gethin Williams, Steph Cornford, and Sarah Shannon for computational support. D.J.B. gratefully acknowledges the award of an ERC advanced grant (CDREG, 322998). Simulations were performed using the facilities of the Advanced Computing Research Centre at the University of Bristol http://www.bris.ac.uk/acrc/. P.O.H. received funding from UK NERC projects NE/I010912/1 and NE/P002536/1. A.B.H. received funding from NERC grant NE/K016016/1. CRU climate forcing data are available from http://www.cru.uea.ac.uk/data; CRUNCEP forcing data are maintained at http://dods.extra.cea.fr/data/p529viov/cruncep/. The Hoelzmann et al. reconstruction is available from ftp://ftp.ncdc.noaa.gov/pub/data/paleo/pollen/africa6k/. The full DGVM outputs are available from http://www.paleo.bris.ac.uk/∼ggpoh/midhol_vegMIP. Selected variables are also included in the supporting information.