Evaluation of five complementary relationship models for estimating actual evapotranspiration during soil freeze-thaw cycles

Evapotranspiration (ET) is an important part of the global hydrothermal balance, consuming about 60% of precipitated water on land back to the atmosphere and more than 50% of the net radiation available on land (Oki & Kanae 2006; Jung et al. 2010; Hirschi et al. 2017; Panwar et al. 2019). It is the core link of surface water, carbon, and energy exchange (Wang & Dickinson 2012) and an important indicator of regional climate change and water cycle research (Wang et al. 2013; Sörensson & Ruscica 2018). However, direct observations of actual evapotranspiration (ET_a) remain unavailable in many locations around the globe (Han et al. 2014; Liljedahl et al. 2017; Kim et al. 2019), and reliable estimates of ET are crucial for a better understanding of current and future water/energy/carbon budgets (Sörensson & Ruscica 2018).

Several methods have been proposed for estimation of ET_a (McMahon et al. 2016), but the models based on the complementary relationship (CR) proposed by Bouchet (1963) are usually preferred because they require only standard meteorological variables and avoid extremely detailed knowledge of the complex processes and interactions between the soil, the vegetation, and the near-surface boundary layer (Xu & Singh 2005; Huntington et al. 2011). A variety of models based on the CR have been developed in recent decades, and examples of widely known models include the advection-aridity (AA) model proposed by Brutsaert & Stricker (1979), the Granger and Gray (GG) model (Granger & Gray 1989), and the CR areal evapotranspiration (CRAE) model derived by Morton (1983). All three models are linear function for CR. Han et al. (2012) developed a nonlinear approach and extended it to a sigmoid generalized complementary function (hereafter referred to as the Han model) (Han & Tian 2018). Inspired by Han et al. (2012), Brutsaert (2015) proposed a generalized nonlinear form of the CR with physical boundary conditions (hereafter referred to as the Brutsaert model).

These five CR models have been used and discussed in various climatic conditions throughout the world. However, compared with warm regions, those five models are rarely discussed in cold regions, where the soil freeze-thaw process affects the hydrothermal processes including ET. Approximately 66 million km² (∼52.5%) of the global land area undergoes a seasonal transition of freezing-thawing conditions every year (Kim et al. 2012; Rowlandson et al. 2018), and a soil freeze-thaw process refers to the physical geological phenomenon that the surface soil temperature drops below 0 °C and then rises above 0 °C (Zhao et al. 2021). The freeze-thaw process of the active layer in the permafrost region could be divided into four stages: summer thawing stage, autumn freezing stage, winter cooling stage, and spring warming stage (Zhao et al. 2000). And in the seasonal frozen region, the soil freeze-thaw process could be divided into four periods: thawing, thawed, freezing, and frozen (Ding et al. 2017). Ma et al. (2015a) evaluated several CR models in the alpine steppe of the Tibetan Plateau, but only during the warm season when the soil had thawed. Hu et al. (2018) evaluated the Brutsaert model across four vegetation types on Mount Gongga in southwestern China. The model evaluation was on a daily scale for a full year, with no respective evaluation of different periods during the soil freeze-thaw cycles. Some researchers have used CR models to calculate ET_a on regional and continental scales (Zhu et al. 2016; Kim & Kaluarachchi 2017; Szilagyi et al. 2017; Liu et al. 2018; Szilagyi & Jozsa 2018; Ma et al. 2019), portions of where certainly must have undergone freeze-thaw cycles. Still, no separate evaluation of the adaptability of the CR models in different periods during soil freeze-thaw cycles has been performed.

Soil freeze-thaw processes are common on the earth's land surface. Permafrost covers approximately one-quarter of the Northern Hemisphere's land area (Obu et al. 2019), and seasonal frozen ground appears throughout the world. Cold regions in which soil freezing-thawing processes are significant and severe are usually located at high elevations and high latitudes. High elevation areas such as the Tibetan Plateau and other mountainous regions are generally important sources of freshwater (Dettinger 2014), and high latitude areas such as the arctic show remarkable change in the water cycle under climate changing (Bring et al. 2016). ET and the water cycle are poorly understood in these cold regions because their extremely harsh natural environmental conditions lead to a scarcity of measurement data (Liljedahl et al. 2017; Yang et al. 2019). Although less ET may occur in the cold season with frozen soil than in the warm season when soil is thawed, it is an integral and indispensable part of the overall ET and water cycle, and it is necessary to evaluate the estimated ET_a throughout the freeze-thaw processes.

