Estimation of surface energy fluxes in the Arctic tundra using the remote sensing thermal-based Two-Source Energy Balance model

The Arctic has become generally a warmer place over the past decades leading to earlier snow melt, permafrost degradation and changing plant communities. Increases in precipitation and local evaporation in the Arctic, known as the acceleration components of the hydrologic cycle, coupled with land cover changes, have resulted in significant changes in the regional surface energy budget. Quantifying spatiotemporal trends in surface energy flux partitioning is a key to 15 forecasting ecological responses to changing climate conditions in the Arctic. An extensive local evaluation of the twosource energy balance model (TSEB) a remote sensing-based model using thermal infrared retrievals of land-surface temperature was performed using tower measurements collected over different tundra types in Alaska in all sky conditions over the full growing season from 2008 to 2012. Based on comparisons with flux tower observations, refinements in the original TSEB net radiation, soil heat flux and canopy transpiration parameterizations were identified for Arctic tundra. In 20 particular, a revised method for estimating soil heat flux based on relationships with soil temperature was developed, resulting in significantly improved performance. These refinements result in mean turbulent flux errors generally less than 50 W·m at half-hourly timesteps, similar to errors typically reported in surface energy balance modelling studies conducted in more temperate climatic regimes. The MODIS leaf area index (LAI) remote sensing product proved to be useful for estimating energy fluxes in Arctic tundra in the absence of field data on local biomass amount. Model refinements found in 25 this work at the local scale build toward a regional implementation of the TSEB model over Arctic tundra ecosystems, using thermal satellite remote sensing to assess response of surface fluxes to changing vegetation and climate conditions.


