Accuracy of five ground heat flux empirical simulation methods in the surface-energy-balance-based remote-sensing evapotranspiration models
Based on the assessment from 230 flux site observations, intra-day and daytime ground heat flux (G) accounted for 19.2 % and 28.8 % of the corresponding net radiation, respectively. This indicates that G plays an important role in remote-sensing (RS) energy-balance-based evapotranspiration (ET) models. The G empirical estimation methods have been evaluated at many individual sites, while there have been relatively few multi-site evaluation studies. The accuracy of the five empirical G simulation methods in the surface-energy-balance-based RS–ET models was evaluated using half-hourly observations. The linear coefficient (LC) method and the two methods embedded with the normalized difference vegetation index (NDVI) were able to accurately simulate a half-hourly G series at most sites. The mean and median Nash–Sutcliffe efficiency (NSE) values of all sites were generally higher than 0.50 in each half-hour period. The accuracy of each method varied significantly at different sites and at half-hour intervals. The highest accuracy was exhibited during 06:00–07:00 LST (all times hereafter are LST), followed by the period of 17:00–18:00. There were 92 % (211/230) sites with an NSE of the LC method greater than 0.50 at 06:30. It showed a slightly higher accuracy during nighttime periods than during daytime periods. The lowest accuracy was observed during the period of 10:00–15:30. The sites with an NSE exceeding 0.50 only accounted for 51 % (118/230) and 43 % (100/230) at 10:30 and 13:30, respectively. The accuracy of the model was generally higher in Northern Hemisphere sites than in Southern Hemisphere sites. In general, the highest and lowest accuracies were observed at the high- and low-latitude sites, respectively. The performance of the LC method and the methods embedded with NDVI were generally satisfactory at the Eurasian and North American sites, with the NSE values of most sites exceeding 0.70. Conversely, it exhibited relatively poor performance at the African, South American, and Oceanian sites, especially the African sites. Both the temporal and spatial distributions of the accuracy of the G simulation were positively correlated with the correlation between G and the net radiation. Although the G simulation methods accurately simulated the G series at most sites and time periods, their performance was poor at some sites and time periods. The application of RS ET datasets covering these sites requires caution. Further improvement of G simulations at these sites and time periods is recommended for the RS ET modelers. In addition, variable parameters are recommended in empirical methods of G simulation to improve accuracy. Instead of the Rn, finding another variable that has a physical connection and strong correlation with G might be a more efficient solution for the improvement, since the weak correlation between G and Rn is the main reason for the poor performance at these regions.
Accurate simulations of evapotranspiration (ET) represent the core of hydrological processes, crop growth, and ecosystem water efficiency simulations (Ponce-Campos et al., 2013). Remote sensing (RS) is the only viable technique that can provide relatively frequent and spatially contiguous ET measurements on global and regional scales (Zhang et al., 2016; Ait Hssaine et al., 2020). Surface energy balance is the main method used in RS ET modeling (Zhang et al., 2016; Chen and Liu, 2020). Ground heat flux (G) accounts for a significant fraction of the surface energy balance (Pauwels and Daly, 2016), but there is insufficient research on these models compared with sensible heat flux (H) (Mohan et al., 2020; Wu et al., 2020). Over bare soils or sparsely vegetated surfaces, G can reach half of the net radiation (Rn) (Heusinkveld et al., 2004). Even under full vegetation cover, G is significant, especially when turbulent processes are less active (Gentine et al., 2012).
Soil heat flux is the heat flux that occurs in a layer of soil. G, which is the soil heat flux at the surface, is difficult to observe directly due to technical limitations (Wang and Bou-Zeid, 2012; Gao et al., 2017), and direct estimation of G using RS data is not possible (Kalma et al., 2008; Allen et al., 2011; Saadi et al., 2018). Soil heat flux is generally measured using heat flux plates near the surface (within a few millimeters of the surface). This measurement is referred to as G′ in this study. However, the difference between G′ and G could be 50 % because of the soil heat storage within the layer from the surface to the flux plate (Heusinkveld et al., 2004; Yue et al., 2011; Wu et al., 2020). A large error is produced if the soil heat storage is ignored in the G calculation (Meyers and Hollinger, 2004; Lu et al., 2018). In addition, the soil above the heat flux plates can lose contact with the rest of the near-surface soil matrix, adversely affecting the water and heat flow (Leuninget al., 2012; Russell et al., 2015). The spatial scale of the G observation is also much smaller than that of the H and latent heat flux (LE) estimates (Shao et al., 2008; Verhoef et al., 2012).
There are numerous schemes for estimating G (Wang and Bou-Zeid, 2012; Gao et al., 2017; Wu et al., 2020); these can be classified into two categories. According to physical mechanisms, G can be calculated from the soil heat conductivity and the vertical gradient of temperature using the so-called gradient approach (Yang and Wang, 2008). A more common approach is to combine G′ at a specific reference depth with the soil heat storage above (Kustas et al., 2000; Wu et al., 2020). G′ can be simulated using the gradient approach or observed by heat flux plates. Soil heat storage in the soil layer above the measured depth can be calculated as the integral over the change in temperature with time multiplied by the volumetric heat capacity of the soil – this is called the calorimetry method (Liebethal and Foken, 2007; Agam et al., 2019). However, applications of these physical mechanism-based approaches are restricted to only a few sites due to the limitations of field observations of soil thermal properties (Mayocchi and Bristowa, 1995; Kustas et al., 2000). Soil thermal properties are affected by soil texture, mineralogical composition, bulk density, and the surrounding environment (e.g., soil moisture and temperature) (Peng et al., 2017; Ju and Hu, 2018). In other words, soil thermal properties vary with time and space. In addition to these physical mechanism-based approaches, G can also be estimated using empirical methods based on Rn, H, or G′ (Cellier et al., 1996; Leuning et al., 2012; Purdy et al., 2016).
To estimate ET in RS models, G is usually obtained from empirical relations with Rn. Choudhury et al. (1987) established an empirical function between G and Rn for bare and vegetated soils. They found that the ratio of G to Rn () was 0.4 for bare soils. For vegetated land cover, the ratio could be described by an exponential relationship with the leaf area index. Kustas and Daughtry (1990) also calibrated the ratio by ground-based measurements for bare soil, alfalfa, and cotton and found the corresponding values to be 0.29 ± 0.05, 0.16 ± 0.035, and 0.27 ± 0.02, respectively. The constant has been used in several RS-based ET models, including the two-source energy balance (TSEB) (Norman et al., 1995), the Atmosphere–Land Exchange Inverse (ALEXI) (Anderson et al., 1997), and the disaggregated ALEXI (DisALEXI) (Norman et al., 2003) models. The ratio values in these models ranged from 0.15 to 0.35. The Global Land Evaporation Amsterdam Model (GLEAM) used 0.05, 0.20, and 0.25 for tall canopy, short vegetation, and bare soil, respectively (Miralles et al., 2011). The is also usually parameterized as an empirical function with the vegetation index in other RS-based ET models, including but not limited to the Surface Energy Balance Algorithm for Land (SEBAL) (Bastiaanssen et al., 1998), Surface Energy Balance System (SEBS) (Su, 2002), Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC) (Allen et al., 2007), and Simplified TSEB (STSEB) (Sánchez et al., 2008) models. The solutions of G in the first two models were also applied to the Simplified Surface Energy Balance Index (S-SEBI) (Roerink et al., 2000) and four-source surface energy balance (SEB-4S) (Merlin et al., 2014) models, respectively. Some studies have modified the parameter values in these empirical relationships, such as in the modified SEBAL (Singh et al., 2008; Faridatul et al., 2020) and TSEB models (Ait Hssaine et al., 2020). In addition, the empirical relations between G and Rn were applied to simulate G in several RS-based global ET datasets. These ET datasets include, but are not limited to, the Breathing Earth System Simulator (BESS) (Jiang and Ryu, 2016), Moderate Resolution Imaging Spectroradiometer (MODIS; MOD16A2) (Mu et al., 2011), GLEAM (Miralles et al., 2011), Numerical Terradynamic Simulation Group (NTSG) (Zhang et al., 2010), and thermal energy balance (Chen et al., 2014; 2021) products. More G empirical estimation methods can be found in Sun et al. (2013) and Bonsoms and Boulet (2022).
Several studies have evaluated the empirical methods in the simulation of G. Liebethal and Foken (2007) evaluated six parameterization approaches for G by using the reference dataset, in which G′ at a depth of −0.2 m and the heat storage in the soil layer above −0.2 m were calculated by the gradient and calorimetry approaches, respectively. They found that the physical mechanism-based calorimetric and simple measurement approaches showed better performance than empirical methods. Similar results were also found in evaluations by Russell et al. (2015) and Gao et al. (2017) for eddy covariance tower measurements in Idaho and Washington, respectively. Saadi et al. (2018) evaluated three empirical methods using flux station measurements and determined that the best results were obtained using the SEBAL method. However, these evaluations were limited to a single-site scale because field observations of soil thermal properties were available only at a few sites. Purdy et al. (2016) evaluated six empirical methods of G simulation against G′ observations at 88 flux sites. This was a very meaningful study on a global scale. However, there is a large difference between G′ and G, which has been described above.
The surface energy balance method provides an alternative solution for assessing the G simulation schemes (van der Tol, 2012). This method could avoid the inconsistent spatial scale of G with that of LE and H in field measurements. The FLUXNET dataset, which contains global flux tower observations, provides a good opportunity to evaluate G simulation methods on a global scale (Knox et al., 2019; Pastorello et al., 2020; Liu, 2021). According to the surface energy balance method, this study uses the FLUXNET dataset to comprehensively evaluate the ability of various schemes to accurately simulate half-hourly G at global flux measurement sites.
This study addresses four key objectives: (1) investigating the temporal and spatial variations and common characteristics of the empirical relationship between G and Rn; (2) evaluating the accuracy of five empirical methods in simulating half-hourly G from Rn; and (3) investigating the performance of five methods at different times during the day and the spatial distribution of simulation accuracy at global flux observation sites. The evaluation results of this study are expected to provide a reference for RS–ET model application and developers.
This study used FLUXNET eddy covariance observations that cover all continents, including FLUXNET2015 (Pastorello et al., 2020) and FLUXNET-CH4 community products (Knox et al., 2019). FLUXNET2015 contains 212 observation sites from 1991 to 2014, while the FLUXNET-CH4 community product contains 81 sites from 2006 to 2018. The longest observational record was 25 years, whereas the shortest was less than one year. Half-hourly data series for LE, H, Rn, and G′ were used. All missing values were eliminated. For example, if there were missing values on a certain day, all data on that day were discarded. Therefore, only days with fully available half-hourly data were used in the analysis. Only sites with a data series longer than 360 d were used. These eliminations ultimately meant that a total of 189 FLUXNET2015 sites and 60 FLUXNET-CH4 sites were used in the analysis because of the lack of observations (Table S1 in the Supplement). There were 19 sites belonging to both FLUXNET2015 and FLUXNET-CH4. G′ was not observed at 63 sites (Table S1). Flux observation data from four sites in Australia were obtained from the TERN OzFlux dataset. These four sites were included in FLUXNET products but with a longer and continuous series up to 2019 (Beringer et al., 2016). In addition, the normalized difference vegetation index (NDVI) data (Vermote et al., 2014), which were derived from surface reflectance data acquired by the advanced very high-resolution radiometer sensor, were used in this study. The dataset had a spatial resolution of 0.05∘ × 0.05∘ and temporal coverage from June 1981 to the present.
2.2.1 The surface energy balance residual method for estimating G
The surface energy balance residual (SEBR) method (Eq. 1) is an effective tool for estimating G using the measurements of other components with a flux tower (van der Tol, 2012). Energy balance is an independent means of assessing G. In the surface energy balance, heat storage in the air between the ground and the height of the eddy covariance system is neglected, as is the horizontal advection of heat and other minor energy sources and sinks (Wilson et al., 2002; Leuning et al., 2012; Stoy et al., 2013). The SEBR method equation is expressed as follows:
where G is ground heat flux, Rn is net radiation, and LE and H are latent and sensible heat flux, respectively. The estimated time series were used as referenced G for evaluating simulations of G in RS–ET models. The mean and standard deviation approach was applied to detect outliers (Liu et al., 2019).
2.2.2 Simulations of G in RS–ET models
Based on the half-hourly series of G simulated by the SEBR method at globally observed sites, the G simulation methods commonly used in five RS ET models were evaluated. The first is the linear coefficient (LC) method. It has been applied in the TSEB, ALEXI, DisALEXI, GLEAM, and other RS ET models to simulate G, but different models use different linear coefficient values. Then the linear models of embedding the NDVI in the form of the power (Bastiaanssen, 1995) and exponential (Choudhury et al., 1987) functions (referenced as LC_NDVI_P and LC_NDVI_E), which were typically applied in the SEBAL (Bastiaanssen et al., 1998) and modified SEBAL (Singh et al., 2008) and SEBS (Chen et al., 2019) ET models, respectively, were evaluated. Finally, the linear models embedded with fractional vegetation coverage, which were applied in the SEBS and STSEB models (referenced as LC_fc_SE and LC_fc_ST), were also evaluated. LC methods embedded with the vegetation index have also been applied in many other RS–ET models, such as the S-SEBI, NTSG, BESS, METRIC, MOD16A2, and SEB-4S models. However, different models may use different parameter values.
In this study, the parameters in these five methods were calibrated rather than the fixed parameter values in the original models. The half-hourly series of the observed Rn and simulated G by the SEBR method at each flux tower site were used for the optimal calibration of the parameters. These five methods are represented as LC, LC_NDVI_P, LC_NDVI_E, LC_fc_SE, and LC_fc_ST. The expressions are as follows:
where α, α1, α2, and b are parameters to be calibrated within the range of [0.01, 1.5], [0.01, 1.5], [0.01, 0.1], and [0.1, 1.5], respectively. At each site, daily series of each half hour were divided into two parts: the first 80 % of the data were used for parameter calibration, and the rest were used for validation. The parameters of these methods were calibrated by the Nash–Sutcliffe efficiency (NSE) at each observation site.
2.2.3 Evaluation criteria
The simulations of G were evaluated using the G series estimated using the SEBR method for each site in each half-hourly period. The criteria tried to evaluate these simulations included the relative error (RE), root mean square error (RMSE), Kling–Gupta efficiency (KGE; Gupta et al., 2009), and NSE (Nash and Sutcliffe, 1970). The RE and RMSE represent the bias deviation from the observed values, whereas the KGE and NSE are indicative of the goodness of fit of the simulated and observed data series. The best-fit value was 1.0, and the goodness of fit deteriorated with increasing deviation from 1.0. The evaluation criteria were calculated as follows:
where Xsi and Xri are the ith values of the simulated and referenced G time series, respectively; n is the time series length; and are the means of the simulated and referenced G, respectively; and SDs and SDr are the standard deviations of the simulated and referenced G, respectively.
3.1 Intra-day distribution of observed surface energy balance items and
The intra-day distribution characteristics of each flux variable were analyzed based on field observation data. Figure 1 shows the intra-day distribution of half-hourly Rn, H, LE, G, and G′. The first four variables and G′ were derived from the mean of 230 and 167 FLUXNET sites (Table S1), respectively. Overall, LE, H, and G accounted for 34.5 %, 46.3 %, and 19.2 % of Rn, respectively. G accounted for 28.8 % of Rn when only daytime periods were considered. This indicates that ignoring G in energy-based models greatly overestimates the ET values. The observed daytime G′ value was only 24.1 % of that of G. Considering the intra-day periods, G′ was only 4.7 % of G. All flux variables were stable and showed little variance from 20:00–06:00 LST (all times hereafter are LST). During this period, the LE was positive and accounted for only 7 % of the total daily LE, whereas other flux variables were negative. It showed a unimodal distribution for all flux variables during the day. The intra-day distribution of H showed the best agreement with the measured Rn (Fig. 1a). However, the intra-day distributions of LE, G, and G′ showed an overall deviation from the measured Rn. The distribution of LE and G′ was generally delayed by half an hour compared to the measured Rn, while that of G was half an hour earlier. The intra-day distribution of each flux variable during the daytime was compared with the sine and Gaussian functions (Fig. 1b–f). The results showed that daytime flux variables were more consistent with the latter than with the sine function, which is commonly used to upscale instantaneous ET to daily values in RS applications. The Gaussian function perfectly matched each flux variable at any time during the day.
The intra-day distribution characteristics of and the ratio of G and H () were also analyzed based on field observation data. The intra-day distributions of and at each site are shown in Fig. 2. During data processing, data points with absolute values greater than 10 in the or daily series of each period were deleted. Outliers in the or series were then removed using the outlier detection method. The varied significantly in different half-hour periods of the intra-day and among the different sites (Fig. 2a). The variation in among the sites was lower during the daytime than that at night. The variation range of among all sites was generally approximately 0.2 in the daytime. This indicates that the of all sites showed high consistency during these periods. At night, the variation mostly ranged between 0.6 and 0.8. In other words, was more consistent across sites during the daytime than at night. The slope and R2 of the linear-fitting curve were −0.012 and 0.92, respectively. The R2 of the polynomial fitting curve reached 0.98 (Fig. 2c).
also varied significantly in different half-hour periods of the intra-day and among the different sites (Fig. 2b). The variation in during each period was greater than that of . In other words, was less consistent across sites than . Like , the variation of among sites during the daytime was significantly lower than that at night. The variation range of among all sites in the daytime was approximately 1.0, while the corresponding range at night was approximately 2.0. In most half-hour periods, the mean and median values of at all sites showed significant differences, with the former generally being approximately 0.5 higher than the latter. During the period of 06:30–16:30, the median and mean values generally showed a unimodal distribution, and the R2 of the polynomial fitting curve for the mean series was 0.95 (Fig. 2d).
3.2 Temporal and spatial analysis of the empirical relationship between G and Rn and between and NDVI
The empirical relationship between G, Rn, and H was analyzed based on the measured data. Overall, G had a strong correlation with Rn but a relatively weak correlation with H. Figure 3 shows the R2 and slope of the linear fitting between G and Rn in each half-hour period of the intra-day. In each period, G and Rn showed a strong linear correlation at most sites, with a fitted R2 generally above 0.4. The mean and median R2 values of all sites were mainly between 0.5 and 0.8. The strong correlation between G and Rn indicates that it is reasonable to use Rn to calculate G in the RS-based energy balance ET models. However, the correlation between G and Rn varied during the different periods. The correlation is relatively high in the periods around 06:00 and 18:00. Especially around 06:00, the R2 of the linear fitting between G and Rn was greater than 0.7 at most sites, and the median R2 of all sites reached 0.8. The correlation between G and Rn in the night periods (20:00–04:30) was slightly stronger than that in the daytime periods (08:00–16:00). During the night periods, the R2 of most sites was generally between 0.45 and 0.70, and the mean and median R2 of all sites were mainly between 0.55 and 0.60. The R2 of most sites generally ranged from 0.40 to 0.65 in the daytime periods, and the mean and median R2 of all sites were concentrated around 0.50. The correlation between G and Rn was relatively low in the period from 10:00 to 15:00, and R2 was lower than 0.65 at most sites, especially in the periods around 14:00, with the mean and median R2 below 0.50.
The slope of the linear fitting of G and Rn in each half-hour intra-day period is shown in Fig. 3b. The fitting slope showed significant differences among the different sites, ranging from 0.1 to 1.1. However, the slopes fitted at most sites also exhibited certain characteristics of a centralized distribution. In each period, the range of the slope at most sites was within 0.3, especially during the daytime periods (08:00–17:30), when the range of the slope was within 0.15. The slopes of most sites were relatively stable during all periods, except for the periods of 06:00–07:00 and 17:00–18:00.
In terms of the seven land cover types, the intra-day performance of each land type was similar to that of all sites except the other type (Fig. 3c and d). The correlation between G and Rn was relatively high during 06:00–07:00 and 17:00–18:00. The correlation in other and wetland types is generally higher than that of other land cover types. In each period, the median R2 of all sites in the two types generally exceeded 0.60, and the highest value even exceeded 0.80. Except for the other type, the difference in the correlation between G and Rn in different land types is mainly reflected in the daytime period. The correlation in the forest and savanna types was significantly lower than that of other types during the daytime, especially for savanna sites, most of which had R2 lower than 0.5 during the daytime. In other-type sites, the correlation between G and Rn in the daytime is stronger than that in the nighttime periods. The slope value of each land cover type in the daytime is lower than that in the night. This intra-day distribution of slope was consistent with that of all sites.
The empirical relationship between G and Rn not only varied significantly at different intra-day periods but also showed great spatial differences among the different sites. The linear fitted R2 between the daily series of G and Rn at each site is shown in Fig. 4. As the median R2 of 48 half-hour periods at each site (Fig. 4a), 91 % of the sites (210/230) showed an R2 > 0.4. The linear-fitting R2 between G and Rn was > 0.6 for 49 % of the sites (114/230). The mean R2 for all sites was 0.58. This indicated that the G in most half-hour periods had a strong correlation with Rn at most observed sites. However, there were also 20 sites where the R2 ranged from 0.2 to 0.4. These sites were mainly located in the middle- and low-latitude regions but were distributed across all observed continents. The correlation was generally stronger in the Northern Hemisphere than in the Southern Hemisphere, with the mean R2 of the northern sites being significantly higher than that of the southern sites. There was a strong correlation between G and Rn at most Eurasian, North American, and Oceanian sites, with a linear fitted R2 generally exceeding 0.4. There was a relatively weak correlation at many African and South American sites, with an R2 value of less than 0.4. At different latitudes, the strongest correlation between G and Rn was found at the middle- and high-latitude (> 45∘) sites with the highest R2 values. The R2 values of these sites exceeded 0.4, with a mean value of 0.65. There was a relatively weak correlation at tropical (< 23.4∘) sites. The R2 values of these sites were relatively low, with a mean value of 0.48.
The spatial distribution of the linear-fitting R2 of G and Rn at 10:30 and 13:30 (Fig. 4b and c) was consistent with that in Fig. 4a, while the R2 value was generally lower than the median of each half-hour period. At 10:30, the mean R2 of all sites was 0.52. A total of 26 % of the sites (60/230) had R2 values lower than 0.4, and 11 sites had R2 values lower than 0.2. Although mainly distributed at low latitudes (< 30∘), these sites were found on all the observed continents. The highest R2 values were at high-latitude (> 60∘) sites, with an average of 0.70. The R2 values of low-latitude sites were significantly lower than those of other sites, with a value of less than 0.4 for most sites and an average of only 0.33. The correlation between G and Rn at 13:30 was weaker than that at 10:30, with R2 values slightly lower than those at 10:30. There was a relatively weak correlation between G and Rn at all African sites (R2 < 0.4). Overall, the results showed a strong correlation between G and Rn at most observed sites. It is reasonable to simulate a daily series of G values from Rn in most areas. However, it is necessary to apply this relationship cautiously in some areas at mid–low latitudes, especially in tropical areas.
For different land cover types, the correlation between G and Rn was the strongest at other and wetland sites. The mean value of the median R2 of statistical sites was 0.71 and 0.67 for these two types, respectively. There was also a strong correlation between G and Rn at cropland, shrubland, and grassland sites. The mean value of corresponding R2 is about 0.60. There was only one site with a median R2 lower than 0.4 in each of the three land cover types. The mean value of corresponding R2 was 0.55 for forest sites. However, the correlation was relatively weak in the savanna-type sites, and the mean value of corresponding R2 is 0.49. The weak correlations between G and Rn (R2 < 0.4) were mainly distributed in the forest (6/97), grassland (7/42), and savanna (3/15) sites.
The empirical correlation between and NDVI was also analyzed. The results showed that the correlation between and NDVI was weak in the daily series but strong in the monthly series. Figure 4d and e show the linear- and exponential-fitted R2 values between the monthly series of and NDVI, respectively, at each observed site. During data processing, only monthly values of observation days greater than 15 d were used, and monthly values affected by frozen soil or snow cover (mean air temperature below −5 ∘C) were excluded. In general, there was a strong correlation between the monthly series and NDVI. The fitted R2 values of the linear and exponential functions were consistent. The median and mean R2 of all sites were 0.71 and 0.65, respectively, and 72 % of the sites (157/218) had an R2 above 0.60. The exponential correlation between and NDVI was stronger than the linear correlation at several sites, including US-Gle and US-Twt. The exponential-fitted R2 of the two sites could be increased from 0.55 and 0.31 of the linear R2 to 0.77 and 0.50, respectively. Conversely, the linear correlation was stronger than the exponential correlation at other sites, such as ES-AMO and CN-QIA. The linear fitted R2 of the two sites could be increased from 0.21 and 0.68 of the exponential R2 to 0.43 and 0.83, respectively. Overall, the spatial distribution of the linear- or exponential-fitted R2 of and NDVI was similar to the linear fitted results of G and Rn. The correlation between and NDVI was stronger in the Northern Hemisphere than in the Southern Hemisphere. The mean R2 value of the northern sites (0.69) was higher than that of the southern sites (0.38). It showed the strongest correlation at the middle- and high-latitude (> 50∘) sites. The R2 values for these sites were generally higher than 0.6, with an average of 0.76. A relatively weak correlation between and NDVI was found at the low-latitude sites, with a mean R2 of 0.38. There were 13 and 15 sites showing weak linear and exponential correlations between and NDVI, respectively, with a fitted R2 lower than 0.2. These sites were mainly located in Australia and southeast Asia.
3.3 Temporal and spatial accuracy of five G simulation methods
In this study, four criteria were tried to evaluate the model. The results showed that only NSE was suitable for the evaluation of different sites and time periods, whereas RE and KGE were not suitable for the evaluation of different sites, and the RMSE was unsuitable for the evaluation of different time periods. The RE, RMSE, and KGE simulated by the LC method at each site and time period are shown in Fig. 5. The RE values tended to be affected by the mean value of the referenced G in Eq. (7). The RE values of all sites were within ± 3 % (Fig. 5a) during the daytime periods from 07:00 to 16:00 because of the relatively large mean value of G during these periods (Fig. 1a). However, when the mean value of the referenced G was low during the period from 17:00 to 06:30, the RE value at each site was generally greater than that during the daytime periods. In particular, if the mean value of the referenced G is much lower (e.g., close to 0) in a half-hour period for a site, even a small simulation bias could result in an extremely large RE value; for example, the RE of some sites exceeded 8000 % at 05:00 and 05:30 (Fig. 5a). The mean RE value of all sites was also too high (e.g., 300 %) during these periods. Therefore, RE is unsuitable for evaluating different time periods or sites. However, the median RE values of all sites may be more robustly used for evaluation than the mean values. Because the RE is included in the formula of the KGE, the shortcomings of the RE are introduced into the KGE. Therefore, KGE is also unsuitable for the evaluation of different time periods or sites. According to Figs. 1a and 5b, there was a positive relationship between the RMSE and G values during half-hour periods in the intra-day period. The RMSE values were directly affected by the G values. Due to the significant variations in G values in each period, the RMSE was not suitable for a comparison evaluation of the simulation accuracy between different periods. The RE, RMSE, and KGE simulated by the LC method in each land cover type are shown in Fig. S1 in the Supplement.
The NSE was used to evaluate the accuracy of the five G simulation methods at different sites and at half-hour intervals. Daily series were randomly assigned to one of two datasets: 80 % were assigned to the calibration dataset, and 20 % were assigned to the validation dataset. The process of random assignment was repeated to generate 100 independent datasets. Results showed that there was little difference between the performance of the models in the calibration and validation datasets. It indicated that these methods are robust. The performance of the LC_NDVI_P in the calibration and validation datasets at some sites can be found in Fig. S2 in the Supplement. Figure 6 shows the NSE values simulated by the LC, LC_NDVI_P, LC_NDVI_E, LC_fc_SE, and LC_fc_ST methods. In general, the performance of each method varies significantly among different sites and time periods. The simulation accuracy of each method showed high consistency among the different half-hour intervals within the intra-day period. It was highest in the period of 06:00–07:00, followed by the period of 17:00–18:00, whereas it was lowest in 10:00–15:30. Regarding the different methods, the accuracy of the first three methods (LC, LC_NDVI_P, and LC_NDVI_E) was significantly higher than that of the last two methods (LC_fc_SE and LC_fc_ST).
The LC method generally demonstrated its ability to accurately simulate the daily series of G at the site scale in each half-hour period, with the NSE of most sites exceeding 0.40. The mean and median NSE values of all sites were generally higher than 0.50 in each time period. It showed the highest accuracy at 06:00 and 06:30, with the NSE of most sites being above 0.70. The accuracy was the lowest from 11:30 to 15:00, with the mean and median NSE values of all stations being between 0.45 and 0.50. The mean and median NSE values of all the sites in the other periods were generally greater than 0.50. This indicates that the simple LC method was able to accurately simulate a half-hourly series of G from Rn at most sites but also lost this ability, with unsatisfactory accuracy, at a few sites.
Although NDVI was embedded in the LC method, the performances of the LC_NDVI_P and LC_NDVI_E methods were similar to those of the LC method. In other words, the consideration of NDVI resulted in limited improvement in the accuracy of the LC method. In addition, the accuracy of the LC_NDVI_E method was significantly lower than that of the first two methods at 14:30 and 15:00, with mean NSE values of only 0.28 and 0.33, respectively. This was because of the low accuracy (NSE < 0.2) of the LC_NDVI_E method at more sites during these two periods.
The accuracy of the LC_fc_SE and LC_fc_ST methods based on fractional vegetation coverage was relatively low (Fig. 6d and e). The two methods showed little difference across sites and periods. The two methods were able to accurately simulate G only at 05:30–07:00 and 17:30–18:00, with the median NSE values of all sites exceeding 0.5. Conversely, the performance was poor in most night and daytime periods, such as 20:00–04:30 and 08:30–15:30. The NSE of most sites was below 0.4, and the mean and median NSE values of all sites were below 0.2. This indicates that the application of these two methods, considering fractional vegetation coverage, requires caution in the G simulation.
Figure 7 shows the NSE simulated by each method in seven land cover types. The intra-day performance of each land cover type was similar to that of all sites, except for the other type, with the highest simulation accuracy in the periods of 06:00–07:00 and 17:00–18:00. The intra-day accuracy varied most greatly at the forest and savanna sites. The median NSE of all sites simulated by the LC_NDVI_E method was close to 0.8 in the period of 06:00–07:00, while the corresponding NSE was only approximately 0.4. It varied little at other land cover types, especially for wetland and shrubland types. The greatest and lowest values of median NSE for all sites simulated by the LC_NDVI_E method were approximately 0.7 and 0.6, respectively. The NSE of the LC, LC_NDVI_P, and LC_NDVI_E methods showed a unimodal distribution in the other-type sites. The NSE was significantly higher in the daytime than during night periods. The highest value was in the morning and noon periods, with the median NSE of all sites exceeding 0.8. The model performance was significantly better than other land cover types. In the other-type sites, the LC_NDVI_E method performed better than other methods, with the median NSE being higher than 0.6 in each time period.
The spatial distribution of the NSE simulated by the LC method at each site is shown in Fig. 8. Overall, there were significant differences in the performance of the LC method among sites, with the lowest and highest NSE values of each site being −0.37 and 0.94, respectively. As for the median NSE of 48 half-hour periods at each site (Fig. 8a), the performance of the LC method was satisfactory at most sites, with the mean NSE of all sites being 0.58. The NSE values of 70 % of the sites (160/230) were higher than 0.5. However, 27 sites had NSE values lower than 0.4, and 5 sites had NSE values lower than 0.2, indicating that the performance was poor at these sites. For different latitudes, the performance was generally satisfactory at the middle and high latitudes, with NSE values above 0.4. The best performance was observed at high latitudes, with a mean NSE value of 0.69. The accuracy of the LC method was generally low at tropical sites, with a mean NSE of 0.47. The performance was generally satisfactory at most sites in Eurasia and North America, with NSE values higher than 0.5. The NSE values of many sites exceeded 0.7 in these regions. Conversely, relatively poor performance was found in the African, South American, and Oceanian sites, especially in the African sites.
The LC method accurately simulated G at 06:30 in most sites (Fig. 8b), with the mean and median NSE values of all sites being 0.73 and 0.78, respectively. The sites with the NSE higher than 0.5 and 0.6 took up 92 % (211/230) and 84 % (193/230), respectively. The NSE was higher than 0.5 at all sites of Eurasia and North America, and the NSE of most sites exceeded 0.7. The method was also able to accurately simulate G at 18:00 in most sites (Fig. 8e), with a mean NSE of all sites being 0.62. A total of 82 % (188/230) and 59 % (135/230) of sites had NSE values exceeding 0.5 and 0.6, respectively.
The LC method performed poorly at many sites at 10:30 (Fig. 8c) and 13:30 (Fig. 8d), where the mean NSE values of all sites were 0.49 and 0.47, respectively. The sites with an NSE exceeding 0.5 only accounted for 51 % (118/230) and 43 % (100/230) during the two half-hour periods, respectively. The method still performed well at high-latitude sites, with mean NSEs of 0.69 and 0.65, respectively. Conversely, it lost the ability to accurately simulate G at tropical sites. The NSE values at most sites were lower than 0.5, with a mean NSE of only 0.29. The performance of the method was poor at the African and South American sites, with NSE values of each site being below 0.5.
For different land cover types, the LC method performed better in the cropland-, wetland-, and other-type sites. The mean values of the median NSE of wetland and other sites were 0.66 and 0.69, respectively. The method was also able to accurately simulate G in the forest-, grassland-, and shrubland-type sites, with the corresponding mean NSE of 0.57 or 0.56. The method performed the worst at the savanna sites, with the corresponding mean NSE being only 0.47. Since the savanna sites are mainly distributed in tropical regions, this is consistent with the relatively poor performance of the tropical-region site, as mentioned above. The performance of the method varied significantly in each land cover type, except in the other type sites. In the wetland-type sites, there were three sites in the United States with an NSE value lower than 0.3. The NSE of the other 35 sites was higher than 0.50, with the highest value being close to 0.90. The grassland sites were distributed in Asia, Europe, North America, and Oceania. The NSE value was greater than 0.5 at each grassland site in Europe. Cropland sites were distributed in Asia, Europe, and the United States. The NSE value was lower than 0.60 at eight sites in the United States, with a mean NSE value of only 0.45. The method was able to accurately simulate G at 11 sites in Europe, except for one site in Mediterranean region, with a mean NSE value of 0.74. The NSEs for the two Asian sites were 0.54 and 0.71, respectively.
According to the evaluation of the five methods mentioned above, LC, LC_NDVI_P, and LC_NDVI_E performed well. Figure 9 shows the optimal values of the parameters of the three methods at each site and at half-hour periods. As for the results of the LC method (Fig. 9a), on the one hand, the optimal parameter values varied significantly in different sites and half-hour periods. The median optimal parameter of all sites was generally around 0.4 in daytime periods and 0.7 at night. The optimal parameter of each site ranged from 0.4 to 1.1 in the night periods. In contrast, the optimal parameters of most sites showed a concentrated distribution. For example, the optimal parameter values of most sites were generally concentrated in the range of 0.25–0.45 during 08:30–16:30, while the corresponding optimal values were generally concentrated in the range of 0.6–0.85 during 19:30–05:30. The median optimal values of all sites were stable at about 0.35 and 0.75 during these two time periods, respectively. The optimal parameter values of the LC_NDVI_P method showed little difference from those of the LC method at each site and half-hour period (Fig. 9b). The optimal parameter (b in Eq. 4) values of the LC_NDVI_E method ranged from 0.01 to 1.4 among different sites and half-hour periods (Fig. 9c). Similar to the LC and LC_NDVI_P methods, the optimal parameter values of the LC_NDVI_E method also showed a concentrated distribution in each period, especially during 7:30–15:00 and 19:30–05:00. Almost all sites with the optimal parameter values ranged from 0.01 to 0.4, and the optimal values of most sites were concentrated in the range of 0.06 to 0.25 during these two periods.
4.1 Limitations and uncertainties
Theoretically, surface energy is balanced. The energy unclosure might be caused mainly by the error of the observed data. The observed G′ instead of G was generally used to investigate the energy balance ratio (Wilson et al., 2002). The energy balance closure problem might be largely caused by the soil heat storage (Foken, 2008). Compared with G, other energy terms can be observed more accurately. The eddy covariance measurements of H and LE are generally considered to be the most accurate observations available and have been widely used as references to evaluate other simulation methods. Eq. (1) makes full use of the surface energy term that can be accurately measured at present. In other words, it assumes that the measurements of Rn, H, and LE are accurate in this study. The uncertainties of the measurements are not considered in this study. However, the uncertainty of G estimated by the SEBR method could be very large at some sites that have a low magnitude of G. Although the majority of sites have G values greater than 10 2 and take up more than 15 % of Rn, there are some sites (20/230) with G values lower than 5 2. There would be great uncertainties in the SEBR method when simulating G values at these sites.
The accuracies of the LC_fc_SE and LC_fc_ST methods, which embed fractional vegetation coverage in the G simulation, were satisfied for monthly simulations but were significantly lower than those of the other three methods in simulating daily values. The weak correlation between and NDVI might be the main reason for the poor performance of these methods. However, the coarse-resolution NDVI data used in this study are not sufficient to represent the scale of flux measurements, especially for sites with heterogeneous land surface. This might be the main reason for this weak correlation. The application of higher-resolution and continuous vegetation index data series is expected to improve the simulation accuracy of these methods. A large error in the G simulation might be induced in the ET modeling process, thereby reducing the accuracy of the ET estimates. In RS–ET models, Rn is generally calculated using radiation balance with RS images and meteorological inputs. However, observed Rn was used for simulating G in this study. In other words, it was assumed that Rn is accurately simulated by the RS–ET models. Therefore, it should be noted that the uncertainty in Rn calculation was also a source of error in G simulations in ET models.
The evaluation results of this study are expected to provide a reference for RS–ET model application and developers. For example, the performance of these methods was good and poor at some sites and time periods and at some land cover types. RS–ET modelers could check the advantage of the models at good-performance regions and find why the models are poor at some other areas, then they should revise the models to improve the accuracy at poor-performance regions. However, how to improve the model to improve the accuracy needs further research.
4.2 Temporal and spatial variations in the simulation of G
Intra-day and daytime G accounted for 19.2 % and 28.8 % of Rn, respectively. This indicates that G plays an important role and cannot be ignored in RS–ET models. Ignoring G would lead to a great deviation in ET estimates, including the term Rn − G. The intra-day and daytime G′ values were only 4.7 % and 24.1 % of the corresponding G, respectively. This indicates that the measured G′ must be carefully used in the application of the surface energy balance equation. If G′ was used instead of G in the equation, the actual G would be largely underestimated. This is consistent with the results of Meyers and Hollinger (2004) and of Lu et al. (2018).
Several methods for simulating G from Rn in energy-balance-based RS–ET applications have been evaluated in this study. , which varied intra-day, was the key to these methods. Santanello and Friedl (2003) provided diurnal covariations in G and Rn. In this study, more concise linear and polynomial functions of the daytime distribution were fitted (Fig. 2) from 230 observation sites worldwide. The linear- and polynomial-fitted R2 values were 0.92 and 0.98, respectively.
The evaluation results show that the accuracy of the G simulation varied significantly in different half-hour intra-day periods. The highest accuracy was exhibited during 06:00–07:00 and 17:00–18:00. It also showed a slightly higher accuracy during night periods than during daytime periods. The lowest accuracy is observed at noon. This is consistent with the correlation between G and Rn, indicating that the accuracy of the G simulation is affected by the correlation between Rn and G. In other words, the stronger the correlation between Rn and G, the higher the accuracy of the simulation of G from Rn. However, RS–ET models are generally applied during daytime periods. For example, MODIS data represent conditions around 10:30 and 13:30, but the simulation accuracy of these two periods is the lowest during intra-day periods. This is an urgent issue to be solved for G simulation in RS–ET applications, which requires the attention of RS–ET modelers.
The performance of G simulation methods also showed significant spatial variation. The accuracy of the G simulation varied significantly among the observation sites, with the corresponding NSE ranging from 0.2 to 0.9. The G simulation of most sites showed high accuracy in most half-hour periods. This verified the reliability of global RS–ET products in these regions because the more accurate G simulation provided a guarantee for an accurate ET simulation, such as in Eurasia and North America. However, there were also some sites with low simulation accuracy, such as most sites in Africa, South America, and Oceania. A large error in the G simulation would be induced in the ET simulation results and reduce the reliability of ET estimates. Therefore, the application of RS–ET estimates and products in these areas needs more caution for its accuracy. The spatial distribution of the model accuracy was also consistent with the correlation between G and Rn. The sites with satisfied model performance were generally characterized by seasonal variations in G and Rn due to the climate in these regions. Conversely, the sites with poor model performance showed little seasonality. Whether G and Rn have seasonal variations was also an important factor affecting the accuracy of empirical methods in simulating G. This was consistent with the evaluation results of the LE simulation accuracy (Liu, 2021; 2022). This study analyzed the performance of the five methods in seven land cover types. The methods performed better in the wetland- and other-type sites than in other ones. RS–ET modellers are recommended to take more cautions in G simulation at the savanna-type sites because the simulation accuracy was generally lowest at this site type. In addition to climate and land cover factors, regional soil and other environmental factors might also be important factors affecting the accuracy of G simulation. More evaluations of G simulation at a regional scale are recommended for further research.
4.3 Applicability of common G simulation methods in RS–ET models
Daytime RS images are generally applied in ET models. The evaluation results of 230 worldwide observation sites showed that the optimal parameter values of most sites were generally concentrated in the range of 0.25–0.45 during daytime periods, with the median of all sites being stable at approximately 0.34. This indicates that the coefficient values applied in most ET models were reasonable, but the coefficient values applied in the GLEAM model were relatively low.
In the RS–ET models, fixed empirical parameters were applied to the global terrestrial G simulation. The fixed parameters might be suitable for some regions but not on a global scale. This study confirmed that the optimal parameter values vary significantly from site to site. Fixed parameter values induced large errors in the G simulations in other regions. Therefore, it is recommended that model developers consider the spatial variations of G simulation parameters in RS–ET modeling on a global scale.
Some RS–ET models embed vegetation indices (e.g., NDVI, leaf area index (LAI) , or fractional vegetation coverage) or land surface temperature (LST) into the coefficient of the LC method, such as the SEBAL and SEBS models. Evaluation of the LC_NDVI_P and LC_NDVI_E methods, which are the LC methods embedded by the NDVI, showed that the improvement in the simulation accuracy was limited by considering the NDVI. The mismatch between flux observations at the site scale and vegetation index data at the grid scale may be one of the reasons for this result. In addition, the term containing the NDVI in these methods could be taken as a whole, which is similar to the coefficient in the LC method. Therefore, the performance of this method is expected to differ slightly from that of the LC method when the parameter is optimally calibrated. Other models embedded by the LAI (Choudhury, et al., 1987; Allen et al., 2011) and LST (Bastiaanssen, 1995; Faridatul et al., 2020), which were not evaluated in this study due to data limitations, may have similar performances, such as the METRIC model embedded by the LAI and the modified SEBAL model (Faridatul et al., 2020) embedded by the LST and NDVI. Saadi et al. (2018) evaluated three such methods using data from a single observation site. Their results showed that the accuracy order was as follows: (1) Bastiaanssen's (1995) method, (2) Jackson et al.'s (1987) method, and (3) Choudhury et al.'s (1987) method, with RMSE values of 0.09, 0.15, and 0.2, respectively. Evaluation of such methods embedded with the LST data is recommended for further research where data is available. The results of this study indicate that the performance of the different methods varied at some sites. However, the differences among the methods were not significant at global sites as a whole.
In general, the LC method and the methods embedded with the NDVI accurately simulated half-hourly G series at most global sites. There was little difference in simulation accuracy between the different models. However, the performance was poor at some sites. Moreover, the optimal values of the model parameters differed among the different sites. This has also verified by Chen et al. (2019). These issues need to be considered in RS–ET models to improve simulation accuracy. The optimal parameter values for each site showed relative stability between different half-hour periods in the daytime, indicating that it was feasible to apply the same coefficient value in different daytime periods.
Intra-day and daytime G accounted for 19.2 % and 28.8 % of Rn, respectively. This indicates that G plays an important role and cannot be ignored in RS–ET models. The accuracy of the five G simulation methods in energy-balance-based RS–ET models was evaluated using half-hourly observations from 230 flux sites. The LC method and the methods embedded with the NDVI could accurately simulate a half-hourly G series at most sites. The mean and median NSE values of all sites were generally higher than 0.50 in each half-hour period. However, the two methods embedded by fractional vegetation coverage in the G simulation showed poor performance in most half-hour periods, except for the periods of 06:00–07:00 and 17:00–18:00, with mean and median NSE values of all sites being below 0.20. The poor performance might be caused mainly by the coarse-resolution vegetation index data, which could not represent the scale of flux measurements. The performance of each method was generally consistent at different sites and time periods.
The accuracy of each method varied significantly at different sites and at half-hour intervals. The highest accuracy was exhibited during 06:00–07:00, followed by the period of 17:00–18:00. A total of 92 % (211/230) of sites exhibited an NSE of the LC method greater than 0.50 at 06:30. It showed a slightly higher accuracy during night periods than during daytime periods. The lowest accuracy was observed at noon periods (10:00–15:30). For example, the sites with an NSE exceeding 0.50 only accounted for 51 % (118/230) and 43 % (100/230) at 10:30 and 13:30, respectively. The NSE values of the different sites ranged from −0.37 to 0.94. The accuracy of the Northern Hemisphere sites was generally higher than that of the Southern Hemisphere sites. In general, it showed the highest accuracy at high-latitude sites and then at middle-latitude sites, while it exhibited the lowest accuracy at low-latitude sites, especially at tropical sites. As for the median NSE of 48 half-hour periods in the LC method, the mean NSE values of the high latitudes and tropical sites were 0.69 and 0.47, respectively. The performance of the LC method and the methods embedded with the NDVI were generally satisfactory at the Eurasian and North American sites, with the NSE values of most sites exceeding 0.70. Conversely, they exhibited relatively poor performance at the African, South American, and Oceanian sites, especially the African sites. Both the temporal and spatial distributions of the accuracy of the G simulation were positively correlated with the correlation between G and Rn. In other words, the sites or periods with stronger correlations between G and Rn have higher simulation accuracy.
Overall, the LC, LC_NDVI_P, and LC_NDVI_E methods accurately simulated the G series at most observation sites and half-hour periods in the intra-day, with an NSE value exceeding 0.50. However, the performance of these methods was poor at some sites and time periods. This negatively affects the accuracy of energy-balance-based RS–ET simulations. The application of RS–ET datasets covering these sites requires caution. The performance was best in the wetland- and other-type sites and was worst at the savanna-type sites. Improvement of G simulation at low-accuracy regions, such as low-latitude regions and savanna-type sites, is recommended for the RS–ET modelers. The weak correlation between the G and Rn is the physical reason for the poor accuracy of G simulation in these regions and sites. Instead of the Rn, finding another variable that has a physical connection and a strong correlation with G might be a more efficient solution to improve the accuracy of the empirical estimation method for G. In addition, the optimal parameter value of each method varied significantly at different sites. Therefore, the fixed parameter values in the G simulation methods do not match the actual situation. Variable parameters are recommended in empirical methods of G simulation to improve accuracy.
The supplement related to this article is available online at: https://doi.org/10.5194/hess-26-6207-2022-supplement.
The author has declared that there are no competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This study was supported and funded by the National Natural Science Foundation of China (grant no. 42171029) and the Strategic Priority Research Program of the Chinese Academy of Sciences (grant nos. XDA2006020202 and XDA23090302). We are grateful to the Data Portal serving the FLUXNET community (https://fluxnet.org/, last access: 8 December 2022) and to OzFlux (http://www.ozflux.org.au/index.html, last access: 8 December 2022) for providing FLUXNET and TERN OzFlux datasets. We would also like to thank Dario Papale, Gilberto Pastorello, and Jason Beringer for their kind assistance with data access. We thank the editor Bob Su and four anonymous reviewers for their constructive comments, which helped to improve the paper.
This research has been supported by the National Natural Science Foundation of China (grant no. 42171029) and the Chinese Academy of Sciences (grant nos. XDA23090302 and XDA2006020202).
This paper was edited by Bob Su and reviewed by four anonymous referees.
Agam, N., Kustas, W. P., Alfieri, J. G., Gao, F., McKee, L. M., Prueger, J. H., and Hipps, L. E: Micro-scale spatial variability in soil heat flux (SHF) in a winegrape vineyard, Irrigation Sci., 37, 253–268, https://doi.org/10.1007/s00271-019-00634-6, 2019.
Ait Hssaine, B., Merlin, O., Ezzahar, J., Ojha, N., Er-Raki, S., and Khabba, S: An evapotranspiration model self-calibrated from remotely sensed surface soil moisture, land surface temperature and vegetation cover fraction: application to disaggregated SMOS and MODIS data, Hydrol. Earth Syst. Sci., 24, 1781–1803, https://doi.org/10.5194/hess-24-1781-2020, 2020.
Allen, R. G., Tasumi, M., and Trezza, R: Satellite-based energy balance for mapping evapotranspiration with internalized calibration (METRIC)-model, J. Irrig. Drain. E., 133, 380–394, https://doi.org/10.1061/(ASCE)0733-9437(2007)133:4(380), 2007.
Allen, R. G., Irmak, A., Trezza, R., Hendrickx, J. M., Bastiaanssen, W., and Kjaersgaard, J: Satellite-based ET estimation in agriculture using SEBAL and METRIC, Hydrol. Process., 25, 4011–4027, https://doi.org/10.1002/hyp.8408, 2011.
Anderson, M. C., Norman, J. M., Diak, G. R., Kustas, W. P., and Mecikalski, J. R: A two-source time-integrated model for estimating surface fluxes using thermal infrared remote sensing, Remote Sens. Environ., 60, 195–216, https://doi.org/10.1016/S0034-4257(96)00215-5, 1997.
Australian Terrestrial Ecosystem Research Network: The OzFlux Data Portal, http://data.ozflux.org.au, last access: 8 December 2022.
Bastiaanssen, W. G. M: Regionalization of surface flux densities and moisture indicators in composite terrain; a remote sensing approach under clear skies in mediterranean climates, SC-DLO, Wageningen, 1995.
Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., and Holtslag, A. A. M: A remote sensing surface energy balance algorithm for land (SEBAL): 1. Formulation, J. Hydrol., 212–213, 198–212, https://doi.org/10.1016/S0022-1694(98)00254-6, 1998.
Beringer, J., Hutley, L. B., McHugh, I., Arndt, S. K., Campbell, D., Cleugh, H. A., Cleverly, J., Resco de Dios, V., Eamus, D., Evans, B., Ewenz, C., Grace, P., Griebel, A., Haverd, V., Hinko-Najera, N., Huete, A., Isaac, P., Kanniah, K., Leuning, R., Liddell, M. J., Macfarlane, C., Meyer, W., Moore, C., Pendall, E., Phillips, A., Phillips, R. L., Prober, S. M., Restrepo-Coupe, N., Rutledge, S., Schroder, I., Silberstein, R., Southall, P., Yee, M. S., Tapper, N. J., van Gorsel, E., Vote, C., Walker, J., and Wardlaw, T: An introduction to the Australian and New Zealand flux tower network – OzFlux, Biogeosciences, 13, 5895–5916, https://doi.org/10.5194/bg-13-5895-2016, 2016.
Bonsoms, J. and Boulet, G: Ensemble machine learning outperforms empirical equations for the ground heat flux estimation with remote sensing data, Remote Sens.-Basel, 14, 1788, https://doi.org/10.3390/rs14081788, 2022.
Cellier, P., Richard, G., and Robin, P: Partition of sensible heat fluxes into bare soil and the atmosphere, Agr. Forest Meteorol., 82, 245–265, https://doi.org/10.1016/0168-1923(95)02328-3, 1996.
Chen, J., and Liu, J: Evolution of evapotranspiration models using thermal and shortwave remote sensing data, Remote Sens. Environ., 237, 111594, https://doi.org/10.1016/j.rse.2019.111594, 2020.
Chen, X., Su, Z., Ma, Y., Liu, S., Yu, Q., and Xu, Z: Development of a 10-year (2001–2010) 0.1∘ data set of land-surface energy balance for mainland China, Atmos. Chem. Phys., 14, 13097–13117, https://doi.org/10.5194/acp-14-13097-2014, 2014.
Chen, X., Su, Z., Ma, Y., and Middleton, E. M: Optimization of a remote sensing energy balance method over different canopy applied at global scale, Agr. Forest Meteorol., 279, 107633, https://doi.org/10.1016/j.agrformet.2019.107633, 2019.
Chen, X., Su, Z., Ma, Y., Trigo, I., and Gentine, P: Remote sensing of global daily evapotranspiration based on a surface energy balance method and reanalysis data, J. Geophys. Res.-Atmos., 126, e2020JD032873, https://doi.org/10.1029/2020JD032873, 2021.
Choudhury, B. J., Idso, S. B., and Reginato, R. J: Analysis of an empirical model for soil heat flux under a growing wheat crop for estimating evaporation by an infrared-temperature based energy balance equation, Agr. Forest Meteorol., 39, 283–297, https://doi.org/10.1016/0168-1923(87)90021-9, 1987.
Faridatul, M. I., Wu, B., Zhu, X., and Wang, S: Improving remote sensing based evapotranspiration modelling in a heterogeneous urban environment, J. Hydrol., 581, 124405. https://doi.org/10.1016/j.jhydrol.2019.124405, 2020.
FLUXNET: The FLUXNET2015 dataset and FLUXNET-CH4 Community Product, https://fluxnet.org/data/download-data/, last access: 8 December 2022.
Foken, T: The energy balance closure problem: An overview, Ecol. Appl., 18, 1351–1367, https://doi.org/10.1890/06-0922.1, 2008.
Gao, Z., Russell, E. S., Missik, J. E. C., Huang, M., Chen, X., Strickland, C. E., Clayton, R., Arntzen, E., Ma, Y., and Liu, H: A novel approach to evaluate soil heat flux calculation: An analytical review of nine methods, J. Geophys. Res.-Atmos., 122, 6934–6949, https://doi.org/10.1002/2017JD027160, 2017.
Gentine, P., Entekhabi, D., and Heusinkveld, B: Systematic errors in ground heat flux estimation and their correction, Water Resour. Res., 48, W09541, https://doi.org/10.1029/2010WR010203, 2012.
Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modeling, J. Hydrol., 377, 80–91, https://doi.org/10.1016/j.jhydrol.2009.08.003, 2009.
Heusinkveld, B. G., Jacobs, A. F. G., Holtslag, A. A. M., and Berkowicz, S. M: Surface energy balance closure in an arid region: role of soil heat flux, Agr. Forest Meteorol., 122, 21–37, https://doi.org/10.1016/j.agrformet.2003.09.005, 2004.
Jackson, R. D., Moran, M. S., Gay, L. W., and Raymond, L. H.: Evaluating evaporation from field crops using airborne radiometry and ground-based meteorological data, Irrigation Sci., 8, 81–90, https://doi.org/10.1007/bf00259473, 1987.
Jiang, C. and Ryu, Y: Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS), Remote Sens. Environ., 186, 528–547, https://doi.org/10.1016/j.rse.2016.08.030, 2016.
Ju, Z. and Hu, C: Experimental warming alters soil hydro-thermal properties and heat flux in a winter wheat field, Arch. Agron. Soil Sci., 64, 718–730, https://doi.org/10.1080/03650340.2017.1361021, 2018.
Kalma, J. D., McVicar, T. R., and McCabe, M. F: Estimating land surface evaporation: A review of methods using remotely sensed surface temperature data, Surv. Geophys., 29, 421–469, https://doi.org/10.1007/s10712-008-9037-z, 2008.
Knox, S. H., Jackson R. B., Poulter B., McNicol G., Fluet-Chouinard E., Zhang Z., Hugelius G., Bousquet, P., Canadell, J. G., Saunois, M., Papale, D., Chu, H., Keenan, T. F., Baldocchi, D., Torn, M. S., Mammarella, I., Trotta, C., Aurela, M., Bohrer, G., Campbell, D. I., Cescatti, A., Chamberlain, S., Chen, J., Chen, W., Dengel, S., Desai, AR., Euskirchen, E., Friborg, T., Gasbarra, D., Goded, I., Goeckede, M., Heimann, M., Helbig, M., Hirano, T., Hollinger, D. Y., Iwata, H., Kang, M., Klatt, J., Krauss, K. W., Kutzbach, L., Lohila, A., Mitra, B., Morin, T. H., Nilsson, M. B., Niu, S., Noormets, A., Oechel, W. C., Peichl, M., Peltola, O., Reba, M. L., Richardson, A. D., Runkle, B. R. K., Ryu, Y., Sachs, T., Schafer, K. V. R., Schmid, H. P., Shurpali, N., Sonnentag, O., Tang, A. C. I., Ueyama, M., Vargas, R., Vesala, T., Ward, E. J., Windham-Myers, L., Wohlfahrt, G., and Zona, D: FLUXNET-CH4 synthesis activity: objectives, observations, and future directions, B. Am. Meteorol. Soc., 100, 2607–2632, https://doi.org/10.1175/BAMS-D-18-0268.1, 2019.
Kustas, W. P. and Daughtry, C. S. T.: Estimation of the soil heat flux/net radiation ratio from spectral data, Agr. Forest Meteorol., 49, 205–223, https://doi.org/10.1016/0168-1923(90)90033-3, 1990.
Kustas, W. P., Prueger, J. H., Hatfield, J. L., Ramalingam, K., and Hipps, L. E: Variability in soil heat flux from a mesquite dune site, Agr. Forest Meteorol., 103, 249–264, https://doi.org/10.1016/S0168-1923(00)00131-3, 2000.
Leuning, R., van Gorsel, E., Massman, W. J., and Isaac, P. R: Reflections on the surface energy imbalance problem, Agr. Forest Meteorol., 156, 65–74, https://doi.org/10.1016/j.agrformet.2011.12.002, 2012.
Liebethal, C. and Foken, T: Evaluation of six parameterization approaches for the ground heat flux, Theor. Appl. Climatol., 88, 43–56, https://doi.org/10.1007/s00704-005-0234-0, 2007.
Liu, Z: The accuracy of temporal upscaling of instantaneous evapotranspiration to daily values with seven upscaling methods, Hydrol. Earth Syst. Sci., 25, 4417–4433, https://doi.org/10.5194/hess-25-4417-2021, 2021.
Liu, Z: Estimating land evapotranspiration from potential evapotranspiration constrained by soil water at daily scale, Sci. Total Environ., 834, 155327, https://doi.org/10.1016/j.scitotenv.2022.155327, 2022.
Liu, Z., Yao, Z., and Wang, R: Evaluation and validation of CryoSat-2-derived water levels using in situ lake data from China, Remote Sens.-Basel, 11, 899, https://doi.org/10.3390/rs11080899, 2019.
Lu, S., Wang, H., Meng, P., Zhang, J., and Zhang, X: Determination of soil ground heat flux through heat pulse and plate methods: Effects of subsurface latent heat on surface energy balance closure, Agr. Forest Meteorol., 260–261, 176–182, https://doi.org/10.1016/j.agrformet.2018.06.008, 2018.
Mayocchi, C. L. and Bristowa, K. L: Soil surface heat flux: some general questions and comments on measurements, Agr. Forest Meteorol., 75, 43–50, https://doi.org/10.1016/0168-1923(94)02198-S, 1995.
Merlin, O., Chirouze, J., Olioso, A., Jarlan, L., Chehbouni, G., and Boulet, G: An image-based four-source surface energy balance model to estimate crop evapotranspiration from solar reflectance/thermal emission data (SEB-4S), Agr. Forest Meteorol., 184, 188–203, https://doi.org/10.1016/j.agrformet.2013.10.002, 2014.
Meyers, T. P. and Hollinger, S. E: An assessment of storage terms in the surface energy balance of maize and soybean, Agr. Forest Meteorol., 125, 105–115, https://doi.org/10.1016/j.agrformet.2004.03.001, 2004.
Miralles, D. G., Holmes, T. R. H., De Jeu, R. A. M., Gash, J. H., Meesters, A. G. C. A., and Dolman, A. J: Global land-surface evaporation estimated from satellite-based observations, Hydrol. Earth Syst. Sci., 15, 453–469, https://doi.org/10.5194/hess-15-453-2011, 2011.
Mohan, M. M. P., Kanchirapuzha, R., and Varma, M. R. R: Review of approaches for the estimation of sensible heat flux in remote sensing-based evapotranspiration models, J. Appl. Remote Sens., 14, 041501, https://doi.org/10.1117/1.JRS.14.041501, 2020.
Mu, Q. Z., Zhao, M. S., and Running, S. W: Improvements to a MODIS global terrestrial evapotranspiration algorithm, Remote Sens. Environ., 115, 1781–1800, https://doi.org/10.1016/j.rse.2011.02.019, 2011.
Nash, J. E. and Sutcliffe, J. V: River flow forecasting through conceptual models part 1-A discussion of principles, J. Hydrol., 10, 282–290, https://doi.org/10.1016/0022-1694(70)90255-6, 1970.
Norman, J. M., Kustas, W. P., and Humes, K. S: A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature, Agr. Forest Meteorol., 77, 263–293, https://doi.org/10.1016/0168-1923(95)02265-Y, 1995.
Norman, J. M., Anderson, M. C., Kustas, W. P., French, A. N., Mecikalski, J. R., Torn, R. D., Diak, G. R, Schmugge, T. J., and Tanner, B. C. W: Remote sensing of surface energy fluxes at 101-m pixel resolutions, Water Resour. Res., 39, 1221, https://doi.org/10.1029/2002WR001775, 2003.
Pastorello, G., Trotta, C., Canfora, E., Chu H., Christianson D., Cheah Y. W., Poindexter, C., Chen, J., Elbashandy, A., Humphrey, M., Isaac, P., Polidori, D., Ribeca, A., van Ingen, C., Zhang, L., Amiro, B., Ammann, C., Arain, M. A., Ardö, J., Arkebauer, T., Arndt, S. K., Arriga, N., Aubinet, M., Aurela, M., Baldocchi, D., Barr, A., Beamesderfer, E., Marchesini, L. B., Bergeron, O., Beringer, J., Bernhofer, C., Berveiller, D., Billesbach, D., Black, T. A., Blanken, P. D., Bohrer, G., Boike, J., Bolstad, P. V., Bonal, D., Bonnefond, J. M., Bowling, D. R., Bracho, R., Brodeur, J., Brümmer, C., Buchmann, N., Burban, B., Burns, S. P., Buysse, P., Cale, P., Cavagna, M., Cellier, P., Chen, S., Chini, I., Christensen, T. R., Cleverly, J., Collalti, A., Consalvo, C., Cook, B. D., Cook, D., Coursolle, C., Cremonese, E., Curtis, P. S., D'Andrea, E., da Rocha, H., Dai, X., Davis, K. J., De Cinti, B., de Grandcourt, A., De Ligne, A., De Oliveira, R. C., Delpierre, N., Desai, A. R., Di Bella, C. M., di Tommasi, P., Dolman, H., Domingo, F., Dong, G., Dore, S., Duce, P., Dufrêne, E., Dunn, A., Dušek, J., Eamus, D., Eichelmann, U., ElKhidir, H. A. M., Eugster, W., Ewenz, C. M., Ewers, B., Famulari, D., Fares, S., Feigenwinter, I., Feitz, A., Fensholt, R., Filippa, G., Fischer, M., Frank, J., Galvagno, M., Gharun, M., Gianelle, D., Gielen, B., Gioli, B., Gitelson, A., Goded, I., Goeckede, M., Goldstein, A. H., Gough, C. M., Goulden, M. L., Graf, A., Griebel, A., Gruening, C., Grünwald, T., Hammerle, A., Han, S., Han, X., Hansen, B. U., Hanson, C., Hatakka, J., He, Y., Hehn, M., Heinesch, B., Hinko-Najera, N., Hörtnagl, L., Hutley, L., Ibrom, A., Ikawa, H., Jackowicz-Korczynski, M., Janouš, D., Jans, W., Jassal, R., Jiang, S., Kato, T., Khomik, M., Klatt, J., Knohl, A., Knox, S., Kobayashi, H., Koerber, G., Kolle, O., Kosugi, Y., Kotani, A., Kowalski, A., Kruijt, B., Kurbatova, J., Kutsch, W. L., Kwon, H., Launiainen, S., Laurila, T., Law, B., Leuning, R., Li, Y., Liddell, M., Limousin, J. M., Lion, M., Liska, A. J., Lohila, A., López-Ballesteros, A., López-Blanco, E., Loubet, B., Loustau, D., Lucas-Moffat, A., Lüers, J., Ma, S., Macfarlane, C., Magliulo, V., Maier, R., Mammarella, I., Manca, G., Marcolla, B., Margolis, H. A., Marras, S., Massman, W., Mastepanov, M., Matamala, R., Matthes, J. H., Mazzenga, F., McCaughey, H., McHugh, I., McMillan, A. M. S., Merbold, L., Meyer, W., Meyers, T., Miller, S. D., Minerbi, S., Moderow, U., Monson, R. K., Montagnani, L., Moore, C. E., Moors, E., Moreaux, V., Moureaux, C., Munger, J. W., Nakai, T., Neirynck, J., Nesic, Z., Nicolini, G., Noormets, A., Northwood, M., Nosetto, M., Nouvellon, Y., Novick, K., Oechel, W., Olesen, J. E., Ourcival, J. M., Papuga, S. A., Parmentier, F. J., Paul-Limoges, E., Pavelka, M., Peichl, M., Pendall, E., Phillips, R. P., Pilegaard, K., Pirk, N., Posse, G., Powell, T., Prasse, H., Prober, S. M., Rambal, S., Rannik,Ü., Raz-Yaseef, N., Reed, D., de Dios, V. R., Restrepo-Coupe, N., Reverter, B. R., Roland, M., Sabbatini, S., Sachs, T., Saleska, S. R., Sánchez-Cañete, E. P., Sanchez-Mejia, Z. M., Schmid, H. P., Schmidt, M., Schneider, K., Schrader, F., Schroder, I., Scott, R. L., Sedlák, P., Serrano-Ortíz, P., Shao, C., Shi, P., Shironya, I., Siebicke, L., Šigut, L., Silberstein, R., Sirca, C., Spano, D., Steinbrecher, R., Stevens, R. M., Sturtevant, C., Suyker, A., Tagesson, T., Takanashi, S., Tang, Y., Tapper, N., Thom, J., Tiedemann, F., Tomassucci, M., Tuovinen, J. P., Urbanski, S., Valentini, R., van der Molen, M., van Gorsel, E., van Huissteden, K., Varlagin, A., Verfaillie, J., Vesala, T., Vincke, C., Vitale, D., Vygodskaya, N., Walker, J. P., Walter-Shea, E., Wang, H., Weber, R., Westermann, S., Wille, C., Wofsy, S., Wohlfahrt, G., Wolf, S., Woodgate, W., Li, Y., Zampedri, R., Zhang, J., Zhou, G., Zona, D., Agarwal, D., Biraud, S., Torn, M., and Papale, D: The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data, Sci. Data, 7, 225, https://doi.org/10.1038/s41597-020-0534-3, 2020.
Pauwels, V. R. N. and Daly, E: Advantages of analytically computing the ground heat flux in land surface models, Hydrol. Earth Syst. Sci., 20, 4689–4706, https://doi.org/10.5194/hess-20-4689-2016, 2016.
Peng, X. Y., Heitman, J., Horton, R., and Ren, T. S: Determining near-surface soil heat flux density using the gradient method: a thermal conductivity model-based approach, J. Hydrometeorol., 18, 2285–2295, https://doi.org/10.1175/JHM-D-16-0290.1, 2017.
Ponce-Campos, G. E., Moran, M. S., Huete, A., Zhang, Y., Breslo, C., Huxman, T. E., Eamus, D., Bosch, D. D., Buda, A. R., Gunter, S. A., Scalley, T. H., Kitchen, S. G., McClaran, M. P., McNab, W. H., Montoya, D. S., Morgan, J. A., Peters, D. P. C., Sadler, E. J., Seyfried, M. S., and Starks, P. J: Ecosystem resilience despite large-scale altered hydroclimatic conditions, Nature, 494, 349–352, https://doi.org/10.1038/nature11836, 2013.
Purdy, A. J., Fisher, J. B., Goulden, M. L., and Famiglietti, J. S: Ground heat flux: an analytical review of 6 models evaluated at 88 sites and globally, J. Geophys. Res.-Biogeo., 121, 3045–3059, https://doi.org/10.1002/2016JG003591, 2016.
Roerink, G. J., Su, Z., and Menenti, M: S-SEBI: a simple remote sensing algorithm to estimate the surface energy balance, Phys. Chem. Earth Pt. B, 25, 147–157, https://doi.org/10.1016/S1464-1909(99)00128-8, 2000.
Russell, E. S., Liu, H., Gao, Z., Finn, D., and Lamb, B: Impacts of soil heat flux calculation methods on the surface energy balance closure, Agr. Forest Meteorol., 214–215, 189–200, https://doi.org/10.1016/j.agrformet.2015.08.255, 2015.
Saadi, S., Boulet, G., Bahir, M., Brut, A., Delogu, É., Fanise, P., Mougenot, B., Simonneaux, V., and Lili Chabaane, Z: Assessment of actual evapotranspiration over a semiarid heterogeneous land surface by means of coupled low-resolution remote sensing data with an energy balance model: comparison to extra-large aperture scintillometer measurements, Hydrol. Earth Syst. Sci., 22, 2187–2209, https://doi.org/10.5194/hess-22-2187-2018, 2018.
Sánchez, J. M., Kustas, W. P., Caselles, V., and Anderson, M. C: Modelling surface energy fluxes over maize using a two-source patch model and radiometric soil and canopy temperature observations, Remote Sens. Environ., 112, 1130–1143, https://doi.org/10.1016/j.rse.2007.07.018, 2008.
Santanello, J. A. and Friedl, M. A: Diurnal covariation in soil heat flux and net radiation, J. Appl. Meteorol., 42, 851–862, https://doi.org/10.1175/1520-0450(2003)042<0851:DCISHF>2.0.CO;2, 2003.
Shao, C., Chen, J., Li, L., Xu, W., Chen, S., Gwen, T., Xu, J., and Zhang, W: Spatial variability in soil heat flux at three Inner Mongolia steppe ecosystems, Agr. Forest Meteorol., 148, 1433–1443, https://doi.org/10.1016/j.agrformet.2008.04.008, 2008.
Singh, R. K., Irmak, A., Irmak, S., and Martin, D. L: Application of SEBAL model for mapping evapotranspiration and estimating surface energy fluxes in south-central Nebraska, J. Irrig. Drain. Eng., 134, 273–285, https://doi.org/10.1061/(ASCE)0733-9437(2008)134:3(273), 2008.
Stoy, P. C., Mauder, M., Foken, T., Marcolla, B., Boegh, E., Ibrom, A., Altaf Arain, M., Arneth, A., Aurela, M., Bernhofer, C., Cescatti, A., Dellwik, E., Duce, P., Gianelle, D., van Gorsel, E., Kiely, G., Knohl, A., Margolis, H., McCaughey, H., Merbold, L., Montagnani, L., Papale, D., Reichstein, M., Saunders, M., Serrano-Ortiz, P., Sottocornola, M., Spano, D., Vaccari, F., and Varlagin, A: A data-driven analysis of energy balance closure across FLUXNET research sites: the role of landscape scale heterogeneity, Agr. Forest Meteorol., 171–172, 137–152, https://doi.org/10.1016/j.agrformet.2012.11.004, 2013.
Su, Z: The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes, Hydrol. Earth Syst. Sci., 6, 85–100, https://doi.org/10.5194/hess-6-85-2002, 2002.
Sun, Z., Gebremichael, M., and Wang, Q: Evaluation of empirical remote sensing-based equations for estimating soil heat flux, J. Meteorol. Soc. Jpn., 91, 627–638, https://doi.org/10.2151/jmsj.2013-505, 2013.
van der Tol, C: Validation of remote sensing of bare soil ground heat flux, Remote Sens. Environ., 121, 275–286, https://doi.org/10.1016/j.rse.2012.02.009, 2012.
Verhoef, A., Ottlé, C., Cappelaere, B., Murray, T., Saux-Picart, S., Zribi, M., Maignan, F., Boulain, N., Demarty, J., and Ramier, D: Spatio-temporal surface soil heat flux estimates from satellite data: Results for the AMMA experiment at the Fakara (Niger) supersite, Agr. Forest Meteorol., 154–155, 55–66, https://doi.org/10.1016/j.agrformet.2011.08.003, 2012.
Vermote, E., Justice, C., Csiszar, I., Eidenshink, J., Myneni, R., Baret, F., Masuoka, E., Wolfe, R., and Claverie, M: NOAA Climate Data Record (CDR) of Normalized Difference Vegetation Index (NDVI), Version 4, NOAA National Centers for Environmental Information, https://doi.org/10.7289/V5PZ56R6, 2014.
Wang, Z. H. and Bou-Zeid E: A novel approach for the estimation of soil ground heat flux, Agr. Forest Meteorol., 154–155, 214–221, https://doi.org/10.1016/j.agrformet.2011.12.001, 2012.
Wilson, K. B., Goldstein, A., Falge, E., Aubinet, M., Baldocchi, D., Berbigier, P., Bernhofer, C., Ceulemans, R., Dolman, H., Field, C., Grelle, A., Ibrom, A., Law, B. E., Kowalski, A., Meyers, T., Moncrieff, J., Monson, R., Oechel, W., Tenhunen, J., Valentini, R., and Verma, S: Energy balance closure at FLUXNET sites, Agr. Forest Meteorol., 113, 223–243, https://doi.org/10.1016/S0168-1923(02)00109-0, 2002.
Wu, B., Oncley, S. P., Yuan, H., and Chen, F: Ground heat flux determination based on near-surface soil hydro-thermodynamics, J. Hydrol., 591, 125578, https://doi.org/10.1016/j.jhydrol.2020.125578, 2020.
Yang, K. and Wang, J. M: A temperature prediction-correction method for estimating surface soil heat flux from soil temperature and moisture data, Sci. China Ser. D, 51, 721–729, https://doi.org/10.1007/s11430-008-0036-1, 2008.
Yue, P., Zhang, Q., Niu, S., Cheng, H., and Wang, X: Effects of the soil heat flux estimates on surface energy balance closure over a semi-arid grassland, Acta Meteorol. Sin., 25, 774–782. https://doi.org/10.1007/s13351-011-0608-4, 2011.
Zhang, K., Kimball, J. S., Nemani, R. R., and Running, S. W: A continuous satellite-derived global record of land surface evapotranspiration from 1983 to 2006, Water Resour. Res., 46, W09522, https://doi.org/10.1029/2009WR008800, 2010.
Zhang, K., Kimball, J. S., and Running, S. W: A review of remote sensing based actual evapotranspiration estimation, WIREs Water, 3, 834–853, https://doi.org/10.1002/wat2.1168, 2016.