Long-term spatial distributions and trends of the latent heat fluxes over the global cropland ecosystem using multiple satellite-based models

Estimating cropland latent heat flux (LE) from continental to global scales is vital to modeling crop production and managing water resources. Over the past several decades, numerous LE models were developed, such as the moderate resolution imaging spectroradiometer LE (MOD16) algorithm, revised remote sensing-based Penman–Monteith LE algorithm (RRS), the Priestley–Taylor LE algorithm of the Jet Propulsion Laboratory (PT-JPL) and the modified satellite-based Priestley-Taylor LE algorithm (MS-PT). However, these LE models have not been directly compared over the global cropland ecosystem using various algorithms. In this study, we evaluated the performances of these four LE models using 34 eddy covariance (EC) sites. The results showed that mean annual LE for cropland varied from 33.49 to 58.97 W/m2 among the four models. The interannual LE slightly increased during 1982–2009 across the global cropland ecosystem. All models had acceptable performances with the coefficient of determination (R2) ranging from 0.4 to 0.7 and a root mean squared error (RMSE) of approximately 35 W/m2. MS-PT had good overall performance across the cropland ecosystem with the highest R2, lowest RMSE and a relatively low bias. The reduced performances of MOD16 and RRS, with R2 ranging from 0.4 to 0.6 and RMSEs from 30 to 39 W/m2, might be attributed to empirical parameters in the structure algorithms and calibrated coefficients.


