**1. Introduction**

Terrestrial evapotranspiration (ET) is an important part of the water cycle and surface energy balance, which is not only the nexus of the water cycle in the soil–vegetation–atmosphere

**Citation:** Shang, C.; Wu, T.; Ma, N.; Wang, J.; Li, X.; Zhu, X.; Wang, T.; Hu, G.; Li, R.; Yang, S.; et al. Assessment of Different Complementary-Relationship-Based Models for Estimating Actual Terrestrial Evapotranspiration in the Frozen Ground Regions of the Qinghai-Tibet Plateau. *Remote Sens.* **2022**, *14*, 2047. https://doi.org/10.3390/rs14092047

Academic Editor: Sayed M. Bateni

Received: 10 March 2022 Accepted: 21 April 2022 Published: 25 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

continuum, but also a crucial link in the mass exchange and energy balance of the land surface [1,2]. Accurate estimation of ET is important to deeply understand how much liquid or solid water from land transforms as vapor in the atmosphere, which has significant implications for regional weather or climate conditions, and water resource management, such as ecological conservation and agricultural management, etc. [3–5].

In practice, ET estimation is often challenging because of complex land–atmosphere interactions, which often involve relevant micrometeorological, hydrological, or ecological methods [6,7]. In recent decades, grea<sup>t</sup> progress has been made on the theories and technologies of actual evapotranspiration (ETa) regardless of observations or simulations, significantly contributing to our knowledge of the ET process [8–12]. Current ET estimation approaches, such as land surface models or remote sensing technologies, have been successfully applied to various climatic or ecosystem zones [13–15]; however, the above approaches generally have relatively higher uncertainty in parameterization schemes at specific zones and require detailed information on vegetation or soil properties [16], which is especially difficult for data-scarce regions. For example, the Penman-Monteith (PM) [17,18] equation is often used to calculate potential evapotranspiration (ETp), which is close to ETa under ample water supply conditions, and with a higher accuracy. However, under water-limited conditions, the PM method assumes typically that ETa is proportional to the ETp (rescaling with soil moisture content), giving rise to practical difficulties such as soil moisture data unavailability in harsh environments where the estimation has a higher uncertainty. In addition, the assumption of the positive relationship between ETa and ETp may be rejected by observational evidence that supports negative correlations between the two variables [19]. Han [20] emphasized that the correlation between ETa and ETp depends mainly on water availability rather than always being positive.

Bouchet [21] first proposed another important hypothesis about ETa estimation, often called the "Bouchet hypothesis". In contrast to the PM method, Bouchet [21] described the negative relationship between ETa and ETp caused by land-atmosphere feedback, which is also referred to as the complementary relationship (CR) approach. The greatest advantage of the CR approach is that it does not require any vegetation and soil information, and only a few model parameters and routine meteorological forcing data are needed [22,23]. Various models based on the CR approach have been proposed over the past decades [24]. The most widely used CR-based models include the AA (advection-aridity) model [25], CRAE (complementary relationship areal evapotranspiration) model [26], and GG (Granger and Gray) model [27], which are all based on the linearly symmetric CR approach. However, some later studies [28,29] have found that the complementary relationships are not symmetric, and they further proposed various asymmetric linear CR-based functions. Nevertheless, Han [30] found that both symmetric and asymmetric linear CR-based functions have poor generality and are applicable only to some mild climate regions. To extend the applicability of the CR function, Han [31] and Brutsaert [32] generalized the CR approach; subsequently, many studies have further refined the CR approach based on the above [33–37]. At present, various CR-based models have been applied to different ecosystems (shrubland, cropland, grassland, etc.) or climate (arid, semiarid, temperate) regions around the world on annual, monthly, daily, and sub-daily timescales [38–42], but very few studies focused on frozen ground regions on the Qinghai-Tibet Plateau (QTP), which may limit the development and application of the CR approach in frozen ground regions.

The QTP accounts for approximately a quarter of China's land area, with an average altitude greater than 4000 m above sea level, known as the "roof of the world" and the "third Pole" [43]. Its vast topography and high altitude have significant impacts on the weather and climate of East Asia and even the world through thermal and dynamic forcing [44,45]. Additionally, frozen ground is widely distributed on the QTP [46], accounting for 96% of the area (including approximately 40% permafrost and 56% seasonally frozen ground). In these regions, freeze-thaw seasonal cycles and permafrost degradation [47,48] have certainly affected the regional hydrothermal balance of the soil–vegetation–atmosphere system; however, one of the current obstacles is that little is known about how seasonal

