A coupled ground heat flux–surface energy balance model of evaporation using thermal remote sensing observations

One of the major undetermined problems in evaporation (ET) retrieval using thermal infrared remote sensing is the lack of a physically based ground heat flux (G) model and its integration within the surface energy balance (SEB) equation. Here, we present a novel approach based on coupling a thermal inertia (TI)-based mechanistic G model with an analytical surface energy balance model, Surface Temperature Initiated Closure (STIC, version STIC1.2). The coupled model is named STIC-TI. The model is driven by noon–night (13:30 and 01:30 local time) land surface temperature, surface albedo, and a vegetation index from MODIS Aqua in conjunction with a clear-sky net radiation sub-model and ancillary meteorological information. SEB flux estimates from STIC-TI were evaluated with respect to the in situ fluxes from eddy covariance measurements in diverse ecosystems of contrasting aridity in both the Northern Hemisphere and Southern Hemisphere. Sensitivity analysis revealed substantial sensitivity of STIC-TI-derived fluxes due to the land surface temperature uncertainty. An evaluation of noontime G (G_i) estimates showed 12 %–21 % error across six flux tower sites, and a comparison between STIC-TI versus empirical G models also revealed the substantially better performance of the former. While the instantaneous noontime net radiation (R_Ni) and latent heat flux (LE_i) were overestimated (15 % and 25 %), sensible heat flux (H_i) was underestimated (22 %). Overestimation (underestimation) of LE_i (H_i) was associated with the overestimation of net available energy (R_Ni−G_i) and use of unclosed surface energy balance flux measurements in LE_i (H_i) validation. The mean percent deviations in G_i and H_i estimates were found to be strongly correlated with satellite day–night view angle difference in parabolic and linear pattern, and a relatively weak correlation was found between day–night view angle difference versus LE_i deviation. Findings from this parameter-sparse coupled G–ET model can make a valuable contribution to mapping and monitoring the spatiotemporal variability of ecosystem water stress and evaporation using noon–night thermal infrared observations from future Earth observation satellite missions such as TRISHNA, LSTM, and SBG.

1 Introduction

One of the outstanding challenges in evaporation (ET) estimation through surface energy balance (SEB) models concerns an accurate characterization of ground heat flux in the open canopy architecture with mixed vegetation such as savanna or in ecosystems with low mean fractional vegetation cover, prevailing water stress, and strong seasonality in soil moisture. Ground heat flux (G) is an intrinsic component of SEB (Sauer and Horton, 2005), affecting the net available energy for ET (the equivalent water depth of latent heat flux, LE) and sensible heat flux (H). It represents an energy flow path that couples the surface with the atmosphere and has important implications for the underlying thermal regime (Sauer and Horton, 2005). Depending on the vegetation fraction and water stress, the magnitude of instantaneous G varies greatly across different ecosystems. In the humid ecosystems with predominantly dense canopies and high mean fractional vegetation cover, G contributes a small proportion to the surface energy balance. Dense canopy cover leads to less transmission of radiative fluxes through multiple layers of canopies, which results in low warming of the soil floor. Due to persistently high soil water content, humid ecosystems generally show low diurnal and seasonal variability in G. In contrast, the magnitude of G is substantially large in arid and semi-arid ecosystems with sparse and open canopy and soil moisture deficits. Due to the prevailing feedback between the physics of ground heat flux, land–atmosphere interactions, and vegetation ecophysiology, evaporation modeling in the complex ecosystems remains a challenging task (Wang et al., 2013; Kiptala et al., 2013). This paper addresses the challenge of simultaneous estimation of G and ET by combining thermal remote sensing observations with a mechanistic G model and a SEB model.

SEB models mainly emphasize estimating sensible heat flux by resolving the aerodynamic conductance (g_A) and computing LE as a residual SEB component as follows.

$$\begin{array}{}\text{(1)}& \mathit{LE}=R_{\mathrm{N}}-G-H\end{array}$$

R_N is the net radiation. The proportion of R_N that is partitioned into G depends upon soil properties like its albedo, soil moisture, and soil thermal properties such as thermal conductivity and heat capacity, which vary with mineral, organic, and soil water fractions. SEB models use land surface temperature (LST or T_S) as an important lower boundary condition for estimating H and LE. Due to the extraordinarily high sensitivity of T_S to evaporative cooling and soil water content variations, thermal infrared (TIR) remote sensing is extensively used in large-scale evaporation diagnostics (Kustas and Anderson, 2009; Mallick et al., 2014, 2015a, 2018a; Cammalleri and Vogt, 2015; Anderson et al., 2012). Evaporation estimation through SEB models commonly employs empirical sub-models of G in a stand-alone mode. Despite the utility of mechanistic G models being demonstrated in different studies (Verhoef, 2004; Murray and Verhoef, 2007; Verhoef et al., 2012), no TIR-based evaporation study attempted to couple a mechanistic G model with a SEB model.

The SEB models for ET estimation driven by remote sensing observations generally use linear and nonlinear relationships for estimating G, and such methods generally employ R_N, T_S, albedo (α_R), and NDVI (e.g., Bastiaanssen et al., 1998; Friedl et al., 2002; Santanello and Friedl, 2003). While the inclusion of T_S and albedo serves as a proxy for soil moisture and surface characteristic effects in G, inclusion of NDVI provides a scaling of the G–R_N ratio for different fractional vegetation covers. Unfortunately, the empirical approaches do not include any information on soil temperature or daily temperature amplitude. These empirical models also lack universal consensus. Setting G as a fraction of R_N does not solve the energy balance equation and disregards the role of thermal inertia of the land surface (Mallick et al., 2015b). This could introduce substantial uncertainty in LE estimation because G effectively couples the surface energy balance with energy transfer processes in the soil thermal regime. It provides physical feedback to LE through the effects of soil moisture, temperature, and conductivity (thermal and hydraulic) (Sauer and Horton, 2005). Such feedbacks are most critical in the arid and semi-arid ecosystems, where LE is significantly constrained by the soil moisture dry-down. The limits imposed on LE by the water stress consequently result in greater partitioning of the net available energy (i.e., R_N−G) into H and G (Castelli et al., 1999).

When LE is reduced due to soil moisture dry-down, both G and T_S tend to show rapid intra-seasonal rise. Therefore, the surface energy balance equation could be linked with the mechanistic G model, T_S harmonics (Verhoef, 2004), and soil moisture availability. Realizing the importance of direct estimates of G in LE, and invigorated by the advent of TIR remote sensing, Verhoef et al. (2012) demonstrated the potential of a TI-based mechanistic model (Murray and Verhoef, 2007) (MV2007 hereafter) for spatiotemporal G estimates in semi-arid ecosystems of Africa. Some studies also emphasized the importance of using noontime and nighttime T_S and R_N for estimating G (Mallick et al., 2015b; Bennet et al., 2008; Tsuang, 2005). The method of MV2007 has so far been tested in a stand-alone mode, and no remote sensing method has so far been attempted to combine such a mechanistic G model (e.g., the MV2007-TI model) with a SEB model for coupled energy–water flux estimation and validation.

By integrating T_S into a combined structure of the Penman–Monteith (PM) and Shuttleworth–Wallace (SW) models, an analytical SEB model was proposed by Mallick et al. (2014, 2015a, 2016). The model, Surface Temperature Initiated Closure (STIC), is based on finding an analytical solution for aerodynamic and canopy–surface conductance (g_A and g_S), where the expressions of the conductances were constrained by an aggregated water stress factor. By physically linking water stress (T_S-derived) with g_A and g_S, STIC established direct feedback between T_S, H, and LE and simultaneously overcame the need for empirical parameterization to estimate the conductances (Mallick et al., 2016, 2018a). Different versions of STIC have been extensively validated in different ecological transects (tropical rainforest to woody savanna) and aridity gradients (humid to arid) (Trebs et al., 2021; Bai et al., 2021; Mallick et al., 2015a, 2016, 2018a, b; Bhattarai et al., 2018, 2019). Based on the conclusions of Verhoef et al. (2012), Mallick et al. (2014, 2015a, b, 2016, 2018a, b, 2022), Bhattarai et al. (2018, 2019), and Bai et al. (2021), there is a need to address some of the challenges in SEB modeling, which are (i) accurate estimation of G and ET in sparse vegetation, (ii) testing of the utility of coupling a TI-based G model with an analytical SEB model for accurately estimating G and ET, and (iii) detailed evaluation of a coupled G–SEB model at the ecosystem scale. Realizing the significance of the mechanistic G model (MV2007) and the advantage of the analytical STIC model and to mitigate some of the overarching gaps in SEB modeling in sparsely vegetated open canopy systems, this study presents the first-ever coupled implementation of MV2007 G with the most recent version of STIC (STIC1.2). We name this new coupled model STIC-TI, and it requires noon–night T_S and associated remotely sensed land surface variables as inputs. We performed subsequent evaluation of STIC-TI in nine terrestrial ecosystems in arid, semi-arid, and sub-humid climates in India, the United States of America (USA) (Northern Hemisphere), and Australia (Southern Hemisphere) at the eddy covariance flux tower sites. The current study addresses the following research questions and objectives.

• i.

  What is the performance of STIC-TI G estimates when compared with conventionally used empirical G models in ecosystems with low mean fractional vegetation cover (f_c) (≤0.5) and with larger soil exposure to radiation, for example, in savanna?

• ii.

  How do the estimates from STIC-TI LE and H fluxes compare with LE and H observations in diverse terrestrial ecosystems that represent a varied range of f_c (0.25–0.5) covering cropland, savanna, and mulga vegetation (woodlands and open forests dominated by the mulga tree – Acacia aneura) spread across arid, semi-arid, sub-humid, and humid climates over a vast range of rainfall (250 to 1730 mm), temperature (−4 to 46 ^∘C), and soil regimes?

• iii.

  What is the seasonal variability of G and evaporative fraction from the STIC-TI model in a wide range of ecosystems having contrasting aridity and vegetation cover?

It is important to mention that assessing the performance of STIC-TI LE and H with respect to other SEB models is not within the scope of the present study. The prime focus of the current study is to assess the sensitivity of STIC-TI, temporal variability of the retrieved SEB fluxes, and cross-site validation of the individual SEB components.

A list of the variables, their symbols, and the corresponding units is given in Appendix A.

2 Study area and datasets

2.1 Study site characteristics

The present study was conducted using data from nine flux tower sites (four sites in India; three sites in Australia; two sites in the USA) equipped with eddy covariance (EC) measurement systems. The distribution of the flux tower sites considered for the present study are shown in Fig. 1 below. The sites cover a wide range of climates, vegetation types, a low fractional vegetation cover (f_c) of around 0.5, and contrasting aridities (Table 1). In India, a network of EC towers was set up under the Indo-UK INCOMPASS (INteraction of Convective Organization and Monsoon Precipitation, Atmosphere, Surface and Sea) program (Turner et al., 2019) at Jaisalmer (IND-Jai) in Rajasthan state, Nawagam (IND-Naw) in Gujarat state, and Samastipur (IND-Sam) in Bihar state and under the Newton Bhabha program (Morisson et al., 2019a, b) at Dharwad (IND-Dha) in Karnataka state. The flux footprint for EC towers in India varied from 500 m to 1 km (Bhat et al., 2019). In the present study, about 90 % of the fluxes came from an area within 500 m to 1 km of the EC tower. Therefore, the relative contribution of vegetated land surface area to the fluxes is close to 90 % (Schmid, 2002; Vesala et al., 2008). The remaining percentage of fluxes originated from an area beyond the flux footprint. The mean annual f_c was found to vary from 0.25 to 0.52, with a standard deviation (SD) ranging from 0.1 to 0.16.

