Study area

The state of Chihuahua is located in the north-central region of Mexico (31° 47' 04" and 25° 33' 32" north latitude, and 103° 18' 24" and 109° 04' 30" west longitude). Chihuahua represents 12.6 % of the national territory, is bordered to the northeast by the United States of America, to the east by Coahuila de Zaragoza and Durango, to the south by Durango and Sinaloa, and to the west by Sonora and Sinaloa (Figure 1). The territory has four types of climates: very dry (40 %), dry and semi-dry (33 %), temperate-sub-humid (24 %) and warm-sub-humid (3 %). The average annual temperature is 17 °C, and the average annual precipitation is 500 mm (^{INEGI, 2017}).

Figure 1 Study area and location of weather stations. 

Meteorological data and processing

Daily temperature and precipitation records were obtained from 33 meteorological stations managed by the Servicio Meteorológico Nacional (SMN) (^{CONAGUA-SMN, 2019}) for the period 1960-2013 (Figure 1, Table 1). This data series includes outlier data cleaning, which is performed by SMN staff. Data were subjected to the double mass curve test, where each station was compared to surrounding stations for consistency, as recommended by ^{Gao et al. (2017)}, ^{Ruíz-Álvarez et al. (2020)} and ^{Searcy and Hardison (1960)}, and all measurements were consistent. A homogeneity test was also performed using the PROCLIMDB package developed for the R program (^{Štěpánek et al., 2009}), which works similarly to CLIMATOL (^{Guijarro, 2016}). This analysis showed that climatic information from all the stations used is homogeneous.

Relative humidity and wind speed records were collected from 20 automatic weather stations (EMAs), managed by the Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias (INIFAP), for the latest available period (1 to 3 years) (^{INIFAP, 2019}). These records were used to validate the relative humidity values required to operate the PMMF method. Finally, solar radiation was calculated using the equation described by ^{Allen et al. (1998)}.

The missing values of databases were estimated with the interpolation method known as inverse distance weighted (IDW), according to Equation 1 (^{Shepard, 1968}):

Where Zp is the estimated value for the point, n is the number of points used in the interpolation, Z _i is the known meteorological record, d _i is the distance between stations (m), and p is the power depending on the degree of terrain undulation: 1 for flat and 2 for steep.

Empirical methods for calculating reference evapotranspiration

The empirical methods used to calculate ET_o were selected based on the availability of meteorological information in the study area. Table 2 shows the empirical methods used to calculate ET_o. The reference method was PMMF, while the methods compared or evaluated were Hargreaves (Hg), Turc (Tc), Blaney Criddle (BC), Thornthwaite (Th), Malmstrom (Mm), Makkink (Mk) and Hargreaves modified by Droogers (Dg).

Statistical analysis

Comparison of means of each empirical method regarding the reference method (PMMF) was performed with an analysis of variance (P ≤ 0.05) and by the Tukey-Kramer least significant difference (LSD). The hypotheses tested were: H_O) the mean values of each empirical method is equal to that of the reference method. H_A) at least one mean value of the empirical methods is different from that of the reference method. A linear regression analysis was performed between the values of each empirical method against those of the reference method. In this case, the hypotheses tested were H_O) the slope of the regression line is equal to 0 and H_A) the slope of the regression line is different from 0. The precision was evaluated with the coefficient of determination (R²), root mean square error (RMSE), mean error (ME) and the Nash-Sutcliffe efficiency (NSE) calculated according to Equations 2 to 5 (^{Alexandris et al., 2006}; ^{Cervantes-Osornio et al., 2012}; ^{Djaman et al., 2018}; ^{Nash & Sutcliffe, 1970}; ^{Tabari, 2010}).

where a _i is the value estimated by the empirical method, t _i is the value obtained with the reference method, ā is the average of the estimates with the empirical method, t-is the average of the reference method values, and N is the total number of measurements.

Spatial and temporal analysis

The geographic information was organized in a Geographic Information System (GIS) with QGIS software ver. 3.8 "Zanzibar". The IDW method was used to obtain regular interpolated surfaces of monthly, seasonal, and average annual values of meteorological variables and ET_o with a 1 km x 1 km pixel size. A spatio-temporal analysis of ET_o values was performed with the maps obtained, including precipitation and temperature spatial and temporal distribution. Due to the volume of information generated in this study, only the annual and seasonal maps: spring (March-May), summer (June-August), autumn (September-November) and winter (December-February), are shown.

