A Comparative Study of GRACE with Continental Evapotranspiration Estimates in Australian Semi-Arid and Arid Basins: Sensitivity to Climate Variability and Extremes

This study examines the dynamics and robustness of large-scale evapotranspiration products in water-limited environments. Four types of ET products are tested against rainfall in two large semi-arid to arid Australian basins from 2003 to 2010: two energy balance ET methods which are forced by optical satellite retrievals from MODIS; a newly developed land surface model (AWRA); and one approach based on observations from the Gravity Recovery and Climate Experiment (GRACE) and rainfall data. The two basins are quasi (Murray-Darling Basin: 1.06 million km2) and completely (Lake Eyre Basin: 1.14 million km2) endorheic. During the study period, two extreme climatic events—the Millennium drought and the strongest La Niña event—were recorded in the basins and are used in our assessment. The two remotely-sensed ET products constrained by the energy balance tended to overestimate ET flux over water-stressed regions. They had low sensitivity to climatic extremes and poor capability to close the water balance. However, these two remotely-sensed and energy balance products demonstrated their superiority in capturing spatial features including over small-scale and complicated landscapes. AWRA and GRACE formulated in the water balance framework were more sensitive to rainfall variability and yielded more realistic ET estimates during climate extremes. GRACE demonstrated its ability to account for seasonal and inter-annual change in water storage for ET evaluation.


Introduction
Evapotranspiration (ET) governs water and energy exchange between the atmosphere and the Earth's surface [1]. Globally, more than half of the solar energy absorbed by the land surface is used for evaporation and transpiration [2]. Annually, approximately 60% of the water precipitated over  [37]): A, equatorial; B, arid; C, warm temperate; W, desert; S, steppe; f, fully humid; s, summer dry; m, monsoonal; w, winter dry; h, hot arid; k, cold arid; a, hot summer; b, warm summer. (b,c) Digital elevation model (DEM) and land use maps were accessed from [38,39], respectively. The Lake Eyre Basin (1.14 million km 2 ) encompasses 81.5% hot dessert (arid climate) and ~15% hot steppe (semi-arid climate). Rainfall within the basin varies widely and is usually unable to meet the atmospheric demand. The annual rainfall over the arid region of the LEB barely reaches 200 mm/year while the north-east edge receives rainfall of up to 700 mm/year under the influence of infrequent tropical storms [40]. The annual potential evaporation (PET) averaged from 1961 to 1990 is of 1453 mm/year (sourced from Australian Bureau of Meteorology data). Floods in the basin are ephemeral and extremely variable [41]. Flows and floodwaters that form during heavy rainfall events carry only 1% of total rainfall; a large fraction (~99%) of rainfall is lost through evaporation and transpiration [35].The Murray-Darling Basin (1.06 million km 2 ) contains a transition from subtropical to dry arid climate. Rainfall distributions within the basin vary greatly, decreasing from the south-eastern and eastern boundaries (between 600 and 800 mm) towards its western and north-western boundaries (between 100 and 300 mm). The average annual PET was 1236 mm/year . The basin consists of three large river systems and ~30,000 wetlands [42]. As the country's food bowl, ~80% of its area is used for agriculture [43]. Nearly half the surface water is redistributed by irrigation [35]. Many areas in the basin have little or no regular runoff [36,44]; the water levels were reduced to historically low levels during the Millennium Drought from mid-1990s to 2009 [45,46].

Rainfall, Potential Evaporation, and Discharge Data
Rainfall data, provided by the Australian Bureau of Meteorology (BoM, Melbourne, Australia, data access: [47]), were used for two purposes: (1) to be used with GRACE's measurements to obtain  [37]): A, equatorial; B, arid; C, warm temperate; W, desert; S, steppe; f, fully humid; s, summer dry; m, monsoonal; w, winter dry; h, hot arid; k, cold arid; a, hot summer; b, warm summer. (b,c) Digital elevation model (DEM) and land use maps were accessed from [38,39], respectively. The Lake Eyre Basin (1.14 million km 2 ) encompasses 81.5% hot dessert (arid climate) and~15% hot steppe (semi-arid climate). Rainfall within the basin varies widely and is usually unable to meet the atmospheric demand. The annual rainfall over the arid region of the LEB barely reaches 200 mm/year while the north-east edge receives rainfall of up to 700 mm/year under the influence of infrequent tropical storms [40]. The annual potential evaporation (PET) averaged from 1961 to 1990 is of 1453 mm/year (sourced from Australian Bureau of Meteorology data). Floods in the basin are ephemeral and extremely variable [41]. Flows and floodwaters that form during heavy rainfall events carry only 1% of total rainfall; a large fraction (~99%) of rainfall is lost through evaporation and transpiration [35].The Murray-Darling Basin (1.06 million km 2 ) contains a transition from subtropical to dry arid climate. Rainfall distributions within the basin vary greatly, decreasing from the south-eastern and eastern boundaries (between 600 and 800 mm) towards its western and north-western boundaries (between 100 and 300 mm). The average annual PET was 1236 mm/year . The basin consists of three large river systems and~30,000 wetlands [42]. As the country's food bowl,~80% of its area is used for agriculture [43]. Nearly half the surface water is redistributed by irrigation [35]. Many areas in the basin have little or no regular runoff [36,44]; the water levels were reduced to historically low levels during the Millennium Drought from mid-1990s to 2009 [45,46].

