Dynamic Scaling of the Generalized Complementary Relationship Improves Long-term Tendency Estimates in Land Evaporation

Most large-scale evapotranspiration (ET) estimation methods require detailed information of land use, land cover, and/or soil type on top of various atmospheric measurements. The complementary relationship of evaporation (CR) takes advantage of the inherent dynamic feedback mechanisms found in the soil-vegetation-atmosphere interface for its estimation of ET rates without the need of such biogeophysical data. ET estimates over the conterminous United States by a new, globally calibrated, static scaling (GCR-stat) of the generalized complementary relationship (GCR) of evaporation were compared to similar estimates of an existing, calibration-free version (GCR-dyn) of the GCR that employs a temporally varying dynamic scaling. Simplified annual water balances of 327 medium and 18 large watersheds served as ground-truth ET values. With long-term monthly mean forcing, GCR-stat (also utilizing precipitation measurements) outperforms GCR-dyn as the latter cannot fully take advantage of its dynamic scaling with such data of reduced temporal variability. However, in a continuous monthly simulation, GCR-dyn is on a par with GCR-stat, and especially excels in reproducing long-term tendencies in annual catchment ET rates even though it does not require precipitation information. The same GCR-dyn estimates were also compared to similar estimates of eight other popular ET products and they generally outperform all of them. For this reason, a dynamic scaling of the GCR is recommended over a static one for modeling long-term behavior of terrestrial ET.


Introduction
Land surface evapotranspiration (ET) is a central component in the Earth's energy, water, and carbon cycles (Wang and Dickinson, 2012;Fisher et al., 2017). Accurate ET information is therefore essential for a better understanding of land−atmosphere interactions (Seneviratne et al., 2006) and the biosphere's water−carbon coupling (Biederman et al., 2016;Feng et al., 2016). It also improves drought predictions (Pendergrass et al., 2020) and helps to find answers to water resources sustainability issues (Condon et al., 2020).
While the globally distributed eddy-covariance flux towers have contributed significantly to our knowledge of ET across a wide range of ecosystems [see a recent review by Baldocchi (2020)], the spatiotemporal variation of global ET and its response to the changing climate remains highly uncertain (Mueller et al., 2011;Liu et al., 2016) because the estimation of long-term, spatially resolved ET is yet laden by difficulties in parameterizing the biophysical processes (e.g., root water uptake, stomatal resistance and its response to CO 2 concentration changes) that control ET in the current land surface models (LSMs) (Ukkola et al., 2016;Ma et al., 2017) and remote sensing algorithms (Vinukollu et al., 2011;Velpuri et al., 2013). In addition to possible model structural errors, the uncertainties in the estimated ET can also arise from errors in their gridded vegetation (Fang et al., 2019) and soil (Zheng and Yang, 2016) parameters due to the large degree of complexity/heterogeneity found in terrestrial ecosystems. For example, most LSMs within NLDAS-2 (the North American Land Data Assimilation System, phase 2) still utilize the NOAA normalized difference vegetation index data developed by Gutman and Ignatov (1998) on a five-year-mean monthly basis without any interannual variation as input (Xia et al., 2012), failing to reasonably capture the impact of vegetation changes on ET. Besides, a recent sensitivity study by Li et al. (2018) demonstrated that the Noah-MP LSM cannot always capture the effect of spatial changes in forest and/or soil types on the simulated ET because of the inherent uncertainties in multiple land cover and soil texture data.
As an alternative, the complementary relationship (CR) (Bouchet, 1963) of evaporation inherently accounts for the dynamic feedback mechanisms found in the soil−vegetation−atmosphere interface, and thus provides a viable, robust alternative for land ET estimation relying solely on standard atmospheric forcing without the need for any soil or vegetation data. The description in the next two paragraphs of the applied CR method parallels that of Ma and Szilagyi (2019).
The generalized nonlinear version of the complementary relationship (GCR) by Brutsaert (2015) relates two scaled variables, x = E w E p −1 and y = E E p −1 as Here, E (mm d −1 ) is the actual ET rate, while E p (mm d −1 ) is the potential ET rate, i.e., the ET rate of a plot-sized wet area in a drying (i.e., not completely wet) environment, typically specified by the Penman (1948)

