Observations and snow model simulations of winter energy balance terms within and between different coniferous forests in southern boreal Finland

The boreal environment is governed by forests and by winter. The coniferous forest canopy has been shown to affect spatial patterns of snow depth and snow quality (e.g. Lundberg & Koivusalo 2003; Stähli & Gustafsson 2006; Veatch et al. 2009). Local differences in the canopy interception of snow lead to significant small-scale variation in the mass balance of boreal forest snow cover (Hardy & Albert 1995; Hedström & Pomeroy 1998; Link & Marks 1999).

The canopy is the governing factor in the snow surface energy exchange in boreal coniferous forests (Davis et al. 1997). Energy transfer to and from the snow include radiation, sensible and latent heat flux, soil heat flux, and heat flux by mass advection. Large spatial variability is seen both in incoming shortwave and incoming longwave radiation fluxes below coniferous forest canopies. The temporal persistence of the spatial distribution of the fluxes depends on the canopy structure and is different for longwave and shortwave fluxes (Essery et al. 2008).

In most of the cases, radiation is the primary source of energy for the snow cover below the forest canopy (Pomeroy & Dion 1996; Link & Marks 1999; Hardy et al. 2004). Observations at alpine sites show that radiation dominates the snow cover energy balance throughout the snow season, and radiation accounts for most of the energy available for snow melt (Marks & Dozier 1992; Davis et al. 1997). Incoming shortwave radiation is effectively attenuated by needles, branches and trunks; intercepted snow also affects the winter-time radiation levels in coniferous forests (Stähli et al. 2009). Spatial variability of the radiation components is important especially during snow melt (Pomeroy et al. 2008; Musselman et al. 2012a, b), most of the variability is due to variability in the direct beam component (Musselman et al. 2013). Litter on the snow surface or dark ground visible through the thinner snow cover often found near tree trunks lowers the albedo of the snow and thus may increase the amount of radiation absorbed (Melloh et al. 2001; Pomeroy et al. 2001).

There is also large spatial variability in the sub-canopy thermal regime and, consequently, in the longwave radiation reaching the snow surface (Rowlands et al. 2002). Shortwave radiation absorbed by the canopy is thermally radiated to the snow surface as longwave radiation. The highest shortwave radiation fluxes are observed in forest openings, but the highest longwave radiation values are observed close to warm, sunlit tree trunks (Essery et al. 2008).

The canopy also affects the snow evaporation in the forest by emitting longwave radiation during the night, keeping the snow surface relatively warm (Bernier & Swanson 1993). The canopy shelters the snow cover from turbulent heat exchanges of latent and sensible heat, which depend on wind speed and below-canopy temperature and humidity gradients (Harding & Pomeroy 1996; Storck et al. 1999). Canopy density (defined as fraction of overlying hemisphere occupied by canopy), and in heterogeneous forests also the size of the forest openings, affect the level of forest snow turbulent heat fluxes. Turbulent heat fluxes are relatively small, but they are important for the snow energy balance during mid-winter and toward the melting period (Koivusalo & Kokkonen 2002; Marks et al. 2008).

Snow is among the most important factors affecting the distribution of plants in alpine and arctic areas (Vestergren 1902; Billings & Bliss 1959). The spatial variability in the distribution of snow affected by the tree and canopy distribution has been shown to play a role in ground vegetation distribution even in a homogeneous Scandinavian pine forest (Rasmus et al. 2011). The effect of tree and canopy distribution on ground vegetation may be direct or indirect through the resulting variability in the snow cover (Økland et al. 1999; Rasmus et al. 2011). Under the snow, the intensity of photosynthetically active radiation (PAR) is low, as solar radiation is effectively absorbed by the snow (Curl et al. 1972; Gerland et al. 1999). Varying incoming shortwave radiation fluxes and snow depths within a forest lead to spatial variability in snow–soil interface temperature and in PAR, both above and below the snow, at the dwarf shrub level.

