Challenges and uncertainty in plot-scale emissivity and surface temperature estimation using ux tower measurements

Land surface temperature (LST) is a preeminent state variable that controls the energy and water exchange between the Earth’s surface and the atmosphere. At the landscape-scale, LST is derived from thermal infrared radiance measured using space-borne radiometers. At the plot-scale, the ﬂux tower recorded longwave radiation components are inverted to retrieve LST. Since the down-welling longwave component was not measured routinely until recently, usually only the up-welling longwave component is used for the plot-scale LST retrieval. However, we found that ignoring reﬂected down-welling longwave radiation for plot-scale LST estimations can lead to substantial error. This also has important implications for estimating the correct surface emissivity using ﬂux tower measurements, which is needed for plot-scale LST retrievals. The present study proposes a new method for plot-scale emissivity and LST estimation and addresses in detail the consequences of omitting down-welling longwave radiation as frequently done in the literature. Our analysis uses ten eddy covariance sites with different land cover types and found that the LST values obtained using both up-welling and down-welling longwave radiation components are 0.5 to 1.5 K lower than estimates using only up-welling longwave radiation. Furthermore, the proposed method helps identify inconsistencies between plot-scale radiometric and aerodynamic measurements, likely due to footprint mismatch between measurement approaches. We also found that such inconsistencies can be removed by slight corrections to the up-welling longwave component and subsequent energy balance closure, resulting in realistic estimates of surface emissivity and consistent relationships between energy ﬂuxes and surface-air temperature differences. Landscape-scale day-time LST obtained from satellite data (MODIS TERRA) was strongly correlated with our plot-scale estimates for most of the sites, but higher by several Kelvin at two sites. We also quantiﬁed the uncertainty in estimated LST and surface emissivity using the different methods and found that the proposed method does not result in increased uncertainty. The results of this work have signiﬁcant implications for the combined use of aerodynamic and radiometric measurements to understand the interactions and feedbacks between LST and surface-atmosphere exchange processes. the observed input variables used for the estimation of plot-scale ε and LST has an associated uncertainty. Here we present exemplary results for Alice Springs, which showed the highest correlation between plot-scale and landscape scale LST


