Regional Daily ET Estimates Based on the Gap-Filling Method of Surface Conductance

Abstract

Remote sensing allows regional evapotranspiration (ET) values to be obtained. Surface conductance is a key variable in estimating ET and controls surface flux interactions between the underlying surface and atmosphere. Limited by the influence of clouds, ET can only be estimated on cloud-free days. In this study, a gap-filling method is proposed to acquire daily surface conductance, which was coupled into a Penman-Monteith (P-M) equation, to estimate the regional daily ET over the Hai River Basin. The gap-filling method is coupled with the canopy conductance, surface conductance and a simple time extension method, which provides more mechanisms and is more comprehensive. Field observations, including eddy covariance (EC) fluxes and meteorological elements from automatic weather station (AWS), were collected from two sites for calibration and validation. One site is located in Guantao County, which is cropped in a circular pattern with winter wheat and summer maize. The other site is located in Miyun County, which has orchard and summer maize crops. The P-M equation was inverted to the computed surface conductance at the field scale, and latent heat fluxes from EC were processed and converted to daily ET. The results show that the surface conductance model used in the gap-filling method performs well compared with the inverted surface conductance, which suggests that the model used here is reasonable. In addition, the relationship between the results estimated from the gap-filling method and EC measurements is more pronounced than that between the other method and the EC measurements. The $R^2$ values improve from 0.68 to 0.75 at the Guantao site and from 0.79 to 0.88 at the Miyun site. The improvement mainly occurs during the growing crop season, according to the temporal variations in the results.

1. Introduction

Energy and water flux exchange is an important physical process in the research field of hydrometeorology and global change. Quantitative research on surface heat flux contributes significantly to water cycle and climatic change interpretations. Evapotranspiration (ET) comprises vegetation transpiration and evaporation from water and soil, which serves as a key component in the interaction between the surface and atmosphere. Global total ET accounts for more than half of the total amount of precipitation. Accurate estimations of ET and its spatiotemporal distribution are meaningful for optimizing water resource management [1,2,3].

Ground observations and ET computations are highly precise in estimating individual sites or at local scales and include observations based on pan evaporation, lysimeter, Bowen ratio, eddy covariance (EC) and large aperture scintillometer (LAS) [4,5,6,7,8,9,10]. Comparably, EC can provide continuously precise measurements of water and heat flux over a homogeneous surface, revealing the character of mass and energy exchange between specific surface features and the atmosphere. EC measurements have been mostly used for parameter calibration and validation in land surface models in previous studies. However, field-scale measurements are limited by their spatial representation and are difficult to promote widely due to the cost. Regarding the underlying surface heterogeneity and turbulence complexity, regional interpolation of ET from sparse measurements is unreliable [11].

Remote sensing enables multi-temporal regional ET monitoring at low costs and multiple scales. Remote sensing-based ET models usually provide model input parameters from remote sensing measurements, such as land surface temperature, albedo and vegetation indexes, to separate sensible heat flux and latent heat flux to solve the surface energy balance equation. Satellites equipped with multi-scale, high-spectrum, multi-view angle sensors are currently designed to provide a more accurate description of surface topography and hydrothermal conditions.

The Penman-Monteith (P-M) [12,13] equation combines aerodynamics with energy balance and considers the surface conductance by mass transfer, which achieves the ET estimation for an unsaturated surface. Surface conductance is a crucial parameter in controlling energy and water exchanges between surface and atmosphere, which is influenced by vegetation-physiology, type of underlying surface and environmental factors and is more complex when mixed pixels exist. As remote sensing techniques develop, leaf area index ($LAI$)-based surface conductance models have been put forward to estimate ET using the P-M equation [14,15,16,17]. However, accurate calibration for model parameters requires field observations that limit the application of models in wider areas.

Due to the difficulty in determining surface conductance, the Surface Energy Balance System (SEBS) [11,18,19] and Surface Energy Balance Algorithm for Land (SEBAL) [1,20,21] were developed to compute sensible heat flux and latent heat flux, and therefore, surface conductance can be acquired as a residual through the energy balance equation. SEBS complemented aerodynamic resistance with residual resistance, which is determined by KB-1, to solve the difference between the land surface temperature and aerodynamic temperature, as well as the difference between mass transfer and heat transfer. SEBS possesses a robust physical basis compared to SEBAL, but SEBS depends on highly accurate meteorological data. In addition, the precision of the model significantly depends on KB-1, which is difficult to determine. SEBAL assumes that the land surface temperature has a linear relationship with the difference between air temperature and aerodynamic temperature, and parameters in the linear relationship are computed by a hot point and cold point selected in the study area. The model has strict requirements for extreme points, which is always ensured empirically. The linear relationship proposed in the model depends on the uniformity of the surface and climatic conditions, which results in a geographic scale problem [22]. Limited by the availability of satellite data, both models only work on cloud-free days.