Rainfall, Potential Evaporation, and Discharge Data
Rainfall data, provided by the Australian Bureau of Meteorology (BoM, Melbourne, Australia, data access: [47]), were used for two purposes: (1) to be used with GRACE's measurements to obtain basin-scale ET estimates through the water balance equation; (2) to evaluate ET datasets and identify their potential deficiencies. The rainfall data is a gridded daily product interpolated from weather station observations [48].
In addition, potential evaporation (hereafter, denoted as E p ) data was obtained from Australian Water Availability Project (AWAP, Canberra, Australia, data access: [49]). Within the Priestley-Taylor framework, E p was forced by the gridded meteorological data obtained from BoM [50]. A comparison of rainfall with potential evaporation gives a first indication of areas that are mainly controlled by the availability of water versus energy. Both the rainfall and E p datasets have a spatial resolution of 0.05 • × 0.05 • , and are available from January 2003 to December 2010.
We took the Lock 1 discharge measurements (see the location in Figure 1c) as the river outflow from the MDB. The data was originally archived as daily records (accessed from: [51]), and were further aggregated to monthly equivalent water height (mm/month) divided by total basin area. The LEB is a completely closed basin and no discharge data were required for this basin.

Model-Based ET Estimates
Three continental modelled ET products were used in this study: PT-CMRS, PM-Mu (a.k.a. MOD16), and AWRA-L. They vary in approach, data inputs, and ground calibration (see Table 1). To make these datasets consistent in time and space, ET results sourced from the three models were converted into monthly values from 2003 to 2010 at 0.05 • (~5 km) spatial resolution. A brief summary of the forcing datasets, mechanism, and the major features for each ET model is given in Table 1.
• PT-CMRS model The model. forced by MODIS retrievals and gridded meteorological datasets, was modified based on the Priestley-Taylor (PT) formula ( [52]; see Equation (1)) to dynamically represent the actual ET variations over the Australian continent (hereafter, denoted as "PT-CMRS").
where R n is the surface net radiation(MJ m −2 d −1 ); G is the soil heat flux (MJ m −2 d −1 ); ∆ is the slope of the curve relating saturation water vapor pressure to temperature (kPa K −1 ); γ is the psychometric constant (kPa K −1 ). The prominent advantage associated with the PT method is that it does have limited input data requirement; wind speed is not compulsory in the model. The developers replaced the original empirical constant α = 1.26 with a flexible scaling factor formulated by the Enhanced Vegetation Index (EVI) and the Global Vegetation Moisture Index (GVMI) [53]. This scheme allows GVMI to separate surface water bodies against bare soil when EVI is low and to detect vegetation water content when EVI is high. Also, the PT-CMRS model accounts for precipitation interception. This model was calibrated at seven flux sites in Australia and validated in 227 catchments. Readers can refer the details about the PT-CMRS ET model from [45]. Data can be downloaded from: [54]. • PM-Mu model ( [10,13]) Based on Cleugh's method [10], the MOD 16 ET data set developed by Mu ([13]; denoted here as 'PM-Mu') -integrated MODIS-retrieved leaf area index (LAI) into the Penman-Monteith (PM) equation ( [55]; see Equation (2)), and improved the estimation of surface resistance and soil evaporation [13].
where ρ a is the mean air density at constant pressure (kg m −3 ); c p is the specific heat of air at constant pressure (MJ kg −1 K −1 ); (e s − e a ) is the water vapor pressure deficit between the saturated air pressure and the actual air pressure (kPa); r s and r a represent the surface/aerodynamic resistance (m s −1 ); λ is latent heat of evaporation (MJ m −2 d −1 ). Further modifications were described in [19]. The dataset is available at global scale. Results have been tested at global and flux-site scales. Data can be downloaded from: [56]. The PM model is process-based and constrained by energy balance. It requires considerably more input data than the PT model; some of this data (especially the wind speed) are barely available over large basins or regions. However, meteorological data have been obtained from atmospheric model reanalysis or gridded meteorological datasets.
• ET estimates from AWRA-L land surface model The AWRA (Australian Water Resources Assessment) ET output (denoted as "AWRA") was produced by a multi-model system simulating hydrological processes and dynamics in landscape, river and groundwater systems, and water use all over Australia [57]. This system is jointly developed by the BoM and the Commonwealth Scientific and Industrial Research Organization (CSIRO). Either the PM or the PT functions were used, depending on the availability of wind speed data [58]. Within each cell, ET is summarized as: where E t is vegetation transpiration; E i , E r , E g , and E s are evaporation from rainfall interception, open water bodies, groundwater, and soil profiles, respectively. AWRA ET estimates balance the requirement between water/energy conservation and data unavailability. However, the model neglects lateral water flows, which leads to an underestimation of ET over areas receiving inflows, such as wetlands, floodplains, and irrigated farmland. Gauged runoff and eddy covariance flux tower observations were used to visually fit some components and parameters of the model.

