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  R_n, 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 R_n is the main reason for the poor performance at these regions.

1 Introduction

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 (R_n) (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 R_n, 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 R_n. Choudhury et al. (1987) established an empirical function between G and R_n for bare and vegetated soils. They found that the ratio of G to R_n ($G/R_{\mathrm{n}}$) 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 $G/R_{\mathrm{n}}$ 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 $G/R_{\mathrm{n}}$ 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 R_n 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 R_n; (2) evaluating the accuracy of five empirical methods in simulating half-hourly G from R_n; 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.

2 Materials and methods

2.1 Data

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, R_n, 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 Methods

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:

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

where G is ground heat flux, R_n 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 R_n 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:

$$\begin{array}{}\text{(2)}& {\displaystyle }G=\mathit{\alpha }\cdot R_{\mathrm{n}},\text{(3)}& {\displaystyle }G=\mathit{\alpha }\cdot [\mathrm{1}-\mathrm{0.98}\cdot {\mathrm{NDVI}}^{\mathrm{4}}]\cdot R_{\mathrm{n}},\text{(4)}& {\displaystyle }G=[\mathit{\alpha }\cdot \exp (-b\cdot \mathrm{NDVI}\left)\right]\cdot R_{\mathrm{n}},\text{(5)}& {\displaystyle }G=[{\mathit{\alpha }}_{\mathrm{1}}+({\mathit{\alpha }}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{1}})\cdot (\mathrm{1}-\mathrm{fc}\left)\right]\cdot R_{\mathrm{n}},\text{(6)}& {\displaystyle }G=\mathit{\alpha }\cdot (\mathrm{1}-\mathrm{fc})\cdot R_{\mathrm{n}},\end{array}$$

where α, α₁, α₂, 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:

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathrm{RE}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{\displaystyle \frac{X_{\text{si}}-X_{\text{ri}}}{X_{\text{ri}}}},\text{(8)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left[\right(X_{\text{si}}-{\stackrel{\mathrm{‾}}{X}}_{\mathrm{s}})-(X_{\text{ri}}-{\stackrel{\mathrm{‾}}{X}}_{\mathrm{r}}){]}^{\mathrm{2}}},\text{(9)}& {\displaystyle }\mathrm{KGE}=\mathrm{1}-\sqrt{(\mathrm{1}-\mathrm{CC}{)}^{\mathrm{2}}+{\mathrm{RE}}^{\mathrm{2}}+{\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{SD}}_{\mathrm{s}}}{{\mathrm{SD}}_{\mathrm{r}}}}\right)}^{\mathrm{2}}},\text{(10)}& {\displaystyle }\mathrm{SD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{\sum }_{i=\mathrm{1}}^n(X_i-\stackrel{\mathrm{‾}}{X}{)}^{\mathrm{2}}},\text{(11)}& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(X_{\text{si}}-X_{\text{ri}}{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(X_{\text{ri}}-{\stackrel{\mathrm{‾}}{X}}_{\mathrm{r}}{)}^{\mathrm{2}}}},\end{array}$$

where X_si and X_ri are the ith values of the simulated and referenced G time series, respectively; n is the time series length; ${\stackrel{\mathrm{‾}}{X}}_{\mathrm{s}}$ and ${\stackrel{\mathrm{‾}}{X}}_{\mathrm{r}}$ are the means of the simulated and referenced G, respectively; and SD_s and SD_r are the standard deviations of the simulated and referenced G, respectively.

Figure 1 Intra-day distribution of raw and normalized R_n, H, LE, G, G^′, and the values from the sine and Gaussian functions (a is raw R_n, H, LE, G, and G^′; b–f are normalized R_n, H, LE, G, and G^′, respectively).

3 Results

3.1 Intra-day distribution of observed surface energy balance items and $G/R_{\mathrm{n}}$

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 R_n, 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 R_n, respectively. G accounted for 28.8 % of R_n 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 R_n (Fig. 1a). However, the intra-day distributions of LE, G, and G^′ showed an overall deviation from the measured R_n. The distribution of LE and G^′ was generally delayed by half an hour compared to the measured R_n, 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.

Figure 2 Intra-day distribution of the ratio of G and R_n ($G/R_{\mathrm{n}}$) and the ratio of G and H ($G/H$) at each site (a and b are intra-day distributions of $G/R_{\mathrm{n}}$ and $G/H$ ratios at each site, respectively; c and d are the fitted lines of the mean values of $G/R_{\mathrm{n}}$ and $G/H$ in the daytime, respectively).

