## Abstract

Rangelands are often characterized by a patchy mosaic of vegetation types, making measurement and modeling of surface energy fluxes particularly challenging. The purpose of this study was to evaluate surface energy fluxes measured using three eddy covariance systems above and within two rangeland vegetation sites and use the data to improve simulations of turbulent energy fluxes in a multilayer plant canopy model: the Simultaneous Heat and Water (SHAW) model. Model modifications included adjustment of the wind profile roughness parameters for sparse canopies, extending the currently used K-theory approach to include influence of the roughness sublayer and stability functions within the canopy, and in a separate version of the model, introducing Lagrangian far-field turbulent transfer equations (L theory) in lieu of the K-theory approach. There was relatively little difference in simulated energy fluxes for the aspen canopy using L-theory versus K-theory turbulent transfer equations, but L theory tracked canopy air temperature profiles better during the growing season. Upward sensible heat flux was observed above aspen trees, within the aspen understory, and above sagebrush throughout the active snowmelt season. Model simulations confirmed the observed upward sensible flux during snowmelt was due to solar heating of the aspen limbs and sagebrush. Thus, the eddy covariance (EC) systems were unable to properly quantify fluxes at the snow surface when vegetation was present. Good agreement between measured and modeled energy fluxes suggest that they can be measured and simulated reliably in these complex environments, but care must be used in the interpretation of the results.

## 1. Introduction

Rangeland ecosystems are often characterized by a patchy mosaic of vegetation types, making measurement and modeling of surface energy fluxes in these environments particularly challenging. Understanding the role of different vegetation communities in modulating energy, water, and carbon fluxes is critical to quantifying the energy flux, carbon storage, and water balances of these ecosystems.

Eddy covariance (EC) systems have gained popularity as a means to measure the surface energy, water, and carbon fluxes. However, landscape variability and complex terrain often necessitate deployment of EC systems at less than ideal sites in order to quantify fluxes in patchy ecosystems. In the past decade, application of EC to measure fluxes of heat, water, and carbon has been extended to more complex sites (Baldocchi et al. 2000a; Turnipseed et al. 2002, 2003; Pomeroy et al. 2003; Kosugi et al. 2007; Hiller et al. 2008; Marks et al. 2008; Reba et al. 2009).

EC applications to understory vegetation are subject to much caution because the underlying hypotheses are generally not valid in conditions typically present in forest understories. Despite these limitations, several studies have found the method capable of accurately quantifying fluxes from the forest understory (Baldocchi et al. 2000b; Marks et al. 2008; Reba et al. 2009) and have partitioned fluxes within the canopy by having separate EC systems above and below the canopy (Blanken et al. 1997; Constantin et al. 1999; Blanken et al. 2001; Wilson et al. 2000; Scott et al. 2003; Jarosz et al. 2008).

EC systems deployed above and within plant canopies can be particularly useful in developing and refining energy balance models. Leuning (2000) used EC systems and measured temperature, humidity, and carbon profiles within a rice canopy to evaluate atmospheric stability corrections and compare with a multilayer canopy model. Baldocchi et al. (2000b) evaluated fluxes of mass and energy near the forest floor pine stands with eddy covariance measurements and a soil–plant–atmosphere exchange model. Wohlfahrt (2004) used EC measurements to compare modeled fluxes and concentrations within and above a mountain meadow canopy for different Lagrangian models and parameterization schemes. Dufrêne et al. (2005) compared a multilayer model of net ecosystem exchange for a beech forest with EC measured fluxes. Marks et al. (2008) used EC measurements beneath a pine canopy to evaluate the turbulent flux calculations of an energy balance snowmelt model. Haverd et al. (2009) used a combination of above-canopy EC measurements, within-canopy chamber measurements, and a soil–vegetation–atmosphere transfer (SVAT) model to determine vertical profiles of the Lagrangian time scale that optimized agreement between measured and modeled temperature and vapor and carbon dioxide profiles.

Combining flux measurements with model simulations helps define knowledge gaps in current understanding of flux processes, enhances our understanding and interpretation of the observed fluxes, and assists in gauging the reliability of measured and modeled fluxes. Objectives of this study were to evaluate observed surface energy fluxes measured above sagebrush, aspen, and the understory beneath an aspen canopy using eddy covariance systems and to use the data to improve simulations of turbulent energy fluxes in a multilayer plant canopy model. The Simultaneous Heat and Water (SHAW) model was used. Previously, the model has been based on K-theory turbulent transfer within the plant canopy. Lagrangian turbulent transfer was introduced into the model used in this study, as well as provisions for computing wind profile parameters for sparse plant canopies.