and long-term ETa rates have changed in the frozen ground regions. In recent decades, several comprehensive observation sites have been established [49,50], improving our understanding of the ET process, but large-scale ETa estimation in the frozen ground regions of the QTP is still sparse. Hence, the above CR approach using only routine meteorological observations seems to be a feasible method. Wang [51] have improved and validated one CR-based model when applied to a permafrost sites on the QTP, but it is still unclear about which CR-based models are best for ETa estimation in this region.

This study utilized four comprehensive field sites, which are all located in the frozen ground regions of the QTP, based on five CR-based models with in situ measurement meteorological and eddy-covariance (EC) flux data, to calculate variables and optimize parameters in CR-based models. The objective of this study was to clarify the applicability of five widely used CR-based models in the frozen ground regions of the QTP and further investigate the uncertainty of the CR approach in ETa estimation. The results of this study could provide a reference for the applicability of CR-based models in the frozen ground regions of the QTP and provide a simple and feasible option for future large-scale terrestrial ETa in such data-scarce regions.

### **2. Materials and Methods**

### *2.1. Site Description*

In this study, in situ measurement data were collected from four comprehensive observation field sites (Figure 1). Tanggula (TGL, 91.86◦E, 32.58◦N, 5100 m asl) and Xidatan (XDT, 94.13◦E, 35.72◦N, 4538 m asl) lie in the south edge and northern limit of permafrost regions of QTP, covered by alpine steppe and alpine meadow, where the annual average air temperature is about −4.7 ◦C and −3.6 ◦C, and the annual accumulated precipitation is about 352 and 384.5 mm, respectively [52]. Another two field sites, Nagqu (BJ, 91.90◦E, 31.37◦N, 4509 m asl) and Nam Co (NAMORS, 90.96◦E, 30.77◦N, 4730 m asl), lie in the seasonally-frozen ground regions of QTP, covered by alpine meadow and alpine steppe, where the annual average air temperature is about 0.54 ◦C and −0.36 ◦C, and the annual accumulated precipitation is about 436 and 462 mm, respectively [53]. All the sites are near the Qinghai-Tibet highway.

**Figure 1.** Locations of four observation field sites on the QTP.

### *2.2. In Situ Measurement Data and Data Processing*

All four selected field sites include eddy covariance flux data, meteorological data, soil temperature, soil moisture content, and soil heat flux data in a shallow soil layer. Eddy

