Incorporation of Net Radiation Model Considering Complex Terrain in Evapotranspiration Determination with Sentinel-2 Data

Evapotranspiration (ET) is the primary mechanism of water transformation between the land surface and atmosphere. Accurate ET estimation given complex terrain conditions is essential to guide water resource management in mountainous areas. This study is based on the ETWatch model driven by Sentinel-2 remote sensing data at a spatial resolution of 10 m incorporating a net radiation model considering the impact of a complex terrain. We tested our model with two years of data in two regions with a high relief near the Huairou (2020) and Baotianman (2019) weather stations. Regarding the validation results of the ET model, the coefficient of determination (R2) reached 0.84 in Huairou and 0.86 in Baotianman, while the root mean square error (RMSE) value reached 0.59 mm in Baotianman and 0.82 mm in Huairou. The validation results indicated that the model is applicable in regions with a complex terrain, and the ET results can capture topographic textures. In terms of the slope aspect, the ET value on south-facing slopes is higher than that on north-facing slopes in both study areas. Accurate ET monitoring in mountainous regions with a high relief yields a profound meaning in obtaining a better understanding of the characteristics of heat and water fluxes at different vegetation growth stages and underlying surface types, which can provide constructive suggestions for water management in mountainous areas.


Introduction
Evapotranspiration (ET), mainly comprising plant transpiration and soil evaporation, is the essential pathway of water transformation between the land surface and atmosphere, linking changes in surface water, carbon cycling, and surface energy [1][2][3]. Factors such as meteorological, vegetation, and radiation conditions and soil moisture influence the ET process [4,5]. In mountainous areas with a complex topography, topographic factors require greater consideration to obtain better ET estimates [6,7].
Regarding meteorological factors, air pressure and temperature and water vapor pressure tend to decrease with increasing elevation, wind speed, and precipitation [8][9][10]. In addition, radiation conditions can vary considerably from slope to slope [11,12]. In general, slopes facing south in the Northern Hemisphere are considered sunny slopes, while slopes facing north are considered shady slopes [13]. Sunny slopes tend to receive more solar meaning of the ET process remains ambiguous. Energy balance residual methods are often limited by the low spatial resolution of the thermal infrared band [57]. Methods based on PM equations are mechanistic approaches but only require meteorological data, surface net radiation, and surface resistance to water vapor transmission as input data. With the continuous development of optical remote sensing, high-spatial resolution satellites, such as Sentinel-2, can facilitate ground observations with a spatial resolution of 10 m, which provides a fine picture of underlying surface conditions. Moreover, the revisit period of Sentinel-2 is shorter than that of Landsat series satellites. With a shorter revisit period, the key parameters in ET calculation, such as the normalized difference vegetation index (NDVI) and surface albedo, which exhibit a limited variability during a short period [58,59], can be extended to the daily scale via data interpolation in day-by-day ET simulations.
Accurate estimation of the spatial distribution of net radiation is essential for ET calculation [60,61]. The topography of the underlying mountainous surface is complex, and the solar incidence conditions vary with the slope gradient. Moreover, the surrounding topography imposes a shading effect on the slope surface, resulting in considerable differences in the received incident radiation between various slope gradients and aspect directions. Therefore, it is necessary to consider the influence of terrain factors on the incident solar radiation in net radiation calculations to obtain better ET estimates. The availability of a high-resolution digital elevation model (DEM), e.g., TanDEM-X, which exhibits a spatial resolution of 0.4 arcsec (approximately 12 m), allows the surface relief characteristics of mountainous areas to be reflected in greater detail, thus facilitating the acquisition of high-spatial resolution surface radiation data. To our knowledge, no studies have incorporated a high-resolution DEM in remote sensing-based ET calculations.
Considering the above reasons, this study aims to estimate ET with the ETWatch model at a high spatial resolution (10 m) based on Sentinel-2 and TanDEM-X data to help decisionmakers and land and water managers develop suitable water resource management strategies in areas with a complex topography. This work follows two specific objectives: (1) to develop a daily net radiation-terrain-based model and (2) to estimate ET, considering remote sensing data, in different agroecosystems in areas with a complex terrain. The key findings could be helpful for land and water managers in similar areas.

Study Area
This study was performed in two regions covering the Huairou and Baotianman stations. Huairou is located in the southern part of the Yan Mountains in northeastern Beijing in northern China, which experiences a warm temperate semihumid continental monsoon climate with four distinct seasons, with both rain and heat during the same period. Consequently, warm and humid conditions prevail in summer, while cold and dry conditions occur in winter. The average annual temperature ranges from 9 • C to 13 • C, and the average annual precipitation ranges from 600 to 700 mm, mainly concentrated in the period from June to August. The Huairou station is located in the southeastern part of Huairou, and the longitudinal and latitudinal coordinates of the site are 116 • 39 35 E and 40 • 25 22 N, respectively, at an altitude of 328 m. Baotianman is located in the eastern part of the Qinling Mountains, on the southern slope of the Funiu Mountains in Neixiang County, Henan Province, Central China, which belongs to the transitional area from the northern subtropical zone to the warm temperate zone, and exhibits the transitional characteristics of the eastern monsoon zone, with a monsoonal continental climate and four distinct seasons. Summer is hot, winter is cold, the temperature quickly rises in spring, the annual average temperature reaches 15.1 • C, and the annual average precipitation reaches 855.6 mm. Baotianman station is located within the Baotianman Nature Reserve, with geographical coordinates of 111 • 56 07 E and 33 • 29 59 N and an altitude of 1410.7 m. The first selected study area encompassed the Huairou station, covering approximately 400 km 2 , and the second selected study area encompassed the Baotianman station, covering with geographical coordinates of 111°56′07″E and 33°29′59″N and an altitude of 1410.7 m. The first selected study area encompassed the Huairou station, covering approximately 400 km 2 , and the second selected study area encompassed the Baotianman station, covering approximately 225 km 2 (Table 1 and Figure 1). We denoted these two study areas, i.e., Huairou and Baotianman, as HR and BTM, respectively.

Data Sources
The data considered in this study mainly include remote sensing data, in-situ EC measurement data, meteorological data, elevation data and certain auxiliary data, like landcover, as described in the following subsections. Table 2 shows the summary of the dataset used in this study.

