Influence of radiative forcing factors on ground–air temperature coupling during the last millennium: implications for borehole climatology

Past climate variations may be uncovered via reconstruction methods that use proxy data as predictors. Among them, borehole reconstruction is a well-established technique to recover the long-term past surface air temperature (SAT) evolution. It is based on the assumption that SAT changes are strongly coupled to ground surface temperature (GST) changes and transferred to the subsurface by thermal conduction. We evaluate the SAT–GST coupling during the last millennium (LM) using simulations from the Community Earth System Model LM Ensemble (CESM-LME). The validity of such a premise is explored by analyzing the structure of the SAT–GST covariance during the LM and also by investigating the evolution of the long-term SAT–GST relationship. The multiple and single-forcing simulations in the CESM-LME are used to analyze the SAT–GST relationship within different regions and spatial scales and to derive the influence of the different forcing factors on producing feedback mechanisms that alter the energy balance at the surface. The results indicate that SAT–GST coupling is strong at global and above multi-decadal timescales in CESM-LME, although a relatively small variation in the long-term SAT– GST relationship is also represented. However, at a global scale such variation does not significantly impact the SAT– GST coupling, at local to regional scales this relationship experiences considerable long-term changes mostly after the end of the 19th century. Land use land cover changes are the main driver for locally and regionally decoupling SAT and GST, as they modify the land surface properties such as albedo, surface roughness and hydrology, which in turn modifies the energy fluxes at the surface. Snow cover feedbacks due to the influence of other external forcing are also important for corrupting the long-term SAT–GST coupling. Our findings suggest that such local and regional SAT–GST decoupling processes may represent a source of bias for SAT reconstructions from borehole measurement, since the thermal signature imprinted in the subsurface over the affected regions is not fully representative of the long-term SAT variations.

