Improving Meteorological Input for Surface Energy Balance System Utilizing Mesoscale Weather Research and Forecasting Model for Estimating Daily Actual Evapotranspiration

: Using Surface Energy Balance System (SEBS) to estimate actual evapotranspiration (ET) on a regional scale generally uses gridded meteorological data by interpolating data from meteorological stations with mathematical interpolation. The heterogeneity of underlying surfaces cannot be effectively considered when interpolating meteorological station measurements to gridded data only by mathematical interpolation. This study aims to highlight the improvement of modeled meteorological data from the Weather Research and Forecasting (WRF) mesoscale numerical model which fully considers the heterogeneity of underlying surfaces over the data from mathematical interpolation method when providing accurate meteorological input for SEBS model. Meteorological data at 1 km resolution in the Hotan Oasis were simulated and then were put into SEBS model to estimate the daily actual ET. The accuracy of WRF simulation was evaluated through comparison with data collected at the meteorological station. Results found that the WRF-simulated wind speed, air temperature, relative humidity and surface pressure correlate well with the meteorological stations measurements (R 2 are 0.628, 0.8242, 0.8089 and 0.8915, respectively). Comparison between ET calculated using the meteorological data simulated from the WRF (ETa-WRF) and meteorological data interpolated from measurements at met stations (ETa-STA) showed that ETa-WRF could better reflect the ET difference between different land cover, and capture the vegetation growing trend, especially in areas with sparse vegetation, where ETa-STA intends to overestimate. In addition, ETa-WRF has less noise in barren areas compared to ETa-STA. Our findings suggest that WRF can provide more reliable meteorological input for SEBS model than mathematical interpolation method. 0.8242, 0.8089 and 0.8915. This evaluation indicated that WRF-simulated meteorological parameters could be used for the estimation of ET.


Introduction
Evapotranspiration (ET) is an important parameter of the water cycle [1,2], which characterizes the transport of moisture from the land surface to the atmosphere. It consists of soil evaporation and vegetation transpiration, and plays an important role in regulating local and regional climate. On the global scale, ET accounts for more than half of the global total precipitation [3][4][5], which means that most of the water transported from the atmosphere to the surface will be not stored. Therefore, accurate estimation of ET and its spatial distribution is important to regional water resources management and the macrocontrol decision making [6].
Traditional ET estimation used meteorological parameters from stations [7][8][9][10][11]. Although it can provide relatively accurate ET over the stations, the heterogeneity of underlying surface and hydrothermal transport on a large scale makes single point ET unable to meet the needs of regional research [12,13].
Surface Energy Balance System (SEBS) is a commonly used model for regional ET estimation based on surface energy balance [14][15][16][17]. Its theoretical basis is heat exchange roughness model, BAS (bulk atmospheric boundary layer (ABL) similarity) and MOS (Monin-Obukhov similarity). First frictional velocity, sensible heat flux and Obukhov stability length are determined and then the actual ET is calculated using surface energy balance index [18,19]. SEBS, using satellite observation, in situ measurements, and meteorological model simulation as input for ET estimation [18], has been validated in sparsely vegetated and dry areas [20][21][22]. To estimate regional ET, SEBS requires accurate gridded meteorological data as input [23]. Earlier studies used gridded meteorological data interpolated from site measurements using mathematical interpolation [24][25][26][27] to provide input for the SEBS model. However, using mathematical interpolation from station measurements cannot effectively consider the heterogeneity of the underlying surface [28]. Most weather stations are scattered, and meteorological elements such as wind speed, humidity and other sensitive variables which are closely related to ET estimation are significantly discrepant due to the influence of the topography and surrounding environment [29,30]. Therefore, high spatial resolution meteorological data are needed for ET estimation using SEBS especially in heterogeneous areas.
The Weather Research and Forecasting (WRF) model is a new generation of fully compressible, non-static and high resolution mesoscale forecasting models, which has been developed by many scientific research institutes in the United States since 1997 [31]. It is a model system that integrates numerical weather forecasting, atmospheric simulation and data assimilation, and adopts an improved physical process scheme with multiple levels of nesting and the ability to easily locate in different geographic locations, and is widely used in mesoscale weather simulation and prediction [32,33]. WRF can simulate high spatio-temporal resolution meteorological elements, while fully considering the underlying surface heterogeneity.
This study will use the WRF model to simulate gridded meteorological data with high spatial resolution, which will be used as input for SEBS model to estimate ET. The accuracy of WRF simulation results is validated using observed meteorological data. Then the daily actual ET of the growing season, calculated using the meteorological data from the WRF (ETa-WRF) and from mathematic interpolation of the stations (ETa-STA) are compared to evaluate the improvement of using WRF simulation as input for SEBS.

