The effect of a dynamic soil scheme on the climate of the mid-Holocene and the Last Glacial Maximum

In order to account for coupled climate–soil processes, we have developed a soil scheme which is asynchronously coupled to a comprehensive climate model with dynamic vegetation. This scheme considers vegetation as the primary control of changes in physical soil characteristics. We test the scheme for a warmer (mid-Holocene) and colder (Last Glacial Maximum) climate relative to the preindustrial climate. We find that the computed changes in physical soil characteristics lead to significant amplification of global climate anomalies, representing a positive feedback. The inclusion of the soil feedback yields an extra surface warming of 0.24 C for the mid-Holocene and an additional global cooling of 1.07 C for the Last Glacial Maximum. Transition zones such as desert–savannah and taiga–tundra exhibit a pronounced response in the model version with dynamic soil properties. Energy balance model analyses reveal that our soil scheme amplifies the temperature anomalies in the mid-to-high northern latitudes via changes in the planetary albedo and the effective longwave emissivity. As a result of the modified soil treatment and the positive feedback to climate, part of the underestimated mid-Holocene temperature response to orbital forcing can be reconciled in the model.


Introduction
The signature of climate variability on orbital timescales (i.e.mid-Holocene, 6 kyr ago, and Last Glacial Maximum, 21 kyr ago) as derived from proxy archives serves as a test bed for a critical evaluation of general circulation models (GCMs; see, e.g., Braconnot et al., 2007Braconnot et al., , 2012;;Hargreaves et al., 2013;Lohmann et al., 2013;Harrison et al., 2014).In this context, terrestrial vegetation dynamics have been identified as a crucial part of physical land surface characteristics within climate models in order to assess climate variability (Kutzbach et al., 1996;Texier et al., 1997Texier et al., , 2000;;Claussen 1997;Claussen and Gayler, 1997;Braconnot et al., 1999Braconnot et al., , 2007)).Subsequently, land surface schemes with interactive vegetation have been implemented into GCMs (e.g.Foley et al., 1998Foley et al., , 2000)).In addition to vegetation dynamics, however, various studies also highlight the impact of physical soil characteristics on climate, as it represents another component of the land surface (Kutzbach et al., 1996;Doherty et al., 2000;Levis et al., 2004;Knorr and Schnitzler, 2006;Shellito and Sloan, 2006;Wanner et al., 2008;Brovkin et al., 2009;Micheels et al., 2009;Knorr et al., 2011;Krapp and Jungclaus, 2011;Vamborg et al., 2011;Lohmann et al., 2015).So far, model studies designed to examine soil dynamics have mainly focused on the sensitivity of individual physical soil characteristics and have therefore primarily addressed the monocausal impact of land albedo (Bonfils et al., 2001;Levis et al., 2004;Knorr and Schnitzler, 2006;Schurgers et al., 2007;Vamborg et al., 2011), soil texture, and maximum water-holding field capacity (Wang, 1999;Levis et al., 2004) on climate.Furthermore, Jiang (2008) investigated the role of physical soil characteristics (soil colour, soil texture) for the Last Glacial Maximum (LGM) by applying an atmosphere GCM, which was asynchronously coupled to a terrestrial biosphere model.In this procedure soil characteristics in a grid cell were iteratively adapted if the dominant plant functional Published by Copernicus Publications on behalf of the European Geosciences Union.type had changed.The integrated vegetation-soil feedback reinforced glacial cooling to some extent (Jiang, 2008).Further progress is shown in Vamborg et al. (2011), in which the authors implemented a dynamic background albedo scheme within a comprehensive fully coupled GCM with vegetation dynamics (earth system model, ESM).In this model, the model additionally accounts for land surface albedo changes based on the calculation of carbon fluxes between terrestrial carbon pools (vegetation, litter and phytomass, and soil).
So far, various studies have focused on the climate effect of physical soil parameter changes.However, the combined impact and feedback of soil albedo, texture, and total waterholding field capacity within the fully coupled atmospherevegetation-ocean system remains to be addressed.The current inability of comprehensive state-of-the-art ESMs to compute the dynamics of physical soil characteristics reveals a potential failure to represent paleoclimate dynamics and quantifying paleoclimate sensitivity.
Here, we investigate the impact of dynamic soils on the climate of the mid-Holocene and the Last Glacial Maximum by means of the comprehensive ESM with vegetation dynamics, COSMOS (Community of Earth System Models; Jungclaus et al., 2006;Raddatz et al., 2007).We develop a simple soil scheme which incorporates the computation of soil albedo, soil texture, and total water-holding field capacity and couple this asynchronously to the ESM.The coupled version of COSMOS allows us to examine the potential soil feedback in the atmosphere-vegetation-ocean system.
The model approach utilizes boundary conditions of the preindustrial period, the mid-Holocene, and the Last Glacial Maximum (LGM).The simulations of the extended model set-up are further compared to the respective standard timeslice studies (Wei and Lohmann, 2012;Zhang et al., 2013) to show that dynamic soils drive significant positive feedbacks in the coupled climate system for past climate states.

Methods
We use the ESM configuration presented by Jungclaus et al. (2006) and describe the design of the soil scheme and the asynchronous coupling procedure.The experimental design includes a short description of the applied model boundary conditions that are representative of the mid-Holocene and LGM (Wei and Lohmann, 2012;Zhang et al., 2013).Furthermore, a forcing factor separation method (Stein and Alpert, 1993;Zhang et al., 2015) has been applied to isolate the potential soil feedback in the climate system.