covariance flux data, including latent heat flux (LE) and sensible heat flux (H), were used to identify the performance of the different CR-based models. All EC systems were installed on towers about 3 m distance from the ground, consisting of a three-dimensional sonic anemometer (CAST3, Campbell Scientific Inc., USA) that measured instantaneous horizontal (u,v), vertical (w) wind speeds and sonic air temperature fluctuations, and an open path infrared gas analyser (Li-7500, LI-cor Inc., USA) that measured the water vapor density and carbon dioxide concentrations fluctuations. The EC instruments were all sampled at a frequency of 10 Hz. The corresponding meteorological data, including roughly 2 m downwards/upwards shortwave and longwave radiation, air temperature, relative humidity, wind speed, and pressure, were employed to drive each CR-based model. Soil temperature, soil moisture content, and soil heat flux data in shallow soil layer (within 10 cm from ground surface) were used to calculate ground heat flux for variables in CR-based models. All above observations were recorded at intervals of 30 min. Data of TGL and XDT were provided from Cryosphere Research Station on the Qinghai-Tibet Plateau, Northwest Institute of Eco-Environment and Resources, CAS (http://new.crs.ac.cn/, 1 September 2021). Data of BJ and NAMORS were provided from National Tibetan Plateau Third Pole Environment Data Center, CAS (https://data.tpdc.ac.cn/en/data/b9ab35b2-81fb-4330-925f-4d9860ac47c3/, 1 September 2021). In addition, two different selected study periods of each field site were determined by data integrity and continuity; observation of first and second year was taken as calibration and validation period, respectively.

The EC flux raw data processing was conducted on EddyPro7.0.6 software, and the main processing steps included spike detection, lag corrections of H2O/CO2 relative to the vertical wind component, sonic virtual temperature correction, coordinate rotation (planar fit rotation), corrections for Webb-Pearman-Leuning density fluctuations, and frequency response correction. The 30 min data output by EddyPro software was screened as follows: (1) the data of instrument error was removed; (2) 30-min record data with more than 10% missing values were excluded; (3) night data with weak turbulence (friction velocity (u\*) less than 0.1 m s<sup>−</sup>1) were eliminated; (4) data with quality control code "2" ("2" indicate poor quality) were filtered out [54]. Then, half-hourly data during data available period accounts for filling the gaps; for gaps within seven consecutive days, they were filled using REddyProcWeb online tool (https://www.bgc-jena.mpg.de/bgi/index.php/ Services/REddyProcWeb, 1 September 2021); for gaps longer than seven consecutive days, they were filled alternatively using aerodynamic method [55] or Bowen ratio energy balance method [56] for 30-min intervals. Finally, 30 min gap-filled flux data were aggregated into daily values. Ground heat flux were computed by linear method and thermal diffusion equation (TDE); detailed calculation processes are described in Yang [57] and Yao [58]. The daily energy balance closure ratios are 0.85, 0.79, 0.95, and 1.1 at TGL, XDT, BJ, and NAMORS, respectively (Figure 2).

It should be noted that on snow cover days, snow sublimation was one of the main evapotranspiration forms for four selected sites. Due to lack of observation of snow cover, this study employed the albedo instead, defined as the ratio of upward short wave radiation to downward shortwave radiation, to identify whether the observation field was snow-free or not. To exclude the influence of solar elevation angle, daily mean albedo was calculated as the average of half-hourly albedo from 10:00 to 14:00, local time. According to the variations of albedo at each observation fields, snow cover days were determined. In addition, snowfall still occurs during warm season at some observation fields, however, snow generally melts quickly to liquid water in few hours; thus, on a daily basis, only sublimation during cold season (from October to April in the next year) were considered. Then, the latent heat of vaporization— *λv* (J Kg−1), for snow free days—and the latent heat of sublimation— *λs* (J Kg−1), for snow cover days—were calculated using Equation (1), and the results were then used to calculate the actual evapotranspiration or sublimation [59]:

$$\begin{cases} \lambda\_{\upsilon} = (2500 - 2.4T\_d) \times 10^3 \\ \lambda\_s = (2834.1 - 0.149T\_d) \times 10^3 \end{cases} \tag{1}$$

where *Ta* (◦C) is the air temperature measured at 2 m high.

**Figure 2.** Comparison of available energy (*Rn* − *G*) and sum of turbulent fluxes (H + LE) for the daily-averaged measurements after gap-filled during study period at four selected sites: (**a**) TGL; (**b**) XDT; (**c**) BJ; (**d**) NAMORS. The black solid line is the 1:1 line, and the black dash line fitted by the linear regression.

### *2.3. Complementary Relationship (CR) Approach*

The CR approach involves three types of evaporation, namely, the actual evaporation (*E*, note: followed by "*E*" instead), the apparent potential evaporation (*Epa*), and the wetenvironment evaporation (*Epo*). The CR approach is expressed as follows: for relatively larger and homogeneous surfaces with minimum advection, under certain net radiation inputs, *E* decreases with the availability of limited water over the underlying surface, and the energy that would be consumed by the latent heat flux thus becomes the sensible heat, thus increasing *Epa*. Hence, one can predict water-limited *E* by gauging how much *Epa* is raised from the hypothetical evaporation rate that should occur under the full wetness (*Epo*). Previous studies [60] hold that a unit decrease in *E* yields a corresponding unit increase in *Epa*, signifying a symmetric CR. This can be expressed as follows:

$$E\_{pa} - E\_{po} = E\_{po} - E \tag{2}$$

*Epa* and *Epo* is calculated as follows:

$$E\_{pa} = E\_{rad} + E\_{aero} = \frac{\Delta(R\_{\text{fl}} - G)}{\Delta + \gamma} + \frac{\gamma f(\text{l}I)(\varepsilon\_{\text{l}} - \varepsilon\_{\text{d}})}{\Delta + \gamma} \tag{3}$$

$$E\_{\rm pot} = \alpha\_{\rm t} E\_{\rm rad} = \alpha\_{\rm t} \frac{\Delta (R\_n - G)}{\Delta + \gamma} \tag{4}$$

where *Erad* is the radiation term, *Eaero* is the aerodynamic term, *αe* is an analog of the dimensionless Priestley-Taylor coefficient, *(Rn* − *G)* is the available energy (mm <sup>d</sup>−1), Δ is the slope of saturation vapor pressure curve at the air temperature *Ta* (kPa ◦C−1), *γ* is the psychrometric constant (kPa ◦C−1), and *eo* and *ea* are the saturation and actual vapor pressure of the air (kPa), respectively. *f*(*U*) (mm d−<sup>1</sup> kPa−1) is the wind function calculated by Penman's original empirical linear equation:

$$f(\mathcal{U}) = 2.6(1 + 0.54\mathcal{U}\_2) \tag{5}$$

where *U*2 is the measured wind speed (m s<sup>−</sup>1) at 2 m height.

Here, five CR-based models were selected to estimate E: the modified AA model [28] represents the linear CR version; the polynomial CR function [32] and sigmoid CR function [36] both represents the generalized nonlinear CR concept, and they represent different improved versions of the original CR approach mainly in physical boundary constraints; in addition, the calibration-free CR function and rescaled CR function are two improved versions based on the polynomial CR function [33–35].

### 2.3.1. Modified AA Model

The modified AA model extended the symmetric CR principle (herein referred to as K2006) and the equation becomes asymmetry as follows:

$$y\_K = \frac{b+1}{b}x\_K - \frac{1}{b} \tag{6}$$

where *b* is a coefficient that depicts the proportion of the sensible heat that increases *Epa*, *xK* = *Epo/Epa*, *yK* = *E/Epa*. The only physical constraint is *E* ≤ *Epo* ≤ *Epa*; when the water availability of the landscape is not limited, *E* is assumed to proceed at *Epa* and *E=Epo = Epa*.

### 2.3.2. Polynomial Generalized Complementary Function

Brutsaert [32] reformulated a new general polynomial complementary function satisfying boundary conditions based on strictly physical constraints, referred to as the B2015 model, and its boundary conditions were as follows:

$$\begin{cases} \begin{aligned} y\_B &= 1, \mathbf{x}\_B = 1 \\ y\_B &= 0, \mathbf{x}\_B = 0 \\ \frac{d y\_B}{d \mathbf{x}\_B} &= 1, \mathbf{x}\_B = 1 \\ \frac{d y\_B}{d \mathbf{x}\_B} &= 0, \mathbf{x}\_B = 0 \end{aligned} \end{cases} \tag{7}$$

where *xB = Epo/Epa*, *yB = E/Epa*. The function form is as follows:

$$y\_B = (2 - c)\mathbf{x}\_B^2 - (1 - 2c)\mathbf{x}\_B^3 - c\mathbf{x}\_B^4 \tag{8}$$

*c* is an adjustable parameter which need to be locally calibrated.

### 2.3.3. Calibration-Free CR Function

Szilagyi [35] utilized the same function forms as the B2015 model (herein referred to as S2017), and the equation is as follows:

$$y\_S = 2X\_S^{\;\;2} - X\_S^{\;\;3} \tag{9}$$

where *yS* is the same as *yB*, *XS* is defined as

$$X\_S = \frac{E\_{p\text{max}} - E\_{\text{pr}}}{E\_{p\text{max}} - E\_{p\text{o\\_sc}}} \frac{E\_{p\text{o\\_sc}}}{E\_{\text{pr}}} \tag{10}$$

Here two variables (*Ep* max and *Epo\_c*) are introduced to improve the performance of S2017 function. *Ep* max is the maximum value that *Epa* can theoretically reach, which may appear when the air loses all moisture, that is,

$$E\_{p\text{max}} = \frac{\Delta\_{Tdry}(R\_n - G)}{\Delta\_{Tdry} + \gamma} + \frac{\gamma f(\text{\textdegree I}) \left(\varepsilon\_{o, Tdry} - 0\right)}{\Delta\_{Tdry} + \gamma} \tag{11}$$

where <sup>Δ</sup>*Tdry* and *eo,Tdry* are the slope of the saturation vapor pressure curve and saturated vapor pressure, respectively, at extreme dry environment air temperature, *Tdry* (◦C). *Tdry* can be estimated from the adiabatic line as follows:

$$T\_{dry} = \frac{\varepsilon\_{o,Twb}(T\_d - T\_{wb})}{\varepsilon\_{o,Twb} - \varepsilon\_d} + T\_{wb} \tag{12}$$

where *Twb* (◦C) is the wet-bulb temperature, and *eo,Twb* (kPa) is the saturated vapor pressure at *Twb*. *Twb* is derived by another iteration of solving for the Bowen ratio during adiabatic changes, that is,

$$
\gamma \frac{T\_{wb} - T\_a}{e\_{o, Twb} - e\_a} = -1\tag{13}
$$

Another one variable *Epo\_c* is defined as

$$E\_{\rm po\\_c} = \alpha\_{\rm \epsilon} E\_{\rm rad} = \alpha\_{\rm \epsilon} \frac{\Delta\_{\rm weak} (R\_n - G)}{\Delta\_{\rm new} + \gamma} \tag{14}$$

where Δ*wea* (kPa ◦C−1) is the slope of the saturation vapor pressure curve at *Twes*. *Epo\_c* is the wet environment evapotranspiration rate calculated with Δ*wea* instead of Δ. Note that S2017 function uses *Epo\_c* instead of *Epo* because the latter variable is based on local air temperature, which is physically improperly employed to calculate wet environment evaporation. Szilagyi and Jozsa (2008) [61] suggested Equation (4) should use wet environment air temperature (*Twea*). *Twea* cannot be measured directly under water-limited conditions but can be approximated by the wet surface temperature (*Twes*). According to Szilagyi [62], *Twes* can be solved iteratively by implementing the Bowen ratio of a small wet patch, that is,

$$\beta\_w = \frac{R\_n - G - E\_{\rm pr}}{E\_{\rm pr}} \approx \gamma \frac{T\_{\rm vvs} - T\_a}{\varepsilon\_{o, \rm Tuves} - \varepsilon\_a} \tag{15}$$

in which *βw* is the Bowen ratio of the wet patch (assuming that the available energy for the wet patch is close to that of the drying surface), and *eo,Twes* is the saturated vapor pressure at *Twes*(≈*Twea*). Note that *Twes* may be larger than *Ta* when air is close to saturation, and in such cases, *Twea* should be replaced by *Ta*.

### 2.3.4. Rescaled Complementary Function

Crago [34] used a similar method to the S2017 model, with *Ep* max introduced to rescale *xB*, and new analytical forms (herein referred to as C2018) have been proposed as follows:

$$\begin{array}{l} \chi\_{\text{min}} = \frac{E\_{ps\\_c}}{E\_{p\text{max}}}\\ X\_{\text{C}} = \frac{x - x\_{\text{min}}}{1 - x\_{\text{min}}}\\ y\_{\text{C}} = X\_{\text{C}} \end{array} \tag{16}$$

where *yC* are consistent with the corresponding variables of *yS* in the S2017 model.

⎧⎪⎨⎪⎩

### 2.3.5. Sigmoid Generalized Complementary Function

Han [36] brought in the minimum and maximum limits of *Erad/Epen* (*<sup>x</sup>*min and *x*max) under an assumed constant *Erad* and rederived four new boundary conditions fitting to two widely accepted assumptions following Penman's combination theory, namely, *dE/dEpen = 0* in extremely arid environments and *E=Epen* in completely wet environments. The boundary conditions are as follows:

$$\begin{cases} \begin{aligned} y\_H &= 0, \mathbf{x}\_H \to \mathbf{x}\_{\min} \\ y\_H &= 1, \mathbf{x}\_H \to \mathbf{x}\_{\max} \\ \frac{d y\_H}{d \mathbf{x}\_H} &= 0, \mathbf{x}\_H \to \mathbf{x}\_{\min} \\ \frac{d y\_H}{d \mathbf{x}\_H} &= 0, \mathbf{x}\_H \to \mathbf{x}\_{\max} \end{aligned} \end{cases} \tag{17}$$

The following new sigmoid function in accordance with the above boundary conditions were also proposed (herein referred to as H2018):

$$y\_H = \frac{1}{1 + m\left(\frac{x\_{\text{max}} - x\_H}{x\_H - x\_{\text{min}}}\right)^H} \tag{18}$$

where *xH = Erad/Epen*, *yH = E/Epen*, and the calculation of *Erad* and *Epen* (equal to *Epa*) refers to Equation (3); *m* and *n* are two unknown parameters that need to be locally calibrated.

$$\begin{cases} \text{ } m = \frac{4a\left(1 + b^{-1}\right)\left(x\_{0,\text{S}} - x\_{\text{min}}\right)\left(x\_{\text{max}} - x\_{0,\text{S}}\right)}{x\_{\text{max}} - x\_{\text{min}}}\\ \quad m = \left(\frac{x\_{0,\text{S}} - x\_{\text{min}}}{x\_{\text{max}} - x\_{0,\text{S}}}\right)^{n} \end{cases} \tag{19}$$

where *x*0.5 is a variable that corresponds to *yH* = 0.5 (simply reflecting a temperate environment). Han [36] conducted that *x*0.5 = 0.5+*b*−<sup>1</sup> *<sup>α</sup>e*(<sup>1</sup>+*b*−<sup>1</sup>), *x*max and *x*min can be simply set as 1 and 0 on a daily time scale, respectively. For longer time scales, *x*max and *x*min needed to be locally calibrated.

### *2.4. Model Parameter Calibration Methods*

For the CR approach, one of the most important steps is to determine unknown parameter values. Generally, "optimal" parameter values could be either derived from transplanting from elsewhere, or are locally calibrated. First, "optimal" parameter values between two different sites are roughly equal under similar climatic or underlying surface conditions. In our study, the vegetation type of all four selected sites was grassland, so we set the "optimal" parameter values from other grassland sites reported in previous literature as the default values in each CR-based model. Note that the unknown parameters in the K2006, B2015, and H2018 models used corresponding values from the US-Fpe flux site [36] as default values, and the unknown parameters in S2017 and C2018 model used corresponding values from the Riggs Creek site at Australia [34] as default values, respectively. Second, unknown parameters of each CR-based model were all optimized by minimizing the root mean square error of the simulated and measured daily *E* during the calibrated period, and then the calibrated parameter values were used for the validation period.

### *2.5. Model Evaluation Criteria*

To assess the model performance, five statistical metrics, including root mean square error (RMSE), mean absolute error (MAE), mean bias error (MBE), Nash-Sutcliffe efficiency (NSE), and correlation coefficient (CC) are used to evaluate the accuracy of simulated *E* with in situ measurements:

$$\text{RMSE} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} (E\_{sim,i} - E\_{obs,i})^2} \tag{20}$$

