Linking the complementary evaporation relationship with the Budyko framework for ungauged areas in Australia

. While the calibration-free complementary relationship (CR) has performed excellently in predicting terrestrial evapotranspiration (ET a ) , how to determine the Priestley–Taylor coefﬁcient ( α e ) is a remaining question. In this work, we evaluated this highly utilizable method, which only requires atmospheric data, with in situ ﬂux observations and basin-scale water-balance estimates (ET wb ) in Australia, proposing how to constrain it with a traditional Budyko equation for ungauged locations. We found that the CR method with a constant α e transferred from fractional wet areas performed poorly in reproducing the mean annual ET wb in un-regulated river basins, and it underperformed advanced physical, machine-learning, and land surface models in closing grid-scale water balance. This problem was remedied by linking the CR method with a traditional Budyko equation that allowed for an upscaling of the optimal α e from gauged basins to ungauged locations. The combined CR–Budyko framework enabled us to reﬂect climate


Introduction
Evapotranspiration (ET a ) plays a pivotal role in water and energy exchanges between the land and the atmosphere. On the global scale, more than 60 % of terrestrial precipitation (P ) returns to the atmosphere through plants' vascular systems and soil pores while consuming over 70 % of surface net radiation (Trenberth et al., 2007(Trenberth et al., , 2009). Since it is tightly coupled with carbon cycles, abnormally low ET a would indicate food insecurity and low ecosystem sustainability (Jasechko, 2018;Kyatengerwa et al., 2020;Pareek et al., 2020;Swann et al., 2016). In severe cases, ET a limited by deficient soil moisture can lead to extreme heatwaves that further propagate the water deficit in space and time (Miralles et al., 2014;Mueller and Seneviratne, 2012;Schumacher et al., 2022).
Despite great community efforts for sharing in situ observations (e.g., Baldocchi, 2020;Novick et al., 2018), ET a gauging networks are unevenly established over land surfaces and often subjected to error sources (e.g., unclosed energy balance) and limited data lengths (Ma et al., 2021). Inevitably, modeling approaches are needed to predict ET a in ungauged or poorly gauged areas or to characterize it on a long timescale in a large area. Hence, various approaches have been proposed, including physical models (e.g., Martens et al., 2017;Zhang et al., 2016), machinelearning techniques (e.g., Jung et al., 2019;Tramontana et al., 2016), and conceptual land surface schemes (e.g., Guimberteau et al., 2018;Haverd et al., 2018).
Published by Copernicus Publications on behalf of the European Geosciences Union. 5956 D. Kim et al.: Linking the complementary evaporation relationship Those modeling approaches typically require P data and land surface information (e.g., remote-sensing vegetation indices) to quantify available soil moisture to the vaporization process. However, due in part to uncertainty associated with P data (Sun et al., 2018) and model structures (Samaniego et al., 2017;Zhang et al., 2019), resulting ET a estimates have shown substantial disparities. In the comprehensive intercomparison by Pan et al. (2020), for example, the 14 advanced land surface models generated the global mean ET a varying widely between 450 and 700 mm a −1 . Such a large incongruity in modeled ET a was also found by the earlier Global Soil Wetness Project (Schlosser and Gao, 2010), suggesting that an alternative method is necessary to circumvent the uncertainty sources.
A practical method to simulate ET a without P data and land surface schemes is the complementary relationship (CR) of evaporation (Bouchet, 1963). It uses the evident fact that the air over a water-limited surface amplifies its vapor pressure deficit (VPD), while this effect disappears when the same surface is amply wet (Chen and Buchberger, 2018;Ramírez et al., 2005;Zhou et al., 2019). Based on the atmospheric self-adjustment, numerous equations have been formulated to predict ET a only using routine meteorological data (e.g., Anayah and Kaluarachchi, 2014;Crago and Crowley, 2005;Crago and Qualls, 2013;Hobbins et al., 2004;Huntington et al., 2011;Kahler and Brutsaert, 2006 among others). In particular, the definitive derivation by Brutsaert (2015) and the following modifications (Crago et al., 2016;Crago and Qualls, 2021;Szilagyi, 2021;Szilagyi et al., 2017) provided strong physical foundations to the early principle of Bouchet (1963). They have excellently predicted ET a at various spatial and temporal scales (e.g., Brutsaert et al., 2017Brutsaert et al., , 2020Crago and Qualls, 2018;Ma et al., , 2021 and allowed users to assess vegetation droughts over national and continental areas (e.g., Kim et al., 2019Kim et al., , 2021Kyatengerwa et al., 2020).
Nevertheless, definitive CRs still require at least some ET a data to calibrate the parameters that determine the hypothetical wet-surface evaporation (ET w ; Qualls and Crago, 2020); thus, they are not fully free of P data or parameterization. For instance, Brutsaert et al. (2020) calibrated the single parameter of the CR of Brutsaert (2015) with flux observations and basin-scale P and runoff (Q) data to estimate mean annual ET a across the globe. For evaluating four definitive CRs from the derivation of Brutseart (2015), Crago et al. (2022) also calibrated their parameters against eddy-covariance flux observations. To date, Szilagyi et al. (2017) have proposed the only CR formulation that purely uses routine meteorological data; however, it depends on a questionable assumption that the parameter for ET w is constant over a large continental area, being counterfactual to experimental studies on the Priestley and Taylor (1972) coefficient (e.g., Assouline et al., 2016;Baldocchi et al., 2016;Parlange and Katul, 1992;Wang et al., 2014). Given the complex space-time links between climate, soil, and vegetation (Hagedorn et al., 2019;Mekonnen et al., 2019;Rodriguez-Iturbe, 2000), the aerodynamic component of ET w is unlikely represented by a fixed fraction of the net radiation.
Owing to the data required for parameter calibration, the state-of-the-art CR formulations might not be applicable in ungauged locations. In part, this problem can be mended by an additional constraint for determining the essential parameters, and the traditional Budyko framework can come into play. A Budyko function (e.g., Fu, 1981;Yang et al., 2008) explains the mean ratio of ET a to P (i.e., surface water balance) simply by climatological aridity and a few implicit parameters, simultaneously closing the surface energy budget (Mianabadi et al., 2020). Although Bouchet's principle has often been linked with the water balance described by Budyko functions (e.g., Carmona et al., 2016;Chen and Buchberger, 2018;Lhomme and Moussa, 2016;Zhang and Brutsaert, 2021), this theoretical link has been ignored when predicting ET a by the definitive CRs. Kim and Chun (2021) explicitly showed that the atmospheric selfadjustment is tightly coupled with the climatological aridity within a Budyko function. This implies that the optimal parameter for a definitive CR should vary with climates rather than staying constant.
In this work, we showed that a Budyko equation could become an important physical constraint when predicting ET a by a definitive CR over a continental area. Here, a practical approach was proposed to determine the parameter reasonably in ungauged locations via a case study for the Australian continent, where the performance of the CR method remained unknown in many parts. Based on the analytical relationship between the CR and the Budyko framework, we showed why the parameter of the CR is not independent of local climate conditions and addressed how to reflect spatially varying climates in its essential parameter.
2 Methodology and data 2.1 The polynomial CR by Szilagyi et al. (2017) For the case study, we employed the calibration-free CR formulated by Szilagyi et al. (2017). It describes the atmospheric self-adjustment to surface moisture conditions using three evaporation rates, namely, ET a , ET w , and the potential evaporation (ET p ). ET a is the actual moisture flux from a land surface to the atmosphere, and ET w is the hypothetical ET a rate that should occur with ample water availability. ET p is the atmospheric capacity to receive water vapor that responds actively to soil moisture conditions. By defining the two dimensionless variables as x ≡ ET w /ET p and y ≡ ET a /ET p , Szilagyi et al. (2017) derived a polynomial function from four definitive boundary conditions. Under ample water conditions, ET p does not deviate from ET w and ET a (i.e., ET p = ET w = ET a ); hence, the corresponding zero-order boundary condition is (i) y = 1 for x = 1. In contrast, ET a must be nil over a desiccated surface (i.e., y = 0), and by energy balance, the surface net radiation should be fully transformed to the sensible heat flux. Then, the atmospheric VPD would be amplified at the maximum level with the same net radiation and wind speed. Defining the maximum ET p rate as E pmax , another zeroorder boundary condition is given as (ii) y = 0 for x = x min ≡ ET w /E pmax . When x = 1 (i.e., ample water), changes in ET a would be controlled by changes in ET w , yielding a first-order boundary condition as (iii) dy/dx = 1 for x = 1. Over a desiccated surface, ET a stays at zero, even when ET w changes; thus, another first-order boundary condition becomes (iv) dy/dx = 0 for x = x min . The simplest polynomial equation satisfying the four boundary conditions is where X rescales the variable x into [0, 1] as Equation (1) allows users to estimate ET a with no land surface information because ET p , ET w , and E pmax are all obtainable from a set of net radiation, air temperature, dew-point temperature, and wind speed data. ET p and E pmax can be estimated by the Penman (1948) equation: where (·) is the slope of the saturation vapor pressure curve (kPa • C −1 ); T a is the mean air temperature ( • C); γ is the psychrometric constant (kPa • C −1 ); R n is the surface net radiation less the soil heat flux (MJ m −2 d −1 ); λ v is the latent heat of vaporization (MJ kg −1 ); f u = 2.6 (1 + 0.54 u 2 ) is the Rome wind function (mm d −1 kPa −1 ), where u 2 is the 2 m wind speed (m s −1 ); and VPD is calculated by e s (T a ) minus e s (T dew ), where e s (·) is the saturation vapor pressure (kPa), and T dew is the dew-point temperature ( • C). T dry in Eq. (4) is the air temperature ( • C) at which the lower atmosphere is devoid of humidity presumably by the adiabatic drying process: where T wb is the wet-bulb temperature ( • C) at which the saturation vapor pressure curve intersects with the adiabatic wetting line. Thus, it is obtained by To estimate ET w in Eq.
(2), the Priestley-Taylor (1972) equation has been a typical choice (e.g., Brutsaert, 2015;Crago et al., 2016;Han and Tian, 2018;Szilagyi et al., 2017): where α e is the Priestley-Taylor coefficient ranging usually within [1.10, 1.32] (Szilagyi et al., 2017), and T w is the wetenvironment air temperature ( • C). T w can be approximated with the wet-surface temperature (T ws ) because the vertical air temperature gradient is negligible under a wet environment. Given its independence on areal extent (Szilagyi and Schepers, 2014), T ws can be approximated by the implicit Bowen ratio (β) of a small wet patch: Equation (8) assumes that the available radiation for the wet patch is close to that of the drying surface (Szilagyi et al., 2017). T ws might be higher than T a when the air is close to saturation. In such a case, T ws should be capped by T a when calculating ET w . The single parameter of the polynomial CR, i.e., α e , is analytically obtainable by inserting the Priestley-Taylor equation into the Bowen ratio of a wet environment (Szilagyi et al., 2017) as where α e must fall within the theoretical limit of [1, 1 + γ / (T a )] (Priestley and Taylor, 1972).

The analytical relationship between the polynomial CR and a Budyko function
Since Eq. (9) is only applicable in a wet environment, Szilagyi et al. (2017) identified wet locations in a continental area based on the fact that the air close to saturation should have high relative humidity (RH) with T ws > T a . Thus, they calculated α e values at locations with RH > 90 % and T ws > T a + 2 • C, and the average value was used to predict ET a for a continental area. However, the spatially constant α e is unlikely suitable in such a large area under diverse climates because the equilibrium between the atmosphere and the underlying surface is intertwined with the partitioning of P to ET a and Q over the surface. Kim and Chun (2021) analytically related Eq. (1) with the traditional Turc-Mezentsev equation and found that the selfadjustment of ET p (i.e., x) is tightly linked with climatological aridity and land properties. For the independence between P and "the possible maximum ET a " of the Budyko framework, Kim and Chun (2021) reformulated the traditional equation with 0 ≡ ET w /P instead of the commonly 5958 D. Kim et al.: Linking the complementary evaporation relationship used aridity index ( ≡ ET p /P ) as where the parameter n implicitly represents the factors affecting the P partitioning other than the climatic drivers. When dividing Eq. (10) by , it is found that the Budyko Eq. (10) is intertwined with the Eq. (1) of the CR: Equation (11) implies that the self-adjustment of ET p (i.e., x) is tightly related with the climatic condition (i.e., ) and the implicit land property (i.e., n).
While the x and n can be achievable from a set of ET a , ET p , E pmax , and P values by inverting Eq. (11), such an approach is not applicable in locations with no ET a data. To quantify x values only using ET p , E pmax , and P , Kim and Chun (2021) developed a regression equation between x and , x min , and n values from 513 gauged river basins around the world. We used the same regression-based regionalization. Considering x min = xET p /E pmax , the nonlinear expression in Eq. (11) can be approximated by a multiple regression as wherex is the approximate ratio of ET w to ET p , and b 0 , b 1 , and b 2 are the intercept and the regression coefficients. Since the implicit parameter n is unavailable in ungauged locations, Eq. (12) needs to be further simplified by neglecting the last term: where c 0 , c 1 , and c 2 are the intercept and the coefficients of the approximated regression. Ifx is known by the regression by Eq. (13), the parameter α e can be estimated using the Priestley-Taylor equation as whereα e is the Priestley-Taylor coefficient that approximately satisfies the CR and the Budyko equations together, and ET eq is the equilibrium ET a (mm d −1 ) at which VPD is nil under a wet environment. It should be noted that P , ET p , E pmax , and ET eq within Eqs. (10)-(13) must be on a timescale where the Turc-Mezentsev equation is valid (typically longer than a year), andα e is still bounded within [1, 1 + γ / (T a )].

Atmospheric forcing, eddy-covariance, and runoff data
We examined the CR-Budyko combined framework in the Australian continent lying within (10-45 • S, 113-155 • E). The required atmospheric forcing data (R n , T a , T dew , and u 2 ) were collected from the advanced ERA5-Land reanalysis archive (Muñoz-Sabater et al., 2021) of the European Centre for Medium-Range Weather Forecasts (https://cds. climate.copernicus.eu; last access: 10 December 2021). The monthly averages of surface latent and sensible heat fluxes, 2 m air temperature, 2 m dew-point temperature, and 10 m U and V wind speed components at 0.1 • × 0.1 • were downloaded for 1981-2020. R n was calculated by summing the two heat fluxes, and the 10 m wind speed components were converted to u 2 using the logarithmic wind profile (Allen et al., 1998).
We also collected the Australian edition of the Catchment Attributes and Meteorology for Large sample Studies (CAMELS; Fowler et al., 2021) series of datasets (https://doi.org/10.1594/PANGAEA.921850). The CAMELS datasets comprise daily time series of 19 hydrometeorological variables at 222 unregulated river basins in Australia up to 2014, and we selected the 71 basins larger than 500 km 2 to contain at least five CR ET a estimates within the boundaries. The water-balance ET a (ET wb ) (i.e., ET wb ≈ P − Q) of each basin was calculated for the two periods of 1981-1997 and 1998-2014. The mean annual ET wb for the former period was used for the regressions with Eqs. (12) and (13), and the predicted ET a was evaluated against the latter.
As a point-scale evaluation dataset, the annual flux observations were taken from the 15 eddy-covariance stations (Table 1) included in the FLUXNET2015 archive (https:// fluxnet.org/; last access: 1 July 2021). We chose the flux towers with two or more annual means and adopted the energybalance-corrected latent heat flux observations with the quality measures "LE_F_MDSQC" higher than 0.70. Given the fine resolution of the ERA5-Land forcing data, we believed that the ET a estimates by CR could be directly compared with the point-scale observations.
In addition, as a grid-scale evaluation reference, the SILO P data at 0.01 • × 0.01 • were collected from the Queensland government (https://www.longpaddock.qld.gov.au/silo/ gridded-data; last access: 1 June 2021) together with the Global RUNoff (GRUN) ENSEMBLE (Ghiggi et al., 2021) (https://doi.org/10.6084/m9.figshare.12794075; last access: 1 October 2021). The global Q data were produced at 0.5 • × 0.5 • using a machine-learning algorithm trained by in situ streamflow observations, and potential biases were reduced by simulations with 21 sets of atmospheric forcing (Ghiggi et al., 2021). The SILO P was used to calculate at each grid of the forcing data. After bilinearly unifying the resolutions of SILO P and GRUN Q data, we also calculated  Against the grid-scale ET wb estimates, performance of the polynomial CR was also compared with three ET a products from a physical, a machine-learning, and a land surface model. The physical model was the Global Land Evaporation Amsterdam Model (GLEAM) v3.2 (Martens et al., 2017; https://www.gleam.eu; last access: 3 June 2020) based on the Priestley-Taylor equation constrained by microwave-derived soil moisture, surface temperature, and vegetation optical depth. The machine-learning ET a product was the FluxCom (http://www.fluxcom.org/; last access: 18 March 2019) that upscaled in situ observations at 224 eddy-covariance towers using 11 algorithms (Jung et al., 2019). We used the version forced by the CRUNCEPv8 that has the longest data length from 1950 to 2016. The land surface model product was the ERA5-Land monthly ET a (https://cds.climate.copernicus.eu; last access: 7 July 2021) simulated by the advanced Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (Balsamo et al., 2015). All the modeled ET a datasets were bilinearly regridded to 0.5 • × 0.5 • for 1998-2014 to be compared with the grid-scale ET wb data. Figure 1a depicts the spatial distribution of the inverted aridity index ( −1 = P /ET p ) that can traditionally categorize climate conditions. The mean ratios between SILO P and ET p for 1998-2014 indicated that 83 % of the Australian land surfaces were under arid (0.05 < −1 < 0.2) and semi-arid climates (0.2 < −1 < 0.5). Semi-humid (0.5 < −1 < 0.65) and humid climates ( −1 > 0.65) were only found in the northern and southeastern coastal areas and the southwestern edge where major cities and agricultural lands have developed. Despite the high aridity, hyper-arid climates ( −1 > 0.05) were not found in Australia.

Performance of the calibration-free CR in Australia
We first examined the calibration-free approach by Szilagyi et al. (2017)    with RH > 90 % and T ws > T a + 2 • C, at which the α e values from Eq. (9) were within 1.15 ± 0.047 (mean ± standard deviation). Though the two conditions were met in some mountainous areas in the southeastern part, we excluded them because unexpectedly high α e values were obtained. The mean α e = 1.15 fell within the theoretical limits and was equal to the value used in the prior studies in China  and the conterminous United States .
Using the CR with α e = 1.15, we predicted ET a over the entire Australian continent (Fig. 1b). The distribution of the resulting mean ET a for 1998-2014 was coherent with that of −1 . The mean CR ET a ranged in 262 ± 85.3 and 547 ± 173 mm a −1 under arid and semi-arid climates, respectively. On the other hand, CR ET a in semi-humid and humid locations was much higher, at 886 ± 187 mm a −1 and 1010 ± 213 mm a −1 , respectively. The calibration-free CR predicted the continental mean ET a to be as high as 489 mm a −1 for 1981-2012, and it was about 11.3 % higher than the estimate for the same period (439 mm a −1 ) by Zhang et al. (2016). The mean fraction of ET a to P for 1998-2014 (97 %) was larger than the typical ET a value in Australia (∼ 90 %; Glenn et al., 2011), implicating that the constant α e = 1.15 seemed to make the CR overrate ET a .
The overestimation of the calibration-free CR was confirmed by the flux observations and the basin-scale ET wb (Fig. 2). The percent bias (p-bias) of CR ET a to the pointscale annual ET a was +10.4 %, while it became more than doubled when compared to the basin-scale ET wb . Though the Pearson correlation coefficients (Pearson r) were significantly high between the CR ET a and the two evaluation references, the low Nash-Sutcliffe efficiency (NSE) to ET wb implies that the CR method could perform poorly in wet river basins. The regression slopes in Fig. 2 also indicate that the calibration-free CR tends to increasingly overestimate as climate becomes wetter. The root mean square error (RMSE) of CR ET a to ET wb was higher than to the point observations. Although it appeared to perform acceptably at the 15 flux towers, the CR method produced considerable biases in the 71 CAMELS basins. The performance measures were not as excellent as the same CR method had shown in the United States (Ma et al., 2021;Kim et al., 2019) and in China .
One may argue that the mean α e derived from fractional wet areas is unlikely representative of the large Australian continent, and this might introduce the biases to CR ET a estimates. Hence, we re-simulated CR ET a with the estimate by Ma et al. (2021) (α e = 1.10) from a global-scale analysis. Figure 3a shows that the predicted ET a became nearly unbiased at the 15 flux tower locations and seemingly suggests that the decreased α e could become a solution to improving the CR method. Nevertheless, the fixed α e still made the CR overestimate ET a in the CAMELS basins under (semi-)humid climates, albeit slightly ameliorated (Fig. 3b).

The empirical relationship betweenx and climate conditions
Figures 2 and 3 imply that the calibration-free CR with a fixed α e was unlikely good at closing the local water balance in (semi-)humid river basins. To resolve this problem with the CR-Budyko framework, first we estimated the climato-logical x and the parameter n of the CAMELS basins using Eq. (11) with the mean annual ET wb , P , ET p , and E pmax for 1981-1997. Figure 4a-c illustrate the scatter plots between the resultant x and corresponding , ET p /E pmax , and n values. Pearson r values between the x and the other three variables were −0.88, −0.59, and 0.44, respectively (significant at the 1 % level), suggesting that the self-adjustment of ET p is not only correlated with climate conditions, but with land surface properties at least in part. By regressing between the x values and the log-transformed , ET p /E pmax and n, we obtained an empirical relationship that enables us to spatially predict the mean annual ratio of ET w to ET p as x =0.949 − 0.204 ln ( ) + 0.231 ln ET p /E pmax +0.0712 ln (n) .
The regression coefficients were all significant at the 1 % level, and the coefficient of determination (R 2 ) was 0.98. The regression equation was further approximated by discarding n from the explanatory variables: x = 1.023 − 0.220 ln ( ) + 0.210 ln ET p /E pmax .
The R 2 value of Eq. (17) declined to 0.93. We found that the simple regression between x and further reduced R 2 to 0.90. While the heterogeneous land properties exert nonnegligible influences, the regression analyses indicate that the climatic condition dominantly explains the spatial variation of the atmospheric self-adjustment. Equation (17) performed excellently in reproducing the x values from CR with and ET p /E pmax (Fig. 4d). The NSE, RMSE, Pearson r, and p-bias between the predictedx and the x from CR were 0.93, 0.03, 0.96, and 0.0 %, respectively.

Evaluation of the CR and the advanced models against the grid ET wb
By multiplyingx by the mean annual ratio between ET p and ET eq , we determinedα e across the Australian land surfaces.
The resultingα e values ranged within 1.13 ± 0.114, and the median value was almost equal to the global estimate (1.10) by Ma et al. (2021). They were relatively high in the northwestern and the northern part while being below the mean in the southern and the eastern parts (Fig. 5a). On 19 % of the surfaces,α e values were unity, and thus they might become below the theoretical limit unless bounded. We again generated CR ET a using the spatially varyingα e values (Fig. 5b). The mean CR ET a for 1998-2014 ranged in 249 ± 78.8 and 530 ± 172.0 mm a −1 under arid and semi-arid climates, while it decreased to 805.2 ± 209 and 932 ± 239 mm a −1 in semi-humid and humid regions, respectively. The flux observations were still acceptably regenerated with the less biases than in the case of α e = 1.15 (Fig. 6a). Theα e based on the Budyko framework significantly reduced the biases introduced by the constant α e in (semi-)humid basins. Albeit some biases remained, the water-balance ET wb for 1998-2014 in the CAMELS basins were better reproduced by the spatially varyingα e (Fig. 6b).
To confirm the improved performance of the combined CR-Budyko method across Australia, we resampled the new CR ET a estimates to 0.5 • × 0.5 • and compared them with the grid ET wb data. The ET a products by GLEAM, Flux-Com, and ERA5-Land were evaluated with the grid evaluation reference. As shown, the CR method with a constant α e = 1.15 overrated the mean annual ET a along the eastern and the northern coastlines (Fig. 7b), underperforming the physical, the machine-learning, and the land surface models (Figure 8a). Although the smaller α e = 1.10 made the CR method perform better, its predictability was still poorer than the three advanced models, and the residual variance was as large as in the case of α e = 1.15 (Fig. 8b).
In contrast, when employing theα e conditioned by local climate conditions, the same CR formulation could alleviate the overestimation along the coastlines (Fig. 7c). The Budyko-function-basedα e led the CR ET a estimates to neatly agree with the grid ET wb , and the residual variance was much smaller than in the case of α e = 1.10 (Fig. 8c). The CR method withα e clearly outperformed the three advanced models in reproducing the grid ET wb estimates ( Fig. 8d-f). Although the referenced grid ET wb has some error sources associated with the upscaling of P and Q, our comparative evaluation suggests that conditioning α e with local climate conditions could substantially reduce the uncertainty of CR ET a estimates in ungauged areas.

Constraining the CR with the Budyko framework for ungauged areas
The CR explains the dynamic equilibrium between the atmospheric ET p and the underlying moisture conditions, while the Budyko framework describes the steady-state water balance with climatic controls (i.e., P and ET w ). The analytical link between the CR and the Budyko equations, hence, implies that the atmospheric self-adjustment needs to be conditioned by the long-term climate conditions. Constraining the Turc-Mezentsev equation by the polynomial CR, Kim and Chun (2021) found that Q changes would be more sensitive to climatic changes than when they were not linked. In the opposite direction, the CR can be constrained by the Budyko equation to determine its essential parameter. In Crago and Qualls (2018), the optimal α e for the linear CR of Crago et al. (2016) varied largely between 1.00 and 1.43. This point-scale experiment has already suggested that a constant α e is unlikely suitable for definitive CRs to predict ET a in Australia. The ratio between the aerodynamic and the radiation components of ET w is evidently affected by the heat entrainment from the top of the boundary layer (Baldocchi et al., 2016), the dissimilarity between heat and water vapor sources (Assouline et al., 2016), the large-scale synoptic changes (Guo et al., 2015), and the horizontal advection of dry-air mass (Jury and Tanner, 1975). More recently, Han et al. (2021) proved the nonlinear dependence of ET w on ET eq , and Yang and Roderick (2019) showed α e changing with R n over ocean surfaces. Hence, the constant α e assumption underpinning the calibration-free CR is counterintuitive to the theoretical and empirical evidence. Although Ma et al. (2021) found some global applicability of the calibrationfree CR, its performance remains unknown in most of the Australian surfaces and in many ungauged basins over the world.
Since ET a plays a pivotal role in the terrestrial water and energy balances, the partitioning of R n into the latent and the sensible heat fluxes cannot be independent of the partitioning of P into ET a and Q. On a mean annual scale, P and ET w are the major determinants of the P partitioning, and thus the parameter α e might not be independent of P. Given the large variability of P, assuming a fixed α e across a continental area may introduce considerable biases to CR ET a estimates. Thus, discarding available P data may not be a good choice when predicting ET a by the CR method in ungauged areas. It is noteworthy that dominantly explained the spatial variation of the mean annual x of the 71 CAMELS basins, and theα e values conditioned by local climates were of a large spatial variation. This suggests that the CR with a constant α e may produce unreliable ET a estimates in ungauged locations.
Nonetheless, the low performance with a constant α e does not indicate that the CR method underperforms the sophisticated ET a models. The simple polynomial CR seemed to outperform the advanced the advanced physical, machinelearning, and land surface models, when its parameter was conditioned by local climates. The proposed CR-Budyko framework enabled us to regionalize the optimal α e for the CR method from gauged basins to ungauged locations in an empirical manner. It should be highlighted that the CR with spatially varyingα e produced much smaller residual variance than the three advanced models.

Remaining issues and caveats
In seven Australian eddy-covariance flux towers, Crago et al. (2022) found that the optimal α e for the polynomial CR was 1.35 for predicting daily ET a in the dimensionless form (i.e., y = ET a /ET p ). However, it was increased to 1.42, 1.45, 1.47, and 1.50 to simulate the dimensional latent heat fluxes at daily, weekly, monthly, and annual timescales, respectively. This implies that the timescale would largely affect the optimal α e for the definitive CRs. Though the stationary Budyko equation can become a constraint at a mean-annual scale, how to capture the scale dependence of α e is a remaining question.
Further questions can arise as to how to quantify ET p and E pmax . For example, the α e values from ET p with the Rome wind function rely upon an unrealistic assumption that the aerodynamic resistance on a vegetated surface is equivalent to that of an open-water surface. It is still unknown if this assumption is practically valid because the Penman equation with the Rome wind function may result in unrealistically high ET p , even on a large wet area (McMahon et al., 2013). Given the importance of the aerodynamic resistance in mod-  ulating surface temperature (Chen et al., 2020), ignoring its variability may become a significant error source for the CR method at both annual and subannual timescales.
In addition, there are some caveats in our case study. We employed the meteorological data different from those used in Ma et al. (2021). The ERA5-Land dataset is a downscaled version of the ERA5 data (Hersbach et al., 2020) by which Ma et al. (2021) predicted ET a globally. Ma et al. (2021) incorporated remotely sensed albedo and emissivity together with a correction factor when calculating R n , whereas we used the sum of the ERA5-Land latent and sensible heat fluxes. Those input differences may lead to differences in CR ET a estimates.
The gridded GRUN Q, too, has some uncertainty sources, though it is the ensemble of many runoff simulations from 21 different atmospheric forcing inputs. In the machine-leaning process by Ghiggi et al. (2021), some Q observations affected by human activities (e.g., dam regulation and return flows from groundwater abstraction) might not be excluded, potentially disrupting the empirical relationship between atmospheric forcing and natural flows. In addition, the uncertainty of SILO P might be non-negligible in areas with limited weather stations and in mountainous areas (Fu et al., 2022). Though we reduced the potential biases of the gridded P and Q datasets by temporal averaging, the grid-scale ET wb estimates should be treated as plausible values rather than exact observations.

Summary
Via a case study in Australia, we showed that the polynomial CR by Szilagyi et al. (2017) is unlikely to perform well when local climate conditions are neglected. The assumption of a constant Priestley-Taylor coefficient cannot reflect the longterm water balance; thereby, CR ET a estimates can be biased. We resolved this problem by conditioning the CR with the traditional Budyko equation, and it allowed for a reasonable determination of the essential parameter in ungauged locations. The following conclusions are worth emphasizing: 1. The constant Priestley-Taylor coefficient transferred from fractional wet locations could make the CR method perform poorly in closing the local water balance. The unrealistic assumption could make the CR method underperform the advanced physical, machinelearning, and land surface models.
2. The Budyko framework can play a role in determining the degree of ET p adjustment at the mean annual scale. It allows for upscaling of the Priestley-Taylor coefficients from gauged to ungauged locations.
3. The Priestley-Taylor coefficients conditioned by local climates made the CR better close the basin-scale water balance. The spatially varying Priestley-Taylor coeffi-cients seemed to make the CR method outperform the advanced ET a models. Author contributions. DK, MC, and JAC organized this study together. DK built the research framework, simulated ET a with the CR method, and drafted the manuscript. JAC processed the modeled ET a datasets and reviewed the draft, and MC actively participated in discussing the results.
Competing interests. The contact author has declared that neither of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Financial support. This research has been jointly supported by the KEITI of MOE (grant no. 2022003640001) and the NRF of MSIT (grant no. NRF-2022R1A2C2010266).
Review statement. This paper was edited by Adriaan J. (Ryan) Teuling and reviewed by two anonymous referees.