Spatiotemporal Estimation of Reference Evapotranspiration for Agricultural Applications in Punjab, Pakistan

Abstract

Estimation of reference evapotranspiration (ET_o) is a key element in water resources management and crop water requirement which, in turn, affects irrigation scheduling. ET_o is subject to the influence of various climatic parameters including minimum temperature (T_min), maximum temperature (T_max), relative humidity (RH), windspeed (WS), and sunshine hours (SH). Usually, the influence of the climatic parameters and a dominating climatic factor influencing ET_o is estimated on yearly basis. However, in diverse climatic regions, ET_o varies with the varying climate. Therefore, this study aims to estimate the spatiotemporal variation in the influence of the climatic parameters on ET_o in Punjab, Pakistan, for the period 1950–2021, specifically focusing on decennial, annual, and monthly patterns. The study area was divided into five agroclimatic zones. The Penman–Monteith model was used to estimate ET_o. The influence was assessed using geographic weighted regression (GWR) and multiscale geographic weighted regression (MGWR) as the primary methods. As per results from MGWR, ET_o in Punjab was highly influenced by the T_min, T_max, and WS. Additionally, annual ET_o exhibited a higher value in southern Punjab in comparison to northern Punjab, with a range of 2975 mm/year in the cotton–wheat zone to 1596 mm/year in the rain-fed zone. Over the course of the past seventy years, Punjab experienced an average increasing slope of 5.18 mm/year in ET_o. T_min was the highest monthly dominant factor throughout the year, whereas WS and SH were recorded to be the dominant factor in the winters, specifically. All in all, accurate estimation of ET_o, which serves as an essential component for crop water requirement, could potentially help improve the irrigation scheduling of crops in the agroclimatic zones.

1. Introduction

The determination of reference evapotranspiration (ET_o) holds significant importance in the realm of remote sensing applications and the regulation of agricultural water usage. In crop water management, it is customary to employ a theoretical grass reference crop possessing specific attributes as a benchmark surface for the purpose of approximating ET_o [1]. The estimation of ET_o holds significant value in the management of hydrological and agricultural water systems. Nonetheless, it is imperative to conduct a thorough evaluation of the quantity of water necessary for the growth of any given crop [1]. Accurately determining ET_o across different land use regions is crucial for efficient management of water resources [2]. Methods for transferring water vapor, such as soil water balance [3,4] and remote sensing techniques [5], have been employed to quantify ET_o. In addition, several hybrid techniques have been used for forecasting ET_o [6,7].The accurate computation of ET_o is a complicated undertaking that entails the consideration of multiple meteorological variables. Several regions around the world have reported the spatiotemporal variation in the climatic parameters, mainly including temperature [8], which ultimately impacts ET_o. Consequently, the development of a precise empirical model that can effectively account for the complexities of the process is a challenging task, as noted in reference [9]. Additionally, climate change plays a significant role in disturbing the surrounding environment in different ways [10], including depleting water levels [11,12], a disturbed hydrological cycle which includes ET_o [13,14] as an important component, and ultimately economic loss [15].

A variety of interpolation techniques are employed in GIS software, each distinguished by the models utilized to forecast unknown values on the basis of sample points. The inverse distance weighted (IDW) and Kriging techniques are among the most frequently employed methods. The inverse distance weighting (IDW) method is a deterministic interpolation technique that relies on nearby reference point measurements or established mathematical formulae to evaluate the texture of the resulting interpolation area [16]. Kriging is a geostatistical method that incorporates the spatial configuration of sample points by assigning weights that differ according to their respective locations. Kriging methodology not only generates an interpolated or predicted surface but also provides a level of precision or confidence in the estimates, which is in contrast to the IDW interpolation technique. The spatiotemporal estimation of ET_o across 41 stations in Pakistan was investigated by Shah et al. [17] during the period of 2006 to 2015. The research discovered that ET_o demonstrated variability and manifested both positive and negative trends across various seasons in multiple regions of Pakistan. Geographically weighted regression (GWR) could potentially serve as a more appropriate technique for addressing this issue with greater precision. The authors Fotheringham et al. (2017) [18] presented a method called multiscale geographically weighted regression (MGWR) to tackle the issues arising from nonstationary connections and heterogeneity. This particular model has the ability to precisely replicate the effects of different independent variables on a dependent variable with regards to spatial variability.