Snow modeling can be an efficient tool in hydrological and ecological applications. The Swiss SNOWPACK model (Bartelt & Lehning 2002; Lehning et al. 2002a, b) is an established and widely used model of snow mass and energy balance. It is one of the few existing snow structure models and can estimate with good temporal resolution the layered structure within the snow cover and the physical properties of the layers (grain size, form and bonding, snow temperature, density and hardness, volumetric fractions of ice, liquid water and air). A canopy sub-model recently included in the model enables estimating the forest snow properties. The use of the SNOWPACK in the forested areas is nevertheless limited because the canopy module has not been extensively tested.

The first aim of this study is to compare the wintertime ranges in incoming shortwave and longwave radiation and temperature regimes between two different types of southern boreal coniferous forest and within a single forest stand. The second aim is to use these observations to validate the ability of the model SNOWPACK to simulate the energy balance terms and temperature regimes of the snow below these different canopies. The effect of the choice of the canopy parameters on simulation quality is also addressed as part of the validation.

Two sites with fairly flat topography in a pine and a spruce forest were chosen for the work. The study sites are located close to each other in Evo, southern Finland (61° 11′ N; 25° 5′ E, 130 m above sea level) with a mean annual temperature of 3.9 °C and mean annual precipitation of 630 mm. Locations for the radiation and temperature measurements within the forests were carefully chosen to ensure the representativeness of the observations when only a limited amount of observational devices – with fixed locations – could be used. Detailed mapping of the pine site was performed, including location of all trees taller than 0.5 m, mean diameter at breast height and mean tree height of the trees taller than 2 m, and basal tree area of the site. The horizontal canopy dimensions of trees taller than 20 m were mapped with a modified Cajanus tube by locating the edge points of the canopy in eight different directions for each tree (Rautiainen et al. 2005). Based on these measurements, indices of tree influence on different points on the forest floor were calculated. The one-sided leaf area index (LAI), and corresponding canopy openness (or sky view fraction (SVF)) of the site were observed at a 1 m resolution using a Licor LAI-2000 Plant Canopy Analyzer along a 40-point transect at the mid-line of the site. Combining the information gained from mapping and LAI/SVF observations, locations representing sparse pine canopy (LAI 1.58 and SVF 0.33) and dense pine canopy (LAI 1.93 and SVF 0.19) were selected. The measurement location below the spruce canopy was based on the LAI and SVF transect observations (LAI 3.80 and SVF 0.15), representing spruce canopy with an average coverage. Both sites were estimated to have a mean canopy height of 14 m. A more detailed description of the spatial distribution of trees, canopy, ground vegetation and snow at the pine site can be found in Rasmus et al. (2011).

Snow depth was observed manually twice a month between January and April (as a mean of a 40-point transect along the mid-line of the sites) at the pine and spruce sites during winter 2008/2009, and monthly at the sparse pine and the spruce locations during winter 2009/2010.

SNOWPACK is a one-dimensional model for snowpack structure, mass and energy balance, developed in the Swiss Federal Institute for Snow and Avalanche Research (SLF) for avalanche warning purposes. SNOWPACK is a predictive model that uses Lagrangian finite elements to solve for heat and mass transfer, stresses, and strains within the snow cover. The model is physically based: energy balance, mass balance, phase changes, water and water vapor movement are included, and the layer calculations are based on snow microstructure (crystal size and form, bond size, number of bonds per crystal). A complete description of the model can be found in Bartelt & Lehning (2002) and Lehning et al. (2002a,,b).

