Effect of Forest Canopy Structure on Wintertime Land Surface Albedo: Evaluating CLM5 Simulations With In‐Situ Measurements

Land Surface Albedo (LSA) of forested environments continues to be a source of uncertainty in land surface modeling, especially across seasonally snow covered domains. Assessment and improvement of global scale model performance has been hampered by the contrasting spatial scales of model resolution and in‐situ LSA measurements. In this study, point‐scale simulations of the Community Land Model 5.0 (CLM5) were evaluated across a large range of forest structures and solar angles at two climatically different locations. LSA measurements, using an uncrewed aerial vehicle with up and down‐looking shortwave radiation sensors, showed canopy structural shading of the snow surface exerted a primary control on LSA. Diurnal patterns of measured LSA revealed strong effects of both azimuth and zenith angles, neither of which were adequately represented in simulations. In sparse forest environments, LSA were overestimated by up to 66%. Further analysis revealed a lack of correlation between Plant Area Index (PAI), the primary canopy descriptor in CLM5, and measured LSA. Instead, measured LSA showed considerable correlation with the fraction of snow visible in the sensor's field of view, a correlation which increased further when only considering the sunlit fraction of visible snow. The use of effective PAI values as a simple first‐order correction for the discrepancy between measured and simulated LSA in sparse forest environments substantially improved model results (64%–76% RMSE reduction). However, the large biases suggest the need for a more generic solution, for example, by introducing a canopy metric that represents canopy gap fraction rather than assuming a spatially homogeneous canopy.

Wintertime LSA across forested regions varies in response to intercepted snow (Bartlett & Verseghy, 2015;Stähli et al., 2009), tree species (Kuusinen et al., 2012), and canopy structure itself (Bright et al., 2018). Webster and Jonas (2018) further highlighted the importance of canopy structural shading for LSA during clear sky conditions. Additionally, as outgoing shortwave radiation is governed by reflections and scattering, sub-canopy radiation processes are vital for overall LSA. Together these factors control LSA, but how well they are represented by LSMs remains uncertain. In order to identify such missing or poorly represented processes in LSMs, there is a strong need to bridge the gap between the contrasting scales of model resolution and in-situ measurements (Williams et al., 2009). This is vital to identify sub-grid processes that must be explicitly accounted for, and those that can be neglected or represented in a simplified manner (Bierkens et al., 2015). For example, a recent modeling study has demonstrated that radiative heterogeneities within a forest stand often remain unaccounted for in coarse-resolution LSMs due to simplified canopy structure representation (Rosati et al., 2020). Mao et al. (2016) also used such an approach to investigate the ability of a LSM to replicate carbon allocation, which further helped to identify unconsidered sub-grid processes of interest for the LSM community.
Global climate models are usually run at coarse horizontal resolutions (∼50-200 km grid cells), leaving satellite data as the only way to assess model performance with regards to simulated LSA over the spatial extents of the simulations (e.g., Loranty et al., 2014;Malmros et al., 2018;Thackerayet al., 2015;Wang et al., 2004). At the resolution of commonly used satellite data (e.g., MODIS MOD10A1 500 m grid), forest heterogeneity in terms of canopy gaps and edges is not resolved. While higher resolution satellite albedo retrievals are becoming more widely available (e.g., 30 m LANDSAT or Worldview4 by Digital Globe), the spatial and temporal resolution is still not high enough for such data to be used as a benchmark to investigate model performance regarding forest-induced variability of wintertime LSA. At the same time, in-situ measurements of LSA are usually at the point-scale and struggle to represent areas greater than the sensor footprint, which would be required to represent the spatially integrated perspective of a coarse grid cell. A possible solution towards mitigating this spatial discontinuity data problem of snow covered forests are airborne platforms (e.g., Lundquist et al., 2018), which when compared to stationary tower measurements, allow for variable measurement heights and the flexibility to obtain data at many locations above a forest canopy. The resulting increased measurement resolution enables impacts of individual spatially and temporally varying factors to be resolved (e.g., the effect of canopy structure on LSA). In order to then close the gap to coarse grid cell model resolution, it is necessary to combine such spatially distributed airborne measurements with point-scale simulation, which allows evaluation of algorithms used within global climate modeling frameworks at the process level, further facilitating LSM development.
This study uses uncrewed aerial vehicle (UAV) measurements of LSA coupled to coincident down-looking hemispherical images over snow covered forests to evaluate the simulation of (a) diurnal, (b) seasonal, and (c) spatial variability of shortwave radiative fluxes by the Community Land Model 5.0 (CLM5). UAV measurements of broadband LSA were made over a large range of canopy structures in two climatically different locations (Switzerland and Finland), across wide ranges of solar zenith and azimuth angles. UAV LSA measurements are first used to investigate the effect of shading on LSA and then used to evaluate CLM5 simulations of LSA, both spatially and temporally as a function of solar angle and canopy structure. These results may motivate the land modeling community to incorporate explicit canopy structure going forward.