Model description
In the current ESM set-up, three model components (atmosphere, sea-ice-ocean, land surface) are coupled to each other within the COSMOS model suite.The atmospheric component of COSMOS is ECHAM5 (Roeckner et al., 2003), with a spectral resolution of T31L19 (3.75 × 3.75 • horizontal grid-space resolution and 19 vertical layers), which enables a direct comparison with previous studies of the mid-Holocene and LGM (Wei and Lohmann, 2012;Zhang et al., 2013).The ocean model MPIOM uses an orthogonal curvilinear grid (3 • × 1.5 • lateral resolution) with a higher resolution at the grid poles (Greenland, Antarctica) and 40 unevenly spaced vertical layers (Marsland et al., 2003).The transfer of fluxes and momentum between atmosphere and ocean is handled by the coupler OASIS3 without any flux corrections (Jungclaus et al., 2006).The modular land surface scheme JSBACH (Raddatz et al., 2007) with vegetation dynamics (Brovkin et al., 2009) is embedded in the ECHAM5 atmosphere model.JSBACH is based on the semiempirical terrestrial ecosystem model BETHY, which incorporates an energy and water balance, photosynthesis, phenology, and a carbon balance compartment (Knorr, 2000).JSBACH employs a subgrid-scale tiling approach in which the actual cover fraction of each plant functional type (PFT) is associated with a tile within each grid cell.In our study, JSBACH uses the standard model configuration of eight PFTs.The land surface energy balance considers soil albedo, leaf area index of the actual PFT, the albedo of stems, snow fraction on the ground and on the canopy, and the masking effect of snow under the canopy.Further background information on COSMOS and its components can be found in Jungclaus et al. (2006).So far, COSMOS has been applied to various palaeo set-ups such as the Holocene (Wei et al., 2012;Wei and Lohmann, 2012;Varma et al., 2012;Lohmann et al., 2013), Last Glacial Maximum and Marine Isotope Stage 3 (Gong et al., 2013;Zhang et al., 2013Zhang et al., , 2014)), last interglacial (Lunt et al., 2013), Pliocene (Stepanek and Lohmann, 2012), and Miocene (Knorr et al., 2011;Knorr and Lohmann, 2014).

Design of the soil scheme
The design of the soil scheme is based on the concept that an organism actively regulates and transforms its abiotic environment (Lovelock, 1979).We therefore specify soil types that correspond to PFTs.Applying this concept, the soil scheme diagnoses soil types by the composition of computed PFTs for each grid cell.The computed PFTs are derived from the land surface model JSBACH (Fig. 1; see Sect.2.4).This procedure does not explicitly capture all pedogenic factors (e.g.parent rock; Targulian and Krasilnikov, 2007); however, some may be indirectly considered by vegetation demands (e.g.plant water accessibility, viable temperature limits).The aim of the soil type classification is to provide constraints of physical characteristics (total water-holding field capacity, albedo, and texture) for the land surface scheme of COSMOS and therefore does not refer to the pedologic nomenclature.
First, the soil scheme takes the mean output of 100 model years, provided by COSMOS, to calculate the physical soil properties (soil albedo in the spectrum of visible and nearinfrared light, soil texture, total water-holding field capacity of soils).Thereafter, the model integrates 600 years in total, including an asynchronous coupling time step of 200 years between COSMOS and the soil scheme.Thus, for a model integration period of 600 years, the soil scheme iterates three times in total.The asynchronous coupling is defined by using the previous 100-model-year mean of COSMOS as input for the soil scheme (Fig. 1).Within this study, the coupled version of the ESM with the soil scheme is called COS-MOS_soil.At the end of the model simulation, the final 100 years of model integration, including the physical soil characteristics of the third iteration are used for analysis.

Deriving physical soil characteristics
The physical soil characteristics of JSBACH comprise snowfree soil albedo (α) in the visible (α vis , 0.3-0.7 µm) and nearinfrared (α nir , 0.7-3 µm) light spectrum, total water-holding field capacity of soil (h cws ), and soil data flags (soil texture) from the United Nations Food and Agriculture Organization (FAO) soil classification (Hagemann et al., 1999;Hagemann, 2002;Rechid et al., 2009).The original data set of the snow-free soil albedo is derived from modified satellite measurements of the Moderate Resolution Imaging Spectroradiometer (MODIS; Rechid et al., 2009).The FAO soil data flags (refer to Gildea and Moore in Henderson- Sellers et al., 1986) represent a general grain size classification ranging from sand to clay.The assessment of the total water-holding field capacity of soils (h cws ) reveals considerable uncertainties and is estimated by a data set of the parameter h pwp (permanent wilting point) and h ava (maximum amount of water that plants may extract from the soil before they start to wilt), which is based on optimised rooting depths (Hagemann et al., 1999).
To derive physical soil characteristics that agree with vegetation changes, we use calculated PFTs from the final 100year mean of a preindustrial model output (Wei et al., 2012;Fig. 2).Physical soil characteristics are specified for all dominant PFTs (defined by a cover fraction > 50 % of the grid cell) by global mean calculation (Table 1).Additionally, we define a soil type that represents global desert areas (defined by the cover fraction that lacks vegetation cover for > 50 model years) and a subgroup specifying the Arabian Peninsula-Sahara desert region (defined by < 50 % vegetation cover within a zonal belt of 13-35.26• N).The desert subgroup is mandatory when considering extremes of h cws and soil albedo values of bare soils (Fig. 2, yellow crosses).Further, we subdivide C3 grass into "warm-tolerance" (defined by the mean annual temperature, MAT, > 0 • C) and "cold-tolerance grass" (MAT < 0 • C).In total, the adjusted classification consists of 10 soil types (Fig. 2, Table 1).The FAO soil data flags constitute a coarse-soil texture classification, ranging from 1 (sand) to 5 (clay), and are rather insensitive to the specified soil types (Table 1).Though calculated by the soil scheme, we do not further consider this variable in our study.

