2.1 Experimental design

The Mid-Holocene (MH) time-slice climate experiment (6000 BP) represents the initial state for the last 6000-year transient simulation with dynamical vegetation. It is thus considered as a reference climate in this study. Because of this, and to save computing time, model adjustments made to set up the model content and the model configuration were mainly done running MH and not preindustrial (PI) simulations (Tables 1 and 2). Only a subset of PI simulations is available for comparison with modern conditions. All the simulations were run long enough (300–1000 years) to reach a radiative equilibrium and be representative of a stabilized MH climate (Fig. 1). They are free of any artificial long-term trends after the adjustment phase, as were IPSL PMIP3 MH simulations (Fig. 1; Kageyama et al., 2013a).

Table 1Tests done to set up the model IPSL version in which we included the dynamical vegetation. For all these simulations the vegetation map is prescribed to the 1860 CE map used in PI-PMIP3. The different columns highlight the name of the test and the initial state to better isolate the different factors contributing to the adjustment curves in Fig. 1. We include in parentheses the tag of the simulation that corresponds to our internal nomenclature for memory.

Download Print Version | Download XLSX

Table 2Simulations run to initialize the dynamical vegetation starting from bare soil or from vegetation map and soil moisture resulting from an off-line ORCHIDEE simulation with dynamical vegetation switched on and using the PI L11 simulated climate as boundary conditions. We include in parentheses the tag of the simulation that corresponds to our internal nomenclature for memory.

* Simulations with an asterisk are considered as references for the model version and the transient simulations.

Download Print Version | Download XLSX

Most tests follow the MH-PMIP3 protocol (Braconnot et al., 2012). This is only due to the fact that this work began before the PMIP4 boundary conditions were available. But the transient simulation (TRHOLV, for TRansient HOLocene simulation with dynamical Vegetation) and the 1000-year-long MH simulations with or without dynamical vegetation that were run to prepare the initial state for it follow the PMIP4-CMIP6 protocol (Otto-Bliesner et al., 2017, Table 1). In all simulations the Earth's orbital parameters are derived from Berger (1978). The MH-PMIP3 protocol uses the trace gases (CO₂, CH₄ and N₂O) reconstruction from ice core data by Joos and Spahni (2008). It has been updated for PMIP4, using new data and a revised chronology that provides a consistent history of the evolution of these gases across the Holocene (Otto-Bliesner et al., 2017). The difference in forcing between PMIP4 and PMIP3 was estimated to be −0.8 W m⁻² by Otto-Bliesner et al. (2017). This is the order of magnitude found for the imbalance in net surface heat flux at the beginning of the MH-FPMIP4 simulation. This simulation started from L11Aer run with the PMIP3 protocol (Fig. 1a). It uses the same model version but follows the PMIP4 protocol. For the subset of PI experiments Earth's orbit and trace gases are prescribed at their 1860 CE values, i.e., the beginning of the industrial era. For the MH and PI time slice experiments, boundary conditions do not vary with time. For the transient simulations the Earth's orbital parameters and trace gases are updated every year.

Figure 1Illustration of the effect of the different adjustments made to produce Mid-Holocene simulations with the modified version of the IPSLCM5A-MR model in which the land surface model ORCHIDEE includes a different soil hydrology and snow models (see text for details). The three panels show the global average of (a) net surface heat flux (W m⁻²), (b) evaporation (kg m⁻²), and (c) 2 m air temperature (^∘C). The different colored lines represent the results for the different simulations reported in Table 1.

Download

In standard versions of the IPSL model, aerosols are accounted for by prescribing the optical distribution of dust, sea salt, sulfate, and particulate organic matter (POM), so as to take into account the coupling between aerosols and radiation (Dufresne et al., 2013). For MH simulations these variables are prescribed at 1860 CE values, for which the level or sulfate and POM is slightly higher than the values found in the Holocene (Kageyama et al., 2013a). Here, except for the first few tests (Table 1), we prescribe only dust and sea salt at their 1860 values and neglect the other aerosols. A fully coupled dust–sea-salt–climate version of the model that does not consider the other aerosols is under development for long transient simulations. For future comparisons it is important to have a similar model setup. Indeed, compared to the version with all aerosols, considering only dust and sea salts imposes a radiative difference of about 2.5 W m⁻² in external climate forcing. Its footprint appears on the net heat flux imbalance at the beginning of L11Aer. It leads to a global air temperature increase of 1.5 ^∘C (Fig. 1c). The largest warming over land is found in the Northern Hemisphere, but the ocean warms almost everywhere by about 1 ^∘C, except in the Antarctic circumpolar current (Fig. 2b). The warmer conditions favor higher precipitation with a global pattern rather similar to what is found in future climate projections (Fig. 2a). This offset affects the mean climate state and is larger than the expected effect of Holocene dusts.

Figure 2Mid-Holocene annual mean precipitation (mm d⁻¹) and 2 m air temperature (^∘C) difference between (a, b) L11Aer and L11, (c, d) L11 and PMIP3, (e, f) PMIP3L11AerEv and L11Aer, and (g, h) FPMIP4 and PMIP3. See Table 1 and text for details on the different simulations.

Download

2.2 The IPSL Earth system model and updated version of the land surface scheme

For these simulations, we use a modified version of the IPSLCM5A model (Dufresne et al., 2013). This model version couples the LMDZ atmospheric model with 144×142 grid points in latitude and longitude (2.5^∘ × 1.27^∘) and 39 vertical levels (Hourdin et al., 2013) to the ORCA2 ocean model at 2^∘ resolution (Madec, 2008). The ocean grid is such that resolution is enhanced around the Equator and in the Arctic due to the grid stretching and pole shifting. The LIM2 sea-ice model is embedded in the ocean model to represent sea-ice dynamics and thermodynamics (Fichefet and Maqueda, 1999). The ocean biogeochemical model PISCES is also coupled to the ocean physics and dynamics to represent the marine biochemistry and the carbon cycle (Aumont and Bopp, 2006). The atmosphere–surface turbulent fluxes are computed taking into account fractional land–sea area in each atmospheric model grid box. The sea fraction in each atmospheric grid box is imposed by the projection of the land–sea mask of the ocean model on the atmospheric grid, allowing for a perfect conservation of energy (Marti et al., 2010). Ocean, sea ice, and atmosphere are coupled once a day through the OASIS coupler (Valcke, 2006). All the simulations keep exactly the same set of adjusted parameters as in Dufresne et al. (2013) for the ocean–atmosphere system.

