Afforestation impact on soil temperature in regional climate model simulations over Europe

In the context of the first phase of the Euro-CORDEX Flagship Plot Study Land Use and Climate Across Scales (LUCAS), we investigate the biophysical impact of afforestation on the seasonal cycle of soil temperature over the European 20 continent with an ensemble of ten regional climate models. For this purpose, each ensemble member performed two idealized land cover experiments in which Europe is covered either by forests or grasslands. The multi-model mean exhibits a reduction of the annual amplitude of soil temperature (AAST) due to afforestation over all European regions, although this is not a robust feature among the models. In Mediterranean, the spread of simulated AAST response to afforestation is between -4 C to +2 C at 1 meter below the ground while in Scandinavia the inter-model spread ranges from -7 C to +1 C. We show that the 25 large range in the simulated AAST response is due to the representation of the summertime climate processes and is largely explained by inter-model differences in leaf area index (LAI), surface albedo, cloud fraction and soil moisture, when all combined into a multiple linear regression. The changes in these drivers essentially determine the ratio between the increased radiative energy at surface (due to lower albedo in forests) and the increased sum of turbulent heat fluxes (due to mixingfacilitating characteristics of forests), and consequently decide the changes in soil heating with afforestation in each model 30 Finally, we pair FLUXNET sites to compare the simulated results with observation-based evidence of the impact of forest on soil temperature. In line with models, observations indicate a summer ground cooling in forested areas compared to open lands. The vast majority of models agree with the sign of the observed reduction in AAST, although with a large variation in the magnitude of changes. Overall, we aspire to emphasize the biophysical effects of afforestation on soil temperature profile with this study, given that changes in the seasonal cycle of soil temperature potentially perturb crucial biochemical processes. 35

continent with an ensemble of ten regional climate models. For this purpose, each ensemble member performed two idealized land cover experiments in which Europe is covered either by forests or grasslands. The multi-model mean exhibits a reduction of the annual amplitude of soil temperature (AAST) due to afforestation over all European regions, although this is not a robust feature among the models. In Mediterranean, the spread of simulated AAST response to afforestation is between -4 o C to +2 o C at 1 meter below the ground while in Scandinavia the inter-model spread ranges from -7 o C to +1 o C. We show that the 25 large range in the simulated AAST response is due to the representation of the summertime climate processes and is largely explained by inter-model differences in leaf area index (LAI), surface albedo, cloud fraction and soil moisture, when all combined into a multiple linear regression. The changes in these drivers essentially determine the ratio between the increased radiative energy at surface (due to lower albedo in forests) and the increased sum of turbulent heat fluxes (due to mixingfacilitating characteristics of forests), and consequently decide the changes in soil heating with afforestation in each model 30 Finally, we pair FLUXNET sites to compare the simulated results with observation-based evidence of the impact of forest on soil temperature. In line with models, observations indicate a summer ground cooling in forested areas compared to open lands.
The vast majority of models agree with the sign of the observed reduction in AAST, although with a large variation in the magnitude of changes. Overall, we aspire to emphasize the biophysical effects of afforestation on soil temperature profile with this study, given that changes in the seasonal cycle of soil temperature potentially perturb crucial biochemical processes.
Robust knowledge on biophysical impacts of afforestation on soil conditions and its feedbacks on local and regional climate is needed in support of effective land-based climate mitigation and adaption policies.

