Estimating Understory Temperatures Using MODIS LST in Mixed Cordilleran Forests

Satellite remote sensing provides a rapid and broad-scale means for monitoring vegetation phenology and its relationship with fluctuations in air temperature. Investigating the response of plant communities to climate change is needed to gain insight into the potentially detrimental effects on ecosystem processes. While many studies have used satellite-derived land surface temperature (LST) as a proxy for air temperature, few studies have attempted to create and validate models of forest understory temperature (Tust), as it is obscured from these space-borne observations. This study worked to predict instantaneous values of Tust using daily Moderate Resolution Imaging Spectroradiometer (MODIS) LST data over a 99,000 km2 study area located in the Rocky Mountains of western Alberta, Canada. Specifically, we aimed to identify the forest characteristics that improve estimates of Tust over using LST alone. Our top model predicted Tust to within a mean absolute error (MAE) of 1.4 ◦C with an overall model fit of R2 = 0.89 over two growing seasons. Canopy closure and the LiDAR-derived standard deviation of canopy height metric were found to significantly improve estimations of Tust over MODIS LST alone. These findings demonstrate that canopy structure and forest stand-type function to differentiate understory air temperatures from ambient canopy temperature as seen by the sensor overhead.


Introduction
There is growing interest in using remote sensing for observing Earth's microclimates in order to acquire a better understanding of ecosystem functions such as vegetation phenology, carbon flux, and the timing of species trophic interactions (e.g., [1]).Surface air temperature (T air ) is a critical driver of these environmental processes, especially at high latitudes where seasonal transitions are more highly pronounced [2].However, obtaining reliable measurements of T air across the landscape is often hindered by missing values due to gaps in the data [3].These voids are typically a result of the limited density of meteorological stations distributed across the landscape [4].Spatial interpolation of these data is prone to error because of the distance and fine-scale temperature aberrations between stations, especially in mountainous areas where their distribution tends to be even sparser [5][6][7].Remote sensing of land surface temperature (LST) helps to resolve this issue by providing an inherently spatialized gridded surface that covers even the most remote regions.The Moderate Resolution Imaging Spectroradiometer (MODIS) has proven invaluable for monitoring year-round global temperature in high temporal detail [8].MODIS LST data is particularly popular because of its high level of post-processing which produces accuracies of ≤1 • C [9].MODIS LST has successfully Remote Sens. 2016, 8, 658 2 of 21 been used to predict T air (e.g., [10,11]), however T air recordings are commonly made in open areas or high above the forest canopy from eddy covariance flux towers (e.g., [12,13]).With observations at heights in excess of 30 m, the relationship between MODIS LST and air temperature in the understory (T ust ) is relatively unknown.
Radiometers do not measure air temperature directly, rather the LST estimations are obtained through a split-window algorithm which translates thermal infrared (TIR) observations into skin temperature of the observable land surfaces (i.e., bare land, urban areas, forest canopy) [14].LST is quite different from T air and they do not generally correlate well [15].This relationship is improved over vegetated areas, since emitted radiance is affected by surface type [16].There is a simple and physical relationship between fractional vegetation cover and thermal properties detected by the sensor [17].The disparity between LST and T air is greater over areas of exposed soil and other surfaces with low emittance [18].As a result, the direct relationship between LST and T air needs to be determined using a supplementary estimation technique.
Variable forest cover types introduce additional complexities by creating irregularities in the latent heat flux and other radiative exchanges [19,20].The prevalent methodologies used to parameterize T air from LST include (i) physical models of the surface energy balance; and (ii) the temperature-vegetation index (TVX) method [21].Energy-balance models are based on physical processes, and require large amounts of radiation and latent heat flux data typically not provided by remote-sensing techniques [15,22].The spectral information obtained by the sensor is not able to perceive many of these complex physical variables, such as wind speed and soil moisture [23].
The challenge in measuring the processes that govern T air above vegetation leave most models under-parameterized, and this is the main reason for the development of the temperature-vegetation index (TVX) method [17].The premise of the TVX method is an established linear and negative correlation between a remotely-sensed vegetation index, often the normalized difference vegetation index (NDVI), and satellite-derived LST [24,25].The TVX assumes that the radiometric temperature of a fully vegetated canopy is in equilibrium with T air [26].The premise is that dense vegetation radiates heat so efficiently that it remains close in temperature to the surrounding ambient air [18].In effect, the assumption is that the top-of-canopy temperature is the same as the canopy itself [11].How this assumption translates to estimating understory temperatures, and the relationship with MODIS LST across different vegetative strata, requires further investigation.
The difference between air temperatures above and below the forest canopy is generally a result of near-surface environmental lapse rates and distance from the ground [1].The perspective of forests from a space-borne sensor is of the top of the canopy, and the thermal information beneath is hidden from view due to the relative opacity of the forest itself.The boundary layer separating the emissive properties of the understory and the supracanopy air mass is occupied by different types and densities of canopy vegetation.Therefore estimations of T ust , without assuming that it is in equilibrium with the canopy, should incorporate physical characteristics of the canopy in addition to MODIS LST.These characteristics can be acquired in a number of ways.Traditionally, forest attributes have been obtained through mensuration techniques developed by the timber industry [27].Metrics such as tree height, diameter, and canopy closure are generally used to estimate the volume of harvestable timber in a stand.These methods have since been adapted to characterize forest landscapes for ecosystem management and studies of species ecology (e.g., [28]).Recently, these conventional mensuration practices have been complemented by the development of airborne laser scanning (LiDAR) data.LiDAR systems procure detailed models of vegetation in three dimensions, with the fundamental benefit of penetrating the canopy to extract additional information from the understory [29,30].Most of the existing literature emphasizes characteristics of the overstory, with relatively few papers using LiDAR to examine the understory within an ecological context [31][32][33].Nijland et al. [34] compared LiDAR-derived canopy metrics with more conventional climate and land cover based variables to model the distribution of various understory species with mixed success.
The importance of estimating understory temperature is to increase our understanding of its cause-effect relationship with other natural processes.For example, the timing of phenological events is closely linked to changes in ambient temperature, where the life histories of many species are explicitly adapted to the seasonal environments in which they live [35].Monitoring microclimate dynamics provides insight into species sensitivities to weather anomalies and global trends in temperature from climate change [36,37].The MODIS LST imagery offers a mechanism to observe this significant driver of ecosystem dynamics over the broad scale, which in turn provides maps as a basis for species conservation, habitat management, or identifying additional environmental responses to climate change [38][39][40][41].The present effort to model T ust is motivated by our involvement in an ongoing grizzly bear (Ursus arctos) research program where we are working to monitor critical habitat for this threatened species in Alberta, Canada [42].These animals strategically exploit a wide variety of understory plants that provide variable nutrition at differential times of the growing season [43].The eventual goal is to model the bears' response to shifting understory phenology by mapping T ust accumulations, or growing degree days (GDD) (e.g., [44][45][46]).
The objective of this study is to predict instantaneous measurements of T ust using daily MODIS LST data and a combination of LiDAR-derived canopy metrics and conventional forest inventory variables.Although MODIS LST has been previously correlated with T air [47][48][49], estimating T ust has not been studied in detail.There is also an absence of in situ meteorological data measured at comparable resolutions to the remotely sensed imagery with which it is correlated [1].Using a multivariate statistical approach, we ranked the comparative importance of LiDAR and conventional forest inventory variables in estimating instantaneous T ust at the MODIS pixel-scale.In situ T ust measurements were collected over two seasons using a widespread network of meteorological stations distributed throughout the foothills and Rocky Mountains of western Alberta, Canada.Although it is normally assumed that the vegetated canopy is the same temperature as the surrounding ambient air [18,26], the understory will experience dissimilar temperatures as a result of the general moderating effect of the canopy [50].We hypothesized that based on the interactions between MODIS LST and forest canopies [16], that supplementary forest canopy and stand variables would significantly improve model predictions of T ust over those made by LST alone.

