Simulation of Evapotranspiration at a 3-Minute Time Interval Based on Remote Sensing Data and SEBAL Model

Featured Application: Our research not only provides a method for estimating evapotranspiration, but also provides the possibility for additional remote sensing models to appear on a “minute” or even “second” time scale. Abstract: Using remote sensing to estimate evapotranspiration minute frequency is the basis for accurately calculating hourly and daily evapotranspiration from the regional scale. However, from the existing research, it is di ﬃ cult to use remote sensing data to estimate evapotranspiration minute frequency. This paper uses GF-4 and moderate-resolution imaging spectroradiometer (MODIS) data in conjunction with the Surface Energy Balance Algorithm for Land (SEBAL) model to estimate ET at a 3-min time interval in part of China and South Korea, and compares those simulation results with that from ﬁeld measured data. According to the spatial distribution of ET derived from GF-4 and MODIS, the texture of ET derived from GF-4 is more obvious than that of MODIS, and GF-4 is able to express the variability of the spatial distribution of ET. Meanwhile, according to the value of ET derived from both GF-4 and MODIS, results from these two satellites have signiﬁcant linear correlation, and ET derived from GF-4 is higher than that from MODIS. Since the temporal resolution of GF-4 is 3 min, the land surface ET at a 3-min time interval could be obtained by utilizing all available meteorological and remote sensing data, which avoids error associated with extrapolating instantaneously from a single image.


Introduction
Evapotranspiration (ET) is an important climatic factor besides solar radiation and atmospheric circulation, which controls the energy and mass exchange between the earth's ecosystem and the atmosphere, thus affecting the water balance of the ecosystem [1][2][3]. The accurate estimation of ET minute frequency at a regional scale is crucial to better understand detailed surface hydrological processes and provide more efficient catchment water management [4,5]. Traditionally ET is estimated for points or patches on the land surface but spatial heterogeneity in land surface characteristics precludes robust upscaling to the regional scale [6][7][8]. However, surface energy balance models using remote sensing data (specifically terrain, soil humidity, air and land surface temperature) enables accurate ET estimates at the regional scale [6,[9][10][11]. For example, in the Haihe River Basin of China, ground verification studies suggested the root mean square error associated with estimates of ET was approximately 0.32 mm/3 h [12]. Further, land-atmosphere interactions in California Delta of USA, 2 of 15 showed that the root mean square error was within 0.98 mm/day when calculating the one-day ET value [13].
Accurate estimation of evapotranspiration in irrigation areas is the basis for optimal management of water resources and water-saving irrigation. Remote sensing evapotranspiration model is of great significance for improving irrigation water use efficiency and simulating crop yield [14,15]. At present, remote sensing evapotranspiration model is mainly combined with meteorological data to evaluate irrigation water demand, irrigation efficiency and crop water monitoring in water resources management in irrigation areas [16][17][18]. In the above research, the input data with high temporal resolution played a very important role in the fine description of the above process. At present, the precipitation and temperature at minute frequency can already be obtained through the measured data of meteorological stations, and then the spatial distribution of temperature and precipitation in the irrigation area can be obtained through spatial interpolation, but the evapotranspiration is "cumulative data", which is different from "instantaneous data" such as temperature and precipitation, and the change at minute frequency is difficult to measure [19][20][21][22]. At present, most of the data measured by meteorological stations are the evaporation measured by the evaporation pan, not the evapotranspiration. Although lysimeter, scintillometer, flux tower and other methods can be used to measure the evapotranspiration with time interval of minutes, the cost of measuring instruments is huge and it is difficult to realize widespread arrangement like meteorological stations [21,23]. Therefore, using remote sensing data to simulate evapotranspiration at minute frequency has very significant potential applications.
However, limited by the current temporal resolution of satellites, the minimum time interval for simulating evapotranspiration using remote sensing data is hourly. Moreover, estimating ET with 1 h intervals is at the expense of expanding the spatial resolution of ET, for example (Liu et al., 2015), using FY meteorological satellite images simulated the evapotranspiration of Tibetan Plateau at 1 h interval, but the spatial resolution of ET was 1000 m, and the result could not accurately reflect the spatial diversity of ET [24]. Ideally, to obtain accurate ET at a regional scale, satellite data with both high temporal (i.e., hourly) and spatial (i.e., <100 m) resolution would be used. However, current satellite imagery either has a high temporal resolution and low spatial resolution (e.g., the time resolution of moderate-resolution imaging spectroradiometer is less than half a day, but its spatial resolution is 250 m to 1000 m.) or has a low temporal resolution and high spatial resolution (e.g., the spatial resolution of Landsat 8 is 15 m to 100 m, but its time resolution is 16 days). Consequently, it is almost impossible to use remote sensing data to obtain ET within 100 m spatial resolution with a temporal resolution of less than 1 h.
The high temporal and spatial resolution, and large width of the GF-4 data, offers significant potential to advance the understanding of ET on a regional scale, in addition to more accurate understanding of land surface energy balances and exchanges. However, due to the lack of thermal infrared band of GF-4, we cannot obtain the land surface temperature at the coincidence time with the visible bands of GF-4. In addition, we find that because there was no remote sensing data with time resolution of minutes, the existing evapotranspiration models do not consider using these remote sensing data, and the existing evapotranspiration models are not suitable for simulating ET with time interval of minutes. Researchers have been unable to take advantage of GF-4 application in simulating ET, and most of the papers were focused on weather monitoring and land target recognition [25,26]. Consequently, the aim of this research was to simulate ET based on GF-4, MODIS data and SEBAL (Surface Energy Balance Algorithm for Land) model, and compared the ET estimates derived from both MODIS and GF-4 imagery. The objectives are: (1) to characterize the variation at 3-min interval in ET estimated based on GF-4, MODIS and meteorological data across different land covers. (2) To compare ET estimates derived from MODIS and GF-4 captured at approximately the same time. (3) To verify MODIS and GF-4 simulated ET using field measured data.