equation as
∆ where (hPa °C −1 ) is the slope of the saturation vapor pressure curve at air temperature, T (°C), and γ is the psychrometric constant (hPa °C −1 ). R n and G are net radiation at the land surface and soil heat flux into the ground, respectively (the latter is typically negligible on a daily or longer time scale), in water equivalent of mm d −1 . The e * term denotes the saturation, while e [= e * (T d )] is the actual vapor pressure of the air [hPa; their difference is called the vapor pressure deficit (VPD)]. T d is the dewpoint temperature, and f u is a wind function, often formulated (e.g., Brutsaert, 1982) as where u 2 (m s −1 ) is the 2-m horizontal wind speed. The so-called wet-environment ET rate, E w (mm d −1 ), of a well-watered land surface having a regionally significant areal extent, is specified by the Priestley and Taylor (1972) equation: The dimensionless Priestley−Taylor (PT) coefficient, α, in Eq. (4), normally attains values in the range of [1.1−1.32] (Morton, 1983). For large-scale model applications of gridded data, Szilagyi et al. (2017) proposed a method of finding a value for α, thus avoiding the need for any calibration.
Very soon after the publication of the GCR, Crago et al. (2016) and Szilagyi et al. (2017) introduced a necessary scaling into Eq. (1) by means of a time-varying wetness index, w = (E p,max − E p )(E p,max − E w ) −1 , to define the dimensionless variable, X, as X = wx, by which Eq. (1) keeps its original nonlinear form, i.e., Note that Eq. (4) in Priestley and Taylor (1972) was designed with measurements under wet environmental conditions; therefore, Δ should be evaluated at the wet-environment air temperature, T w (°C), instead of the typical dryingenvironment air temperature, T (Szilagyi and Jozsa, 2008;Szilagyi, 2014). By making use of a mild vertical air temperature gradient (de Vries, 1959;Szilagyi and Jozsa, 2009;Szilagyi, 2014) observable in wet environments (as R n is consumed predominantly by the latent heat flux at the expense of the sensible one, and water representing an unusually high latent heat of the vaporization value found in nature), T w can be approximated by the wet surface temperature, T ws (°C). Note that T ws may still be larger than the drying-environment air temperature, T, when the air is close to saturation, but the same is not true for T w , due to the cooling effect of evaporation. In such cases, the estimated value of T w should be capped by the actual air temperature, T (Szilagyi, 2014;Szilagyi and Jozsa, 2018). Szilagyi and Schepers (2014) demonstrated that T ws is independent of the size of the wet area. Thus, T ws can be obtained through iterations from the Bowen ratio (β) of the sensible and latent heat fluxes (Bowen, 1926) when applied over a small, plot-sized, wet patch (by the necessary assumption that the available energy for the wet patch is close to that for the drying surface) the Penman equation is valid for, i.e., .
(6) E p,max in the definition of X within Eq. (5) is the maximum value that E p can achieve (under unchanging available energy for the surface) during a complete dry-out (i.e., when e becomes close to zero) of the environment, i.e., in which T dry (°C) is the so-achieved dry-environment air temperature. The latter can be estimated from the (isoenthalp) adiabat of an air layer in contact with the drying surface (Szilagyi, 2018a), i.e., where T wb (°C) is the wet-bulb temperature. T wb can be obtained with the help of another iteration of writing out the Bowen ratio for adiabatic changes (e.g., Szilagyi, 2014), such as For a graphical illustration of the saturation vapor pressure curve, the different temperatures and the related ET rates defined, please refer to Ma and Szilagyi (2019). The same source also includes a brief description of how the CR evolved into Eq. (5) over the past 40 years. Additionally, it plots selected historical CR functions over sample data, and explains how assigning a value of α is performed without resorting to any calibration. A sensitivity analysis of the ET rates in Eq. (5) to their atmospheric forcing is found in Ma et al. (2019).
While Brutsaert et al. (2020)  (2019) performed one (and the same one) via a dynamic wetness index. Whereas the wetness index assigns increasing values to wetter environmental conditions, the aridity index does the same to drier ones. Brutsaert et al. (2020) did not include this dynamic wetness index method in their study, and therefore the present work was initiated to fill this gap.

