Comparison of Shuttleworth–Wallace and Dual Crop Coefficient Method for Estimating Evapotranspiration of a Tea Field in Southeast China

Abstract

Determination of evaporation (E) and transpiration (T) in tea fields separately is important in developing precise irrigation scheduling and enhancing water use efficiency. In this study, the Shuttleworth–Wallace (S-W) model was applied to simulate the variations of E and T based on the data from 2015 to 2018 in a tea field in southeast China. The dual crop coefficient (D-K) method recommended by FAO-56 was also applied to calculate E and T, using the same data set to compare with the S-W model. The measured crop coefficient (K_c) ranged from 0.43 to 1.44 with the average value was 0.90 during 1–150 DOY (days of year), and the measured K_c tended to be stable with the average value of 0.83 during 151–365 DOY in 2015. The S-W model estimated ET_c with root mean square error (RMSE) and R² of 0.45 mm d⁻¹ and 0.97, while for the D-K method the values were 0.61 mm d⁻¹ and 0.95. Therefore, both approaches could estimate the E and T separately in tea fields in southeast China, however, the D-K method had a slightly poorer accuracy compared to the S-W model in the estimation of ET_c.

1. Introduction

Tea, as one of the oldest (since 3000 BCE) and most popular nonalcoholic caffeine-containing beverages, has been integrated into the social and daily life in China. At present, tea is commercially cultivated in more than 3.80 million hectares of land on a continent-wide scale, and 5.56 million metric tons of tea world-wide were produced annually in 2014 [1]. Irrigation is essential to ensure tea production in southeast China, and the appropriate amount of irrigation water at the right time directly increases tea quality. However, most of the tea farmers in southeast China are still using an empirically determined irrigation quota due to lack of accurate irrigation basis for tea plants. Therefore, accurate determination of tea evapotranspiration (ET_c) is urgently needed to develop precise irrigation scheduling and enhance water use efficiency in this region [2,3].

Single source models, such as the Penman–Monteith (PM) method [4,5] pan evaporation method [6], Takakura method [7] and Priestley–Taylor method [8], are often used to determine ET_c for grain crop fields [9]. However, unlike grain crops, tea fields always have wide row space (120 cm) between tea plants so that large soil surface is uncovered by the tea plants. Hence, the above single source ET_c models that do not consider independently soil evaporation (E) are especially problematic because the soil evaporation ratio (E/ET_c) is large [10]. Many researchers have reported that the soil evaporation (E) should not be neglected in estimating ET_c and it accounts for 20–30% of the ET_c for a cherry orchard [10], 25% of ET_c for greenhouse tomato [11], 17.8–26.4% of ET_c for a mulched greenhouse hot pepper [9], and 7.75–21.87% of ET_c for greenhouse cucumber [12]. Accordingly, the two source Shuttleworth–Wallace (S-W) model [13] and the dual crop coefficient (D-K) method [14], which can estimate transpiration (T) and E separately, have been validated to estimate ET_c in different ecosystems [11,15,16,17,18,19].

The S-W model is the first analytical model, combining T and E by formulating the different media through which evaporative flux travels as resistances [20]. Three aerodynamic resistances among different interface (bulk boundary layer resistance: $r_a^a$, aerodynamic resistance between mean canopy flow and reference height: $r_a^c$, aerodynamic resistance between soil surface and mean canopy flow: $r_a^s$) are introduced simultaneously to regulate the transfer between these surfaces and the atmosphere [11,21]. Canopy and soil surface resistances ($r_s^c$ and $r_s^s$) are introduced simultaneously into the S-W model to regulate the transfer of energy from plants and soil, respectively [11]. Ortega-Farias et al. [22] and Zhao et al. [19] parameterized the S-W model by adopting a Jarvis type $r_s^c$ formulation and power function between $r_s^s$ and soil water content at the top layer to estimate vineyard ET_c. Zhu et al. [17] parameterized the S-W model by adopting a Leuning type $r_s^c$ formulation and exponential function between $r_s^s$ and soil water content to estimate alpine grassland ET_c on the Qinghai–Tibetan plateau. Huang et al. [12] parameterized the S-W model based on the measured stomatal resistance and E, by adopting an exponential function between $r_s^c$ and net radiation and a power function between $r_s^s$ and soil water content to estimate cucumber ET_c in a Venlo-type greenhouse. Previous studies indicated that the S-W model performed well in estimating ET_c of vineyard [19,22,23], paddy field [24], orchard [10,25], greenhouse tomato [11], and maize [26,27], but it is challenging to parameterize the two resistances.