Consequently, using in situ long sequence daily artificial lysimeters data and automatic meteorological observations in Qilian Mountains, the aims of this study are (1) to quantify ET_a in a cold high elevation mountainous area and (2) to evaluate the ability of five CR models to estimate daily ET_a during soil freeze-thaw cycles in cold regions.

The experiment was conducted at the Qilian Alpine Ecology and Hydrology Research Station, located in the headwaters region of Heihe River in the Qilian Mountains (Figure 1). A standard meteorological observation field (99°52.9′E, 38°16.1′N, 2,980 m) was constructed at the station on 20 June 2009, and an automatic meteorological tower was installed to record half-hourly measurements of meteorological variables including air temperature, relative humidity, wind speed and direction, radiations, soil heat flux, and multiple layers of soil temperature and water content. According to the meteorological data from the observation field and automatic meteorological tower, from 2010 to 2017, the mean air temperature and relative humidity were approximately 0.9 °C and 54.0%, respectively. The mean annual precipitation was 487.1 mm; more than 90% of which occurred during the plant growing season (May to September). Rainfall was the main form of precipitation, and snowfall accounted for about 9.6% of the total precipitation. According to field observation, even if there was snowfall, snow was difficult to accumulate, and most of them melt in 1 day. The existence of snow cover was rare, so the effect of snow was not considered in this study. The research site was classified as seasonal frozen soil, and the maximum frozen soil depth was about 240 cm. The dominant species in the area were Kobresia capillifolia (Decne.) C.B. Clarke and Carex moorcroftii (Falc. ex Boott). The canopy cover exceeded 98%, and the vegetation root depth was nearly 35 cm. Following the method recommended by Allen et al. (1998), the half-hour meteorological data were processed into daily mean values for this study.

In this study, x_min = 0 and x_max = 1 because the complementary curve is insensitive to the value of x_min and x_max (Han & Tian 2018). More detailed information on the procedure was given by Han & Tian (2018).

The five models can be divided into two groups: (1) uncalibrated models (GG and CRAE models), all of the parameters were used according to the values measured or recommended by the existing models; (2) calibrated models (AA, Han, and Brutsaert models), key parameters need to be calibrated. The Priestley–Taylor parameter α is an important parameter for the AA, the Han, and the Brutsaert models; and it is not quite the Priestley–Taylor parameter but is merely a weak analog of it (Crago & Qualls 2018; Hu et al. 2018). It is usually has a default value of 1.26 (Priestley & Taylor 1972; Ma et al. 2015a; Zhang et al. 2017), and the values of α range between 0.87 and 1.44 in various climatic conditions in China (Liu et al. 2016). Since the original AA model version was selected in this study, the α was an only parameter need to be calibrated. In addition to Priestley–Taylor parameter α, b of the Han model in Equations (30) and (32) and c of the Brutsaert model in Equation (34) need to be calibrated, respectively. Considering the seasonal variability of α in the CR (Yang et al. 2013; Liu et al. 2016), the monthly variation in the mean parameters were calibrated, with parameters optimized by minimizing the mean absolute error (MAE) of the estimated actual evaporation (Legates & McCabe 1999; Han & Tian 2018).

The measured ET_a and the corresponding meteorological data from 2010 to 2011 were selected to calibrate the model parameters, and the data from 2013 to 2017 were selected to evaluate the five CR models (the data for 2012 were excluded because of gaps). The model parameters were calibrated daily and then processed into monthly mean values. During the calibration, the parameter b of Han and c of Brutsaert were unconstrained, and the parameters of α of AA, Han, and Brutsaert models were bounded between 0.87 and 1.44. The calibrated monthly parameters values are shown in Table 1.

Note that, to avoid the possible influence of precipitation events and snow sublimation on evaluation, ET_a measurements from days with precipitation were not considered in the evaluation of the five CR models.

