Reference evapotranspiration estimate with missing climatic data and multiple linear regression models

Study area

The Adana Plain is located in the Mediterranean Region and consists of two parts, Çukurova Plain and Upper Plain. This study selected the Adana plain because It is Turkey’s most extensive and fertile delta plain. Moreover, several crops can be planted in the same crop year in Adana, where fertile cultivated land occupies a large area due to the suitable climate. Cereal, fruit, vegetable, and citrus are grown in the region, and their production efficiency is high (KafalıYılmaz, 2019). In Turkey, Adana Plain ranks first in citrus, watermelon, soybean, and peanut production (Kadeş, 2019). Rainfall amounts and distribution are insufficient during the growing season, so irrigation is indispensable in crop production in Adana Plain. This study includes April to October, covering the growing season of the main crops in the Adana Plain. Long-dated daily climate variables covering the period from 1998 to 2019 were used in the study.

This study used the Adana weather station in the center of Adana, operated by the Turkish State Meteorological Service (TSMS) (Fig. 1). Adana has a hot-summer Mediterranean climate. According to the Köppen-Geiger classification, Adana has a warm temperate climate with dry summer and hot summer (classified as Csa). (Öztürk, Çetinkaya & Aydın, 2017). The annual average (1929-2021) precipitation is 668.7 mm, of which about 50% of rainfall occurs in December, January, and February. The yearly average (1929–2021) sunshine hours and temperature are 7.5 h and 19.2 °C, respectively (Turkish State Meteorological Service, 2022).

FAO56-Penman-Monteith (PM) method

FAO recommends the FAO56-PM method as the sole standard approach to estimating ETo and validating other models, and It is expressed as (Allen et al., 1998): (1)${\displaystyle ETo=\frac{0.408\times \Delta \times \left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}\right)+\gamma \times \frac{900}{{T_m}_c+273}\times {\mathrm{U}}_2\times \left({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}\right)}{\Delta +\gamma \times \left(1+0.34\times {\mathrm{U}}_2\right)}}$ (2)${\displaystyle e_{\mathrm{s}}=\frac{0.6108\times exp\left[\frac{17.27\times Tma{x}_c}{Tma{x}_c+237.3}\right]+0.6108\times exp\left[\frac{17.27\times Tmi{n}_c}{Tmi{n}_c+237.3}\right]}{2}}$ (3)${\displaystyle e_{\mathrm{a}}=\frac{RH}{100}\times e_{\mathrm{s}}}$ (4)${\displaystyle \Delta =\frac{4098\times \left[0.6108\times exp\left(\frac{17.27\times {T_m}_c}{{T_m}_c+237.3}\right)\right]}{{\left({T_m}_c+237.3\right)}^2}}$

The wind speed data measured at 10 m was converted to a standard height of 2 m using Eq. (5) (Allen et al., 1998) (Eq. (5)). (5)${\displaystyle U_2=U_Z\times \left[\frac{4.87}{ln\left(67.8\times Z-5.42\right)}\right].}$ Where ETo = reference evapotranspiration (mm d⁻¹); R_n = net radiation (MJ m⁻² d⁻¹); G = soil heat flux density (MJ m⁻² d⁻¹), considered zero for daily estimates; Δ = slope of vapor pressure curve (kPa. °C⁻¹); γ = psychrometric constant (0.0672 kPa K⁻¹); Tm_c = mean air temperature (°C); Tmax_C = maximum air temperature (°C); Tmin_C = minimum air temperature (°C); U₂ = wind speed at 2 m height (m s⁻¹); Z = height of the wind speed measurement (10 m); e_s = saturation vapor pressure (kPa); e_a = actual vapor pressure (kPa); (e_s − e_a) = saturation vapor pressure deficit (kPa); RH = mean relative humidity (%).