The IND-Jai site represents the arid western zone over desert plains of the natural grassland ecosystem. The region receives very low rainfall (100–300 mm) during monsoon and experiences a wide range in air temperature, high solar radiation, wind speed, and high evaporative demand (Raja et al., 2015). The IND-Naw site represents a semi-arid agroecosystem in the middle Gujarat agroclimatic zone of northwestern India and has a pre-dominant rice–wheat cropping system. The IND-Sam site has a sub-humid climate of the northwestern alluvial plain zone in the Indo-Gangetic Plain (IGP) situated in eastern India, and this site also follows rice–wheat crop rotation. IND-Dha represents the humid sub-tropical climate of the transition zone in southern India, and this site is comprised of crops.

Figure 1 Locations of the flux tower sites in India, Australia, and the USA overlaid on a climate-type map. Image source: Peel et al. (2017).

In the USA, two EC tower sites were located at Tonzi Ranch (US-Ton) and Vaira Ranch (US-Var), in the lower foothills of the Sierra Nevada. Both EC stations are part of the AMERIFLUX Management Project (https://ameriflux.lbl.gov/, last access: 3 December 2022). US-Ton is classified as an oak savanna woodland. While the overstorey is dominated by blue oak trees (40 % of the total vegetation) with intermittent grey pine trees (3 trees ha⁻¹), the understorey species include a variety of grasses and herbs. The mean annual rainfall at this site is 559 mm. US-Var is a grassland-dominated site, and the growing season is confined to the wet season only, typically from October to early May. The mean annual rainfall at this site is 559 mm. The mean annual f_c was found to vary from 0.18 to 0.26, and SD is of the order of 0.06 to 0.07.

In Australia, three EC tower sites were located at Howard Springs (AU-How), Alice Springs Mulga (AU-ASM), and Adelaide River (AU-Ade) in the Northern Territory as part of the OzFlux network (Beringer et al., 2016) and the Terrestrial Ecosystem Research Network (TERN), which is supported by the National Collaborative Infrastructure Strategy (NCRIS) (http://www.ozflux.org.au/monitoringsites/index.html, last access: 3 December 2022). AU-How is situated in the Black Jungle Conservation Reserve representing an open woodland savanna, and the mean annual rainfall is 1750 mm. The AU-ASM is located at Pine Hill cattle station near Alice Springs. The woodland is characterized by a mulga canopy, and mean annual rainfall is 306 mm. AU-Ade represents savanna with a mean annual rainfall of 1730 mm. The mean annual f_c varied from 0.21 to 0.48 with an SD range of 0.08–0.17. A description of Australian flux sites is given in Beringer et al. (2016). Average heights of vegetation are 1.15 m at IND-Naw, 1 m at IND-Jai, 1.23 m at IND-Sam, 1.5 m at IND-Dha, 6.5 m at AU-ASM, 15 m at AU-How, 7 m at AU-Ade, 10 m at US-Ton, and ≤0.5 m at US-Var.

Table 1 An overview of the eddy covariance flux tower site characteristics used in the present study.

T_A: air temperature during the observation period. P: rainfall (mm) measured using the rain gauge at the flux tower site during the study period. IND is for India, AU is for Australia, and US is for the United States. SD is the standard deviation of the annual mean f_c which is computed from the NDVI as mentioned in Sect. 3.1.

2.2 Datasets

2.2.1 Micrometeorological data at flux tower sites

Standardized, controlled, and harmonized SEB flux and meteorological data from nine EC towers were used in the present analysis. In Australia, H and LE were measured through the EC systems and R_N was measured through net radiometers at varying heights of 15 m (AU-Ade), 23 m (AU-How), and 11.6 m (AU-ASM), respectively. In India, the EC measurement height was maintained approximately at 8 m above the surface, except at IND-Dha, where it was installed at a height of 4.2 m. In the USA, the SEB measurements were carried out at tower heights of 23 m at US-Ton and 2 m US-Var. A summary of the instrumentation is given in Table A2 of Appendix A. All the flux tower sites were equipped with a range of meteorological instrumentation which measured diurnal air temperature (T_A) and relative humidity (R_H), four components of the net radiation (R_N, consisting of downwelling and upwelling shortwave and longwave radiation: SW ↓, SW ↑, LW ↑ and LW ↓, respectively) above the vegetated canopy. In addition, the diurnal soil heat flux (G) and soil temperature (T_ST) were measured at all three Australian sites and both US sites. In India, the diurnal soil heat flux was measured only at IND-Dha.

For the Indian sites, the raw EC measurements of the turbulent wind vectors (u, v, and w, for horizontal, meridional, and vertical, respectively), sonic temperature (T), and CO₂ and water vapor mass density were recorded at a sampling rate of 20 Hz. Raw EC data were post-processed to obtain level-3 quality-controlled and harmonized surface fluxes at 30 min flux-averaging intervals using the EddyPRO^® Flux Calculation software (LI-COR Biosciences, Lincoln, Nebraska, USA) using the data-handling protocol described by Bhat et al. (2019). The EC data from the OzFlux sites were averaged over 30 min recorded by the logger and processed through levels using the PyFluxPro standard software processing scripts as mentioned in Isaac et al. (2017). Level 3 (L3) used in this paper was produced using PyFluxPro (Isaac et al., 2017) employing Dynamic INtegrated Gap filling and partitioning for the Ozflux (DINGO) system as described in Donohue et al. (2014) and Beringer et al. (2016). The quality-checked EC data at 30 min intervals for the two AMERIFLUX sites US-Ton and US-Var were acquired from https://doi.org/10.17190/AMF/1245971 and https://doi.org/10.17190/AMF/1245984, respectively.

2.2.2 Remote sensing data

Optical and thermal remote sensing observations available from the Moderate Resolution Imaging Spectroradiometer (MODIS) (Didan, 2015) onboard the Aqua platform were used in the present study (Table 2) for estimating G and its associated SEB fluxes. These include 8 d land surface temperature (LST or T_S) at 13:30 and 01:30, surface emissivity (ε_s) (MYD11A2), daily surface albedo (α_R) (MCD43A3), and 16 d NDVI (MYD13A2). The overpass times of MODIS Aqua are at 13:30 and 01:30. The 8 d average values of clear-sky T_S available from MYD11A2 data were used (source: https://vip.arizona.edu/documents/viplab/MYD11A2.pdf, last access: 3 December 2022) for all nine flux tower sites. Since the MYD21A2 LST product is known to provide better accuracy (1–1.5 K) (Hulley et al., 2016) as compared with the MYD11A2 LST over semi-arid and arid ecosystems, the former was also used in STIC-TI to compare LE and H estimates (see Table 5 in Sect. 4.4) with the estimates of the MYD11A2 LST over the arid and semi-arid sites (IND-Jai, IND-Naw, US-Ton). The noon–night pair of thermal remote sensing observations from Aqua are close to the times of occurrences of maximum and minimum soil surface temperature (see Fig. 2) and are therefore ideal for soil heat flux modeling using thermal inertia. The MODIS Terra overpass times are at 11:00 and 23:00 and are far from the times of occurrences of minimum–maximum soil temperatures. Therefore, MODIS Aqua acquisition times were used in the present study.

Table 2 A summary of MODIS Aqua optical and thermal remote sensing products used in the present study.

3 Methodology

3.1 Coupled soil heat flux–SEB model

In this paper, we modified a thermal inertia (TI)-based soil heat flux (G) model using noon–night thermal remote sensing observations and thereafter coupled the TI-based G with STIC1.2. A clear-sky net radiation (R_N) model was also introduced into this coupled model, and the R_N estimation algorithm is described in Appendix B. The estimation of G by modifying the MV2007-TI approach and its coupling with STIC1.2 is the most novel component of the modeling scheme, and it is therefore described in the main body of the paper (Sect. 3.1.1). Such a coupling enabled the implementation of a mechanistic G model along with an analytical SEB model using optical–thermal remote sensing data. The coupled model is hereafter referred to as STIC-TI.

3.1.1 MV2007 soil heat flux model based on TI

The functional form for estimating instantaneous G (G_i hereafter) (Eq. 2 below) is based on the harmonic analysis of soil surface temperature and is described in detail by Murray and Verhoef (2007) and Maltese et al. (2013).

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{G}_{\mathrm{i}}& =\mathrm{\Gamma }\left[\left(\mathrm{1}-\mathrm{0.5}{f}_{\mathrm{C}}\right)\right(\sum _{n=\mathrm{1}}^{k}A\sqrt{n\mathit{\omega }}\sin (n\mathit{\omega }t\\ & +{\mathit{\varphi }}_n^{\prime }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathit{\pi }\mathrm{\Delta }t}{\mathrm{12}}}\left)\right)]=\mathrm{\Gamma }{J}_{\mathrm{S}}\end{array}\end{array}$$

G_i is the soil heat flux at the surface at a particular instance (W m⁻²); Γ is the soil thermal inertia (J m⁻² K⁻¹ s^−0.5), k is the total number of harmonics used, A is the amplitude (^∘C) of the nth soil surface temperature (T_ST) harmonic, ω is the angular frequency (rad s⁻¹), t is the time (s), ${\mathit{\varphi }}_n^{\prime }$ is the phase shift of the nth soil surface temperature harmonic (rad), J_S is the summation of harmonic terms of soil surface temperature (K), and Δt (s) is the time offset between the canopy composite temperature and the below-canopy soil surface temperature. Here, we represent G_i and A as ecosystem-scale (≤1 km) soil heat flux and surface soil temperature amplitude (averaged from the soil surface to 10 cm depth), respectively, and assume it to be valid for a different vegetated landscape.

Since we have considered a single pair (noon–night corresponding to 13:30 and 01:30) of MODIS Aqua LST data in the present study, the phase shift (${\mathit{\varphi }}_n^{\prime }$) is taken as 0, and the number of harmonics is taken as 1 (k=1) for estimating G_i. Thus, Eq. (2) is modified as follows,

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{G}_{\mathrm{i}}& =\mathrm{\Gamma }\left[\left(\mathrm{1}-\mathrm{0.5}{f}_{\mathrm{C}}\right)\left(A\sqrt{\mathit{\omega }}\sin \left(\mathit{\omega }t+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathit{\pi }\mathrm{\Delta }t}{\mathrm{12}}}\right)\right)\right]\\ & =\mathrm{\Gamma }{J}_{\mathrm{S}}\end{array}\end{array}$$

