A new Snow-SVAT to simulate the accumulation and ablation of seasonal snow cover beneath a forest canopy

Abstract

A new snow—soil—vegetation—atmosphere transfer (Snow-SVAT) scheme, which simulates the accumulation and ablation of the snow cover beneath a forest canopy, is presented. The model was formulated by coupling a canopy optical and thermal radiation model to a physically based multi-layer snow model. This canopy radiation model is physically based yet requires few parameters, so can be used when extensive in situ field measurements are not available. Other forest effects such as the reduction of wind speed, interception of snow on the canopy and the deposition of litter were incorporated within this combined model, SNOWCAN, which was tested with data taken as part of the Boreal Ecosystem—Atmosphere Study (BOREAS) international collaborative experiment. Snow depths beneath four different canopy types and at an open site were simulated. Agreement between observed and simulated snow depths was generally good, with correlation coefficients ranging between r² = 0.94 and r² = 0.98 for all sites where automatic measurements were available. However, the simulated date of total snowpack ablation generally occurred later than the observed date. A comparison between simulated solar radiation and limited measurements of sub-canopy radiation at one site indicates that the model simulates the sub-canopy downwelling solar radiation early in the season to within measurement uncertainty.

List of Symbols

Introduction

Snowmelt is a vital source of water for many parts of North America, northern Eurasia, Austria and Switzerland and many other mountainous areas, and also carries nutrients in solution and suspension. Prediction of the timing and magnitude of meltwater runoff is crucial for efficient management of this natural resource. The effects of the forest cover on the timing of snowmelt are complex, and elements of these effects have been studied in detail over many decades (see Murray and Buttle, 2003, Table 2, for a summary of the results of studies of accumulation differences between open sites and forest stands).

Snowmelt beneath a forest canopy is generally observed to occur later or more slowly than at an open site (Federer and Leonard, 1971; Hardy and Hansen-Bristow, 1990; Metcalfe and Buttle, 1995; Koivusalo and Kokkonen, 2002). Rouse (1984) noted more rapid snowmelt beneath the canopy from enhanced thermal emission, but later ablation than at an open site as a result of increased accumulation in the forest. Suzuki and others (1999) found that whilst net radiation was greater beneath a forest canopy, the contribution of sensible heat enhanced snowmelt rates at an open site. Others have found that the relative timing of snowmelt at forested and open sites varies between years, between sites or with different atmospheric conditions (Yamazaki and Kondo, 1992; Murray and Buttle, 2003). Further investigation of the energetic and mass effects of the forest on the snow cover, through both field measurements and the use of physically based models, is needed to determine the conditions under which a particular forest will generally enhance or inhibit the rate of snowmelt. These studies are also important for the representation of the forest-snow ecosystem in larger-scale models.

Recent studies have shown that within general circulation models (GCMs) the representation of forests with underlying snow cover has a crucial impact on the model results (Viterbo and Betts, 1999; Betts, 2000). This is particularly important for boreal forests, which cover 15% of the Earth’s land surface and have snow on the ground for 68 months of the year.

The interaction of radiation between the forest canopy and snow cover is complex, as the radiation balance affects the snowpack properties including snow grain-size, which in turn affect the sub-canopy radiation balance through albedo change. However, this interaction is not simulated in GCMs, and has not been studied in detail, even at the point scale.

Current snow-soil-vegetation-atmosphere transfer (Snow-SVAT) schemes have awide range of complexity, from single- ordual-layer snow models such as UEB (Tarboton and Luce, 1996) or SNOBAL (Marks and Winstral, 2001) to the multi-layer representation of snow by Hardy and others (1998), who linked the Geometric-Optic Radiative Transfer model, GORT (Li and others, 1995), with the physically based multi-layer snow model SNTHERM (Jordan, 1991).

A common approach is to represent the canopy radiative effect with a two-parameter Beer’s law formulation that requires calibration or extensive field measurements. Whilst this method may give impressive results at individual sites, use of a physically based model applicable to all sites would be an improvement. The approach taken by Hardy and others (1998) has led to probably the most advanced Snow- SVAT to date, with a geometrical-optical model for solar radiation between tree crowns and a radiative transfer model for solar radiation within tree crowns (with down- welling thermal radiation estimated by a sky-view factor weighted average of the atmospheric and canopy thermal radiation). However, a large number of parameters are required to initialize GORT, and these are difficult to measure in the field. In addition, the two models are not coupled, so all solar radiation data are processed in bulk and are then used to drive SNTHERM independently. Hence multiple scattering of solar radiation between the canopy and snow surface is not considered.