Introduction
The effects of global change are reflected in land surface temperature (LST) anomalies and their interannual variability 1 . It controls the magnitude and variability of the surface energy balance (SEB) components and simultaneously gets modulated by the SEB partitioning 2,3 . LST contains imprints of surface moisture and is extremely sensitive to evaporative cooling, which makes it a preeminent variable for studying evaporation and surface-atmosphere exchange [4][5][6] . It directly affects the amount of emitted longwave radiation and influences the saturation vapor pressure at the surface that drives latent heat flux. Thus, the ecohydrological functioning and carbon-water coupling are largely controlled by the surface temperature of the soil-vegetation system 7 . The availability of an extensive network of eddy covariance measurements (FLUXNET) allows us to understand the interactions and feedbacks between the surface-atmosphere exchange processes such as evaporation, transpiration, and its control by the atmosphere and vegetation at the diurnal time scale. However, the unavailability of direct LST measurements at the same scale hinders a detailed understanding of the interactions and feedbacks between LST and surface-atmosphere exchange processes, which is of utmost importance to the climate modeling community 8 In the last two decades, plot-scale radiometric data collected at eddy covariance sites (ECS) have gained popularity for in-situ LST retrieval due to its wide availability and high temporal resolution 9,10 . In addition to this, the LST estimates at plot-scale originate from a relatively homogeneous footprint in comparison to the satellite-derived LST (MODIS pixels). ECS measurements are primarily used to assess the impacts and feedbacks of climate change on key ecosystem fluxes 11 . By definition, LST is a thermodynamic temperature that can be felt or measured by an accurate thermometer at the land surface-atmosphere point-of-contact and is independent of wavelength 12 . The instantaneous value of LST is the result of interplay between the net radiation at the surface, ground heat flux (G), sensible heat flux (H) and latent heat flux (LE) 13 . Thus, LST can also be used for the estimation of H 14 and LE 15 between the surface and the atmosphere. LST provides the lower-boundary condition in SEB models for diagnostic estimates of LE and is highly relevant for drought monitoring 2,5,16 . Inversion of the longwave radiation in FLUXNET data to obtain LST has been found to strongly depend on the emissivity of the underlying surface 17 , which is not available as routine measurement. Therefore, estimating in-situ LST is not straightforward due to the involvement of two unknowns (LST and emissivity) inside one measurement variable (up-welling longwave radiation). To circumvent this challenge, we conducted simultaneous retrievals of LST and emissivity by exploiting the longwave radiation components in conjunction with associated SEB flux measurements.
The SEB components can be sub-divided into radiative components (often lumped in net radiation, R net ) and thermodynamic components, including sensible, latent and ground heat flux (H, LE, G respectively): As the surface-to-air temperature difference drives the exchange of sensible heat between surface and atmosphere, all components of Eq. (1) depend on the LST. Net radiation (R net ) can be sub-divided into down-welling and up-welling components. Only a fraction of solar top-of-the-atmosphere radiation reaches the Earth's surface, as some is reflected back to space by clouds, some is absorbed by the atmosphere and emitted later as longwave radiation. The emitted longwave radiation as a function of surface temperature (T s , K) and surface emissivity (ε) is given by Stefan-Boltzmann (SB) equation: where σ (W m −2 K −4 ) is the SB constant, ε is the surface emissivity ranging between 0 and 1, and T s (K) is the LST. Emissivity is defined as efficiency of a surface to emit thermal energy relative to a perfect black body. For a land surface, it depends on soil type, vegetation cover, soil moisture, soil chemistry, roughness, spectral wavelength, temperature and view angle 18 Putting the radiative components together, we can sub-divide R net into: Reflected shortwave in Eq. (3) is expressed as R sre f = αR sdown , where α is the surface albedo. Considering Kirchhoff's law, whereby the emissivity of a surface equals its absorptivity, emissivity values below unity result in reflected longwave radiation, expressed as: Radiometric temperature or LST is the "ensemble directional radiometric surface temperature" 18 , and can be estimated from the infrared radiance emanating from a given surface with known emissivity 19 . The emitted and down-welling longwave radiance are measured at given angle within its instantaneous field of view (fov) by a downward facing sensor relatively close to the surface (a few meters for an eddy covariance tower). The radiation received by a pyrgeometer or infrared sensor is a combination of the radiation emitted and reflected by the surfaces in its fov.
Substitution of Eqs. (4 and 2) into Eq. (5) yields R lup as a function of emissivity, surface temperature and down-welling longwave radiation: Eq. (6) is then solved for LST as a function of measured longwave and known surface emissivity: In order to invert LST as shown in Eq. (7), ε values are required. However, radiometers at ECS do not measure spectral bands separately to deduce emissivity directly. Therefore, we will deduce ε from observations of sensible heat flux (H), which is defined as the heat transfer driven by a surface-to-air temperature difference. It can be expressed mathematically in analogy to Ohm's law as: where T a (K) is the temperature of the air measured at a reference height above the surface, C p (J kg −1 K −1 ) is the specific heat capacity of air, ρ (kg m −3 ) is the air-density, and r a (s m −1 ) is the total resistance to heat transport from surface to the atmosphere. In simplified form, we write: where m (m s −1 ) is a proportionality constant (defined as m = ρC p /r a and broadly referred to as heat transfer coefficient) and depends on surface characteristics and micro-meteorology 20 . It is evident from Eq. (9) that for T s − T a = 0, H will be zero. This boundary condition and the linear relationship between H and ∆T has been used in the past to estimate ε at the plot-scale from observed H, T a and estimated T s using measured longwave radiation 21,22 . Another approach for plot-scale ε estimation filters the data where H is close to zero, substitutes T s in Eq. (6) by T a and solves for ε 23 . However, due to surface heterogeneity, sparse canopies are prone to footprint mismatch between the aerodynamic (flux tower) footprint and radiometric (hemispherical) footprint [24][25][26] , where the aerodynamic footprint represents the area contributing to measured sensible heat flux (10s to 100s of meters fetch), while the radiometric footprint is dominated by the surface below the sensor at a 60 o viewing angle, contributing to the measured longwave radiation (used for T s estimation). This can result in a different boundary condition i.e. at ∆T = 0, H = 0 as expressed in Eq. (10): where H is representative of the sensible heat flux from the eddy covariance tower footprint, T s is representative of all the radiating surfaces in the radiometric sensor's view, and c is interpreted as the H from surfaces in the aerodynamic footprint that are not seen by the radiometer. Plot-scale estimation of ε and LST using observed H, T a , R lup and R ldw as described above and in the Methods section, may be prone to substantial uncertainty. It is unclear how uncertainties in observed fluxes propagate into the uncertainty of estimated LST and ε. By design, infrared thermal (IRT) sensors only measure up-welling infrared radiance and therefore cannot explicitly account for the amount of reflected down-welling infrared radiation in the signal. For a long time, down-welling longwave R ldw was not routinely observed at ECS 27 and was also considered to be the most poorly quantified component of the radiation budget 28 . Therefore, the second term in Eq. (6) is commonly omitted, arguing that ε ≈ 1, and therefore Eq. (6) is simplified to Eq. (2) 29 : Eq. (11) can be solved for T s to yield what we will term the "short equation" (seq) for T s : Note that the above derivation is actually flawed, as the second term of Eq. (6) was omitted arguing that ε ≈ 1, and yet ε was retained in the first part of the equation. Nevertheless, even with the availability of down-welling longwave measurements 30 , the use of Eq. (11) is still a common practice 17,29 . This gives rise to the question if the short equation (Eq. (12)) is adequate to estimate LST from ground-based measurements. In the remainder of this paper, we will refer to LST obtained using the long equation (Eq. (7)) as T leq and to LST obtained using the short equation (Eq. (12)) as T seq .
To better understand and improve approaches of plot-scale LST estimation, the present study addresses the following research questions: 1. Can we obtain an adequate estimate of plot-scale LST while neglecting the reflected down-welling longwave radiation? 2. Does the estimation of plot-scale ε based on observed sensible heat flux (H) have an advantage over satellite-derived ε for plot-scale LST estimation?
3. How much uncertainty is introduced in plot-scale LST and ε due to uncertainty in measured EC fluxes?
To answer these questions, we analysed data for ten eddy covariance sites in different biomes and climates (see Table 2). Plot-scale broadband monthly emissivity was derived using observed H and estimated ∆T as proposed by Holmes et al. 21 . Plot-scale LST was estimated using plot-scale or landscape scale emissivity with (Eq. 7) and (Eq. 12). Plot-scale LST was compared with MODIS LST (TERRA satellite-sensed) for the times of satellite overpass. Uncertainty in ε and LST due to uncertainty in observed fluxes was calculated using SOBOL based uncertainty analysis (SAlib) 31 . See the Methods section for more details.