## 2. Materials and methods

The study area is the Reynolds Mountain East (RME) catchment located in the southwestern portion of the Reynolds Creek Experimental Watershed (RCEW) operated by the U.S. Department of Agriculture (USDA) Agricultural Research Service, Northwest Watershed Research Center. RCEW is located in the Owyhee Mountains of southwestern Idaho. RME is a 39.0-ha headwater catchment that ranges in elevation from 2024 to 2139 m MSL (Pierson et al. 2001). The catchment is dominated by low and mountain big sagebrush (*Artemesia arbuscula* and *Artemesia tridentada vaseyana*) and rocky ground covering 69% of the catchment; patches of aspen (*Populus tremuloides*) and willow (*Salix* spp.) cover 26% and fir (*Abies* spp.) covers the remaining 5% of the catchment area.

Three EC systems described by Flerchinger et al. (2010) were established to monitor fluxes across the RME catchment as part of a long-term study to characterize the hydrology of this mountainous headwater catchment. One system was located over a wind-exposed sagebrush site, another within the understory of the aspen grove, and a third situated above the aspen canopy. Observations from these three systems overlapped during 2007, which is the focus of the current study.

Vegetation at the sagebrush site consists of about half sagebrush with the remainder consisting of equal amounts of native grasses and forbs. Sagebrush at the site is approximately 60 cm in height with a leaf area index (LAI) of 0.77 based on point frame measurements. The site is a gently rolling ridge top with slope varying from 1% to 3%. The aspen site consists of an aspen grove with an understory of grasses and forbs. Average tree height was 9.5 m and the maximum tree height was 15 m (Flerchinger et al. 2010). Stem area index (SAI) of the trunks and limbs prior to the growing season was 0.45; maximum LAI of the aspen measured during the growing season was 1.35 in August. Understory vegetation measurements indicated an LAI of approximately 1.2 for the grasses and forbs, which were about 75 cm in height.

EC systems used to measure turbulent fluxes consisted of a three-dimensional sonic anemometer (model CSAT3, Campbell Scientific, Inc., Logan, Utah) and an open path infrared gas analyzer (IRGA; model LI-7500, LI-COR, Inc., Lincoln, Nebraska) sampled at 10 Hz. The EC systems were located at 5 m above the ground surface at the sagebrush site and 4.5 and 19.25 m at the aspen site. Shortwave and longwave radiation, air temperature, and humidity were collected every 30 min using a four-component net radiometer (CNR-1, Kipp & Zonen, Delft, The Netherlands), and a temperature/humidity probe (HMP45C, Viasala, Helsinki, Finland). Ground heat flux was measured with up to six heat flux sensors (HFP01, Hukseflux, The Netherlands) installed 0.08 m deep within the soil and three sets of self-averaging thermocouples installed at 0.02 and 0.06 m deep. A single set of soil heat flux sensors were shared by the understory and above aspen sites. Soil moisture used to compute volumetric heat capacity of the soil was measured hourly at 0.03 m using Hydra Probe II soil moisture sensors. Wind speed, humidity, and air temperature were measured using 2D sonic anemometers and temperature/humidity probes at heights of 3, 9, and 15 m within the aspen canopy. Dual-gauge precipitation systems especially designed for the windy and snow-dominated conditions prevalent in the area were used to measure precipitation (Hanson 1989; Hanson et al. 1999, 2004).

Processing of the eddy covariance data is described by Flerchinger et al. (2010). Postprocessing of the 30-min EC data consisted of sonic temperature correction (Schotanus et al. 1983), density correction (Webb et al. 1980), and coordinate rotation (Kaimal and Finnigan 1994). Soil heat flux measured at 0.08 m was corrected for heat storage above the heat flux plates. Heat stored within the canopy was computed empirically (Arain et al. 2003; Wu et al. 2007) based on the change in temperature over the 30-min time period (Blanken et al. 1997, 1998). Quality of the EC data assessed by energy balance closure was limited to periods without snow cover because snow temperatures and energy stored within the snowpack were not available.

EC systems are well known for their inability to close the energy balance (Twine et al. 2000; Wilson et al. 2002). Flerchinger et al. (2010) assessed the energy balance closure by regressing the turbulent fluxes *H + LE* versus net radiation (*R _{n}*), canopy storage (

*S*), and ground heat flux (

*G*) during periods when snow was not present. They reported a slope of 0.84 for the sagebrush and 0.74 above the aspen. However, a slope of 0.38 was reported for the aspen understory site because of gaps in the canopy that exposed the net radiometer in the understory to direct radiation for extended periods during midafternoon hours. Upon excluding these hours from the energy balance closure analysis, the slope of the regression line improved to 0.70.