The land surface scheme is the ORCHIDEE model (Krinner et al., 2005). It is coupled to the atmosphere at each atmospheric model 30 min physical time steps and includes a river runoff scheme to route runoff to the river mouths or to coastal areas (d'Orgeval et al., 2008). Over the ice sheet, water is also routed to the ocean and distributed over wide areas so as to mimic iceberg melting and to close the water budget (Marti et al., 2010). This model accounts for a mosaic vegetation representation in each grid box, considering 13 (including two crops) plant functional types (PFTs) and interactive carbon cycle (Krinner et al., 2005).

We made several changes in the land surface model (Table 1). The first one concerns the inclusion of the 11-layer physically based hydrological scheme (de Rosnay et al., 2002) that replaces the 2-layer bucket-type hydrology (Ducoudré et al., 1993). The 11-layer hydrological model had never been tested in the fully coupled mode before this study. We paid particular attention to closing the water budget of the land surface model to ensure that O (1000 years) simulations will not exhibit spurious drift in sea level and salinity. In addition, the new prognostic snow model was included (Wang et al., 2013). The scheme describes snow with three layers that are distributed so that the diurnal cycle and the interaction between snowmelt and runoff are properly represented. In order to avoid snow accumulation on a few grid points, snow depth is not allowed to exceed 3 m. The excess snow is melted and added to soil moisture and runoff while conserving water and energy (Sylvie Charbit and Christophe Dumas, personal communication, 2017). Because of a large cold bias at high latitudes in the first tests, we reduced the bare soil albedo that is used to combine fresh snow and vegetation in the snow-aging parameterization. Other changes concern the adjustments of some of the parameterizations. The way the mosaic vegetation is constructed in ORCHIDEE favors bare soil too much when leaf area index (LAI) is low (Guimberteau et al., 2018). To overcome this problem, an artificial factor of 0.70 was implemented in front of bare soil evaporation to reduce it (Table 1). This factor is compatible with the order of magnitude of the reduction brought by the implementation of a new evaporation parametrization for bare soil in the current IPSLCM6A version of the model (Philippe Peylin, personal communication, 2017). For all the other surface types the evaporation is computed as in L11. The last adjustment concerns the combination of snow albedo with the vegetation albedo. The procedure was different when vegetation was interactive or prescribed. Now the combination of snow and vegetation albedo is based on the effective vegetation cover in the grid box in both cases. It leads to a larger albedo than with the IPSL-CM5A-LR reference version when vegetation is prescribed. It counteracts the effect of the fresh snow albedo reduction.

2.3 Impact of the different changes on model climatology and performances

Figures 1 and 2 highlight how the changes discussed in Sect. 2.2 affect the model adjustment and climatology. The hydrological model (L11) produces about 1.25 mm d⁻¹ higher global annual mean evaporative rates than MH-PMIP3. The water cycle is more active in L11. Precipitation is enhanced in the midlatitudes and over the tropical lands (Fig. 2c) where larger evapotranspiration and cloud cover both contribute to cooling the land surface (Fig. 2d). A higher evaporative rate should lead to a colder global mean temperature (Fig. 1c). This is not the case. The large-scale cooling over land is compensated for by warming over the ocean (Fig. 2d), caused by reduced ocean evaporation and changes in the ocean–land heat transport. The radiative equilibrium is achieved at the top of the atmosphere with the same global mean long-wave and shortwave radiation budget in the two simulations (L11 and MH_PMIP3).

The factor introduced to reduce bare soil evaporation did not lead to the expected reduction in evaporation (Fig. 1b). Indeed, when evaporation is reduced, soil temperature increases and the regional climate becomes warmer, allowing for more moisture in the atmosphere and thereby more evaporation where soil can supply water (Figs. 2e, f and 1). Therefore, differences resulting from bare soil evaporation do not show up on the precipitation map (Fig. 2e) but in the increased temperature over land in the Northern Hemisphere (Fig. 1f). This is consistent with similar findings when analyzing land use feedback (Boisier et al., 2012). Note that in Fig. 2f about 0.1 ^∘C of the 0.4 ^∘C global warming in L11AerEv is still a footprint of the warming induced by the aerosol effect described in Sect. 2.1 but that it does not alter our conclusions on the regional temperature–evaporation feedback.

The difference between MH-FPMIP4 and MH-PMIP3 represents the sum of all the changes in the land surface model and forcing discussed above (Fig. 2g, h). A PI simulation performed with the new model version (PI-FPMIP4, Table 1) allows us to assess how they affect the model performances. A rapid overview of model performances is provided by a simple set of metrics derived from the PCMDI Metric Package (Gleckler et al., 2016, see Appendix 1). Figure 3 highlights that temperature biases are reduced in PI-FPMIP4 at almost all model levels but that biases are enhanced for precipitation and total precipitable water compared to PI-PMIP3 (comparison of blue and black lines in Fig. 3). Taken together all the changes we made have little effect on the bias pattern (Fig. 3a). The model performs quite well compared to the CMIP5 ensemble of PI simulations, except for cloud radiative effect (Fig. 3). The effect of cloud in the IPSLCM5A-LR simulations results mainly from low level clouds over the ocean (Vial et al., 2013; Braconnot and Kageyama, 2015). The atmospheric tuning is exactly the same as in the default IPSLCM5A-LR version, and the introduction of all the changes described above has almost no effect on the cloud radiative effect. Overall, the model version with the 11-layer hydrology has similar skill to the IPSLCM5A reference (Dufresne et al., 2013), and we are confident that the version is sufficiently realistic to serve as a basis on which we can include the dynamical vegetation.

Figure 3(a) Spatiotemporal root mean square differences (rms_xyt) and (b) annual mean global model bias (bias_xy) computed on the annual cycle (12 climatological months) over the globe for the different preindustrial simulations considered in this paper (colored lines) and individual simulations of the CMIP5 multi-model ensemble (gray lines). The metrics for the different variables are presented as parallel coordinates, each of them having their own vertical axis with corresponding values. In these plots, “ta” stands for temperature (^∘C) with “s” for surface (850 and 200 for 850 and 200 hPa, respectively), “prw” for total water content (g kg⁻¹), “pr” for precipitation (mm d⁻¹), “rlut” for outgoing long-wave radiation (W m⁻²), and “rltcre” and “rltcre” for the cloud radiative effect at the top of the atmosphere in the shortwave and long-wave radiation, respectively (W m⁻²). See Appendix A for details on the metrics.

Download

2.4 Initialization of the Mid-Holocene dynamical vegetation for the transient simulation

We added the vegetation dynamics by switching on the dynamical vegetation model described in Zhu et al. (2015). Compared to the original scheme (Krinner et al., 2005), it produces more realistic vegetation distribution in mid- and high-latitude regions when compared with present-day observations.

Two different strategies have been tested to initialize the dynamical vegetation (Table 2). In the first case (MH-Vmap), the initial vegetation distribution was obtained from an off-line simulation with the land surface model forced by the CRU-NCEP 1901–1910 climatology. In the second case (MH-Vnone), the model restarted from bare soil with the dynamical vegetation switched on, using the same initial state as MH-Vmap for the atmosphere, the ocean, the sea ice, and the land ice. As expected, the evolution of bare soil, grass, and trees is very different in MH-Vmap and MH-Vnone during the first adjustment phase (black and blue curves in Fig. 4a–c). Vegetation adjusts in less than 100 years (1200 months) in MH-Vmap (blue curve). This short-term adjustment indicates that the climate–vegetation feedback has a limited impact on vegetation when the initial state is already consistent with the characteristics of the simulated climate. In MH-Vnone, which starts from bare soil (black curve), the adjustment has a first rapid phase of 50 years for bare soil and about 100 years for grass and trees, followed by a longer phase of about 200 years. The latter corresponds to a long-term oscillation that has been induced by the initial coupling shock between climate and land surface. Note that PMIP4 instead of PMIP3 MH boundary conditions were used to run the last part of these simulations (red and yellow curves in Fig. 4a–c). In the coupled system, most of the vegetation adjustment takes about 300 years, which is longer than results of off-line ORCHIDEE simulations (less than 200 years). Since MH-Vnone started from a coupled ocean–atmosphere–ice state at equilibrium, this result also indicates that the land–sea–atmosphere interactions do not alter the global energetics of the IPSL model much in this simulation where atmospheric CO₂ is prescribed. The two simulations converge to very similar global vegetation cover. Figure 4 suggests that there is only one global mean stable state for the Mid-Holocene with the IPSL model, irrespective of the initial vegetation distribution (see also Table A1, Appendix A1).

Figure 4Long-term adjustment of vegetation for the Mid-Holocene (MH)(a, b, c) and preindustrial (PI) climate(c, d, e), when starting from bare soil (Vnone) or from a vegetation map (Vmap). The 13 ORCHIDEE PFTs have been gathered as bare soil, grass, trees, and land use. When the dynamical vegetation is active, only natural vegetation is considered. Land use is thus only present in one simulation, corresponding to a preindustrial map used as a reference in the IPSL model (Dufresne et al., 2013). The corresponding vegetation is referred to as PI_prescribed. The x axis is in months, starting from 0, which allows us to plot all the simulations that have their own internal calendar on the same axis.

Download

For the transient simulation, we decided to use the results of the MH-VNone simulation as initial state (Table 2) because it leads to more realistic forest in the PI-Vnone simulation (see discussion in Sect. 4). We performed a preindustrial simulation (PI-Vnone) using MH-Vnone as an initial state and switching on the orbital parameters and trace gases to their PI values. Figure 3 indicates that the vegetation feedback slightly degrades the global performances for PI temperature and brings the model performance close to the IPSLCM5A-LR CMIP5 version. It also contributes to reducing the mean bias in precipitable water, evaporation, precipitation, and long-wave radiation, but it has no effect on the bias pattern (assessed by the rmst (spatiotemporal root mean square error) in Fig. 3; see also Appendix). Vegetation thus has an impact on climate, but its effect is of smaller magnitude than those of the different model and forcing adjustments done to set up the model version we use here. Section 4 provides a more in-depth discussion of the vegetation state.

3.1 Long-term forcing

Starting from the MH-Vnone simulation, the transient simulation of the last 6000 years (TRHOLV) allows us to test the response of climate and vegetation to atmospheric trace gases and Earth's orbit (see Sect. 2.1). The atmospheric CO₂ concentration slowly rises throughout the Holocene from 264 ppm 6000 years ago to 280 for the preindustrial climate at around −100 BP (1850 CE) and then experiences a rapid increase from −100 BP to 0 BP (1950 CE) (Fig. 5). The methane curve shows a slight decrease and then follows the same evolution as CO₂, whereas NO₂ remains around 290 ppb throughout the period. The radiative forcing of these trace gases is small over most of the Holocene (Joos and Spahni, 2008). The largest changes occurred with the industrial revolution. The rapid increase in the last 100 years of the simulation has an imprint of about 1.28 W m⁻².

Figure 5Evolution of trace gases CO₂ (ppm), CH₄ (ppb), and N₂O (ppb) and seasonal amplitude (maximum annual – minimum annual monthly values) of the incoming solar radiation at the top of the atmosphere (W m⁻²) averaged over the Northern Hemisphere (black line) and the Southern Hemisphere (red line). These forcing factors correspond to the PMIP4 experimental design discussed by Otto-Bliesner et al. (2017).

Download

The major forcing is caused by the slow variations in the Earth's orbital parameters that induce a long-term evolution of the magnitude of the incoming solar radiation seasonal cycle at the top of the atmosphere (Fig. 5). It corresponds to decreasing seasonality (difference between the maximum and minimum monthly values for each year) in the Northern Hemisphere and increasing seasonality in the Southern Hemisphere (Fig. 5). It results from the combination of the changes in summer and winter insolation in both hemispheres (Fig. 6). These seasonal changes are larger at the beginning of the Holocene (about −8 W m⁻² per millennium in the Northern Hemisphere, NH, and +5 W m⁻² per millennium in the Southern Hemisphere, SH) and then the rate of change linearly decreases in the NH (increases in the SH) from 4500 to about 1000 BP (Fig. 5). There is almost no change in seasonality in the NH over the last 1000 years, whereas in the SH seasonality starts to decrease again by 2000 BP. The shape of insolation changes is thus different in both hemispheres and so is the relative magnitude of the seasonal cycle between the two hemispheres. This would be seen whatever the calendar we used to compute the monthly means because of the seasonal asymmetry induced by precession in the MH (see Joussaume and Braconnot, 1997; Otto-Bliesner et al., 2017).

Figure 6Long-term evolution of incoming solar radiation at the top of the atmosphere (TOA) (W m⁻²; upper section of panels) and associated response of temperature (^∘C) and precipitation (mm yr⁻¹) expressed as the difference to the 6000 BP initial state and smoothed by a 100-year running mean) for (a) NH summer, (b) NH winter, (c) SH summer, and (d) SH winter. Temperatures are plotted in red and precipitation in blue for summer, and they are, respectively, plotted in orange and green for winter. NH summer and SH winter correspond to June to September averages, whereas NH winter and SH summer correspond to December to March averages. All curves, except insolation, have been smoothed by a 100-year running mean.

