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

The actual evapotranspiration (ETa) estimated models based on the complementary relationship (CR) theory have been applied in various climatic conditions around the world. However, in cold regions, the evaluation of the adaptability of the CR models was performed through complete freeze-thaw cycles, and the adaptability during various periods of soil freeze-thaw cycles has not been evaluated separately. Daily ETa was measured by lysimeters on alpine grassland in the Qilian Mountains from 2010 to 2017, and the measurements were used to evaluate five CR models during the thawing, thawed, freezing, and frozen periods, respectively. The five models comprised the advection-aridity (AA) model of Brutsaert and Stricker, the GG model proposed by Granger and Gray, Morton’s CR areal evapotranspiration (CRAE) model, the Han model, and Brutsaert model. The results show that all five CR models were only able to estimate daily ETa during the thawed period. None of the models could estimate the daily ETa during the thawing, freezing, or frozen periods. The basic assumptions of the CR may not be suitable for non-thawed periods with complex energy processes, and no complementary behavior was shown in the non-thawed periods. The CR models must be applied with caution during freeze-thaw cycles in cold regions.

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 2 (∼52.5%) of the global land area undergoes a seasonal transition of freezing-thawing conditions every year (Kim et al. ; Rowlandson et al. ), 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. ).
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. ). And in the seasonal frozen region, the soil freeze-thaw process could be divided into four periods: thawing, thawed, freezing, and frozen 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 ), and high latitude areas such as the arctic show remarkable change in the water cycle under climate changing (Bring et al. ). 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. ; Yang et al. ). 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.

Study area
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 0 E, 38 16.1 0 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 half-hour meteorological data were processed into daily mean values for this study.

Background of the CR
The complementary principle involves the relationship between ET a , potential ET (ET po ), and 'apparent' potential ET (ET pa ) (Brutsaert ). The CR theory proposed first as a symmetric linear function by Bouchet (), and Brutsaert & Parlange () broadened it to an asymmetric linear function: where ET pa is the 'apparent' potential evaporation, which describes the ET that occurs from a small saturated surface; ET po is the potential ET or the wet environment ET, which describes the ET that takes place from the extensive wellwatered surface; and ET a is the actual evapotranspiration.
Definitions of and differences between the three ET were given by Brutsaert (). b is the proportionality.
By normalizing Equation (1) with ET pa , Equation (1) can be expressed as follows: where f is a function of the quantity inside the brackets. The differences in the CR models are (1) how to define function f and (2) how to calculate ET pa and ET po .
AA model. Based on the consideration of 'advection-aridity', Brutsaert & Stricker () derived the AA model, which can be expressed as follows: where ET AA a is ET a estimated by the AA model and ET AA pa is calculated by combining information from the energy budget and water vapor transfer in the equation proposed by Penman (): where Δ is the slope of the saturation vapor pressure curve at the air temperature (kPa/ C), γ is the psychometric constant (kPa/ C), R n is the net radiation (MJ/m 2 /day), G is the soil heat flux (MJ/m 2 /day), λ is the latent heat of vaporization (MJ/kg), and e s and e a are the saturation and actual vapor pressure at air temperature (kPa), respectively. f (U) is the wind function, which includes the 2-m wind speed (U 2 , m/s), and should be derived via the Monin-Obukhov similarity theory as follows: where α is the parameter and other symbols are as defined previously.
GG model. For the GG method, Granger & Gray () introduced the concepts of relative ET (G*) and showed that an equation similar to that of Penman could also be derived following the CR approach: where ET GG a is ET a estimated by the GG model, E a is the drying power of the air, which in general can be written as follows: Based on daily estimated values of ET a from the water balance, Granger & Gray () showed that a unique relationship exists between G* and a parameter that they called the relative drying power D, given as follows: Later, Granger () modified Equation (13) to CRAE model. For the CRAE method, ET a can be estimated using Equation (15): where ET CRAE a is ET a estimated by the CRAE model; however, the calculation of ET pa and ET po is different. To calculate ET pa , Morton () decomposed the Penman equation into two separate parts to describe the energy balance and the vapor transfer process. The most prominent advantage of CRAE is that it does not require wind speed as an input (Ma et al. (b)). The 'equilibrium temperature', T p , which is defined as the temperature at which the energy balance equation and the vapor transfer equation for a moist surface give the same result, was introduced to calculate the potential ET.
where f v is a vapor transfer coefficient (W/m 2 /mbar), ε is the surface emissivity, p is the atmospheric pressure (mbar), σ is the Stefan-Boltzmann constant (W/m 2 /K 4 ), T p and T a are the equilibrium and air temperatures ( C), e p is the saturation vapor pressure (mbar) at T p , e d is the saturation vapor pressure (mbar) at dew point temperature, λ is the latent heat of vaporization (W day/kg), γ is the psychometric constant (mbar/ C), and R n and G are as defined previously but in units of W/m 2 . T p can be obtained via iterations from Equations (16) and (17). Note that the units used in the CRAE model are adopted here to ensure that the correct values of the empirical constants are included.
f v is given by: where p s is the sea-level atmospheric pressure (mbar), f z is a constant, and ξ is a dimensionless stability factor estimated from: where e a is the saturation vapor pressure (mbar) at air temperature, Δ is the slope of the saturation vapor pressure curve (mbar/ C) at air temperature, and b 0 ¼ 1.0 for the CRAE model.
To estimate ET po , Morton () also modified Equation (9) by Priestley & Taylor () to account for the equilibrium temperature dependence of both the available energy and the slope of the saturation vapor pressure curve: where (R n À G) p is the net available energy for the surface at ET rad is the proportion of the radiation term in ET Han pa : Considering comparison with other models, Equation (22) was derived a modified form: In Equations (22) and (27), m and n are empirical coefficients: m ¼ In this study, x min ¼ 0 and   (33)) and a supposedly more flexible quartic formulation with an additional adjustable parameter (Equation (34)): where ET B a is ET a estimated by the Han model and c is an adjustable parameter. When c ¼ 0, Equation (34)