The MGWR method has garnered considerable interest from diverse fields owing to its efficacy in modeling the spatial associations among multiple dependent and independent variables [18]. The effectiveness of ET_o estimation models is inherently connected to the specific meteorological conditions of a given location, as stated in reference [19]. Therefore, it is crucial to ascertain appropriate models that are contingent on the accessibility of meteorological data. The significance of ET_o on hydrological processes and water resources is particularly noteworthy in predominantly arid regions of Pakistan, surpassing its impact in other climatic zones. The exacerbation of droughts and aridity presents a formidable challenge to the sustainable development of agriculture and water resources in Pakistan, a nation heavily reliant on agriculture. Accurate estimation of ET_o is essential for efficient monitoring and management of droughts and aridity. The inadequacy of meteorological information has posed a difficulty in accurately determining ET_o across various regions of Pakistan. Therefore, it is crucial to investigate suitable ET_o models that can be suggested for ET_o estimation in diverse areas. The weather in Pakistan varies between different seasons on a monthly level, which makes it difficult to understand the influence of climatic parameters and the dominant factor influencing ET_o in each month. One of the major difficulties in water resources management is the accurate estimation and forecasting of ET_o. However, the most recent estimation of ET_o and the dominant factor influencing it on a monthly level are missing, which negatively impacts water resources management ultimately impacting the end-user.

According to the studied literature, it can be safely assumed that up-to-date estimations of the reference evapotranspiration and, especially, of the monthly influence of a dominant factor (out of minimum temperature, maximum temperature, relative humidity, windspeed, and sunshine hours) on ET_o for Punjab, Pakistan, are missing. Therefore, this study investigates the spatiotemporal variation of reference evapotranspiration across Punjab, Pakistan. Furthermore, spatial variation in the monthly impact of influencing parameters on ET_o and the dominant factor have been estimated.

2. Materials and Methods

2.1. Study Area

The selection of Punjab, Pakistanm for the analysis of the impact of influencing parameters on ET_o is justified by the fact that agriculture constitutes a significant proportion of the country’s overall exports [20], with Punjab contributing to 60% of this portion [20], as depicted in Figure 1. Moreover, Punjab has attempted to tackle the problems of food security for many years; therefore, the assessment of the influence of climatic conditions on crop production in such an agriculture-rich region is imperative. Furthermore, accurate ET_o estimation in this region can aid in decision making about irrigation scheduling, water allocation for different crops, and the scheduling of other agricultural activities [21]. Moreover, this study assesses the behavior of ET_o in five different agroclimatic zones of Punjab, with 36 selected study points, i.e., districts of the province. Different approaches can be used for the classification of Punjab districts into regions, such as dividing Punjab into three regions, north, central, and south, or disintegrating rural Punjab into north and south. In ref. [22], Amjad et al. used the methodology of Pinckney [23] to classify rural areas of Punjab province into five agroclimatic or crop zones. The agroclimatic zones of the Punjab were classified according to the published report by the International Food Policy Research Institute.

The categorization of various zones is primarily determined by the primary kharif season crop, as wheat dominates the rabi season in most regions. The kharif season in Pakistan encompasses the time frame from April to September, during which monsoon rains are prevalent. In contrast, the rabi season occurs from October to March, featuring cooler temperatures and reduced precipitation. Agricultural activities during the rabi season heavily rely on irrigation methods for crop cultivation. In regions where irrigation is practiced, rice and cotton are the predominant kharif crops of utmost importance. Rice cultivation exhibits a propensity for prevalence in regions characterized by elevated water tables, dense soils, and copious precipitation, whereas cotton production is commonly observed in arid locales. Therefore, a significant differentiation can be observed between regions that are suitable for rice cultivation and those that are suitable for cotton cultivation. Based on these two crops, i.e., cotton and rice, alternating with the dominant crop, wheat, two zones, namely, the cotton–wheat zone and rice–wheat zone, were created. The cotton–wheat zone, located in the southern and southwestern regions of Punjab, exhibits a climatic pattern characterized by hot summers and cool winters. The precipitation in this region is comparatively lower than that in the rice–wheat zone, rendering it conducive for the cultivation of cotton during the kharif season and wheat during the rabi season. The rice–wheat zone, predominantly situated in the central and eastern regions of Punjab, exhibits a climate distinguished by warm summers and mild winters. The monsoon season is characterized by a substantial precipitation pattern that occurs during the kharif season, creating favorable conditions for the cultivation of rice. In contrast, the rabi season is characterized by lower levels of precipitation, making it conducive for the cultivation of wheat.

A distinct geographical region, with Faisalabad at its center, exhibits favorable conditions for cultivating multiple kharif crops, without any single crop exerting dominance. This region has been categorized as the mixed-crop zone. The region exhibits a confluence of climatic conditions observed in different zones.