Download

3.2 Long-term climatic trends

Changes in temperature and precipitation follow the long-term insolation changes in each hemisphere during summer until about 2000 to 1500 BP (Fig. 6). Then trace gases and insolation forcing become equivalent in magnitude and small compared to MH insolation forcing, until the last period where trace gases lead to a rapid warming. The NH summer cooling reaches about 0.8 ^∘C and is achieved in 4000 years. The last 100-year warming reaches 0.6 ^∘C and almost counteracts, for this hemisphere and season, the insolation cooling. SH summer (JJAS) and NH winter (NDJF) temperatures are both characterized by a first 2000-year warming. It reaches about 0.4 ^∘C. It is followed by a plateau of about 3000 years before the last rapid increase of about 0.6 ^∘C that reinforces the effect of the Holocene insolation forcing. During SH winter, temperature does not seem to be driven by the insolation forcing (Fig. 6d). In both hemispheres. summer precipitation trends correlate well with temperature trends, as is expected from a hemispheric first-order response driven by a Clausius–Clapeyron relationship (Held and Soden, 2006). This is not the case for winter conditions because one needs to take into account the changes in the large-scale circulation that redistribute heat and energy between regions and hemispheres (Braconnot et al., 1997; Saint-Lu et al., 2016).

Figure 7Fraction of the evolution of the seasonal amplitude of temperature (red) and precipitation (blue) represented by the projection of these climate variables on the evolution of the seasonal amplitude of insolation as a function of latitude. The solid line stands for the raw signal and the dotted line for the signal after a 100-year smoothing.

