Parameterization of Four Models to Estimate Crop Evapotranspiration in a Solar Greenhouse

Working to simplify mechanistic models on the basis of reliability for estimating crop evapotranspiration (ET) in a greenhouse is still worthwhile for horticulturists. In this study, four ET models (Penman–Monteith, Priestley–Taylor, and Shuttleworth–Wallace models, and the Crop coefficient method) were parameterized after taking the restriction effect of resistance parameters in these models on ET into account, named as PA-PM, PA-PT, PA-CC, and PA-SW, respectively. The performance of these four parameterized models was compared at different growth stages, as well as the entire growing season. Tomatoes that were ET-grown in a solar greenhouse without a heating device were measured using weighting lysimeters during 2016–2017 and 2019–2021, in which data from 2016 were used to adjust the model parameters, and data from the other four study years were used to examine the model performance. The results indicated that the PA-PT and PA-CC models have a better performance in estimating tomato ET at four growth stages, while the PA-PM and PA-SW performed well only at the development and middle stages. Compared to the ET that was measured with the weighting lysimeters, the ET that was predicted using the PA-PM model was 27.0% lower at the initial stage, and 8.7% higher at the late stage; the ET that was computed using the PA-SW model was 19.5% and 13.6% higher at the initial and late stages, respectively. The PA-PT model yielded the lowest root mean square error and the highest index of agreement against the other models over the entire growing season, indicating that the PA-PT model is the best recommended model for estimating tomato ET in a solar greenhouse.