Plot-scale ε using long and short equation
Following the method proposed by Holmes et al. 21,22 , plot-scale monthly ε was estimated at the study site by fitting ε to minimise the root mean square error (RMSE) of the regression between H and T s − T a (SI Fig. 12). In Fig. 1a, c, and d, we reproduced the original data of Figs. 2a, 3C, and 3Q from Holmes et al. (2009) 21 to validate our interpretation of their approach using the short equation (Eq. (12)). We noted only marginal differences between the two results based on the short equation, which are likely due to different fitting algorithms. The replication of the H(∆T ) plot using the long equation (Eq. (7)) with the same data is given in Fig. 1b and the monthly ε values are shown in Fig. 1c, d, indicated by blue stars. The retrieved LST values were slightly higher when using Eq. (7) (compare a and b in Fig. 1). The use of the long equation (Eq. (7)) resulted in substantially (10%) lower values of ε as compared to the values estimated by Holmes et al. 21 for the common study sites (Brookings, Fig. 1c and Yatir, Fig. 1d).
Another approach for plot-scale ε estimation (Maes et al. (2019) 23 ) in combination with Eq. 7) resulted in even lower ε values for Brookings, as shown in Fig. 1c (red stars), whereas at Yatir, this approach gave an ε value higher than 1 (red star in Fig. 1d). Note that the long equation also yielded H(∆T ) relationship for many more months at Yatir Forest (blue star) than the short equation (black dots) as shown in Fig. 1d, as it resulted in achieving a strong correlation between H and ∆T (section 3 for details). The pattern of lower ε and higher LST using the long equation compared to the short equation was confirmed for all the ten sites used in the present study (Table SI 2

Landscape scale vs plot-scale estimates of ε and LST
At each site, LST was estimated using both the short equation (T seq , Eq. 12) and the long equation (T leq , Eq. 7). In the first step, tower-based longwave radiation and landscape scale broadband ε from MODIS spectral ε (ε MODIS , Eq. 14) was used.

4/19
The yearly daytime surface-to-air temperature difference for each study site is estimated and shown in Fig. 2. At all sites, Eq. (12) resulted in higher day-time plot-scale T s estimates as compared to Eq. (7), when using ε MODIS , with the medians of surface-to-air temperature differences (∆T ) differing by 0.8 to 1.5K (Fig. 2). The difference in ∆T using the two equations is highest at the water limited sites, e.g. AS and YA. Note that for two sites ( Table 2 for site abbreviations. revealed strong correlations between plot-scale and landscape scale LST estimates but systematically lower plot-scale LST ( Fig.  3a, b). Use of plot-scale ε plot for LST estimation (T seq and T leq ) resulted in substantial reduction of the bias as shown in Fig. 3c, d. This trend in bias reduction was similar at other sites (Table SI2 for details). The minimum bias is found at TUM, a closed canopy (eucalypt forest) and the highest bias was obtained at LF and HS, heterogeneous ecosystems with sparse canopies (woodland savanna). However, for some sites, weak correlation between satellite-derived and local LST estimates were also evident (at DU, R 2 was reduced from 0.8 to 0.4, see Table SI2). Also, using plot-scale ε for LST estimation resulted in positive T s − T a at LF and HS as shown in SI3, Fig

Plot-scale ε estimation using long equation with intercept
In order to account for the possibility of bias between radiometric and aerodynamic measurements (e.g. due to footprint mismatch of measuring devices or instrument bias) we also fitted Eq. 10, i.e. a relationship allowing for an intercept in the linear fit between H and ∆T (instead of forcing it through zero as in Fig. 1) for plot-scale ε estimation. As shown in Fig. 4, the plot-scale ε values resulting from this approach (H = m∆T + c) were substantially closer to the landscape-scale ε values compared with the approach without intercept (H = m∆T ), as shown in Table ( (Fig. 5a) without change in other regression paramaters (m, RMSE, R 2 ). In this study, we did not apply any energy balance closure scheme, as a Bowen ratio closure, although resulting in higher R 2 values at HS, also led to even greater intercept (c) (Fig. 5b). Interestingly, adding 40 W m −2 to the measured up-welling longwave radiation and subsequent energy balance closure largely removes the intercept and at the same time increases R 2 , as shown in (SI Fig.  13). Also the bias between MODIS and plot-scale LST is reduced from -10.66 K (Table.1 . In a first step, we used satellite-derived landscape scale broadband emissivity from MODIS (ε MODIS , Eq. 14) for estimating plot-scale LST from tower-based longwave measurements, and compared these with landscape-scale LST extracted from MODIS LST dataset (T MODIS ) .
estimations ( Table 1). The uncertainty in plot-scale ε estimated using Eq. 7 ('leq') and Eq. 9 (i.e. without intercept in H(∆T )) was mainly in the range of ±0.02 to ±0.05, with a maximum of ±0.2 if outliers are included (blue color in Fig. 6a). The short equation (Eq.12, 'seq') resulted in a vary narrow range of ε values between 0.94 and 0.99 throughout the year, with very small uncertainty (around ±0.01, black boxes in Fig. 6a). Interestingly, the differences in ε uncertainty did not propagate into differences in LST uncertainty, which were around ±0.2 K at the hourly scale for each equation if plot-scale emissivity was used (blue boxes in Fig. 6b and black boxes in Fig. 6c). In fact, if landscape-scale values of ε were used, the LST uncertainty was even bigger (±0.5 K, orange boxes in Fig. 6b and c). The comparison of uncertainty in epsilon and T s − T a using the long equation with and without intercept is shown in SI Fig. 14. The uncertainty in hourly T s − T a more than doubles if an intercept is allowed during the estimation of ε (SI Fig. 14b).

Discussion
Our analysis revealed a fundamental flaw in the commonly used short equation (Eq. (12)) for estimating plot-scale LST and ε, as it does not produce the same results as the long equation (Eq. (7)) even with high values of ε (MODIS ε). In fact, the short equation strongly over-estimates the sensitivity of LST to ε (SI Fig. 11), as it neglects the fact that low emissivity results in a greater fraction of reflected longwave in the sensor signal (compare Eq. (12) and (7)). The sensitivity of the long equation (Eq. (7)) to ε is driven by the contrast between R lup and R ldwn , whereas for the short equation (Eq. (12)), it is only driven by observed R lup (SI Fig. 9). For instance, an error of 0.01 in ε at a water-limited site (e.g. AS) can cause an error of 0.17 K using Eq. (7) and 0.79 K using Eq. (12) respectively (SI Fig. 11). This means that small errors in ε can result in large differences in LST when using the short equation, or conversely, unrealistic LST values can conveniently be rectified by slightly changing the  Same analysis and legends as in Fig. 4c), but (a) After adding 40 (W m −2 ) to measured R lup , and (b) after closing the energy imbalance using a Bowen ratio closure scheme.
ε value. This is illustrated e.g. in Fig. 6, where estimation of plot-scale ε resulted in similar LST values between the short and long equations, but with vastly different ε values and much greater uncertainty in estimated ε using the long equation compared to the short equation. Considering that the short equation ignores an important component of longwave radiation, it must be concluded that in this case, it achieves seemingly the right results for the wrong reasons. The reduced sensitivity of the long equation (Eq. (7)) to ε is of advantage for plot-scale LST estimation, since plot-scale ε is usually unknown and therefore used as an approximate value 17 Table 2.
linearfit was allowed to have an intercept (Table 1, plot-scale ε). The intercept (i.e. ∆T = 0 at H = 0) could be caused by combining measurements coming from instruments (radiometer, eddy covariance system) with different footprints 25 . The mismatch of source areas becomes important if the surface underlying the instruments has heterogeneous land cover. Although "footprint awareness" is often omitted at ECS under the assumption of homogeneity 24 , in patchy vegetation, the radiometer can be "seeing" a different vegetation fraction than that contributing to EC measurements, meaning that H = 0 at ∆T = 0. This problem was not detected by Holmes et al. 21 , as the short equation (Eq. (12)) was used, and due to its high sensitivity to ε (SI Fig. 11a) even a small reduction in ε corrected the offset in H(∆T ) (Fig. 1a). In contrast, when repeating the same analysis using the long equation (Eq. (7)), a larger reduction in ε is required to remove the intercept, resulting in lower ε (Fig. 1b). By allowing an intercept in the H(∆T ) linear fit, we implicitly account for the possibility of a footprint mismatch or instrument bias in the data. This small change in methodology enables us to detect such problems by inspecting the value of the intercept (c). Considering the aerodynamic footprint to be larger than the radiometric footprint 24, 25 , a positive intercept can be interpreted as the H from the aerodynamic footprint which is not seen by the radiometer. The intercept was very high for the sites HS and LF (Table. 1). A close inspection of the H(∆T ) plots at these sites (SI Fig.  8) revealed negative day-time T s − T a (Fig. 2), which may suggest an underestimation of R lup . While testing this hypothesis at HS (having the highest intercept, Fig. 4c) we found that adding roughly 40 W m −2 (approx. 8% of observed R lup , Fig. 5) in observed R lup led to significant reduction in the intercept from 294 W m −2 (Fig. 4c) to 17 W m −2 and positive day-time T s − T a (Fig. (5a)). The other linear regression parameters (m, R 2 , RMSE) were not affected ( compare Fig. (5a and 4c). The hemispherical view of the radiometers looking down at a heterogeneous canopy makes it possible that they "see" more tree crowns and less soil than the area contributing to the eddy covariance footprint. This could lead to an underestimation of R lup , and an underestimation by 30 − 40 W m −2 would be equivalent to approximately 5-10% of the observed flux, which is within the range of a typical energy imbalance found at this site (SI Fig. 13). Previous studies have found a dependence of footprint mismatch on wind direction 24-26 , but we did not find a significant relation between monthly intercept and dominant wind direction at Howard Springs.
Surface heterogeneity has also been recognized as one of the potential causes for the lack of energy balance closure observed at most ECS at diurnal scales 9, 34 . However, in our analysis the use of an energy balance closure scheme (based on the Bowen ratio) led to much lower values of plot-scale ε using Holmes approach with the long equation and without intercept. In contrast, if an intercept was allowed, energy balance closure led to an increase in positive intercept (Fig. 5b). Perhaps this is the reason why other studies on plot-scale ε estimation have also used the observed fluxes without correction [21][22][23] .
However, our analysis suggests that the footprint mismatch may cause a small bias in the up-welling longwave radiation measurements that is not accounted for in conventional energy balance closure approaches. When we added 35 W m −2 (instead of 40 W m −2 , see Fig. 5a) to the measured up-welling longwave radiation and subsequently closed the energy balance at the HS site (which had the largest H(∆T ) intercept), we largely removed the intercept and at the same time obtained realistic ε values and an increased R 2 (SI Fig. 13a). In addition, the bias between MODIS LST and plot-scale LST at HS was reduced by 6.4 K (SI Fig. 13b) compared to using up-welling longwave without correction. Uncertainty in monthly plot-scale emissivity due to uncertainty in H, R lup , R ldw and T a , using Eq. 6 ('leq', blue) and Eq. 11 ('seq', black). (b) Hourly uncertainty in T s − T a on 15 July based on Eq. 6, due to uncertainty in R lup , R ldw and T a when landscape-scale emissivity is used (ε MODIS , orange) or due to uncertainty in H, R lup , R ldw and T a when plot-scale emissivity is used (ε plot , blue). (c) Same as Panel b, but based on Eq. 11.
When estimating plot-scale LST using landscape-scale ε values, we found at many sites with a sparse canopy strongly negative bias in comparison to MODIS LST, which is in agreement with previous studies where the bias for sparse canopies reached up to 12 K 35 . Plot-scale LST estimates based on plot-scale ε using a linear H(∆T ) fit without an intercept largely reduced this bias between plot-scale and MODIS LST (Table 1) and also reduced the uncertainty in diurnal LST (Fig. 6b,c) in comparison to the use of landscape-scale ε. However, the resulting plot-scale ε values were unrealistically low at some sites (Table 1, center). In contrast, allowing an intercept (H = m * ∆T + c) in plot-scale ε estimation resulted in more realistic ε values at these sites, but very large intercept values (over 200 W m −2 at some sites), indicating that the plot-scale LST values cannot be used in combination with the observed aerodynamic fluxes at these sites, as strongly positive H at 0 surface-air temperature difference is physically inconsistent (Fig.  4c). In addition, this approach increased the bias between plot-scale and MODIS LST at most of the study sites (Table 1). At the sites with the largest intercept values, we found that an assumed bias in up-welling longwave radiation by only 6-9% would largely remove the intercept and also reduce the bias between MODIS and plot-scale LST (Fig. 5a, Fig. 13b). A detailed analysis of such bias and potential correction approaches is beyond the scope of this study. Given that the fit of a linear model without intercept is statistically questionable in general 36 , and the fact that such a fit resulted in unrealistically low values of ε at some sites, we conclude that fitting a model with intercept is the more robust approach, and that a significant intercept should be used as a red flag for the utility of the data for estimation of plot-scale LST.
Note that the fluxes observed at ECS are representative of the composite signal from both, soil and vegetation, which typically have a different range of surface temperatures and emissivities 37 . The ε of soil strongly depends on soil moisture 9/19 content 38 , whereas the ε of a canopy depends on its structural attributes and leaf area index, the latter of which can vary strongly at the seasonal scale 39 . For example, the laboratory-measured directional ε for various canopy elements (bark, leaf and its arrangement, stem wood) ranged between 0.9 to 1 40 . Laboratory measurements of thermal infrared reflectance spectra suggest that the ε uncertainty due to structural unknowns, such as leaf orientation, is more significant than the differences in leaf component emissivity among plant species 33 . Consequently, it is clear that the ε of a surface is a function of many factors and detailed study of all these factors is out of scope of the present study. Derivation of landscape-scale broadband ε land from narrowband spectral emissivity is a first-order approximation for capturing the integrated effects of land cover from MODIS spectral bands 37 , whereas the derivation of plot-scale ε from EC flux data, as analyzed in this study, provides an independent, alternative estimate of ε and at the same time provides a quality check on the correspondence between radiometric and aerodynamic measurements at a given ECS.
In summary, our results reveal that the short equation (Eq. (12), neglecting down-welling longwave radiation) leads to biased estimates of LST and hugely over-estimated sensitivity of LST to surface emissivity. Therefore, the use of Eq. (12) is not recommended and should be replaced by Eq. (7) if down-welling longwave radiation measurements are available. At some sites, the use of Eq. (7) resulted in plot-scale LST estimates that were far below satellite-derived landscape-scale LST values, and also inconsistent with plot-scale flux data (negative surface-air temperature difference when sensible heat flux is strongly positive). In many previous studies, such bias would have been removed by slightly lowering surface emissivity (ε), but the reduced sensitivity of Eq. (7) to ε would require unrealistically low values of ε to remove the low-bias in LST. When estimating plot-scale ε values, realistic estimates based on Eq. (7) are only possible at these sites if we include an intercept in the H(∆T ) relationship, but this again results in very high intercept values (over 200 W m −2 ). However, assuming that the intercept is a consequence of a foot-print mismatch between the aerodynamic and radiometric measurements, a small correction in up-welling longwave (6-9%) and subsequent energy balance closure largely removed the intercept and and produced realistic ε values and self-consistent H(∆T ) plots. This approach also reduced the bias between plot-scale LST and MODIS LST, although it did not improve the weak correlation between them (SI Fig. 13b). The intercept value can be used as a criterion to check the consistency of observed data (radiometric and aerodynamic measurements) before using them in combination, as a strong intercept indicates inconsistency between observed sensible heat flux and surface-to-air temperature difference. Therefore, the proposed method of fitting a linear relation with intercept to H and ∆T data has the potential to provide more robust benchmark data sets for model evaluation and validation at the plot-scale. Overall, the implications of our study are of particular relevance for the research community interested in understanding the diurnal and seasonal feedback in soil-vegetation systems based on fluxes, emissivity and LST estimated at the plot-scale.

Research data
Tower data: ECS collect micro-meteorological measurements above the surface (vegetation canopy) using towers (flux tower) following common measurement protocols 11 . The towers are generally equipped with an instrument made up of pyrgeometers or radiometers to measure up-welling and down-welling shortwave and longwave radiation, which is further used to calculate net radiation (Eq. 3). Besides radiative fluxes, measurement at ECS also include sensible and latent heat fluxes, net carbon-dioxide exchange and a range of meteorological variables, such as air temperature, humidity and wind speed. T a is the air temperature measured at a reference height above the canopy. Each flux measurement is accompanied by a flagging system based on the second CarboEurope-IP QA/QC workshop 41 . In our current work we use only high quality (flag 0) data. For the analysis, ten sites were selected to represent a variety of land cover types and climates (Table 2). Eight sites belong to the North Australian Tropical Transect (NATT) and two sites (Yatir Forest, Brookings) are chosen to replicate results from Holmes et.al 21 . Eddy covariance level 3 data is obtained from http://data.ozflux.org.au/portal/pub/listPubCollections.jspx for Australian sites. The data for Brookings was obtained from Ameriflux whereas the data for Yatir Forest was obtained through personal communication with Professor Yakir's lab in order to obtain the older version of the data, which was used by Holmes et al. 21 MODIS data: Landscape-scale emissivity and LST data (MODIS product MOD11A1) was downloaded from NASA earth data . It is a level 3 daily LST product gridded in the sinusoidal projection at a spatial resolution of 0.928 km by 0.928 km. The daily LST pixel values in each granules (tile contains 1200 x 1200 grids in 1200 rows and 1200 columns) is retrieved by the generalized split-window algorithm under clear-sky conditions and MODIS LST values are averaged by overlapping pixels in each grid with overlapping areas as weight 42 . The downloaded data in hierarchical data format (hdf), were converted into tagged image file format (tiff) using a python package called PyModis 43 . From tiff the files are converted into csv format. The details of data extraction and conversion can be found at https://renkulab.io/projects/gitanjali.thakur/ modis_lst_fpar/. Thermal remote sensing (MODIS) is used to measure spectral emissivity through four channels (28,29,30,31) at wavelengths ranging between 8-12 µm 37 and the system of equations is iteratively solved for a given range of wavelengths (9 -12µm) to obtain ε and LST using radiative transfer models 27,37,44 In the current study, dataset columns used to compare plot-scale LST are: day time daily LST and local view time. In order to obtain landscape-scale ε the emissivity from bands 31 and 32 are used. These bands have more stable emissivities ranging from 0.92-1, thus used to derive broadband emissivity 42 . Plot-scale emissivity estimation: This approach was initially proposed by Holmes 21 to estimate plot-scale ε using the short equation (Eq. 12) with H, R lup , T a . In the present work, we have used both the long equation (Eq. 7) and short equation (Eq. 12) to estimate ε. The prime variables used in the study are H, R lup , R ldwn , and, T a , whereas the ancillary variables R n and wind speed (W s) are used to filter the data for analysis. The data filtering criteria are (R n > 25 W m −2 ) and wind speed (W s > 2ms −1 ) 21 . Data (R lup , R ldwn , H, T a ) satisfying the filtering criteria is segregated monthly. Assuming an initial ε value, observed longwave radiation (R lup , R ldwn ) is used to estimate T s for each valid data point within a month. For each month, a linear regression (with or without intercept, see main text) is performed between sensible heat (H) and surface-air temperature difference (T s − T a ) (Fig. 7b) using Scipy https://docs.scipy.org/doc/scipy-0.16.0/ reference/generated/scipy.stats.linregress.html. This procedure is repeated using a predefined range of ε from 0.4 until 0.998 with a step size of 0.002. From all the values of ε that resulted in strong correlation between H and T s − T a (R 2 > 0.5), the value that achieved the smallest RMSE was chosen as the suitable ε for the given month. An illustration plot for RMSE and ε is shown in SI Fig. 12. The ε is obtained for each month using both, long (Eq. (7)) and short equation (Eq. (12)) and termed as ε leq and ε seq respectively, as shown in Fig. 7b. For sites like HS and LF resulting in a high value of intercept, we hypothesized an underestimation of R lup and corrected it by adding 6 − 8% of the observed day time observed R lup to the measured data before closing the energy balance using Bowen ratio closure, and estimating monthly plot-scale ε.
Recently, another approach for plot-scale ε estimation using Eq. (5) was discussed by Maes et.al (2019) 23 . For each month, the observed R lup , R ldwn , T a are filtered for non rainy days using the filtering criteria (−2 < H < 2 , and α < 0.4). The ε is estimated by substituting T s = T a in Eq. (5) as shown in Eq. (13). The monthly ε was represented as the median of ε obtained by substituting filtered data in Eq. (13).
LST comparison: LST estimates based on optimum ε for each month using a linear regression with and without intercept H(∆T ) is compared to the MODIS daily LST. For each month the T s estimates using plot-scale ε and MODIS ε corresponding to the exact TERRA day-time of pass was obtained (30 minute tower data was upsampled to each minute using linear interpolation). TERRA (satellite) overpasses at local solar time between 10:30 am to 12 pm in ascending mode 12 . Plot-scale LST at local satellite overpass ( T leq and T seq ) was obtained using measured longwave radiation that was interpolated to match the exact time of overpass and corresponding monthly plot-scale ε. Plot-scale daily LST is compared to MODIS LST in terms of mean, bias, RMSE and R 2 using a robust linear regression model (scipy stat model,https: //github.com/scipy/scipy/blob/v1.7.1/scipy/stats/_stats_mstats_common.py as shown in Fig.  7a. The goodness of fit between plot-scale and landscape-scale LST was determined by looking at R 2 (Fig. 7b). The bias is estimated as the mean of deviation between daily MODIS LST and ground-based T s (T leq,seq − T MODIS ). See SI table for data sources and acronyms in SI Table 3.
General approach: We estimate landscape-scale broadband ε using MODIS spectral ε as shown in Bahir et. al (2017) 45 .
ε MODIS = 0.4587ε 31 + 0.5414ε 32 (14) Tower-based longwave radiation measurement (R lup , R ldwn ) passing the filtering criteria (as mentioned in plot-scale emissivity estimation) along with MODIS based ε was used to invert LST using Eq. (7) and Eq. (12). The obtained plot-scale LST was compared to landscape-scale MODIS LST using a robust linear regression as mentioned above and shown in Fig. 7a. Uncertainty estimation: Uncertainty in plot-scale ε and LST was quantified based on an assumed systematic error (caused by a potential bias in measurement devices) at the study sites. In a first step, based on the literature 28, 46 the error bounds of each input variable (H, R lup , R ldw , T a ) used for plot-scale ε estimation were defined. The error bounds for R lup and R ldwn are -5 to 5 W m −228 , for H, we used -20 to 20 W m −2 and for T a we used -1 to 1 K 46 . The error samples (perturbation) within these bounds were generated using the Saltelli sampling scheme (using python based package SALIB 47 ). Each error sample is added to the monthly segregated measured fluxes as explained above. Observed fluxes combined with perturbed fluxes are used to estimate T s using Eq. (12) and Eq. (7). The obtained range of diurnal T s and observed T a with perturbation is used to calculate the uncertainity in ∆T ; an example for July 15 is shown in Fig. 6c. Perturbed sensible heat flux(H+sample error) and perturbed ∆T is used to obtain ε as described in plot-scale ε estimation. The range of monthly plot-scale ε obtained is reported as uncertainty in monthly ε.

Author contributions statement
GT, SJS, KM, and IT designed the study. GT conducted the analysis and prepared the manuscript together with SJS. KM, IT and MS provided textual input, corrections, suggestions and critical discussion of the manuscript.

Code and data availability
The data and code used for this study is available on renkulab https://renkulab.io/projects/gitanjali. thakur/plot-scale-lst-emissivity. Before final publication, static version of the repository will be uploaded to zenodo.org and receive a separate DOI  Figure 7. Schematic representation of steps followed for plot-scale LST retrieval using eddy covariance measurement (a) Landscape emissivity and longwave measurement and compared to Landscape-scale LST (T MODIS ) (b) Plot-scale emissivity estimation using observed H, R ldwn , R lup and plot-scale LST is estimated using plot-scale ε is compared to landscape-scale emissivity. The R 2 , RMSE, Bias are mentioned in Fig. (3).