Biomass heat storage dampens diurnal temperature variations in forests

Observational evidence suggests that compared to non-forested areas, forests have a cooling effect on daytime land surface temperature (LST) and a warming effect on nighttime LST in many regions of the world, thus implying that forests dampen the diurnal temperature range. This feature is not captured by current climate models. Using the Community Land Model 5.0 (CLM5.0), we show that this diurnal behavior can be captured when accounting for biomass heat storage (BHS). The nighttime release of energy absorbed by the vegetation biomass during the day increases both nighttime LST and ambient air temperature in forested regions by more than 1 K. The daytime cooling is weaker than the nighttime warming effect, because the energy uptake by the biomass is compensated by a reduction in the turbulent heat fluxes during day. This diurnal asymmetry of the temperature response to BHS leads to a warming of daily mean temperatures, which is amplified during boreal summer warm extremes. Compared to MODIS, CLM5.0 overestimates the diurnal LST range over forested areas. The inclusion of BHS reduces this bias due to its dampening effect on diurnal LST variations. Further, BHS attenuates the negative bias in the nighttime LST difference of forest minus grassland and cropland, when compared to MODIS observations. These results indicate that it is essential to consider BHS when examining the influence of forests on diurnal temperature variations. BHS should thus be included in land surface models used to assess the climatic consequences of land use changes such as deforestation or afforestation.

A coding error was found and corrected in the implementation of the leaf energy storage, causing an underestimation of this storage. A corrected version of the code can be found at: https://github.com/ RonnyMeier/ctsm/tree/ERL-106754_corrected. For any future work, please use the code from this repository instead of the repository mentioned at the end of section 2.1 on page 3.
After this correction, the diurnal range of the biomass heat flux in figure 1 is increased by 3.98 W m −2 on average (area-weighted mean, figure C1). The impact on temperature, however, is very limited and the main conclusions of the study remain unchanged

