# Spatial-Temporal Interpolation of Reference Evapotranspiration for Pakistan

Authors: Syed Muhammad Saeed Shah; M. El-Morshedy; Wahidullah Mansoor
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The unequal water distribution in the universe has resulted in more than 2 billion people living in water-stressed areas. Globally, withdrawals of water resources are still below the critical level. Pakistan is also affected by this serious problem. In the context of environmental evaluation, irrigation scheduling, and water resource management, evapotranspiration predictions are crucial. The evapotranspiration estimation methods may need to be evaluated daily or monthly to understand the climate change effects in the local areas. This paper investigates spatial-temporal interpolation of evapotranspiration data from 41 Pakistani meteorological stations. To estimate evapotranspiration, we have used the average time series data between 2006 and 2015 on temperature, relative humidity, wind speed, and solar radiation. We developed a new modified Hargreaves equation to estimate evapotranspiration and evaluate the performance of different models. We found that our modified Hargreaves model could perform better than the original Hargreaves model. To analyze trends in different seasons, such as rainy, dry, and annual, FAO-56 Penman-Monteith, Blaney-Criddle, and Hargreaves-Samani models are used. We used Ordinary Kriging for Functional Data to map and predict evapotranspiration spatially and temporally in various locations.


---

## Body

## 1. Introduction

It is necessary to mention that evaporation and evapotranspiration are two different concepts. The concept of evapotranspiration (ET0) refers to the elimination of liquid water from a vegetation-rich region and its transport to the atmosphere through transpiration and evaporation. The evaporation process involves converting liquid water into water vapor and evaporating from surfaces such as wet vegetation, lakes, and soil into the atmosphere. In transpiration, water in the plant, tissues are eliminated to the atmosphere, most significantly through stomata, which are the small opening in the leaves of plants and grasses.It is undeniable that human beings are experiencing a water crisis due to rapid population growth, rising temperatures, droughts, and land depletion. Water shortage questions center on how much water we use in daily life. In terms of water usage, agriculture consumes 70 percent of it, as outlined by Montgomery and Elimelech [1].Burn and Hesch [2] determined the changes in precipitation, ET0, and water deficits due to public involvement with climate change.Donatelli et al. [3] showed that accurate ET0 prediction for agricultural reasons is crucial to meet the challenges of effective water management.Using data from 41 Australian sites from 1975 to 2004, Roderick et al. [4] examined wind speed, humidity, temperature, and radiation and concluded that humidity and temperature establishment trends have little influence on the pan evaporation patterns. A study of ET0 provides a basis for agricultural management and planning, especially in Mediterranean-controlled ecosystems, where rainfall changes with temperatures and irrigation are someway accomplished by this water supply; see Bates et al. [5].Generally, the direct impact of weather change on water resources has been studied throughET0. Hydrological changes in tropical areas have significant impacts on weather change. Consequently, climate change will lead to an increase in temperature and variations in rainfall prototypes. A higher temperature will stimulate the hydrological system and water assets, causing an increase in ET0; see Shahid [6].The irrigation system of Pakistan is the most significant among the globe. Hence, an appropriate assessment ofET0 and its possible variation due to the weather change is very critical for the scheduling and management of water assets.We chose four methods for predicting the referenceET0 in this paper. First, a Penman-Monteith is often used when climatological information is available or may be reliably predicted from the empirical equations. A second method is Hargreaves, which is a temperature and radiation-based process that can be applied when climatological data is inadequate. It estimates reference ET0 on time intervals equal to or longer to one month. The third procedure considered is BCR. It is a temperature-based method that gives a prediction of ET0 on a monthly basis. The use of this method has been questioned, but due to its widespread application around the globe, we use it to determine climate change. The fourth procedure, Modified Hargreaves, is a temperature and radiation-based technique that contains much empiricism in its present form.Many researchers have examined trends in climate variation to determine whether climate variation is significant; see Kousari and Ahani [7]. Among the four weather variables, humidity is the most significant factor influenced by solar radiation, wind speed, temperature, and the order of the four changes with region and season, affecting wind-speed patterns at the regional level. Moreover, one can not deny the importance of the quality of weather data and estimation techniques for climatic changes [8].There are a wide variety of techniques available for analyzing the geo-statistical data; see Cressie and Wikle [9]. Many software macros and packages have been developed in order to carry out Spatio-temporal modeling. For example, spate package gives functionality by Stochastic Potential Difference Equation (SPDE); see Sigrist et al. [10]. The RandomFields package is appropriate for analyzing and simulating spatial and spatiotemporal random fields; see Schlather et al. [11]. ComPRandFld gave a set of procedures for the analysis of spatiotemporal random fields and suggested most spatial-temporal covariance models, the pairwise composite likelihood for estimation, and spatiotemporal tapering for kriging. Bourotte et al. [12] developed a package that provides functionality for analyzing and estimating multivariate spatial-temporal records. In this respect, this package gives an edible class for nonseparable models and provides an interesting starting point.Space-time kriging relies on the data set and can be computationally expensive. Hence, local kriging is an attractive substitute. It raises the question of choosing the “closest” neighbors from the space-time space. An obvious choice in such a situation would be to choose the space-time regions possessing the highest correlation with the unobserved regions. The relationship between spatial and temporal distance in finding out the extent of correlation will vary, depending on the space-time covariance model. We also considered a specialR package geofd for spatially distributed data with respect to time for predicting the functional data [13]. The k-nearest neighbors in this metric spatial-temporal space S×T are chosen to use the R package FNN [14].The difference between techniques based on the actual data is significant [15]. However, from the estimation point of view, the technique based on ordinary kriging for function-valued spatial data is the simplest one. Further, the estimation of linear coregionalization models for modeling the spatial dependence among curves is significant, conditional upon the number of fundamental functions used for soothing the more extensive set of data. In this context, R package geofd is applied within the statistical limits, which only include functions for carrying out the spatial forecast of functional data by “ordinary kriging” for function-valued spatial data; see [16].This paper considers different evapotranspiration methods to analyze a long-term time series of weather data from various weather stations located in Pakistan from 2006 to 2015. The aim of this paper is to model and characterize monthly, seasonal, and annual distribution patterns, trends, and spatiotemporal variability in the reference evapotranspiration of Pakistan.


## 2. Material and Methods

### 2.1. Area of Study

We consider 41 meteorological locations of Pakistan. These 41 stations are controlled by the Pakistan Meteorological department (PMD) and are located in the Northern and Southern areas of Pakistan. Evapotranspiration estimates are computed using four evapotranspirationET0 models for the monthly data set related to minimum temperature (Tmin), maximum temperature (Tmax), minimum relative humidity (RH min), maximum relative humidity (RH max), average monthly Wind speed (U2), and monthly solar radiation (Rs) of each station for 10 years (from 2006 to 2015).


### 2.2. Models Used for Reference Evapotranspiration

This section reviews the various techniques analyzed in this paper for predicting the referenceET0. It is essential to know the input needs of the methods to determine, which is the most appropriate and optimal model, aligned with the precision.


#### 2.2.1. FAO-56 Penman Monteith (PM)

In this paper, we consider the FAO-56(PM) model of the United Nations Food Agriculture Organization as a standard reference model. The Penman models family is known as the most preciseET0 model in almost every climate [17]. The main drawback of the Penman family is that they often need climatic data inputs, thereby restraining their usefulness in data-restricted sites. Furthermore, there are other significant variations of the FPM model, each with its own set of contradictions regarding clear and vivid climate regions, crops, etc. The FAO-56(PM) model is a substitute adjusted form of actual 1948 Penman model, but FAO-56(PM) model merit for bulk surface resistance (rs) and aerodynamic (ra). With the combination of these two terms, the original Penman method can accurately predict the actual procedures of ET0 from the surface of vegetation. In addition, one can modify the model to fit a specific kind of vegetation, such as a particular crop. The mathematical representation of the FAO−56PM equation is given by(1)λET0PM=δRn−G+ραcpes−ea/raδ+λ1+rs/ra,whereλ  = the latent heat of vaporization Jkg−1,ET0= the evapotranspiration Jm−2,δ= the slope of saturation vapor pressure temperature curve PaK−1,γ= the psychometric constant PaK−1,Rn= the Net Radiation Jm−2s−1,G = soil heat flux Jm−2s−1,ρ= mean air density kgm−3,Cp= specific heat of air Jkg−1k−1,ra and rs are the aerodynamic and surface resistances s−1m.In theFAO56PM manual, we can find formulae for calculating each component, along with illustrations of the units. Allen et al. [18] adjusted the fundamental form of FAO-56(PM) to involve these variables and to form a simplified PM equation FAO-56(PM). Mathematically, we have(2)ET0PM=0.408δRn−G+λ900/T+273U2eS−eaδ+λ1+0.34u2.


#### 2.2.2. Hargreaves-Samani (HGR)

The Hargreaves-Samani model is very convenient, appropriate, and computationally precise to apply to meteorological temperature data in various climate locations. The mathematical representation of the Hargreaves-Samani equation presented inFAO−56PM is as follows:(3)ETHGR=0.00233∗Ra∗Tmax−Tmin0.5Tmax+Tmim2+17.80.


#### 2.2.3. Blaney-Criddle (BCR)

The major drawback of theFAO−56PM method is that it requires the information of the large number of meteorological parameters. These parameters are not always available for all cites. To predict ET0, Blaney-Criddle equations have been developed. For example, the BCR model would be a great asset when only focusing on temperature observations. The BCR model is given by(4)ET0=P∗0.46∗Tmean+8.13,whereTmean  = Average temperature 0C.P = Average daily Percentage of annual daytime hours according to the latitude.


#### 2.2.4. Development of Modified Hargreaves-Samani Model

Significant variables such as solar radiation, relative humidity, wind speed, and temperature are present in the available data. It is also worth knowing that almost all meteorological stations in Pakistan report monthly average relative humidity (RH) and air temperature (T) data. Allen et al. [18] proposed that Hargreaves is an accurate option to use when enough data is not available to predict ET0 by FAO−56PM. Hargreaves’s method requires only minimum and maximum extraterrestrial radiation and air temperature. Subburayan et al. [19] used a similar method that we are also using to calibrate the Hargreaves equation in our study. The new modified Hargreaves equation can be written as follows:(5)Y=CXd,whereY =  FAO−56PM,X=Tmax−Tmin,andC=Tmax+Tmim/2+17.8∗0.0023∗Ra.We have measured the values ofC, X and d (the exponent) by using the data set of ten years of average monthly maximum and minimum relative humidity, maximum and minimum air temperature, wind speed, and solar radiation for each of 41 weather stations of Pakistan. Hence, the value of exponent d=0.382 is found to be appropriate for all stations. The resulting modified Hargreaves equation becomes(6)ETHarg=0.0023∗Ra∗Tmax−Tmin0.382Tmax+Tmim2+17.8.


### 2.3. Performance Evaluation Criteria

The BCR, HGR, and modified Hargreaves conclusions are compared with theFAO−56PM prediction technique used in this analysis as the reference ET0 technique. A variety of statistical parameters have been used, including the Standard Error of Estimate (SEE), the Mean Absolute Error (MAE), the Coefficient of Residual Mass (CRM), and Modeling Efficiency (ME), to evaluate the performance of four models. Moreover, the range of statistical indicators is shown in Table 1.Table 1
Range of statistical indicators.
Statistical indicatorRangeCRM(−∞ to +∞)SEE(0 to∞)ME(−∞ to 1)MAE(0 to+∞)