Temporal analysis of precipitation and temperature

Figure 2 shows the maximum and minimum average annual temperature and the annual cumulative precipitation for 1960-2013. The maximum temperature in this period ranged from 28.5 °C (2011) to 25 °C, while the minimum temperature was between 10.7 and 7.7 °C (1973). On the other hand, the average annual cumulative precipitation was 270 mm; the highest value was in 2004 (493 mm) and the lowest in 1960 (55 mm).

Figure 2 Mean annual temperature and annual cumulative precipitation in Chihuahua, Mexico (1960-2013). 

Chihuahua shows a monomodal cycle concerning temperature and precipitation (Figure 3). The average maximum temperature ranged from 18 °C in January to 34.4 °C in June, while the minimum temperature varied from 0.3 °C in January to 17 °C in July. The highest monthly average precipitation was in July (66.6 mm), and the lowest value (4.8 mm) was recorded in March. The figure shows that the high-temperature season coincides with the rainy season (June-September). The dry season coincides with low temperatures (October-March); therefore, the planning of sowing and harvesting is necessary.

Figure 3 Mean monthly temperature (maximum and minimum) and mean monthly precipitation (1960-2013). 

Spatial distribution of temperature and precipitation

Chihuahua is the largest state in Mexico, which generates significant differences in the distribution of temperature and precipitation. Figure 4 shows the spatial distribution of temperatures (maximum and minimum) and average seasonal rainfall in the state. The maximum temperature (36.5 °C) occurs in summer in Ojinaga (800 m a. s. l.), while the minimum (14.2 °C) happens in winter in Maderas (2 110 m a. s. l.). Therefore, agriculture is mainly carried out in the spring-summer cycle (May-October).

Figure 4 Spatial distribution of seasonal temperature (maximum and minimum) and precipitation in Chihuahua (1960-2013). 

In Chihuahua, temperature decreases at a rate of 6 °C per 1 000 m ascent, according to the linear regression model in Figure 5a (R² = 0.74). Temperature varies from 34 °C at 169 m a. s. l. to 15 °C at 3 286 m a. s. l.

Figure 5 Linear regression between elevation (m a. s. l.) versus a) mean annual temperature (°C) and b) cumulative annual precipitation (mm) in weather stations of the Comisión Nacional del Agua (CONAGUA) in Chihuahua, Mexico.  

Regarding precipitation, the maximum value was for the municipality of Bachiniva (99 mm, 2 017 m a. s. l.), while the minimum was in spring in Ojinaga (1.7 mm, 770 m a. s. l.). It is estimated that in 75 % of the state's surface, the annual cumulative precipitation is less than 300 mm, with a gradient that increases from east to west: from 73 mm in the region adjacent to the United States to a maximum of 552 mm in the Sierra Madre Occidental.

The precipitation-height relationship is complex. Figure 5b shows the model obtained by linear regression of this relationship. The R² value was 0.28 (P ≤ 0.05), with a precipitation gradient increasing at a rate of 9.7 mm per 1 000 m of ascent. This relationship, although low, indicates that precipitation could be influenced by other variables, such as terrain relief and slope, especially in mountainous regions (^{Geng et al., 2017}).

Reference evapotranspiration values

The statistical parameters (mean, standard deviation, quartiles) and outliers of ET_o allows comparing the methods used and studying their frequency distribution. The diagrams in Figure 6 show the variation of the average daily ET_o calculated with the eight empirical methods, which indicate a mean ranging from 2.5 to 5.3 mm·day^-1 and a standard deviation between 1 and 2.2 mm·day^-1. The PMMF method reached an average value of 5.3 mm·day^-1, with a standard deviation of 2.2 mm·day^-1. The Dg and Hg methods are the closest to PMMF, with a mean of 5.2 and 4.5 mm·day^-1, respectively; this represents underestimating ET_o values of 2 and 15 % concerning PMMF. This coincides with that reported by ^{Allen et al. (1998)}. They indicate that the methods that come closest to PMMF consider temperature and radiation as the primary sources of energy that promote evapotranspiration, like Hg and Dg.

