Surface Energy Flux Estimation in Two Boreal Settings in Alaska Using a Thermal-Based Remote Sensing Model

: Recent Arctic warming has led to changes in the hydrological cycle. Circum-Arctic and circumboreal ecosystems are showing evidence of “greening” and “browning” due to temperature warming leading to shrub encroachment, tree mortality and deciduousness. Increases in latent heat ﬂux from increased evapotranspiration rates associated with deciduous-dominated ecosystems may be signiﬁcant, because deciduous vegetation has extremely high-water use and water storage capacity compared to coniferous and herbaceous plant species. Thus, the impact of vegetation change in boreal ecosystems on regional surface energy balance is a signiﬁcant knowledge gap that must be addressed to better understand observed trends in water use / availability and tree mortality. To this end, output from a two-source energy balance model (TSEB) with modiﬁcations for high latitude boreal ecosystems was evaluated using ﬂux tower measurements and Terra / Aqua MODIS remote sensing data collected over the two largest boreal forest types in Alaska (birch and black spruce). Data under clear and overcast days and from leaf-out to senescence from 2012 to 2016 were used for validation. Using ﬂux tower observations and local model inputs, modiﬁcations to the model formulation for soil heat ﬂux, net radiation partitioning, and canopy transpiration were required for the boreal forest. These improvements resulted in a mean absolute percent di ﬀ erence of around 23% for turbulent daytime ﬂuxes when surface temperature from the ﬂux towers was used, similar to errors reported in other studies conducted in warmer climates. Results when surface temperature from Terra / Aqua MODIS estimates were used as model input suggested that these model improvements are pertinent for regional scale applications. Vegetation indices and LAI time-series from the Terra / Aqua MODIS products were conﬁrmed to be appropriate for energy ﬂux estimation in the boreal forest to describe vegetation properties (LAI and green fraction) when ﬁeld observations are not available. Model improvements for boreal settings identiﬁed in this study will be implemented operationally over North America to map surface energy ﬂuxes at regional scales using long time series of remote sensing estimates as part of NOAA’s GOES Evapotranspiration and Drought Information System.