$$\text{MAE} = \frac{1}{N} \sum\_{i=1}^{N} \left| E\_{\text{sim},i} - E\_{\text{obs},i} \right| \tag{21}$$

$$\text{MBE} = \frac{1}{N} \sum\_{i=1}^{N} (E\_{\text{sim},i} - E\_{\text{obs},i}) \tag{22}$$

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{N} (E\_{sim,i} - E\_{obs,i})^2}{\sum\_{i=1}^{N} \left(E\_{obs,i} - \overline{E\_{obs}}\right)^2} \tag{23}$$

$$\text{CC} = \frac{\sum\_{i=1}^{N} \left( E\_{\text{sim},i} - \overline{E\_{\text{sim}}} \right) \left( E\_{\text{obs},i} - \overline{E\_{\text{obs}}} \right)}{\sqrt{\sum\_{i=1}^{N} \left( E\_{\text{sim},i} - \overline{E\_{\text{sim}}} \right)^{2}} \sqrt{\sum\_{i=1}^{N} \left( E\_{\text{obs},i} - \overline{E\_{\text{obs}}} \right)^{2}}} \tag{24}$$

where *N* is the number of observation, *Esim,i* (*i* = 1, 2, 3, ... , *N*) is the simulated actual evapotranspiration value, *Eobs,i* is the observed actual evapotranspiration value, *Esim* and *Eobs* are the mean values of *Esim,i* and *Eobs,i*, respectively.