The approach of the model, SNOWCAN, described in this paper, is to use a minimum number of parameters to represent the absorption and scattering of optical and thermal radiation by the canopy, with the requirement that these parameters may easily be obtained in the field or by remote sensing. The canopy radiation model is physically based and accounts for multiple scattering of radiation, yet is computationally efficient. Here the radiation model has been coupled to a complex multi-layered snow model to simulate the evolution of snowpack stratigraphy, but the canopy radiation model could be coupled to a simpler snow model for incorporation within larger-scale models, where computational efficiency is vital.

The three-parameter canopy radiation model has been coupled to a multi-layer physically based snow model in such a way that the radiation balance affects and is affected by the change in snow structure. This offers a distinct advantage over current Snow-SVATs, which pre-process seasonal solar radiation data for canopy effects, and therefore cannot model the snow albedo feedback process effectively. In addition, the absorption, emission and scattering of thermal radiation by the canopy is considered explicitly in SNOWCAN. Other forest effects, such as the interception of snow by the canopy and reduction of wind speed, have also been incorporated in the coupled model.

SNOWCAN was used to simulate snow accumulation and ablation beneath four different canopy types (aspen, jack pine, black spruce and mixed) and at an open site within the Canadian boreal forest, which were studied as part of the Boreal Ecosystem-Atmosphere Study (BOREAS) experiment (Sellers and others, 1997). Simulated snow depths at these sites generally showed good agreement with automated or manual measurements of snow depth, although the final rate of snowmelt was underestimated by the model.

Sub-canopy solar radiation beneath a mature jack pine canopy was simulated to within the uncertainty of the rather limited measurements made. However, this study highlights the need for an extensive multi-site, multi-season experiment with a carefully collected dataset that includes sub-canopy solar and thermal radiation measurements that are spatially and temporally concurrent with turbulent flux and snowpack measurements.

Model Description

Snow model

The snow model used in SNOWCAN is a multi-layer, onedimensional model, based on the physical principles behind SNTHERM (Jordan, 1991). SNOWCAN is based on a La- grangian adaptive grid that deforms with the compaction of the snow cover. Layers are added to the model as precipitation falls, and mass and energy transfer is calculated between the layers and across the model boundaries. Principles of energy, mass and momentum conservation are applied to the snowpack to derive equations to simulate temperature and mass changes over the season. A full description of SNOWCAN is given in Tribbeck (2002) and is summarized below. Fluxes within the snowpack are defined as positive upwards, whereas surface fluxes are defined as positive downwards by meteorological convention. A control volume schematic of the mass and energy fluxes is shown in Figure 1.

Fig. 1. Control volume schematic of mass and energyfluxes.

Conservation of energy

The energy balance for the snowpack is given by Equation (1), where the enthalpy is represented in terms of the temperature (Equation (2)), and Equation (1) is solved to determine the temperature of each snow layer.

(1)

(2)

At the snow surface, a Neumann boundary condition (specification of the derivative of the solution) is used to specify the transfer of heat between the snowpack and atmosphere, as described in Equation (3).

(3)

A Dirichlet boundary condition is used to define the energy transfer at the lower boundary, where a constant temperature is assumed at some depth in the soil. Errors associated with this assumption are discussed later.

Conservation of mass

Conservation of mass is applied to the total medium and to the constituents separately (ice, liquid water and water vapour). The last term of the mass-conservation equation (Equation (4)) vanishes for the total medium since M _k,k′ = -M _k′.k.

(4)

At the surface, mass is added to the snowpack from precipitation:

(5)

where the density of new snow is determined from the wet- bulb temperature. Mass is also added to, or removed from, the surface by turbulent exchange:

(6)

The turbulent exchange coefficients were first defined in Jordan (1992), although these were subsequently improved by Jordan and others (1999). Windless exchange coefficients are used to maintain minimal turbulent transfer at zero wind speeds to ensure model stability (Jordan and others, 1999).

Liquid-water and water-vapour flow within the soil are not currently modelled in SNOWCAN. Liquid water and water vapour are assumed to run off horizontally at the snow-soil interface. Temperature and phase changes within the soil are simulated in SNOWCAN, hence soil water content is required for model initialization. The assumption of horizontal runoff at this interface is discussed later.

Conservation of momentum

Liquid water drains through the snowpack relative to the ice matrix. In SNOWCAN, liquid water flow is assumed to be laminar and unsaturated. Inertial, convective and phase- change terms in the momentum balance and capillary forces are assumed to be negligible, so the flow of liquid water may be represented by Darcy’s law:

(7)