#### 2.3.1. Model Efficiency

If the model efficiency (ME) is close to 1.0, all predicted (other models) values would tend to match better with the observedFAO−56PM variable. The ME value (near zero) indicates the model’s poor performance, and the ME value (negative) indicates the predicted values are not as good as the simple methods using observed averages. It can be given mathematically as follows:(7)ME=∑i=1ny0,i−y¯2−∑i=1nyp,i−y0,i2∑i=1ny0,i−y¯2,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model FAO−56PM.y¯  = Average of observed variable.


#### 2.3.2. Coefficient of Residual Mass

The coefficient of residual mass (CRM) is widely used to assess the performance of different models. Moreover, it also shows the overestimation or underestimation of the models. A zero value of CRM shows that the estimation is better, and positive and negative values of CRM indicate that the technique is underestimating and overestimating the observed variableFAO−56PM, respectively. It can be shown mathematically as(8)CRM=∑i=1ny0,i−∑i=1nyp,i∑i=1ny0,i,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model. FAO−56PM.


#### 2.3.3. Mean Absolute Error

Mean absolute error is widely used as a performance measure for comparing two or more models. A small value of MAE shows good performance of the model. It can be shown mathematically as(9)MAE=∑i=1NETPM−ETOMN,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


#### 2.3.4. Standard Error of Estimate

The standard error of estimate (SEE) can also be used to compare the performance of different models. It is usually used to relate the good performance of fit between the estimated and the observed variables. Thus, SEE can be adopted in the performance computation of models to know about the precision and accuracy of the involvedET0 technique by contrasting it with the standard reference procedure. A general value of SEE preferably near to zero shows good performance of the model. A large value deviating from zero indicates a poor fit. It can be shown mathematically as(10)SEE=∑i=1NETPM−ETOM2N−10.5,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


### 2.4. Seasonal Data

In this study, 41 meteorological stations of Pakistan are chosen to evaluate the season-based performance of two models based on radiation, and one model based on temperature for predicting for three seasons, i.e., rainy, dry, and annual, on annual time stage. The major objective of the current study is to examine the efficiency of these models. The seasonal data of Bahawalnagar is presented in Table2.Table 2
Summary of Mean Weather Parameters for 10 years of data for Bahawalnagar.
MonthTminTmaxTmeanHumidityWind speedRaJan3.8518.4911.1766.93.20521.1Feb7.942214.9763.54.637525.8Mar13.6628.1420.953.75.372531.4Apr19.2134.8327.0239.555.89536.8May24.439.6932.04533.66.5740June26.6939.9533.3243.16.86541.2July26.5836.0831.3364.455.94540.6Aug25.9734.9130.4469.055.2338Sep23.7234.0128.86567.34.26533.4Oct18.7432.4825.6158.253.1627.6Nov11.4327.119.26562.151.97522.2Dec5.520.713.167.752.3319.8


### 2.5. Trend Analysis

The mean decadeET0 values for the time period 2006–2015 were evaluated for 41 meteorological stations of Pakistan study regions using three variable models. A statistical pair test is used to contrast the monthly and decade ET0 for the three models. Initially, a comparison is made between the two temperature-based and radiation-based methods and then between each method based on the temperature and the FAO−56PM standard model. The standard error of estimate (SEE) and the coefficients of determination R2 were computed to know about the best technique. The standard error of estimate (SEE) is the prediction of the average deviation of the regression from FAO−56PM predicted values. The equations for standard error of estimate (SEE) and coefficients of determination R2 are given as follows:(11)SEE=∑i=1Nxi−x′N−22,R2=∑i=1Nxi−x¯xi′−x¯′∑i=1Nxi−x¯2∑i=1Nxi′−x¯i′22.Here,xi represents the FAO−56PM observed ET0, xi′ is the regression estimated ET0 for the ith value, xi and xi′ represent the mean values of the corresponding variable, and N represents the number of observations taken into account. A simple linear regression x=0+1y with zero intercepts and unit slope and x and y as the endogenous and exogenous variable, respectively.


### 2.6. Ordinary Kriging for Functional Data (OKFD)

Increasing amounts of space-time data are being gathered and processed due to technological development and the requirement of society to study variables that change in both spatial and temporal domains, for example, climate, air quality variations, and crop production. Analysis of space and time correlations is vital to study the causes and character of variability but is also required for estimating the observations at points from the nearing values.The spatiotemporal interpolation is easy and effective to consider because it can give more appropriate estimations than the space interpolation ignoring the time points. Delicado et al. [20] suggested a classical usual kriging predictor but considered the curves into account instead of unidimensional data; that is, every curve is measured by a scalar technique. This method is called “Ordinary Kriging for function-valued spatial data.” Delicado et al. [20] also sorted out the difficulty of spatial forecasting of functioned observations by measuring each examined curve through a functional technique.Ferraty and Vieu [21] elaborated a functional variate as a usual variable X attaining the observations in an unlimited dimensional space. A functional set of data x1,……,xn is the observation of n function variables Y1;….;Yn distributed as Y. Let T=c;d⊆R. We implement the functional data constituted as(12)L2T=Y:T⇒R,∣∫TYt2dt<∞.Note thatL2T with the inner product x,y=∫Tytxtxdt explains a Euclidean space. Following Giraldo et al. [13], we elaborate a usual functional process as Yst:s∈D⇒Rd;t∈T⇒R, where d is usually taken as 2, and Yst is a function variable for any s∈D. Let s1…..sn be the hypothetical points in D and predict that we can examine a recognition of the functional usual process Yst at these “n” sites, ys1t;…..;ysnt. The Ordinary Kriging for the functional data is a geostatistical method for forecasting Ys0t, the functional usual process at s0, where s0 is an unsampled region.It is commonly predicted that the functional usual process is second-order isotropic and stationary; that is, the average and the variance function are constant, and the covariance relies only upon the extension between modeling points (however, the technique could also be carried out without predicting these parameters). Formally, we may suppose that(1)
EYst=mt and VYst=σ2t for all s∈D and all t∈T.(2)
COVYsit,Ysjt=Csi−sjt=Cijh,t;si,sj∈D,t∈T,whereh=si−sj.(3)
1/2V=Ysit−Ysjt=γsi−sjt=ci,jh,t;si,sj∈D,t∈T,      h=si−sj.These assumptions imply thatVYsit−Ysjt  =  EYsit−Ysjt2 and γsi−sjt=σ2t−Csi−sjt. According to Giraldo et al. [15], the OKFD predictor is defined by(13)Y^sot=∑i=1nλiYsit,λ1,…..,λn∈R.The forecaster in equation (12) is identical to the simple usual kriging forecaster of Cressie and Wikle [9] but takes curves into account rather than variables. The curve is a linear aggregation of examined curves. We assume that each calculated curve is an entire datum, so we account for the whole curve as the only quantity. The weights λ in equation (12) give the effect of the curves encircling the unsampled region where we like to carry out our estimation. The curves from the regions nearer to the estimation value will automatically have more effect than the others, which are away. These weights are predicted in a way that the forecaster in equation (12) becomes excellent Best Linear Unbiased Predictor (BLUP). We predict that every examined function can be such in terms of K basis function, B1t,….,BKt, as(14)Ysit=∑l=1KailBlt,i=1,2,..,n,where ai=ai1,…..,aiK, Bt=B1t,…..,BKt.These expressions are short versions of the Fourier series (for periodic functions) or B-splines expansions.In order to find BLUP, we first contemplate the unbiasedness. From the average constant situation above, we need the condition∑i=1nλi=1. In a typical geostatistical setting, we suppose that the examinations are recognition of a usual field Ys:s∈D,D∈Rd. The kriging predictor is defined as Y^so=∑i=1nλiYsi and the BLUP is achieved by reducing(15)σso2=VY^si−Yso,subject to ∑i=1nλi=1. Furthermore, in multivariate geo-statistics, the data consists of Ys1,…,Ysn; that is, we have observations of a spatial vector-valued process Ys:s∈D, where Ys=Ys1,….,Ysm and D∈Rd. In this context, VY^si−Yso is a matrix, and the BLUP of m variables at an unsampled region s0 can be achieved by reducing(16)σso2=∑i=1mVY^si−Yso.Bounded by restraints that guarantee unbiasedness, that is, reducing the proof of the mean-squared estimated error matrix is exposed to some constraints. Extending the technique to the functional framework by substituting the summation with an integral, the “n” parameters in equation (13) are achieved by calculating the given constrained optimization difficulty [15].(17)minλi,….,λn∫TVY^si−Ysodt|∑i=1nλi=1,which after some simplification can be given as(18)∑i=1n∑j=1nλiλj∫TCijh,tdt+∫Tσ2tdt−2∑i=1n∫TCi0h,tdt+2μ∑i=1nλi−1,where μ is the Lagrange multiplier used to consider the unbiasedness limitation. Reducing equation (18) in context with λ1,….,λn and μ, we carry out the given linear system of equations, which allows us to predict the model, with these systems of equations being shown in matrix form:(19)λ1∫Tγs1−s1tdt+…..+λn∫Tγs1−sntdt−μ=∫Tγso−s1tdt,λ1∫Tγsn−s1tdt+…..+λn∫Tγsn−sntdt−μ=∫Tγso−sntdt,λ1+λ2+……+λn=1,The function(20)γh=∫Tγsi−sjtdt,is called the trace-variogram. To deal with the system in equation (20), a predictor of the trace-variogram is required. Given that, we are supposing that Yst has a constant average function mt over D.(21)VYsit−Ysjt=EYsit−Ysjt2.By using Fubini’s algorithm,(22)γh=12E∫TYsit−Ysjt2dt,where si,sjϵD with h=si−sj. Then, adherence of the typical method-of-moments for this predicted value shows the following predictor:(23)γ^h=12Nh∑i,j∈Nh∫TYsit−Ysjt2dt,where Nh=si,sj:si−sj=h, and Nh is the number of distinct elements in N(h). For random spaced data, there are not normally adequate observations parted accurately by a distance h. Then, N(h) is changed to si,sj:si−sj∈h−ε;h+ε, with ε>0 being a small value.Once we have predicted the trace-variogram for K valueshk,k=1,2,……,K, a parametric technique h;θ like “Exponential,” “Gaussian,” “Spherical,” or “Mathern” must be equipped [22]. The estimation of trace-variance of the functional data, the usual kriging depending on the trace-variogram, is given by(24)σ^so2=∫TVY^sot−Ysotdt=∑i=1nλi∫Tγsi−sotdt−μ.This technique should be regarded as a global perplexity measure because it assimilates the typical point-wise estimation variance of the usual kriging. Spatial functional data and spatiotemporal parameters have a causal relation in the context that we have developed a spatial process by time or by any additional feature. The space and the time models contemplate the development of a spatial process by time and model the spatiotemporal dependence. Thus, we have one variableYs,t, and we need to forecast a variable at an unsampled region. In the case of a space function, Yst is a function, so we are estimating a function.


## 2.1. Area of Study

We consider 41 meteorological locations of Pakistan. These 41 stations are controlled by the Pakistan Meteorological department (PMD) and are located in the Northern and Southern areas of Pakistan. Evapotranspiration estimates are computed using four evapotranspirationET0 models for the monthly data set related to minimum temperature (Tmin), maximum temperature (Tmax), minimum relative humidity (RH min), maximum relative humidity (RH max), average monthly Wind speed (U2), and monthly solar radiation (Rs) of each station for 10 years (from 2006 to 2015).


## 2.2. Models Used for Reference Evapotranspiration