Introduction
The rapid development of facility agriculture around the world provides high-quality vegetables, which greatly meet the demand for human beings for a quality life.In China, the greenhouse industry has created great economic value for local farmers while improving planting efficiency [1,2].Greenhouse production, unlike field cultivation, requires rigorous and correct irrigation scheduling and climate conditions to maintain high production levels and crop quality [3].Improving the reliability and precision of the irrigation amount and frequency for greenhouse crops over the entire growing season is a challenge for horticulturists and agriculturalists.Evapotranspiration (ET), as the outgoing vapor flux at the boundary between the soil, plant, and atmosphere, is a key variable for optimizing irrigation scheduling, and has been supported by many studies [3][4][5][6][7][8].
Measuring ET, or simulating it numerically, is a common way for mastering its pattern changes.In greenhouse climates with semi-enclosed structures, instruments based on field measurements require relatively complex physical principles and techniques, such as soil water balance [9][10][11][12], lysimeters [13][14][15], Bowen ratio [16], and sap flow plus microlysimeters [17,18].In addition, instrumental measurement is expensive, time-consuming, and complex; further, it is difficult to avoid errors when using instrumental measurement [5,16,19,20].Given the limitation of field measurements, therefore, obtaining ET by establishing a reliable mathematical model provides an effective alternative method that many studies give priority to [7].
At present, models frequently used in greenhouse for computing crop ET include the Penman-Monteith (PM), Priestley-Taylor (PT), and Shuttleworth-Wallace (SW) models, and the Crop coefficient method (CC).The single-layer standard PM model has been widely used for estimating ET in a variety of greenhouse types.The canopy cover is schematized as a single "big leaf" in the PM model and placed at a certain height within the canopy.The successful use of the PM model relies upon two key parameters, e.g., surface resistance (r s ) and aerodynamic resistance (r a ).The r s is a bulk resistance parameter that acts in series between the canopy surface and soil surface, while the r a can be defined as a function of the wind speed and canopy height [11].However, the PM "single-layer" model faces the problem of application with simple and easily accessible parameters, especially for r s and r a .In the past few years, some empirical and semi-empirical models have been developed to calculate the r s by introducing terms or functions to make these models more suitable for special environmental conditions and improve the accuracy of ET estimation [21][22][23][24][25][26][27].Among them, the Jarvis-Stewart model has been widely applied due to its physical representation of the environmental conditions controlling ET [4,13,14,28].Another approach to estimating the r a parameter in greenhouses is the heat transfer coefficient approach [29], which can be used due to the low wind speed (<0.5 m/s at most of the days).However, calculating the heat transfer coefficient requires the accurate real-time acquisition of the canopy temperature, which has high requirements on measuring instruments, giving rise to uncertainties in the estimation of r a , especially under a sparse canopy space.In general, the complex computation process of r a and r s has slowed down the routine application of the PM model by farmers.Allen et al. [11] suggested that the canopy resistance can be considered as a constant, independent of climatic conditions; this simplifies the calculation procedure of the PM equation, but has only been verified on irrigated grass.Fernández et al. [30] introduced r a as a fixed value of 295 s/m into the PM model for computing crop ET in greenhouses.These results provide guidance for the parameterization of the PM model in greenhouses.
The PT model is a simplification of the Penman equation but it neglects the influence of vapor pressure deficit on crop ET [31].The PT model assumes that ET is the product of the equilibrium evaporation and PT coefficient (α PT ).The former depends only on solar radiation and air temperature, while the latter was an empirical parameter, influenced by the complexity of the natural environment [32,33].Therefore, using the PT model to calculate ET in a greenhouse with low wind speed is convenient.Previous research reported that the successful use of the PT model depends on the accurate determination of α PT [34][35][36][37].A well-known α PT was determined to be 1.26 by using daily measurements of ET on many vegetated areas, while subsequent studies and theories pointed to a nearconstant α PT ∈ [1.2, 1.3] [31,34,38,39].Conversely, other studies indicated that the value of the α PT had a large fluctuation over the whole growth period, and was mainly controlled by the mulching styles, leaf area index, air temperature, soil water content, relative humidity, etc. [35,[40][41][42][43].The improved α PT equation, which was developed by considering the above factors, generally has a higher simulation accuracy, but the simple PT model becomes more complicated due to the many parameters involved.Therefore, developing a simple α PT equation without losing the precision of the PT model in a greenhouse climate is still worthwhile work.
The crop coefficient method was a common method for computing ET since the publication of the FAO56 guidelines in 1998.In this method, crop ET was generated using the product of the specific crop coefficient (K c ) and reference evapotranspiration (ET o ).Here, ET o is the potential evaporating power of the atmosphere at a uniform underlying surface and can be calculated from weather data in reference crops, such as clipped grass or alfalfa [11,19].Our previous study reported that taking the value of r a as a constant of Plants 2024, 13, 1579 3 of 22 308 s/m in a solar greenhouse was very close to the measurement [17].The K c synthesizes physical and physiological differences between crop-specific surfaces and reference surfaces, and varies with management factors (e.g., crop type, climate, irrigation methods, etc.); hence, the K c plays an important role in practical applications [44].The segmented curves of the K c were described in the FAO 56 guidelines, which is an approximate representation of the change in crop vegetation from planting to harvesting.At the initial and middle stages, the K c can be represented by two horizontal lines, since the variation in the crop growth in these two stages with time is small.The K c at the development and late stages is characterized by rising and falling line segments, respectively, which is closely related to the canopy cover of crop vegetables.The K c values for most crops were recommended in the FAO56 guidelines, which provided a convenient method for users [19].Using the same K c values in different climates to estimate ET, however, was subject to significant uncertainties, as the K c varies greatly with the environment [45,46].As a result, the K c needs to be amended based on local standards to guide irrigation scheduling more accurately.Therefore, investigating the variation patterns in the K c and modifying them according to special climate conditions, are important for developing ET-based irrigation scheduling.
The dual-layer SW model consists of two standard PM equations used to calculate soil evaporation and plant transpiration, respectively, and weighted by a set of coefficients representing the combination of soil and canopy resistance [47].In the SW model, two surface resistances are employed to regulate the transfer of heat and mass from soil (r s s ) and plants (r s c ), and three aerodynamic resistances (r a a , r a c , r a s ) are used to regulate fluxes to the atmosphere [24].Although the SW model was widely applied to calculate ET for most crop types, it has a large number of parameters that are hard to parameterize, resulting in more cumbersome calculation procedures [48].The performance of the SW model has been examined in solar greenhouses for tomatoes [13] and in Venlo-type greenhouses for cucumbers [49]; also, several resistance parameters of the SW model were analyzed over the whole growth stage.In addition, Zhao et al. [50] conducted a sensitivity analysis to uncertainties in input parameters and resistance for the SW model, and found that sensitivity of compute ET in the r a a , r a c , and r a s are unapparent, which provides a valuable reference for the parameterization of the SW model in greenhouses.
In conclusion, an effective approach for significantly improving the accuracy of ET estimation is to recalibrate the parameters in the model, which has been widely raised by many scholars.However, the greater number of parameters that need to be calibrated, the greater the risk of overfitting [48].Therefore, parameterizing the PM, PT, CC, and SW models in solar greenhouses according to the characteristic of a semi-enclosed climate is still meaningful work, so that the precision is acceptable and model structure is simplified.The purpose of our study is (1) to parameterize the Penman-Monteith, Priestley-Taylor, and Shuttleworth-Wallace models, and the Crop coefficient method, for estimating crop ET in solar greenhouses in North China, and (2) to evaluate the performance of these models at different growth stages and the entire growth period.

Penman-Monteith Model
The Penman-Monteith (PM) model which is used for calculating ET describes the canopy cover as a "big leaf" placed at a reference height above the crop; the structure diagrams of the PM model can be found in Allen et al. [11].The formula is given as follows: In the PM model, physically based parameters of aerodynamic resistance (r a ) and surface resistance (r s ) are two key factors that affect the accuracy of the PM model.Three approaches for calculating r a were compared and evaluated in our previous study [17]; in this work, the r a variations were analyzed over the whole growth stage.To our surprise, setting the r a as a constant of 308 s/m is acceptable to calculate the ET.So, the measured r s can be calculated by inverting Equation (1) with the measured ET from three weighing lysimeters, and by taking r a as 308 s/m in this study.To simplify the calculation method of the r s , the statistical analysis revealed that a power function can be fitted for the relationship between the r s and LAI (see Figure 1).this work, the ra variations were analyzed over the whole growth stage.To our surprise, setting the ra as a constant of 308 s/m is acceptable to calculate the ET.So, the measured rs can be calculated by inverting Equation (1) with the measured ET from three weighing lysimeters, and by taking ra as 308 s/m in this study.To simplify the calculation method of the rs, the statistical analysis revealed that a power function can be fitted for the relationship between the rs and LAI (see Figure 1).

Priestley-Taylor Model
The Priestley-Taylor (PT) model, based on equilibrium evaporation, is a simplified form of the Penman model [31]; the formula is expressed as follows: The successful use of the PT model lies primarily in determining the PT coefficient (αPT) [34,35], although some studies have shown that the coefficient αPT could be used as a constant (1.26) on some vegetated surfaces [33,36].Conversely, many studies believe that the αPT was variable over the whole growing stage [34,42].In this study, the coefficient αPT is computed by re-arranging the PT model, where the data of λET is from the weighing lysimeters.Thereafter, the relationship between the αPT and the days after transplanting (DAT) was established by using the data from the 2016 study years (see Figure 2).

Priestley-Taylor Model
The Priestley-Taylor (PT) model, based on equilibrium evaporation, is a simplified form of the Penman model [31]; the formula is expressed as follows: The successful use of the PT model lies primarily in determining the PT coefficient (α PT ) [34,35], although some studies have shown that the coefficient α PT could be used as a constant (1.26) on some vegetated surfaces [33,36].Conversely, many studies believe that the α PT was variable over the whole growing stage [34,42].In this study, the coefficient α PT is computed by re-arranging the PT model, where the data of λET is from the weighing lysimeters.Thereafter, the relationship between the α PT and the days after transplanting (DAT) was established by using the data from the 2016 study years (see Figure 2).Here α values were calculated by inverting the Priestley-Taylor model by using the 2016 ET data.The black line is the fitting regression line.

Crop Coefficient Method
In the FAO56 manual, the crop coefficients (Kc) can be computed by the ratio of the measured ET to the reference evapotranspiration (ETo).

Crop Coefficient Method
In the FAO56 manual, the crop coefficients (K c ) can be computed by the ratio of the measured ET to the reference evapotranspiration (ET o ).
The ET o was calculated by using the FAO 56 Penman-Monteith formula with a fixed r a of 308 s/m [17], and is given as follows: Four growth stages (initial, development, middle, and late) were divided according to FAO Irrigation and Drainage Paper No. 56; the initial stage runs from the planting date to approximately a 10% ground cover, the development stage runs from a 10% ground cover to an effective full cover, the middle stage runs from the effective full cover to the start of maturity, and the late stage runs from the start of maturity to the harvest [11].The observed length of the growth seasons for the greenhouse tomato is given in Table 1.The K c is defined as K cini , K cmid, and K cend at the initial, middle, and late stages, respectively.Furthermore, the values of K cini , K cmid, and K cend in the 2016 study year were obtained by using Equation (3) as 0.49, 1.03, and 0.86 (Figure 3 shows each respective segment of the linear relationship).

Shuttleworth-Wallace Model
In the Shuttleworth-Wallace (SW) model, soil evaporation and plant transpiration as one-dimensional models were combined.The surface resistance involved occurs at the

Shuttleworth-Wallace Model
In the Shuttleworth-Wallace (SW) model, soil evaporation and plant transpiration as one-dimensional models were combined.The surface resistance involved occurs at the interfaces between the soil and canopy and is further regulated by the aerodynamic resistance, which governs the transfer of water and heat between the soil and plants, plants and canopy, as well as soil and canopy; the structure diagrams of the SW model can be found in [47].
The meaning of the letters in the above equations can be found in the list of symbols and in Shutleworth and Wallace's paper [47].A and A s are given as follows: The extinction coefficient of light attenuation (C) was monitored by using two radiation sensors; one was placed 30 cm above the plant canopy and the other was placed on the soil surface.Based on the data from the 2016 study year, we found that the maximum value of C was 0.85 when the LAI reached its peak, and a value of 0 was taken for bare soil [13].Thereafter, linear equations between C and LAI were used to compute C on each day.
In the SW model, five resistances should be precisely calculated, including two aerodynamic resistances (r a a and r a s ), two surface resistance (r s c and r s s ), and one boundary layer resistance (r a c ); the detailed calculations can be found elsewhere.The SW model improvement in this study is mainly developed to consider how to deal with these resistance parameters.Zhao et al. [50] indicated that the fluctuation of r a a and r a s was particularly insensitive to the SW model; our study attempts to define them as constants at different growth stages.The values of r a a and r a s in 2016 can be taken as 1124 and 576, 670 and 1005, 390 and 900, and 430 and 930 s/m at the initial, development, middle, and late stages, respectively (Figure 4a).So, these values were also used in other study years.To simplify the method of obtaining r s c and r s s , we tested the relationship between r s c , r s s , and environmental factors in the 2016 study year, and found that a power function can be effectively fitted for the relationship between r s c and LAI (y = 228.18x−0.898 , R 2 = 0.98) (Figure 4b), while an exponential function can be effectively fitted for the relationship between r s s and θ s (y = 1225.6e−9.677x , R 2 = 1.0) (Figure 4d).These two relationships are also applied to other study years.The r a c was also insensitive to the SW model; however, we found that r a c and LAI have a good power function relationship (y = 112.06x−1.011 , R 2 = 0.96) (Figure 4c).These parameters required by the Penman-Monteith, Priestley-Taylor, and Shuttleworth -Wallace models, and the Crop coefficient method, after parameterization are shown in Table 2.

Parameterized Models Meteorological-Based Parameters Crop Parameters Soil Parameters
Penman The meaning of the letters in Table 2 can be found in the list of symbols.

Evaluation of Model Performance
Ten indexes were utilized to assess the performance of these improved models, e.g., the regression coefficient (b 0 ), determination coefficient (R 2 ), mean absolute error (MAE), root mean square error (RSME), rank sum ratio (RSR), average absolute error (AAE), average relative error (ARE), percent bias (PBIAS), modeling efficiency (EF), and index of agreement (d IA ).Their formulas are given as follows [13,51]: Plants 2024, 13, 1579 where O i and P i (i = 1, 2, . .., n) are the observed and predicted values, and O and P represent the mean values, respectively.The MAE can be used to reflect the actual prediction error.
The RMSE is the variance in the estimated error.The RSR can be used to compare the mean, median, variance, and other indicators of two or more samples.The AAE can express the mean size of the estimation error.The size of the error in relative terms can be expressed by the ARE.The PBIAS can express the average tendency of the predicted values to their corresponding observations.The ratio of the mean square error to the variance in the observed data can be defined as the EF.The d IA can be used to represent the largest relative value that can occur from each observation-model simulation pair of values.In these indexes, the b 0 , R 2 , EF, and d IA are close to 1, while the MAE, RMSE, RSR, AAE, and ARE are close to 0, indicating that the improved model performs better.All of the figures, tables, and data processing in this paper are completed by using the software of Microsoft Excel 2016.

Experiment Site and Description
The experiment was conducted from March to July 2016-2017 and 2019-2021 at the Xinxiang Comprehensive Experimental Base of the Chinese Academy of Agricultural Sciences located in Henan Province (35 • 86 ′ N, 113 • 68 ′ E, altitude 78.7 m).This site had a temperate continental climate with an average annual temperature of 14.2 • C, of which July was the hottest month, with a mean temperature of 27.1 • C, while January was the coldest month, with a mean temperature of 0.7 • C, and an average annual humidity of 68%.The average annual rainfall was 548 mm and the annual evaporation was 1908 mm.The most precipitation occurred from June to September, with a total rainfall of 409.7 mm, accounting for 72% of the annual precipitation.A local steel-framed solar greenhouse was employed to carry out this study, which was topped with a 0.02 mm thick polyethylene film and fitted with vents (20 m in length × 0.5 m in width), with brick walls on the north, east, and west sides.The greenhouse was located in the east-west direction and covered an area of 510 m 2 .Tomato seedlings with varieties of Solanum lycopersicum and c.v. Jinding were transplanted into the greenhouse in early March, with a planting density of 5.7 plants per square meter, and a plant spacing of 30 cm; the division of the tomato growth stages can be found in Table 1.Irrigation events were controlled by referencing the 20 cm standard evaporation pan; irrigation was performed when the accumulated pan evaporation reached 20 ± 2 mm.In our study, drip irrigation was employed with a discharge rate of 1.1 L/h, and the irrigation amount was controlled by water meters.

Measurements
Three categories of content were measured, including meteorological factors, crop growth parameters, and soil characteristics.
Meteorological factors in the central area of the greenhouse, 2 m above the ground, were monitored by using an automatic weather station, including air temperature (T a ), relative humidity (RH), net radiation (R n ), canopy temperature (T c ), soil heat flux (G), and wind speed (u 2 ).A CS215 (Campbell Scientifc Inc., Logan, UT, USA) combined air temperature and relative humidity sensor was used for measuring T a and RH.R n was monitored by using an NR LITE2 net radiometer (Kipp & Zonen, Delft, The Netherlands), with a sensitivity of 10 µV/(W•m 2 ).The factor of T c was measured by using an SI-111 infrared radiation pyrometer (Campbell Scientifc Inc., Logan, UT, USA), installed 0.3 m above the canopy surface.Below 5 cm of the soil surface, two HFP01 soil heat flux plates (Hukseflux, Delft, The Netherlands) were employed to measure and adjust the G [40].In addition, a high-precision air velocity meter (Wind Sonic, Gill, Nottingham, UK) with a precision of ±0.02 m/s was used to measure u 2 inside the greenhouse.All meteorological data were collected in the data logger (Campbell Scientifc Inc., Logan, UT, USA) every 30 min.
The plant height and leaf area were measured manually at 7 d intervals; ten replications were set in this study.The missing data were obtained by using the 'piecewise cubic Hermite interpolation polynomial" method in MATLAB 2020 software.The leaf area index (LAI) can be defined as the ratio of the plant leaf area to the unit soil surface area, which was determined by summing the rectangular area (leaf length × maximum width) of each leaf multiplied by a reduction factor of 0.64 [6].
The soil water content on the surface (0-5 cm) was continuously monitored using five ECH2O sensors (5TE, Decagon Devices, Inc., Pullman, WA, USA), which were distributed 15 cm away from the drippers along the line.The data were collected at 30 min intervals using an EM50 data logger (Decagon Devices, Inc., Pullman, WA, USA).
Crop evapotranspiration (ET) was monitored by using three weighing lysimeters, located in the middle of the greenhouse, with the specification of 1.0 m in length × 1.0 m in width × 1.2 m in depth.The lysimeter was filled from top to bottom with field soil for the first 80 cm depth, followed by a 10 cm thick fine aggregate and a 10 cm thick coarse aggregate to facilitate percolation.Six tomato plants were transplanted in each lysimeter with the same planting pattern as the field; bamboo poles were used to support plants to avoid lodging.The data of ET were reported using a micro-computer at 1.0 h intervals.

Variations of Daytime Meteorological Factors, Leaf Area Index, and Soil Moisture in Five Study Years
The variation patterns of meteorological factors (R n , VPD, T a ), reference evapotranspiration (ET o ), LAI, and surface soil water content (θ s ) during the whole growth period in 2016-2017 and 2019-2021 were investigated primarily.Figure 5a shows the variation in the R n in the daytime during the five study years.The daytime R n ranged from 0.1 to 16.2 MJ/(m 2 •d), with an average of 8.2 MJ/(m 2 •d), and the maximum R n often occurred on sunny days in June or July, while minimum values were frequently observed in overcast or rain days in every year.Figure 5b depicts that the daytime VPD varied from 0.1 to 4.0 kPa with an average of 1.6 kPa for the five study years.The variations in VPD in the greenhouse seem to be irregular, which is possibly due to the frequent occurrence of a high T a and high RH inside the greenhouse.The variations in the T a were consistent with the R n , which increases as the growth period progresses.The T a during the five years was between 11.5 and 37.0 • C with the average of daytime T a at seedling, development, middle, and late stages of 20.9, 23.8, 28.8, and 31.7 • C, respectively (Figure 5c).
The variations in the ET o showed a very similar trend with the R n , which fluctuated from 0.1 to 5.9 mm with a mean of 3.1 kPa.Because R n was the main control factor of the ET o .The average of the ET o at the seedling, development, middle, and late stages was 2.3, 2.9, 3.6, and 3.6 mm, respectively (Figure 5d).The tomato LAI presented a regular parabolic curve, and the LAI at the initial stage was relatively low, generally below 0.5 m 2 /m 2 .The maximum value of the LAI for every year differed remarkably (Figure 5e), with a maximum value of 3.2, 3.8, 2.8, 3.2, and 3.2 m 2 /m 2 in the five years.The difference among the five years was mainly caused by crop growth and development.The variations in the θ s are displayed in Figure 5f; it varied from 0.17 to 0.35.Before the development stage, the soil moisture is relatively low, which is related to local cultivation practices.Generally, 20 mm of irrigation water was injected after transplantation; then, the tomato seedlings are squatted until they enter the development stage.Thereafter, irrigation events were performed again.spiration (ETo), LAI, and surface soil water content (θs) during the whole growth period in 2016-2017 and 2019-2021 were investigated primarily.Figure 5a shows the variation in the Rn in the daytime during the five study years.The daytime Rn ranged from 0.1 to 16.2 MJ/(m 2 •d), with an average of 8.2 MJ/(m 2 •d), and the maximum Rn often occurred on sunny days in June or July, while minimum values were frequently observed in overcast or rain days in every year.Figure 5b depicts that the daytime VPD varied from 0.1 to 4.0 kPa with an average of 1.6 kPa for the five study years.The variations in VPD in the greenhouse seem to be irregular, which is possibly due to the frequent occurrence of a high Ta and high RH inside the greenhouse.The variations in the Ta were consistent with the Rn, which increases as the growth period progresses.The Ta during the five years was between 11.5 and 37.0 °C with the average of daytime Ta at seedling, development, middle, and late stages of 20.9, 23.8, 28.8, and 31.7 °C, respectively (Figure 5c).
The variations in the ETo showed a very similar trend with the Rn, which fluctuated from 0.1 to 5.9 mm with a mean of 3.1 kPa.Because Rn was the main control factor of the ETo.The average of the ETo at the seedling, development, middle, and late stages was 2.3, 2.9, 3.6, and 3.6 mm, respectively (Figure 5d).The tomato LAI presented a regular parabolic curve, and the LAI at the initial stage was relatively low, generally below 0.5 m 2 /m 2 .The maximum value of the LAI for every year differed remarkably (Figure 5e), with a maximum value of 3.2, 3.8, 2.8, 3.2, and 3.2 m 2 /m 2 in the five years.The difference among the five years was mainly caused by crop growth and development.The variations in the θs are displayed in Figure 5f; it varied from 0.17 to 0.35.Before the development stage, the soil moisture is relatively low, which is related to local cultivation practices.Generally, 20 mm of irrigation water was injected after transplantation; then, the tomato seedlings are squatted until they enter the development stage.Thereafter, irrigation events were performed again.

Performance of Parameterized Models for Simulating Daily ET at Four Growth Stages
The performance of the parameterized Penman-Monteith (PA-PM), Priestley-Taylor (PA-PT), Shuttleworth-Wallace (PA-SW) models, and the Crop coefficient method (PA-CC) for estimating greenhouse tomato ET was evaluated.Tables 3-6 and Figure 6 show the comparison of daytime ET estimated using the four parameterized models and measured using the weighting lysimeters at the initial, development, middle, and late stages in 2017 and 2019-2021.

Performance of Parameterized Models for Simulating Daily ET at Four Growth Stages
The performance of the parameterized Penman-Monteith (PA-PM), Priestley-Taylor (PA-PT), Shuttleworth-Wallace (PA-SW) models, and the Crop coefficient method (PA-CC) for estimating greenhouse tomato ET was evaluated.Tables 3-6 and Figure 6 show the comparison of daytime ET estimated using the four parameterized models and measured using the weighting lysimeters at the initial, development, middle, and late stages in 2017 and 2019-2021.6).The PA-PT and PA-CC models show a better performance for estimating the daily ET at the initial stage; the b 0 was close to 1.0, the MAE ranged from 0.21 to 0.29 mm/d, and the RMSE ranged from 0.28 to 0.38 mm/d in four years (Tables 4 and 5).