Data Sources
The data considered in this study mainly include remote sensing data, in-situ EC measurement data, meteorological data, elevation data and certain auxiliary data, like landcover, as described in the following subsections. Table 2 shows the summary of the dataset used in this study. The remote sensing data utilized in this study largely include Sentinel-2, Moderate Resolution Imaging Spectroradiometer (MODIS), Fengyun-2F (FY-2F), Atmospheric Infrared Sounder (AIRS), National Centers for Environmental Prediction (NCEP) and Soil Moisture Active Passive (SMAP) data. Sentinel-2 comprises two multispectral imaging satellites, 2A and 2B, launched successively, carrying Multispectral Imager (MSI) sensors with an orbital altitude of 786 km and an amplitude of 290 km. These satellite sensors cover 13 spectral bands with a ground resolution from 10 to 60 m in each data band, ranging from visible and near-infrared to shortwave infrared (SWIR), with different spatial resolutions. The revisit period is 10 days for one satellite and 5 days for two satellites, considered together. In this study, Sentinel-2 1C-level data were downloaded from the European Space Agency (ESA) Copernicus Data Centre website (https://scihub.copernicus.eu/dhus/#/home (accessed on 2 December 2021)). Atmospheric, radiometric, and geometric corrections of the downloaded 1C-level data were performed with the Sen2Cor plug-in in Sentinel Application Platform (SNAP) software with terrain correction and BRDF correction options. The images were mosaicked and clipped to obtain the surface reflectance in the two study areas. Surface reflectance data in the near-infrared and red bands were obtained to calculate the NDVI, which was considered to simulate the vegetation cover and leaf area index (LAI) and subsequently applied in canopy conductance and soil evaporation calculations. We chose Sentinel-2 rather than Satellite products such as AVHRR, MODIS and Himawari because these satellite products cannot meet the needs of ET estimation in regions with complex terrain due to relatively low spatial resolution. As for Landsat, the temporal resolution is not satisfactory, which cannot provide enough monitoring during the main growing season. The surface albedo was calculated through multiband fitting and applied in subsequent net radiation calculations. Combined with a cloud mask, the Savitzky-Golay (S-G) filtering method was chosen to temporally extend the NDVI and albedo on cloud-free days to daily scales, which has been widely applied in relevant studies [62,63]. Cloud classification data based on FY-2F geostationary satellite data acquired from the China National Satellite Meteorological Center (http://satellite.nsmc.org.cn/PortalSite/Default.aspx (accessed on 5 December 2021)) were employed to estimate sunshine hours [64]. MODIS land surface temperature data products (MOD11A1) at a 1-km spatial resolution were collected from the Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 5 December 2021)). Observations provided by the AIRS installed on the Aqua satellite of the Goddard Earth Sciences Data and Information Services Center (GES DISC) (https://disc.gsfc.nasa.gov/ (accessed on 7 December 2021)) provided vertical distribution data of the air temperature and humidity, and the height of the atmospheric boundary layer was determined according to a previous study [65]. The air temperature and humidity at the height of the atmospheric boundary layer were retrieved from AIRS data, while the wind velocity was determined with NCEP reanalysis data obtained from the Physical Sciences Laboratory of the National Oceanic and Atmospheric Administration (NOAA) (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html (accessed on 5 December 2021)) [66]. National Aeronautics and Space Administration (NASA)-US Department of Agriculture (USDA) enhanced SMAP global soil moisture data provided global information at a 10-km spatial resolution and were mainly employed for daily surface conductance reconstruction [67].

In-Situ Tower Observation Data
The in-situ observations obtained at flux towers mainly included radiation component observation and EC data. The radiation component observation data were considered to validate the net radiation-terrain-based model, while the EC data were employed to validate the resulting ET model estimates. In this study, we collected data in HR in 2020 and BTM in 2019. Observation data pertaining to HR were collected from the field observatory, while BTM data were provided by ChinaFLUX (http://www.chinaflux.org/ (accessed on 10 December 2021)). The adopted EC observation instruments largely comprised ultrasonic anemometers and CO 2 /H 2 O infrared gas analyzers. The observation tower in HR is 40 m high, while the observation tower in BTM is 36 m high. The EC instruments were set up at a height of 30 m in HR and 29 m in BTM, while the radiation component observation instruments were placed at 20 m in HR and 22 m in BTM. The data sampling frequency of the EC observations at the two sites was 10 Hz, which were averaged and stored for 30 min. The underlying surface surrounding the HR observation tower is dominated by arborvitae, while oak trees dominate the underlying surface surrounding the BTM observation tower. Since high-resolution ET values are calculated in this study, it was necessary to average the model calculation results based on the footprint area of each EC instrument to match the flux observation results during model validation. This study applied a flux footprint prediction (FFP) model to calculate the footprint area [68], and the results are shown in Figure 2. The collected EC data were processed following the literature [44], including the exclusion of data outliers, data before and after precipitation, and nighttime data under extremely low-turbulence conditions (the friction velocity is lower than 0.2 m/s). It should be noted that approximately 88 days of EC data from May to August 2020 were missing in HR because of instrument or data logger issues.
The in-situ observations obtained at flux towers mainly included radiation component observation and EC data. The radiation component observation data were considered to validate the net radiation-terrain-based model, while the EC data were employed to validate the resulting ET model estimates. In this study, we collected data in HR in 2020 and BTM in 2019. Observation data pertaining to HR were collected from the field observatory, while BTM data were provided by ChinaFLUX (http://www.chinaflux.org/ (accessed on 10 December 2021)). The adopted EC observation instruments largely comprised ultrasonic anemometers and CO2/H2O infrared gas analyzers. The observation tower in HR is 40 m high, while the observation tower in BTM is 36 m high. The EC instruments were set up at a height of 30 m in HR and 29 m in BTM, while the radiation component observation instruments were placed at 20 m in HR and 22 m in BTM. The data sampling frequency of the EC observations at the two sites was 10 Hz, which were averaged and stored for 30 min. The underlying surface surrounding the HR observation tower is dominated by arborvitae, while oak trees dominate the underlying surface surrounding the BTM observation tower. Since high-resolution ET values are calculated in this study, it was necessary to average the model calculation results based on the footprint area of each EC instrument to match the flux observation results during model validation. This study applied a flux footprint prediction (FFP) model to calculate the footprint area [68], and the results are shown in Figure 2. The collected EC data were processed following the literature [44], including the exclusion of data outliers, data before and after precipitation, and nighttime data under extremely low-turbulence conditions (the friction velocity is lower than 0.2 m/s). It should be noted that approximately 88 days of EC data from May to August 2020 were missing in HR because of instrument or data logger issues.

Meteorological Data
The first part of the meteorological data encompassed daily meteorological data obtained from the China Meteorological Data Service Center (https://data.cma.cn/ (accessed on 28 November 2021)). The data at each station mainly include the relative air humidity (RH), wind speed (Vwind), atmospheric pressure (PRS), maximum air temperature (Tmax), minimum air temperature (Tmin), mean air temperature (Tmean), and sunshine hours (sunt). These parameters were extrapolated from the point scale to the spatial scale with the kriging interpolation method [69]. Sea-level values of four parameters, namely, Tmax, Tmin, Tmean, and PRS, were first calculated based on the site elevation using empirical relationships between these parameters and elevation [51,70]. After interpolation, the extrapo-

Meteorological Data
The first part of the meteorological data encompassed daily meteorological data obtained from the China Meteorological Data Service Center (https://data.cma.cn/ (accessed on 28 November 2021)). The data at each station mainly include the relative air humidity (RH), wind speed (V wind ), atmospheric pressure (PRS), maximum air temperature (T max ), minimum air temperature (T min ), mean air temperature (T mean ), and sunshine hours (sunt). These parameters were extrapolated from the point scale to the spatial scale with the kriging interpolation method [69]. Sea-level values of four parameters, namely, T max , T min , T mean , and PRS, were first calculated based on the site elevation using empirical relationships between these parameters and elevation [51,70]. After interpolation, the extrapolated sea-level values in all pixels were transformed into values at the corresponding elevation with DEM data. All the interpolated meteorological data were processed to follow the same spatial reference system and pixel size as the remote sensing data. The second part of the meteorological data originated from China Meteorological Radiation Data International Exchange Stations, also obtained from the China Meteorological Data Service Center (https://data.cma.cn/ (accessed on 28 November 2021)). These stations provided daily values of radiation observations, mainly including the total radiation, direct beam radiation, reflected radiation, and diffuse radiation. Observation data were collected at two stations (nearest to HR and BTM) from 2000 to 2020 to calibrate the solar radiation model in the net radiation calculation process.

Elevation Data
The elevation data, namely, TanDEM-X global DEM data, was mainly applied in the net radiation calculation process and provided by the German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt or DLR) (https://tandemx-science.dlr.de/ (accessed on 15 November 2021)) with a spatial resolution of 0.4 arcsec (~12 m) [71].

Auxiliary Data
Auxiliary data, including land cover, was obtained in this study. ChinaCover, a land cover dataset for China, was employed in this study and was mostly involved in the analysis of ET model results and parameterization processes such as canopy conductance, which was provided by the Aerospace Information Research Institute, Chinese Academy of Sciences (AIRCAS). ChinaCover contains 38 secondary classes and 6 primary classes (forestland, grassland, cropland, water bodies, built-up land, and bare land) with a spatial resolution of 30 m [72]. Figure 3 shows a summary of the input data in this study and provides the proposed slope-scale ET model framework.

Net Radiation Calculation
In this study, the net radiation (R n ) was calculated with Equation (1) [73].
where R s↓ is the downward shortwave radiation, α is the surface albedo, and R nl is the net longwave radiation. The downward shortwave radiation (R s↓ , incident solar radiation) was calculated as the sum of the direct solar radiation (R s↓_dir ), sky diffuse radiation (R s↓_di f ) and reflected radiation in adjacent regions (R s↓_adj ). Hence, the total incident solar radiation R s↓ at the surface, considering a complex terrain, can be expressed as Equation (2). The direct solar radiation (R s↓_dir ) was calculated with Equation (3a-e) [74,75].
where R s↓_dir0 (W/m 2 ) is the direct solar radiation on a horizontal surface, R b denotes the ratio of the direct solar radiation on an inclined surface to that on a horizontal surface, α day (radians) is the daily average solar elevation angle calculated based on the latitude (ϕ, radians) and day of the year (DOY) [75], θ z (radians) is the daily solar zenith angle complementary to α day , φ s (radians) is the daily solar azimuth angle, δ is the declination angle (radians) and can be calculated with the DOY [73], β (radians) is the slope, and A (radians) is the slope orientation. R s↓_dir0 can be estimated with the daily sunshine hours sunt and daily maximum possible sunshine hours sunt max with validated linear, quadratic, or cubic empirical regression models of R s↓_dir0 /R a and sunt/sunt max , expressed as Equation (4a-c) [76,77].
where R a is the daily extraterrestrial solar radiation (MJ·m −2 ·day −1 ) calculated according to [73], and the units can be converted into W/m 2 through multiplication by 10 6 24×3600 . Moreover, sunt max can be determined based on δ and ϕ, and lowercase letter with numeric subscripts such as a 1 , b 2 , and c 3 are empirical regression coefficients.
The sky diffuse radiation (R s↓_di f ) can be calculated by employing the sky view factor (Φ sky ) to correct the diffuse sky radiation on a horizontal surface (R s↓_di f 0 ), expressed as Equation (5a,b) [78,79].
where V sky is the sky view angle (radians). In a single pixel, the minimum sky view angle in eight directions (namely, N, NE, N, NW, W, SW, S, and SE) was calculated, and the average minimum sky view angle was regarded as the sky view angle in this pixel [80]. In each image pixel (P 0 ), the elevation angle from every pixel P within radius of L (L here was set to 10 [81]) along the above eight directions to pixel P 0 (θ) was calculated with Equation (6a,b). where X and Y are the image coordinates of pixel P, X 0 and Y 0 are the image coordinates of pixel P 0 , H and H 0 are the altitudes (m) of pixel P and P 0 , respectively, and S p is the pixel size (m). Thus, the corresponding zenith angle (θ 0 ) can be expressed as Equation (7).
The average minimum zenith angle in eight directions is defined as the sky view angle V sky . R s↓_di f 0 can also be estimated in a similar manner to the daily sunshine hours sunt and daily maximum possible sunshine hours sunt max with validated linear, quadratic or cubic empirical regression models [76,77].
The reflected radiation in adjacent regions (R s↓_adj ) was calculated with Equation (8a,b) and the approximation method of Dozier and Frew (1990) [82].
The net longwave radiation (R nl ) was calculated mainly based on the meteorological parameters, as expressed in Equation (9) [73].
where σ is the Stefan-Boltzmann constant (4.903 × 10 −9 MJ·K −4 ·m −2 ·day −1 ), e a is the actual vapor pressure (KPa) and can be calculated with the relative humidity and air temperature [83], T max and T min are the daily maximum and minimum air temperatures (K), respectively, and R s↓0 is the clear-sky solar radiation (MJ·m −2 ·day −1 ) and can be calculated with R a and altitude [73]. The surface shortwave broadband albedo (α), determining the incident solar radiation amount absorbed by the underlying surface, was calculated with the linear combination of the Sentinel-2 surface reflectance data using the method reported by Li et al. (2018) [84], thus realizing narrow-to-broadband conversion via Equation (10).
where i is the band number of the Sentinel-2 images, x i is the surface reflectance in band i, and a i and c are empirical regression coefficients, as listed in Table 3.

Evapotranspiration Calculation
In this study, the net radiation model considering the influence of terrain, as introduced in Section 2.3.1, was incorporated into the ETWatch model to estimate the ET value at the slope scale in the two study areas.
It is difficult to perform ET calculations on cloudy days. The ETWatch model mainly captures key information on sunny days and combines multisource remote sensing and meteorological data to calculate energy balance components. First, based on the concept of the residual method, the instantaneous energy balance components (net radiation [85], soil heat flux [86], and sensible heat flux [87]) were calculated to determine the instantaneous ET. Second, following the assumption that the evaporation fraction remains constant during the day, the instantaneous ET value on sunny days was extended to the daily ET value on a certain time scale [88]. Then, the surface conductance on sunny days was obtained with an inverted PM equation. Next, the surface conductance on sunny days was extended to the daily scale, combined with meteorological, biophysical, soil and radiation factors [89]. Finally, the daily ET was calculated through the PM equation. To conclude, the ETWatch model integrates a series of submodels, namely, FY sunshine hour (FYsunt) model [90], atmospheric boundary layer (ABL) model [65], aerodynamic roughness length (A z0m ) model [91], vapor pressure deficit (vpd) model, net radiation (R n ) model [85], sensible heat flux (H) model [87], soil heat flux (G) model [86], and gap-filling model of the surface conductance (r s ) [89], which jointly realize the fine description of the temporal and spatial distribution patterns of the key parameters in the ET estimation process, thus facilitating improvement of the regional ET estimation reliability [69]. Figure 4 shows the ETWatch model architecture.

Net Radiation Results
In the net radiation calculation process, this study focused on the influence of terrain factors on solar radiation and divided solar radiation into three parts: direct solar radiation, sky diffuse radiation, and reflected radiation in adjacent regions. Figure 5 shows the spatial distribution of the slope and sky view factor in the two study areas. The sky view factor is also an important indicator reflecting the characteristics of the terrain, and the value ranges from 0 to 1. A sky view factor value close to 1 suggests a good view of the sky in the surrounding hemisphere, which generally occurs in plain areas, ridges, or mountain peaks, while a sky view factor value close to 0 often occurs in low-lying or valley bottom areas. Concerning the situation surrounding HR, the southeast side of the study area is a plain area, and the main land cover types include farmlands and built-up surfaces. When the terrain is relatively flat, the sky view factor is high, while the sky view factor is low under complex mountainous terrain conditions. In terms of the situation surrounding BTM, most of the study area is located in mountainous areas, and only some narrow topographically flat areas exhibit high sky view factor values.
Remote Sens. 2022, 14, x 11 of 29 spatial distribution of the slope and sky view factor in the two study areas. The sky view factor is also an important indicator reflecting the characteristics of the terrain, and the value ranges from 0 to 1. A sky view factor value close to 1 suggests a good view of the sky in the surrounding hemisphere, which generally occurs in plain areas, ridges, or mountain peaks, while a sky view factor value close to 0 often occurs in low-lying or valley bottom areas. Concerning the situation surrounding HR, the southeast side of the study area is a plain area, and the main land cover types include farmlands and built-up surfaces. When the terrain is relatively flat, the sky view factor is high, while the sky view factor is low under complex mountainous terrain conditions. In terms of the situation surrounding BTM, most of the study area is located in mountainous areas, and only some narrow topographically flat areas exhibit high sky view factor values. Figure 5. Spatial distribution of the slope (a1,b1) and sky view factor (a2,b2) in HR (a1,a2) and BTM (b1,b2). Figure 6 shows scatter plots of the sky view factor versus the elevation and slope in the two study areas. The elevation in HR is clearly concentrated within the range from 50-1100 m, and the sky view factor ranges from 0.7-0.875, with the highest sky view factor values in a large plain area at a low elevation (lower than 250 m), as shown in Figure 6(a1). The elevation in BTM largely ranges from 450-1300 m, and the sky view factor is mainly concentrated within the range from 0.7-0.8 (please refer to Figure 6(b1)). As shown in Figure 6(a2,b2), the sky view factor in both study areas tends to decrease with increasing slope. Figure 5. Spatial distribution of the slope (a1,b1) and sky view factor (a2,b2) in HR (a1,a2) and BTM (b1,b2). Figure 6 shows scatter plots of the sky view factor versus the elevation and slope in the two study areas. The elevation in HR is clearly concentrated within the range from 50-1100 m, and the sky view factor ranges from 0.7-0.875, with the highest sky view factor values in a large plain area at a low elevation (lower than 250 m), as shown in Figure 6a1. The elevation in BTM largely ranges from 450-1300 m, and the sky view factor is mainly concentrated within the range from 0.7-0.8 (please refer to Figure 6b1). As shown in Figure 6a2,b2, the sky view factor in both study areas tends to decrease with increasing slope.  Figure 7, based on radiation observations at the stations nearest the two study areas, and the equations and related statistics are summarized in Table  4. In regard to the horizontal direct radiation ( ↓_ ), the linear fit was the worst, the quadratic and cubic polynomials were very similar, and the cubic polynomial yielded the best fit. The adjusted coefficient of determination (R 2 ) values were 0.760 and 0.708 for HR and BTM, respectively, while the RMSE values were 0.115 and 0.111, respectively. Regarding the horizontal sky diffuse radiation ( ↓_ ), the best fit was still obtained with the cubic polynomial, with adjusted R 2 values of 0.588 and 0.350 for HR and BTM, respectively, and RMSE values of 0.052 and 0.048, respectively.  Calibration results of the relationship between R s↓_dir0 /R a and R s↓_di f 0 /R a and sunt/ sunt max are shown in Figure 7, based on radiation observations at the stations nearest the two study areas, and the equations and related statistics are summarized in Table 4. In regard to the horizontal direct radiation (R s↓_dir0 ), the linear fit was the worst, the quadratic and cubic polynomials were very similar, and the cubic polynomial yielded the best fit. The adjusted coefficient of determination (R 2 ) values were 0.760 and 0.708 for HR and BTM, respectively, while the RMSE values were 0.115 and 0.111, respectively. Regarding the horizontal sky diffuse radiation (R s↓_di f 0 ), the best fit was still obtained with the cubic polynomial, with adjusted R 2 values of 0.588 and 0.350 for HR and BTM, respectively, and RMSE values of 0.052 and 0.048, respectively.
The abovementioned calibrated equations with the cubic polynomial were substituted into Equations (3) and (6) to obtain solar and net radiation results, respectively. In the validation process of the net radiation calculation model, the model calculation results were compared to the obtained ground observations (Figure 8  The abovementioned calibrated equations with the cubic polynomial were substituted into Equations (3) and (6) to obtain solar and net radiation results, respectively. In the validation process of the net radiation calculation model, the model calculation results were compared to the obtained ground observations (Figure 8), and the accuracy of the net radiation model calculation results in both HR and BTM was satisfactory, with R 2 values of 0.87 and 0.85, respectively, and RMSE values of 15.64 and 25.90 W/m 2 , respectively. With reference to the mean error (ME), the HR results were relatively close to the overall average value (4.45 W/m 2 ), while the BTM results were overestimated (10.54 W/m 2 ). The validation results indicated that the accuracy of the net radiation model calculation results is satisfactory, and the validated model can be adopted in the subsequent ET calculations.   Figure 9 shows the day-by-day variation in the calculation results of the net radiation model compared to the in-situ observations, including the downward shortwave radiation. The variation patterns of the net radiation model calculation results were consistent with the in-situ observations in both study areas. As an essential component of the net radiation balance at the surface, the downward shortwave radiation is a major energy source for processes such as plant growth and atmospheric circulation. Figure 9 clearly shows that the downward shortwave radiation greatly influences the net radiation.  Figure 9 shows the day-by-day variation in the calculation results of the net radiation model compared to the in-situ observations, including the downward shortwave radiation. The variation patterns of the net radiation model calculation results were consistent with the in-situ observations in both study areas. As an essential component of the net radiation balance at the surface, the downward shortwave radiation is a major energy source for processes such as plant growth and atmospheric circulation. Figure 9 clearly shows that the downward shortwave radiation greatly influences the net radiation.  Figure 10 shows a comparison of the net radiation model performance in the different months, characterized by the coefficient of determination (R 2 ), relative mean error (rME), and relative root mean square error (rRMSE), between the net radiation model calculation   Figure 9 shows the day-by-day variation in the calculation results of the net radiation model compared to the in-situ observations, including the downward shortwave radiation. The variation patterns of the net radiation model calculation results were consistent with the in-situ observations in both study areas. As an essential component of the net radiation balance at the surface, the downward shortwave radiation is a major energy source for processes such as plant growth and atmospheric circulation. Figure 9 clearly shows that the downward shortwave radiation greatly influences the net radiation.  Figure 10 shows a comparison of the net radiation model performance in the different months, characterized by the coefficient of determination (R 2 ), relative mean error (rME), and relative root mean square error (rRMSE), between the net radiation model calculation  Figure 10 shows a comparison of the net radiation model performance in the different months, characterized by the coefficient of determination (R 2 ), relative mean error (rME), and relative root mean square error (rRMSE), between the net radiation model calculation results and in-situ observations. As shown in Figure 10a, the net radiation models in both study areas performed better in summer and autumn than in winter and spring, with the highest correlation coefficient value of 0.9 in August in HR, and exceeding 0.7 in May, June, and September. The highest correlation coefficient value in BTM was 0.85 in April, exceeding 0.7 in all months except January, October, November and December. The net radiation results for HR were overestimated in all months except for March, April, and October, while the net radiation results for BTM were slightly underestimated in November. In terms of rME (Figure 10b), the deviation in BTM was large. Overall, the deviation in February, July, August, September and December exceeded 15%. The largest deviation in HR occurred in August at 21%, while the absolute values of rME in March, April, October, November and December were less than 10%. As shown in Figure 10c, the variation in rRMSE in HR ranged from 14% (December) to 25% (February), and the variation in rRMSE in BTM ranged from 17% (March) to 44% (February).
June, and September. The highest correlation coefficient value in BTM was 0.85 in April, exceeding 0.7 in all months except January, October, November and December. The net radiation results for HR were overestimated in all months except for March, April, and October, while the net radiation results for BTM were slightly underestimated in November. In terms of rME (Figure 10b), the deviation in BTM was large. Overall, the deviation in February, July, August, September and December exceeded 15%. The largest deviation in HR occurred in August at 21%, while the absolute values of rME in March, April, October, November and December were less than 10%. As shown in Figure 10c, the variation in rRMSE in HR ranged from 14% (December) to 25% (February), and the variation in rRMSE in BTM ranged from 17% (March) to 44% (February). Based on the spatial distribution of the annual average net radiation (Figure 11), the net radiation on built-up surfaces in the plain area to the southeast of HR was low, and the net radiation over water bodies was the highest. Moreover, BTM is mainly mountainous, and the variation in net radiation in this area was low. In both study areas, the variation in net radiation in the vegetation-covered mountainous area could indicate certain textures with the topography, reflecting the influence of the topography on the net radiation to a certain extent. Based on the spatial distribution of the annual average net radiation (Figure 11), the net radiation on built-up surfaces in the plain area to the southeast of HR was low, and the net radiation over water bodies was the highest. Moreover, BTM is mainly mountainous, and the variation in net radiation in this area was low. In both study areas, the variation in net radiation in the vegetation-covered mountainous area could indicate certain textures with the topography, reflecting the influence of the topography on the net radiation to a certain extent.

ET Results
Validation of the ET model was accomplished by comparing the model calculation results to EC data, and the validation results are shown in Figure 12a,b. The R 2 value was 0.84 for HR and 0.86 for BTM, and from the perspective of deviation, the model calculation results were slightly higher in both study areas, with ME values of −0.39 and 0.11 mm for HR and BTM, respectively. The RMSE was lower in BTM (0.59 mm) than that in HR (0.82 mm). The ET model performed relatively well in both topographically complex study areas.  Figure 13 compares the day-by-day variation in the ET model results to the in-situ observations, including the net radiation. Again, the day-by-day variation in the ET model results for both study areas suitably agreed with that in the in-situ observations. The ET variation pattern was similar to that of the net radiation. Net radiation is an essential component of the surface energy balance, and as the energy source of the latent heat flux (ET), this parameter highly influences ET.

ET Results
Validation of the ET model was accomplished by comparing the model calculation results to EC data, and the validation results are shown in Figure 12a,b. The R 2 value was 0.84 for HR and 0.86 for BTM, and from the perspective of deviation, the model calculation results were slightly higher in both study areas, with ME values of −0.39 and 0.11 mm for HR and BTM, respectively. The RMSE was lower in BTM (0.59 mm) than that in HR (0.82 mm). The ET model performed relatively well in both topographically complex study areas.

ET Results
Validation of the ET model was accomplished by comparing the model calculation results to EC data, and the validation results are shown in Figure 12a,b. The R 2 value was 0.84 for HR and 0.86 for BTM, and from the perspective of deviation, the model calculation results were slightly higher in both study areas, with ME values of −0.39 and 0.11 mm for HR and BTM, respectively. The RMSE was lower in BTM (0.59 mm) than that in HR (0.82 mm). The ET model performed relatively well in both topographically complex study areas.  Figure 13 compares the day-by-day variation in the ET model results to the in-situ observations, including the net radiation. Again, the day-by-day variation in the ET model results for both study areas suitably agreed with that in the in-situ observations. The ET variation pattern was similar to that of the net radiation. Net radiation is an essential component of the surface energy balance, and as the energy source of the latent heat flux (ET), this parameter highly influences ET.  Figure 13 compares the day-by-day variation in the ET model results to the in-situ observations, including the net radiation. Again, the day-by-day variation in the ET model results for both study areas suitably agreed with that in the in-situ observations. The ET variation pattern was similar to that of the net radiation. Net radiation is an essential component of the surface energy balance, and as the energy source of the latent heat flux (ET), this parameter highly influences ET.  Figure 14 shows a comparison of the ET model performance in the different months between the ET calculation results and in-situ observations. As shown in Figure 14a, both study areas exhibited a good performance. The highest R 2 value in HR was 0.84 in April, while only January, February, March, November and December yielded R 2 values lower than 0.5. The R 2 value in August was the highest, namely 0.87 in BTM, and there were only four months with R 2 values below 0.4 (March, April, November and December). As shown in Figure 13b, the ET model results were overestimated in February, March, April, and June in HR, while in BTM, the ET model results were underestimated in January, February, May, June, and July. As shown in Figure 14b, the largest deviation (rME) occurred in March (exceeding 70%), while the absolute values from May to October were less than 10%. The model results were all overestimated in HR. The largest deviation (rME) in HR occurred in June (51%), while the absolute values remained within 10% in January, February and September. Regarding rRMSE (Figure 14c), the seasonal variation characteristics were similar to those of rME. The rRMSE values varied between 33% (September) and 82% (March) in HR and between 15% (June) and 102% (March) in BTM. Overall, the accuracy of the ET model results was satisfactory during the main vegetation growing season.  Figure 14 shows a comparison of the ET model performance in the different months between the ET calculation results and in-situ observations. As shown in Figure 14a, both study areas exhibited a good performance. The highest R 2 value in HR was 0.84 in April, while only January, February, March, November and December yielded R 2 values lower than 0.5. The R 2 value in August was the highest, namely 0.87 in BTM, and there were only four months with R 2 values below 0.4 (March, April, November and December). As shown in Figure 13b, the ET model results were overestimated in February, March, April, and June in HR, while in BTM, the ET model results were underestimated in January, February, May, June, and July. As shown in Figure 14b, the largest deviation (rME) occurred in March (exceeding 70%), while the absolute values from May to October were less than 10%. The model results were all overestimated in HR. The largest deviation (rME) in HR occurred in June (51%), while the absolute values remained within 10% in January, February and September. Regarding rRMSE (Figure 14c), the seasonal variation characteristics were similar to those of rME. The rRMSE values varied between 33% (September) and 82% (March) in HR and between 15% (June) and 102% (March) in BTM. Overall, the accuracy of the ET model results was satisfactory during the main vegetation growing season. Figures 15 and 16 show month-by-month spatial distributions of the ET model calculation results in 2020 in HR and 2019 in BTM. The monthly ET changes in the two study areas reflected a pronounced seasonality, gradually increasing with vegetation growth, peaking in summer, and gradually decreasing thereafter. Since the net radiation calculation considered the topography influence, the ET results could also reflect the differences between the various topographic conditions to a certain extent.
ote Sens. 2022, 14, x 18 of Figure 14. Monthly performance of the ET model expressed by the coefficient of determination (a), relative mean error, rME (b), and relative root mean square error, rRMSE (c) in HR (blue b and BTM (orange bar). Figures 15 and 16 show month-by-month spatial distributions of the ET model cal lation results in 2020 in HR and 2019 in BTM. The monthly ET changes in the two stu areas reflected a pronounced seasonality, gradually increasing with vegetation grow peaking in summer, and gradually decreasing thereafter. Since the net radiation calcu tion considered the topography influence, the ET results could also reflect the differen between the various topographic conditions to a certain extent.

Model Performance
In this study, topographic factors, such as the elevation and sky view factor, were introduced to calculate the downward shortwave radiation, and a net radiation model considering topographic factors was constructed. Then, the daily net radiation calculation results were incorporated into a remote sensing based ETWatch model with a high spatial resolution under complex topographic conditions. The model performance was satisfactory, and a similar accuracy was achieved to that reported in previous studies [6,[92][93][94].
The incorporation of high-resolution DEM data, combined with high-resolution Sentinel 2 remote sensing data, provided net radiation results with a high spatial resolution and accuracy. To explore whether the involvement of topographic factors yields more accurate net radiation simulation results, we also calculated the net radiation with the ETWatch model (Model_ETW), which does not account for topographic factors. Comparing the net radiation calculation results of the present model to those of Model_ETW and those based on the Global Land Data Assimilation System (GLDAS) [95], Climate Forecast System version 2 (CFSv2) [96], ERA5 [97], and Clouds and the Earth's Radiant Energy System (CERES) [98] datasets considering the in-situ Rn observations, both study areas obtained the highest accuracy with the model proposed in this study in regard to the correlation coefficient (R) and RMSE ( Figure 17). Consideration of topographic factors yielded Rn model results with a higher accuracy than that obtained with Model_ETW, and all Rn datasets in both HR and BTM. The net radiation data mentioned above were incorporated into the ET calculation process, and the ET calculation results were compared in a similar manner, as shown in Figure 18. The results revealed that the ET calculation results based on the Rn model data suitably agreed with the observation data in both study areas, thus generating the best performance. Furthermore, the accuracy of the ET model results was slightly improved with the help of the terrain-considered Rn results.

Model Performance
In this study, topographic factors, such as the elevation and sky view factor, were introduced to calculate the downward shortwave radiation, and a net radiation model considering topographic factors was constructed. Then, the daily net radiation calculation results were incorporated into a remote sensing based ETWatch model with a high spatial resolution under complex topographic conditions. The model performance was satisfactory, and a similar accuracy was achieved to that reported in previous studies [6,[92][93][94].
The incorporation of high-resolution DEM data, combined with high-resolution Sentinel 2 remote sensing data, provided net radiation results with a high spatial resolution and accuracy. To explore whether the involvement of topographic factors yields more accurate net radiation simulation results, we also calculated the net radiation with the ETWatch model (Model_ETW), which does not account for topographic factors. Comparing the net radiation calculation results of the present model to those of Model_ETW and those based on the Global Land Data Assimilation System (GLDAS) [95], Climate Forecast System version 2 (CFSv2) [96], ERA5 [97], and Clouds and the Earth's Radiant Energy System (CERES) [98] datasets considering the in-situ R n observations, both study areas obtained the highest accuracy with the model proposed in this study in regard to the correlation coefficient (R) and RMSE ( Figure 17). Consideration of topographic factors yielded R n model results with a higher accuracy than that obtained with Model_ETW, and all R n datasets in both HR and BTM. The net radiation data mentioned above were incorporated into the ET calculation process, and the ET calculation results were compared in a similar manner, as shown in Figure 18. The results revealed that the ET calculation results based on the R n model data suitably agreed with the observation data in both study areas, thus generating the best performance. Furthermore, the accuracy of the ET model results was slightly improved with the help of the terrain-considered R n results.  As shown in Figures 11, 15 and 16, the Rn results and final ET results in this study revealed certain topographic characteristics in terms of the spatial distribution. Rn in both HR and BTM decreased with increasing slope (Figure 19(a1,b1) Figure  19(a2,b2), respectively). In both study areas, the mean net radiation values were higher on both sunny and semisunny slopes than on shady and semishady slopes, respectively. Rn, as the main source of energy in the ET process, also incorporated these topographic features into the final ET results. Figure 20 shows that ET in both study areas also decreased with increasing slope. The mean annual ET value on sunny and semisunny slopes was 616.93 mm, and the mean annual ET value on shady and semishady slopes was 596.22 mm in HR. In contrast, the mean annual ET value on sunny and semisunny slopes reached 761.58 mm, and the mean annual ET value on shady and semishady slopes reached 655.53 mm in BTM. ET on south-facing slopes was higher than that on north-facing slopes in both study areas.  As shown in Figures 11, 15 and 16, the Rn results and final ET results in this study revealed certain topographic characteristics in terms of the spatial distribution. Rn in both HR and BTM decreased with increasing slope (Figure 19(a1,b1) Figure  19(a2,b2), respectively). In both study areas, the mean net radiation values were higher on both sunny and semisunny slopes than on shady and semishady slopes, respectively. Rn, as the main source of energy in the ET process, also incorporated these topographic features into the final ET results. Figure 20 shows that ET in both study areas also decreased with increasing slope. The mean annual ET value on sunny and semisunny slopes was 616.93 mm, and the mean annual ET value on shady and semishady slopes was 596.22 mm in HR. In contrast, the mean annual ET value on sunny and semisunny slopes reached 761.58 mm, and the mean annual ET value on shady and semishady slopes reached 655.53 mm in BTM. ET on south-facing slopes was higher than that on north-facing slopes in both study areas. As shown in Figures 11, 15 and 16, the R n results and final ET results in this study revealed certain topographic characteristics in terms of the spatial distribution. R n in both HR and BTM decreased with increasing slope (Figure 19(a1,b1) (Figure 19(a2,b2), respectively). In both study areas, the mean net radiation values were higher on both sunny and semisunny slopes than on shady and semishady slopes, respectively. R n , as the main source of energy in the ET process, also incorporated these topographic features into the final ET results. Figure 20 shows that ET in both study areas also decreased with increasing slope. The mean annual ET value on sunny and semisunny slopes was 616.93 mm, and the mean annual ET value on shady and semishady slopes was 596.22 mm in HR. In contrast, the mean annual ET value on sunny and semisunny slopes reached 761.58 mm, and the mean annual ET value on shady and semishady slopes reached 655.53 mm in BTM. ET on south-facing slopes was higher than that on north-facing slopes in both study areas. Figure 19. Average daily Rn at the different slopes (a1,b1) and aspect ranges (a2,b2) in HR (a1,a2) and BTM (b1,b2). Aspect range abbreviations such as N, NW and W are explained in Table 5. Figure 20. Average yearly ET at the different slopes (a1,b1) and aspect ranges (a2,b2) in HR (a1,a2) and BTM (b1,b2). Aspect range abbreviations such as N, NW and W are explained in Table 5. Figure 19. Average daily R n at the different slopes (a1,b1) and aspect ranges (a2,b2) in HR (a1,a2) and BTM (b1,b2). Aspect range abbreviations such as N, NW and W are explained in Table 5.  Figure 19. Average daily Rn at the different slopes (a1,b1) and aspect ranges (a2,b2) in HR (a1,a2) and BTM (b1,b2). Aspect range abbreviations such as N, NW and W are explained in Table 5.

Uncertainties and Future Research Directions
Although the validation results indicated that the model performs reliably, there remain certain shortcomings in the proposed framework that should be further investigated, as suggested below.
For the estimation of direct and diffuse radiation for horizontal surface in this study, we used ground-based observations to empirically regress the fraction of direct radiation and diffuse radiation with relative sunshine duration hours, respectively, and selected the model with the best performances for further estimation. The empirical regression method was chosen for this study because of the sufficient observation data. However, in the actual case, when the radiation from the sun is transmitted to the surface, it is partially absorbed by gases such as water vapor and ozone. Besides, the dust, molecules and cloud droplets in the atmosphere also have a scattering effect on the solar radiation. In future studies, models that take into account the actual transmission should be used to estimate direct solar radiation and sky diffuse radiation [99][100][101].
Many existing ET models are based on the energy balance residual concept [102][103][104], such as the ETWatch model. These models require the land surface temperature (LST) as input, which can be limited by the accuracy and spatial resolution of thermal infrared remote sensing data. LST data with a 1-km spatial resolution retrieved from MOD11A1 instead of LST Landsat data were obtained in this study due to the better data availability. Recent studies, such as gap-filling of fine-spatial resolution LST data [105] and spatial sharpening of coarse-resolution LST data [106], have improved the availability and quality of LST data. These key findings should be considered to enhance the prospects of energy balance-based ET models.
In the net radiation calculation model proposed in this study, only the effect of the topography on the downward shortwave radiation was considered. In regard to longwave radiation, the main influencing factors included the cloud cover, surface emissivity, relative humidity, surface temperature, and local atmospheric circulation [107,108]. Considering that the topography also influences longwave radiation to a certain degree, the study of Yan et al. (2020) [109] can be incorporated in future model development to consider the influence of the topography on the directionality of longwave radiation and anisotropy in thermal radiation to obtain terrain-corrected net longwave radiation data, which can then be applied in net radiation calculation.
Soil moisture is an important factor in the process of daily surface conductance reconstruction. The SMAP data considered in this study, with a coarse resolution (10 km), can hardly reflect the soil moisture variation among different underlying surface types and topographic environments. Studies have been reported combining Sentinel-1 and Sentinel-2 images to simulate soil moisture with a high spatial resolution [110][111][112]. Combining these studies should be considered in future studies to introduce high-spatial resolution soil moisture information into ET calculations.
The NDVI is an important parameter in the ET estimation process. Optical remote sensing data can only provide cloud-free daily observations, and the daily NDVI must be obtained through interpolation. The revisit period is generally long for remote sensing data with a high spatial resolution, but the available cloud-free data are often minimal. As such, the daily NDVI obtained via interpolation can overlook information on vegetation status changes on cloudy and rainy days, resulting in the situation in which the interpolated NDVI variations differ from the actual changes. Microwave remote sensing data are generally not affected by atmospheric conditions (clouds and aerosols) and can achieve full-time, all-weather monitoring of vegetation conditions. In future studies, the radar vegetation index obtained from Sentinel-1 [113] should be considered in ET calculations instead of the NDVI. We applied Sentinel-2 only to acquire vegetation information, while the potentials of the Harmonized Landsat and Sentinel-2 (HLS) dataset could be explored in future studies. The HLS dataset incorporates surface reflectance data from both Sentinel-2 Multi-Spectral Instrument (MSI) and Landsat-8 Operational Land Imager (OLI) and can achieve land surface monitoring at the spatial resolution of 30 m every 2 to 3 days [114]. More details on vegetation growth conditions could be found based on HLS datasets, which are beneficial to ET studies over complex terrain.

Conclusions
In this study, a net radiation model considering complex topographic conditions was incorporated into the ETWatch model, and net radiation and ET results with a 10-m spatial resolution were obtained by combining high-spatial resolution remote sensing Sentinel-2 images with meteorological data in the two study areas, i.e., HR and BTM. The ET results were validated against situ observations. The R 2 value was 0.84 in HR and 0.86 in BTM, while the RMSE was lower in BTM (0.59 mm) than that in HR (0.82 mm). The model accuracy was satisfactory, better than ET results with a no-terrain-considered Rn model, and the day-to-day variation in the model results also suitably agreed with that in the in-situ observations. The inclusion of topographic factors, such as the slope and slope direction, ensured that the ET results reflected certain topographic characteristics. The mean annual ET value on sunny and semisunny slopes was 616.93 mm, and the mean annual ET value on shady and semishady slopes was 596.22 mm in HR. In contrast, the mean annual ET value on sunny and semisunny slopes reached 761.58 mm, and the mean annual ET value on shady and semishady slopes reached 655.53 mm in BTM. High spatial and temporal resolutions of ET data in mountainous areas with a complex terrain could contribute to a greater understanding of the characteristics of heat and water fluxes at different vegetation growth stages and substrate types and guide ecological protection and rational allocation of water resources in mountainous areas.
Author Contributions: L.W. was responsible for the experimental design, manuscript preparation and data processing and presentation. B.W. contributed to the conceptual design, manuscript review, funding acquisition and project administration. A.E., Z.M., S.L., X.N., N.Y. and W.Z. contributed to data processing and manuscript review. All authors have read and agreed to the published version of the manuscript.