Study Area
The study area extends over 700 km along the eastern slopes of the Alberta Rocky Mountains (49 • N to 55 • N) in western Canada (Figure 1).T ust observation plots were placed at near-regular latitudinal intervals along elevational transects covering a gradient from 800 m to 1800 m.The forests of the region are typified by deciduous trembling aspen (Populus tremuloides) and balsam poplar (Populus balsamifera) in mixed stands with white spruce (Picea glauca) at lower elevations.The mountains are dominated by mature conifer forests of lodgepole pine (Pinus contorta) and Engelmann spruce (Picea engelmannii).Stand composition of the observation plots ranged from 40% being purely coniferous, 40% mixed, and 20% entirely deciduous.

In Situ Tust Observations
Within the study area, twenty-nine broadly dispersed Tust observation plots were selected to obtain a range of temperature regimes, forest composition types, and variability in canopy closure.Each plot was 250 × 250 m (62,500 m 2 ) in size and located within comparatively homogenous forest stands that extended a minimum of 250 m beyond the edges of the plot.This spatial buffer retained representative temperatures for that particular stand type, minimizing the influence from surrounding divergent cover types.This stand size also corresponds with the areal unit scale of a 1 km MODIS LST pixel (Figure 2).Within each plot, four meteorological sensors recorded 8-bit hourly temperature from early spring (pre green-up) until after autumn senescence for two growing seasons: 2011 and 2012.The sensors were enclosed in solar radiation shields 1 m above ground-level, approximating the typical height of shrubby understory vegetation.

In Situ T ust Observations
Within the study area, twenty-nine broadly dispersed T ust observation plots were selected to obtain a range of temperature regimes, forest composition types, and variability in canopy closure.Each plot was 250 × 250 m (62,500 m 2 ) in size and located within comparatively homogenous forest stands that extended a minimum of 250 m beyond the edges of the plot.This spatial buffer retained representative temperatures for that particular stand type, minimizing the influence from surrounding divergent cover types.This stand size also corresponds with the areal unit scale of a 1 km MODIS LST pixel (Figure 2).Within each plot, four meteorological sensors recorded 8-bit hourly temperature from early spring (pre green-up) until after autumn senescence for two growing seasons: 2011 and 2012.The sensors were enclosed in solar radiation shields 1 m above ground-level, approximating the typical height of shrubby understory vegetation.

In Situ Tust Observations
Within the study area, twenty-nine broadly dispersed Tust observation plots were selected to obtain a range of temperature regimes, forest composition types, and variability in canopy closure.Each plot was 250 × 250 m (62,500 m 2 ) in size and located within comparatively homogenous forest stands that extended a minimum of 250 m beyond the edges of the plot.This spatial buffer retained representative temperatures for that particular stand type, minimizing the influence from surrounding divergent cover types.This stand size also corresponds with the areal unit scale of a 1 km MODIS LST pixel (Figure 2).Within each plot, four meteorological sensors recorded 8-bit hourly temperature from early spring (pre green-up) until after autumn senescence for two growing seasons: 2011 and 2012.The sensors were enclosed in solar radiation shields 1 m above ground-level, approximating the typical height of shrubby understory vegetation.

MODIS Land Surface Temperature Observations
The daily LST product is collected by the sun-synchronous, near-polar orbiting Terra (MOD11A1) and Aqua (MYD11A1) satellites, with Aqua in an ascending orbit and Terra in a descending orbit, having two daily equatorial crossings each at 13:30, 01:30, and 10:30, 22:30, respectively (local solar time) [51].LST is derived from the MODIS thermal emissivity band channels 31 (10.78-11.28µm) and 32 (11.77-12.27 µm).Atmospheric effects are corrected using a split-window algorithm that compares the differential thermal radiation absorption between these two bands [52].Version v005 was the most up-to-date MODIS LST product available, with an accuracy of ≤1 • K, and provided significantly improved spatial coverage, stability, and accuracy compared to previous versions [9,47].Data download and processing was automated using R programming language [53].The tiled daily LST data were obtained from the NASA Land Processes Distributed Active Archive Center (LP DAAC) for all dates corresponding with the in situ T ust observations.The MODIS Reprojection Tool [54] was used to reproject the MOD11A1, MYD11A1, and coincident local solar view time products from Sinusoidal to 1 km gridded geotiff's in Universal Transverse Mercator projection (UTM Zone 11N, NAD83 datum) using nearest-neighbor resampling.Pixel values were then rescaled and converted to degrees Celsius.An abridged methodology is shown in Figure 3.

