Modeling and Partitioning of Regional Evapotranspiration Using a Satellite-Driven Water-Carbon Coupling Model

The modeling and partitioning of regional evapotranspiration (ET) are key issues in global hydrological and ecological research. We incorporated a stomatal conductance model and a light-use efficiency-based gross primary productivity (GPP) model into the Shuttleworth–Wallace model to develop a simplified carbon-water coupling model, SWH, for estimating ET using meteorological and remote sensing data. To enable regional application of the SWH model, we optimized key parameters with measurements from global eddy covariance (EC) tower sites. In addition, we estimated soil water content with the principle of the bucket system. The model prediction of ET agreed well with the estimates obtained with the EC measurements, with an average R2 of 0.77 and a root mean square error of 0.72 mm·day−1. The model performance was generally better for woody ecosystems than herbaceous ecosystems. Finally, the spatial patterns of ET and relevant model outputs (i.e., GPP, water-use efficiency and the ratio of soil water evaporation to ET) in China with the model simulations were assessed.


Introduction
Evapotranspiration (ET) is a key process of the ecosystem water cycle and energy balance and is closely linked to ecosystem productivity [1,2].Thus, detailed and precise knowledge of regional ET is important to obtain a better understanding of the global carbon and water cycles [2][3][4].Simulating ET at regional scales with high accuracy remains challenging despite the development of many relevant models [5,6].For example, a comparison of 15 model simulations from the Global Soil Wetness Project-2 (GSWP-2) revealed that the annual mean global land surface ET ranged from 272 to 441 mm•year −1 , indicating large discrepancies among the models [7].A comparison by Chen et al. [6] of eight ET models revealed a range of annual mean ET values from 535 to 852 mm in China, confirming large model-to-model differences.Thus, further improvements in regional ET simulation and partitioning are needed.
The Penman-Monteith model (P-M model) [8] and the Shuttleworth-Wallace model (S-W model) [9] are favored because they capture the physical process understanding to the best of our current knowledge.At the site scale, since most input variables can be measured directly or estimated from meteorological and biotic measurements, these two models are widely used [10][11][12][13][14].At the regional scale, however, the P-M and S-W models are seldom used because some key variables, particularly the canopy stomatal conductance, are difficult to derive using remote sensing data.Cleugh et al. [15] developed a remote sensing ET algorithm using the P-M scheme in which canopy stomatal conductance was estimated as a function of the Normalized Difference Vegetation Index (NDVI) and Leaf Area Index (LAI).Based on Cleugh's algorithm, Mu et al. [16,17] developed a new version of a global remote sensing ET model.Compared with measurements from eddy covariance (EC) systems, Mu's model generally illustrated good performance for many FLUXNET sites [16].
The S-W model is a two-source model developed from the P-M model.This model simulates the water vapor flux from plants (i.e., plant transpiration) and soil (i.e., soil evaporation) separately.The S-W model has been widely used because of its simple and accurate consideration of hydrological processes and good performance [11][12][13][14]18,19].Most studies indicate that the performance of S-W for diverse ecosystems is superior to that of other ET models, particularly in terms of partitioning [10,20].If the variables and parameters are available, the S-W model could be promising for modeling and partitioning ET at the regional scale.
Using the Ball-Berry model [21] to estimate canopy stomatal resistance and a soil moisture function to estimate soil surface conductance, we previously deployed the S-W model to simulate ET with high accuracy in four grassland ecosystems [14].By further adopting a light-use efficiency-based gross primary productivity (GPP) model, we developed a model, SWH, for estimating and partitioning ET using meteorological data and remote sensing products [22].Compared to EC measurements, the SWH model produced satisfactory estimates of ET and GPP in a temperate forest and alpine grassland.However, two features of the current version of the SWH model prevent its application at the regional scale.First, key parameters that were optimized at individual sites prior to model development are not available for regional application.Second, soil water content, an input variable, is measured directly with soil moisture sensors, which are difficult to access at the regional scale.Therefore, to apply the SWH at the regional scale, parameterization and soil water content estimation must be addressed.Accordingly, the objectives of this study are: (1) to optimize key parameters of the model by using the in situ measurements from FLUXNET towers across the globe; (2) to estimate soil water content with a bucket system model; (3) to evaluate the model performances on ET and GPP simulation in different biome types with the FLUXNET measurements; (4) using China as a case study region, to assess the spatial patterns of ET, GPP, the fraction of soil evaporation (E/ET) and water use efficiency (WUE, GPP/ET) based on the outputs of the new model.