The mean values of ET_a measured by two lysimeters were taken as the ET_a at the research site, and Figure 3 shows the variations in ET_a and precipitation in various periods of the soil freeze-thaw cycles from 2010 to 2017. The mean annual ET_a was 441.5 mm. ET_a accounted for 91.7% of precipitation in the corresponding period and represented the largest loss of water cycle on the research site. In the years with low precipitation, such as 2010 and 2012, ET_a even exceeded precipitation because soil freezing could store precipitation in abundant years for ET in the next year (Sugimoto et al. 2003; Ohta et al. 2008). The vast majority of ET_a and precipitation both occurred during the thawed period. The mean of the ratio of ET_a during the thawed period and for the whole year was 91.4%, and that value for precipitation was 93.5%. The mean daily ET_a values for the thawed, freezing, frozen, and thawing periods were 2.19, 0.28, 0.16, and 0.41 mm/day, respectively.

Figures 4 and 5 show the performance of the estimation of daily ET_a using five CR models for the various freeze-thaw periods, and Table 2 shows the values of the evaluation criteria. The performance of the five models was significantly better during the thawed period than during the other three periods, with significant more R² and NSE. In summary, during the thawed period, the calibrated models, with mean MAE, RMSE, and NSE of 0.67 mm/day, 0.83 mm/day, and 0.76, respectively, showed the better performances than uncalibrated models, with mean MAE, RMSE, and NSE of 0.79 mm/day, 1.02 mm/day, and 0.65, respectively. There were slight differences in the model performance of the calibration models. The Han and Brutsaert models with two calibrated parameters showed the better performance than the AA model with one calibrated parameter, and the Han model showed slight better performance than the Brutsaert model. The performance order of the five models during the thawed period was Han model, Brutsaert model, AA model, GG model, and CRAE model. In the freezing, frozen, and thawing periods, none of the five models were able to estimate the daily ET_a at the research site.

Based on the fact that the calibration models showed better performances than the uncalibrated models, the CRs of different periods of the freeze-thaw cycle were analyzed by the calibration models. Figure 6 shows the dimensionless form of the CR for the four freeze-thaw periods. The observed daily ET_a showed significant complementary behavior only during the thawed period; no linear or nonlinear CR was seen in the other three periods.

Figure 7 shows the scaled daily rates of actual evapotranspiration (ET_a/ET_po) and scaled apparent potential evapotranspiration (ET_pa/ET_po) against the moisture index (ET_a /ET_pa) for the four freeze-thaw periods. Still, complementary behavior was only shown during the thawed period, not in the other three periods.

As the Priestley–Taylor parameter α, an important parameter for the AA model, Han model, and Brutsaert model, was calibrated using local meteorological data, it is necessary to analyze the sensitivity of this parameter. The calibrated monthly α in this study ranged from 1.02 to 1.08, which is similar to the values reported in other studies in which this formulation was used in high elevation regions (Ma et al. 2015a; Hu et al. 2018) and ranged from 0.87 and 1.44 in various climatic conditions in China (Liu et al. 2016). To analyze the influence of α on the model's performances, a series of values from 0.87 to 1.44, at intervals of 0.01, were introduced to Equations (9), (29), and (38), and ET_a was then estimated using the AA, Han, and Brutsaert models. Figure 8 shows the variation in NSE produced by the AA, Han, and Brutsaert models relative to changes in α for various freeze-thaw periods. During the thawed period, the NSE from all three models could approach 1 with a reasonable α. During the other three periods, with α values from 0.87 to 1.44, the NSE from all three models were less than 0. This sensitivity analysis of α indicated that the CR models could not estimate ET_a during the freezing, frozen, and thawing periods for any value of a.