MODIS Land Surface Temperature Observations
The daily LST product is collected by the sun-synchronous, near-polar orbiting Terra (MOD11A1) and Aqua (MYD11A1) satellites, with Aqua in an ascending orbit and Terra in a descending orbit, having two daily equatorial crossings each at 13:30, 01:30, and 10:30, 22:30, respectively (local solar time) [51].LST is derived from the MODIS thermal emissivity band channels 31 (10.78-11.28μm) and 32 (11.77-12.27 μm).Atmospheric effects are corrected using a split-window algorithm that compares the differential thermal radiation absorption between these two bands [52].Version v005 was the most up-to-date MODIS LST product available, with an accuracy of ≤1°K, and provided significantly improved spatial coverage, stability, and accuracy compared to previous versions [9,47].Data download and processing was automated using R programming language [53].The tiled daily LST data were obtained from the NASA Land Processes Distributed Active Archive Center (LP DAAC) for all dates corresponding with the in situ Tust observations.The MODIS Reprojection Tool [54] was used to reproject the MOD11A1, MYD11A1, and coincident local solar view time products from Sinusoidal to 1 km gridded geotiff's in Universal Transverse Mercator projection (UTM Zone 11N, NAD83 datum) using nearest-neighbor resampling.Pixel values were then rescaled and converted to degrees Celsius.An abridged methodology is shown in Figure 3. Accurate LST measurements are highly dependent on the amount of cloud cover at the time of overpass, therefore a two-step filter process was adapted from Neteler [6] and Zorer et al. [46] to obtain only the highest quality pixels for each observation plot.The first step used the integrated image bitword quality assessment (QA) flags to mask out cloud or corrupted image pixels.Data were selected with "0" flags in each quality control field, e.g., 00 in "Mandatory QA flags", 00 in "Data quality flag", average emissivity error ≤0.01, average LST error ≤1°K, and view-angle less than 40 degrees from nadir [55].The second step involved a 1 km proximity search around each observation plot pixel, considering that cloud-contaminated pixels not detected by the QA flags most commonly occur along the edges of clouds [6,56].Any adjacent pixel with more than a 5°K difference from an observation plot pixel was considered corrupted (i.e., from vapor or cloud shadow) and not used in the model.Overall, 1920 MODIS LST images were processed.
To obtain instantaneous measurements of in situ Tust, temporally linking observations with coincident MODIS LST view times was critical.This procedure can be problematic, considering MOD11A1 and MYD11A1 tile data can contain a mosaic of non-sequential and temporally unrelated surface temperatures resulting from multiple daily imaging from adjacent swaths [48].Off-nadir viewing is performed to maximize the amount of cloud-free ground surface visible to the sensor during each orbit.To separate the image in to spatially corrected view times, a longitudinal correction Accurate LST measurements are highly dependent on the amount of cloud cover at the time of overpass, therefore a two-step filter process was adapted from Neteler [6] and Zorer et al. [46] to obtain only the highest quality pixels for each observation plot.The first step used the integrated image bitword quality assessment (QA) flags to mask out cloud or corrupted image pixels.Data were selected with "0" flags in each quality control field, e.g., 00 in "Mandatory QA flags", 00 in "Data quality flag", average emissivity error ≤0.01, average LST error ≤1 • K, and view-angle less than 40 degrees from nadir [55].The second step involved a 1 km proximity search around each observation plot pixel, considering that cloud-contaminated pixels not detected by the QA flags most commonly occur along the edges of clouds [6,56].Any adjacent pixel with more than a 5 • K difference from an observation plot pixel was considered corrupted (i.e., from vapor or cloud shadow) and not used in the model.Overall, 1920 MODIS LST images were processed.
To obtain instantaneous measurements of in situ T ust , temporally linking observations with coincident MODIS LST view times was critical.This procedure can be problematic, considering MOD11A1 and MYD11A1 tile data can contain a mosaic of non-sequential and temporally unrelated surface temperatures resulting from multiple daily imaging from adjacent swaths [48].Off-nadir viewing is performed to maximize the amount of cloud-free ground surface visible to the sensor during each orbit.To separate the image in to spatially corrected view times, a longitudinal correction was made for each plot and the equation of time (1) was applied to the image solar view times to yield local view times (±30 s) with corrections for daylight savings (2) [57].
where the angle x is defined as a function of the day of year and produces the difference between mean solar time and true solar time on a given date.
where t s is solar time, EOT is the equation of time, and LC is the site longitudinal correction from the nearest standard meridian.Once calculated, local LST view times were then rounded to the nearest hour and coupled with the corresponding hourly T ust plot measurement for every clear sky observation during the study period.

LiDAR-Derived Canopy Metrics
Airborne laser scanning data for the study area were provided by the Government of Alberta's Forest Management Branch.The data were collected between 2003 and 2010 (although the majority is circa 2007 during predominantly leaf-on conditions) using sensors capable of detecting four returns per pulse.Data densities were typically 1-2 returns/m 2 .The LiDAR data were normalized to height above ground-level and processed into a suite of height and structure metrics using standard processing routines available in FUSION [58].The majority of metrics used first returns in accordance with recommendations in Bater et al. [59], who found that first-return vegetation-height metrics were more stable than those calculated from all-pulse returns.The canopy metrics were averaged over the entire 250 m study plot to approach the areal-unit scale of the MODIS imagery.
Individual LiDAR-derived metrics were selected for use in the models based on their most plausible and intrinsic role in explaining the T ust /LST interface, and also on their predictive success in previous structural understory LiDAR models (e.g., [34]) (Table 1).For example, the ratio pa14 is considered a proxy for canopy closure (1.4 m being the understory height threshold [60]), while standard deviation of canopy height (stdd) can be thought of as absolute vegetation roughness [61].A number of indices were included because they incorporate all pulse returns.Although being less stable, they potentially provide more information on subcanopies and understory vegetation than metrics derived from first returns alone [62].This was an effort to maximize the utility of the LiDAR data, as it only provides information on vegetation structure, with no additional nuanced ecological information on species assemblies and abundance.