This section reviews the various techniques analyzed in this paper for predicting the referenceET0. It is essential to know the input needs of the methods to determine, which is the most appropriate and optimal model, aligned with the precision.


### 2.2.1. FAO-56 Penman Monteith (PM)

In this paper, we consider the FAO-56(PM) model of the United Nations Food Agriculture Organization as a standard reference model. The Penman models family is known as the most preciseET0 model in almost every climate [17]. The main drawback of the Penman family is that they often need climatic data inputs, thereby restraining their usefulness in data-restricted sites. Furthermore, there are other significant variations of the FPM model, each with its own set of contradictions regarding clear and vivid climate regions, crops, etc. The FAO-56(PM) model is a substitute adjusted form of actual 1948 Penman model, but FAO-56(PM) model merit for bulk surface resistance (rs) and aerodynamic (ra). With the combination of these two terms, the original Penman method can accurately predict the actual procedures of ET0 from the surface of vegetation. In addition, one can modify the model to fit a specific kind of vegetation, such as a particular crop. The mathematical representation of the FAO−56PM equation is given by(1)λET0PM=δRn−G+ραcpes−ea/raδ+λ1+rs/ra,whereλ  = the latent heat of vaporization Jkg−1,ET0= the evapotranspiration Jm−2,δ= the slope of saturation vapor pressure temperature curve PaK−1,γ= the psychometric constant PaK−1,Rn= the Net Radiation Jm−2s−1,G = soil heat flux Jm−2s−1,ρ= mean air density kgm−3,Cp= specific heat of air Jkg−1k−1,ra and rs are the aerodynamic and surface resistances s−1m.In theFAO56PM manual, we can find formulae for calculating each component, along with illustrations of the units. Allen et al. [18] adjusted the fundamental form of FAO-56(PM) to involve these variables and to form a simplified PM equation FAO-56(PM). Mathematically, we have(2)ET0PM=0.408δRn−G+λ900/T+273U2eS−eaδ+λ1+0.34u2.


### 2.2.2. Hargreaves-Samani (HGR)

The Hargreaves-Samani model is very convenient, appropriate, and computationally precise to apply to meteorological temperature data in various climate locations. The mathematical representation of the Hargreaves-Samani equation presented inFAO−56PM is as follows:(3)ETHGR=0.00233∗Ra∗Tmax−Tmin0.5Tmax+Tmim2+17.80.


### 2.2.3. Blaney-Criddle (BCR)

The major drawback of theFAO−56PM method is that it requires the information of the large number of meteorological parameters. These parameters are not always available for all cites. To predict ET0, Blaney-Criddle equations have been developed. For example, the BCR model would be a great asset when only focusing on temperature observations. The BCR model is given by(4)ET0=P∗0.46∗Tmean+8.13,whereTmean  = Average temperature 0C.P = Average daily Percentage of annual daytime hours according to the latitude.


### 2.2.4. Development of Modified Hargreaves-Samani Model

Significant variables such as solar radiation, relative humidity, wind speed, and temperature are present in the available data. It is also worth knowing that almost all meteorological stations in Pakistan report monthly average relative humidity (RH) and air temperature (T) data. Allen et al. [18] proposed that Hargreaves is an accurate option to use when enough data is not available to predict ET0 by FAO−56PM. Hargreaves’s method requires only minimum and maximum extraterrestrial radiation and air temperature. Subburayan et al. [19] used a similar method that we are also using to calibrate the Hargreaves equation in our study. The new modified Hargreaves equation can be written as follows:(5)Y=CXd,whereY =  FAO−56PM,X=Tmax−Tmin,andC=Tmax+Tmim/2+17.8∗0.0023∗Ra.We have measured the values ofC, X and d (the exponent) by using the data set of ten years of average monthly maximum and minimum relative humidity, maximum and minimum air temperature, wind speed, and solar radiation for each of 41 weather stations of Pakistan. Hence, the value of exponent d=0.382 is found to be appropriate for all stations. The resulting modified Hargreaves equation becomes(6)ETHarg=0.0023∗Ra∗Tmax−Tmin0.382Tmax+Tmim2+17.8.


## 2.2.1. FAO-56 Penman Monteith (PM)

In this paper, we consider the FAO-56(PM) model of the United Nations Food Agriculture Organization as a standard reference model. The Penman models family is known as the most preciseET0 model in almost every climate [17]. The main drawback of the Penman family is that they often need climatic data inputs, thereby restraining their usefulness in data-restricted sites. Furthermore, there are other significant variations of the FPM model, each with its own set of contradictions regarding clear and vivid climate regions, crops, etc. The FAO-56(PM) model is a substitute adjusted form of actual 1948 Penman model, but FAO-56(PM) model merit for bulk surface resistance (rs) and aerodynamic (ra). With the combination of these two terms, the original Penman method can accurately predict the actual procedures of ET0 from the surface of vegetation. In addition, one can modify the model to fit a specific kind of vegetation, such as a particular crop. The mathematical representation of the FAO−56PM equation is given by(1)λET0PM=δRn−G+ραcpes−ea/raδ+λ1+rs/ra,whereλ  = the latent heat of vaporization Jkg−1,ET0= the evapotranspiration Jm−2,δ= the slope of saturation vapor pressure temperature curve PaK−1,γ= the psychometric constant PaK−1,Rn= the Net Radiation Jm−2s−1,G = soil heat flux Jm−2s−1,ρ= mean air density kgm−3,Cp= specific heat of air Jkg−1k−1,ra and rs are the aerodynamic and surface resistances s−1m.In theFAO56PM manual, we can find formulae for calculating each component, along with illustrations of the units. Allen et al. [18] adjusted the fundamental form of FAO-56(PM) to involve these variables and to form a simplified PM equation FAO-56(PM). Mathematically, we have(2)ET0PM=0.408δRn−G+λ900/T+273U2eS−eaδ+λ1+0.34u2.


## 2.2.2. Hargreaves-Samani (HGR)

The Hargreaves-Samani model is very convenient, appropriate, and computationally precise to apply to meteorological temperature data in various climate locations. The mathematical representation of the Hargreaves-Samani equation presented inFAO−56PM is as follows:(3)ETHGR=0.00233∗Ra∗Tmax−Tmin0.5Tmax+Tmim2+17.80.


## 2.2.3. Blaney-Criddle (BCR)

The major drawback of theFAO−56PM method is that it requires the information of the large number of meteorological parameters. These parameters are not always available for all cites. To predict ET0, Blaney-Criddle equations have been developed. For example, the BCR model would be a great asset when only focusing on temperature observations. The BCR model is given by(4)ET0=P∗0.46∗Tmean+8.13,whereTmean  = Average temperature 0C.P = Average daily Percentage of annual daytime hours according to the latitude.


## 2.2.4. Development of Modified Hargreaves-Samani Model

Significant variables such as solar radiation, relative humidity, wind speed, and temperature are present in the available data. It is also worth knowing that almost all meteorological stations in Pakistan report monthly average relative humidity (RH) and air temperature (T) data. Allen et al. [18] proposed that Hargreaves is an accurate option to use when enough data is not available to predict ET0 by FAO−56PM. Hargreaves’s method requires only minimum and maximum extraterrestrial radiation and air temperature. Subburayan et al. [19] used a similar method that we are also using to calibrate the Hargreaves equation in our study. The new modified Hargreaves equation can be written as follows:(5)Y=CXd,whereY =  FAO−56PM,X=Tmax−Tmin,andC=Tmax+Tmim/2+17.8∗0.0023∗Ra.We have measured the values ofC, X and d (the exponent) by using the data set of ten years of average monthly maximum and minimum relative humidity, maximum and minimum air temperature, wind speed, and solar radiation for each of 41 weather stations of Pakistan. Hence, the value of exponent d=0.382 is found to be appropriate for all stations. The resulting modified Hargreaves equation becomes(6)ETHarg=0.0023∗Ra∗Tmax−Tmin0.382Tmax+Tmim2+17.8.


## 2.3. Performance Evaluation Criteria

The BCR, HGR, and modified Hargreaves conclusions are compared with theFAO−56PM prediction technique used in this analysis as the reference ET0 technique. A variety of statistical parameters have been used, including the Standard Error of Estimate (SEE), the Mean Absolute Error (MAE), the Coefficient of Residual Mass (CRM), and Modeling Efficiency (ME), to evaluate the performance of four models. Moreover, the range of statistical indicators is shown in Table 1.Table 1
Range of statistical indicators.
Statistical indicatorRangeCRM(−∞ to +∞)SEE(0 to∞)ME(−∞ to 1)MAE(0 to+∞)


### 2.3.1. Model Efficiency

If the model efficiency (ME) is close to 1.0, all predicted (other models) values would tend to match better with the observedFAO−56PM variable. The ME value (near zero) indicates the model’s poor performance, and the ME value (negative) indicates the predicted values are not as good as the simple methods using observed averages. It can be given mathematically as follows:(7)ME=∑i=1ny0,i−y¯2−∑i=1nyp,i−y0,i2∑i=1ny0,i−y¯2,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model FAO−56PM.y¯  = Average of observed variable.


### 2.3.2. Coefficient of Residual Mass

The coefficient of residual mass (CRM) is widely used to assess the performance of different models. Moreover, it also shows the overestimation or underestimation of the models. A zero value of CRM shows that the estimation is better, and positive and negative values of CRM indicate that the technique is underestimating and overestimating the observed variableFAO−56PM, respectively. It can be shown mathematically as(8)CRM=∑i=1ny0,i−∑i=1nyp,i∑i=1ny0,i,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model. FAO−56PM.


### 2.3.3. Mean Absolute Error

Mean absolute error is widely used as a performance measure for comparing two or more models. A small value of MAE shows good performance of the model. It can be shown mathematically as(9)MAE=∑i=1NETPM−ETOMN,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


### 2.3.4. Standard Error of Estimate

The standard error of estimate (SEE) can also be used to compare the performance of different models. It is usually used to relate the good performance of fit between the estimated and the observed variables. Thus, SEE can be adopted in the performance computation of models to know about the precision and accuracy of the involvedET0 technique by contrasting it with the standard reference procedure. A general value of SEE preferably near to zero shows good performance of the model. A large value deviating from zero indicates a poor fit. It can be shown mathematically as(10)SEE=∑i=1NETPM−ETOM2N−10.5,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


## 2.3.1. Model Efficiency

If the model efficiency (ME) is close to 1.0, all predicted (other models) values would tend to match better with the observedFAO−56PM variable. The ME value (near zero) indicates the model’s poor performance, and the ME value (negative) indicates the predicted values are not as good as the simple methods using observed averages. It can be given mathematically as follows:(7)ME=∑i=1ny0,i−y¯2−∑i=1nyp,i−y0,i2∑i=1ny0,i−y¯2,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model FAO−56PM.y¯  = Average of observed variable.


## 2.3.2. Coefficient of Residual Mass

The coefficient of residual mass (CRM) is widely used to assess the performance of different models. Moreover, it also shows the overestimation or underestimation of the models. A zero value of CRM shows that the estimation is better, and positive and negative values of CRM indicate that the technique is underestimating and overestimating the observed variableFAO−56PM, respectively. It can be shown mathematically as(8)CRM=∑i=1ny0,i−∑i=1nyp,i∑i=1ny0,i,wheren = number of observations.yp,i  = is the estimated models [HGR, BCR and M-H].y0,i  = observed variable or model. FAO−56PM.


## 2.3.3. Mean Absolute Error