Site Description
Field experiments for this study were conducted in both subalpine (Swiss Alps) and boreal (Finnish Lapland) forest environments, covering two different climatic conditions. Locations of field sites ( Figure 1c) were selected to cover a range of tree species with varying canopy structures and a gradient of canopy density at each site. A spruce site in Davos Laret, Switzerland (Figure 1d), and a pine site in Sodankylä (Figure 1e) served as intensive observation sites, where the bulk of the field measurements and corresponding model experiments was conducted. In order to widen the range of observed tree species, additional measurements were performed at three further subalpine sites (larch, beech, and pine site) as well as one further boreal site (birch site). The larch site was located in Flin, Switzerland (blue star in Figure 1c), the beech site in Maienfeld, Switzerland (green star in Figure 1c) and the pine site on the Ofenpass, Switzerland (pink star in Figure 1c), while the birch site was located in proximity to the pine site in Sodankylä. Since these additional sites were each only visited on one clear sky day during the 2018/19 winter season, these data were used to test the impact of sun-view fraction on LSA but not for broader model evaluation purposes.

Data Collection
At each site, we performed in-situ measurements of LSA, which were used for the two main objectives of this study: to investigate the effect of shading on LSA, and, to assess CLM5 model performance at the process level. Albedo was measured using a DJI S1000 octocopter equipped with up-and down-looking Kipp and Zonen CMP3 shortwave radiation sensors. CMP3 pyranometers measure shortwave broadband albedo in the spectral range of 300-2,800 nm. A full description of the UAV platform and sensors can be found in Webster and Jonas (2018). Programmed flight plans for each site allowed exact repetition of measurement surveys at varying solar angles throughout a day, and across a winter season. Flying height of the UAV was programmed to be between 15 and 30 m above the top of canopy. Each flight plan consisted of between 6 MALLE ET AL.  and 12 waypoints (WPs); the UAV stopped and hovered for 15 s at each of the WPs. As an example, Figure 1b shows a transect through a 3D lidar point cloud along the flight plan for the spruce site, where WP locations are indicated by green crosses. Table 1 summarizes field site characteristics and amount of data collected per site, including number of flights, WPs per flight, and captured zenith/azimuth angle ranges.
Radiation data for each WP was only used if at least 6 s of data were within a 3 m buffer in x, y, and z direction around a defined WP coordinate, and a filter was applied to omit any data exceeding a sensor tilt threshold of 5° or more, as suggested by Bogren et al. (2016). A down-looking Sony Alpha NEX6 16.1 MP camera with a Yasuhara Madoka f/4 7.3 mm 180° fish-eye lens was attached to the UAV and set to trigger every second during each flight. Each down-looking hemispherical photograph reflected the 180° view-field of the down-looking pyranometer. These hemispherical photographs were each binarized, discretized into 10 concentric analysis rings and analyzed following the methodology in Malle et al. (2019): First, the land cover descriptors tree-view fraction (VF Tree ) and snow-view fraction (VF Snow ) were obtained from the ratio between black and white pixels within each analysis ring, which was further weighted by the sine of the elevation angle. VF Snow was then split into sunlit snow-view fraction (VF Sun ) and shaded snow-view fraction (VF Shade ) based on a secondary binarization and analysis. Figure 1a shows an example down-looking hemispherical photograph as well as the resulting categorized image, demonstrating the spatial extents of the corresponding view-fractions.
For each WP location at the intensive observation sites, canopy coverage (CC), defined as the ratio of area covered by the vertical projection of the canopy relative to ground area (e.g., Mazzotti et al., 2019), was derived from a canopy height model (CHM). A CHM was generated based on terrestrial lidar data for the pine site in Sodankylä and based on airborne lidar data for the spruce site in Laret, following Mazzotti et al. (2020) and Khosravipour et al. (2014). The CHMs were binarized based on a 2 m threshold and 10 m circular fetches were used to compute local CC around each point.
In order to obtain Plant Area Index (PAI, sometimes interchangeably referred to as Vegetation Area Index or Leaf Area Index (LAI) (Jonckheere et al., 2004)) at each WP location, in-situ up-looking hemispherical images were taken at all 9 WPs in Laret and at 6 out of the 10 in Sodankylä, while at the remaining 4 boreal sites synthetic hemispherical images were generated from terrestrial lidar data following Webster et al. (2020). Hemispherical images were first binarized into tree and sky, using a threshold algorithm provided by the model HPEval . PAI was then calculated using the image analysis software "Hemisfer" (Schleppi et al., 2007)  ote. Same colors are used as in Figure 1c, note that the pine and birch site in Finland have the same color due their close proximity.