The rain-fed zone or the barani region is deemed distinctive owing to its dependence on precipitation for agricultural purposes. The rain-fed zone situated in the northern and northwestern regions of Punjab exhibits a climate characterized by higher aridity. The precipitation levels in this particular region are constrained.

Furthermore, the region adjacent to the left bank of the Indus River in Punjab, characterized by limited irrigation infrastructure and lower levels of agricultural activity, is recognized as a distinct area referred to as the low-intensity zone. The climate in this zone is characterized by aridity, resulting in minimal precipitation. It is imperative to acknowledge that the aforementioned descriptions offer a comprehensive outline of the climatic conditions prevalent in the agroclimatic zones of Punjab. However, it is crucial to recognize that there may be deviations within each zone due to localized factors, including but not limited to elevation, proximity to water bodies, and microclimatic circumstances. The delineation of the five agroclimatic zones, namely, the low-intensity zone, the rice–wheat zone, the mixed-crop zone, the cotton–wheat zone, and the rain-fed zone, is attributed to the variations in growing climatic conditions and agricultural practices observed from north to south.

Table 1. Summary of the climatic parameters (1950–2021) for the five agroclimatic zones in Punjab, Pakistan.

Agroclimatic Zone	Maximum Temperature Tmax (°C)	Minimum Temperature Tmin (°C)	Relative Humidity RH (%)	Sunshine Hours SH (h)	Windspeed WS (m/s)	Reference Evapotranspiration ETo (mm/day)
Cotton–wheat zone	32.97	18.44	40.4	11.99	2.85	6.49
Rice–wheat zone	30.47	17.35	67.49	11.99	3.47	5.47
Low-intensity zone	32.35	17.86	44.99	11.99	3.29	6.52
Mixed-crop zone	31.8	18.02	57.37	11.99	3.62	6.15
Rain-fed zone	28.89	15.45	74.17	12	4.04	5.03

presents the average values of the five climatic parameters influencing ET_o for the period 1950–2021 for the five agroclimatic zones.

The low-intensity zone includes the districts Bhakkar, Dera Ghazi Khan, Layyah, Mianwali, Muzaffargarh, and Rajanpur, with 32.3 °C, 17.8 °C, 44.9%, 11.9 h, 3.2 m/s, and 6.5 mm/day of T_max, T_min, RH, SH, WS, and ET_o, respectively, with an average for the period of 1950–2021. The rice–wheat zone includes districts Gujranwala, Gujrat, Hafizabad, Kasur, Lahore, Mandi Bahauddin, Nankana Sahib, Narowal, Sheikhupura, and Sialkot, with 30.4 °C, 17.3 °C, 67.4%, 11.9 h, 3.4 m/s, and 5.4 mm/day of T_max, T_min, RH, SH, WS, and ET_o, respectively. The mixed-crop zone includes the districts Chiniot, Faisalabad, Jhang, Khushab, Okara, Sargodha, and Toba Tek Singh, with 31.8 °C, 18.0 °C, 57.3%, 11.9 h, 3.6 m/s, and 6.1 mm/day of T_max, T_min, RH, SH, WS, and ET_o, respectively. The cotton–wheat zone includes the districts Bahawalnagar, Bahawalpur, Khanewal, Lodhran, Multan, Pakpattan, Rahim Yar Khan, Sahiwal, and Vehari, with 32.9 °C, 18.4 °C, 40.4%, 11.9 h, 2.8 m/s, and 6.4 mm/day of T_max, T_min, RH, SH, WS, and ET_o, respectively. The rain-fed zone includes the districts Attock, Chakwal, Jhelum, and Rawalpindi, with 28.8 °C, 15.4 °C, 74.1%, 12.0 h, 4.0 m/s, and 5.0 mm/day of T_max, T_min, RH, SH, WS, and ET_o, respectively.

2.2. Data Collection