(6)${\displaystyle Rn=Rns-Rnl}$ (7)${\displaystyle Rns=0.77\times Rs}$ (8)${\displaystyle Rnl=\left[\sigma \times \right.\left(\frac{Tma{x}_K^4+Tmi{n}_K^4}{2}\right)\times \left(0.34-0.14\times \sqrt{e_a}\right)\times \left(1.35\times \frac{R_s}{R_{so}}-0.35\right)}$ (9)${\displaystyle Rso=0.75xRa}$ where Rns = net solar radiation (MJ m⁻² d⁻¹); Rnl = net longwave radiation (MJ m⁻² d⁻¹); Rs = solar radiation (MJ m⁻² d⁻¹); Rso = clear-sky solar radiation; Ra = extraterrestrial radiation (MJ m⁻² d⁻¹); Tmax_k = maximum absolute temperature (K); Tmin_k = minimum absolute temperature (K).

Daily values of Ra were computed using the equations given in the FAO56 Paper by Allen et al. (1998).

Procedures for estimating missing climatic data

The FAO56-PM method can only be used when a full climatic dataset is available. However, Allen et al. (1998) suggested the following procedures for estimating absent Rs, U, Δe data:

Estimating missing solar radiation (Rs) data with the Angstrom formula

Rs is often estimated from sunshine data using the Angstrom formula (Eq. (10)) (10)${\displaystyle R_s=\left(a_s+b_s\times \frac{n}{N}\right)\times R_a}$

where N = daylight hours (hour); n = actual sunshine duration (hour). a_s = 0.25; b_s = 0.50 are recommended where no actual Rs data are available (Allen et al., 1998). This study used as a_s = 0.25 and b_s = 0.50. In this study, Rs estimated from sunshine data using the Angstrom formula is denoted by the symbol Rs_n.

Estimating missing solar radiation (Rs) data with the Hargreaves’ radiation formula

Where sunshine data are lacking, the Hargreaves’ radiation formula can be used to estimate Rs. (Eq. (11)) (Allen et al., 1998). (11)${\displaystyle R_s=K_{rs}\times \sqrt{\left(Tma{x}_C-Tmi{n}_c\right)}\times R_a}$ Where Krs = adjustment coefficient (≅0.16 for inland areas dominated by landmass; ≅ 0.19 for coastal areas adjacent to the coast of a large land mass). This study used 0.16 for Krs. In this study, Rs estimated from air temperature differences using the Hargreaves’ radiation formula is denoted by the symbol Rs_T.

Procedure for estimating wind speed (U)

Where wind speed data are not measured, mean wind speed data observed at a close location within the same homogenous region can be used instead of actual wind speed data, as proposed by Allen et al. (1998). This study used the closest weather station’s average wind speed (U₂) of 0.95 m s⁻¹ (1998–2019) instead of observed data.

Procedure for estimating vapor pressure deficit (es − ea =Δe)

Where air humidity data are missing, vapor pressure deficit ( Δe) can be calculated based on temperature data suggested by Allen et al. (1998). e_a is obtained by assuming that the dew point temperature (Tdew) is close to the daily minimum temperature (Tmin). Then, e_a iscalculated by Eq. (12). (12)${\displaystyle e_a=0.6108\times exp\left(\frac{17.27\times Tmi{n}_c}{Tmi{n}_c+237.3}\right)}$

Hargreaves-Samani (HS) and HS (A & G) equations

FAO suggests that when Rs, RH, and U data are missing, ETo can be estimated with the Hargreaves-Samani equation only using Tmax and Tmin (Hargreaves & Samani, 1985) and is expressed by Eq. (13): (13)${\displaystyle ETo=0.0023\times Ra\times {\left(Tma{x}_c-Tmi{n}_c\right)}^{0.5}\times \left({T_m}_c+17.8\right).}$

In this study, a modified HS equation developed by Almorox & Grieser (2016) for a warm temperate climate with dry summer (as classified Cs) was also used to estimate ETo (Eq. (14)). (14)${\displaystyle ETo=0.00419\times Ra\times {\left(Tma{x}_c-Tmi{n}_c\right)}^{0.2342}\times \left({T_m}_c+17.8\right)-0.0208.}$

Multiple linear regression (MLR) models

MLR is a statistical model that includes multiple predictive variables. The general form of the MLR model is: (15)${\displaystyle Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_n{X}_n}$ where Y is the dependent variable; X₁, …, and X_n are independent variables; ß₀ is the unknown intercept; ß₁, …, ß_n are partial regression coefficients of the function.

Multiple regression analysis was performed to develop linear models to estimate daily ETo values from meteorological variables, using ETo results based on calculations made by the FAO-56 PM equation (Tmax, Tmin, RH, U, Rs). The FAO56-PM ETo was employed as the dependent variable, while meteorological parameters were used as independent variables. The t-stat and p-values were examined to determine whether the best-fit line adequately fits the data. In several previous studies (Irmak et al., 2003; Noi, Degener & Kappas, 2017; Mattar & Alazba, 2019), one of the most common validation methods is dividing sample data into calibration and a validation dataset (for example, 70% and 30%, respectively). This study developed the model parameters using fifteen years (70%) of daily climatic data (April 1, 1998–October 31, 2012). Seven years (30%) of daily climatic data (April 1, 2013–October 31, 2019) were used to validation of the models.

Data analysis

The ETo estimated using the MLR models and FAO56-PM with missing data were compared with the ETo data computed using the FAO56-PM equation with complete datasets. The performance of the methods and models was determined by regression analysis. The performances of ETo estimates were also evaluated using the following equations below, suggested by Willmott (1981), Karunanithi et al. (1994), and Jacovides & Kontoyiannis (1995). All the statistical indices were calculated daily.

(16)${\displaystyle d=1-\left[\frac{\sum _{i=1}^n{\left(P_{i-}{O}_i\right)}^2}{\sum _{i=1}^n{\left(\left|P_i-\overline{O}\right|+\left|O_i-\overline{O}\right|\right)}^2}\right],\left(0\le d\le 1\right)}$ (17)${\displaystyle RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(P_i-O_i\right)}^2}}$ (18)${\displaystyle RE=\frac{RMSE}{\overline{O}}}$ (19)${\displaystyle MBE=\frac{1}{n}\sum _{i=1}^n\left(P_i-O_i\right)}$ (20)${\displaystyle MAE=\frac{1}{n}\sum _{i=1}^n\left|P_i-O_i\right|}$ where n = the total number of data; P_i = estimated ETo values (mm d⁻¹); O_i = FAO56-PM ETo values (mm d⁻¹); $\overline{O}$ = the mean of FAO56-PM ETo values (mm d⁻¹); d = the index of agreement; RMSE = root mean square errror; RE = relative error; MBE = mean bias error; MAE = mean absolute error.