Soil scheme
The soil scheme computes grid-cell-specific soil characteristics sol n (n refers to h cws , fao, α s,vis , α s,nir ) according to the soil classification described in Sect.2.3.The soil characteristic sol n is calculated by a soil index.The soil index includes the cover fraction of the soil types (f i ), which are generally derived by PFT-specific cover fractions (Table 1).The soil index weights the soil characteristic of a specific soil type (f i sol i ) by the cumulative cover fraction of all soil types 10 i=1 f i per grid cell: (1) If 10 i=1 f i does not occupy the full grid cell space, the residual is interpreted as desert (see Sect. 2.3, Table 1).The weighted soil characteristics are summed to a grid cell average of soil characteristics and thereafter transferred to the atmospheric component of COSMOS.
gases as compared to the preindustrial period (Wei et al., 2012).Adjustments of LGM boundary conditions include orbital parameter settings, greenhouse gases, sea-level and ice sheet changes, and an initial salinity change in the glacial ocean by +1 ‰ (Zhang et al., 2013).The length of the simulations PI_ctl, HOL_ctl, and LGM_ctl amounts to at least 3000 years for each model run (Wei et al., 2012;Zhang et al., 2013).The model simulations, including dynamic soils, refer to the respective time-slice model simulation (PI_ctl, HOL_ctl, LGM_ctl) and continue the integration for another 600 years in total.All model output shows mean climate conditions for the final 100 years of the simulations.

Analysis of the soil feedback to the climate system
The respective soil feedback to the climate of the mid-Holocene and the LGM is analysed by using a factor sepa-ration method for identifying synergies in numerical models (Stein and Alpert, 1993;Lunt et al., 2012;Zhang et al., 2015).The method calculates the effect between applied changes to a reference state with respect to an unperturbed reference state (f 0,0 ), e.g. the calculated single effect of changing two forcing factors f 1,0 and f 0,1 separately, as well as the combination of both factors f 1,1 .Thereafter, the synergetic effect f1,1 between both factors can be calculated (Stein and Alpert, 1993).To isolate the synergetic effect f1,1 between the two factors, the anomaly term of their linear combinations is calculated: (2) Assuming a linear (climate) system between two factors, the combination of both factors is identical to the sum of each single factor; thus, the synergy term equals 0. Devi-  By means of the factor separation method, the model simulations are factorised representing the background climate (f climate,0 ) and dynamic soils (f 0,soil ).Thereafter, we disentangle the joint effect of climate and dynamic soils (f climate,soil ), which we call soil feedback.The soil feedback is indicated by fclimate,soil : The factors f 0,0 and f 0,soil refer to the preindustrial model run without (PI_ctl) and with dynamic soils (PI_sol), respectively.The defined background climate of the model run f climate,0 represents several changes in the boundary conditions with respect to the preindustrial period that are applied for each time slice (mid-Holocene, LGM; see Sect.2.5), and the "dynamic soils" factor f 0,soil refers to model runs, which have been performed including the interactive soil scheme (COSMOS_soil).The combination of both factors is given by f climate,soil .
The forcing of the mid-Holocene model run is mainly driven by insolation anomalies as compared to the preindustrial period; therefore, the soil feedback is largely attributed to changes in orbital forcing.However, the multitude of changes in the LGM boundary conditions that is representative of the LGM background climate (see Sect. 2.5) does not allow the identification of a monocausal relationship.The calculated soil feedback of the respective model studies is presented in Sects.3.1 and 3.2.

Results
This section gives a brief presentation on the climate anomalies, which are introduced by coupling of the soil scheme and the ESM for the preindustrial period.In the next step, the simulated mid-Holocene and LGM climate that are derived from the extended model set-up including dynamic soils (COSMOS_soil) are described, and, thereafter the soil feedback to the respective climate state is shown.

Preindustrial simulation (PI_sol -PI_ctl)
The preindustrial mean surface air temperature (SAT), precipitation, and evaporation anomalies between COS-MOS_soil and COSMOS (PI_sol -PI_ctl) show the effect of coupling the soil scheme to the ESM (Fig. 3).In general, the global mean SAT increases by 0.27 • C. High northern latitude SATs locally increase by about 1 • C at the Canadian Arctic Archipelago, and Baffin Bay is affected by reduced sea-ice cover as a consequence of changes in the soil characteristics.Regional warming in remote areas like the Southern Ocean is sensitive to sea-ice changes as a function of anomalous wind fields and ocean current shifts.
In parallel with boreal sea-ice retreat, the soil scheme particularly favours the expansion of forests in the high northern  steady-state conditions.Therefore, soil characteristics in a final state may not reflect the same status of present soil characteristics.For a more detailed discussion, the reader is referred to the beginning of Sect. 4.