Mean absolute error is widely used as a performance measure for comparing two or more models. A small value of MAE shows good performance of the model. It can be shown mathematically as(9)MAE=∑i=1NETPM−ETOMN,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


## 2.3.4. Standard Error of Estimate

The standard error of estimate (SEE) can also be used to compare the performance of different models. It is usually used to relate the good performance of fit between the estimated and the observed variables. Thus, SEE can be adopted in the performance computation of models to know about the precision and accuracy of the involvedET0 technique by contrasting it with the standard reference procedure. A general value of SEE preferably near to zero shows good performance of the model. A large value deviating from zero indicates a poor fit. It can be shown mathematically as(10)SEE=∑i=1NETPM−ETOM2N−10.5,whereN = total number of observations.ETPM  =  ET0 calculated with FAO-56(PM).ETOM  =  ET0 calculated with [HGR, BCR and M-H].


## 2.4. Seasonal Data

In this study, 41 meteorological stations of Pakistan are chosen to evaluate the season-based performance of two models based on radiation, and one model based on temperature for predicting for three seasons, i.e., rainy, dry, and annual, on annual time stage. The major objective of the current study is to examine the efficiency of these models. The seasonal data of Bahawalnagar is presented in Table2.Table 2
Summary of Mean Weather Parameters for 10 years of data for Bahawalnagar.
MonthTminTmaxTmeanHumidityWind speedRaJan3.8518.4911.1766.93.20521.1Feb7.942214.9763.54.637525.8Mar13.6628.1420.953.75.372531.4Apr19.2134.8327.0239.555.89536.8May24.439.6932.04533.66.5740June26.6939.9533.3243.16.86541.2July26.5836.0831.3364.455.94540.6Aug25.9734.9130.4469.055.2338Sep23.7234.0128.86567.34.26533.4Oct18.7432.4825.6158.253.1627.6Nov11.4327.119.26562.151.97522.2Dec5.520.713.167.752.3319.8


## 2.5. Trend Analysis

The mean decadeET0 values for the time period 2006–2015 were evaluated for 41 meteorological stations of Pakistan study regions using three variable models. A statistical pair test is used to contrast the monthly and decade ET0 for the three models. Initially, a comparison is made between the two temperature-based and radiation-based methods and then between each method based on the temperature and the FAO−56PM standard model. The standard error of estimate (SEE) and the coefficients of determination R2 were computed to know about the best technique. The standard error of estimate (SEE) is the prediction of the average deviation of the regression from FAO−56PM predicted values. The equations for standard error of estimate (SEE) and coefficients of determination R2 are given as follows:(11)SEE=∑i=1Nxi−x′N−22,R2=∑i=1Nxi−x¯xi′−x¯′∑i=1Nxi−x¯2∑i=1Nxi′−x¯i′22.Here,xi represents the FAO−56PM observed ET0, xi′ is the regression estimated ET0 for the ith value, xi and xi′ represent the mean values of the corresponding variable, and N represents the number of observations taken into account. A simple linear regression x=0+1y with zero intercepts and unit slope and x and y as the endogenous and exogenous variable, respectively.


## 2.6. Ordinary Kriging for Functional Data (OKFD)

Increasing amounts of space-time data are being gathered and processed due to technological development and the requirement of society to study variables that change in both spatial and temporal domains, for example, climate, air quality variations, and crop production. Analysis of space and time correlations is vital to study the causes and character of variability but is also required for estimating the observations at points from the nearing values.The spatiotemporal interpolation is easy and effective to consider because it can give more appropriate estimations than the space interpolation ignoring the time points. Delicado et al. [20] suggested a classical usual kriging predictor but considered the curves into account instead of unidimensional data; that is, every curve is measured by a scalar technique. This method is called “Ordinary Kriging for function-valued spatial data.” Delicado et al. [20] also sorted out the difficulty of spatial forecasting of functioned observations by measuring each examined curve through a functional technique.Ferraty and Vieu [21] elaborated a functional variate as a usual variable X attaining the observations in an unlimited dimensional space. A functional set of data x1,……,xn is the observation of n function variables Y1;….;Yn distributed as Y. Let T=c;d⊆R. We implement the functional data constituted as(12)L2T=Y:T⇒R,∣∫TYt2dt<∞.Note thatL2T with the inner product x,y=∫Tytxtxdt explains a Euclidean space. Following Giraldo et al. [13], we elaborate a usual functional process as Yst:s∈D⇒Rd;t∈T⇒R, where d is usually taken as 2, and Yst is a function variable for any s∈D. Let s1…..sn be the hypothetical points in D and predict that we can examine a recognition of the functional usual process Yst at these “n” sites, ys1t;…..;ysnt. The Ordinary Kriging for the functional data is a geostatistical method for forecasting Ys0t, the functional usual process at s0, where s0 is an unsampled region.It is commonly predicted that the functional usual process is second-order isotropic and stationary; that is, the average and the variance function are constant, and the covariance relies only upon the extension between modeling points (however, the technique could also be carried out without predicting these parameters). Formally, we may suppose that(1)
EYst=mt and VYst=σ2t for all s∈D and all t∈T.(2)
COVYsit,Ysjt=Csi−sjt=Cijh,t;si,sj∈D,t∈T,whereh=si−sj.(3)
1/2V=Ysit−Ysjt=γsi−sjt=ci,jh,t;si,sj∈D,t∈T,      h=si−sj.These assumptions imply thatVYsit−Ysjt  =  EYsit−Ysjt2 and γsi−sjt=σ2t−Csi−sjt. According to Giraldo et al. [15], the OKFD predictor is defined by(13)Y^sot=∑i=1nλiYsit,λ1,…..,λn∈R.The forecaster in equation (12) is identical to the simple usual kriging forecaster of Cressie and Wikle [9] but takes curves into account rather than variables. The curve is a linear aggregation of examined curves. We assume that each calculated curve is an entire datum, so we account for the whole curve as the only quantity. The weights λ in equation (12) give the effect of the curves encircling the unsampled region where we like to carry out our estimation. The curves from the regions nearer to the estimation value will automatically have more effect than the others, which are away. These weights are predicted in a way that the forecaster in equation (12) becomes excellent Best Linear Unbiased Predictor (BLUP). We predict that every examined function can be such in terms of K basis function, B1t,….,BKt, as(14)Ysit=∑l=1KailBlt,i=1,2,..,n,where ai=ai1,…..,aiK, Bt=B1t,…..,BKt.These expressions are short versions of the Fourier series (for periodic functions) or B-splines expansions.In order to find BLUP, we first contemplate the unbiasedness. From the average constant situation above, we need the condition∑i=1nλi=1. In a typical geostatistical setting, we suppose that the examinations are recognition of a usual field Ys:s∈D,D∈Rd. The kriging predictor is defined as Y^so=∑i=1nλiYsi and the BLUP is achieved by reducing(15)σso2=VY^si−Yso,subject to ∑i=1nλi=1. Furthermore, in multivariate geo-statistics, the data consists of Ys1,…,Ysn; that is, we have observations of a spatial vector-valued process Ys:s∈D, where Ys=Ys1,….,Ysm and D∈Rd. In this context, VY^si−Yso is a matrix, and the BLUP of m variables at an unsampled region s0 can be achieved by reducing(16)σso2=∑i=1mVY^si−Yso.Bounded by restraints that guarantee unbiasedness, that is, reducing the proof of the mean-squared estimated error matrix is exposed to some constraints. Extending the technique to the functional framework by substituting the summation with an integral, the “n” parameters in equation (13) are achieved by calculating the given constrained optimization difficulty [15].(17)minλi,….,λn∫TVY^si−Ysodt|∑i=1nλi=1,which after some simplification can be given as(18)∑i=1n∑j=1nλiλj∫TCijh,tdt+∫Tσ2tdt−2∑i=1n∫TCi0h,tdt+2μ∑i=1nλi−1,where μ is the Lagrange multiplier used to consider the unbiasedness limitation. Reducing equation (18) in context with λ1,….,λn and μ, we carry out the given linear system of equations, which allows us to predict the model, with these systems of equations being shown in matrix form:(19)λ1∫Tγs1−s1tdt+…..+λn∫Tγs1−sntdt−μ=∫Tγso−s1tdt,λ1∫Tγsn−s1tdt+…..+λn∫Tγsn−sntdt−μ=∫Tγso−sntdt,λ1+λ2+……+λn=1,The function(20)γh=∫Tγsi−sjtdt,is called the trace-variogram. To deal with the system in equation (20), a predictor of the trace-variogram is required. Given that, we are supposing that Yst has a constant average function mt over D.(21)VYsit−Ysjt=EYsit−Ysjt2.By using Fubini’s algorithm,(22)γh=12E∫TYsit−Ysjt2dt,where si,sjϵD with h=si−sj. Then, adherence of the typical method-of-moments for this predicted value shows the following predictor:(23)γ^h=12Nh∑i,j∈Nh∫TYsit−Ysjt2dt,where Nh=si,sj:si−sj=h, and Nh is the number of distinct elements in N(h). For random spaced data, there are not normally adequate observations parted accurately by a distance h. Then, N(h) is changed to si,sj:si−sj∈h−ε;h+ε, with ε>0 being a small value.Once we have predicted the trace-variogram for K valueshk,k=1,2,……,K, a parametric technique h;θ like “Exponential,” “Gaussian,” “Spherical,” or “Mathern” must be equipped [22]. The estimation of trace-variance of the functional data, the usual kriging depending on the trace-variogram, is given by(24)σ^so2=∫TVY^sot−Ysotdt=∑i=1nλi∫Tγsi−sotdt−μ.This technique should be regarded as a global perplexity measure because it assimilates the typical point-wise estimation variance of the usual kriging. Spatial functional data and spatiotemporal parameters have a causal relation in the context that we have developed a spatial process by time or by any additional feature. The space and the time models contemplate the development of a spatial process by time and model the spatiotemporal dependence. Thus, we have one variableYs,t, and we need to forecast a variable at an unsampled region. In the case of a space function, Yst is a function, so we are estimating a function.


## 3. Results and Discussions

### 3.1. Estimation of Evapotranspiration Models