Download

We further estimate the links between the long-term climate response and the insolation forcing for the different latitudinal bands by projecting the zonal mean temperature and precipitation seasonal evolution on the seasonal evolution of insolation. We define the seasonal amplitude for each year as the difference between the maximum and minimum monthly values. We consider for each latitude the unit vector S:

where $\mathit{S}=\frac{\mathit{S}{\mathit{W}}_{\mathrm{is}-\mathrm{TOA}}}{\Vert \mathit{S}{\mathit{W}}_{\mathrm{is}-\mathrm{TOA}}\Vert }$ represents the norm of the seasonal magnitude of the incoming solar radiation at TOA over all time steps (SW_is−TOA(t); $t=-\mathrm{6000}$ years to 0, with an annual time step). Any climatic variable (V) can then be expressed as

with

where b is the unit vector orthogonal to S. The ratio $\mathit{\alpha }/(\mathit{\alpha }+\mathit{\beta }{)}^{\mathrm{1}/\mathrm{2}}$ measures in which proportion a signal projects onto the insolation (Fig. 7). Figure 7 confirms that the projection of temperature and precipitation onto the insolation curve is larger in the Northern than in the Southern Hemisphere. The best match is obtained between 10 and 40^∘ N where about 80 % of the temperature signal is a direct response to the insolation forcing. The projections are only 50 % in the tropics in the Southern Hemisphere. These numbers go up to 90 % if a 100-year smoothing is applied to temperature. Seasonal precipitation is also projected to be as high as 90 % when considering the filtered signal, confirming the strong links between temperature and precipitation in the NH over the long timescale. The projection is poorer, but not zero, when the raw precipitation signal is considered. At the Equator and at high latitudes in both hemispheres the projection is poor or zero. At the Equator, the MH insolation forcing favors a larger north–south seasonal march of the ITCZ (Intertropical Convergence Zone) over the ocean and the inland penetration of Afro-Asian monsoon precipitation during boreal summer. Surface temperature is reduced in regions where precipitation is enhanced due to the combination of increased cloud cover and increased surface evaporation (Joussaume et al., 1999; Braconnot et al., 2007a). When the monsoon retreats to its modern position, surface temperature in these regions increases, thereby enhancing the monsoon's seasonal cycle. It is thus out of phase compared to insolation forcing. This is also true over SH continents where temperatures in regions affected by monsoons do not follow the local insolation and have a similar seasonal evolution to the Northern Hemisphere. This out-of-phase relationship is consistent with glacier reconstructions (Jomelli et al., 2011). At a higher latitude the projection of the raw signal is not good because of the large decadal variability. North of 40^∘ N the mixed layer depth is also larger (about 200 m) than in the tropics (about 70 m), which contribute to damping the seasonal change over the ocean. Thus the seasonal temperature response is flatter than the shape of the seasonal insolation forcing, which leads to a poor projection over midlatitude and high-latitude ocean especially in the ocean-dominated SH (Fig. 7). Sea-ice cover has also little change north of 80^∘ N, which also damps the changes in seasonality (Fig. 8). These changes are, however, amplified by the increase in sea ice during summer in the Arctic resulting from cooler conditions with time and by the reduction in the winter sea-ice cover in the Labrador and the Greenland–Iceland–Nordic seas (Fig. 8a, b). For snow cover the conditions are contrasted depending on the regions (Fig. 8a, b, d), with an decrease in the maximum cover over Eurasia related to a long-term rise in minimum temperature (Fig. 8d)