years
1  At the development and middle stages, four parameterized models averagely produced the error of daily ET, ranging from 1.1% to 8.2% (Figure 6); these models yielded a b 0 close to 1.0, and most of the R 2 values were above 0.8, with the MAE changing from 0.39 to 0.57 mm/d, and the RMSE ranging from 0.51 to 0.71 mm/d.Furthermore, most of the EF values were near or above 0.8, and the d IA values were higher than 0.9.The statistical analysis also shows that the RSR values were lower than 0.5, and the AAE values were below 0.6 mm/d; the ARE varied from 16.5 to 31.5, and the PBIAS was within a positive and negative 10 value (Tables 3-6).
At the late stage, the PA-PM and PA-SW models averagely overestimated the daily ET by 8.7% and 13.6%, respectively, while the daily ET estimated by PA-PT and PA-CC models was very close to the measurements in the four study years (Figure 6).The statistical analysis shows that the MAE and RMSE yielded by the PA-PT and PA-CC models were lower than 0.5 and 0.7 mm/d, whereas the MAE and RMSE produced by the PA-PM and PA-SW models were higher than 0.55 and 0.75 mm/d in 2017 and 2019-2021 (Tables 3-6).
In addition, compared with the PA-PM and PA-SW models, the PA-PT and PA-CC models have a lower AAE, ARE, and PBIAS, and a higher EF and d IA .
In summary, the precision of the PA-PM and PA-SW models in simulating the daily greenhouse tomato ET at the initial and late stages was unsatisfactory, while the PA-PT and PA-CC models showed a better performance.At the development and middle stages, the daily greenhouse tomato ET calculated by the four parameterized models agreed well with the measured ET from the weighting lysimeters.Among them, the PA-PT model was preferred in the greenhouse for estimating daily tomato ET due to its brief structure and fewer parameters.