The Atmosphere-Land Exchange Inverse (ALEXI) model was developed as a robust regional framework for the two-source energy balance (TSEB) land surface model [23], which computes the local partitioning of the surface fluxes and radiometric temperature between the soil and vegetation components of a model grid cell [24]. The ALEXI model coupled the TSEB with an atmospheric boundary layer (ABL) model to simulate the land-atmosphere feedback of air temperature, which is used as a boundary condition for the sensible heat flux [25]. Given morning surface temperature rise measurements from geostationary satellites, ALEXI related the time-integrated sensible heat from the land surface during the morning hours to the rise in ABL temperature without requiring the near-surface air temperature measurement [25,26]. An associated flux disaggregation algorithm (DisALEXI), which is also built on the TSEB land surface representation, can be used to downscale the GEOS-based ALEXI flux estimates to a finer spatial scale using land surface temperature data from moderate- to high-resolution thermal infrared imaging systems [27,28,29,30,31]. Although regional daily ET can be obtained from ALEXI, the resolution depends on the geostationary satellite data, which are always in coarse resolution. Moreover, the models are limited by complexity during the computational process.

The ETWatch coupled residual method with the P-M equation is used to acquire regional daily ET. The residual method estimates sensible heat flux on cloud-free days and latent heat flux can be obtained as a residual from the energy balance equation. Subsequently, surface conductance for cloud-free days can be inverted using the P-M equation and extended to daily scale through several environmental variables. Daily ET can thus be estimated using the P-M equation with daily surface conductance [32,33,34,35]. In this study, an enhanced method for gap-filling of surface conductance is proposed and applied in Hai River Basin. Daily surface conductance can be estimated from those on clear-sky days with the method. The enhanced method couples a surface conductance model with a time extension method, which guarantees a strong, mechanistic physical process and is more comprehensive in considering impact factors. In addition, the gap-filling method can avoid calibration of empirical parameters and direct computation of surface conductance thus to ensuring precision and simplicity in estimating daily regional surface conductance. When coupled with ETWatch, daily regional ET can be obtained without complex computational procedures and are not limited by specific data. Multi-source remote sensing and meteorological data were collected for input. Field observations from two sites were used for calibration and validation.

2. Study Area and Data

2.1. Study Area

The Hai River Basin contains the largest water system in North China, which flows through eight provinces (112°–120°E, 35°–45°N) and covers approximately 320,600 km². The basin has the semi-humid and monsoonal climate typical of mid-latitude regions. The climate is cold during the winter but recovers quickly in the spring, with a mean annual temperature range of 1.5–14 °C. The mean annual precipitation in the basin is 539 mm, and the mean annual ET is 512 mm, most of which occurs during the summer. Figure 1 shows the position and land cover map of the whole basin using multi-temporal MODIS data and an object-based classification method. The underlying surfaces of the basin are mainly mountains and plateaus in northern and western regions and plains in the southern and eastern regions. Winter wheat and summer maize are the two main crops grown on arable land and are rotated with the change of season. Figure 2 shows the phenology of the two crops.

Two ground observation sites within the basin, which have typical crop types, are selected for validation. These sites are Miyun station (MY) and Guantao station (GT). The positions of the two sites are displayed in Figure 1. Site MY is in a valley of Miyun County in Beijing, which is in the north of the basin (40°37′N, 117°19′E). MY has less rainfall and lower temperatures than the basin, and it is covered by 2.5 km² of orchard and summer maize and is surrounded by mountains. Site GT is in the North China plain (36°30′N, 115°07′E), which has a homogeneous surface and is planted with two crops (winter wheat from October to June and summer maize from June to September). At both sites, an observation system consisting of EC and an automatic weather station (AWS) was installed to monitor the flux measurements and meteorological elements.

2.2. Site Measurements

The EC flux data are transferred automatically and simultaneously from the stations to the server. All data are sampled at a 10 Hz frequency and averaged every 30 min with the post-processing software EdiRe (University of Edinburgh) and EddyPro (Li-COR Biosciences) equipped instruments. Sensible and latent heat fluxes for 2016 were collected from both stations, and energy closure calibration was completed based on the Bowen ratio closure method, due to the data energy imbalance. To guarantee the flux data reliability, abnormal values were removed because of instrument malfunction and bad weather.