Remote Sensing Data
The GF-4 satellite (launched on 29 December 2015) has a geosynchronous orbit fixed at 105.6 • E, an altitude of 36,000 km and a swath width of 400 km [27]. GF-4 can acquire remote sensing images in the extent of 2 • S-68 • N and 20 • E-180 • E. The GF-4 satellite is an important addition to the civil series of Chinese high-resolution earth observing satellites, as it provides both high spatial resolution (<100 m) and temporal (20 s) resolution visible, and near and medium wave infrared imagery with a pixel resolution of down to 50 m (Table 1). GF-4 images were obtained from (http://data.cma.cn/). Land surface temperature (LST) products and surface reflectance products of MODIS were obtained from (https://ladsweb.modaps.eosdis. nasa.gov). Digital elevation model (DEM) was obtained from (https://search.earthdata.nasa.gov/). Considering that weather should be sunny or cloudless at the imaging time, the imaging time should be close enough, and the imaging area of two images should be overlapped as much as possible. The imaging time of MODIS and GF-4 were shown in Table 2. LST and DEM data were resized or resampled, using the nearest neighbor method, to either 50 or 250 m to match the spatial resolution of the GF-4 and MODIS imagery, respectively. 1 Although the temporal resolution of GF-4 is 20 s, these data are only provided to specific institutions. The GF-4 data for the 3-min time interval used in this paper is open to all researchers, but the number of images and imaging time in each region will not be the same. Therefore, the number of GF-4 images in Table 2 is different between China and South Korea.

Meteorological and Flux Data
The time of sunrise, sunset was obtained from (http://www.sunrisesunset.com/). Daily maximum LST, maximum wind speed, mean wind speed, minimum air temperature, and maximum air temperature at 6 km resolution were obtained from the China Meteorological Administration (http://new-cdc.cma.gov.cn). Meteorological data were resized, using the nearest neighbor method, Appl. Sci. 2020, 10, 4919 4 of 15 to either 50 or 250 m to match the spatial resolution of the GF-4 and MODIS imagery, respectively. The flux data of GDK and GCK sites were provided by Korea Flux Monitoring Network (KoFlux), they were used to verify simulated ET, with a 30 min measuring interval. The position of GDK is 37 • 44 56 N, 127 • 08 57 E, and the position of GCK is 37 • 44 54 N, 127 • 09 44 E (Figure 1). More details about GDK and GCK sites can be obtained from Asia Flux website (http://asiaflux.net/). details about GDK and GCK sites can be obtained from Asia Flux website (http://asiaflux.net/).

Location and Study Area
The study area is located in China and South Korea (Figure 1). The study area ranges from an altitude of 100 m -1100 m, comprising a hilly area and regions of flat terrain. Climatologically, it belongs to the north temperate continental monsoon climate, and experiences distinct seasonal variation as well as impacts from the monsoon advancing and retreating. To minimize the influences of undulating terrain and hill shade, six study areas were selected: (CB) mainly cultivated land and built-up areas, (CG) mainly of garden and urban construction land, (CF) mainly cultivated land and forest land, (KB) mainly cultivated land and village construction land, (KE) mainly forest land and river wetland (GDK and GCK Flux sites were in KE region), and (KD) mainly cultivated land and forest land (Figure 1).

SEBAL Model
The SEBAL model is a physically based land surface energy balance model that uses remotely sensed input and has been widely used to calculate ET [28][29][30][31][32]. Input data to the SEBAL model include: ground elevation, visible light and near infrared bands from remote sensing images (used to calculate the Normalized Difference Vegetation Index, NDVI value), the thermal infrared band (used to calculate the Land Surface Temperature, LST), and temperature and wind speed at the imaging time. The two papers were referred for the sensible heat calculation, hot and cold pixel selection, net radiation and soil head flux calculation [33,34]. Then, we calculated the ET of MODIS and GF-4 at imaging time in Table 2 using SEBAL model.

Calculation of NDVI
The MODIS NDVI was calculated by [35]: in Equation (1), the NDVI denotes the normalized vegetation index, and and denote surface reflectance in the near infrared band and the red band, respectively.

Location and Study Area
The study area is located in China and South Korea ( Figure 1). The study area ranges from an altitude of 100 m-1100 m, comprising a hilly area and regions of flat terrain. Climatologically, it belongs to the north temperate continental monsoon climate, and experiences distinct seasonal variation as well as impacts from the monsoon advancing and retreating. To minimize the influences of undulating terrain and hill shade, six study areas were selected: (CB) mainly cultivated land and built-up areas, (CG) mainly of garden and urban construction land, (CF) mainly cultivated land and forest land, (KB) mainly cultivated land and village construction land, (KE) mainly forest land and river wetland (GDK and GCK Flux sites were in KE region), and (KD) mainly cultivated land and forest land ( Figure 1).

SEBAL Model
The SEBAL model is a physically based land surface energy balance model that uses remotely sensed input and has been widely used to calculate ET [28][29][30][31][32]. Input data to the SEBAL model include: ground elevation, visible light and near infrared bands from remote sensing images (used to calculate the Normalized Difference Vegetation Index, NDVI value), the thermal infrared band (used to calculate the Land Surface Temperature, LST), and temperature and wind speed at the imaging time. The two papers were referred for the sensible heat calculation, hot and cold pixel selection, net radiation and soil head flux calculation [33,34]. Then, we calculated the ET of MODIS and GF-4 at imaging time in Table 2 using SEBAL model.

Calculation of NDVI
The MODIS NDVI was calculated by [35]: in Equation (1), the NDVI denotes the normalized vegetation index, and ρ nir and ρ red denote surface reflectance in the near infrared band and the red band, respectively.

of 15
The GF-4 NDVI was calculated as for the MODIS NDVI, shown in Equation (1), but the apparent reflectance of the near infrared (the red band) was derived by [36]: in Equation (2), where ρ λ is the at-satellite reflectance, ESUN λ is the mean solar exoatmospheric irradiance in W·(m 2 ·ster·µm) −1 , θ s is the solar zenith angle in degrees, and d is the earth-sun distance in astronomical units. L λ is the at-satellite spectral radiance in W·(m 2 ·ster·µm) −1 . θ s can be obtained from the head file. The ESUN value for each GF-4 band has not been published yet, but can be calculated using the spectral response function of the solar spectrum curve and the sensor, as shown in Equation (3) [37].
In Equation (3), ESUN denotes the band average solar radiation outside of the atmosphere; λ 1 and λ 2 denote the upper and the lower integration limit of wavelength of band range; E(λ) denotes the solar spectrum radiation of the remote sensor out of the atmosphere at band λ, where different solar spectrum curves have different E(λ) values. Current studies show that the World Radiation Center (WRC) solar spectrum curve is the most favorable for the calculation of the ESUN using this sensor at a medium resolution in China [38]. According to the WRC solar spectrum curve, the E(λ) value could be in the range of λ 1 and λ 2 ; and S(λ) denotes the spectral response function of the remote sensor at band λ. The ESUN values for all GF-4 bands calculated using Equation (3) are shown in Table 3. The equation converting the DN value of satellite images to radiance images using the absolute calibration coefficient is [39]: where GF-4 provides gain values (A) under 5 statuses, for example, status 2-6-4-6-6 denotes an integration time for full color, blue, green, red, and near infrared of 2, 6, 4, 6, and 6 ms, respectively, the other statuses please see Table 4. Users need to determine the status of image bands according to the appropriate XML file parameters of GF-4 images, and then select the corresponding calibration coefficient. simulate LST at any time using following equations [40][41][42]. In this paper, two daily LST images after sunrise were as input from MODIS (Table 1).
Daytime LST variation after sunrise was derived by: in which: T(t) denotes the LST at the time t; T B,set denotes the LST at sunset; T B,max denotes the daily maximum LST; W 2 = π/(DL − 2p) denotes the angular frequency of the sinusoid in the second stage; DL denotes the daytime length; p = TIME x − NOON, TIME x denotes the time when the maximum LST appears, TIME x can be attained from meteorological stations, NOON denotes the time when the largest solar altitude appears, which is usually selected as 12.0; and ϕ 2 = π/2 − W 2 × TIME x denotes the phase angle.
This paper used the MODIS LST at the overpass time as inputs to calculate T B,max , T B,min , and T B,set , and then calculated the LST at GF-4 imaging time using Equations (5) and (6).

Calculation of Daytime Air Temperature
Since GF-4 do not take images at night, the general formula to calculate daytime air temperature is [42]: in which T a denotes daytime air temperature, and T min and T max denote daily minimum and daily maximum air temperature. In Equation (11), S(t) is a function of time (t) with data range of 0-1, represented as: in which t denotes any time, NOON denotes the time when the largest solar altitude appears, DL denotes the daytime length, and P denotes the time difference between the highest air temperature and the largest solar altitude. The time difference results from the intrinsic difference in heat storage of soil and air. When air temperature rises, T min denotes the lowest air temperature at present day in Equation (11); when air temperature decreases, T min denotes the lowest air temperature on the next day.

Calculation of Wind Speed
Daily variation in wind speed has the following features: the wind speed is low from nighttime to a time (t 1 ) in the morning, at which point the wind speed gradually increases to the maximum value until a time (t 2 ) in the afternoon, then gradually decreases to the minimum value until a time (t 3 ) at night. The t 1 , t 2 , and t 3 vary with location: in the study area, t 1 = 1.0, t 2 = 14.0, and t 3 = 0.0, and the variation of wind speed with time is expressed by the following equation [42]: in which: W a denotes wind speed at any time (m/s), W max and W min denote daily maximum and daily minimum wind speed (m/s), respectively, t a denotes any time, The daily minimum wind speed can be calculated by using the daily air movement distance and the daily maximum wind speed as input data: in Equation (10): TOT w denotes daily air movement distance (km/d), which can be calculated from the daily average wind speed (m/s) available from meteorological stations and the length of a day (24 h × 3600 s/h).

Verification and Evaluation of Simulation Results
We estimated ET from SEBAL model using MODIS and GF-4 at imaging time in Table 2 per minute by using the methods in Section 2.2. The following methods were used for verifying ET: (1) cross-validation ET obtained from different remote sensing data; (2) verify with ET data measured in field. In this paper, we first cross-verified ET of GF-4at 11:36 (T 8 ) using ET of MODIS at 11:35 (MT 0 ). Then ET at 13:50 (T 5 ) of GF-4 and 13:45 (MT 1 ) of MODIS were cross-verified. Finally, we used the field measured ET data to verify them from the Flux sites in South Korea, which was close to the imaging time of remote sensing images. Since the measuring time interval at KoFlux sites was 30 min, we chose the measured ET at 11:30 and 14:00 to verify the simulated ET at 11:22 (T 4 ) and 13:58 (T 7 ) of GF-4 respectively. Then, the field measure ET at 12:00 (MT 0 ) was used to verify the simulated ET of MODIS at 12:10 (MT 0 ).
In order to evaluate the difference between measured ET and remote sensing data, we use root mean square error (RMSE) and mean relative error (MRE) to analyze the error.
in Equations (11) and (12), Y i is the measured value, Y i is the simulated value obtained from the model, and n is the number of sample points [43,44]. The smaller the values of RMSE and MRE, the higher the simulation accuracy of the model [43,44].

Results
When the NDVI of GF-4 is taken as the input of SEBAL model, the spatial resolution of ET is 50 m, while when NDVI of MODIS is taken as the input of SEBAL model, the spatial resolution of ET is 250 m. in Equation (10): TOTw denotes daily air movement distance (km/d), which can be calculated from the daily average wind speed (m/s) available from meteorological stations and the length of a day (24 h×3600 s/h).

Verification and Evaluation of Simulation Results
We estimated ET from SEBAL model using MODIS and GF-4 at imaging time in Table 2 per minute by using the methods in Section 2.2. The following methods were used for verifying ET: (1) cross-validation ET obtained from different remote sensing data; (2) verify with ET data measured in field. In this paper, we first cross-verified ET of GF-4at 11:36 (T8) using ET of MODIS at 11:35 (MT0). Then ET at 13:50 (T ) of GF-4 and 13:45 (MT ) of MODIS were cross-verified. Finally, we used the field measured ET data to verify them from the Flux sites in South Korea, which was close to the imaging time of remote sensing images. Since the measuring time interval at KoFlux sites was 30 min, we chose the measured ET at 11:30 and 14:00 to verify the simulated ET at 11:22 (T ) and 13:58 (T ) of GF-4 respectively. Then, the field measure ET at 12:00 (MT ) was used to verify the simulated ET of MODIS at 12:10 (MT ).
In order to evaluate the difference between measured ET and remote sensing data, we use root mean square error (RMSE) and mean relative error (MRE) to analyze the error.
in Equations (11) and (12), Yi is the measured value, Yi' is the simulated value obtained from the model, and n is the number of sample points [43,44]. The smaller the values of RMSE and MRE, the higher the simulation accuracy of the model [43,44].

Results
When the NDVI of GF-4 is taken as the input of SEBAL model, the spatial resolution of ET is 50 m, while when NDVI of MODIS is taken as the input of SEBAL model, the spatial resolution of ET is 250 m.  As seen in Figures 2-4, since the spatial resolution of NDVI of GF-4 is higher than that of MODIS, the derived texture of ET is more obvious than that from MODIS. In order to cross-verify ET obtained from MODIS and GF-4, the bilinear method was used and the GF-4 ET data were resampled to 250 m, which was identical to the spatial resolution of MODIS. As shown in Figure 5, 420 total sampling points were randomly selected in areas CB, CG and CF (Figure 5a), and 870 total sampling points were randomly selected in areas KB, KE and KD (Figure 5b) to find variability between the two ETs. Meanwhile we verified the simulated ET with the measured ET data of KoFlux as shown in Figure 6.  As seen in Figures 2-4, since the spatial resolution of NDVI of GF-4 is higher than that of MODIS, the derived texture of ET is more obvious than that from MODIS. In order to cross-verify ET obtained from MODIS and GF-4, the bilinear method was used and the GF-4 ET data were resampled to 250 m, which was identical to the spatial resolution of MODIS. As shown in Figure 5, 420 total sampling points were randomly selected in areas CB, CG and CF (Figure 5a), and 870 total sampling points were randomly selected in areas KB, KE and KD (Figure 5b) to find variability between the two ETs. Meanwhile we verified the simulated ET with the measured ET data of KoFlux as shown in Figure 6.  As seen in Figures 2-4, since the spatial resolution of NDVI of GF-4 is higher than that of MODIS, the derived texture of ET is more obvious than that from MODIS. In order to cross-verify ET obtained from MODIS and GF-4, the bilinear method was used and the GF-4 ET data were resampled to 250 m, which was identical to the spatial resolution of MODIS. As shown in Figure 5, 420 total sampling points were randomly selected in areas CB, CG and CF (Figure 5a), and 870 total sampling points were randomly selected in areas KB, KE and KD (Figure 5b) to find variability between the two ETs. Meanwhile we verified the simulated ET with the measured ET data of KoFlux as shown in Figure 6. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 16  As can be seen from Figure 5, most ET simulation values of MODIS are less than that of GF-4's, but there is a strong correlation between the ET derived from GF-4 and MODIS, R = 0.581 (p <0.001) in areas CB, CG and CF, and R = 0.810 (p <0.001) in areas KB, KE and KD. While the difference in imaging time of MODIS and that of GF-4 is 1 to 5 min, it is certain that the ETs derived from both satellites have a significant linear correlation at the same imaging time. In addition, as can be seen from Figure 6, the difference between the simulated and measured values of ET was not obvious.
The following two methods were used to verify the simulated ET. Firstly, we used the field measured data of KoFlux to verify the simulated ET of MODIS and GF-4. Then, we took the ET calculated by MODIS as the measured value, and the ET calculated by GF-4 as the simulated value. The verification results were shown in Table 5.  As can be seen from Figure 5, most ET simulation values of MODIS are less than that of GF-4's, but there is a strong correlation between the ET derived from GF-4 and MODIS, R = 0.581 (p <0.001) in areas CB, CG and CF, and R = 0.810 (p <0.001) in areas KB, KE and KD. While the difference in imaging time of MODIS and that of GF-4 is 1 to 5 min, it is certain that the ETs derived from both satellites have a significant linear correlation at the same imaging time. In addition, as can be seen from Figure 6, the difference between the simulated and measured values of ET was not obvious.
The following two methods were used to verify the simulated ET. Firstly, we used the field measured data of KoFlux to verify the simulated ET of MODIS and GF-4. Then, we took the ET calculated by MODIS as the measured value, and the ET calculated by GF-4 as the simulated value. The verification results were shown in Table 5. As can be seen from Figure 5, most ET simulation values of MODIS are less than that of GF-4's, but there is a strong correlation between the ET derived from GF-4 and MODIS, R = 0.581 (p < 0.001) in areas CB, CG and CF, and R = 0.810 (p < 0.001) in areas KB, KE and KD. While the difference in imaging time of MODIS and that of GF-4 is 1 to 5 min, it is certain that the ETs derived from both satellites have a significant linear correlation at the same imaging time. In addition, as can be seen from Figure 6, the difference between the simulated and measured values of ET was not obvious.
The following two methods were used to verify the simulated ET. Firstly, we used the field measured data of KoFlux to verify the simulated ET of MODIS and GF-4. Then, we took the ET calculated by MODIS as the measured value, and the ET calculated by GF-4 as the simulated value. The verification results were shown in Table 5. As can be seen from Figure 6 and Table 5, RMSE = 1.04 (mm/day) and the MRE was 16.95% when validating simulated ET values using KoFlux data, which showed that it was feasible to simulate ET using the methods in this paper, while most GF-4 simulated ET values were higher than those of MODIS's, the MRE of ET when validating GF-4 using MODIS was less than 50%.
In order to further analyze the variation of ET in all study areas at time T 1 -T 15 in areas CB, CG and CF, T 1 − T 7 in areas KB, KE and KD, the minimum, maximum, and average values, and the standard error of ET at above time was calculated by ArcMAP 10.3 for all pixels, as shown in Figures 7 and 8.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 16 As can be seen from Figure 6 and Table 5, RMSE = 1.04 (mm/day) and the MRE was 16.95% when validating simulated ET values using KoFlux data, which showed that it was feasible to simulate ET using the methods in this paper, while most GF-4 simulated ET values were higher than those of MODIS's, the MRE of ET when validating GF-4 using MODIS was less than 50%.
In order to further analyze the variation of ET in all study areas at time T1-T15 in areas CB, CG and CF, T − T in areas KB, KE and KD, the minimum, maximum, and average values, and the standard error of ET at above time was calculated by ArcMAP 10.3 for all pixels, as shown in Figures  7 and 8.  As shown in Figure 7, the trends of the average, minimum, maximum and average ET in different study areas have significant differences, but the minimum ET in areas CF at time T1-T15 does not. The minimum ET does not vary in study area at time T1-T15, and the pixel with ET = 0 always exists. As shown in Figures 7 and 8, the maximum value trend in all areas is the same as that of the average values, while the fluctuation in trend of the maximum and minimum value in all areas is different. This is mainly due to the difference among soil types and meteorological conditions in the six areas, causing the impacts on ET to vary [3,45]. This validates that due to the different of surface types and meteorological conditions, even two MODIS remote sensing images show daily variation in remote sensing pixel ETs. This daily variation does not follow the constant Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 16 As can be seen from Figure 6 and Table 5, RMSE = 1.04 (mm/day) and the MRE was 16.95% when validating simulated ET values using KoFlux data, which showed that it was feasible to simulate ET using the methods in this paper, while most GF-4 simulated ET values were higher than those of MODIS's, the MRE of ET when validating GF-4 using MODIS was less than 50%.
In order to further analyze the variation of ET in all study areas at time T1-T15 in areas CB, CG and CF, T − T in areas KB, KE and KD, the minimum, maximum, and average values, and the standard error of ET at above time was calculated by ArcMAP 10.3 for all pixels, as shown in Figures  7 and 8.  As shown in Figure 7, the trends of the average, minimum, maximum and average ET in different study areas have significant differences, but the minimum ET in areas CF at time T1-T15 does not. The minimum ET does not vary in study area at time T1-T15, and the pixel with ET = 0 always exists. As shown in Figures 7 and 8, the maximum value trend in all areas is the same as that of the average values, while the fluctuation in trend of the maximum and minimum value in all areas is different. This is mainly due to the difference among soil types and meteorological conditions in the six areas, causing the impacts on ET to vary [3,45]. This validates that due to the different of surface types and meteorological conditions, even two MODIS remote sensing images show daily variation in remote sensing pixel ETs. This daily variation does not follow the constant As shown in Figure 7, the trends of the average, minimum, maximum and average ET in different study areas have significant differences, but the minimum ET in areas CF at time T 1 -T 15 does not. The minimum ET does not vary in study area at time T 1 -T 15 , and the pixel with ET = 0 always exists. As shown in Figures 7 and 8, the maximum value trend in all areas is the same as that of the average values, while the fluctuation in trend of the maximum and minimum value in all areas is different. This is mainly due to the difference among soil types and meteorological conditions in the six areas, causing the impacts on ET to vary [3,45]. This validates that due to the different of surface types and meteorological conditions, even two MODIS remote sensing images show daily variation in remote sensing pixel ETs. This daily variation does not follow the constant linear rule, meaning that extrapolating instantaneous ET at the imaging time to a hourly or daily time scale will cause relatively large errors [46,47]. Although GF-4 has no thermal infrared band and the LST at imaging time could not be obtained, its temporal resolution is relatively high and the land surface ET at a three-minute time interval can be attained by using all available meteorological and MODIS data. By using this method, the error associated with extrapolating instantaneous ET from one remote sensing image could be avoided and the real spatial diversity of ET at various imaging time could be obtained.

Discussion
While this paper compared ETs derived from GF-4 and MODIS using the SEBAL model, the extent of the differences from SEBAL model ET input data needs further analysis. The reason for the difference in ET simulation results between the two sensors needs further discussion, and the validation of ET simulation results requires careful consideration as well: (1) Errors in SEBAL model input data lead to a decrease in ET simulation accuracy. However, in many areas, simulation is limited by available data, and often, the wind speed, LST, air temperature, and other data cannot be attained. While simulating the above parameters using the present model, various kinds of error must exist. If the above meteorological data could be replaced by observed data, the model accuracy could be significantly improved. In addition, the ESUN of GF-4 is not published yet, and the subsequent ESUN simulation that was used has errors which affected both the calculation of NDVI and the simulated values of ET.
(2) Besides the errors associated with SEBAL input data, the dependence on remote sensing data can lead to different ETs derived from the GF-4 and MODIS sensors. Many researchers have pointed out that the local variables of the SEBAL model would change with the scope of the remote sensing image, which could lead to additional uncertainty associated with the domain dependence of the remote sensing model accuracy [7,28,34,48]. On the other hand, remote sensing images with different spatial resolutions would also cause differences in local variables, and thus could lead to completely different dry-wet conditions, dependent upon the resolution of the remote sensing models. For all established models based on remote sensing images, the model outputs may have uncertainties if the determination of variables or input parameters depends on the scope and resolution of those images.
(3) The objective of this paper was not to validate whether the GF-4 or MODIS remote sensor has higher accuracy in ET retrieval, but to quantify the similarities and differences between ETs derived from both satellites by comparing their retrievals. It is clear that more efficient ground observations are essential to validate model results. Currently, widely used ET remote sensing model validation techniques include the lysimeter method, the field water balance method, the Bowen ratio method, the eddy covariance method and the large aperture scintillometer method [23,[49][50][51][52]. However, in practical operation, the spatio-temporal scale of the observed data does not always match with that of the model retrievals. Therefore, we can only use KoFlux data with 2-10 minutes difference in remote sensing images for verification. The GF-4 is only the beginning of high orbit and high temporal-spatial resolution satellite technology, and as more research and data are gathered, new observational approaches will emerge.
(4)It is also necessary to consider cross-validation using remote sensing data with higher spatial resolution. We tried to verify GF-4 imaged at 11:14 in regions CB, CG and CF by using Landsat 8 imaged at 10:29. Meanwhile in regions KB, KE and KD, Landsat 8 imaged at 10:04 was used to verify GF-4 imaged at 11:14, the verification results of ET were shown in Table 6. As can be seen from Table 6, when the simulated ET of GF-4 was verified by Landsat 8, the MRE of ET of GF-4 was very large. The imaging time of Landsat 8 was 44 min earlier than that of GF-4 in regions CB, CG and CF, and the imaging time of Landsat 8 was 70 min earlier than that of GF-4 in region KB, KE and KD. With such a long time difference, we think that the results of Table 6 can not accurately represent the actual accuracy of GF-4. However there is still a strong correlation between the ET derived from GF-4 and Landsat 8, R = 0.548 (p < 0.001).

Conclusions
(1) The spatial distribution of ET derived from MODIS and GF-4 showed that the texture of ET derived from GF-4 was more obvious than that from MODIS, and GF-4 was able to express the variability of ET spatial distribution. The correlation between ETs derived from two sensors showed significant linear correlation; and the ET values derived from GF-4 were higher than that of MODIS.
(2) Even in the same study area, trends in the maximum value, the average value, and the standard error of ET at different times were not the same. This validates that even if in two MODIS remote sensing images scope, due to different surface types and meteorological conditions, the daily variation of ET in the remote sensing pixel scale does not follow the constant linear rule.
(3) Since the temporal resolution of GF-4 is 3 minutes in our paper, the land surface ET every three-minutes could be obtained by using all available meteorological data and other remote sensing data. At this scale, the extrapolation of instantaneous ET from a single image was avoided, and the spatial diversity at each imaging time was attained.
(4) This paper has validated ET derived from GF-4 and MODIS, and the verification results also showed that the error was within the normal range, but more observational methods to validate ET estimation results of GF-4 are needed, which should be a major goal for future studies.

Perspectives
The GF-4 satellite can not only collect images with a large scope and at high temporal and spatial resolutions, but can also carry out "staring observation" in a specific region according to user instructions, which will play an important role in disaster reduction, meteorology, earthquake, forestry, and obtaining precise measurements in other fields. For example, GF-4 could practically observe a typhoon, since middle and low orbit satellites have high spatial resolution, but the retrogression time is relatively long, and the continuous observation of interesting areas is difficult. Compared with present high orbit continuous observation meteorological satellites, GF-4 could not only continuously provide information on the development of a typhoon, but could also take observations and measurements of details like typhoon textures.
In addition, GF-4 will further accelerate the development of advanced models supported by remote sensing data. Currently, models in disparate fields such as ecology, hydrology, geochemistry, and meteorology are driven by remote sensing data, but are limited by the temporal and spatial resolution of available remote sensing data, and face issues of fulfilling simultaneous goals such as having "short time", "high accuracy", and "large space" [53][54][55]. Most models driven by remote sensing are on time scales of days, months, or years. Even if a part of these models could realize hourly time scales from meteorological satellites such as those available from NOAA, the spatial resolution of their simulation results is over 1 km, which could not accurately estimate the spatial diversity of observation pixels [56][57][58]. The comprehensive advantages of GF-4 of increased temporal and spatial resolution and large widths would facilitate the appearance of various kinds of remote sensing models on fine time scales at minutes or even seconds.