Performance of Parameterized Models for Simulating Daily ET over the Whole Growth Stage
To evaluate the performance of four parameterized models for simulating daily ET systematically, 4 year weighting lysimeter data from 2017 and 2019-2021 were used to test the model precision over the whole growth stage (Table 7 and Figure 7).To our surprise, minor errors were found in all four parameterized models; the PA-PM and PA-SW models slightly overestimated the daily ET by 5.4% and 3.0%, respectively, while the daily ET was slightly underestimated by the PA-PT and PA-CC models, with values of 0.7% and 2.7% in the four study years (Figure 7).The statistical analysis shows that all models yielded the R 2 higher than 0.8; particularly, the R 2 of the PA-PT model was higher than 0.9 (Table 7).The MAE and RMSE produced from the PA-SW model were higher than other models, while the PA-PT model has the lowest MAE and RMSE, which were 0.38 and 0.51 mm/d (Table 7), indicating that the PA-PT model has the better performance among the four parameterized models.The indexes of the RSR, AAE, and PBIAS have a similar trend to the MAE and RMSE, indicating that the PA-PT model has a lower estimation error compared to other models.In addition, a smaller ARE was found in the PA-CC model (Table 7), indicating that the estimation error amplitude of the PA-CC model was small over the whole growth period, that is, the stability of the PA-CC model in greenhouses was better than the PA-PM and PA-SW models.