Figure 8(a) Total change in snow cover (kg m⁻²) and sea-ice fraction (%) integrated over the last 6000 years, and evolution from the Mid-Holocene of annual mean maximum summer and minimum winter values for (b) sea-ice averages over the Northern Hemisphere, (c) snow (solid lines) and 2 m air temperature (dotted lines) averaged for all regions north of 60^∘ N, and (d) snow and 2 m air temperature over Eurasia. In (b, c), and (d) black, dark blue and light blue stand, respectively, for the annual mean, maximum, and minimum annual monthly values for sea-ice or snow cover, and black, green, and red stand for annual mean, minimum, and maximum air temperature.

Download

3.3 Long-term vegetation trends

These long-term climate evolutions have a counterpart on the long-term evolution of vegetation (Fig. 9). At the global and hemispheric scale, the long-term vegetation trends correspond to reductions or increases in the area covered by vegetation, reaching 2 % to 4 % of the total land area depending on vegetation type (Fig. 9). The global vegetation averages reflect the Northern Hemisphere changes where most of the vegetated continental masses are located. As expected from the different long-term trends in insolation, the long-term evolution of tree and grass cover is the opposite in the two hemispheres. Note, however, that bare soil slightly increases in both hemispheres.

Figure 9Long-term evolution of the simulated (a) bare soil, (b) grass, and (c) tree covers, expressed as the percentage (%) of global, NH, or SH continental areas, and (d) GPP (PgC yr⁻¹) over the same regions. Annual mean values are smoothed by a 100-year running mean.

Download

In the Northern Hemisphere the changes follow the changes in summer temperature, with the best match obtained for grass, which increases almost linearly until 2000 BP and then remains quite stable. In the Southern Hemisphere the phasing between vegetation change and temperature is not as good, again because this hemisphere is dominated by ocean conditions rather than land conditions. However, the tree expansion reaches a maximum between 2000 and 1000 BP, and then the tree cover slightly decreases, which corresponds to the slight cooling in SH summer temperature. The gross primary productivity (GPP; Fig. 9d) is driven in both hemispheres by the changes in tree cover. It accounts for a reduction of about 5 PgC yr⁻¹. It is, however, possible that the GPP change is misestimated in this simulation because CO₂ is prescribed in the atmosphere, which implies that the carbon cycle is not fully interactive. Figure 10 compares the vegetation map obtained for the preindustrial period in TRHOLV (150–100 BP) with MH vegetation. It shows that bare soil increases in semiarid regions in Africa and Asia as well as in southern Africa and Australia. The reduction in tree PFTs is maximal north of 60^∘ N, in south and southeast Asia, the Sahel, and most of North America. Tree PFTs are replaced by grass PFTs. In the Southern Hemisphere, forest cover increases in southern Africa, southeast South America, and part of Australia.

In the last 100 years the effect of trace gases and in particular the rapid increase in the atmospheric CO₂ concentration leads to a rapid vegetation change characterized by tree regrowth, which is dominant in the NH (Figs. 9 and 10). This tree recovery counteracts the reduction from the Mid-Holocene in mid- and high NH latitudes (Fig. 10b, e, and h). This effect is consistent with the observed historical growth in gross primary production discussed by Campbell et al. (2017). The GPP increase in the last 100 years results from increased atmospheric CO₂. It suggests that the CO₂ effect counteracts the tree decline induced by insolation. When reaching 0 BP (1950 CE), bare soil remains close to PI, grass is reduced by 3 %, and trees are increased by about 3 %.

Figure 10Vegetation map comparing (a, d, g) the Mid-Holocene (first 50 years) and the preindustrial (50 years around 1850 CE (150 to 100 BP)) periods of the transient simulation, (b, d, h) the differences between the historical period (last 50 years) and the preindustrial period of the transient simulation, and (c, f, i) the difference between preindustrial climate for the transient simulation and the PI-Vnone simulations. For simplicity we only consider bare soil (a–c), grass (d–f), and trees (g–i).

Download

3.4 Regional trends

Figure 11 highlights the long-term vegetation trends for three regions that, respectively, represent climate conditions north of 60^∘ N, over the Eurasian continent, and in the West African monsoon Sahel–Sahara region. These are regions for which there are large differences in MH–PI climate and vegetation cover (Figs. 9 and 10). They have also been chosen because they are widely discussed in the literature and are also considered as tipping points for future climate change (Lenton et al., 2008). They are well suited to provide an idea of different characteristics between regions.

Figure 11Long-term evolution of bare soil, grass, and trees, expressed as the percentage of land cover north of 60^∘ N, over Eurasia, and over West Africa. The different values are plotted as differences with the first 100-year averages. A 100-year running mean is applied to the curves before plotting.