### **3. Results**

### *3.1. Variations in ET Rates at Four Observation Sites*

The daily variations in actual evapotranspiration (*E*), sublimation, radiation term (*Erad*), and apparent potential evapotranspiration (*Epa*) rate at four observation field sites are shown in Figure 3. Note that relatively larger errors of observation may exist at NAMORS during some periods, especially in the warm season when *E* is larger than *Epa*, which is discrepant with Penman's evaporation theory. Visibly, among four observation field sites–whether located in permafrost regions or seasonally frozen ground regions–are both clear seasonal variations of each ET variable: higher *E* and *Epa* usually occurs during the warm season, maximum *E* could reach approximately 4 mm d−<sup>1</sup> at TGL and XDT, and 5 mm d−<sup>1</sup> or more at BJ and NAMORS. For *Epa*, which exhibited greater daily fluctuations compared with *E* and *Erad*, the maximum *Epa* ranges from 5 mm d−<sup>1</sup> to 7 mm d−<sup>1</sup> at four selected sites: low *E* occurred during cold seasons and was generally close to zero when soils were frozen, except at NAMORS, and *Epa* and *Erad* was relatively lower, which was mainly controlled by net radiation. Additionally, some obvious differences between permafrost regions and seasonally frozen ground regions were that greater variations in *E* occurred at seasonally frozen ground regions than permafrost regions, and the onset time of increasing *E* after the frozen period in seasonally frozen ground regions (usually in late March) was earlier than that of permafrost regions (usually in late April), which was due to the onset time of thawing being earlier in the former. As the ground ice thaws, surface soil moisture increases, and *E* also increases.