LSA in CLM5
CLM5, the land component of the Community Earth System Model, is a state of the art, process-based LSM that simulates carbon, nitrogen, and energy exchange between the atmosphere and the terrestrial Earth (Lawrence et al., 2018(Lawrence et al., , 2019. In CLM5, vegetation is conceptually represented by a homogeneous layer of leaves covering the land surface (also referred to as the big leaf approach). Shortwave radiative transfer through this single layer is calculated via the two-stream approximation (Dickinson, 1983), consisting of a coupled pair of inhomogeneous differential equations for the upward and downward radiative flux. The two-stream approximation is solved separately for the visible and near infrared spectrum, as well as for incident direct and diffuse radiation, ultimately determining how forest canopies control reflection, absorption, and transmission of solar energy. CLM5 uses five spectral bands for the shortwave radiation calculations, one visible (300-700 nm) and four near infrared (700-1,000 nm, 1,000-1,200 nm, 1,200-1,500 nm, 1,500-5,000 nm). Note that corresponding CMP3 pyranometer measurements cover all five bands, four of which entirely and the last one partly. Leaves are assumed to be randomly distributed in space and treated as flat surfaces with identical optical properties on top and bottom. The canopy is horizontally uniform and all diffuse radiative fluxes are viewed as isotropic in both upward and downward direction. The full solution to the two-stream approximation and its limitations are presented in detail in Dickinson (1983) and Sellers (1985).
Overall LSA in CLM5 is further modulated by forest floor albedo which depends on the interplay of snow albedo, soil albedo, and snow cover fraction. Forest floor snow albedo in CLM5 is computed via the SNow ICe and Aerosol Radiation model (SNICAR), which yields overall snow albedo (Flanner & Zender, 2005). Soil albedo, also needed as a boundary condition for the SNICAR snow radiative transfer calculation, is dependent on soil color class which is prescribed for each grid cell and soil water saturation level. Fractional snow covered area calculations follow Swenson and Lawrence (2012), with different parameterizations used for snow accumulation and snow depletion periods. Compared to previous model versions, CLM5 includes an updated snow interception scheme, which allows temperature and wind induced snow unloading, further providing the fraction of the canopy covered in snow used for LSA calculations. However, this study only focuses on days without snow in the forest canopy. A technical description of CLM5 including new model features in the latest version can be found in Lawrence et al. (2019).

CLM5 Model Setup
Spatial heterogeneity of the land surface is accounted for in CLM5 by a sub-grid hierarchy system. Each grid cell is split into different land units (vegetation, glacier, lake, urban, crop), and vegetated land-units further consist of up to 15 Plant Functional Types (PFTs) plus bare ground. For the purpose of this study, we apply the global-scale modeling algorithms of CLM5 to the point-scale for which CLM5 features a dedicated point mode (PTCLM). Each PTCLM simulation corresponds to a WP on a flight path. These point-mode CLM5 simulations of snowpack mass and energy fluxes were run at 21 WPs, nine at the subalpine spruce site (Laret, CH), and 12 at the boreal pine site (Sodankylä, FIN), using site-specific meteorological forcing and waypoint-specific forest canopy data. While the two sites contained different tree species (spruce, pine), they were both classified as the PFT Needle-leaf Evergreen Trees (NET). Hence, each CLM5 run was set up to simulate a grid cell being fully (100%) covered by NET of variable densities. Meteorological driving data (incident short and long-wave radiation, air temperature, relative humidity, wind speed, pressure, and precipitation) at 1 h temporal resolution was derived from automatic weather stations (AWS) in unforested, open locations within 1 km of each forest site. CLM5 simulations ran between January 2017 and August 2018 at Laret and between January 2018 and August 2019 at Sodankylä to coincide with evaluation data while also allowing one year of spin-up. A spin-up was necessary to ensure soil moisture and soil temperature were in equilibrium and not affecting temporal dynamics and physical properties of the simulated snow cover evolution.
In CLM5, vegetation structure for each PFT is described by monthly varying LAI and Stem Area Index (SAI), as well as heights of canopy top and canopy bottom. For global simulations, percent PFT are derived from MODIS satellite data (Lawrence and Chase, 2007), as are monthly LAI and SAI values. The combination of LAI and SAI results in the PAI. Figure 2 shows the global data set of wintertime NET PAI (MODIS-derived, 0.25°) commonly used in large-scale CLM5 simulations. For our CLM5 simulations in point-mode, however, these canopy descriptors were calculated from a combination of above-mentioned terrestrial lidar data (boreal sites; canopy heights), airborne lidar data (subalpine sites; canopy heights), and hemispherical photography (all sites; PAIs). PAI was split into LAI and SAI by applying the proportions given in the 0.25° MODIS-derived LAI and SAI standard CLM5 input data set (0.28 for Sodankylä, 0.33 for Laret). In addition, synthetic experiments were set up for both Sodankylä and Laret: PAI was varied from 0 to 5 at 0.1 increments, while the meteorological and remaining point-specific canopy data remained the same. This allowed us to specifically investigate the sensitivity of simulated LSA to PAI.

Model Adaptions for This Study
CLM5 simulations at unforested open sites adjacent to forest measurement locations allowed parameter optimization, for example, to account for gauge under-catch of measured precipitation. Snow under-catch factors for solid precipitation were determined as 1.3 and 1.5 for Laret and Sodankylä, respectively, both sites featuring unshielded gauges. Initial simulations at all forested sites showed that snow melted too quickly when compared to observations. The problem was attributed to an exaggerated sensible heat exchange between the canopy air space and the snow surface. This was mitigated by reducing the turbulent transfer coefficient under canopy (csoil cn ). Additionally, the fraction of the canopy covered by snow (fcansno) used for LSA calculations was adapted to be less sensitive to the amount of intercepted snow. In the original code, canopy albedo was unrealistically responsive to even tiny amounts of snow in the canopy. Fractional snow covered area in CLM5 is calculated following Swenson and Lawrence (2012), which uses a snow covered area shape function. This function depends on the standard deviation of topography within a grid cell. As the impact of topography was assumed to be small for all sites, a comparatively small standard deviation value (25) was chosen, which delayed the onset of fractional snow cover and showed a rapid decline in snow cover extent toward the end of the melt season. Collectively, these model adaptions and parameter choices led to more realistic snowpack and LSA developments at all waypoint locations. Table 2

Empirical Findings: LSA and Sunlit Snow-View Fraction
LSA measurements from a UAV allowed a seasonal and diurnal analysis of LSA in relation to solar position and canopy structure at needleleaf and broadleaf forests. Figure 3 shows all LSA measurements taken during clear sky, full snow cover, and no interception conditions. Each point represents data from UAV measurements at one unique WP, contrasting LSA with coincidently measured VF Sun . The analysis is integrated across the spatial (different WPs within each site) and temporal (repeated measurements throughout the day) dimension. Strong positive relationships between LSA and VF Sun were found at all sites, all of which were statistically significant. Relationships were stronger for needleleaf sites (R 2 = 0.79-0.94) compared to deciduous sites (R 2 = 0.59-0.88). A linear regression model was fitted to data points from each field site, indicating y-intercepts between 0.07 and 0.18 as well as gradients between 0.38 and 0.72 for sites with deciduous environments and gradients between 0.38 and 0.61 for sites in needleleaf forest environments. The strong positive relationships suggest that canopy structural shading of the snow surface is exerting a primary control on LSA on clear sky days. Note again that each site encompassed several WPs, hence, in this evaluation, variation in VF Sun is a function of the zenith and azimuth angle of the solar position as well as the local VF Snow , which is explaining some of the spread.

Effect of Solar Zenith Angle on Observed and Simulated LSA
To disentangle the effects of zenith angle and VF Snow on LSA, we show data from a series of flights at the spruce site in Davos Laret over the course of one day (Figure 4a). At waypoints where the forest canopy is sparse (VF Snow > 0.3), LSA decreases by up to 11% as the zenith angle increases. This occurs due to the fact that an increasing solar zenith angle corresponds to an increased canopy shaded area of the snow surface, as long as the canopy is sparse enough to permit direct insolation of the forest floor at all.
Most of the decline in LSA happens at low to medium zenith angles (35°-55°), while there is only a slight further decrease toward very large zenith angles (>75°) following sunrise and preceding sunset. In contrast, at waypoints with denser canopy (VF Snow < 0.25), LSA was seen to increase by a few percent as zenith angle increases. In denser forest the proportion of sunlit snow is low at all times, and sun flecks become rare even in the middle of the day. As the solar zenith angle increases, the proportion of diffuse radiation in incoming solar radiation also increases. Hence, below a certain VF Snow , this increase in diffuse radiation and the corresponding increase in multiple scattering leads to increased LSA with increased solar zenith angle.
Corresponding CLM5 simulations (Figure 4b) show only a minimal response to diurnal changes of the zenith angle. Contrary to our measurements, simulated LSA increases slightly with zenith angle even for very sparse points (VF Snow > 0.4), but only by less than 5%. These small differences are a function of the interplay between forest floor snow albedo (dependent on grain size and the angle of incidence of the solar MALLE ET AL.  rays) and canopy radiative transfer calculations in CLM5. At WPs with denser canopy, CLM5 simulations replicate magnitudes of measured LSA fairly well. Most notably, however, CLM5 simulations feature a large positive bias of up to 40% in LSA estimates for sparse canopy. This magnitude bias is investigated further in the following sections.

Effect of Azimuth Angle
In addition to zenith angle, the azimuth can also be important for LSA, as it controls the orientation of shadows on the snow surface. As an illustrative example, we present data from a boreal pine site in Sodankylä ( Figure 5). The location of trees and canopy gaps in relation to the solar position resulted in a large proportion of the snow surface being sunlit in the morning, while the site had a smaller sunlit-view fraction at the same zenith angle in the afternoon. This resulted in a non-symmetrical diurnal signal of MALLE ET AL.   CLM5 simulations show only a minimal variation of the simulated LSA with azimuth angle (Figure 5c) and the measured decrease in LSA at WPs with sparse canopy during the second half of the solar cycle is not represented in the CLM5 output. CLM5 does not consider the azimuth angle, so any corresponding variation of LSA is of indirect nature only. Figure 4d further shows a strong response of simulated LSA to VF Sun (R 2 = 0.63), however, only between WPs but not within each WP. Consistent with our so-far results, for this site and day, simulated LSAs for sparse canopy was consistently overestimated, with a high bias of up to 20% (i.e., simulated LSA of 0.6 compared to a measured albedo of 0.4).

Seasonal LSA Evolution
Next, we compared simulated LSA time series with measurements taken in clear sky conditions between February and May 2018 at the spruce site in Davos Laret. CLM5 consistently overestimated LSA in very sparse environments (CC < 20%), while LSA bias was mixed across dense canopies (CC > 50%). Specifically, compared to the individual measurements throughout the season at both dense and sparse WPs, CLM5 simulations overestimated LSA by 5%-66% in 55% of the cases and underestimated LSA by 5%-12% in 24% of the cases, while the remaining 21% were within 5% of the measurements. Underestimation of LSA predominantly occurred during low zenith angles around solar noon. Forest floor snow albedo also modulates overall LSA, and it is at the very sparse areas where it has the strongest effects. In most cases, however, the effect of the forest canopy prevails. A clear response to interception events can be noted as spikes in simulated LSA coincide with snowfall events and corresponding snow in the canopy. Furthermore, during springtime starting in mid-April, simulated LSA strongly declines responding to melt-out of snow on the ground, especially in sparse canopies. Simulated LSA in Figure 6 also demonstrates that, as CC decreases from dense to sparse, the diurnal shape of simulated albedo changes from a parabola shape with the vertex at the lowest point to an upside-down parabola, with this transition occurring at approximately 30% CC.

Effect of Canopy and Land Surface Descriptors
Measured LSA is strongly correlated to canopy shading of the snow surface, which is itself a function of the position of the Sun relative to the structure of the canopy surrounding each WP. Here, we focus on the effect of canopy and land surface descriptors, rather than a moving solar angle. Therefore, in an attempt to isolate the relationship between LSA and various canopy and land surface descriptors from the effect of the moving Sun, we limited the following analysis to data from a confined range of zenith angles (57°-67°), MALLE ET AL.  which was chosen to include the maximum amount of measurements at both the pine site in Sodankylä and the spruce site in Davos Laret. Figure 7 shows the correlation between either measured (row one in Figure 7) or simulated LSA (row two in Figure 7) and four different metrics that relate to canopy structure and radiation transfer (PAI, 1-CC, VF Snow , VF Sun ). Although CLM5 is built around a radiation transfer model that emphasizes an exponential relationship between LSA and PAI, measurements show comparatively little correlation between those variables. On the contrary, measurements suggest a strong correspondence between LSA and VF Sun , a relationship that is currently unaccounted for in CLM5. Correlations for 1-CC and VF Snow constitute a transition between the above extreme cases. The four descriptors represent a sequence of canopy metrics from spatially generalized (PAI), to local environment (CC), to specifically account for the spatial arrangement of surrounding canopy elements independent of the position of the Sun (VF Snow ), or in direct relation to the position of the Sun (VF Sun ). Along this sequence, correlation to LSA increases for observational data (Figures 7a-7d), but decreases for CLM5 simulations (Figures 7e-7h). The increasing correlation for observed data (Figures 7a-7d) is more pronounced for the pine than for the spruce site. For both sites, however, VF Sun had the best correlation with measured LSA (0.89 and 0.96 for the pine and the spruce site, respectively) and the weakest with simulated LSA (0.48 and 0.6 for the pine and the spruce site, respectively).

LSA Sensitivity to PAI
Data from WPs over sparse canopy environments have demonstrated CLM5 considerably overestimates LSA in these environments. In an attempt to explain the high bias, we performed a sensitivity analysis based on the CLM5 LSA simulations with synthetically varying PAI, whose values were selected to lie between 0 and 5 at 0.1 increments. For these simulations, input data and model settings remained unchanged to our previous simulations. Focusing on one clear sky day in Sodankylä (April 26, 2019), Figure 8 highlights that most of the reduction in simulated LSA (approximately 60%) with an increasing PAI happens between PAI  0 and PAI 1 while simulated LSA only decays marginally with PAI when PAI exceeds 2. These synthetic modeling results show that in sparse forest environments, even small changes in PAI value can have large effects on simulated LSA, with the caviat that once a LSA saturation is reached (∼PAI 2), the input PAI only has marginal effects.

Effective PAI Values
The high LSA sensitivity to PAI in sparse forest environments (Figure 8) and the systematic high biases in CLM5 compared to measurements (Figures 4-6) suggest that the canonical exponential relationship between LSA and PAI (Figure 7e) may exaggerate LSA increases with decreasing canopy density (Figure 8). This begs the question, what effective PAI values would be needed in the radiation transfer equations to replicate measured LSA? To this end, we derived an optimal PAI for each WP location individually ("best fit PAI") by minimizing the error of each WP simulation. Specifically, we minimized the sum of the absolute differences between the WP LSA simulation and all UAV-based LSA measurements performed at this WP location. Best fit PAI values were mostly higher than those estimated from point specific hemispherical images, particularly for the spruce site. We further fitted a logarithmic function to those best fit PAI values, whereby each data point was weighted dependent on the resulting improvement of simulated LSA when using the best fit PAI values. This approach provided a parameterization of effective PAI as a function of HP-derived PAI, the results of which are summarized in Figure 9. Largest improvements of simulated LSA by using the best fit PAI across all measurement points were found for PAI values <1, with larger effects for the spruce compared to the pine site ( Figure 9). The shown best fit PAI equations were applied to compute effective PAIs at each WP location, which were subsequently used to re-run CLM5 simulations. CLM5 LSA simulations using these effective PAIs showed a substantial improvement for both the spruce (Figures 10a and 10b, 76% RMSE reduction) and the pine site (Figures 10c and 10d, 64% RMSE reduction).

Discussion
UAV-based measurements of LSA during clear sky conditions showed a strong correlation with VF Sun across a large range of canopy structures in two climatically different locations. This highlights the strong influence canopy shading of the snow surface has on LSA during clear sky conditions, which was demonstrated to be valid across a much wider range of canopy structures and PFTs than initially shown by Webster and Jonas (2018). In addition, results in this study strengthen the spatial transferability of shading as a key control on outgoing shortwave radiation below the forest canopy, shown at a single site by Malle et al. (2019). Overall, sub-and above-canopy measurements from these three studies demonstrate the importance of canopy shading for radiative processes in snow-dominated environments, yet modeling results presented in this study further demonstrate that the resulting spatial and temporal LSA heterogeneities are unaccounted for in LSMs due to the simplified representation of canopy structure.
Assessment of LSM performance is complicated due to the often contrasting spatial scales of model resolution and in-situ measurements as well as the concurrence of parametric and structural uncertainties in modeling frameworks (Keenan et al., 2011). In this study, the combination of spatiotemporal LSA measurements with corresponding point-scale CLM5 simulations allowed identification of deficiencies in model process representation, for example, canopy structural shading of the snow surface. Canopy structural shading of the snow surface is controlled by the geometrical arrangement of trees in relation to the solar zenith and azimuth angles. UAV measurements showed both zenith and azimuth angles control diurnal variability of LSA. However, corresponding CLM5 simulations were not able to replicate the observed diurnal or spatial variability in LSA. The two-stream approximation used in CLM5 (Sellers, 1985) only takes zenith angle into account, which changes simulated LSA by marginal amounts in early mornings and late afternoons. Note that we did not investigate the effect of topography on LSA as all field sites were flat. Largest discrepancies between simulations and measurements of LSA were identified in sparse forest environments, consistent with Yuan et al. (2017) who found the two-stream approximation performed best in dense forests with limited spatial heterogeneity. Gaps and edges, predominant in sparse forest environments, increase the potential for direct insolation of the snow surface and result in LSA which are highly variable in space and time. It is in these environments where CLM5 simulations showed the most substantial deficiencies and failed to replicate the measured diurnal and spatial variability in LSA.
CLM5 uses the big leaf approach which simplifies the canopy structure into a homogeneous layer of absorbent matter. Moreover, absorption is parameterized as a function of PAI, which is a spatially integrated metric that neglects any structural heterogeneity within the canopy of interest. In consequence, CLM5 is unable to reflect the complex interplay between solar position and gaps in the surrounding canopy as evidenced by subcanopy radiation measurements Mazzotti et al., 2019), which in turn influence LSA. The particularly striking biases for sparse canopies suggest that CLM5 not only over-represents the penetration of direct light through sparse canopy, but further does not account for the spatial and temporal heterogeneity in canopy snow surface shading. Simulations further showed how in sparse forest environments even small changes in PAI can have large effects on simulated LSA. This is in line with findings by Manninen and Jääskeläinen (2018), who albeit only looking at diffuse radiation, found that CLM4 LSA decays more steeply with increasing PAI compared to two other albedo schemes, the Canadian Land Surface Scheme (Verseghy et al., 1993) and the forest reflectance model PARAS (Rautiainen & Stenberg, 2005). Wang et al. (2016) also found that most of the LSA reduction in CLM4 happened as PAI is increased from 0 to 1, and that an LSA saturation occurred for PAI values exceeding 3. The global data set of wintertime NET PAI (MODIS-derived, 0.25°) commonly used in CLM5 simulations (see Figure 2), shows that very sparse forest canopies (PAI < 1) are found in approximately 8% of the NET covered pixels. These pixels are mostly located around the circumpolar taiga to tundra transition, an area particularly sensitive to climate change (Serreze et al., 2000;Swann et al., 2010). Dense forest canopies (PAI > 2), on the other hand, account for approximately 73% of the NET area. Given that our model sensitivity experiment showed that CLM5-simulated LSA was insensitive to PAI values greater than 2, but we did not see this saturation in our measurements, it is likely that CLM5 is underperforming in LSA estimates across almost three quarters of the Northern Hemisphere NET area. We were able to mitigate part of the problem by using effective PAI values. Effective PAIs could be a first step to remove the overall bias in LSA for sparse canopies without the need of structural MALLE ET AL. changes to the radiation transfer scheme built into CLM5. However, such an approach would not help to fix the limited diurnal variation of simulated LSA, which would necessitate a more detailed representation of radiation transfer through canopy gaps. It is worth noting here that as CLM forms the foundation of several CMIP5 models, for example, the Norwegian Earth System model NorESM1 (Thackeray et al., 2018), the identified weaknesses in CLM5 are hence likely to extend into other modeling environments.
Investigation into the relationship between model error and land surface descriptors showed that VF Sun has the best correlation with measured LSA. While VF Sun can be straightforward to measure in this instance, incorporating this variable into LSMs would be much more challenging. In recent years, airborne lidar data has become a valuable source of detailed canopy structure information, which has motivated the development of models that can characterize radiation transfer through forest canopy at much greater detail Musselman et al., 2013;Zellweger et al., 2019). Recent efforts have resulted in methods that are efficient enough to resolve insolation patterns down to meter and minute resolution even over fairly large areas (∼1 km 2 , Webster et al., 2020). The benefit of such a detailed representation of radiation transfer processes within a snow model was first demonstrated by Musselman et al. (2012). Currently, there are at least two forest snow models capable of fully distributed simulations at very high spatial resolution, that include similarly detailed radiation transfer models and reproduce observed snow accumulation and melt dynamics more accurately than conventional models (Broxton et al., 2015;Mazzotti et al., 2020). While these approaches may be too computationally expensive for inclusion in LSMs intended for hemisphere-scale applications, they demonstrate the increasing potential for spatially distributed radiation transfer schemes capable of explicitly resolving radiation transfer through canopy gaps. In a coarse scale LSM like CLM5, structural heterogeneity needs to be better described in order to explicitly account for diurnal variability of LSA. We suggest including a metric which accounts for canopy structure in 2D, combined with a radiative transfer approach that allows for temporal interaction between solar position and such 2D canopy structures. Such an approach would allow the temporal variability of canopy shading of the snow surface to be more explicitly represented and should be a focus in the development of the next generation of LSMs.