Introduction
In the recent past, a significant increase in air temperature in the Arctic [1], also known as the Arctic amplification, is leading to an accelerated increase in air temperatures compared to other parts of the globe inducing broad sea-ice retreat, snow and ice melting and increases in sea level [2,3]. Arctic and sub-Arctic ecosystems in Alaska, dominated by tundra and boreal forest land covers, are witnessing singular changes due to climate warming including widespread permafrost degradation, increases in the area burned and the severity of wildfires, decreased thickness and duration of winter snow cover, acceleration of the hydrological cycle and, rises in river discharge and important vegetation changes in both structure and distribution [4][5][6][7][8][9][10][11]. Furthermore, land surface Arctic hydrologic feedbacks to changing climate are actively coupled to the energy balance of these ecosystems [12]; and in turn, the partitioning of this energy balance plays a critical role in the hydrologic cycle regulation [13].
Circum-Arctic and circumboreal ecosystems are showing evidence of "greening" and "browning" [14][15][16][17], respectively, especially in Alaska and Western Canada. Due to increased temperatures there has been a rise in deciduous shrubs in the Circum-Arctic, which is responsible for the "greening" of the Arctic tundra [16,18] resulting in nearly a 14% increase in vegetation cover [4]. On the other hand, some boreal forest regions are showing "browning" caused by higher coniferous tree mortality because of warmer and drier conditions. Moreover, a higher frequency and intensity of wildfires is changing the boreal forest composition, increasing deciduous tree and shrub cover composition due to their higher drought tolerance and their ability to establish quickly after disturbance [15,[19][20][21]. Hence, the transition zone between forest and tundra ecosystems in the northern boundary is expanding both latitudinally and in elevation and is leading to an increase of tree heights and shrub growth leading to denser and taller canopies [16,22].
The impact of vegetation change in boreal ecosystems on the regional surface energy balance is a significant knowledge gap that must be addressed. Filling this knowledge gap is critical for quantifying feedbacks to, and responses and vulnerability of, this landscape to continued climate warming, shifts in hydrology, and increased disturbance from wildfires [13,23]. Increases in latent heat flux associated with deciduous dominated ecosystems may be significant because deciduous vegetation have extremely high water use and storage compared to native coniferous and herbaceous plant species [24,25].
Currently, networks such as FLUXNET allow quantifying local surface energy fluxes by collecting data from different ecosystems worldwide, including the Arctic. However, measurements taken from the available flux towers are not necessarily representative of the variation across the surrounding landscape. Additionally, in Arctic and boreal regions there are very few FLUXNET towers from 56 • N to 71 • N compared to other latitudes [26]. Consequently, the flux network in high latitudes cannot capture the critical changes in vegetation and surface energy balance occurring with climate warming. Finally, with a significantly lower spatial distribution of flux towers in the Arctic and boreal regions compared to mid-latitude regions, combined with their remoteness, the harsh environment and the maintenance and travel costs, data from existing towers have significant temporal gaps [27].
Due to the lack of the spatiotemporal information on the water and energy balance exchange in boreal ecosystems and the role of vegetation changes in driving these budgets, there is a critical need to apply remote sensing-based energy balance models to quantify these impacts across the landscape. Therefore, this study focuses on refining and evaluating a diagnostic remote sensing-based energy balance model for estimating seasonal dynamics of surface energy fluxes over the boreal forest, using measurements of land surface temperature retrieved from thermal infrared sensors on satellite platforms as a key boundary condition. Specifically, an improved version of the Two-Source Energy Balance model (TSEB, Norman et al. [28]), not yet examined for high latitude boreal ecosystems, is evaluated with two eddy covariance flux towers over birch and black spruce, which are vegetation types representative of the boreal forest in interior Alaska. The model is applied for all-sky conditions from 2012 to 2016 using critical inputs of surface temperature derived from (a) local flux tower observations and (b) satellite remote sensing estimates. Remote sensing estimates of vegetation properties (leaf area index (LAI), NDVI and, EVI) and surface temperature from the Aqua and Terra MODIS are used as inputs in TSEB to evaluate further implementation for a regional Atmosphere-Land Exchange Inverse (ALEXI) modelling system [29], implemented operationally over North America as part of NOAA's GOES Evapotranspiration and Drought Information System [30].

Overview of the Two-Source Energy Balance (TSEB) Model
To estimate surface energy fluxes, the Two-Source Energy Balance (TSEB) model in its series version [28,31,32] was used with later modifications for Arctic tundra environments and applied in clear sky and cloudy conditions [27]. Table 1 lists all TSEB modifications for the boreal forest applied and evaluated in this study. The model considers the surface radiometric temperature (T RAD ), either measured from field-based radiometers or from satellite thermal sensors, to be a composite of both soil (T S ) and canopy (T C ) temperatures weighted by the fraction of vegetation (f C ) observed at a certain radiometer or thermal sensor viewing angle (ω): Estimation of f C (ω) in canopies with a spherical leaf angle distribution can be approximated by: where LAI is the leaf area index and Ω is a clumping factor that considers the degree to which the vegetation is non-randomly distributed, for example as in a row crop or in a clumped sparsely vegetated canopy. The surface energy balance equation for both the canopy ( C ) and soil ( S ) components of the combined soil-canopy-atmosphere system is formulated as: where R N is the net radiation for both the soil and canopy components estimated considering the short-wave and long-wave radiation divergence. Due to a high cloud cover presence in Arctic and boreal ecosystems, atmospheric emissivity (ε a ) in the net longwave radiation configuration needs to be modified when applied for all-sky conditions (including both clear sky and overcast conditions) [33][34][35] as follows: where s is the solar radiation to the potential solar radiation ratio, e is vapor pressure and T A is air temperature.
In the TSEB formulation [28], the total R N and both the sensible (H) and latent heat (LE) fluxes are summed from the soil and canopy component fluxes: A series-resistance version of TSEB allow interaction between soil and canopy sensible heat fluxes and temperature differences as follows: Here, T AC is the air temperature in the canopy-air space (K), R S is the resistance to the heat flow in the boundary layer above the soil (s·m −1 ), R X is the total boundary layer resistance of the vegetation leaf layer (s·m −1 ), and R A is the aerodynamic resistance computed from the stability-corrected temperature profile equations (s·m −1 ). Further information on the resistance formulation is given in Appendix A.
The canopy latent heat (LE C ) is estimated following the Priestley-Taylor formula [36] initially assuming a potential rate for LE C, and soil latent heat (LEs) is solved as a residual as follows: where f G is green vegetation fraction (-), γ is the psychrometric constant (kPa K −1 ), ∆ is the slope of the saturation vapor pressure versus temperature (kPa K −1 ), α PTC is the Priestley-Taylor coefficient initially determined for the canopy component with a value of 1.26 for general conditions tested during the growing season in rangelands and croplands. However, TSEB internally modifies its initial value under conditions of canopy stress (see Section 2.2.1 for further details). Finally, in the TSEB scheme, G is computed as portion of R NS , accounting for a daily cycle through a phase shift between G and R NS suggested by [37] following: Here, c G is the maximum value of the G/R NS ratio that, based on experimental data for several conditions in rangelands and croplands, can be assumed to be around~0.3 [38], although it may vary depending on soil class and soil moisture. However, for boreal forest applications, the c G value was modified to 0.07 (see Section 2.2.2 for further details). B is chosen to minimize the deviation of c G from Equation (14), t is time in seconds relative to solar noon, and S is the phase shift between G and R NS in seconds. To estimate the LE flux from the canopy component, TSEB internally changes its initial α PTC values of 1.26 to allow an acceptable partitioning between LE C and LE S under stressed vegetation conditions. The initial value of α PTC is down-adjusted when TSEB results in negative values for LE S , given that condensation on the soil is unlikely to occur during the day [28]. This condition typically occurs if the potential transpiration flux is too high to be consistent with the observed surface temperature. The canopy temperature estimate in is too low, resulting in a high soil temperature derived from Equation (1). This leads to high soil sensible heat (Equation (10)) exceeding the available energy at the soil surface, and hence LE S < 0 result from Equation (13). When this non-physical condition is encountered midday, it is assumed that the canopy is stressed and α PTC is iteratively reduced until LE S values are higher than 0. This iterative scheme works well in ecosystems where canopy values of α PTC are relatively conservative under unstressed conditions [39], while for stressed canopies where the soil surface is usually dry, a LE S value close to zero is a reasonable assumption [40].
In croplands, rangelands and temperate forests where the TSEB model was previously evaluated, the assumptions that (a) the vegetation is fully green (f G = 1); and (b) that the maximum potential (unstressed) transpiration rate is obtained with α PTC = 1.26 [36]; yielded successful results. However, Ref. [41] suggested that for canopies that are more conservative in their water use (not transpiring at the potential rate) or not completely green, α PTC and f G need to be adjusted to yield a reasonable LE partitioning. Moreover, refs. [27,42] and also found acceptable TSEB results for natural vegetation and Arctic tundra using initially a lower α PTC , suggesting that lower values for natural ecosystems with vegetation adapted to water-limited environments would be more realistic and yield more accurate LE values.
In many boreal forests, tree growth is limited by low precipitation and low temperatures that in turn restrict photosynthetic capacity and reduce root hydraulic conductivity and stomatal conductance, resulting in low-LAI canopies that exert a significant resistance to transpiration [43]. Moreover, conifer forests (such as black spruce forests) growing in upland regions of the boreal zone, evaporate at rates between 25 and 75% of equilibrium evaporation defined as a special case when α PTC = 1. On the other hand, evaporation rates from deciduous forests (such as birch forests) may approach equilibrium rates [43]. The standard TSEB parameterization of α PTC , assuming an initial value of 1.26, may not reach the low values resulting the controls of climate and vegetation on the energy exchange of boreal forests.
Although specific values of α PTC for the boreal forest are not found in the literature, α PT measurements for the whole system (canopy and soil) and hereinafter referred to as α PTS, are available for deciduous forest and boreal conifer evergreen systems [43][44][45][46][47][48][49][50][51] with an average value of 0.6 ± 0.17 and 0.9 ± 0.2, respectively, and representative for black spruce and birch forests. Moreover, α PTS also may show seasonal variations that can vary significantly with LAI, vapour pressure deficit (VPD), and soil moisture [52].
For modelling purposes, deciduous forest and boreal conifer evergreen averaged α PTS values together with the original TSEB α PTC value of 1.26 were used to estimate the most appropriate initial α PTC value for TSEB. Moreover, to capture seasonal phenology in green vegetation cover, a simple vegetation index ratio proposed by [43] was applied: Additionally, to evaluate the effect of α PTC and f G on the latent heat estimation, a sensitivity analysis on α PTC and f G was conducted [42]. This consists in running TSEB with α PTC values increasing at 0.05 intervals from 0.4 to 1.3 for both black spruce and birch, respectively, with (a) the f G derived from the EVI and NDVI ratio and (b) a f G = 1 assuming active transpiration from the canopy .
Finally, α PTS measurements were computed to (1) analyze its seasonal behavior at both flux towers as a function of soil moisture, VPD and phenology, and (2) compare its value with results from the sensitivity analysis. Values of α PTS were computed using the flux tower data using the following equation: where E and E eq are the measured evaporation from the flux towers and the equilibrium evaporation rate, respectively, with E eq computed as follows:

Modifications and Evaluation for Soil Heat Flux
For soil heat flux estimation, a c G mean value of 0.3 is not likely appropriate for boreal forest settings due to forest floor residue and subcanopy moss and vegetation For tundra vegetation, large errors using a c G of 0.3 were also found [27]. Although there is a lack of specific c G values for spruce or deciduous forests in the literature, data derived from different studies in summer months and in similar boreal settings [43,44,47,50,[53][54][55][56][57] allowed computing a mean value of c G for spruce and deciduous forests of 0.07 ± 0.03.
Although this approach was successfully applied for croplands, a new approach in Arctic tundra [27] to estimate G using the relationship between G and T RAD in Equation (18) yielded better agreement than using both constant or phase shifted c G values in Equation (14). This methodology also accounted for a phase shift by using the maximum value of c GT in Equation (14) instead of c G as follows: Thus, to estimate G, a modified c G value for the boreal forest of 0.07 was applied to the original TSEB formulation (Equation (14)). Moreover, the approach using the G-T RAD relationship (Equations (14) and (18)) was also fitted and evaluated with 60% and 40% of the flux tower dataset, respectively, and then applied to the full dataset at the black spruce and birch forest sites.

Study Area and Instrumentation
Two experimental sites were established from 2011 to 2012 in two of the main covers of the boreal forest, birch and black spruce forests, for the TSEB calibration and evaluation in boreal settings ( Figure 1). The first setting, located at the North Campus of the University of Alaska Fairbanks (UAF; 64 • 51 56.803"N and 147 • 50 34.154"W), was installed in a needle-leaf black spruce forest with discontinuous permafrost. The overstory was predominantly black spruce (Picea mariana) with a canopy cover around 60% and a mean tree height of about 5 m. The understory is covered by shrubs (Betula nana, Ledum palustre, Alnus incana) and mosses (Sphagnum spp.). The flux tower was installed in 2011 with a total height of 23 m. a.g.l. Unfortunately, in 2014, the thaw of discontinuous permafrost around the flux towers caused tower integrity issues and tower height was reduced to 15 m. This led to a data gap of around 1 year, from July 2014 to June 2015.
The second setting was located at the Caribou Poker Creeks Research Watershed (CPCRW; 65 • 10 17.962"N and 147 • 28 17.137"W) in a deciduous paper birch (Betula neoalaskana) forest with a summer canopy cover around 95% and a mean tree height of about 16 m and a tower height of 23 m a.g.l. In this case, the understory is covered by shrubs (Betula nana, Ledum palustre) and mosses (Sphagnum spp.).Both field sites were equipped with a sonic anemometer, a gas analyzer operating at 20 Hz sampling rate, four-component net radiometer sensor and air temperature sensors at different heights (for further details on tower instrumentation see [58]). Ground heat was monitored at both sites by temperature and soil moisture probes and heat flux plates installed in the subsurface soil layers. Precipitation was measured at black spruce flux tower using a rainfall gauge. At the birch flux tower rainfall data from another meteorological station operated by the CPCRW long-term ecological research (LTER) network were used. Both towers were operated year-round, although some data gaps were present in deep winter due to power shortage or flux tower access limitations due to harsh weather conditions.

Model Input and Evaluation
The TSEB model was run in two modes: (a) using local measurements of surface temperature as input, and (b) using remote sensing of surface temperature estimates, and in both cases with meteorological data from the flux towers (air temperature, atmospheric pressure, wind speed, vapor pressure, and solar radiation). For both modes, vegetation properties such as fG and LAI were derived using remote sensing estimates. As a local source of surface temperature, TRAD, upwelling longwave data from the four component net radiation sensor were converted to TRAD [58]. For model runs using satellite-based TRAD, Terra/Aqua MODIS TRAD from Terra/Aqua MODIS LST product (MOD11) was used. To ensure the best quality data, only high quality LST data with an average error ≤ 1 K and with a view-angle up to 35 degrees according to quality flags of the product were considered. Pixel selection was done by a footprint analysis based on [59] methodology accounting for 90% of the flux tower footprint cover.

Model Input and Evaluation
The TSEB model was run in two modes: (a) using local measurements of surface temperature as input, and (b) using remote sensing of surface temperature estimates, and in both cases with meteorological data from the flux towers (air temperature, atmospheric pressure, wind speed, vapor pressure, and solar radiation). For both modes, vegetation properties such as f G and LAI were derived using remote sensing estimates. As a local source of surface temperature, T RAD , upwelling longwave data from the four component net radiation sensor were converted to T RAD [58]. For model runs using satellite-based T RAD , Terra/Aqua MODIS T RAD from Terra/Aqua MODIS LST product (MOD11) was used. To ensure the best quality data, only high quality LST data with an average error ≤1 K and with a view-angle up to 35 degrees according to quality flags of the product were considered. Pixel selection was done by a footprint analysis based on [59] methodology accounting for 90% of the flux tower footprint cover.
To evaluate R N , LE, H, and G TSEB model output, a conservative approach was taken to screen out the 30-min time-step flux tower data to ensure inclusion of only high-quality data. To cover the whole growing season at both sites, only snow-free data between May to September from 2012 to 2016 data were considered, identified using Terra/Aqua MODIS snow cover products. Then, daily rainfall events were also removed. Finally, two more filters to ensure data quality and daytime conditions were applied: (a) a 30-min timescale energy closure > 70% and, (b) R N > 100 W·m −2 . This resulted in a total of 289 days at UAF and 392 at CPCRW. It is important to note that the final evaluation dataset was substantially reduced due to the large amount of rainy days, which are common in the boreal forest.

Micrometeorological Data Processing
High-frequency (20 Hz) eddy covariance raw data collected at both sites were screened to identify and remove nonphysical values and data spikes [60]. Half-hour mean wind speed and direction were then computed after wind velocity components of the coordinate system were rotated to align with prevailing wind direction [61,62]. Sonic temperature data was then converted to air temperature by adjusting for humidity effects [63] and corrected for sensor displacement and frequency response attenuation [64,65]. The moisture and carbon dioxide fluxes were then corrected for the effects of buoyancy and water vapor density [66] and as a result, 30-min turbulent fluxes were calculated. The daily average energy balance closure for the selected period at both black spruce and birch sites was 0.87 and 0.91, respectively, which is in agreement with many other eddy covariance studies [67]. To apply the surface energy balance expression, residuals were allocated to latent heat flux [68], ensuring the closure of the energy balance and further comparability with modelled fluxes. Corrections to account for soil heat storage were applied on the soil heat flux plates measurements together with soil bulk density sampled at each site (1244 and 1208 kg·m −3 for birch and black spruce, respectively) as well as soil moisture and soil temperature probes [69,70].

Remote Sensing Estimates of Vegetation Properties
The TSEB configuration used in this study requires estimates of LAI and f G derived from EVI and NDVI. LAI field measurements were obtained intermittently at both sites during the study time period with a Decagon LP-80 Ceptometer using 50 regular samples within a 200 m radius nearby the flux towers. However, in-situ measurements did not cover the growing season to senescence, and field measurements of NDVI or EVI were unavailable. To produce temporally smoothed NDVI, EVI and LAI data, TIMESAT [71] was used together with MODIS 4-day LAI product (MCD15A3H), both MODIS daily reflectance product at 250 and 500 m (MOD09GQ and MOD09GA, respectively) and their quality flags (see [27] for more details). The clumping factor Ω (Equation (2)) was set to 0.8 and 0.7 for birch and black spruce sites, respectively, based on forest inventories and image photointerpretation.

Results and Discussion
In this section, we present an evaluation of daytime instantaneous surface energy flux estimates from TSEB over all-sky conditions using local and remotely sensed inputs of T RAD , and applying different values of α PTC and f G .

Model Performance Using In-Situ T RAD Measurements
For both flux datasets used here, measured G was a relatively small term with a daytime average value between 9 and 17 and W·m −2 for the black spruce and birch sites, respectively, in comparison to daytime average R N of around 321 W·m −2 at both sites. Fitting Equation (14) for c G for both spruce and deciduous sites resulted in values for c G , B and S of 0.07, −7200 and 250,000, respectively. When the variable c TG computed from Equation (18) replaced c G in Equation (14), fitted values for c TG, B and S were 0.9, −7200 and 200,000, respectively. In both cases, a 2 h phase shift after the maximum T RAD at noon showed a negligible influence on the results when applying a B variation of ±15,000 s.
Using c G = 0.07 in Equation (14) for G estimation, model estimates showed lower agreement with observed soil heat fluxes at the black spruce site, with R 2 , the mean absolute percent difference (MAPD), the root mean square error (RMSE) and mean bias error (MBE) values of 0.1, 95%, 19 W·m −2 and 17 W·m −2 , respectively, and with R 2 , MAPD, RMSE and MBE values of 0.05, 98%, 26 W·m −2 and 24 W·m −2 for the birch site (results not shown in Table 2). When G was estimated replacing c G in Equation (14) with c GT estimated from Equation (18) a better agreement was obtained, reducing by 50% discrepancies with measured G compared to the standard TSEB approach using a constant c G ( Table 2). This is also in agreement with [27] who found better performance using the G-T RAD relationship for Arctic tundra. In this case, modeled G had negligible bias compared to observations, with R 2 , MAPD, RMSE and MBE values of 0.47, 44%, 5 W·m −2 and 1 W·m −2 , the for black spruce site, and with R 2 , MAPD, RMSE and MBE values of 0.65, 47% 7 W·m −2 and −3 W·m −2 , for the birch site ( Figure 2 and Table 2). Although RMSE is similar at both sites, the lower observed soil heat flux at the black spruce site, due a thicker moss layer that increases soil insulation, led to a higher MAPD value for the birch site.
Modeled R N yielded comparable results at both sites with a high correlation and almost negligible bias, with values of MAPD of around 5%, an RMSE from 18 W·m −2 to 22 W·m −2 , an MBE from 0 W·m −2 to 4 W·m −2 and an R 2 of 0.98 ( Figure 2 and Table 2). These results are similar to those found in Arctic tundra using the same R N modelling configuration [27]. This suggests that this model configuration may be also applied for the whole growing season in all-sky conditions for boreal forest settings. Moreover, they are also aligned with previous TSEB model findings in different cover types and for clear sky in which a 5% of error percentage (MAPD) was described [31,[39][40][41][72][73][74]. Table 2. Results of the model agreement and error estimation compared to observed flux tower data using an initial TSEB α PTC value of 1.26 at both settings, 0.6 for black spruce setting, and 0.9 for birch setting. n is the number of 30-min periods selected. Error estimators (RMSE, mean absolute difference (MAD), MBE) are in W·m −2 while MAPD is in %. G was modelled using c GT in Equation (18) applied to Equation (14). X are the mean values for the estimated energy balance components in W·m −2. Turbulent heat fluxes, H and LE, estimated with model modifications in G and R N , yielded reasonable agreement with 30-min observed fluxes when α PTC was adjusted for boreal vegetation (0.6 and 0.9 for black spruce and birch settings, respectively), and poorer agreement using the original α PTC of 1.26, with an RMSE difference between both configurations of 23 W·m −2 to 28 W·m −2 (Table 2 and Figure 2).

Black Spruce|α
For birch and black spruce, RMSE for LE and H for the original α PTC configuration ranged from 64 to 77 W·m −2 , while for the adjusted α PTC configuration LE and H ranged from 41 to 49 W·m −2 , errors comparable to those described in other studies [75], and according to [76] after daily integration of instantaneous daytime fluxes the errors will tend to be reduced on the order of 10-15%. The error percentage (MAPD) in LE and H using the adjusted α PTC values were also reduced from 32% to 39% to less than 25%. For both original and adjusted α PTC values, the model tended to overestimate LE and underestimate H. Model LE and H under-and overestimation may be also explained due to the lack of instrument closure and methodological uncertainties, insufficient estimation of storage terms, unmeasured advective fluxes, landscape scale heterogeneity or instrument spatial representativeness, among others [67,68,70,77,78].
In comparison, Sanchez et al. [79] reported similar RMSE results of around 50 W·m −2 for both H and LE in a Finnish boreal forest over a two-month summer validation study with a simplified version the TSEB. These results suggest that the adjusted lower α PTC values for deciduous and coniferous forest covers better describe controls of climate and vegetation on the energy exchange of boreal forest under current conditions. This is in line with other studies in deciduous and coniferous forests [27,41,42,49] suggesting that vegetation-type adjusted values of α PTC could yield more accurate H and LE values in cases where the natural vegetation is adapted to the local climate conditions. Although a vegetation class-dependent value of α PTC might seem disadvantageous as ancillary land use cover or vegetation cover information is needed, current initiatives in mapping the boreal forest vegetation types [80,81] should enable use of adjusted values of α PTC , thus, decreasing the bias associated with using a universal value of α PTC of 1.26 for regional applications.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 24 lack of instrument closure and methodological uncertainties, insufficient estimation of storage terms, unmeasured advective fluxes, landscape scale heterogeneity or instrument spatial representativeness, among others [67,68,70,77,78]. In comparison, Sanchez et al. [79] reported similar RMSE results of around 50 W·m −2 for both H and LE in a Finnish boreal forest over a two-month summer validation study with a simplified version the TSEB. These results suggest that the adjusted lower αPTC values for deciduous and coniferous forest covers better describe controls of climate and vegetation on the energy exchange of boreal forest under current conditions. This is in line with other studies in deciduous and coniferous forests [27,41,42,49] suggesting that vegetation-type adjusted values of αPTC could yield more accurate H and LE values in cases where the natural vegetation is adapted to the local climate conditions. Although a vegetation class-dependent value of αPTC might seem disadvantageous as ancillary land use cover or vegetation cover information is needed, current initiatives in mapping the boreal forest vegetation types [80,81] should enable use of adjusted values of αPTC, thus, decreasing the bias associated with using a universal value of αPTC of 1.26 for regional applications.

Evaluation of Remote Sensing Vegetation Properties Estimates
Time series of LAI from the MODIS product used to partition soil and canopy temperatures (Equation (2)) captured seasonal dynamics from green-up to senescence for the whole period at both flux towers (Figure 3). Moreover, comparison with in-situ field estimates of LAI (not originally intended for LAI MODIS product evaluation) yielded a reasonable error with a RMSE of 0.7 and 0.6 and an MBE of −0.7 and 0.2 at the black spruce and birch sites, respectively. For black spruce, fG computed using Equation (15) a good correspondence with LAI dynamics for the whole growing

Evaluation of Remote Sensing Vegetation Properties Estimates
Time series of LAI from the MODIS product used to partition soil and canopy temperatures (Equation (2)) captured seasonal dynamics from green-up to senescence for the whole period at both flux towers (Figure 3). Moreover, comparison with in-situ field estimates of LAI (not originally intended for LAI MODIS product evaluation) yielded a reasonable error with a RMSE of 0.7 and 0.6 and an MBE of −0.7 and 0.2 at the black spruce and birch sites, respectively. For black spruce, f G computed using Equation (15) a good correspondence with LAI dynamics for the whole growing period (May to September) was found. The birch site showed a similar behavior from June to August. However, in May and September the pattern was different due to an already green understory with low LAI in leaf-out and senescence (Figure 3

Seasonal Dynamics in Surface Energy Fluxes
Monthly estimation of RN, LE, H, and G with adjusted αPTC and the modified G configuration showed an acceptable overall agreement with observations from June through August (months with the greatest vegetation activity in the boreal forest), with RMSE values lower than 50 W·m −2 and an mean MAPD around 23% (Tables 3 and 4, and Figure 4). However, model performance deteriorated to some extent at leaf-out and at senescence. During these periods, LE was overestimated by the model, and H underestimated. This may be related to unreliable estimates of fG or αPTC at the start of leaf-out and at the end of the senescence periods.

Seasonal Dynamics in Surface Energy Fluxes
Monthly estimation of R N , LE, H, and G with adjusted α PTC and the modified G configuration showed an acceptable overall agreement with observations from June through August (months with the greatest vegetation activity in the boreal forest), with RMSE values lower than 50 W·m −2 and an mean MAPD around 23% (Tables 3 and 4, and Figure 4). However, model performance deteriorated to some extent at leaf-out and at senescence. During these periods, LE was overestimated by the model, and H underestimated. This may be related to unreliable estimates of f G or α PTC at the start of leaf-out and at the end of the senescence periods.    Due to the lack of α PTC values in the literature, a first approach to understand its seasonal behavior is assuming α PTS computed in Equations (16) and (17) with flux tower data as representative of the α PTC behavior. Results showed that from June to August α PTS values were similar to other studies and used in this study to model turbulent fluxes with a value of 0.89 ± 0.05 for birch and 0.56 ± 0.03 for black spruce ( Figure 5). However, in the leaf-out and senescence (May and September) α PTS yielded lower values. For black spruce, values were 0.49 and 0.56, similar to the reference α PTC value of 0.6 used to estimate LE. However, for birch values were almost half of the reference α PTC value of 0.9 that was used, 0.47 and 0.53, respectively, similar to those reported in other boreal deciduous forests [82]. In both cases, α PTS behaved according to VPD. Moreover, results from the sensitivity analysis (Figures 6 and 7) for αPTC showed similar temporal behavior to that from αPTS values ( Figure 5). For black spruce, the lowest error (RMSE) in LE estimation was found with an αPTC value of around 0.55 together with a variable fG using the EVI and NDVI relationship (Equation (15)). The RMSE ranged from 40 W·m −2 to 50 W·m −2 for the whole growing season. Using a constant fG value of 1.0, the error increased for all LE estimations for the whole period suggesting that a reduction in αPTC alone was not enough to properly capture LE dynamics.
For birch, lower errors were found from June to August when a αPTC value of 0.8 and a variable fG were applied ( Figure 6). However, in May and September, there was an improvement in LE when lower αPTC values of around 0.5 and a variable fG were applied of around 15 W·m −2 and 2 W·m −2 , respectively. This behavior is similar to αPTS for the whole system (soil + canopy) seasonal dynamics ( Figure 5). When fG was considered to be 1, the overall error was higher. Moreover, lower values of αPTC of around 0.6 were needed to obtain similar values to those using a variable fG. This led to more unrealistic αPTS values compared to those found in Figure 5 and reported in the literature.
To further evaluate the sensitivity of choice of αPTC on LE estimation, differences in MAPD obtained using a range in αPTC from that obtained using the reference values of 0.6 and 0.9 for black spruce and birch, respectively, are shown in Figure 7. MAPD for black spruce was almost insensitive to use of lower values of αPTC, and errors increase steadily for αPTC above the reference value. This suggests that an initial αPTC value of 0.6 should be used for modeling black spruce water use. For all values of αPTC, applying a variable fG improved the LE estimation as already reported by [41]. Birch showed similar pattern from June to August, with an αPTC of around 0.8. However, the error difference within a αPTC interval from 0.8 to 1 in absolute value was around 1%, steadily increasing when lower and higher values were applied. This suggests that an initial αPTC mean value of 0.9 found in the literature may be applied to model LE regionally for birch. However, in May, αPTC lower values of 0.5 yielded a 7% improvement while in September of around 2%, suggesting that a value of 0.5 would be more appropriate to estimate LE. This moderate improvement in September could be related to the issues in capturing fG dynamics at the end of the season due to an already green understory but with low evaporation rates, although more research is needed. Moreover, results from the sensitivity analysis (Figures 6 and 7) for α PTC showed similar temporal behavior to that from α PTS values ( Figure 5). For black spruce, the lowest error (RMSE) in LE estimation was found with an α PTC value of around 0.55 together with a variable f G using the EVI and NDVI relationship (Equation (15)). The RMSE ranged from 40 W·m −2 to 50 W·m −2 for the whole growing season. Using a constant f G value of 1.0, the error increased for all LE estimations for the whole period suggesting that a reduction in α PTC alone was not enough to properly capture LE dynamics.
For birch, lower errors were found from June to August when a α PTC value of 0.8 and a variable f G were applied ( Figure 6). However, in May and September, there was an improvement in LE when lower α PTC values of around 0.5 and a variable f G were applied of around 15 W·m −2 and 2 W·m −2 , respectively. This behavior is similar to α PTS for the whole system (soil + canopy) seasonal dynamics ( Figure 5). When f G was considered to be 1, the overall error was higher. Moreover, lower values of α PTC of around 0.6 were needed to obtain similar values to those using a variable f G . This led to more unrealistic α PTS values compared to those found in Figure 5 and reported in the literature.
To further evaluate the sensitivity of choice of α PTC on LE estimation, differences in MAPD obtained using a range in α PTC from that obtained using the reference values of 0.6 and 0.9 for black spruce and birch, respectively, are shown in Figure 7. MAPD for black spruce was almost insensitive to use of lower values of α PTC , and errors increase steadily for α PTC above the reference value. This suggests that an initial α PTC value of 0.6 should be used for modeling black spruce water use. For all values of α PTC , applying a variable f G improved the LE estimation as already reported by [41]. Birch showed similar pattern from June to August, with an α PTC of around 0.8. However, the error difference within a α PTC interval from 0.8 to 1 in absolute value was around 1%, steadily increasing when lower and higher values were applied. This suggests that an initial α PTC mean value of 0.9 found in the literature may be applied to model LE regionally for birch. However, in May, α PTC lower values of 0.5 yielded a 7% improvement while in September of around 2%, suggesting that a value of 0.5 would be more appropriate to estimate LE. This moderate improvement in September could be related to the issues in capturing f G dynamics at the end of the season due to an already green understory but with low evaporation rates, although more research is needed.      Model performance in estimating R N was comparable at both sites with a high correlation (R 2 from 0.95 to 0.99) and a low MAPD of around 5% for most of the season. For G, seasonal values at both sites had lower values and narrow range compared to the remaining surface energy fluxes which makes it more difficult for the model to be able to capture the seasonal dynamics. However, the proposed method to estimate G was able to capture the seasonal dynamics at both sites from June to August, corresponding to the growing season peak, yielding lower MAPD results.

Model Performance Using Remote Sensing T RAD Estimates
As a first step towards regional implementation of the modified TSEB model over the boreal region, T RAD from satellite data was used as input. In addition boreal-specific modifications to α PTC and G were implemented using c GT (Equation (18)) instead of c G in Equation (14), thus allowing for a phase shift. Modifications in α PTC included an α PTC value for the whole season of 0.6 for black spruce, and an α PTC value of 0.5 in May and September, and 0.9 from June to August for birch, were applied to evaluate model performance. Overall results ( Table 5 and Figure 8) showed commensurate agreement for R N and G than those using T RAD from the pyrgeometer at each site, yielding LE and H slightly higher RMSE values but similar error percentage. These results suggest that regional implementation of TSEB for the boreal forest with the model modifications applied in this study are adequate to retrieve surface energy fluxes. R N yielded similar results compared to model application with pyrgeometer T RAD with an RMSE and an MAPD of around 20 W·m −2 and 4% for both black spruce and birch sites, respectively. G also showed similar results with an RMSE and an MAPD of around 7 W·m −2 and 40% for both settings.
For both boreal sites, LE and H modeled with remotely sensed T RAD yielded slightly higher RMSE of around 65 and 55 W·m −2 , respectively, although the MAPD was consistent to those results found with pyrgeometer data of around 30 and 20%, respectively. Thus, modifications in α PTC and f G together with remote sensing estimates of T RAD suggests that TSEB can be applied regionally to estimate turbulent fluxes. However, an error and bias increase compared to local estimates of T RAD using flux tower data might be explained in part by biases in T RAD -T A , due to atmospheric and emissivity corrections, sensor biases or unrepresentative T RAD at 1 km resolution relative to the flux tower footprint and meteorological forcing data. Thus, when applied regionally, TSEB is typically performed in time-differential mode because time-differences in T RAD reduce errors caused by uncertainty in atmospheric correction and emissivity for determining the absolute T RAD value.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 24 Model performance in estimating RN was comparable at both sites with a high correlation (R 2 from 0.95 to 0.99) and a low MAPD of around 5% for most of the season. For G, seasonal values at both sites had lower values and narrow range compared to the remaining surface energy fluxes which makes it more difficult for the model to be able to capture the seasonal dynamics. However, the proposed method to estimate G was able to capture the seasonal dynamics at both sites from June to August, corresponding to the growing season peak, yielding lower MAPD results.

Model Performance Using Remote Sensing TRAD Estimates
As a first step towards regional implementation of the modified TSEB model over the boreal region, TRAD from satellite data was used as input. In addition boreal-specific modifications to αPTC and G were implemented using cGT (Equation (18)) instead of cG in Equation (14), thus allowing for a phase shift. Modifications in αPTC included an αPTC value for the whole season of 0.6 for black spruce, and an αPTC value of 0.5 in May and September, and 0.9 from June to August for birch, were applied to evaluate model performance. Overall results ( Table 5 and Figure 8) showed commensurate agreement for RN and G than those using TRAD from the pyrgeometer at each site, yielding LE and H slightly higher RMSE values but similar error percentage. These results suggest that regional implementation of TSEB for the boreal forest with the model modifications applied in this study are adequate to retrieve surface energy fluxes. RN yielded similar results compared to model application with pyrgeometer TRAD with an RMSE and an MAPD of around 20 W·m −2 and 4% for both black spruce and birch sites, respectively. G also showed similar results with an RMSE and an MAPD of around 7 W·m −2 and 40% for both settings.
For both boreal sites, LE and H modeled with remotely sensed TRAD yielded slightly higher RMSE of around 65 and 55 W·m −2 , respectively, although the MAPD was consistent to those results found with pyrgeometer data of around 30 and 20%, respectively. Thus, modifications in αPTC and fG together with remote sensing estimates of TRAD suggests that TSEB can be applied regionally to estimate turbulent fluxes. However, an error and bias increase compared to local estimates of TRAD using flux tower data might be explained in part by biases in TRAD-TA, due to atmospheric and emissivity corrections, sensor biases or unrepresentative TRAD at 1 km resolution relative to the flux tower footprint and meteorological forcing data. Thus, when applied regionally, TSEB is typically performed in time-differential mode because time-differences in TRAD reduce errors caused by uncertainty in atmospheric correction and emissivity for determining the absolute TRAD value. Comparison of MODIS TRAD estimates with flux tower TRAD observations yielded a RMSE and a MBE of 1.9 K and −1.3 K at the birch site, respectively, and a RMSE and a MBE of 2.2 K and −1.3 K, respectively, at the black spruce site (Figure 9), thus, underestimating ground estimates of TRAD. This resulted in an increased bias and error in LE and H estimation lowering the final agreement. Therefore, to avoid an increase in bias on TRAD, for circumpolar boreal regions an alternative time- Comparison of MODIS T RAD estimates with flux tower T RAD observations yielded a RMSE and a MBE of 1.9 K and −1.3 K at the birch site, respectively, and a RMSE and a MBE of 2.2 K and −1.3 K, respectively, at the black spruce site (Figure 9), thus, underestimating ground estimates of T RAD . This resulted in an increased bias and error in LE and H estimation lowering the final agreement. Therefore, to avoid an increase in bias on T RAD , for circumpolar boreal regions an alternative time-differential TSEB technique using a dual-time difference (DTD model) in observed T RAD and air temperature [41,83] will be implemented adding the model modifications found in this paper. differential TSEB technique using a dual-time difference (DTD model) in observed TRAD and air temperature [41,83] will be implemented adding the model modifications found in this paper.  Moreover, it is worth noting that for the whole study period, there was a limited amount of TRAD data from MODIS. Although this is in part due to the constraints of data selection applied for the TSEB evaluation, clouds are an important source of contamination that hamper satellite applications using optical and thermal wavebands in the Arctic and boreal ecosystems. For instance, only 5% of the selected dataset was considered to be completely cloud-free. Thus, for regional applications of surface energy fluxes retrieval using remote sensing thermal infrared data as TRAD, a multiplatform approach from several polar satellites may be needed to increase the frequency of useable data, particularly during the growing season. Guzinski et al. [41] demonstrated good performance of a MODIS-driven DTD implementation over monitoring sites in Denmark, which have frequent cloud cover and are not well-sampled by polar orbits using combined Terra and Aqua satellite overpasses. Fortunately, near the poles, moderate resolution polar-orbiting thermal imaging systems, including Terra/Aqua MODIS, NOAA AVHRR, NPOES VIIRS, Sentinel 3 or microwave systems such as Aqua AMSR-E Ku-band [84] provide multiple image acquisitions per day, that allow mimicking a coarse temporal resolution geostationary system but at higher spatial resolutions increasing the likelihood of thermal infrared cloud-free pixels.
For regional implementation, meteorological inputs (air temperature, atmospheric pressure, wind speed, vapor pressure, and solar radiation) can be obtained from the NCEP North American Regional Reanalysis dataset, available at a spatial resolution of 32 km and a temporal resolution of 3 h. Land cover maps required to assign initial αPT values can be obtained from current land cover Moreover, it is worth noting that for the whole study period, there was a limited amount of T RAD data from MODIS. Although this is in part due to the constraints of data selection applied for the TSEB evaluation, clouds are an important source of contamination that hamper satellite applications using optical and thermal wavebands in the Arctic and boreal ecosystems. For instance, only 5% of the selected dataset was considered to be completely cloud-free. Thus, for regional applications of surface energy fluxes retrieval using remote sensing thermal infrared data as T RAD , a multiplatform approach from several polar satellites may be needed to increase the frequency of useable data, particularly during the growing season. Guzinski et al. [41] demonstrated good performance of a MODIS-driven DTD implementation over monitoring sites in Denmark, which have frequent cloud cover and are not well-sampled by polar orbits using combined Terra and Aqua satellite overpasses. Fortunately, near the poles, moderate resolution polar-orbiting thermal imaging systems, including Terra/Aqua MODIS, NOAA AVHRR, NPOES VIIRS, Sentinel 3 or microwave systems such as Aqua AMSR-E Ku-band [84] provide multiple image acquisitions per day, that allow mimicking a coarse temporal resolution geostationary system but at higher spatial resolutions increasing the likelihood of thermal infrared cloud-free pixels.
For regional implementation, meteorological inputs (air temperature, atmospheric pressure, wind speed, vapor pressure, and solar radiation) can be obtained from the NCEP North American Regional Reanalysis dataset, available at a spatial resolution of 32 km and a temporal resolution of 3 h. Land cover maps required to assign initial α PT values can be obtained from current land cover initiatives [80,81], or using the Terra/Aqua MODIS land cover product or the USGS National Land Cover Database (https://www.usgs.gov/centers/eros/science/national-land-cover-database). Although black spruce and birch are the two main covers of the boreal forest of interior Alaska, other minor covers representative of wetland ecosystems such as fen are also present in the boreal forest and can be parameterized by [27]. Finally, vegetation properties such as f G or LAI can be derived by temporally smoothing Terra/Aqua NDVI, EVI and LAI product data.

Conclusions
Modifications of a two-source energy balance (TSEB) model were applied to estimate daytime surface energy fluxes on two main vegetation types of the boreal forest of interior Alaska, black spruce and birch. An extensive model evaluation from leaf-out to senescence from 2012 to 2016 using local thermal data from the flux towers and Terra/Aqua MODIS remote sensing estimates as inputs was conducted. Model modifications included: (a) R N estimation for all-sky conditions (overcast and cloud-free), (b) a refined model for soil heat flux (G) previously applied for Arctic tundra and based on the T RAD -G relationship, and, (c) a α PTC parameterization for estimating canopy transpiration adopting an initial value of 1.26 and adjusted values of 0.6 and 0.9 for black spruce and birch, respectively.
Results showed that TSEB can appropriately model surface energy fluxes in the boreal forest from leaf-out to senescence. The modified model for soil heat flux estimation (G) yielded lower errors half that from the standard TSEB formulation. The R N model configuration for all sky conditions yielded similar errors to previous studies only for clear-sky conditions. TSEB modifications for boreal settings with adjusted α PTC provided turbulent heat flux estimates of H and LE with a mean RMSE value less than 50 W·m −2 and error percentage of around 23% in comparison with flux tower data that is comparable with errors typically described by other studies modelling surface energy fluxes. Results with the original α PTC configuration showed higher RMSE suggesting that vegetation and climate-type adjusted values of α PTC indicating different transpiration response to atmospheric forcing for black spruce, and birch would yield more accurate H and LE values.
A sensitivity analysis on α PTC and f G showed that for black spruce, an initial α PTC value of 0.6 together with a variable f G can be applied to estimate LE regionally for the whole growing season. For birch, although an initial α PTC value of 0.9 was successfully applied from June to September, seasonally adjusted values were also needed in the leaf-out and senescence periods with a α PTC value of 0.5 due to the influence of understory phenology.
When remote sensing estimates of T RAD were used as input, TSEB yielded slightly higher RMSE of around 60 W·m −2 for turbulent heat fluxes (H and LE), which was mainly attributed to an increase in bias compared to tower-based T RAD measurements. However, similar errors for G and R N were found compared with tower-based T RAD measurements.
For vegetation properties, Terra/Aqua MODIS vegetation indices (NDVI and EVI) and LAI time series fitted with TIMESAT provided reasonable estimates of LAI and f G , resulting in the model reproducing with good fidelity the temporal trends of energy partitioning together with boreal forest seasonality from leaf-out to senesce. This fact is particularly important for regional applications as LAI and f G field observations are not always available.
For regional applications, future work will be focused on incorporating the TSEB model improvements in boreal forest settings described in this study within the Atmospheric Land EXchange Inverse (ALEXI) surface energy balance regional modelling framework by using long time series of remote sensing estimates. It is also worth noting that for the whole study period, limited cloud-free MODIS thermal data were available. As clear-sky sky data are needed for regional applications, a multiplatform approach from several satellites polar satellites such as Terra/Aqua MODIS, NOAA AVHRR, NPOES VIIRS. Alternatively, microwave-derived T RAD using Aqua AMSR-E will be also needed to increase the likelihood of cloud-free data during non clear-sky conditions.

Appendix A. Summary of Equations Used to Estimate Aerodynamic Resistances in TSEB
In this appendix, the resistances involved in Equation (9), Equations (10) and (11) are summarized. R A is the aerodynamic resistance in the surface layer computed from the stability-corrected temperature profile equations (s·m −1 ) defined as: where d O is the displacement height, U is the wind speed measured at height z U , k is von Karman's constant (≈0.4), z T is the height of the T A measurement, Ψ M and Ψ H are the Monin-Obukhov stability functions for momentum and heat, respectively, and z OM is the aerodynamic roughness length. Both d O and z OM were estimated using LAI and canopy height (h C ) as follows [85,86]: where n is the within-canopy wind speed profile extinction coefficient parameterized as: where Cd is the drag coefficient of the foliage elements with a value of 0.2, and the ratio u * /u(h) is parameterized as: u * u(h) = 0.360 − 0.264 e (−15.1 * C d * LAI) (A4) Finally, z OM is defined as: R S is the resistance to the heat flow in the boundary layer above the soil (s·m −1 ) defined as: where a ≈ 0.004 m·s −1 , a ≈ 0.012 m·s −1 , U S is the wind speed (m·s −1 ) at the height above the soil surface where the effect of the soil surface roughness is minimal (~5 cm) and estimated from the [87] wind attenuation model through the canopy layer (see Equation (A8)). Finally, R X is the total boundary layer resistance of the vegetation leaf layer (s·m −1 ) defined as: where C is derived from weighting a coefficient in the equation for leaf boundary layer resistance (assumed a constant~90 s 1/2 m −1 ) and S is leaf size (m). U d+Zm is given by: where U C is the wind speed at canopy height (h C ) and a is the wind attenuation coefficient of the [87] wind attenuation model which is a function of LAI, h C and leaf size. For a full description of the TSEB resistance formulations see [28,31,32].