The summary of averagedET0 estimation in (mm/day) by using FAO−56PM, BCR, and HGR for Bahawalnagar station is given in Table 3. These results are given in Tables 4 and 5.Table 3
Summary of annual evapotranspiration in (mm/day) by usingFAO−56PM, BCR, and HGR for Bahawalnagar.
JanFebMarAprMayJuneJulyAugSepOctNovDecFAO−56PM2.873.724.856.026.596.65.735.194.754.163.252.69BCR3.183.754.795.967.097.516.996.645.995.184.083.26HGR5.387.2910.6414.9917.9317.6414.1412.6111.510.217.495.49Table 4
Selected Meteorological stations of Pakistan with their coordinates and elevation.
Sr#Name of stationsLatitudeLongitudeElevation1Balakot34.5572.35995.42CHerat33.8271.5513723D.i.khan31.9270.86172.34Dir35.271.8513755Kalam35.8472.982103.16Chilas35.4274.112507Bunji35.6674.613728Gilgit32.9274.3424609Astore35.3374.9216810Chhor29.8869.72511Hyderabad25.3868.422812Jacobebad28.368.475513Naawabshah26.2568.373714Baddin24.6368.9915Karachi24.966.932216Shhukr27.768.8668.5217Rohri27.6668.96618Mirpur khas25.52691519Jiwani25.0661.815620Pajngur26.9764.196821Pasni25.2663.48922Gawadar25.1362.3329. 8623Lesbella26.2366.178724Ormara25.264.66225Turbat25.9863.0615526Bahawalnagar29.3373.85161.0527Bahawalpur29.3371.7811028D g khan30.0570.63148.129Faislabad31.4370.63148.130Jhang31.2672.3119031Jhelum32.9873.75287.1932Khanpore28.6570.6888.4133Lahore apt31.5874.4216.1534Lahore pbo31.5574.3321435Mangla33.0673.63283.336Mbdin32.9673.8252.9737Multan30.271.43121.9539Okara30.873.4318040Ry khan28.4370.3282.9341Sargodha32.0572.6618742Sialkot32.5174.53255.1Table 5
Extraterrestrial radiation for selected meteorological stations of Pakistan.
StationsJanFebMarAprMayJuneJulyAugSepOctNovDecBalakot18.723.729.936.14041.640.837.632.125.619.917.3CHerat18.723.729.936.14041.640.837.632.125.619.917.3Di khan19.924.830.736.54041.440.737.932.826.621.118.5Dir17.522.62935.74041.740.837.431.524.618.716.1Kalam17.522.62935.74041.740.837.431.524.618.716.1Chilas17.522.62935.74041.740.837.431.524.618.716.1Bunji17.522.62935.74041.740.837.431.524.618.716.1Gilgit19.924.830.736.54041.440.737.932.826.621.118.5Astore17.522.62935.74041.740.837.431.524.618.716.1Chhor21.125.831.436.84041.240.63833.427.622.219.8Hyderabad23.427.832.835.439.940.640.238.334.529.324.522.2Jacobebad22.325.832.237.14040.940.438.233.928.523.321Nawabshah23.427.832.835.439.940.640.238.334.529.324.522.2Baddin24.628.833.535.639.740.339.938.334.930.225.523.3Karachi24.628.833.535.639.740.339.938.334.930.225.523.3Shhukr22.325.832.237.14040.940.438.233.928.523.321Rohri22.325.832.237.14040.940.438.233.928.523.321MP khas23.427.832.835.439.940.640.238.334.529.324.522.2Jiwani23.427.832.835.439.940.640.238.334.529.324.522.2Panjgur23.427.832.835.439.940.640.238.334.529.324.522.2Pasni23.427.832.835.439.940.640.238.334.529.324.522.2Gawadr23.427.832.835.439.940.640.238.334.529.324.522.2Lasbella23.427.832.835.439.940.640.238.334.529.324.522.2Ormara23.427.832.835.439.940.640.238.334.529.324.522.2Turbat23.427.832.835.439.940.640.238.334.529.324.522.2Bahawalngr21.125.831.436.84041.240.63833.427.622.219.8Bahawalpur21.125.831.436.84041.240.63833.427.622.219.8D g khan21.125.831.436.84041.240.63833.427.622.219.8Faislabad19.924.830.736.54041.440.737.932.826.621.118.5Jhang19.924.830.736.54041.440.737.932.826.621.118.5Jhelum19.924.830.736.54041.440.737.932.826.621.118.5Khanpore22.325.832.237.14040.940.438.233.928.523.321Lahore apt19.924.830.736.54041.440.737.932.826.621.118.5Lahore Pbo19.924.830.736.54041.440.737.932.826.621.118.5Mangla18.723.729.936.14041.640.837.632.125.619.917.3Mbdin19.924.830.736.54041.440.737.932.826.621.118.5Multan21.125.831.436.84041.240.63833.427.622.219.8Okara21.125.831.436.84041.240.63833.427.622.219.8R Y khan22.325.832.237.14040.940.438.233.928.523.321Sargodha19.924.830.736.54041.440.737.932.826.621.118.5Sialkot19.924.830.736.54041.440.737.932.826.621.118.5


### 3.2. Performance Evaluation of Evapotranspiration Models

The results of Model Efficiency (ME) and Coefficient of Residual Mass (CRM) for different models are given in Table6(a). The results of other performance measures, i.e., Mean Absolute Error (MAE) and Standard Error of Estimate (SEE), are given in Table 6(b). From the results of these performance indicators, we can see that the HGR model has better performance as compared to the BCR and FPM models in nonseasonal data.Table 6
Performance Evaluation of Selected Models based on ME and CRM (Left panel), Performance Evaluation of Selected Models based on MAE and SEE (Right panel).
(a) Performance based on ME and CRMMECRMStationsBCRM-HO-HBCRM-HO-HBahawalnagar2.11.974.343.3221.7147.77Bahawalpur1.621.216.04213.3466.54DG Khan1.421.225.491.9813.4660.49Faislabad1.971.545.162.6816.8356Jhang1.981.646.162.6915.8458.1Khanpore1.390.996.591.810.9772.59Lahore APT2.162.034.173.3322.3545.9Lahore PBO1.621.035.271.8711.457.98Mangla1.841.285.82.4114.1463.8M B din1.510.926.511.6310.1971.64Multan2.252.032.02322.2951.04Okara1.340.767.111.6912.6883.3R Y Khan1.611.216.582.213.3672.43Sargodha1.50.896.751.639.8874.31Sialkot1.610.986.451.7610.7770.99Jehalum1.450.916.491.6510.0371.43Badin4.745.011.785.7855.219.59Chhor5.956.21.997.2568.2321.9Gawadar3.253.62.593.5639.728.49Hyderabad8.588.832.479.6897.1727.18Jecobabad2.62.524.593.6227.7750.56Jiwini1.962.582.352.3428.4625.92Karachi4.264.761.134.6652.3812.44Lasbella5.455.442.86.3759.9230.89Mirpur khas2.162.045.323.0722.5458.59Nawabshah6.166.461.997.4871.0621.98Ormara1.852.183.722.1324.0440.92Panjgur6.026.041.776.8466.4519.52Pasni2.913.422.253.2937.7224.78Rohri2.422.634.253.0428.9546.79Shukkur3.433.424.054.6637.744.59Turbot9.479.452.1510.44103.923.71Astore0.810.823.230.999.1235.53Balakot1.060.615.891.186.7264.88Bunji0.860.514.861.045.5853.53CHerat0.840.783.561.118.6439.23Chilas1.681.633.832.0917.9742.18D I Khan1.150.76.781.437.77874.59Dir0.770.365.70.874.02762.76Gilgit1.3915.591.7711.0361.57Kalam1.40.424.511.964.71149.7(b) Performance based on MAE and SEEMAESEEStationsBCRM-HO-HBCRM-HO-HBahawalnagar2.11.974.343.3221.7147.77Bahawalpur1.621.216.04213.3466.54DG Khan1.421.225.491.9813.4660.49Faisalabad1.971.545.162.6816.8356Jhang1.981.646.162.6915.8458.1Khanpore1.390.996.591.810.9772.59Lahore APT2.162.034.173.3322.3545.9Lahore PBO1.621.035.271.8711.457.98Mangla1.841.285.82.4114.1463.8M B din1.510.926.511.6310.1971.64Multan2.252.032.02322.2951.04Okara1.340.767.111.6912.6883.3R Y Khan1.611.216.582.213.3672.43Sargodha1.50.896.751.639.8874.31Sialkot1.610.986.451.7610.7770.99Jehalum1.450.916.491.6510.0371.43Badin4.745.011.785.7855.219.59Chhor5.956.21.997.2568.2321.9Gawadar3.253.62.593.5639.728.49Hyderabad8.588.832.479.6897.1727.18Jecobabad2.62.524.593.6227.7750.56Jiwini1.962.582.352.3428.4625.92Karachi4.264.761.134.6652.3812.44Lasbella5.455.442.86.3759.9230.89Mirpur khas2.162.045.323.0722.5458.59Nawabshah6.166.461.997.4871.0621.98Ormara1.852.183.722.1324.0440.92Panjgur6.026.041.776.8466.4519.52Pasni2.913.422.253.2937.7224.78Rohri2.422.634.253.0428.9546.79Shukkur3.433.424.054.6637.744.59Turbot9.479.452.1510.44103.923.71Astore0.810.823.230.999.1235.53Balakot1.060.615.891.186.7264.88Bunji0.860.514.861.045.5853.53CHerat0.840.783.561.118.6439.23Chilas1.681.633.832.0917.9742.18D I Khan1.150.76.781.437.77874.59Dir0.770.365.70.874.02762.76Gilgit1.3915.591.7711.0361.57Kalam1.40.424.511.964.71149.7


### 3.3. Estimation of Trend Analysis

Different statistical indicators detect the presence of trends in ET0 time series estimates. The statistical indicators for detecting trend and seasonal variation in the time series data includeR2, regression, and standard error of the estimate. The ET0 estimated by these temperature-dependent models are significantly different from those of FAO−56PM. According to R2, ET0 linear expressions can be applied between the ET0 of the FAO−56PM and HGR, which is considered for the dry, rainy, and annual seasons for all locations. Similarly, the expressions can be fitted between the ET0 of FAO−56PM and BCR, in the research areas for the dry, rainy, and annual seasons for 41 cities of Pakistan. The temperature-dependent methods were selected based on their evaluation during the season time step as displayed in Table 7, and so on for the remaining stations. However, regarding ET0 estimation by using both temperature-dependent methods, BCR model performs better as compared to HGR; see Figure 1.Table 7
Models performance according to decade (2006–2015) and SEE for Blaney-Criddle and Hargreaves models, where (∗) indicates statistically significant trends.
SeasonModelInterceptSlopeR2SEEBahawalnagarDryBCR2.460.720.680.82HGR3.373.140.977RainyBCR1.370.70.840.3∗HGR−1.052.390.932.62AnnualBCR0.421.050.880.88HGR−3.383.120.987.49BahawalpurDryBCR2.260.780.640.89HGR3.193.170.977.53RainyBCR1.790.630.840.4∗HGR−0.972.340.942.91AnnualBCR0.541.050.880.98HGR−3.173.120.988.04DG KhanDryBCR2.160.790.760.83HGR−3.373.140.977.63RainyBCR1.80.660.840.49HGR−0.492.320.952.77AnnualBCR0.820.990.920.96HGR−2.943.090.988.11FaisalabadDryBCR4.620.370.441.32HGR−0.342.640.976.9RainyBCR2.660.310.430.64HGR0.481.890.912.79AnnualBCR2.180.690.461.47HGR−1.732.750.97.44JhangDryBCR4.80.330.391.29HGR0.052.560.987.03RainyBCR2.730.280.410.69HGR0.551.880.912.91AnnualBCR2.350.640.431.47HGR−1.462.680.97.61KhanporeDryBCR4.460.410.431.19HGR−0.92.780.987.46RainyBCR2.590.320.360.6∗HGR0.411.970.913.27AnnualBCR1.750.750.441.36HGR−2.532.920.898.14Lahore APTDryBCR4.630.360.471.34HGR−0.562.680.976.79RainyBCR2.640.340.410.66HGR0.41.930.93.84AnnualBCR2.350.650.441.5HGR−1.662.750.97.36Lahore PBODryBCR4.660.380.531.53HGR−0.592.680.986.36RainyBCR2.80.330.440.6∗HGR0.451.940.922.67AnnualBCR2.460.690.481.66HGR−1.492.750.936.9BCR: Blaney Criddle, HGR: HargreavesFigure 1
Scatter diagram between monthlyET0 estimates by FAO−56PM against Original Hargreaves, BCR, and Modified Hargreaves in Bahawalnagar and Bahawalpur weather stations of Pakistan.


### 3.4. Seasons Based Model Performance