Study Area
The Hotan Oasis (79.10°-80.71°E, 36.81°-37.81°N) is located in the southern edge of the Taklimakan Desert in Xinjiang, China (Figure 1a). The total area of the Hotan oasis area is 8420 Km 2 . The south of the oasis with high elevation is adjacent to the mountainous area, and the north of the oasis is an alluvial-diluvial inclined plain with relatively low elevation. This is a semi-arid area with annual average precipitation less than 50 mm, annual evaporation above 2600 mm, and annual average relative humidity lower than 39%. The crops of the Hotan Oasis are mainly wheat, maize, cotton and rice. Usually, walnuts and wheat are intercropped. May to September is the growing season of crops in the Hotan Oasis, and the actual ET of the study area from May to September of 2015 is calculated in this paper.

Data
Data used for WRF simulation and ET calculation in this paper are listed in Table 1. The 300 m 2015 Global Land Cover Data (CCI-LC2015) released by the Climate Change Initiative (CCI) project of European Space Agency (ESA) is used to update the old underlying surface in WRF. The reason for choosing this map is the high classification accuracy for land covers with distinct boundaries. For example, the classification accuracy of dry farmland (such as wheat and maize) can reach 89%-92%, and that of irrigated farmland (such as rice) can reach 83%-89% [34]. Comparison with the field survey data found that CCI-LC2015 has good classification accuracy in the Hotan Oasis.
The data used for the required initial and boundary conditions for driving the WRF model are Final Operational Global Analysis (FNL), produced jointly by the National Center for Environmental Forecasting (NCEP) and the National Center for Atmospheric Research (NCAR) of the United States, with spatial resolution of 1° × 1°, and temporal resolution of 6 h. NCEP claimed that it may be the best choice for archiving and analyzing data of long-term business models [35]. It has advantages in small-and medium-scale weather analysis and can compensate for the shortcomings of conventional observation data in weather analysis [36].
The required satellite imagery is MODIS Terra 500 m resolution surface reflectance products (MOD09A1) and 1 km resolution land surface temperature and specific emissivity products (MOD11A2). 30 m ASTER GDEM V2 data developed by METI and NASA are used to provide topographic information. The Digital Elevation Model (DEM) and derived slope and aspect data are utilized to calculate the solar radiation received at the surface, considering the shading effect of terrain [37].
Meteorological measurements of wind speed, air temperature, relative humidity, and surface pressure were collected from Hotan meteorological station and three surrounding stations: Pishan, Yutian and Minfeng (see Figure 1). Meteorological measurements were used to obtain gridded meteorological data. The method of inverse distance and the influence of elevation has been taken into account when conducting interpolation.
A 30 m land cover map was derived from Landsat-8/OLI imagery, and evaluated through field survey in 2015.

Methods
To improve the ET estimated by the SEBS model, the WRF-simulated wind speed, air temperature, relative humidity, and surface pressure are used as input for SEBS. First, wind speed, air temperature, relative humidity, and surface pressure are simulated using WRF and evaluated using site measurements. Further, comparison is performed between ET estimates using WRF simulation as input and using data interpolated from site measurement as input. The flowchart of calculating ET is shown in Figure 2.