Mid-Holocene climate (HOL_sol -PI_sol)
The application of Holocene boundary conditions to COS-MOS_soil (HOL_sol) leads to a global surface air temperature (SAT) warming of 0.34 • C with respect to PI_sol.The anomalous SAT warming of the standard Holocene model simulation (HOL_ctl -PI_ctl: +0.1 • C) undergoes a substantial temperature rise by including the soil feedback ( fHOL,soil = +0.24• C).High northern latitudes exhibit largest amplitudes, i.e. the Barents Sea is affected by sea-ice retreats (Fig. 4).Hatched areas shown in Fig. 4 and the subsequent Figs. 5 and 6 mark anomalies in regions with no significant climate change (significance level: 0.05).A general feature of the mid-Holocene climate pattern is given by negative temperature anomalies, accompanied by significant precipitation (Fig. 5) and evaporation anomalies (Fig. 6), which arise especially in regions that are sensitive to the African and Asian monsoon.In principle, the mid-Holocene land surface is basically characterised by higher forest (+2.3 %) and grass abundance (+1.2 %) at the expense of desert fraction (−3.9 %).More specifically, the forest cover expands in northern Canada, northern Siberia and at the Sahara-Sahel transition; this expansion is accompanied by anomalous warming (Fig. 4) and increasing precipitation (Fig. 5), respectively, or a combination of these two factors.Simultaneously, decreased forest fraction is simulated in regions exhibiting negative precipitation anomalies with respect to PI_sol, as can be exemplarily seen in the South America tropical rainforests (Fig. 7).
Figure 8 shows the simulation of soil albedo (α) and total water-holding field capacity (h cws ) under different climate background conditions of the mid-Holocene and LGM.For the mid-Holocene, the global mean of h cws increases by 2.38 cm, which is most pronounced in northern Africa at the Sahara-Sahel transition.As a result, the global soil moisture budget increases by 1.85 cm as compared to PI_sol (Fig. 9); this increase is mainly attributed to the soil feedback (+68 %).A change in the physical soil characteristics (soil feedback) enhances the atmospheric hydrological transport from the ocean towards the continent (+1.29 kg m −2 yr −1 ), which in turn contributes to elevated land surface runoff and drainage (+3.35 kg m −2 yr −1 ).
Globally, the land surface albedo decreases by 0.011, with maximum temperature anomalies at the transition zone of the Sahara-Sahel region.This is accompanied by changes in soil characteristics and a northward migration of vegetation towards the Sahara desert (Fig. 4).The high northern latitudes, in particular, undergo a "darkening" of the land surface due to the integrated effect of lowered soil albedo (Fig. 8), albedo effects caused by the expansion of forest cover (Fig. 7), and regional snow depth changes (given in millimetres of water equivalent) in the continental realm of the Barents Sea (Fig. 10).The replacement of C3 grass by forest in these locations increases vegetation canopy and governs associated changes in soil characteristics: a higher abundance of boreal evergreen forests in the high northern latitudes leads to a decreased soil albedo and to an increase in the total waterholding field capacity relative to soils that are associated with C3 grasses.

LGM climate (LGM_sol -PI_sol)
The LGM climate, as calculated in COSMOS_soil (LGM_sol), reveals an overall global SAT cooling of 7.03 • C. The impact of the implemented soil scheme ( fLGM,soil ) yields an additional cooling of 1.07 • C with maximum anomalies in the high northern latitudes and localised cooling anomalies in the low-latitude desert regions (Fig. 4).Like the overall temperature decrease, the water vapour capacity in the atmosphere -as a function of temperature that is given by the Clausius-Clapeyron relation -is strongly reduced, which therefore dampens precipitation and evaporation fluxes in a glacial climate (Fig. 5, Fig. 6).Total mean precipitation decreases by 0.42 mm day −1 with pronounced amplitudes at the inner tropical convergence zone near the equator and at latitudes > 40 • .The equatorward-expanding sea-ice cover reduces precipitation in the high latitudes and governs increased rainout at the sea-ice edge (Fig. 5).Precipitation and snow depth decrease especially in the western part of northern Siberia, whereas regional cooling governs increased snow depth in eastern Siberia and close to the Fennoscandian ice sheet (Fig. 10).The global forest cover decrease (−14.9 %) is compensated for by the expansion of grass (+10 %) and desert area (+12.7 %).Given that subaerial shelf areas developed due to the sea-level drop, forest (+6 %) and grass (+7 %) can evolve there.Major retreats in boreal forest cover arise in Eurasia in a zonal belt between 50-65 • N and in proximity to the Scandinavian and Laurentide ice sheets (Fig. 11).Cold desert areas as well as snow depth increase in the high northern latitudes (Fig. 10).Interestingly, forest cover increases slightly along ∼ 30 • S most likely as a response to greater precipitation there.Tropical forest cover is drastically reduced by 59 % compared to PI_sol.In most areas grasses replace forests, and deserts (e.g. in the Arabian Peninsula, Sahara-Sahel borderline, and northern Siberia) expand at the expense of grasses.
The computed change in physical soil characteristics shows a general decrease in h cws (−8.3 cm; Fig. 8b) as a response to climate and vegetation change.However, this effect is mainly compensated for (global mean h cws , including subaerial shelf regions: −0.4 cm) by considering additional soil processes in subaerial regions, as an effect of the sea level drop (Fig. 8b).The strongest decline in total water-holding field capacities of soils is localised at the southern borderline of the Eurasian and Laurentide ice sheets and the Sahara desert.This responds to a decrease in surface temperatures and precipitation, which is in line with a reduction in the forest fraction and a parallel expansion of the desert area (Fig. 11a, e).
Comparing Figs.8b and 9b, the total water-holding field capacity of soils constitutes a direct soil wetness control.Thereby, elevated soil wetness along the ∼ 30 • S latitudinal band (Fig. 9b) is not directly linked to h cws but rather reflects a change in atmospheric circulation that is associated with in-creased precipitation, evaporation, and forest cover (Figs. 5c,6c,11a).
The global change in soil albedo (Fig. 8b) yields a remarkable brightening of the land surface, with pronounced amplifications in the high latitudes (Fig. 9d).This generally accompanies desert areas expanding at the expense of grass and forest cover (Fig. 11) and thus gives a rise to the overall shortwave reflectivity of the land surface (Fig. 9d).Especially the Siberian region provides maximum changes in soil and land albedo (and h cws ) that are associated with a southward shift of the taiga-tundra transition.Amplified by changes in the physical soil characteristics at the transition zone, the latitudinal shift is well represented by a general retreat of Siberian forests and an additional southward shift of the boreal treeline (Fig. 11b).North of the boreal treeline (> 60 • N), the incorporated soil feedback promotes an expansion of the polar desert (Fig. 11f) and an increase in the snow cover (Fig. 10d) that parallels a demise of grass cover there (Fig. 11d).The integrated feedback effect of boreal soil changes that interacts with vegetation snow cover and sea ice is reflected by additional regional surface air temperature cooling of more than 3 • C (Fig 4d)

