Simulated single-layer forest canopies delay Northern Hemisphere snowmelt

Abstract. Single-layer vegetation schemes in modern land surface models have been found to overestimate diurnal cycles in longwave radiation beneath forest canopies. This study introduces an empirical correction, based on forest-stand-scale simulations, which reduces diurnal cycles of sub-canopy longwave radiation. The correction is subsequently implemented in land-only simulations of the Community Land Model version 4.5 (CLM4.5) in order to assess the impact on snow cover. Nighttime underestimations of sub-canopy longwave radiation outweigh daytime overestimations, which leads to underestimated averages over the snow cover season. As a result, snow temperatures are underestimated and snowmelt is delayed in CLM4.5 across evergreen boreal forests. Comparison with global observations confirms this delay and its reduction by correction of sub-canopy longwave radiation. Increasing insolation and day length change the impact of overestimated diurnal cycles on daily average sub-canopy longwave radiation throughout the snowmelt season. Consequently, delay of snowmelt in land-only simulations is more substantial where snowmelt occurs early.



Introduction
Forest canopy cover modulates longwave radiation received by the ground, which consequently differs from atmospheric longwave radiation. This process is called longwave en-25 hancement and has been shown to result in substantial postive net longwave radiation of the surface when snow cover is prevalent, especially under clear skies and during snowmelt (Webster et al., 2016). In contrast, net longwave radiation fluxes are typically negative for snow under clear-sky con-30 ditions in unforested areas as has been observed for ever-green Canadian boreal forests (Harding and Pomeroy, 1996;Ellis et al., 2010). Moreover, forest cover has been reported to enhance snowmelt for a subarctic open woodland during overcast days and early in the snowmelt season (Woo and 35 Giesbrecht, 2000). However, the impact of forest coverage on snowmelt varies regionally as a function of forest density and meteorological conditions, with the importance of shortwave and longwave radiation changing throughout the snowmelt season (Sicart et al., 2004;Lundquist et al., 2013). 40 Meteorological conditions control longwave enhancement as clear skies increase insolation and thereby vegetation temperature while radiative temperature of the sky is reduced (Sicart et al., 2004;Lundquist et al., 2013;Todt et al., 2018). Therefore, values for longwave enhancement, i.e. the ratio 45 of below-canopy to above-canopy longwave radiation, are higher under clear skies but close to 1 or even smaller for overcast conditions due to similar vegetation temperature and radiative temperature of the sky. Vegetation density impacts longwave enhancement by scaling the respective contribu-50 tions of vegetation and atmosphere to sub-canopy longwave radiation as well as by governing the impact of meteorological forcing on vegetation temperatures (Todt et al., 2018). While observations have shown trunks heating up due to insolation and emission of longwave radiation consequently in-55 creasing (Rowlands et al., 2002;Pomeroy et al., 2009), diurnal variations in tree temperatures depend on exposure to insolation and thus vegetation density (Webster et al., 2016).
About a fifth of seasonally snow-covered land over the Northern Hemisphere is covered by boreal forests (Rutter 60 et al., 2009), indicating the process of longwave enhancement affects a substantial fraction of global snow cover. Considerable challenges persist in the representation of snow cover and snowmelt in the current generation of climate models as historical simulations from Climate Model In-65 tercomparison Project's fifth phase (CMIP5) underestimate observed trends and interannual variability of spring snow cover extent (SCE) (Derksen and Brown, 2012;Brutel-Vuilmet et al., 2013;Rupp et al., 2013;Mudryk et al., 2014;Thackeray et al., 2016)