with the two boundary values (i.e., Δt=1.5 h for f_c=1 and Δt=0 for f_c=0, f_c being the vegetation fraction), and a linear approach is proposed here to describe the time offset Δt as a function of f_c (Maltese et al., 2013). For a given day, f_c was derived by normalizing NDVI with the upper–lower limits of the annual NDVI cycle.

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }t=\mathrm{1.5}{f}_{\mathrm{c}}\end{array}$$

Scaling function for estimating ecosystem-scale surface soil temperature amplitude (A)

Estimating ecosystem-scale A involves two steps, (a) computing point-scale soil surface temperature amplitude (from the surface to 0.1 m depth) (T_STA hereafter) from the available measurements of soil surface temperature and (b) linking T_STA with remote sensing variables to develop scaling functions for A. Point-scale soil temperature measured at different depths within the top 10 cm soil layer at AU-ASM, US-Ton, and US-Var was averaged and considered to be a representative surface soil temperature (0–10 cm). For Ind-Dha and AU-Ade, single-depth (10 cm) soil temperature measurement was used. Studies also showed that LST carries some signal beneath the skin of the surface (Johnston et al., 2022).

Several studies suggested a theoretical sinusoidal trajectory of soil surface and sub-surface temperatures (Gao et al., 2010), where the amplitude is maximum at the surface, and it gradually decreases with depth to become 37 % of the surface amplitude until the damping depth (Hillel, 1982). However, at deeper depths, soil temperatures remain constant with time and do not show many fluctuations as compared with surface or near-surface soil temperatures. This invariant soil temperature is called deep soil temperature (Mihailovic et al., 1999). However, the diurnal surface soil temperature measurements (within the top 0.1 m depth) across different flux tower sites showed a sinusoidal–exponential behavior, i.e., a sinusoidal pattern from sunrise until the afternoon and an exponential pattern from afternoon through sunset to the next sunrise. An illustrative example of the theoretical and observed trajectories of surface soil temperature is shown in Fig. 2. This diurnal surface soil temperature variation has a single harmonic component (Gao et al., 2010). To compute T_STA, a theoretical half-curve of the sinusoidal pattern is assumed and was derived from measurements as exemplified in Fig. 2.

Figure 2 An illustrative example of the typical diurnal variation of observed soil temperature (T_ST) (from the surface to 0.1 m depth) at OzFlux sites and timings of MODIS AQUA observations. Here, T_STmax and T_STmin are maximum and minimum point-scale observed soil surface temperatures.

It is evident from Fig. 2 that T_STmin represents the minimum surface soil temperature occurring 1–1.5 h after sunrise, and T_STmax occurs during 12:30–15:00 local time. The in situ measured T_ST on completely clear-sky days at OzFlux sites was used to extract T_STmax and T_STmin, and T_STA was derived as T_STmax–T_STmin from the theoretical half-curve of the sinusoidal pattern.

T_STA is generally influenced by several land surface characteristics such as surface temperature and surface albedo of the soil–canopy complex, surface heat capacities, fractional canopy cover, and thermal conductivity. T_S and α_R are the major thermal and reflective land surface properties that have strong synergy with surface soil temperature dynamics. Hence, we have used bivariate regression analysis to develop a scaling function for estimating ecosystem-scale T_STA (top to 0.1 m depth). The bivariate regression is based on the difference of noon (d) and night (n) T_S data and α_R (Duan et al., 2013; Tian et al., 2014) from MODIS Aqua. The scaling function given in Eq. (5) estimates ecosystem-scale T_STA (symbolized as “A” in Eq. 5) from the surface to 0.1 m soil depth:

$$\begin{array}{}\text{(5)}& A=B_{\mathrm{1}}\left(T_{\mathrm{Sd}}-T_{\mathrm{Sn}}\right)+B_{\mathrm{2}}\left({\mathit{\alpha }}_{\mathrm{R}}\right)+B_{\mathrm{3}}\end{array}$$

Here, B₁, B₂, and B₃ are coefficients of the regression model. T_Sd and T_Sn are noontime and nighttime LST, respectively. The results of this regression analysis are elaborated on in Sect. 4.1.

Estimating Γ

Γ is the key variable for estimating G_i using Eq. (2). MV2007 adopted the concept of normalized thermal conductivity (Johansen, 1975) and developed a physical method to estimate Γ as follows:

$$\begin{array}{}\text{(6)}& \mathrm{\Gamma }=e^{\left[{\mathrm{{\rm Y}}}^{\prime }\left(\mathrm{1}-S_{\mathrm{r}}^{({\mathrm{{\rm Y}}}^{\prime }-\mathit{\delta })}\right)\right]}\left({\mathit{\tau }}_{\ast }-{\mathit{\tau }}_{\mathrm{0}}\right)+{\mathit{\tau }}_{\mathrm{0}}\end{array}$$

where τ_∗ and τ₀ are the thermal inertia for saturated and air-dry soil (J m⁻² K⁻¹ s^−0.5), ${\mathit{\tau }}_{\mathrm{0}}=D_{\mathrm{1}}{\mathit{\theta }}_{\ast }+D_{\mathrm{2}}$; ${\mathit{\tau }}_{\ast }=D_{\mathrm{3}}\left({\mathit{\theta }}_{\ast }^{-\mathrm{1.29}}\right)$, Υ^′ (unitless) is a parameter depending on the soil texture (Murray and Verhoef, 2007; Minasny and Hartemink, 2011; Anderson et al., 2007), S_r (m³ m⁻³) is relative saturation and is equal to ($\mathit{\theta }/{\mathit{\theta }}_{\ast }$), and δ (unitless) is the shape parameter which is dependent on the soil texture. θ_∗ (m³ m⁻³) is the soil porosity (equal to the saturated soil moisture content when the soil moisture suction is zero), θ (cm³ cm⁻³) is the volumetric soil moisture, and D₁, D₂, and D₃ are coefficients which were derived from a large number of experimental data. The reported global values of D₁, D₂, and D₃ were taken as −1062.4, 1010.8, and 788.2, respectively (Maltese et al., 2013). The values for θ_∗ and the shape parameter for soil textures across the study sites were specified according to Van Genuchten (1980). The details are mentioned in Table E1 of Appendix E.

In the present study, the relative soil moisture saturation, $S_{\mathrm{r}}(\mathit{\theta }/{\mathit{\theta }}_{\ast })$, is represented in terms of an aggregated moisture availability (M) of the canopy–soil complex through a linear function (Eq. 12). In the case of zero canopy cover, M represents the soil moisture availability from the surface to 0.1 m depth. In the sparse and open canopy, rates of moisture availability from soil to root and root to canopy were assumed to be the same.

Theoretically, M is expressed as the available soil moisture fraction between field capacity (θ_fc) and the permanent wilting (θ_wp) point as given in Eq. (7) below.

$$\begin{array}{}\text{(7)}& M={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{wp}}}{{\mathit{\theta }}_{\mathrm{fc}}-{\mathit{\theta }}_{\mathrm{wp}}}}\end{array}$$

where θ_fc (m³ m⁻³) is the volumetric soil moisture at the field capacity (at a suction of 330 hPa), and θ_wp (m³ m⁻³) is the volumetric soil moisture at the permanent wilting point (at a suction of 15 000 hPa). Since θ_fc, θ_∗, and θ_wp are characteristic volumetric soil moisture contents corresponding to specific suctions and depend on the soil texture, dividing the numerator and denominator in Eq. (7) by θ_∗ gives the following expression:

$$\begin{array}{}\text{(8)}& M={\displaystyle \frac{\frac{\mathit{\theta }}{{\mathit{\theta }}_{\ast }}-\frac{{\mathit{\theta }}_{\mathrm{wp}}}{{\mathit{\theta }}_{\ast }}}{\frac{{\mathit{\theta }}_{\mathrm{fc}}}{{\mathit{\theta }}_{\ast }}-\frac{{\mathit{\theta }}_{\mathrm{wp}}}{{\mathit{\theta }}_{\ast }}}}\end{array}$$

Due to their dependence on soil texture, the ratios (${\mathit{\theta }}_{\mathrm{fc}}/{\mathit{\theta }}_{\ast }$) and (${\mathit{\theta }}_{\mathrm{wp}}/{\mathit{\theta }}_{\ast }$) are treated as constants. These are represented as C and C^′ in the later equations (Eqs. 9, 10, and 11). The constants C and C^′ vary from 0.3 to 0.8 and from 0.1 to 0.4 (Murray and Verhoef, 2007; Minasny and Hartemink, 2011; Anderson et al., 2007), respectively, over different soil textures.

$$\begin{array}{}\text{(9)}& {\displaystyle }M={\displaystyle \frac{\frac{\mathit{\theta }}{{\mathit{\theta }}_{\ast }}-C^{\prime }}{C-C^{\prime }}}\text{(10)}& {\displaystyle }M\left(C-C^{\prime }\right)=\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\ast }}}\right)-C^{\prime }\end{array}$$

By replacing S_r in Eq. (6) as $\mathit{\theta }/{\mathit{\theta }}_{\ast }$ and by rearranging Eq. (10), the following linear function is obtained.

$$\begin{array}{}\text{(11)}& S_{\mathrm{r}}={\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\ast }}}=M\left(C-C^{\prime }\right)+C^{\prime }=M^{\prime }\end{array}$$

Thus, the modified equation to calculate Γ is given by Eq. (12) as follows:

$$\begin{array}{}\text{(12)}& \mathrm{\Gamma }=e^{\left[{\mathrm{{\rm Y}}}^{\prime }\left(\mathrm{1}-M^{{\prime }^{({\mathrm{{\rm Y}}}^{\prime }-\mathit{\delta })}}\right)\right]}\left({\mathit{\tau }}_{\ast }-{\mathit{\tau }}_{\mathrm{0}}\right)+{\mathit{\tau }}_{\mathrm{0}}\end{array}$$

By substituting the values obtained from Eqs. (4), (5), and (12) into Eq. (3), we obtained the instantaneous ecosystem-scale G_i corresponding to the MODIS Aqua noontime overpass. The intrinsic link between G_i estimates through MV2007-TI and the SEB scheme in STIC1.2 is made through M, where the computation of M follows the procedure as described in Mallick et al. (2016, 2018a, b) and Bhattarai et al. (2018) (description in Appendix C).

Estimating M

In STIC1.2, an aggregated moisture availability (M) of the canopy–soil complex is expressed as the ratio of the “vapor pressure difference” between the aerodynamic roughness height of the canopy (i.e., source/sink height) and air to the “vapor pressure deficit” between aerodynamic roughness height and the atmosphere:

$$\begin{array}{}\text{(13)}& M={\displaystyle \frac{\left(e_{\mathrm{0}}-e_{\mathrm{A}}\right)}{\left(e_{\mathrm{0}}^{\ast }-e_A\right)}}={\displaystyle \frac{\left(e_{\mathrm{0}}-e_{\mathrm{A}}\right)}{\mathit{\kappa }\left(e_{\mathrm{S}}^{\ast }-e_{\mathrm{A}}\right)}}={\displaystyle \frac{s_{\mathrm{1}}\left(T_{\mathrm{0}\mathrm{D}}-T_{\mathrm{D}}\right)}{\mathit{\kappa }{s}_{\mathrm{2}}\left(T_{\mathrm{S}}-T_{\mathrm{D}}\right)}}\end{array}$$