The accurate evapotranspirationET0 estimation plays a significant role in water management for agricultural purposes in Pakistan. Dividing the mean monthly and decade data according to dry and rainy seasons increases the performance of BCR and HGR in 41 stations of Pakistan. The assessment of HGR and BCR has been determined on the basis FAO−56 model for 41 stations of Pakistan.The study explores that the determination of evapotranspirationET0 on a seasonal basis, i.e., dry and rainy season, increases the efficiency of two temperature-dependent methods. In addition, wind speed and radiation seem to be significant variables affecting the ET0 estimation in the considered condition. The results indicate that the wind speed change between the dry and rainy seasons may affect ET0. In the semiarid climate of the selected stations, the current study results indicate that a seasonal approach must be used when estimating ET0 based on temperature-dependent models.It is found that the BCR model performs better than HGR in seasons as shown in Table8 (a) but Modified Hargreaves performs better in nonseasonal data as shown in Table 8 (b). Moreover, we need more attention for collecting the wind data since ET0 is sensitive to its effect. Finally, the estimated seasonal temperature-dependent models are suggested to be used in the estimation procedure of reference ET0 rather than other evapotranspiration models, which are determined by the FAO−56 model. The newly proposed technique can open a great wide horizon for resolving the problems of current irrigation management by accurate and timely estimation of ET0 in Pakistan.Table 8
Performance of models for selected stations of Pakistan for seasonal data (right side) and for nonseasonal data (left side).
(a) Performance of models for selected stations of Pakistan (nonseasonal data)ModelsStationsBCRM-HO-HBahawalnagarAverageGoodNormalBahawalpurAverageGoodNormalDG KhanGoodAverageNormalFaisalabadAverageGoodNormalJhangAverageGoodNormalKhanporeGoodAverageNormalLahore APTAverageGoodNormalLahore PBOGoodAverageNormalManglaGoodAverageNormalM B dinAverageGoodNormalMultanAverageGoodNormalOkaraAverageGoodNormalR Y KhanAverageGoodNormalSargodhaAverageGoodNormalSialkotGoodAverageNormalJehalumAverageGoodNormalBadinAverageGoodNormalChhorAverageGoodNormalGawadarAverageGoodNormalHyderabadAverageGoodNormalJecobabadGoodAverageNormalJiwiniAverageGoodNormalKarachiAverageGoodNormalLasbellaAverageGoodNormalMirpur khasNormalGoodAverageNawabshahAverageGoodNormalOrmaraGoodAverageNormalPanjgurAverageGoodNormalPasniNormalGoodAverageRohriAverageGoodNormalShukkurAverageGoodNormalTurbotAverageGoodNormalAstoreAverageGoodNormalBalakotGoodAverageNormalBunjiNormalGoodAverageCHeratAverageGoodNormalChilasAverageGoodNormalD I KhanAverageGoodNormalDirNormalGoodAverageGilgitAverageGoodNormalKalamAverageGoodNormal(b) Performance of models for selected stations of Pakistan (seasonal data)BCRHGRStationsDryRainyAnnualDryRainyAnnualBahawalnagarGoodGoodAverageAverageAverageNormalBahawalpurGoodGoodAverageAverageAverageNormalDG KhanAverageAverageGoodGoodGoodNormalFaislabadGoodGoodAverageAverageAverageNormalJhangGoodGoodAverageAverageAverageNormalKhanporeAverageAverageGoodGoodGoodNormalLahore APTGoodGoodAverageAverageAverageNormalLahore PBOAverageAverageGoodGoodGoodNormalManglaAverageAverageGoodGoodGoodNormalM B dinGoodGoodAverageAverageAverageNormalMultanGoodGoodAverageAverageAverageNormalOkaraGoodGoodAverageAverageAverageNormalR Y KhanGoodGoodAverageAverageAverageNormalSargodhaGoodGoodAverageAverageAverageNormalSialkotAverageNormalGoodGoodGoodNormalJehalumGoodGoodAverageAverageAverageNormalBadinGoodGoodAverageAverageAverageNormalChhorGoodGoodAverageAverageAverageNormalGawadarGoodGoodAverageAverageAverageNormalHyderabadGoodGoodAverageAverageAverageNormalJecobabadGoodAverageGoodGoodGoodNormalJiwiniGoodGoodAverageAverageAverageNormalKarachiGoodGoodAverageAverageAverageNormalLasbellaGoodGoodAverageAverageAverageNormalMirpur khasGoodGoodNormalNormalNormalAverageNawabshahGoodGoodAverageAverageAverageNormalOrmaraAverageNormalGoodGoodGoodNormalPanjgurGoodGoodAverageAverageAverageNormalPasniGoodGoodNormalNormalNormalAverageRohriGoodGoodAverageAverageAverageNormalShukkurGoodGoodAverageAverageAverageNormalTurbotGoodGoodAverageAverageAverageNormalAstoreGoodGoodAverageAverageAverageNormalBalakotNormalAverageGoodGoodGoodNormalBunjiGoodGoodNormalNormalNormalAverageCHeratGoodGoodAverageAverageAverageNormalChilasGoodGoodAverageAverageAverageNormalD I KhanGoodGoodAverageAverageAverageNormalDirGoodGoodNormalNormalNormalAverageGilgitGoodGoodAverageAverageAverageNormalKalamGoodGoodAverageGoodAverageNormal


### 3.5. Spatial Temporal Interpolation reference Evapotranspiration

#### 3.5.1. Area of Study

We observe the monthly evapotranspirationET0 data from all four provinces of Pakistan for the duration of 2006–2015 displayed in Figure 2. The data were collected from the Meteorological Department of Pakistan with their coordinates displayed in Figure 2(b). We use R package libraries geofd for flowing the data (by Fourier or B-splines basis) and geoR for fitting a variogram model to the predicted trace-variogram functions.Figure 2
Evapotranspiration estimate for 41 locations and grid points of observed data set according to decade 2006–2015 (120 months).Since the evapotranspirationET0 data set used for analysis is periodic, therefore, Fourier basis functions are the most accurate selection [23]. Nevertheless, for elaborative intention, we also use B-spline basis functions. We can estimate at a single location or multiple locations. Our study covers both scenarios. We ease the ET0 data by applying a B-splines basis and carry out the estimation for an unattended area (Arifwala, Leiah, Mianwali, Sargodha, and Gojra). Also, by applying a Fourier basis, we forecast evapotranspiration ET0 curves at ten randomly selected sites.


#### 3.5.2. B-Splines Basis Method

By applying the B-splines basis technique, we can estimate one space (location) at a time by applying various values of basis (n basis). These values were applied for the estimation of the space-time curve. Firstly, we create the B-spline basis in this technique for plotting the data. Four covariance models are used, that is, “spherical,” “exponential,” “Gaussian,” and “Mathern.” By applying B-splines to smooth the data, Figure 3 illustrates the results. We select only one station for estimation space (Arifwala), which is situated in the south Punjab of Pakistan. Figure 3(b) shows the sight of the selected location in a grid map.Figure 3
Smoothing data by using B-Spline Basis according to decade 2006–2015 (120 months).After essential computation, the last step in the space-time estimation by applying a B-spline basis is to check the manner for the chosen location at Arifwala, plotted in Figure4. The change of evapotranspiration for Arifwala is spotted by a line. Also, predictions of evapotranspiration for Mianwali, Leiah, Gojra, Sargodha, and Chakwal are also shown in Figure 4.Figure 4
Prediction of unobserved locations by using B-Spline Basis according to decade 2006–2015 (120 months).


#### 3.5.3. Fourier Basis Method

In addition to the B-spline method, the Fourier basis method can also be used to flow the space-time estimated curve. This method provides a more accurate result when the data is periodic. Several stations can be estimated at the same time using this procedure. Figure5 shows the variogram cloud for four chosen covariance models by applying the Fourier basis technique. We also apply the previous four covariance techniques, and the exponential covariance still has the least SEE. The estimation of ten casually chosen locations by this technique is shown in Figure 5(b).Figure 5
Smoothing data by using Fourier Basis and Prediction of unobserved Locations according to decade 2006–2015 (120 months).TheR package geofd is used for weather data to predict the unknown locations. This package contains functions for modeling the trace-variogram function and for carrying out spatial-temporal prediction using the method of ordinary kriging for functional data. The functional data analysis package provides techniques for smoothing data by using different basis functions.


## 3.1. Estimation of Evapotranspiration Models

The summary of averagedET0 estimation in (mm/day) by using FAO−56PM, BCR, and HGR for Bahawalnagar station is given in Table 3. These results are given in Tables 4 and 5.Table 3
Summary of annual evapotranspiration in (mm/day) by usingFAO−56PM, BCR, and HGR for Bahawalnagar.
JanFebMarAprMayJuneJulyAugSepOctNovDecFAO−56PM2.873.724.856.026.596.65.735.194.754.163.252.69BCR3.183.754.795.967.097.516.996.645.995.184.083.26HGR5.387.2910.6414.9917.9317.6414.1412.6111.510.217.495.49Table 4
Selected Meteorological stations of Pakistan with their coordinates and elevation.
Sr#Name of stationsLatitudeLongitudeElevation1Balakot34.5572.35995.42CHerat33.8271.5513723D.i.khan31.9270.86172.34Dir35.271.8513755Kalam35.8472.982103.16Chilas35.4274.112507Bunji35.6674.613728Gilgit32.9274.3424609Astore35.3374.9216810Chhor29.8869.72511Hyderabad25.3868.422812Jacobebad28.368.475513Naawabshah26.2568.373714Baddin24.6368.9915Karachi24.966.932216Shhukr27.768.8668.5217Rohri27.6668.96618Mirpur khas25.52691519Jiwani25.0661.815620Pajngur26.9764.196821Pasni25.2663.48922Gawadar25.1362.3329. 8623Lesbella26.2366.178724Ormara25.264.66225Turbat25.9863.0615526Bahawalnagar29.3373.85161.0527Bahawalpur29.3371.7811028D g khan30.0570.63148.129Faislabad31.4370.63148.130Jhang31.2672.3119031Jhelum32.9873.75287.1932Khanpore28.6570.6888.4133Lahore apt31.5874.4216.1534Lahore pbo31.5574.3321435Mangla33.0673.63283.336Mbdin32.9673.8252.9737Multan30.271.43121.9539Okara30.873.4318040Ry khan28.4370.3282.9341Sargodha32.0572.6618742Sialkot32.5174.53255.1Table 5
Extraterrestrial radiation for selected meteorological stations of Pakistan.
StationsJanFebMarAprMayJuneJulyAugSepOctNovDecBalakot18.723.729.936.14041.640.837.632.125.619.917.3CHerat18.723.729.936.14041.640.837.632.125.619.917.3Di khan19.924.830.736.54041.440.737.932.826.621.118.5Dir17.522.62935.74041.740.837.431.524.618.716.1Kalam17.522.62935.74041.740.837.431.524.618.716.1Chilas17.522.62935.74041.740.837.431.524.618.716.1Bunji17.522.62935.74041.740.837.431.524.618.716.1Gilgit19.924.830.736.54041.440.737.932.826.621.118.5Astore17.522.62935.74041.740.837.431.524.618.716.1Chhor21.125.831.436.84041.240.63833.427.622.219.8Hyderabad23.427.832.835.439.940.640.238.334.529.324.522.2Jacobebad22.325.832.237.14040.940.438.233.928.523.321Nawabshah23.427.832.835.439.940.640.238.334.529.324.522.2Baddin24.628.833.535.639.740.339.938.334.930.225.523.3Karachi24.628.833.535.639.740.339.938.334.930.225.523.3Shhukr22.325.832.237.14040.940.438.233.928.523.321Rohri22.325.832.237.14040.940.438.233.928.523.321MP khas23.427.832.835.439.940.640.238.334.529.324.522.2Jiwani23.427.832.835.439.940.640.238.334.529.324.522.2Panjgur23.427.832.835.439.940.640.238.334.529.324.522.2Pasni23.427.832.835.439.940.640.238.334.529.324.522.2Gawadr23.427.832.835.439.940.640.238.334.529.324.522.2Lasbella23.427.832.835.439.940.640.238.334.529.324.522.2Ormara23.427.832.835.439.940.640.238.334.529.324.522.2Turbat23.427.832.835.439.940.640.238.334.529.324.522.2Bahawalngr21.125.831.436.84041.240.63833.427.622.219.8Bahawalpur21.125.831.436.84041.240.63833.427.622.219.8D g khan21.125.831.436.84041.240.63833.427.622.219.8Faislabad19.924.830.736.54041.440.737.932.826.621.118.5Jhang19.924.830.736.54041.440.737.932.826.621.118.5Jhelum19.924.830.736.54041.440.737.932.826.621.118.5Khanpore22.325.832.237.14040.940.438.233.928.523.321Lahore apt19.924.830.736.54041.440.737.932.826.621.118.5Lahore Pbo19.924.830.736.54041.440.737.932.826.621.118.5Mangla18.723.729.936.14041.640.837.632.125.619.917.3Mbdin19.924.830.736.54041.440.737.932.826.621.118.5Multan21.125.831.436.84041.240.63833.427.622.219.8Okara21.125.831.436.84041.240.63833.427.622.219.8R Y khan22.325.832.237.14040.940.438.233.928.523.321Sargodha19.924.830.736.54041.440.737.932.826.621.118.5Sialkot19.924.830.736.54041.440.737.932.826.621.118.5