Model Description
The SWH model simulates plant transpiration and soil water evaporation separately with algorithms based on the Penman-Monteith equation [8].A full description of the SWH model is available in Appendix A and in Hu et al. [14].Five resistances are necessary to estimate ET with the SWH model (Figure 1, Appendix A), i.e., soil surface resistance (r ss ), canopy stomatal resistance (r sc , the reverse of canopy stomatal conductance), the aerodynamic resistances encountered by the water flux leaving leaf lamina (r ac ) or soil surface (r as ) before being incorporated into the mean canopy flow and the transfer resistance between the hypothetical mean canopy flow and the reference height (r aa ).r ss and r sc are the most critical for model performance among the five resistances.r ss was estimated as a function of soil water content [23,24]: where SWC and SWC s are the soil water content and saturated water content in the surface soil (m 3 •m −3 ) and b 1 (s•m −1 ), b 2 and b 3 (s•m −1 ) are empirical constants.b 1 is fixed as 3.5 s•m −1 [23].
Remote Sens. 2017, 9, 54 3 of 20 where SWC and SWCs are the soil water content and saturated water content in the surface soil (m 3 •m −3 ) and b1 (s•m −1 ), b2 and b3 (s•m −1 ) are empirical constants.b1 is fixed as 3.5 s•m −1 [23].Using the scheme of the Biome-Biogeochemical(Biome-BGC model) [24], SWC was estimated using a one-layer bucket model in which soil water storage was estimated as the balance of water input (i.e., precipitation and snowmelt) and output (i.e., ET, runoff and sublimation).Runoff was calculated as the surplus precipitation overpassing the water storage at the point of the saturated and field-holding capacity.Snowmelt and sublimation were calculated as functions of air temperature and solar incident radiation.Thereafter, SWC was calculated using soil water storage and the parameter 'effective soil depth' (d), which is determined by parameterization, explained below.Due to the insensitivity of the model output, the same empirical values of 0.45 and 0.35 m 3 •m −3 were used for saturated water content and field-holding capacity, respectively, for all ecosystems.Details for calculating the SWC are available in the Biome-BGC handbook (http://www.ntsg.umt.edu).
We estimated rsc by introducing the Ball-Berry model in our study [21]: where g0 and a1 are empirical parameters, Pn (µmol•m 2 •s −1 ) is the photosynthetic rate and hs is the leaf surface relative humidity, which was approximated using the air relative humidity.CS is the leaf surface CO2 content.CS should vary seasonally and inter-annually, considering it is insensitive to model simulation; hence, we simplified it as a constant (390 ppm).Note that a varying CO2 must be used if long-term trend analysis is to be conducted.
Pn is a key driving variable in the estimation of rsc.We used GPP, which was estimated using a satellite-based light-use efficiency model, to replace this variable [22]: Using the scheme of the Biome-Biogeochemical(Biome-BGC model) [24], SWC was estimated using a one-layer bucket model in which soil water storage was estimated as the balance of water input (i.e., precipitation and snowmelt) and output (i.e., ET, runoff and sublimation).Runoff was calculated as the surplus precipitation overpassing the water storage at the point of the saturated and field-holding capacity.Snowmelt and sublimation were calculated as functions of air temperature and solar incident radiation.Thereafter, SWC was calculated using soil water storage and the parameter 'effective soil depth' (d), which is determined by parameterization, explained below.Due to the insensitivity of the model output, the same empirical values of 0.45 and 0.35 m 3 •m −3 were used for saturated water content and field-holding capacity, respectively, for all ecosystems.Details for calculating the SWC are available in the Biome-BGC handbook (http://www.ntsg.umt.edu).
We estimated r sc by introducing the Ball-Berry model in our study [21]: where g 0 and a 1 are empirical parameters, P n (µmol•m 2 •s −1 ) is the photosynthetic rate and h s is the leaf surface relative humidity, which was approximated using the air relative humidity.C S is the leaf surface CO 2 content.C S should vary seasonally and inter-annually, considering it is insensitive to model simulation; hence, we simplified it as a constant (390 ppm).Note that a varying CO 2 must be used if long-term trend analysis is to be conducted.
P n is a key driving variable in the estimation of r sc .We used GPP, which was estimated using a satellite-based light-use efficiency model, to replace this variable [22]: where PAR is the incident photosynthetically-active radiation (µmol•m −2 •s −1 ) and FPAR is the fraction of PAR absorbed by the canopy, which was estimated as a function of NDVI [25] (FPAR = 1.24NDVI − 0.168).
At the very beginning, we used the MODIS LAI/FPAR product directly, finding many outliers or no data available at many sites, which hampered the job of parameterization.Therefore, we used NDVI to estimate FPAR and LAI.ε is the light-use efficiency (µmol•CO 2 •µmol −1 PPFD, photosynthetic photon flux density(PPFD)) and is down-regulated by air temperature and the vapor pressure deficit (VPD): 5) where ε max is the apparent quantum yield or maximum light-use efficiency and f (T) and f (VPD) are the downward-regulation scalars for the effects of temperature and VPD (kPa) on the light-use efficiency of vegetation, respectively.T min , T max and T opt are the minimum, maximum and optimum air temperatures ( • C) for photosynthetic activity, respectively.If the air temperature falls below T min or increases above T max , f (T) is set to zero.In this study, T min , T opt and T max were set to 0, 20 and 40 • C, respectively [26].A minimum and maximum VPD, which indicates the start and full constraint of VPD was set as 0.5 and 3.5 kPa, respectively [27].
Note that various schemes are used to down-regulate ε in different light-use efficiency-based GPP models [26,[28][29][30].We compared most schemes and their combinations to yield the solution used in this study, which has a relatively high accuracy and is convenient for use at the regional scale [31].
In the SWH model, LAI and the light extinction coefficient (λ) are needed to estimate the downward solar radiation reaching the soil surface.LAI was estimated as a function of FPAR and λ using Beer's law (Equation (A6) in Appendix A).A plant functional type-specific λ was obtained based on a meta-analysis of global terrestrial ecosystems [32].Soil heat flux (G) is needed to estimate the net radiation to the soil surface (Equation (A5) in Appendix A), which was estimated as a function of net radiation and NDVI (G = 0.1R n (1 − 0.98NDVI 4 ) [33].The model time step was set to 8 days because the satellite products were calculated as 8-day composites.

Site Information
The EC flux data and in situ meteorological measurements from 68 FLUXNET sites (The LaThuile free fair-use Dataset), were extracted from the Level 4 product, which was downloaded directly from the FLUXNET website (Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), 2013).Data for 8 additional sites were obtained from the Chinaflux network.After removing the sites for which there was less than one year of measurements and no net radiation measurements or NDVI data were available, 63 sites remained for analysis in this study (Table S1), including 7 cropland sites, 11 grassland sites, 5 savanna sites, 3 wetland sites, 3 shrubland sites, 7 deciduous broadleaf forest sites (DBF), 5 evergreen broadleaf forest sites (EBF), 19 evergreen needleleaf forest sites (ENF) and 3 mixed forest sites (MF).

Data
Continuous records of half-hourly GPP were calculated based on the net ecosystem exchange (NEE) and ecosystem respiration (Re), with Re determined with nocturnal thermal sensitivity curves [34].The data gaps of NEE associated with equipment failures and unsuitable micrometeorological conditions were filled with the temperature response approach or light response approach [34].A continuous one-day GPP and ET were summed to obtain daily values.The daily values were then averaged to 8-day intervals in accordance with the temporal resolution of the MODIS data.These EC-derived GPP and ET were then used to calibrate and test the model.Meteorological data (i.e., air temperature, precipitation, relative humidity, wind speed, R n and PAR) were the input driving variables for the modified SWH model.In situ measurements of these variables at the flux tower sites were used when deriving the site-specific parameter sets.In addition to the flux data from FLUXNET towers, we also collected data of annual values of GPP and ET at 43 non-Chinaflux sites in China, which was extracted from the published literature (site information is available in Yu et al. [35]).Since this dataset included only the annual values, it was not used for parameterization, but only for testing the model performance.
To estimate regional ET in China based on in situ measurements at 750 nationwide meteorological stations, air temperature, precipitation, relative humidity and wind speed were spatially interpolated at a resolution of 10 km with the method of a bi-variate thin plate with elevation as a covariate using the ANUSPLIN software package [36].R n was calculated as the balance of incident downward and upward shortwave and longwave radiation, in which the MODIS albedo product was used to estimate the net downward shortwave radiation, and the model of Allen et al. [37] was used to estimate the net upward radiation.Incident shortwave radiation and PAR were estimated with the Ångstrøm equation [38].Detailed information of this method is available in Zhu et al. [39].
NDVI is another key driving variable in the SWH model.MODIS NDVI (MOD13Q1) data are satellite products acquired from the website of ORNL DAAC (1 km, http://daac.ornl.gov).These MODIS products contain some cloud-contaminated or missing data.Therefore, before inputting the data into the model, these products were processed using the software package TIMESAT3.0(an asymmetric Gaussian method) to exclude noise and fill gaps [40].

Parameterization
We used in situ measurements of ET and GPP with EC systems at the 63 sites across the globe to optimize b 2 , b 3 , a 1 , g 0 , ε max and d.We first optimized the parameters site by site, yielding 63 site-specific parameter sets.The parameters were optimized by Monte Carlo simulation, which was described in detail in our previous study [14].Briefly, we performed 10,000 Monte Carlo simulations to select ten top-performance parameter sets, and the mean of the ten top-performance parameter sets was regarded as the best-fit parameter set.As expected, the site-specific parameters resulted in desirable model performance at most sites.For regional application of the model, we next calculated the bulk average values of the parameters for each biome based on the site-specific parameters; i.e., all parameter sets of the sites belonging to the same biome were averaged as the parameter set for that biome.Note that the flux data and meteorological data at some sites were gap-filled due to diverse technique problems causing the data missing, which might cause large uncertainties if the gaps were very big, and this in return may cause unreliable parameterization.To minimize this adverse effect we excluded the sites with poor calibration performance, i.e., if the R 2 between EC measurements and model predictions was less than 0.6, the parameter estimated at this site was not used to estimate the biome-average parameters.Finally, the estimated parameter sets at three sites (i.e., ID-Pag, US-Ha1, and DHS in Table S1) were excluded.
At the very beginning, we calibrated the model using the measurements from half of the sites, and the measurements from the remaining sites were used to test the model performance (Figure S1).Comparing with the results with all sites being used for calibration, we yielded very close parameter estimates for each biome, and the simulated ET was nearly identical (data not shown).Assuming that using data of more sites would make the estimated parameters more convincing, we therefore used the parameters with all sites for calibration in this study.All results of the model simulation in this paper were based on the biome-specific parameter sets.

Model Performance and Sensitivity Analysis
With the calculated biome-specific parameters and soil water content being estimated with the bucket scheme, we evaluated the model with two independent datasets of ET and GPP.The first dataset was the in situ measurements at the 63 flux tower sites, and the second was the GPP and ET in ecosystems in China extracted from published literature (see the details in Section 2.3).The driving meteorological variables for this first dataset were the in situ measurements at the tower sites and were the 10-km interpolated products for the second dataset.
The determinant coefficients of the relationship between model and observation (R 2 ) and the root mean square error (RMSE) were used to quantify the difference between the model predictions and observations: where O i and M i are the observed and modeled values, respectively.N is the total number of observations at one site.Sensitivity analysis is a critical step in the identification of the potential parameters and driving variables that cause model uncertainty.Some advanced methodologies for performing sensitivity analysis have been developed in recent years [41].Here, we conducted a conventional one-at-a-time sensitivity analysis.This method involves varying the model parameters individually while fixing the remaining parameters [42].We ran the SWH model by increasing or decreasing each of the model parameters or driving variables by 10% while fixing the others.The model output sensitivity (in %) was calculated as the relative differences in mean daily ET between the model runs (run +10% or run −10% ) and the reference run (run base ): The maximum value between VR +10% and VR −10% was regarded as the sensitivity of the modeled ET to the selected parameter or the driving variable.In the sensitivity analysis, we mainly focused on the six estimated parameters with large variability among the sites (Table 1).In addition, because the soil water content was estimated and the NDVI product may have large uncertainty under cloudy weather conditions, these two driving variables were also included in the sensitivity analysis.* The number of flux sites of savanna, wetland and shrubland is too small to estimate their parameters with acceptable confidence.Considering these biomes types are physically close to grassland and the values of optimized parameters are also close, we combines the four biome types and use the same parameter set.PPFD: photosynthetic photon flux density (µmol•m −2 •s −1 ).

Parameterization and Model Sensitivity
The optimized parameters for each biome type are presented in Table 1.The slope coefficient in the Ball-Berry model, a 1 (Equation ( 2)), was estimated as 7.5 to 10.3.The maximum light-use efficiency, ε max , was estimated as 0.0011 mg•CO 2 •µmol −1 PPFD for forest and grassland and 0.0022 mg•CO 2 •µmol −1 PPFD for cropland, indicating that croplands have the maximum light-use efficiency, possibly twice as high as that of forests and grasslands.However, a large site-to-site variability of the parameters was detected within each biome type.For example, the magnitude of variability (standard deviations) of g 0 is comparable to the mean values.In comparison, the within-biome variability was least for ε max , which was relatively conserved across sites.
Sensitivity analysis results indicated that the modeled ET was most sensitive to NDVI; 10% changes in NDVI resulted in approximately 4% changes in ET (Table 2).NDVI was followed by a 1 (or ε max , ca. 3%), g 0 (1.0% to 1.8%), b 2 (0.5% to 2.5%) and b 3 (0.5% to 1.3%).The sensitivity to soil water content and soil depth was minor.Ten percent changes in each of these two parameters introduced changes of less than 1% in ET (except a 1.7% change for SWC in grassland).Comparing the parameter sensitivities among the vegetation types revealed that the model was significantly more sensitive to the parameters when estimating r ss in grassland than in cropland and forest (Table 2).For example, the sensitivity of b 2 in grassland (2.5%) was three-fold greater than that in forest (0.8%) and five-fold greater than that in cropland (0.5%).In addition, the model sensitivity of SWC in grassland (1.7%) was obviously much larger than in forest (0.5%) and cropland (0.4%).These results imply that the accurate estimation of soil surface resistance is more important for grassland than for the other two biome types.

Model Performance
The SWH model performed well for both ET and GPP simulations.Collectively, the SWH model explained 65% of the variation in eight-day ET and 68% in GPP at the 63 sites (p < 0.01, Figure 2).By comparison, the MODIS ET algorithm explained 53% of the variations in ET (p < 0.01).In general, both the SWH model and the MODIS ET algorithm underestimated ET at the sites or in the season with high ET.
Compared with the MODIS algorithm, the SWH model yielded a higher R 2 at 58 of the 63 sites and lower RMSE at 36 of the 63 sites.The bulk average R 2 of the SWH model for all 63 sites was 0.13 higher (0.77 for SWH and 0.64 for MODIS), and the average RMSE was 0.11 mm•day −1 lower (0.72 for SWH and 0.83 for MODIS).The performance of the two models was comparable for DBF, savanna and shrubland, whereas SWH exhibited superior performance for EBF, wetland and MF (Figure 3).Remote Sens. 2017, 9, 54 8 of 20      When the interpolated climate data were used as the input variables, the model predictions were also consistent with the observations (i.e., the annual values of GPP and ET at 43 non-Chinaflux sites in China published in the literature).The model explained 77% and 84% of the spatial variability of GPP and ET, respectively (Figure 4).Note that similar to the validation for flux sites at the eight-day scale, certain underestimations of both GPP and ET were observed for the sites (or years) with high GPP and ET values.In addition, the regionalized model (i.e., driving with interpolated climate data) also illustrated good performance at the eight-day timescale.The model explained 62% and 77% of the variations in ET and GPP, respectively (Figure S2).When the interpolated climate data were used as the input variables, the model predictions were also consistent with the observations (i.e., the annual values of GPP and ET at 43 non-Chinaflux sites in China published in the literature).The model explained 77% and 84% of the spatial variability of GPP and ET, respectively (Figure 4).Note that similar to the validation for flux sites at the eight-day scale, certain underestimations of both GPP and ET were observed for the sites (or years) with high GPP and ET values.In addition, the regionalized model (i.e., driving with interpolated climate data) also illustrated good performance at the eight-day timescale.The model explained 62% and 77% of the variations in ET and GPP, respectively (Figure S2).  1) are not included in the illustration.

Regional GPP, ET, WUE and E/ET in China
According to the model simulations performed at the regional scale, GPP and ET exhibit similar spatial patterns in China (Figure 5).Ecosystems with high GPP and ET were mainly distributed in southern and eastern China, where forest and cropland are the dominant ecosystem types.The sites used for model parameterization (Table 1) are not included in the illustration.

Regional GPP, ET, WUE and E/ET in China
According to the model simulations performed at the regional scale, GPP and ET exhibit similar spatial patterns in China (Figure 5).Ecosystems with high GPP and ET were mainly distributed in southern and eastern China, where forest and cropland are the dominant ecosystem types.Ecosystems with low GPP and ET were mainly distributed in northern and western China (including the Tibetan Plateau), where grassland is the dominant ecosystem type.The average ET of the ecosystems in China in 2000 to 2010 was 408.7 mm•year −1 .Similarly, the mean value of GPP was 809 g•C•m −2 •year −1 , WUE was 1.68 g•C•mm −1 •H 2 O and E/ET was 55%.Spatially, WUE was high in the southern and eastern areas, where forest and cropland are the major ecosystem types, but low in the northern and western areas, where grassland is the dominant ecosystem type.For E/ET, however, the spatial pattern was the opposite.
the Tibetan Plateau), where grassland is the dominant ecosystem type.The average ET of the ecosystems in China in 2000 to 2010 was 408.7 mm•year −1 .Similarly, the mean value of GPP was 809 g•C•m −2 •year −1 , WUE was 1.68 g•C•mm −1 •H2O and E/ET was 55%.Spatially, WUE was high in the southern and eastern areas, where forest and cropland are the major ecosystem types, but low in the northern and western areas, where grassland is the dominant ecosystem type.For E/ET, however, the spatial pattern was the opposite.