To simplify the parameterization process for modeling T and E separately, the Food and Agricultural Organization of the United Nations (FAO) developed an indirect method, named dual crop coefficient (D-K) method, to estimate ET_c [14]. The D-K method, allows the separation of E and T and divides crop coefficient (K_c) into basal crop coefficient (K_cb) and soil evaporation coefficient (K_e) [11,19]. Due to its practical simplicity, the D-K method is widely adopted to estimate ET_c in different types of sparse crops and climatic regions [26]. Despite the increasing relevance of using the D-K method for modeling crop ET_c, determination of K_cb and K_e is challenging [28]. The straightforward adoption of generalized K_cb recommended by FAO-56 can lead to errors in the estimation of ET_c and its components, because the dividing of crop growth period and associated crop coefficients are closely related to local climate and crop condition [19]. Zhao et al. [29] observed K_cb for maize were 0.20, 1.10, and 0.45 (at the initial, mid, and end stage), respectively, in the semiarid to sub-humid climate of the North China Plain. Contrastingly, Miao et al. [30] reported that observed K_cb for maize were 0.10, 1.15, and 0.25 (at the initial, mid, and end stage), respectively, in the desert climate of the Hetao irrigation district. Rosa et al. [31] also reported the observed K_cb for maize were 0.07, 1.15, and 0.20 (at the initial, mid, and end stage), respectively, in the dry sub-humid climate of southern Portugal.

While previous studies have reported that the S-W and D-K methods had different performances for various crops grown in different regions, no study has been conducted to assess and compare the performances of the S-W and D-K methods for estimating ET_c for perennial plants like tea in southeast China. Therefore, the purpose of this study was to (1) parameterize the resistances and coefficients of soil-atmosphere and plant-atmosphere interface of the S-W and D-K model based on detailed meteorological data from a tea field in southeast China; (2) compare and evaluate the performances of the S-W and D-K method in estimating ET_c based on the measured ET_c by Bowen ratio energy balance method, so as to provide accurate estimation of ET_c for tea fields in southeast China.

2. Materials and Methods

2.1. Field Observation

The experiment was conducted in a tea field located at Jiangsu province, China (31°65′ N, 119°23′ E, 23 m a.s.l) from 2015 to 2018. The experiment site is located in a subtropical monsoon climate, average annual precipitation is 1058.8 mm, and about 750 mm occurs in spring and summer; mean, maximum, and minimum air temperatures through the year are 15.2 °C, 36.6 °C, and −5.7 °C, respectively. The soil text in the study area is medium loam and the field capacity is 28% [32].

Anji White tea, which is one of the cultivars of tea in China, was used for the experiment. The size of the observation field was 55 × 58 m and was surrounded by other tea fields, with the total area larger than 7000 m² [32]. The tea plants were transplanted into the field on 5 March 2014, with 120 cm spacing between the rows, and 52 cm spacing within a row. According to the FAO-56 and the actual measurements of the tea plants, the period from transplanting to the end of 2018 was seen as the initial growing stage of the tea plants, and the average leaf area index (LAI) was around 2.5.

The latent heat flux of tea field was measured by a Bowen ratio energy balance system installed in the center of field. The system consists of two layers of high-accuracy temperature and humidity sensors, radiometers, a soil heat flux plate, a three-cup anemometer, a solar power supply equipment, and a data logger. Two layers (1.5 and 2.0 m above the ground surface) of air temperature (T_a) and relative humidity (RH) were measured with sensors HMP155A (Vaisala, Finland) with high measurement accuracies (±0.2 °C from −10 to 40 °C for T_a, ± 1 % from 0 to 90 % for RH). The absolute error between the two HMP155A sensors was calibrated by setting them at the same height before the field observation. The radiometers (CNR-4, Kipp and Zonen, Amsterdam, The Netherland), which can measure downward shortwave/longwave and upward shortwave/longwave radiation separately, were installed at 2.5 m above the ground. The calculations of longwave radiation components were corrected as a thermal effect caused by instrument heating. Soil heat flux was measured at 2 cm depth with a soil heat plate HFP01-L10 (Campbell Scientific, Logan, UT, USA). Wind speed was measured with a three-cup anemometer A100L2 (MetOne, New York, NY, USA) at 2.0 m above the ground. Soil volumetric water content and soil temperature were measured by Hydra Probe sensors by setting the sensors at five different depths within the soil layer (5, 10, 20, 50, and 70 cm). All the data were obtained every 1 s and recorded by a data logger CR3000-NB (Campbell Scientific, USA). Due to some technical problems, the data 335-356 days of year (DOY) in 2016 and 193-300 DOY in 2017 were lost.

In this study, the direction of the prevailing winds during the growing season were westerly, the influence of fetch was not considered due to the similar coverage, and the irrigation intensity for 200 m of upwind direction of the observation field [32].

2.2. Bowen Ratio Energy Balance Method

The latent heat flux can be determined by the Bowen ratio energy balance method as follows:

$$\mathit{LE}=\frac{R_n - G}{1+\beta }$$

(1)

$$\beta =\gamma \frac{\mathsf{\Delta }T}{\mathsf{\Delta }e}$$

(2)

where LE is the latent heat flux (W m⁻²), R_n the net radiation (W m⁻²), G is the ground heat flux (W m⁻²), the ΔT and Δe are the temperature (°C) and vapor pressure (kPa) difference between the two measurement layers, respectively, and γ is the psychrometric constant (kPa °C⁻¹).

2.3. Shuttleworth–Wallace (S-W) Model

The evapotranspiration (ET_c), tea plants transpiration (T), and soil evaporation (E) can be calculated by the following expressions based on the S-W model [13]:

$${\mathit{ET}}_c=T+E$$

(3)

$${T=C}_c\frac{{\mathsf{\Delta }A+\left(\right(\rho }_a{c}_{p}D-{\mathsf{\Delta }r}_a^c{A}_s{)/(r}_a^a{+r}_a^c\left)\right)}{{\mathsf{\Delta }+\gamma (1+(r}_s^c{/(r}_a^a{+r}_a^c\left)\right))}$$

(4)

$${E=C}_s\frac{{\mathsf{\Delta }A+\left(\right(\rho }_a{c}_{p}D-{\mathsf{\Delta }r}_a^c(A-A_s{\left)\right)/(r}_a^a{+r}_a^c\left)\right)}{{\mathsf{\Delta }+\gamma (1+(r}_s^s{/(r}_a^a{+r}_a^s\left)\right))}$$

(5)

$$C_c={\left. [1+\frac{R_c{R}_a}{R_s{(R}_c{+R}_a)}\right]}^{-1}$$

(6)

$$C_s={\left. [1+\frac{R_s{R}_a}{R_c{(R}_s{+R}_a)}\right]}^{-1}$$

(7)

where C_c is canopy resistance coefficient; C_s is soil surface resistance coefficient; Δ is the slope of the saturation vapor pressure curve at temperature (kPa °C⁻¹); c_p is specific heat of the air at constant pressure (= 1013 J kg⁻¹ °C⁻¹); ρ_a is air density (kg m⁻³); D is water vapor pressure deficit (kPa); $r_s^c$ is the canopy resistance (s m⁻¹); $r_a^c$ is bulk boundary layer resistance of the vegetative elements in the canopy (s m⁻¹); $r_a^a$ is aerodynamic resistance between mean canopy flow and reference height (s m⁻¹); $r_a^s$ is aerodynamic resistance between soil surface and mean canopy flow (s m⁻¹); and $r_s^s$ is soil surface resistance (s m⁻¹). The S-W model adopts the concept of a bulk boundary layer resistance, $r_a^c$, which controls transfer between the surface of vegetation and the canopy air steam. Vertical transport is controlled by two further aerodynamic resistances ($r_a^a$ and $r_a^s$). $r_a^a$ is the transfer resistance between the hypothetical mean canopy flow and the reference height above the crop. $r_a^s$ is the aerodynamic resistance encountered by the energy fluxes leaving the substrate before they are incorporated into the mean canopy flow [13]. For simplicity, it is assumed that the various aerodynamic resistances are identical for sensible and latent heat. More detailed approaches that parameterize the three aerodynamic resistances ($r_a^c$, $r_a^a$, and $r_a^s$) can be found in [13].

Values of R_a, R_s, and R_c were calculated as follows:

$$R_a{=(\mathsf{\Delta }+\gamma )r}_a^a$$

(8)

$$R_s=\left(\mathsf{\Delta }+\gamma \right)r_a^s{+\gamma r}_s^s$$

(9)

$$R_c=\left(\mathsf{\Delta }+\gamma \right)r_a^c{+\gamma r}_s^c$$

(10)

A and A_s (W m⁻²) are the available energy leaving the canopy and soil surface, respectively, and were calculated as

$${A=R}_n-G$$

(11)

$$A_s{=R}_{\mathrm{ns}}-G$$

(12)

R_ns is net radiation absorbed by soil surface and can be calculated using Beer’s law

$$R_{\mathit{ns}}{=R}_n \mathit{exp}(-C \mathrm{LAI})$$

(13)

where C is extinction coefficient of the crop for R_n, which was set to 0.7 in this study according to the Lambert–Beer Law [33].

Estimation of Resistances

The canopy resistance ($r_s^c$) computed from the inversed Equation (4) can be written as Equation (14).

$$r_s^c=\frac{C_C (\mathsf{\Delta }A\left(r_a^a{+r}_a^c\right){+(\rho }_a{c}_{p}D-{ \mathsf{\Delta }r}_a^c{A}_s\left)\right)}{T \gamma }-\frac{\mathsf{\Delta }\left(r_a^a{+r}_a^c\right)}{\gamma }-r_a^a-{ r}_a^c$$

(14)

One approach for modeling $r_s^c$, suggested by Katerji and Perrier [34], was established by a relationship between two ratios ($r_s^c/r_a$ and $r^{*}/r_a$) [32].

By defining a climatic resistance given as

$$r^{*}=\frac{\mathsf{\Delta }+\gamma }{\gamma }\frac{{\rho }_a{c}_p{(e}_s-{ e}_a)}{{\mathsf{\Delta }(R}_{n }- G)}$$

(15)

A climatic resistance is related to the isothermal resistance and first introduced by Monteith [35] and represents the surface resistance for equilibrium evaporation. The value mainly depends on climatic characteristics, although R_n and G are also influenced by the characteristics of the vegetative surface [36]. Katerji and Perrier [34] presented a linear link between $r_s^c/r_a$ and $r^{*}/r_a$. In this study, we adopted a non-linear functional relationship between $r_s^c/r_a$ and $r^{*}/r_a$ presented by Yan et al. [37]