For the estimation of ET_o, satellite data were used. The data of 5 different climatic parameters were obtained from 2 online sources (summarized details of the data sources are presented in

Table 2. Summary of the data sources used in this study.

Dataset Name	Variables	Temporal Extent	Spatial Resolution	Ref.
CRU–TS	Tmin, Tmax	1901–2021	All land areas (excluding Antarctica) at 0.5° resolution	[24]
NCEP–NCAR Reanalysis 1 (Data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov (accessed on 28 January 2023))	RH and WS	1948–2021	All land areas (excluding Antarctica) at 2.5° resolution	[25]
Sunshine hours (calculated)	SH	–	–	–

). For minimum temperature (T_min) and maximum temperature (T_max), the satellite data were obtained from the Climatic Research Unit gridded Time Series (CRU–TS). However, for relative humidity (RH) and windspeed (WS), the satellite data were obtained from National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) (NCEP–NCAR) Reanalysis 1 for Punjab. For further analyses, data points were extracted for the 36 districts (categorized into 5 agroclimatic zones as per [22,23]) in the Punjab from 1950 to 2021. Figure 2 illustrates the workflow of the presented study. The present study assessed the decadal patterns of ET_o spanning from 1950 to 2021. In order to analyze the trend of reference evapotranspiration, the available data were divided into seven distinct decades spanning 1950 to 2021. Additionally, the data were further categorized into five distinct zones, namely, the cotton–wheat zone, the low-intensity zone, the rice–wheat zone, the mixed-crop zone, and the rain-fed zone, all of which are located within the Punjab province of Pakistan.

2.2.1. ET_o Estimation

In this study, the standard Penman–Monteith (PM) model was utilized to calculate the reference evapotranspiration in mm/day. It was developed in the mid-1960s by Allen and Tideman, and it combines both energy and mass balance principles [26]. The PM ET_o model is given as Equation (1). Details of the parameters are available in the literature [26].

$$\mathrm{E}{\mathrm{T}}_{\mathrm{o}}=\frac{0.408∆\left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}\right)+\mathsf{\gamma }\frac{900}{\mathrm{T}+273}{\mathrm{u}}_2({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}})}{∆+\mathsf{\gamma }(1+0.34{\mathrm{u}}_2)}$$

(1)

where ET_o is the reference evapotranspiration in mm day⁻¹; Δ is the slope of the vapor pressure curve in kPa °C ⁻¹; R_n is the net radiation at the crop surface in MJ m² day⁻¹; G is the soil heat flux density in MJm².day⁻¹, which is considered as 0 [1]; γ is the psychrometric constant in kPa °C ⁻¹; T is the mean air temperature, which is determined from minimum and maximum temperatures in °C; u₂ is the wind speed at a height of 2 m height in ms⁻¹; e_s is the saturation vapor pressure in kPa; and e_a is the actual vapor pressure in kPa.

2.2.2. Trend Analysis

The present investigation employed a prevalent technique for trend analysis, namely, linear regression, due to its ability to detect and measure patterns in data across temporal dimensions. Through the application of a linear model to the dataset, researchers are able to approximate the pace and orientation of variation in a given variable that is of significance. The aforementioned approach proves to be especially advantageous in the examination of extended patterns in climatic parameters, such as ET_o. In this study, the Mann–Kendall (MK) test was employed to identify fluctuations in the timeseries data [27]. The Sen’s slope technique has been utilized in combination with the MK test to examine patterns in ET_o data. The technique, formulated by Sen in 1968, offers a computation of the pace of alteration in a given variable across a period [28]. This study employs linear trend analysis to examine the Punjab region, utilizing the linear regression method for the period spanning 1950 to 2021.

2.2.3. Data Quality and Integrity Assessment

In the current scenario, data were analyzed through basic statistical methods such as minimum and maximum values, range, mean, standard deviation, skewness, and kurtosis.

Table 3. Descriptive statistics of the climatic parameters for the estimation of ET_o for Punjab, Pakistan, for the period of 1950–2021.

Variable	Mean	St. Dev.	Minimum	Maximum	Range	Skewness	Kurtosis
Tmin (°C)	17.631	0.516	16.426	18.774	2.348	–0.14	–0.47
Tmax (°C)	31.496	0.499	30.168	32.412	2.243	–0.25	–0.46
RH (%)	55.747	2.529	49.407	61.944	12.537	–0.33	0.02
WS (m/s)	3.3819	0.5137	2.3519	4.3133	1.9615	–0.32	–1.19
ETo (mm/year)	2339.7	149.1	2094.1	2691.3	597.1	0.63	–0.34

presents an overview of the descriptive statistics of the climatic parameters for the estimation of ET_o for Punjab, Pakistan, for the period 1950–2021. The analysis conducted and presented in Table 3 demonstrates that the studied variables exhibited varying skewness values, encompassing both positive and negative values. This suggests the presence of asymmetry in the distributions of these variables. Moreover, the kurtosis values indicate that the climatic variables demonstrate either a distribution with lower peak and heavier tails (platykurtic) or a distribution with higher peak and lighter tails (leptokurtic).