Introduction
Latent heat flux (LE) plays a key role in the energy and water cycles in agricultural ecosystems [1,2]. a large number of studies have shown that LE is a vital variable for developing precise irrigation scheduling and enhancing water use efficiency in agricultural production due to the close relationship between soil water depletion and the rate of evapotranspiration (equivalent production [22]. Therefore, long-term spatial distributions and trends of LE based on different LE models in cropland ecosystems at the global scale under various climatic conditions are still needed. In the present study, we evaluated four LE models, including MODIS Algorithm (MOD16), Revised Remote Sensing-Based Penman-Monteith LE Algorithm (RRS), Priestley-Taylor Algorithm (PT-JPL) and Modified Satellite-Based Priestley-Taylor Algorithm (MS-PT), using the LE measurements from thirty-four eddy covariance sites over cropland ecosystems, and results from temporal trend and spatial distribution of the estimated cropland LE by satellitebased observations. The objectives of this study were as follows: (1) compare LE algorithm performances over cropland ecosystems based on 34 flux tower sites provided by the EC towers network and (2) analyze the differences in the temporal and spatial patterns of LE over the global cropland ecosystem during 1982~2010.

Ground EC data
Ground data collected from 34 EC towers of FLUXNET, Asiaflux, Ameriflux were used for the analysis. We used the FLUXNET data provided by FLUXNET community (http://fluxnet. fluxdata.org/data/la-thuile-dataset/). The radiation data included continuous series of halfhourly incident solar radiation (Rs, W/m 2 ), surface net radiation (Rn, W/m 2 ), sensible heat flux (H, W/m 2 ), latent heat flux (LE, W/m 2 ) and ground heat flux (G, W/m 2 ). The meteorological data included continuous series of half-hourly incident relative humidity (RH), air temperature (Ta, K), diurnal air temperature range (ΔT, K), wind speed (Ws, m/s) and vapor pressure (e, Pa). The meteorological data and radiation data provided input variables, while groundmeasured LE were used as the validation [24,25]. We converted the half-hourly data into daily data. Quality control and gap-filling techniques were utilized to strengthen the quality of the data sets. The tower sites covered main cropland ecosystems over Europe, Asia and America. The data were selected to cover the period of 2000-2009, with at least one year of data available (Table 1).
Although it is ideal to measure LE using the EC technique, the sum of H and LE was overall less than the available energy. EC sites underestimated LE, and an H flux was generally observed over most EC sites in Europe and North America [25]. Twine et al. (2000) [26] also reported the energy closure problem based on the Southern Great Plains 1997 Hydrology Experiment; the closure issue became more significant upon consideration of the long-term water balance. Therefore, we used the method provided by Twine et al. to eliminate the effect of the unclosed energy problem.
where LE c is the corrected LE flux; H o and LE o are the uncorrected sensible heat flux and LE flux, respectively. NDVI and FPAR values for the EC sites were derived from AVHRR products. The AVHRR NDVI data were downloaded using the ftp link (http://ecocast.arc.nasa.gov/data/pub/gimms/). LAI values were obtained from the GLASS product. GIMMS NDVI data at a 0.833-degree spatial resolution were used for model verification at the EC sites. FPAR data were calculated using the method proposed by Myneni and Williams (1994) [27]; this method was verified in various ecosystems [28]. The GLASS data can be found on the web site.

Satellite and reanalysis datasets
When estimating LE at regional scale for the physically based LE algorithms, we obtained AVHRR products, including bimonthly GIMMS NDVI and UMD Land Cover data with 1 km spatial resolution. We used the MODIS products with 0.05-degree spatial resolution, including annual land-cover type (MCD12C1). We also used the monthly LAI product (GLASS01B01), with 0.05-degree spatial resolution, from the GLASS datasets. The digital elevation model (DEM) product was obtained from the SRTM30 dataset. The DEM data can be found on the web site (https://dds.cr.usgs.gov/srtm/version2_1/SRTM30/). We used MERRA datasets for 1982-2010 to obtain daily radiation and meteorological data, including Rs, Rn, RH, Ta, ΔT and Ws. The MERRA dataset is a NASA reanalysis for the satellite era that uses a new version of the Goddard Earth Observing System Data Assimilation System Version 5 (GEOS-5). It has been reported that the MERRA dataset produced comparable results for global and energy cycle research [29]. We obtained MERRA data from the Goddard Earth Science Data and Information Services Center web site (https://gmao.gsfc.nasa.gov/reanalysis/MERRA/). Considering the different temporal and spatial resolutions of the input data, including daily MERRA (½ degree latitude × degree longitude spatial resolution), bimonthly GIMMS NDVI (0.833-degree spatial resolution), UMD Land Cover data (1-km spatial resolution), and monthly MODIS (0.05 degrees), we averaged the daily MERRA data and bimonthly GIMMS NDVI to yield monthly data sets. For the land cover classification, we used the UMD Land Cover data before 2001 and MCD12C1 after that. We spatially interpolated coarse resolution MERRA data and GIMMS NDVI to the resolution of 0.05 degrees using a spatial linear interpolation scheme. Adjustments were made to ensure all inputs parameters in LE models had consistent temporal and spatial resolutions.  [12] as follows: (1) improving the calculation of the vegetation cover fraction with FPAR; (2) modifying the estimation of aerodynamic, boundary-layer, and canopy resistance; (3) calculating the soil heat flux using available energy; (4) calculating LE based on daytime and nighttime LE; (5) calculating evaporation from the wet canopy surface and transpiration from the dry surface; (6) calculating evaporation from the saturated soil surface and evaporation from the moist soil surface. Previous validation results showed that the MOD16 algorithm successfully improved LE estimation using 46 AmeriFlux sites [31]. The MOD16 algorithm was applied to produce MODIS terrestrial LE using the MODIS products and GMAO daily meteorological dataset.

Revised remote sensing-based Penman-Monteith LE algorithm (RRS).
Previous studies showed that high air temperature significantly decreases leaf stomata, causing structural defects [32]. Therefore, Yuan et al. (2010) [33] revised the equations dealing with the temperature constraint by accounting for stomata conductance and energy allocation between the vegetation canopy and soil surface [34]. The Beer-Lambert law was also used to exponentially partition the net radiation between the canopy and soil surface. Model parameters across different vegetation types were set to be consistent to reduce the effects of misclassification. Validation results based on 23 EC flux tower sites in China provided by Chen et al. (2014) [35] showed that the RRS algorithm improved LE estimation compared to the MOD16 algorithm.  [14] created a bio-meteorological method to calculate actual ET (ET, equivalent to LE). The inputs of this model included net radiation (Rn), normalized difference vegetation index (NDVI), soil adjusted vegetation index (SAVI), maximum air temperature (Tmax), and water vapor pressure (ea). The total ET includes the canopy transpiration and evaporation from soil and rain water intercepted by the canopy. Each component flux is determined by eco-physiological factors or conditions to reduce potential LE to actual LE using the Priestley-Taylor equation.

Modified satellite-based Priestley-Taylor algorithm (MS-PT).
To overcome the computational complexities of aerodynamic resistance parameters, Yao et al. (2013) [36] revised the Priestley and Taylor equation using an index of soil water deficit calculated using the apparent thermal inertia and diurnal temperature range to constrain soil evaporation and using the revised linear two-source model (N95) to calculate vegetation transpiration. For canopy interception, the relative surface wetness index was modified to the fourth power of the soil moisture constraint. The model only required four inputs, including Rn, Ta, DT, and NDVI. Validation by Yao et al. (2015) [37] showed that the MS-PT algorithm had good performance over global forest biomes. Feng et al. (2015) [38] found that the MS-PT algorithm improved LE estimation over global grassland, shrubland and savanna.

Statistical analysis
We compared these models at the global scale using the modern-era retrospective analysis for research and applications (MERRA) dataset, the Normalized Difference Vegetation Index (NDVI) and Fractional Photosynthetically Active Radiation (FPAR) from the moderate resolution imaging spectroradiometer (MODIS) products, the leaf area index (LAI) from the global land surface satellite (GLASS) products, and the digital elevation model (DEM) from the shuttle radar topography mission (SRTM30) dataset to generate LE at the regional scale during 1982~2010.
The four models were first validated at 34 EC towers. We calculated the daily LE using the meteorological variables collected from EC sites and the NDVI and LAI from the satellite data (GIMMS NDVI and GLASS01A01). The vegetation index at the EC sites was calculated by computing the average of the surrounding pixels. Before the comparison, we calculated the coefficient of determination (R 2 ), RMSE (root mean square error), and average bias for estimated LE to evaluate the model performances. We used Taylor diagrams [39] to analyze the performance of the four LE models. The Taylor diagrams provided a summary of how well patterns matched observations in terms of statistical indexes. These statistical indexes included standard deviations, centered root mean square error and their correlation. The radial distance from the origin indicated the standard deviations; the cosine of the azimuth angle showed correlation. The standard deviation of simulated patterns is shown to be proportional to the radial distance from the observed point to the patterns.
For regional estimates of LE, we applied simple spatially interpolated method to the MERRA data and AVHRR data to match the MODIS spatial resolution. Daily MERRA data and half monthly AVHRR NDVI data were combined to yield monthly data. After the interpolation process, we calculated LE over a cropland biome from 1982 to 2010 based on four models. The bias computed by the LE values based on one model minus that of the others was also calculated.
To estimate the long trend variability of LE, we calculated annual LE using four models and applied a linear trend analysis [40] using a linear model (y = ax + b), where t, a and b are the time, slope and intercept of the regression line, respectively. Student's t-test was used to analyze and classify the significance of the trend at the 1% and 5% level.

Models evaluation
All four models explained 50%~76% of the LE variability over all sites (Fig 1). The MS-PT showed the highest R 2 (0.76). The averaged bias for all four models ranged from -14.9 to 2.38 W/m 2 . Generally, MOD16 showed a slight positive bias and the remaining methods had negative biases. RRS had the highest negative bias (14.43 W/m 2 ) compared with the other models. The RMSEs for the four models ranged from 23 to 35 W/m 2 . The MS-PT showed the lowest RMSE (23.74 W/m 2 ) compared with the other methods.
Generally, MS-PT had good overall performance for cropland (Fig 2). This result was consistent with the previous study [36]. In

Spatial patterns of LE
All models generally correlated well with ground observations of daily LE at most sites [17], with R 2 ranging from 0.4 to 0.9, biases ranging from 15 to 50 W/m 2 , and RMSE ranging from  Generally, the four models showed similar spatial patterns of LE, with high LE in warm and humid areas and low LE in arid and cold areas (Fig 4). However, some substantial differences in LE estimates existed among the models. The average annual LE from MS-PT was 36.26 W/ m 2 , higher than RRS (33.49 W/m 2 ) and MOD16 (35.63 W/m2) and lower than PT-JPL (39.73 W/m 2 ). East of the Ural Mountains near Chelyabinsk, all models showed low LE estimates ranging from 20~50 W/m 2 (Fig 5). In the northwest part of India, the LE was relatively low (30~40 W/m 2 ), but in the southeast part, the LE was approximately 50~80 W/m 2 . This might be attributed to the spatial heterogeneity of available energy and humidity. The largest difference in LE estimates occurred between the eastern and western parts of Europe. The LE values in the western part of Europe (40~60 W/m 2 ) were relatively higher than those in the eastern part (20~40 W/m 2 ). This was likely because the water availability with large-scale atmospheric circulation influenced LE distribution. All models showed a low LE estimate (20~40 W/m 2 ) in the northern part of North America and high LE estimate (40~70 W/m 2 ) in the southern part.
We further analyzed the spatial variability of magnitude of inter-model differences in the LE estimates in different parts of the world (Fig 6). MOD16 showed higher LE than MS-PT in the northwest part of North America (approximately 0~10 W/m 2 ). However, MOD16   in Northern India, the LE estimates by MS-PT were higher than those by RRS, with the difference ranging from 30 to 70 W/m 2 . These results indicated that the differences of four LE models are not spatial homogeneous. This might be attributed to different sensitivity of input parameters of LE models in different climatic regimes [41].

Annual mean and long-term trends
All models showed strong seasonality in LE, and they captured the magnitudes and seasonal variations of LE (Fig 7). The results indicated that the LE produced by the four models had similar monthly patterns, with maxima in the summer (62~85 W/m 2 ) and minima in the winter (11~42 W/m 2 ). However, there were also substantial differences among models. The PT-JPL predicted the highest LE compared with the other models in spring and summer. The difference between MS-PT and PT-JPL mainly occurred in spring, when the temperature changed significantly (Fig 8). MOD16 and RRS almost showed similar magnitudes, with a slight difference occurring during summer months.  Latent heat fluxes estimation over the global cropland ecosystem using multiple satellite-based models Similar slightly annual increasing trends in LE estimates from 1982 to 2010 were produced by the four models (Fig 9). The LE estimates predicated by MOD16, RRS, MS-PT and PT-JPL were in the range of 31~42 W/m 2 .
The long-term changes in LE differed among the four models and had different spatial patterns among the four models (Fig 10). All models showed strongly increasing LE in Asia. Specifically, approximately 46, 49 and 50% of area in Asia were dominated by the increasing trend in RRS, MOD16 and MS-PT, respectively. However, only 29% of the area was dominated by the increasing trend in PT-JPL. Generally, all models showed an increasing trend in Northern India (Fig 11). However, in the northeast part of China and east of the Ural Mountains, MOD16, RRS and MS-PT showed an increasing trend, while the opposite was found for PT-JPL, with a slightly decreasing trend. In the Far East Areas of Russia, MOD16, RRS, and PT-JPL simulated strongly positive LE trends; however, MS-PT showed negative LE trends. In Europe, increasing trends in LE accounted for approximately 61, 66 and 74%, according to MOD16 and MS-PT. However, only 45% and 14% of the area was dominated by the increasing trends in RRS and PT-JPL. Consistent LE increasing trends were found among the four models for Western Europe. PT-JPL simulated a slightly decreasing trend in Eastern Europe, while increasing ET trends were indicated by the other four models. In North America, increasing trends in LE accounted for approximately 16, 27, 33 and 44%, according to the PT-JPL, RRS, MOD16 and MS-PT, respectively. In the central south of Canada, all models showed an increasing trend. In Iowa, RRS and PT-JPL showed a decreasing trend. In Australia, no substantial changes were found by the four models. In South America, only MOD16 showed increasing trends in LE, which covers approximately 47% of the total area in South America. There is an obvious Latent heat fluxes estimation over the global cropland ecosystem using multiple satellite-based models difference of annual trend between PT-JPL and the other three models. These differences have connection with climate regime and local climate condition via meteorological parameters such as relative humidity and diurnal temperature range. These meteorological parameters have closely relationship with constraint parameters in these algorithms. The MOD16 and RRS generally adopt these meteorological parameters to calculate canopy resistance and soil resistance but PT-JPL and MS-PT adopt these parameters to calculate scale factors to downscale potential LE to actual LE. Relative humidity is used as input in PT-JPL but no used in MS-PT. So the variations of meteorological parameters related to climate regime cause the difference of annual trend cause the difference of annual trend among the models.

Model performance
Ground measurements and estimated LE were in good agreement over cropland sites, suggesting that it can be considered as a suitable method for evaluating the cropland LE over a variety Latent heat fluxes estimation over the global cropland ecosystem using multiple satellite-based models of climate conditions. The reason might be reduced the number of input parameters, which only included net radiation, NDVI, air temperature and diurnal temperature range. The fewer input parameters eliminated the effect caused by input uncertainties. The MS-PT also avoided the complex calculations of aerodynamic and surface resistance. The good performance of MS-PT was similar to that of our previous studies [36][37][38].   [42] also found that the MS-PT model improved the LE estimation compared with the 40 EC sites data.
Both MS-PT and PT-JPL were developed based on a simple LE equation for radiationand temperature-based equilibrium by Priestley and Taylor [13]. However, the different performances of the two models might be caused by the different calculations of the ecophysiological constraint functions. Fisher et al. (2008) [14] calculated relative surface wetness by squaring relative humidity (RH) twice. However, for MS-PT, the relative surface wetness was the fourth power of the soil moisture constraint. The soil moisture constraint in MS-PT was the exponential algorithm of the Apparent Thermal Inertia (ATI). The ATI was simply estimated as the inverted diurnal temperature range. For PT-JPL, the soil moisture constraint was the exponential algorithm of RH and saturation VPD. Therefore, appreciable differences might occur due to the calculation of actual soil surface evaporation and actual vegetation transpiration.
The results of MOD16 were generally concordant with the results of RRS. However, there were small discrepancies in the LE estimates by the two models. This difference might be attributed to the difference in model structures related to the calculations of surface resistance. Both LE models used the same controlling factors from VPD and minimum temperature (Tmin) on stomatal conductance. Mu et al. (2011) [12] calibrated VPD and Tmin for each biome. However, these variables were different in RRS. Both models used different ways to partition the net radiation between the canopy and soil surface. For MOD16, the fraction of absorbed photosynthetically active radiation was a proxy of the vegetation cover fraction. However, Fc was calculated using the Beer-Lambert law with LAI. All these factors led to different results of LE predicted by the two models.

Spatial and temproal differences
The satellite estimates of LE predicted by these LE models generally displayed a similar distribution, with large values at low latitudes but smaller values at high latitudes. This was probably due to uneven distribution of the solar radiation. In North America, LE values in the northern part were lower than those in the southern part. A similar pattern occurred in East Asia and North India. Vinukollu et al. (2011) [43] calculated LE using three models, including the Surface Energy Balance System (SEBS), MOD16 and a modified Priestley-Taylor approach (PT-JPL), and found a similar distribution pattern. Zhang et al. (2010) [44] quantified global terrestrial ET using the modified Penman-Monteith approach from 1983 to 2006 and found strong regional variations and latitudinal gradients corresponding to global climate patterns. The results also showed that coastal areas had higher LE than land areas, such as Europe and South Australia.
Inter-model differences generally ranged from -30 to 10 W/m 2 based on MERRA data for 1982~2010. Although these LE models were developed based on a physical mechanistic equation, different results might be attributed to the model structure and sensitivity of input variables, which were related to the regional climatic environment. Previous studies also supported our conclusion and highlighted the model structure and the uncertainty of input data. In this study, global decadal cropland LE was calculated using four process-based LE models. We found an increasing trend of global average annual LE and different variations of LE in specific regions over the cropland ecosystem. Numerous studies showed increasing transpiration due to the increase of the growing seasons in recent decades [46][47][48][49][50]. The increasing trend of cropland LE might be attributed to increasing irrigation water diversions. McLaughlin et al. found an increasing trend of irrigation water from 1960 to 2010. The regional long-term trend of LE generally varied with changes of anthropogenic activities and regional climate. Sacks et al. (2011) [51] analyzed corn and soybean data over the US corn belt from 1981~2005 and found that the annual-average LE increase was 0.50 W/m 2 , or 6.3 mm/year, over the 25 years. Riediger et al. (2014) [52] used a framework of the comprehensive research cluster Nachhaltiges Land management in Norddeutschen Tiefland and predicted increasing scenarios in ET at Uelzen in Germany.

Uncertianties analysis
Uncertainties in predicting LE were considered to be largely induced by the forcing data [53][54][55]. Liu et al. (2015) [54] analyzed the impact of input data uncertainties on (ET) using the Climate Research Unit (CRU) TS3.1 data, the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim Reanalysis, the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis data, the Global Meteorological Forcing Data set for land surface modeling from Princeton University, and MERRA data for 1979-2008 in Northern China. They found that MERRA displays the lowest correlation with measurements for temperature, wind and relative humidity. These uncertainties might affect the accuracy of the final calculation. Zhao et al. (2013) [56] also reported that MERRA surface solar radiation has an average bias error of +20.2 W/m 2 at different time scales based on American FLUXNET sites. NDVI or LAI were often used as a surrogate of vegetation moisture and for the calculation of vegetation coverage. However, the uncertainties in NDVI or LAI led to erroneous subsequent calculation of canopy conductance.
In addition to the uncertainties in forcing data, another source of uncertainty might stem from the scale mismatch between the remote sensing data and the EC observations. McCabe et al. (2006) [57] used three different remote sensing data from Landsat (60 m), ASTER (90 m), and MODIS (1020 m) to study the effects of sensor resolution on land surface flux estimates. These scholars found that MODIS-based flux retrievals could not effectively characterize the different latent heat flux patterns between soybean and corn. Kustas et al. (2004) [58] also found a lack of spatial variability in the coarse resolution data-based LE retrievals compared with different pixel resolutions of remote sensing inputs.

Conclusions
In this study, four LE satellite-based LE algorithms, including MOD16, RRS, PT-JPL and MS-PT, were evaluated over a cropland ecosystem based on 34 EC flux towers sites. We compared these LE model performances with ground LE measurements and satellite forcing data. All models had different assumptions, data requirements, empirical model parameters and model structures ranging from complex Penman-Monteith and energy balance schemes to simpler Priestley-Taylor approaches. The results showed that all the LE models produced acceptable results for cropland ecosystems. MS-PT produced the best performance with respect to LE prediction, with a higher R 2 , a lower bias and a lower RMSE. MOD16 and RRS showed the lowest LE retrievals compared with the other models.
All algorithms had high LE at low latitudes and low LE at high latitudes. The results showed that the magnitudes and seasonal variations of LE were well characterized by the four models.
A slightly increasing trend in the annual average LE from 1982 to 2010 was also observed. Considering that uncertainties still existed, the limitations associated with the model comparison, as analyzed above, might also cause modeling errors. Intercomparison of models over a range of diverse bioclimates and longer periods are recommended for future studies.