The goal of this study was to evaluate the use of CR models to estimate daily ET_a in cold regions, and the results show that the CR models were applicable only during the thawed period; in other words, the CR was only shown during the thawed period. The CR theory was first proposed by Bouchet (1963), and later broadened by Brutsaert & Parlange (1998). The core assumption is that when water is limited, the energy difference between ET_pa and ET_po, has a linear or nonlinear relationship with the energy difference between ET_pa and ET_a. However, this assumption was proposed on intuitive grounds, and its validity has never been rigorously justified (Brutsaert 2015). When the soil is thawed, the ET_a process in cold regions may be similar to that in warm regions, and the assumption is suitable for this condition. During the freezing and thawing periods, the prominent characteristic of the cold region is the intraday freeze-thaw cycles (Guo et al. 2011); thus the energy was used not only to evaporate the water into vapor but also in the phase change of ice-water. When the soil was frozen, unfrozen liquid water content in soil was very low, and the sublimation of solid water, whose physical and energy consumption processes differ from those of evaporation, was the main composition of ET (Knowles et al. 2012). In general, the energy consumption pattern during the non-thawed periods differs from that during the thawed period because in addition to soil evaporation, energy is used for the phase change of ice-water and sublimation. The altered plant transpiration affected by freeze-thaw cycles can also change the energy balance in cold regions (Cable et al. 2014). These complex energy processes may be inconsistent with the basic assumptions of the CR, which may explain why none of the five CR models could estimate the daily ET_a during the freezing, frozen, and thawing periods.

Although the applicability of the CR models during various periods of the freeze-thaw cycles in cold regions has not been evaluated separately, some studies have compared the estimated and measured daily ET_a values for the complete freeze-thaw cycles, and the results showed good model performance (Zhu et al. 2016; Hu et al. 2018; Liu et al. 2018). However, a statistics trap may apply to these results because they evaluated the models across the complete freeze-thaw cycle. The results of this study show that the CR models were not able to estimate the daily ET_a during the non-thawed periods (Figure 5; Table 2). When the statistics were extended to all of the data, all five models could estimate the daily ET_a at the research site very well. The model performance for the complete freeze-thaw cycle was even better than that for the thawed period, with a larger R², lower MAE and RMSE values, and a larger NSE values (Figure 9; Table 3). The reason for this statistical trap is that the ET that occurred during the thawed period accounted for the vast majority of the total ET for the entire freeze-thaw cycle (Figure 3). If the model performance was tested on a monthly scale, all five models could estimate ET_a for the thawed months (May, June, July, August, and September), but not the non-thawed months (January, February, March, April, October, November, and December) (Figure 10; Table 4). Like the daily scale, a statistical trap was found when the statistics were extended to all months, and the model performance for the complete freeze-thaw cycle was better than that for the thawed period. The evaluation of the model performance for the complete freeze-thaw cycle does not indicate that the CR models can calculate the daily or monthly ET_a values separately for different periods.

Daily ET_a during the thawed period is significantly higher than that during the freezing, frozen, and thawing periods, and ET_a that occurs during the thawed period accounts for the vast majority of the total ET_a. All five CR models can estimate daily ET_a for the thawed period. The performance order of the five models were Han model, Brutsaert model, AA model, GG model, and CRAE model The basic assumptions of the CR may not be suitable for the non-thawed periods with complex energy processes. None of the five models could estimate the daily ET_a value for the freezing, frozen, and thawing periods, and no complementary behavior was seen during the non-thawed periods.

The environment in cold areas is extremely harsh, and ET_a measurements in such regions are seriously lacking. Although CR models cannot estimate daily ET_a during the non-thawed periods, estimated ET_a remains an alternative means of understanding ET and water cycle in cold regions because CR models can effectively calculate ET_a during the thawed period, which accounts for most of the total ET_a in those regions. It should be noted, however, that estimations of ET_a from CR models during the non-thawed periods must be treated with caution.

This work was carried out with financial support from the National Key Research and Development Project (2019YFC1510505), National Natural Sciences Foundation of China (41690141 and 41671029), and Joint Research Project of Three-River Headwaters National Park, Chinese Academy of Sciences and the People's Government of Qinghai Province (LHZX-2020-11). The measurement data used in this study are provided by Qilian Alpine Ecology and Hydrology Research Station (http://hhsy.casnw.net/) or can be downloaded freely from the figshare website (https://doi.org/10.6084/m9.figshare.9901151.v2). The authors thank Mr Baoshan Ma for field lysimeters measurement in the harsh natural environmental condition.

All relevant data are available from an online repository or repositories (https://doi.org/10.6084/m9.figshare.9901151.v2).