WRF Model
Land heterogeneity not only determines the microclimate but also affects mesoscale atmospheric circulation [38,39]. The accurate land cover will conduct reliable and accurate numerical modelling [40][41][42]. Therefore, the old land cover data in WRF is first updated with a 300 m global land cover data of 2015.
The coarse-resolution reanalysis data was needed to obtained initial and boundary conditions which were used to drive the model. The mechanism of WRF is finite difference theory and is to simulate weather conditions by applying dynamics and thermodynamics theory in an atmospheric physical environment [32,33,43,44].
The simulation is performed using WRF v3.9.1, which was set up in three nested domains (D01, D02 and D03, see Figure 3, the upper half) with horizontal grid spacing of 25 km, 5 km and 1 km, respectively. The outermost mesh (D01) center has a latitude and longitude of 37.735°N and 82.914°E. Table 2 shows a summary of WRF configuration and physical parameterization. The simulation time is 22 May 2015-7 October 2015, which is the growing season of the crops in Hotan Oasis. To stabilize the model, the first 48 h were considered the spin-up time. The mode was set to output results every hour.  Cover map which has been reclassified and integrated according to USGS classification system; the lower half displays the distribution of test sites, as the figure shows more refined details than the above, and the color and settings of the legend are slightly different.

SEBS Model
The SEBS model is based on the surface energy balance equation in which the surface energy at a certain moment can be expressed as: where is the net radiation, is the soil heat flux, is the turbulent sensible heat flux, and is the turbulent latent heat flux ( is the latent heat of vaporization and E is the actual ET).
The equation to calculate the net radiation is given by: where is the albedo, is the downward solar radiation, is the downward longwave radiation, is the emissivity of the surface, is the Stefan-Bolzmann constant, and is the surface temperature.
The equation to calculate soil heat flux is parameterized as: in which it is assumed that the ratio of soil heat flux to net radiation is Г = 0.05 for full vegetation canopy [45] and Г = 0.315 for barren [46]. An interpolation is then performed between these limiting cases using the fractional canopy coverage, .
In SEBS algorithm, the energy balance computation at limiting conditions was conducted to derive the relative evaporation (Λ ). Under the dry-limit, the latent heat becomes zero due to the limitation of soil moisture and the sensible heat flux is at its maximum value. The expression follows as: Under the wet-limit, where the evaporation takes place at potential rate (λ ), the sensible heat flux takes its minimum value, , i.e., The relative evaporation can be derived as: Finally, the evaporative fraction (Λ) which is the energy used for the ET process is given by: The difference between the SEBS model and others lies in putting forward the concept of evaporative fraction. After the evaporative fraction is determined, the formula of daily average ET can be deduced as follows: where Λ is the daily evaporative fraction and is the density of water (kg/ ).

Evaluation of WRF Simulation and ET
First, WRF-simulated wind speed, air temperature, relative humidity, and surface pressure are evaluated by comparing this data with Hotan meteorological station measurements. Since the station's daily observation is obtained by averaging the observation of four time slots (02:00, 08:00, 14:00, 20:00) of the day, the simulated values is processed in the same way. Root mean square error (RMSE), BIAS , and coefficient of determination ( ) between simulated values and station measurements are adopted to evaluate the WRF simulation [47]. The calculation of these statistical indexes are as follows: where is No. day simulation and is No. day measurement. and are the mean of simulations and measurements respectively, and is the number of days. The ETa-WRF (ET calculated using the meteorological data simulated from the WRF) is analyzed and evaluated in the following two aspects: (i) the temporal variation of ET in different land covers; and (ii) the improvement of ETa-WRF compared to ETa-STA (ET calculated from meteorological data interpolated from measurements at met stations) for different underlying surfaces.
The temporal variation of ETa-WRF is analyzed over homogeneous pixels which is chosen based on a field-proven high-resolution (30 m) land cover map. As the WRF simulation has a spatial resolution of 1 km, the 30 m land cover map is further transformed into a land-cover diversity map with pixels of 1 km × 1 km according to the fraction of the dominant land cover [48]. Then one homogenous 1 km pixel of each land cover is selected as a test site (see Figure 3, the lower half). In order to explore the WRF simulation in estimation of ET over complex terrain, test sites are selected both in a mountainous region and a plain area. The test site of the mountainous region is displayed in the DEM map of Figure 1. All test sites are selected at or near the field survey points to ensure accuracy.