Parameters of the SWH Model
The values of the optimized parameters were consistent with previously-published results.For example, the slope coefficient in the Ball-Berry model (Equation ( 2)) was estimated as 7.5 to 10.3, consistent with the values in Community Land Model(CLM)4.0[43].The estimated maximum lightuse efficiency values (0.0011 mg•CO2•µmol −1 PPFD for forest and grassland and 0.0022 mg•CO2•µmol −1 PPFD for cropland) are located in the middle range of εmax variations in the meta-analyses of Kergoat et al. [44] and Garbulsky et al. [45].Note that we derived εmax using the NDVI data to estimate FPAR (Equation (3)).A higher εmax would have been obtained if the Enhanced Vegetation Index (EVI) had been used to derive this parameter due to the lower magnitude of EVI compared to NDVI.Our results indicated that the maximum light-use efficiency is two-fold higher in cropland than in forest and grassland, most likely due to the intensive management of cropland through techniques such as irrigation and nutrient manipulation.

Parameters of the SWH Model
The values of the optimized parameters were consistent with previously-published results.For example, the slope coefficient in the Ball-Berry model (Equation ( 2)) was estimated as 7.5 to 10.3, consistent with the values in Community Land Model(CLM)4.0[43].The estimated maximum light-use efficiency values (0.0011 mg•CO 2 •µmol −1 PPFD for forest and grassland and 0.0022 mg•CO 2 •µmol −1 PPFD for cropland) are located in the middle range of ε max variations in the meta-analyses of Kergoat et al. [44] and Garbulsky et al. [45].Note that we derived ε max using the NDVI data to estimate FPAR (Equation (3)).A higher ε max would have been obtained if the Enhanced Vegetation Index (EVI) had been used to derive this parameter due to the lower magnitude of EVI compared to NDVI.Our results indicated that the maximum light-use efficiency is two-fold higher in cropland than in forest and grassland, most likely due to the intensive management of cropland through techniques such as irrigation and nutrient manipulation.
Our results indicated large site-to-site variability for most of the estimated parameters, which implies that the use of the look-up table approach for each biome type will introduce uncertainties into the model.For better region-scale performance of the model, we originally expected to identify some empirical relationships between the parameters and climate or biotic variables for application to the estimation of parameters at other sites.Unfortunately, we failed to identify any regression correlation between the parameters and climatic-or vegetation-related variables.The factors affecting the spatial variations in the parameters may be too complex to use a simple regression method to depict the spatial pattern.Alternatively, the number of sites in this study may not be sufficient to clarify the mechanism that controls the spatial patterns of the parameters.Data from more sites may clarify the spatial patterns of the parameters.In contrast to most parameters, the estimated maximum light-use efficiency, ε max , was relatively constant across sites (Table 2).This finding is consistent with a recent study by Yuan et al. [46], who reported that the use of a constant ε max introduces only a minor bias in the estimation of GPP.
The data of observations may introduce some uncertainties for parameterization due to the effect of gap-filling at some flux sites.However, this uncertainty is minor and has quite little impact on the estimate of biome-specific parameters.First, the quality of the FLUXNET dataset had been comprehensively assessed and controlled before release.The very bad data were assigned the value of −9999, which was excluded for parameterization in our study.Second, we have excluded the sites with poor calibration performance, which could be possibly due to the uncertainty from observations (see the Methods).In addition, assuming that the gap-filled observation data exert large uncertainties on model simulation, one can expect that the model performance would be very poor, and the parameterization results with data from different groups of sites would be much different and the simulated ET inconsistent.However, our results indicate that the model illustrated reliable performance, implying that the estimated parameter sets were reliable.Furthermore, using the data of 32 sites and 64 sites to estimate the parameters, respectively, we found that the two methods yielded very close values of biome-specific parameters, and the simulated ETs with the two parameter sets were nearly identical (data not shown).