Model applications
The time-varying (and thus dynamic) scaling of x (Crago et al., 2016;Szilagyi et al., 2017) (5), is necessary because the GCR of Brutsaert (2015) unrealistically predicts near-zero land evaporation only when E w in x itself approaches zero. This is because the potential evaporation rate, E p , in the denominator of x always assumes well-bounded values due to physical limits on the range of its constituents, i.e., net radiation, soil heat flux, air temperature, wind speed, and VPD.
An alternative, static scaling of x by Brutsaert et al. (2020) takes place via an adjustable parameter, α c , that acts as the PT-α value for wet environments. Since Eq. (4) can also be written as E w = αE e , where E e is the equilibrium evaporation rate of Slatyer and Mcllroy (1961), thus the scaled variable, X, becomes X = α c E e E p −1 = α c xα −1 . The spatially variable (but constant through time at a given location) value of α c was then related to a long-term aridity index by Brutsaert et al. (2020), with the latter defined as the ratio of the mean annual E p and rain depth, and globally calibrated with the help of additional water-balance data, requiring altogether seven parameters in highly nonlinear equations.
Note that the X = wx scaling by Crago et al. (2016) and Szilagyi et al. (2017) requires only the forcing variables (R n , G, T, u 2 , and VPD), without the need for external precipitation/rain data, which is significant as precipitation is possibly the most uncertain meteorological variable to predict in climate models. It is important to mention that w changes with each value of x, unlike α c . As Szilagyi et al. (2017) demonstrated, a (temporally and spatially) constant value of the PT α, necessary for x, can be set by the sole use of the forcing variables, without resorting to additional water-balance data of precipitation and stream discharge, thus making the approach calibration-free on a large scale (Szilagyi, 2018b;Ma et al., 2019;Ma and Szilagyi, 2019) where wetenvironmental conditions, necessary for setting the value of α, can likely be found. Note that setting a constant value of α is also necessary for Brutsaert et al. (2020) in order to force their spatially variable but temporally constant α c values to reach a predetermined value of about 1.3 under wet conditions. Despite almost half a century of research following publication of the Priestley and Taylor (1972) equation, there is still no consensus about what environmental variables (atmospheric, radiative, and/or surficial properties), and exactly how their spatial and temporal averaging, influence the value of the PT α. Until compelling information is available on these variables, a spatially and temporally constant α value may suffice for modeling purposes.
As was found by Szilagyi (2018b), the value of the PT α depends slightly on the temporal averaging of the forcing data, i.e., whether or not the monthly values are long-term averages [yielding α = 1.13 (Szilagyi et al., 2017) and 1.15 (Szilagyi, 2018b), respectively]. Therefore, here, it is tested if such is the case for the globally calibrated model of Brutsaert et al. (2020). Namely, if its performance is affected by similar changes (from long-term mean monthly values to monthly values) in the input/forcing variables, then some caution must be exercised during its routine future application, and recalibration of its seven parameters may be necessary. Note that besides the different scaling of x, everything is the same (including data requirements) in the two GCR model versions applied here, except that Δ in E e is evaluated at the measured air temperature in Brutsaert et al. (2020) while the same in E w (= αE e ) is evaluated at an estimated wet-environment air temperature (Szilagyi et al., 2017) explained above.
Both model versions (denoted for brevity by GCR-stat and GCR-dyn, respectively) were tested over the conterminous United States, first with long-term averages (1981−2010) of monthly 32-km resolution North American Regional Reanalysis (NARR) (Mesinger et al., 2006) radiation and 10-m wind (u 10 ) data [reduced to 2-m values via u 2 = u 10 (2/10) 1/7 (Brutsaert, 1982)], as well as with 4.2-km PRISM air, and dewpoint temperature values (Daly et al., 1994) followed by a continuous 37-year simulation of monthly values over the 1979−2015 period. The NARR data were resampled to the PRISM grid by the nearest neighbor method. Monthly soil heat fluxes were considered negligible. Evaluation of the model estimates were performed by water-balance estimates of basin-representative evaporation rates (E wb ) with the help of United States Geological Survey two-and six-digit Hydrologic Unit Code (HUC2 and HUC6) basin ( Fig. 1) discharge data (Q) together with basin-averaged PRISM precipitation (P) values as E wb = P − Q, either on an annual (for trend analysis) or long-term mean annual basis. The application of a simplified water balance is justifiable as soil-moisture and groundwater-storage changes are typically negligible over an annual (or longer) scale (Senay et al., 2011) for catchments with no significant trend in the groundwater-table elevation values.