Performance of WRF-Simulated Hydro-Meteorological Variables
The correlation, RMSE, and BIAS between WRF-simulated wind speed, air temperature, relative humidity and surface pressure compared with ground measured are shown in Figure 4. The WRFsimulated air temperature, relative humidity and surface pressure had good correlation with the measured value, with coefficient of determination (R 2 ) of 0.8242, 0.8089 and 0.8915, respectively. However, the WRF-simulated wind speed had a lower coefficient of determination (0.628), and was overestimated compared to site measured with bias of 1.0299 m/s. Earlier studies found that wind speed is one of the most difficult variables to simulate, and the simulated values are usually higher than the observation [49][50][51]. The simulated relative humidity is underestimated compared to site measured with bias of −10.1912%.

Spatial Distribution
The 1 km daily ET of the study area was estimated using the SEBS model with WRF-simulated gridded meteorological data as input, covering the whole process from sowing to harvesting in 2015. The ET of each month from May to September of 2015 is shown in Figure 5.
In general, ETa-WRF captured the dynamic variation of ET in the study area. At the beginning of the growing season (May), ET is relatively low and shows uneven distribution (Figure 5a), as the crops were generally in the sowing period. The water bodies have the highest ET in May. In June, as the leaf is spreading and temperature is increasing, ET (Figure 5b) increased obviously compared with that in May. By July, the vegetation was in full growth, therefore ET reached the maximum (Figure 5c), showing a homogeneous distribution over the vegetation area. ET began to decrease in August (Figure 5d) due to the leaf being out of vegetation and the successive ripening of crops, and then decreased with certain areas close to 0 as the crops were being harvested in September ( Figure  5e), with the water bodies still having the highest ET.  Figure 6 shows the temporal variation of ETa-WRF for each land cover over the growing season of 2015. It is clear that the barren plain has the lowest ET (close to 0) in all land cover throughout the growing season. The barren mountain has slightly higher ET than the barren plain, and with small temporal variation. ET of open shrublands and grasslands has similar temporal variation and values to that of urban and built-up lands. One of the possible reasons for this phenomenon could be that open shrublands and grasslands have sparse vegetation coverage, and so do the urban and built-up lands. ET of croplands has much dynamic temporal variation during the growing season; two possible reasons are phenology variation of dense vegetation coverages and irrigation. Among all the land covers, paddy and water bodies have the highest ET. At the beginning of the growing season, ET of paddy was even higher than that of water bodies, and one possible cause may be that the paddy was in its transplanting stage and the field was fully covered by water, so ET of paddy may be the combined contribution of crops and water. The ET of paddy began to decline gradually and was eventually lower than that of the water bodies at the beginning of August, and one possible reason is the drainage of the paddy fields. Water bodies had high ET during the whole period and had small temporal variation with a slight decreasing trend in the late summer.

Temporal Variation
In summary, ETa-WRF could reflect the temporal dynamics during the growing season of vegetation, also can capture the difference of ET in different land cover.

Discussion
Improvement of ET estimation using WRF simulation compared to using station measurement will be discussed in this section.