where e₀ and $e_{\mathrm{0}}^{\ast }$ are the actual and saturation vapor pressure at the source/sink height, e_A is the atmospheric vapor pressure, $e_{\mathrm{S}}^{\ast }$ is the saturation vapor pressure at the surface, T_0D is the dew-point temperature at the source/sink height, T_S is the LST, T_D is the air dew-point temperature, s₁ and s₂ are the psychrometric slopes of the saturation vapor pressure and temperature between the (T_0D−T_D) versus (e₀−e_A) and (T_S−T_D) versus ($e_{\mathrm{S}}^{\ast }-e_{\mathrm{A}}$) relationship, and κ is the ratio between ($e_{\mathrm{0}}^{\ast }-e_{\mathrm{A}}$) and ($e_{\mathrm{S}}^{\ast }-e_{\mathrm{A}}$). To solve Eq. (13), estimation of T_0D is necessary. An initial estimate of $T_{\mathrm{0}\mathrm{D}}(T_{\mathrm{0}\mathrm{D}}=[(e_{\mathrm{S}}^{\ast }-e_{\mathrm{A}})-s_{\mathrm{3}}{T}_{\mathrm{S}}+s_{\mathrm{1}}{T}_{\mathrm{D}}]/(s_{\mathrm{1}}-s_{\mathrm{3}}\left)\right)$ and M was obtained following Venturini et al. (2008), where s₁ and s₃ were approximated in T_D and T_S, respectively. However, Eq. (13) cannot be solved directly because there are two unknowns in one equation. However, since T_0D also depends on LE (Mallick et al., 2016, 2018a), an iterative update of T_0D (and M) was carried out by expressing T_0D as a function of $\mathit{LE}(T_{\mathrm{0}\mathrm{D}}=T_{\mathrm{D}}+(\mathit{\gamma }\mathit{LE}/\mathit{\rho }{c}_p{g}_A{s}_{\mathrm{1}}\left)\right)$, which is described in detail by Mallick et al. (2016, 2018a) and Bhattarai et al. (2018). In the numerical iteration, s₁ was not updated to avoid numerical instability, and it was expressed as a function of T_D.

3.1.2 STIC-TI: coupling the modified MV2007-TI and STIC 1.2

The initiation of the coupling between MV2007-TI and STIC1.2 was executed by linking G_i estimates from the modified MV2007-TI with M estimates from STIC1.2. Having the initial estimates of M (through Eq. 13), an initial estimation of G_i was made from Eq. (2), where S_r in Eq. (11) was replaced with the initial estimates of M^′. From the initial estimates of G_i (Eq. 2) and R_Ni (equations in Appendix B), initial estimates of LE_i and H_i were obtained through the Penman–Monteith energy balance (PMEB) equation. Analytical expressions of the conductances for estimating H and LE through the PMEB equation were obtained by solving the state equations as described in the Appendix. The process was then iterated by updating $T_{\mathrm{0}\mathrm{D}}[T_{\mathrm{0}\mathrm{D}}=T_{\mathrm{D}}+(\mathit{\gamma }\mathit{LE}/\mathit{\rho }{c}_{\mathrm{p}}{g}_{\mathrm{A}}{s}_{\mathrm{1}}\left)\right]$ and M in every time step (as mentioned in Mallick et al., 2016, 2018a) and re-estimating G_i (using Eq. 3), net available energy (R_Ni−G_i), conductances, LE_i and H_i until stable estimates of LE_i were obtained. The conceptual block diagram and algorithm flow of STIC-TI are shown in Fig. 3a and b, respectively.

Figure 3 (a) Conceptual diagram of the STIC-TI model showing different input variables and model outputs and (b) algorithmic flow for estimating G and associated SEB fluxes through STIC-TI.

Examples of the iterative stabilization of G_i and LE_i for the Indian, Australian, and US ecosystems of India are shown in Fig. 4. The iterative stabilization of G_i and LE_i was obtained between 8 and 25 iterations for all the sites.

Figure 4 Illustrative examples of iterative stabilization of STIC-TI G_i (yellow marker line) and LE_i (grey marker line) in (a) IND-Jai, (b) AU-ASM, and (c) US-Ton.

The noteworthy features of STIC-TI are (1) estimating G by modifying the mechanistic MV2007-TI model using noon and midnight T_S information from thermal remote sensing observations available through a polar-orbiting satellite platform (e.g., MODIS Aqua), (2) coupling the mechanistic MV2007-TI G model with STIC1.2 to simultaneously estimate surface moisture availability (M), G, and SEB fluxes, (3) introducing water stress information in G (through M) to better constrain the aerodynamic and canopy–surface conductances as well as the SEB fluxes, and (4) derivation of the amplitude of ecosystem-scale surface soil temperature (from the topsoil to 0.1 m soil depth).

3.1.3 Generation of remote sensing inputs

Two of the key variables in SEB modeling are T_S and ε_s. These two variables were retrieved at 923 m spatial resolution from MODIS Aqua noon–night TIR observations (MYD11A2) in bands 11.03 and 12.02 µm using a generalized split-window algorithm (Wan and Li, 1997). For optimal retrieval, tractable sub-ranges of atmospheric column water vapor and lower-boundary air surface temperature were used. Land surface emissivity was estimated from land cover types and anisotropy factors. The MYD21A2 LST product was generated using the temperature–emissivity separation (TES) algorithm (Hulley et al., 2016) and an improved water vapor scaling method to remove the atmospheric effects. Albedo was estimated from the MODIS (MCD43A2 Version 6.0) Bidirectional Reflectance Distribution Function and Albedo (BRDF/Albedo) daily dataset (Schaaf et al., 2002) at 462 m spatial resolution. Actual albedo is a value which is interpolated between white-sky and black-sky albedo as a function of fractional diffuse skylight (which is a function of aerosol optical depth). NDVI was obtained from the level-3 MODIS vegetation index product (MYD13Q1, version 6.1), which is generated every 16 d at 250 m (m) spatial resolution. All the input remote sensing variables mentioned in Table 2 were resampled to the spatial resolution of the MYD11A2 product (923 m).

3.2 Sensitivity and statistical analysis

The accuracy of STIC-TI heavily depends on the accuracy of T_S, NDVI, and α_R due to the dual role of T_S in estimating M and G_i, the role of NDVI in G_i, and the combined role of T_S and α_R in estimating R_Ni. Therefore, one-dimensional sensitivity analysis was conducted to assess the impacts of uncertainty in T_S, NDVI and α_R on G_i, H_i and LE_i. The sensitivity was assessed by varying noontime T_S by ±0.5, ±1.0, and ±1.5 K (keeping nighttime T_S constant so that amplitude can vary automatically), varying NDVI by ±0.05, ±0.10, and ±0.15, and varying albedo by ±0.02, ±0.05, and ±0.10, respectively. SEB fluxes were computed by using T_S, NDVI, and α_R for three different periods of the year in all eight ecosystems. Sensitivity analyses were conducted by increasing and decreasing systematically T_S, NDVI, and α_R from its central value while keeping the other variables and parameters constant. This procedure was selected because the fluxes and intermediate outputs of the STIC-TI model reflect an integrated effect due to uncertainty in T_S. In the first run, SEB fluxes were computed using in situ T_S measurements obtained from the flux tower outgoing longwave radiation measurements. Then T_S was increased and decreased at a constant interval and a new set of fluxes was estimated. In the similar way, α_R and NDVI were increased and decreased at constant intervals, and a new set of fluxes was computed. The sensitivity of STIC-TI was assessed by Eq. (14).

$$\begin{array}{}\text{(14)}& \mathrm{Sensitivity}={\displaystyle \frac{E_{\mathrm{i}\mathrm{0}}-E_{\mathrm{iM}}}{O_{\mathrm{i}}}}\cdot \mathrm{100}\end{array}$$

E_i0 is the estimated (original) model output, and E_iM is the estimated (modified) output obtained by changing the variable whose sensitivity is to be tested. O_i is the actual measurements. Apart from the sensitivity analysis, the following set of statistical metrics was used to assess model performances.

$$\begin{array}{}\text{(15)}& {\displaystyle }{R}^{\mathrm{2}}={\left({\displaystyle \frac{\sum _{i=\mathrm{1}}^n(E_{\mathrm{i}}-\stackrel{\mathrm{‾}}{E})(O_{\mathrm{i}}-\stackrel{\mathrm{‾}}{O})}{\sqrt{\sum _{i=\mathrm{1}}^n(E_{\mathrm{i}}-\stackrel{\mathrm{‾}}{E}{)}^{\mathrm{2}}}\sqrt{\sum _{i=\mathrm{1}}^n(O_{\mathrm{i}}-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}}}\right)}^{\mathrm{2}}\text{(16)}& {\displaystyle }\mathrm{RMSE}=\sqrt{\sum _{i=\mathrm{1}}^n{{\displaystyle \frac{(E_{\mathrm{i}}-O_{\mathrm{i}})}{n}}}^{\mathrm{2}}}\text{(17)}& {\displaystyle }\mathrm{BIAS}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n(E_{\mathrm{i}}-O_{\mathrm{i}})}{n}}\text{(18)}& {\displaystyle }\mathrm{MAPD}={\displaystyle \frac{\mathrm{100}}{n}}\sum _{i=\mathrm{1}}^n\left|{\displaystyle \frac{E_{\mathrm{i}}-O_{\mathrm{i}}}{O_{\mathrm{i}}}}\right|\text{(19)}& {\displaystyle }\mathrm{KGE}=\mathrm{1}-\sqrt{{\left(r-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{E}}}{{\mathit{\sigma }}_{\mathrm{o}}}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{E}}{\stackrel{\mathrm{‾}}{O}}}-\mathrm{1}\right)}^{\mathrm{2}}}\end{array}$$

R² is the coefficient of determination, RMSE is the root-mean-square error, BIAS is the mean bias, MAPD is the mean absolute percent deviation, KGE is the Kling–Gupta efficiency, n is the total number of data pairs, and the bar indicates the mean value of the measured variable and the model estimates of the same variable. E_i and O_i are the model-estimated and model-measured SEB fluxes, r is the Pearson correlation coefficient, $\stackrel{\mathrm{‾}}{O}$ is the average of the measured values, $\stackrel{\mathrm{‾}}{E}$ is the average of the estimated values, σ_o is the standard deviation of observation values, and σ_E is the standard deviation of the estimated values. The KGE has been widely used for calibration and evaluation of hydrological models in recent years, and it combines the three components of the Nash–Sutcliffe efficiency (NSE) of model errors (i.e., correlation, bias, ratio of variances or coefficients of variation) in a more balanced way, but it has not been widely used for analyzing the ET model performances. KGE = 1 indicates perfect agreement between modeled estimates and observations. The performance of a model is considered “poor” for KGE between 0 and 0.5, and models with negative KGE values are considered “not satisfactory”.

4 Results

4.1 Ecosystem-scale surface soil temperature amplitude (A)