Introduction
Forests play a critical role in the climate system by regulating land-atmosphere exchanges of greenhouse gases, energy, and water (e.g. Pongratz et al 2010, de Noblet-Ducoudré et al 2012. The influence of forests on climate can be divided into two categories of processes: (1) the globally acting biogeochemical processes, representing the release or sequestration of greenhouse gases, and (2) the biogeophysical processes, representing direct alterations of the local energy and/or water budget. While the first category of processes appears to be dominant on a global scale (Pongratz et al 2010), the latter alter the climate considerably at local to regional scale and are thus essential in determining the influence of forests on local climate ( Lee et al 2011, Li et al 2015, Bright et al 2017, Duveiller et al 2018.
In the absence of snow, forests are observed to have a local cooling impact on both land surface temperature (LST) and above-canopy air temperature during the day compared to open land (i.e. grassland and cropland; Lee et al 2011, Li et al 2015, Alkama and Cescatti 2016, Duveiller et al 2018. During nighttime, forests tend to be only slightly colder than open land in tropical and subtropical areas and even warmer than open land in the mid-latitudes. Hence, forests tend to dampen the diurnal temperature range (DTR) everywhere but the boreal regions (Lee et al 2011, Duveiller et al 2018. In the mid-latitudes, the sign of the temperature difference of forest minus open land even changes during the diurnal cycle, from a negative Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. daytime difference to a positive one at night (Lee et al 2011, Zhang et al 2014. This feature is missing in all the climate models analyzed in the LUCID, CMIP5 and LUCAS intercomparisons (Lejeune et al 2017, Davin et al 2019. A considerable fraction of the models even exhibits a warming effect of forests on daily maximum 2 m air temperature (T2M) and a cooling effect on daily minimum T2M, opposing observations completely. Such biases could be partly related to the calculation of T2M in models, which does not necessarily correspond the temperature 2 m above the canopy in forests (Meier et al 2018, Winckler et al 2018. However, the nighttime warming by forests is still not represented in CLM4.5 when considering LST, which is more directly comparable to remote sensing observations (Meier et al 2018). This systematic model bias suggests that land surface models are either missing or are poorly representing an important process affecting the energy redistribution at the land surface.
There is a general consensus about the processes leading to the daytime cooling effect of forests compared to open land. Namely, forests have higher surface roughness, which results in stronger turbulent heat fluxes, as well as higher evaporative fraction, which is associated with more evaporative cooling (Davin and de Noblet-Ducoudré 2010, Lee et al 2011, Vanden Broucke et al 2015, Schultz et al 2017. Indeed, reducing the negative bias in the ET difference of forest minus open land in CLM4.5 also attenuated the positive bias in the daytime LST difference between these land cover types (Meier et al 2018). On the other hand, little is known about the mechanism behind the nighttime warming effect of forests in the mid-latitudes. Lee et al (2011) hypothesized that under stable conditions warm air from the planetary boundary layer is transported more effectively towards the land surface over forests, due to their higher surface roughness compared to the one of open land. Vanden Broucke et al (2015) observed that the nighttime warming of forests is not only related to sensible heat fluxes but also to more incoming longwave radiation. They hypothesized that the latter flux could be explained by (1) a stronger greenhouse gas effect over forests due to a moister boundary layer, (2) higher aerosol loading over forests, or (3) the presence of warmer air due to more turbulent mixing over forests (i.e. the higher incoming longwave radiation is a side effect of the stronger turbulent mixing over forests). An additional explanation for the nighttime warming effect of forests could be that energy accumulated in the plant biomass during the day is released into the canopy space during the night (Schultz et al 2017, Meier et al 2018. Due to the high amounts of biomass in forests compared to grasslands or croplands, Biomass heat storage (BHS, i.e. heat storage in the biomass itself) is likely much larger over forests and could thus be a driver behind the nighttime warming effect of forests. Indeed, in situ observational studies indicate that the diurnal amplitude of the biomass heat flux (BHF; i.e. the change of energy stored in the biomass) can be substantial (Aston 1985, Moore and Fisch 1986, McCaughey and Saxton 1988, Meesters and Vugts 1996, Vogt et al 1996, Silberstein et al 2001, Meyers and Hollinger 2004, Oliphant et al 2004, dos Michiles and Gielow 2008, Garai et al 2010, Lindroth et al 2010, Kilinc et al 2012, Burns et al 2015. While the BHF is negligible when integrated over longer time scales, it typically exhibits diurnal amplitudes of 15-75 W m −2 in mature forests (e.g. Oliphant et al 2004, Haverd et al 2007, Lindroth et al 2010, Kilinc et al 2012, Burns et al 2015. A flux of this magnitude appears sufficient to considerably alter diurnal temperature variation over forests. Hence, several land surface model have already included the process BHS (e.g. Verseghy et al 1993, Samuelsson et al 2011, Boone et al 2017.
Still, the relevance of BHS for the local climate has to our knowledge not been assessed at a global scale. In this study, we therefore investigate, how BHS affects diurnal LST and T2M variations over forested regions using the Community Land Model 5.0 (CLM5.0). Swenson et al (2019) recently introduced a scheme simulating energy storage in leaves and stems in CLM5.0. Here, we extend the BHS scheme of Swenson et al (2019), by coupling BHS to the biomass carbon stocks simulated by the model with an active carbon and nitrogen cycle. As a consequence, the leaf area index (LAI), leaf biomass, stem biomass, and vegetation height are diagnosed directly from the simulated vegetation carbon pools instead of relying on uncertain parameter values, thereby improving the internal consistency within the model. We then present a first global estimate of the diurnal amplitude of the BHF and compare the simulated BHF to in situ studies to assess the realism of the BHS scheme. In a second step, we assess how BHS affects temperatures at grid cell level and evaluate the modeled DTR with remote sensing observations from the MODerate resolution Imaging Spectroradiometer (MODIS) system. Finally, we investigate how BHS affects the sensitivity of LST to land cover at sub-grid scale, which we compare to two MODIS-based data sets (Li et   CLM5.0 simulates thermodynamic processes, such as absorption, reflection, and emission of shortwave and infrared radiation, sensible and latent heat fluxes from the vegetation and soil, and heat storage in the soil column. The hydrology of the land surface considers infiltration, runoff, canopy interception, and evapotranspiration, distinguishing between the water, snow, and ice phase. Further, CLM5.0 can represent the exchange of carbon and nitrogen via processes such as photosynthesis, autotrophic and heterotrophic respiration, litterfall, nitrogen deposition, nitrogen fixation, and nitrogen mineralization, denitrification, and fire. Recently, Swenson et al (2019) implemented a BHS scheme to a post-release version of CLM5.0 which allows the representation of heat storage in stems and leaves. An additional vegetation temperature, the stem temperature, was introduced in this scheme, besides the already existing bulk canopy temperature (representing the temperature of the leaves). As a consequence, the energy balance is solved both for the leaf and the stem.
The leaf energy balance is closed at each time step, under the assumption that leaf temperature (T leaf ) is in balance with the leaf energy fluxes given the high surface area/volume ratio (i.e. the leaf temperature is iteratively adjusted until the energy balance is closed): where S leaf  and L leaf  are the net solar and the net longwave radiation absorbed by the leaves, H leaf and λE leaf are the sensible and latent heat fluxes from the leaves, and c leaf is the heat capacity of the leaves. The energy balance of the stem is computed after calculating T leaf as follows: where S stem  and L stem  are the net solar and the net longwave radiation absorbed by the stem, H stem is the sensible heat flux from the stem, c stem is the heat capacity of the stem, and T stem is the stem temperature. Hence, the energy surplus (deficit) of the stem (i.e. the left-hand side of equation (2)) is then taken up (released) by the stem and the stem temperature is adjusted accordingly. The BHF represents the change in energy content of the biomass and is calculated as follows: We introduce a few modifications and different parameter choices compared to the original implementation of Swenson et al (2019). In addition to the LAI, the BHS module requires a number of plant functional type (PFT)-specific parameters. We assume a dry wood density of 500 kg m −3 (Swenson et al 2019). For the fraction of the biomass that is water ( f w ) we use a value of 0.7 for the leaves (Bonan et al 2018) and a value of 0.5 for the stem, which is representative for the outer section of the stem (Herrington 1969). For the tree number density, the value of the closest biome of Crowther et al (2015) is chosen for each PFT. Further, we couple the BHS scheme to the prognostic carbon and nitrogen module in CLM5.0 (referred to as the Bgc-mode in the model community). This allows to calculate the aboveground leaf biomass (ALB, kg m −2 ) and aboveground stem biomass (ASB, kg m −2 ) based on the simulated carbon pools: where C leaf and C stem are the leaf carbon pool and stem carbon pools, which are prognostically simulated by CLM5.0, and b carbon is the ratio of dry biomass to carbon with a constant value of 2 kg/kgC (Bonan et al 2018). This approach enables the calculation of a breast-height diameter (DBH) from the ASB instead of assuming a globally constant value. As in Swenson et al (2019), the resistance to heat transfer between the interior of the tree and the tree surface, r int , is linked to the DBH: where r w is resistance to heat transfer per meter of stem diameter. We use r w as a tuning parameter to achieve realistic BHF values and arrive at the same value of 1000 s m −2 as in Swenson et al (2019). In this study however, r int may vary spatially and temporally as the DBH is linked to the stem carbon pool. Additionally, the Bgc-module computes the LAI, litter area index (SAI), and vegetation height based on the vegetation carbon pools instead of using values estimated from remote sensing observations (Lawrence et al 2018).

Simulation setup
The simulations presented in this study are analyzed over the period of 2002-2010 (since the observational data described in the section below start from 2002 and the atmospheric forcing is available only until 2010). They were run in offline mode, forced by the Global Soil Wetness Project (GSWP3) data set, at 0.5°s patial resolution (Kim 2017). The percentage of the different PFTs within each grid cell is derived from MODIS observations, as described in Lawrence and Chase (2007), and is kept constant throughout the simulation at the coverage of the year 2000 ( figure A1). We use separated soil columns in our simulations to suppress unrealistically large lateral ground heat fluxes between different PFTs (Meier et al 2018, Schultz et al 2016. To isolate the effect of BHS, we run a control simulation with no BHS (CLM-CTL). The second simulation, CLM-BHS, is run in the exact same configuration but with an active BHS scheme. As in Swenson et al (2019), we raise the upper cap of the Monin-Obukhov stability parameter in the surface layer from 0.5 to 100 to assess the full impact of BHS in both simulations (thereby making the atmospheric stability virtually unconstrained). Both the CLM-CTL and CLM-BHS simulations start from the same initial state and are run for the years 1997 to 2010 (note that this experimental design mutes feedbacks arising from BHS on the simulated carbon stocks themselves). There is therefore an additional 5 years of spin up before the analysis period of 2002-2010. The initial state in 1997 is retrieved by running the model for 146 years at 0.5°spatial resolution, starting from an already spun-up, pre-industrial state, which was interpolated from the original resolution of 0.9°×1.25°to 0.5°×0.5°resolution. For this additional spin up, we first cycle five times through the atmospheric forcing data of the years 1901-1910 for the period of 1851-1900. After, we force the model with the reanalysis data of the years 1901-1996. During this spin up, BHS remains inactive. Figure A2 illustrates the resulting aboveground biomass (AGB) of different land cover types.

Observational data
We evaluate the modeled LST with MODIS observations both at grid cell level and in terms of the sub-grid difference between forest minus open land. The MODIS instruments aboard the satellites Aqua and Terra provide LST measurements at approximately 01:30 and 13:30 solar time at 0.05°×0.05°resolution. For the evaluation at grid cell level we use the data of the monthly MYD11C3 data product from July 2002 to December 2010. After masking out observations with a reported LST error estimate larger than 1 K and/or an emissivity error estimate larger than 0.01 (as in Li et al 2015) and discarding pixels with a land fraction lower than 80%, the original data is averaged to the model resolution of 0.5°×0.5°. The resulting monthly values are then used to derive a multi-year monthly average. From this multi-year monthly average we then calculate the multi-year total average, excluding pixels with no valid data for at least month. The model performance in terms of its representation of the local LST difference of forest minus open land is evaluated using the two MODIS-based data sets of Li et al (2015) and Duveiller et al (2018), subsequently called MOD-Li15 and MOD-Du18. These data sets compare LST over different land covers within moving windows of 9 by 5 pixels or 5 by 5 pixels to infer a LST difference of forest minus open land at roughly 01:30 and 13:30 solar time. Given the relatively high resolution of the two data sets, it is a reasonable assumption that the different pixels within the moving window experience similar climatic conditions. Thus, the data sets are assumed to capture the local impact of a conversion from forest to open land.
MOD-Li15 and MOD-Du18 differ in methodology, time frame, and the MODIS products utilized. MOD-Li15 uses the MYD11C2 product extracting data from the period 2002-2012. On the other hand, MOD -Du18 employs MYD11C3 data from 2008 to 2012. The MYD11C2 and MYD11C3 products are both aggregated from daily LST observations (MYD11C1) to 8 d average values in the case of MYD11C2 and monthly average values in the case of MYD11C3 (Wan and Hulley 2015a, 2015b). Although the two data sets exhibit a comparably high spatial resolution, the observed pixel usually still comprises different land cover types. Thus, it is challenging to isolate a LST over forest or open land only. MOD-Li15 defined a threshold of 80% in the coverage of a certain pixel by forest or open land for the pixel to be classified as forest or open land, respectively, and calculated the difference between the forest minus the open land pixels within search windows of 9 by 5 pixels. The LST of the different pixels was corrected for elevation differences using a lapse rate which was inferred from the MODIS observations. Further, the comparison samples were masked for an elevation difference of less then 500 m. The more recent MOD-Du18 data set on the other hand used a multi linear regression model within windows of 5 by 5 pixels to establish a relationship between the fraction of different vegetation types and LST. Additionally, Duveiller et al (2018) applied stricter criteria to mask out pixels with too high elevation variability but did not apply an elevation adjustment to the LST data.

Model evaluation
MOD-Li15 and MOD-DU18 depict the local LST difference of forest minus open land under similar atmospheric forcing. We extract a comparable signal from our simulations by calculating the sub-grid difference of forest minus open land in the variables of interest following the approach of Malyshev et al (2015). In this approach the average in a certain variable over the forest tiles within a grid cell is subtracted from the average over the open land tiles within the same grid cell (i.e. under the same atmospheric forcing). Note that this difference is comparable to the signal of an open-land-to-forest transition (i.e. afforestation or reforestation) rather than a deforestation signal. A more detailed description of this methodology can be found in Meier et al (2018).
The LST is calculated as the weighted average between the leaf temperature (T leaf ) and the ground temperature (T grnd ): The vegetation emissivity, e v , is calculated as in equation (4.20) in Lawrence et al (2018): here L and S are the LAI and stem area index, respectively, and m is the average inverse optical depth for longwave radiation (with a value of 1). T grnd is inferred from the snow temperature (T snow ), the temperature of the top soil layer (T soil ), and the surface water temperature (T H O 2 ) as follows:  (2018), we assume that the LST difference at 01:30 and 13:30 reasonably represent the daily minimum and maximum LST difference (we find that this assumption is reasonable in the context of CLM5.0 output). Hence, ΔLST min (f-o) and ΔLST max (f-o) refer to the LST difference of forest minus open land at 01:30 and 13:30, respectively. In addition, we analyze the effect of forests on the DTR, which can be calculated as the difference of ΔLST max (f-o) minus ΔLST min (f-o) (Duveiller et al 2018). Accordingly, this variable is called ΔDTR(f-o). While the advantage of LST is that there exist observations with a global coverage, T2M is used more frequently as a temperature metric and is thus more readily interpreted. We therefore also analyze the impact of BHS on T2M at grid cell level (equation (5.58) in Lawrence et al 2018). T2M is a diagnostic variable in CLM5.0, representing the air temperature 2 m above the apparent sink of sensible heat (i.e. 2 m above the sum of the roughness length for sensible heat and the displacement height). Note that this height normally lies below the canopy for forests but above the canopy for shorter vegetation, such as grassland or cropland. However, the spatial variability in the observations is considerable and much higher than in the model, which can be related to several factors. First, the observational studies utilize differing methodologies to estimate BHS. Second, those studies are conducted for a limited time period and are thus impacted strongly by the meteorological conditions during the measurement campaign. Finally, the geometry of the vegetation at the site has an impact on BHS. An example for the latter factor is the study of Kilinc et al (2012), which was conducted at a site with relatively few but enormous trees (average mass of 29Mg per tree), thereby resulting in a relatively small BHF when compared to the large amount of AGB. The opposite is the case in Burns et al (2015), where the exceptionally high tree number density of 0.4 m −2 is contributing to a comparably high BHF range. As can be seen in figure 1(b), CLM-BHS can not capture such exceptionally large BHFs, which is related to two reasons: first, the land surface characterization in the simulation is not necessarily representative for a specific site.

Impact of BHS at grid cell level
The comparison of the CLM-BHS and the CLM-CTL simulations during boreal summer (JJA) shows that BHS has a considerable impact on T2M in areas with a large fraction covered by forests (see figures 2 and A1(a)). The diurnal T2M impact induced by BHS is asymmetric: while the daily minimum T2M often increases by more than 1 K, the daily maximum cooling effect of BHS rarely exceeds 0.1 K (figures 2(a) and (b)). As a consequence, daily average T2M tends to be higher in CLM-BHS. This warming frequently exceeds 0.4 K both for daily mean T2M and daily mean LST (see figure A3 for LST). BHS also has important implications regarding extreme conditions. In fact, the warming effect of BHS during the 5% hottest days of JJA is stronger in most regions than the effect on the mean, frequently exceeding 1 K over forested areas (figures 2(d) and A4(c)). In addition, the implications of BHS vary with season. For example, the daytime cooling effect of BHS is more pronounced during winter in the northern latitudes, resulting in a weaker increase of daily mean T2M and daily mean LST due to BHS during this season (figures A5 and A6). The temperature impact of BHS displayed in figures 2, A3, A5, and A6 originates mostly from the forest PFTs. When considering only the forest PFTs, BHS often increases the daily minimum T2M by more than 2 K ( figure 3(b)). On the other hand, the impact of BHS over the open land PFTs is marginal, due to the negligible storage capacity of leaves ( figure 3(d)).
The increase in nighttime temperatures and decrease in daytime temperatures resulting from BHS,   4(a)). In comparison to MODIS, CLM-CTL tends to exhibit negative bias in DTR over regions with low forest coverage, and a positive bias over regions with higher forest coverage both in the tropics and subtropics as well as in the mid-latitudes (figures 4(b) and 5). This finding is supported by a pattern correlation coefficient of 0.48 between the DTR bias of CLM-CTL ( figure 4(b)) and the forest fraction in the land cover data ( figure A1(a)). BHS reduces part   Points with a mean which is insignificantly different from zero in a two-sided t-test at 95% confidence level are marked with a black dot. All data from the 2002 to 2010 analysis period corresponding to a given latitude and a given month are pooled to derive the sample set for the test. The numbers next to the titles are the area-weighted spatiotemporal root-mean-squared deviation of the respective simulation against the MOD-Li15 data set (Meier et al 2018). Panel (d) shows the zonal annual mean of MOD-Li15 (green, range between the 10th and 90th percentiles in gray), CLM-CTL (blue, range between the 10th and 90th percentiles in blue), and CLM-BHS (red, range between the 10th and 90th percentiles in orange). Note that on this subfigure results have been smoothed latitudinally with a simple moving average over 4°. The same with MOD-Du18 in panels (e)-(h). For values over different latitudinal bands see figure A7. of the positive DTR bias over forested regions, by dampening the DTR (figures 4(c)-(e)). The relative reduction of the DTR bias in forested regions is most distinct in the northern mid-latitudes, where CLM-CTL typically overestimates the DTR by 1-5 K. The relative improvement is more moderate in tropical regions, as CLM-CTL frequently overestimates the DTR by more than 5 K, which is not alleviated completely by BHS. Therefore, the overall tendency to overestimate the DTR in forested regions and underestimate the DTR in sparsely vegetated regions partly persists in CLM-BHS, despite these improvements (figure 5). ΔLST min (f-o) is on average in the order of 1 K higher in CLM-BHS as compared to CLM-CTL, resulting in a considerable reduction of the global root-meansquared deviation (RMSD) in comparison to both observational data sets (figure 6). In contrast to the clearly negative ΔLST min (f-o) in CLM-CTL, ΔLST min (f-o) of CLM-BHS is close to zero in tropical regions as observed in the MODIS-based products. However, there is still some disagreement on the nighttime signal of forests between CLM-BHS and the observations in the extra-tropics. The nighttime warming by forests between 20°N and 55°N is sometimes too weak in the BHS simulation, especially in comparison to MOD-Li15 (figures 6(d) and (h)). However, the nighttime warming by forests at these latitudes is less pronounced in MOD-Du18 than in MOD-Li15, highlighting a certain degree of uncertainty in the observations. Similarly, CLM does not capture the distinctly positive ΔLST min (f-o) between 25°S and 40°S which is observed in MOD-Li15, whereas the BHS simulation agrees with MOD-Du18 overall at these latitudes. At high-latitudes on the other hand, BHS impairs the already positively-biased ΔLST min (f-o). During the day, BHS consistently decreases ΔLST max (f-o) by about 0.1 K, which causes only a small change in RMSD (figure A8). In offline simulations, BHS appears therefore mostly relevant for the land cover sensitivity of the model during night, whereas its effect is less distinct during the day.

Diurnal asymmetry of BHS impact on temperature
We found that BHS affects daytime temperature much less than nighttime temperature, although roughly the same amount of energy absorbed by the vegetation during the day is released during night (i.e. the daily average BHF is close to zero). We hypothesize that this diurnal asymmetry of the temperature impact is related to a different structure of the surface layer during daytime compared to the nighttime. During the day, the surface layer tends to be unstable and the turbulent heat fluxes are directed from the land surface towards the atmosphere. An initial reduction of the surface temperature due to the energy uptake by the vegetation decreases the surface layer instability and thus also the turbulent heat fluxes. The lower turbulent heat fluxes lead to a higher surface temperature (because less energy is transported away from the surface) and form therefore a negative feedback to the initial temperature decrease ( figure 7). In other words, the energy uptake by BHS during day is mainly compensated by a reduction of the turbulent heat fluxes. During the late night on the other hand, the surface layer is often stable. Therefore, the energy released from the biomass is compensated less efficiently by increased turbulent heat fluxes. Further, the sensible heat flux can be directed towards the land surface during the late night. In this case, the increase in surface temperature due to BHS reduces the surface layer stability and the suppression of turbulent mixing, resulting in more intrusion of relatively warm boundary layer towards the surface. Consequently, the feedback from the surface layer to the initial surface temperature increase can even be an additional warming of the surface ( figure 7).
This behavior can be seen in point scale simulations in the tropics at the site of dos Michiles and Gielow (2008). As long as there is solar irradiance, the energy taken up/released by the vegetation is compensated for to a large extent by a reduction/increase in turbulent heat fluxes (figure 8(a), compensation calculated from the difference of the respective flux between a CLM-CTL-like simulation and that of a CLM-BHSlike simulation). Figure 8(b) displays the fraction of the BHF that is compensated by decreased/increased emission of longwave radiation by the land surface ( f rad ) against the gradient between the atmospheric temperature and the vegetation temperature (as an indicator for near-surface stability). f rad remains low and surprisingly constant throughout the day, indicating an only weak reduction of the LST when BHS is positive. Once the Sun sets, the surface layer grows more stable, thereby inhibiting the compensation by the turbulent heat fluxes. Therefore, f rad gradually increases during early night up to roughly 0.6 and remains relatively constant during the second half of the night (accompanied by a strong nighttime increase in LST). As hypothesized before, sensible heat even provides a small positive feedback during the late night, due to the lower surface stability resulting from BHS (The dark red area in figure 8 a from 2:00 to 7:00 am indicates that sensible heat has a warming effect compensating part of the emitted longwave radiation effect). A similar diurnal cycle of f rad occurs in the mid-latitudes at the site of Lindroth et al (2010) during the summer months ( figure A9).
The strong dependency of f rad on the surface layer stability can also explain why the daytime cooling effect of BHS in the mid-latitudes is more pronounced during winter. During this season the daytime compensation of BHF by reduced turbulent heat fluxes is less efficient as stable conditions occur more frequently than during summer (Chan and Wood 2013). Hence, a larger proportion of the BHF is compensated by reduced emission of longwave radiation. At the site of Lindroth et al (2010), f rad often exceeds 0.1 during the mid-latitudinal winter even when the BHF is positive (i.e. the biomass takes up energy), which is normally not the case during summer (see figures A9 and A10). Hence, the characteristic diurnal cycle of f rad which emerges in the tropics and during summer in the mid-latitudes, does not develop as clearly during DJF.
According to our modeling results, the diurnally asymmetric temperature difference between forest and nearby open land, which was found in observational studies (Lee et al 2011, Li et al 2015, Vanden Broucke et al 2015, Schultz et al 2017, Duveiller et al 2018, can thus be explained by differing relevant processes during day and night. In various studies it was found that the differences in the evaporative fraction and albedo are more relevant for the temperature difference of forest minus open land during day (Liu et al 2005, Vanden Broucke et al 2015, Meier et al 2018, whereas it appears that BHS mainly affects nighttime temperatures (Swenson et al 2019 and this study). However, another part from this asymmetry is likely related to the higher turbulent heat fluxes over forests due to their high surface roughness (Schultz et al 2017). The sensible heat flux can thus have a cooling effect over forests during day and a warming effect during night compared to open land, if this flux

Limitations and knowledge gaps
While this study demonstrates the importance of BHS on a global scale, it could not address several aspects. The evaluation of the DTR in the CLM5.0 simulations with MODIS at grid-cell-level is to some degree an unfair comparison, as part of the biases could be the result of discrepancies between the atmospheric forcing used in our simulations and the actual conditions during the MODIS observations. As a consequence, even a perfect land surface model could not achieve a perfect agreement with the observations. Unfortunately, it is not straight-forward to assess the disagreement between MODIS and GSWP3 as LST is not directly comparable to atmospheric temperature at the lowest level. Yet, we find that the DTR in the MODIS LST data often considerably exceeds the DTR in the GSWP3 forcing, at locations where CLM5.0 underestimates the DTR in LST (see figures 4 and A11). Also, the MODIS observations can be retrieved only under cloud-free conditions. Unfortunately, the GSWP3 forcing does not contain any cloud cover information, hence making it difficult to mask out cloudy days in the model output. Further, the disagreement in DTR between CLM5.0 and MODIS could also be related to bare soil or open land PFTs, where BHS is unimportant due to the small (or nonexistent) biomass of these land cover types. In fact, our analysis revealed that the DTR is often underestimated by CLM5.0 in sub-tropical and mid-latitudinal regions, which are mostly covered by bare soil or open land (see figures 4(b) and A1). Hence, the missing process of BHS is likely not the only source for the difference between our CLM5.0 simulations and the MODIS DTR.
At least two potential feedbacks to BHS cannot be addressed with our simulation set up. First, land-only simulations inherently mute the feedbacks from the atmosphere. We found that approximately 90% of the energy-uptake by the biomass is compensated by a reduction of the turbulent heat fluxes during the day. BHS therefore reduces the daytime energy input from the land surface into the boundary layer. This reduced energy input could result in dynamical feedbacks in the boundary layer. It seems therefore necessary to analyze simulations coupled to the atmosphere to capture the full effect of BHS. Another feedback of BHS that we do not assess in this study are the changes in vegetation structure induced by BHS. The alteration of the local climate due to BHS affects primary productivity and respiration (not shown). This can not only alter the vegetation phenology (e.g. vegetation height, LAI, biomass) but also the carbon budget of the land surface. We therefore encourage further research in this direction.
Finally, there remain several aspects of energy storage at the land surface that could be further improved in the CLM modeling framework. First, the model we use restricts the energy storage at the land surface to the ground and biomass energy storage, neglecting the sensible and latent heat stored in the canopy air space, which can be of comparable magnitude in observations as BHS (e.g. dos Michiles and Gielow 2008, Lindroth et al 2010, Kilinc et al 2012). Bonan et al (2018) find that introducing a multi-layer canopy model, which represents the roughness sublayer and thus the canopy air space, did reduce but not alleviate completely the positive bias in the DTR in forests. However, this model version did not account for BHS in the stem. Hence, the combination of a multi-layer canopy model and heat storage in the stem could potentially alleviate the overestimation of the DTR in forests completely. Second, we assume a uniform stem temperature in our model, which is of course not the case in reality. The existing BHS parameterization could thus be refined with the 'force-restore' method, which represents two stem layers (Haverd et al 2007), or even by including multiple stem layers, as proposed in the 'analog model' by Herrington (1969). Further, a number of parameters such as the tree number density are assumed to be globally constant for the different PFTs. In reality however, such parameters can vary considerably within a PFT. Hence, it could be beneficial to include spatially explicit data if available from observations.

Implications for LUC impacts on climate
Recently, the local and regional climate impact of landuse changes such as deforestation has been a heavily discussed topic in the literature. These climatic impacts are often assessed by contrasting climate model simulations that differ in land cover. However, all of the CMIP5, LUCID, and LUCAS regional climate models miss or underestimate the observed nighttime cooling effect and daytime warming effect associated with deforestation (Lejeune et al 2017, Davin et al 2019. This indicates, that current climate models miss an important part of the climate signal of forest-related land-use changes. Up to now, BHS has often been neglected in global climate models, although observed BHFs are of non-negligible magnitude across the diurnal cycle. Given our results of including BHS in CLM5.0, we reach the following three main conclusions regarding the role of BHS in the climate system: (1) BHS dampens the DTR in forested regions, thereby improving the model agreement with MODIS observations. (2) BHS affects LST stronger under stable surface layer conditions as compared to unstable conditions. As a consequence, the BHS-induced temperature increase during the night tends to exceed the temperature reduction during the day.
(3) The effect of BHS is especially pronounced during warm extremes. These conclusions indicate that the representation of BHFs is important, especially in the case of diurnal temperature variations and/or temperature extremes in forests. Climate models used to assess the biogeophysical impact of forest-related landuse changes should thus consider including BHS in their land surface component.          Table 1. Comparison of aboveground biomass (AGB) (kg m −2 ) and diurnal range of the biomass heat flux (BHF) (W m −2 ) in in situ observations and the CLM-BHS simulation. The months during which measurements were made go from January (1) to December (12). Values marked with a asterisk were estimated by reading from graphs displaying the diurnal cycle of biomass heat fluxes. For some studies the respective plant functional type was not available at the location of the observational study. Therefore, the value of the closest PFT was used instead when comparing to the observational studies of Haverd et al (2007), Silberstein et al (2001), Meesters and Vugts (1996), and Garai et al (2010