## 3.2. Performance Evaluation of Evapotranspiration Models

The results of Model Efficiency (ME) and Coefficient of Residual Mass (CRM) for different models are given in Table6(a). The results of other performance measures, i.e., Mean Absolute Error (MAE) and Standard Error of Estimate (SEE), are given in Table 6(b). From the results of these performance indicators, we can see that the HGR model has better performance as compared to the BCR and FPM models in nonseasonal data.Table 6
Performance Evaluation of Selected Models based on ME and CRM (Left panel), Performance Evaluation of Selected Models based on MAE and SEE (Right panel).
(a) Performance based on ME and CRMMECRMStationsBCRM-HO-HBCRM-HO-HBahawalnagar2.11.974.343.3221.7147.77Bahawalpur1.621.216.04213.3466.54DG Khan1.421.225.491.9813.4660.49Faislabad1.971.545.162.6816.8356Jhang1.981.646.162.6915.8458.1Khanpore1.390.996.591.810.9772.59Lahore APT2.162.034.173.3322.3545.9Lahore PBO1.621.035.271.8711.457.98Mangla1.841.285.82.4114.1463.8M B din1.510.926.511.6310.1971.64Multan2.252.032.02322.2951.04Okara1.340.767.111.6912.6883.3R Y Khan1.611.216.582.213.3672.43Sargodha1.50.896.751.639.8874.31Sialkot1.610.986.451.7610.7770.99Jehalum1.450.916.491.6510.0371.43Badin4.745.011.785.7855.219.59Chhor5.956.21.997.2568.2321.9Gawadar3.253.62.593.5639.728.49Hyderabad8.588.832.479.6897.1727.18Jecobabad2.62.524.593.6227.7750.56Jiwini1.962.582.352.3428.4625.92Karachi4.264.761.134.6652.3812.44Lasbella5.455.442.86.3759.9230.89Mirpur khas2.162.045.323.0722.5458.59Nawabshah6.166.461.997.4871.0621.98Ormara1.852.183.722.1324.0440.92Panjgur6.026.041.776.8466.4519.52Pasni2.913.422.253.2937.7224.78Rohri2.422.634.253.0428.9546.79Shukkur3.433.424.054.6637.744.59Turbot9.479.452.1510.44103.923.71Astore0.810.823.230.999.1235.53Balakot1.060.615.891.186.7264.88Bunji0.860.514.861.045.5853.53CHerat0.840.783.561.118.6439.23Chilas1.681.633.832.0917.9742.18D I Khan1.150.76.781.437.77874.59Dir0.770.365.70.874.02762.76Gilgit1.3915.591.7711.0361.57Kalam1.40.424.511.964.71149.7(b) Performance based on MAE and SEEMAESEEStationsBCRM-HO-HBCRM-HO-HBahawalnagar2.11.974.343.3221.7147.77Bahawalpur1.621.216.04213.3466.54DG Khan1.421.225.491.9813.4660.49Faisalabad1.971.545.162.6816.8356Jhang1.981.646.162.6915.8458.1Khanpore1.390.996.591.810.9772.59Lahore APT2.162.034.173.3322.3545.9Lahore PBO1.621.035.271.8711.457.98Mangla1.841.285.82.4114.1463.8M B din1.510.926.511.6310.1971.64Multan2.252.032.02322.2951.04Okara1.340.767.111.6912.6883.3R Y Khan1.611.216.582.213.3672.43Sargodha1.50.896.751.639.8874.31Sialkot1.610.986.451.7610.7770.99Jehalum1.450.916.491.6510.0371.43Badin4.745.011.785.7855.219.59Chhor5.956.21.997.2568.2321.9Gawadar3.253.62.593.5639.728.49Hyderabad8.588.832.479.6897.1727.18Jecobabad2.62.524.593.6227.7750.56Jiwini1.962.582.352.3428.4625.92Karachi4.264.761.134.6652.3812.44Lasbella5.455.442.86.3759.9230.89Mirpur khas2.162.045.323.0722.5458.59Nawabshah6.166.461.997.4871.0621.98Ormara1.852.183.722.1324.0440.92Panjgur6.026.041.776.8466.4519.52Pasni2.913.422.253.2937.7224.78Rohri2.422.634.253.0428.9546.79Shukkur3.433.424.054.6637.744.59Turbot9.479.452.1510.44103.923.71Astore0.810.823.230.999.1235.53Balakot1.060.615.891.186.7264.88Bunji0.860.514.861.045.5853.53CHerat0.840.783.561.118.6439.23Chilas1.681.633.832.0917.9742.18D I Khan1.150.76.781.437.77874.59Dir0.770.365.70.874.02762.76Gilgit1.3915.591.7711.0361.57Kalam1.40.424.511.964.71149.7


## 3.3. Estimation of Trend Analysis

Different statistical indicators detect the presence of trends in ET0 time series estimates. The statistical indicators for detecting trend and seasonal variation in the time series data includeR2, regression, and standard error of the estimate. The ET0 estimated by these temperature-dependent models are significantly different from those of FAO−56PM. According to R2, ET0 linear expressions can be applied between the ET0 of the FAO−56PM and HGR, which is considered for the dry, rainy, and annual seasons for all locations. Similarly, the expressions can be fitted between the ET0 of FAO−56PM and BCR, in the research areas for the dry, rainy, and annual seasons for 41 cities of Pakistan. The temperature-dependent methods were selected based on their evaluation during the season time step as displayed in Table 7, and so on for the remaining stations. However, regarding ET0 estimation by using both temperature-dependent methods, BCR model performs better as compared to HGR; see Figure 1.Table 7
Models performance according to decade (2006–2015) and SEE for Blaney-Criddle and Hargreaves models, where (∗) indicates statistically significant trends.
SeasonModelInterceptSlopeR2SEEBahawalnagarDryBCR2.460.720.680.82HGR3.373.140.977RainyBCR1.370.70.840.3∗HGR−1.052.390.932.62AnnualBCR0.421.050.880.88HGR−3.383.120.987.49BahawalpurDryBCR2.260.780.640.89HGR3.193.170.977.53RainyBCR1.790.630.840.4∗HGR−0.972.340.942.91AnnualBCR0.541.050.880.98HGR−3.173.120.988.04DG KhanDryBCR2.160.790.760.83HGR−3.373.140.977.63RainyBCR1.80.660.840.49HGR−0.492.320.952.77AnnualBCR0.820.990.920.96HGR−2.943.090.988.11FaisalabadDryBCR4.620.370.441.32HGR−0.342.640.976.9RainyBCR2.660.310.430.64HGR0.481.890.912.79AnnualBCR2.180.690.461.47HGR−1.732.750.97.44JhangDryBCR4.80.330.391.29HGR0.052.560.987.03RainyBCR2.730.280.410.69HGR0.551.880.912.91AnnualBCR2.350.640.431.47HGR−1.462.680.97.61KhanporeDryBCR4.460.410.431.19HGR−0.92.780.987.46RainyBCR2.590.320.360.6∗HGR0.411.970.913.27AnnualBCR1.750.750.441.36HGR−2.532.920.898.14Lahore APTDryBCR4.630.360.471.34HGR−0.562.680.976.79RainyBCR2.640.340.410.66HGR0.41.930.93.84AnnualBCR2.350.650.441.5HGR−1.662.750.97.36Lahore PBODryBCR4.660.380.531.53HGR−0.592.680.986.36RainyBCR2.80.330.440.6∗HGR0.451.940.922.67AnnualBCR2.460.690.481.66HGR−1.492.750.936.9BCR: Blaney Criddle, HGR: HargreavesFigure 1
Scatter diagram between monthlyET0 estimates by FAO−56PM against Original Hargreaves, BCR, and Modified Hargreaves in Bahawalnagar and Bahawalpur weather stations of Pakistan.


## 3.4. Seasons Based Model Performance