Actual evapotranspiration measurement
The lysimeter is the standard instrument used to measure where ΔW is the weight difference measured by the lysimeter (g/day), ρ is the water density (1 g/cm 3 ), S is the surface area of the soil within it (cm 2 ), P is the precipitation (mm/day), and I is the infiltration (g/day), measured by weighing daily, as with the lysimeters.

Period division
The complete soil freeze-thaw cycle can be divided into four periods: thawing, thawed, freezing, and frozen. The soil temperatures at depths of 0 and 40 cm were chosen to divide the periods. The soil temperature at 0 cm was used to indicate whether the surface soil was frozen or thawed.
The soil temperature at 40 cm was selected because ET basically occurs near the surface and because of the vegetation root depth and the size of the lysimeters. The criteria for division were as follows: T 0 >0 and T 40 >0 Thawed T 0 0 and T 40 >0 Freezing T 0 0 and T 40 0 Frozen T 0 >0 and T 40 0 Thawing where T 0 and T 40 represent the soil temperature at depths of 0 and 40 cm ( C), respectively. Figure 2 shows an example of the four-period division in 2017.

Estimation of parameters
The five models can be divided into two groups: (1) Table 1.

Evaluation criteria
Statistical indices were used for the quantitative analysis of the ET a modeling performance. The ET a values measured by the lysimeters and estimated by five CR models were compared using a series of statistics. The errors were calculated as follows: where R 2 , MAE, RMSE, and NSE are the coefficient of determination, the mean absolute error (mm/day), the root-mean-square error (mm/day), and the Nash-Sutcliffe efficiency, respectively. n is the number of statistical days, and cov and σ are the covariance and standard deviation, respectively, E is the estimated ET a values by CR models (mm/day), and M represents the ET a values measured by the lysimeters (mm/day).
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.   frozen, and thawing periods, none of the five models were able to estimate the daily ET a at the research site.

CR during soil freeze-thaw cycles
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.  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 ). 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. ); thus the energy was used not only to evaporate the water into vapor but also in the phase change of icewater. 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. ). 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. ). 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.

Statistics traps in model performance
Although  (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 nonthawed 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

CONCLUSIONS
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 nonthawed 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.