Introduction
There is currently a strong policy focus on afforestation as a possible greenhouse gases mitigation strategy to meet ambitious climate targets (Grassi et al., 2017). The biogeochemical effects of afforestation or reforestation are mostly related to increased 40 carbon stocks stored in vegetation and soil, as the total carbon stored in forests is nearly three times larger than carbon stored in croplands (Devaraju et al., 2015). However, understanding the full climate consequences of the large-scale deployment of such a strategy requires to consider also the biophysical effects of afforestation arising from changes in evapotranspiration efficiency, rooting depths and soil water holding capacity, surface roughness and surface albedo (Betts, 2000;Bonan, 2008;Davin and de Noblet-Ducoudre, 2010;Perugini et al., 2017;Duveiller et al., 2018).
Previous studies have attempted to quantify the biophysical impact of land-use changes (LUC) on global scale, employing either an ensemble of earth system models (ESMs) (Pitman et al., 2009;Noblet-Ducoudré et al., 2012;Boisier et al., 2012;Lejeune et al., 2018) or applying a single ESM individually (Claussen et al., 2001;Davin et al., 2007;Li et al., 2016). Davin and de Noblet-Ducoudre, 2010 analysed an ESM's sensitivity to idealized global deforestation, indicating that the net biophysical impact results from the balance between radiative and non-radiative processes. In the same study, deforestation 50 induced a warming over the tropical zone owing to a reduction in evapotranspiration rate and surface roughness, whereas a deforestation-induced cooling simulated over the temperate and boreal zones, because an albedo increase provided the dominant influence in these regions. In the context of Land-Use and Climate, IDentification of Robust Impacts model intercomparison project, Noblet-Ducoudré et al., 2012 diagnosed the LUC effects over North America and Eurasia between the present and the pre-industrial era. They found that deforestation caused a systematic surface albedo increase across all 55 seasons, leading to a reduction in available energy accompanied by a decrease in the sum of turbulent fluxes. Furthermore, Lejeune et al., 2018 using a suite of simulations from Coupled Model Intercomparison Project Phase 5 concluded that moderate deforestation over Eurasia and North America has substantially led to a local warming of present-day hot extremes since preindustrial time.
Regional Climate Models (RCMs) constitute dynamical downscaling techniques applied over limited-area domains with 60 boundary conditions either from global reanalysis or global climate model (GCM) output (Katragkou et al., 2015;Giorgi, 2019;Rummukainen, 2016). RCMs operate on higher resolutions than GCMs adding value in regions with complex orography and capturing exreme events (Soares et al., 2012;Cardoso et al., 2013;Warrach-Sagi et al., 2013). RCMs have been also used individually to address the LUC effects on regional scale (Gálos et al., 2013;Tölle et al., 2018;Cherubini et al., 2018;Belušić et al., 2019). Lejeune et al., 2015 used a state-of-the-art RCM to explore the biophysical impacts of possible future deforestation on Amazonian climate. They demonstrated that the projected land cover changes for 2100 could slightly increase the mean annual surface temperature by 0.5 o C and decrease the mean annual rainfall by -0.17 mm day -1 compared to present conditions. Flagship Pilot Study (FPS). It was initiated jointly by the European branch of the Coordinated Downscaling Experiments EURO-CORDEX (Jacob et al., , 2020 and the global model intercomparison study "Land-Use and Climate, IDentification of robust impacts" LUCID (Noblet-Ducoudré et al., 2012). In the first phase of LUCAS, for the first time multimodel and multi-physics simulations were performed under a common experimental protocol to address the RCMs sensitivity to idealized land use changes in Europe. The first experiment assumed a maximum forest coverage while the second a 80 maximum grass coverage over Europe.
Contrasting these two idealized LUC experiments, Davin et al., 2020 analyzed the robustness of RCMs responses to afforestation and according to their results, afforestation implied an albedo-induced warming over northern Europe during winter and spring. Furthermore, the summer near-surface temperature response to afforestation was subject to large uncertainty, strongly related with disagreement among models in land-atmosphere interactions. Analyzing a part of RCM ensemble established within LUCAS FPS, Breil et al., 2020 examined the impact of afforestation on the diurnal temperature cycle in summer. Their results revealed that afforestation dampened the diurnal surface temperature cycle, while the opposite was true for the temperature in the lowest atmospheric model level. Afforestation could also enhance snow melt and modify the land-atmosphere interactions in sub-polar and alpine climates through changes in snow-albedo effect in winter and spring (Mooney et al., 2021).

90
The responses of atmospheric processes to afforestation have been extensively discussed in previous studies. However, the changes in soil temperature profile following the afforestation remain unexplored up to now in LUCAS community.
MacDougall and Beltrami, 2017 suggested that deforestation may has led to a long-term warming of the ground, associated with a reduction of heat fluxes towards the atmosphere. Here, we investigate the biophysical impact of afforestation on soil temperature across Europe, as simulated by a suite of ten RCMs established within the frame of the first phase of FPS LUCAS.

95
The comparison between two extreme LUC scenarios, representing the Europe entirely covered by forest and grass respectively, enable us to gain insights into the biophysical impacts of theoretical afforestation on soil temperature variations (Sect. 3.1). In order to explain the inter-model spread in annual amplitude of soil temperature (Sect. 3.4), we examine the changes in surface energy balance components with respect to differences in land-use parameters across RCMs (Sect. 3.2) and the response of soil moisture content to afforestation in summer (Sect. 3.3). In addition, we compare the simulated impact on 100 AAST with observational evidence based on FLUXNET paired sites, classified as forest or open land (Sect. 3.5).

Regional Climate Model ensemble
Two idealized LUC experiments are carried out using an ensemble of ten RCMs. Table 1 provides a brief description of the RCM ensemble characteristics, while more information about the land and atmospheric setups can be found in Davin et al., 105 2020 and in Table S1 in the supplementary material. Compared to Davin et al., 2020 the current model ensemble includes simulations from two additional RCMs (CCLM-CLM5.0 and WRFc-NoahMP) while one of the RCMs (RCA) is not included here because the necessary variables for the analysis were missing. Compared to CCLM-CLM4.5, CCLM-CLM5.0 is coupled with a modified version of CLM 5.0  that includes biomass heat storage Meier et al., 2019). WRFc-NoahMP shares the same land component as WRFb-NoahMP but differs in the atmospheric set-up.