The accurate evapotranspirationET0 estimation plays a significant role in water management for agricultural purposes in Pakistan. Dividing the mean monthly and decade data according to dry and rainy seasons increases the performance of BCR and HGR in 41 stations of Pakistan. The assessment of HGR and BCR has been determined on the basis FAO−56 model for 41 stations of Pakistan.The study explores that the determination of evapotranspirationET0 on a seasonal basis, i.e., dry and rainy season, increases the efficiency of two temperature-dependent methods. In addition, wind speed and radiation seem to be significant variables affecting the ET0 estimation in the considered condition. The results indicate that the wind speed change between the dry and rainy seasons may affect ET0. In the semiarid climate of the selected stations, the current study results indicate that a seasonal approach must be used when estimating ET0 based on temperature-dependent models.It is found that the BCR model performs better than HGR in seasons as shown in Table8 (a) but Modified Hargreaves performs better in nonseasonal data as shown in Table 8 (b). Moreover, we need more attention for collecting the wind data since ET0 is sensitive to its effect. Finally, the estimated seasonal temperature-dependent models are suggested to be used in the estimation procedure of reference ET0 rather than other evapotranspiration models, which are determined by the FAO−56 model. The newly proposed technique can open a great wide horizon for resolving the problems of current irrigation management by accurate and timely estimation of ET0 in Pakistan.Table 8
Performance of models for selected stations of Pakistan for seasonal data (right side) and for nonseasonal data (left side).
(a) Performance of models for selected stations of Pakistan (nonseasonal data)ModelsStationsBCRM-HO-HBahawalnagarAverageGoodNormalBahawalpurAverageGoodNormalDG KhanGoodAverageNormalFaisalabadAverageGoodNormalJhangAverageGoodNormalKhanporeGoodAverageNormalLahore APTAverageGoodNormalLahore PBOGoodAverageNormalManglaGoodAverageNormalM B dinAverageGoodNormalMultanAverageGoodNormalOkaraAverageGoodNormalR Y KhanAverageGoodNormalSargodhaAverageGoodNormalSialkotGoodAverageNormalJehalumAverageGoodNormalBadinAverageGoodNormalChhorAverageGoodNormalGawadarAverageGoodNormalHyderabadAverageGoodNormalJecobabadGoodAverageNormalJiwiniAverageGoodNormalKarachiAverageGoodNormalLasbellaAverageGoodNormalMirpur khasNormalGoodAverageNawabshahAverageGoodNormalOrmaraGoodAverageNormalPanjgurAverageGoodNormalPasniNormalGoodAverageRohriAverageGoodNormalShukkurAverageGoodNormalTurbotAverageGoodNormalAstoreAverageGoodNormalBalakotGoodAverageNormalBunjiNormalGoodAverageCHeratAverageGoodNormalChilasAverageGoodNormalD I KhanAverageGoodNormalDirNormalGoodAverageGilgitAverageGoodNormalKalamAverageGoodNormal(b) Performance of models for selected stations of Pakistan (seasonal data)BCRHGRStationsDryRainyAnnualDryRainyAnnualBahawalnagarGoodGoodAverageAverageAverageNormalBahawalpurGoodGoodAverageAverageAverageNormalDG KhanAverageAverageGoodGoodGoodNormalFaislabadGoodGoodAverageAverageAverageNormalJhangGoodGoodAverageAverageAverageNormalKhanporeAverageAverageGoodGoodGoodNormalLahore APTGoodGoodAverageAverageAverageNormalLahore PBOAverageAverageGoodGoodGoodNormalManglaAverageAverageGoodGoodGoodNormalM B dinGoodGoodAverageAverageAverageNormalMultanGoodGoodAverageAverageAverageNormalOkaraGoodGoodAverageAverageAverageNormalR Y KhanGoodGoodAverageAverageAverageNormalSargodhaGoodGoodAverageAverageAverageNormalSialkotAverageNormalGoodGoodGoodNormalJehalumGoodGoodAverageAverageAverageNormalBadinGoodGoodAverageAverageAverageNormalChhorGoodGoodAverageAverageAverageNormalGawadarGoodGoodAverageAverageAverageNormalHyderabadGoodGoodAverageAverageAverageNormalJecobabadGoodAverageGoodGoodGoodNormalJiwiniGoodGoodAverageAverageAverageNormalKarachiGoodGoodAverageAverageAverageNormalLasbellaGoodGoodAverageAverageAverageNormalMirpur khasGoodGoodNormalNormalNormalAverageNawabshahGoodGoodAverageAverageAverageNormalOrmaraAverageNormalGoodGoodGoodNormalPanjgurGoodGoodAverageAverageAverageNormalPasniGoodGoodNormalNormalNormalAverageRohriGoodGoodAverageAverageAverageNormalShukkurGoodGoodAverageAverageAverageNormalTurbotGoodGoodAverageAverageAverageNormalAstoreGoodGoodAverageAverageAverageNormalBalakotNormalAverageGoodGoodGoodNormalBunjiGoodGoodNormalNormalNormalAverageCHeratGoodGoodAverageAverageAverageNormalChilasGoodGoodAverageAverageAverageNormalD I KhanGoodGoodAverageAverageAverageNormalDirGoodGoodNormalNormalNormalAverageGilgitGoodGoodAverageAverageAverageNormalKalamGoodGoodAverageGoodAverageNormal


## 3.5. Spatial Temporal Interpolation reference Evapotranspiration

### 3.5.1. Area of Study

We observe the monthly evapotranspirationET0 data from all four provinces of Pakistan for the duration of 2006–2015 displayed in Figure 2. The data were collected from the Meteorological Department of Pakistan with their coordinates displayed in Figure 2(b). We use R package libraries geofd for flowing the data (by Fourier or B-splines basis) and geoR for fitting a variogram model to the predicted trace-variogram functions.Figure 2
Evapotranspiration estimate for 41 locations and grid points of observed data set according to decade 2006–2015 (120 months).Since the evapotranspirationET0 data set used for analysis is periodic, therefore, Fourier basis functions are the most accurate selection [23]. Nevertheless, for elaborative intention, we also use B-spline basis functions. We can estimate at a single location or multiple locations. Our study covers both scenarios. We ease the ET0 data by applying a B-splines basis and carry out the estimation for an unattended area (Arifwala, Leiah, Mianwali, Sargodha, and Gojra). Also, by applying a Fourier basis, we forecast evapotranspiration ET0 curves at ten randomly selected sites.


### 3.5.2. B-Splines Basis Method

By applying the B-splines basis technique, we can estimate one space (location) at a time by applying various values of basis (n basis). These values were applied for the estimation of the space-time curve. Firstly, we create the B-spline basis in this technique for plotting the data. Four covariance models are used, that is, “spherical,” “exponential,” “Gaussian,” and “Mathern.” By applying B-splines to smooth the data, Figure 3 illustrates the results. We select only one station for estimation space (Arifwala), which is situated in the south Punjab of Pakistan. Figure 3(b) shows the sight of the selected location in a grid map.Figure 3
Smoothing data by using B-Spline Basis according to decade 2006–2015 (120 months).After essential computation, the last step in the space-time estimation by applying a B-spline basis is to check the manner for the chosen location at Arifwala, plotted in Figure4. The change of evapotranspiration for Arifwala is spotted by a line. Also, predictions of evapotranspiration for Mianwali, Leiah, Gojra, Sargodha, and Chakwal are also shown in Figure 4.Figure 4
Prediction of unobserved locations by using B-Spline Basis according to decade 2006–2015 (120 months).


### 3.5.3. Fourier Basis Method

In addition to the B-spline method, the Fourier basis method can also be used to flow the space-time estimated curve. This method provides a more accurate result when the data is periodic. Several stations can be estimated at the same time using this procedure. Figure5 shows the variogram cloud for four chosen covariance models by applying the Fourier basis technique. We also apply the previous four covariance techniques, and the exponential covariance still has the least SEE. The estimation of ten casually chosen locations by this technique is shown in Figure 5(b).Figure 5
Smoothing data by using Fourier Basis and Prediction of unobserved Locations according to decade 2006–2015 (120 months).TheR package geofd is used for weather data to predict the unknown locations. This package contains functions for modeling the trace-variogram function and for carrying out spatial-temporal prediction using the method of ordinary kriging for functional data. The functional data analysis package provides techniques for smoothing data by using different basis functions.


## 3.5.1. Area of Study

We observe the monthly evapotranspirationET0 data from all four provinces of Pakistan for the duration of 2006–2015 displayed in Figure 2. The data were collected from the Meteorological Department of Pakistan with their coordinates displayed in Figure 2(b). We use R package libraries geofd for flowing the data (by Fourier or B-splines basis) and geoR for fitting a variogram model to the predicted trace-variogram functions.Figure 2
Evapotranspiration estimate for 41 locations and grid points of observed data set according to decade 2006–2015 (120 months).Since the evapotranspirationET0 data set used for analysis is periodic, therefore, Fourier basis functions are the most accurate selection [23]. Nevertheless, for elaborative intention, we also use B-spline basis functions. We can estimate at a single location or multiple locations. Our study covers both scenarios. We ease the ET0 data by applying a B-splines basis and carry out the estimation for an unattended area (Arifwala, Leiah, Mianwali, Sargodha, and Gojra). Also, by applying a Fourier basis, we forecast evapotranspiration ET0 curves at ten randomly selected sites.


## 3.5.2. B-Splines Basis Method

By applying the B-splines basis technique, we can estimate one space (location) at a time by applying various values of basis (n basis). These values were applied for the estimation of the space-time curve. Firstly, we create the B-spline basis in this technique for plotting the data. Four covariance models are used, that is, “spherical,” “exponential,” “Gaussian,” and “Mathern.” By applying B-splines to smooth the data, Figure 3 illustrates the results. We select only one station for estimation space (Arifwala), which is situated in the south Punjab of Pakistan. Figure 3(b) shows the sight of the selected location in a grid map.Figure 3
Smoothing data by using B-Spline Basis according to decade 2006–2015 (120 months).After essential computation, the last step in the space-time estimation by applying a B-spline basis is to check the manner for the chosen location at Arifwala, plotted in Figure4. The change of evapotranspiration for Arifwala is spotted by a line. Also, predictions of evapotranspiration for Mianwali, Leiah, Gojra, Sargodha, and Chakwal are also shown in Figure 4.Figure 4
Prediction of unobserved locations by using B-Spline Basis according to decade 2006–2015 (120 months).


## 3.5.3. Fourier Basis Method

In addition to the B-spline method, the Fourier basis method can also be used to flow the space-time estimated curve. This method provides a more accurate result when the data is periodic. Several stations can be estimated at the same time using this procedure. Figure5 shows the variogram cloud for four chosen covariance models by applying the Fourier basis technique. We also apply the previous four covariance techniques, and the exponential covariance still has the least SEE. The estimation of ten casually chosen locations by this technique is shown in Figure 5(b).Figure 5
Smoothing data by using Fourier Basis and Prediction of unobserved Locations according to decade 2006–2015 (120 months).TheR package geofd is used for weather data to predict the unknown locations. This package contains functions for modeling the trace-variogram function and for carrying out spatial-temporal prediction using the method of ordinary kriging for functional data. The functional data analysis package provides techniques for smoothing data by using different basis functions.


## 4. Conclusion

In this study, the primary aim was to explain spatiotemporal interpolation and estimate referenceET0 techniques. Based on the analysis of our study, the air temperature modified Hargreaves model was used instead of the Hargreaves-Samani model since it mainly underestimates ET0 in tropical humid and hot areas. A value of 0.382 for the Hargreaves exponent is suggested instead of its actual value, i.e., 0.5, in modified Hargreaves for Pakistan. The performance of models is measured using modeling efficiency, coefficient of residual mass, mean absolute error, and standard error of the estimate. In this study, we inferred that purely temperature-dependent models such as BCR worked really well also if data of significant meteorological parameters such as Relative humidity, Wind speed, and Solar radiation is limited. Moreover, our suggested modified Hargreaves model is convenient to use ET0 assessment. We used radiation-based techniques in all Pakistan’s stations, i.e., FAO−56PM, and our suggested method for Pakistan gives accurate estimated results. From several models contrary to the typical reference model FAO−56PM, our suggested modified Hargreaves model is found to be on higher value for the coefficient of determination R2, while the actual Hargreaves-Samani model is found to be on the lower value. The ET0 changes in the average seasonal time magnitude from 2006 to 2015 worked out by applying BCR and HGR models in Pakistan. The upward and downward changes are observed in the seasonal ET0 in selected regions of Pakistan.Cressie and Wikle [9] suggested kriging estimator(s) for carrying out spatial forecasting of functional data. In this way, a balanced step, through B-splines or Fourier basis functions, is implemented at first. A technique for setting up spatial dependence between functions is then suggested, along with a method for achieving a spatial forecast for an unusual place. We predict one unsampled location at once, that is, Arifwala, Leiah, Gojra, Mianwali, Chakwal, and Sargodha, by using the technique of B-Spline basis because this predicts only one location at a time. As our data show periodic patterns, we use Fourier basis instead of B-spline based for predicting more than one unsampled location and smoothing the trends, as well as Ordinary Kriging for Functional Data (OKFD) for spatiotemporal prediction of functional data.



---

*Source: 5488725-2022-04-28.xml*