The 2016 meteorological data collected from the AWS contain six parameters: wind velocity, air temperature, relative humidity, pressure, net radiation and soil moisture. The data are recorded by loggers at an interval of 10 min, and then, they are averaged to acquire daily values. A small amount of data, which obviously deviated from physical meaning, were deleted and gap-filled with a linear interpolation method [36].

The flux data were collected for nearly the whole year from both stations and covered the main crop growing season. The meteorological data for station GT were adequate, whereas the last four months of data from station MY were lost due to an instrument malfunction. The detailed information from collected ground observations are listed in

Table 1. Information from collected ground observations.

Station	Land Cover	Measurements	Period
Guantao	Winter wheat Summer maize	EC	January–December 2016
		AWS
Miyun	Orchard Summer maize	EC	January–December 2016
		AWS	January–Auguat 2016

.

2.3. Remote Sensing Data

In this study, several satellite datasets were collected and include MODIS, NCEP, FY-2D and AMSR-2. Four MODIS images were needed to cover the basin, and the preprocessing included radiometric, atmospheric and geometric corrections, which were already completed by the data supplier. MYD09GA and MYD11A1 are surface reflectance and land surface temperature products, respectively. The data were converted to image files from hierarchical data format (HDF) files, and then, they were mosaicked and clipped to the study area using interactive data language (IDL) programs. Normalized difference vegetation index ($NDVI$) can be calculated using the MODIS 1,2-band, and land surface albedo was converted to broadband from narrowband [37]. Cloud pixels were detected using a method proposed by Ackerman [38]. Time series datasets for $NDVI$ and albedo were computed by linear interpolation and smoothed by Savitzky-Golay filter [39]. MOD07 provided the atmospheric temperature and humidity profiles, and a planetary boundary layer (PBL) mixing height derived method was applied to extract information for a PBL [40]. All the MODIS products are 500 m and 1 km resolution and can be downloaded from the website of LAADS DAAC (https://ladsweb.modaps.eosdis.nasa.gov/). The wind velocity of a PBL can be extracted from NCEP reanalysis data, which has a spatial resolution of 2.5° and is provided by NOAA Earth System Research Laboratory (https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html), according to the height of PBL. FY-2D is the second operational vehicle of the first-generation geostationary meteorological satellite system operated by the Chinese Meteorological Administration (CMA). The sensor aboard the satellite consists of one visible channel with a resolution of 1.25 km and four infrared channels with resolutions of 5 km. Products are distributed on data services website in HDF format (http://satellite.nsmc.org.cn/portalsite/default.aspx). Cloud data from the FY-2D satellite were used to compute the sunshine hours, instead of the interpolated meteorological sunshine hours [41,42]. AMSR-2 is a sensor aboard GCOM-W1, which can provide daily global soil moisture of the top soil with resolution of 10 km (https://gcom-w1.jaxa.jp/auth.html). All the above data were projected to Albers conical equal area projections and resampled to 1 km resolution using an earth resource development assessment system (ERDAS), which is a software package for processing remote sensing images.

In addition, daily meteorological elements, including wind velocity, air temperature, relative humidity and air pressure, which were recorded at 49 national weather stations over the basin, were collected from the China Meteorological Data Service Center. All parameters were corrected with a digital elevation model (DEM) and then interpolated into daily maps at 1 km resolution using the inverse distance weighting method [43]. The interpolated results were projected to Albers conical equal area projections as satellite data.

3. Methodology

3.1. Surface Conductance Model

3.1.1. Canopy Conductance

The canopy conductance ($G_c$) depends both on the atmospheric factors and upon available water content in the soil [44,45,46]. The model used in this study to calculate the canopy conductance ($G_c$) is as follows [47]:

$$G_c=g_{c,max}·LAI·F_1·F_2·F_3·F_4^4,$$

(1)

$LAI$ was obtained using the following empirical equation:

$$LAI=\sqrt{\frac{1+NDVI}{1-NDVI}NDVI},$$

(2)

The factor $F_1$ represents the effect of photosynthetically active radiation and can be defined as follows [48]:

$$F_1=\frac{g_{c,min}/g_{c,max}+f}{f+1},$$

(3)

$$f=0.55\frac{R_s}{R_{gl}}\frac{2}{LAI},$$

(4)

$g_{c,min}$ and $g_{c,max}$ are minimum and maximum values of canopy conductance, respectively. $g_{c,min}$ is relatively insensitive in the calculation above and always set as a small constant ($2\times {10}^{-4}\mathrm{m}/\mathrm{s}$) [46]. $R_s$ is the downward shortwave radiation, and $R_{gl}$ is a limiting parameter associated with the actual underlying surface. $LAI$ is the leaf area index [49].

The second factor $F_2$ characterizes the soil water content influence on the root zone. Due to the difficulty obtaining the soil water content for the root zone in a region, the soil water content product of the top soil from ASMR-2 was used in combination with a transformed equation as follows [50]:

$$F_2=0.1LAI+\left(1-0.1LAI\right)\left\{1-exp\left[\theta \left(-0.5LAI-1\right)\right]\right\},$$

(5)

$$\theta =\frac{{\theta }_g-{\theta }_{wilt}}{{\theta }_{cr}-{\theta }_{wilt}},$$

(6)

${\theta }_g$ is the soil water content, and ${\theta }_{cr}$ is the field capacity, which is approximately 0.75 times the soil’s porosity [51]. ${\theta }_{wilt}$ is the wilting point of the soil water content and is common for many herbaceous species [52]. Values recommended for ${\theta }_{cr}$ and ${\theta }_{wilt}$ for different soil types were adopted [53].

The factor $F_3$ addresses the effects of the water vapour pressure deficit $D_a$ over the canopy and can be described as follows [48,49]. $C_v$ is a species-dependent empirical parameter. A larger $D_a$ leads to a smaller $F_3$ and $G_c$, which increases the resistance of the canopy evaporation.

$$F_3=1-C_v·D_a,$$

(7)

The factor $F_4$ takes into consideration the air temperature constraint and is proposed as follows [54]:

$$F_4=1-0.0016{\left(T_{opt}-T_a\right)}^2,$$

(8)

where $T_a$ is the air temperature above the canopy. $T_{opt}$ is the optimal temperature for the growth of many crops and always set as 298 K. $F_4$ decreases when $T_a$ departs from the optimal temperature.

Several empirical parameters in the equations above are related to the actual surface and should be calibrated using ground observations. In this study, the calibrated results from other studies [53,55] were adopted, which are listed in

Table 2. Surface-dependent parameters settings.

Underlying Surface	${\mathit{g}}_{\mathit{c},\mathit{m}\mathit{a}\mathit{x}} \left(\mathbf{m}/\mathbf{s}\right)$	${\mathit{C}}_{\mathit{v}} \left(h/Pa\right)$	${\mathit{R}}_{\mathit{g}\mathit{l}} \left(W/{\mathbf{m}}^2\right)$
Deciduous broad-leaf forest	1/100	0.025	30
Evergreen needle-leaf forest	1/150	0.025	30
Evergreen broad-leaf forest	1/150	0.025	30
Cropland	1/40	0.023	100
Grassland	1/40	0.0155	100
Wetland	1/150	0.0155	100
Bare soil	1/50	--	--

.

3.1.2. Surface Conductance

Total ET is the sum of transpiration from the plant canopy $E_c$ and evaporation from the soil $E_s$ [15,16].

$$E=E_c+E_s,$$

(9)

The P-M equation has been widely used in previous studies to calculate the ET of the land surface [14,16,17,33,44,45,50,56,57]. Therefore, the P-M equation is used for $E$ and $E_c$ in this study as described below:

$$\lambda E=\frac{\Delta A+\rho c_p{D}_a{G}_a}{\Delta +\gamma \left(1+G_a/G_s\right)},$$

(10)

$$\lambda E_c=\frac{\Delta A_c+\rho c_p{D}_a{G}_a}{\Delta +\gamma \left(1+G_a/G_c\right)},$$

(11)

where $\lambda$ is the latent heat of evaporation, $\Delta$ is the slope of the curve relating saturation water vapour pressure to temperature, $\gamma$ is the psychrometric constant, $A$ and $A_c$ are the available energy absorbed by the land surface and canopy surface, respectively, $\rho$ is the air density, $c_p$ is the specific heat of air at a constant pressure, $G_a$ is the aerodynamic conductance, and $G_s$ and $G_c$ are surface conductance and canopy conductance, respectively. Five variables used here, i.e., $\Delta$, $\gamma$, $\rho$, $c_p$, $D_a$, are assumed to be known or can be estimated using meteorological results. This study adopted $G_a$ without considering the stability corrections because the P-M equation is relatively insensitive to aerodynamic conductance, especially when $G_a$ >> $G_s$ [14]. Here $G_a$ was calculated as Equation (12) described, where $k$ is von Karman’s constant (0.41), $u$ is wind velocity, $z$ is the height of wind velocity, $d$ is zero displacement height, $z_{om}$ and $z_{ov}$ are the roughness lengths to momentum water vapor transfer. The quantities $d$, $z_{om}$ and $z_{ov}$ were estimated using $d=2h/3$, $z_{om}=0.123h$ and $z_{ov}=0.1z_{om}$, where $h$ is canopy height [58]. The estimation of $h$ can be described as Equation (13), and the corresponding parameters including maximum height ($h_{max}$), minimum height ($h_{min}$), maximum $NDVI$ ($NDV{I}_{max}$) and minimum $NDVI$ ($NDV{I}_{min}$) of various underlying surface can refer to field measurements and exist achievements of other researchers [59]. $A$ equals the result of daily net radiation $R_n$ minus the daily soil heat flux $G$. Daily $G$ can be neglected here due to its opposing character in the daytime versus nighttime. Daily $R_n$ can be estimated using sunshine hours from FY-2D, albedo from MODIS products and meteorological results [41].

Soil evaporation has an important role in the total ET, especially in sparse vegetation. The assumption that soil evaporation occurs at some fraction $f$ of the equilibrium rate at the soil surface is as follows [14,15]:

$$G_a=\frac{k^{2}u}{\ln \left(\left(z-d\right)/z_{om}\right)\ln \left(\left(z-d\right)/z_{ov}\right)},$$

(12)

$$h=h_{min}+\frac{h_{max}-h_{min}}{NDV{I}_{max}-NDV{I}_{min}}\left(NDVI-NDV{I}_{min}\right),$$

(13)

$$\lambda E_s=\frac{f\Delta A_s}{\Delta +\gamma },$$

(14)

In the equations above, the total available energy $A$ is partitioned into that absorbed by the canopy surface $A_c$ and by the soil surface $A_s$. The vegetation fraction $f_c$, which is estimated from the $NDVI$, is applied in the ratio of total available energy absorbed by the canopy and soil:

$$A_c=A{f}_c,$$

(15)

$$A_s=A\left(1-f_c\right),$$

(16)

$$f_c={\left(\frac{NDVI-NDV{I}_{min}}{NDVI-NDV{I}_{max}}\right)}^2,$$

(17)

The soil evaporation fraction $f$ used in this study is relevant to the top soil moisture status where evaporation occurs. According to the complementary relationship hypothesis, water vapour pressure deficit and humidity response to the variation in moisture in the top soil via land-atmosphere interactions [60,61], $f$ can be described as follows:

$$f=a{\left(RH\right)}^{D_a/\beta },$$

(18)

where $RH$ is the air humidity, $\beta$ is the empirical parameter, which is set to 100, and $a$ is the added factor for adjustment. Combining Equations (9)–(16), the surface conductance can be obtained with the following:

$$G_s=G_c\left[\frac{1+\frac{\left(1-f_c\right)G_a}{\left(\Delta /\gamma +1\right)G_c}\left[f-\frac{\left(\Delta /\gamma +1\right)\left(1-f\right)G_c}{G_a}\right]+\frac{\rho c_p{D}_a{G}_a}{\Delta R_n}}{1-\left(1-f_c\right)\left[f-\frac{\left(\Delta /\gamma +1\right)\left(1-f\right)G_c}{G_a}\right]+\frac{\rho c_p{D}_a{G}_a}{\Delta R_n}}\right],$$

(19)

3.2. Gap-Filling Method

ETWatch is an integration of the “residue approach” and P-M equation [25,32,33,34,35,43]. The residue approach is an energy balance model combined with the mass transfer method, which can compute ET on cloud-free days based on SEBS. ETWatch consists of several subsystems: net radiation model, soil heat flux model, sensible heat flux mode, aerodynamic roughness length model and PBL model. ETWatch computes ET on cloud-free days with multi-source remote sensing data. The basic files are used as auxiliary data in the computation of subsystems, which include a digital terrain model, latitude file, aspect file and slope file. The sensible heat flux $H$ is computed first, and then, the latent heat flux $\lambda E$ can be obtained as the residue. The latent heat flux can be converted to actual ET in diurnal scales through the latent heat of evaporation and evaporative fraction (EF) method [62,63]. Figure 3 shows the estimation flowchart for ET on cloud-free days.

Daily ET data are difficult to acquire because of cloud cover over the basin. A gap-filling method of surface conductance extended from clear-sky surface conductance, which can be computed with clear-sky ET results from ETWatch and an inversion from the P-M equation, is proposed in this study to generate daily surface conductance and to compute the daily ET using the P-M equation.

$$G_{s,clear}=\frac{\lambda E_{clear}{G}_a\gamma }{\Delta A-\left(\Delta +\gamma \right)\lambda E_{clear}+\rho C_p{D}_a{G}_a},$$

(20)

$$Func=G_c\left[\frac{1+\frac{\left(1-f_c\right)G_a}{\left(\Delta /\gamma +1\right)G_c}\left[f-\frac{\left(\Delta /\gamma +1\right)\left(1-f\right)G_c}{G_a}\right]+\frac{\rho c_p{D}_a{G}_a}{\Delta R_n}}{1-\left(1-f_c\right)\left[f-\frac{\left(\Delta /\gamma +1\right)\left(1-f\right)G_c}{G_a}\right]+\frac{\rho c_p{D}_a{G}_a}{\Delta R_n}}\right],$$

(21)

$$\frac{G_{s,daily}}{G_{s,clear}}=\frac{Fun{c}_{daily}}{Fun{c}_{clear}},$$

(22)

where $G_{s,clear}$ is the surface conductance on clear-sky days, $\lambda E_{clear}$ is the diurnal latent heat flux on clear-sky days resulted from ETWatch, Equation (20) is the inversion of P-M equation. Equation (19) is an expression for surface conductance proposed in Section 3.1.2 and the right part of the equation can be marked as $Func$ shown in Equation (21). The subscript $daily$ and $clear$ represent daily scale and clear-sky days respectively. Once the $G_{s,clear}$ is acquired, daily surface conductance ($G_{s,daily}$) can be estimated using the gap-filling method expressed as Equation (22). $Fun{c}_{daily}$ and $Fun{c}_{clear}$ have the same form according to Equation (21) and consist of a series of empirical parameters including $g_{c,max}$, $R_{gl}$, $C_v$, which remain unchanged for the same kind of land cover over time and several external variables such as vegetation fraction, temperature, radiation, wind velocity, etc., which change significantly. Therefore, the Equation (20) suggests that the daily surface conductance estimation is mainly based on the temporal variation in environmental factors between unclear-sky days and clear-sky days, which also indicates that surface-dependent parameters including minimum canopy conductance, empirical parameters related to the radiation and temperature factors have less important roles in the estimation because of their stabilities over time, although they have different values among land covers.

3.3. Evaluation Index

For the quantitative analysis of model performance, a set of evaluation indexes was selected. The coefficient of determination ($R^2$) is used to provide a measure of consistency between estimation and observation. The mean bias error ($BIAS$) and the root-mean-square error ($RMSE$) are usually used to measure the average difference between estimation and observation. The mean relative error ($MRE$) describes the relative extent of deviation of estimation from observation. The calculation methods of these evaluation indexes are listed as follows:

$$R^2=\frac{{\left[{\sum }_{i=1}^n\left(P_i-\overline{P}\right)\left(O_i-\overline{O}\right)\right]}^2}{{\sum }_{i=1}^n{\left(P_i-\overline{P}\right)}^2{\sum }_{i=1}^n{\left(O_i-\overline{O}\right)}^2},$$

(23)

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{n}},$$