Download

North of 60^∘ N a substantial reduction in trees at the expense of grass starts at 5000 BP (Fig. 11a, b). Vegetation almost reaches its preindustrial conditions around 2500 BP. The largest trends are found between 5000 and 2500 BP in this region, and this reflects the timing of the NH summer cooling well. The change in total forest in Eurasia is small. A first change is followed by a second one around 3000 BP. Despite the 100-year smoothing applied to all the curves, they exhibit large decadal to multi-centennial variability. Over West Africa (Fig. 11c), the largest trends start slightly later (4500–5000 BP) and are more gradual until 500 BP. The vegetation trends are also punctuated by several centennial events that do not alter the long-term evolution as much as some of these events do in the other two boxes.

The variability found for vegetation is also found in temperature and precipitation at the hemispheric scale (Fig. 6). It is even higher at the regional scale in midlatitudes and high latitudes (Fig. 8). This variability is not present in the imposed forcing. It results from internal noise. Because of this it is difficult, for example, to say whether the NH winter temperature trend was rapid until 4000 BP and then temperature remains stable or whether the event impacting temperature and precipitation around 4800 to 4500 BP masks a more gradual increase until 3000 BP, as is the case for NH summer where the magnitude of the temperature trend is larger than its variability (Fig. 6). Note that some of these internal fluctuations reach half of the total amplitude of the regional vegetation trends (Fig. 11) and that they are a dominant signal over Eurasia, where the long-term mean change in the total tree cover is small (Figs. 10 and 11). Temperature and precipitation are well correlated at this centennial timescale (Fig. 6).

Despite the dry bias over the Sahel region in this version of the model, the timing of the vegetation changes over West Africa reported in Fig. 11 is consistent with the major features discussed for the end of the African humid period (Liu et al., 2007; Hély et al., 2014). In particular, the replacement of grass by bare soil starts earlier than the reduction in the tree cover located further south (Fig. 11). At the scale of the Sahel region, we do not have abrupt changes in vegetation but gradual ones. It is, however, abrupt at the grid cell level. These changes are associated with the long-term decline in precipitation, as well as the southward shift of the tropical rain belt associated with the African monsoon (Fig. 12). The location (latitude) of the rain belt is estimated here as the location of the maximum summer precipitation zonally averaged between 10^∘ W and 20^∘ E over West Africa. Most of the southward shift of the rain belt occurs between MH and 3500 BP and corresponds to a difference of about 1.8^∘ N of latitude over this period. Then the southward shift is smaller, with a total shift of 2.5^∘ N of latitude diagnosed in this simulation. The comparison of Figs. 11 and 12 clearly shows that the rapid decrease in vegetation occurs after the rapid southward shift of the rain belt. An interesting point is that the amount of precipitation is also shifted in time compared to the location of the rain belt. It suggests that the vegetation feedback on precipitation is still effective during the first period of precipitation decline and that it might have amplified the reduction in precipitation when vegetation is reduced over the Sahel region.

As seen in Fig. 10, the NH decrease in forest cover is mainly driven by the changes that occur north of 60^∘ N (Figs. 8 and 10). These trends reflect more or less what is expected from observations (Bigelow et al., 2003; Jansen et al., 2007; Wanner et al., 2008). It results from the summer cooling that affects both the summer sea-ice cover in the Arctic, the summer snow cover over the adjacent continent, and the amplification of the insolation forcing south of 70^∘ N by snow–vegetation albedo feedback. Further south over Eurasia, Fig. 11 suggests that there are only marginal changes in Eurasia in terms of vegetation. Figure 13 shows that the total tree cover over this region does not reflect the mosaic vegetation and forest composition well. Indeed, the long-term decrease in forest is dominated by the decrease in temperate and boreal deciduous trees. Boreal needleleaf evergreen trees do not change, whereas the temperate ones increase. This figure also highlights that the long-term change in Eurasian tree composition throughout the Mid- to late Holocene is punctuated by centennial variability. The different trees also have different timing and variability. Boreal forests are more sensitive to variability during the first 3000 years of the simulation, whereas temperate broadleaf trees exhibit larger variability in the second half. The large events have a climatic counterpart (Fig. 8), so that the composition of the vegetation is a result of a combined response to the long-term climatic change and to variability. These two effects can lead to different vegetation composition depending on stable or unstable vegetation states (Scheffer et al., 2012). Decadal vegetation changes have been discussed for recent climate in these regions (Abis and Brovkin, 2017), which suggests that despite the fact that our dynamical vegetation model might underestimate vegetation resilience, the rapid changes in vegetation mosaic are a key signal over Eurasia. Future model–data comparisons should consider composition changes and variability to properly discuss vegetation changes over this region.

Figure 12Evolution of (a) the location of the West African monsoon annual mean (black) and maximum (red) rain belt in degrees of latitude and (b) annual mean (black), minimum (green), and maximum (red) monthly precipitation (mm d⁻¹) averages over the Sahel region. The first 100 years have been removed and a 100-year running mean applied before plotting.

Download

4.1 Simulated versus reconstructed vegetation

Section 3 shows how climate and vegetation respond to insolation and trace gases. The simulated changes are in broad agreement with what is expected from various sources of data. However, Sect. 2 mentions model adjustments and biases. They all contribute to the difficulty of producing the right vegetation changes in the right place, at the right time, and for the right reasons. It is thus important to fully understand what we can expect in terms of realism from this simulation. We investigate it for the Mid-Holocene and modern climate, for which we can use the BIOME6000 vegetation reconstruction (Harrison, 2017).