Improvement of ETa-WRF Among Different Land Covers Over Homogeneous Pixel
Comparison between ETa-WRF and the corresponding ETa-STA is shown in Figure 7. Figure 8 displays the temporal variation of normalized difference vegetation index (NDVI) of each land cover. Earlier studies found that in arid or semi-arid areas without irrigation, the ET is mainly caused by the transpiration of vegetation, and basically maintains a positive correlation with the vegetation index, such as NDVI [52][53][54][55][56].  The comparison of the two figures shows that the ETa-WRF of the barren plain is lower than that of ETa-STA in the whole growth season, and ETa-STA shows an abnormally high peak in July compared with other months, while the change of ETa-WRF is relatively small. It can be seen from Figure 8 that the NDVI measures of the barren plain and barren mountain are close to 0 and have relatively small temporal variation from May to September, significantly lower than that of open shrublands and grasslands and urban and built-up lands. However, ETa-STA results showed that the ET of the barren mountain was higher than that of the open shrublands and grasslands, and even higher than that of the urban and built-up lands and croplands in certain months (Figure 7a). The ETa-STA of the barren mountain was generally several orders of magnitude higher than ETa-WRF over the whole growing season. However, the selected barren mountain test site has rather low vegetation coverage and is not covered by perennial snow or glaciers, which are unlikely to have high ET. Figure 7a shows that the ETa-STA of open shrublands and grasslands fluctuated greatly in the growing season, except that it is higher than ETa-WRF in July, while it is significantly lower than ETa-WRF in other months. As can be seen from Figure 8, the NDVI of open shrublands and grasslands is always close to or even higher than that of built-up lands, while the ETa-STA of open shrublands and grasslands is much smaller than that of built-up lands except for July, which obviously does not conform to the changing characteristics of NDVI, while the ETa-WRF performs better in line with the general trend of NDVI. The ETa-STA results of the barren mountain, open shrublands and grasslands test sites did not agree with the conclusions drawn by earlier studies [52][53][54][55][56]. We acknowledge that the comparison between the simulation of WRF and the measurements does show some deviation and uncertainty, but compared with the pure mathematical interpolation based on the site data, WRF simulation can consider the heterogeneity of the underlying surface better. In future research, we will adjust the physical parameterization scheme related to the WRF model to explore how to obtain more accurate simulation results.

Macroscopical Improvement of ETa-WRF Over the Underlying Surface
The 7.20-7.27/2015 is the time when vegetation was in full growth and had relatively higher ET, and it was selected as the date for the comparison between ETa-WRF and ETa-STA. It was obvious that in barren areas with low NDVI (less than 0.08) in the south of the study area (Figure 9), ETa-STA (Figure 10a) still has relatively high values (mostly 4-5 mm/day), which is even close to the ET of open shrublands and grasslands and urban and built-up lands. Especially in the barren mountainous area where the NDVI level is also low, the ET was even higher than that of croplands. While ETa-WRF (Figure 10b) improved, the barren area has ET close to 0. This is because the meteorological parameters simulated by WRF can fully consider the heterogeneity of the underlying surface, therefore providing better ET estimation over complex underlying surfaces.

Conclusions
Accurate grid meteorological data are crucial for the accurate estimation of ET using the SEBS model. This paper used the meteorological model WRF to simulate high-spatial resolution gridded meteorological data with full consideration of the underlying surface heterogeneity, which was further used as input for SEBS model to generate ET of the Hotan Oasis in China, from May to September 2015.
The WRF-simulated wind speed, temperature, relative humidity and air pressure correlated well with site-measured values with a RMSE of 1.1244 m/s, 1.9552 ℃, 11.8113%, 1.5337 hPa and R 2 of 0.628, 0.8242, 0.8089 and 0.8915. This evaluation indicated that WRF-simulated meteorological parameters could be used for the estimation of ET.
Compared with ET estimated using meteorological data interpolated from station measurements as input, ETa-WRF could better reflect the ET difference between different land covers and capture the vegetation growing trend, especially in areas with sparse vegetation.
Reliable gridded meteorological data at a regional scale can be obtained considering the underlying surface heterogeneity using a WRF model, and regional ET which can better capture the underlying surface can be generated, which can possibly provide a reliable reference for precise management of water resources and rational formulation of macrocontrol policies. Future research will focus on improving WRF-simulated wind speed and relative humidity variables for ET estimation, and conducting quantitatively evaluation of ET.