The scaling functions developed to estimate ecosystem-scale (1 km) surface soil temperature amplitude (A) from point-scale T_STA were used to estimate G_i. However, before the development of the scaling functions, analysis was carried out to investigate the relationship of soil temperature amplitude between the two different spatial scales. The scatterplot (Fig. 5a) of the noon–night LST difference (ΔTs) versus T_STA for different albedo classes showed a linear increase in ΔTs with increasing T_STA. However, some divergence of data points within the cluster was also noticed, which could be associated with different albedo (α_R) levels. A bivariate linear function was fitted between T_STA as a predictand (Y) versus ΔTs (T_sd−T_sn) and α_R as predictors (X₁ and X₂, respectively). The function was found to be $Y=\mathrm{0.59}{X}_{\mathrm{1}}-\mathrm{51.3}{X}_{\mathrm{2}}+\mathrm{8.66}$ by combining the data of nine ecosystems (r=0.86). The coefficients in the above expressions correspond to B₁ (0.59), B₂ (51.3), and B₃ (8.66) of Eq. (5). The estimated amplitude from these ecosystem-scale predictors and scaling functions was treated as an ecosystem-scale surface soil temperature amplitude (A).

Figure 5 (a) Two-dimensional scatterplots between (ΔTs) versus T_STA at different α_R levels over different ecosystems. Here T_STA on the y axis is the observed soil temperature amplitude that is used to develop the scaling function, and ΔTs is the noon–night LST difference of MODIS AQUA. (b) Validation of ecosystem-scale estimates of A from the above functions over different sites.

Validation of ecosystem-scale estimates of A from the above functions over different sites is shown in Fig. 5b with respect to T_STA for the independent datasets. The estimated A was found to have an MAPD of 19.9 %, negative bias, and R²=0.90 over different ecosystems. The temporal variation of estimated A and T_STA is shown in Fig. D1 in Appendix D. Further analysis was carried out to investigate the bias in A at three fractional vegetation cover (f_c) classes (f_c<0.3; $\mathrm{0.3}\le f{f}_{\mathrm{c}}\le \mathrm{0.5}$; f_c>0.5) representing bare soil (class 1), 30 %–50 % canopy cover (class 2), and more than 50 % canopy cover (class 3), respectively. While negative bias was noted for class 1 and class 3 (−0.54 and −0.83 ^∘C), the bias was positive (0.49 ^∘C) in the intermediate f_c, which represents sparse and patchy canopy cover. The signals of surface albedo, emissivity, and temperature of soil surface and canopy are relatively pure in class 1 and class 3 as compared with class 2, where the surface signal carries more heterogeneity. Given that T_STA is computed from the in situ measurements, it is likely to carry more heterogeneity in class 2 as compared with the other two classes. The land surface emissivity in MYD11A2 was estimated from land cover types and the anisotropy factor, which have differential impacts on T_ST and T_S, leading to such an opposite bias in class 2 as compared with class 1 and class 3.

4.2 Sensitivity analysis of STIC-TI G_i, LE_i and H_i to land surface variables

4.2.1 Sensitivity of G_i to land surface variables

The average sensitivity of G_i to three land surface variables (T_S, NDVI, α_R) by combining the estimates of wet and dry periods is shown in Fig. 6. G_i was found to be substantially sensitive to T_S, with error magnitude ranging from 2 % to 18 % due to T_S uncertainties of ±0.5–2.5 K (Fig. 6a), with greater sensitivity to T_S during the summer season. The median sensitivity of G_i due to ±5 %–10 % uncertainty in α_R varied from 5 % to 12 % in all the ecosystems (Fig. 6b). The uncertainties in NDVI revealed 2 % to 15 % error in G_i estimates (Fig. 6c), and no significant difference in the mean sensitivity due to NDVI uncertainties was noted between the ecosystems. The sensitivity of G_i decreased with increasing values of NDVI.

4.2.2 Sensitivity of LE_i and H_i to land surface variables

Both LE_i and H_i were sensitive to T_S on the order of 2 %–29 % (LE_i) and 5 %–35 % (H_i) for a T_S uncertainty of ±0.5–2.5 K from its mean values (Table 3). Interestingly, LE_i was more sensitive to T_S uncertainties as compared with H_i in the rainfed ecosystems. The highest mean sensitivity of LE_i to T_S was found in arid (IND-Jai: 2 %–28 %), semi-arid (AU-ASM: 5 %–21 %), tropical savanna (IND-Dha: 3 %–26 %), savanna (US-Ton: 4 %–29 %), and arid (US-Var: 3 %–26 %) ecosystems. The mean sensitivity of H_i to T_S was maximum in sub-humid (IND-Sam: 2 %–32 %), semi-arid (IND-Naw: 2 %–28 %), and savanna (AU-Ade: 8 %–17 %) (Table 3). A greater sensitivity of the SEB fluxes due to α_R uncertainties was found than due to NDVI. The median sensitivity of LE_i and H_i due to 10 % uncertainty from mean α_R varied within 2 %–16 % in all the ecosystems (Table 3). By contrast, errors in the two SEB fluxes were substantially low (2 %–13 %) due to ±0.05–0.15 uncertainty from mean NDVI (Table 3).

Figure 6 Sensitivity of STIC-TI G_i due to uncertainties in T_S (a), α_R (b), and NDVI (c) for nine flux tower sites in India, the United States, and Australia. The uncertainties were introduced by taking the mean values of these variables during three different periods (summer, rainy, and winter) of a year. Mean uncertainties of the three periods are presented in the figure.

Table 3 Sensitivity (in percent) of LE_i and H_i due to T_S, NDVI, and α_R uncertainties.

4.3 Comparative evaluation of the STIC-TI and contemporary G_i models

The performances of the STIC-TI and existing G_i models were evaluated and compared with respect to in situ G_i measurements. The existing models reported by Moran et al. (1989), Bastiaanssen et al. (1998), Su (2002), and Boegh et al. (2004) have been considered for comparing with the TI-based model. These four existing models are referred to here as MOR89, BAS98, SU02, and BO04, respectively. While the models MOR89, SU02, and BO04 are based on linear regression between G versus NDVI, BAS98 is based on multivariate regression of G with NDVI, LST, and α_R. The performance of the STIC-TI was substantially better as compared with MOR89, SU02, and BO04 with respect to MAPD (19 %), RMSE (22 W m⁻²), and the coefficient of determination (R²=0.8) when compared with in situ measurements over one Indian, three Australian and two US flux tower sites (Table 4) and comparable with the BAS98 G_i model. The validation plot of retrieved noontime G_i from STIC-TI is shown in Fig. 7.

Figure 7 Validation of noontime (13:30) G_i estimates with respect to in situ measurements in different ecosystems. The regression between the two sources of G_i is G_i (STIC-TI) = 0.90 G_i (tower) −0.10.

Table 4 A comparison of error statistics of G_i estimates from STIC-TI and existing G_i models over different ecosystems.

The RMSE varied from 9 to 20 W m⁻², with MAPD ranging from 12 % to 21 % across individual flux tower sites. A high magnitude of G_i was predicted in the arid and semi-arid systems (120–240 W m⁻²) as compared with the humid systems (20–90 W m⁻²), which was in close correspondence to the observations. The model also captured the range of G_i that is generally found in different biomes (20–140 W m⁻² for grasslands, 20–90 W m⁻² for cropland) (Purdy et al., 2016). Due to the paucity of G_i measurements, direct validation of G_i was only possible for 32 d (concurrent to MODIS overpass) at the IND-Dha site. Overall, STIC-TI tends to provide reasonable G estimates for the terrestrial ecosystems, with a soil temperature amplitude above 5 ^∘C.

4.4 Evaluation of STIC-TI LE_i, H_i, and EF

The modeled versus measured LE_i and H_i showed good agreement in all nine ecosystems with RMSE in LE_i and H_i estimates using the MYD11 LST product on the order of 29–62 W m⁻² and 26–61 W m⁻², MAPD of 9 %–31 % and 20 %–36 %, BIAS of −29 to 38 W m⁻² and −44 to 32 W m⁻² (Fig. 8a, b; Table 5), and a high R² of 0.8.

Figure 8 (a) Validation of STIC-TI LE_i estimates with respect to in situ measurements in different ecosystems. (b) Validation of STIC-TI H_i estimates with respect to in situ measurements in different ecosystems.

Table 5 Error statistics of STIC-TI LE_i and H_i estimates with respect to EC measurements in different ecosystems of India, the USA, and Australia using the MYD11A2 LST product for all nine sites and using the MYD21A2 LST product for three semi-arid and arid sites. The statistics obtained by using MYD21A2 LST are shown in the parentheses.

Arid ecosystems in India (IND-Jai) and the US (Ton and Var) and the semi-arid ecosystem in Australia (AU-ASM) revealed relatively high MAPD (31 %, 25 %, 27 %, and 28 %) (Table 5). In general, STIC-TI was able to produce the dominant convective heat fluxes with respect to the EC measurements, as evident through the low RMSE for H_i and high RMSE for LE_i in IND-Jai, US-Ton, US-Var, and AU-Ade, where LE_i is inherently low except for a few rainy days. A uniform distribution of data points around the 1 : 1 validation line (Fig. 8a) indicated an overall low BIAS in LE_i estimates. However, modeled H_i was consistently lower than the observations (negative BIAS) in the tropical savanna (IND-Dha and AU-How) and semi-arid (IND-Naw) ecosystems (−44)–(−25) and −26 W m⁻²), while a consistent positive BIAS was observed in AU-ASM (semi-arid), AU-Ade (savanna), and US-Var (arid) (Fig. 8b; Table 5). This consequently led to an overall low negative BIAS (−10 W m⁻²) and a relatively low R² in H_i (R²=0.8) as compared with the errors in LE_i (BIAS = 15 W m⁻², R²=0.9). The regression between the modeled and tower measurements of LE_i is LE_i (STIC-TI) = 0.98LE_i (tower) − 0.266. The regression between the modeled and tower measurements of H_i is H_i (STIC-TI) = 0.93H_i (tower) + 4.90. The KGE statistics varied in the range of 0.71–0.95 for LE_i and in the range of 0.64–0.91 for H_i, respectively, across all nine flux tower sites, thus revealing a reasonably high efficiency of the model in capturing the magnitude and variability of SEB fluxes.

The impact of the MODIS Aqua day–night view angle difference (δVZA) on STIC-TI fluxes was investigated. Estimated errors in terms of the mean percent deviation in LE_i, H_i, and G_i with respect to measurements for each 10^∘ bin over 16 angular bins within ±80^∘ were analyzed in response to the mean δVZA of each angular bin. G_i errors (X) were found to be significantly correlated with δVZA (Y) in a parabolic ($Y=\mathrm{0.0027}{X}^{\mathrm{2}}-\mathrm{0.0025}X+\mathrm{1.4919}$; r=0.73) pattern (refer to Appendix F, Fig. F1a). Errors in G_i on the order of −5 % to 10 %, 10 %–15 %, and > 15 % were largely found to be within ±30, ±45, and >45 to −80^∘ δVZA, respectively. The errors in H_i were found to have a strong linear ($Y=-\mathrm{0.1452}X+\mathrm{1.1146}$, r=0.77) dependence on δVZA (refer to Appendix F, Fig. F1b). However, a weak dependence of LE_i errors ($Y=-\mathrm{0.0878}X+\mathrm{2.0314}$, r=0.5) on δVZA (refer to Appendix F, Fig. F1c) was found, as the majority of the errors were within ±10 %, which corresponded to ±60^∘ δVZA. The nature of the relations and the degree of dependency of the model flux errors on δVZA in this study would be helpful for minimizing the error budget in surface energy balance fluxes from future thermal infrared missions with day–night observations.