ET Estimates from GRACE Rainfall and Discharge Observations
A number of studies have derived ET estimates from GRACE observations in large basins, e.g., [11,33,34].
Here we computed ET estimates using regional GRACE solutions which are characterized by reduced north-to-south striping (using constrained regularization) and no contamination from other parts of the world [59,60]. Several studies over different continents-South America [61], Australia [62], and Africa [63]-demonstrated that this regional approach offers a reduction of both north-south striping due to the distribution of GRACE satellite tracks and temporal aliasing of correcting models that are present in the global GRACE solutions. GRACE regional solutions were available at a spatial resolution of 2 • × 2 • from 4 July 2003 to 3 December 2010 at intervals of 10 days [62]. Change in the basin water storage (∆S) is calculated as the difference between two successive GRACE terrestrial water storage anomalies (TWSA) against the average (∆TWS): Using the water balance equation at basin-scale, ET was obtained as the difference between the total amount of rainfall over the basin (P), the river discharge at the outlet (Q) and the change in ∆S over a specific time period ∆t: To be consistent with the above modelled ET datasets, calculations were performed over 30-day time period. Missing data from 20 January 2004 to 29 January 2004 were linearly interpolated. Combining GRACE ∆TWS observations and BoM rainfall datasets, monthly basin-average ET estimates over the MDB were computed using Equation (5) with all variables set. However, due to the fact that Q is unavailable at gridded-cell scale, GRACE ET maps were approximated as (P − ∆S) with Q assumed to be negligible. This operation was also applied to the LEB.