2.2.4. Spatial Distribution of ET_o

The utilization of spatial visualization techniques for ET_o distribution analysis facilitates the identification of regions with elevated water demand and potential depletion or overuse of water resources [1]. Furthermore, it has the potential to identify regions with favorable temperature [29] and humidity conditions for cultivating various crops, thereby facilitating crop selection and land-use management. The utilization of Inverse Distance Weighting (IDW) for the purpose of interpolating ET_o data and creating spatial representations of ET_o across various regions of Pakistan has been extensively documented in numerous studies [30]. As a result, the ET_o data underwent spatial interpolation using the IDW interpolation technique in the ArcGIS software, as opposed to more intricate methods such as Kriging. This process yielded evapotranspiration maps for various climatic periods spanning from 1950 to 2021, which are presented in Figure 3 for the purpose of change analysis.

2.2.5. Spatial Relationship between Climatic Variables and ET_o

The present study employed Pearson’s correlation coefficient to establish the correlation between the influencing parameters and the reference evapotranspiration. This coefficient was preferred over other methods of evaluating linear relationships due to its distinct advantages. The utilization of Pearson’s correlation coefficient (R) is widespread in the assessment of models in studies related to ET_o, as evidenced by several studies [31]. The requisite quantity of data for establishing a correlation is comparatively lower than in other techniques of data assessment, typically ranging from three to seven variable pairs [32]. Furthermore, it has the capability to conduct an analysis on a considerable number of variables without causing substantial overstatement of outcomes [32]. A sensitivity analysis was carried out to assess the impact of ambient weather variables on the computed ET_o. Bois et al. [33] have discussed several techniques for performing sensitivity analyses, which have been documented in the literature. In this study, a monthly analysis was performed to establish a correlation between the influencing parameters and ET_o. The mathematical expression denoting the Pearson correlation is represented by Equation (2):

$$r=\frac{\sum _{i=1}^N(x-\stackrel{-}{x})(y-\stackrel{-}{y})}{\sqrt{\sum _{i=1}^N(x-\stackrel{-}{x}{)}^2\sum _{i=1}^N(y-\stackrel{-}{y}{)}^2}}$$

(2)

In this study, the observed climatic variable was denoted as x and the estimated ET_o was represented as y; $\stackrel{-}{x}$ and $\stackrel{-}{y}$ are their subsequent mean values. Similarly, the coefficients were used to determine the most influential climatic variables on ET_o in the study areas.

3. Results

3.1. Trend of ET_o

Figure 4 depicts the results of the Mann–Kendall trend analysis applied to annual ET_o in the five agroclimatic zones of Punjab for analyzing the heterogeneity of the timeseries ET_o. In Figure 4a, a heterogeneity in the year 2000 was observed in ET_o for Punjab during the period of 1950 to 2021. The sudden fluctuation in the trend of ET_o could be attributed to the changing climate, which ultimately impacts the climatic parameters which influence ET_o. Figure 4b shows cumulative trend of ET_o for the period of 1950–2021, which is increasing through the time period. A lowest ET_o of 2094.12 mm/year was observed in the year 1963 for the Punjab, whereas the highest ET_o, i.e., 2691.26 mm/year, was observed in 2021, which clearly illustrates the impact of rising temperatures. Similarly, trends of the ET_o for the five agroclimatic zones are presented in Figure 4; all of the zones show rapidly increasing trends of ET_o except in the low-intensity zone, which could be attributed to lower precipitation in the area. The agroclimatic zones comprising the respective districts are also presented in Figure 4.

3.2. Decennial ET_o