Conventional Forest and Topographic Variables
Before the widespread availability of LiDAR and optical remote sensing, forest structure and compositional characteristics were (and still are) traditionally collected by personnel on the ground [27].
These ground data, combined with LST and derived topographic variables, have provided suitable model results for estimating T air .For example, Parmentier et al. [20] used a number of covariates such as forest land cover type, elevation, and aspect to predict monthly average T air to approximately 2.5 • C. Variables representing land cover and vegetation type are particularly important since they affect surface emissivity [14].Differences between AM and PM LST values have been found to vary between 0.3 • C to 3.2 • C depending on cover type [63].Taking into consideration this sensitivity of LST to land cover and canopy density, forest tree species composition (percent conifer) in each of the 29 observation plots was measured every 50 m along two 250 m stratified-random linear transects (14 subplots per plot).A representative sample of tree species was obtained using a 360 • sweep of each subplot using a wedge prism relascope (Figure 4).Canopy closure was measured during leaf-on conditions using five hemispherical photographs at each subplot and later processed with WinSCANOPY software [64] to calculate percent closure.Increased closure was assumed to have insulative properties that provide thermal cover for the ambient understory air mass [16,65].
Remote Sens. 2016, 8, 658 7 of 21 approximately 2.5 °C.Variables representing land cover and vegetation type are particularly important since they affect surface emissivity [14].Differences between AM and PM LST values have been found to vary between 0.3 °C to 3.2 °C depending on cover type [63].Taking into consideration this sensitivity of LST to land cover and canopy density, forest tree species composition (percent conifer) in each of the 29 observation plots was measured every 50 m along two 250 m stratified-random linear transects (14 subplots per plot).A representative sample of tree species was obtained using a 360° sweep of each subplot using a wedge prism relascope (Figure 4).Canopy closure was measured during leaf-on conditions using five hemispherical photographs at each subplot and later processed with WinSCANOPY software [64] to calculate percent closure.Increased closure was assumed to have insulative properties that provide thermal cover for the ambient understory air mass [16,65].A 10 m photogrammetrically compiled digital elevation model was used to extract topographic variables including elevation, slope, aspect, and annual terrain solar radiation in watt-hours per square meter (WH/m 2 ) using ArcGIS 10.3 software [66].Throughout the year the difference between same-day AM and PM LST observations can vary up to 5 °C [63].Because of this, Julian day was included in the models to account for this variability in the surface energy balance at different times of the growing season.Benali et al. [47] found Julian day to have little predictive effect in their Tair models, so the quadratic polynomial of Julian day was used in this study to better represent the seasonal temperature curve (Table 2).A 10 m photogrammetrically compiled digital elevation model was used to extract topographic variables including elevation, slope, aspect, and annual terrain solar radiation in watt-hours per square meter (WH/m 2 ) using ArcGIS 10.3 software [66].Throughout the year the difference between same-day AM and PM LST observations can vary up to 5 • C [63].Because of this, Julian day was included in the models to account for this variability in the surface energy balance at different times of the growing season.Benali et al. [47] found Julian day to have little predictive effect in their T air models, so the quadratic polynomial of Julian day was used in this study to better represent the seasonal temperature curve (Table 2).

Statistical Methods and Model Ranking
Regression approaches have previously demonstrated accurate predictions of T air using solely LST or LST associated with ancillary variables [4,13,46,49,67].In this application, models were divided into three groupings similar to the organization used by Benali et al. [47], where they modeled instantaneous MODIS LST observations during the (i) daytime; (ii) night; and (iii) using both datasets combined (hereafter referred to as the daytime, night, and day-night models).Separating the dataset in this manner was intended to isolate any effect from differing diurnal near-surface temperature gradients, insolation variability during the day, and outgoing radiation at night.Within the view-time groupings, candidate models were organized into subgroupings established on explanatory variable type (conventional variables, LiDAR-derived metrics, and both types combined) to assess any marked improvement in model performance between them.Generalized linear models were derived in STATA 13 statistical software [68] to estimate instantaneous values of T ust .Random effects (GLS estimator) were used to offset bias from using multiple temperature observations within single sites.A cross-correlation between explanatory variables tested for both linear and monotonic relationships to ensure no collinearity within models.Akaike information criterion (AIC) was used to perform rankings of the candidate models.This procedure seeks to find a balance between model fit (log likelihood) and the number of variables in the model, penalizing those with too much complexity, ranking the best-fit and most parsimonious model as the highest [69].A total of 4153 day and night clear-sky observations were made during the study period (Figure 5); 20% of these data were randomly withheld for validation [70].

Statistical Methods and Model Ranking
Regression approaches have previously demonstrated accurate predictions of Tair using solely LST or LST associated with ancillary variables [4,13,46,49,67].In this application, models were divided into three groupings similar to the organization used by Benali et al. [47], where they modeled instantaneous MODIS LST observations during the (i) daytime; (ii) night; and (iii) using both datasets combined (hereafter referred to as the daytime, night, and day-night models).Separating the dataset in this manner was intended to isolate any effect from differing diurnal nearsurface temperature gradients, insolation variability during the day, and outgoing radiation at night.Within the view-time groupings, candidate models were organized into subgroupings established on explanatory variable type (conventional variables, LiDAR-derived metrics, and both types combined) to assess any marked improvement in model performance between them.Generalized linear models were derived in STATA 13 statistical software [68] to estimate instantaneous values of Tust.Random effects (GLS estimator) were used to offset bias from using multiple temperature observations within single sites.A cross-correlation between explanatory variables tested for both linear and monotonic relationships to ensure no collinearity within models.Akaike information criterion (AIC) was used to perform rankings of the candidate models.This procedure seeks to find a balance between model fit (log likelihood) and the number of variables in the model, penalizing those with too much complexity, ranking the best-fit and most parsimonious model as the highest [69].A total of 4153 day and night clear-sky observations were made during the study period (Figure 5); 20% of these data were randomly withheld for validation [70].

