Simulation of air and ground temperatures in PMIP3/CMIP5 last millennium simulations: implications for climate reconstructions from borehole temperature profiles

For climate models to simulate the continental energy storage of the Earth’s energy budget they must capture the processes that partition energy across the land-atmosphere boundary. We evaluate herein the thermal consequences of these processes as simulated by models in the third phase of the paleoclimate modelling intercomparison project and the fifth phase of the coupled model intercomparison project (PMIP3/CMIP5). We examine air and ground temperature tracking at decadal and centennial time-scales within PMIP3 last-millennium simulations concatenated to historical simulations from the CMIP5 archive. We find a strong coupling between air and ground temperatures during the summer from 850 to 2005 CE. During the winter, the insulating effect of snow and latent heat exchanges produce a decoupling between the two temperatures in the northern high latitudes. Additionally, we use the simulated ground surface temperatures as an upper boundary condition to drive a one-dimensional conductive model in order to derive synthetic temperature-depth profiles for each PMIP3/CMIP5 simulation. Inversion of these subsurface profiles yields temperature trends that retain the low-frequency variations in surface air temperatures over the last millennium for all the PMIP3/CMIP5 simulations regardless of the presence of seasonal decoupling in the simulations. These results demonstrate the robustness of surface temperature reconstructions from terrestrial borehole data and their interpretation as indicators of past surface air temperature trends and continental energy storage.