In this study, the criteria given below, suggested by Stockle, Kjelgaard & Bellocchi (2004), were used to interpret the performance statistics of the methods in conjunction with the RMSE, MAE, MBE, b (slope), and a (intercept). The score is very good if d ≥ 0.95 and RE ≤ 0.10. The score is good if d ≥ 0.95 and 0.10 <RE ≤ 0.15. The score is acceptable if d ≥ 0.95 and 0.15 <RE ≤ 0.20. The score is marginal if d ≥ 0.95 and 0.20 <RE ≤ 0.25. Other combinations of d and RE values and R² <0.85 indicate poor performance.

Estimating daily ETo with missing data, using the FAO56-PM method, HS and HS(A&G) equations

The relationships between the ETo values estimated by the FAO56-PM method with complete data and those with the FAO56-PM method with missing data, HS, HS(A&G) equations are given in Fig. 2. It can be seen from Fig. 2 that the FAO56-PM ETo with missing U data correlated well with the FAO56-PM method with complete climatic data. Concerning regression equations, FAO56-PM with missing Rs data (Rs_n) gave the best-predicted values: a slope close to unity (0.9825) and intercept close to zero (0.25), and coefficient of determination close to one (R² = 0.9224). The FAO56-PM with missing e_a data gave the second-best values: (b = 0.9941, a = 0.2573, and R² = 0.8475). It is clear from Fig. 2 that all methods underestimated and overestimated ETo values during the growing season. However, the FAO56-PMT method using only a set of temperature data referred to in the literature as a reduced set FAO56-PM equation generally overestimated ETo for more days. HS(A&G) equations gave similar results, as well. The FAO56-PMT method overestimated the ETo during 2,907 days and underestimated the ETo during 1,801 days. These overestimates and underestimates were 3,696 days and 1,012 days for the HS equation, whereas 3,356 and 1,352 days for HS(A&G) equation, respectively.