Reliability and Merits of the SWH Model
Compared with the EC observations, the SWH model yielded an average R 2 of ca.0.8 and a RMSE of ca.0.7 mm•day −1 at all 63 sites (Table 3).This performance is comparable to the performance achieved using models driven by climatic and remote sensing data.Mostly, models yielded an R 2 close to or less than that of SWH, e.g., 0.71 and 0.68 for the Remote Sensing Penman Monteith (RS-PM) model and EC-LUE model [30], 0.64 for the modified PM model [47] and 0.75 for the Breathing Earth System Simulator (BESS) model [48].In addition, Vinukollu et al. [49] evaluated the performances of three process-based models, including a single source energy budget model, a Penman-Monteith-based approach and a Priestley-Taylor-based approach at the global scale.Their intercomparison of model estimates with 12 EC towers yielded a mean correlation of 0.57.Chen et al. [6] reported mean R 2 and RMSE values of 0.5 to 0.8 and 0.45 to 0.75 mm•day −1 , respectively, for eight ET models at 23 EC sites.
Compared with previous model reports, our estimated mean annual ET in China (408.7 mm) is within the range of the global average, i.e., 310 to 547 mm [5].This value is slightly higher than the estimate of Yao et al. [50], who reported a mean annual ET in China of 364.9 mm.Recently, Chen et al. [6] estimated ET in China using eight ET models; the annual ET ranged from 535 to 852 mm among the models.The variation in these results may be the result of differences among model structures [6].In addition, the different sources and spatial resolutions of the input climate data used in the different models may also underlie this variation [51].
Our estimated E/ET (0.55) is consistent with the estimated global average based on ten process models, which yielded values ranging from 0.36 to 0.75, with an average of 0.58 [5].However, a recent meta-analysis by Schlesinger and Jasechko [52] indicated that global E/ET is 39% (±15% SD), suggesting that the average E/ET in ecosystems in China, when estimated using the SWH model, may be higher than the global average.This higher estimate may be partly due to the high proportion of grassland in China (ca.40% of the country).In addition, the model's underestimation of GPP may also contribute to the higher E/ET.
The SWH model simulated that the ecosystems with high ET were mainly distributed in southern and eastern China, whereas low-ET ecosystems were mainly distributed in northern and western China.Similarly, Yao et al. [50] and Chen et al. [6] also reported that ET was highest in the humid tropics and sub-tropics, intermediate in temperate regions and lowest in both cold and arid regions.Water-use efficiency (WUE) exhibited a spatial pattern similar to the patterns observed for GPP and ET, with high values in the forest ecosystems in southern and eastern China and low values in the grasslands in northern China.This result is consistent with the results of inter-site comparisons based on EC measurements [14,53].This spatial pattern may be primarily due to the lower E/ET in ecosystems with dense canopies, which may allow these ecosystems to use water more efficiently [54].
As a simplified water-carbon coupling model, the SWH model represents a compromise between biogeochemical models that involve many detailed processes and conventional remote sensing ET models that only involve several empirical functions.By using the scheme of the S-W model, the SWH model estimates the water vapor fluxes in ecosystems with a mechanistic approach.However, by introducing the stomatal conductance model and the light-use efficiency-based GPP model, the SWH model largely reduces the processes and parameters required to simulate the photosynthesis process (Table 4).In addition to ET partitioning, the SWH model also estimates GPP, another key variable in ecological processes, which helps to address the coupling relationship between carbon and hydrological processes (e.g., ecosystem WUE) at multiple spatiotemporal scales [54].The third merit of the SWH model is that relatively few input variables are required.One remote sensing product (i.e., NDVI) and six meteorological variables (i.e., air temperature, precipitation, relative humidity, wind speed, R n and PAR) are the input driving variables for the model (Figure 1).Most of the variables are accessible from public data sources, facilitating regional and global applications.