In reality, liquid-water percolation is complex, and the assumption of horizontal heterogeneity applied here is certainly not valid. Heterogeneities such as rapid drainage channels or ice lenses in the snowpack can enhance or delay melt percolation. Marsh (1991) discusses these effects and their implications in depth.

Snow compaction and grain growth

The algorithms used to represent snow densification in SNOWCAN were developed by Jordan (1991), based on the work of Anderson (1976). Parameterization of these algorithms applied in SNOWCAN differs slightly from that used in SNTHERM (Jordan, 1991), as described later.

Snowpack settling and metamorphism under a small or negligible temperature gradient tend to lead to spherical grains that pack more efficiently. Densification of the snow- pack by this mechanism is represented by grid compression, according to Equation (8), where the parameter ξ is calibrated for each site.

(8)

This process is limited by increasing snow density. Above a critical density, the rate of metamorphism compaction is reduced by a factor e^{−0 046(γi—γcrit)}, where the critical density is determined by calibration at each site but is typically in the range 100−150 kg m^-3. Wet-snow metamorphism is faster, and the rate of change of layer width given by Equation (8) is multiplied by a factor 2.

Densification of the snow cover also results from the overburden pressure of the snow, P _s. This much slower process is described by:

(9)

where η₀ is the viscosity coefficient at the melting point and at the limit of zero snow density and p _s is the snow density. The viscosity coefficient is determined by calibration for each site.

Beneath a forest canopy, snow densification is also affected by the impact of snow, ice and liquid water unloaded from the canopy. This process is spatially and temporally heterogeneous, and because of its complexity has not been incorporated within this model. The value of the densifica- tion parameters determined by calibration may reflect additional compaction by this process.

Grain growth occurs by the hand-to-hand vapour diffusion method described by Yoshida (1963), and is controlled by the size of grains, the vapour flux and the liquid-water content of the snowpack. The rate of grain growth in SNOWCAN is represented by the following equations:

(10)

(11)

(12)

Canopy radiation model

The optical and thermal radiation canopy model, RM (Pearson and others, 1999), approximates the canopy as a turbid medium with the statistics of a spherical leaf-angle distribution function. Radiation incident on the canopy elements is absorbed or scattered isotropically, and multiple scattering is simulated between the canopy and snow surface. Thermal radiation is treated in the same way as solar radiation. This canopy radiation model is physically based, yet is simple in its implementation and as a result is computationally efficient, and requires minimal extra parameters that are commonly measured or well known. These are its advantages over existing models that may represent some processes more directly but are unusable without comprehensive field measurements.

Radiative transfer theory is not applicable to discontinuous canopies (Li and others, 1995) as assumptions of horizontal homogeneity break down, so this model is only appropriate for continuous canopies. Roujean (1996) reported that the effect of vegetation clumping was to enhance the probability of light penetration relative to a random leaf distribution. However, there are good reasons for using such a simplified representation. Firstly, RM’s level of detail and realism is concomitant with that of the snow model. Secondly, a more detailed radiation model cannot be initialized or parameterized adequately from data that are generally available or can reasonably be collected. The positions, shapes and optical properties of nearby trees and their components are essentially unknowable without extensive, timeconsuming field measurements, which are not practicable for most field efforts. Thirdly, RM is less sensitive to the structure of the vegetation than one might expect (Pearson and others, 1999), and the simplifications embodied in the model leave it, nevertheless, as a physically based description of the thermal and optical properties of a two-layer system, containing adjustable parameters that enable the response of a dense canopy to be modelled successfully. Given the model limitations, RM would not be appropriate for canopies with high heterogeneity, such as large canopy gaps, or with very low canopy density (leaf area index (LAI) < 1).

Figure 2 is a schematic of the canopy radiative fluxes, and demonstrates the multiple scattering of radiation between the canopy and the underlying surface. Radiative transfer theory is used to determine the transmission and reflection coefficients for isotropic radiation (defined in Equation (18) and (19)). Here, solar and thermal radiative fluxes are treated similarly, but with different reflection and transmission coefficients. Given these coefficients, the total sub-canopy downwelling and upwelling isotropic radiative fluxes can be determined by Equations (13) and (14). This formulation is simple yet enables multiple scattering between the canopy and underlying surface to be considered with negligible increase in computational efficiency.

(13)

(14)

Fig. 2. RMmodel canopy fluxes.

Direct solar radiation is approximated as a delta-function beam in the direction of the sun. The interception of radiation depends on its path length through the canopy and hence its direction.The penetration of direct solar radiation is modelled with change in sun angle, as given by Equation (15), where the vegetation optical depth is a function of the LAI as discussed later in this section.

(15)