Introduction
Recent increases in global surface temperatures (Field et al 2014) are the result of the Earth's energy imbalance arising mainly from an increase in the atmospheric concentration of greenhouse gases (Hansen et al 2011, von Schuckmann et al 2016. These changes affect all of Earth's climate subsystems: atmosphere, ocean, land and cryosphere (Levitus et al 2001, Huang 2006, Hansen et al 2010, Lyman et al 2010, Zwally et al 2011. Although most of the Earth's recent energy gains have been stored in the ocean, the land surface has stored the second largest amount and is a critical factor in the evolution of processes such as permafrost melting, which is potentially an important positive feedback within the climate system (Knorr et al 2005, Lawrence et al 2008. Correctly simulating the energy storage in each climate subsystem, including the continental component of the Earth's energy budget, is therefore important for general circulation models (GCMs) to produce robust future climate change projections under different emission scenarios, and to assess the future evolution of potentially important climate feedback processes in the shallow subsurface (Hansen et al 2005, Koven et al 2013.
Estimates of heat gain in the continental subsurface prior to widespread instrumental observations are those obtained from geothermal data by direct interpretation of the data at the local scale. At global scales, the subsurface heat gain is indirectly inferred from ground surface temperature (GST) histories from borehole temperature profiles. Interpretations of GST histories from geothermal data as indicators of climate change nevertheless assume a long-term thermal coupling between the lower atmosphere and the continental subsurface (Beltrami 2002. Multiple studies have shown evidence that supports this assumption at decades and longer time scales (González-Rouco et al 2003, Chapman et al 2004. In contrast, some studies have pointed to the decoupling between high-latitude air and ground temperatures in winter, because of long-term changes in snow cover , Mann and Schmidt 2003, Stieglitz et al 2003. Such impacts could induce a potential bias in surface-temperature reconstructions from geothermal data Huang 2000, González-Rouco et al 2006), but the extent and impact of decoupling on decadal and longer timescales is a matter of debate (González-Rouco et al 2003, 2006, Chapman et al 2004.
In this paper, we use the set of last-millennium (LM) model experiments from the third phase of the paleoclimate modelling intercomparison project (PMIP3) and the fifth phase of the coupled model intercomparison project (CMIP5), to test a fundamental assumption of borehole paleoclimatology: near-surface air and ground temperature changes are strongly coupled over decades to centuries, making inversions of terrestrial borehole temperature profiles representative of surface air temperature changes over such timescales. This assumption is tested across the ensemble of PMIP3/CMIP5 GCM LM simulations, which, in contrast to earlier experiments that used a single model, includes a larger set of state-of-the-art GCMs with different configurations and initial conditions, as well as a variety of land-surface model components responsible for simulating the energy partitioning at the air-ground interface. These simulations additionally include a larger representation of last millennium forcings (Schmidt et al 2011(Schmidt et al , 2012 than used in the GKSS ECHO-g LM simulations (Crowley 2000, González-Rouco et al 2003, 2006. The PMIP3/CMIP5 ensemble therefore includes a broader representation of the effect of different climate forcings on long-term air-ground coupling. For instance, the CMIP5 simulations include the effect of aerosol forcing, a factor that was not taken into account in the ECHO-g past millennium simulations. Our analysis therefore is the first to use a state-of-the-art ensemble of LM climate model simulations to evaluate the effects of the seasonal relationships between air and ground temperatures and the ability of temperature-profile inversion schemes to recover changes in surface air temperatures over decades and centuries.

Data and methods
We use LM PMIP3 simulations (Braconnot et al 2012) and historical CMIP5 (Taylor et al 2012) simulations from the models listed in table 1. The PMIP3 simulations are performed with the same configurations and resolutions as the corresponding GCMs in the CMIP5 experiments, allowing direct evaluation of the paleoclimate simulations in conjunction with the historical runs. We do not use two GCM simulations in the PMIP3/CMIP5 archive: the MIROC-ESM model has a long-term temperature drift (Sueyoshi et al 2013) and the complete ground temperature output for the IPSL-CM5A-LR LM simulation was not available on the Earth System Grid Federation's server.
The temporal interval of our analysis is defined by the LM simulations that extend from 850 to 1850 common era (CE), concatenated to the CMIP5 Historical simulations (Mieville et al 2010), spanning the period from 1850 to 2005 CE. Although LM and historical simulations are not a single continuous simulation, the discontinuity for the variables employed at 1850 CE falls within the range of the simulated climate variability (e.g. Coats et al 2014). We use the monthly surface air temperature at 2 m (SAT, tas in the CMIP5 variable catalogue), as well as the monthly GST, linearly interpolated to a depth of 1 m (GST, tsl in the CMIP5 variable catalogue) in order to have a reference depth for all models (e.g. Koven et al 2013). We also use the deepest subsurface temperature layer (DST) of each model.
In order to verify the methodology for reconstructing GST histories from borehole temperature profiles, we also employ a purely conductive 1D forward model (Beltrami andMareschal 1992, Beltrami et al 1992), which uses the anomalies of the PMIP3/ CMIP5 simulated air and ground temperatures as upper boundary conditions on a semi-infinite homogeneous half space, with a zero-flux bottom boundary condition. This forward model is used to generate synthetic profiles to a depth of 500 m. Subsequently, the synthetic profiles are inverted to retrieve the surface-temperature histories, following the singular value decomposition method of inversion to reconstruct surface-temperature histories from geothermal data .

Results
The mean annual global SAT, continental SAT, GST and DST anomalies, from 850 to 2005 CE with respect to their 1900-1990 CE means are highly correlated at multi-decadal time scales (figure 1 and table 2). Global and continental SAT anomalies (figure 1) display some differences in their respective global means for the five Table 1. General circulation models (GCMs) used in this analysis, the land-surface component of the GCMs, the total simulated subsurface depth, the depth of the deepest subsurface temperature (DST), the number of subsurface layers and the land-surface component references.  simulations. The evolution of these anomalies, however, is highly correlated (>0.9) at decadal to century times scales. If the last century is not taken into account, the correlation values are smaller, but still highly significant (see tables 2 and 3). As expected, the global SAT anomalies are warmer than the continental SAT anomalies for most of the simulation period in all models, due to the well-known differences between ground and ocean thermal properties. Over the continental areas, global GST and DST are warmer than SAT for all the PMIP3/CMIP5 simulations, except for the MPI-ESM-P, in which the continental SAT, GST and DST anomalies are similar throughout the simulation. In the subsurface, the global GST and DST anomalies are identical for the shallow GCMs, showing correlations 0.9 in the BCC-CSM1.1, GISS-E2-R and MPI-ESM-P models, while correlations are slightly reduced (∼0.7) for the CCSM4 (bottom boundary at 43.74 m) and MRI-CGCM3 (bottom boundary at 10 m) simulations. These differences between shallow and deep model temperature anomalies are due to the larger attenuation and phase-shift of the subsurface temperatures at the deepest layer for the deep GCMs in relation to the shallow bottom boundary GCMs (Smerdon and Stieglitz 2006). These results suggest the existence of a strong coupling between GST and SAT over continental areas at time-scales from decades to centuries, although some studies have noted seasonal differences in the airground thermal relationship (González-Rouco et al 2003, Lin et al 2003, Mann and Schmidt 2003, Stieglitz et al 2003, Chapman et al 2004. To study the seasonal effect of winter and summer on the air and ground temperature relationships within the PMIP3/CMIP5 past millennium simulations, we evaluated the continental SAT, GST and DST anomalies in the north hemisphere (NH) extratropical areas (27°N-90°N, 180°W-180°E) in winter (December, January and February) and summer (June, July and August) separately (figure 2). In NH summer, the continental SAT, GST and DST anomalies are similar for all of the simulations, showing a stable temporal coupling and high correlation between air and subsurface temperatures (figure 2, right panel; table 3). Air and subsurface temperature anomalies in NH winter (figure 2, left panel) show some decoupling, with GST and DST warmer than continental SAT for four of the simulations. The one exception is the MPI-ESM-P simulation, in which the coupling remains nearly constant in both seasons. Winter decoupling across the ensemble of models is reflected in the correlation coefficients; all the models show lower correlation coefficients between continental SAT and GST in NH winter than in NH summer, except the MPI-ESM-P model that shows correlations of 0.9 independent of the season.
To explore the possible causes of winter decoupling within the PMIP3/CMIP5 simulations we compute the correlation coefficients in each grid-point between SAT and GST from 850 to 2005 CE seasonly for the globe. The correlation in NH summer (June, July and August, annual mean) between SAT and GST (figure 3, right panels) shows high correlation coefficients globally, except in areas where freezing phenomena persists during the summer. This behaviour appears in four of the GCMs, but the MPI-ESM-P simulation yields correlation coefficients that are high globally. In the NH winter (December, January and February, annual mean) (figure 3, left panels), the correlation coefficients between SAT and GST decrease for the five simulations in the northern hemisphere. High correlation coefficients remain in the southern hemisphere. For the BCC-CSM1.1, GISS-E2-R, MRI-CGCM3 and CCSM4 simulations, the correlation coefficients decrease to less than 0.5 in part of North America and Eurasia, where the presence of snow and Table 3. Seasonal temporal correlation coefficients between temperatures in NH extratropical areas (27°N −90°N, 180°W −180°E), from 850 to 2005 CE (from 850 to 1900 CE in brackets) for the annual and the 21-year filtered time series. Employed variables are the SAT over the continental areas (Continental SAT), ground surface temperature at 1 m (GST) and annual subsurface temperature at the deepest layer (DST) from each PMIP3/CMIP5 simulation. All of the correlation coefficients are significant 95% level using a phase-randomizing bootstrapping technique with 1000 Monte Carlo runs (Ebisuzaki 1997 freezing phenomena are expected to cause decoupling on annual timescales (e.g. Kane et al 2001, Sokratov and Barry 2002, Beltrami and Kellman 2003, Stieglitz et al 2003, Smerdon and Stieglitz 2006. Although the correlation coefficients in NH winter decrease with respect to those in NH summer in the MPI-ESM-P simulation, this decrease is weak and the correlation coefficients remain high globally in both seasons. The decrease of the correlation coefficients in snow covered areas for the PMIP3/ CMIP5 simulations reinforces the importance of the insulating effect of snow cover on air-ground coupling. The strong air-ground coupling in NH winter and NH summer for the MPI-ESM-P simulation may be due to the use of the same value of the thermal conductivity for frozen and thawed subsurface, as well as  4). The 21-year filtered spatial correlation patterns remain similar to the annual correlation patterns, although the correlation coefficients are universally larger for the filtered quantities. This behaviour is present for the five model simulations, indicating that at low frequencies the coupling between air and ground temperatures is enhanced. The correlation coefficients also display spatial variability over areas where snow cover is not present Figure 2. The 21-year running mean of simulated anomalies relative to 1900-1990 CE in NH extratropical areas (27°N-90°N, 180°W -180°E), using continental SAT, GST and DST in NH winter (December, January and February annual mean, left column) and NH summer (June, July and August annual mean, right column).
(figures 3 and 4), likely due to factors such as the spatial variation of moisture (Dirmeyer 2011) and vegetation cover, the latter of which increases latent heat fluxes through evapotranspiration (Bonan 2001). Local correlation coefficients between SAT and GST are enhanced through the tropics in the MRI and GISS models, but areas of low vegetation such as the Sahara desert also have large correlation coefficients in multiple models during both the NH summer and winter. Multiple non-cryogenic processes are therefore relevant to the coupling between SAT and GST temperatures that are simulated differently across the model ensemble. It is beyond the scope of this investigation to diagnose each of those processes and their impact on the coupling between air and ground temperatures, rather it is sufficient to note that they impact the coupling by differing degrees within each of the models and have different temporal signatures based on the annual and 21-year filtered correlation patterns.
Regardless of the underlying seasonal processes that cause short-term coupling variations, we can asses the possibility that they may alter the surfacetemperature histories reconstructed from geothermal data, by generating synthetic subsurface temperature anomaly profiles from the PMIP3/CMIP5 continental SAT, GST and DST anomalies. We do so using an idealized semi-infinite homogenous subsurface model as described in the data and methods section to generate profiles at all locations on the continents. Although the network of borehole measurements samples a subset of global land grids, González-Rouco et al (2006) demonstrated that the network sampling was sufficient for estimating results from a complete representation of the global land grid. We therefore use a complete sampling of continental grid points as representative of results from a subsampling of grid cells based on the actual distribution of the borehole network. The synthetic temperature anomaly profiles are then inverted to obtain the corresponding ground surfacetemperature histories, following standard borehole climatology methods for reconstructing the past millennium temperature changes (e.g. Mareschal and Beltrami 1992). The synthetic subsurface temperature profiles (figure 5) derived from the SAT, GST and the DST anomalies as upper boundary conditions, show similar profiles across the five simulations. In all cases, each profile shows a different subsurface temperature anomaly in the upper 100 m, as a response to each time history of the simulated SAT anomalies. Most importantly, however, the GST histories recovered by inversion of each subsurface temperature anomaly profile retain the low-frequency changes in the air and subsurface temperatures of each model simulation. These results support the ability of borehole climatology methods to reconstruct surface-temperature histories from geothermal data, in spite of seasonal decoupling processes. With regard to snow cover specifically, the absence of a low-frequency influence is likely due to the fact that simulated changes in the snow cover thickness and length of snow cover season do not change significantly over the simulation period (Goodrich 1982).

Conclusions
Results support the assumption that air-ground thermal coupling is strong on global and regional spatial scales and over decades and centuries within the PMIP3/CMIP5 last millennium simulations. In NH winter particularly, the PMIP3/CMIP5 simulations show a short-term decoupling in the northern high latitudes, due to the insulating effect of snow cover. Because these seasonal variations of temperature do not penetrate more than a few metres into the subsurface, however, the short-term decoupling does not affect the GST histories retrieved from geothermal data. Furthermore, inversions recover the GST histories from the temperature trends of the PMIP3/ CMIP5 temperature simulations.
The above conclusions are consistent with those of (González- Rouco et al 2003Rouco et al , 2006) despite significant differences between the earlier GKSS ECHO-g model and the state-of-the-art GCMs from the PMIP3/CMIP5 ensemble used herein. In particular, these conclusions stand across a collection of models with different landsurface model components. The earlier ECHO-g simulations also use different external forcings and did not include aerosol forcings. Despite each of these differences and the diversity of models across the PMIP3/CMIP5 ensemble, the thermal coupling between the lower atmosphere and the continental subsurface is maintained over decades to centuries within the five PMIP3/ CMIP5 GCM simulations, validating the interpretation of underground temperatures as indicators of climatic change on decades to centuries. These findings also add to the now large body of work supporting the use of borehole reconstructions in the ensemble of proxy-estimated temperatures of the last millennium and the importance of the borehole estimates as benchmarks for low-frequency temperature change over the last 500-1000 years.