Recently a canopy sub-model has been added to the SNOWPACK model to simulate the impact of vegetation on the upper boundary conditions of the underlying snow cover. The canopy sub-model simulates in particular the transmission of shortwave and longwave radiation, turbulent heat exchange and precipitation interception and throughfall. The details of the canopy model can be found in previous publications (Lehning et al. 2006; Stähli et al. 2009; Musselman et al. 2012a, b; Rasmus et al. 2013) and only a summary is given here with special attention to the representation of radiation and interception processes and the influence of LAI and SVF. The canopy sub-model simulates the temperature and the storage of intercepted precipitation (frozen and unfrozen) in a single-leaf vegetation layer, characterized by three structural parameters: mean vegetation height, LAI and SVF. The canopy temperature is calculated by solving a canopy energy balance equation where the canopy net radiation is balanced by the sum of sensible and latent heat fluxes between the canopy and the atmosphere (the model version used in this study does not take into account any heat storage in the canopy). The basis for the radiation transmission model is the absorption of radiation as an exponential function of LAI following Beer's Law. The absorption coefficient for shortwave radiation includes the impact of solar elevation angle for direct beam path length through the canopy (Chen et al. 1997) and SVF (Gryning et al. 2001), multiple reflections between the canopy layer and the snow or soil surface below (Taconet et al. 1986), as well as the impact of intercepted snow on the canopy albedo as described in Stähli et al. (2009) and Musselman et al. (2012a, b). The transmission model for longwave radiation follows the same principles as for diffuse shortwave radiation, with constant SVF and radiation absorption as a function of LAI in the canopy covered fraction (1-SVF). The downward longwave radiation at the snow surface is a sum of radiation from the canopy as function of canopy temperature, and atmospheric radiation either transmitted through the canopy or directly through the gaps in the canopy corresponding to the fraction SVF. The parameter values for LAI and SVF are further used in the calculation of snow and rain interception and for the turbulent heat exchange between the snow surface and the atmosphere. The direct throughfall fraction is assumed equal to the SVF. The remaining precipitation is intercepted following the snow interception model suggested by Pomeroy et al. (1998) where the interception capacity is linearly related to LAI. The intercepted precipitation is evaporated as part of the canopy latent heat flux. Intercepted snow will be unloaded when air temperature is increased above 0 °C, when the snow interception capacity is assumed to be drastically reduced. The turbulent heat exchange is calculated using aerodynamic resistances between the atmosphere and the canopy air layer following the usual bulk-formulation based on displacement height and surface roughness length. Additional resistances are added for the fluxes from the canopy air layer to the canopy layer and the snow/soil layer, respectively, following the simplified two-layer model suggested by Blyth et al. (1999). The resistance between the snow/soil layer and the canopy layer is further increased as an exponential function of LAI (details available in Musselman et al. 2012a, b). SNOWPACK has been validated in several studies in varying climatic conditions (e.g. Lehning et al. 1998; Lundy et al. 2001; Rasmus et al. 2007); the canopy radiation transmission sub-model has been evaluated by Stähli et al. (2009) and Musselman et al. (2012a, b). In the present study, the snow cover part of the model was applied without calibration to local conditions in order to detect how well the radiation and temperature conditions could be reproduced by running the model with various canopy parameterizations.

The SNOWPACK model uses air temperature (°C), relative humidity (%), wind velocity (ms⁻¹) and direction (°), incoming shortwave and longwave radiation (Wm⁻²) and precipitation (mm), either above the canopy or from an open area as meteorological input data. For the model runs in the present study, input data were obtained from the Finnish Meteorological Institute (FMI) station near the open field site, as the study data for the open site were limited to the winter 2007–2008. This observational data were nevertheless used for input data sensitivity tests. Synoptic meteorological observations with 3 h time resolution were used together with the FMI grid estimates of the incoming shortwave radiation of the site. Incoming longwave radiation was estimated using the difference between potential (theoretical value of incoming shortwave radiation not taking cloudiness into account) and observed incoming shortwave radiation, air temperature and relative humidity of each 3 h time step (the method is described in Konzelmann et al. (1994)).

The model was run for the open, pine (sparse and dense) and spruce sites for the study winters between October and April. The model outputs analyzed included snow–soil interface temperature and energy fluxes at the snow surface.

Even though the LAI values of the sites are observed using an optical method (Licor LAI-2000 Plant Canopy Analyzer), the commonly used correction factors may increase the value of LAI significantly. It is also possible to estimate LAI from hemispheric photos, using biomass based estimates or with remote sensing techniques. Further, SVF can be either determined optically (using Licor LAI-2000 Plant Canopy Analyzer) or by empirical relationships between LAI and SVF. There is no agreement yet among snow modelers on which method yields the best estimates of commonly used canopy parameters. It has been shown that the model output is sensitive to the canopy parameters (Rasmus et al. 2013), and furthermore, it is not known which parameter combination will reproduce the below-canopy radiation conditions most reliably.