Turbulent fluxes during the growing season were adjusted to force energy balance closure while maintaining the Bowen ratio (Twine et al. 2000). This was problematic when the Bowen ratio approached −1.0. Therefore, whenever the magnitude of *H + LE* was greater than the error in the energy balance, *H* and *LE* were adjusted equally to compensate for the error and to force energy balance closure. Unrepresentative net solar radiation measurements during the afternoon hours within the understory also created problems. The feasibility of closing the energy balance by substituting the understory solar radiation measurements with simulated values was explored. Forcing energy balance closure was not possible when snow was present.

Fetch is not a concern for most wind directions at the sagebrush site; however, the aspen sites are limited to approximately 160 m of fetch for acceptable wind directions. Wind direction had very little effect on energy balance closure for the aspen understory. The above aspen tower was located near the northeast edge of the aspen grove to maximize fetch. Therefore, periods when the wind did not originate from the direction of the aspen (170°–290° from north) were removed from analysis for the above aspen site. Analysis of the sagebrush site found that the slope of the energy balance closure for north winds was significantly different from the slopes for other wind directions (Flerchinger et al. 2010); however, data for the sagebrush site was not screened for wind direction, as energy balance closure was reasonable for all directions.

Footprint analysis conducted by Flerchinger et al. (2010) for the three sites indicated that the region exerting the maximum influence on measured fluxes under neutral conditions was 54 m from the EC system for the sagebrush site, 46 m for the aspen understory, and 36 m for the above aspen site. Approximately 65%–75% of the flux upwind of the aspen towers originated within the 160-m fetch.

Flux data were filtered for spikes, instrument malfunctions, and out-of-range signals. Two separate study periods representing the winter season (days 55–80) and growing season (days 150–225) were selected for analysis. These periods were selected because they had generally complete EC and meteorological data. Even so, approximately 22% of the hourly flux values for these periods were rejected for the sagebrush and aspen understory, and approximately 52% for the above aspen site because of either data quality or unfavorable wind direction. Therefore, the simulated and measured diurnal response of the surface energy fluxes averaged over 5- or 10-day periods, computed from composite hourly averages, were used to compare with simulated and measured fluxes.

## 3. Model description

Version 2.5 of the SHAW model (Flerchinger et al. 2009b) was used in the present study. The SHAW model has been tested and applied extensively over a range of vegetation types in semiarid and arid environments, particularly in the surrounding RCEW (Flerchinger et al. 1996, 1998), and Link et al. (2004) previously validated the model for fir forest canopies. The model simulates the surface energy balance, evapotranspiration, and fluxes within a multispecies plant canopy using detailed physics of heat and water transfer through the soil–plant–atmosphere continuum, making it ideal for use in this study.

Solar and longwave radiation exchange between canopy layers, residue layers, and the snow or soil surface are computed by considering downward direct and upward and downward diffuse radiation being transmitted, reflected, and adsorbed by each layer. Based on the analysis of Flerchinger et al. (2009a), the approach of Dilley and O’Brien (1998) is used for clear-sky estimates of incoming longwave radiation and adjusted for cloud cover estimated from input solar radiation. Flerchinger et al. (2009b) describe improvements to the model for within-canopy radiation transfer. The proportion of net solar and longwave radiation balance for each canopy layer absorbed by each vegetation type is computed based on albedo and emissivity of the vegetation, weighted by the extinction coefficient and LAI of each vegetation type within the layer. Nontranspiring plant parts (stems, trunks, etc.) are simulated as a separate vegetation type with their own properties, including an input stem area index.

The model accumulates snow depending on air temperature at the time of precipitation. Energy and mass transfer calculations for snow within the SHAW model are patterned after the point energy and mass balance model developed by Anderson (1976). Albedo is estimated from grain size, which in turn is estimated from snow density. Snow density is adjusted for compaction, metamorphosis, and settling processes. Vegetation characteristics (height, leaf area, etc.) are adjusted for that portion protruding above the snowpack; the portion of vegetation covered by snow is not simulated by the model.

### a. SHAW model version 2.5

Turbulent flux above canopy is computed in version 2.5 of the model based on the Monin–Obukov stability theory. Sensible heat is computed from

where *ρ _{a}*,

*c*, and

_{a}*T*are the density (kg m

_{a}^{−3}), specific heat (J kg

^{−1}C

^{−1}), and temperature (°C) of air at the measurement reference height

*z*

_{ref};