Experimental design
In LUCAS, each participating RCM undertook two different simulations, applying the same experimental design. In the first 115 experiment, called FOREST, models are forced with a vegetation map representing a Europe fully covered by trees, where they can realistically grow. Bare lands and water bodies were conserved as in original model maps. In the second experiment, called GRASS, the models integrate the same vegetation map, with the only difference that trees are entirely replaced by grasslands. These maps are shown in Figure S1 and detailed description about the creation of maps and the way they are implemented into the respective RCMs can be found in Davin et al., 2020. All simulations are performed over the Euro-120 CORDEX domain (Jacob et al., 2020) with a spatial resolution of 0.44 o (~50 km), forced by ERA-Interim reanalysis data (Dee et al., 2011) at their lateral boundaries and at the lower boundary over sea. Our analysis covers the 30-year period 1986-2015 and focuses on the following eight European sub-regions as described in Christensen and Christensen, 2007: Alps , British Isles , Eastern Europe , France , Iberian Peninsula , Mediterranean , Mid-Europe and Scandinavia (Figure 1).
We consider the FOREST minus GRASS differences, implying the impact of theoretical maximum afforestation on soil 125 temperature in Europe. Fourier's second law of heat conduction is widely used by LSMs to update temperature in each soil layer (Eq. 1): where is the time rate of soil temperature (K s -1 ) and is the spatial gradient of soil temperature (K m -1 ) in the vertical direction z (m). The quantity k represents the thermal diffusivity (m 2 s -1 ) defined at the layer node depth z(m) and is equal to organic matter content), on bulk density and soil wetness. In our experiments, soil texture remains unchanged and RCMs do not account for possible occurrence of heat sources or sinks (such as organic matter or carbon decomposition) in the realm where soil heat flow takes place. Thus, the potential changes in soil wetness with afforestation constitute the main driver of 135 differences in soil thermal diffusivity in our experiments. In this way, we use soil moisture response to afforestation as a potentially explanatory variable of soil temperature variations in RCMs.
Similar to Breil et al., 2020, we employ the residual of energy balance at land surface in order to express the surface energy input into the ground. Specifically, we define as energy input into ground the residual energy amount resulting from available radiative energy (net shortwave + incoming longwave radiation) minus the sum of turbulent heat fluxes (latent and sensible 140 heat flux), without accounting for likely deviation of surface energy budget from assumed balance in models (Constantinidou et al., 2020b). Our analysis on the changes of surface energy balance components due to afforestation is carried out for summer season, when models disagree both on the sign and magnitude of soil temperature response. Thus, the land surface is assumed to be snow-free. Also, the current RCMs do not account for heat storage into biomass over land surface, apart from CCLM-CLM5.0. A detailed description on the structure of land-atmosphere exchange in the different LSMs is provided in Breil et al., 145 2020.