The previous analyses support estimations of global/hemispheric past temperature obtained from borehole temperature inversions. First, they support the overall performance of the methodology in retrieving past GST histories from borehole profiles. 5 Second, the use of up-to-date forcing representations in CMIP5 model ensembles also ensures that long-term alterations of surface properties like those induced by LULC changes, the effect of anthropogenic aerosols cooling and the potential long-term snow cover feedbacks induced by both forcings, do not seem to bias inversion results at global/hemispheric scales. In spite of these positive results, the analysis of last generation model experiments including the complete set of agreed CMIP5 forcings (Schmidt et al., 2011(Schmidt et al., , 2012 in the so-called all-forcing experiments does not allow for insights into the individual effect of 10 each of the forcing and the related feedbacks, both at global/hemispherical or at continental/regional scales. This would be desirable in order to quantify the effect of the individual forcings on SAT-GST coupling, and to obtain estimates of its particular temporal evolution that would allow for disentangling its contribution from that of other forcings. Therefore, the present work considers all-forcings and single-forcing types of experiments in order to address how state-of-the-art climate models simulate SAT-GST coupling from global to regional scales and to evaluate the potential influence of the external forcings on 15 the SAT-GST relationship. This will also provide information about where and when a decoupling of SAT-GST may exist, with implications for borehole inversion practices at different spatio-temporal scales. For this purpose we use the Community Earth System Model-Last Millennium Ensemble (CESM-LME; Otto-Bliesner et al., 2016), which to date is the largest existing ensemble of LM simulations with a single model and has not been used in previous assessments of this kind. The CESM-LME includes all-and single-forcing experiments that consider jointly or separately the transient evolution of solar variability, vol-20 canic activity, orbital changes, GHGs, anthropogenic aerosols and LULC (see details in Sect. 2) Therefore, using all-forcing simulations (ALL-F hereafter) allows for evaluating the SAT-GST relationship throughout the LM in a realistic representation of the real world conditions. Additionally, the single-forcing experiments are suitable to identify the specific role that each forcing might play and its 'fingerprint' (Hegerl et al., 2011) on the presence of biases in the SAT-GST coupling. Most external natural (solar and volcanic variability) and anthropogenic (GHGs, LULC and aerosols) have the poten-25 tial to indirectly affect SAT-GST thermodynamics through snow cover feedbacks (Hartmann et al., 2013). Additionally, LULC changes can also have direct influence in SAT-GST. During the LM the Earth's land cover has been substantially modified by the replacement of natural ecosystems with agricultural land, especially since the industrial period (Pontgratz et al., 2008;Hurtt et al., 2009). Such land use changes have the potential of altering the surface energy balance by the modification of energy fluxes and moisture budgets (Betts, 2001;Betts et al., 2007;Brovkin et al., 2006;Anderson et al., 2011;Zhao and Jackson, 30 2014) driving to direct impacts on the atmosphere and ground-surface temperature relationship (MacDougall and Beltrami, 2017). Furthermore, vegetation changes may additionally lead to indirect non-linear effects that are relevant in the borehole context. The most important one is the fact that deforestation at high latitudes leads to an increase in snow cover since low vegetation accumulates continuous snow cover more readily than forests in early winter favoring its permanence longer in spring (Myhre et al., 2013). The local, regional and large scale implications of these interactions on long-term climate variability and 35 used in Landrum et al. (2012). It is composed by a total of 30 simulations including a subset of 10 simulations that incorporate all the LM external forcings and smaller subsets of single-forcing simulations considering each forcing individually (see Table   2 for details and references therein). The simulations used in this study and how they trace back to the original experiment names produced by Otto-Bliesner et al. (2016) are shown in Table 3.

3 Methods
In the present work the LM refers to the period from 850 to 2005 CE while the periods from 850 to 1850 and 1851 to 2005 CE refer to pre-industrial and industrial, respectively. 2 m air temperature is used for SAT and the first land model level (Table   1) represents GST. Additionally, the ST at different depths are used for addressing different aspects of the subsurface heat transport within the CESM-LME. The simulated ST does not include data over the Antarctic region therefore this region is 10 excluded from the analysis.
The relationship between SAT and GST is analyzed from two complementary perspectives. First, the covariance structure during the LM of the global SAT and GST anomalies relative to 1851-2005 CE is assessed at interannual scale. Likewise, the spatial patterns of the correlation between SAT and GST are also analyzed. This spatial analysis is extended to the differences between SAT and GST LM mean values as this gives an additional measure of the energy exchange across the air-ground inter- 15 face (Bartlett et al., 2004) that may inform about a potential air-ground temperature decoupling. In addition, the same spatial analysis is done for the differences between GST and ST L8 in order to get insights into processes affecting the pure conductive subsurface transport assumption within the shallow subsurface. For specific cases where SAT and GST exhibit some kind of decoupling, a description of the main processes that lead to this decoupling at local and regional scales is presented. For this purpose, the evolution of SAT, GST and ST at different depths is assessed during the last 105 years of the simulation period 20 (1900( -2005 for the cases of interest. In some particular examples additional variables (e.g. latent and sensible heat fluxes or snow cover) are also included for a better description of the physical processes. The purpose of such analyses is providing some examples rather than detailing the minutiae of the processes that lead to different SAT-GST responses at many different locations. For the analysis described above the ALL-F ensemble will be used, thus, providing from the point of view of the forcings considered a more complete and realistic representation of real world conditions than the single forcing runs. This 25 allows for using the surrogate reality of the CEMS-LME model as a test bed for detecting potential sources of deviations in the SAT-GST relationship during the LM.
Second, we focus on the assessment of the long-term trend in the SAT-GST relationship. Non-stationarity in such association at long timescales (multi-decadal to multi-centennial) might imply that GST variations would not be representative of the SAT variations, thus rendering unreliable inferences of past climate change (Bartlett et al., 2004). A two-phase regression model 30 (Solow, 1987) is applied to address this issue, that allows for detecting changes in the long-term trend of the SAT-GST relationship. The year-of-change is identified, as well as the trend before and after the change. The timing of change and the magnitude of trends are suggestive of long-term changes in surface coupling. Initially, this assessment is performed at a global scale by 6 Clim. Past Discuss., https://doi.org/10.5194/cp-2018-72 Manuscript under review for journal Clim. Past Discussion started: 4 July 2018 c Author(s) 2018. CC BY 4.0 License.
analyzing the LM evolution of global SAT minus GST anomalies. The analysis of the long-term trend is then extended to the spatial domain starting from the evaluation of SAT and GST linear trends independently during the industrial period since the largest changes are expected to occur within this period. In this part besides the ALL-F ensemble we have also considered the GHG-only, LULC-only and OZ/AER-only ensembles because they have the highest potential of altering the land surface characteristics. This allows for identifying the influence of each forcing on the LM SAT-GST coupling, including changes in 5 snow cover, soil moisture, latent and sensible heat variations. The linear trend analysis of SAT and GST gives a general view of their evolution through the industrial times and a first glance at the most important areas where long-term SAT-GST decoupling may exist. Finally, considering that the timing of variations in the land surface properties may have started before the industrial period, the two-phase regression analysis is also applied to the complete LM.

Results
CESM-LME supports the assumption that SAT is tightly coupled with GST at global scales and above multi-decadal scales. the low pass filter influence of the heat conduction below the surface, whereas for SAT, GST and ST L8 31-yr running mean low pass filter outputs are shown. Subsurface temperature anomalies closely track SAT anomalies with relative small differences between them, thus indicating that air and soil temperature are coupled above multi-decadal time scales. The correlation coefficients for the 31-yr filtered series in Table 4 indicate high correlation (p<0.05) for the soil layers close to the surface that diminishes slightly with depth due to the phase shift of the signal. Table 4 also indicates that the correlation considering the 20 high-frequency variations (yearly; left column) is only high at the levels close to the surface, whereas at the deepest layer the correlation is low since the high-frequency variations are progressively filtered out and phase shifted as depth increases.
Despite the strong coupling between air and subsurface temperatures at global scales, the existence of a relative small offset between SAT and GST that grows backwards in time is evident for the annual (Fig. 1a) averages. This indicates a slight longterm decoupling between SAT and GST. Previous works have argued that the nature of this offset arises from the changes in 25 snow cover and its influence insulating the ground from very low air temperatures (e.g. Mann and Schmidt, 2003;Bartlett et al., 2004). González-Rouco et al. (2009) suggested that beside snow cover, soil moisture has also the potential of accounting for the thermal difference between SAT and GST. In order to explore the influence of such processes on the global SAT-GST relationship this analysis is extended considering winter and summer seasons (DJF and JJA hereafter) independently (Fig. 1c,e).
The offset between SAT and GST observed in the annual plot is apparent only in the DJF season. Indeed, the differences are 30 larger during this season than in the annual data whereas the JJA evolution of SAT and GST anomalies is virtually identical.
In addition, the correlation coefficients (Table 4) suggest slightly lower correlation for DJF than for JJA, even if SAT and GST are highly correlated and significant for both seasons, especially for the filtered values. The largest offset occurring in DJF 7 Clim. Past Discuss., https://doi.org/10.5194/cp-2018-72 Manuscript under review for journal Clim. Past Discussion started: 4 July 2018 c Author(s) 2018. CC BY 4.0 License. reinforces also the important role of snow cover in decoupling air and soil temperatures. Figure 1b,d,f further illustrates the strong relation between SAT-GST offset and the snow cover by displaying the 31-yr low pass filter outputs of the LM evolution of global snow cover and SAT-GST differences. Note that on the NH winter (DJF) and, due to its influence, also on annual averages, the decrease (increase) in snow cover lead to a decrease (increase) in the SAT-GST offset. This is due to the insulating effect of snow that keeps GST close to zero while SAT can reach large negative values. Thus, an increase of snow cover leads 5 to larger negative SAT-GST differences. For JJA on the contrary, an in-phase relationship is found at all time scales. Long-term trends change both in snow cover and in SAT-GST after the end of the 19th century. During the boreal summer increases of snow also enhance SAT-GST differences due to insulation of the ground from the warmer summer SAT; and conversely for snow cover decreases. This effect dominates the global average over that of the JJA austral winter during which SAT-GST and snow cover changes experience an anti-phase relationship as described above. Therefore, the NH influences anti-phase 10 covariability of snow cover and SAT-GST during DJF (detrended correlations, r=-0.52, p<0.05) and annual (r=-0.67, p<0.05) and in-phase covariability during JJA (r=0.62, p<0.05).
Dashed lines in Fig. 1b,d,f indicate the long-term trend in both SAT-GST offset and snow cover. Interestingly, the SAT-GST offset remains constant within the pre-industrial period but its long-term trend experienced a change during the industrial period. The behavior is similar for the annual and seasonal cases with some difference in the timing. Such variations in the 15 long-term relationship should be taken into account since they play an important role in decoupling SAT and GST at time scales relevant for the borehole theory. Nonetheless, this relative small change in the long-term trend has limited impact for the global SAT-GST coupling (Sect. 4.1).
The spatial variability of the relationship between SAT and GST gives further insights about the role of different processes on the SAT-GST coupling. Figure 2 shows the spatial distribution of the differences between mean LM SAT and GST as well 20 as the correlation coefficients for annual, DJF and JJA averages in the ALL-F ensemble. Maps have been obtained from the ALL-F 2 member; dots indicate that 80% of the ensemble members agree in delivering significant differences or correlation for a given grid-point. In the case of the annual temperatures ( Fig. 2a), GST is generally warmer than SAT being the differences low in most of the globe (less than 2 • C) except in the NH mid-and high-latitudes areas where the differences are higher (up to 15 • C). The correlation maps (Fig. 2b) provide a similar pattern with high and significant values in most the globe (>0.8 in 25 regions located below 45 • north) and lower correlations over NH mid-and high-latitudes; especially over the east of Siberia.
A similar behavior is seen in the DJF season although more pronounced. During these months GST is much warmer than SAT reaching differences up to 30 • C at the northernmost part of North America and Eurasia (Fig. 2c). Differences are smaller at mid-and low-latitudes. The influence of the ocean over coastal areas providing larger SAT relative to GST is noticeable.
Likewise, the correlation is lower over northern snow covered areas (Fig.2d) while in the rest of the globe it remains high. On 30 the contrary, Fig. 2e,f shows that during JJA, when the snow cover is scarce, the SAT-GST coupling is strong globally with temperature differences lower than 2 • C, and high correlation coefficients (>0.9). Consequently, the role of the snow cover in decoupling SAT and GST is highlighted. Positive correlation values are low over borderline areas where snow cover is more variable, thus producing variability in the SAT-GST offset and thereby altering the covariance structure. Close to these areas in central Asia the high negative correlation within the Tibetan Plateau is noteworthy (see discussion below).  Figure 3 illustrates the SAT-GST decoupling due to the snow cover at the local scale for a particular grid-point as an example.
The grid-point is located over a region with considerable snow cover during the boreal cold season (northeastern Russia).
The DJF and JJA evolutions of SAT, GST, ST at different depths and the snow cover for the last 105 years from the ALL-F 2 simulation are shown. Note that in DJF the snow covers 100% of the grid-cell during almost the whole period. Thus the soil is insulated and the difference between SAT and GST is~-15 • C on average. The temperature of the deeper layers are presented 5 in order to illustrate the amplitude attenuation and phase shift with depth of the temperature signal. Note that during DJF, GST is only slightly below 0 • C and the agreement of its variations with those of SAT is only noticeable in the largest changes of both while ST L6 and deeper STs are above 0 • C. In JJA, SAT and GST are very similar and their low frequency variability propagates to deeper levels, all above 0 • C.
Some aspects of the SAT-GST spatial distribution deserves further attention. One of the most noteworthy is that over the Ti-10 betan Plateau region temperature differences between SAT and GST are as large as in the mid and high NH latitudes for both annual and DJF (Fig. 2a, c). However SAT and GST are negative correlated (Fig.2b,d), being the only region of the globe where this occurs. For JJA the correlation nonetheless is positive and high (Fig. 2e). The nature of this opposite phase arises from discontinuous snow cover over this region during DJF. Usually the snow cover insulates the soil from the colder air, avoiding the heat exchange with the atmosphere (as shown in Fig. 3). Nevertheless, discontinuous snow cover insulates the soil only 15 partially, leading to this particular SAT-GST interaction. Figure 4 displays this behavior for a grid-point located over the Tibetan Plateau. SAT, GST, ST at different depths, snow cover and surface sensible heat flux (SHFLX) are shown. In DJF during periods of low snow cover the fraction of surface expose to the atmosphere allows for exchanging energy from the warmer soil to the colder air. Conversely, when snow cover is high, the large fraction of insulated soil reduces almost completely the heat transfer from the soil to the atmosphere. Therefore with lower (higher) fractions of snow cover, higher (lower) heat transfer 20 takes place with GST decreasing (increasing) and SAT increasing (decreasing). Indeed, there is a high negative correlation (-0.84, p<0,05) between the snow cover fraction and SHFLX which is the main way for energy dissipation within this region since latent heat fluxes (LHFLX) are negligible. Higher/lower albedo due to variations in snow cover fraction also contribute to the negative SAT-GST correlation over this region (not shown). During JJA on the other hand, snow cover is negligible thus SAT-GST coupling is strong. Note that GST and all ST are above zero and GST is higher than SAT as it is warmed by radiative 25 gain and transferring heat to the atmosphere, thus the positive correlation (0.5, p<0,05) between SAT-GST and SHFLX.
The spatial SAT-GST differences during JJA (Fig. 2e) depict other relevant aspects that have an influence on the SAT-GST relation at relative short-time scales. SAT is generally colder than GST globally however for JJA there are large areas inland, mainly located at the south-east of the United States, some parts of central and eastern Europe and eastern Asia with warmer SAT relative to GST. Variations in LHFLX from DJF to JJA are driving this effect. Figure 5 shows that the areas where LHFLX 30 increase in JJA relative to DJF are related to the same areas where SAT is higher than GST in JJA (Fig. 2e). Therefore there is a direct relation between the increase of evapotranspiration in JJA and the ground temperature response at these locations.
The timeseries in Fig. 5 display this behavior for a grid-point located over southeastern of China. During JJA the surplus of energy due to higher solar radiation reaching the surface is mostly dissipated as latent heat leading to a net heat loss at the ground surface. Note the anti-correlation between surface soil moisture and LHFLX (-0.46; p<0,05) during JJA as well as the large SHFLX values that take place when soil moisture is lowest and evapotranspiration is limited. Therefore, the high rate of LHFLX contribute to cool the surface and GST tends to be lower relative to SAT. Soil water content also exhibits large changes during JJA in consistence with the large evapotranspiration (Fig. 5, top), thus providing the source of moisture that contributes to temperate SAT and cool GST. Large variations in evapotranspiration from DJF to JJA are also present at mid-latitudes of the SH summer continents (America, Africa and Australia) although only a very limited impact on Fig. 2e is perceived; especially 5 over the western coast. Over these regions high incoming energy impinging the surface during SH summer (not shown) supports high rate of latent heat fluxes, however as soil water becomes a limiting factor, more energy is dissipated as sensible heat and the ground surface is warmed. Therefore GST experiences higher temperature than SAT on average. Figure 5 also shows large evapotranspiration over the tropical rain forest in America and Africa, both in DJF and JJA that does not translate to positive SAT-GST differences. Over the rain forest, the energy fluxes at the surface do not vary significantly from DJF to JJA since 10 incoming radiation is relatively constant throughout the year, and transferred to evaporation and evapotranspiration within the canopy while soils are well watered by precipitation to support the large amounts of evapotranspiration. This situation leads to a small range of variation both in SAT and GST with very small differences between them and higher GST relative to SAT. Figure 2a,c,e shows similar higher SAT relative to GST also at some coastal areas both in JJA and in DJF. During DJF the effect is mainly present in the NH mid-latitudes as, for instance, in most of the coasts of Europe, both the east and the west coasts 15 of North America, as well as in Japan and the east coast of China. During JJA on the other hand, this behavior is seen mostly in the SH mid-latitudes, as in southern South America, South Africa and the south of Australia. Therefore there is a clear relationship between the winter season for each of the hemispheres and this SAT-GST temperature response in the CESM-LME simulations. Interestingly, higher SAT relative to GST is also evident at some coasts over tropical regions both in JJA and in DJF mainly over southeastern Asia, the south of the Indian subcontinent, the Gulf of Guinea and some areas of South and 20 Central America. In the CESM-LME the atmospheric grid-box of the coastal areas is partitioned into land and ocean fractions, For the areas with sea ice formation an additional sea ice fraction is considered (Neale et al., 2012). This configuration of the coastal grid-points leads to a partition of the energy fluxes at the surface into those of the land fraction and those of the ocean part. During the cold season such partition is determining the higher SAT warming relative to GST since the relatively low net radiation that impinges the surface at mid-and high-latitudes limits the ground surface heating as well as the energy fluxes out 25 of the land surface fraction. In contrast, the energy fluxes from the ocean surface to the air above are large mainly as a result of the temperature difference between the water and the air above, since the water is relatively warmer than the continental air in winter. Strong winds blowing predominantly onshore over most of the coast also contribute with advection of latent heat from the relatively warmer ocean. The dissipation of energy from the ocean fraction to the atmosphere warms the air so the net effect is the higher SAT relative to GST in the winter season. Over the tropical coasts that exhibit the same behavior the energy 30 fluxes out of the ocean fraction of each grid-point also contribute to the higher warming of the air relative to the ground surface.
Nevertheless, although at these locations high rates of evapotranspiration all year long play also an important role since they generate evaporative cooling of the ground surface as in the example describe in Fig. 5.
The different examples used to illustrate the most important processes that may influence the air and soil temperature relationship at short time scales depict also relevant information about the propagation with depth of the annual cycle. For instance, at Comparable pictures with some differences in the magnitude of gradients can be seen in the previous figures.
Therefore, it is interesting to understand also the propagation below the surface. GCMs simulate purely conductive regimes, and the temperature variations that propagate to deeper soil layers are established at or near the ground surface (Smerdon 5 et al., 2003). Thus it is important to asses the propagation of the temperature signal within the shallow subsurface. This issue is addressed by analyzing the relationship between GST and ST L8 . where frozen ground persists during JJA. Indeed, the nature of the large departure in the temperature response at the shallow subsurface at these locations arises from non-conductive processes related to latent heat release/uptake of freezing and thawing of the water content above 1 m depth that may account for the subsurface heat transfer (Kane. et al., 2001).
To illustrate this mechanism, Fig. 7 shows the temperature evolution of SAT, GST and ST at different depths as well as the soil 20 ice content (SIC) in the upper soil layers for a grid-point located at the north of Canada. Over these areas the SIC in the upper 1 m depth increases (decrease) during the cold (warm) season. During JJA SAT and GST increase/decrease at the same rate since no ice is present at the ground surface so it is warmed by radiative gain and heat is transferred to the atmosphere. However, for deeper soil layers the energy available is employed in melting the SIC (note the lower SIC in L6 during JJA relative to DJF) and latent heat is required so these layers do not experience a temperature increase as the shallowest ones. Therefore the tempera-25 ture at L8 (~1 m depth) is kept below/near 0 • C during the warm season due to the the zero-curtain effect (Outcalt et al., 1990), while GST is centered around 12 • C, leading to differences of ca. 15 • C between GST and ST L8 and low correlation coefficient (0.28 for this grid-point) during JJA. In turn, during DJF SAT distributes around -35 • C and the frozen ground experiences skin temperatures of about -8 • C and of about -4 • C at 1 m depth. Consequently, the temperature offset between GST and ST L8 is largest in summer. According to this, there are some non-conductive processes associated with permanent frozen soils in the 30 shallow subsurface that are included in the CLM4 parametrization (Lawrence et al., 2011) and that play an important role on the heat transport. The higher SIC in L8 during JJA relative to DJF (Fig. 7) arises from the fact that freezing of the active layer begins in late autumn and, due to the release of the latent heat of fusion, the freezing front is inhibited (Hinkel et al., 2001) and only reach L8 in spring when SIC at this layer peaks. Then, as the thaw front penetrates downward in summer, ice in the shallowest soil is melted, but it does not reach L8 until autumn, when SIC at this layer is lowest. Therefore due to the thawing/freezing processes, seasonal changes at the upper and deeper subsurface levels are phase-shifted.

SAT-GST long-term changes
The mechanisms that have been described have an impact on the coupling between SAT and GST at short-time scales but they do not affect the long-term SAT-GST association if they are stationary since its influence would be constant at long time 5 scales. However if such mechanisms experience variations with time, the SAT-GST relationship would also change with time, thus the thermal signature imprinted in the subsurface would not be representative of the long-term SAT variations (Bartlett et al., 2005). Figure 1b illustrates the existence of a constant offset between SAT and GST within the preindustrial period that changes during the industrial times indicating a variation in the long-term SAT-GST relationship. This may be relevant in the interpretation of borehole climate reconstructions since it may induce a long-term decoupling between SAT and GST in the 10 CESM-LME. At global scale the change in the long-term SAT-GST offset have an impact of about 0.05 • C (Fig. 1a,b) and thus, do not seem to be, at these scales, very relevant. However, the impacts could be larger for other GCMs with higher climate sensitivity or a different representation of surface processes that may contribute to decouple GST from SAT (e.g. snow cover).
Likewise, within the CESM-LME simulations impacts on decoupling may be important at regional or local scales.
To examine the spatial distribution of the long-term SAT and GST evolution during industrial times, we evaluate the linear 15 trends of both of them independently during this period at every land model grid-point. Besides the ALL-F ensemble we have also considered in this analysis the anthropogenic single-forcing ensembles (Fig. 8), bearing in mind their potential influence on the processes that modulate the relationship between SAT and GST, as for instance: variations in snow cover, soil moisture and albedo, among others. The results in this section are shown considering information from all members of the ensembles (see details in Table 3 the global warming pattern due to the influence of GHGs (Hartmann et al., 2013). Indeed, if only the contribution of GHGs is considered the warming would be higher and globally distributed. Figure 8 also indicates that the cooling in the ALL-F ensemble is mainly driven by the contribution of the LULC and OZ/AER external forcings. For instance, the cooling trends over the Baltic sea and the north of the Black and Caspian seas that dominate the SAT and GST cooling trends during the industrial period in the ALL-F ensemble are the result of the influence of LULC changes (Fig. 8e,f) with additional contribu-30 tions of OZ/AER (Fig. 8g,h). In addition, the negative trends of both SAT and GST over some areas of Africa as well as over the north east of Brazil are also detectable in the LULC-only ensemble (Fig. 8e,f). Likewise, the OZ/AER-only ensemble also contributes to the cooling over Brazil, and the strong negative trends observed at south-east China in the ALL-F ensemble are clearly identifiable in this ensemble (Fig. 8g,h).
Although the general pattern of cooling/warming during industrial times is broadly similar for SAT and GST, it differs substantially in some regions. Note that there are considerable differences in the amplitude of warming trends in SAT and GST over Fennoscandia as well as at the northernmost part of North America in the ALL-F ensemble (Fig. 8a,b). Similarly, over some areas of central and eastern Europe SAT and GST industrial trends have different sign. There are also considerable differences in the amplitude of the cooling in SAT and GST over northeastern Brazil. Such dissimilar behaviors of SAT and GST during 5 the industrial period are connected to variations in the energy fluxes at the surface in response to changes in the land surface characteristics due to the influence of the external forcings during this period.
For a more detailed analysis of SAT and GST long-term relationship, we applied a two-phase regression model (see Sect. 3) during the complete LM at every land model grid-point to the differences between SAT and GST in the ALL-F ensemble as well as in the GHG-only and LULC-only ensembles (Fig. 9). In the case of the OZ/AER-only ensemble the linear trends of the 10 differences between SAT and GST during the industrial period are presented in Fig. 10 since this set of simulations span solely from 1850 to 2005 CE. Figure 9a,b,c shows the year-of-change of the long-term SAT-GST trends for the ALL-F, GHG-only and LULC-only ensembles respectively. Figure 9d,e,f (9g,h,i) represents the magnitude of the change before (after) it takes place. Note that except for some areas, most continental regions experience statistically significant dates of change in phase with the start of the industrial period, during the 19th century (darker blue). For all of them SAT-GST offset trends before the 15 change are negligible or slightly negative in the ALL-F ensemble (Fig. 9d), indicating a stable relationship between SAT and GST globally before industrial times. This characteristic is also evident in the GHG and LULC-only ensembles (Fig. 9e,f). The strong positive or negative trends before the change (Fig. 9d,e,f) correspond to early dates of change around the 10th and 11th centuries (light green colors in the year-of-change maps) as a result of a change in the simulated differences between SAT and GST early in the LM and a stable behavior in GST-SAT differences thereafter. Note that for all of these regions the trend after 20 the change is negligible and none of them is statistically significant (Fig. 9c,f,i).
The trends after the change on the other hand, indicate significant regional variations in the long-term SAT-GST relationship.
In general annual SAT minus GST yields negative values (as shown in Fig. 2a) with the exception of the coastal areas as explained in Sect. 4. Thus, for continental areas (SAT-GST<0), positive (negative) trends indicate that differences tend to get smaller (larger) in absolute values, and conversely for the limited coastal areas where SAT-GST differences are positive. Figure   25 1b allows for visualization of this behavior. Note the positive trend after the change when the difference between SAT and GST anomalies becomes smaller with time. Regionally, several circumstances account for impacting SAT-GST long-term coupling.
On the one hand, decreasing SAT-GST differences over land may emerge from two conditions. Firstly, when there is a higher warming rate of SAT relative to GST as depicted in Fig. 8a,b over the northernmost part of North America, Fennoscandia, northeast Russia and some areas of central Eurasia. Secondly, when there is a cooling of both SAT and GST but the later 30 decreases at a higher pace as described in Fig. 8a,b for the region of northeastern Brazil and some areas of Africa. Note that these two scenarios are represented in Fig. 9g with positive trends over these regions after the change. On the other hand, the increase in SAT-GST difference arises either from the effect of rising GST in the presence of stable/decreasing SAT or due to the higher warming rate of GST relative to SAT. The former case is appreciated in Fig. 8a,b for the areas of central and eastern Europe as well as eastern U.S whereas the latter is found over the Indian subcontinent and southeastern Asia. Note that both cases are represented in Fig. 9g with negative trends.
Trends of differences between SAT and GST from the GHG-only and the LULC-only ensembles help understanding the relative contributions to the long-term variations seen in the ALL-F simulations. For instance, the GHG-only ensemble (Fig. 9h) shows similar positive trends as Fig. 9g over northern North America, Fennoscandia, northeast Russia and central Eurasia, though with a much larger magnitude and geographical extension. Correspondingly, negative trends after the change in the 5 LULC-only ensemble (Fig. 9i) are comparable to those in Fig. 9g for the regions of central and eastern Europe, eastern U.S., the Indian subcontinent and southeastern Asia. Additionally, the positive values over Brazil as well as over central and southern Africa in Fig. 9g are also depicted in Fig. 9i.
Interestingly, the two-phase regression analysis does not expose any variation in SAT-GST long-term relationship over southeastern China, where the linear trends during the industrial period show a relative strong decrease in both SAT and GST in the 10 ALL-F and the OZ/AER-only ensembles (Fig. 8a,b,g,h). Furthermore, the linear trend of SAT-GST differences during industrial times for the OZ/AER-only simulations (Fig. 10) does not exhibit any SAT-GST decoupling over this region either. This suggests that the dominant effect of OZ/AER forcing on the SAT and GST responses over this region is not affecting their long-term coupling. Nonetheless, Fig. 10 illustrates some interesting aspects of the SAT-GST relationship in the OZ/AER-only ensemble such as the negative contribution to the SAT-GST trends over North America, northern Europe, the Tibetan Plateau 15 and central Asia. Additionally, the positive trends over northern Siberia are also notable as well as the positive values over some relative small areas of central and eastern Africa, the coast of Angola and eastern Brazil. Although the bulk of these SAT-GST responses depicted in Fig. 10 does not translate to SAT-GST long-term decoupling in the ALL-F ensemble, they play an important role either in counteracting the influence of other external forcings or contributing to decoupling-related processes over some regions. 20 The following paragraphs aim at providing an insight into the relative contribution of the individual forcings and the associated physical mechanisms to the variations of the long-term SAT-GST association detected in Figs. 9 and 10.
In the cases of the long-term variations due to the LULC influence, changes in vegetation cover can be considered the main driving factor since they alter the radiative fluxes and water cycling at the surface due to the modification of the physical properties such as albedo, roughness and evapotranspiration (Pongratz et al., 2010). Figure 11 gives an example of how long-term 25 changes in the energy fluxes at the surface due to LULC changes do impact the SAT-GST coupling at long time scales. It shows the 31-yr low pass filter outputs of SAT, GST, reflected shortwave radiation (RSW) and SHFLX evolution for a characteristic grid-point over the Great Lakes region (U.S.) where a warming of GST relative to SAT is simulated in CESM-LME during the industrial period. Results are shown for one of the members of the ALL-F, LULC-only, GHG-only and OZ/AER-only ensembles. Around 1800 CE SAT tends to decrease whereas GST tends to increase in both the ALL-F and LULC-only simulations, 30 while GHG-only and OZ/AER-only simulations do not display the same behavior that produces larger differences between SAT and GST represented by negative trends in Fig. 9d. At the same time RSW and SHFLX exhibit large long-term variations in the ALL-F and LULC-only simulations. Therefore, this modification of the long-term SAT-GST relationship is clearly the response to LULC changes. Such variations in the surface energy fluxes over this region are likely a response of vegetation replacement from forested areas to grassland or croplands. Forested landscapes dissipate SHFLX more efficiently to the atmo- sphere due to a higher surface roughness than open fields (Jackson et al., 2008). In addition, lower vegetation types have higher reflectivity than forests. All the previous contribute to SAT decreases over these regions, especially in winter. Furthermore, deforestation process at mid-and high-latitudes tends to positively feedback with increases in snow cover (Anderson et al., 2011). This type of changes in LULC contribute to increasing albedo which is reinforced by changes in snow cover at these latitudes. Additionally, higher DJF snow cover tends to increase insulation of the soil from the cold overlying air. Therefore, the 5 different temperature response of SAT and GST is related to higher albedo that leads to cooling the atmosphere, the reduction of SHFLX from the ground surface to the atmosphere and the additional snow cover positive feedback. This particular LULC process is important for corrupting the SAT-GST coupling at time scales relevant for the borehole theory (centennial) since the thermal signature recorded in GST during the industrial period would not be representative of the past long-term SAT variations in regions where this effect is dominant. For the areas of central and eastern Europe where a GST warming relative to SAT 10 is also observed in Fig. 9c the mechanisms are similar to those described in Fig. 11 because these areas were also subject to intense transformation from forested areas to cropland at the beginning of the industrial period according to the LULC forcings considered in the CESM-LME (Pongratz et al., 2008;Hurtt et al., 2009).
Changes in vegetation cover are also important for long-term SAT and GST temperature differential response over tropical regions although the driver mechanisms are different than those at mid-and high-latitudes (Lee et al., 2011), and they deserve 15 being considered. At the northeast of Brazil both SAT and GST have negative trends in both the ALL-F and the LULC-only simulations during the industrial period (Fig. 8a,b,e,f). However, the decrease of GST is much larger than the one of SAT as represented in Fig. 9c,i with positive trends after the change. The changes observed at this location in the energy fluxes likely correspond to transitions from open lands to a forested area (reforestation or afforestation) leading to lower albedo and higher evapotranspiration rate as is shown in Fig. 12. This situation leads to an apparent long-term cooling of GST relative to SAT at this location. The temperature response over this area is 25 influenced by different mechanisms. First, the conversion from lower-to higher-type vegetation reduces the solar radiation that impinges on the surface and GST decreases by a radiative effect. Second, forested lands usually have lower albedo thus absorb more shortwave radiation (Zhao and Jackson, 2014). This surplus of energy is balanced by the increase in transpiration, and consequently GST decreases also by a non-radiative process. The latter is especially important in humid climates (von Randow et al., 2004) as the one in this example. According to this, there is a net heat loss at the ground surface with a higher decrease 30 in GST relative to the overlying air. The SAT-GST coupling becomes strong again after some period when the new vegetation cover reaches a stable estate by mid 20th century. For the areas of Africa where there is also a cooling of GST relative to SAT (Fig. 9g,i), the mechanisms are comparable to those describe for Fig. 12.
Over some of these tropical regions there is also a contribution from the OZ/AER forcing to the SAT-GST response. Note that Fig. 10 shows positive trends over Uganda, the coast of Angola and over the east of Brazil, also noticeable in Fig. 9c. At these regions the incoming solar radiation is reduced due to the effect of aerosols as they are an important element for cloud formation, thus contributing to a higher reflectivity of solar radiation (Tao et al., 2012). The description of the specific processes related to aerosols-cloud interaction goes beyond the scope of this study, therefore only the influence on the energy balance at the surface is addressed. As lower shortwave radiation impinges the surface, the energy gain decreases and therefore the ground surface heating is lower. The reduction in the energy gain at the surface is compensated by a lower dissipation 5 via sensible heat whereas the fluxes of latent heat remains relatively constant or even increases in some areas due to a higher moisture as a result of increased precipitation. To illustrate this mechanisms Fig. 13 shows the temporal evolution of SAT, GST, SHFLX, LHFLX and surface incoming shortwave radiation (SSW) for a grid-point located at the region of Uganda. Note in the ALL-F simulation the reduction in SSW after 1900 CE that leads to a decrease in SHFLX. Interestingly, for this location LHFLX experiences an increase at the same time as a result of increase precipitation (not shown) thus providing the source of 10 moisture for evapotranspiration. This situation leads to a higher decrease in GST relative to SAT due to a net loss of energy at the ground surface. Similar SAT and GST responses are also observed in the OZ/AER-only and LULC-only simulations, whereas in the GHG-only, SAT and GST increase.
The SAT-GST decoupling processes described above for individual grid-points are important also at larger spatial scales. Figure   14(left) shows an extension of the mechanism depicted in Fig. 12 including a larger area over the northeast of Brazil (between 15 1 • S-11 • S and 47 • W-35 • W). The negative trend since the 18th century is less accentuated for SAT than for GST (-0,014 and -0,053 • C/decade respectively in the All-F simulation and -0.007 and -0.033 in the LULC-only one) indicating strong contribution of past LULC changes. The same air-subsurface temperature response occurs in other tropical and subtropical areas such as the east of Africa.
The regional analysis is extended to the area of southeastern China in order to illustrate additional information about the in-20 fluence of different external forcings on the SAT-GST relationship. As discussed before, the negative trends for both SAT and GST over this region during the industrial period represented in the ALL-F and OZ/AER-only ensembles (Fig. 8a,b,g,h) do not entail a corruption of the SAT-GST long-term coupling. Figure 14 Although the majority of the important long-term variations in the SAT-GST relationship at regional and local scales observed in the ALL-F ensemble (Fig. 9c)  driver for long-term SAT-GST decoupling. The positive trends after the change over Fennoscandia, northeast Russia and north of North America observed in the ALL-F ensemble (Fig. 9c) can be to a large extent explained by the influence of GHG-only ensemble inasmuch as a broadly similar picture is also portrayed over these regions in Fig. 9f. In the case of the GHG-only ensemble, the strong warming of SAT relative to GST over these regions is driven by the increasing air temperature during industrial times due to the the positive radiative forcing of GHGs in the presence of a considerable long-term reduction in 5 simulated snow cover. González-Rouco et al. (2009) showed that such an scenario would lead to a higher exposure of soil to cold winter air, therefore the soil would partially record colder temperatures, previously prevented by the snow cover insulating effect. In the ALL-F ensemble this effect is damped for the fact of considering additional forcings that keep the snow cover relatively constant during industrial times. For instance, the contribution of the OZ/AER forcing is particularly important for counteracting the effect of GHGs since it leads to colder climate conditions due to its negative radiative forcing. Note the strong 10 negative trends in Fig. 10 over North America, northern Europe and the Tibetan Plateau that partially balance the effect of the GHGs over these regions. Nonetheless, in the ALL-F there is still an overall SAT warming relative to GST since the relatively stable snow cover is insulating the soil from a warmer SAT reducing the overall offset between them (Bartlett et al., 2005). leading to a small increase in the global SAT-GST offset. In the same way the OZ/AER-only ensemble shows an increase in global snow cover while the SAT-GST offset in industrial times exhibits a relatively slight increase. The interaction between different external forcings in the ALL-F ensemble leads to an stable snow cover during the industrial period since the sharp 25 decrease in snow induced by the GHG forcing is partially compensated by the counteracting effect of the LULC and OZ/AER forcings. Additional forcings such as volcanic may also contribute to counteracting GHGs effect at multi-decadal timescales (not shown). Consequently, in the presence of a warmer climate there is a difference in the warming rate of SAT and GST in industrial times at global scale in the ALL-F ensemble member (0,025 and 0,018 • C/decade respectively). This scenario leads to a net effect of long-term decrease in the SAT-GST differences starting around 1800 CE as discussed in Fig. 1 and also 30 evident in Fig. 15.
This work evaluates the stationarity of the coupling between SAT and GST temperatures as simulated by the CESM in an ensemble of experiments spanning the LM. The initial motivation for this work roots on previous literature , 2009García-García et al., 2016) that addresses the realism of the borehole hypothesis for climate reconstruction, namely, that SAT and GST vary synchronously and that reconstructing past GST changes from borehole temperature 5 profiles is a good proxy for past SAT variations. The use of the CESM-LME allows for analyzing the influence of forcing changes on the SAT-GST covariability, both individually and as a group, by considering the different all-forcing and singleforcing ensembles. Additionally, having several experiment ensemble members for each given forcing type, allows for disentangling the effects of internal variability from that of the forcing response. Ultimately, the coupling between SAT and GST is assessed at global and also at regional/local scales. This last purpose needs considering different mechanisms that contribute 10 to SAT-GST variability within different climate types. In doing so, a variety of factors and conditions that contribute to the surface energy balance with different mechanisms is provided.
The CESM-LME shows that at global scale the SAT-GST coupling is strong above multi-decadal timescales since GST tracks SAT throughout the LM, as found in previous works. However in spite of the strong coupling, the CESM-LME also reflects that the SAT-GST relationship has not remained constant along the whole LM at these spatio-temporal scales. Hence, the nature of 15 such variation is evaluated.
Globally, snow cover is the most important agent in modulating the connection between SAT and GST. Therefore the variation of the SAT-GST relationship described by the CESM-LME simulations should be in principle driven by variations of global snow cover. Nevertheless, the simulated snow cover remains relatively stable at the time when SAT-GST coupling varies, thus this change can not be solely explained by the influence of the snow cover. With this in mind, we explored with some de-20 tail different processes that may influence the SAT-GST relationship at different spatio-temporal scales. Firstly, we address processes acting at seasonal time scales that were identified from an spatial analysis of the SAT-GST differences and correlations. Secondly, the long-term evolution of the SAT-GST relationship is evaluated in the ALL-F, LULC-only, GHG-only and OZ/AER-only ensembles.
Several processes over different regions either relevant during winter/summer play an important role on impacting the SAT 25 and GST coupling. Some examples are: snow cover over mid and high latitudes, discontinuous snow cover over the Tibetan Plateau region and seasonal variations in the energy fluxes at the surface. Although these processes are important for disrupting the SAT-GST relationship at seasonal scales, they have no implications on the long-term coupling if they are stationary.
Nonetheless, if they experience variations with time the SAT-GST long-term relation may be impacted.
As discussed in Sect. 4.1 some of the anthropogenic external forcings have the potential to impose long-term variations on 30 processes that regulate the relationship between SAT and GST. Among them, the LULC changes are the most important ones since they modify the energy fluxes at the ground-air interface, and consequently corrupt the SAT-GST coupling locally and regionally at various time scales. An example is the response to the deforestation processes triggered by the expansion of agriculture mainly during the industrial period at mid and high latitudes, where SAT and GST long-term coupling is impacted due to the variations in the albedo, surface roughness and hydrology. Similar decoupling processes related to LULC changes are found over different regions around the globe, as the one described in Sect. 4.1 over northeastern Brazil, some areas over Africa and over the Indian subcontinent. All of them driven by the long-term modifications of the energy fluxes at the surface either from increased evapotranspiration, reduced energy dissipation via sensible heat and others. Besides these kind of decoupling processes induced by individual forcings, the interactions of a variety of mechanisms and feedbacks from different 5 external forcings can also exert an influence on the long-term SAT-GST relationship at different spatial scales. For instance, the effect of GHGs leads to a reduction of the snow cover during industrial times that is counterbalanced by the opposite effect of both LULC and OZ/AER forcings. As a consequence, the snow cover remains relatively stable over some regions during the industrial period in the presence of a warmer climate. This situation leads to a difference in the SAT and GST long-term evolution during the industrial period, since the snow cover is insulating the soil from a warmer SAT. This effect is present 10 over the NH high-latitudes of North America, Fennoscandia and northeastern Russia. Indeed, at global scale the combination of a steady snow cover under warmer climate conditions is the dominant effect for explaining the variations in the long-term SAT-GST relationship.
Our findings indicate that the assumption of a strong relation between SAT and GST may be impacted from local to regional scales by different mechanisms especially by the influence of LULC changes due to the modification of the energy balance at 15 the surface. Therefore, the interpretation of temperature reconstructions from borehole measurements at this spatial scales must consider LULC changes as a source of possible bias. The effects of additional external forcings may also exert some influence on processes such as variations in the snow cover, hydrology and other land surface properties, that may in turn feedback on the SAT-GST long-term coupling. At global scale the influence of such local and regional decoupling processes is only ca. 0.05 • C, hence the SAT-GST coupling at this spatial scale is supported by the CESM-LME. Res. Lett., 30, 1721, doi:10.1029/2003GL017 337, 2003.        (Table 3). Bottom. Same as middle for SHFLX, LHFLX and water content in the top 10 cm of the soil (W-10cm).