Evaluation of ET Estimates Using the Budyko Diagram Scheme
Water balance at basin-scale is governed by Equation (5). In the Budyko framework, the nature of the annual water balance is determined by the ratio E/P (evaporation efficiency) as a function of E p /P (drought index or climatic aridity) accounting for the partition of rainfall between evaporation and runoff [64]: In extremely dry cases, if a basin is provided with sufficient evaporative energy but limited precipitation, then E will approximate to P; conversely, E will be mainly determined by E p in a wet basin. These form two asymptotes to constrain the Budyko curve in a boundary. A Budyko diagram is usually used at long-term average scale. Since GRACE could provide the annual basin water storage change information, we assumed that the Budyko diagram could be valid at inter-annual time step. Annual E and P were presented as Ea and Pa, respectively.
Noticing that the indices of Ea/Pa and E p /Pa could be sensitive to the different sources of Pa and E p , we assumed that E p /Pa would not vary widely as E p in MDB and LEB is much larger than Pa; Pa is based on measurements from rainfall stations. Figure 2 shows the spatial distribution across the two basins of average annual ET for each dataset and average annual rainfall over the 2003-2010 period. The estimated annual ET maps are spatially consistent with corresponding rainfall maps. The highest ET rates were found in south-eastern MDB (>800 mm/year) where the climate is more humid and intense canopy transpiration occurs. The northern parts of MDB and LEB also present a high ET rate of 400-600 mm/year, which is attributed to the strong ET fluxes caused by tropical rainstorms during the summer. Controlled by an arid climate, 100-200 mm/year of rainfall, the central part of the LEB and the western MDB show the lowest ET rates of all four datasets. However, we identified some differences between ET estimates. Over the extremely dry regions of the central LEB (rainfall: 100-200 mm/year), ET estimates from AWRA and GRACE were more reasonable than the values provided by two energy balance constrained methods (PM-Mu and PT-CMRS) which predicted 200-300 mm/year of water lost via ET processes over the central LEB. Over humid regions, PM-Mu estimated an annual ET flux lower than the other three methods. Figure 3 zooms in three smaller regions (boxes A, B, and C shown in Figure 2) that are prone to inundation (A and B) or irrigation (C). The red patches in contrast to the dry background in PT-CMRS represent the high-ET patterns over irrigated areas as well as flood plains composed of open water bodies and dense forests, demonstrating the dataset's ability to capture relatively small ET features and heterogeneities over those landscapes. This ability is mostly attributed to the integration of MODIS-based EVI and GVMI indices in PT-CMRS [53]. Although it was also driven by MODIS optical products, PM-Mu did not clearly present those small-scale ET features in space. AWRA also considered ET flux from canopy cover and open waters as well as saturated soil. Nevertheless, its coarse gridding system and fractional index parameterization for different land covers make the model homogenize those landscape features.  However, we identified some differences between ET estimates. Over the extremely dry regions of the central LEB (rainfall: 100-200 mm/year), ET estimates from AWRA and GRACE were more reasonable than the values provided by two energy balance constrained methods (PM-Mu and PT-CMRS) which predicted 200-300 mm/year of water lost via ET processes over the central LEB. Over humid regions, PM-Mu estimated an annual ET flux lower than the other three methods. Figure 3 zooms in three smaller regions (boxes A, B, and C shown in Figure 2) that are prone to inundation (A and B) or irrigation (C). The red patches in contrast to the dry background in PT-CMRS represent the high-ET patterns over irrigated areas as well as flood plains composed of open water bodies and dense forests, demonstrating the dataset's ability to capture relatively small ET features and heterogeneities over those landscapes. This ability is mostly attributed to the integration of MODIS-based EVI and GVMI indices in PT-CMRS [53]. Although it was also driven by MODIS optical products, PM-Mu did not clearly present those small-scale ET features in space. AWRA also considered ET flux from canopy cover and open waters as well as saturated soil. Nevertheless, its coarse gridding system and fractional index parameterization for different land covers make the model homogenize those landscape features. However, we identified some differences between ET estimates. Over the extremely dry regions of the central LEB (rainfall: 100-200 mm/year), ET estimates from AWRA and GRACE were more reasonable than the values provided by two energy balance constrained methods (PM-Mu and PT-CMRS) which predicted 200-300 mm/year of water lost via ET processes over the central LEB. Over humid regions, PM-Mu estimated an annual ET flux lower than the other three methods. Figure 3 zooms in three smaller regions (boxes A, B, and C shown in Figure 2) that are prone to inundation (A and B) or irrigation (C). The red patches in contrast to the dry background in PT-CMRS represent the high-ET patterns over irrigated areas as well as flood plains composed of open water bodies and dense forests, demonstrating the dataset's ability to capture relatively small ET features and heterogeneities over those landscapes. This ability is mostly attributed to the integration of MODIS-based EVI and GVMI indices in PT-CMRS [53]. Although it was also driven by MODIS optical products, PM-Mu did not clearly present those small-scale ET features in space. AWRA also considered ET flux from canopy cover and open waters as well as saturated soil. Nevertheless, its coarse gridding system and fractional index parameterization for different land covers make the model homogenize those landscape features.   Figure 2) that are prone to inundation or irrigation. Table 2 shows the mean seasonal ET estimates and corresponding rainfall averaged over the period 2003-2010 and the fractions of ET to rainfall. The four methods used to estimate ET in this study were able to capture the seasonal ET patterns in each basin. LEB is much drier (has little rainfall against high potential evaporation requirement) than the MDB. Therefore, it is likely that the ET fluxes in the LEB is lower than in the MDB through the four seasons. However, in a comparison of the magnitude of seasonal ET to rainfall, the four methods behaved inconsistently. Again, seasonal PM-Mu estimates were constantly lower than other ET datasets in MDB, with the lowest ratio of 0.4 found in winter. The ratios of ET/P are frequently above 1.0 in PT-CMRS and PM-Mu during the dry seasons in LEB; the largest ratios (1.5 for PT-CMRS and 1.6 for PM-Mu) are observed in autumn, indicating that estimated ET is greater than the rainfall by 20-30 mm on average. This is also clearly shown by the large gap in Figure 4b. Both methods PM-Mu and PT-CMRS actually did not use the input data directly related to water balance (as rainfall or groundwater storage variations) or to water stress (as thermal infrared data) but mainly relied on the evolution of vegetation as monitor from remote sensing data. Inclusion of thermal data could have picked some low and high moisture conditions. LAI alone may not capture the variation on soil moisture, typically when the bare soil evaporation becomes significant after rainfall or under extreme dry conditions. By comparison, constrained by rainfall in its algorithm, mean seasonal ETs estimated by AWRA and GRACE were closer to the rainfall levels, with their ratios falling into a range of 0.7-1.2.  Table 2 shows the mean seasonal ET estimates and corresponding rainfall averaged over the period 2003-2010 and the fractions of ET to rainfall. The four methods used to estimate ET in this study were able to capture the seasonal ET patterns in each basin. LEB is much drier (has little rainfall against high potential evaporation requirement) than the MDB. Therefore, it is likely that the ET fluxes in the LEB is lower than in the MDB through the four seasons. However, in a comparison of the magnitude of seasonal ET to rainfall, the four methods behaved inconsistently. Again, seasonal PM-Mu estimates were constantly lower than other ET datasets in MDB, with the lowest ratio of 0.4 found in winter. The ratios of ET/P are frequently above 1.0 in PT-CMRS and PM-Mu during the dry seasons in LEB; the largest ratios (1.5 for PT-CMRS and 1.6 for PM-Mu) are observed in autumn, indicating that estimated ET is greater than the rainfall by 20-30 mm on average. This is also clearly shown by the large gap in Figure 4b. Both methods PM-Mu and PT-CMRS actually did not use the input data directly related to water balance (as rainfall or groundwater storage variations) or to water stress (as thermal infrared data) but mainly relied on the evolution of vegetation as monitor from remote sensing data. Inclusion of thermal data could have picked some low and high moisture conditions. LAI alone may not capture the variation on soil moisture, typically when the bare soil evaporation becomes significant after rainfall or under extreme dry conditions. By comparison, constrained by rainfall in its algorithm, mean seasonal ETs estimated by AWRA and GRACE were closer to the rainfall levels, with their ratios falling into a range of 0.7-1.2.   In Figure  . Four ET methods responded to the climate extremes dynamically. The smaller the anomalies were, the more similar the modelled ET estimates were to their multi-year average, and vice versa. Small anomalies were found in PT-CMRS and PM-Mu, which implies that ET estimates reproduced by these two energy balance constrained models do not apparently vary during the extremely wet or dry spells. AWRA and GRACE performed sensitively to climatic extremes. In spite of some spikes, the GRACE-based ET product tended to provide a relatively high/low ET with a magnitude consistent with rainfall.

