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.
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, , and water deficits due to public involvement with climate change.
Donatelli et al. [3] showed that accurate 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 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 through . 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 ; see Shahid [6].
The irrigation system of Pakistan is the most significant among the globe. Hence, an appropriate assessment of 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 reference 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 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 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 special R 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 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 evapotranspiration 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 reference . 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 precise 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 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 equation is given bywhere = the latent heat of vaporization , the evapotranspiration , the slope of saturation vapor pressure temperature curve , the psychometric constant , the Net Radiation , G = soil heat flux , mean air density , specific heat of air , and are the aerodynamic and surface resistances .
In the 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.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 in is as follows:
2.2.3. Blaney-Criddle (BCR)
The major drawback of the 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 , 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 bywhere = Average temperature . 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 by . 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:where Y = , , and .
We have measured the values of , and (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 is found to be appropriate for all stations. The resulting modified Hargreaves equation becomes
2.3. Performance Evaluation Criteria
The BCR, HGR, and modified Hargreaves conclusions are compared with the prediction technique used in this analysis as the reference 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.
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 observed 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:where n = number of observations. = is the estimated models [HGR, BCR and M-H]. = observed variable or model . = 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 variable , respectively. It can be shown mathematically aswhere n = number of observations. = is the estimated models [HGR, BCR and M-H]. = observed variable or model. .
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 aswhere N = total number of observations. = calculated with FAO-56(PM). = 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 involved 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 aswhere N = total number of observations. = calculated with FAO-56(PM). = 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 Table 2.
2.5. Trend Analysis
The mean decade 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 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 standard model. The standard error of estimate (SEE) and the coefficients of determination 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 predicted values. The equations for standard error of estimate (SEE) and coefficients of determination are given as follows:
Here, represents the observed , is the regression estimated for the value, and represent the mean values of the corresponding variable, and represents the number of observations taken into account. A simple linear regression with zero intercepts and unit slope and and 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 is the observation of function variables distributed as Y. Let . We implement the functional data constituted as
Note that with the inner product explains a Euclidean space. Following Giraldo et al. [13], we elaborate a usual functional process as , where d is usually taken as 2, and is a function variable for any . Let be the hypothetical points in and predict that we can examine a recognition of the functional usual process at these “n” sites, . The Ordinary Kriging for the functional data is a geostatistical method for forecasting , the functional usual process at , where 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) and for all and all .(2), where .(3), .
These assumptions imply that = and . According to Giraldo et al. [15], the OKFD predictor is defined by
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, , aswhere , .
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 . In a typical geostatistical setting, we suppose that the examinations are recognition of a usual field . The kriging predictor is defined as and the BLUP is achieved by reducingsubject to . Furthermore, in multivariate geo-statistics, the data consists of ; that is, we have observations of a spatial vector-valued process , where and . In this context, is a matrix, and the BLUP of m variables at an unsampled region can be achieved by reducing
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].which after some simplification can be given aswhere is the Lagrange multiplier used to consider the unbiasedness limitation. Reducing equation (18) in context with 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:
The functionis 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 has a constant average function over .
By using Fubini’s algorithm,where with . Then, adherence of the typical method-of-moments for this predicted value shows the following predictor:where , and 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 , with being a small value.
Once we have predicted the trace-variogram for K values , a parametric technique 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
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 variable , and we need to forecast a variable at an unsampled region. In the case of a space function, is a function, so we are estimating a function.
3. Results and Discussions
3.1. Estimation of Evapotranspiration Models
The summary of averaged estimation in (mm/day) by using , BCR, and HGR for Bahawalnagar station is given in Table 3. These results are given in Tables 4 and 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 Table 6(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.
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 include , regression, and standard error of the estimate. The estimated by these temperature-dependent models are significantly different from those of . According to , linear expressions can be applied between the of the and HGR, which is considered for the dry, rainy, and annual seasons for all locations. Similarly, the expressions can be fitted between the of 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 estimation by using both temperature-dependent methods, BCR model performs better as compared to HGR; see Figure 1.
3.4. Seasons Based Model Performance
The accurate evapotranspiration 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 model for 41 stations of Pakistan.
The study explores that the determination of evapotranspiration 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 estimation in the considered condition. The results indicate that the wind speed change between the dry and rainy seasons may affect . In the semiarid climate of the selected stations, the current study results indicate that a seasonal approach must be used when estimating based on temperature-dependent models.
It is found that the BCR model performs better than HGR in seasons as shown in Table 8 (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 is sensitive to its effect. Finally, the estimated seasonal temperature-dependent models are suggested to be used in the estimation procedure of reference rather than other evapotranspiration models, which are determined by the 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 in Pakistan.
3.5. Spatial Temporal Interpolation reference Evapotranspiration
3.5.1. Area of Study
We observe the monthly evapotranspiration 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.
Since the evapotranspiration 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 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 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 ( 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.
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 Figure 4. 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.
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. Figure 5 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).
The R 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 reference 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 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 assessment. We used radiation-based techniques in all Pakistan’s stations, i.e., , and our suggested method for Pakistan gives accurate estimated results. From several models contrary to the typical reference model , our suggested modified Hargreaves model is found to be on higher value for the coefficient of determination , while the actual Hargreaves-Samani model is found to be on the lower value. The 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 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.
Data Availability
The data are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.