On the Canadian Prairies, agricultural practices result in millions of hectares of standing crop stubble that gradually emerges during snowmelt. The importance of stubble in trapping wind-blown snow and retaining winter snowfall has been well demonstrated. However, stubble is not explicitly accounted for in hydrological or energy balance snowmelt models. This paper relates measurable stubble parameters (height, width, areal density, and albedo) to the snowpack energy balance and snowmelt with the new, physically based Stubble–Snow–Atmosphere Model (SSAM). Novel process representations of SSAM quantify the attenuation of shortwave radiation by exposed stubble, the sky and vegetation view factors needed to solve longwave radiation terms, and a resistance scheme for stubble–snow–atmosphere fluxes to solve for surface temperatures and turbulent fluxes. SSAM results were compared to observations of radiometric snow-surface temperature, stubble temperature, snow-surface solar irradiance, areal-average turbulent fluxes, and snow water equivalent from two intensive field campaigns during snowmelt in 2015 and 2016 over wheat and canola stubble in Saskatchewan, Canada. Uncalibrated SSAM simulations compared well with these observations, providing confidence in the model structure and parameterization. A sensitivity analysis conducted using SSAM revealed compensatory relationships in energy balance terms that result in a small increase in net snowpack energy as stubble exposure increases.
Snow meltwater is important for crop germination and early season growth in no-till nonirrigated farming systems commonly found in cold, semiarid agricultural regions. Stubble, the standing residue of cultivated grain and oilseed crops, is characterized by stalks that remain erect throughout snow accumulation and ablation. Enhanced snow accumulation with increasing stubble height has been well demonstrated and is due to deposition of blowing snow and suppression of wind erosion of snowfall (Pomeroy and Gray 1995). However, snowmelt models either ignore short vegetation (Gray and Landine 1987; Marks et al. 1998), assume protruding vegetation can be represented by modifying surface albedo (Liston and Hiemstra 2011), or simulate the bending over and burial of grasses and shrubs by snow (Ménard et al. 2014; Liston and Hiemstra 2011). The interactions between stubble stalks and snow occur over large areas with regional implications for hydrology and climate. Cold-region no-till crop production systems, characterized by standing stubble, are found throughout the American Midwest, Eurasian Steppe (Derpsch and Friedrich 2009), and Canadian Prairies. In the Canadian Prairie provinces (Manitoba, Saskatchewan, and Alberta), the area of no-till crop production increased from 1.7 million hectares in 1990 to 17.3 million hectares in 2016 (Statistics Canada 2016). Despite the large-scale conversion to no-till systems, the importance of snow meltwater, and the snow-trapping characteristic of stubble (Kort et al. 2011), a quantitative understanding of how the snowpack energy balance changes with the gradual exposure of stubble is lacking. How stubble may or may not influence snowmelt rates, and therefore runoff, infiltration, and land–atmosphere interactions, is unknown. The large extent of no-till crop production means that even small changes in snowmelt or land–atmosphere forcing may have regional implications.
Forest and short vegetation canopies are analogous to stubble, and their influences on the snowmelt energy balance have been the subject of substantial research. Canopies attenuate the transmission of shortwave radiation (Bewley et al. 2005; Ellis and Pomeroy 2007; Musselman et al. 2015; Pomeroy et al. 2009; Reid et al. 2014) and enhance subcanopy longwave irradiance to the snow surface (Essery et al. 2008a; Sicart et al. 2006; Webster et al. 2016). Approaches to estimate shortwave attenuation vary between simple Beer’s law methods (Sicart et al. 2003; Mahat and Tarboton 2012; Pomeroy and Dion 1996) to more complex methods that use either two-stream solutions (Mahat and Tarboton 2012), consider sky view factors estimated from hemispherical photography (Musselman et al. 2012), or implement computationally expensive ray tracing (Essery et al. 2008b; Musselman et al. 2015). Longwave radiation contributions are often estimated using sky view factors (Rowlands et al. 2002) in conjunction with observed or modeled canopy temperatures (Musselman and Pomeroy 2017; Pomeroy et al. 2009; Webster et al. 2016). The major difference between stubble and forest interactions on radiation transfer behavior relate to the relative sizes of the elements. Stubble height, unlike a forest, is on the same order as snow depth, and thus stubble will transition over the melt period from being buried to becoming fully exposed. The shortwave radiation attenuation and longwave emittance from exposed stubble is therefore dynamic. In contrast, forest canopy contributions are generally static as the bulk of the canopy is well above the surface.
In short and sparse canopies, turbulent transfer is often estimated by local gradient diffusion approaches (K theory; Wallace 1991). The K theory predicts that increased stubble exposure over melting snow will increase surface roughness, thereby increasing the ability of the snow and stubble surface to absorb momentum, leading to increasing turbulence and turbulent transfer (Prueger and Kustas 2005). In contrast, exchange specific to the snow surface below the exposed stubble, the surface of interest in this study, is a function of the stubble exposure and does not reflect the areal average increase in turbulent transfer as predicted by K theory (Bewley et al. 2010). Alternate resistance parameterizations are required to account for observations of suppressed turbulent transfer due to the decoupling of the surface from the atmosphere by stubble influencing wind velocity profiles (Mahat et al. 2013), displacing the airflow from the surface (Brun et al. 1984; Burt et al. 2005; Cutforth and McConkey 1997) and ultimately reducing wind speeds at the surface (Aase and Siddoway 1980).
In the absence of relevant previous research, the extent to which stubble exposure will attenuate shortwave radiation, enhance longwave irradiance, and modify turbulent fluxes is unclear. It is important to understand how these relative changes will manifest themselves in terms of the net snowpack energy balance over the snowmelt period. The overall objective of this study is to understand how exposed stubble modifies the snowpack energy balance. Specifically, its purpose is to 1) develop and validate a model to simulate the snowpack energy balance as a function of the exposed stubble characteristics and 2) use this model to develop a quantitative understanding of the compensatory relationships between stubble characteristics and the snowpack energy balance.
2. Stubble–snow–atmosphere snowmelt model development
a. Snowpack energy balance
The role of stubble in modifying energy transfer to the underlying snow surface is manifested through the snowpack energy balance, given as (Gray and Male 1981)
where all terms have units of watts per square meter, and is net shortwave irradiance, is net longwave irradiance, is sensible heat flux, is latent heat flux, and is the flux of energy advected to the snow through precipitation. The ground heat flux is negligible during snowmelt periods on the Canadian Prairies and is hereafter neglected (Pomeroy and Goodison 1997; Granger and Male 1978). The sum of the left-hand side of Eq. (1) is the net energy flux for the snowpack that either changes the internal energy of the snowpack or melts snow . The sign convention is for positive fluxes to be directed toward the snow. The energy balance interactions and mass fluxes for the snow–stubble–atmosphere are visualized in Fig. 1. The model developed is called the Stubble–Snow–Atmosphere Model (SSAM).
1) Shortwave radiation
Beer’s law relates light attenuation to the properties of a scattering medium and represents the transmission of shortwave radiation through a stubble canopy. The transmittance of stubble can be expressed as
where (dimensionless) is the extinction coefficient and PAI (m2 m−2) is the plant area index. The plant area index is conceptualized as the one-sided surface area of vegetation per square meter (Campbell and Norman 1998; Eagleson 2002; Shaw and Pereira 1982). Stubble may be idealized as a collection of uniform vertical cylinders; therefore, PAI is proposed to be estimated as half the surface area of a collection of vertical cylinders, which simplifies to
where r (m) is the radius of an individual stubble stalk, (m) is the height of the exposed stubble, and is the areal density of stubble stalks (number of stalks per square meter). This ignores the surface area of stubble ends, which is minor. A minimum is set to 0.001 m to avoid numerical instabilities throughout the model. Shortwave radiation is composed of diffuse (subscript d) and direct beam (subscript b) components. For direct radiation, comprises the ratio of the shadow area to surface area, which for an opaque vertical cylinder simplifies to (Eagleson 2002)
where (rad) is the solar elevation angle. This is modified to account for transmission associated with forward scattering through a simple parameterization from Goudriaan (1977):
where (dimensionless) is a scattering coefficient that is the sum of a material’s reflectance and transmittance (Goudriaan 1977; Wang 2003). Stubble stalks are opaque; therefore, it is assumed that ≈ stubble albedo [ (dimensionless)]. For diffuse radiation, where the radiation is incident upon a surface from all directions, canopies with a vertical leaf angle distribution have = 0.55 (Eagleson 2002). In the absence of observations of direct and diffuse shortwave components, the incoming above-canopy shortwave radiation (W m−2) is partitioned by calculation of the diffuse fraction of incoming shortwave radiation fd (dimensionless) with the empirically based model of Reindl et al. (1990). The shortwave radiation incident at the snow surface (W m−2) is represented as
where from Eq. (2) is calculated separately for direct and diffuse (both dimensionless) as
The bulk transmittance of the stubble (dimensionless) is given as
The net shortwave radiation at the snow surface must account for the albedo of the snow surface (dimensionless), which gives
The decreases over time as snow ages and undergoes metamorphosis and is estimated prognostically from the application of Verseghy (1991) as
where is the albedo decay coefficient (dimensionless), (0.75) is the minimum snow albedo, (s) is the interval duration, and identifies the time interval. In the event of snowfall, is refreshed to 0.90. The Verseghy (1991) snow albedo model was selected for its simplicity, and more physically detailed snow albedo models, such as Snow, Ice, and Aerosol Radiation (SNICAR; Flanner et al. 2007) or physically based snow albedo model (PBSAM; Aoki et al. 2011), could be implemented in future applications.
2) Longwave radiation
Net longwave radiation at the snow surface is the sum of incoming (W m−2) and outgoing (W m−2) longwave radiation. Exposed stubble modifies through changing the stubble surface temperature (K) and the sky view factor (dimensionless) from the perspective of the snow as (Pomeroy et al. 2009)
where (W m−2) is the incoming longwave radiation above the canopy, (dimensionless) is the emissivity of the stubble, and (5.67 × 10−8 W m−2 K−4) is the Stefan–Boltzmann constant. A novel simulation solution for the term is expressed in appendix A. The term is expressed via the Stefan–Boltzmann law and depends upon the radiometric snow-surface skin temperature (K) as
where (dimensionless) is the snow emissivity. The first term on the right represents emission from the snow, and the rightmost term represents reflectance of the incoming radiation. Emissivity for all sources, snow and stubble, is assumed to be ~1.0 to account for longwave reflections (Mahat and Tarboton 2012), which removes the reflectance term from Eq. (13).
(i) Stubble temperature
The is estimated by solving the surface energy balance of a stubble stalk, adapted from Musselman and Pomeroy’s (2017) bulk approach to simulate the surface temperature of a tree trunk. The surface energy balance of a single stubble stalk is therefore
where is net shortwave radiation, is the net longwave radiation, is sensible heat flux, and is the change in stubble stalk energy storage with respect to time. Fluxes are normalized to the stubble surface area, or volume for , for a single stalk so all units are in watts. Stubble is senescent, so latent heat flux is neglected. Calculation of uses the transmittance parameterization [Eq. (9)] as
where (dimensionless) is the albedo of stubble. The first term in the square brackets represents incoming shortwave radiation from the atmosphere, and the second term represents a first-order reflectance from the snow surface. The residual of is the portion of the incident upon the stubble and dividing by provides the mean incident shortwave radiation for each stalk. This parameterization accounts for the decrease in transmissivity with and conserves energy. The term is quantified as
where (dimensionless) is the view factor composed of other stubble from the perspective of a single stalk. The terms consider (from left to right) incoming longwave radiation from the atmosphere, snow, and adjacent stubble, respectively, and emittance from the stalk. An analytical simulation solution to quantify is shown in appendix B. The residual of is divided between the underlying snow surface and atmosphere to estimate contributions of each source. The term is estimated as
where (K) is the air temperature and (s m−1) is the resistance to sensible heat transfer between the stubble stalk and the atmosphere. The term calculates heat storage as a change from the previous interval as (Gouttevin et al. 2015)
where (J K−1) is the stubble stalk heat capacity multiplied by the temperature difference in (K) between time steps. The term is defined as
where (m3) is the stubble stalk volume, (kg m−3) is the volumetric mass density of stubble, and (J kg−1 K−1) is the specific heat capacity of the stubble. For wheat stubble is 121 kg m−3, and is 1630 J kg−1 K−1 (Ahn et al. 2009). A sensitivity analysis, not shown, varied and by several orders of magnitude and demonstrated that these parameters have a negligible influence on , due to the small volume of the individual stalks, and are hereafter assumed transferable between stubble types. The is a function of forced and free convection coefficients for sensible heat transfer between a stubble stalk and the surrounding air as
Following Monteith and Unsworth (2008), is calculated as
where (20.2 × 10−6 m2 s−1 at 10°C; Denny 1993) is the molecular diffusivity for heat in air. The Nusselt number (dimensionless) estimates the degree of turbulent transfer due to forced convection as a function of the Reynolds number (Monteith and Unsworth 2008):
where Re is the Reynolds number of the stubble stalk cylinder that quantifies the ratio of inertial forces to viscous forces in a fluid:
where (m s−1) is wind speed at the top of the stubble canopy, (kg m−1 s−1) is the viscosity of air, and the stubble stalk radius is taken as the length scale. The calculation of utilizes the logarithmic wind profile assumption to give
where (m s−1) is the wind speed at the measurement height (m), (m) is the canopy displacement height [Eq. (35)], and (m) is the aerodynamic roughness length [Eq. (34)]. Sutherland (1893) calculates as
where and are coefficients with values of 1.458 × 10−6 kg m−1 s−1 K−1/2 and 110.4 K, respectively. The term is calculated as
The free convection Nusselt number follows Monteith and Unsworth (2008) and is estimated as a function of the Grashof number (Gr):
(ii) Radiometric snow-surface temperature
The used to simulate is estimated from the energy balance of an infinitesimally small snow surface layer in equilibrium with the atmosphere and thermally decoupled from the underlying snowpack (Pomeroy et al. 2016). The internal change of temperature in a snowpack is supplemented by absorption of incoming shortwave radiation and buffered by the thermal heat capacity of the snow mass while the skin surface cools because of longwave exitance and latent heat losses from sublimation. Thus, the skin surface energy balance is
where (0.05; Pomeroy et al. 2016) is the shortwave absorption factor that represents the amount of shortwave radiation absorbed in the nearly transparent near surface snow layer; (kg m−3) is the air density; (1005 J kg−1 K−1) is the specific heat capacity of air; (s m−1) is the aerodynamic resistances for exchange between the surface the atmosphere [Eq. (32)]; (2.835 × 106 J kg−1) is the latent heat of sublimation; and and (both kg kg−1) are the specific humidities of the unsaturated air and saturated snow surface, respectively.
3) Turbulent fluxes
The total sensible and latent turbulent heat exchanges between snow, stubble, and the atmosphere are solved in a simple parallel resistance bulk transfer approach. In bulk gradient form
where (K) is the snow temperature at the turbulent exchange interface. The parameterization of with respect to an emerging stubble assumes an exponential wind profile within and a logarithmic wind profile above the stubble canopy, as (Mahat et al. 2013)
where (0.4) is the von Kármán constant, and (dimensionless; 2 for wheat; Brutsaert and Parlange 1992) is the exponential wind decay coefficient. The eddy diffusion coefficient (m2 s−1) at , which assumes is equivalent for momentum, sensible heat, and water vapor exchange (Mahat et al. 2013), is calculated as
Choudhury and Monteith (1988) provide parameterizations for and sensitive to vegetation features that were successfully applied to wind profiles observations within stubble canopies (Aiken et al. 2003) as
where is the drag coefficient of the stubble element [dimensionless; 0.5 from Aiken et al. (2003)] and is the snow-surface roughness (0.005 m; Pomeroy et al. 2016). The silhouette area index [SAI (m2 m−2)], the vertical cross section of roughness elements per unit area, is (Aiken et al. 2003)
Bulk latent energy exchange is restricted between the snow and atmosphere while sensible heat exchange occurs between the stubble, snow, and atmosphere as
The bulk sensible heat flux contribution from stubble (W m−2) is adapted from Ménard et al. (2014) as
Turbulent exchange occurs within the snow, not at the infinitesimally thin snow surface represented by , therefore represents a distinct snowpack temperature. The to calculate and comes from the relation to the snowpack internal energy that is detailed in section 2a(5). Atmospheric stability is assumed to be neutral to avoid making unrealistic and unsupported assumptions of effective surface parameters or stability corrections derived over a homogenous surface. It would be physically unrealistic for a parallel resistance scheme of snow and stubble elements to have both stable and unstable stability corrections in the same atmospheric column. Exposed stubble over a snow surface has not been the focus of detailed turbulent transfer studies, so no publications examining stability schemes over this unique surface are available to inform the model choice of stability scheme.
4) Energy advected by precipitation
The energy advected by precipitation is positive toward the snowpack in the case of rainfall and negative in the case of snowfall. A simple representation is
where Rain (kg m−2) is rainfall, (4184 kJ kg−1 K−1) is the specific heat capacity of water, is temperature of the precipitation (K), Snow (kg m−2) is snowfall, (2010 kJ kg−1 K−1) is the specific heat capacity of ice, and is the melting temperature of water (273.15 K). The may be assumed from or estimation of the hydrometeor temperature (Harder and Pomeroy 2013).
5) Internal energy change and melt energy
During snowmelt, often exhibits a diurnal pattern with a negative flux at night, from longwave emittance, and a positive flux during the day, from shortwave and enhanced turbulent terms. Shallow snowpacks have a limited capacity to moderate this variation in energy exchange, and to accurately simulate the diurnal pattern of the snowpack warming, ripening, melting, refreezing, and cooling, the internal energy of a snowpack (J) needs to be explicitly tracked. This is done in a prognostic manner as
where is initialized at as a function of initial SWE and , as
The dependence of on relates the term to the change in over time as
The partitioning of into and is a function of whether the snow is isothermal, which can be inferred from as
To avoid excessive cooling of the snowpack, due to numerical instabilities when snow is very shallow, a minimum internal energy (J) is defined as
where (K) is the minimum temperature for the preceding 24 h (Gray and Landine 1988). In any case where is less than , the associated is reset to that associated with . A full description of includes energy associated with ice, water, and vapor states. This simple model only tracks the energy of the ice portion as vapor contributions are minimal (Pomeroy and Goodison 1997), and the liquid water component is a distinct storage term in the mass balance.
6) Energy balance solution
SSAM’s solution is complicated by the interdependence of , , and . In addition, Eqs. (14) and (40) include internal energy storage terms for a stubble stalk and the snowpack, respectively. A prognostic solution is employed within a nonlinear equation solver in R (nleqslv; Hasselman 2017) that uses a Jacobian optimization scheme with a Broyden update to simultaneously solve the coupled energy balances found in Eqs. (14), (29), and (1). The solution is initialized with during a nighttime interval when temperature differences are minimal. Because of the phase change of water from solid to liquid, and must be ≤273.15 K. In situations where the solution results in and/or > 273.15 K, and/or are set to 273.15 K, SSAM is rerun, and the resulting positive term is available to ripen or melt snow.
b. Snowpack mass balance
The influence of stubble characteristics upon snow depletion is simulated with coupling of the snowpack energy balance with a simple single-layer snowpack mass balance model. The snowpack mass balance allows for simulation of both accumulation and depletion of snow as
where (kg m−2) is meltwater discharge, (kg m−2) is sublimation or deposition, (kg m−2) is snow erosion or deposition from blowing snow, and (kg m−2) is blowing snow sublimation. In cases of blowing snow, changes in SWE will correspond to changes in per the ratio , assuming is unaffected. The energy balance is coupled to the snowpack mass balance through conversion of to melt in terms of equivalent SWE as
where (334 kJ kg−1) is the latent heat of fusion. The term also represents a mass exchange of either sublimation or deposition and is put into terms of equivalent SWE as
where (2835 kJ kg−1) is the latent heat of sublimation. Only once snow becomes isothermal, , will snow begin to melt. As snow is a porous media, the initial melt increases the liquid water content [LWC (kg m−2)] of the snow rather than immediately discharging meltwater from the snowpack. The maximum liquid-water-holding capacity of snow [ (kg m−2)], the amount of water that can be held by snow without draining, defines how much snowmelt is needed to fully ripen a snowpack. This is estimated as (Essery 2015)
where is the density of water (1000 kg m−3), (m) is the snow depth, (dimensionless) is the liquid water capacity of snow (0.001–0.08; Pomeroy and Brun 2001), and (dimensionless) is snow porosity as,
where is the density of ice (917 kg m−3). Once the snowpack is ripe, any additional depletes the SWE and is discharged from the snowpack. Meltwater flow and retention is complicated by heterogeneous snow structure and is shown to have large implications for meltwater discharge from a snowpack (Leroux and Pomeroy 2017; Marsh and Woo 1984) but is ignored hereafter as snow structure is not represented in the single-layer snow model. A negative refreezes the liquid water prior to the cooling of the snowpack. The conditional algorithm used to partition into or and track changes in , LWC, and SWE is detailed in Fig. 2.
3. Data and methods
The field site near Rosthern, Saskatchewan, Canada, is representative of a no-till agricultural region on the northern Canadian Prairies, where agricultural practices control physical characteristics of the vegetation cover. The landscape has little relief and is interspersed with woodlands and wetlands. Snow depth accumulation is typically less than 0.5 m. Pomeroy et al. (1993, 1998) described the snow accumulation and melt energetics of similar environments.
To assess SSAM, snowmelt field campaigns in 2015 and 2016 collected observations to test SSAM components for a selection of stubble treatments. Site characteristics are summarized in Table 1.
1) Shortwave radiation
Direct observation of employed an array of Apogee SP110 pyranometers. Snow depth is dynamic, therefore sensors were mounted on threaded rods to allow vertical adjustment (Fig. 3). Sensors were cleaned and adjusted to the snow surface after snow accumulation or ablation events. Periods when sensors were buried were removed from the analysis. At each observation site, two pyranometers were placed, one within the stubble row and the other between the stubble rows, to account for variability in . Herein, the observations reported from each site are the average of the two sensors.
2) Radiometric snow-surface temperature
Snow-surface temperature was observed using Apogee SI-111 infrared radiometers. At each observation site, a SI-111 was fixed to a mobile platform that was shifted as needed to restrict observation to snow surfaces.
3) Stubble temperature
Stubble temperature is challenging to measure as stubble elements are very small. Two approaches were taken. First, thermocouples (30-gauge Type T) were inserted into the stalks through a small incision. At each site eight thermocouples were inserted over the vertical extent of . A challenge is that thermocouples sample the interior temperature of the stubble, which may differ from the stubble surface. The second approach involved intermittent thermography of the stubble using a FLIR T650 thermal camera. Significant challenges exist with thermography, specifically, the ability to focus on the stubble elements, the smearing of the emitted radiation over the coarse pixels, the variability and uncertainty of surface emissivity, and the environmental conditions (Muniz et al. 2014; Shea and Jamieson 2011). Depending on input uncertainty, the standard error of thermography can be up to 3°C (Muniz et al. 2014). Bias was corrected with the difference between the snow temperature observed by the FLIR and an adjacent SI-111.
4) Eddy covariance
Eddy covariance (EC) instrumentation was deployed during the 2015 observation campaign to Short15 and Tall15 treatments to observe the areal average LE and . The identical setups used LI-COR 7500A open-path infrared gas analyzers in conjunction with Campbell Scientific CSAT3 sonic anemometers. Sensor heights were 1.8 m on both sites to ensure the flux footprints remained within the stubble treatment domains while also sampling a representative span of eddy sizes. Data were logged at 20 Hz and postprocessed with default settings in Eddy Pro Software (v.6.2.0) to give 30-min average flux observations. The Mauder and Foken (2006) procedure assessed data quality, and only zero-flag data were used in this study; from 59% (Short15) to 77% (Tall15) of the time series data were retained.
5) Meteorological data
A permanent meteorological reference station adjacent to the instrumented stubble treatments observed and with a Kipp and Zonen CNR1 net radiometer, and wind direction with an RM Young 05103 Wind Monitor, precipitation with an alter-shielded Geonor TB-200 weighing gauge, and and RH with a Campbell Scientific HMP45C212. Precipitation phase was estimated by applying Harder and Pomeroy’s (2013) psychrometric approach, and snow undercatch was corrected using the correction of Smith (2009).
6) Stubble characteristics
Information on stubble characteristics—, , , , and row spacing (row)—is required to parameterize SSAM. These parameters were sampled for the canola and wheat stubble treatments during the observation campaigns. The respective treatments’ values are given in Table 1 with remaining characteristics summarized by crop type in Table 2. The exposure of over the course of melt is estimated from the difference between the measured snow-free and the measured over snowmelt from the snow surveying described in section 3b(7).
Plant area index
Independent observations of PAI relative to were obtained with a Decagon AccuPAR LP-80 Ceptometer. The LP-80 is a portable line quantum sensor, 0.84-m probe length with 80 sensors, that measures canopy photosynthetically active radiation (PAR) transmissivity. Leaf area index (LAI), analogous to PAI in this situation, is related to PAR transmissivity with an empirical extinction coefficient (Welles and Cohen 1996). The LP-80 is appropriate for PAI observations in a short discontinuous stubble canopy and provides an independent validation of the proposed PAI parameterization (Decagon Devices 2016). Vertical PAI profiles were obtained by sampling at 0.05-m intervals.
7) Snow surveys
Snow surveys provided regular observations of , measured by a snow probe, and , measured by snow coring with an ESC-30 snow tube. The snow courses for each treatment consisted of 130 and 17 observations in 2015 and 90 and 12 observations in 2016. Snow surveys were conducted on 1–3-day intervals during the melt period. Areal average SWE was estimated from multiplication of observed and mean . The sampling uncertainty of SWE estimates is quantified with bootstrapping. The and observations were resampled 10 000 times to develop robust estimates of the respective 95% confidence interval (CI). The 95% confidence for SWE () is calculated as (Steppuhn and Dyck 1974)
c. Model validation
The microscale nature of the stubble and snow environment and limitations of available instrumentation prevents direct validation of all SSAM energy balance terms. The PAI parameterization was assessed by comparing observed PAI profiles to estimates calculated with the observed stubble properties. The assessment spanned the snowmelt period from 8 to 29 March for both 2015 (at four sites) and 2016 (at two sites). The PAI for calculation of used as the difference between the radiometer height and stubble height and the observed mean for the respective stubble types. The reference station supplied the observations. SSAM was run with input data of 15-min and and 2-m , RH, and from the meteorological reference station from 8 to 29 March for the respective years to test , , areal average LE and , and SWE performance. Snow properties (SWE, , and ) were initialized from the snow survey and radiometer observations at the time of maximum observed SWE. The variable in the term and were adjusted based on relative field conditions: and in 2015, as the presence of significant ice layers led to faster decline and lower , and and in 2016, as the relatively clean snowpack led to slower decline and greater . Model performance was assessed by comparing against the mean temperature of exposed in situ stubble thermocouples and thermography estimates, against the infrared radiometer observations, and LE against the EC observations, and SWE against snow survey observations.
Model performance was assessed with the root-mean-square error (RMSE) and model bias (MB). Each test provides a different perspective on model performance: RMSE is a weighted measure of the difference between the observation and model (Legates and McCabe 1999), and MB indicates the mean over or underprediction of the model versus observations (Fang and Pomeroy 2007). All error metrics are rounded to two decimal places, so any MB values reported as 0 are actually <0.0049.
d. Model sensitivity
The overall influence stubble exposure has on the terms of the snow energy balance was explored with a sensitivity analysis of SSAM. Canola and wheat stubble, defined by Table 2 parameters, was simulated with varying from 0 to 0.5 m at 0.02-m increments. A consistent baseline for the meteorological inputs, corresponding to typical snowmelt conditions in the middle of March in this region, used = 3°C, RH = 75%, = 4 m s−1, = 500 W m−2, and = 266 W m−2. While keeping these variables constant, each meteorological input was varied across a range consistent with daytime snowmelt conditions on the Canadian Prairies (Table 3). A sun angle of 37° and an value of 0.75 were specified. Steady-state conditions, no change in with time, and a ripe snowpack, = 0°C, are assumed.
4. Results and discussion
a. Model performance
1) PAI parameterization performance
The PAI parameterization of SSAM controls shortwave radiation interception and turbulent exchange processes. Observations of PAI, as they vary with , are plotted in Fig. 4a and demonstrate the clear differences between a sparse canola stubble and a dense wheat stubble. The PAI parameterization [Eq. (3)] is plotted against observations in Fig. 4b and demonstrates that the proposed model accounts for differences in stubble characteristics. The model successfully estimates PAI with low RMSE and MB.
2) Shortwave radiation performance
The predicted struggles to account for the randomness of the discontinuous stubble that causes the surface to vary between fully exposed and shaded. With only two sensors at each site, the observed is more variable than the areal average behavior on an hourly interval. The effect of stubble gaps and limited sensors is apparent in the relatively large values of RMSE and the scatter of modeled and observed hourly subcanopy radiation in Fig. 5. Cumulative radiation is of greater concern for snowmelt modeling than instantaneous values, and the model and observations are much more appropriate to capture this dynamic. There is little difference in cumulative values as demonstrated with low MB values, ±0.03 W m−2, and little difference between cumulative observed and modeled subcanopy radiation (Fig. 5).
3) Stubble temperature performance
The comparison of estimated to observations reveals systematic differences as SSAM consistently overestimates relative to the observations from thermocouples during peak sun angles (Fig. 6). This is an expected consequence of the core stalk temperature being buffered by the stubble heat capacity and conductivity. The higher FLIR observations, relative to thermocouples, are comparable to SSAM values during the daytime, confirming that surface temperatures exceed the core stalk thermocouple measurements. Considering thermography uncertainty, SSAM is performing well; RMSE of thermography estimates are near the expected FLIR uncertainty of 3°C (not shown).
4) Radiometric snow-surface temperature performance
The estimated showed similar accuracies to those reported by Pomeroy et al. (2016) for homogenous snow cover with RMSEs between 1.4° and 1.8°C (Fig. 7). The temporal relationship between observed and modeled is strong, though SSAM tends to overestimate the nighttime cooling on cold calm nights.
5) Turbulent fluxes
Assessment of the turbulent fluxes is limited to Tall15 and Short15 sites. Therefore, transferability of the SSAM resistance scheme to canola is untested. Over wheat treatments, LE showed excellent temporal agreement with limited scatter (Fig. 8) and low errors while (Fig. 9) had weaker performance. The overprediction of relates to deficiencies in process representation and differences between what SSAM estimates and EC observations represent. First, the parallel resistance conceptualization does not allow turbulent exchange between snow and stubble that would moderate and ; this will lead to overprediction of positive and negative . In multiple source models, is often calculated with a serial resistance scheme that couples source temperatures through calculation of a canopy air temperature (Ménard et al. 2014; Norman et al. 1995; Best et al. 2011). Unfortunately, the physical process understanding of turbulent exchange between the dynamic exposure of short and sparse stubble elements and snow during melt is insufficient to represent conceptually as a serial resistance scheme. Turbulent exchanges in open environments are intermittent (Helgason and Pomeroy 2012a,b), and there is no consistent establishment of a distinct within-canopy air temperature needed for a serial resistance scheme to be valid. Second, there is a mismatch between the surface characteristics at the point scale of SSAM and the EC flux footprint. The EC observations reflect a flux footprint where the combined spatial variability of the stubble and snow depth will lead to a range of that differs from the model-simulated uniform . Snow cover and heterogeneity in the flux footprint are dynamic due to the variability of wind direction in each interval; therefore, noise, rather than a systematic bias, is expected and observed in scatterplots of Figs. 8 and 9. Despite differences in assumptions and scales, the representation of areal average and LE does capture the temporal behavior and magnitude of the observations. The errors are reasonable with respect to the literature as all turbulent transfer parameterizations are challenged by issues of energy balance closure and heterogeneous flux footprints (Andreas 2002; Helgason and Pomeroy 2012b; Marks et al. 2008).
6) Snow water equivalent and exposed stubble height
The SWE depletion and emergence simulations performed well with respect to the areal average SWE and observations from snow surveys (Fig. 10). Both melt seasons were characterized by an initial rapid melt in early March that was suspended by cooler weather and a moderate snowfall event, with final melt occurring rapidly at the end of March. The main difference was that 2015 had greater premelt SWE. The modeled SWE depletion and for both seasons capture these dynamics well. The model tends to underestimate through RMSE as <5.6 cm. In contrast, for SWE the model overestimated the initial melt for Wheat16 treatments and underestimated the initial melt for Tall15. Wheat16 and Canola16 simulations captured the accumulation of the midmelt snowfall event well while Short15 and Tall15 overestimated SWE after the midmelt snowfall. The timing and final melt-out sequence was captured for all sites. The SWE RMSEs were greater in 2015 than 2016 but were reasonable at less than 14 mm. These results demonstrate that SSAM’s energy balance approach, with respect to independent observations, successfully represented the melting snowpack energy balance and can predict realistic SWE depletion and emergence. A control run is also presented that demonstrates limited differences in SWE depletion with and without stubble despite an obvious difference in . It is important to reiterate that no calibration was used to optimize model performance, which provides confidence in the model process representations.
7) Validation summary
A challenge of validating SSAM is that it is extremely difficult to obtain direct observations of the processes represented in the model, due to the small-scale nature of the stubble elements and their dynamic emergence from snow during melt. The successful representation of PAI, , , , , LE, and SWE, observations of state variables, or fluxes, provides independent validation of each energy balance term. This gives confidence that SSAM can successfully reproduce the complex and dynamic interactions evident between the stubble, snow, and atmosphere.
b. Snow energy balance compensation
The sensitivity analysis of SSAM with respect to variations in stubble properties and meteorological inputs articulates the nonlinear interactions that lead to energy balance compensation. Generally, there is a limited change in with respect to increased (Fig. 11). Increasing increases , as there is more warm stubble in the view factor of the snow surface, while decreases due to the reflection and absorption of in the stubble. Turbulent fluxes, and , show a complex response to through its relationship with . As increases, the above-canopy exchange increases due to increases, and the within-canopy exchange decreases, due to increases. The surface exchange coefficient, in Fig. 12, initially decreases as there is a larger decrease of within-canopy exchange than above-canopy exchange as stubble emerges from the snow. An inflection point occurs near 0.04-m for both canola and wheat stubble when the rate of change of the within canopy exchange decreases. Greater exchange within the sparse canola canopy leads to an increase in exchange with increased after the inflection point. The denser wheat stubble exchange has limited change with below the inflection point. The influence of this relationship with exposed stubble is evident in all simulations of .
The behavior of the response to stubble exposure with respect to variations in , RH, and is a consequence of the turbulent snow–atmosphere interactions (Fig. 13). As is expected for °C, increases in °C will increase and as temperature and humidity gradients increase. A positive relationship between RH and is also demonstrated and, with the conditions simulated, two behaviors are expressed. With RH = 25% the large negative forces a negative that will initiate refreezing and possible cooling of the snowpack rather than melt. The term is quite sensitive to RH, and the values for RH > 25% show large increases that translate into large increases in . There is no response of to RH as RH cannot influence the temperature of a wet isothermal snowpack. Turbulent fluxes are directly related to and show a clear positive relationship with , , and . The denser wheat stubble, relative to the more open canola canopy, shows less response in , , and to variations in , , RH, and .
Variations of and have large impacts upon that are functions of and , respectively (Fig. 14). The is that is transmitted to and absorbed by the snow; therefore, increased will increase as and are independent terms. There is a small increase in with respect to increased as shortwave radiation is absorbed by the stubble and reemitted as longwave radiation. There is minimal response, not shown, in the and terms to variation. There is an inverse relationship between and ; as increases, decreases, which shifts the source of incoming longwave at the snow surface from the sky to the stubble. When the longwave flux from the stubble is less than from the sky, in the case of a high , will decrease with increasing and vice versa. The response of and to and is greater for wheat, relative to canola, as there are greater decreases in and with .
The cumulative difference in the net energy balance terms over the course of snowmelt will depend upon the dynamic and interacting response of meteorological conditions and stubble exposure. From the sensitivity analysis it is expected that canola or tall stubble will have greater , and therefore melt rates, than wheat or short stubble. The differences are subtle because of the complex energy balance interactions that act to compensate each other, which is very clear when comparing the SWE simulation with and without stubble in Fig. 10. To describe fully the melt patterns in this region, differences in snow accumulation due to stubble differences need to be considered.
The outcome of this work has two main implications. The first is that the lack of representation of stubble emergence in regions of seasonal snow cover and no-till agriculture is a clear deficiency in current land–atmosphere models, which affects their ability to represent snowmelt processes in agricultural regions. To improve understanding of land–atmosphere feedbacks, land surface schemes need to include the dynamics of stubble emergence. From a turbulent transfer perspective, the dynamic change in with respect to , a decline of 15% and 11% from = 0 m to = 0.04 m followed by a 5% and 21% increase from = 0.04 m to = 0.5 m for wheat and canola stubble, respectively, are unmodeled dynamics in current land surface schemes. The strength of the land–atmosphere coupling has consequences for the radiative terms of the energy balance as changes with snowmelt. The exchange between the snow and stubble surface with the atmosphere demonstrates large energy balance differences with and without stubble (Fig. 15). Exposed stubble increases shortwave absorption, which translates into greater longwave emittance and sensible heat transfer from the surface to the atmosphere. Wheat stubble has a greater stubble surface area and greater turbulent fluxes than canola, such that canola stubble is an energy sink that cools the atmosphere and wheat stubble is an energy source that warms the atmosphere. Since land–atmosphere models provide the lower-boundary conditions for numerical weather prediction and climate models, the predictive capacity of weather and climate from short to long time scales might be enhanced with these additional physics. The second is that the clear increase in premelt SWE with increasing stubble height (Pomeroy and Gray 1995) may not necessarily translate into greater frozen soil infiltration. Greater melt rates associated with this greater may lead to lower infiltration efficiency and greater runoff efficiency. This potential compensatory relationship requires further research to quantify the hydrological implications of stubble management.
Quantification of the snow energy balance response to stubble exposure improves the understanding of land–atmosphere interactions and the role of stubble management upon snowmelt processes in semiarid cold agricultural production regions. The proposed SSAM model represents the snow energy balance underlying exposed stubble and is validated successfully against subcanopy , , , , LE, and SWE observations. A sensitivity analysis of SSAM with respect to stubble exposure and meteorological forcing shows that stubble exposure demonstrates complex energy balance interactions. Generally, stubble will decrease , , and and increase , resulting in a subtle increase in . The change in with may have snowmelt infiltration and runoff implications but requires further study that addresses stubble influence upon accumulation processes. In contrast, the ability for stubble to influence individual energy balance terms will have land–atmosphere modeling implications and provides a compelling argument to include these additional physics into land–atmosphere models. The ability to quantify these stubble–snow–atmosphere energy balance interactions in small and large-scale models is important to understanding the impact of stubble management on snowmelt processes, land–atmosphere interactions, and hydrology in cold agriculture regions.
Funding comes from the Natural Sciences and Engineering Research Council of Canada through Discovery Grants, Research Tools and Instruments, the Changing Cold Regions Network, and the Canada Research Chairs programme. Field and technical assistance from Bruce Johnson, Chris Marsh, Kevin Shook, and Michael Schirmer, and the flexibility of Nathan Janzen and Robert Regehr, the farmers of the study area, are gratefully acknowledged.
Snow-Surface Sky View Factor Parameterization
The sky view factor is calculated from the perspective of a snow surface underlying exposed stubble. The view factor from the perspective of a surface plane composed of a perpendicular cylinder is solved with a solution provided by Sparrow et al. (1962). With assumed to be the residual of the cylinder view factor, the solution for a single cylinder is
where (m) is the radius of the cylinder, (m) is the distance of the surface of interest to the origin of the cylinder, and (m) is the height of the cylinder. The areal average of for a collection of stubble stalks is simulated as follows:
Simulate the stubble locations randomly per , row, and stubble row width.
Identify a representative sample area within the simulation domain.
Remove stalks located behind other stalks and calculate distance between remaining stalks and each selected coordinate of the sample area.
Sum individual view factors as per Eq. (A1) at each sample coordinate.
Calculate mean from values at each sample coordinate.
The simulation domain used in this analysis is 4 m × 4 m and the sample area (1 m × 0.5 row spacing) is in the middle of the domain. The sample coordinates are spaced every 0.02 m in both and dimensions.
Stubble-to-Stubble View Factor Parameterization
The stubble-to-stubble view factor is the portion of the view factor from the perspective of a single stubble stalk composed of other stubble. For two parallel cylinders a solution is given by Juul (1982) as
where (m) and (m) are the radii of parallel cylinders 1 and 2, (m) is the distance between the origins of the cylinders, and (m) is the height of the stubble. Simulation for case of mean is as follows:
Simulate the stubble locations randomly per , row, and stubble row width.
Identify stalks that will provide a representative estimate of .
Remove stalks located behind other stalks and calculate distance between each remaining stubble stalk and each selected stalk.
Sum individual view factors per Eq. (B1) to give at each stalk.
Calculate mean for selected stalks.
The simulation domain used is 4 m × 4 m, and the stubble stalks of interest are a 1-m row in the middle of the domain.