The first term in Equation (15) accounts for the proportion of radiative flux that is not intercepted by the canopy, whilst the second term calculates the proportion of intercepted radiative flux that is scattered downwards. Intercepted direct solar radiation that is scattered from the canopy, non- intercepted radiation that is subsequently reflected from the underlying surface, diffuse solar radiation and thermal radiation are assumed to be isotropic. The transmission and reflection coefficients required to solve Equations (13) and (14) arise from the solution of the total radiant flux (M _iso) exiting from a planar surface, as described by Pearson and others (1999). With the approximations that diffuse radiative flux is isotropic and that the flux is exponentially attenuated in the canopy, the total flux exiting the bottom of the canopy is determined as:

(16)

(17)

The incident flux is Iiso, hence the intercepted flux is determined as πI _iso [1-2E₃(τV)]. The reflected radiative flux is assumed to be scattered equally upwards and downwards, so the reflection coefficients are:

(18)

(19)

The function f(τV) = [1-2E ₃(τV)] has the following properties:

(20)

The second derivative condition ensures that changes in the scattering properties with canopy structure are rapid at low vegetation concentrations, but slower at higher vegetation concentrations. This corresponds physically to the higher probability of leaves overlapping at higher leaf concentrations, and therefore having a lesser effect on the radiation transfer.

The dependence of the transmission and reflection coefficients on the canopy structure arises from the relation between the vegetation optical depth and the canopy LAI. The LAI was chosen to represent the canopy because it is a standard measurement, easily obtained from ground-based measurements or from remote sensing. The optimal method of LAI measurement is determined by the areal extent of the system represented in the simulation and the available validation data.

For a spherical leaf distribution, it can be shown that the vegetation optical depth is related to the LAI as τv = A _L/2 (Monteith and Unsworth, 1990). However, this assumption may not be appropriate for canopies with a high degree of clumping, so other distribution functions may be used to determine the optical depth dependence on LAI, or other empirical formulae used where the radiative interception efficiency of the canopy is known.

As described previously, diffuse solar radiation is assumed to be isotropic from the sky, and direct solar radiation is approximated as a delta-function beam in the direction of the sun. These assumptions were made for simplicity of calculation, but may be a source of error, since diffuse shortwave flux is observed to be non-isotropic, particularly at low sun angles. Atmospheric thermal radiative flux is required to drive the model, whereas canopy and snow surface thermal radiative fluxes are calculated using Stefan’s law. The relative contributions of the atmosphere, canopy and snow to the sub-canopy thermal radiative flux depend on the canopy structure.

Additional canopy effects

In addition to the radiative effects of a canopy, other forest processes affect the energy and mass of an underlying snow- pack. SNOWCAN also accounts for the reduction of wind speed by the canopy, interception of snow by the canopy, deposition of litter and reduction of surface albedo by litter. These additional effects are considered to be less important than the radiative canopy effects in the regions studied, as described in the following sections, and therefore are not treated with as much detail in the model.

Wind speeds beneath a canopy are less than those in the open, and therefore the turbulent exchange is also reduced (Price and Dunne, 1976). In SNOWCAN, the turbulent transfer mechanisms developed by Jordan (1992) for an open snow cover are also applied beneath the forest canopy, but the reduction of wind speed beneath the canopy is modelled using the empirical algorithm for the boreal canopy (w _f = w _o/5), developed by Link and Marks (1999) derived from BOREAS data.The minimum sub-canopy wind speed is limited to 0.2 ms^-1 in order to maintain model stability. This algorithm is consistent with measurements by Price and Dunne (1976) and Koivusalo and Kokkonen (2002). A fully coupled atmosphere-canopy model (including simulation of canopy temperature) will be included in the future.

Snow intercepted by the forest canopy may be retained for a period of time before throughfall occurs as a result of wind displacement, additional loading or partial melt, or is lost to the atmosphere by sublimation. Detailed measurements of snow loading on and sublimation from forest canopies have been carried out by Hedstrom and Pomeroy (1998) and Lundberg and others (1998), and a survey of the available literature by Pomeroy and others (1998) suggested that up to 60% of annual snowfall is intercepted by the canopy, and 25−45% of the annual snowfall is lost to the atmosphere by sublimation.

Here, measurements of throughfall are used to drive SNOWCAN, but when these are not available, a constant proportion of the precipitation is assumed to be intercepted by the canopy, and it is assumed that all intercepted snow is lost by sublimation. This is a simplification, but a detailed physical model of interception processes, such as that of Par viainen and Pomeroy (2000), is beyond the scope of SNOW- CAN at present. The effect of intercepted snow on canopy single-scatter albedo is assumed to be negligible.