Evaluation of Four Parameterized Models to Compute ET in Greenhouse
In our study, four models were parameterized to compute the ET of greenhousegrown tomatoes, including the Penman-Monteith, Priestley-Taylor, and Shuttleworth-Wallace models, and the Crop coefficient method.The results indicated that the PA-PM and PA-SW models severely under-and overestimated the tomato ET at the initial stage, respectively, while they slightly overestimated the tomato ET at the late stage.The performance of the PA-PM and PA-SW models, however, was better at the development and middle stages (Figure 6 and Tables 3 and 6).The PA-PT and PA-CC models performed well at the four growth stages.The performance of the PA-PM model at the initial stage was somewhat different from the results obtained by Palmer et al. [52]; these authors pointed out that the PM model tended to over-predict the ET when the LAI becomes small.Its imperfection is not surprising; the surface resistance of the PA-PM model in our study was fitted using the LAI as a power function, resulting in the fitted rs increasing rapidly with the decrease in the LAI (LAI < 0.5 m 2 /m 2 at the initial stage).The larger the rs, the smaller the ET simulated by the PM model, which has been confirmed in many studies [4,14,17,53,54].Conversely, the rs was underestimated at the late stage due to a higher LAI, The consistency index evaluation of the four parameterized models shows that the EF and d IA yielded from PA-PT model were the highest, indicating that the mean square error of the observed data was very close to the variance, and a small gap in the maximum relative value of the simulation results was obtained by using the PA-PT model.Although the consistency index of the PA-PM, PA-CC, and PA-SW models were slightly lower compared with the PA-PT model, their performance was still satisfactory, with the EF and d IA above 0.8 and 0.95 (Table 7), respectively.In a word, the PA-PT model was recommended due to a high simulation accuracy and easier access to model parameters.