Results
The quadratic function of Julian day and LST was set as the standard null model, as these two variables were found to explain more than 80% of the variability in Tust.Julian day explained 20%

Results
The quadratic function of Julian day and LST was set as the standard null model, as these two variables were found to explain more than 80% of the variability in T ust .Julian day explained 20% more variability than LST alone, which is considered good for ecological phenomena [71].Using all clear-sky observations, separate baseline model rankings of the two variable types (conventional and LiDAR-derived) was performed to get an overview of their relative importance in predicting T ust (Table 3).The delta AIC (∆ i ) is a measure of each model relative to the top-ranked model, where the criterion suggests substantial evidence for models within ∆ i < 2. Values between 3 and 7 indicate that the model has considerably less support, and ∆ i > 10 indicates that the model is very unlikely [72].Akaike weights (w i ) can be interpreted as the probability that the top-ranked model is the best in the set of candidate models [69,73].Of the conventional metrics, the best model was explained (in addition to the null) by the interaction between LST and percent conifer, which describes the proportion of coniferous trees within a sample unit.The effect of the interaction term in the model can be interpreted as a linkage between these two variables where the influence of forest composition on T ust depends on the coincident LST value at the time of observation.Without the interaction, we are simply trying to predict T ust as solely the unique effect of percent conifer.The top two conventional metric models lead the rankings, each having Akaike weights of 0.50 (percent conifer and the interaction of LST with percent conifer).There was weak support for the addition of elevation, since it is equally ranked with the percent conifer interaction model despite increasing model complexity.By and large, employing conventional variables imparts a major improvement in estimating T ust over using LST alone, which was ranked very poorly.
There was no observed improvement in the prediction of T ust using exclusively LiDAR-derived metrics over the null model (Table 3).The best-performing was skewness, although there is some optimism in the strength of the other LiDAR models since they are all within 2 AIC values of the top ranked model (the null).Any LiDAR-derived metric model containing an interaction term had significantly reduced performance which was reflected by a large increase in AIC value.The overall top-ranked models from the final groupings (conventional, LiDAR-derived, and both combined) were directly compared to gauge which supplementary variable type best predicts T ust .This comparison was performed for each overpass dataset (daytime, night, and day-night).The combination of conventional variables and LiDAR-derived metrics clearly produced the best models, regardless of view time, as reflected by the large Akaike weights (Table 4).Percent canopy closure (cc) and the LiDAR-derived metric standard deviation of canopy height (stdd) both appear in all of the top-ranked models.They are often in the form of an interaction term with LST, which is a strong indication of an interrelationship with the forest canopy.At night, the standardized coefficients show that for every percent increase in canopy closure, there will be a 0.06 • C increase in T ust (Table 5).During the day, for every unit increase in closure, T ust decreases 0.12 • C. Also in the top daytime models was the percent returns above 1.4 m (pa14) metric, which is the LiDAR-derived equivalent to canopy closure.This metric had a similar negative daytime influence as closure, where an increase in this variable produces a slight drop in T ust .For both the daytime and day-night models, the interaction term of the LiDAR-derived standard deviation of canopy height was positive, and implies that when LST interacts with increasingly rough canopies, T ust increases, but very little.At night this relationship is inverted, and a higher standard deviation in height results in a cooling of T ust .The percent conifer coefficient was negative, suggesting that as the percentage of conifers increases, T ust will decrease.Overall, canopy closure and the standard deviation of canopy height were the most influential ancillary variables for estimating T ust .The daytime model was the best performing overall, with a model fit of R 2 = 0.89, and the lowest validation error of mean absolute error (MAE) = 1.4 • C.This prediction accuracy parallels estimations of T air [4,13,49], which is noteworthy considering the predictions are of temperatures beneath the forest canopy (Table 6; Figure 6).
The daytime model was the best performing overall, with a model fit of R 2 = 0.89, and the lowest validation error of mean absolute error (MAE) = 1.4 °C.This prediction accuracy parallels estimations of Tair [4,13,49], which is noteworthy considering the predictions are of temperatures beneath the forest canopy (Table 6; Figure 6).

Top Model Model Fit MAE RMSE BIAS
The best day-night model had the same model fit as the daytime model (R 2 = 0.89), but had a 0.5 • C increase in error (MAE of 1.9 • C).Although this model had no overall bias, the validation data show a slight T ust underprediction in cooler temperatures that occur during the very start and end of the growing season (Figure 7).The top model derived from night observations had the lowest model fit of the three MODIS view times (R 2 = 0.77) and a small error increase over the daytime model in estimation error (MAE = 2.1 • C).The model appears to underestimate higher temperatures and overestimate colder temperatures which seem to balance out, producing no predictive bias on average.The best day-night model had the same model fit as the daytime model (R 2 = 0.89), but had a 0.5 °C increase in error (MAE of 1.9 °C).Although this model had no overall bias, the validation data show a slight Tust underprediction in cooler temperatures that occur during the very start and end of the growing season (Figure 7).The top model derived from night observations had the lowest model fit of the three MODIS view times (R 2 = 0.77) and a small error increase over the daytime model in estimation error (MAE = 2.1 °C).The model appears to underestimate higher temperatures and overestimate colder temperatures which seem to balance out, producing no predictive bias on average.