Seasonal Variations
In Figure 5 Four ET methods responded to the climate extremes dynamically. The smaller the anomalies were, the more similar the modelled ET estimates were to their multi-year average, and vice versa. Small anomalies were found in PT-CMRS and PM-Mu, which implies that ET estimates reproduced by these two energy balance constrained models do not apparently vary during the extremely wet or dry spells. AWRA and GRACE performed sensitively to climatic extremes. In spite of some spikes, the GRACE-based ET product tended to provide a relatively high/low ET with a magnitude consistent with rainfall.

Inter-Annual Variations
Inter-annual comparisons were performed at basin-scale to test the response of the four ET datasets. Figure 4 shows monthly ET estimated by the four ET methods over the MDB and LEB from 2003 to 2010 (GRACE ET starts from August 2003), against the rainfall and potential evapotranspiration (Ep) for the same period. At monthly time steps, the correlation coefficients (R) were computed to measure the agreement between each ET pair as well as their sensitivity to rainfall.
Monthly Ep rates (ranging from 50 to 250 mm/month in MDB and from 80 to 250 mm/month in LEB) were always higher than the rainfall levels throughout the study period, indicating water-limited conditions. In that case, ET fluxes in our study areas is predominantly controlled by water availability rather than by energy. Two energy-based models show a quite low linear correlation with rainfall; RPT-CMRS = 0.35 and RPM-Mu = 0.41 for the MDB, and RPT-CMRS = 0.51 and RPM-Mu