The data presented in Figure 5 indicate a decadal pattern that demonstrates a rise in ET_o throughout the duration of the study. The cartographic representation of the 1950s (depicted in Figure 5a) revealed a diminished ET_o in the rain-fed and rice–wheat zones (i.e., 1328.03 mm/year), whereas the southwestern region exhibited relatively elevated values (i.e., 2370.56 mm/year). Figure 5b illustrates the fluctuation of difference in ET_o during the second decade (1960–1969) throughout Punjab with respect to the ET_o base year 1950s. The study noted that the highest difference in the ET_o values across the northern and southern regions were found in the cotton–wheat zone, with Rahim Yar Khan particularly exhibiting a 10.19 mm/year elevated ET_o compared to the base year. The mixed crop zone, also known as the central region, typically exhibited a moderate-to-low level of ET_o. In majority of the cotton–wheat zone, the difference ET_o was negative (up to −91.31 mm/year) compared to the base year (illustrated in Figure 5b). The results obtained during the time frame of 1970–1979 indicate a rise in the difference of ET_o in the rice–wheat zone and rain-fed zone (i.e., up to 116.26 mm/year) when compared to the base year. This is illustrated in Figure 5c. During the fourth decade spanning 1980 to 1989, there was an observed increase in the difference of ET_o in comparison to the base year 1950s, as illustrated in Figure 5d. The highest positive difference in observed ET_o was 84.77 mm/year in the rice–wheat zone, rain-fed zone, and somewhat in the mixed-crop zone compared to the base year. As depicted in Figure 5e, the rate of difference in ET_o exhibited an increasing trend from the northern region to the southern region (i.e., from 115.74 mm/year to −79.35 mm/year, respectively), or from the rain-fed zone to the cotton–wheat and low-intensity zone during the period 1990–1999. Drastically elevated trends in the difference of ET_o were observed in the year 2000–2009 compared to the base year 1950s, as shown in Figure 5f. Typically, rain-fed zone was the most affected zone (i.e., 234.73 mm/year), followed by the rice–wheat zone (i.e., 200.8 mm/year), compared to the base year 1950s. According to Figure 5g, the highest difference in ET_o was observed in the northern region (i.e., 734.57 mm/year) compared to the southern region (320.29 mm/year) in the year 2010–2021 relative to the base year 1950s. Finally, Figure 5h illustrates the ET_o levels in the period of 2010–2021. i.e., 2690.85 mm/year in the southern region of the area, whereas 2042.86 mm/year was recorded in the northern regions of the area. Over the course of the past decade, the cotton–wheat and low-intensity zones have exhibited the highest ET_o values, whereas the mixed-crop zone has demonstrated a moderate rate (i.e., 2400 mm/year to 2500 mm/year) with slight anomaly in the Jhang district. Conversely, the rice–wheat zone in the northeastern region and the rain-fed zone in the northwestern region have displayed relatively low ET_o rates; however, they were the most affected zones when compared to the base year 1950s due to climate change. The rate of reference evapotranspiration is directly proportional to temperature, resulting in a significant fluctuation in ET_o rate during the month of September as supported by the findings of previous research [34]. As per the data presented in Figure 5, the highest evapotranspiration rate recorded was 2690.8 mm/year in the year 2010–2021. Prior research has identified comparable patterns that align with the results of the current investigation.

Table 4. Decennial trend of ET_o against each of the five agroclimatic zones from 1950 to 2021 in Punjab, Pakistan.

Agroclimatic Zone	ETo (mm/year) per Decade							ETo Trend (mm/year)
	1950–1950	1960–1969	1970–1979	1980–1989	1990–1999	2000–2009	2010–2021	1950–2021
Low-intensity zone	2243	2230	2221	2220	2287	2227	2237
Cotton–wheat zone	2157	2096	2068	2072	2123	2293	2323
Rice–wheat zone	1999	1990	2100	2097	2105	2070	2046
Rain-fed zone	1801	1805	1904	1850	1873	1956	1965
Mixed-crop zone	2064	2044	2117	2114	2173	2164	2149

presents the decennial pattern of ET_o in relation to the five agroclimatic zones in Punjab, Pakistan, spanning from 1950 to 2021.

3.3. Sensitivity Analysis

3.3.1. Geographic Weighted Regression (GWR)