Uncertainties of the SWH Model
Our results illustrated different performances across sites and PFTs.Further analysis indicated that the model performance for ET prediction was significantly correlated with the performance for GPP prediction (p < 0.01, Figure S3), indicating that accurate estimate of GPP is critical for ET simulation in the SWH model.In this study, GPP was estimated using a light-use efficiency model.The remotely-sensed NDVI product was a key driving variable in GPP estimation, which may introduce large uncertainty in some regions (discussed below).Uncertainties may also be incurred by the inaccurate estimate of light-use efficiency.Many efforts have been made to solve this issue, such as improving drought constraint [51], accounting for diffusive radiation [57,58] and updating the empirical values in the equations of environmental constraints (e.g., T min , T opt , T max , VPD max ) in the GPP model [31].Therefore, improving GPP estimation would be anticipated with the outcome of these efforts, which thus will be a benefit for improving ET modeling with the SWH model.
Our results illustrated that the model performances were poorer in EBF and grasslands.Monitoring of tropical forests using shortwave radiation from satellite-based observations is difficult due to atmospheric effects, particularly severe cloud contamination in tropical regions with EBF ecosystems [6,51].Therefore, the lower-quality NDVI product likely underlies the poorer model performance in EBF (Figure 3).For grasslands, the most sensitive parameters are likely responsible for the model uncertainty.As our sensitivity analysis indicated, the model was much more sensitive to the parameters used to estimate the soil surface resistance in grassland than in cropland and forest (Table 2).In addition, inaccurate estimates of GPP may be another important source of the poorer estimates of ET in grassland.For example, we recently compared four widely-used light-use efficiency-based GPP models.The results indicated that most models yielded poorer estimates of GPP in grassland than other ecosystem types [31].
The different model performances among biome types may also partly be determined by the scheme used to estimate ET components.In the community of the two-source ET models, three different vegetation parameterization schemes are adapted, i.e., layer approach, patch approach and hybrid approach.The SW model uses the layer approach, and it does not distinguish soil evaporation from inter-canopy soil surface and under-canopy soil surface.One limitation of the layer approach is that it suits uniformly-vegetated surfaces better [59].This may be an additional reason why the SWH model performed better in biomes with a more uniform vegetation distribution.
Another uncertainty associated with the SWH model is its lack of consideration of intercepted canopy water evaporation, which may introduce biased estimations of ET in forest ecosystems with high and frequent rainfall events and, consequently, potentially reduce the model performance at EBF sites.However, the MODIS algorithm, which estimates wet canopy evaporation independently, also clearly exhibited poorer predictive ability in EBF (Figure 3).This result suggests that canopy interception is not an important source of the inconsistencies between the SWH model and observations at EBF sites.
Soil water content is a necessary driving variable in the SWH model and was estimated using a one-layer bucket model.We expected that inaccurate estimates of soil water content would introduce large uncertainty into the modeled ET.However, the sensitivity analysis indicated that the modeled ET does not depend on soil water as strongly as anticipated.A 10% change in soil water content introduced a less than 1% change in ET (Table 2).This result is reasonable considering the time step of the model.It is acknowledged that ET would depend highly on soil water content at diurnal and daily time scales.At daily and lower time scales, the effect of soil water content is largely through impacts on plant stomatal conductance.In our study, ET was estimated at an eight-day time step, in which the effect of soil water content would largely ne via impacts on Leaf Area Index (which can be quantified with NDVI), instead of stomata.This might be the reason why NDVI is the most sensitive variable in the model.In addition, soil moisture was only used to estimate the soil surface resistance of evaporation, r ss , in the SWH model (Equation (1)), which makes the ET simulation depend on soil water content very little.
The SWH model simulates ET, as well as the major components of soil water evaporation and plant transpiration.However, we failed to quantify the model performance for estimating each component in this study.We recently analyzed eddy flux tower measurements from three Chinaflux forest sites at which flux towers were installed both above and below the canopy, allowing us to partition ET with direct measurements.The results indicated that the model outputs were generally consistent with the observations [53].However, due to the scarcity of such observations, and the possible overestimation of E/ET mentioned above, we acknowledge that further work is needed to fully evaluate the partitioning performance of the model.

