Mid-Holocene climate change over China: model–data discrepancy

The mid-Holocene period (MH) has long been an ideal target for the validation of general circulation model (GCM) results against reconstructions gathered in global datasets. These studies aim to test GCM sensitivity, mainly to seasonal changes induced by the orbital parameters (longitude of the perihelion). Despite widespread agreement between model results and data on the MH climate, some important differences still exist. There is no consensus on the continental size (the area of the temperature anomaly) of the MH thermal climate response, which makes regional quantitative reconstruction critical to obtain a comprehensive understanding of the MH climate patterns. Here, we compare the annual and seasonal outputs from the most recent Paleoclimate Modelling Intercomparison Project Phase 3 (PMIP3) models with an updated synthesis of climate reconstruction over China, including, for the first time, a seasonal cycle of temperature and precipitation. Our results indicate that the main discrepancies between model and data for the MH climate are the annual and winter mean temperature. A warmer-than-present climate condition is derived from pollen data for both annual mean temperature (∼0.7 K on average) and winter mean temperature (∼1 K on average), while most of the models provide both colder-than-present annual and winter mean temperature and a relatively warmer summer, showing a linear response driven by the seasonal forcing. By conducting simulations in BIOME4 and CESM, we show that surface processes are the key factors creating the uncertainties between models and data. These results pinpoint the crucial importance of including the non-linear responses of the surface water and energy balance to vegetation changes.

1 Introduction

Much attention from paleoclimate studies has been focused on the current interglacial (the Holocene), especially the mid-Holocene (MH; 6±0.5 ka). The major difference in the experimental configuration between the MH and pre-industrial (PI) arises from the orbital parameters, which brings about an increase in the amplitude of the seasonal cycle of insolation in the Northern Hemisphere and a decrease in the Southern Hemisphere (Berger, 1978). Thus, the MH provides an excellent case study on which to base an evaluation of the climate response to changes in the distribution of insolation. Great efforts have been devoted by the modelling community to designing MH common experiments using similar boundary conditions (Joussaume and Taylor, 1995; Harrison et al., 2002; Braconnot et al., 2007a, b). In addition, much work has been done to reconstruct the paleoclimate change based on different proxies at the global and continental scale (Guiot et al., 1993; Kohfeld and Harrison, 2000; Prentice et al., 2000; Bartlein et al., 2011). The greatest progress in understanding MH climate change and variability has consistently been made by comparing large-scale analyses of data with simulations from global climate models (Joussaume et al., 1999; Liu et al., 2004; Harrison et al., 2014).

However, the source of discrepancies between model results and data is still an open and stimulating question. Two types of inconsistencies have been identified: (1) the model and data showing opposite signs; for instance, paleoclimate evidence from data records indicates an increase of about 0.5 K in global annual mean temperature during the MH compared with PI (Shakun et al.. 2012; Marcott et al., 2013), while there is a cooling trend in model simulations (Liu et al., 2014). (2) The same trend being displayed by both model and data but with different magnitudes. Previous studies have shown that while climate models can successfully reproduce the direction and large-scale patterns of past climate changes, they tend to consistently underestimate the magnitude of change in the monsoons of the Northern Hemisphere as well as the amount of MH precipitation over northern Africa (Braconnot et al., 2012; Harrison et al., 2015). Moreover, significant spatial variability has been noted in both observations and simulations (Peyron et al., 2000; Davis et al., 2003; Braconnot et al., 2007a; Wu et al., 2007; Bartlein et al., 2011). For instance, Marcott et al. (2013) reconstructed a cooling trend of global temperature during the Holocene, mainly from marine records (∼80 %). Based on 642 sub-fossil pollen data points, Marsicek et al. (2018) show that long-term warming defined the Holocene until around 2000 years ago for the European and North American continents. The different trends of pollen- and marine-based reconstructions indicate the spatial variability of annual temperature change during the MH over the globe, which has already been investigated by Bartlein et al. (2010). That makes regional quantitative reconstruction (Davis et al., 2003; Mauri et al., 2015) essential to obtain a comprehensive understanding of the MH climate patterns and to act as a benchmark to evaluate climate models (Fischer and Jungclaus, 2011; Harrison et al., 2014).

China offers two advantages with respect to these issues. The sheer expanse of the country means that the continental response to insolation changes over a large region can be investigated. Moreover, the quantitative reconstruction of seasonal climate changes during the MH, based on the new pollen dataset, provides a unique opportunity to compare the seasonal cycles for models and data. Previous studies indicate that warmer- and wetter-than-present conditions prevailed over China during the MH and that the magnitude of the annual temperature increases varied from 2.4 to 5.8 K spatially, with an annual precipitation increase in the range of 34–267 mm (e.g. Sun et al., 1996; Jiang et al., 2010; Lu et al., 2012; Chen et al., 2015). However, Jiang et al. (2012) clearly show a mismatch between multi-proxy reconstructions and model simulations. In terms of climate anomalies (MH–PI), besides the ∼1 K increase in summer temperature, 35 out of 36 Paleoclimate Modelling Intercomparison Project (PMIP) models reproduce annual (∼0.4 K) and winter temperatures (∼1.4 K) that are colder than the baseline, and a drier-than-baseline climate in some western and middle regions over China is depicted in models (Jiang et al., 2013). Jiang et al. (2012) point out the model–data discrepancy over China during the MH, but the lack of seasonal reconstructions in their study limits comparisons with simulations.