Table 1 shows the performance statistics between ETo calculated by the FAO56-PM with a complete dataset and estimated by the FAO56-PM with missing data and estimated by the HS and HS(A&G) equations. The daily methods’ ranking was made, considering all statistical indices in Table 1. As shown in Table 1, when U data were missing, the closest weather station’s average wind speed of 0.95 m s⁻¹ was a perfect option for the FAO56-PM method because of minor errors. This method’s score was accepted as very good according to the criteria given by Stockle, Kjelgaard & Bellocchi (2004), and had the lowest RMSE, MAE, RE, and MBE values among the methods and ranked first (RMSE = 0.372 mm d⁻¹, MAE = 0.248 mm d⁻¹, RE = 0.088, and MBE = 0.009 mm d⁻¹). Sentelhas, Gillespie & Santos (2010) found similar results in Southern Ontario, Canada. When just U data were missing, the standard U data from a nearby weather station was a perfect option for the FAO56-PM method due to minor errors. Allen et al. (1998) also suggest using a wind speed value of 2 m s⁻¹ when wind speed data are unavailable. However, in our research, 98% of data from measurements showed U values below 2 m s⁻¹. The studies by Popova, Kercheva & Pereira (2006) in South Bulgaria and Córdova et al. (2015) in the high Andes of southern Ecuador presented that the results of RMSE and MBE were near zero when using a wind speed value of 2 m.s⁻¹. Djaman et al. (2016) reported that the performance of FAO56-PM was unsuitable in dry conditions when U was considered as 2 m s⁻¹. However, when using the daily average U in the same conditions, FAO56-PM estimated ETo accurately with MBE values between −0.05 to 0.04.

When e_a, usually obtained from relative humidity, was missing; in the case of estimating e_a by assuming that dew point temperature (Tdew) is close to the daily minimum temperature (Tmin), the agreement between ETo calculated with the complete data and ETo estimated with missing data indicated good performance due to low errors (RMSE = 0.564 mm d⁻¹, MAE = 0.421 mm d⁻¹, RE = 0.133, and MBE = 0.232 mm d⁻¹). Several locations presented similar results when e_a was estimated by Tmin (Jabloun & Sahli, 2008; Popova, Kercheva & Pereira, 2006; Djaman et al., 2016). The study by Jabloun & Sahli (2008) obtained good results for sub-humid and semi-arid locations of Tunisia when e_a was estimated by Tmin. Excluding Tunis, for all other locations, the RMSE and MBE values ranged from 0.233 to 0.303 mm d⁻¹ and from −0.212 to −0.078 mm d⁻¹. In the study performed by Popova, Kercheva & Pereira (2006) in the Trace plain of south Bulgaria, when e_a was estimated by Tmin, regression coefficients (R²) were found close to 1, the slope of simple linear regression (b) ranged from 0.93 to 0.97, and the values of the standard error of estimates (SEE) ranged from 0.15 to 0.27 mm d⁻¹. Similarly, Djaman et al. (2016) obtained good results in different agro-ecological zones of Burkina Faso. In the study, the values of b varied from 0.90 to 1.14, R² ranged from 0.77 to 0.92, and RMSE and MBE varied from 0.07 to 0.22 mm d⁻¹ and from −0.45 to 0.22 mm d⁻¹, respectively.

When Rs data were missing, using the Angstrom formula to estimate Rs was a perfect option, with very low errors (RMSE = 0.390 mm d⁻¹, MAE = 0.272 mm d⁻¹, RE = 0.092, and MBE = 0.176 mm d⁻¹). Its score can be accepted as very good. If sunshine hours were not measured, using Hargreaves’ radiation formula to estimate Rs was unacceptable due to d < 0.95 and R² <0.85. So, Its score can be accepted as poor. In addition, the slopes of regression and intercept between ETo calculated with the complete data and those estimated with missing Rs data, b = 0.7398 and a = 0.9435, are not good values. However, when Hargreaves’ radiation formula calculates the Rs (Rs_T), the ETo estimates were more accurate than the FAO56-PMT method, HS, and HS(A&G) equations according to all statistical indices (Table 1). Different studies using sunshine and temperature-based models to estimate Rs observed similar results to our research. The study by Jahani, Dinpashoh & Nafchi (2017) in Iran showed that sunshine-based models had more accurate estimations than temperature-based ones. The study by Trnka et al. (2005) in Central European lowlands presented that the Angstrom-Prescott model was the best method for estimating Rs. Aladenola & Madramootoo (2014) recommended the Hargreaves’ radiation formula for estimating Rs in Canada’s absence of solar radiation and sunshine data, unlike our research. More research is needed to find the best model for estimating Rs for different locations.