Conclusions
The SWH model prediction of ET and GPP agreed well with the observations recorded at global flux tower sites.In general, the model performance was better for woody ecosystems than for herbaceous ecosystems.Using interpolated climate data for China as the input climate forcing, the model simulations were also consistent with the measurements.Based on the simulation, high ET, GPP, WUE and E/ET values were observed in southern and eastern China, where forest and cropland are the dominant ecosystem types.The opposite pattern was observed in northern and western China, where grassland is the dominant ecosystem type.Notably, the model performance was poorer where ε max is the apparent quantum yield or maximum light-use efficiency and f (T), f (W) and f (VPD) are the downward regulation scalars for the effects of temperature and VPD on light-use efficiency of vegetation, respectively.T min , T max and T opt are the minimum, maximum and optimum air temperature ( • C) for photosynthetic activity, respectively.If the air temperature falls below T min or increases above T max , f (T) is set to zero.In this study, T min , T opt and T max were set to 0, 20 and 40 • C, respectively.If the estimated f (VPD) was less than zero, f (VPD) was set to zero [27] VPD max was set to 3.5 kPa.

Figure 1 .
Figure 1.Flowchart of evapotranspiration (ET) estimation with the SWH model.The meanings of the acronyms in this figure are available in Appendix B.

Figure 1 .
Figure 1.Flowchart of evapotranspiration (ET) estimation with the SWH model.The meanings of the acronyms in this figure are available in Abbreviation.