Inter-Annual Variations
Inter-annual comparisons were performed at basin-scale to test the response of the four ET datasets. Figure 4 shows monthly ET estimated by the four ET methods over the MDB and LEB from 2003 to 2010 (GRACE ET starts from August 2003), against the rainfall and potential evapotranspiration (E p ) for the same period. At monthly time steps, the correlation coefficients (R) were computed to measure the agreement between each ET pair as well as their sensitivity to rainfall.
Monthly E p rates (ranging from 50 to 250 mm/month in MDB and from 80 to 250 mm/month in LEB) were always higher than the rainfall levels throughout the study period, indicating water-limited conditions. In that case, ET fluxes in our study areas is predominantly controlled by water availability rather than by energy. Two energy-based models show a quite low linear correlation with rainfall; R PT-CMRS = 0. 35  Remarkable differences can be found during wet (ET PT-CMRS and ET PM-Mu < P) or dry spells (ET PT-CMRS and ET PM-Mu > P; particularly in LEB). By comparison, AWRA has high correlation with rainfall in both basins, 0.7 < R < 0.8. AWRA and GRACE ET estimates were closer to rainfall levels (especially GRACE) due to the fact that they were forced by rainfall data. They also had a good consistency with each other: R MDB = 0.67 and R LEB = 0.69. Although PT-CMRS had a higher linear relationship with AWRA (0.77 ≤ R MDB/LEB ≤ 0.79), the amplitude of the PT-CMRS was significantly lower than AWRA during the extreme seasons. In general, the two energy-balanced products exhibited similar patterns at monthly time steps, and displayed large differences with AWRA and GRACE-based ET in dry and wet periods. PM-Mu was systematically lower than other ET sources in MDB. Figure 6 presents the capability of each ET dataset in capturing ET variations in different hydrological years. PM-Mu constantly provided the lowest ET estimates over the MDB; a common short of~100 mm/year was found when compared to PT-CMRS. Compared to AWRA and GRACE estimates, PT-CMRS provided larger values over the MDB in dry hydrological years but lowest values during wet years (e.g., August 2009-July 2010; impacted by La Niña event). ET overestimation by the two energy balance constrained models was commonly observed in the LEB, with annual ET flux twice or larger than the rainfall during the extremely dry hydrological year of August 2004-July 2005 and August 2007-July 2008. Moreover, these two estimates did not exhibit any response to rain variations from year to year and provided almost constant annual values. Constrained by the water balance equation, AWRA and GRACE showed a good sensitivity to annual rainfall variations both in timing and magnitude. GRACE-derived ET estimates at annual scales showed larger inter-annual variations than AWRA over the MDB (with ET larger than rain for some years), while very similar results were obtained over the LEB. = 0.37 for the LEB. Remarkable differences can be found during wet (ETPT-CMRS and ETPM-Mu < P) or dry spells (ETPT-CMRS and ETPM-Mu > P; particularly in LEB). By comparison, AWRA has high correlation with rainfall in both basins, 0.7 < R < 0.8. AWRA and GRACE ET estimates were closer to rainfall levels (especially GRACE) due to the fact that they were forced by rainfall data. They also had a good consistency with each other: RMDB = 0.67 and RLEB = 0.69. Although PT-CMRS had a higher linear relationship with AWRA (0.77 ≤ RMDB/LEB ≤ 0.79), the amplitude of the PT-CMRS was significantly lower than AWRA during the extreme seasons. In general, the two energy-balanced products exhibited similar patterns at monthly time steps, and displayed large differences with AWRA and GRACE-based ET in dry and wet periods. PM-Mu was systematically lower than other ET sources in MDB. Figure 6 presents the capability of each ET dataset in capturing ET variations in different hydrological years. PM-Mu constantly provided the lowest ET estimates over the MDB; a common short of ~100 mm/year was found when compared to PT-CMRS. Compared to AWRA and GRACE estimates, PT-CMRS provided larger values over the MDB in dry hydrological years but lowest values during wet years (e.g., August 2009-July 2010; impacted by La Niña event). ET overestimation by the two energy balance constrained models was commonly observed in the LEB, with annual ET flux twice or larger than the rainfall during the extremely dry hydrological year of August 2004-July 2005 and August 2007-July 2008. Moreover, these two estimates did not exhibit any response to rain variations from year to year and provided almost constant annual values. Constrained by the water balance equation, AWRA and GRACE showed a good sensitivity to annual rainfall variations both in timing and magnitude. GRACE-derived ET estimates at annual scales showed larger inter-annual variations than AWRA over the MDB (with ET larger than rain for some years), while very similar results were obtained over the LEB.  A Budyko diagram was used to confirm that ET fluxes in our two basins were determined by water availability rather than the energy factor. Estimates were plotted in the Budyko diagrams to examine the variation of the water closure property across hydrological years (Section 2.2.4).
The aridity indices in the Budyko diagrams (Figure 7) confirm that LEB is more arid than MDB (LEB: 4.0 < E p /Pa < 11.5; MDB: 2.7 < E p /Pa < 5.0). In MDB, Ea/Pa ratios estimated by PM-Mu were sitting in a range 0.5-0.8, demonstrating a systematically low percentage of rainfall transforming into ET. For the same basin, PT-CMRS had some Ea/Pa values lying beyond y = 1 in some dry years (larger E p /Pa), manifesting a slight outperformance of PT-CMRS ET during drought conditions. Such a phenomenon was more apparent in the LEB. The dryer the year in the LEB (E p /Pa > 10), the more unrealistically the PM-Mu and PT-CMRS estimates; Ea/Pa ratios reached 1.8 (PM-Mu) and 1.7 (PT-CMRS) when LEB received the lowest annual rainfall (158.6 mm/year) from August 2007 to July 2008. By comparison, Ea/Pa predicted by the two water-balance-constrained methods (GRACE and AWRA) were much closer to the asymptotic curve. AWRA and GRACE predicted 85.1% and 97.5% of rainfall lost via evaporation in both basins during the wettest year August 2009-July 2010, respectively; that is to say, 15% (from AWRA) and 2.5% (from GRACE) of rainfall would be converted into runoff and exit the basin or be stored in the aquifer. A Budyko diagram was used to confirm that ET fluxes in our two basins were determined by water availability rather than the energy factor. Estimates were plotted in the Budyko diagrams to examine the variation of the water closure property across hydrological years (Section 2.2.4).
The aridity indices in the Budyko diagrams (Figure 7) confirm that LEB is more arid than MDB (LEB: 4.0 < Ep/Pa < 11.5; MDB: 2.7 < Ep/Pa < 5.0). In MDB, Ea/Pa ratios estimated by PM-Mu were sitting in a range 0.5-0.8, demonstrating a systematically low percentage of rainfall transforming into ET. For the same basin, PT-CMRS had some Ea/Pa values lying beyond y = 1 in some dry years (larger Ep/Pa), manifesting a slight outperformance of PT-CMRS ET during drought conditions. Such a phenomenon was more apparent in the LEB. The dryer the year in the LEB (Ep/Pa > 10), the more unrealistically the PM-Mu and PT-CMRS estimates; Ea/Pa ratios reached 1.8 (PM-Mu) and 1.7 (PT-CMRS) when LEB received the lowest annual rainfall (158.6 mm/year) from August 2007 to July 2008. By comparison, Ea/Pa predicted by the two water-balance-constrained methods (GRACE and AWRA) were much closer to the asymptotic curve. AWRA and GRACE predicted 85.1% and 97.5% of rainfall lost via evaporation in both basins during the wettest year August 2009-July 2010, respectively; that is to say, 15% (from AWRA) and 2.5% (from GRACE) of rainfall would be converted into runoff and exit the basin or be stored in the aquifer.