Leaf litter and other forest debris are deposited onto the snow surface from the canopy. The debris may become buried beneath subsequent snowfall, and become reexposed as the snow melts. Forest litter within the snowpack affects the radiation balance of the snowpack. Litter absorbs more solar radiation than does the snow, and emits thermal radiation accordingly. The spectral effect of litter on snow surface albedo was examined by Melloh and others (2002).

Hardy and others (2000) simulated the deposition of litter within SNTHERM with the following equation:

(21)

where litter deposited on the surface is assumed to be retained within the snow layers as they are subsequently buried beneath fresh precipitation. This algorithm has also been implemented within SNOWCAN. The form of Equation (21) represents an increase in the probability that additional litter will overlap that already present. For a constant litter-deposition rate, the rate of change of litter coverage decreases with time. The litter is assumed to change the surface albedo by:

(22)

Snow-SVAT solution

The meteorological data required to drive the snow model are available, typically, at hourly intervals. SNOWCAN has an adaptive time-step (typically varying between 5 s and 15 min), which is chosen according to the accuracy of the solution. The forcing data are interpolated accordingly, and sub-canopy radiation fluxes are determined. The mass-balance equation is solved first, to determine the effective saturation of the snowpack. Subsequently the energy balance is solved to determine the snow temperature. Snow albedo is determined from the physical properties of the snowpack (snow grain-size) and meteorological fluxes. The new snow albedo is used to determine the sub-canopy radiation balance at the next time-step. A control volume finite- element approach is used and the governing equations are linearized, lending to the solution of the tridiagonal matrix equation.

The forcing data required to drive SNOWCAN are: incident solar (direct and diffuse) and thermal radiation, wind speed, relative humidity and air temperature, measured above the canopy or at a nearby open site, and also canopy temperature and precipitation water equivalent. Only three additional parameters are required to describe the canopy: LAI, canopy emissivity and canopy singlescatter albedo.

Site and Data Description

SNOWCAN was tested with data taken as part of BOREAS, a collaborative multidiscipline experiment designed to examine processes that affect climate change (Sellers and others, 1997). A mature jack pine (Pinus banksiana) site (SSA-OJP) within the southern study area (SSA) was selected for an initial test of SNOWCAN, and the period January−April 1994 was chosen because of the range of data collected and because previous studies have focused on this year and site (Hardy and others, 1997).

Link and Marks (1999) simulated the snowpack at this site and five additional sites for the winter period 1994/95. The two additional sites within the northern study area (NSA) have a mature jack pine canopy (NSA-OJP) and a mixed black spruce (Picea mariana) /poplar (Populus balsami- fera) canopy (NSA-YTH). Within the SSA, the three additional sites were an open site (SSA-OPEN), and sites with a mature aspen (Populus tremuloides) canopy (SSA-OA) and a mature black spruce canopy (SSA-OBS). Here, meteorological data used to drive SNOBAL (Link and Marks, 1999) were also used to test SNOWCAN for a range of canopy characteristics.

For the initial test of SNOWCAN at the SSA-OJP site in winter 1993/94, 15 min averaged meteorological data were used to drive the model. These data were measured by instruments mounted 10 m above the canopy on a tower, and precipitation throughfall was measured beneath the canopy. Data missing for 1 hour or less were linearly interpolated. Missing data over longer periods were estimated from data at the SSA-OA site situated approximately 21 km southwest of the SSA-OJP site. Snow depth at the SSA-OJP, used to test SNOWCAN, was measured automatically with an ultrasonic snow-depth gauge, which was situated in a small gap between canopy crowns. Downwelling sub-canopy solar radiation was measured over three clear days in February 1994 by nine radiometers that were moved and relevelled daily, as described by Hardy and others (1997).

Available data at the six different sites used to drive the multi-site simulations for winter 1994/95 are summarized in Table 1 (Link and Marks, 1999). Hourly meteorological data from flux towers above the NSA-OJP, NSA-YTH, SSA-OA and SSA-OJP canopies were used to drive the SNOWCAN simulations at these sites. Meteorological data from the SSA-OJP site were also used to drive simulations at the SSA-OPEN and SSA-OBS sites.

Table 1. Model parameterization differences between sites with different canopy characterizations

Model initialization and calibration

For the initial test of SNOWCAN at the SSA-OJP site in January−April 1994, the density profile on 1 January 1994 was assumed to be logarithmic with an integrated density of 220 kg m^-3. At the lower boundary, the temperature of the lowest soil layer, which was 1 m below the soil surface, was assumed to be at 1°C for the duration of the simulation. The initial snow temperature profile was assumed to be linear between the soil and air temperatures, the average grain-size of each snow layer was assumed to be 1 mm and the initial litter content of the snowpack was assumed to be negligible in this simulation.