The intra-day distribution characteristics of $G/R_{\mathrm{n}}$ and the ratio of G and H ($G/H$) were also analyzed based on field observation data. The intra-day distributions of $G/R_{\mathrm{n}}$ and $G/H$ at each site are shown in Fig. 2. During data processing, data points with absolute values greater than 10 in the $G/R_{\mathrm{n}}$ or $G/H$ daily series of each period were deleted. Outliers in the $G/R_{\mathrm{n}}$ or $G/H$ series were then removed using the outlier detection method. The $G/R_{\mathrm{n}}$ varied significantly in different half-hour periods of the intra-day and among the different sites (Fig. 2a). The variation in $G/R_{\mathrm{n}}$ among the sites was lower during the daytime than that at night. The variation range of $G/R_{\mathrm{n}}$ among all sites was generally approximately 0.2 in the daytime. This indicates that the $G/R_{\mathrm{n}}$ of all sites showed high consistency during these periods. At night, the variation mostly ranged between 0.6 and 0.8. In other words, $G/R_{\mathrm{n}}$ was more consistent across sites during the daytime than at night. The slope and R² of the linear-fitting curve were −0.012 and 0.92, respectively. The R² of the polynomial fitting curve reached 0.98 (Fig. 2c).

$G/H$ also varied significantly in different half-hour periods of the intra-day and among the different sites (Fig. 2b). The variation in $G/H$ during each period was greater than that of $G/R_{\mathrm{n}}$. In other words, $G/H$ was less consistent across sites than $G/R_{\mathrm{n}}$. Like $G/R_{\mathrm{n}}$, the variation of $G/H$ among sites during the daytime was significantly lower than that at night. The variation range of $G/H$ 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 $G/H$ 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 R² of the polynomial fitting curve for the mean series was 0.95 (Fig. 2d).

Figure 3 Intra-day distribution of the linear fitted R² (a, c) and slope (b, d) between G and R_n (c and d are classified by land cover types).

3.2 Temporal and spatial analysis of the empirical relationship between G and R_n and between $G/R_{\mathrm{n}}$ and NDVI

The empirical relationship between G, R_n, and H was analyzed based on the measured data. Overall, G had a strong correlation with R_n but a relatively weak correlation with H. Figure 3 shows the R² and slope of the linear fitting between G and R_n in each half-hour period of the intra-day. In each period, G and R_n showed a strong linear correlation at most sites, with a fitted R² generally above 0.4. The mean and median R² values of all sites were mainly between 0.5 and 0.8. The strong correlation between G and R_n indicates that it is reasonable to use R_n to calculate G in the RS-based energy balance ET models. However, the correlation between G and R_n varied during the different periods. The correlation is relatively high in the periods around 06:00 and 18:00. Especially around 06:00, the R² of the linear fitting between G and R_n was greater than 0.7 at most sites, and the median R² of all sites reached 0.8. The correlation between G and R_n 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 R² of most sites was generally between 0.45 and 0.70, and the mean and median R² of all sites were mainly between 0.55 and 0.60. The R² of most sites generally ranged from 0.40 to 0.65 in the daytime periods, and the mean and median R² of all sites were concentrated around 0.50. The correlation between G and R_n was relatively low in the period from 10:00 to 15:00, and R² was lower than 0.65 at most sites, especially in the periods around 14:00, with the mean and median R² below 0.50.

The slope of the linear fitting of G and R_n 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 R_n 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 R² 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 R_n 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 R² lower than 0.5 during the daytime. In other-type sites, the correlation between G and R_n 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.

Figure 4 The linear- or exponential-fitting R² between daily series of G and R_n and that between monthly series of $G/R_{\mathrm{n}}$ and NDVI at each observed site (a is the median R² of 48 half-hour periods at each site; b and c are the R² values between G and R_n at 10:30 and 13:30, respectively; d and e are linear and exponential functions fitting R² between $G/R_{\mathrm{n}}$ and NDVI, respectively).

The empirical relationship between G and R_n not only varied significantly at different intra-day periods but also showed great spatial differences among the different sites. The linear fitted R² between the daily series of G and R_n at each site is shown in Fig. 4. As the median R² of 48 half-hour periods at each site (Fig. 4a), 91 % of the sites (210/230) showed an R² > 0.4. The linear-fitting R² between G and R_n was > 0.6 for 49 % of the sites (114/230). The mean R² for all sites was 0.58. This indicated that the G in most half-hour periods had a strong correlation with R_n at most observed sites. However, there were also 20 sites where the R² 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 R² of the northern sites being significantly higher than that of the southern sites. There was a strong correlation between G and R_n at most Eurasian, North American, and Oceanian sites, with a linear fitted R² generally exceeding 0.4. There was a relatively weak correlation at many African and South American sites, with an R² value of less than 0.4. At different latitudes, the strongest correlation between G and R_n was found at the middle- and high-latitude (> 45^∘) sites with the highest R² values. The R² 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 R² values of these sites were relatively low, with a mean value of 0.48.