Evaluation of Four Parameterized Models to Compute ET in Greenhouse
In our study, four models were parameterized to compute the ET of greenhouse-grown tomatoes, including the Penman-Monteith, Priestley-Taylor, and Shuttleworth-Wallace models, and the Crop coefficient method.The results indicated that the PA-PM and PA-SW models severely under-and overestimated the tomato ET at the initial stage, respectively, while they slightly overestimated the tomato ET at the late stage.The performance of the PA-PM and PA-SW models, however, was better at the development and middle stages (Figure 6 and Tables 3 and 6).The PA-PT and PA-CC models performed well at the four growth stages.The performance of the PA-PM model at the initial stage was somewhat different from the results obtained by Palmer et al. [52]; these authors pointed out that the PM model tended to over-predict the ET when the LAI becomes small.Its imperfection is not surprising; the surface resistance of the PA-PM model in our study was fitted using the LAI as a power function, resulting in the fitted r s increasing rapidly with the decrease in the LAI (LAI < 0.5 m 2 /m 2 at the initial stage).The larger the r s , the smaller the ET simulated by the PM model, which has been confirmed in many studies [4,14,17,53,54].Conversely, the r s was underestimated at the late stage due to a higher LAI, which was because infinite types of tomato varieties were used in our study, resulting in less leaf aging at the late stage.But the performance of the PA-PM model was not greatly affected, as its overestimation was only 8.7% (Figure 6).During the development and middle periods, the LAI of tomatoes ranged from 1.0 to 3.5 m 2 /m 2 ; the PA-PM model can effectively predict the tomato ET, indicating that calculating the r s was feasible by using the LAI as a power function.So, the assumption that the underlying surface is treated as a "big leaf" in the PM model was the key to ensure its accuracy [7,14,22,33,48,52,54].In addition, such results also illustrated that surface resistance played a critical role in the performance of the PA-PM model, especially during the sparse canopy stage.
The PA-SW model over-predicted the ET by 19.5% and 13.6% at the initial and late stages, respectively (Figure 6).These results may be due to the lower tomato LAI, the overestimation of available energy intercepted by the canopy in the PA-SW model, as well as the restriction of soil water stress on water transport.Referring to local agronomic practices, tomato seedlings transplanted in greenhouses were irrigated with only 20 mm water to enhance the expansion of root system, and no water was applied before entering the development stage; additionally, a small amount of water was supplied at the late stage to improve the tomato quality.Thus, soil water stress frequently occurred at the initial and late stages.A number of studies have also reported that crop ET will be over-predicted by the SW model under water stress and low LAI conditions [13,50,[54][55][56][57][58].The PA-SW model performs well during the development and middle stages, which have been confirmed in our previous work [13].Additionally, the PA-SW model still has errors in estimating the tomato ET, which is related to the large number of parameters in the model.A large number of parameters may increase the uncertainty of the model [53].Notwithstanding its limitation, the PA-SW model indeed provides a simple method to predict tomato ET in greenhouse where the LAI > 1.0 m 2 /m 2 and without water stress.
The PA-PT and PA-CC models performed well at all four growth stages, and fewer parameters were embedded in these two models, among which two important coefficients, α PT and K c , determined the successful computation of tomato ET in the greenhouse.Herein, the α PT was obtained by inverting Equation ( 2) and was then fitted with days after transplanting; the K c was obtained by inverting Equation (3) and was calculated using the linear interpolation method.Both the α PT and K c were corrected using measured data from the 2016 study year.The PT model was a simplification based on the available energy penman method, in which the radiation factor was the main component [31].Yan et al. [59] indicated that the latent heat flux was the primary component of the net radiation over the whole growing season in the greenhouse and was principally constrained by net radiation.Similar results were also found elsewhere [3,[5][6][7]16,[59][60][61].So, the performance of the PA-PT model was conceivable.In our study, tomatoes were grown in the same greenhouse using a similar time, pattern, and management practice; thus, the K c fluctuated smoothly, although the greenhouse environment varies slightly in each year.Therefore, taking the K c as 0.49, 1.03, and 0.86 at the initial, middle, and late stages, respectively, was perfect for computing the tomato ET in the greenhouse.The values of K c in this study were lower than the values recommended by the FAO56 guidelines (0.6, 1.15, and 0.9), especially at the middle stage, which may be due to the higher humidity and lower air speed in the greenhouse.In addition, the walkway installed in the greenhouse reduces the canopy coverage, which was another reason for the decrease in the K c .Hanson and May [62] concluded that the K c of greenhouse tomatoes at the middle stage ranged from 0.99 to 1.08 based on the data from three study years.Qiu et al. [63] also reported that the tomato K c in the greenhouse at the mid-season mark was much lower than the values recommended by the FAO56 guidelines and fluctuated between 0.77 and 0.97.
The four parameterized models can effectively compute the daily ET of greenhousegrown tomatoes over the entire growth stage, although slight over-or underestimates are inevitable, but the errors were within the allowable range (MAE ranged from 0.38 to 0.51 mm/d) (Table 7 and Figure 7).The successful use of the models was mainly attributed to the consideration of the integrated restriction effect of the greenhouse microenvironment, plant morphology, and soil moisture on the model parameters.Therefore, the parameterized method of the PM, PT, CC, and SW models to estimate the ET for greenhouse-grown tomatoes with drip irrigation was acceptable.

Advantages and Disadvantages of Four Parameterized Models
The PA-PT and PA-CC models can provide an accurate estimation of ET for tomatoes in greenhouses at different growth stages, especially at the initial and late stages with low coverage and water deficit.Additionally, the PA-PT and PA-CC models were simpler and easier to use than the PA-PM and PA-SW models, due to fewer parameters.Of course, some errors still occur when estimating ET using these two models.Firstly, the discrepancy of R 2 and EF between ET estimated by the PA-PT model and measurements was larger than PA-PM and PA-CC models at the initial stage, which may be attribute to the lower available energy received by the canopy in sparse LAI stage.Secondly, the ET o should be confirmed accurately in the PA-CC model; however, the PM model fixed a constant r a of 308 s/m, which was acceptable.Therefore, putting forward the exact formula for calculating the ET o in a solar greenhouse was quite clearly necessary.
As for the PA-PM model, it provided an accurate ET estimate when the soil subsurface was fully covered; in this case, the underlying surface can be treated as a whole, completely satisfying the hypothesis of the "big leaf".Additionally, the performance of the PA-PM model also decreases when soil water was limited [64].Regrettably, the PA-PM model does not take into account the contribution of soil layers to evaporation; hence, it is not suitable for partial or sparse canopy conditions.Nyolei et al. [4] indicated that the application of the surface resistance model was not appropriate in the case where soil water vapor was the main source of the ET, especially during the early stages of crop development, as the approach underestimates the ET during this stage.Therefore, the LAI was used as a power function to calculate the r s at the initial stage, which may produce a great error.
The PA-SW model is probably the most cumbersome of these four parameterized models because many parameters should be adjusted before application, which may increase the uncertainty of the model.Even though five resistances in the PA-SW model were parameterized in our study, we still did not improve the accuracy when the LAI was below 1.0 m 2 /m 2 .Nyolei et al. [4] considered that developing a segmented method of surface resistance calibration can improve the performance of models to make them better adaptable to environmental changes.However, at present, it is still unclear how all environmental variables in greenhouses correlate with surface resistance, because the surface resistance was simultaneously affected by several environmental variables.Li et al. [28] also indicated that unsegmented models decrease the canopy resistance, thereby overestimating the ET.Notwithstanding its limitation, the PA-SW model does serve as a simple approach for users to choose in estimating the tomato ET under the dense canopy stage.Additionally, the PA-SW model can also be used to estimate soil evaporation and plant transpiration separately, which provides an effective tool for reducing "ineffective evaporation" from the soil layer.
In conclusion, our results can provide a simple method and theoretical basis for estimating crop ET in greenhouses, especially by parameterizing the Priestley-Taylor model, whose parameters are more readily available and more acceptable to farmers.

Conclusions
The application of mechanistic evapotranspiration models is often limited due to the large number of parameters and the need for calibration according to environmental variables.In this study, the Penman-Monteith, Priestley-Taylor, and Shuttleworth-Wallace models, and the crop coefficient method were parameterized, and their performance in the estimation of tomato ET in a greenhouse was estimated.The PA-PT model estimated the ET more accurately than the PA-PM, PA-CC, and PA-SW models, both at different stages and the entire growing season.Furthermore, the parameters in the PA-PT model were minimal and easier to obtain, thereby reducing the model uncertainty.The second recommended model was the PA-CC model; although the accuracy was slightly lower than the PA-PT model, its convenience may be more acceptable to farmers.The discrepancy was found between the ET estimated by the PA-PM and PA-SW models and measurements at the initial and late stages, where the underlying surface was sparse and crops frequently suffered from water stress.Considering the large number of parameters in the Shuttleworth-Wallace model, the PA-SW model provides a simple approach to predict crop ET in greenhouses, which is helpful for its practical applications.By comparing the structure and accuracy of these four parameterized models, the PA-PT model is recommended to estimate crop ET in a greenhouse.These parameterized models can be used in an easy and relatively accurate way for determining greenhouse tomato water consumption and developing appropriate irrigation scheduling.

Figure 1 .
Figure 1.Relationship between surface resistance (rs) in Penman-Monteith model and leaf area index (LAI) by using 2016 study year.The black line is the fitting regression line.

Figure 1 .
Figure 1.Relationship between surface resistance (r s ) in Penman-Monteith model and leaf area index (LAI) by using 2016 study year.The black line is the fitting regression line.

Figure 2 .
Figure 2. Relationship between Priestley-Taylor coefficient (α) and days after transplanting (DAT).Here α values were calculated by inverting the Priestley-Taylor model by using the 2016 ET data.The black line is the fitting regression line.
calculated by inverting weighing lysimeter Kc used in the CC model I: Initial II: Development III: Middle IV: Late

Figure 3 .
Figure 3. Variations of crop coefficient (K c ) at four growth stages in 2016 study year.Here blue circles o are K c values estimated by inverting the crop coefficient model using the 2016 ET data.The solid black line is the average value of K c at different growth stages.

Figure 4 .
Figure 4. Resistance parameters in the Shuttleworth-Wallace model and its relationship with the influence factors.(a) Variations in ra a and ra s during the growing period, (b) relationship between rs c and LAI, (c) relationship between ra c and LAI, (d) and relationship between rs s and θs in 2016 study year.The black line is the fitting regression line.

Figure 4 .
Figure 4. Resistance parameters in the Shuttleworth-Wallace model and its relationship with the influence factors.(a) Variations in r a a and r a s during the growing period, (b) relationship between r s c and LAI, (c) relationship between r a c and LAI, (d) and relationship between r s s and θ s in 2016 study year.The black line is the fitting regression line.

Figure 7 .
Figure 7.Comparison of daily ET estimated by the parameterized Penman-Monteith model (ET PA-PM , (a)), Priestley-Taylor model (ET PA-PT , (b)), crop coefficient method (ET PA-CC , (c)), Shuttleworth-Wallace model (ET PA-SW , (d)) and that measured by weighting lysimeter (ET WL ) over the whole growth stage in 2017 and 2019-2021.The black line is the linear fitting regression line, and the blue line is the 1:1 line.

Table 1 .
Information of growing stages and irrigation evens for greenhouse tomato during 2016-2017 and 2019-2021.Variations of crop coefficient (Kc) at four growth stages in 2016 study year.Here blue circles ο are Kc values estimated by inverting the crop coefficient model using the 2016 ET data.The solid black line is the average value of Kc at different growth stages.

Table 2 .
The different parameters required by four parameterized models.

Table 3 .
Statistical results of daily ET estimated using the parameterized Penman-Monteith model at four growth stages in 2017 and 2019-2021.

Table 4 .
Statistical results of daily ET estimated using the parameterized Priestley-Taylor model at four growth stages in 2017 and 2019-2021.

Table 5 .
Statistical results of daily ET estimated using the parameterized crop coefficient model at four growth stages in 2017 and 2019-2021.

Table 6 .
Statistical results of daily ET estimated using the parameterized Shuttleworth-Wallace model at four growth stages in 2017 and 2019-2021.At the initial stage, the PA-PM model presented a trend of underestimating the daily ET, whereas daily ET values were averagely overestimated by the PA-SW model in four study years.The statistical analysis shows that the PA-PM model seriously underestimated the ET by 27.0% (Figure 6(a1)), with a MAE of 0.27 mm/d, RMSE of 0.35 mm/d, and a positive PBIAS occurring in the PA-PM model (Table3).However, the PA-SW model significantly overestimated the ET by 19.5% (Figure6(a4)), with a MAE of 0.39 mm/d, RMSE of 0.56 mm/d, and a negative PBIAS occurring in the PA-SW model (Table

Table 7 .
Statistical results of daily ET estimated using the parameterized Penman-Monteith model (PA-PM), Priestley-Taylor model (PA-PT), crop coefficient model (PA-CC), Shuttleworth-Wallace model (PA-SW) over the whole growth stage in 2017 and 2019-2021.