FLUXNET observational data
We use measured or high-quality gap-filled data of soil temperature on monthly scale from the FLUXNET2015 Tier 2 dataset to complement the model-based analysis. Detailed documentation on data and processing methods can be found in Pastorello et al., 2020.
In order to extract the potential effect of afforestation from observations, we employ a space-for-time analogy by searching for pairs of neighboring flux towers located in forest (deciduous, evergreen or mixed trees) and open land (grasslands or croplands), respectively. This approach has been used in previous studies aiming to investigate biophysical impacts of local LUC and evaluate LSM performance (Broucke et al., 2015;Chen et al., 2018). In search for site pairs, the following criteria were defined: the two sites have to 1) be located in the Euro-CORDEX domain, 2) differ in the type of vegetation, one site 155 being forested and the other one being either cropland or grassland, 3) have a linear distance within the horizontal resolution of the performed simulations (less than 50 km), 4) have a common measurement period of at least two years, and 5) provide measurements at common depth below the ground surface. In total, we found 14 sites that met our criteria and combined in ten pairs. Their locations are depicted in Figure 1 and their characteristics are reported in Table 2. The median linear distance between the paired sites is 11.4 km and their median elevation difference is 125 m.

160
The close proximity between the flux towers of paired sites ensures almost similar atmospheric conditions, so that differences can be primarily attributed to the different vegetation cover. Applying a simple linear correlation test, the differences either in elevation or separation between the flux towers of paired sites are not the dominant factors in determining the changes in AAST (r = -0.2 and r = -0.3, respectively).
For comparison with the RCMs, we consider the observed mean monthly soil temperature differences (forest minus open land) averaged over all paired sites. This is then compared with the mean of the grid cells matching the locations of the observational pairs in the various RCMs (FOREST minus GRASS). Modelled soil temperature was linearly interpolated to the common measurement depth that is available for each pair site and averaged over the time period 2003-2014 which covers the observational time span.
Last but not least, the observational setup does not fully resemble the experimental design applied in RCM ensemble. The 170 spatial scale of afforestation applied in models is significantly different from the small forest patches the flux towers are located in. The theoretical maximum afforestation in RCMs has the potential to induce changes in large-scale atmospheric circulation, which can create teleconnections (Swann et al., 2012) that modify the regional cloud cover (Laguë and Swann, 2016) and thus the regional climate conditions. Such feedbacks are not realistic in observations, where most forest measurement locations are located in relatively small forest patches surrounded by open land and is almost unlikely to alter the climate conditions on 175 regional scale.

Soil temperature response
The afforestation (FOREST minus GRASS) effect on the annual amplitude of soil temperature (AAST) at 1 meter below the ground surface is shown in In winter, the soil temperature differences due to afforestation are small in the majority of simulations and with a tendency for an increase.

Surface energy availability
As reported in previous section, the simulated AAST response exhibits great variability during the summer season, when models disagree both on the sign and magnitude of changes. For this reason, it is essential to examine the changes in the 230 available energy to warm the ground across RCMs in summer. The heterogeneity in the changes of surface energy availability with afforestation is largely consistent with the disagreement in the changes of AAST among RCMs. Thus, it is crucial to explore the origin of large inter-model spread in changes of surface energy balance in summer. Below, we examine the afforestation impact on the different components of surface energy balance for each RCM over Mediterranean (Figure 6) and Scandinavia (Figure 7). Similar figures can be found for the rest European 245 sub-regions in the supplementary material (Figures S12-S17). The analysis of differences in surface energy balance components is performed with respect to changes in land-use characteristics in each RCM, such as leaf area index (LAI), surface roughness and surface albedo. Positive (negative) values indicate an increase (decrease) due to afforestation.
In both regions, all models (except from CCLM-TERRA) consistently show an increase in net shortwave radiation at the surface due to afforestation, which is a result of lower albedo in FOREST compared to GRASS experiment. The changes vary 250 across RCMs from -5 W m -2 to 25 W m -2 over Mediterranean and from -15 W m -2 to 35 W m -2 over Scandinavia. In Scandinavia, the changes in net shortwave radiation are stronger than Mediterranean. This is attributed to the fact that the forests in Scandinavia consist of needleleaf trees, which have lower albedo values compared to broadleaf trees dominating in the rest regions of Europe. Furthermore, the WRF configurations exhibit more pronounced increases in net shortwave radiation with respect to other RCMs, which is linked to stronger reductions in albedo values in these simulations (Figure 6f, Figure   255 7f). Moreover, the albedo effect is further intensified by a reduction in cloud fraction with afforestation over Scandinavia in WRF configurations (Figure 7c). In CCLM-TERRA, the reduced net shortwave radiation is due to a pronounced increase in cloud fraction with afforestation triggered by a strong and widespread increase in evaporation rates . Cloud fraction is also increased with afforestation in other CCLM configurations, however the reduced incoming shortwave radiation is offset by the albedo effect and thus the changes in net shortwave radiation have positive sign in these simulations.
The increase in available radiative energy at the surface with afforestation is followed by an increase in sensible heat flux, which is another robust feature among simulations. According to Breil et al., 2020, the  To sum up, all RCMs respond to afforestation in the same way. That is, afforestation leads to increased available radiative energy at the surface due to lower albedo values in FOREST experiment compared to GRASS. In parallel, a large part of this additional radiative energy is transformed into turbulent heat energy due to the mixing-facilitating forest characteristics, such as the high LAI and roughness values, which enhance the heat exchange between the ground and upper atmosphere. The