*T*is the temperature (°C) of the exchange surface; and

_{s}*r*is the resistance to surface heat transfer (s m

_{H}^{−1}) corrected for atmospheric stability. The exchange surface can be either the top of the canopy, the residue layer, the snow surface, or the soil surface. Latent heat flux is given by

where *λ* is latent heat of vaporization (J kg^{−1}), *E* is flux of water vapor from the atmosphere to the exchange surface (kg m^{−3}), *ρ*_{υs} and *ρ*_{υa} are vapor density (kg m^{−3}) of the exchange surface and at the reference height *z*_{ref}, and the resistance value for vapor transfer, *r*_{υ}, is taken to be equal to *r _{H}*. Resistance to turbulent heat transfer,

*r*, is computed from

_{H}where *u _{*}* is friction velocity, taken as

*k* is von Kármán constant; *d* is the zero-plane displacement; *u*_{ref} is wind speed at reference height; *z _{H}* and

*z*are the surface roughness parameters for the temperature and momentum profiles; and

_{m}*ψ*and

_{H}*ψ*are the diabatic correction factors for heat and momentum, computed as a function of atmospheric stability. If a plant canopy is present, the surface roughness parameter for the momentum profile,

_{m}*z*, is assumed 0.13

_{m}*h*and the zero-plane displacement,

_{c}*d*, is 0.77

*h*, where

_{c}*h*is plant canopy height. Surface roughness parameter for the temperature profile,

_{c}*z*, is assumed to be 0.2

_{H}*z*(Campbell 1997).

_{m}The SHAW model uses K theory for within-canopy turbulent transfer. Transfer of heat and vapor within the canopy is dependent upon location within the canopy according to the following expression (Flerchinger et al. 1998):

where *K* is the dispersion coefficient for turbulent heat and mass transfer within the canopy (m^{2} s^{−1}), and *ψ _{H}* is the diabatic correction factor dependent on stability computed from

*H*above the canopy. Assuming

*K*constant for

*z < d*and relating

*ψ*within the canopy to stability above the canopy are recognized limitations within this version of the model.

_{H}Although the SHAW model can simulate energy and mass fluxes quite well in a variety of situations, it is known that the SHAW model tends to melt snow prematurely in forested canopies (Link et al. 2004; Smith 2011). This has been attributed to difficulties in computing turbulent fluxes at the snow surface beneath a forest canopy. Therefore, modifications to the turbulent transfer routines were explored herein.

### b. Model modifications

SHAW version 2.5 had no provision for sparse canopies, but these are often present on rangelands and in developing crops. Estimation of effective zero-plane displacement, *d _{e}*, and roughness,

*z*, for sparse canopies was adopted from Zeng and Wang (2007):

_{me}Here, *d _{g}* and

*z*are the displacement height and surface roughness beneath the plant canopy;

_{mg}*z*is input to the model while

_{mg}*d*is set equal to snow depth or taken as zero for a residue or soil surface;

_{g}*V*is a function leaf area index:

where *L _{T}* is total leaf area index and

*L*

_{cr}is a critical value assumed for full canopy cover, taken as 2.0. However, as subsequently discussed in the results, the following linear interpolation for

*z*similar to the form of Eq. (6) seemed to work better than the logarithmic interpolation in Eq. (7):

_{me}Wind speed within the canopy was estimated using the approach of Nikolov and Zeller (2003) that assumes an exponential decay computed from cumulative leaf area index; Nikolov and Zeller (2003) based their equations on analysis and data presented by Massman (1987, 1997).

Two approaches were explored to improve model turbulent transfer above and within the canopy. Extension of the K-theory approach currently used in the model was explored by adopting approaches described by Oleson et al. (2010) and Harman and Finnigan (2008). The second approach was to adopt the Lagrangian far-field dispersion (Raupach 1989) with the atmospheric stability corrections given by Leuning (2000), or L theory. A limitation of K theory is that it requires local gradients to be in equilibrium with fluxes, which is not always the case in plant canopies. Raupach (1989) developed localized near-field (LNF) theory as an alternative to K theory, recognizing that scalar concentrations at any point are the result of advection from near-field sources and diffusion from sources farther away. While implementing the full LNF theory is beyond the scope of this study, we recognize that accuracy of simulated concentrations within the plant canopy may be compromised by using only the Lagrangian far-field dispersion.

#### 1) Extension of K theory

Harman and Finnigan (2008) describe flux-gradient theory within a plant canopy and the roughness sublayer. A generalized expression for friction velocity including the effects of the roughness sublayer is

where *L* is the Obukov roughness length. Herein, *z _{me}* and