In this study, we used six possible combinations for LAI and SVF (Table 1). In the ‘optical effective’ parameterization the values obtained by using the LAI-2000 were used for LAI and SVF without correction. In the ‘optical corrected’ parameterization the LAI values were corrected for needle clumping by a correction factor of 1.65 (Stenberg 1996). Two experimental functions, the ‘Pomeroy SVF’ and ‘Canadian land surface scheme (CLASS) SVF’ parameterizations, were also used to estimate the SVF from either effective or corrected optical LAI values, respectively (Verseghy et al. 1993; Pomeroy et al. 2002).

Equation (1) is used, for example, in the CLASS for global circulation models, and it was originally developed for estimating long wave radiation in coniferous forests (Verseghy et al. 1993). The land surface scheme of the Swedish Rossby Centre regional climate model uses Equation (1) for all types of forest and for the estimation of both shortwave and longwave radiation (Samuelsson et al. 2006).

The absolute differences among the sites increased toward the spring, especially in the monthly mean incoming shortwave but also in the longwave radiation (Table 2). The pattern is less clear in the relative differences. Maximum differences between the study environments are seen during the afternoon hours, both in absolute and relative terms (Table 3).

Comparison of the absolute differences between the simulated and observed mean monthly snow–soil interface temperature and incoming radiation fluxes also indicated that the use of the corrected LAI values better approximates the observations than the use of the uncorrected effective LAI (data not shown). SVF could be either optically estimated or based on empirical equation. In most of the cases the error in incoming shortwave radiation increased toward the spring. If simulated values were higher than the observed ones during the winter months, the situation may be reversed during April. These seasonal patterns indicate a model deficiency related to solar elevation angle. The solar elevation angle was taken into account in the model, both in the direct beam absorption coefficient and in the modification of the direct beam SVF. In both cases, the model parameterization should increase the shortwave transmission as a function of increasing solar angle. It should also be noted that most of the seasonal deviation in downward shortwave radiation was reduced by using the optical corrected LAI and the higher SVF estimates in the model. The impact of local gaps in the forest canopies and the corresponding sensitivity to solar elevation angle and time of the year could however be further taken into account, either by calibration of the canopy transmissivity parameters or by detailed canopy gap fraction modeling using hemispherical photos following Musselman et al. (2012a, b). The impact of local gaps was also pointed out by Stähli et al. (2009) using below-canopy radiation measurements along a 10 m rail. Error in incoming longwave radiation was more constant over time, with a systematic underestimation of about 10–20 Wm⁻² independent of forest stand but with a clearly better results for simulations with the optical corrected LAI and as small SVF as possible. Snow–soil interface temperature was highly sensitive to the choice of parameterization only in the spruce forest case; quality of these simulations was linked to the quality of snow depth simulations (Figure 5). A general decrease was seen in the differences toward the spring. When comparing the sparse pine forest cases during winters 2008–2009 and 2009–2010, differences between the simulations and observations were generally larger during the winter 2009–2010 with thicker snow cover.

The observations on near surface radiation and temperature conditions in southern boreal forests with different canopy conditions reinforce previous findings on the importance of canopy and canopy distribution on winter time snow surface mass and energy balance (Hedström & Pomeroy 1998; Lundberg & Koivusalo 2003; Stähli et al. 2009).

The effects of tree species were clearly seen in the observations. Incoming shortwave radiation values were higher in the pine forest compared to the spruce forest at all times, and because of the dense canopy structure, the spruce forest had higher incoming longwave radiation levels below the canopy than the pine forest. Differences in incoming longwave radiation levels between the pine and spruce sites are caused by notable differences in LAI and SVF between the sites, but they may also may be due partly to structural differences, such as branch orientation, in the canopies of different tree species. The forest snow–soil interface temperatures were higher below the sparse pine canopy during the snow season, but the ground warmed more efficiently below the dense canopy after snow melt. This can be explained by the slightly thicker snow cover (and slightly delayed snow melt) below the sparse canopy. Higher levels of incoming shortwave radiation were observed below the sparse pine canopy, and a slight tendency toward lower incoming longwave radiation flux was also seen.

Differences in observed radiation flux values, both absolute and relative, were highest during the afternoon hours of the day and increased toward the spring. This suggests that a relatively small error in the canopy cover estimate for a hydrological or climate model may lead to large uncertainties in the surface flux estimates during those times.