The simulations for the following winter (1994/95) began before the start of the snow season. The initial soil temperature profile was assumed to be cold, and the model allowed to spin up to warmer soil temperatures. The model spin-up time is of the order of days, whereas the first snow deposition event occurred several weeks later. The parameters used in the simulations at each of the six sites are given inTable 1.

The LAIs were derived from measurements detailed by Chen and others (1997), with the exception of the young mixed canopy site (NSA-YTH), where no information was available. LAI for this site was assumed to be similar to that of a nearby young jack pine canopy. The canopy single-scatter albedo and emissivity were assumed to be 0.25 and 0.977 respectively. Litter deposition rates were derived from Hardy and others (1998). For turbulent transfer calculations, the measurement heights used were relative to the ground (26.8 m), and the roughness length used in these simulations was 5 mm.

Site calibration of the densification parameters was carried out by model comparison with measured snow depth for the longest period between precipitation events. These parameters were derived from the January−April 1994 SSA-OJP simulation, and were also applied to the winter 1994/95 simulations. For the simulations presented here, ξ = 6.942 × 10^-6. This parameter is greater than the value of ξ = 2.778 × 10^-6 used in SNTHERM, and reflects greater initial compaction. The viscosity coefficient is equal to 9 × 10⁷ N s m² for the data presented in this paper, which is greater than the values suggested by Kojima (1955) (1.7 × 10⁵Nsm² < η ₀ < 1.9 × 10⁶Nsm²), Anderson (1976) (η₀ = 3.6 × 10⁶ N sm²) and Brun and others (1989) (η0 = 1.0 × 10⁵N sm²), and reflects much slower densification as a result of pressure. Hardy and others (1997) also adjusted the compaction algorithms used by Jordan (1991) to correct for overestimation of snow densifica- tion, but did not state what changes were made.

Results

Simulated snow depth beneath the boreal jack pine canopy is shown for the 1994 accumulation and ablation period in Figure 3. Simulated and observed snow depths show reasonable agreement, having a regression coefficient of r ² = 0.90. Figure 3 shows differences between simulated and observed snow depths on days when significant precipitation has occurred. These differences arise from the location of the precipitation gauge (under a canopy crown, hence model input data are affected by interception processes) and the snow-depth gauge (in a small gap between crowns, hence validation data are less affected by interception processes). If the snow-depth data are used to estimate the timing of precipitation events, the agreement between measured and simulated snow depth is improved, with a regression coefficient of r² = 0.94 (simulation not shown).

Fig. 3. Simulated snow depth beneath a jack pine boreal forest canopy.

SNOWCAN simulates significant melt events on days 60−64, 70−73, 88−90 and 97−108, and captures the change in depth during these melt periods, as shown in Figure 3. However, at the end of the ablation period, the rate of melt is underestimated by the model, with persistence of the simulated snowpack several days after the observed snow- pack has disappeared.

Figure 4 shows a comparison between simulated subcanopy downwelling solar radiation and hourly averages of measured sub-canopy solar radiation. Also shown on this graph are the forcing data used to drive the model. The simulated reduction in solar radiation as a result of canopy shading by SNOWCAN is of the correct order of magnitude, although simulated sub-canopy radiation is greater than the averaged radiometer measurements, particularly at low sun angles. The average difference between simulated and measured solar radiation is 3.0 W m^-2 for this 3 day period.

Fig. 4. Simulated sub-canopy solar radiation compared to areal averaged measurements.

At the lower boundary, simulated soil temperatures were compared to measurements at 10, 20 and 50 cm depth within the soil (r² = 0.89, 0.89, 0.76 respectively), as shown in Figure 5. The general temperature trend, driven by the upper boundary conditions, is followed at all depths. Simulation of surface temperature is generally good, except in the first 20 simulation days, where the simulated soil temperature is approximately 4°C colder than the measured temperature. At greater depths within the soil, however, the simulated temperatures diverge from, and are less variable than, the observations. Peaks in the simulated soil temperature on days 63, 74 and 91 correspond to simulated melt periods shown in the snow-depth graph (Fig. 3). Propagation of the warm front downwards after melt periods was simulated with lag times similar to the observations.

Fig. 5. Comparison between measured temperatures at fixed depths within soil and simulated soil temperatures driven by a constant lower-boundary temperature.