*d*are substituted for surface roughness and displacement height in the above expression. The influence of the roughness sublayer on momentum transfer, , is quantified by

_{e}Harman and Finnigan (2008) give the following expression:

where *l _{m}* is the mixing length for momentum,

*β*=

*u**/

*u*,

_{h}*u*is wind speed at the top of the canopy, and

_{h}*c*and

_{m}*c*

_{2m}are proportionality constants whose values are computed from stability and canopy characteristics (Harman and Finnigan 2008). To simplify,

*β*is taken herein as 0.3, which is a typical value under neutral conditions. Expressions analogous to Eqs. (11) and (12) can be written for scalar transfer (i.e., heat and water vapor). From Oleson et al. (2010), it follows that the resistance to turbulent heat transfer can be expressed as

In lieu of Eq. (5), transfer within the canopy was taken as

Wind speed at the top of canopy was computed from Eq. (10). Recognizing that a logarithmic wind profile typically exists in the bottom 20%–25% of the canopy (Norman 1979), flux between the bottom canopy layer and the snow or soil surface was computed from Eqs. (1) through (4); surface roughness parameters for the snow or soil surface and wind speed in the bottom canopy layer were used to compute friction velocity, stability, and sensible heat at the snow or soil surface.

#### 2) Lagrangian far-field theory

The Lagrangian far-field dispersion coefficient (m^{2} s^{−1}) is given by

where *σ _{w}*(

*ζ*) and

*τ*(

_{L}*ζ*) are the standard deviation of the vertical velocity and the Lagrangian time scale as functions of stability. Leuning (2000) presented equations to ensure a smooth transition for

*τ*and

_{L}*σ*above and through the canopy. The equations and parameters provided by Leuning (2000) were used for

_{w}*z*> 0.8

*h*. Within the canopy, the Lagrangian time scale can be approximated by

_{c}and *σ _{w}* can be approximated by

Here, *u _{*g}* is friction velocity computed near the ground surface, using wind speed within the canopy and surface roughness of the snow or soil surface. The stability functions

*ϕ*and

_{H}*ϕ*were limited to the range of −2 to 1 and are computed as

_{w}and

The stability parameter applicable within the canopy was computed as

where *g* is acceleration due to gravity. Near the surface (*z* < 0.25*h _{c}*), the gradient Richardson number was used:

Here *T*_{c,i}, *u*_{c,i}, and *z _{i}* are temperature wind speed and height of canopy layer

*i*, and

*u*is wind speed at the ground surface. Leuning (2000) used the relation between

_{g}*ζ*and

*R*presented by Kaimal and Finnigan (1994):

_{i}The value of *R _{i}* is limited to 0.175 to minimize the high sensitivity of ζ on

*R*as

_{i}*R*approaches 0.2.

_{i}### c. Model application

Versions of the model were initialized with measured soil temperature and soil water profiles on day 51 for the winter simulation and on day 150 for the growing season. The aspen site was run with 8 soil layers to a depth of 1.0 m and the sagebrush site with 11 layers to a depth of 1.8 m. Free drainage by gravity was assumed for the lower soil water boundary of the sagebrush site; observed soil water content at 1.0 m was used for the aspen site because of its proximity to shallow groundwater. Simulations were driven with observed air temperature, wind speed, humidity, precipitation, and solar radiation measured above the aspen canopy and the sagebrush site. The snowpack was initialized based on snow depth and density measurements taken on day 51. Surface roughness values were set to typical values of 0.15 and 1.0 cm for the snow and soil surface, respectively. Plant parameters for stomatal resistance, root resistance, and leaf resistance of the aspen, sagebrush, and aspen understory were taken from Flerchinger et al. (1996). Nodes within the aspen canopy were specified at heights of 0.3, 3, 6, 9, and 15 m; simulated turbulent fluxes between the 3- and 6-m heights were compared to fluxes measured at the 4.5-m aspen understory EC system. Node spacing within the sagebrush canopy was selected by the model and limited to 0.5 LAI. LAI of the aspen, aspen understory, and sagebrush were assumed uniformly distributed with height. Root-mean-square deviation (RMSD) and mean bias error (MBE) were used to compare model results with observed surface fluxes and the aspen canopy temperatures profiles.

## 4. Results

### a. Energy balance closure of the understory