Some studies found that *E* is larger than rainfall at the annual scale at observation sites on the QTP [59,63]; the same phenomenon was also found in our study (see Table 1). Due to missing values of rainfall during the study period at BJ and NAMORS, we instead used data from the nearby Naqu and Dangxiong country meteorological station provided from China Meteorological Administration. The exception of this was when there were observation errors about precipitation or ET; this phenomenon is often explained by surface soil moisture providing available water for evaporation. Although limited by land surface water supply, soil evaporation is one of the main forms of water consumption for the four selected field sites.

**Figure 3.** The daily variations of *E* (blue circles), sublimation (red circles), *Erad* (pink line), and *Epa* (gray area) during the study period at the observation fields: (**a**) TGL; (**b**) XDT; (**c**) BJ; (**d**) NAMORS. "W" represents warm season and "C" represents cold season.



Note: N1 represents available samples of first observational year, N2 represents available samples of second observational year. a Rainfall data on second year from Naqu country meteorological station; b Rainfall data on second year from Dangxiong country meteorological station.

### *3.2. Evaluating Model Performance*

3.2.1. Model Performance with Default and Calibrated Parameter Values on a Daily Timescale

This study first employed the default and calibrated parameter values to simulate the daily actual evapotranspiration (*Esim*) and then compared with in situ measurements, respectively (Figure 4).