(24)

$$BIAS=\frac{{\sum }_{i=1}^n\left(P_i-O_i\right)}{n},$$

(25)

$$MRE=\frac{100}{n}{\sum }_{i=1}^n\frac{P_i-O_i}{\overline{O}},$$

(26)

where $P_i$ is No.i time estimation and $O_i$ is No.i time measurement. $\overline{P}$ and $\overline{O}$ are mean of estimations and measurements respectively, and $n$ is the number of measurements.

4. Results

4.1. Estimation of Surface Conductance

Surface conductance can be estimated as the inversion of the P-M equation with EC and meteorological measurements. Therefore, the surface conductance model described in Section 3.1 was calibrated using inverted values calculated with part of site measurements and key parameters, such as $g_{c,max}$ and $a$, which were determined at sites GT and MY. The surface conductance estimations using satellite and meteorological data at both sites were compared to the inverted results using site measurements. The validation results are shown in Figure 4. The sample numbers indicate that both satellite data and EC measurements are of good quality.

As shown in the linear fitting results in Figure 4, GT shows a high coefficient of determination value ($R^2$ = 0.82) and low root mean square error value ($RMSE$ = 0.0025 m/s). MY shows values similar to those of GT. Both sites have smaller fluctuations compared to the inverted results. In addition, the plots exhibit an obvious seasonal change with the crop calendar. The surface conductance estimations are relatively high values (scattered points greater than 0.015 m/s) from April to May and July to August at GT and from July to August at MY, which are consistent with the growth of the winter wheat and summer maize. A significant correlation between the estimation and inverted results indicated that the surface conductance model mentioned in this study can effectively estimate the surface conductance, when the empirical parameters are accurately calibrated.