Gaps in the canopy exposed the net radiometer in the understory to direct radiation for extended periods during the midafternoon hours, resulting in poor energy balance closure for the understory (slope = 0.38 and *r*^{2} = 0.76) reported by Flerchinger et al. (2010). Excluding these hours resulted in a slope of 0.70 with an *r*^{2} of 0.80. The SHAW model has been shown to simulate within-canopy radiation fluxes quite well (Flerchinger et al. 2009b), so simulated net radiation within the aspen understory was substituted for measured values to obtain values more representative of average conditions within the understory. Doing so resulted in a slope for the energy balance closure of 0.67 and an *r*^{2} of 0.89. Although the slope is slightly less than when the midafternoon hours were excluded, simulated net radiation yielded a somewhat tighter correlation in the energy balance closure. Thus, simulated net radiation was used to force energy balance closure for the understory during the growing season. Further corroboration of the understory heat fluxes was explored by comparing directly with sensible heat fluxes computed between the 3- and 9-m heights using the dispersion coefficient in Eq. (14), the gradient Richardson number, and observed canopy temperature and wind speed. Correlation between these gradient-computed fluxes and the EC fluxes was 0.83 with a slope of 0.76. Thus, fluxes measured in the understory were deemed quite reasonable.

### b. Interpolation approach for z_{me}

Adjustment of the zero-plane displacement for the sparse leafless canopy (SAI = 0.45) based on Eq. (6) resulted in a decrease of *d _{e}* from 11.5 to 4.9 m. Logarithmic interpolation of

*z*from Eq. (7) resulted in a value of 0.03 m while linear interpolation resulted in 0.82 m. From Eq. (10) and neglecting stability effects, the normalized wind speed above the zero-plane displacement can be expressed as

_{me}To minimize the influence of stability, wind data during the winter period at the 19.25-m height were screened for all data above 5 m s^{−1} coming from an appropriate fetch direction; stability effects on were ignored. The normalized wind speeds for the 15- and 9-m heights, taken as the left-hand side of Eq. (23), were computed, averaged, and plotted against the normalized height, taken as the right-hand side of Eq. (23). The resulting plot presented in Fig. 1 shows that linear interpolation of *z _{me}* yields a much closer match to the 1:1 line compared to logarithmic interpolation.

The interpolation approach for *z _{me}* had relatively little influence on simulated wintertime energy fluxes for the aspen. RMSD for the aspen understory sensible heat flux was within 0.1 W m

^{−2}whether logarithmic or linear interpolation was used; the interpolation approach changed RMSD for sensible heat flux above the aspen by less than 0.5 W m

^{−2}. However, it did impact simulated temperature profiles within the canopy. The simulated canopy air temperature profile for both K theory and L theory resulted in much poorer simulations using logarithmic interpolation. RMSD for simulated temperature for the 3-m height was 2.8°C for L theory with logarithmic interpolation compared to 0.9°C using linear interpolation. K theory simulated canopy air temperature using logarithmic interpolation resulted in RMSD of 1.2°C compared to 0.7°C using linear interpolation.

The interpolation approach had only a small influence on simulated fluxes for the sagebrush site, but logarithmic interpolation actually tended to be slightly better. RMSD for K theory using logarithmic interpolation versus linear interpolation was 22.9 versus 23.6 W m^{−2} during the winter and 47.3 versus 47.4 W m^{−2} during the growing season. RMSD for wintertime sensible heat fluxes using L theory were 11.6 W m^{−2} using logarithmic interpolation versus 11.5 W m^{−2} using linear; growing season RMSD was 26.2 W m^{−2} for logarithmic and 26.4 W m^{−2} for linear. Leaf area index of the aspen during the growing season was sufficiently high that there was no difference between linear and logarithmic interpolation of *z _{me}*. Based on the combined results for the aspen and sagebrush sites, a linear interpolation for the surface roughness parameters was adopted for the current study, but further study is warranted.

### c. Energy fluxes for the aspen canopy

Above-canopy energy fluxes during the snow-covered period are plotted in Fig. 2. Interestingly, even though positive air temperatures exist above a 0° snowpack during active snowmelt (after day 70), measured midday sensible heat fluxes were actually upward, indicating unstable conditions at the top of the canopy. This is because solar heating of the aspen trees created upward sensible heat flux in the upper layers of the canopy. Obviously somewhere within the canopy, heat flux must be downward toward the snow surface during snowmelt. However, even at the 4.5-m height within the understory, there was sufficient solar warming of the aspen trees that midday measured sensible heat flux was positive, or upward, as plotted in Fig. 3. This suggests the existence of an extremely stable layer immediately above the snow surface. Measuring flux across this thin layer is somewhat challenging as measurements would need to made much closer to the snow surface.