The spatial distribution of the linear-fitting R² of G and R_n at 10:30 and 13:30 (Fig. 4b and c) was consistent with that in Fig. 4a, while the R² value was generally lower than the median of each half-hour period. At 10:30, the mean R² of all sites was 0.52. A total of 26 % of the sites (60/230) had R² values lower than 0.4, and 11 sites had R² values lower than 0.2. Although mainly distributed at low latitudes (< 30^∘), these sites were found on all the observed continents. The highest R² values were at high-latitude (> 60^∘) sites, with an average of 0.70. The R² 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 R_n at 13:30 was weaker than that at 10:30, with R² values slightly lower than those at 10:30. There was a relatively weak correlation between G and R_n at all African sites (R² < 0.4). Overall, the results showed a strong correlation between G and R_n at most observed sites. It is reasonable to simulate a daily series of G values from R_n 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 R_n was the strongest at other and wetland sites. The mean value of the median R² of statistical sites was 0.71 and 0.67 for these two types, respectively. There was also a strong correlation between G and R_n at cropland, shrubland, and grassland sites. The mean value of corresponding R² is about 0.60. There was only one site with a median R² lower than 0.4 in each of the three land cover types. The mean value of corresponding R² was 0.55 for forest sites. However, the correlation was relatively weak in the savanna-type sites, and the mean value of corresponding R² is 0.49. The weak correlations between G and R_n (R² < 0.4) were mainly distributed in the forest (6/97), grassland (7/42), and savanna (3/15) sites.

Figure 5 The RE, RMSE, and KGE simulated by the LC method at each site and at half-hour intervals.

The empirical correlation between $G/R_{\mathrm{n}}$ and NDVI was also analyzed. The results showed that the correlation between $G/R_{\mathrm{n}}$ and NDVI was weak in the daily series but strong in the monthly series. Figure 4d and e show the linear- and exponential-fitted R² values between the monthly series of $G/R_{\mathrm{n}}$ 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 $G/R_{\mathrm{n}}$ series and NDVI. The fitted R² values of the linear and exponential functions were consistent. The median and mean R² of all sites were 0.71 and 0.65, respectively, and 72 % of the sites (157/218) had an R² above 0.60. The exponential correlation between $G/R_{\mathrm{n}}$ and NDVI was stronger than the linear correlation at several sites, including US-Gle and US-Twt. The exponential-fitted R² of the two sites could be increased from 0.55 and 0.31 of the linear R² 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 R² of the two sites could be increased from 0.21 and 0.68 of the exponential R² to 0.43 and 0.83, respectively. Overall, the spatial distribution of the linear- or exponential-fitted R² of $G/R_{\mathrm{n}}$ and NDVI was similar to the linear fitted results of G and R_n. The correlation between $G/R_{\mathrm{n}}$ and NDVI was stronger in the Northern Hemisphere than in the Southern Hemisphere. The mean R² 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 R² values for these sites were generally higher than 0.6, with an average of 0.76. A relatively weak correlation between $G/R_{\mathrm{n}}$ and NDVI was found at the low-latitude sites, with a mean R² of 0.38. There were 13 and 15 sites showing weak linear and exponential correlations between $G/R_{\mathrm{n}}$ and NDVI, respectively, with a fitted R² 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.

Figure 6 The NSE simulated by the (a) LC, (b) LC_NDVI_P, (c) LC_NDVI_E, (d) LC_fc_SE, and (e) LC_fc_ST methods based on optimized parameters in each site and half-hour intervals.

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 R_n 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 The NSE simulated by the (a) LC, (b) LC_NDVI_P, (c) LC_NDVI_E, (d) LC_fc_SE, and (e) LC_fc_ST methods in each land cover type.

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.

Figure 8 Spatial distribution of NSE values simulated by the LC method (a is the median NSE of 48 half-hour values; b–e represent the NSE values at 06:30, 10:30, 13:30, and 18:00, respectively).

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.

Figure 9 The optimal parameters of the (a) LC, (b) LC_NDVI_P, and (c) LC_NDVI_E methods in each site and half-hour period.

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 Discussion

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 R_n, 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 $\mathrm{W}/\mathrm{m}$² and take up more than 15 % of R_n, there are some sites (20/230) with G values lower than 5 $\mathrm{W}/\mathrm{m}$². 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 $G/R_{\mathrm{n}}$ 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, R_n is generally calculated using radiation balance with RS images and meteorological inputs. However, observed R_n was used for simulating G in this study. In other words, it was assumed that R_n is accurately simulated by the RS–ET models. Therefore, it should be noted that the uncertainty in R_n 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 R_n, 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 R_n − 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 R_n in energy-balance-based RS–ET applications have been evaluated in this study. $G/R_{\mathrm{n}}$, which varied intra-day, was the key to these methods. Santanello and Friedl (2003) provided diurnal covariations in G and R_n. In this study, more concise linear and polynomial functions of the $G/R_{\mathrm{n}}$ daytime distribution were fitted (Fig. 2) from 230 observation sites worldwide. The linear- and polynomial-fitted R² 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 R_n, indicating that the accuracy of the G simulation is affected by the correlation between R_n and G. In other words, the stronger the correlation between R_n and G, the higher the accuracy of the simulation of G from R_n. 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 R_n. The sites with satisfied model performance were generally characterized by seasonal variations in G and R_n due to the climate in these regions. Conversely, the sites with poor model performance showed little seasonality. Whether G and R_n 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.

5 Conclusions

Intra-day and daytime G accounted for 19.2 % and 28.8 % of R_n, 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 R_n. In other words, the sites or periods with stronger correlations between G and R_n 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 R_n is the physical reason for the poor accuracy of G simulation in these regions and sites. Instead of the R_n, 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.