Figure 6 Box plot for reference evapotranspiration (ET_o) values estimated from eight empirical methods. 

The methods with the lowest dispersion were Mm (1.5 to 4.5 mm·day^-1) and BC (2.2 to 5.3 mm·day^-1), which coincides with that reported by ^{Le et al. (2012)}, who indicate that in regions with high variation of climate, height and land use conditions, empirical methods do not always provide satisfactory results. PMMF and Dg methods had the highest values of standard deviation (S) (2.2 and 1.9 mm·day^-1, respectively), while Mm and BC had the lowest values of S (1.0 mm·day^-1, in both cases).

Statistical analysis

The comparative statistical analysis between the empirical methods identified the closest methods to the reference method (PMMF). Table 3 and Figure 7 show the analysis of variance for the comparison of means and linear regression between ET_o values estimated by empirical methods and the PMMF reference method. The comparison of means shows that, except for the Dg method, the empirical methods show statistically significant differences (P > 0.05) concerning PMMF. In the case of linear regression, all methods showed a statistically significant linear relationship with the PMMF method.

Figure 7 Linear regression between the seven empirical methods and Penman-Monteith method modified by the FAO (PMMF). Hg = Hargreaves; Tc = Turc; BC = Blaney Criddle; Th = Thornthwaite; Mm = Malmstrom; Mk = Makkink; Dg = Hargreaves modified by Droogers. 

All methods had R² values above the minimum value (0.70) reported by ^{Moriasi et al. (2007)}, indicating a good fit concerning the reference method. RMSE values ranged from 1.2 to 2.9 mm·day^-1, while ME ranged from -0.7 to -2.7 mm·day^-1and NSE went from -0.55 to 0.71; this suggests that all methods underestimated ET_o calculated with PMMF. The results obtained agree with those reported by ^{Lozoya et al. (2014)}. The method closest to PMMF was Hg with R² = 0.83, RMSE = 1.1632 mm·day^-1, ME = -0.695 mm·day^-1, NSE = 0.71. It can be an alternative to calculate ET_o in regions where only temperature records are available, as in Chihuahua, since the value of extraterrestrial radiation can be obtained from tables or with the equation recommended by ^{Allen et al. (1998)}. However, if solar radiation measurements, wind speed, temperature and humidity are available, it would be feasible to use PMMF (^{Lozoya et al., 2014}).

The comparative analysis shows that the methods that consider temperature and radiation (Hg, Tc, BC, Th and Mk) performed better regarding PMMF than those using temperature and other variables such as relative humidity and precipitation (Mn and Dg).

Analysis of monthly reference evapotranspiration

A comparative analysis of the average daily values of ET_o on a monthly scale, estimated with the different empirical methods, shows that all the methods studied have similar results throughout the year, although with significant differences depending on the season (Figure 8). According to the results, the maximum ET_o values occur between June and July (4.5 to 8.3 mm·day^-1), and the minimum between January and December (0.5 to 2.5 mm·day^-1). Six methods underestimate PMMF most of the year with values ranging from 0.1 to 4.1 mm·day^-1 (Hg, Tc, BC, Th, Mm and Mk). The months with the most considerable differences concerning PMMF (between 2 and 4 mm·day^-1) are in the rainy spring and summer seasons (February and October), while the minor differences (< 2 mm·day^-1) correspond to late autumn and early winter (November, December and January) (Figure 8), which is in agreement with that reported by ^{Ayyoub et al. (2017)}. The Dg method is the only one that overestimates PMMF almost throughout the year, except in July, August, and September.