Discussion
Catchment/basin water balance method provides the simplest way to solve ET estimation by meeting only three water budget components as its input variables. However, the component of ∆S is usually the hardest part to access or measure. In most cases, the common way to address the absence of ∆S is to regard it as negligible over a long-term period [31,32,65]. But when it comes to the annual or inter-annual scales, neglecting of ∆S would lead to an imbalance of water budget equation [66][67][68]. In fact, short-term ∆S can be a crucial component in water budget, particularly when a basin meets climate extremes. Time series of ∆S anomaly over the MDB and the LEB exhibits large annual (several tenths of mm of equivalent water height) and inter-annual variations (with a decrease from 2003 to 2009 and an increase since) between 2003 and 2010 ( Figure 8). Table 3

Discussion
Catchment/basin water balance method provides the simplest way to solve ET estimation by meeting only three water budget components as its input variables. However, the component of ∆S is usually the hardest part to access or measure. In most cases, the common way to address the absence of ∆S is to regard it as negligible over a long-term period [31,32,65]. But when it comes to the annual or inter-annual scales, neglecting of ∆S would lead to an imbalance of water budget equation [66][67][68]. water loss or gain in ∆S would take up 10-20% of total annual rainfall for both basins. By including GRACE-estimated ∆S, the closure of the basin water balance equation is improved; what's more, the phase and amplitude of the annual and seasonal cycle of ET can be ascertained [69].
Water 2017, 9,614 14 of 19 equation is improved; what's more, the phase and amplitude of the annual and seasonal cycle of ET can be ascertained [69].  In terms of the uncertainty of GRACE-based ET estimates, they are largely dependent on the quality of GRACE TWSA and thus computed ∆S. A previous study using regional GRACE solutions over Australia found an uncertainty of 19.1 mm when computing the standard deviation of the TWSA over a xeric region [52]. Another major source of uncertainty associated with GRACE ET estimates is attributed to the quality of rainfall data. Usually, arid and semi-arid basins have poor rainfall monitoring networks. In our case, MDB has much denser monitoring networks than LEB does; the latter has a large unpopulated region without rainfall stations [63]. At basin scales, the RMSE calculated between BoM and TRMM 3B43 rainfall datasets are 3.0 and 27.4 mm/month for MDB and LEB, respectively. By blending gauge records with ancillary rainfall data, such as radar and satellite measurements or an ensemble method, may improve the accuracy of rainfall input data [63,64].
Optical satellite imagery allows landscape conditions such as vegetated surface, biome types, and surface standing waters to be distinguishable at unprecedentedly high spatiotemporal resolution (e.g., 0.05° and 8/16 days for MODIS), which facilitates an enhancement in simulating large-scale ET processes. For example, an integration of RS vegetation index (VIs, including LAI, NDVI, and EVI) allows for a separation between plant transpiration and soil evaporation, as well as for adding the canopy rainfall interception component into ET modelling. In our study, PM-Mu employs a complex and process-based scheme to address vegetation transpiration and soil evaporation separately; the new version sets up numerous thresholds to define canopy stomatal conditions over biomes and wet canopy [19]. PT-CMRS combines the EVI and GVMI indexes into a simple and dynamic scaling strategy of the Priestley-Taylor framework to deal with transpiration or evaporation from vegetation, open water bodies, and bare soil; it does consider canopy rainfall interception but merely as a scaled precipitation component. In addition, PT-CMRS ET estimates  In terms of the uncertainty of GRACE-based ET estimates, they are largely dependent on the quality of GRACE TWSA and thus computed ∆S. A previous study using regional GRACE solutions over Australia found an uncertainty of 19.1 mm when computing the standard deviation of the TWSA over a xeric region [52]. Another major source of uncertainty associated with GRACE ET estimates is attributed to the quality of rainfall data. Usually, arid and semi-arid basins have poor rainfall monitoring networks. In our case, MDB has much denser monitoring networks than LEB does; the latter has a large unpopulated region without rainfall stations [63]. At basin scales, the RMSE calculated between BoM and TRMM 3B43 rainfall datasets are 3.0 and 27.4 mm/month for MDB and LEB, respectively. By blending gauge records with ancillary rainfall data, such as radar and satellite measurements or an ensemble method, may improve the accuracy of rainfall input data [63,64].
Optical satellite imagery allows landscape conditions such as vegetated surface, biome types, and surface standing waters to be distinguishable at unprecedentedly high spatiotemporal resolution (e.g., 0.05 • and 8/16 days for MODIS), which facilitates an enhancement in simulating large-scale ET processes. For example, an integration of RS vegetation index (VIs, including LAI, NDVI, and EVI) allows for a separation between plant transpiration and soil evaporation, as well as for adding the canopy rainfall interception component into ET modelling. In our study, PM-Mu employs a complex and process-based scheme to address vegetation transpiration and soil evaporation separately; the new version sets up numerous thresholds to define canopy stomatal conditions over biomes and wet canopy [19]. PT-CMRS combines the EVI and GVMI indexes into a simple and dynamic scaling strategy of the Priestley-Taylor framework to deal with transpiration or evaporation from vegetation, open water bodies, and bare soil; it does consider canopy rainfall interception but merely as a scaled precipitation component. In addition, PT-CMRS ET estimates heavily rely on ground calibration (none of the flux tower sites are distributed in inland dry areas) [53]. Obviously, there is a tradeoff between accuracy and parsimoniousness in ET modelling, as well as a consideration of data availability from ground observations.
In arid and semiarid environments, soil evaporation often makes up the majority of the total ET due to low vegetation coverage. This means that ET models should be able to reflect the relation with soil moisture conditions. Land surface models like AWRA partition precipitation to soil moisture and groundwater systems (other LSMs may not), and thus put a soil moisture constraint onto ET. However, the soil water constraint of energy-based ET products, such as PM-Mu and PT-CMRS, tends to be weak; most models tend to determine ET fluxes by energy factor, in particular based on net radiation estimation rather than by water constraint. Long et al. [16] mentioned that such RS-based ET models only make soil moisture implicitly linked to VIs and atmospheric variables. PM-Mu has soil moisture information indirectly linked to the LAI/NDVI and vapor pressure deficit (VDP); so does PT-CMRS, replaced by EVI and GVMI. These models implicitly assume that vegetation develops in agreement with water availability, which may be true over a long term or "standard" climatic period, but can fail in conditions when drought or rainfall present situations far from average. These models also rely on parameters (e.g., describing stomatal response) that are difficult to assess at large scale, as well as on the accuracy of land use maps. In semi-arid and arid regions, such as MDB and LEB, rain may fall and evaporate directly from bare soil, hence without major influence on vegetation indices and derived ET estimates. That is a potential explanation as to why satellite ET datasets like PT-CMRS and PM-Mu are less sensitive to rainfall/soil moisture and cannot balance well the water budget in these arid and semi-arid environments. Other RS-based ET models rely on the use of thermal infrared data that provide a strong link to water stress and moisture availability and which may improve significantly energy balance estimates. However, today there are no available operational ET products based on thermal infrared, in particular because such models are difficult to implement over large areas such as a continent or even a country. They either rely on a highly accurate characterization of spatial variation of climatic variables (in particular air temperature) or require a very homogeneous climatic zone for implementation.