Figure 9 Relationship between the KGE of STIC-TI (G_i and LE_i) with $\mathrm{\Delta }{T}_{\mathrm{S}}/f_{\mathrm{c}}$ in different terrestrial ecosystems.

Further investigation was done on whether the KGE for STIC-TI G_i and LE_i follows any systematic pattern, and the ratios ΔT_S and f_c were used as a proxy for surface heterogeneity and dryness. The plot of the KGE of G_i and LE_i with this ratio is shown in Fig. 9. KGE–G_i was found to show a systematic decrease with an increase in the ΔT_S–f_c ratio of up to 40, after which it remained unchanged with an increase in the ratio. Although the KGE of LE_i also decreased (20 % reduction) with an increase in the ΔTs–f_c ratio, KGE-LE_i was found to increase beyond ΔTs−f_c 40. This revealed that the model efficiency remained high (>0.8) within certain dryness limits (ΔTs−f_c ratio <20 and >50) and that the efficiency reduced moderately (within 0.7–0.8) for intermediate dryness. Interestingly, the use of MYD21A2 LST in STIC-TI showed improvements (see the parentheses in the different columns in Table 5) in LE_i and H_i error statistics as compared with using MYD11A2 LST in terms of higher R² and KGE and lower RMSE in LE_i (3 %–8 % less) and H_i (2 %–3 % less) for semi-arid and arid sites such as IND-Jai, IND-Naw, and US-Ton.

An independent evaluation of multi-temporal heat fluxes over two US flux sites for the years 2016–2018 is shown in Figs. 10 and 11. STIC-TI G_i estimates with the MYD11A2 LST product showed a close match with in situ measurements with respect to intra- and inter-annual variability in G_i, followed by LE_i and H_i. This further demonstrates the merit of the coupled model for reproducing ecosystem-scale G_i estimates, especially for shorter and open canopies.

Figure 10 Illustrative examples of the temporal evolution of STIC-TI-derived fluxes using the MYD11A2 LST product versus observed SEB fluxes for 3 consecutive years from 2016 to 2018 in a grassland ecosystem in the United States (e.g., US-Var).

Figure 11 Illustrative examples of the temporal evolution of STIC-TI-derived fluxes using the MYD11A2 LST product versus observed SEB fluxes for 3 consecutive years from 2016 to 2018 in a woody savanna ecosystem in the United States (e.g., US-Ton).

The temporal behavior of STIC-TI and the observed evaporative fraction (EF) (ratio of LE and R_N–G) (Fig. 12) along with observed monthly rainfall (P) distinctly captured the substantial temporal variability in EF during the dry-to-wet transition in the Indian study sites, which also corresponded to low (high) θ and P. In IND-Naw and IND-Sam, a marked rise (>0.4) in STIC-TI EF was noted during days of the year (DOYs) 25 to 75, when wheat is grown under assured irrigation. The impact of irrigation is thus captured by the substantial increase in EF in the absence of P. In contrast, the rainfed grassland system (IND-Jai) showed a peak EF (∼ 0.8), which corresponded to southwestern monsoon rainfall during June to September and a progressive decline in EF during the dry-down period in October to April corresponding to the post southwestern monsoon phase. Some intermittent spikes in EF were also noted during the dry-down phase in both STIC-TI and observations. The intermittent EF spikes during the soil moisture dry-down phase could be due to enhanced LE through moisture advection from the surrounding vegetation causing a greater enhancement of evaporation than expected. This is known as the “clothesline effect”, which frequently occurs in semi-arid and arid ecosystems. In addition to IND-Jai, the response of both modeled and measured EF to wet and dry spells was also noted during the southwestern monsoon period at all other flux tower sites of India.

Figure 12 Illustrative examples of the temporal variation of STIC-TI-derived EF using the MYD11A2 LST product with respect to measured EF and P in (a) IND-Naw, (b) IND-Jai, (c) IND-Sam, and (d) IND-Dha.

The temporal behavior of EF from the STIC-TI using the MYD11A2 LST product and EC measurements along with measured θ and P at the OzFlux and AmeriFlux sites also revealed (Fig. 13) a close correspondence of STIC-TI to EC observations. Low EF (0.05–0.40) during the dry season around DOY 100–250 and high EF (>0.4) during the wet season (DOY 1–120 and 300–360) in AU-ASM, US-Ton, and US-Var was observed. The analysis showed that STIC-TI EF can capture the annual variability of observed EF and its responses across different ecosystems during the wet and dry seasons. The plots of STIC-TI EF versus measured θ (in the inset of Fig. 13) revealed triangular scatter close to the right-angled triangle with a positive slope of the hypotenuse in three ecosystems: AU-ASM, US-Var, and US-Ton. This showed, in the water-controlled ecosystems, that distinct wet–dry seasons exist and that the positive EF–θ relationship is an outcome of the soil moisture controls on transpiration during the dry season.

Figure 13 Comparison of the temporal variation of STIC-TI-derived EF using the MYD11A2 LST with respect to measured EF, θ, and P in (a) AU-ASM, (b) US-Var, and (c) US-Ton. The scatterplots in the inset show the relationship between STIC-TI EF with respect to measured θ.

5 Discussion

5.1 Interaction of flux and internal SEB metrices

From Sect. 4.1 we found a relatively reduced sensitivity of G_i to T_S uncertainties. In any given condition, if an over(under)estimation of M due to noontime T_S uncertainties (through Eq. 13) leads to an over(under)estimation of Γ, the effects of such over(under)estimation of Γ (due to noontime T_S uncertainties) tend to be compensated by under(over)estimation of amplitude A (in Eq. 5), ultimately leading to a reduction in the sensitivity of G_i to T_S. While the scatter between G versus A for a wide range of Γ (Fig. 14a) revealed large scatter with increasing amplitude under the dry conditions (low Γ), the scatter between Γ versus T_S for different M (Fig. 14b) revealed an exponential reduction in Γ with increasing T_S and dryness and almost no significant change in Γ with increasing T_S at a constantly high dryness (M<0.25). Thus, the confounding effects of Γ, A, and M through Eqs. (3), (5), (12), and (13) led to a reduction in sensitivity of G to T_S, as exemplified in Fig. 14.

Figure 14 Response plots among parameters of the TI-based G_i model, such as (a) G_i versus amplitude (A) for varying Γ and (b) noontime T_S versus Γ with varying M.

Concerning LE_i and H_i, dual uncertainties could be propagated in both the fluxes through daytime T_S (through M and G_i), leading to high sensitivity of these two SEB fluxes due to T_S perturbations. The relatively high sensitivity of LE_i to T_S (as compared with H_i) in the non-irrigated ecosystems could be due to partial compensation of $g_{\mathrm{A}}/g_{\mathrm{S}}$ in both the numerator and denominator of the PMEB equation for H (Eq. C7 of Appendix C). A recent study (Fig. 10 in Mallick et al., 2018a) showed high sensitivity of g_S due to T_S (a 1 % change in T_S led to a 5.2 %–7.5 % change in g_S) as compared with g_A sensitivity to T_S (a 1 % change in T_S led to a 1.6 %–2 % change in g_A), suggesting that errors in g_S due to T_S uncertainty tend to be larger than errors in g_A. Partial cancellation of the conductance errors in the numerator of Eq. (C7) might have resulted in compensation of H_i errors in the water-limited ecosystems. In this environment, the variability of LE_i is mainly dominated by $g_{\mathrm{A}}/g_{\mathrm{S}}$, which makes LE_i highly sensitive due to T_S uncertainties. Combined uncertainty due to $g_{\mathrm{A}}/g_{\mathrm{S}}$ in the denominator and g_A in the numerator of Eq. (C6) resulted in greater sensitivity in LE_i to T_S in the arid and tropical savanna ecosystems (Mallick et al., 2015, 2018a; Winter and Eltahir, 2010). The very low sensitivity of LE_i and H_i due to uncertainties in NDVI is because NDVI was not used in the conductance parameterizations and effects due to NDVI in STIC-TI only being propagated through G_i. The sensitivity of LE_i and H_i to albedo was mainly due to the dependence of net radiation (R_Ni) on albedo, and any resultant uncertainty in R_Ni (due to albedo) tends to be reflected in the sensitivity of LE_i and H_i to albedo.

5.2 Possible sources of errors in SEB flux evaluation

In STIC-TI, underestimation and overestimation errors in G_i in different ecosystems (Fig. 7) could originate due to the errors in the MYD11A2 LST product. A host of studies previously reported that the standard deviations of errors in retrieved emissivity in bands 31 and 32 are 0.009, and the maximum error in retrieved T_S of the MOD11A1 LST falls within 2–3 K, which is mainly due to the errors in surface emissivity correction (Duan et al., 2017; Wan, 2014). In the present analysis, we found an overestimation error of MODIS T_S in the range of 0.5–1.5 K when compared with in situ infrared temperature measurements at the tropical savanna site. As mentioned in Sect. 3.1, a positive (negative) bias in T_S would tend to an overestimation (underestimation) of amplitude (A) in Eq. (5), underestimation (overestimation) of M in Eq. (13), and consequent underestimation (overestimation) of Γ (Eq. 12) and G_i, respectively. Furthermore, the standard deviation of NDVI surrounding the tower sites varied from 0.01 to 0.05 when compared with the ground measurements, which could be another source of error in the STIC-TI model. In addition, NDVI saturates at LAI > 3. However, STIC-TI provides direct estimates of ecosystem G and is independent of R_N.

Despite the comparable accuracy of current G estimates with the G model of Bastiaanssen et al. (1998), the foundation of STIC-TI lies in the use of soil moisture characteristics with varying soil textural types which are known to influence the soil heat conductance and thereby G. Thus, the control of soil moisture on evaporation is explicitly included in STIC-TI as opposed to the semi-empirical G function of Bastiaanssen et al. (1998). The higher accuracies of the TI-based thermal diffusion model as compared with R_N-based empirical G models were also reported by Purdy et al. (2016) at daily or longer timescales in cropland and grassland. All these G model estimates may at times differ from in situ measurements due to not accounting for leaf litter presence or the layer on the soil floor in the remote-sensing-based G model.