Figure 6 demonstrates simulation of snow depth at six different sites in the boreal forest, each with different canopy characteristics. As in the 1994 simulation shown in Figure 3, the rate of ablation at the end of the season is underestimated by the model. Within the SSA, the snowpack at the open site forms before the snowpack beneath the canopies. Small periods of snowmelt occur on days 43−44,122−123 and 144−145 at the forest sites but not the open site, and the snowmelt on days 162−166 is more pronounced at the forested sites than at the open site.

Fig. 6. Snow depths simulated at six different sites with different canopy characteristics and compared with measurements.

The agreement between simulated and observed snow depth is demonstrated by the regression coefficients given in Table 2. The correlation coefficients in italics apply to periods when snow depth was measured manually on a biweekly basis. A significance test at these sites shows significant positive correlation between measured and simulated snow depths at the 1 % confidence level. Also shown in Table 2 is the first date when the ground is snow-free at the end of the season. The simulated snow melt-out date is later than the observed melt-out date in all cases, with the possible exception of the SSA-OBS site, where insufficient data are available to determine the model accuracy in this respect. For the sites where automatic snow-depth measurements are available, the maximum time difference between simulated and observed meltout date is 1 week.

Table 2. Correlation coefficient agreement between observed and simulated snow depths for winter 1994/95 simulations. Also shown are the day of year on which the actual and the simulated ground surfaces are snow-free

Despite the errors in the simulation of the final melt period, the general trend in simulated snow melt-out dates between sites with different canopy characteristics follows the observations. Within the NSA, the snow beneath the higher LAI canopy (NSA-OJP) melts before the snow beneath the lower-density canopy (NSA-YTH). This is in contrast to the SSA, where the snow beneath the increasingly higher LAI is prolonged relative to the lower canopy densities and the open site. As SSA-OPEN, SSA-OJP and SSA- OBS were driven by identical forcing data, differences between these three sites with very different canopy structure are examined in closer detail.

Changes between the simulated sub-canopy net solar and net thermal radiation contributions to the energy balance under different canopy conditions are demonstrated in Figure 7. As the canopy density increases, net solar radiative energy decreases but net thermal radiative energy input increases. At the SSA-OPEN site, thermal radiation has a mostly cooling effect, in contrast to the SSA-OBS site where thermal radiation may contribute more energy than solar energy to the energy balance. On days 25-26, 144 and 162-163, the SSA-OJP and SSA-OBS sites have strong thermal contributions to the energy balance, which are not apparent at the SSA-OA (not shown) and SSA-OPEN sites. From Figure 6 these periods correspond to the early development of the snowpack at the open site (from day 22) and the melt events at the forested sites (days 144-145 and 162-166).

Fig. 7. Contrast between simulated net thermal and net solar radiation at multiple sites within the borealforest.

The contributions of turbulent energy transfer to the energy balance of the simulated snowpack beneath different canopies are demonstrated in Figure 8. In general, the simulated turbulent transfer of energy between the air and snow surface beneath the canopy is small: less than 20 W m^-2 in magnitude, with small differences between sites. Peaks in sensible-heat transfer occur during periods of snowmelt, particularly evident at the SSA-OJP and SSA- OBS forested sites on days 144−145 and 162−166. Evaporative cooling is more pronounced during snowmelt periods when the relative humidity is not close to saturation.

Fig. 8. Simulated reduction in turbulent energy exchange beneath forest canopies.

Simulated snow water equivalents (SWEs) and liquid- water content of the snowpack at these three sites are presented in Figure 9. Early development of the snowpack at the open site contributed to greater accumulated SWE at the open site compared to the forested sites. The difference in SWE between the open and forested sites was enhanced by melt on days 43−44, 82, 122, 144−145 and 162−166 at the forested sites. Melt on days 170−178 was greatest at the SSA- OBS site, but also greater at the SSA-OPEN than the SSA- OJP site, which led to a smaller SWE at the SSA-OPEN site than the forested SSA-OJP site for the first time in the season. From day 180 onwards, the melt rate was faster at the SSA-OPEN site than at the SSA-OBS, which resulted in smaller SWE and earlier complete snow ablation at the SSA-OPEN site. The increased accumulated SWE and faster rate of snowmelt at the end of the season at the SSA- OJP site compared to the SSA-OBS site resulted in nearsimultaneous total ablation of the snowpack at both these forested sites.

Fig. 9. Change in simulated snowpack SWE under three different canopy conditions, driven by identicalforcing data.

Radiative differences between these three sites during late ablation are highlighted in Figure 10. During the day, net radiation input to the snowpack decreases with increasing canopy density. However, significant cooling of the open snowpack occurred at nighttime. This nocturnal cooling was suppressed at both canopy sites.