Energy balance model
To separate the effect of radiative components of the isolated soil feedback on the respective climate, we additionally present results of a zero-dimensional energy balance model.The radiative balance of the soil feedback is quantified by applying a zero-dimensional energy balance model (0-dim EBM) for a grey atmosphere (Sellers, 1969): The 0-dim EBM comprises the total irradiance (S = 1367 Wm −2 ), the Stefan-Boltzmann constant (σ = 5.67 × 10 −8 Wm −2 K −4 ), the planetary albedo α, and the effective longwave emissivity ε.The planetary albedo is calculated by the ratio of incoming and outgoing shortwave radiation at the top of the atmosphere, and the effective longwave emissivity is given by the fraction of outgoing longwave radiation at the top of the atmosphere and the surface -all radiative fluxes are provided by the globally and regionally averaged model output.Thereafter, the radiative equilibrium temperature T is derived by the 0-dim EBM.
For the preindustrial period (PI_ctl) the 0-dim EBM constitutes a planetary albedo of 0.32 and an effective longwave emissivity of 0.584, which equals a radiative equilibrium temperature of 289.4 K, in line with previous EBM considerations (Heinemann et al., 2009).To compare the radiative contributions in the prevailing model studies, α and ε are expressed in terms of temperature anomalies with respect to the preindustrial period (Fig. 12).Motivated by the 0-dim EBM, we apply radiation diagnostics at the surface on a regional scale comprising the tropics (30 • S-30 • N) and mid-to-high northern latitudes (30-90 • N).
For the mid-Holocene and LGM (without dynamic soils), the 0-dim EBM calculates a similar warming (+0.albedo and effective longwave emissivity contribute to an amplification of the original temperature anomaly, taking dynamic soils into account.The anomaly of COSMOS_soil and COSMOS with respect to the preindustrial period (PI_sol -PI_ctl: +0.1 • C) shows the error of procedure by including the soil scheme in the model set-up (see Sect. 3.1).This anomaly accounts for the offset between model runs with and without dynamic soils and should be considered in the comparison of the two model versions.As a reference, the regional and global deviations regarding the error of procedure (PI sol ) are displayed as well (Fig. 12).Nevertheless, we state that anomalies between model studies with (COSMOS_soil) and without the interactive dynamic soil component (COS-MOS) consistently reveal stronger soil effects as compared to the error of procedure (PI sol ).Interestingly, based on the 0-dim EBM analysis the Holocene temperature rise without dynamic soils (+0.1 • C) is solely attributed to a reduction in the effective longwave emissivity, corroborated by an increase in the water vapour content in the atmosphere (+0.28 kg m −2 ).Dynamic soils in the mid-Holocene simulation additionally amplify the Holocene warming by changes in both α (+0.14 • C) and ε (+0.19 • C).The EBM, which is applied to low latitudes (30 • S-30 • N) of the Holocene, calculates minor temperature reductions (HOL_ctl -PI_ctl: −0.1 • C).Moderate cooling in the low latitudes is largely associated with orbital forcing during the mid-Holocene, by insolation redistribution from low to the high latitudes (Berger and Loutre, 1991;Laepple and Lohmann, 2009).Instead, including dynamic soils in the climate system, the EBM calculates a change in temperature that overcompensates for the initial cooling of the standard model set-up (HOL_sol -PI_ctl: +0.3 • C, Fig. 12b).In the mid-to-high northern latitudes (30-90 • N), the surface radiation diagnostics exerts amplified temperature changes throughout all of the model studies as compared to the global climate signature (Fig. 12c).The boreal temperature change (30-90 • N) is further corroborated by dynamic soils reinforcing the trend of the original signal in HOL_ctl and LGM_ctl by changes in α (HOL: +27 %; LGM: +21 %) and ε (HOL: Generally, changes in the planetary albedo of the model are caused by the integrated albedo changes in cloud cover, snow, sea ice, vegetation, and the effect of dynamic soils shown.Interestingly, additional boreal planetary albedo changes (Fig. 12c) that are associated with the soil feedback are relatively weak for the mid-Holocene but much stronger for the Last Glacial Maximum.This fact may be best explained by the limited changes in the vegetation canopy for the Holocene (Fig. 7) as compared to the combined effect of vegetation reductions and the exposition of relatively high albedo soils during the LGM (Fig. 11).
A change in the effective longwave emissivity is fundamentally affected by changes in radiatively active gas concentrations, water vapour and clouds in the atmosphere.Especially the dynamic changes in the total water-holding field capacity of soils affect net evaporation and precipitation fluxes between the atmosphere and the land surface, hence controlling part of the water vapour content in the atmo- sphere.In turn, additional temperature anomalies that are related to soil albedo changes feed back to the total capacity of water vapour in the atmosphere.However, the present set of model studies does not allow a quantification of the separated feedbacks between the atmosphere and the total waterholding field capacity of soils and soil albedo.For the LGM temperature, the effective longwave emissivity change that is induced by dynamic soils is weaker than compared to the effect of soil albedo changes.A likely explanation for the dominance of the albedo effect is given by additional regional soil formation in subaerial shelf seas that dampens the overall reduction in total water-holding field capacities of soils (Fig. 8b).Instead, the additional mid-Holocene temperature rise that is related to a reduced effective longwave emissivity as a consequence of soil processes is stronger than the temperature effect which is associated with planetary albedo changes.Especially in the low latitudes, dynamic soils may even reverse a conservative regional cooling trend towards a warming signal via radiative effects of increased atmospheric water vapour as shown by changes in the effective longwave emissivity (Fig. 12b).

Discussion
The design of the soil scheme, similar to equilibrium terrestrial vegetation models (e.g.Prentice et al., 1992), does not account for pedogenesis (soil evolution) and rather reflects the potential soil feedback to climate.As a consequence, the soil scheme may overestimate the soil evolution in regions such North America (not shown) in which present soil formation still lags behind as a consequence of Northern Hemisphere ice sheet retreat during the last deglaciation 19-7 kyr ago (Fairbanks et al., 1992;Harden et al., 1992).As defined by the nature of equilibrium models, the soil scheme calculates a specific soil type which is expected under steady-state climate conditions.The global SAT warming of PI_sol with respect to the static soil properties of PI_ctl may be related to the modelled steady state of the soil properties.Especially in the north-eastern part of North America, prescribed soil properties (PI_ctl) as derived by the MODIS satellite measurement product (Rechid et al., 2009) are characterised by relatively high soil albedos and low total water-holding field capacities (Hagemann et al., 1999).This is typical for soils at an early stage of development.Harden et al. (1992) show that over the past 18 kyr, the age and development of exposed soils in North America have been directly controlled by the Laurentide ice sheet retreat.Such time-dependent lagging processes are not captured by the present equilibrium soil scheme; therefore, anomalous SAT warming in northern North America reflects the temperature response of more developed soil types.Nevertheless, as shown for the mid-Holocene model runs, the soil scheme is not sensitive at locations which have been covered by the Laurentide ice sheet (Fig. 4b).
Given the timescales of various pedogenic processes, specific soil processes can be categorised as rapid (10 1−2 years, e.g.litter formation), medium-rate (10 3−4 years, e.g.mollic, umbric humification), and slow (10 5−6 years, e.g.ferrallitisation, saprolitisation) (Targulian and Krasilnikov, 2007).Zones in which soils are associated with a low climatic potential of pedogenesis (e.g. in cold and hot deserts) reveal rapid pedogenic processes and are considered to be young soils.More developed soils can be found in humid boreal and temperate regions, whereas the humid tropics with the most developed soils have the highest climate potential of pedogenesis and are characterised by slow specific pedogenic processes (Targulian and Krasilnikov, 2007).
So far, the scheme does not account for the potential evolution of soil types, and therefore any transient effects, including lag effects, are currently not considered in the simulated climate-soil feedback.Future transient GCM studies may utilise more sophisticated soil schemes, which should have the aim of deriving time-dependent variables in terms of nonlinear progressive and regressive soil development, acting on different timescales (Johnson et al., 1990;Hoosbeek and Bryant, 1992).

Dynamic soil feedback to the mid-Holocene climate
Compared to the preindustrial period, the mid-Holocene climate is predominantly forced by seasonal insolation changes, with insolation governing warming in high latitudes and cooling in the low latitudes.So far, most of the ESMs cannot reconcile the amplitude of regional temperature changes as suggested by the geological data, raising the question of lacking or impeded feedbacks in present state-of-the-art ESMs (O'ishi and Abe-Ouchi, 2011; Hargreaves et al., 2013;Lohmann et al., 2013).When including the dynamics of physical soil characteristics, our results indicate a warmer mid-Holocene climate than shown in most previous model studies with dynamic vegetation (Gallimore et al., 2005;Braconnot et al., 2007;Otto et al., 2009a, b).Our mid-Holocene study with incorporated dynamic soils exhibits a global SAT rise of 0.34 • C, similar to the study of O'ishi and Abe-Ouchi (2011) showing a global SAT anomaly of +0.36 • C by applying a general climate model with slabocean-atmosphere-vegetation dynamics (dynamic soils not included).O'ishi and Abe-Ouchi (2011) find a mid-Holocene warming in their model which is mainly attributed to vegetation changes (+0.23 • C).However, our model simulations (with vegetation and without dynamic soils) result in a moderate SAT rise (HOL_ctl -PI_ctl: +0.10 • C), whereas the synergy analysis suggests a relatively strong soil feedback (+0.24 • C) that leads to most of the mid-Holocene warming (+71 % of the overall temperature change).Further, the EBM analysis suggests that mid-Holocene warming (HOL_ctl; Wei et al., 2012) is likely attributed to a decrease in the effective longwave emissivity (Fig. 12) as a result of increased water vapour in the atmosphere (+0.28 kg m −2 ).However, the additional consideration of dynamic soils in the model can significantly corroborate the global temperature response by changes in both the planetary albedo and the effective longwave emissivity (Fig. 12).Sundqvist et al. (2010a, b) reconstruct the mid-Holocene climate north of 60 • N as being 2 • C warmer than the preindustrial climate, with seasonal variations of +1 • C in summer and +1.7 • C in winter.They conclude that the year-round warming in the high northern latitudes is mainly attributed to warming during spring and/or autumn.We find that north of 60 • N, polar amplification leads to +1.15 and +1.83 • C of warming without and with dynamic soils, respectively.The highest seasonal warming peaks occur in winter (+2.37 • C) and autumn (+2.35 • C), which can be attributed to changes in insolation; however, these changes are much stronger in magnitude than other studies have estimated (e.g.Sundqvist et al., 2010a, b).Due to a lack of spatial sampling heterogeneity in proxy sample availability (a majority of the data are found in terrestrial European climate archives; Otto et al., 2009a, b), the multiproxy approach of Sundqvist et al. (2010a, b) may be biased (Sundqvist et al., 2010a, b;Otto et al., 2009a, b).Nevertheless, Holocene winter warming in high latitudes seems to be essential for the northward migration of trees and the expansion of temperate deciduous forests in Europe by reducing the limiting effect of the temperature of the coldest month (Prentice et al., 2000).
Climate reconstructions of the mid-Holocene that encompass the high northern latitudes indicate an anomalous spring warming (Sundqvist et al., 2010a, b), though insolation during this season was reduced.Otto et al. (2009a, b) address this discrepancy by means of an atmospherevegetation model study, in which they test the sensitivity of increased forest canopy masking snow cover.The atmosphere-vegetation feedback reveals a spring warming of 0.34 • C together with a forest cover expansion of 13 %, which helps to reduce part of the model-proxy data mismatch (Otto et al., 2009a, b).Despite lower spring insolation, HOL_sol warms up by about 0.98 • C, strongly amplified by the soil feedback (+0.72 • C).HOL_sol exhibits an increase in forest cover (+11%) in high latitudes (60-90 • N) but does not show strong forest cover changes (+0.7 %) in midlatitudes (35-60 • N), as also indicated by the pollen record (Prentice et al., 2000).Midlatitude temperature anomalies during spring are close to PI_sol (−0.07 • C), compensated for by the soil feedback (+0.4 • C) and reducing the effect of low insolation forcing during spring.
In principle, physical soil characteristics of the mid-Holocene climate provide year-round elevated soil moisture, as a first-order effect of increased total water-holding field capacities of soils, and lower soil albedo as compared to the preindustrial period.However, depending on the seasonal variation, both soil characteristics may respond differently: for example, the local insolation increase at the transition from boreal winter to spring triggers increased snowmelt due to the presence of relatively dark soils in the high northern latitudes (60-90 • N), which accounts for an amplified spring warming.The model simulations show that the effect of darker soils, in particular, increases spring snowmelt (+0.15 mm day −1 ) in the high northern latitudes (60-90 • N) during the Holocene as compared to the standard model study without dynamic soils (+0.05 mm day −1 ).
In addition to this, the general effect of increased total water-holding field capacity of soils in the midlatitudes (35-60 • N) develops in parallel with increased water vapour in the atmosphere as a result of higher soil moisture throughout the year.Nevertheless, especially for the comparison of the mid-Holocene and preindustrial model simulations with dynamic soils (HOL_sol -PI_sol), atmospheric water vapour shows a lesser decrease during spring (−0.07 kg m −2 ) and a stronger summer increase (+3.04 kg m −2 ) compared to the spring (−0.44 kg m −2 ) and summer response (+2.22 kg m −2 ) of the standard model simulations (HOL_ctl -PI_ctl).This response reflects the interaction between higher availability of soil moisture, evaporation processes, and anomalies in seasonal insolation changes in the mid-Holocene as compared to the preindustrial period.For the mid-Holocene, elevated soil moisture provides additional evaporation fluxes, which in turn increases water vapour in the atmosphere.This process counterbalances lower local spring insolation as compared to the preindustrial period by reducing the effective longwave emissivity (see EBM analysis, Sect.3.4).
Though spring insolation is reduced compared to the preindustrial period, the soil and vegetation feedback (Wohlfahrt et al., 2004) potentially provide a fundamental mechanism to explain this anomalous spring warming (Sundqvist et al., 2010a, b;Otto et al., 2009a, b).

Dynamic soil feedback to the Last Glacial Maximum climate
In comparison to the coupled GCM studies of PMIP2, the overall SAT cooling in this study (−7.03  2005) quantified the single effects of lowering atmospheric CO 2 , ice sheets, and vegetation dynamics for the LGM time period.They constrained additional global mean cooling to 0.6 • C with a pronounced response in the high latitudes, which is induced by the consideration of vegetation dynamics.The amplified polar cooling, as seen in their model, can be partly inferred from a southward shift of the boreal treeline, i.e. the taiga-tundra transition, which accounts for 2 • C of cooling over Eurasia.More specifically, the southward migration of the taiga-tundra transition reinforces an initial cooling by replacing comparably low-albedo forests with high-albedo grass cover.Especially during the spring season, the retreat of forest cover exposes the snowcovered ground that results in a late start of the snowmelt as a consequence of the radiative effects (e.g.Jahn et al., 2005;Brovkin et al., 2003).The climate response of the taiga-tundra feedback largely agrees with pollen-based biome reconstructions for the LGM; this requires strong amplified high-latitude cooling in order to infer an extensive equatorward regression of the boreal treeline to ∼ 55 • N (e.g.Prentice et al., 2000;Bigelow et al., 2003;Kaplan et al., 2003).Referring to the combination of vegetation and soil dynamics within our model simulations, an overall decline in forest cover is shown in the mid-to-high northern latitudes as a response of amplified cooling (Fig. 11a).The soil feedback analysis suggests the interplay between soil processes and the taiga-tundra feedback that results in an additional southward retreat of the boreal forest cover from ∼ 61.2 to ∼ 57.5 • N (Fig. 11b).Moreover, the temperature response to changes in physical soil characteristics north of the taiga-tundra transition also reinforces the replacement of cold deserts at the expense of grass cover (Fig. 11d, f).As suggested by the EBM analysis, the dominant effect of soil processes, cooling the mid-to-high latitudes, is the result of (soil) albedo changes.However, reduced atmospheric water vapour that is associated with changes in the total water-holding field capacity of soils is also important (see EBM analysis, Sect.3.4).The additive interaction of climate feedbacks, e.g.vegetation dynamics, with dynamic soil changes in the high northern latitudes provides part of an explanation for an increased amplification of the temperature response (more than −3 • C).This is stronger as compared to the effect of vegetation dynamics without soil dynamics (−2 • C), shown by Jahn et al. (2005).The integrative and mutual effect of, in particular, sea-ice, snow, vegetation, and soil changes may provide sufficient cooling in the high northern latitudes to explain the extensive southward regression of the boreal treeline that has been reconstructed for the LGM (Prentice et al., 2000;Bigelow et al., 2003).

Conclusions
This study highlights the importance of soil feedbacks within climate models in order to capture the sensitivity of past climates as suggested by the proxy record.We developed a soil scheme simplified by the assumption that vegetation actively transforms and regulates the abiotic environment (Lovelock, 1979).Therefore, the physical soil characteristics are tightly linked to bioclimatic factors like vegetation and thus climate change.The soil scheme is asynchronously coupled to a comprehensive climate model with vegetation dynamics.The scheme is tested for time slices of the Last Glacial Maximum and mid-Holocene.We find that dynamic soils reinforce the original climate signal for climates warmer and colder than the preindustrial period with a significant amplification in the mid-to-high northern latitudes.This could -at least partly -reconcile the overly weak response in climate models to external forcing during the mid-Holocene (Braconnot et al., 2012;Hargreaves et al., 2013;Lohmann et al., 2013).
Various studies provide ample evidence that terrestrial vegetation dynamics positively feed back to climate changes.As shown herein, the consideration of dynamic soils on climate change undergoes the same trend of direction.Claussen et al. (2006) show that only the fully integrated interaction of atmosphere, ocean, and vegetation dynamics provides the strongest amplitude of climate variation by precessional forcing.In addition to and by analogy with vegetation dynamics (Claussen, 2009), we have shown that the soil feedback may also reinforce the climate response during the late Quaternary.
For the pre-Quaternary, climate reconstructions reveal particularly strong warming in the polar regions (e.g.Jenkyns et al., 2004;Moran et al., 2006;Huber and Caballero, 2011).The simulation of this polar amplification is a key challenge for current model approaches (e.g.Knorr et al., 2011;Krapp and Jungclaus, 2011;Fedorov et al., 2013;Salzmann et al., 2013).To provide a potential solution for the model-data discrepancy, numerous studies provoke processes for polar amplification, dealing with energy balance and heat transport changes (e.g.Sloan et al., 1995;Otto-Bliesner and Upchurch, 1997;Sloan and Pollard, 1998).Alternatively, we suggest that the positive climate feedback due to soil dynamics offers a testable alternative contribution towards understanding polar amplification of greenhouse climate states in the geological past.

Figure 2 .
Figure 2. Physical soil properties (abscissa) used in the dynamical vegetation module JSBACH associated with plant functional types along latitudes.Albedo values are derived from modified satellite measurements of the Moderate Resolution Imaging Spectroradiometer (MODIS;Rechid et al., 2009), and total water-holding field capacity of soils are taken fromHagemann et al. (1999).

Figure 3 .Figure 4 .
Figure 3. Differences in 100-year annual mean of the preindustrial climate state with dynamic soils and fixed soil characteristics (PI_sol -PI_ctl) included for (a) surface air temperature in degrees Celsius, (b) total precipitation in millimetres per day, (c) evaporation in kilograms per square metre per day, (d) forest fraction, (e) grass fraction, and (f) desert fraction.Hatched areas identify non-significant values, calculated by means of a Student t test with 0.05 levels of significance.

Figure 7 .
Figure 7. Vegetation differences in 100-year mean of the mid-Holocene and preindustrial climate state with dynamic soils (HOL_sol -PI_sol) and enclosed soil feedback ( fHOL,soil ) included for (a) forest fraction and (b) soil impact on forest fraction, (c) grass fraction and (d) soil impact on grass fraction, (e) desert fraction and (f) soil impact on desert fraction.

Figure 11 .
Figure 11.As Fig. 7 but for the Last Glacial Maximum.

Figure 12 .
Figure 12.Application of a zero-dimensional energy balance model (EBM) on model studies of the preindustrial period (including dynamic soils, PI_sol), mid-Holocene (HOL) and Last Glacial Maximum (LGM).Panel (a) Global EBM and surface radiation diagnostics that are used for (b) low latitudes (30 • S-30 • N) and (c) high northern latitudes (30-90 • N).Investigations of the EBM and surface radiation diagnostics refer to model results with (blank bars) and without (hatched bars) dynamic soils with respect to the preindustrial period (PI_ctl).The effect of planetary albedo and effective longwave emissivity is expressed in terms of temperature changes.