**Figure 4.** Comparison of simulated daily actual evapotranspiration (*Esim*) rates by five CR-based models with default parameter against measurements made using EC system at (**<sup>a</sup>**–**<sup>e</sup>**) TGL; (**f**–**j**) XDT; (**k**–**<sup>o</sup>**) BJ; (**p**–**<sup>t</sup>**) NAMORS, respectively. The blue scattered dots and dash lines represent *E* estimated by default parameter values and its fitting line; the red scattered dots and dash lines represent *E* estimated by calibrated parameter values and its fitting line; the black solid line is the 1:1 line.

The results indicate that *Esim* was underestimated to different degrees by each CRbased model at all four selected sites, whether the default or calibrated parameter values were employed; the highest NSE value was 0.92 at TGL and the lowest NSE value was 0.58 at NAMORS. The average RMSE value was 0.4 mm d−<sup>1</sup> at TGL and at XDT, 0.58 mm d−<sup>1</sup> at BJ, and 0.63 mm d−<sup>1</sup> at NAMORS when default values were employed (Table 2). After being locally calibrated, the highest NSE value was 0.92 at TGL and lowest NSE value was 0.69 at NAMORS. The average RMSE value was 0.36 mm d−<sup>1</sup> at TGL, 0.38 mm d−<sup>1</sup> at XDT, 0.56 mm d−<sup>1</sup> at BJ, and 0.6 mm d−<sup>1</sup> at NAMORS (Table 2), with only a small improvement in accuracy when compared with the results simulated by default parameter values. In addition, the CC values of five CR-based models all exceeded 0.87. Therefore, a better performance of the CR approach was found for daily *E* estimation in the frozen regions of the QTP, and unknown parameter values from other similar climatic or underlying surfaces could also be applicable without local calibration due to a lack of in situ measured *E*, such as in data-scarce regions.


**Table 2.** Statistical results of CR-based models with default and calibrated (in parentheses) parameter values against daily actual evapotranspiration measurements from four observation sites. Bold font indicates the best model based on maximum NSE values among five CR-based models.

However, a relatively larger bias occurred when there were lower *E* values for each model at all four sites, which indicated poor applicability of the CR approach. There were even some negative *Esim* values by the K2006 model when *E* was low at all four observation sites, probably because there were no strict boundary conditions constrained in the K2006 model, which led to abnormal values near the border of the boundary. In general, the CR approach may be inappropriate during the cold season in the frozen regions of the QTP because lower *E* values mainly occurred during the cold season. The S2017 model performed poorly among the five CR-based models, with relatively lower NSE values and higher RMSE and MBE values (Table 2), especially at NAMORS, where the MAE values reached 0.57 mm d−<sup>1</sup> and 0.52 mm d−<sup>1</sup> with the default and calibrated parameter values, respectively, signifying the largest bias during the study period. The C2018 model performed best at TGL, BJ, and NAMORS, and performed very well at XDT, and the NSE values of the C2018 model were all above 0.75. The overall performances of the B2015 and H2018 models were slightly inferior to that of the C2018 model, and the NSE and RMSE values of the above two models were close to those of the C2018 model.

Note that for the C2018 model, although there was only one unknown parameter to calibrate, there were more variables, such as wet and dry environment temperature, which meant a relatively more complex calculation process. For the B2015 model, previous studies have pointed out its deficiency in physical boundaries, which may be inappropriate for extremely dry or wet environments and periods. Therefore, the H2018 model, with relatively robust physical boundary constraints and a simpler calculation process, may be the most appropriate CR-based model for daily *E* estimation in the frozen regions of the QTP.

### 3.2.2. Performance of Different CR-Based Functions against Relationships among Three Evapotranspiration Variables

To determine relationships among three evapotranspiration variables (*E*, *Epo* and *Epa*) during warm and cold seasons and to further elucidate the applicability and optimal form of CR-based functions in the frozen ground regions of the QTP, the performance of different function forms using calibrated parameter values combined with observation data is presented in Figure 5. Here, all data at each field site were divided into three parts: daily ET variables that fulfill the preconditions of each CR-based model, judged by " po (*Epo\_c*) ≤ *Epa*" for the K2006, B2015, S2017, and C2018 models, "*E* ≤ *Epo* and *E* ≤ *Epen*" for the H2018 model, and the above data can be used to further divide based on the warm and cold season, respectively. The remaining data are daily evapotranspiration variables that do not fulfill the CR-based models. Note that although parameter calibration was performed, there were still some data that did not satisfy the CR-based models, which may have been due to observation errors or inapplicability of the CR approach to that time.