The overestimation (underestimation) of LE_i (H_i) is also due to the effects of the spatial resolutions of different input variables on these two SEB fluxes and conducted statistical evaluation with respect to the measured SEB fluxes. Eswar et al. (2017) demonstrated the need for spatial disaggregation models for monitoring LE_i at field scale using contextual models by disaggregation of the evaporative fraction (Λ) and downwelling shortwave radiation ratio (R_G). Using different disaggregation models, they estimated LE_i at 250 m spatial resolution and reported an RMSE of 30–32 W m⁻² as compared with LE_i obtained at a 1000 m spatial resolution with an RMSE of 40–70 W m⁻² over different sites in India. Anderson et al. (2007) reviewed different validation experiments conducted in diverse agricultural landscapes (Anderson et al., 2004, 2005; Norman et al., 2003) and reported RMSE in LE_i in the range of 35–40 W m⁻² (15 %) at 30–120 m disaggregated spatial resolution. Current analysis also brought out the need for noon–night thermal imaging with a spatial resolution finer than 1000 m to adequately capture the magnitude and variability of LE_i in the terrestrial ecosystems, especially agroecosystems where average field sizes are lower (<0.5 ha) and fragmented, such as in India and other sub-continents.

As seen in Fig. 8a and Table 5, there is a gross overestimation of LE_i with respect to the tower observations when MYD11A2 LST is used. The consistent positive BIAS in STIC-TI LE_i in five out of nine sites is presumably due to the overestimation of R_Ni (Fig. B1 of Appendix B) and the underestimation of G_i. Figure 7 shows an overestimation of G_i for three OzFlux sites and US sites and an underestimation of G_i for the Indian site with G_i (STIC-TI) = 0.90 G_i (tower) − 0.10 and an overestimation of R_Ni at the ecosystem scale, with R_Ni (STIC-TI) = 0.78R_Ni (tower) + 58.92 (Appendix B2). This means that a systematic overestimation of net available energy (R_Ni−G_i) will be obvious in cases where STIC-TI shows an underestimation of G_i. Since available energy is an important component for estimating LE through the PMEB equation, an overestimation of net available energy leads to an overestimation of LE by STIC-TI. Sensible heat flux will be consequently underestimated due to the complementary nature of the PMEB equation. It may also be noted that the use of MYD21A2 LST led to a relatively better accuracy in LE_i (3 %–8 %) and H_i (2 %–3 %) as compared with using MYD11A2 LST in semi-arid and arid ecosystems. The higher retrieval accuracy of MYD21A2 LST using the TES algorithm over MYD11A2 LST that uses a split-window algorithm (Wan and Li, 1997) is the main reason for obtaining higher accuracy in LE_i and H_i estimates.

The standard deviations of the MODIS Aqua day–night overpass time over the study sites were found to be within 30–45 min (Sharifnezhadazizi et al., 2019), and the expected deviation in LST from the mean local time would be around ±0.75 K (Sharifnezhadazizi et al., 2019). Sensitivity analysis showed that resultant uncertainties in STIC-TI heat flux estimates would be on the order of ±5 %–7 % due to this LST uncertainty.

5.3 Effects of SEB closure

Given that there is a widespread lack of SEB closure ($H+\mathit{LE}\ne R_N-G$) or residual energy balance, knowledge of the impact of different vegetation types and climatic variables on SEB “non-closure” is essential. A recent study by Dare-Idowu et al. (2021) covering eight growing seasons and three crops (maize, wheat, and rapeseed) in two sites of southwestern France showed that the systematic effect of each site on SEB closure was stronger than the influence of crop type and stage. The same study revealed a greater percentage of SEB closure under unstable atmospheric conditions and in the prevailing wind directions, and sensible heat advection accounted for more than half of the imbalance at both the sites.

In our study, using unclosed SEB observations for Indian sites in the absence of in situ G_i observations also added to the consistent positive BIAS in the statistical evaluation of LE_i. A ubiquitous lack of energy balance closure on the order of 10 %–20 % worldwide at most of the EC sites is reported in different literatures (Stoy et al., 2013; Wilson et al., 2002), which implies a systematic underestimation (overestimation) of LE_i (EC tower) (and/or the H_i EC tower). Accommodating an average 15 % imbalance in LE_i (EC tower) would tend to diminish the positive BIAS in STIC-TI. Therefore, the pooled gain (0.98) and positive BIAS between the STIC-TI and tower LE_i are determined by the overestimation of (R_Ni−G_i) combined with the underestimation of measured LE_i from the EC towers. An underestimation of H_i (negative BIAS) is associated with two reasons: (a) ignoring the two-side aerodynamic conductance of the leaves (Jarvis and McNaughton, 1986; Monteith and Unsworth, 2013; Schymanski et al., 2017), which could lead to substantial underestimation of H_i, and (b) the complementary nature of the PMEB equation; if LE_i is overestimated, H_i will be underestimated. In addition, frequent micro-advection fluxes alter measured in situ H and LE fluxes, but these advection conditions are not explicitly accounted for in the current STIC-TI model. At the EC tower sites, the fraction of the residual energy balance to R_N can be quantified with respect to vegetation/crop growth characteristics or biophysical properties. However, where G observations are lacking, such as in many Indian EC tower sites, the TI-based G model can be used to fill up the missing G observations to quantify residual energy balance and to correct the SEB non-closure.

6 Summary and conclusions

This study addressed one of the long-term uncertainties in thermal remote sensing of evaporation modeling in open-canopy heterogeneous ecosystems, which is associated with empirical methods of estimating ground heat flux. We demonstrated for the first time physical integration and coupling of a mechanistic ground heat flux model with an analytical evaporation model (Surface Temperature Initiated Closure, STIC) within the surface energy balance equation. The model is called STIC-TI, which uses satellite-based land surface temperature from MODIS Aqua and associated biophysical variables, and it has minimal independence of in situ measurements. The estimation of evaporation through STIC-TI does not require any empirical function for inferring the biophysical conductances. STIC-TI is independent of uncertain parameterizations of surface roughness and atmospheric stability and does not involve any look-up table for biome or plant functional attributes. By linking noon–night land surface temperature with a harmonics equation of thermal inertia and soil moisture availability, STIC-TI derived the ground heat flux and subsequently coupled it with evaporation. Independent validation of STIC-TI with respect to eddy covariance flux measurements from nine terrestrial ecosystems in arid, semi-arid, and sub-humid climates in India, the USA, and Australia led us to the following conclusions.

• i.

  While the MODIS Aqua day–night view angle difference showed a strong impact on ground heat flux and sensible heat flux estimated deviations of STIC-TI (with respect to measurements), a relatively weak dependence of latent heat flux errors on the day–night view angle difference was noted.

• ii.

  The most notable advantages of STIC-TI are firstly that it provides direct estimates of ground heat flux while simultaneously integrating the effects of soil water stress on ground heat flux and evaporation through the inclusion of noon–night land surface temperature information. Secondly, the ecosystem-scale surface soil temperature amplitude used in the ground heat flux model can advance our understanding of associated terrestrial ecosystem processes.

The requirement of few input variables in STIC-TI generates promise for surface–atmosphere exchange studies using readily available data from the current-generation remote sensing satellites (e.g., MODIS, VIIRS) that have noon–night land surface temperature observations. STIC-TI can also be potentially used for distributed ET mapping from future high spatial resolution (∼ 50–60 m) TIR missions with noon–night observations with high revisits, such as the Indo-French mission, TRISHNA (Thermal infrared Imaging Satellite for High-resolution Natural Resource Assessment) (Lagouarde et al., 2018, 2019), ESA's LSTM (Land Surface Temperature Monitoring), and NASA SBG (Surface Biology and Geology), respectively.

Appendix A

Table A1 A list of symbols, their descriptions, and units used in the present study.

Table A2 Summary of instruments used, height or depth and period of measurements, and measured variables at nine EC flux tower sites.

Appendix B

B1 Clear-sky instantaneous net radiation (R_Ni) model

Net radiation (R_N) is defined as the difference between the incoming and outgoing radiation, which includes both longwave and shortwave radiation at the Earth's surface.

Terrestrial R_N has four components: downwelling and upwelling shortwave radiation (R_G and R_R) and downwelling and upwelling longwave radiation (R_L↓ and R_L↑), respectively.

$$\begin{array}{}\text{(B1)}& R_{\mathrm{N}}=\left(R_{\mathrm{G}}-R_{\mathrm{R}}\right)+\left(R_{\mathrm{L}\downarrow }-R_{\mathrm{L}\uparrow }\right)\end{array}$$

Out of these four terms mentioned in Eq. (B1), R_G and R_L↓ are dependent on various factors such as geographic location, season, cloudiness, aerosol loading, and atmospheric water vapor content and less on surface properties. On the other hand, the upwelling radiations in Eq. (B1) strongly depend on the surface properties such as surface reflectance and emittance, land surface temperature, and soil water content (Zerefos and Bais, 2013).

Instantaneous net radiation (R_Ni) can be derived using Eq. (B2) as follows (Mallick et al., 2007):

$$\begin{array}{}\text{(B2)}& {\displaystyle }{R}_{\mathrm{Ni}}=R_{\mathrm{ns}}-R_{\mathrm{nl}}\text{(B3)}& {\displaystyle }{R}_{\mathrm{ns}}=(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{R}})R_{\mathrm{G}}\text{(B4)}& {\displaystyle }{R}_{\mathrm{nl}}=R_{\mathrm{L}\downarrow }-R_{\mathrm{L}\uparrow }\end{array}$$

where R_ns is net shortwave radiation (W m⁻²), R_nl is net longwave radiation (W m⁻²), and α_R is the broadband surface albedo shortwave spectrum.

A WMO (World Meteorological Organization) shortwave radiation model (Cano et al., 1986) calibrated over Indian conditions (Mallick et al., 2007, 2009) was used to compute R_G using the following equation:

$$\begin{array}{}\text{(B5)}& R_{\mathrm{G}}={\mathit{\tau }}_{\mathrm{sw}}{G}_{\mathrm{sc}}E(\sin {\mathit{\beta }}_{\mathrm{e}}{)}^{\mathrm{1.15}}\end{array}$$

where τ_sw is the global clear-sky transmissivity for the shortwave radiation (0.7), G_sc is the solar constant (1367 W m⁻²), ε is the eccentricity correction factor due to variation in the Sun–Earth distance, and β_e is the Sun elevation in degrees.

R_L↓ at any instance was calculated as follows:

$$\begin{array}{}\text{(B6)}& R_{\mathrm{L}\downarrow }={\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }(\mathrm{273.14}+T_{\mathrm{A}}{)}^{\mathrm{4}}\end{array}$$

where σ is the Stefan–Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴), T_A is the air temperature (^∘C), and ε_a is the atmospheric emissivity.

Atmospheric emissivity (ε_a) was computed using the following equation (Bastiaanssen et al., 1998):

$$\begin{array}{}\text{(B7)}& {\mathit{\epsilon }}_{\mathrm{a}}=\mathrm{0.85}-\ln {\mathit{\tau }}_{\mathrm{sw}}^{\mathrm{0.09}}\end{array}$$

R_L↑ at any particular instance was calculated as follows:

$$\begin{array}{}\text{(B8)}& R_{\mathrm{L}\uparrow }={\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }(\mathrm{273.14}+T_{\mathrm{s}}{)}^{\mathrm{4}}\end{array}$$

where ε_s is the surface emissivity in the thermal infrared (8–14 µm) spectrum and T_S is the land surface temperature (^∘C).

B2 Evaluation of STIC-TI R_Ni