4.2. Comparison between Different Methods with EC Measurements

The method proposed in this study (hereafter referred to as the enhanced method) and another method (hereafter referred to as original method) were coupled into the P-M equation to estimate the daily basin ET, and the comparison was executed to evaluate the performance of these two methods. The original method estimates the daily surface conductance by multiplication and temporal variation of limited variables including $LAI$, net radiation ($R_n$), soil moisture ($SM$) and wind velocity ($U$) [33]. Equation (27) gives the expression of the original method. Figure 5 reports the estimation of daily ET plotted against the EC measurements at the GT and MY sites.

The figure above clearly demonstrates that two methods can estimate ET to a reasonable extent, while the enhanced method shows less dispersed results and a higher correlation with the measurements than the original method. For site GT, the $R^2$ values increase from 0.68 to 0.75, and the $RMSE$ decreased from 0.90 to 0.79. The improvement mainly occurred during the crop growing season, which can be proven because the high ET values in the scatter plots displayed various distributions among the two methods. This finding indicates that the enhanced method can better describe the surface conductance variation in the crop growth. The same analysis was used to evaluate the differences between the two methods at site MY, and comparable results were obtained. The enhanced method also performed better than the original method, with $R^2$ values that improved from 0.79 to 0.88 and an $RMSE$ that reduced from 0.80 to 0.62.