$$\frac{r_s^c}{r_a}=a \times \frac{r^{*}}{r_a}+\mathit{b} \times \sqrt{\frac{r^{*}}{r_a}}+c$$

(16)

The aerodynamic resistance ($r_a$) can be calculated from Perrier [38,39]:

$$r_a=\frac{\ln \left. [\left(\mathit{x}-\mathit{d}\right) / \left({\mathit{h}}_c-\mathit{d}\right)\right] \ln \left. [\left(\mathit{x}-\mathit{d}\right){ / \mathit{z}}_0\right]}{{\mathit{uk}}^2}$$

(17)

where k is the Karman constant (=0.40); x is the reference height (=2 m); h_c is the mean crop height (m); d is the zero plane displacement (m); u is the wind speed at the reference height (m/s); and z₀ is the roughness length of the crop relative to momentum transfer (m). The z₀ and d are defined as 0.63 and 0.13 of the canopy heights, respectively [32].

The soil surface resistance ($r_s^s$) was calculated using the Ortega-Farias [22] power model, which is expressed as follows:

$$r_s^s=19 (\frac{{\mathit{\theta }}_{\mathit{sat}}}{{\mathit{\theta }}_{\mathit{swc}}}{)}^{3.5}$$

(18)

where $r_s^s$ is soil surface resistance (s m⁻¹); θ_swc is volumetric soil water content in the top layer of soil at 10 cm depth; θ_sat is saturated volumetric soil water content at 10 cm depth (0.40 m³ m⁻³).

In this study, the E was assumed to be negligible during spring shoots with low moisture content in 2015, so the measured T was replaced by the measured LE from March to May in 2015 to parameterize $r_s^c$ by solving Equation (14). The $r_s^c$ sub-model was integrated into S-W model for predicting LE and the model’s accuracy was validated by comparing the predicted and measured LE based on the data from 2016 to 2018.

2.4. Dual Crop Coefficient (D-K) Method

In the dual crop coefficient method, the ET_c is defined as the product of crop coefficient (K_c) and reference evapotranspiration (ET₀), and K_c is divided into soil evaporation coefficient (K_e) and basal crop coefficient (K_cb) [14].

$${\mathit{ET}}_c{=\mathit{K}}_c{ \times \mathit{ET}}_0$$

(19)

$${T=\mathit{K}}_{cb}{ \times \mathit{ET}}_0$$

(20)

$${E=\mathit{K}}_e{ \times \mathit{ET}}_0$$

(21)

2.4.1. Reference Evapotranspiration

The ASCE-EWRI [40] standardized the PM method for grass reference ET₀ with a condensed, simplified form from the original PM method:

$${\mathit{ET}}_0=\frac{0.408\mathsf{\Delta }\left(R_n-G\right)+\gamma \frac{\mathrm{900}}{T_a+\mathrm{273}}{\mathit{u}(e}_s-{ e}_a)}{\mathsf{\Delta }+\gamma (1+0.34\mathit{u})}$$

(22)

where ET₀ is reference evapotranspiration in mm d⁻¹ for daily time steps, R_n and G are net radiation at the crop surface and soil heat flux density at the soil surface in MJ m⁻² d⁻¹ for daily time steps, and T_a is the daily air temperature at 2.0 m height (°C). The meteorological data in Equation (22) was obtained in the Bowen ratio energy balance system installed in the center of the tea field.

2.4.2. Basal Crop Coefficient (K_cb) and Soil Evaporation Coefficient (K_e)

Based on the FAO-56 [14], K_cb can be expressed as:

$${\mathit{K}}_{\mathit{cb}}={\mathit{K}}_{c \mathit{min}}{+(\mathit{K}}_{\mathit{cb} \mathit{full}}-{ \mathit{K}}_{c \mathit{min} }) \times \left. [1 - \exp (-0.7 \times \mathrm{LAI})\right]$$

(23)

$${\mathit{K}}_{\mathit{cb} \mathit{full}}=\min \left(1.00+0.1\mathit{h}, 1.2\right)+\left. [0.04\left(\mathit{u} -2\right)-0.004\left({\mathit{RH}}_{\mathit{min}}-45\right)\right]{\left(\frac{\mathit{h}}{3}\right)}^{0.3}$$

(24)

where K_{c min} is the minimum value of basal crop coefficient for bare soil (= 0.1), K_{cb full} is the basal crop coefficient when crops have nearly full ground cover, h is the tea plant height.

K_e can be expressed as:

$${\mathit{K}}_e{=\mathit{K}}_r\left({\mathit{K}}_{c \mathit{max}}-{\mathit{K}}_{\mathit{cb}}\right) \le { \mathit{f}}_{\mathit{ew}}{\mathit{K}}_{c \mathit{max}}$$

(25)

$${\mathit{K}}_{c \mathit{max}}=\max \left(\left\{1.2+\left. [0.04\left(\mathit{u}-2\right)-0.004\left({\mathit{RH}}_{\mathit{min}}-45\right)\right]{\left(\frac{\mathit{h}}{3}\right)}^{0.3}\right\},\left\{{\mathit{K}}_{\mathit{cb}}+0.05\right\}\right)$$