In order to examine the relationship between the distribution of ET_o and its influencing factors, a correlation analysis was conducted using Geographically Weighted Regression (GWR) at 36 district points. Subsequently, maps were generated using Inverse Distance Weighting (IDW) gridding. Figure 6 presents the resultant maps that illustrate the correlation between ET_o and various climatic parameters for the time period 1950–2021. A correlation ranging from 0.29 to 0.89 with standard deviation of 0.1 was observed between ET_o and T_min reported in Figure 6a. According to Figure 6a, low correlation (i.e., 0.29) was observed in the southern region of the area covering most of the cotton–wheat zone and the low-intensity zone. This could be attributed to the gradual change in elevation from north to south, which is known to impact ET_o. In Figure 6b, a good correlation ranging 0.78 to 0.91 with standard deviation of 0.03 was observed between ET_o and T_max. According to Figure 6b, the correlation trend follows a similar direction as that of ET_o and T_min; however, the minimum value of the correlation is raised to 0.78, covering most of the cotton–wheat zone, the low-intensity zone, and the mixed-crop zone. This could be attributed to the changing climate and rising temperatures across the country. A correlation ranging from 0.67 to 0.89 with standard deviation of 0.05 was observed between ET_o and RH reported in Figure 6c. A minimum correlation of 0.67 was observed between ET_o and RH in the central region of the mixed-crop zone, whereas a correlation of 0.7 to 0.8 was observed between ET_o and RH in most of the cotton–wheat zone, the low-intensity zone, and the mixed-crop zone. In the rain-fed and rice–wheat zones, a correlation of 0.8 to 0.89 was observed between ET_o and RH. This trend could be attributed to the frequent rainfall events in the northern region of the area leading to increased RH compared to the infrequent rainfall events in the southern region leading to drier conditions and lesser RH. According to Figure 6d, a correlation ranging from 0.29 to 0.63, with standard deviation 0.06, was observed between ET_o and WS. Minimum correlation, i.e., 0.29, was observed in the southern region of the cotton–wheat zone and the low-intensity zone due to lesser windspeed in these southern regions of the area. A moderate correlation of 0.6 was observed between ET_o and WS in the central Punjab region, comprising some districts of the cotton–wheat, low-intensity, and rice–wheat zones. This could be attributed to typically higher windspeeds in these regions due to their geographical location and the location of the Himalayan ranges. In the northern region of the area, namely, the rain-fed zone and the rice–wheat zone, a correlation of 0.3 to 0.4 was observed between ET_o and WS.

3.3.2. Multiscale Geographic Weighted Regression (MGWR)

The MGWR model was utilized to estimate ET_o in Punjab from 1950–2021, and Figure 7 displays the monthly distribution of the relative contribution of minimum temperature (T_min) in this estimation. The findings presented in Figure 7a demonstrate the varying levels of contribution of T_min during the month of January across different regions, ranging from 0.17% to 0.75%. The study revealed that ET_o was significantly impacted by the rain-fed and low-intensity zones, whereas the rice–wheat and mixed-crop zones exhibited a relatively lower level of influence. Figure 7b demonstrated a comparable pattern where the minimum temperature (T_min) had an impact on ET_o during February. The range of this impact was between 0.18% and 0.76%, indicating a significant level of contribution in both the rain-fed and cotton–wheat zones. During the month of March, the contribution of T_min varied between 0.15 to 0.74%. The rain-fed zone was observed to be the most affected, as depicted in Figure 7c. The impact of T_min during the month of April is depicted in Figure 7d, revealing a significant influence on the rain-fed, low-intensity, and cotton–wheat zones, with a range of 0.17% to 0.71%. The rice–wheat zone, rain-fed zone, low-intensity zone, and mixed-crop zone experienced significant impact in May, with an overall range of 0.22% to 0.76%, as illustrated in Figure 7e. During the month of June, the T_min had an impact on the value, which varied between 0.31% and 0.85%. The rice–wheat zone, low-intensity zone, and rain-fed zone were the areas that were most affected, as shown in Figure 7f. The data presented in Figure 7g indicate that the values observed for July ranged between 0.25% and 0.80%. The rice–wheat zone and mixed-crop zones exhibited comparatively lower levels of influence, whereas the remaining zones were subject to significant impact. According to Figure 7h, the contribution of T_min during the month of August varied across different regions, with values ranging from 0.32% to 0.85%. The study revealed that the rain-fed zone had a significant impact on ET_o, whereas the remaining zones exhibited a comparatively lower level of contribution. Figure 7i illustrates the impact of T_min on ET_o for the month of September. The observed values exhibited a significant level of contribution within the rain-fed zone, ranging from 0.27% to 0.86%. The rain-fed zone in October exhibited a similar pattern, with the highest impact ranging from 0.22% to 0.72%, as shown in Figure 7j. The influence of T_min on ET_o was observed in November, with values ranging from 0.16% to 0.72%, as depicted in Figure 7k. The rice–wheat zone, cotton–wheat zone, low-intensity zone, and mixed-crop zone were found to have experienced the least amount of impact. Finally, it should be noted that the impact of T_min on ET_o in the rain-fed region remained consistently significant throughout the month of December, exhibiting the greatest level of influence. The observed values exhibited a range between 0.16% and 0.75%, with the mixed-crop zone displaying the least amount of impact, as depicted in Figure 7l.

The findings derived from Figure 7 demonstrate that the impact of climatic variables on ET_o exhibits variability on a monthly basis. Figure 8 depicts the dominant factor that exerted the highest influence on ET_o on a monthly basis from 1950 to 2021 in the region of Punjab, Pakistan.