280
balance between the increased available radiative energy and the increased sum of turbulent heat fluxes will determine if the surface energy input into the soil will be increased or decreased with afforestation in each RCM. Since these processes are differently weighted in each modelling system depending on land-use characteristics, the resulting energy input into the soil varies within the model ensemble in terms of the sign and magnitude of changes. In CCLM-TERRA, CCLM-VEG3D, CCLM-CLM4.5, CCLM-CLM5.0 and RegCM-CLM4.5 the soil heating is decreased with afforestation in summer over Mediterranean

285
and Scandinavia, because the increased available radiative energy is compensated by the increased sum of turbulent heat fluxes.
On the other hand, REMO-iMOVE and the sub-ensemble built around NoahMP exhibit an increase in soil heating with afforestation, since the increase in the sum of turbulent heat fluxes is not enough to compensate their pronounced increase in net shortwave radiation.

Soil moisture 290
The changes in soil moisture could also have key role in describing the simulated soil temperature response to afforestation, because they affect the thermal diffusivity within the soil column. It is expected that a drier (wetter) soil column would lead to a larger (smaller) AAST owing to its smaller (larger) heat capacity, when considering equal soil heat fluxes between the two experiments.
In Figure 8 we map the mean summer differences in soil moisture content ( The surface water balance (P-E), defined as the difference between precipitation (P) and total evapotranspiration (E), decreases 305 with afforestation during summer in the majority of models over the whole Europe ( Figure S18). In most simulations, the decrease in the terrestrial water budget originates from increased evapotranspiration rates with afforestation. In summer, high LAI values do not allow solar radiation to reach the ground surface, as a result soil evaporation is limited and transpiration dominates overall evapotranspiration (Bonan, 2008). Specific characteristics, such as the big leaf area, the deep roots, the great available energy due to low albedo and the mixing of the upper atmospheric boundary layer because of the high surface This is probably linked with low atmospheric demands for hydrates in FOREST experiment of CCLM-VEG3D . In WRFa-NoahMP, the use of Grell-Freitas as convection scheme, exploits the transpiration facilitating features of forests causing extreme soil drying from very early in summer. Therefore, the evapotranspiration rate lowers with afforestation, 315 because the dry soil is not able to satisfy the atmospheric needs for hydrates.
The soil moisture changes with depth would indirectly reveal the afforestation effect on the evapotranspiration process during summer. The water uptake for transpiration occurs in different depths within the soil column for grasslands and forests. In grasslands, the soil water needed for transpiration is extracted from shallow layers, because the large fraction of their roots is located there, depleting the moisture of upper soil. On the other hand, forests have a deeper root distribution, thus consuming 320 water from a bigger soil water reservoir. In Figure 9 we show the afforestation-induced soil moisture changes within the top 1 meter of the soil over Mediterranean and Scandinavia. Similar plots for the other sub-regions can be found in Figure S19 of the supplementary material. The heterogeneity of SMC changes with depth is evident in most models, especially in Mediterranean. In Scandinavia, distinct drying of the uppermost soil layers is shown by some models, especially CCLM-CLM4.5 and CCLM-CLM5.0, which is related to changes in water amounts from snow melt. The different structures of land models and the various descriptions of physiological characteristics of plants in LSMs, such as the root distributions, differentiate the pattern of SMC changes with depth among the simulations. Also, possible biases in the representation of surface fluxes potentially affect the afforestation effect on soil moisture. For example, in CCLM-TERRA the latent heat fluxes are strongly increased with afforestation, as discussed in previous studies Breil et al., 2020), inducing intense drying of the soil column.