The dynamical vegetation module simulates fractional cover of 13 PFTs. These PFTs cannot be directly compared with the reconstructed biome types based on pollen and plant macrofossil data from the BIOME 6000 dataset (Harrison, 2017). In order to facilitate the comparison, we converted the simulated PFTs into eight mega-biomes, using the biomization method algorithm proposed by Prentice et al. (2011). The algorithm uses a mixture of simulated climate and vegetation characteristics (see Appendix and Fig. A2). Alternative thresholds as proposed in previous studies (Joos et al., 2004; Prentice et al., 2011) were tested to account for the uncertainties in the biomization method (see Fig. A2). At a first glance MH-Vnone reproduces the large-scale pattern found in the BIOME6000 reconstruction (Fig. 14a). The comparison, however, indicates that the boreal forest tree line is located too far south. It results from a combination of a cold bias in temperature in these regions and a systematic underestimation of forest biomass in Siberia with ORCHIDEE when forced by observed present-day climate (Guimberteau et al., 2018). Such underestimation of tree biomass could lead to too low a tree height in ORCHIDEE and thus to the replacement of boreal forest by dry woodland according to the biomization algorithm (Fig. A2a). Also, vegetation is underestimated in West Africa, consistent with a dry bias (not shown). The underestimation of the African monsoon precipitation is present in several simulations with the IPSL model (Braconnot and Kageyama, 2015) and is slightly enhanced in summer when the dynamical vegetation is active. With interactive vegetation, however, equatorial Africa is more humid (Fig. 15a). Figure 14c provides an idea of the major mismatches between simulated vegetation and the BIOME6000 reconstructions. In particular, the simulation produces too much desert where we should find grass and shrubs. It also produces too much tundra instead of boreal forest and too much savannah and dry woodland in several places that should be covered by temperate trees, boreal trees, or tundra, confirming the visual map comparison (Fig. 14c). Similar results are found when considering the preindustrial climate in TRHOLV compared to the BIOME6000 preindustrial biome reconstruction (Fig. A2d). These are systematic biases. These systematic biases are confirmed when comparing the simulated PFTs for PI with those of the 1860 CE map estimated from observations and used in simulations with prescribed vegetation (see Table A1 for regions without land use).

Figure 13Evolution of the different tree PFTs in Eurasia, expressed as the percentage change compared to their 6000 BP initial state. Each colored line stands for a different PFT. Values have been smoothed by a 100-year running mean.

Download

Figure 14(a) Simulated mega-biome distribution by MH-Vnone, converted from the modeled PFT properties using the default algorithm described in Fig. A1. Panels (b) and (c): reconstructions in BIOME 6000 DB version 1 for the MH and PI periods (Harrison, 2017). (d) Number of pixels where reconstruction is available and the model matches (or does not match) the data. Note that multiple reconstruction sites may be located in the same model grid cell, in which case we did not group them so that each site was counted once. Numbers in parentheses on the x axis in (d) represent the number of sites for each biome type. Panel (e) same as (c) but for the number of matches between (e) the BIOME6000 MH (6 k) and PI (0 k) reconstructions at pollen sites and (f) the simulated mega-biomes for MH and PI at each model grid cell.

Download

It is not possible to estimate the vegetation feedback on the long-term climate evolution from the transient simulation. It is, however, possible to infer how the dynamical vegetation affects the mean climatology for the MH, a period for which simulations with prescribed and dynamic vegetation are available. Metrics discussed in Sect. 2.4 (Fig. 3) show that the introduction of the dynamical vegetation in the model reduces the amount of precipitation and that the climate is dryer. The simulations with dynamical vegetation only consider natural vegetation, whereas the 1860 CE map we prescribe when vegetation is fixed includes land use. In regions affected by land use, all MH simulations produce less bare soil (3 %), more tropical trees (5 %), similar temperate tree cover, increased boreal tree cover (10 %), and a different distribution between C3 versus C4 grass (see Table A1). In Eurasia where croplands are replaced by forest, the lower forest albedo induces warmer surface conditions (Fig. 15b). Also, when snow combines with forest instead of grasses, the snow–vegetation albedo is lower, leading to the positive snow–forest feedback widely discussed for the last glacial inception (de Noblet et al., 1996; Kutzbach et al., 1996). Figure 15a also highlights that precipitation is increased over the African tropical forest and reduced over South America. In most regions the impact of vegetation is much smaller than the impact of the changes in the land surface hydrology and forcing strategy discussed in Sect. 2.3 (Fig. 2).

Figure 15Impact of the dynamical vegetation and initialization of vegetation on the simulated climate. Differences for annual mean (a, c, e) precipitation (mm d⁻¹) and (b, d, f) 2 m air temperature (^∘C) between (a) the MH in the TRHOLV simulation and (b) the MH simulation without dynamical vegetation (MH FPMIP4), (c) the Mid-Holocene and (d) the preindustrial simulations in the TRHOLV simulation, and (e) and (f) the two preindustrial simulations initialized from bare soil (PI-Vnone) or a vegetation map for vegetation (PI-Vmap). See Table 2 and text for details on the simulations.

Download

The differences between the MH simulated vegetation map and the 1860 CE map reflect both systematic model biases and vegetation changes related to the MH climate differences with PI. We can infer from Figs. 15 and 16 that vegetation has a positive warming feedback in the high latitudes during MH. Part of the differences between the MH and the PI conditions in Fig. 15c, d is dominated by the impact of vegetation. Similar patterns as those obtained for the impact of vegetation are found over Eurasia for temperature, or southeast Asia and North America for precipitation. For the grid points where BIOM6000 data are available for both MH and PI (0 k – 1860 CE), the major simulated biome changes occur for savannah and wood and grass and shrubs (Fig. 14e). Differences are also found for trees and tundra, but to a lesser extent. The comparison with similar estimates from BIOME6000 reconstructions indicates that grass and shrubs exhibit major changes and that trees show larger differences compared to the simulation. The model shift between savannah and wood and grass and shrubs is consistent with the noted bias for savannah and the fact that the tree cover is underestimated in NH latitudes (Fig. 14).