As seen in Table 1, for conditions in which U and e_a data were missing (-U, -e_a), the ETo estimates were very similar to those estimating e_a from minimum temperature (−e_a). However, when U and e_a data were missing, the slope was not so close to the unity ( b = 0.7345), and the intercept was not so close to zero ( a = 1.3922). In addition, due to d < 0.95, this combination indicated poor performance. The other combinations were unacceptable according to statistical indices and presented poor performance. As shown in Table 1, when RE is expressed as a percentage, these values range from 8.8% to 22.6% according to methods. The HS method showed the highest error at 22.6%, and the FAO56-PM method with missing U data showed the lowest error at 8.8% in estimating daily ETo. HS equation gave the highest MBE (0.550 mm d⁻¹), whereas the FAO56-PM method with missing U data gave the lowest MBE (0.0009 mm d⁻¹).

A study in a semi-arid highland environment by Benli et al. (2010) evaluated the performance of the FAO56-PM method for cases where RH, U, Rs, or all three parameters would be missing, comparing with weighing lysimeter measurements. They estimated Rs using Hargreaves’ radiation formula and estimated e_a accepting Tdew = Tmin. They used the average wind speed of U₂ = 2 m s⁻¹ instead of actual wind speed data because this value is the average wind speed of 2000 weather stations around the globe. The FAO56-PM with complete data ranked first in the study (RMSE = 1.32 mm d⁻¹), FAO56-PM with missing U data ranked second (RMSE = 1.36 mm d⁻¹), and FAO56-PM with missing e_a data ranked third (RMSE = 1.34 mm d⁻¹). Similar to our research, the FAO56-PM method with missing U data estimated ETo accurately in their study. However, the HS equation (RMSE = 1.40 mm d⁻¹; MAE = 54 mm) estimated ETo slightly better than the FAO56-PMT method (RMSE = 1.40 mm d⁻¹; MAE = 55 mm) in their study, unlike our study.

A study in Southern Ontario, Canada, by Sentelhas, Gillespie & Santos (2010) evaluated the performance of FAO56-PM and alternative methods when data were missing. They used mean wind speed data observed in a close area. Using the procedure given in FAO56 Paper, they estimated Rs as a function of air temperature and also estimated vapor pressure deficit based on temperature data. The results proved that the FAO56-PM method estimated daily ETo accurately when wind speed and relative humidity data or all two parameters were missing. However, estimating Rs from temperature differences was not accurate enough in Southern Ontario because this approach resulted in poor ETo estimation. These results are in agreement with our study.

In our study, the FAO56-PMT method, using only a set of temperature data, could not accurately estimate daily ETo (Table 1). However, It provided more accurate results than HS and HS (A&G) equations. Similar results have been found in the studies performed by Trajkovic (2005) in Serbia, Popova, Kercheva & Pereira (2006) in the Trace plain area of south Bulgaria, Jabloun & Sahli (2008) in Tunisia, and López-Moreno, Hess & White (2009) in Pyrenees. However, Martinez & Thepadia (2010) reported that the HS equation produced more minor overestimation errors than the FAO56-PMT method. In their study, the FAO56-PMT method indicated the greatest errors in coastal areas, while the HS equation showed the greatest errors inland and island locations in Florida. Todorović, Karic & Pereira (2013) reported that where aridity prevails, the HS equation’s results were likely better than those of PMT. In contrast, PMT results were better for less arid climates, from semi-arid to humid. Martinez-Cob & Tejero-Juste (2004) determined that in windless locations of a semi-arid region in Spain, the HS equation gave unreliable estimates for daily ETo, and they recommended that the HS equation be used for 10-day periods. Torres, Walker & McKee (2011) reported that simple methods give less accurate results for daily ETo estimates than weekly and monthly ETo estimates.

Estimating daily ETo using MLR models