Figure 2 .
Figure 2. Performances of the SWH model in estimating: (a) evapotranspiration (ET); and (b) gross primary productivity (GPP); and the performance of the MODIS ET algorithm in estimating ET (c).All of the sites used in this study use an eight-day timescale (n = 9012).ET_EC and GPP_EC are the ET and GPP derived from the measurements using eddy covariance (EC) systems.

Figure 3 .
Figure 3. Mean difference between the model predictions (including the SWH model and the MODIS ET algorithm) and observations in terms of (a) the determinant coefficient of the model-observation relationship (R 2 ); and (b) the root mean square error (RMSE) for each plant functional type.The numbers in parentheses indicate the number of sites used in this study.The error bars indicate standard deviations.

Figure 2 .of 20 Figure 2 .
Figure 2. Performances of the SWH model in estimating: (a) evapotranspiration (ET); and (b) gross primary productivity (GPP); and the performance of the MODIS ET algorithm in estimating ET (c).All of the sites used in this study use an eight-day timescale (n = 9012).ET_EC and GPP_EC are the ET and GPP derived from the measurements using eddy covariance (EC) systems.

Figure 3 .
Figure 3. Mean difference between the model predictions (including the SWH model and the MODIS ET algorithm) and observations in terms of (a) the determinant coefficient of the model-observation relationship (R 2 ); and (b) the root mean square error (RMSE) for each plant functional type.The numbers in parentheses indicate the number of sites used in this study.The error bars indicate standard deviations.