Discussion
Forest type, density, and canopy structure all have an impact on the relationship between MODIS LST and understory air temperature.In this paper, we examined which metrics coupled with MODIS LST improve estimates of Tust.This is, to our knowledge, the first study of its kind to specifically model air temperatures in the understory using MODIS LST.We statistically predicted instantaneous values of Tust at a 250 m plot-scale to within MAE 1.4 °C-2.1 °C (depending on the time of day).This is comparable to previous canopy-top Tair models that had error ranging from 1.3 °C to 3.5 °C and similar model fit [4,10,13,47].The fit slightly exceeded that obtained by Hanes and Schwartz [1], who modeled maximum daily Tair at the base of a flux tower (1.5 m) (R 2 = 0.86).
The Tust estimates were consistently accurate over a range of forest types and densities.This may be a result of the overall effectiveness of the model variables at representing the forest canopy characteristics.This variability could not be explained by LST alone (null model) as demonstrated in the AIC rankings.Despite the cost of increased model complexity, those which contained canopy metrics consistently ranked higher than the LST null model.A likelihood ratio test was performed to determine if these top-ranked models estimated Tust significantly better than the null.The results found that the addition of canopy metrics produce a significant improvement over the LST null in the day-night combined model (p < 0.001), and marginally significant in the daytime and night models (p = 0.08; p = 0.07).
The improvement of the models through the addition of the canopy metrics implies that there is a distinction between canopy-top (Tair) and understory temperatures.The difference is likely a result of moderation from the canopy [50,74].Hanes and Schwartz found that as canopy vegetation becomes increasingly dense during spring leaf out, MODIS LST values begin to approach air temperatures at 1.5 m above ground.In areas where canopy density is reduced, the influence of the overlying forest is diminished and therefore LST should deviate from Tair.However, the Tust estimates across the study area are consistently accurate despite the broad range of forest densities in the plots (10%-70% closure).This finding implies that, unlike Tair, the relationship between LST and Tust is likely consistent regardless of canopy closure.
Modeling day and night LST as separate datasets highlighted daytime as the optimal diurnal period to observe Tust.The daytime model was best overall, with a 10% increase in model fit and 0.7 °C lower error than at night.The majority of previous Tair studies (e.g., [12]) have found increased accuracy during the night, which is attributed to less interference from reflected daytime radiation [75].The improved daytime accuracy in the Tust models may be owed to modeling entirely within