Table 2 presents the MLR models developed for the thirteen combinations. The regression analysis (t-stat, p-value) was made to determine independent variables significantly affecting ETo at a 95% confidence level (Table 3). The values of the t-stat should be greater than 1.96 or less than −1.96 to verify the statistical goodness of the regression coefficients. The absolute values of the t-stat for all MLR models were greater than 1.96, except for Tmin in the MLR-2 model and Tmax in the MLR-9 model; the p-values confirmed this. The most significant variable in the MLR-4, MLR-8, and MLR-11 models is Rs, with the highest t-stat values. The most significant variable in the MLR-6, MLR-10, and MLR-13 models is Rs_n, with the highest t-stat values. This indicates that Rs and n variables have more effect on estimating ETo than other variables in Adana Plain conditions. Similarly, Pereira et al. (2015) reported that Rs could dominate ETo estimation in mostly humid and sub-humid climates during summer, where the relative power of the vapor pressure deficit and wind term of the FAO56-PM is small compared to the radiation term. Several studies observed similar results to our research (Dimitriadou & Nikolakopoulos, 2022; Yirga, 2019; Perugu, Singam & Kamasani, 2013; Choi et al., 2018). In a study in Peloponnese, Greece, by Dimitriadou & Nikolakopoulos (2022), sunshine hours (n) and net radiation (Rn) significantly affected ETo more than the other variables. Among one-input models, the MLR model with n as the input variable exhibited the best performance (RMSE = 0.409 mm d⁻¹, MAE = 0.320 mm d⁻¹, R² = 0.960). In contrast, the MLR model with U as the sole explanatory variable had the poorest performance in their research (RMSE = 2.001 mm d⁻¹, MAE = 1.953 mm d⁻¹, R² = 0.030). In a study by Yirga (2019), for the Megecha catchment in Ethiopia to model ETo with MLR, the results of the multiple correlations showed that n was a strong positive correlation with ETo (r = 0.82), and the effect of n on ETo was more than other variables. Perugu, Singam & Kamasani (2013) performed a study in India, and the correlation analysis indicated that n, T, and U influenced the ETo more than other variables.

RH data in the five models (MLR-3, MLR-7, MLR-11, MLR-12, and MLR-13) had less impact on the estimation of ETo with small t-stat values (Table 3). In a study in South Korea by Choi et al. (2018), Tmean, Tmin, Tmax, and n had strong positive correlations with ETo. In contrast, RH negatively correlated and had less effect on ETo, similar to our research. According to a study in China by Gong et al. (2006), FAO56-PM ETo was most sensitive to RH among the four climate variables, which vary by season and region and are affected by regional wind speed patterns, and followed this by Rs, T, and U. In the study by Nandagiri & Kovoor (2006), factor analysis indicated that U appears to be an essential variable in the arid climate. In contrast, n seems more dominant in subhumid and humid climates. Mohsin & Lone (2021) observed in the temperate Kashmir Valley for all the stations the strongest correlation of ETo with Tmin. In Egypt, Mattar & Alazba (2019) modeled ETo from climatic variables using GEP and MLR models. The results showed that using RH, U, or both in the models increased the performance of GEP and MLR models. Niaghi, Hassanijalilian & Shiri (2021) applied MLR models to Red River Valley with a subhumid climate to estimate ETo from climatic variables. They grouped input combinations into air temperature-based (Tmax, Tmin), mass-transfer-based (Tmax, Tmin, U), and radiation-based (Rs, Tmax, Tmin). Similar to our research, the radiation-based MLR model performed better than other combinations in their study (RMSE = 0. 68 mm d⁻¹, MAE = 0.51 mm d⁻¹, and R² = 88%). In our study, as shown in Tables 4 and 5, the MLR-4 model (Tmax, Tmin, Rs) showed better performance than the MLR1 model (Tmax, Tmin) and MLR2 model (Tmax, Tmin, U).

In Fig. 3, MLR model estimates were compared to the FAO56- PM method via scatter plots for the calibration and validation processes. The b (slope) values are close to one, and a (intercept) values are close to zero in MLR-8 and MLR-11 models with R² of 0.9528, 0.9530 in the calibration process, and R² of 0.9514, 0.9515 in the validation process, respectively. Also, MLR-4, MLR-10, and MLR-13 models performed well according to their b, a, and R² values.

The values of RMSE, MAE, RE, R², and d for MLR models during the calibration and validation process are presented in Tables 4 and 5. Among all models, the MLR1 model, estimating ETo with only Tmax and Tmin, showed the lowest performance with the values of RMSE of 0.915, 1.117 mm d⁻¹, MAE of 0.750, 0.940 mm d⁻¹, RE (%) of 22.5, 24.2, R² of 0.3498, 0.4234 in the calibration and validation process, respectively. When Rs data were added to the MLR-1 model, the ETo prediction accuracy of the new model (MLR-4) increased by more than 52% for RMSE, MAE, and RE in the calibration and validation processes. In the MLR-6 model, Rs was estimated by n data using the Angstrom formula, and the MLR-6 model was designed with Tmax, Tmin, and Rs_n. MLR-6 model estimated daily ETo with more than 32% accuracy than the MLR-1 model for RMSE, MAE, and RE values in the calibration and validation processes. Similarly, in a study by Yamaç (2021) in a semi-arid highland environment, when Rs data were added to the ANN1 model, which has Tmax and Tmin, the ETo estimation accuracy of the new model (ANN2) for MSE, RMSE, and MAE values increased by more than 30% in training and testing subsets.