Figure 3 .
Figure 3. Mean difference between the model predictions (including the SWH model and the MODIS ET algorithm) and observations in terms of (a) the determinant coefficient of the model-observation relationship (R 2 ); and (b) the root mean square error (RMSE) for each plant functional type.The numbers in parentheses indicate the number of sites used in this study.The error bars indicate standard deviations.

Figure 4 .
Figure 4. Comparisons of the yearly values of modeled GPP (a); and ET (b) with observations from EC systems for ecosystems in China (n = 49).The sites used for model parameterization (Table1) are not included in the illustration.

Figure 4 .
Figure 4. Comparisons of the yearly values of modeled GPP (a); and ET (b) with observations from EC systems for ecosystems in China (n = 49).The sites used for model parameterization (Table1) are not included in the illustration.

Figure 5 .
Figure 5. Spatial patterns of ET, GPP, water-use efficiency (WUE: GPP/ET) and the ratio of soil water evaporation to ET (E/ET) in China from 2001 to 2010.

Figure 5 .
Figure 5. Spatial patterns of ET, GPP, water-use efficiency (WUE: GPP/ET) and the ratio of soil water evaporation to ET (E/ET) in China from 2001 to 2010.

Table 1 .
Look-up table of model parameters for each biome type.The values in parentheses are standard deviations.

Table 3 .
Performances of the SWH and MODIS algorithms in ET and GPP prediction.R 2 _ET, R 2 _GPP and R 2 _MODIS are the R 2 of the relationships between ET modeled by SWH, GPP by SWH, MODIS ET and the observations, respectively.RMSE_ET (mm•day −1 ), RMSE_GPP (g•C•m −2 •day −1 ) and RMSE_MODIS (mm•day −1 ) are the root mean square errors between ET modeled by SWH, GPP by SWH, MODIS ET and the observations, respectively.PFT, plant functional type.

Table 4 .
List of ET models driven by remote sensing and climate data.Note that the models estimating ET indirectly with the energy balance approach are not included.PML: Penman-Monteith-Leuning, PT-Fi: Priestley-Taylor-Fisher. RS-PM: Remote Sensing Penman-Monteith.BESSS: Breathing Earth System Simulator.