(26)

where K_{c max} is the maximum value of K_c following rain or irrigation, f_ew is the fraction of the soil that is wetted (=0.5) for irrigation, K_r is the evaporation reduction coefficient dependent on the cumulative depth of water depleted from the soil surface, which is expressed as follows [14,41]:

$${\mathit{K}}_r=\frac{\mathit{TEW}-D_e}{\mathit{TEW}-\mathit{REW}}=\frac{\mathrm{1000}({\mathit{\theta }}_{\mathit{SWC}} -0.5{\mathit{\theta }}_{\mathit{wp}}{)\mathit{Z}}_e}{\mathit{TEW}-\mathit{REW}}$$

(27)

where TEW is total evaporable water (mm), which is the maximum depth of water that can be evaporated from the soil when the soil surface has been initially completely wetted; D_e is the cumulative depth of evaporation (depletion) from the soil surface layer (mm); REW is the readily evaporable water, which is the maximum depth of water that can be evaporated from the soil surface without restriction; θ_swc is the actual surface volumetric soil water content and θ_wp is the surface soil water content at wilting point (=0.12 m³ m⁻³ in this study); Z_e is the depth of the surface soil layer that is subject to drying by way of evaporation (=0.10 m in this study).

2.5. Evaluation of Models’ Performance

Statistical indices included a liner regression with 0 interception between observed and simulated values, root mean square error (RMSE), mean absolute error (MAE), index of agreement (d), and Bias were calculated for validating the accuracy of the models:

$$\mathit{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(E_i-O_i\right)}^2}{N}}$$

(28)

$$MAE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\left|E_i-O_i\right|$$

(29)