. Snow Model Intercomparison
Project's second phase (SnowMIP2) identified less skill in modelling snow for forested than for open sites, which was attributed to complex processes between atmosphere, snow, and vegetation (Essery et al., 2009;Rutter et al., 2009).
Among models displaying deficiencies in simulating snow 10 cover evolution across boreal forests is the Community Land Model (CLM) version 4 and its parent Community Climate System Model version 4 (Thackeray et al., 2014(Thackeray et al., , 2015. CLM uses a one-layer vegetation scheme and CLM version 4.5 (CLM4.5) has been found to show deficiencies in simulation 15 of sub-canopy longwave radiation and longwave enhancement with overestimated diurnal cycles under clear-sky conditions (Todt et al., 2018). Similar issues have been mitigated in CLM4.5 by subdividing the roughness layer (Bonan et al., 2018) and in SNOWPACK, a one-dimensional snow cover 20 model, by partitioning the canopy into two layers with separate energy balances and consequently separate vegetation temperatures (Gouttevin et al., 2015), which results in different longwave radiation fluxes emitted upward and downward from the vegetation.

25
In order to avoid implementing multiple canopy layers in a global land model and associated computational costs, this study presents an alternative guided by the effect of separate vegetation layers on sub-canopy longwave radiation. A correction to sub-canopy longwave radiation is implemented 30 in CLM4.5 to reduce overestimated diurnal cycles, damping variations in longwave radiation emitted downward and, consequently, increasing variations in longwave radiation emitted upward. While simulation of sub-canopy longwave radiation and longwave enhancement by land surface models had 35 so far been assessed using forest stand-scale forcing and evaluation data, this study uses land-only simulations of CLM4.5 and snow-off dates derived from global observations of snow water equivalent (SWE) to assess the impact of overestimated diurnal cycles in sub-canopy longwave radiation on 40 simulated global snow cover and snowmelt. Therefore, this study has three objectives: i. To develop a correction of sub-canopy longwave radiation simulated by single-layer vegetation in CLM4.5; ii. To evaluate the effect of this correction on simulated di-45 urnal cycles and daily averages of sub-canopy longwave radiation; iii. To quantify the impact of corrected sub-canopy longwave radiation on snow cover and snowmelt across the Northern Hemisphere.

50
Section 2 presents methodological details about treatment of sub-canopy longwave radiation in CLM4.5, the physical basis for the empirical correction, configuration of global land-only simulations, and calculation of snow-off date from SWE observations. Calculation of correction factors is de-55 tailed in Sect. 3, and their impacts on the simulated energy balance and seasonal cycle of snow cover are presented in Sect. 4. We conclude with a brief discussion in Sect. 5.

Methods
2.1 Sub-canopy longwave radiation in CLM4.5 60 Vegetation in CLM4.5 is parameterized as a single layer using a "big-leaf" approach (Oleson et al., 2013). Sub-canopy longwave radiation is calculated as the sum of atmospheric longwave radiation LW atm and longwave radiation emitted by vegetation LW veg , weighted by vegetation emissivity ε v : 65 using the Stefan-Boltzmann law with Stefan-Boltzmann constant σ = 5.67 10 −8 Wm −2 K −4 and vegetation temperature T v . Vegetation temperature is calculated based on an energy 70 balance, net radiation minus turbulent heat fluxes. Radiative transfer of direct and diffuse shortwave radiation is calculated via a two-stream approximation (Sellers, 1985) considering one reflection from ground to canopy. Net longwave radiation is calculated from atmospheric longwave ra-75 diation, vegetation temperature, and (ground) surface temperature and determined by vegetation emissivity and emissivity of the ground. Calculation of turbulent heat fluxes in CLM4.5 is based on Monin-Obukhov similarity theory and described by Oleson et al. (2013). Vegetation emissivity de-80 pends on Leaf Area Index (LAI) and Stem Area Index (SAI) and is calculated as This parameter is a combination of emissivity in the physical sense and a weighing parameter based on vegetation density, 85 however, we will stick to this denomination here for consistency with the nomenclature of the technical description of CLM4.5 (Oleson et al., 2013). CLM4.5 subdivides grid cells based on land units, e.g. vegetated, glacier, or urban, and vegetated land units based 90 on Plant Functional Types (PFTs), with up to 16 possible PFTs as well as bare soil. Sub-canopy longwave radiation is calculated for each PFT present in a grid cell, with separate values of LAI, SAI, and vegetation temperature for each PFT. All PFTs within one vegetated land unit share a single 95 column of snow and soil, so that fluxes from vegetation to the ground are weighted averages over all PFTs. Consequently, changes in fluxes from an individual PFT affect snow cover beneath every PFT in a particular vegetated land unit.

Correction of sub-canopy longwave radiation in
CLM4.5 For this study, a correction factor f corr was implemented in CLM4.5 to reduce unphysical diurnal variations in subcanopy longwave radiation. As atmospheric longwave radiation is an input variable to CLM4.5, from either forcing datasets or the atmospheric component of CESM (Community Atmosphere Model, CAM), correction factors were used to scale longwave radiation emitted from vegetation: Conceptually, correction factors represent a vegetation structure consisting of multiple individual layers, so that longwave radiation fluxes emitted upward and downward from the vegetation are no longer equal by design. In a multi-layer canopy configuration, the uppermost layer contributes most 15 to longwave radiation emitted upward to the atmosphere and directly absorbs incoming longwave and shortwave radiation fluxes. Conversely, the lowest layer contributes most to longwave radiation emitted downward to the surface, but is insulated from atmospheric fluxes by the canopy layers above.

20
Using this multi-layer canopy configuration as a guideline, longwave radiation emitted by the canopy was partitioned asymmetrically upward and downward in CLM4.5. The resulting above-canopy longwave radiation flux to the atmosphere was calculated as g with emissivity of the ground ε g and temperature of the 30 ground T g . Ground emissivity in CLM4.5 is calculated as a weighted sum of emissivities of snow (0.97) and soil (0.96), weighted by the fraction of snow covering a grid cell. In Eq. (4), the first term represents atmospheric longwave radiation transmitted through the vegetation, reflected by the 35 ground, and transmitted through the vegetation to the atmosphere; the second term represents longwave radiation emitted from the vegetation reaching the atmosphere; and the third term represents longwave radiation emitted by the ground and transmitted through the vegetation to the atmo-40 sphere. The second term combines longwave radiation emitted by the vegetation directly to the atmosphere (first term in parentheses) and longwave radiation emitted downward from the vegetation, reflected by the ground, and transmitted through the vegetation to the atmosphere (second term in 45 parentheses). For f corr > 1, LW above decreases as reduction in the first term in parentheses (2−f corr ) outweighs increase in the second term in parentheses (1 − ε v ) (1 − ε g ) f corr , while LW sub in Eq. (3) increases. For f corr < 1, LW sub decreases and LW above increases. Note that the sum of LW sub 50 and LW above was not changed by the introduction of f corr , which guaranteed conservation of energy. The calculation of vegetation temperature in CLM4.5 was not altered by this approach, so that the temperature of the single vegetation layer represented an average of multiple (theoretical) layers sug-55 gested by asymmetrical upward and downward longwave radiation fluxes.

Global offline simulations with CLM4.5
Offline simulations of CLM4.5 were forced by prescribed atmospheric data, using the CRUNCEP version 7 data set 60 covering 1981 to 2016 and thus snow seasons 1981/82 to 2015/16 (Viovy, 2018). The impact of correction factors on longwave enhancement, snow cover, and snowmelt was assessed by comparing two simulations, a control run (henceforth, CTRL) and a run in which correction factors were im-65 plemented (henceforth, CORR). Correction factors were applied to evergreen needleleaf trees in CLM4.5 as given in Eq.
(3) and Eq. (4). Two PFTs, Needleleaf Evergreen Boreal Trees (NEBTs) and Needleleaf Evergreen Temperate Trees (NETTs), represent evergreen forests across snow-covered 70 areas in CLM4.5 and grid-cell coverage by these two PFTs is shown in Fig. 1a. Plant Area Index (PAI), the sum of LAI and SAI, is shown in Fig. 1b as a weighted average of NEBTs and NETTs.

75
A blended data set of five global observation-based SWE products (henceforth, Blended-5) covering the period 1981 to 2010 (Mudryk et al., 2015) was used to estimate snow-off dates across the Northern Hemisphere and evaluate simulation of snowmelt in CTRL and CORR. In contrast to simula-80 tions, observations display snow persisting for physically unrealistical durations, which necessitates a SWE threshold to estimate snow-off dates (Krinner et al., 2018). While Mudryk et al. (2017) and Krinner et al. (2018) used thresholds of 4mm and 5mm, respectively, for estimates of spatial snow 85 cover extent, a smaller SWE value was necessary to represent the precise timing of meltout within individual grid cells. A threshold of 1mm was used in this study to define meltout for the Blended-5 mean, and snow-off date was defined as the first day of a year for which SWE did not exceed this 90 threshold. Sensitivity of snow-off dates to threshold values was tested for the range 0.5mm to 4mm, however, the overall conclusions of this study are unchanged for different thresholds.

Calculation of correction factors
95 Todt et al. (2018) created a "toy model", which utilized forest stand-scale forcing data to evaluate sub-canopy longwave radiation in CLM4.5 and revealed systematic simulation errors that depend on meteorological conditions. These me-teorological conditions were categorized via insolation and cloudiness represented by effective emissivity of the sky, which is calculated as using air temperature T air . Based on those stand-scale sim-5 ulations, this study calculated correction factor f corr from ε sky and insolation SW in as Coefficients b 0,...,3 relate to the intercept of the equation, ε sky , insolation, and interaction of ε sky and insolation, re-10 spectively, and were calculated via multiple linear regression from stand-scale simulation errors expressed as ratios ( Fig. 2) and observations of ε sky and insolation at forest stands listed in Table 1. Consequently, correction factors were calculated as inverses of these ratios to scale longwave 15 radiation in CLM4.5. For example, if stand-scale simulations revealed an overestimation of longwave radiation by 25% for particular values of ε sky and SW in , correction factors in global simulations would be 1.25 −1 = 0.8 for the same meteorological conditions.

20
As CLM4.5 only simulates longwave radiation emitted from vegetation, simulation errors were calculated for LW veg that was derived from sub-canopy longwave radiation via Eq. (1) and Eq. (2) using measurements of atmospheric longwave radiation and PAI given in Table 1. Error 25 ratios as a function of ε sky and insolation as well as estimates based on regression coefficients are shown in Fig. 2. Nighttime estimates are a linear function of ε sky as SW in = 0, while daytime estimates include potential non-linear interactions of ε sky and SW in . Both daytime and nighttime sim-30 ulation errors generally increase in magnitude with clearer skies.
Regression coefficients as outlined in Eq. (6) are shown in Fig. 3 for every site and season, differentiated for day and night. Intercept b 0 and regression coefficient for ε sky b 1 agree 35 in sign for all sites and agree in magnitude for all sites except Abisko (panels a, b, d, e), with little interannual variability for the two sites with multiple years of data (Alptal and Seehornwald). In contrast to Abisko, b 0 and b 1 for the deciduous forest at Cherskiy are similar to those for evergreen sites 40 Alptal, Seehornwald, and Sodankylä despite featuring different vegetation type, structure, and density. Regression coefficients involving insolation agree in sign but differ in magnitude among Alptal, Cherskiy, Seehornwald, and Sodankylä (panels c and f), with similar values for the latter two sites 45 due to little interannual variability for Seehornwald. In contrast, interannual variability is large for Alptal with higher magnitudes for all four years combined compared to Seehornwald and Sodankylä, while magnitudes are smallest for Cherskiy. For Abisko, five out of six regression coefficients 50 display smallest magnitudes due to deciduous vegetation and consequently low vegetation density as well as smaller simulation errors compared to other sites (Todt et al., 2018). Overall, uncertainties are largest for Abisko due to a short evaluation period, with no regression coefficient being significantly 55 different from zero (or one as in the case of intercept b 0 ).
For implementation in global simulation CORR, regression coefficients were calculated based on one season each of Alptal, Cherskiy, Seehornwald, and Sodankylä in order to balance dense and sparser sites. Despite featuring a decid-60 uous PFT, Cherskiy was included as regression coefficients are similar to evergreen sites. Individual seasons for Alptal, 2005, andSeehornwald, 2009, were chosen based on similarity of regression coefficients to those for all years combined of the respective site. Regression coefficients for these four 65 sites combined are shown as red lines in Fig. 3. Estimates of simulation errors based on these regression coefficients are shown in Fig. 2 and explain 60% of variance in nighttime errors and 59% of variance in daytime errors.
4 Effect of correction in global simulations of CLM4.5 70

Sub-canopy longwave radiation -case study Alptal, Switzerland
For the location of Alptal, in contrast to other forest stands used in this study, forest stand and CLM4.5 grid cell feature similarly high vegetation densities (PAIs of 4.1 m 2 m −2 75 and 3.7 m 2 m −2 , respectively) and thus similar vegetation emissivities ε v (0.983 and 0.975, respectively). This allows for a comparison of diurnal cycles of sub-canopy longwave radiation as well as longwave enhancement between standscale measurements and offline simulations. Implementation 80 of correction factors in CLM4.5 results in decreased subcanopy longwave radiation during daytime and increased sub-canopy longwave radiation during nighttime, thereby reducing diurnal cycles. For the grid cell representing Alptal, diurnal ranges decrease from about 70 Wm −2 to about 30 85 Wm −2 during snowmelt season ( Fig. 4a and Fig. 4b).
Observations at the forest stand show an average diurnal range of about 15 Wm −2 during snowmelt season. Simulations and observations display a similar range of intraseasonal variability but do not agree in evolution and daily average of 90 sub-canopy longwave radiation. Implementation of correction factors increases average sub-canopy longwave radiation, seen in Fig. 4b, for two reasons. Firstly, daytime correction depends on insolation, which changes throughout the snow cover season so that daytime correction varies to a 95 higher degree than nighttime correction. Secondly, nights are longer than days prior to the boreal spring equinox, which leads to nighttime increases outweighing daytime decreases. Consequently, correction results in increased average subcanopy longwave radiation even for equal magnitudes of 100 daytime overestimation and nighttime underestimation.
Comparison of simulated and measured longwave enhancement is shown in Fig. 4c and Fig. 4d for Alptal. As for sub-canopy longwave radiation, the diurnal cycle of simulated longwave enhancement is reduced by implementation of correction factors with increased enhancement at night and 5 decreased enhancement at daytime. Reduction of daytime longwave enhancement increases throughout the snowmelt season, which is due to increasing insolation and thus increasing reduction of sub-canopy longwave radiation during daytime. Longwave enhancement values vary between 1.1 and 1.4 in CTRL, which is predominately driven by diurnal cycles. The diurnal cycle of longwave enhancement is reduced by more than 50% in CORR, resulting in a diurnal range similar to observations and increased daily average longwave enhancement. Simulated longwave enhancement 15 displays little intraseasonal variability, with variations mostly due to the overestimated diurnal cycle. This indicates that intraseasonal variability in sub-canopy longwave radiation is largely due to variations in atmospheric longwave radiation. In contrast, measured longwave enhancement values range 20 from less than 1 to more than 1.6 and display little diurnal variability but high variability on synoptic timescales, which results in a different daily average of longwave enhancement compared to simulations. Moreover, lower average longwave enhancement for observations indicates more overcast con-25 ditions, which lead to smaller diurnal cycles in sub-canopy longwave radiation compared to simulations. Therefore, correction factors improve the realism of diurnal cycles of subcanopy longwave radiation and longwave enhancement, encouraging usage for evaluation of impact on snow cover.

30
The contrast in variability between simulated and observed longwave enhancement can be seen in Fig. 5. Observations show a large range of longwave enhancement values that are closely tied to effective emissivity of the sky, which represents clear-sky (low ε sky ) and overcast (high ε sky ) con-35 ditions. Observed longwave enhancement increases for decreasing ε sky as the contrast between vegetation temperatures, increasing due to higher insolation, and effective temperature of the sky increases. Spread in observed longwave enhancement is small throughout the range of ε sky , indi-40 cating little diurnal variability and the process of longwave enhancement depending on meteorological conditions. Simulations display a narrow range of ε sky , which causes the lack of intraseasonal variability seen in Fig. 4c. The spread in simulated longwave enhancement values is substantially 45 larger compared to observations for the respective range in ε sky representing overestimated diurnal cycles. Implementation of correction factors reduces the spread in longwave enhancement values and increases average longwave enhancement (see Fig. 4d), however, spread in longwave enhance-50 ment is still overestimated and average longwave enhancement is underestimated in CORR compared to observations for the respective range in ε sky .

55
Lack of variability in simulated ε sky , as seen for the grid cell of Alptal, across the Northern Hemisphere results in spatially similar correction factors that largely dependent on insolation. However, variability in both insolation and diurnal ranges of atmospheric longwave radiation indicate small 60 variations in meteorological forcing that are not represented by ε sky . Therefore, ε sky in simulations may indicate clearsky conditions even when insolation and atmospheric longwave radiation suggest more overcast conditions, resulting in overestimated correction factors and overcorrection of sub-65 canopy longwave radiation. This overcorrection results in larger nighttime than daytime values of sub-canopy longwave radiation in contrast to atmospheric longwave radiation and occurs mostly along continental coasts (Fig. 6). Consequently, a contour line is used in the following to denote an 70 overcorrection for 10% of days.
To demonstrate the impact of correction factors spatially, maps of longwave enhancement beneath evergreen needleleaf forests in CLM4.5 are shown in Fig. 7a and Fig. 7b. Averages over boreal winter and spring show an enhance-75 ment of longwave radiation beneath canopies by about 20% to 30% and display little differences across boreal forests, which is due to small spatial variability in both ε sky and vegetation density (Fig. 1). CORR displays increased average longwave enhancement north of 40 • N with an additional en-80 hancement of longwave radiation of up to 5% beneath dense boreal forests. Changes in longwave enhancement generally increase with latitude as daytime correction factors vary with insolation while nighttime correction factors are independent of latitude. A higher increase in longwave enhancement can 85 be seen for higher vegetation density within regions covered by boreal forests (Fig. 1b) due to weighting of contributions to subcanopy longwave radiation (Eq. (3)).

Snow cover and snowmelt
Changes in sub-canopy longwave radiation induced by the 90 correction increase the net energy flux to the surface, which can be seen for grid cell-averaged snow surface temperature ( Fig. 7c and Fig. 7d). Simulated average snow surface temperatures are determined by latitude, topography, and continentality, reaching values of less than -40 • C in the moun-95 tainous regions of northeastern Siberia (Fig. 1c), and range between -20 • C and -15 • C for boreal forests, the outlines of which can be seen in central Siberia and central North America. The impact of correction factors is limited to grid cells for which vegetation is dominated by evergreen needleleaf 100 trees and implementation results in an increase in average snow surface temperature of up to 2 • C. The lack of spatial variability is caused by little spatial variability in meteorological conditions, high vegetation density, and similarly high PFT coverage across boreal forests ( Fig. 1a and Fig. 1b). 105 Cold content, the energy required to raise snow temperatures to 0 • C, is used to quantify the impact of correction factors on the entire snow column. Average cold content simulated by CLM4.5 mostly reaches values of up to 4 MJm −2 and exceeds 5 MJm −2 only in glaciated grid cells (Fig. 7e).
In CTRL, simulated average cold content ranges between 1.5 MJm −2 and 3 MJm −2 across boreal forests, with lowest values in northeastern Europe and highest values in eastern Siberia, western Canada, and Quebec. Relative changes in cold content from CTRL to CORR display spatial differ-10 ences with cold content generally decreasing across boreal forests (Fig. 7f). Reductions in average cold content reach up to 30% in northeastern Europe and western North America and up to 20% in central North America. Across Siberian boreal forests, relative reductions decrease from west to east 15 from more than 20% to about 10%. Spatial differences in relative reductions correspond to spatial differences in average cold content, with higher relative reductions for smaller averages, representing a more even spatial pattern of absolute reductions in cold content as indicated by changes in snow 20 surface temperature (Fig. 7d).
Spatial patterns in snow-off date are similar to those in cold content with higher cold content corresponding to later meltout ( Fig. 7g and Fig. 7h). Changes in snow-off date from CTRL to CORR display stark spatial contrasts with meltout 25 happening up to 10 days earlier in central Europe and on the western coast of North America. Meltout is advanced by about 5 days for boreal forests in northeastern Europe and western Siberia and slightly less for boreal forests in central North America. In contrast, meltout is delayed in mountains 30 of southeastern Siberia (Fig. 1c), where meltout occurs late among boreal forests.
As offline simulations lack spatial variability in ε sky , latitude (through insolation) and duration of snow on the ground (through day length) control spatial differences in impact 35 of correction of sub-canopy longwave radiation on snow-off date. Changes in longwave enhancement due to correction of sub-canopy longwave radiation before and after the boreal spring equinox, approximated by averages over February/March and April/May, display opposite signs across the 40 Northern Hemisphere (Fig. 8), with shorter (longer) days than nights before (after) the equinox resulting in an increase (decrease) in daily average longwave enhancement. Generally, lower insolation at higher latitudes leads to a more positive impact of correction on daily average longwave enhance-45 ment, increasing (decreasing) positive (negative) changes in longwave enhancement with increasing latitude before (after) the boreal spring equinox. Across mid-latitudes, increase in daily average longwave enhancement over February and March is roughly similar to decrease in daily average long-50 wave enhancement over April and May, while increase over February and March outweighs decrease over April and May across high latitudes including most of the regions covered by boreal forests.
Reasons for spatial differences in changes of snow-off date 55 across Siberian boreal forests are explored in Fig. 9. Snowoff dates are similar spatially in CTRL, likely caused by higher elevations in southeastern Siberia compensating for less cold content, and meltout generally occurs past the boreal spring equinox in northwestern and southeastern Siberia. 60 However, higher insolation for southeastern Siberia results in higher reductions of daytime sub-canopy longwave radiation by correction factors and consequently smaller increases in daily average sub-canopy longwave radiation prior to the boreal spring equinox compared to northwestern Siberia. Al-65 though changes in sub-canopy longwave radiation are still positive in southeastern Siberia accumulated over the snow season, causing a decrease in cold content, reduction in daily average sub-canopy longwave radiation by correction factors past the boreal spring equinox cancels out the previous in-70 crease and consequently, snowmelt is slightly delayed. In contrast to southeastern Siberia, meltout is slightly accelerated in central North America although both latitude and meltout date are similar, as relative reductions in cold content are generally higher. However, differences in changes in 75 meltout date between central North America and southeastern Siberia are minor.

Snow-off date in comparison to observations
Simulated and observed snow-off dates are compared in Fig. 10 for grid cells with consistent snow cover throughout 80 preceding December and coverage by evergreen needleleaf trees of at least 50%. Simulations CTRL and CORR generally feature a narrower probability density function (PDF) of snow-off dates, indicating a shorter snowmelt season, and later meltout compared to observations across the en-85 tire Northern Hemisphere (Fig. 10a). While shapes of observed PDFs are well represented by simulations over Eurasia (Fig. 10b, d), observations show a clearer, shorter peak of meltout compared to simulations over mountainous western North America (Fig. 10c). Correction of sub-canopy long-90 wave radiation displays little impact when accumulated over the entire Northern Hemisphere, however, it systematically reduces the delay of simulated snow-off dates throughout the snowmelt season. PDFs of snow-off dates for regional subsets reflect spatial patterns seen in Fig. 7h, with minor 95 differences between CTRL and CORR over most of western North America (Fig. 10c) and eastern Siberia (Fig. 10d) but substantial acceleration of snow-off dates over western Siberia and eastern Europe (Fig. 10b) due to correction of sub-canopy longwave radiation. 100 The regionally limited impact of corrected sub-canopy longwave radiation is highlighted by filtering PDFs of snowoff date for grid cells with average differences in snow-off date between CORR and CTRL of at least 3 days (Fig. 10e,  f). Correction of sub-canopy longwave radiation improves 105 timing of meltout in filtered grid cells, especially over western Siberia and eastern Europe where the filtered PDF for Todt et al.: Simulated single-layer forest canopies delay Northern Hemisphere snowmelt 7 CORR, in contrast to CTRL, closely resembles observations. PDFs of snow-off dates derived from Blended-5 SWE display sensitivity to threshold choices, however, this uncertainty is generally smaller than differences between simulations and observations. 5 5 Discussion Todt et al. (2018) found roughly similar magnitudes for daytime overestimations and nighttime underestimations of subcanopy longwave radiation in CLM4.5; however, this study shows that different durations of day and night over the snow 10 cover season result in a net positive impact of correction on daily averages of sub-canopy longwave radiation. Correction factors change throughout the snowmelt season due to increasing insolation and length of day. Consequently, net impact on daily averages of sub-canopy longwave radiation 15 varies resulting in spatial differences in impact on cold content over the snow cover season and meltout date. Net increase in sub-canopy longwave radiation during the snow cover season is highest for regions of early snowmelt where snow is already comparatively warm, which results in accel-20 erated snowmelt. Lundquist et al. (2013) showed that forests enhance snowmelt compared to open area in regions where winters are warm and mid-winter melt events happen, during which longwave enhancement outweighs shading. Spatial differences in change of meltout date broadly agree with 25 this pattern as the highest acceleration of melt occurs for regions with warmer winters as indicated by snow surface temperatures (Fig. 7c), suggesting that mid-winter melt events could be underestimated by CLM4.5. Conversely, correction of sub-canopy longwave radiation results in slightly delayed 30 snowmelt in southeastern Siberia albeit average cold content over the entire snow cover season being reduced. This delay is due to meltout happening substantially later than the boreal spring equinox and high insolation during the snowmelt period, which result in reduction in daytime sub-canopy long-35 wave radiation outweighing increased sub-canopy longwave radiation during night. Consequently, overestimated diurnal cycles of sub-canopy longwave radiation in CLM4.5 lead to spatial differences in impact on snowmelt timing across boreal forests in offline simulations.

40
Previous comparison between offline simulations of CLM4 and observations have shown CLM4 failing to accurately simulate the timing of both snow ablation and snow accumulation across boreal forests, with snowmelt compressed into the period March to May (Thackeray et al., 2014(Thackeray et al., , 2015.

70
Changing seasonality in a warming climate may have implications for snowmelt and longwave enhancement. Future warming will lead to earlier snowmelt, when less energy from insolation is available for melt, which will likely result in lower melt rates (Musselman et al., 2017). A shortened 75 snow season indicates more asymmetrical lengths of day and night during snowmelt and consequently, overestimated diurnal cycles of sub-canopy longwave radiation in CLM4.5 could result in even higher underestimations in daily averages. Moreover, underestimated sub-canopy longwave radia-80 tion suggests that CLM4.5 underestimates melt rates in general. In turn, future projections are complex, as corrected and thus increased sub-canopy longwave radiation might cancel out reduced energy from insolation due to earlier snowmelt. Nonetheless, the contribution of longwave enhancement to 85 snowmelt is likely to increase in the future, further necessitating accurate simulation of sub-canopy longwave radiation.
Implementation of correction factors resulted in realistic average diurnal ranges of sub-canopy longwave radiation and longwave enhancement, but more substantial underestima-90 tion than overestimation of longwave enhancement seen in Fig. 5 suggests that the impact of shortcomings in CLM4.5 on snow cover and snowmelt might still be underestimated by this study. Gouttevin et al. (2015) and Todt et al. (2018) have shown the implementation of biomass heat storage to result 95 in a net positive impact on sub-canopy longwave radiation as well as a slight reduction of diurnal cycles. This suggests that heat storage by biomass could further reduce nighttime underestimation in CLM4.5 and improve the simulation of sub-canopy longwave radiation and longwave enhancement. 100

Conclusions
This study assessed the impact of deficiencies in simulated longwave enhancement by forest canopies on snow cover in CLM4.5. Sub-canopy longwave radiation simulated by 8 Todt et al.: Simulated single-layer forest canopies delay Northern Hemisphere snowmelt CLM4.5's single-layer vegetation was corrected based on the damping effect of multiple canopy layers. Correction factors were derived from forest stand-scale simulations and subsequently implemented for evergreen needleleaf trees in global land-only simulations of CLM4.5. Correction reduces overestimated diurnal cycles of sub-canopy longwave radiation by decreasing daytime overestimations and nighttime underestimations. This results in a net increase of sub-canopy longwave radiation over the entire snow cover season, due to longer nights than days. Consequently, correction results in 10 increasing average snow temperatures and earlier meltout, indicating that CLM4.5 underestimates snow temperatures and delays snowmelt due to overestimated diurnal cycles of sub-canopy longwave radiation. Comparison with observations confirmed a delay of meltout in land-only simula-15 tions of CLM4.5 across boreal forests, which is decreased by correction of sub-canopy longwave radiation. While landonly simulations exhibit a spatially uniform underestimation of snow temperatures by CLM4.5 across evergreen boreal forests, the impact of correction on meltout timing displays 20 spatial differences that depend on insolation and duration of snow on the ground. The effect of overestimated diurnal cycles on daily average sub-canopy longwave radiation changes throughout the snowmelt season as insolation and length of day increase. Consequently, CLM4.5 delays snowmelt more 25 in regions of warmer snow cover and earlier meltout. However, spatial variability in impact on snow cover is limited in land-only simulations of CLM4.5 due to a lack of variability in meteorological conditions.
Author contributions. MT, NR, and CF designed the experiments. MT and LW performed the simulations and MT analyzed them. MT prepared the manuscript with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict 50 of interest. Table 1. Forest stands used for calculation of correction factors based on simulations by Todt et al. (2018). Vegetation density is given as Plant Area Index (PAI), the one-sided area of plant components per unit ground surface area including stems, branches, and leaves or needles. Abisko and Cherskiy feature deciduous vegetation, so that trees were leafless throughout the simulation periods and PAI values do not consider leaves or needles.       Figure 5. Longwave enhancement measured (green) at the forest stand of Alptal, Switzerland and simulated in CTRL (black) and CORR (red) for boreal evergreen needleelaf trees in the corresponding gridcell of Alptal, Switzerland as a function of effective emissivity of the sky. Each data point represents an hourly average seen in Fig. 4c.  . Averages in CTRL (a, c, e, g) and differences between CORR and CTRL (b, d, f, h) for longwave enhancement beneath evergreen needleleaf trees (a, b), snow surface temperature (c, d), cold content (e, f), and snow-off date (g, h). Longwave enhancement is averaged over December to May while snow surface temperature and cold content are averaged over entire snow cover seasons. Differences CORR -CTRL are calculated as averages of differences between each individual snow cover season. For panels c-h, a mask is applied to filter out grid cells that are not perennially snow-covered. Black lines demarcate continental areas with less than 10% of overcorrected days. Green lines demarcate areas with coverage by evergreen needleelaf trees of at least 50%.    N). Panels e and f are as panels a and b, respectively, but only for grid cells with average changes in snow-off dates of at least 3 days (as seen in Fig. 7h).