Comparison of the clear-sky R_Ni estimates with respect to in situ measurements revealed RMSE in R_Ni on the order of 27–72 W m⁻², MAPD 8 %–24 %, BIAS (−67)–50 W m⁻², and R² varying from 0.62 to 0.90 across all the sites (Fig. B2, Table B2). Among the nine sites, a consistent underestimation of R_Ni was noted in IND-Dha, US-Ton, and US-Var (with BIAS of −23, −61, and −67 W m⁻²), whereas substantial overestimation of R_Ni was found in IND-Sam, IND-Naw, and AU-ASM with BIAS of 50, 37, and 43 W m⁻², respectively (Table B2).

Figure B1 Validation of STIC-TI-derived R_Ni estimates with respect to in situ measurements in different ecosystems. The regression equation between modeled versus in situ R_Ni is R_Ni (STIC-TI) = 0.78R_Ni (tower) + 58.92.

Table B1 Performance evaluation statistics of clear-sky R_Ni estimates in nine different ecosystems.

Appendix C

C1 Estimating SEB fluxes using the STIC1.2 analytical model and thermal remote sensing data

STIC1.2 (Mallick et al., 2014, 2015a, b, 2016, 2018a) is a one-dimensional physically based SEB model and is based on the integration of satellite LST observations into the Penman–Monteith energy balance (PMEB) equation (Monteith, 1965). In STIC1.2, the vegetation–substrate complex is considered as a single unit. Therefore, the aerodynamic conductances from individual air–canopy and canopy–substrate components is regarded as an “effective” aerodynamic conductance (g_A), and surface conductances from individual canopy (stomatal) and substrate complexes are regarded as an “effective” canopy–surface conductance (g_S) which simultaneously regulates the exchanges of sensible and latent heat fluxes (H and LE) between the surface and atmosphere. One of the fundamental assumptions in STIC1.2 is the first-order dependence of these two critical conductances on M through T_S. Such an assumption enabled an integration of satellite LST into the PMEB model (Mallick et al., 2016). The common expression for LE and H according to the PMEB equation is as follows.

$$\begin{array}{}\text{(C1)}& {\displaystyle }\mathrm{LE}={\displaystyle \frac{s\mathit{\varphi }+\mathit{\rho }{c}_{\mathrm{p}}{g}_{\mathrm{A}}{D}_{\mathrm{A}}}{s+\mathit{\gamma }\left(\mathrm{1}+\frac{g_{\mathrm{A}}}{g_{\mathrm{S}}}\right)}}\text{(C2)}& {\displaystyle }H={\displaystyle \frac{\mathit{\gamma }\mathit{\varphi }\left(\mathrm{1}+\frac{g_{\mathrm{A}}}{g_{\mathrm{S}}}\right)-\mathit{\rho }{c}_{\mathrm{p}}{g}_{\mathrm{A}}{D}_{\mathrm{A}}}{s+\mathit{\gamma }\left(\mathrm{1}+\frac{g_{\mathrm{A}}}{g_{\mathrm{S}}}\right)}}\end{array}$$

In the above equations, the two biophysical conductances (g_A and g_S) are unknown, and the STIC1.2 methodology is based on finding analytical solutions for the two unknown conductances to directly estimate LE (Mallick et al., 2016, 2018a). The need for such analytical estimation of these conductances is motivated by the fact that g_A and g_S can neither be measured at the canopy nor at larger spatial scales, and there is no universally agreed appropriate model of g_A and g_S that currently exists (Matheny et al., 2014; van Dijk et al., 2015). By integrating T_S with standard SEB theory and vegetation biophysical principles, STIC1.2 formulates multiple state equations in order to eliminate the need to use the empirical parameterizations of the g_A and g_S and also to bypass the scaling uncertainties of the leaf-scale conductance functions to represent the canopy-scale attributes. The state equations for the conductances are expressed as a function of those variables that are mostly available as remote sensing observations and weather forecasting models. In the state equations, a direct connection to T_S is established by estimating M as a function of T_S. The information of M is subsequently used in the state equations of conductances, aerodynamic variables (aerodynamic temperature, aerodynamic vapor pressure), and evaporative fraction, which is eventually propagated into their analytical solutions. M is a unitless quantity, which describes the relative wetness (or dryness) of a surface and also controls the transition from potential to actual evaporation, which implies M→1 under saturated surface conditions and M→0 under extremely dry conditions. Therefore, M is critical for providing a constraint against which the conductances are estimated. Since T_S is extremely sensitive to the surface moisture variations, it is extensively used for estimating M in a physical retrieval scheme (detail in Appendix A3 of Bhattarai et al., 2018; Mallick et al., 2016, 2018a). It is hypothesized that linking M with the conductances will simultaneously integrate the information of T_S into the PMEB model. To illustrate, we express the state equations by symbols, ${\mathrm{sv}}_{\mathrm{1}}=f\mathit{\{}{c}_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},v_{\mathrm{1}},v_{\mathrm{2}},v_{\mathrm{3}},v_{\mathrm{4}},{\mathrm{sv}}_{\mathrm{3}},{\mathrm{sv}}_{\mathrm{5}}\mathit{\}}$; ${\mathrm{sv}}_{\mathrm{2}}=f\mathit{\{}{v}_{\mathrm{4}},{\mathrm{sv}}_{\mathrm{1}},{\mathrm{sv}}_{\mathrm{5}},{\mathrm{sv}}_{\mathrm{6}}\mathit{\}}$; ${\mathrm{sv}}_{\mathrm{3}}=f\mathit{\{}{c}_{\mathrm{3}},v_{\mathrm{3}},v_{\mathrm{4}},{\mathrm{sv}}_{\mathrm{4}},{\mathrm{sv}}_{\mathrm{5}}\mathit{\}}$; ${\mathrm{sv}}_{\mathrm{4}}=f\mathit{\{}{c}_{\mathrm{3}},v_{\mathrm{3}},{\mathrm{sv}}_{\mathrm{1}},{\mathrm{sv}}_{\mathrm{2}},{\mathrm{sv}}_{\mathrm{7}},{\mathrm{sv}}_{\mathrm{8}}\mathit{\}}$. Here, f, sv, v, and c denote the function, state variables, input variables (five input variables; radiative and meteorological), and constants (three constants), respectively. Here sv₁ to sv₄ are g_A, g_S, aerodynamic temperature (T₀), evaporative fraction (Λ), and sv₈ is M. Given the estimates of M, net radiative energy (R_Ni−G_i), T_A, R_H, the four state equations are solved simultaneously to derive analytical solutions for the four state variables and to produce a surface energy balance “closure” that is independent of empirical parameterizations for g_A, g_S, T₀, and Λ. The state equations of STIC are given below.

$$\begin{array}{}\text{(C3)}& {\displaystyle }{g}_{\mathrm{A}}={\displaystyle \frac{\mathit{\varphi }}{\mathit{\rho }{c}_{\mathrm{p}}\left[(T_{\mathrm{0}}-T_{\mathrm{A}})+\left(\frac{e_{\mathrm{0}}-e_{\mathrm{A}}}{\mathit{\gamma }}\right)\right]}}\text{(C4)}& {\displaystyle }{g}_{\mathrm{S}}=g_{\mathrm{A}}{\displaystyle \frac{\left(e_{\mathrm{0}}-e_{\mathrm{A}}\right)}{\left(e_{\mathrm{0}}^{\ast }-e_{\mathrm{0}}\right)}}\text{(C5)}& {\displaystyle }{T}_{\mathrm{0}}=T_{\mathrm{A}}+\left({\displaystyle \frac{e_{\mathrm{0}}-e_{\mathrm{A}}}{\mathit{\gamma }}}\right)\left({\displaystyle \frac{\mathrm{1}-\mathrm{\Lambda }}{\mathrm{\Lambda }}}\right)\text{(C6)}& {\displaystyle }\mathrm{\Lambda }={\displaystyle \frac{\mathrm{2}\mathit{\alpha }s}{\mathrm{2}s+\mathrm{2}\mathit{\gamma }+\mathit{\gamma }\frac{g_{\mathrm{A}}}{g_{\mathrm{S}}}\left(\mathrm{1}+M\right)}}\end{array}$$

Detailed derivations of these four state equations are given in Mallick et al. (2016, 2018a). However, the analytical solutions to the four state equations contain three accompanying unknowns; e₀ (vapor pressure at the source/sink height), $e_{\mathrm{0}}^{\ast }$ (saturation vapor pressure at the source/sink height), and Priestley–Taylor coefficient (α), and as a result there are four equations with seven unknowns. Consequently, an iterative solution was needed to determine the three unknown variables (as described in Appendix A2 in Mallick et al., 2016). Once the analytical solutions of g_A and g_S are obtained, both variables are returned into Eq. (13) to directly estimate LE.

In STIC-TI, an initial value of α was assigned as 1.26; initial estimates of $e_{\mathrm{0}}^{\ast }$ were obtained from T_S through temperature–saturation vapor pressure relationship, and initial estimates of e₀ were obtained from M as, $e_{\mathrm{0}}=e_{\mathrm{A}}+M(e_{\mathrm{0}}^{\ast }-e_{\mathrm{A}})$. Initial T_0D and M were estimated according to Venturini et al. (2008) as described in Sect. 3.2, and initial estimation of G was performed from initial M using the equation sets Eqs. (2)–(11). With the initial estimates of these variables; first estimate of the conductances, T₀, Λ, H, and LE were obtained. The process was then iterated by updating $e_{\mathrm{0}}^{\ast }$, D₀, e₀, T_0D, M, and α (using Eqs. A9, A10, A11, A17, A16, and A15 in Mallick et al., 2016), with the first estimates of g_S, g_A, T₀, and LE, and re-computing G, φ, g_S, g_A, T₀, Λ, H, and LE in the subsequent iterations with the previous estimates of $e_{\mathrm{0}}^{\ast }$, e₀, T_0D, M, and α until the convergence of LE was achieved. Stable values of G, conductances, LE, H, T₀, $e_{\mathrm{0}}^{\ast }$, e₀, T_0D, M, and α were obtained within ∼ 25 iterations. The inputs needed for computation of LE_i (Eq. C6) are air temperature (T_A), land surface temperature (T_S), relative humidity (R_H), net radiation (R_Ni) and soil heat flux (G_i).

Appendix D

The temporal variation of estimated A and T_STA is shown in Fig. D1. The annual variations of T_STA in different ecosystem was found to be within the ranges of 1–4 ^∘C.

Figure D1 Temporal variation of A and T_STA in (a) AU-ASM (2013), (b) US-Ton (2014), and (c) US-Var (2014).

Appendix E

Table E1 Soil textural properties and their values used in the present study (Murray and Verhoef, 2007; Minasny and Hartemink, 2011; Anderson et al., 2007).

Appendix F

The day–night view angle effect on errors of STIC-TI heat flux estimates from measurements is shown in Fig. F1.

Figure F1 Dependence of STIC-TI model flux error in terms of mean percent deviation from measurements of day–night view angle differences of MODIS Aqua expressed as a mean of the 10^∘ bin for (a) G_i, (b) H_i, and (c) LE_i.