Discussion
Forest type, density, and canopy structure all have an impact on the relationship between MODIS LST and understory air temperature.In this paper, we examined which metrics coupled with MODIS LST improve estimates of T ust .This is, to our knowledge, the first study of its kind to specifically model air temperatures in the understory using MODIS LST.We statistically predicted instantaneous values of T ust at a 250 m plot-scale to within MAE 1.4 • C-2.1 • C (depending on the time of day).This is comparable to previous canopy-top T air models that had error ranging from 1.3 • C to 3.5 • C and similar model fit [4,10,13,47].The fit slightly exceeded that obtained by Hanes and Schwartz [1], who modeled maximum daily T air at the base of a flux tower (1.5 m) (R 2 = 0.86).
The T ust estimates were consistently accurate over a range of forest types and densities.This may be a result of the overall effectiveness of the model variables at representing the forest canopy characteristics.This variability could not be explained by LST alone (null model) as demonstrated in the AIC rankings.Despite the cost of increased model complexity, those which contained canopy metrics consistently ranked higher than the LST null model.A likelihood ratio test was performed to determine if these top-ranked models estimated T ust significantly better than the null.The results found that the addition of canopy metrics produce a significant improvement over the LST null in the day-night combined model (p < 0.001), and marginally significant in the daytime and night models (p = 0.08; p = 0.07).
The improvement of the models through the addition of the canopy metrics implies that there is a distinction between canopy-top (T air ) and understory temperatures.The difference is likely a result of moderation from the canopy [50,74].Hanes and Schwartz found that as canopy vegetation becomes increasingly dense during spring leaf out, MODIS LST values begin to approach air temperatures at 1.5 m above ground.In areas where canopy density is reduced, the influence of the overlying forest is diminished and therefore LST should deviate from T air .However, the T ust estimates across the study area are consistently accurate despite the broad range of forest densities in the plots (10%-70% closure).This finding implies that, unlike T air , the relationship between LST and T ust is likely consistent regardless of canopy closure.
Modeling day and night LST as separate datasets highlighted daytime as the optimal diurnal period to observe T ust .The daytime model was best overall, with a 10% increase in model fit and 0.7 • C lower error than at night.The majority of previous T air studies (e.g., [12]) have found increased accuracy during the night, which is attributed to less interference from reflected daytime radiation [75].The improved daytime accuracy in the T ust models may be owed to modeling entirely within forested sites and the inherent difference of estimating T ust as opposed to T air .Essentially, the radiative exchanges in the understory are sufficiently unique to those above the canopy resulting in improved daytime T ust observations.In addition, the predominance of conifer trees in the study area may have augmented the performance of the daytime model, since low albedo coniferous canopy absorbs and emits energy more effectively than deciduous [76]; green-needle forest is considered a near-perfect emitter of thermal radiation [77].This is supported in a daytime model by the positive interaction between LST and percent conifer.This same interaction term is not present in the night model, perhaps because there is no potential for interaction with sunlight.
Conifer forests are also drier than deciduous, ultimately affecting the latent heat flux between them throughout the diurnal cycle.Williamson et al. [48] and Niclòs et al. [4] found that soil moisture was a significant predictor of T air as it affects surface emissivity.Considering at-site soil moisture was not measured in this study, GIS-derived variables could have been used as proxies for moisture such as a wet area map [34].Cresswell et al. [78] explain that since the Earth's skin surface temperature is generally higher than T air during the day and cooler at night, this will lead to overestimation and underestimation of air temperature in diurnal T air models, respectively.However, this estimation bias was not observed between the T ust models, perhaps from the moderating effect of the canopy.At night, there was a positive LST interaction with canopy closure which in turn increased T ust .During the day, this relationship was negative and for every percent increase in closure, T ust decreased 0.12 • C.These results corroborate that the canopy shades the understory during the day and insulates it at night [33,79].This also implies that the understory microclimate is unique to the overlying canopy and therefore emphasizes the distinction between T air and T ust .
Regarding the performance of the LiDAR-derived metrics, they contributed to all of the top-ranked models, indicating some overall utility.The standard deviation of canopy height was the most predictive either as a solitary variable (at night) or interacting with LST (daytime).It would be more intuitive that average height or canopy density (percent returns above 1.4 m) should explain more variability in T ust ; rather it was the roughness of the forest as expressed by the standard deviation of height [33].During the daytime T ust rises when forest structure becomes increasingly rough, while at night the opposite occurs, where a higher standard deviation results in lower T ust .A rough surface would be characteristic of forest with a strong variation in age structure or disturbance patterns such as trails or cut lines.Although not in the top-ranked models, LiDAR-derived skewness in the distribution of tree heights appeared strongly in many of the models, having a significant yet subtle influence on T ust .Skewness is a similar metric to standard deviation and signifies canopy complexities resulting from variable age structure or directional trends in height across the plot due to soil conditions or slope [80].The correlation of this metric with T ust was at all times negative, suggesting that as skewness in height increases, T ust drops slightly.The relative contribution to model improvement using the LiDAR metrics raises the question of the advantage of these data over the conventional and GIS-derived variables.Ultimately, it depends on the scale and objectives of the research; in this T ust analysis, the benefits of the LiDAR metrics were not as substantial as the conventional variables.Conversely, the cost of traditional ground measurements are considerably more than LiDAR surveys [81], especially in contrast to the increasing availability of LiDAR technology and platforms such as unmanned aerial systems [82,83].
Downscaling the 2-pulse-per-meter LiDAR metrics to the 250 m plot-scale could have impacted their effectiveness as explanatory variables of T ust .This was done to match the areal unit scale of the observation plots, however it is unlikely that the effect was significant as averaging would only remove the extreme and outlying values.Tompalski et al. [84] compared stand-level LiDAR estimates of tree volume to individual volume measurements and found them to be very similar.The areal scale of 250 m was conceivably an improvement over point-source temperature measurements made by lone meteorological stations or flux towers.The relatively low root mean square error (RMSE) in the T ust models may be indicative of measuring in situ temperature at a broader 250 m scale instead of discrete points on the landscape.The observation plots were located in the center of more expansive areas of homogeneous forest as a solution to the disparity between the resolutions of the larger MODIS LST pixel and the T ust plots.Olsson and Jonsson [85] found that different scales of gridded temperature source-data (e.g., MODIS LST) had no significant influence on the accuracy of their forest phenology models.
For this mitigation approach to operate effectively, it was assumed that temperature across the broader homogenous stand was uniform.This assumption was validated over two seasons, by recording in situ air temperature within a 250 m plot placed at the center of a broad expanse (2700 km 2 ) of perfectly homogenous land cover (native grassland).Employing the identical methodologies as in the forested plots, these data were used to train the LST null model.Air temperature estimates in the validation plot had an average RMSE increase of 2 • C over the equivalent null model estimates in the forested plots.Theoretically, the error should have been lower since the validation plot had completely uniform land cover relative to the principal observation plots.However, this result substantiates that forest cover improves statistical estimates of air temperature using LST [22].
The tiled and projected M*D11A1 products were chosen for this study based on their ease of use and successful application in previous models (e.g., [6,46]).The lower-level LST swath data are less modified and unaffected by projection error or temporal adjustments.However, opting for the higher-level product was a tradeoff between these errors versus the benefit of having the obvious cloud-contaminated LST values removed.During pre-processing the LST imagery selected for this study was reprojected and resampled, which introduces geometric error: chiefly the potential to alter the LST pixel value with values from adjacent cover types.We attempted to mitigate this issue by working in forest stands that were much larger than the plots in which observations were made.Some error may have originated from intrinsic uncertainties in the radiometric or geometric precision of the original MODIS LST data [6], yet this was mitigated by using only the highest quality QA filtered pixels for the analysis.The local time correction made during pre-processing produced times that were truly instantaneous or within 30 min of a MODIS overpass.It is unlikely that this period was sufficient to allow any significant change in temperature on the ground, but it cannot be ruled out completely.
Terrain complexity within the study area introduced an assortment topographic and microclimatic factors, but these had little effect on the model results.This was expected since this study was a within-pixel analysis with no broader spatial interpolation.Observation plots with limited topography were intentionally selected as "targets" for MODIS in an attempt to constrain the variability of T ust solely to forest type and structure.Elevation was the only significant topographic variable in the models because small differences in elevation produce considerable changes in the environmental lapse rate (e.g., [7]).The use of aspirated solar radiation shields worked to limit in situ T ust measurement bias.However, variability in canopy cover between plots may have exposed some sensors to slightly more radiation than others, potentially altering the results [86].
Clear-sky overpasses occurred predominantly during the summer and autumn months, however a considerable proportion occurred in autumn.The canopy is mostly intact for the first half of September, but leaf senescence may have introduced some bias in the small proportion of plots that were deciduous.Given the enduring issue of cloud over the mountains, there was a tradeoff between increased observations during a reliably clear time of year and the transition in canopy attributes.An extended observation period was required to observe the full phenological progression of the understory and to obtain a range of seasonal temperatures.Other sources of model uncertainty most likely occurred during exceedingly cold times of year.The MODIS LST product is accurate to ±1 • C between the temperatures of −10 • C and 50 • C [9], there were early and late season temperatures during the study period that were approximately −20 • C potentially reducing model accuracy.This evidence of seasonal variability affecting model estimates is substantiated by Hanes and Schwartz [1] who found that differences between T air and MODIS LST are increasingly divergent during colder parts of the year in temperate climates, mainly due to snow cover.However, for the purposes of using T ust for monitoring vegetation phenology, the accuracy of T ust models is not critical during these colder periods since the plants are still dormant and GDD accumulations have not yet begun.
Overall, the models demonstrate that estimates of T ust using MODIS LST can be refined with additional variables which characterize the interface between the top of the forest canopy and the understory.These results facilitate the next step in this research to monitor understory phenology over the entire study area.This requires spatially extending the models beyond the observation plots across the gridded MODIS LST surface to produce maps of daily T ust .GDD maps derived from T ust should provide improved accuracy in predictions of understory plant phenology than those derived from canopy-top T air .Hanes and Schwartz [1] found that top-of-flux tower GDD accumulations underpredicted plant phenology opposed to measurements made closer to the ground since lapse rates produce cooler temperatures at height.The instantaneous T ust models derived in this study will ultimately facilitate the development of improved land surface phenology maps.