Conclusions
For the first time, this study used temporally and spatially distributed in-situ measurements of LSA to evaluate the performance of the radiative transfer scheme implemented in a LSM over seasonally snow covered forest stands. An extensive data set of airborne LSA measurements using an UAV co-registered to canopy structure and snow surface shading information, captured a large range of canopy structures and solar angles in subalpine (Switzerland) and boreal (Finland) locations. These measurements were used to evaluate corresponding point-scale simulations with the CLM5. Measurements revealed a strong correlation between LSA and sunlit-snow across a range of tree species, solar angles, and canopy structures, further emphasizing the dominant control canopy structural shading has over LSA. CLM5 simulations did not adequately capture the measured spatial and temporal variability in LSA. This is partly due to the generalized canopy structure in the radiative transfer scheme, which over-simplifies the interaction between solar position and canopy heterogeneity. Additional errors then resulted in the amount of solar radiation incident at the snow surface, particularly in sparse forest environments. Analysis of diurnal patterns of measured LSA further revealed strong effects of both azimuth and zenith angles on LSA, which were not replicated by CLM5 simulations. Decay in simulated LSA with increasing PAI (the main canopy descriptor in CLM5) between 0 and 2 was too steep. As simulations further showed strong overestimations of LSA in sparse forest environments, we propose the use of effective PAI values as a simple first-order correction for this discrepancy. However, more complex canopy descriptors are needed which account for gaps and capture the structural heterogeneity of forest stands. Such model developments would help decrease uncertainty in LSA simulations across seasonally snow covered forest environments, with profound implications for snow albedo feedbacks.

Data Availability Statement
LSA data presented in this study are available from the WSL data repository Envidat at their website (www.envidat.ch/dataset/lsa_forest_snow).