Canopy parameter values had a great effect on the quality of sub-canopy snow mass and energy balance simulations by the SNOWPACK model. The model output is especially sensitive to the canopy density and the tree species (characterized in the model by LAI and SVF); aspects like leaf orientation or canopy structure, varying between the species, are not taken into account. In most of the cases, the simulation results of radiative terms, snow depth or snow–soil interface temperature had the best agreement with the measurements when using a combination of optically obtained LAI values with needle clumping correction (Stenberg 1996) and either optical or empirical SVF values as a canopy parameterization.

The discrepancies seen in some of the snow–soil interface temperature simulations are most probably linked with the occasional difficulties to correctly simulate the snow accumulation and canopy interception. Owing to unstable weather during autumn, it is normal at the study sites for the snow cover to form and melt several times, before settling rather late during the early winter. Problems in the snow cover formation timing can cause discrepancies between simulated and observed snow–soil interface temperatures. Comparison with snow accumulation data from a nearby FMI station showed that the model captured quite well the snowless and snow covered periods during the early winter, although the snow depths may not be precise. There could be couple of days mismatch in the final formation date. Canopy and canopy parameterization affected the melt rate and melt date of the simulated snow cover, but although the canopy estimates had effects on snow accumulation, formation date of the simulated snow cover was not highly sensitive to the canopy. The bottom boundary conditions of the SNOWPACK model, including the rough parameterization of the soil column properties, especially of heat conductivity, and perhaps an unrealistic initial temperature profile may also have affected the simulation results. There is also naturally uncertainty in manual snow depth observations.

Pomeroy et al. (2009) showed that the level of downwelling longwave radiation below the canopy, observed in their study, was enhanced because the canopy, warmed by absorbed shortwave radiation, re-radiated the energy absorbed as longwave radiation. Thus, using air temperature as an approximation of the canopy temperature leads to an underestimation of the energy balance term related to longwave radiation below coniferous canopies. In a previous study using the SNOWPACK model, which simulates the canopy temperature, a problem with the longwave transmissivity simulations was shown to be associated with too high canopy temperatures and overestimation of the longwave radiation (Stähli et al. 2009). In contrast to the latter study, the results presented here showed an underestimation of longwave radiation which was rather constant in time, around 10–20 Wm⁻² on a monthly basis, without any clear seasonal dependence. The use of corrected LAI values (higher values compared to the effective LAI) and the smallest SVF values reduces the relative importance of the atmospheric longwave radiation and increases the importance of the longwave radiation emitted by the canopy. The simulated longwave radiation below the canopy was systematically improved using this parameterization, without increasing seasonal deviations. This suggests that the underestimation of below-canopy longwave radiation could also be a function of errors in the estimation of atmospheric longwave radiation, and not only related to seasonal or diurnal dynamics in the canopy longwave radiation emission.

For snow–soil interface temperatures, the differences in model performance were rather small compared to the magnitude of the underestimation during the cold winter months. There was a tendency that the simulations with higher SVF had a larger underestimation of snow–soil interface temperatures as well as slightly larger underestimation of longwave radiation, suggesting that the longwave radiation balance is governed more by the SVF than by the LAI, and that the best SVF estimate with regard to this particular longwave model was given by the Pomeroy function, followed second by the optical estimates and thirdly by the CLASS model. However, for shortwave radiation, the optical measured SVF gave the best results followed by the CLASS function and the Pomeroy function last. The snow depth simulations were more influenced by the SVF parameterizations in years with less snow, and more by the LAI estimates in years with more snow. This is logical since the direct throughfall and interception capacity are directly linked to SVF and LAI, respectively, in the model. In snow-rich years, simulated snow depths were closer to observations when using the optical effective LAI, which is in line with Hedstrom & Pomeroy (1998) who suggested that effective LAI should be more related to snow interception capacity, whereas the double-sided LAI would be more relevant for liquid water interception. In years and sites where SVF were more important for the similarity between observed and simulated snow depth, the best simulations were obtained with the higher SVF values, suggesting that direct throughfall is exceeding the SVF estimated by optical measurements. In summary, these results suggest that LAI corrected for clumping is the most appropriate for radiation transmission modeling, whereas effective LAI is most appropriate for snow interception modeling. It is also clear that the seasonal behavior of canopy shortwave radiation transmission, most likely linked to solar elevation and impact of local gaps in the canopy, need to be better described to improve modeling of local conditions, such as those represented by the observations in this study.