Results and discussion
With the long-term mean monthly data, GCR-stat performed slightly but consistently better than GCR-dyn (Fig. 2), reflected best in the Nash−Sutcliffe model efficiency (NSE) and root-mean-square error (RMSE) values, both models providing unbiased, basin-averaged mean annual ET estim-ates. This outcome is unsurprising, as GCR-stat takes advantage of measured precipitation while GCR-dyn does not.
However, the picture changes when switching from long-term mean monthly forcing values to monthly values in a continuous simulation (Fig. 3). GCR-dyn, with a slightly changed PT-α value [from 1.13 to 1.15, using the procedure of Szilagyi et al. (2017)] continues to produce unbiased estimates of basin-averaged mean annual evaporation values. However, the globally calibrated GCR-stat model underestimates the water-balance-derived values by about 10% [i.e., relative bias (RB) of −0.09 for both basin scales] and produces reduced interannual variability (see the horizontally elongated "crosses " for the HUC2 basins in Fig. 3) in comparison with GCR-dyn. Reduced model performance of GCR-stat is also apparent in the long-term linear tendencies (obtained as least-squares fitted linear trends) of the basin-averaged annual evaporation values (Fig. 4) by being less effective than GCR-dyn in reproducing the observed linear trends in the water-balance data.
As on a mean-annual basis GCR-stat performs better than (with mean monthly values) or about equal to (in a continuous simulation) GCR-dyn by exploiting precipitation data (which on the watershed scale naturally serves as an upper bound for land ET), its weakened performance in trends can only be explained by the same reliance on the long-term means of the precipitation (and E p ) rates in the (therefore) static α c values that will hinder its response to slow (decadal) changes in aridity. The same problem cannot occur in GCR-dyn, since its wetness index (w) is updated in each step of calculations.
The current GCR-dyn model has already been shown to (a) yield correlation coefficient values in excess of 0.9 with local measurements of latent heat fluxes across diverse climates in China (Ma et al., 2019) in spite of large differences in spatial representativeness (i.e., grid resolution vs flux measurement footprint), and (b) outperform several popu- lar large-scale ET products over the conterminous United States (Ma and Szilagyi, 2019). These products include three LSMs-namely, Noah (Chen and Dudhia, 2001), VIC (Liang et al., 1994), and Mosaic (Koster and Suarez, 1996); two reanalysis products-namely, NCEP-II (Kanamitsu et al., 2002), and ERA-Interim (Dee et al., 2011); another two remote-sensing based approaches-namely, GLEAM (Miralles et al., 2011;Martens et al., 2017) and PML (Zhang et al., 2017;Leuning et al., 2008); and a spatially upscaled eddy-covariance measurement product, FLUXNET-MTE (Jung et al., 2011). In a comparison with water-balance data, GCR-dyn turns out to produce even better statistics on the HUC2 scale than the spatially interpolated eddy-covari-ance measured ET fluxes (Fig. 5), which is remarkable from a calibration-free approach. GCR-dyn especially excels in the long-term linear tendency estimates of the HUC2 ET rates, demonstrated by Figs. 6 and 7. As FLUXNET-MTE employs several temporally static variables for its spatial interpolation method, its ability to capture long-term trends is somewhat limited (Jung et al., 2011). On the contrary, the dynamic scaling inherent in GCR-dyn automatically adapts to such trends and identifies them rather accurately.
Among the different popular large-scale ET models, GCR-dyn produces multi-year mean annual ET rates closest in its spatial distribution to those of FLUXNET-MTE (Fig. 8), with a spatially averaged ET value almost identical (both in its spatial average and standard deviation) to that of GLEAM (Fig. 8), which is a remote-sensing based approach. Note that all models of the comparison (except GCR-dyn) rely on precipitation data as input, which greatly aids ET estimations, as on a regional scale and long-term basis precipitation forms an upper bound for terrestrial ET rates; plus, it may drive any required soil-moisture water-balance calculations.
The spatial distribution of the modeled multi-year linear ET trends is displayed in Fig. 9. Again, GCR-dyn produces results closest in spatial distribution to FLUXNET-MTE in terms of the statistically significant trends and to GLEAM in general. For a more detailed discussion of model comparisons (including additional model descriptions), please refer to Ma and Szilagyi (2019).
In conclusion, it can be stated that the GCR of evaporation (Brutsaert, 2015) is a very effective tool in land ET modeling, as it requires only a few, largely meteorological forcing variables, and avoids the need for detailed soil-moisture and land-cover information. Although attractive, as its (altogether seven) parameters have already been globally precalibrated, the GCR model version (GCR-stat) of Brutsaert et al. (2020) may, however, not perform optimally in estimating long-term tendencies in basin-wide ET rates. This is par- ticularly the case in comparison to an earlier, calibrationfree version (GCR-dyn), having no precalibrated parameter values but requiring that its sole, temporally−and spatially−constant parameter (i.e., the PT coefficient) be set using the actual forcing data through a largely automated method, described in Szilagyi et al. (2017). Since in a continuous monthly simulation both models performed about the same, while the GCR-dyn produced better long-term tendencies, a dynamic scaling of E w E p −1 is recommended over a static one in future applications of the GCR of evaporation.
As has been recommended before (Szilagyi, 2018b;Szilagyi and Jozsa, 2018;Ma and Szilagyi, 2019), it is reiterated here that GCR-dyn, due to its minimal data requirement, calibration-free nature and dynamic scaling, may continue to serve as a diagnostic and benchmarking tool for more complex and data-intensive models of terrestrial ET rates, including calibration/verification (for past values) and reality checking (for future scenario values) of the LSM-pre-     dicted ET rates of existing regional and general circulation models.
Open Access this article is distributed under the terms of the creative commons attribution 4.0 international license (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the creative commons license, and indicate if changes were made.