Conclusions
This study, examines the dynamics of continental ET products in water-limited environments. Four ET datasets derived from PT-CMRS, PM-Mu, AWRA, and GRACE were compared in the Murray-Daring and Lake Eyre Basins, against rainfall variations and during climate extremes. Two energy balance constrained ET methods, PT-CMRS and PM-Mu, which are forced by optical satellite retrievals, have poor response to high water variability; they provided unrealistically high ET values (beyond the rainfall) during the drought and low values when rainfall became exceptionally high. This problem can be attributed to the lack of water constraint following water balance. In contrast, AWRA and GRACE are forced by rainfall data, demonstrating their dynamics in addressing the spatial/temporal rainfall variability over water-limited environments. Absolute error associated with each ET product is not directly measurable but a measure of error may be established via a basin-scale comparison with the ensemble mean.
Our results imply that: (1) current ET models based on the energy balance and derived using optical data are not accurate over water-limited areas and that they may need to include further water constraint(s) for ET estimation over arid and semi-arid regions; (2) GRACE observations provide a valuable tool for quantifying basin-scale water storage change at annual or sub-annual scales, as well as an independent source for large-scale ET mapping and validation; (3) a promising way to enhance ET estimation over large water-limited basins/regions may rely on merging high-resolution (but poor in basin-water-balance-constrained) optical RS-based ET products with well water-balance-constrained (but coarse in spatial resolution) GRACE ET estimates [20,70].