The seasonal differences in simulation and observation of shortwave radiation transmission were to a large degree reduced by using the LAI corrected for needle clumping and the higher estimates of SVF. The diurnal differences in simulated and observed shortwave radiation were more likely an effect of unsymmetrical distribution of local gaps in the forest canopies, perhaps related to adaption of the trees to the dominating direction of shortwave radiation coming from the south at low solar angles at these high latitude Finnish sites, or changes in the branch orientation due to stem moisture conditions. The differences in diurnal longwave radiation between simulations and observations could instead be related to the lack of canopy heat storage in the model, since a canopy heat storage would decrease and delay the diurnal canopy temperature cycle, and the corresponding emission of longwave radiation from the canopy. The model version used in the current paper did not include canopy heat storage, and this hypothesis could not be tested. However, in a recent study by Gouttevin et al. (2015), canopy heat storage was introduced in the SNOWPACK canopy model, which improved both canopy trunk temperatures and below-canopy longwave radiation.

Effect of time resolution was very small in our results. Also the use of FMI grid estimates of the incoming shortwave radiation did not seem to play a significant role in interpretation of the simulation outputs below the canopy. The situation was different with the longwave radiation. The use of observed radiation from the nearby open field significantly improved the below-canopy longwave radiation simulations compared to the use of the FMI grid estimate of the shortwave radiation, and longwave radiation derived from it. All of the simulation outputs underestimated the observed below-canopy radiation levels; the underestimation was less in the simulations that used the observed radiation as input.

The effect of the forest canopy is, at its best, described in simple ways in snow, hydrological and climate models (e.g. Essery 1998; Rutter et al. 2009). Variation between different forest types is seldom taken into account, and consideration of variation inside a forest is normally omitted. Additional variations among the results of simulations are introduced by the different possible values for the commonly used canopy parameters LAI and canopy openness (SVF). It is also clear that the study was limited (a single set of instruments within each canopy). Nevertheless, detailed mapping of the site, including mapping of LAI and SVF (presented in Rasmus et al. (2011)) facilitated finding representative locations for the instruments, and even though the results cannot acknowledge the small-scale variability of the energy balance terms, the study gives a realistic estimate of the range of the radiation levels experienced in these southern boreal forests. One way forward is an ensemble-type approach, where a set of LAI and SVF is given as model input, and a corresponding set of model outputs is produced, expressing the natural variability in radiation (and snow) conditions within the forest.

Spatial variability of the boreal forest canopy, leading to variability in seasonal snow cover through spatial differences in snow mass and energy balance, may have significant ecological consequences. Evergreen dwarf shrubs such as lingonberry have been shown to retain their photosynthetic capacity through the winter, although the photosynthesis below the snow is strongly light-limited (Lütz et al. 2005; Lundell et al. 2008, 2010; Saarinen & Lundell 2010). Depending on the snow and light conditions, subnivean photosynthesis has been estimated to compensate for a significant part of the daily respiratory losses of these plants even at the low PAR levels prevailing under the snow (Starr & Oberbauer 2003). Spatial variation in snow cover, leading to variation in temperature and PAR regimes, can thus affect species composition on a range of scales, and be a factor in shaping vegetation both at the landscape and at the stand level.

Hanna Leisti, Olli-Kalle Kauppinen and Tuomas Niskanen assisted us with the fieldwork. The Lammi biological station and the HAMK University of Applied Sciences provided logistic support. Heikki Hänninen and Helena Åström from the PECC-group are acknowledged for their valuable support on this project. Our research was funded through the Academy of Finland (Project numbers 115532 and 122194), the Environmental Research Center of the University of Helsinki HERC, the Finnish Cultural Foundation, the Maj and Tor Nessling Foundation, the Jenny and Antti Wihuri Foundation and Svenska kulturfonden.