Impact of deforestation on subsurface temperature profiles: implications for the borehole paleoclimate record

Subsurface temperature profiles measured in boreholes are one of the important archives of paleoclimate data for reconstructing the climate of the past 2000 years. Subsurface temperatures are a function of past ground surface temperatures (GST), however GSTs are influenced both by changes in land-use and changes in regional climate. Thus the history of deforestation at borehole sampling locations represents a potential uncertainty in the reconstructed temperature history at the site. Here a fully coupled Earth system model is used estimate the magnitude of the subsurface temperature anomaly from deforestation events from a global perspective. The model simulations suggest that warming of the ground surface is the dominant response to deforestation, consistent with the limited field data that exist. The magnitude of the temperature anomaly varies by environment with a global average anomaly of 0.85 °C with a range of −0.48 °C to 1.78 °C. The warming originates from a reduction in the efficiency of turbulent energy flux to the atmosphere overcompensating an increase in albedo. Overall our simulations suggest that deforestation has a large impact on subsurface temperatures for centuries following deforestation and thus GST reconstructions should take into account previous deforestation events.


Introduction
Energy exchange between the atmosphere and the land surface determines the ground surface temperature [1,2]. Across most of the land surface this energy exchange is mediated by living vegetation which influences the roughness of the surface, the albedo of the surface, and the flux of water and latent heat to the atmosphere [2]. Removal of vegetation through deforestation alters the surface energy balance through an assortment of mechanisms. The removal of tall vegetation (trees) reduces the roughness length of the surface which decreases the turbulence driven sensible heat exchange between the land surface and the atmosphere. Except over melting snow or ice, the sensible heat flux on average transfers heat from the surface to the atmosphere, thus the removal of tall trees should make this mechanism less efficient and create a surface warming effect [1,2]. Both coniferous and deciduous trees tend to have higher albedo than grasses [3]. By shedding snow and exposing their needles during winter coniferous trees decrease the albedo of the surface in regions where they occur in conjunction with snow [4]. Thus removal of trees tends to increase the albedo of the surface creating a cooling effect. In addition, removal of trees decreases shading of the surface which will have a warming effect on the surface [4]. Therefore deforestation alters the energy balance of the surface in counteracting ways which contingent on the specifics of the environment could cause either warming or cooling of the surface.
Warming or cooling of the surface can also be caused by local climate change. No matter the origin of the temperature change the temperature anomaly is propagated into the ground creating a temperature anomaly at depth [5][6][7][8][9]. Thus the temperature profile of the subsurface is a function of past ground surface temperatures. This observation has led to the field of borehole paleoclimatology, which uses thermal physics to infer past ground surface temperatures from geothermal data [6,10]. Such geothermal data have been widely used at local, regional and global scales to reconstruct past ground surface temperature changes [8,9,[11][12][13][14], to estimate the magnitude of the continental heat storage [15,16] and to assess the fidelity of land surface components of climate models [17].
Fundamental to retrieving climatic information from subsurface temperature data is the assumption that the Earth's upper crust is in thermal equilibrium, and that the temperature distribution in its upper few kilometres is determined by the long-term (< 1000 years) surface air temperature and the heat flow from the interior of the Earth (constant at the > 1 Ma time scale). Under these conditions, surface energy balance changes are recorded in the subsurface as temperature anomalies to the background quasi-steady state geothermal field. The contribution from the interior of the Earth can be estimated and removed [7], and the subsurface thermal anomalies can be interpreted as indicators of climate-induced changes at the ground surface. Changes in surface energy balance propagate typically such that, the upper ∼ 500 m of the subsurface yields information on the climate of the past millennium. There is evidence of long-term coupling between the lower atmosphere and the ground from direct calculations [18,19], and also from analyses of CMIP5/PMIP3 last millennium general circulation models' simulations [20]. Thus geothermal data offer, not only a complementary view of long-term temperature changes at local, regional, hemispheric and global scales, but also a record of the changes in subsurface energy storage over the last millennia [21].
Temperatures measured in boreholes must pass a quality criteria selection to ensure that these data are not contaminated by non-climatic factors [22][23][24][25][26][27][28][29]. One of the most serious problems affecting geothermal data is that many of the data-sets have been collected at sites affected by past deforestation [30][31][32]. Therefore it is of interest to examine from a theoretical perspective what is the expected magnitude of the impact of deforestation on subsurface temperatures, the time after which the signal can be detected in subsurface data, and how the effect of deforestation varies across many different climates and environments. Here, we performed a series of model experiments to explore the effect of deforestation on subsurface temperatures using a fully coupled spatially resolved intermediate complexity Earth system model. We find that for most environments deforestation causes warming of the surface and that a global mean maximum signal of 0.15°C with depth remains in the subsurface temperatures profile 500 years after deforestation.

Model description
The University of Victoria Earth system climate model (UVic ESCM) is a climate model of intermediate complexity composed of the full three-dimensional ocean general circulation model coupled to a simplified moisture and energy balance atmosphere, thermodynamic-dynamic sea-ice module [33], and land surface scheme [34,35]. The UVic ESCM is a deterministic model that conserves mass and energy to machine precision [36]. The land surface scheme is a highly modified variant on the Hadley Centre's Moses scheme [37]. The land-atmosphere interface consists of five Plant Function Types (PFTs) and bare ground. The five PFTs are broadleaf trees, needle-leaf trees, C3 grass, C4 grass and shrubs. The PFTs are governed by the Triffid dynamic vegetation model [37]. Each PFT competes for space on each grid-cell with dominance determined by a competition hierarchy, suitability to local climate, and the resources each PFTs has available to invest in expansion [38]. Human land use is prescribed by designating a fraction of a grid-cell that can only be covered by C3 or C4 grasses, or bare ground. When human land-use increases the shrub, broadleaf, and needle-leaf PFTs are removed and replaced with bare ground. Grass PFTs are allowed to over-take the agricultural area and typically within a year will replace bare ground. All PFTs compete over the non-agricultural fraction of the grid-cell. The surface energy balance is solved for each PFTand bareground and takes the form: where SW n is the net shortwave radiation, LW n is net longwave radiation, S H is sensible heat, L H is latent heat flux and G is ground heat flux. The net shortwave radiation is calculated from incoming shortwave radiation (SW ↓ ) and the albedo of the surface (a): The albedo of the surface varies depending on the presence or absence of snow, the PFT and the PFTs state. The snow free albedo is: where a s is the albedo of bare soil, a c is the maximum canopy albedo (unique for each PFT), k is a constant, and L is the leaf-area-index for the PFT, which varies in time [38]. Accounting for the effect of leaf-area-index on albedo is how the model handles shading by vegetation and contributes to the gradient between surface temperature and air temperature. The net longwave radiation is calculated from the incoming longwave (LW ↓ ), and the skin temperature (T Ã ) of the surface: Environ. Res. Lett. 12 (2017) 074014 where s is the Stephan-Boltzmann's constant. The skin temperature is a free variable which is used to close the energy balance [35]. The sensible and latent heat fluxes are solved using the bulk-aerodynamic approach [38]. Sensible heat flux is computes as: and latent heat flux as: where r a is the density of air, C p is the volumetric heat capacity of air T a is 2 m air temperature, r a is aerodynamic resistance, q sat is the saturated specific humidity (taken as a function of skin temperature), q a is the specific humidity, and L v is the latent heat of vapourization or sublimation contingent on the presence or absence of snow. For bare-ground r a is a function of the soil moisture availability and each PFT. r s is stomatal resistance, a complex function of the rate of photosynthesis and soil water availability [39]. As windspeed and vegetation height increase r a decreases reflecting the enhanced sensible and latent heat flux with more turbulent atmospheric conditions. G is calculated by the subsurface scheme described below. The version of the UVic ESCM used here features a deep surface scheme extending to 250 m depth. This deep land-surface was first implemented by [40] to explore the effects of climate change on permafrost thaw. The scheme has subsequently been updated by refs. [41,42] to better account for the permafrost carbon cycle feedback to climate change. The deep subsurface consists of 14 soil layers exponentially increasing in thickness with depth, beginning with a 10 cm thick surface layer. The top eight layers (10 m) are active in the water-cycle, and top six layers (3.35 m) in the carbon-cycle [35]. Water enters the subsurface from the rainfall and is transported through the subsurface based on the porosity of the soil. The deepest hydraulically active layer is gravity drained, with extracted water entering the appropriate river basin. Plant roots extract water from the subsurface to provide for the physiological needs of each PFT [35]. The model includes freeze-thaw dynamics by-which water can exist as solid or liquid based on equations the minimize Gibbs free energy. The equations account for the effect of valence forces on freezing-point such that liquid water and ice can exist simultaneously over a range of temperatures in the simulated subsurface [35]. Heat flow in the subsurface is determined by the temperature gradient of the subsurface, the flow of water, and latent heat exchange when water freezes or thaws. The bottom boundary condition of the subsurface module is a spatially resolved but temporally invariant geothermal heat flux [35].

Experiment Design
The idea of the present study was to explore the effect of deforestation on deep subsurface temperature with a global perspective. The challenge with using an Earth system model for this purpose is that by altering the local climate deforestation also alters the global climate. Therefore if we simply deforest the entire planet we will be unable to distinguish the effect of local deforestation and global deforestation on subsurface temperatures. The strategy we have employed to avoid the interference is to only deforest 1% of the land-surface in each branch simulation. For the initial state of the experiment human land use was set to zero globally and the UVic ESCM was spun-up under otherwise year 1850 Common Era (C.E.) forcing for 2000 years. From this stable initial state 101 branch simulations were conducted. In all but one branch simulation 1% of the model land grid-cells were set to 95% agricultural land use for one year, then allowed to recover. In the final branch simulation a control run is conducted wherein no land-use-change is imposed. The grid-cells where randomly selected without replacement such that each land surface gridcell was deforested in one simulation. The simulations were conducted for 500 years under constant year 1850 C.E. forcing (except for human land use), with model output saved yearly for the first 100 years and at 25 year averages for the final 400 years of the simulation. Comparison between each deforestation simulation and the control run shows that deforesting 1% of the global land surface has a negligible effect on global temperature with maximum deviations from the average temperature of the control simulation of −5 to 3 mK, consistent with [43].
Land grid cells in the UVic ESCM contain a mixture of plant function types and bare ground, therefore what constitutes a 'forest' within the modelling framework needs to be defined. We chose to consider forests as environments where broad-leaf trees, needle-leaf trees and shrubs cover 67% of the grid cell before deforestation.  regions. By 488 years after deforestation only the northern-most regions of the boreal forest and part of the high Himalayas have not recovered fully. The Broad-leaf trees recover slowest. The central forests of Africa and South America recover first with more marginal forests in the Sahel, South East Asia and East Asia recovering slowest. In some grid cells the broadleaf tree plant function type never recovers. This outcome appears to be a manifestation of the tropical forest-grassland instability characteristic of the Triffid dynamic vegetation model [37]. That is, the model has a bi-stability in some regions between grassland and tropical forests, when the grid is 95% deforested the grass plant function types are stably dominant and the forest never returns. The effect occurs at 68 grid cells at the margins of the tropical forests.

Results
The evolution of the temperature anomaly created by deforestation is shown in figure 2. The figure shows the maximum temperature anomaly created by deforestation and the depth of the anomaly at time slices 5, 55, 95, 188 and 488 years after deforestation. Warming dominates most of the deforested region with only a small region of cooling at the fringe of the temperate forests (see supplementary figure S1). The area weighted mean warming at the surface across deforested grid cells in the first year after deforestation is 0.85°C with a range of −0.48°C to 1.78°C at individual grid cells. As the forest recovers at the surface the maximum temperature anomaly becomes smaller in magnitude and positioned deeper in the subsurface (figure 2). Generally the temperature anomaly is smaller and diminishes faster in the tropical forest relative to the temperate forest, except for the tropical grid cells where the forests does not recover (figure 1). Figure 3 shows the mean maximum temperature anomaly with depth over time averaged over the grid points where the forest recovers. The mean maximum temperature anomaly plateaus at ∼ 0.7°C for the first 50 years of the simulation then  [44]. The change in the surface energy balance is shown for net shortwave radiation, net longwave radiation, sensible and latent heat fluxes in figure 4 for 5, 25, and 65 years after deforestation. The change in shortwave radiation causes a cooling effect consistent with an increase in albedo from switching to grass plant function types and from the effect of the removal of trees on the albedo of snow-covered ground. The change in longwave radiation cools the surface (there is an increase in outgoing longwave radiation). The change in sensible heat flux has a warming effect on the surface as one would expect from a reduction in roughness length of the surface caused by removal of tall trees. The change in latent heat flux initially has a warming effect on the surface, consistent with a decline in the roughness length, but by 65 years after deforestation the change in latent heat relative to the control simulation is having a cooling effect on the surface in the boreal forest but continues to have a warming effect in the tropics. The change from a warming to a cooling effect coincides with the regrowth of the boreal forest and is consistent with higher transpiration from a recovering forest.The rate of evapotranspiration is linked to the rate of photosynthesis and gross primary productivity in the model. Thus the cooling effect of the change in latent heat flux in the boreal forest by year 65 is likely linked to higher growth rates from forest recovery creating higher evapotranspiration and overcompensating for the reduction in turbulent fluxes from the change in roughness length. The net effect of the change in surface energy balance caused by deforestation in the model is warming of the surface driven by changes in sensible heat flux, an effect that applies almost world-wide.

Discussion
Reference [30] comprises one of the only studies to examine borehole temperature and micro-meteorological measurements from nearby forested and deforested study sites. The study was conducted with data from Vancouver Island (temperate rainforest) and Yukon (boreal forest) of Canada. The study showed a ground surface temperature warming of 1.8°C for deforested sites relative to forested sites on Vancouver island, and a warming of 1.2°C at deforested sites in Yukon relative to forested sites. Reference [30] deduced that the change in temperature was from a reduction in evapotranspiration being compensated by and increase in outgoing longwave radiation. The results of the present study are generally consistent with the field study of [30] suggesting a warming effect of deforestation resulting change in radiative fluxes being counteracted by a change in turbulent fluxes. The picture painted by our modelling methods suggest a more complicated change in energy balance terms with large changes in albedo and sensible heat flux, variables which [30] did not measure.
Reference [45] examined the consistency of surface air temperature and subsurface temperature trends in the province of Alberta, Canada. The study found that subsurface temperature warmed faster than air temperatures during the latter half of the 20th century in the boreal forest environment in the northern fraction of the province, but not in the prairie ecosystem of southern Alberta. The study examined several possible sources for the difference in temperature trends and concluded that widespread deforestation from land-development and forest fires was the source of the warming [45]. The findings of [45] are generally consistent with the present study which shows a much larger effect of warming from deforestation in boreal forest environments relative to grass and shrub dominated prairie environments. The global and regional effects of large-scale deforestation are complicated and must account for the effects of the change in surface albedo, conversion of biomass into CO 2 , and feedbacks between evapotranspiration, cloud formation and cloud albedo [4,46]. The present study was specifically designed to avoid the global effects of deforestation by deforesting only 1% of the land surface in 101 parallel iterations of the experiment. Thus the local warming effect on ground surface temperatures from deforestation and global cooling from the biogeophysical effect of deforestation [4] seen in other modelling experiments utilizing Earth system models should not be seen as inherently contradictory.
Our results support the caution conventionally applied in borehole paleo-climatological studies where-by borehole temperature records that are suspected to have been affected by deforestation or forest-fires are eliminated from databases before inversion to retrieve surface temperature change or heat flux [29]. The historical and archeological record of human land-use-change is poor in much of the world, even for the past 500 years [47]. Our study suggests that even 500 years following deforestation the signal of deforestation may remain in the subsurface temperature record. Therefore there may be an opportunity to augment historical and archaeological records by using borehole temperature records that have been disturbed by deforestation to help reconstruct land-use-change history using the data rejected for inferring historical climate variability.    global perspective. The simulations suggest that warming of the surface is the dominant response to deforestation with an average temperature change of 0.85°C for model grid cells which were once forested. The temperature anomaly propagates through the subsurface diminishing and becoming deeper with time. By ∼ 500 years after deforestation a signal with an average strength of 0.15°C remains present. Deforesting the surface alters the energy balance at the surface. In our simulations deforestation reduced the absorption of net radiation and reduces the transport of heat away from the surface from the turbulent energy fluxes, suggesting the change in the sensible heat flux as the source of surface warming. Our results suggest that large scale ground surface temperature reconstructions from subsurface temperatures must take into account possible past deforestation events, especially for poorly documented older data, when estimating ground temperature changes. A conclusion consistent with earlier work e.g. [30,45]. This result is of particular importance when estimating the magnitude of subsurface energy storage because of its role in driving near-surface processes with strong potential positive feedbacks on the future Earth system evolution [41]. Overall our simulations suggest that deforestation has a large impact on subsurface temperatures for centuries following deforestation.