$$\mathit{d}= 1-\frac{{{\displaystyle \sum }}_{i=1}^N{\left(E_i-O_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(\left|E_i-\overline{O}\right|+\left|O_i-\overline{O}\right|\right)}^2}$$

(30)

$$\mathrm{Bias}=\frac{{{\displaystyle \sum }}_{\mathit{i}=1}^{\mathit{N}}{(E}_{\mathit{i}}-{\mathit{O}}_{\mathit{i}})}{{{\displaystyle \sum }}_{\mathit{i}=1}^{\mathit{N}}{\mathit{O}}_{\mathit{i}}}$$

(31)

where N is number of observations, E_i and O_i are estimated and observed values, and $\overline{O}$ is the mean observed value, respectively. A perfect model fit will have d = 1 and Bias = 0; positive values of Bias indicate model overestimation, and negative values of Bias indicates model underestimation [9,42].

3. Results

3.1. Interannual Variability of Climatic Factors at the Tea Field

The variations of annual climatic factors (i.e., R_n, T_a, VPD, and u) at the tea field during 2015–2018 are shown in Figure 1. The daytime R_n varied from 0 to 590 W m⁻² with an average value of 243, 230, 238, and 290 W m⁻² in 2015, 2016, 2017, and 2018, respectively. The daytime air temperature (T_a) changed from −7.1 to 36.1 °C with a mean of 17.9 °C during the four years, while the T_a of 2018 was the highest with a value of 19.6 °C and that of 2017 was the lowest with a value of 15.3 °C. Compared with T_a, VPD showed a similar interannual trend. It had a highest average value of 0.9 kPa in 2018, a minimum average value of 0.7 kPa in 2017 among the four years. The daytime wind speed (u) at height of 2 m varied from 0 to 7.5 m s⁻¹, with an average value of 2.7, 2.4, 2.5, and 2.2 m s⁻¹ in 2015, 2016, 2017, and 2018, respectively.

3.2. Parameterization of $r_s^c$

Diurnal variation of experimental values of $r_s^c/r_a$ and $r^{*}/r_a$ is shown in Figure 2, with a quadratic polynomial relationship being found between $r_s^c/r_a$ and ($r^{*}/r_a$)^0.5. The coefficients a, b, and c in Equation (16) in this study were 1.37, −0.18, and −0.17 (with R² = 0.84, RMSE = 0.30) for tea plants. Previous studies have been conducted to commonly express $r_s^c/r_a$ and $r^{*}/r_a$ as a linear function in alfalfa, grass, maize, and canola fields [34,43,44], but in this study we found a best-fit nonlinear relationship corresponded to a dependence of $r_s^c/r_a$ on the square root of $r^{*}/r_a$. The same function type was also reported in several studies, but the coefficients a, b, and c were different from the presented values (He et al. [45]: wheat field, 0.88, 0.82, and −1.95 in arid regions; Yan et al. [37]: maize field, 2.74, −5.9, and 7.04 in semiarid regions, and buckwheat field, 0.73, 1.25, and −0.28 in humid regions) as shown in Figure 1. As pointed out by Rana et al. [46], the coefficients of those relationships depended on the type of crop, its phenological state, and soil water status. Even for the same crop, the differences in coefficients among studies may exist due to different climatic regions [32,36].

3.3. Crop Coefficient (K_c), Basal Crop (K_cb), and Soil Evaporation Coefficient (K_e)

Figure 3a shows the variations of local measured K_c and the volumetric soil water content (θ_swc). Both K_c and θ_swc followed almost the same pattern along with DOY, for which the K_c values ranged from a minimum value of 0.43 in the 107 DOY to a maximum value of 1.44 in the 26 DOY. The measured K_c varied greatly and ranged from 0.43 to 1.44 with an average value of 0.90 during 1–150 DOY, while the K_c tended to be stable with an average value of 0.83 during 151–365 DOY. This may be due to a number of factors, for example; (1) the LAI dramatically changed for collecting tea during the spring and summer seasons; (2) the θ_swc varied more in the rainy season (Figure 3a), which usually occurred during the spring and summer seasons. Can the local K_c values estimate by the θ_swc for perennial plants like tea plants? Figure 3b shows the correlation between the local K_c and θ_swc in 2015, it was surprisingly found that there was a good consistency between the local K_c and θ_swc (R² = 0.92, Pearson correlation coefficient = 0.96 and p value = 0.000476 < 0.01).

The calculated basal crop coefficients (K_cb) of tea plants by Equation (23) based on the data from 2015 to 2018 are shown in Figure 4, together with the recommended K_cb by FAO-56 and the soil evaporation coefficient (K_e) calculated by Equation (25). The FAO-56 recommended K_cb for tea plants were 0.90 at initial stages, respectively [14], while the corresponding values of calculated K_cb were 0.89 at the initial stage (2015–2018). Figure 4 shows that the average calculated K_cb was close to, but lower than, the FAO-56 recommended values in 2015, while the average calculated K_cb was 0.95, 0.98, and 0.96 in 2016, 2017, and 2018, close to but more than the FAO-56 recommended values of 0.95. According to Allen et al. [14], K_cb are related to the local conditions, cultural practices, or crop varieties, but local values of K_cb should not be expected to deviate by more than 0.2 from the recommended values. For K_e, it was generally affected by surface soil water content and canopy coverage ratio [47]. The average value of K_e was 0.32 for 2015, higher than the value of 0.28 for 2016, 0.25 for 2017, and 0.23 for 2018, mainly due to the differences in canopy coverage ratio. From the Figure 4, values of K_e were close to 0 during DOY 240–260 in 2016 and DOY 275-285 in 2018, while values of θ_swc were close to 0 with an average value of 0.04 and 0.07 cm³ cm⁻³ at the corresponding dates.

3.4. The Performance of Two Methods in ET_c Simulation

Figure 5 showed the variations of estimated daily ET_c by the S-W and D-K method and the measured values by Bowen ratio energy balance method in the tea field in 2016-2018. Results showed that both approaches had good performance in the estimation of ET_c, except some overestimation of the D-K method at some unspecified date, e.g., 130 DOY of 2017 and 129 DOY of 2018 (Figure 5). From the regression analysis based on four years data (Figure 6), the average a and R² were 1.07 and 0.97 for the S-W model, while the corresponding values were 0.97 and 0.95 for the D-K method, indicating that both methods had high correlations with the measured ET_c.

Table 1. Statistical analysis of measured and estimated daily evapotranspiration (ET_c, mm d⁻¹) by the S-W and D-K method in 2016–2018.

Year	Model	ETc-Mimulated	ETc-Measured	a	R2	RMSE	MAE	d	Bias
2016	S-W	2.33	2.22	1.07	0.98	0.42	0.27	0.98	0.05
	D-K	2.17		1.00	0.95	0.59	0.43	0.96	−0.02
2017	S-W	2.16	2.01	1.08	0.97	0.51	0.34	0.97	0.07
	D-K	1.95		0.96	0.93	0.67	0.47	0.94	−0.03
2018	S-W	2.58	2.40	1.06	0.98	0.43	0.28	0.98	0.07
	D-K	2.31		0.95	0.96	0.58	0.39	0.97	−0.04
Average	S-W	2.36	2.21	1.07	0.97	0.45	0.30	0.98	0.06
	D-K	2.14		0.97	0.95	0.61	0.43	0.96	−0.03

Note: a is the coefficients of regression; R² is the coefficient of determination; RMSE is the root mean square error (mm d⁻¹); MAE is the mean absolute error (mm d⁻¹); d is the index of agreement.

showed the measured ET_c varied from 2.01 to 2.40 mm d⁻¹ with a mean of 2.21 mm d⁻¹ during the four years, while the measured ET_c of 2017 was the lowest with an average value of 2.01 mm d⁻¹. It can be explained that the most data (193–300 DOY) was lost in 2017. More statistical indices for evaluation of the accuracies of two approaches are shown in Table 1. The RMSE, MAE, and d were 0.45, 0.30 mm d⁻¹, and 0.98 for the S-W model, while the corresponding values were 0.61, 0.43 mm d⁻¹, and 0.96 for the D-K method. Moreover, the biases were 0.06 and −0.03 for the S-W model and D-K method, respectively. The statistical indexes showed both approaches could accurately estimate the daily ET_c in tea fields in southeast China, but the D-K method had a slightly poorer accuracy compared to the S-W model.

To further explore the differences of the S-W and D-K method in a separate estimation of T and E, the ratio of E to ET_c simulated based on two methods was compared (

Table 2. Models simulated values of tea plants transpiration (T) and soil evaporation (E) in 2016–2018.

Year	Model	T (mm d−1)	E (mm d−1)	E/ETc
2016	S-W	1.89	0.59	23.79%
	D-K	2.47	0.71	22.33%
2017	S-W	1.69	0.47	21.76%
	D-K	2.35	0.64	21.40%
2018	S-W	2.05	0.52	20.23%
	D-K	2.59	0.6	18.81%
Average	S-W	1.88	0.53	21.93%
	D-K	2.47	0.65	20.85%

). The simulated T were 1.89 mm d⁻¹ in 2016, 1.69 mm d⁻¹ in 2017, and 2.05 mm d⁻¹ in 2018 by S-W model, while the corresponding E were 0.59, 0.47, and 0.52 mm d⁻¹ for the three years, respectively. The E/ET_c simulated by the S-W and D-K method were 21.93 % and 20.85 %, respectively. The E/ET_c was declining year by year, and the highest value was 23.79 % in 2016. It can be explained that soil evaporation was related to the fraction of ground coverage when a smaller fraction of the soil surface was covered by the tea plants during the initial stage, which created a large wetted soil surface area that was exposed to sun radiation, and consequently higher soil evaporation. It is noteworthy that the difference between the E/ET_c simulated by the S-W and D-K method was small (around 1%).

4. Discussion

4.1. Parametrization of S-W Model

Previous studies have highlighted the fact that the canopy resistance is the most sensitive factor compared to the other resistances in the S-W model due to its comprehensive consideration of meteorological factors, crop characteristics, and soil moisture conditions [10,11,22,48,49]. Therefore, accurate determination of the canopy resistance is particularly important in accurate estimation of ET_c. Bao et al. [49] optimized a canopy resistance ($r_s^c$) sub-model of the S-W model and used a modified Ball–Berry model, including the empirical parameters estimated by the Monte Carlo algorithm for mobile dunes in China’s Horqin Sandy Land, with R² and RMSE between the measurements and simulations as 0.83, 0.32 mm d⁻¹, respectively. Liu et al. [50] improved the Jarvis-type $r_s^c$ sub-model by incorporating a term of effective leaf area index and a function to reflect the influence of the specific soil moisture, to improve the accuracy of the S-W model for rice ET_c in the Taihu lake region of China, with the R² and RMSE were 0.945, 0.934 mm d⁻¹, respectively. In this study, we optimized KP-type $r_s^c$ sub-model by a quadratic polynomial relationship based on the measured daily ET_c in a tea field, with the R² and RMSE at 0.97 and 0.45 mm d⁻¹, respectively. Compared to the previous studies, the S-W model parameterized by the present study achieved a higher accuracy. The main reason is that the KP-type $r_s^c$ sub-model in this study was optimized with a best-fit quadratic polynomial relationship, not a linear function, which is more adaptive to the tea field.

4.2. Prediction of Crop Coefficients

As for perennial plants, the average K_c values (0.86) found for the tea plants in this study were lower than that reported by Pinho Sousa et al. [51] for acai palm (=1.08), Meijide et al. [52] for palm oil (=1.03), and Flumignan et al. [53] for coffee trees (=1.28 for ET₀ < 3 mm d⁻¹, 0.98 for ET₀ > 3 mm d⁻¹). The possible reason for this is that the K_c in this study was determined using the dataset in 2015, while the tea plants were transplanted into the field in March 2014, which means the values of K_c were at the initial stage of the tea plants. Many researchers found that the K_c was significantly related to LAI [54,55,56]. For example, Wang et al. [56] presented a linear relationship between the K_c and LAI with a good result (R² ≥ 0.78, p < 0.01) in a grapevine ecosystem in the Nanhu Oasis of northwestern China. Singh Rawat et al. [55] reported that the K_c and LAI were in a strong relationship (second order polynomia, R² = 0.98) in a semi-arid environment. Guo et al. [54] presented a cubic polynomial function, which was the best for simulating the relationship between the K_c and LAI for spring maize in the arid region of Northwest China. In this study, due to the slight variations of daily LAI of tea plants and the difficulties in the measurement, the LAI was not considered to simulate the changes of the K_c. Instead, it was found that the θ_swc was significant correlated to the daily K_c (R² = 0.92, p < 0.01).

4.3. Comparison of Model Performance

The calculated statistical indices showed that the S-W and D-K method can estimate daily ET_c with good accuracy, except some big discrepancies of the D-K method some days. It can be explained that there are some irregular days for collecting tea, which may affect the daily K_c. Gharsallah et al. [26] compared the S-W model and D-K method in a surface irrigated maize agro-ecosystem in Northern Italy, and showed that the S-W model provided good results for the entire agricultural season, including initial growing stage, while the D-K method with generalized crop coefficients overestimated the ET_c, especially during the middle growth stage. Jiang et al. [27] compared the D-K and the S-W model in predicting daily ET_c and its components of maize, and reported that daily ET_c estimated by the S-W model was closer than by the D-K method to observed ET_c by the eddy covariance system at whole growing stage, while the D-K method overestimated daily ET_c at the initial and development stage in the arid region of northwest China. The principal reason for the poorer performance of the D-K method is that the tabulated K_cb adjusted with local data used for the simulations overestimated the actual, site-specific K_cb. Indeed, Gharsallah et al. [26] found that the K_c in mid-stage, even if adjusted to consider local conditions, was approximately 15% larger than the value derived from eddy covariance observations for maize grown in Northern Italy. However, Zhao et al. [19] showed that an opposite result to the D-K method performed better in in partitioning ET_c than the S-W model, while the S-W model significantly underestimated E, especially around wetting events of a vineyard in an arid region of northwest China. Gong et al. [11] also showed that the D-K method performed better compared to the S-W model in estimating daily tomato ET_c in greenhouses, while the S-W model overestimated ET_c by 17.9% at initial stage and underestimated ET_c by 16.6% at mid-stage in Henan province, China. The possible reasons for the S-W model having performed poorly are: (1) the $r_s^c$ in the S-W model was normally parameterized by the Jarvis approach or the modified Jarvis approach, the effect of water stress was not considered in the $r_s^c$ sub-model; (2) the $r_s^s$ in the S-W model was parameterized using the θ_swc at certain range (20–40%), while θ_swc deviated from the range especially before or just after the wetting events, the $r_s^s$ would not be parameterized by θ_swc accurately. Gong et al. [11] reported that the overestimation of S-W model at the initial stage was mainly due to the underestimation of $r_s^c$ under the soil water stress condition. Additionally, many studies have proved that the $r_s^c$ was the most sensitive variable compared to the other resistance in the S-W model [10,11,19,22].

4.4. Implications of the Modeling

As stated above, the establishment of the parameter sub-model played a crucial role in determining the performance of the S-W and D-K method [11]. The parameter sub-model should combine effects of local meteorological elements, crop growth conditions, and soil water status [19,22]. In this study, the parameter of $r_s^c$ and K_c sub-models were recalibrated by combining the influence of the meteorological elements and soil water status in the tea field. By integrating the $r_s^c$ and K_c sub-models into the S-W and D-K method, both methods had good performances in estimating ET_c of the tea field. However, how the effect of tea plants growth on the performance of the S-W and D-K methods can be integrated still needs to be further investigated.

In agricultural ecosystems, the fractions of E and T in ET_c are affected by soil and water management practices, and they have been used as indicators of crop water use efficiency (WUE) [57,58]. In this study, the range of E/ET_c in the tea field simulated by the S-W and D-K method was 22.33–23.79%, 21.40–21.76%, 18.81–20.23% in 2016, 2017, 2018, respectively. Some researchers have used in situ techniques (e.g., lysimeter method, eddy covariance, heat pulse sensors, and sap flow methods) to obtain the fraction of E in ET_c, and Wang et al. [58] obtained an E/ET_c value of 19% using sap flow gauges (to measure T) and the weighing lysimeter (to measure ET_c) in a maize field; Sauer et al. [59] observed a E/ET_c range of 8 %–12 % in a narrow-row soybean field using the eddy covariance system (to measure ET_c) and sap flow stem gauges (to measure T); Wagle et al. [60] collected high frequency (=10 Hz) time series eddy covariance system observations over a rainfed alfalfa, and reported that the average E/ET_c was approximately 18–23% in central Oklahoma, USA. Compared to the reported values of E/ET_c from previous literatures, there is still space to minimize E loses in the tea fields. Hence, based on the results simulated by the constructed S-W and D-K methods in this study, reducing a proper amount of irrigation water to decrease the soil evaporation at the growing season of tea plants is an advisable way for improving water use efficiency, especially at initial stages in 2016.

5. Conclusions

In this study, we evaluated the S-W and D-K method in estimating daily evapotranspiration (ET_c) and its components of a tea field, using the measurements of ET_c by Bowen ratio energy balance system for four years (2015–2018) in southeast China. The canopy resistance ($r_s^c$) in the S-W model was parameterized by a climatic resistance $r^{*}$, with a quadratic polynomial relationship between $r_s^c/r_a$ and square root of $r^{*}/r_a$. The soil surface resistance ($r_s^s$) in the S-W model was parameterized by the soil water content (θ_swc) at the top layer. The measured crop coefficient (K_c) was defined as the ratio of measured ET_c to ET₀, basal crop coefficient (K_cb), and evaporation coefficient (K_e) in the D-K method was determined based on the measured meteorological data, LAI, and θ_swc data in the tea field.

Both the S-W and D-K method had good performances in estimating ET_c of the tea field, with an average RMSE and R² of 0.53 mm d⁻¹ and 0.96, while the S-W model performed slightly better than the D-K method. The measured and simulated average daily ET_c of the tea field was 2.21 mm d⁻¹ (Bowen ratio energy balance), 2.36 mm d⁻¹, and 2.14 mm d⁻¹ (S-W model and D-K method, respectively). The average ratios of E to ET_c simulated by the S-W and D-K method were 21.93 % and 20.85 %, respectively. The above results indicated that the E should be estimated independently for tea fields to improve the accuracy of ET_c modeling.