According to Figure 8, the dominant factor varied across space and time. In terms of monthly variation, WS was the dominant factor across most of the Punjab in January, as shown in Figure 8a. However, in February, SH was observed to be the dominant factor across majority of the 36 districts in Punjab, as shown in Figure 8b. Similar to January, WS dominated the influence on ET_o across Punjab in March, as shown in Figure 8c. According to Figure 8d, the WS and SH were the most prominent dominant factors out of the five climatic parameters influencing ET_o in April across Punjab. In May, T_min, SH, and WS were the most dominant factors spread across Punjab (shown in Figure 8e). This could be linked to longer solar days, higher wind speed across the plains, and temperature variation due to start of summer months. In June, T_min was the most dominant factor influencing ET_o across Punjab, shown in Figure 8f. Similarly, in July, T_min and T_max were the driving factors influencing reference evapotranspiration across Punjab as shown in Figure 8g. In August, T_min and WS were the most dominant factors shown in Figure 8h. According to Figure 8i, in September, T_min, T_max, and WS were the most dominant factors influencing different agroclimatic zones of the Punjab. According to Figure 8j, in October, in the low intensity zone, SH was the driving factor, whereas in the upper Punjab, RH was the driving factor influencing ET_o. In other zones, T_min and SH were observed to be influencing ET_o. In November, a mix–up of T_min, T_max, WS, and SH influenced ET_o at different stations across Punjab, as shown in Figure 8k; however, WS and SH were the most common dominant factors at majority of the stations. According to Figure 8l, WS and SH were the dominant factors impacting ET_o across Punjab in December.

4. Discussion

T_min series showed more significant trends, which is in accordance with the global trend over the previous decades. The overall pattern assessment revealed that annual T_min had a considerable influence on ET_o in the rain-fed zone. ET_o is impacted in the heavy rainfall [35,36] seasons due to variations in relative humidity. To address the varying influence, the GWR technique was used, which can be effective in allowing model parameters to vary across space [18]. This approach can provide valuable insights into the complex spatial patterns of ET_o and its influencing factors, which can aid in developing more accurate and localized water management strategies. The correlation of climatic factors to ET_o findings demonstrates that all parameters have a substantial influence on ET_o [37]. According to the findings by Hupet et al. [38], it has been observed that radiation has a significant influence on ET_o. Since minimum temperature remains unchanged as much during the wet season as it does during the period of dryness, it has a greater influence on ET_o during the dry season. Therefore, the impact of T_min was observed on ET_o via MGWR and its subsequent influence from the month of January to December for the study period 1950–2021. According to Zhang et al. [39], ET_o was highly influenced by T_min and T_max across China on an annual basis. The study provided a detailed dominant factor analysis throughout the study period (1970–2014) in China. However, the study failed to correlate the impact of monthly variation of weather conditions/seasons on the influencing ability of the climatic parameters. The present study assesses this impact on a monthly basis. Furthermore, the estimated reference evapotranspiration could be correlated with measured ET_o from any of the stations for reference. Additionally, in future research directions, the estimated spatiotemporal ET_o (which serves as an essential component in the crop water requirement estimation) could help in the estimation of zone- and crop-specific irrigation scheduling across Punjab. Ultimately, based on the study, the estimation of the amount of water required for irrigation purpose would help the relevant policymakers to better allocate the water resources in the country.

5. Conclusions

Maximum T_max, T_min, RH, WS, and ET_o (i.e., 33.9 °C, 18.7 °C, 76.5%, 4.2 m/s, and 7.6 mm/day, respectively) were observed across the study area. The highest ET_o was observed in the southern regions of the area. The range of ET_o was also higher in the southern region (i.e., up to 5.57 mm/day). As per GWR, T_max had the highest correlation (i.e., 0.91) with ET_o in the mixed-crop zone. As per the MGWR, T_min, T_max, and WS had the highest influence on ET_o throughout different regions in Punjab. Throughout the study period, ET_o showed an increasing trend. In addition, a heterogeneity was observed in the ET_o trend in year 2001. Out of the five agroclimatic zones, the cotton–wheat zone observed the highest ET_o due to the higher temperatures in the region. In terms of a monthly dominant factor, the highest overall dominant factor was T_min, whereas WS and SH were the highest observed dominant factors in the winters, specifically. It was concluded that T_min and T_max had the highest overall influence on ET_o in Punjab. In addition to these parameters, WS was also observed to be the highest dominant factor influencing ET_o on a monthly basis. In terms of future research directions, the ultimate goal of the proposed study was to investigate the zone- and crop-specific irrigation scheduling which incorporates ET_o as an essential component.