Conclusions
Fluctuations in surface air temperature drive a variety of environmental processes which include radiative energy and gas fluxes, species interactions, and vegetation growth.Recently there has been growing importance on collecting thermal data from satellite platforms to estimate T air at regional and global extents [16,87].However, the use of remote sensing for modeling forest understory temperatures has yet to be explored.In this paper we described a methodology for predicting plot-scale T ust using MODIS LST and a combination of conventional forest inventory variables and LiDAR-derived canopy metrics.The objective aimed to identify which of these ancillary characteristics played a critical role in modifying estimates of the underlying temperatures.T ust was predicted to within 1.4 • C across a variety of forest types and elevations.Canopy closure and the LiDAR-derived standard deviation of canopy height metric were found to significantly improve estimations of T ust over MODIS LST alone.These findings confirm that canopy structure and forest stand-type can function to differentiate understory air temperatures from ambient canopy temperature as observed by the sensor overhead.
This study developed a remote sensing procedure that works to uncover detailed microclimatic conditions beneath the forest canopy.This type of mechanistic data is particularly useful for ecological modeling [88], but is normally obtained from solitary weather stations or gridded top-of-canopy temperature surfaces.The T ust models were trained using in situ temperature data collected over a broader ground surface to approximate the pixel size of the MODIS imagery.This effort to align the spatial scale of observations likely improved the results and is recommended for any remote sensing investigations collecting surface temperature at moderate resolutions.
These models provide useful estimations of T ust at the plot-scale, but to fully benefit from the spatial nature of the remotely-sensed datasets it would be essential to develop broad-extent daily maps of T ust .There is a need for such a dataset, but foreseeable challenges include producing maps without spatial or temporal punctuations in the time series.Current trends in global climate necessitate predictions of how ecosystem processes will respond to these changes [89,90], and continuous maps of air temperature are instrumental in answering this question.Future research in this area promises to be rewarding as space-borne technologies continuously advance alongside ongoing LST product improvements (e.g., the recent release of MODIS collection version 6).These invaluable and freely available data will provide rapid assessments of the earth surface with perpetually increasing accuracy and precision.

Figure 1 .
Figure 1.Study area and understory temperature observation plots.

Figure 1 .
Figure 1.Study area and understory temperature observation plots.

Figure 1 .
Figure 1.Study area and understory temperature observation plots.

Figure 3 .
Figure 3. Simplified flow chart of the understory temperature modeling process.

Figure 3 .
Figure 3. Simplified flow chart of the understory temperature modeling process.

Figure 4 .
Figure 4. Example observation plot sample transects of tree species composition prism sweep sites (percent confer) and hemispherical photography locations (percent closure).

Figure 4 .
Figure 4. Example observation plot sample transects of tree species composition prism sweep sites (percent confer) and hemispherical photography locations (percent closure).

Figure 5 .
Figure 5. Spatial and temporal characteristics of the valid Moderate Resolution Imaging Spectroradiometer (MODIS) overpass dataset including spatial distribution of total clear-sky MODIS overpasses for each observation plot.Also, the proportional percentage of total clear-sky overpasses by month, and the proportional percentage of valid clear-sky observations versus all potential MODIS observations (cloud contamination) for the duration of the study period.

Figure 5 .
Figure 5. Spatial and temporal characteristics of the valid Moderate Resolution Imaging Spectroradiometer (MODIS) overpass dataset including spatial distribution of total clear-sky MODIS overpasses for each observation plot.Also, the proportional percentage of total clear-sky overpasses by month, and the proportional percentage of valid clear-sky observations versus all potential MODIS observations (cloud contamination) for the duration of the study period.

Figure 6 .
Figure 6.In situ understory temperatures vs. predicted Tust values estimated using the top-ranked models for (a) day and night; (b) daytime; and (c) night.These validation data were collected during clear-sky MODIS overpasses within a nearly-deciduous (6% conifer) observation plot with moderate canopy closure (60%) at an elevation of 1330 m.

Figure 6 .
Figure 6.In situ understory temperatures vs. predicted T ust values estimated using the top-ranked models for (a) day and night; (b) daytime; and (c) night.These validation data were collected during clear-sky MODIS overpasses within a nearly-deciduous (6% conifer) observation plot with moderate canopy closure (60%) at an elevation of 1330 m.

Figure 7 .
Figure 7.In situ understory temperatures vs. predicted Tust values estimated using the top-ranked models for (a) day and night; (b) daytime; and (c) night.Data are the randomly withheld clear-sky MODIS overpasses of the study area over the duration of the study period.

Figure 7 .
Figure 7.In situ understory temperatures vs. predicted T ust values estimated using the top-ranked models for (a) day and night; (b) daytime; and (c) night.Data are the randomly withheld clear-sky MODIS overpasses of the study area over the duration of the study period.

Table 2 .
Conventional explanatory variables used for modeling understory temperature.

Table 2 .
Conventional explanatory variables used for modeling understory temperature.

Table 3 .
Baseline Akaike information criterion (AIC) model rankings (top 10) showing relative performance of conventional and LiDAR-derived variable types, as well as individual variable contributions in estimating T ust using all clear-sky (day-night) observations.Model rank was assessed through difference in AIC values (∆i) and weights (w i ) reflecting model likelihood.Model complexity is characterized by the number of model parameters (K i ).

Table 4 .
AIC ranking of the top day-night, night, and daytime models within the candidate groupings of conventional, LiDAR-derived, and the combined metrics to predict T ust .

Table 5 .
Estimated parameters for the top-ranked models from each overpass dataset.Standardized coefficients and standard errors (in parentheses) are presented by model.

Table 6 .
Comprehensive results of the top-ranked Tust models.Coefficient of determination (R 2 ), mean absolute error (MAE), root mean square error (RMSE).

Table 6 .
Comprehensive results of the top-ranked T ust models.Coefficient of determination (R 2 ), mean absolute error (MAE), root mean square error (RMSE).