Figure 16Difference between Vegetation maps obtained with the two different initial states for (a, c, e) Mid-Holocene simulations and (b, d, f) preindustrial simulations. Vmap stands for MH and PI simulations where the Mid-Holocene vegetation has been initialized from a vegetation map and Vnone for MH and PI simulations where the Mid-Holocene has been initialized from bare soil. For simplicity we only consider fractions of (a, b) bare soil, (c, d) grass, and (e, f) trees.

Download

Note that the vegetation differences found between the historical period and the PI period in TRHOLV are not negligible. We can estimate from Fig. 15a, b that neglecting land use leads to an underestimation of about 1 ^∘C in Eurasia between the MH and PI in this TRHOLV simulation. Depending on whether PI or the historical period is used as a reference, the magnitude of the MH changes in vegetation and climate are different. Also, land use has regional impacts and should be considered in PI or in the historical period. This stresses that a quantitative model–data comparison must be treated with caution, knowing that both the reference period (PI or historical) and the complexity of the land surface model (prescribed vegetation, natural dynamical vegetation, land use … ) can easily lead to 1 ^∘C difference in some regions.

4.2 Multiple vegetation states for the preindustrial climate

Another source of uncertainty concerns the stability of the simulated vegetation maps. Several studies suggest that the initial state only has a minor impact on the final climate because there is almost no change in the thermohaline circulation over this period and models do not exhibit major climate bifurcations (e.g., Bathiany et al., 2012). This is the main argument used by Singarayer and Valdes (2010) to justify that their suite of snapshot experiments may provide a reasonable transient climate vision when put together. Is this the case in the TRHOLV simulation when vegetation is fully interactive? This transient simulation does not exhibit much change in the indices of thermohaline circulation, which remains close to 16–18 Sv (sverdrup) (1 Sv =10⁶ m³ s⁻¹) throughout the period. The global metrics (Fig. 3) show that at the global scale the results of the TRHOLV simulations for PI (1860 CE) are similar to those of PI-Vnone. It is also the case for seasonal and extratropical or tropical values (Fig. A1). We can therefore conclude that there is no difference in mean surface climate characteristics between the snapshot PI-Vnone experiments and the PI period simulated in the transient TRHOLV simulation.

Is the vegetation then also similar to the one simulated in PI-VNone? The PI vegetation simulated in TRHOLV shows little differences to the one found for PI-Vnone (Fig. 10c, f, and i). The relative percentages of land covered by the different vegetation classes correspond to 15 % for bare soil, 41 % for grass, and 43 % for trees, respectively. These values are similar to the one found for PI-VNone (15 %, 40 %, and 44 %, respectively) within 1 % error bar. This does not necessarily hold at the regional scale where regional differences are also found between PI-TRHOLV and PI-Vnone. Indeed, Fig. 10 indicates differences in tree and grass cover in Eurasia around 60^∘ N and different geographical coverage between bare soil, grass, and trees over southern Africa and Australia. These differences are very small compared to the differences between MH and PI in TRHOLV, but are as large as the difference between historical (i.e., the last 50 years of the simulation) and PI in a few places in Eurasia. As seen in previous section, these are regions where variability is large and vegetation instable.

We also tested whether the PI vegetation and climate would also be similar when starting from MH-Vmap instead of MH-Vnone (dark pink and orange lines in Fig. 4d–f). This is also a way to have a better idea of the range of response one would expect from ensemble simulations, knowing that we only ran one full transient simulation. For the PI-Vmap simulation, the orbital parameters and trace gases were first prescribed to preindustrial conditions for 15 years while maintaining the vegetation PFTs in each grid cell to those obtained in MH-Vmap (Table 2, Fig. 4). Then, the dynamical vegetation was switched on. It induces a rapid transition of the major PFTs that takes about 10 years before a new global equilibrium is reached (Fig. 4d–f). For PI-VNone presented in Sect. 2.4 the same procedure was applied, but the dynamical vegetation was switched on after 5 years (Table 2 and Fig. 4), and the new equilibrium state is reached without any relaxation or rapid transition.

PI-Vnone and PI-Vmap converge to different global vegetation states (Fig. 4). Compared to the values listed above for PI-Vnone and PI-TRHOLV the respective covers of bare soil, grass, and trees for PI-Vmap are 20 %, 37 %, and 43 %. In particular, PI-Vmap produces a larger bare soil cover than PI-Vnone (Fig. 4d). It is even larger than the total bare soil cover found in the 1860 CE map used in PI simulations when vegetation is prescribed (Fig. 4). Interestingly part of these differences between Vmap and Vnone is found in the Southern Hemisphere and the northern edge of the African and Indian monsoon regions. Since there is almost no difference in MH vegetation between Vmap and Vnone, these differences in PI vegetation drive the vegetation differences between MH and PI (Fig. 16). The MH simulated changes seem larger with Vmap. A previous assessment of model results against vegetation and paleoclimate reconstructions (e.g., Harrison et al., 1998, 2014) suggests that MH–PI vegetation for Vmap would be in better agreement with reconstructed changes from observations in terms of forest expansion in the Northern Hemisphere or grasses in the Sahel (Fig. 16). However, the modern vegetation map for this PI-Vmap simulation has even less forest than PI-Vnone north of 55^∘ N (Fig. 4e, f, and i), for which forest is already underestimated (Fig. A2). These differences in PI vegetation only have a small counterpart in climate. These biases correspond to cooler condition in the midlatitudes and high northern latitude (Fig. 15). In the annual mean, there is almost no impact on precipitation (Fig. 15). In terms of climate these two simulations are very similar and closer to each other than to other simulations, whatever the season or the latitudinal band (Fig. A1). The small differences in climate listed above are thus too small to be captured by global metrics. It suggests that there is no direct relationship between the different vegetation maps and model performances. The different vegetation maps are obtained with a similar climate, which indicates that in this model multiple global and vegetation states are possible under preindustrial climate or that tiny climate differences can lead to different vegetation.