Introduction
Air temperatures in the Alaskan Arctic have shown a significant increase, especially in past decade (Serreze and Barry, 2011).Results from models forced with a range of climate scenarios from the Intergovernmental Panel on Climate Change (IPCC) indicate that by the mid-21st century the permafrost area in the Northern Hemisphere is likely to decrease by 37-81 % (IPCC, 2014).In general, the Arctic has become a warmer place, leading to an acceleration of the hydrologic cycle, earlier snow melt and drier soils due to permafrost degradation (AMAP, 2012;Elmendorf et al., 2012;Rawlins et al., 2010;Sturm et al., 2001;Overduin and Kane, 2006).Furthermore, the hydrologic response of the Arctic land surface to changing climate is dynamically coupled to the region's surface energy balance (Vörösmarty et al., 2001), and the partitioning of energy fluxes plays an important role in modulating the hydrologic cycle of Arctic basins (Rawlins et al., 2010).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Evapotranspiration (ET, in units of mass, kg s −1 m −2 or mm d −1 ) or equivalently, latent heat flux (LE, in energy units, W m −2 ), is an important component of both the land surface hydrologic cycle and surface energy balance.As an example, Kane et al. (2004) reported water loss due to ET in the Imnavait Creek basin in Alaska is about 74 % of summer precipitation or 50 % of annual precipitation, as estimated from water balance computations.Even though ET is a significant component of the hydrologic cycle in Arctic regions, it is poorly quantified in Arctic basins, and the bulk estimates do not accurately account for spatial and temporal variability due to vegetation type and topography (Kane et al., 2000).In the Arctic, values of ET or LE are usually either derived from field estimates (Kane et al., 1990;Mendez et al., 1998) or calculated purely from empirical or quasi-physical models such as those described by Zhang et al. (2000) and Shutov et al. (2006) using meteorological station forcing data.However, due to remoteness, harsh winter conditions and the high costs of maintaining ground-based measurement networks, the data currently collected are also both temporally and spatially sparse.
In Arctic tundra ecosystems, several factors have contributed to the vegetation change, such as increased extent of severe fires, increased extent in deciduous vegetation or shrub encroachment in tundra ecosystems (Myers-Smith et al., 2011;Sturm et al., 2001), among others.Over at least the past three decades, Arctic ecosystems have shown evidence of "greening" (Xu et al., 2013;Bhatt et al., 2010), with about a 17 % increase in peak vegetation greenness for the Arctic tundra biome (Jia et al., 2003).Moreover, the forest-tundra transition zone is observed to be moving further north, tree heights are increasing and shrubs are becoming denser and taller (ACIA, 2004;AMAP, 2012).These changes in vegetation will have an important impact on the surface energy balance, especially in areas where shrubs have made their appearance in former tundra vegetation.This increase in leaf area index (LAI), together with canopy height and changes in the distribution of canopy elements, will augment the multiple scattering and absorption of radiation, likely resulting in a lower albedo (Beringer et al., 2005), although more detailed observations and measurements, particularly for the beginning of the snow-free period and peak growing season, are needed (Williamson et al., 2016).Also, according to Beringer et al. (2005), Bowen ratio increases from tundra to forested sites will result in an increasing dominance of sensible heat (H ) as the primary energy source heating the atmosphere.In the case of a transition from tundra to tall shrub and then to forest, H would likely increase during the growing season from ∼ 15 to nearly 30 %, respectively.This will have an important impact in the tundra energy partitioning, resulting in a positive feedback to the atmosphere that further warms the Arctic climate.However, the magnitude of changes in surface energy partitioning due to vegetation changes and resulting impact on local Arctic climate is still unclear and more research is needed to better under-stand these vegetation-change-atmosphere dynamics (Eugster et al., 2000;Jung et al., 2010).
In the last two decades, surface energy balance methods have demonstrated their utility in modeling water availability using diagnostic retrievals of energy fluxes from in situ or remote sensing data, especially data acquired in the thermal infrared (TIR) region (Kalma et al., 2008).While remote sensing estimates of ET over the Arctic exist from global modeling systems (Mu et al., 2009;Zhang et al., 2010), these modeling systems typically do not compute the full energy balance.To estimate energy fluxes at local scales, on the order of hundreds of meters, initiatives such as FLUXNET (http://fluxnet.ornl.gov/)provide eddy covariance flux measurements at discrete sites situated in different ecosystems across the US and globally.Unfortunately, there are few measurements sites in the Arctic (Mu et al., 2009), making the existing instrument network insufficient to capture pertinent details of the changing Arctic climate and landscape (ACIA, 2004;AMAP, 2012;Serreze and Barry, 2011;Vörösmarty et al., 2001).Detailed process-based (prognostic) land surface models can be also used to estimate coupled water and energy fluxes over landscapes (Duursma and Medlyn, 2012;Ek et al., 2003;Falge et al., 2005;Haverd et al., 2013;Smith et al., 2001;Vinukollu et al., 2012, among others); however, they may neglect important processes that are not known a priori.For example, Hain et al. (2015) demonstrated the value of comparing prognostic and TIR-based diagnostic latent heat flux estimates over the continental US to diagnose moisture sources and sinks that were not well-represented in the prognostic modeling system.
Given the critical need to better understand the water and energy balance over tundra ecosystems, and the role of changing climate and vegetation cover in driving these budgets, the aim of this work is to evaluate and refine a diagnostic TIR remote-sensing-based model for estimating seasonal dynamics of surface energy fluxes (LE, H , net radiation (R n ) and soil heat flux (G)) as well as energy partitioning in the Arctic tundra growing season from 2008 to 2012.Specifically, a refined version of the Two-Source Energy Balance model (TSEB; Norman et al., 1995) for Arctic tundra is evaluated with in situ forcing data from three eddy covariance flux towers in all sky conditions, and remote sensing estimates of vegetation properties (LAI, NDVI (normalized difference vegetation index) and EVI (enhanced vegetation index)) from the Moderate Resolution Imaging Spectroradiometer (MODIS).Model refinements include a new G parameterization and new configurations to retrieve R n (effective atmospheric emissivity), H and LE (two different Priestley-Taylor configurations).The TSEB serves as the land surface scheme in a regional Atmosphere-Land Exchange Inverse (ALEXI) modeling system (Anderson et al., 2011), implemented operationally over North America as part of NOAA's GOES Evapotranspiration and Drought Information System.Although the TSEB has been demonstrated to work well over a range of vegetation and climate conditions at midlatitudes (Anderson et al., 2007(Anderson et al., , 2011;;Choi et al., 2009;Sánchez et al., 2009;Tang, et al., 2011;Timmermans et al., 2007, among others), it has not yet been examined for tundra ecosystems characteristic of high latitudes.
2 Two-Source Energy Balance model: an overview Evapotranspiration (ET) can be estimated by surface energy balance models that partition the energy available at the land surface (R n − G, where R n is net radiation and G is the soil heat flux, both in W m −2 ) into turbulent fluxes of sensible and latent heating (H and LE, respectively, in W m −2 ): where L is the latent heat of vaporization (J kg −1 ) and E is ET (kg s −1 m −2 or mm s −1 ).
The model used in this study is the series version of the TSEB scheme originally proposed by Norman et al. (1995), which has been revised to improve shortwave and longwave radiation exchange within the soil-canopy system and the soil-canopy energy exchange (Kustas andNorman, 1999, 2000).A list of the TSEB inputs can be found in Table 1.TSEB has been successfully applied over rain-fed and irrigated crops and grasslands in temperate and semi-arid climates (Anderson et al., 2004(Anderson et al., , 2012;;Cammalleri et al., 2010Cammalleri et al., , 2012) ) but has not been previously applied over the Arctic tundra.
In the TSEB, directional surface radiometric temperature derived from satellite or a ground-based radiometer, T RAD (θ ) (K), is considered to be a composite of the soil and canopy temperatures, expressed as where T c is canopy temperature (K), T s is soil temperature (K) and f c (θ ) is the fractional vegetation cover observed at the radiometer view angle θ .For a canopy with a spherical leaf angle distribution and LAI, f c (θ ) can be estimated as where the factor indicates the degree to which vegetation is clumped as in row crops or sparsely vegetated shrubland canopies (Kustas andNorman, 1999, 2000).The composite soil and canopy temperatures are used to compute the surface energy balance for the canopy and soil components of the combined land surface system: where R N s is net radiation at the soil surface, R N c is net radiation divergence in the vegetated canopy layer, H c and H s are canopy and soil sensible heat flux, respectively, LE c is the canopy transpiration rate, LE s is soil evaporation and G is the soil heat flux.The net shortwave radiation is calculated from the measured incoming solar radiation and the surface albedo, while net longwave radiation is estimated from the observed air and land surface temperatures, using the Stefan-Boltzmann equation with atmospheric emissivity from the Brutsaert (1975) method.By permitting the soil and vegetated canopy fluxes to interact with each other, Norman et al. (1995) derived expressions for H s and H c expressed as a function of temperature differences where and with the total sensible heat flux H = H c + H s expressed as where ρ is air density (kg m −3 ), C p is the specific heat of air (kJ kg −1 K −1 ), T A c is air temperature in the canopy air layer (K), T a is the air temperature in the surface layer measured at some height above the canopy (K), r X is the total boundary layer resistance of the complete canopy of leaves (s m −1 ), r s is the resistance to sensible heat exchange from the soil surface (s m −1 ) and r a is aerodynamic resistance (s m −1 ) defined by In Eq. ( 9), 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.
The original resistance formulations are described in more detail in Norman et al. (1995) with revisions described in Kustas andNorman (1999, 2000).Weighting of the heat flux contributions from the canopy and soil components is performed indirectly by the partitioning of the R n between soil and canopy and via the impact on resistance values from the fractional amount and type of canopy cover (see Kustas and Norman, 1999).
For the latent heat flux from the canopy, the Priestley-Taylor formula is used to initially estimate a potential rate for LE c : where α PTC is a variable quantity related to the Priestley-Taylor coefficient (Priestley and Taylor, 1972), but in this  case defined exclusively for the canopy component, which was suggested for row crops by Tanner and Jury (1976) and normally set to an initial value of 1.2, is the slope of the saturation vapor pressure versus temperature curve and γ is the psychrometric constant (∼ 0.066 kPa • C −1 ).f g is the fraction of green vegetation that according to Guzinski et al. (2013) and Fisher et al. (2008) can be estimated through the NDVI and the EVI: Under stress conditions, TSEB iteratively reduces α PTC from its initial value.The TSEB model requires both a solution to the radiative temperature partitioning (Eq.2) and the energy balance (Eqs.6 and 7), with physically plausible model solutions for soil and vegetation temperatures and fluxes.Nonphysical solutions, such as daytime condensation at the soil surface (i.e., LE s < 0), can be obtained under conditions of moisture deficiency.This happens because LE c is overestimated in these cases by the Priestley-Taylor parameterization, which describes potential transpiration.The higher LE c leads to a cooler T c and T s must be accordingly larger to satisfy Eq. ( 7).This drives H s higher, and the residual LE s from Eq. ( 12) can become negative.If this condition is encountered by the TSEB scheme, α PTC is iteratively reduced until LE s ∼ 0 (expected for a dry soil/substrate surface).However, there are instances where the vegetation is not transpiring at the potential rate but is not stressed due to its adaption to water and climate conditions (Agam et al., 2010) or the fact that not all the vegetation is green or actively transpiring (Guzinski et al., 2013) (a thorough discussion of conditions that force a reduction in α PTC , can be also found in Anderson et al., 2005 andLi et al., 2005).
The latent heat flux from the soil surface is solved as a residual in the energy balance equation: with G estimated as a fraction of the net radiation at the soil surface (c G ): During the midmorning to midday period, the value of c G can be typically assumed to be constant (Kustas and Daughtry, 1990;Santanello and Friedl, 2003).In this case, a typical value of ∼ 0.3 can be assumed for c G based on experimental data from several sources (Daughtry et al., 1990).However, c G value varies with soil type and moisture conditions as well as time, due to the phase shift between G and R N s over a diurnal cycle (Santanello and Friedl, 2003).
3 TSEB formulation refinements for Arctic tundra

Downwelling longwave radiation estimation: effective atmospheric emissivity for all sky conditions
The original TSEB formulation estimates the downwelling longwave radiation component of R n using the effective atmospheric emissivity (ε) method described in Brutsaert (1975) for clear sky conditions: where e is the water pressure in millibars and T a in K, and C is 1.24 as in the original Brutsaert (1975) formulation.
However, in this study, TSEB is applied for all sky conditions, including clear sky, partially cloudy and overcast conditions.To estimate ε for all sky conditions, Crawford and Duchon (1999) proposed a methodology that incorporated the Brutsaert (1975) clear sky parameterization and the Deardorff (1978) cloudiness correction using a simple cloud modification introducing a cloud fraction term (clf) according to the following equation: The clf is defined as where s is the ratio of the measured solar irradiance to the clear sky irradiance.Shortwave clear sky irradiance used in Eq. ( 16) may be obtained through the methodology proposed by Pons and Ninyerola (2008), where incident clear sky irradiance is calculated through a digital elevation model at a specific point during a particular day of the year, taking into account the position of the Sun, the angles of incidence, the projected shadows, the atmospheric extinction and the distance from the Earth to the Sun.For Arctic areas, Jin et al. ( 2006) suggested an improved formulation of C for clear sky conditions that can also be applied in Eq. ( 15) for all sky conditions, defined as In order to evaluate if the Jin et al. ( 2006) method offered more accurate estimates of ε for Arctic conditions, this method was compared to Brutsaert (1975) formulation used in TSEB in both cases for all sky conditions using Eq. ( 15).

Refinements in soil heat flux parameterization: c G coefficient and definition of a new coefficient based on T RAD
In the Arctic tundra, the propagation of the thawing front in the soil active layer consumes a large proportion (around 18 %) of the energy input from the positive net radiation (Boike et al., 2008a;Rouse, 1984).Moreover, the presence of permafrost in tundra areas may contribute to the large tundra soil heat flux by creating a strong thermal gradient between the ground surface and depth, offsetting the influence of the highly insulative moss cover which would otherwise have been expected to reduce soil heat flux (Myers- Smith et al., 2011;Sturm et al., 2001).Therefore, previous formulations of soil heat flux used in TSEB applications, mainly representative of cropped and sparse-vegetated areas in the US, need to be adjusted and validated for Arctic tundra.
Currently there are several methodologies that allow estimating soil heat flux from tenths of centimeters to meters in depth in the Arctic tundra by using modeling or instrumentation at several depths (Lynch et al., 1999;Ekici et al., 2015;Jiang et al., 2015;Romanovsky et al., 1997;Yao et al., 2011;Zhuang et al., 2001;Hinzman et al., 1998).However, in this study, a simple approach based on the relationship between G and R N s (Eq.13) was used to estimate the soil heat flux in the near-surface soil layer (around 10 cm depth).This approach has less complexity and requires less input data than the methods mentioned above and allows estimating G at regional scales.
In early TSEB implementation, a constant value of c G value around 0.3 was used to estimate G for the midmorning to midday period (Eq.13 based on findings by Kustas and Daughtry (1990) for US study sites.However, this assumption can result in significant errors if applied out of this time range.For diurnal hourly timescales, Kustas et al. (1998) developed a method to estimate c G based on time differences with the local solar noon quantified by a non-dimensional time parameter.Although this approach does not consider the phase shift between G and R N s over a diurnal cycle, a phase shift was included in the model proposed by Santanello and Friedl (2003) in the following form: where A represents the maximum value of c G , B is chosen to minimize the deviation of c G from Eq. ( 13), t is time in seconds relative to solar noon and S is the phase shift between G and R N s in seconds.Values fitted for A, S and B were 0.31, 10 800 and 74 000, respectively.Although c G values for Arctic tundra were not found in the literature, several studies (Beringer et al., 2005;Eugster et al., 2000Eugster et al., , 2005;;Boike et al., 2008b;Eaton et al., 2001;Kodama et al., 2007;Langer et al., 2011;Soegaard et al., 2001;Westermann et al., 2009;Mendez et al., 1998;Lund et al., 2014) present the relationship between R N s and G during the summer months in similar tundra areas.According to these studies, a mean value of 0.14, as a maximum value of c G in Eq. ( 18), can be derived from different analyses of R N s and G over the Arctic tundra.
An alternative parameterization for G was suggested by Santanello and Friedl (2003) for several types of soils with crops and by Jacobsen and Hansen (1999) for Arctic tundra that links the soil heat flux to the diurnal variations in surface radiometric temperature.This approach can also be applied for Arctic tundra as follows: where c TG is a coefficient that represents the relationship between the diurnal variation of T RAD and G.For diurnal hourly timescales, c TG can be also estimated using the phase shift proposed in Eq. ( 18), where, in this case, S is the phase shift between G and T RAD in seconds.This new approach avoids using R N s , which is more difficult to define in tundra systems given the influence of the surface moss layer above the mineral soil.Moreover, A, S and B in Eq. ( 18) can be fitted by using direct measurements of T RAD from thermal field sensors, commonly available on flux towers (pyrgeometer), or thermal data from geostationary or polar satellites.Thus, to evaluate soil heat flux for diurnal hourly timescales, the approaches of Kustas et al. (1998) and Santanello and Friedl (2003) were compared using the original c G value of 0.30 and a new value for Arctic tundra of 0.14, both as maximum values of c G in Eq. ( 18).A, B and S values for the new c TG approach were fitted and tested using an extended evaluation dataset and then compared to these radiation-based methods (see Sect. 4.2).

Priestley-Taylor coefficient
In the original TSEB formulation, the Priestley-Taylor approach for the canopy component of LE is used.In this case, α PTC is normally set to an initial value of 1.26 for the general conditions tested during the growing season in rangelands and croplands.For stressed canopies, TSEB internally modifies α PTC to yield reasonable partitioning between LE c and LE s .
As with the c G coefficient, specific α PTC values for tundra were not found in the literature.Alternatively, measurements of bulk (soil plus canopy) for Arctic tundra systems are available (Beringer et al., 2005;Eaton et al., 2001;Eugster et al., 2005;Engstrom et al., 2002;Mendez et al., 1998;Lund et al., 2014), suggesting a mean value of around 0.92.This bulk value might suggest that α PTC could also be lower for Alaska tundra summer conditions.For natural vegetation, Agam et al. (2010) also suggested that a lower α PTC value might yield better results.Therefore, for modeling purposes, two different values of α PTC values, 0.92 and 1.26, were applied to evaluate which nominal α PTC input to TSEB was more appropriate for Arctic tundra.
4 Study area and data description

Study area
To refine and evaluate the TSEB model for Alaska's Arctic tundra summer conditions, three eddy covariance flux towers (referred to as Fen, Tussock and Heath; see Fig. 1) were selected.These are located across the Imnavait Watershed (∼ 904 m a.s.l.) with eddy covariance and associated meteorological data collection beginning in 2007 (Euskirchen et al., 2012;Kade et al., 2012).A brief description of instrumentation at the tower sites is provided in  4.2 Model inputs, evaluation datasets and metrics

Micrometeorological input data
Data incorporated in this study spanned the period from May to September in 2008 to 2012.These included eddy covariance data for latent and sensible collected at 10 Hz and processed to 30 min means (described below) as well as meteorological data collected at 30 min intervals (Tables 1 and 2).These data, from under all sky conditions, were used to refine and evaluate the model performance (Table 1).This dataset was considered to be representative of the short Arctic tundra vegetative cycle from early growing to senescence as well as used to capture inter-and intra-annual vegetation dynamics.Meteorological inputs for TSEB include wind speed, air temperature, vapor pressure, atmospheric pressure, longwave incoming radiation and solar radiation, all of which were collected at the three measurement sites (see Tables 1 and 2).The surface radiometric temperature T RAD inputs were obtained from the pyrgeometer sensor at the Tussock station and from infrared radiometer sensors at both Fen and Heath stations.

Remote sensing input data: vegetation properties
In addition, TSEB also requires estimates of LAI and the fraction of vegetation that is green to specify f c in Eq. ( 2) and to estimate LE c in Eq. ( 10).While in situ measurements of LAI were not available at the tower sites for the length of this study, the 500 m combined Terra/Aqua MODIS 4day LAI product (MCD15A3H) was available for the study area.This product has been successfully applied in other applications of the TSEB (Guzinski et al., 2013), where sites are considered homogeneous over several kilometers, and serves here as a proxy for local observations.The fraction of vegetation that is green (f g ) in Eq. ( 10) was estimated using NDVI and EVI from MODIS imagery using the daily 250 m reflectance product (MOD09GQ), and using the blue band in the daily 500 m reflectance product (MOD09GA) to correct for residual atmospheric effects, with negligible spatial artifacts.Because MODIS time series contains occasional lower-quality data, gaps from persistent clouds, cloud contamination and other gaps (Gao et al., 2008), a program for analyzing time series of remote sensing imagery, TIMESAT (Jönsson and Eklundh, 2004), was used to produce temporally smoothed NDVI, EVI and LAI by selecting the best estimates through these products' quality flags.Gao et al. (2008) found a good agreement with field measurements when smoothing MODIS LAI data using this distribution and several weights (w) based on the product quality flags (w = 1.0 for LAI retrievals from the radiative-transfer model (high quality) or for LAI retrieval that reaches saturation, w = 0.25 for retrievals from an empirical model and w = 0.0 for all invalid and fill values).Beck et al. (2006) also reported that an asymmetric Gaussian distribution was appropriate for describing vegetation dynamics using NDVI at high latitudes and several weights (w) based on the product quality flags (highest quality/clear, mixed and cloudy were assigned weights of 1, 0.5 and 0, respectively).For NDVI, EVI and LAI time series smoothing, the weights and quality flags proposed by Beck et al. (2006) and Gao et al. (2008) were used.Vegetation height, used to define roughness parameters d o and z om , was assigned based on measurements made in the vicinity of the flux towers (Kade et al., 2012), and the clumping factor was set to 1 for all sites based on the knowledge that Arctic tundra has a variable moss layer with little bare ground.Variability regarding these inputs for the studied periods is shown in Table 1.Moreover, to ensure that only snow-free periods were analyzed, Terra/Aqua MODIS snow cover products (MOD10A1 and MYD10A1) were used to screen days with snow cover at the beginning and end of the growing season.

Micrometeorological flux data for model evaluation
The eddy covariance data used in TSEB evaluation, including latent and sensible heat, were processed with EddyPro ® (2014) software.Changes in mass flow caused by changes in air density (Webb et al., 1980), corrections for frequency attenuation of eddy covariance fluxes following Massman (2000) and Rannik (2001) and storage corrections for calm periods (friction velocity (u * ) was less than 0.1 m s −1 , as suggested by Rocha and Shaver, 2011) were accounted for.The automatic gain control (AGC) value (which represents optical impedance by precipitation) was computed for the infrared gas analyzer (IRGA) and used as a QA/QC variable for both flux and radiation data, with 60 as the maximum threshold value (LI- COR Inc, 2004).Rejection angles of 10 • were also used when the eddy covariance instruments were downwind of a tower to remove flow distortions.In addition, corrections for stationarity, lags and step change, among others, were performed by the flux processing software (for further information on micrometeorological data processing, see Euskirchen et al., 2012, and http: //aon.iab.uaf.edu/data_info).To select the best data available, the above criteria were used to flag the micrometeorological dataset, and no gap-filled data were used.
In addition, soil heat flux plate measurements were corrected to account for soil heat storage above the plate according to the calorimetric methodology proposed by Domingo et al. (2000) and Lund et al. (2014) using existing field measurements of soil bulk density for each site (758, 989 and 1038 kg m −3 for Fen, Tussock and Heath flux stations, respectively), soil moisture from the water content reflectometer and thermocouple averaging soil temperature probes (TCAVs) placed at two depths in the soil (see Table 2).
To evaluate the new c TG approach, a total of 41 068 halfhourly time steps of T RAD and G from 04:00 to 21:00 LST were selected (11 593, 14 454 and 15 021 for Fen, Tussock and Heath flux stations, respectively).Coefficients A, B and S were fitted using 60 % of all available data (fitting subset) aggregated in 30 min time steps for the whole summer period.The remaining 40 % of the data were reserved for model testing (test subset) (see Table 4 for flux stations distribution).To evaluate the TSEB model, including G retrieve from Kustas et al. (1998) and Santanello and Friedl (2003) approaches, a total of 5178 half-hourly time steps (1558, 1273 and 2347 for Fen, Tussock and Heath flux stations, respectively) were subset from the previous selection by imposing three criteria: (a) energy closure at the half-hourly timescale exceeded 70 %, (b) R n was higher than 100 W m −2 in order to ensure daylight conditions and (c) no precipitation was present.20) to (24), respectively.

Evaluation metrics
where e i refers to the estimated value of the variable in question (R n , H , LE or G), o i is the observed value (in situ measurement provided by the flux station), n is the number of data points, and O and E are the average of the o i and e i values, respectively.5 Results and discussion

Evaluation of soil heat flux model refinements for tundra
Both the Kustas et al. (1998;K98) and the Santanello and Friedl (2003;SF03) soil heat flux models used to estimate G at the study sites yielded high errors when a value of c G = 0.3 was used, with MAPD ranging from 90 to 159 %.In this case, the SF03 approach provided better results (Table 3).It is important to note that G is a relatively small term with a maximum value on the order of 50 W m −2 .Both models generally overestimated G with a MBE from 3 to 40 W m −2 , with the SF03 model generating lower biases.Results improved when a c G value of 0.14 was used, with MAPD ranging from 48 to 76 % and with lower RMSE values from 15 to 21 W m −2 and MBE from −4 to −14 W m −2 .With the lower value of c G , the K98 approach provided better results (Table 3).Similar to the original c G , c TG can be also estimated using the Santanello and Friedl (2003) method in Eq. ( 18).Mean diurnal profiles in T RAD and G, averaged over all tundra sites (see Sect. 4.2 and Table 3) showed a phase shift between these variables (Fig. 2).The mean G value for the summer period peaked around 15:00 LT, with a phase shift around 4 h after the maximum T RAD at noon.Using T RAD and G observations at half-hourly time steps from the fitting subset, diurnal c TG curves were derived for the growing season for each of the tower sites, showing reasonable agreement (Fig. 3).A fit to the mean curve yielded parameter values of S = −14 400 s, A = 1.55 and B = 160 000 s.As in the case of Santanello and Friedl (2003), a B variation of ±15 000 s had no significant influence on the results.Statistical comparisons between observed fluxes from the test subset and simulations using the fitted parameters show good agreement and negligible bias (Table 4) with R 2 , MAPD, RMSE and MBE values of 0.68, 37 %, 6 and 0 W m −2 , respectively.In addition, the new model was also evaluated using the same flux subset used in Table 3 to assess the K98 and SF03 configurations, demonstrating improved performance with roughly half the MAPD than K98 and SF03 configurations (Table 4).
The performance of the G parameterization for Arctic tundra reported here is comparable or superior to previous studies reported in the literature using the Santanello and Friedl (2003) or Kustas et al. (1998) approaches for other ecosystems.In shrub-grass-dominated areas and boreal forest, several studies (Anderson et al., 2008;Kustas et al., 1998;Li et al., 2008;Sánchez et al., 2009;Timmermans et al., 2007) reported MAPD and RMSE values ranging from 19 to 59 % and from 15 to 35 W m −2 , respectively.Studies in corn and soybean crops (Anderson et al., 2005;Choi et al., 2009;Li et al., 2005;Santanello and Friedl, 2003) reported MAPD and RMSE values ranging from 19 to 34 % and from 10 to 41 W m −2 , respectively.

Net radiation evaluation: effective atmospheric emissivity
Effective atmospheric emissivity estimated using the Brutsaert ( 1975 Estimated R n for all sky conditions yielded strong agreement with observed values for all flux towers (see Fig. 4 and Table 5) with a mean R 2 , MAPD, MAD and RMSE of 0.99, 7 %, 18 W m −2 and 23 W m −2 , with a tendency to overestimate R n with a MBE of 7 W m −2 .In terms of RMSE and MAPD, all study sites behaved similarly (see Fig. 4).These results are in line with previous TSEB model applications for other cover types and clear sky conditions where a MAPD of around 5 % was reported (Anderson et al., 2000(Anderson et al., , 2005(Anderson et al., , 2008;;Li et al., 2005Li et al., , 2008;;Kustas and Norman, 1999;Guzinski et al., 2013).This suggest that R n estimation using this scheme can be applied regionally under summer all sky conditions in Arctic tundra when a source of solar radiation (METEOSAT or GOES; Cristóbal and Anderson, 2013), air temperature (Cristóbal et al., 2008) and T RAD (MODIS land surface temperature and emissivity product) are available.

Latent and sensible heat fluxes evaluation: α PTC configuration for Arctic tundra
The average energy balance closure using half-hour periods for the evaluation dataset was 88 %, which is in agreement Hydrol.Earth Syst.Sci., 21, 1339-1358, 2017 www.hydrol-earth-syst-sci.net/21/1339/2017/ with the average closure of 90 % for these flux stations (Euskirchen et al., 2012).Lack of closure may be explained by instrument and methodological uncertainties, insufficient estimation of storage terms, unmeasured advective fluxes, landscape-scale heterogeneity or instrument spatial representativeness, among others (Lund et al., 2014;Stoy et al., 2013;Foken et al., 2011;Foken, 2008;Wilson et al., 2001).More recently, there is evidence that non-orthogonal sonics underestimate vertical velocity causing undermeasurement of H and LE on the order of 10 % (Kochendorfer et al., 2012a;Frank et al., 2013), although this is still being debated (Kochendorfer et al., 2012b).While, currently, there is no uniform answer on how to deal with non-closure of the energy balance in eddy covariance datasets, and methods for analyzing the reasons for the lack of closure are still under discussion (Foken et al., 2011), in this study, TSEB output is primarily compared with eddy covariance fluxes as observed, without closure corrections.However, to facilitate comparisons to numerous studies in the literature imposing energy conservation to eddy covariance data when evaluating surface energy balance models (Courault et al., 2005;Kalma et al., 2008), and given that strong evidence is presented in the literature that both H and LE are undermeasured by the eddy covariance technique, additional comparisons with closed fluxes using the Bowen ratio (H BR and LE BR ) approach suggested by Twine et al. (2000) and LE recalculated as the residual (LE RES , e.g., Li et al., 2008) are provided for completeness.Results with closed fluxes are presented to provide bounds on the range in probable model performance and to demonstrate the impact of closure corrections on model evaluation metrics.LE and H estimated through both the new proposed soil heat flux methodology and the all sky R n methodology scheme yielded reasonable agreement with observed halfhourly unclosed turbulent fluxes for both α PTC parameterizations of 0.92 and 1.26 (see Tables 4 and 5, and Fig. 4), although α PTC = 0.92 yielded marginally lower errors for H and LE.Relative errors (MAPDs) were 40 and 25 % for LE and H , respectively, for all combined sites using α PTC = 0.92, and 45 and 27 % using the standard value of α PTC = 1.26, respectively.A slight improvement in H and LE estimates using α PTC ∼ 0.9 also agrees with Agam et al. (2010) who also found better results with lower α PTC for natural vegetation in water-limited environments.
When energy balance closure is imposed, model performance is mainly improved for LE (up to 10 % decrease in MAPD) in part due to the fact that the measured turbulent fluxes are adjusted to achieve energy conservation as required by surface energy balance models.Indeed, the relatively small errors between modeled and measured R n (Table 5) and relatively small and unbiased magnitude in modeled/measured G (Table 5 and Fig. 4) suggest the unclosed turbulent fluxes contribute to the error statistics comparing Tables 5 and 6 with closed H and LE observed fluxes.Nevertheless, since the mean RMSE for all fluxes compared to unclosed and closed turbulent fluxes and for all parameterizations and sites was around 50 W m −2 (Tables 5 and  6), which is commensurate with errors typically reported in other surface energy balance studies (Kalma et al., 2008), these results suggest that a generalized α PTC value of 1.26 in global TSEB applications may adequately reproduce energy fluxes in Arctic tundra during the growing season, from leafout until senescence, while also capturing inter-and intraannual dynamics.However, biases in regional applications  may be reduced by using a land cover class-dependent value of α PTC .
Currently, there is limited research published on application of energy balance models to estimate energy fluxes for Arctic tundra.Mu et al. (2009) reported year-round errors from 20 to 40 % in two Arctic tundra sites in Barrow (Alaska, USA) at daily periods based on a modified aerodynamic resistance-surface energy balance model where the required surface conductance is estimated from remotely sensed LAI based on the Cleugh et al. (2007) formulation.TSEB results, however, were evaluated with half-hourly data in summer conditions and, although they cannot be directly compared with results in this previous study, they show similar errors.As in the case of R n , LE and H results are also in line with previous works for other cover types using in situ data as input to TSEB (Anderson et al., 2000(Anderson et al., , 2008;;Li et al., 2005).

Seasonal dynamics of surface energy fluxes and energy partitioning
In general, monthly estimation of surface energy fluxes showed a good agreement with observations during the growing season.Because the model yielded similar results with both α PTC parameterizations of 0.92 and 1.26, this section only shows the seasonal dynamics with α PTC of 0.92.Because of the undermeasurement issues with eddy covariance data and greater uncertainty and error associated with LE measurements (Wolf et al., 2008), seasonal dynamics of turbulent fluxes were compared with residual LE.Estimated R n yielded a low MAPD around 6 %, increasing up to 12 % at the end of the growing season (Table 7 and Fig. 5).The proposed new method to estimate G yielded better MAPD results from June to August, which coincides with the peak of the growing season in July.A similar pattern was found for LE and H , where the best MAPD results occurred also in the middle of the growing season (June and July).MAPD for LE, H and G tended to be higher in May and September, thus coinciding with earlier plant growth or the senesce periods, respectively.The MODIS LAI product, used to estimate the fractional vegetation cover (Eq. 3) to partition soil and canopy temperatures, performed as a good proxy to capture inter-and intra-annual vegetation dynamics (Fig. 6).Mean seasonal MODIS LAI from May to September for all flux stations was 1.2 ± 0.5 m 2 m −2 .Previous studies close to the study area, in Toolik Lake and Imnavait Creek (Shaver and Chapin, 1991;Shippert et al., 1995;Williams et al., 2001;Williams et al., 2006), reported LAI field estimates ranging from 0.2 to 1.4 m 2 m −2 for different tundra types  around mid-July to mid-August, suggesting LAI overestimation from the MODIS product.Loranty et al. (2010) also reported LAI overestimation when using this product in similar tundra types, finding better agreement using a NDVI-LAI relationship (Shaver et al., 2007;Street et al., 2007), although the nonlinearity in the NDVI-LAI conversion is prone to averaging errors when scaled with remote sensing data (Stoy et al., 2009).Despite MODIS LAI overestimation, it performed well for the Arctic tundra, suggesting utility for regional applications, although NDVI-LAI methods might be considered for future applications.f g estimated through NDVI and EVI also captured interand intra-annual vegetation dynamics (Fig. 6), with a mean seasonal value from May to September for all flux stations of 0.82 ± 0.7.From May to August (from the beginning and almost to the end of the growing period), f g showed a good agreement with LAI dynamics.However, while f g showed a steady increase at the beginning of the growing season, it did not follow MODIS LAI dynamics in September.This caused the model to overestimate LE and underestimate H during this time period, degrading agreement with observed data.The underperformance of the f g methods near the end of the growing season might be related to the presence of a variable moss layer, which can exert strong controls on understory water and heat fluxes in Arctic tundra ecosystems (Blok et al., 2011) and may have masked the actual vegetation dynamics (Fig. 6).Further research is needed to confirm this hypothesis.The pattern of daily estimated surface energy fluxes also compared well to observed fluxes for all sky conditions.As an example, time series of modeled and measured surface energy fluxes are segmented in Fig. 7 8).Observed and estimated Bowen ratio (β) yielded mean values of 0.60 and 0.67, respectively.In all cases, observed and estimated results are in line with previous studies for Arctic tundra (Lynch et al., 1999;Eugster et al., 2000).It is worth noting that the difference between observed and estimated values of LE / R n and H / R n partitions was only around 3 % and for G / R n was almost negligible.From June to August, mean absolute difference values between observed and estimated values for LE / R n and H / R n were around 4 %, increasing up to 15 % in September due to model over-and underestimation, while the G / R n difference was only less than 1 %.
These results suggest that the model is able to reproduce accurately temporal trends of energy partition in concert with tundra vegetation dynamics in the growing vegetation peak from June to August and could be used to monitor changes in surface energy fluxes concurrently with vegetation dynamics.

Conclusions and future work
Parameterizations for R n , G and α PTC used in the TSEB model were evaluated and refined for applications in different tundra types in Alaska over the full Arctic tundra growing season.Results showed that TSEB may adequately reproduce energy fluxes in Arctic tundra during the growing season, from leaf-out until senescence.The modified TSEB provided turbulent heat flux estimates with a mean RMSE value on the order of 50 W m −2 in comparison with unclosed eddy covariance measurements of H and LE collected at three flux towers -commensurate with errors typically reported in other surface energy balance studies.Moreover, as in many other studies using eddy covariance flux data for evaluating model performance, imposing energy balance closure yielded better agreement between measured and modeled turbulent fluxes.The all sky R n estimation scheme tested here yielded similar errors to those from other studies for only clear sky conditions.This demonstrates potential for regional-scale applications when reliable sources of solar radiation, air temperature and T RAD are available.A refined model for soil heat flux (G), based on the soil temperature-G relationship, was evaluated from the green-up phase to senescence using data from multiple years, and yielded errors half the magnitude of the standard TSEB formulation based on the relationship between R N s and G.The TSEB α PTC parameterization for estimating canopy transpiration (LE c ) was tested using the standard TSEB value of 1.26 and a value of 0.92 suggested in the literature for Arctic tundra, and both parameterizations yield similar flux errors suggesting tundraspecific values of α PTC are not needed.
In the absence of in situ measurements of LAI within the vicinity of the tower sites, the MODIS LAI product provided reasonable inputs for localized model testing.The model was able to reproduce accurately temporal trends of energy partitioning in concert with tundra vegetation dynamics in the peak growing season.Moreover, it also has potential to monitor changes in surface energy fluxes in Arctic tundra due to changes in vegetation composition (e.g., shrub encroachment).This is particularly crucial in the Arctic where there is a sparse network of meteorological and flux observations.Further research is needed regarding the specific role of the moss layer in modifying remote sensing estimates of green vegetation cover fraction and soil heat conduction within tundra ecosystems.
Future work will incorporate the TSEB model refinements identified here for Arctic tundra into regional and global applications of the ALEXI surface energy balance modeling system.Model performance within a fully satellite-based remote sensing framework will be compared to the local evaluations reported here at these tundra flux sites.In addition, the diagnostic assessments of ET and surface energy fluxes will be compared with regional output from process-based prognostic land surface models to better understand the strengths and weaknesses of both types of modeling systems.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Location of the Fen, Tussock and Heath flux towers at the Imnavait watershed.Right panel map is in UTM-6N NAD83 with coordinates in kilometers.
For model evaluation, surface energy fluxes (R n , LE, H and G) from the flux datasets (observed values) were compared to TSEB outputs (estimated values) using five metrics describing model errors and biases: the coefficient of determination (R 2 ) was used to indicate the precision of the estimates in relation to observed surface energy fluxes; the root mean square error (RMSE) was used as a measure of accuracy to measure differences between values estimated by the TSEB model and values actually observed by the flux towers; the mean bias error (MBE) was used to indicate cumulative offsets between measured and observed values; the mean absolute difference (MAD) was used to indicate the magnitude of the average absolute difference of observed and estimated values; and finally, the mean absolute percent difference (MAPD) was used to express the magnitude of absolute difference between observed and estimated values relative to the observed mean value, from Eqs. (

Figure 2 .
Figure 2. Mean daytime cycle for G and T RAD in the study area computed using all data available from the Fen, Tussock and Heath flux towers from 2008 to 2012.

Figure 3 .
Figure 3.Time series of modeled c TG and observed c TG values from the Fen, Tussock and Heath flux stations as well as mean values for summer conditions (bars represent standard deviation of the mean).

Figure 4 .
Figure 4. Comparison of modeled vs. measured half-hourly fluxes using α PTC of 0.92.The 1 : 1 line represents perfect agreement with observations.Columns represent results with unclosed observed turbulent fluxes (left), Bowen ratio closure (middle) and residual closure (right) for the Fen, Tussock, Heath sites and all sites combined (rows).
) and Jin et al. (2006) methodologies yielded similar errors in simulated downwelling longwave radiation results, with a R 2 of 0.58 and a RMSE of 26 and 27 W m −2 , respectively.The C coefficient computed through Jin et al. (2006) yielded a value of 1.25 ± 0.009, very close to Brutsaert (1975) C value of 1.24.This suggests that the simpler Brutsaert (1975) C coefficient can be used efficiently to model effective atmospheric emissivity in all sky conditions when combined with Crawford and Duchon (1999) and Pons and Ninyerola (2008) methods for summer Arctic tundra.

Figure 5 .
Figure 5.Comparison of modeled vs. observed half-hourly surface fluxes (using LE from residual closure) by month using α PTC of 0.92 and G estimated by the new model.The 1 : 1 line represents perfect agreement with observations.

Figure 6 .
Figure 6.Mean MODIS LAI and fraction of green vegetation (f g ) temporal dynamics for all flux stations from 2008 to 2012.

Figure 7 .
Figure 7.Comparison of hourly flux tower RN, LE, H and G observations (using LE from residual closure) ( o ) (from 06:00 to 21:00 LST) at the Heath flux tower with model estimates ( e ) using α PTC of 0.92.Each diurnal segment represents flux data averaged by hour over 5-day intervals from 2008 to 2012.

Figure 8 .
Figure 8. Monthly mean observed ( o ) and estimated ( e ) energy partitioning (LE / RN, H / RN and G / RN) and Bowen ratio (β) for all flux stations from 2008 to 2012 using α PTC of 0.92.
for the Heath flux station, with each diurnal segment representing flux data averaged by hour over 5-day intervals from 2008 to 2012.Observed and estimated R n exhibited an excellent agreement showing almost the same daily temporal pattern for the full growing season while LE, H and G yielded a good daily agreement being underestimated in May and September, especially in the case of LE.In terms of observed ( o ) and estimated ( e ) mean season energy flux partitioning, LE o / R No , H o / R No and G o / R No yielded mean values of 0.55, 0.37 and 0.08, respectively, and LE e / RN e , H e / R Ne and G e / R Ne yielded mean values of 0.58, 0.34 and 0.08, respectively (Fig.

Table 1 .
Flux station name and location, period of model evaluation and list of inputs required by the TSEB.Average and standard deviation for the input values were computed for the full period of model evaluation for each site.

Table 2 .
The Fen tower, located at the valley bottom in a wet sedge ecosystem, includes Eriophorum angustifolium and dwarf shrubs such as Betula nana and Salix spp., and vegetation types around the tower are comprised of 52 % wet sedge and 47 % tussock tundra.The Tussock tower, located at the midslope in a moist acidic tussock tundra ecosystem, is dominated by the tussock-forming sedge Eriophorum vaginatum, Sphagnum spp.and dwarf shrubs such as Betula nana and Salix spp.In this case, vegetation types around the flux tower Hydrol.Earth Syst.Sci., 21, 1339Sci., 21,  -1358Sci.,21,, 2017www.hydrol-earth-syst-sci.net/21/1339/2017/

Table 2 .
Trochim et al. (2016)12)Fen, Tussock and Heath flux sites instrumentation (more information is available at http://aon.iab.uaf.edu/imnavait).Apogee infrared radiometers were oriented 45 • off nadir at the three flux stations.Negative measurement of height indicate measurements of depth.An asterisk ( * ) indicates that this instrument is only available at the Tussock flux station.PAR indicates photosynthetically active radiation.are95%tussocktundra.The Heath tower sits atop a broad dry ridge at the top edge of the eastern watershed boundary in a heath tundra ecosystem dominated by dwarf shrubs and lichen.The vegetation here is 20 % heath, but also included 72 % tussock tundra, with the balance made up of sedge meadow and bare soil.Further detailed information about the study is provided inEuskirchen et al. (2012)andTrochim et al. (2016).

Table 3 .
Santanello and Friedl (2003)he soil heat flux estimation usingSantanello and Friedl (2003)(SF03) and Kustas et al. (1998) (K98) methodologies and two values for the maximum c G value.RMSE, MBE and MAD are in W m −2 and MAPD is in %.

Table 4 .
Accuracy statistic for the new c TG approach for the fit and the test.RMSE, MBE and MAD are in W m −2 , MAPD is in % and n is the number of half-hour intervals.

Table 5 .
Accuracy and error statistics from the comparison of modeled vs. observed unclosed and closed surface fluxes using α PTC of 0.92.n is the number of half-hour periods analyzed.RMSE, MAD and MBE are in W m −2 and MAPD is in %.

Table 6 .
Accuracy and error statistics from the comparison of modeled vs. observed unclosed and closed surface fluxes using α PTC of 1.26.n is the number of half-hour periods analyzed.RMSE, MAD and MBE are in W m −2 and MAPD is in %.

Table 7 .
Mean monthly accuracy and error statistics from the comparison of modeled vs. observed surface fluxes (using LE from residual closure) using α PTC of 0.92.n is the number of half-hour periods analyzed.RMSE, MAD and MBE are in W m −2 and MAPD is in %.