The results indicated that observation data were unevenly scattered around CRfitting lines at each site, and all CR-based function forms basically captured relationships among the three types of evaporation during the warm season in the frozen ground regions of the QTP, similar to the simulated results from some other sites around the world [36,42,64]. However, when variable y converges to 0 (that is, *E* is low), which often occurs in cold seasons, the dimensionless variable *x* (or *X*) has a relatively wide range, even above 0.6, because strong radiation (*Erad*) on the QTP leads to a high proportion of radiation to apparent potential evapotranspiration (*Epa*). Usually, a high proportion of radiation indicates there is strong evaporation energy available and higher actual evaporation; however, for water-limited frozen ground regions, water availability is another important factor affecting evaporation. Figure 5a–t shows relatively small deviations during warm seasons and the CR-fitting lines cross near the center of the observation data; five different CR-based functions both captured the daily *E* values well without evident differences, however, differences were clear in the lower *E* values. Because most studies about CR approach focus on warm seasons or growing seasons when *E* is relatively high, it is not enough to take into account the conditions of small *E* values during cold seasons; difficulties lie in high uncertainty of observation data during cold seasons, which leads to a lack of reliable data to validate whether the hypothesis of the CR approach is applied to cold and high-altitude areas or cold seasons from the aspect of the mechanism.

#### 3.2.3. Model Performance with Calibrated Parameter Values on a Monthly Timescale

We also investigated the performance of five CR-based models with calibrated parameter values on a monthly timescale. Figure 6 illustrates that the monthly *Esim* go<sup>t</sup> acceptable accuracy at TGL, XDT, and BJ, but relatively larger deviations occurred at NAMORS. The simulated performance between each CR-based model was very close, which indicated more stable performance to the CR approach on longer timescales. Compared with variations in precipitation, the peak of *E* was not always synchronized with the peak of precipitation on a monthly timescale. Taking TGL (Figure 6a) for an example, the peak of *E* occurred in August and was higher than concurrent precipitation; higher *E* may be contributed by water stored in the soil layer by antecedent precipitation, due to the peak of precipitation in July.

**Figure 5.** Scatter plot of y versus different rescaling *x* for five CR-based functions at four observation sites (dots, **<sup>a</sup>**–**<sup>t</sup>**), blue scattered dots represent data during cold season and red scattered dots represent data during warm season, green scattered dots represent data that could not meet the preconditions of CR approach. The solid lines fitted by CR-based functions with calibrated parameter values.

**Figure 6.** Comparison of monthly *E* values estimated by five CR functions with calibrated parameters against measurements by eddy-covariance system at (**a**) TGL; (**b**) XDT; (**c**) BJ; (**d**) NAMORS, respectively. Number in parentheses are available days during corresponding month.

To further determine the contributions of monthly simulated bias to annual total bias at each field site, the monthly absolute bias of each CR-based model was calculated (Figure 7). Although the performance of the five CR-based models varied from one site to another, overall, the H2018 model seemed to have a smaller annual bias than the other models. According to Figure 7e,i,o,t, the H2018 model underestimated annual *E* at each site, and the bias was within 15 mm. The K2006, B2015, and C2018 models performed slightly inferior to the H2018 model; larger biases all occurred at NAMORS, and the rest of the sites performed relatively well. The S2017 model performed worst among the five CR-based models at each site, and the largest bias occurred at NAMORS, which was close to 100 mm on an annual timescale. Considering the contribution of monthly or seasonal bias to annual total bias (see Table 3), negative annual bias was mainly contributed by the negative bias of the warm season at TGL; however, the negative bias of the cold season contributed more to the total bias at XDT, BJ, and NAMORS. The larger negative monthly bias frequently occurred in October or November, and the maximum bias was −15.3 mm (in October) at TGL, −5.42 mm (in July) at XDT, −11.72 mm (in July) at BJ, and −17.58 mm (in October) at NAMORS.

**Figure 7.** Bias of monthly *E* estimated by five CR-based models at (**<sup>a</sup>**–**<sup>e</sup>**) TGL; (**f**–**j**) XDT; (**k**–**<sup>o</sup>**) BJ; (**p**–**<sup>t</sup>**) NAMORS, respectively. Dash lines represent annual total bias. The grey filled area represents the warm season.

**Table 3.** Bias of simulated *E* by CR-based models with calibrated parameter values against measurements by EC on different timescales (monthly, seasonal, annual) at TGL, XDT, BJ and NAMORS. Grey-filled area is warm season period (May to September), bold font stands for model with minimum absolute bias during corresponding period.