During the winter season, there was very little difference in simulated energy fluxes between the K-theory and L-theory approaches either above the canopy or within the understory. RMSD and MBE were slightly better for above-canopy fluxes from K theory (Table 1), but results are mixed for the understory. Both approaches were able to mimic the midday unstable air conditions at the 4.5-m height during active snowmelt caused by the solar warming of the aspen limbs, yet captured the downward sensible heat flux at the snow surface, as shown in Fig. 3. Simulations with trees having an albedo of 1.0 yielded downward sensible heat fluxes throughout the canopy, confirming the hypothesis that solar warming was the cause for the upward heat flux (data not shown).

Canopy air temperature profiles were represented well by both approaches. Statistics comparing simulated and measured air temperatures are summarized in Table 2. MBE shows slightly better results for K theory, but RMSD values show mixed results for the two approaches.

Snow depth measurements indicate that the snowpack melted out on day 88. Both L-theory and K-theory approaches melted the snowpack a few days prematurely. K theory simulated complete melt on day 83 and L theory on day 84. However, both approaches are a considerable improvement over the original SHAW 2.5 version of the model that used above-canopy stability corrections within the canopy; it simulated complete melt on day 75 with RMSD and MBE for below-canopy sensible heat flux of 66.5 and −50.8 W m^{−2}. MBE for 3-m canopy temperature of −2.6°C was overly influenced by the biased sensible heat flux. Clearly melt dynamics at the site are radiation dominated (Table 1), as observed average turbulent fluxes are small and tend to cancel each other. Simulated sensible heat in the original SHAW 2.5 played far too much of a role in the melt, with downward simulated sensible heat flux throughout the canopy. Replacement of Eq. (5) in the original SHAW 2.5 was the main source of improvement; results using Eqs. (1)–(4) for above-canopy fluxes yielded nearly identical results as using Eqs. (10)–(14) for the extended K-theory approach.

There was very little difference between the K-theory and L-theory approaches for the simulated energy balance above the aspen canopy during the growing season, as shown in Fig. 4 and Table 3. Either approach did well, except RMSD for above-canopy latent heat flux tended to be high for both—around 60 W m^{−2}—and the MBE of around +25 W m^{−2} indicated a tendency to overpredict evapotranspiration. Most of this bias can be attributed to aspen transpiration being overpredicted in the afternoon hours, perhaps because of stress-related stomatal closure not simulated by the model. Stomatal closure within the model is controlled solely by leaf water potential and has no provision for stress due to elevated temperature or vapor deficit. These observations suggest that these provisions should be incorporated into the model, but are beyond the focus of this paper.

Overprediction of afternoon evapotranspiration is apparent in the understory as well, but not as significant as for the aspen (not shown). MBE for the understory was +11.5 W m^{−2} for K-theory approach and +12.1 W m^{−2} for the L-theory approach. As the soil profile dried, turbulent heat transfer shifted from latent heat to sensible heat flux, particularly for the understory, which has shallower roots than the aspen trees.

Unlike the winter season, the canopy air temperature profile was consistently represented better using L-theory turbulent transfer compared to K theory. RMSD of the 3-m simulated canopy air temperature was 2.5°C using K theory and 1.4°C using L theory (Table 4). K theory, in particular, overpredicted midday temperatures when unstable conditions tended to prevail within the canopy (Fig. 5); the reason for this is not clear.

### d. Energy fluxes for the sagebrush site

Simulating net solar radiation above the sagebrush site was more challenging during the winter season compared to the growing season (Table 5) because of brush progressively protruding above the snow as melt progressed, drastically changing the albedo. RMSD for net radiation was 18.3 W m^{−2} during the melt period compared to 8.7 W m^{−2} during the growing season. Unfortunately, measured snow depth was not available for the sagebrush, but based on the six soil heat flux plates, the ground became progressively exposed between late afternoon (hour 17) on day 77 and day 84. Complete melt was simulated on day 83 using the K theory and day 90 using L theory. The original SHAW 2.5 version of the model that used above-canopy stability corrections within the canopy performed much better at the sagebrush than at the aspen site; complete melt at the sagebrush site was simulated on day 84. The reason for this is likely because within-canopy stability effects had less influence at this windy exposed site with sparse cover, particularly before a significant amount of the sagebrush canopy was exposed.