An important issue raised by Liu et al. (2014) is that the discrepancy at the annual level could be due to incorrect reconstructions of the seasonal cycle, a key objective in our paper. Moreover, it has been suggested that the vegetation change can strengthen the temperature response at high latitudes (O'Ishi et al., 2009; Otto et al., 2009), as well as alter the hydrological conditions in the tropics (Z. Liu et al., 2007). However, compared to the substantial land cover changes in the MH derived from pollen datasets (Ni et al., 2010; Yu et al., 2000), the changes in vegetation have not yet been fully quantified and discussed in PMIP3 (Taylor et al., 2012).

In this study, we first present new reconstructions. We use a quantitative biomization method to reconstruct vegetation types during the MH based on a new synthesis of pollen datasets and then use an inverse vegetation model (Guiot et al., 2000; Wu et al., 2007) to obtain the mean annual temperature, the mean temperature of the warmest month (MTWA), the mean temperature of the coldest month (MTCO), the mean annual precipitation, July precipitation and January precipitation over China for the MH. Furthermore, we present a comprehensive evaluation of the PMIP3 simulations performed with state-of-art climate models based on our reconstructions of temperature and precipitation. This is the first time that such progress towards a quantitative seasonal climate comparison for the MH over China has been made. This point is crucial because the MH PMIP3 experiment is essentially one that looks at the response of the models to changes in the seasonality of insolation, and the attempt to derive reconstructions of both summer and winter climate to compare with the simulations will thus enable us to answer the question posed by Liu et al. (2014) on the importance of seasonal reconstructions.

Figure 1 Distribution of pollen sites during the mid-Holocene period in China. The black circle is the original China Quaternary Pollen Database, red circles are digitized ones from published papers and green circles represent the three original pollen data points used in this study. The area in green represents the Tibetan Plateau, and yellow represents the Loess Plateau.

2 Data and methodology

2.1 Data

In this study, we collected 159 pollen records covering most of China for the MH period (6000±500 ¹⁴C yr BP) (Fig. 1). Notably, according to IntCal13 (Reimer et al., 2013), the MH time slice of 6000±500 ¹⁴C yr BP is about 6800 cal BP (the average value), which is not totally consistent with the “mid-Holocene” used in the CMIP5–PMIP3 experiment (6000 cal BP). But for a better comparison with BIOME6000 (in which the MH is defined as 6000±500 ¹⁴C yr BP), we decided to choose the pollen data at 6000±500 ¹⁴C yr BP in our study. In the 159 records, 65 were from the China Quaternary Pollen Database (CQPD, 2000), 3 were original datasets obtained for our study and the others were digitized from pollen diagrams in published papers with a recalculation of pollen percentages based on the total number of terrestrial pollen types. These digitized 91 pollen records were selected according to three criteria: (1) clearly readable pollen diagrams with a reliable chronology and a minimum of three independent age control points since the Last Glacial Maximum (LGM); (2) including the pollen taxa during the 6000±500 ¹⁴C yr BP period with a minimum sampling resolution of 1000 years per sample; and (3) abandoning the pollen records if the published paper mentions the influence of human activity. For the age control, different dating methods are utilized in the collected pollen records, and we applied CalPal 2007 (Weninger et al., 2007) to correct ¹⁴C age into calendar age so that they can be contrasted with each other. For lacustrine records, if the specific carbon pool age is mentioned in the literature, the calendar age is corrected after deducting the carbon pool. Otherwise, the influence of the carbon pool is not considered. The age–depth model for the pollen records was estimated by linear interpolation between adjacent available dates or by regression. Using ranking schemes from the Cooperative Holocene Mapping Project, the quality of dating control for the mid-Holocene was assessed by assigning a rank from 1 to 7; 70 % of the records used in our study fell into the first and second classes (see Table 1 for detailed information) according to the Webb 1–7 standards (Webb III, 1985). Vegetation type was quantitatively reconstructed using biomization (Prentice et al., 1996), following the classification of plant functional types (PFTs) and biome assignment in China by the Members of the China Quaternary Pollen Database (CQPD, 2000), which has been widely tested in surface sediment. The new sites (91 digitized data points and 3 original data points) added to our database improved the spatial coverage of pollen records, especially in the northwest, the Tibetan Plateau, the Loess Plateau and southern regions, where the data in the previous databases are very limited.

Table 1 Basic information on the pollen dataset used in this study.

Modern monthly mean climate variables investigated in this study, including temperature, precipitation and cloudiness (total cloud fraction), have been collected for each modern pollen site based on the datasets (1951–2001) from 657 meteorological observation stations over China (China Climate Bureau, China Ground Meteorological Record Monthly Report, 1951–2001). The MH soil properties and characteristics used in the inverse vegetation model were kept the same with PI conditions, which are derived from the digital world soil map produced by the Food and Agricultural Organization (FAO, 1995). The atmospheric CO₂ concentration for the MH was taken from ice core records (EPICA community members, 2004) and was set at 270 ppmv.

A three-layer back-propagation (BP) artificial neural network (ANN) technique was used for interpolation on each pollen site (Caudill and Butler, 1992). Five input variables (latitude, longitude, elevation, annual precipitation, annual temperature) and one output variable (biome scores) have been chosen in ANN for modern vegetation. The ANN has been calibrated on the training set, and its performance has been evaluated on the verification set (20 %, randomly extracted from the total sets). After a series of training runs, the lowest verification error is obtained with five neurons in the hidden layer after 10 000 iterations. In our study, at each pollen site, we first applied the biomization method to get the biome scores for both the present day and the MH. The anomalies between past (6 ka) and modern vegetation indices (biome scores) were then interpolated to the 0.2×0.2^∘ grid resolution by applying the ANN. After that, the modern grid values were added to the values of the grid of paleo-anomalies to provide gridded paleo-biome indices. Finally, the biome with the highest index was attributed to each grid point. This ANN method is more efficient than many other techniques on the condition that the results are validated by independent datasets, and therefore it has been widely applied in paleoclimatology (Guiot et al., 1996; Peyron et al., 1998). The schematic diagram of ANN (Fig. S1) can be found in the Supplement.

Table 2 Earth's orbital parameters and trace gases as recommended by the PMIP3 project.

Table 3 PMIP3 model characteristics and references.

2.2 Climate models

PMIP, a long-standing initiative, is a climate model evaluation project which provides an efficient mechanism for using global climate models to simulate climate anomalies for past periods and to understand the role of climate feedbacks. In its third phase (PMIP3; Braconnot et al., 2011), the models were identical to those used in the Climate Modelling Intercomparison Project Phase 5 (CMIP5) experiments. The experimental set-up for the mid-Holocene simulations in PMIP3 followed the PMIP protocol (Taylor et al., 2012; https://wiki.lsce.ipsl.fr/pmip3/doku.php/pmip3:design:6k: final, last access: 20 June 2019). The main forcing difference between the MH and PI in PMIP3 is the change in the orbital configuration. More precisely, the orbital configuration in the MH climate has an increased summer insolation and a decreased winter insolation in the Northern Hemisphere compared to the PI climate (Berger, 1978). In addition, the CH₄ concentration is prescribed at 650 ppbv in the MH, while it is set at 760 ppbv in PI (Table 2).

All 13 models (Table 3) from PMIP3 with the MH simulation have been included in our study, including eight atmosphere–ocean (AO) models and five atmosphere–ocean–vegetation (AOV) models. Means for the last 30 years were calculated from the archived time series data on individual model grids for climate variables. Then the near-surface temperature and precipitation flux were bi-linearly interpolated to a uniform 2.5^∘ grid in order to compute the bioclimatic variables (e.g. MAT, MAP, MTWM, MTCO, July precipitation) onto a common grid for comparison with the reconstruction results.

2.3 Vegetation model

The vegetation model BIOME4 is a coupled bio-geography and biogeochemistry model developed by Kaplan et al. (2003). Monthly mean temperature, precipitation, sunshine percentage (an inverse measure of cloud area fraction), absolute minimum temperature, atmospheric CO₂ concentration and subsidiary information about the soil's physical properties like water retention capacity and percolation rates are the main input variables. It represents 13 plant functional types (PFTs), which have different bioclimatic limits. The PFTs are based on physiological attributes and bioclimatic tolerance limits such as heat, moisture and chilling requirements as well as resistance of plants to cold. These limits determine the areas where the PFTs can grow in a given climate. A viable combination of these PFTs defines a particular biome among 28 potential options. These 28 biomes can be further classified into eight megabiomes (Table S1). BIOME4 has been widely utilized to analyse past, present and potential future vegetation patterns (e.g. Bigelow et al., 2003; Diffenbaugh et al., 2003; Song et al., 2005). In this study, we conducted 13 PI and associated MH biome simulations using PMIP3 climate fields (temperature, precipitation and sunshine) as inputs. The climate fields, obtained from PMIP3, are the monthly mean data of the last 30 model years.

Table 4 Important values for each plant life form used in the ΔV statistical calculation as assigned to the megabiomes.

Table 5 Attribute values and the weights for plant life forms used by the ΔV statistic.

2.4 Statistics and interpolation for vegetation distribution

To quantify the differences between simulated (based on BIOME4 forced by the climate model output) and reconstructed (from pollen) megabiomes, a map-based statistic (point-to-point comparison with observations) called ΔV (Sykes et al., 1999; Ni et al., 2000) was applied to our study. ΔV is based on the relative abundance of different plant life forms (e.g. trees, grass, bare ground) and a series of attributes (e. g. evergreen, needle-leaf, tropical, boreal) for each vegetation class. The definitions and attributes of each plant form follow naturally from the BIOME4 structure, and the vegetation attribute values in the ΔV computation were defined for BIOME4 in the same way as for BIOME1 (Sykes et al., 1999). The abundance and attribute values are given in Tables 4 and 5, which describe the typical floristic composition of the biomes. Weighting the attributes is subjective because there is no obvious theoretical basis for assigning relative significance. Transitions between highly dissimilar megabiomes have a weighting of close to 1, whereas transitions between less dissimilar megabiomes are assigned smaller values. The overall dissimilarity between model and data megabiome maps was calculated by averaging ΔV for the grids with pollen data, while the value was set at 0 for any grid without data. ΔV values <0.15 can be considered to point to very good agreement between simulated and actual distributions; 0.15–0.30 is good, 0.30–0.45 is fair, 0.45–0.60 is poor and >0.80 is very poor (adjusted from Zhang et al., 2010).

2.5 Inverse vegetation model

The inverse vegetation model (Guiot et al., 2000; Wu et al., 2007), highly dependent on the BIOME4 model, is applied to our reconstruction. The key concept of this model can be summarized in two points: firstly, a set of transfer functions able to transform the model output into values directly comparable with pollen data is defined. There is no full compatibility between the biome typology of BIOME4 and the biome typology of pollen data. A transfer matrix (Table S2) was defined in our study whereby each BIOME4 vegetation type is assigned a vector of values, one for each pollen vegetation type, ranging from 0 (representing an incompatibility between the BIOME4 type and pollen biome type) to 15 (corresponding to maximum compatibility). Secondly, using an iterative approach, a representative set of climate scenarios compatible with the vegetation records is identified within the climate space, constructed by systematically perturbing the input variables (e.g. ΔT, ΔP) of the model (Table S3).

Table 6 Regression coefficients between the reconstructed climates by inverse vegetation models and observed meteorological values.

The climatic parameters used for regression are the actual values (data source: China Climate Bureau, China Ground Meteorological Record Monthly Report, 1951–2001). MAT is the annual mean temperature, MTCO is the mean temperature of the coldest month, MTWA is the mean temperature of the warmest month, MAP is the annual precipitation, RMSE is the root mean square error of the residuals and ME is the mean error of the residuals. P_jan: precipitation for January; P_jul: precipitation for July.R is the correlation coefficient ± standard error.

The inverse vegetation model (IVM) provides a possibility, for the first time, to reconstruct both annual and seasonal climates for the MH over China. Moreover, it offers a way to consider the impact of CO₂ concentration on competition between PFTs as well as on the relative abundance of taxa and thus makes reconstructions from pollen records more reliable. More detailed information about the IVM can be found in Wu et al. (2007).

We applied the inverse model to modern pollen samples to validate the approach by reconstructing the modern climate at each site and comparing it with the observed values. The high correlation coefficients (R=0.75–0.95), with intercepts close to 0 (except for the mean temperature of the warmest month) and slopes close to 1 (except for the July precipitation), demonstrated that the inversion method worked well for most variables in China (see Table 6).

2.6 Sensitivity test for vegetation feedback

To quantify the vegetation feedback on climate change during the mid-Holocene over China, we performed a sensitivity test using CESM version 1.0.5. This version, developed at the National Center for Atmospheric Research, is a widely used coupled model with dynamic atmosphere (CAM4), land (CLM4), ocean (POP2) and sea-ice (CICE4) components (Gent et al., 2011). Here, we use ∼ 2^∘ resolution for CAM4 (∼ 1.9^∘ for latitude ×2.5^∘ for longitude) in the horizontal direction and 26 layers in the vertical direction. The POP2 adopts a finer grid, with a nominal 1^∘ horizontal resolution and 60 layers in the vertical direction. The land and sea-ice components have the same horizontal grids as the atmosphere and ocean components, respectively.

Two experiments were conducted, including a mid-Holocene (MH) experiment (6 ka) with the original vegetation setting (prescribed as PI vegetation for the MH) and an MH experiment with reconstructed vegetation (6 ka_VEG). In detail, experiment 6 ka used the MH orbital parameters (eccentricity: 0.018682; obliquity: 24.105^∘; longitude of the perihelion: 0.87^∘) and modern vegetation (Salzmann et al., 2008). Compared to experiment 6 ka, experiment 6 ka_VEG used our reconstructed vegetation in China. Except for the changed vegetation, all other boundary conditions were kept unchanged in these two experiments, including the solar constant (1365 W m⁻²), modern topography and ice sheet, and pre-industrial greenhouse gases (CO₂=280 ppmv; CH₄=760 ppbv; N₂O=270 ppbv). Experiment 6 ka was initiated from the default pre-industrial simulation and run for 500 model years. Experiment 6 ka_VEG was initiated from model year 301 of experiment 6 ka and run for another 200 model years. We analysed the computed climatological means of the last 50 model years from each experiment here.

3 Results

3.1 Comparison of annual and seasonal climate changes at the MH

In this study, we collected 159 pollen records, broadly covering the whole of China (Fig. 1). To check the reliability of the collected data, we first categorized our pollen records into megabiomes in line with the standard tables developed for BIOME6000 (Table S1) and compared them with the BIOME6000 dataset (Fig. 2). The match between the collected data and BIOME6000 is more than 90 % (145 out of 159 sites) for both the MH and PI.

Figure 2 Comparison of megabiomes for PI (a, c) and the MH (b, d): (a, b) BIOME6000 and (c, d) pollen data collected in this study.

Based on pollen records, the spatial pattern of climate changes over China during the MH, deduced from IVM, are presented in Fig. 3a–c (points), alongside the results from PMIP3 models (shaded in Fig. 3). For temperature, a warmer-than-present annual climate condition (∼0.7 K on average) is derived from pollen data (the points in Fig. 3a), with the largest increase occurring in the northeast (3–5 K) and a decrease in the northwest and on Tibetan Plateau. On the other hand, the results from a multi-model ensemble (MME) indicate a generally colder annual temperature ($\sim -\mathrm{0.4}$ K on average), with significant cooling in the south and slight warming in the northeast (shaded in Fig. 3a). Of the 13 models, 11 simulate a cooler annual temperature compared with PI than MME. However, two models (HadeGEM2-ES and CNRM-CM5) present the same warmer condition as was found in the reconstruction (Fig. 3d). Compared to the reconstruction, the annual mean temperature during the MH is largely underestimated by most PMIP3 models, which depict an anomaly ranging from $\sim -$1 to ∼0.5 K. Concerning seasonal change, during the MH, MTWA from the data is ∼0.5 K higher than PI, with the largest increase in the northeast and a decrease in the northwest. From model outputs, an average increase of ∼1.2 K is reproduced by MME, with more pronounced warming at high latitudes, which is consistent with the insolation change (Berger, 1978). Figure 3e shows that all 13 models reproduce the same warmer summer temperatures as the data and that HadGEM2-ES and CNRM-CM5 reproduce the largest increases among the models. Although models simulate warmer MTWA, which is consistent with reconstructions, there is a discrepancy between them on MTCO. In Fig. 3c, the data show an overall increase of ∼1 K, with the largest increase occurring in the northeast and a decrease of opposite magnitude on the Tibetan Plateau. Inversely, MME reproduces a decreased MTCO with an average amplitude of $\sim -\mathrm{1.3}$ K, the areas with strongest cooling being the southeast, the Loess Plateau and the northwest. Similarly to the MME, all 13 models simulate a colder-than-present climate with amplitudes ranging from $\sim -\mathrm{2.0}$ K (CCSM4 and FGOALS-g2) to $\sim -\mathrm{0.7}$ K (HadGEM2-ES and CNRM-CM5).

Figure 3 Model–data comparison for annual and seasonal (MTWA and MTCO) temperature (K). For (a–c), points represent the reconstruction from IVM, and shading shows the last 30-year mean simulation results of a multi-model ensemble (MME) for 13 PMIP3 models. The box-and-whisker plots (d–f) show the changes as indicated by each PMIP3 model and the reconstruction. (d) Changes in annual temperature, (e) changes in MTWA and (f) changes in MTCO. The lines in each box show the median value from each set of measurements, the box shows the 25 %–75 % range and the whiskers show the 90 % interval (5th to 95th percentile).

Concerning the annual change in precipitation, the reconstruction shows wetter conditions during the MH across almost the whole of China with the exception of part of the northwest. The southeast presents the largest increase in annual precipitation. All but two models (MIROC-ESM and FGOALS-g2) depict wetter conditions with an amplitude of ∼10 to ∼50 mm. The reconstruction and MME results also indicate an increased annual precipitation during the MH (Fig. 4a), with a much larger magnitude visible in the reconstruction (∼30 mm and ∼230 mm, respectively). The main discrepancy in annual precipitation between simulations and reconstruction occurs in the northeast, which is depicted as drier by the models and wetter by the data. With regard to seasonal change, the reconstruction shows an overall increase in July rainfall (∼50 mm on average), with a decrease in the northwestern regions and East Asian monsoon region in the Yangtze River valley. In line with the reconstruction, the MME also shows an overall increase in rainfall (∼7 mm on average), with a decrease in the northwest for July (Fig. 4b). Notably, a much larger increase is simulated for the south and the Tibetan Plateau by the models, while the opposite pattern emerges along the eastern margin from both models and data. For January precipitation, the reconstruction shows an overall increase in most regions (∼15 mm), except for the northwestern region, while MME indicates a slight decrease (∼3 mm on average). More detailed information about the geographic distribution of simulated temperature and precipitation for each model can be found in Figs. S2–S7.

Figure 4 Model–data comparison for annual, July and January precipitation (mm). For (a, b), points represent the reconstruction from IVM, and shading shows the last 30-year mean simulation results of a multi-model ensemble (MME) for 13 PMIP3 models. The box-and-whisker plots (d–f) show the changes as indicated by each PMIP3 model and the reconstruction. (d) Changes in annual precipitation, (e) changes in July precipitation and (f) changes in January precipitation. The lines in each box show the median value from each set of measurements, the box shows the 25 %–75 % range and the whiskers show the 90 % interval (5th to 95th percentile).

Table S4 provides the biome scores from IVM for pollen data collected from published papers. The reconstructed climate change derived from IVM at each pollen site can be found in Table S5, in which the columns show the median and the 90 % interval (5th to 95th percentage) for feasible climate values produced with the IVM approach. The simulated values for each of the climate variables as shown in the box plots (Figs. 3 and 4) are given in Tables S6 and S7. All the values mentioned above are the mean of the values at 159 pollen sites.

3.2 Comparison of vegetation change at the MH

The use of the PMIP3 database is clearly limited by the different vegetation inputs among the models for the MH period (Table S8). Only HadGEM2-ES and HadGEM2-CC use dynamic vegetation for the MH, and the vegetation of the other 11 models is prescribed to PI with or without interactive leaf area index (LAI), which can introduce a bias to the role of vegetation–atmosphere interaction in the MH climates. To evaluate the model results against the reconstruction for the MH vegetation, we conducted 13 biome simulations with BIOME4 using the PMIP3 climate fields, and the megabiome distribution for each model during the MH is displayed in Fig. 5 (see Fig. S8 for a comparison of PI biomes). To quantify the model–data dissimilarity between megabiomes, a map-based statistic called ΔV (Sykes et al., 1999; Ni et al., 2000) was applied here (see Sect. 2.4).

Figure S9 shows the dissimilarity between simulations and observations for megabiomes during the MH, with the overall values for ΔV ranging from 0.43 (HadGEM2-ES) to 0.55 (IPSL-CM5A-LR). According to the classification of ΔV (see Sect. 2.4) for the 13 models, 12 (all except HadGEM2-ES) showed poor agreement with the observed vegetation distribution. Most models poorly simulate the desert, grassland and tropical forest areas for both periods but perform better for warm mixed forest, tundra and temperate forest. However, this statistic is based on a point-to-point comparison, so the ΔV calculated here cannot represent an estimation of full vegetation simulation due to the uneven distribution of pollen data and the potentially huge difference in the area of each megabiome. For instance, tundra in our data for PI is represented by only 4 points, which counts for a small contribution to ΔV since we averaged it over a total of 159 points, but this calculation could induce a significant bias if these 4 points are representative of a large area of China.

Figure 5 Comparison of interpolated megabiome distributions (plot in the red rectangle) with the simulated spatial pattern from BIOME4 for each model during the mid-Holocene.

So, we used the biome scores based on the artificial neural network technique as described by Guiot et al. (1996) for interpolation (the plots in the red rectangle in Fig. 5) and compared the simulated vegetation distribution from BIOME4 for each model with the interpolated pattern. During the MH, most models are able to capture the tundra on the Tibetan Plateau as well as the combination of warm mixed forest and temperate forest in the southeast. However, all models fail to simulate or underestimate the desert area in the northwest compared to reconstructed data. The main model–data inconsistency in the MH vegetation distribution occurs in the northeast, where data show a mix of grassland and temperate forest, and the models show a mix of grassland and boreal forest.

Figure 6 Changes in the extent of each megabiome as a consequence of simulated climate changes for each model, both expressed as change relative to the PI extent of the same megabiome.

The area statistics carried out for simulated vegetation changes (Fig. 6) reveal that the main difference during the MH, compared with PI, is that grassland replaced boreal forest in large tracts of the northeast (Figs. 5, S8). No other significant difference in vegetation distribution between the two periods was derived from the models. Unlike in the models, three main changes in megabiomes during the MH are depicted by the data. Firstly, the megabiomes were converted from grassland to temperate forest in the northeast. Secondly, a large area of temperate forest was replaced in the southeast by a northward expansion of warm mixed forest. Thirdly, in the northwest and at the northern margin of the Tibetan Plateau, part of the desert area changed into grassland. However, none of the models succeed in capturing these features, especially the transition from grassland into forest in the northeast during the MH. Therefore, this failure to capture vegetation changes between the two periods will lead to a cumulating inconsistency in the model–data comparison for climate anomalies if this computed vegetation was used as a boundary condition in MH climate simulations.

4 Discussion

4.1 Validation and uncertainties of the reconstructions

To investigate the discrepancy between models and reconstructions for the MH climate change over China, the reliability of our reconstruction should first be considered. We therefore compared our reconstruction with previous studies concerning the MH climate change over China based on multiple proxies (including pollen, lake core, paleosol, ice core, peat and stalagmite); the related references and detailed information are listed in the Supplement (Tables S9 and S10). In comparison with PI conditions, most reconstructions reproduced warmer and wetter annual conditions during the MH (Fig. 7), as in our study. In other words, this model–data discrepancy for climate change over China during the MH is common and robust with respect to reconstructions derived from different proxies. Our study reinforces the picture given by the discrepancies between PMIP simulations and pollen data based on a synthesis of the literature.

Figure 7 Comparison between the climate reconstruction and previous reconstruction over China. (a) Previous temperature results. The diamond is the qualitative reconstruction, red is the temperature increase and green is the temperature decrease; the circle is the quantitative reconstruction. (b) Mean annual temperature reconstruction in this study; (c) previous precipitation results. The diamond is the qualitative reconstruction, red is the precipitation decrease and green is the precipitation increase; the circle is the quantitative reconstruction. (d) Mean annual precipitation reconstruction in this study.

However, there could still be some biases due to the reconstruction method. Estimated climates for the present day from IVM were compared with observed climates (Table 6). The slopes and intercepts show a slight bias for annual and January precipitation, while there is a larger bias between the IVM reconstruction and observation for temperature and July precipitation. IVM relies heavily on BIOME4, and since BIOME4 is a global vegetation model, it is possible that the spatial robustness of regional reconstruction could be less than that of global reconstruction due to the failure to simulate local features (Bartlein et al., 2011). Moreover, the output of the model cannot be directly compared to the pollen data, and the conversion of BIOME4 biomes to pollen biomes by the transfer matrix may add a source of uncertainty in the reconstruction. All these biases in reconstruction should be considered in the evaluation of the discrepancy between models and data for climate change during the MH over China.

4.2 Uncertainties for simulations

The discrepancies between models and data for MH climate change can also result from uncertainties in simulation and/or model characteristics. First, the coarse spatial resolution of models can be a factor of discrepancy: previous studies show that the GCMs from PMIP3 are reliable in simulating the geographical distribution of temperature and precipitation over China for the present day. However, compared with observations, most models have topography-related cold biases (Jiang et al., 2016). The climate fields, directly used from the model output without downscaling, will not contain the spatial variability of modern climate in topographically complex areas. Thus, it is necessary to determine to what degree model–data mismatch is related to the rough topography used in the climate models. In PMIP3, MRI-CGCM3 has the highest resolution (atmosphere: 320 × 160 × L48; ocean: 364 × 368 × L51), while IPSL-CM5A-LR has the lowest (atmosphere: 96 × 96 × L39; ocean: 182 × 149 × L31). In Fig. 8, we give the actual modern topography and the interpolated topography used in MRI-CGCM3 and IPSL-CM5A-LR. For MRI-CGCM3, the topography is very close to the observation, so for this model, the model–data discrepancy during the MH over China is not related to the resolution. However, for the model with coarse resolution (IPSL-CM5A-LR), it is true that the coarse version of the model will lead to a bias in topography when regional diversity is discussed. The spatial variations in topography could influence the vegetation and hence the simulated climate. To quantify this impact, we compare the topography and PI climate results of IPSL-CM5A-LR and IPSL-CM5A-MR. Figure 9 shows that the difference in topography caused by model resolution does have an impact on small scales (e.g. southern region of the Tibetan Plateau) but not on the overall pattern. However, small- or regional-scale variations in climate can have a large impact on vegetation and hence reconstructed climate. For better comparison in the future work, downscaled climate variables should be considered.

Figure 8 The topography comparison between models and observations.

Secondly, besides the qualitative consistency among models caused by the protocol of PMIP3 experiments (Table 2), variability in the magnitude of anomalies between models is clearly illustrated by the box plots (Figs. 3 and 4), especially for the temperature anomaly. Figure S10 demonstrates that there is no clear relationship between PI temperature and temperature anomaly (MH–PI). In other words, these disparities in value or even pattern among models are not related to the difference in PI simulations in a simple manner; instead, they reflect the obvious differences in the response by the climate models to the MH forcing, which raises the question of the magnitude of feedbacks among models.

Figure 9 The pre-industrial climate comparison between simulations and observations. “Tas” means temperature 2 m above the surface, and “pr” means precipitation.

Figure 10 Scatter plots showing temperature, cloud cover feedback and surface albedo feedback changes during the MH. The values shown are the simulated 30-year mean anomaly (MH–PI) for the 13 models. (a) Annual mean temperature relative to the annual mean cloud cover feedback and (d) annual surface albedo feedback. (b) Summer (JJA) mean temperature relative to the summer mean cloud cover feedback and (e) summer surface albedo feedback. (c) Winter (DJF) mean temperature relative to the summer mean cloud cover feedback and (f) winter surface albedo feedback. The horizontal and vertical lines in the plots represent the value of 0.

As positive feedbacks between climate and vegetation are important to explain regional climate changes, the failure of the models to represent the amplitude and pattern of the observed vegetation differences (see Sect. 3.2) could amplify and partly account for the model–data disparities in climate change, mainly due to variations in the albedo. Because HadGEM2-ES and HadGEM2-CC are the only two models in PMIP3 with a dynamic vegetation simulation for the MH, we focused on these models to examine the variations in vegetation fraction in the simulations. The main vegetation changes during the MH demonstrated by HadGEM2-ES are an increased tree coverage (∼15 %) and a decreased bare soil fraction (∼6 %), while HadGEM2-CC depicts a ∼3 % decrease in tree fraction and a ∼1 % increase in bare soil (Fig. S11). We made a rough calculation of albedo variance caused solely by vegetation change for both models and for our reconstruction based on the area fraction and albedo value of each vegetation type (Betts, 2000; Bonfils et al., 2001; Oguntunde et al., 2006; Bonan, 2008). The overall albedo change from the vegetation reconstruction during the MH shows a ∼1.8 % decrease when snow free, with a much larger impact (∼4.2 % decrease) when snow covered. The results from HadGEM2-ES are highly consistent with the albedo changes from the reconstruction, featuring a ∼1.4 % (∼6.5 %) decrease without (with) snow, while HadGEM2-CC produces an increased albedo value during the MH (∼0.22 % for snow free, ∼1.9 % with snow cover), depending on its vegetation simulation. Two ideas could be inferred from this calculation: (1) HadGEM2-ES is much better in simulating the MH vegetation changes than HadGEM2-CC, and (2) the failure by models to capture these vegetation changes will result in a much larger impact on winter albedo (with snow) than summer albedo (without snow). In conclusion, there is an obvious advantage of using AOVGCM instead of AOGCM when discussing the MH climate, but the premise is that the AOVGCM can simulate an accurate vegetation distribution.

These surface albedo changes due to vegetation changes could have a cumulative effect on the regional climate by modifying the radiative fluxes. For instance, the spread of trees into the grassland biome in the northeast during the MH, revealed by the reconstruction in our study, should act as a positive feedback to climate warming by increasing the surface net shortwave radiation associated with reductions in albedo due to taller and darker canopies (Chapin et al., 2005). Previous studies show that cloud and surface albedo feedbacks on radiation are major drivers of differences between model outputs for past climates. Moreover, the land surface feedback shows large disparities among models (Braconnot and Kageyama, 2015).

We used a simplified approach (Taylor et al., 2007) to quantify the feedbacks and to compare model behaviour for the MH, thus justifying the focus on surface albedo and atmospheric scattering (mainly accounting for cloud change). Surface albedo and cloud change are calculated using the simulated incoming and outgoing radiative fluxes at the Earth's surface and at the top of the atmosphere (TOA) based on data for the last 30 years averaged from all models. Using this framework, we quantified the effect of changes in albedo on the net shortwave flux at TOA (Braconnot and Kageyama, 2015) and further investigated the relationship between these changes and temperature change. Figure 10 shows that most models produced a negative cloud cover and surface albedo feedback on the annual mean shortwave radiative forcing. Concerning seasonal change, the shortwave cloud and surface feedback in most models tend to counteract the insolation forcing during the boreal summer, while they enhance the solar forcing during winter. A strong positive correlation between albedo feedback and temperature change is depicted, with a large spread in the models owing to the difference in albedo in the 13 models. In particular, CNRM-CM5 and HadGEM2-ES capture higher values of cloud and surface albedo feedback, which could be the reason for the reversal of the decreased annual temperature seen in other models (Fig. 3d).

Figure 11 Climate anomalies between the two experiments (6 ka and 6 ka_VEG) conducted in CESM version 1.0.5. The anomalies (6 ka_VEG-6 ka) of temperature and precipitation at both the annual and seasonal scale are presented, and all these climate variables are calculated as the last 50-year means from two simulations.

However, the vegetation patterns produced by BIOME4 in Fig. 5 are not used in the PMIP3 experimental set-up. They are only determined by the input variables from models. Therefore, the disagreements on the MH vegetation pattern are possibly inherited from the PI. To better quantify the vegetation–climate feedback, two experiments were conducted in CESM version 1.0.5, including a mid-Holocene (MH) experiment (6 ka) with the original vegetation setting (prescribed as PI vegetation for the MH) and an MH experiment with reconstructed vegetation (6 ka_VEG). Figure 11 shows the climate anomalies (6 ka_VEG minus 6 ka) between two simulations for both the annual and seasonal scale. For temperature, it is clear that the 6 ka_VEG simulation reproduces a warmer annual mean climate (∼0.3 K on average) as well as an obviously warmer winter (∼0.6 K on average). For precipitation, the reconstructed vegetation leads to more annual and seasonal precipitation, which can also reconcile the model–data discrepancy of an increased amplitude for precipitation during the MH (the data reproduced a larger amplitude than the model, as revealed by our study). So the mismatch between models and data in terms of MH vegetation could partly account for the discrepancy in climate due to the interaction between vegetation and climate through radiative and hydrological forcing with albedo. These results highlight the value of building a new generation of models able to capture not only the atmosphere and ocean response, but also the non-linear responses of vegetation and hydrology to climate change.

5 Conclusion

In this study, we compare the annual and seasonal outputs from the PMIP3 models with an updated synthesis of climate reconstructions over China, including, for the first time, the seasonal cycle of temperature and precipitation. In response to the seasonal insolation change prescribed in PMIP3 for the MH, all models produce similar large-scale patterns for seasonal temperature and precipitation (higher than present July precipitation and MTWA, lower than present MTCO). The main discrepancy emerging from the model–data comparison occurs for the mean annual temperature and MTCO; data show an increased value and most models simulate the opposite except CNRM-CM5 and HadGEM2-ES, which reproduced the higher-than-present MH annual temperature. By conducting simulations with BIOME4 and CESM, we show that surface processes are the key factors explaining the discrepancies between models and data. These results show the importance of including the non-linear responses of the surface water and energy balance related to vegetation changes. However, it should also be noted that prescribing the vegetation with reconstructed biomes would reduce the power of the biome-based climate reconstruction owing to the potential circularity (prescribe the vegetation to get the vegetation). Moreover, besides the vegetation influence, the impact of rough topography, soil type and other possible factors on model–data discrepancy remains to be investigated in future work.