Fig. 10. Simulated effect of canopy density on net solar and thermal radiation during the ablation period.

Discussion

A new Snow-SVAT has been formulated to simulate the evolution of snow beneath a forest canopy. A physically based snow model and a canopy radiation model were coupled to form SNOWCAN, which uses a simplified parameterization of other forest effects, such as the sublimation of snow intercepted on the canopy, reduction of wind speed beneath the canopy and the deposition of litter on the snow surface. The approach of coupled models allows the radiation feedback mechanisms to be represented, where the sub-canopy radiation both affects and is affected by the physical state of the snowpack.

SNOWCAN is driven by meteorological variables measured either above the canopy, or at a nearby open site. The model was tested with data from the BOREAS experiment, measured above different canopies with the exception of precipitation, which was measured beneath the canopy. Simulated snow depth generally showed good agreement with the measured snow depth, although the final rate of snowmelt was underestimated by SNOWCAN. This error could result from the assumption of gravitational liquid- water flow, since rapid drainage channels are observed to form within snow. Liquid water retained in the simulated snowpack may be refrozen with nocturnal cooling, and energy input during the day must first remelt this liquid water before drainage can occur, so the simulated snowpack is prolonged. Advected heat from bare soil patches, not simulated in SNOWCAN, may also contribute to this error at the end of the season.

A comparison between simulated downwelling subcanopy solar radiation and 3 days of measurements at one site indicated that SNOWCAN simulated the sub-canopy solar radiation for this site and time of year to within measurement error, although this highlights the necessity for further data at different sites and times of year. The partition of radiation may be important in multi-layer snow models such as SNOWCAN, particularly for remote sensing applications, since the distribution of energy absorption within the snowpack affects the structure of the snowpack.

At the lower boundary, soil temperature over the January−April 1994 season at the SSA-OJP site was simulated with a good degree of accuracy at the soil surface, which demonstrates that the assumption of a constant temperature for the lower-boundary condition is viable, although the extent of temperature variation at greater depths in the soil was not simulated accurately by SNOW-CAN. The simulated temperature variations are driven by the energy transfer at the snow surface, and constrained by the specified temperature at the lower boundary. Application of this boundary at 1m depth within the soil may not be deep enough where accurate simulation of the soil thermal structure is required.

The effect of the canopy on the radiative energy input to the snowpack was clearly shown, particularly when simulations of the snow beneath three different canopy conditions driven by identical forcing data were compared. Thermal radiation from the canopy prevented early pack development, and induced or enhanced early- and mid-season snowmelt, which led to increased snow accumulation at an open site. Also greater accumulation was simulated under a thinner canopy than the more dense canopy. Interception of snow and its subsequent removal from the forest canopy has not been considered in these simulations, but could exacerbate the difference as a result of sublimation losses from the intercepted snow.

Late in the season, high sun angles and decreased snow albedo increases the radiative energy input to the open snowpack and results in a faster rate of snowmelt, despite significant nocturnal cooling of the snow surface at the open site, which was suppressed beneath the forest canopy. In general, an increase in the forest canopy density increased melt early and in the middle of the season, but decreased the melt rate late in the season, which led to persistence of the snowpack at the forested sites relative to that at an open site. The snowmelt rate beneath the forest canopies could be reduced further where litter deposition is negligible.

The data available from the BOREAS experiment are insufficient to validate SNOWCAN in full. Evaluation of the model here was governed by the availability of data. Ideally, sub-canopy upwelling and downwelling solar and thermal radiation measurements are required to test the radiation interaction between the forest canopy and snow cover. Measurements of the snowpack temperature, density and grain-size profiles are needed to investigate how well the model represents the physical processes within the snowpack. However, point simulations are not representative of areal averages because the spatial variability is large, particularly in areas where interception processes dominate. An evaluation of SNOWCAN over larger areas also needs to be carried out, and this requires extensive and continuous areal measurements of both sub-canopy radiation and snowpack properties such as depth and SWE.

This study highlights the benefits that could be provided by model validation from intensive field campaign data. Nevertheless, use of existing data in this paper has shown that this Snow-SVAT has the capacity to provide good simulations of snow depth beneath several different canopy types. This approach offers an improvement over existing models, as this model does not require intensive field measurements of canopy parameters, and the interaction of solar and thermal radiation between the snow and canopy is simulated. In addition, this approach allows simulation of the detailed structure of the snowpack, which is useful for remote-sensing applications. A Snow-SVATsuch as SNOW- CAN could also be used to improve parameterization of the forested snowpacks within larger-scale land surface models and GCMs.