During most of the melt period, sensible heat flux was downward (negative) above the sagebrush; however, late in the melt season, as the sagebrush became more exposed and absorbed more solar radiation, simulated sensible heat flux above the sagebrush was upward (positive), even as the snow melted. RMSD of sensible heat flux was lower using L theory (11.5 W m^{−2} compared to 23.6 W m^{−2}; Table 5), but MBE was worse (+3.7 W m^{−2} compared to +0.6 W m^{−2}). Figure 6 presents L-theory sensible heat flux during the period that the sagebrush became exposed as the snow melted away. The top of the sagebrush became exposed in the simulations on day 75; measured sensible heat flux became positive on day 77 and simulated flux approached zero. High wind speeds from the afternoon of day 78 through the morning of day 79 resulted in strong negative sensible heat fluxes, overcoming the warming influence of the absorbed solar radiation by the exposed sagebrush. Simulated flux at the surface (Fig. 6) illustrates that the heat fluxes measured above the sagebrush cannot be used to infer energy balance and melt processes at the snow surface, even for a short-statured sparse canopy such as the sagebrush.

Growing season fluxes presented in Fig. 7 show good agreement between measured and simulated fluxes. As with the aspen simulation, afternoon evapotranspiration was overpredicted. MBE was around +16 W m^{−2} for both turbulent approaches and RMSD was around 34 W m^{−2}. Even so, reduction in evapotranspiration as the soil profile dried (soil water data not shown) was captured well by the model. While latent heat flux was simulated almost identically for the two approaches, sensible heat flux was simulated much better using L theory; RMSD was 26.4 W m^{−2} for L-theory sensible heat flux and 47.4 W m^{−2} for K theory. Turbulent flux at the ground surface was lower using K theory, resulting in a positive bias in sensible heat flux above the canopy (+17.7 W m^{−2}); this was offset by a negative bias in net longwave radiation (−28.1 W m^{−2}) as the simulated plant litter on the soil surface was too warm.

## 5. Summary and conclusions

Intermingling of plant communities is a common feature in rangeland and forest landscapes. The resulting mosaic of vegetation communities creates a challenge for measuring and modeling the energy and mass fluxes in these complex environments, as adequate fetch requirements and complex topography make interpretation of EC results difficult.

The purpose of this study was to define areas for improvement in the SHAW multilayer plant canopy model and to evaluate observed energy fluxes above and within two diverse rangeland sites. Model modifications included incorporating equations for computing wind profile roughness parameters for sparse canopies, extending the K-theory approach to include influence of the roughness sublayer and stability functions within the canopy, and introducing Lagrangian far-field turbulent transfer equations (L theory) in lieu of the K-theory approach used previously. Observed wind profiles and simulation results suggest that the surface roughness parameter computed using the logarithmic interpolation presented in the literature was too small and a linear interpolation was adopted instead. Even so, further study to investigate the appropriate approach for estimating the effective surface roughness parameters for sparse canopies is warranted. The most pronounced improvement from the modifications came from changes in the within-canopy turbulent transfer and stability corrections. There was little difference in simulated energy fluxes for the aspen canopy using L-theory versus K-theory turbulent transfer equations, but L theory tracked canopy air temperature profiles better during the growing season. Sensible heat flux above the sagebrush was simulated markedly better using L theory compared to K theory.

Upward sensible heat flux was observed during active snowmelt both above the aspen trees and within the understory—even within 4 m of the snow surface. EC applications to understory vegetation are subject to much caution; however, model simulations confirmed that observed upward sensible heat flux during snowmelt was due to solar heating of the aspen limbs. Thus, unstable atmospheric conditions existed throughout most of the aspen canopy during melt; although not directly measured, very stable air conditions must have been present near the snow surface. Measuring fluxes within this small stable air layer would be challenging as measurements would need to be made much closer to the snow surface.

Upward sensible heat flux was observed during snowmelt even for the sparse low-statured sagebrush during periods of low wind speed. The influence of the sagebrush in producing these upward fluxes was captured by the model and downward sensible heat flux was simulated at the snow surface. Thus, the three EC systems did not quantify fluxes at the snow surface when vegetation was present, illustrating the difficulty in using EC systems to quantify energetics of snowmelt processes.

Growing season fluxes were simulated reasonably well using either K theory or L theory, but both tended to overestimate afternoon evapotranspiration, particularly above the aspen canopy. Provisions for stress-related stomatal closure beyond the current approach using leaf water potential are recommended for future model development.

Combining flux measurements with model simulations can help define knowledge gaps in current understanding of flux processes and lend credence to measured and/or simulated fluxes. Good agreement between the measured and modeled energy fluxes suggest that they can be measured and simulated reliably in these complex environments, but care must be used in the interpretation of the results.

## Acknowledgments

USDA is an equal opportunity provider and employer.

## REFERENCES

_{2}, H

_{2}O and sensible heat within and above a mountain meadow canopy: A comparison of three Lagrangian models and three parameterisation options for the Lagrangian time scale