$$G_{s,daily}=\frac{G_{s,clear}\times {\left(LAI\times R_n\times SM\times U\right)}_{daily}}{{\left(LAI\times R_n\times SM\times U\right)}_{clear}},$$

(27)

Due to the gaps in time series EC measurements and for convenient comparability, monthly averaged ET was calculated from the field observations and remote sensing estimations. Figure 6 shows a comparison of the results.

As shown in Figure 2, winter wheat and summer maize reach their growth peaks in May and August, respectively. For site GT, the ET temporal variations correspond to the crop calendar, which fluctuates throughout the whole year and exhibits a relatively low value because of the winter wheat harvest. ET estimated from the enhanced method shows a characteristic advance compared to the original method, especially during the growing season. The former has less distinction with measurements in May and August than the latter, which suggests that the surface conductance model used in this paper can better reflect the physiological changes in the crops. For site MY, ET has high values from June to August, following the orchard and summer maize growth periods. The enhanced method also better estimates the peak growth of summer maize in August, compared to the original method.

Several statistical indices were selected to evaluate the estimation accuracy and consistency of measurements. The results are shown in

Table 3. Comparison statistics of the daily and monthly averaged ET for both methods at the Guantao and Miyun sites.

Sites		Original Method				Enhanced Method
		${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{B}\mathit{I}\mathit{A}\mathit{S}$	$\mathit{M}\mathit{R}\mathit{E}$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{B}\mathit{I}\mathit{A}\mathit{S}$	$\mathit{M}\mathit{R}\mathit{E}$
Daily	Guantao	0.68	0.90	−0.42	−22.1%	0.75	0.79	−0.26	−12.5%
	Miyun	0.79	0.80	−0.07	−4.2%	0.88	0.62	−0.05	−3.3%
Monthly Averaged	Guantao	0.96	0.29	−0.38	−22.5%	0.97	0.27	−0.27	−14.6%
	Miyun	0.97	0.26	−0.06	−3.1%	0.98	0.23	−0.05	−2.8%

. The $R^2$ and $RMSE$ values of monthly averaged ET show the same trend as the daily results for both sites. The $MRE$ values range from −14.6% to −12.5% for the GT site and from −3.3% to −2.8% for the MY site, which indicates a smaller estimation deviation of the enhanced method from the measurements. The $BIAS$ values also show notable differences among the results of both methods and show a preferable performance from the enhanced method. Both methods slightly underestimated the ET compared with field observations.

4.3. Spatiotemporal Distribution of Results

The gap-filling method proposed in this study was applied combined with P-M equation to estimate regional daily ET for Hai River Basin. As shown in Figure 7, the ET spatial patterns showed a decrease from southeast to northwest and were consistent with the precipitation pattern. The precipitation in the east and southeast regions was significantly greater than in other regions, and accordingly, the ET values were high in this area. Further analysis of the land cover map (Figure 1) revealed that crops in the southeast region consumed more water than the forest in the northwest region and urban areas. In addition, the area around Bohai Bay showed the highest ET, where the surface was moist, and the rainfall was abundant. Crop land and forest were further classified into dry farmland and paddy field and into arbour forest and shrub forest, respectively, based on Figure 1. The annual ET values for different land covers were counted and are displayed in Figure 8. Water surface was the highest, with more than 700 mm. Paddy field was slightly higher than dry farmland, due to the different crops. For the forest, arbour evaporated more than shrub, while the ET of forest was less than the cropland. Grassland and artificial surface were the lowest among these classes. The difference between ET and precipitation in the artificial surface resulted from many urban areas being covered by impervious surfaces, such as roads and buildings. In addition, Figure 9 shows the temporal variation in the monthly ET over the whole basin. The values increased gradually from January to May and decreased in June, which was influenced by the variation in meteorological elements and phenology of winter wheat. A rapid recover in July and August was shown because of the growth of summer maize.

5. Discussion

Surface conductance is a key parameter in the P-M equation and is difficult to estimate, especially at the regional scale because of its complex underlying surface conditions. In this study, daily surface conductance is estimated from inverted surface conductance for a clear day to meet the requirements of the P-M equation to acquire the regional ET. The gap-filling method used here contains the estimation methods of the surface conductance and temporal change in surface variations, considering both the surface conductance mechanism and environmental factors.

Surface conductance has many impacting factors and cannot be displayed well by the model. To decrease the ET error when surface conductance is used in this study, a surface conductance simulation from the ET inversion is obtained through the P-M equation. Otherwise, ET is separated into two parts: vegetation transpiration calculated from the P-M equation and soil evaporation calculated from the equilibrium ET. Correspondingly, net radiation is also separated into two parts through the vegetation fraction, which simplifies the partition process of net radiation compared with others but matters little in the ET accuracy. In addition, the soil EF calculated through the land-atmosphere interactions replaces the soil moisture usage, can improve the applicability of the soil evaporation model and avoids soil conductance. For the canopy conductance, the Jarvis model is applied, because it performed well in many previous studies [44,46,47]. According to the results shown in Figure 4, the estimations from the surface conductance model have a near-linear relationship with the inversion results from the in situ measurements.

Accurate estimation of surface conductance requires calibration of several key parameters, including maximum canopy conductance, empirical parameters related to radiation, temperature, vapour pressure deficit and the adjustment factor in the soil evaporation calculation. All parameters should be calibrated with in situ measurements, meaning the field measurements for each type of land cover are needed for calibration, which is hard to accomplish. The surface conductance gap-filling method can avoid the calibration and adopt the reference values from other studies, as shown in Table 2. In this study, the surface conductance model is not used to calculate conductance values; it is only used for expression. Daily surface conductance is estimated from clear day surface conductance, which mainly depends on the environmental factor variations. The parameters mentioned above stay constant over time and therefore do not influence the results.

ET differs significantly in various land cover types. The method discussed here is applied only in two sites for two crops. The results are limited by the field datasets for specific underlying surfaces. Further research should be focus on the applicability of the method to other crops and land cover types, as well as in arid and semi-arid regions. In this study, two sites are surrounded by homogeneous and flat arable land. Therefore, field measurements are representative and matched for the MODIS pixels. As shown in Figure 5 and Figure 6, the method proposed in this study performed better than the original method. The latter contains only three environmental factors and adopted a multiplication form. Comparably, the former is more considerate and has a better mechanism. Surface conductance is strongly determined by the surface itself and the environmental factors, especially the meteorological elements. In the future, the contribution of variables for surface conductance should be studied to select the main variables, which will benefit the surface conductance simulation.

6. Conclusions

This study proposed a gap-filling method to acquire daily surface conductance from clear day surface conductance estimation. The method combines a time extension method with the canopy and a surface conductance model and is carried out with multi-source remote sensing data, coupled with the P-M equation to estimate the regional daily ET. The field observations over the winter wheat and summer maize for the two sites are obtained to validate the method’s performance. The results show that the method performs well at both sites, with the $R^2$ values reaching 0.75 and 0.88 at the Guantao and Miyun sites, respectively. Compared to other methods, the method proposed in this study increases the accuracy of the ET estimation to different degrees. The reason attributed to improvements might be the comprehensiveness of the impact factors and mechanism of the latter method. However, limited by the availability of in situ data in this study, this method is applied to only two crops. Further applications are necessary in different land cover types and climatic conditions.