The origin of inter-model spread in AAST
The widespread and homogeneous soil drying with afforestation, mentioned in previous section, is not consistent with the mixed AAST response. On the other hand, it is noted higher agreement between the pattern of changes in soil heating and in AAST. In section 3.2, we showed that the afforestation impact on radiative processes, such as the decrease in surface albedo, increase the available radiative energy at the surface. In parallel, the afforestation effect on non-radiative processes, removes 335 a large part of thermal energy from surface to atmosphere in the form of sensible heat flux. The balance between these processes will determine if the surface energy input into the soil will be increased or decreased with afforestation in each RCM. However, the above biophysical processes are differently weighted across RCMs depending on land-use characteristics, like surface roughness, albedo and LAI, which affect the turbulent mixing and the amount of the absorbed solar energy at the surface.
Furthermore, the response of cloud fraction to afforestation is another important factor, which affects the soil heating, because 340 of its impact on the incoming shortwave radiation at the surface.
With the aim to quantify the effect of changes in above mentioned quantities on the simulated AAST response to afforestation, we conduct a linear regression analysis over all the European sub-regions. More specifically, we use the mean summer changes in albedo, LAI, cloud fraction and soil moisture content as explanatory (independent) variables, to determine to what extent they influence the changes in AAST (dependent variable). When we regress all the explanatory variables against the simulated

345
AAST response, we find that the coefficient of multiple determination (R 2 ) is above 80% in all regions, which indicates the key role of the selected drivers in shaping the effect of afforestation on soil temperature (Figure 10). In southern regions, Mediterranean and Iberian Peninsula, the albedo effect predicts the largest part of the inter-model spread in AAST response.
Over regions of central Europe (Mid-Europe, eastern Europe, France, British Isles) the predictive ability of albedo changes remains strong, although the cloud fraction is the dominating factor which effectively explains the inter-model variance over 350 these regions. Soil moisture also contributes to the explanation of the inter-model spread in AAST over the regions of central Europe, although is not a dominating driver. In Scandinavia, the simulated AAST response is largely explained by differences in LAI across RCMs, with cloud fraction also substantially contributing to the prediction of the inter-model spread. The changes in LAI are potentially connected with the simulated cloud fraction response, since higher LAI values could facilitate the evaporation rates triggering an increase in cloud cover. This interaction effect between two or more physical processes which are used as explanatory variables constitutes a caveat of the used statistical approach, with result to reducing the effectiveness of the corresponding drivers in predicting the response of the dependent variable.