Figure 8 Temporal monthly distribution of reference evapotranspiration (ET_o) estimated with empirical methods and differences concerning the reference method (Penman-Monteith, modified by the FAO (PMMF). 

The most significant differences concerning PMMF occur in June with values ranging from 1.3 to 4.1 mm·day^-1, while the minor differences (0.2 to 1.8 mm·day^-1) occur in December and January. The closest to PMMF is Hg and Dg, with differences ranging from 0.1 to 1.8 mm·day^-1. These extremes occur in December and August, respectively, indicating that the most significant differences are concentrated in the months when the temperature begins to decrease, coinciding with the rainy seasons, similar to that obtained by ^{Lozoya et al. (2014)}. These results indicate that estimating ET_o closest to the reference method (PMMF) on a monthly scale is Hg. It was used throughout the state to obtain homogeneous surfaces of monthly, seasonal and annual values.

Spatial distribution of seasonal reference evapotranspiration

The spatial interpolation of point values of ET_o in areas where there is no climate information allows the elaboration of maps to analyze spatial variation. Figure 9 shows seasonal average maps generated from the spatial interpolation of daily ET_o calculated with the Hg method. The highest ET_o values were reported in the spring (5.6 to 7.1 mm·day^-1) and summer (4.9 to 6.5 mm·day^-1), increasing from higher to lower elevations in a west-east direction. The highest ET_o value (7 mm·day^-1) occurred in the spring at the Colina station (key 8 311, at 1 214 m a. s. l.). The lowest ET_o values occurred in the fall (2.5 to 3.6 mm·day^-1) and winter (2.6 to 3.9 mm·day^-1), with a predominant increase from high to low elevations from north to south. The lowest ET_o value (3.5 mm·day^-1) occurred in the fall at the station La Huerta, Topia (key 10037, at 670 m).

The behavior mentioned above indicates a strong correlation between ET_o, time of the year and elevation, especially during spring-summer and in the direction of the Sierra Madre Occidental. It coincides with that reported by ^{Sun et al. (2020)}. They indicate that ET_o gradually increases with increasing elevation up to a threshold altitude (approximately 1 900 m a. s. l. in the study region) and then decreases. It may be because, during the warmer months, solar radiation is the primary source of energy for evapotranspiration, a dependence that intensifies with increasing altitude (^{Allen et al., 1998}). In contrast, lower radiation availability in fall-winter generates lower ET_o values, causing the spatial gradient to have a predominant relationship with latitude rather than elevation, consistent with that reported by ^{Tang et al. (2019)}. Similar results to those of the PMMF model were obtained by ^{Monterroso-Rivas et al. (2016)}, ^{Ruíz-Álvarez et al. (2014)} and ^{Sheffield et al. (2010)}.

Figure 9 Spatial distribution of reference evapotranspiration (ET_o) estimated by the Hargreaves method (1960-2013). 

The seasonal analysis shows significant differences in the spatial distribution of ET_o in Chihuahua, about the time of the year and altitude. On the other hand, for agricultural purposes, it is essential to represent the spatial distribution of the daily averages and annual cumulative values obtained with the Hg method. A standard deviation map of the annual ET_o is also added to identify the areas with significant variability (Figure 10), which is helpful for water resource planning and management.

Figure 10 Spatial distribution of reference evapotranspiration (ET_o) (a) annual cumulative (b) daily average, (c) standard deviation estimated with the Hargreaves method (period 1960-2013). 

Monthly or annual cumulative ET_o values are crucial for irrigated agriculture when a water balance is required to know the soil moisture content and determine periods of deficit or excess. Annual values ranged from 1 435 to 1 862 mm, while daily values ranged from 3.9 to 5.1 mm·day^-1. These values are similar to those obtained by ^{Orozco-Corral (2010)} and ^{Segura-Castruita and Ortiz-Solorio (2017)}.

Figures 10a and 10b show that the annual cumulative and daily average ET_o values have similar spatial behavior. It can be observed that the minimum values are concentrated in the western and central regions (near the Peñitas and Majalca stations). A radial increase occurs towards the southern area (near the La Huerta station) and north (near the Dublan station) in the direction of the Sierra Madre Occidental. This behavior indicates that elevation plays an essential role in the distribution of ET_o. An increase in height leads to a rise in ET_o, especially in the warm seasons (spring-summer). It is because radiation is the primary energy source for ET_o, as reported by ^{Sun et al. (2020)} and ^{Tang et al. (2019)}.

The standard deviation is a measure of the dispersion of the data concerning the mean. In this sense, Figure 10c shows that the spatial variation of the daily mean values of ET_o ranges between 1.5 and 2.0 mm·day^-1 in the direction of the Sierra Madre Occidental. It confirms the effect of altitude on ET_o estimation since it is observed that the higher elevations have lower variation compared with the great plains of the state had more significant variation during the period analyzed.