In the MLR-5 model, Rs was estimated by Tmax and Tmin data using the Hargreaves radiation formula, and the MLR-5 model was designed with Tmax, Tmin, and Rs_T. The results indicated that the MLR-5 model predicted daily ETo with more than 21% accuracy than the MLR-1 model, based on RMSE, MAE, RE, and R² values in the calibration and validation processes. When RH data were added to the MLR-1 model, the ETo estimation accuracy of the new model (MLR-3) slightly improved with a range of 0.8% to 1.7% in the calibration and validation processes for RMSE, MAE, RE, and R² values.

Only five models, MLR-4, MLR-8, MLR-10, MLR-11, and MLR-13, accurately predicted the daily ETo in the calibration and validation processes (Tables 4 and 5). The results showed that the MLR-11 model had the smallest RMSE (0.246−0.288mm d⁻¹), MAE (0.183−0.210mm d⁻¹), and RE values (6.1−6.2%) and the highest R² (0.9530−0.9515) and d (0.999−0.999) for the calibration and validation processes, respectively. The statistical results of the MLR-8 model were nearly the same as the MLR-11 model (RMSE = 0.246−0.289 mm d⁻¹; MAE = 0.184−0.210 mm d⁻¹; RE = 6.1−6.3%; R² = 0.9528−0.9514; d = 0.999−0.999). According to this, MLR-8 and MLR-11 models’ score was very good according to the criteria given by Stockle, Kjelgaard & Bellocchi (2004) for the calibration and validation processes. Therefore, the MLR-8 model is recommended in this study because it uses lesser climate data than the MLR-11 model. The MLR-8 model, which does not use RH data, can be successfully used to predict daily ETo in the study area. MLR-10 and MLR-13 models’ scores for the validation process according to mentioned criteria were very good, and the MLR-4 model’s score was good. Generally, the results of the models agree with the literature (Tabari et al., 2012; Dai et al., 2008; Yamaç, 2021), concluding that more climatic variables commonly increase model estimation accuracy. Tabari et al. (2012) evaluated MLR models’ different input combinations in a semi-arid highland environment in Iran. Similar to our study, the lowest performance was observed in the MLR model that only used Tmean (RMSE = 0.765 mm d⁻¹, MAE = 0.621 mm d⁻¹, R² = 0.92) while the MLR model that used Tmean, Rs, RH, and U had the best performance (RMSE = 0.552 mm d⁻¹, MAE = 0.442 mm d⁻¹, R² = 0.96) in their study. Dai et al. (2008) assessed the MLR models in northern China’s arid, semi-arid, and sub-humid areas. In the sub-humid area, the MLR models with four inputs (Tmean, RH, n, and U) and three inputs (Tmean, RH, n) gave similar results to the ANN model with the same inputs. MLR models predicted ETo accurately with RMSE = 0.202−0.273mm d⁻¹, RE = 5.2−6.9%, R² = 0.961−0.930, respectively. Yamaç (2021), in a semi-arid highland environment, evaluated ANN and kNN models. The results showed that more climatic variables increased the models’ predictive abilities, and the models with 6 inputs gave the best results in the study.

The main implication of our study is that MLR models developed for Adana Plain for estimating ETo have very satisfactory performance. The findings from our study are significant because agriculture accounts for 69% of global water withdrawals, mainly used for irrigation. This ratio can reach 95% in some developing countries (FAO, 2011). The FAO estimates that the world will need about 60% more food by 2050 and that irrigated food production will increase by more than 50% over the same period (FAO, 2017). To ensure sustainable development and water supply, particularly in arid environments, irrigation professionals need tools to estimate ETo at a large scale (Khoshravesh, Sefidkouhi & Valipour, 2017). Our study’s findings are significant for the agricultural managers and irrigation engineers in a region with similar climatic conditions to estimate ETo.