FLUXNET paired sites
In this section, we compare the simulated impact on AAST with observational evidence of afforestation effect on soil temperature, based on ten FLUXNET paired sites. In winter, simulations and observations illustrate insignificant changes in 360 soil temperature with afforestation (Figure 11). The magnitude of afforestation effect in the observations is amplified during summer, revealing a strong cooling up to -3 o C. The majority of models captures the seasonal pattern of changes in soil temperature and particularly the observed summer cooling, albeit with considerable variation in the magnitude of changes.
CCLM-TERRA shows the largest changes in summer soil temperature (-5

Discussion & Conclusions
In this study, we employed the experimental design established within LUCAS FPS, to investigate the afforestation impact on soil temperature over the Euro-CORDEX domain. Two idealized land cover change experiments performed by an ensemble of ten RCMs, in which the European land surface is represented as fully covered by forest and grass, respectively. The majority of simulations showed a dampening of the annual soil temperature cycle with afforestation, owing to changes in summer soil temperature. A large inter-model spread produced, ranging from -7 o C to +2 o C depending on model and region.
The changes in AAST with afforestation found to be consistent with summer changes in available energy to warm the ground 380 across models and regions. In other words, RCMs which showed a ground cooling following afforestation, tend to simulate a reduction in surface energy input into the ground, and vice versa. What differentiates the sign of changes in soil heating across models, is the balance between two biophysical processes, which are greatly affected by afforestation. First, it is the increased available radiative energy at the surface, due to lower albedo in forests, and second it is the increased sum of turbulent heat fluxes (mostly sensible heat flux), owing to mixing-facilitating characteristics in forests, such as high LAI and surface 385 roughness values, which enhance the heat exchange between ground and atmosphere. Although these physical processes are differently weighted in LSMs depending on land-use characteristics, such as surface albedo, surface roughness and LAI, while subsequent atmospheric feedbacks, such as the cloud cover changes, can influence the surface fluxes. Thus, the magnitude of afforestation effect on net shortwave radiation and on turbulent heat fluxes is differently pronounced across models. In six out of ten RCMs of the ensemble, the increased available radiative energy is compensated by the increased sum of turbulent heat fluxes, thus simulating a decrease in soil heating with afforestation and finally a reduction in soil temperature, while the opposite is true for the other four modelling systems. Finally, the changes in albedo, LAI, cloud fraction and soil moisture found to explain more than 80% of inter-model variance in AAST response in all sub-regions.
Previous studies which addressed the effects of LUC on soil temperature have reported similar results with the present work.
Ni et al., 2019 employed field monitoring on a landscape consisted of tree and grass covered ground, to investigate the soil 395 temperature effects on root water uptake for a time period from July to November. They found that soil temperature under the grass-covered ground had larger fluctuations and slightly higher values compared to tree-covered ground in summer.

405
In line with recent findings from observations and model-based studies (Jia et al., 2017;Ren et al., 2018;Zhang et al., 2018;, we found that afforestation induced a widespread soil moisture reduction in summer, implying smaller soil heat capacity. This was also a robust feature among the models, albeit with a considerable range in the magnitude of changes.
Soil moisture decrease with afforestation resulted from large drying of deep layers, related to the fact that forests and grasslands extract soil water for transpiration process from different soil depths. Although, the homogeneous soil drying and thus the 410 smaller soil heat capacity is not consistent with the afforestation-induced decrease of soil temperature in the majority of models, explaining only a small part of inter-model variance in AAST response in regions of central Europe.
Based on paired observations from FLUXNET dataset, we evaluated the simulated soil temperature response to afforestation.
The vast majority of models agreed with the observational evidence that showed a summer ground cooling in forested areas compared to open land. The paired sites exhibited a mean reduction of -3 o C in AAST, while the simulated response varied 415 from -5 o C to 1 o C.
The current ensemble enables us to address the role of atmospheric and land processes in the representation of biophysical forcing of land cover change, since it involves simulations which share the same atmospheric model coupled to different land components, or share the same LSM with different atmospheric set-ups. The switch from CCLM to RegCM when both coupled to CLM4.5 did not induce important changes in model results, implying the dominance of land processes in these simulations.

420
Among the suite of models which share the NoahMP LSM, the atmospheric configuration selected for WRFb-NoahMP and WRFc-NoahMP significantly refined the afforestation effect on soil temperature, compared to WRFa-NoahMP. Future studies should focus on the evaluation of model performances, similar to Katragkou et al., 2015 andConstantinidou et al., 2020a, in order to identify the origins of systematic biases and improve the representation of climate processes in simulations. Moreover, our results stress the crucial role of LSM in the simulation of the biophysical effects of afforestation on soil conditions. Among 425 the LSMs coupled to the CCLM model, the choice of CLM significantly improves the representation of afforestation impact on AAST. Also, WRF coupled to CLM4.0 agreed better with observations than WRF coupled to NoahMP. Another issue is the problematic behaviors in model performances stemming from unrealistic descriptions of the physical plant functioning in LSMs. Meier et al., 2018 improved the representation of the evapotranspiration with land cover change in CLM4.5, modifying parameters related to transpiration process, such as the root distribution and water uptake formulation.

430
Research has accounted for the contribution of historical deforestation to present climate conditions. Last years, governments and non-governmental organizations are planning re/afforestation programs around the world with the purpose to mitigate the negative effects of anthropogenic activities on climate. With our study, we aspire to contribute to the deeper understanding of the scientific community on the biophysical effects of afforestation on soil conditions. Future studies focused on the consequences of afforestation from biological or chemical aspect, are encouraged to consider our results, in order to draw 435 comprehensive conclusions on important climate processes in which afforestation is involved, such as the carbon sequestration and microbial respiration.