Design and analysis of statistical probability distribution and non-parametric trend analysis for reference evapotranspiration

Accurate estimates of reference evapotranspiration are critical for water-resource management strategies such as irrigation scheduling and operation. Therefore, knowledge of events such as spatial and temporal reference evapotranspiration (ETo) and their related principle of statistical probability theory plays a vital role in amplifying sustainable irrigation planning. Spatiotemporal statistical probability distribution and its associated trends have not yet has explored in Pakistan. In this study, we have two objectives: (1) to determine the most appropriate statistical probability distribution that better describes ETo on mean monthly and seasons wise estimates for the design of irrigation system in Khyber Pakhtunkhwa, and (2) to check the trends in ETo on a monthly, seasonal, and annual basis. To check the ETo trends, we used the modified version of the Mann-Kendall and Sen Slope. We used Bayesian Kriging for spatial interpolation and propose a practical approach to the design and study of statistical probability distributions for the irrigation system and water supplies management. Also, the scheme preeminent explains ETo, on a monthly and seasonal basis. The statistical distribution that showed the best fit ETo result occupying 58.33% and 25% performance for the design of irrigation scheme in the entire study region on the monthly level was Johnson SB and Generalized Pareto, respectively. However, according to the Anderson-Darling (AD) and Kolmogorov–Smirnov (KS) goodness of fit measure, seasonal ETo estimates were preferably suited to the Burr, Johnson SB & Generalized Extreme Value. More research work must be conduct to assess the significance of this study to other fields. In conclusion, these findings might be helpful for water resource management and policymaker in future operations.


INTRODUCTION
Reference evapotranspiration is one of the most crucial functions in the Hydrological progression, and the climate is affected by this phenomenon (Allen et al., 1998). The term evapotranspiration involves evaporation from the soil, and plant surfaces include transpiration from plant leaves. Evaporation involves converting liquid water into water vapor and removing it from evaporation's surface (Allen et al., 1998). Recently, water crisis management planning, agriculture water management, and water resource projects were scheduled based on ET o . An estimate of irrigation water requirement is vital for water planning and operation of a sophisticated irrigation system worldwide. Therefore, ET o estimates have frequently been using in hydrology, agriculture water management, irrigation, drainage engineering designing, water crisis supervision, scheduling, and operational supervision (Valipour, 2014a;Valipour, 2014b;Valipour, 2015b;Valipour, 2015c).
Besides, several researchers showed that an estimate of irrigation water requirement is significant for the function of an intricate irrigation system in many parts of the world (Tabari, Grismer & Trajkovic, 2013;Tabari & Talaee, 2011;Thepadia & Martinez, 2012). An accurate estimation of the reference or crop evapotranspiration (ET o , ET c ) can utilize the accessible water resources such as stored water for the arid period. To estimate ET c in the first stage, ET o is calculated and then multiplied by the harvest coefficient (K C ) (Allen et al., 1998;Doorenbos, Pruitt & Aboukhaled, 1997). Further, Allen et al. (1998) define crop evapotranspiration as an evapotranspiration rate with an inferred crop height from the hypothetical crop and a fixed resistance to the canopy (70 s per meter). The albedo would cautiously acknowledge the proximity to evapotranspiration from a widespread green grass surface of equal height, aggressive growth, and strictly shading the soil. Likewise, many well-known researchers discussed and strongly suggested the FAO-56 (PM) technique as a standard procedure for accurate assessment and estimation of ET o , if lysimeter data is unavailable or limited on reference evapotranspiration (Gavilán et al., 2006;Irmak, Allen & Whitty, 2003;Nandagiri & Kovoor, 2006;Trajkovic, 2007;Utset et al., 2004).
On the other hand, numerous studies have shown that the temperature-based Hargreaves-Samani formula could provide accurate measures of ET o for five or longer days in various corners of the world (Droogers & Allen, 2002;Hargreaves & Allen, 2003;Jensen et al., 1997;Landeras, Ortiz-Barredo & López, 2008;Martí et al., 2015). Ultimately, the Hargreaves-Samani system's implied when partial data such as relative humidity, solar radiation, and wind speed is sufficient. Extensive research has been carried out in different parts of the world for estimating ET o and comparison models by using limited climatic data. They compared the Makkink, Turc, Hargreaves, and Priestley-Taylor methodologies for estimation. The studies revealed that the Radiation-based approach was more precise and proper than the Hargreaves method (Bois et al., 2008;Tabari, Grismer & Trajkovic, 2013;Valipour, 2015a;Xu & Singh, 2002). However, some researchers also observed that the suggested ET o values showed overestimation in damp conditions (Itenfisu et al., 2003;Temesgen et al., 2005;Trajkovic, 2005;Valipour, 2017). Similarly, several research studies also showed underestimated ET o in arid and semi-arid regions in diverse locations (Azhar & Perera, 2011;Benli et al., 2010;Droogers & Allen, 2002;Khoob, 2008). Consequently, in arid and semi-arid environments, local calibration and justification of such concepts are considerable (Da Cunha, Magalhaes & De Castro, 2013).
Moreover, several research findings associated with environmental events have required also evaluating probability distributions and improved clustering procedure based on energy consumption, prediction, and risk assessment in networks (Blain, 2012;Blain & Meschiatti, 2014;Huang et al., 2017;Nam & Choi, 2014;Rahman et al., 2013;Wu, Xiong & Wu, 2017;Zhang et al., 2015). As a result, the elementary probability theory of extreme events has been developing and used extensively. Thus, extreme events' statistical probability distribution modeling remains of great concern and plays a vital part in water resource supervision in different conditions (Katz, 2010).
In contrast, climate change has been detected not only in individual parameters such as temperature and precipitation but also in integrated parameters such as ET o (Ma et al., 2017). The most influential factor for climate change and the water cycle is evapotranspiration, as it is the only linkage between the water balance and the surface energy balance of the land (Wang et al., 2014). Some well-known authors worldwide suggested that knowledge about best fit statistical probability model in site ease is a correct hydrological prediction in the dimensioning drainage and irrigation planning (Fernandes, Wolff & Folegatti, 2019;Nam & Choi, 2014;Uliana et al., 2017;Yoo, Choi & Jang, 2008).
Besides, climate change is almost proven globally; the pattern of evapotranspiration is not apparent. ET o trends can increase or decrease depending on climatic conditions and regions. However, shifts in hydrological cycle components, such as ET o , display multiple increasing and reducing patterns (Donohue, McVicar & Roderick, 2010;Gao et al., 2006;Liu et al., 2010;Shenbin, Yunfeng & Thomas, 2006;Xu et al., 2006;Zhang et al., 2009). Also, numerous researches in several parts of the globe have shown a decline in evapotranspiration over the past decades. However, some have identified the contrary phenomenon, which is a rise in the trend of evapotranspiration (Croitoru et al., 2013;Mohsin & Lone, 2020;Pandey & Pandey, 2013;Shan et al., 2015;Shi et al., 2017;Talaee, Some'e & Ardakani, 2014;Vicente-Serrano et al., 2014a;Vicente-Serrano et al., 2014b;Wang et al., 2014). Therefore, most of the works conducted to identify ET o trend fluctuations checked by the traditional non-parametric Mann-Kendall and Sen Slope estimator tools. In the meantime, those experiments have not measured the effect of coefficients of autocorrelation with different lags. However, only in specific studies is the first order (lag-1) autocorrelation regarded (Amirataee, Montaseri & Sanikhani, 2016;Palizdan et al., 2014;Sayemuzzaman, Jha & Mekonnen, 2015;Tabari, Nikbakht & Talaee, 2012). Consequently, it may be crucial to exclude autocorrelation's effect to ascertain the changes to reference evapotranspiration.
Hence, in the present study, we use the updated version of the Mann-Kendall (Hamed & Rao, 1998) to check monthly, seasonal, and annual ET o pattern variations in the Khyber Pakhtunkhwa. Therefore, it is essential to note reference evapotranspiration patterns and their position in regional dry and wet environments. It could provide a scientific foundation for managing and distributing regional water resources and technical judgment on flood and drought disaster prevention (Ma et al., 2017).
In South Asia, Pakistan is situated at the intersection of Central Asia and the Middle East, enabling its place of tremendous value in the world with the most incredible scenery, stretching from the Arabian Sea, its southern boundary, to the stunning and highest mountain ranges of the Himalayas-Karakoram-Hindukush (HKH) in the north, also recognized as the third pole of the world (Yao et al., 2012). According to the latest census estimate, its total population is (207 774 520) (Statistics, 2017). On the other hand, many working classes are related either directly or indirectly to the agriculture sector (Rehman & Chandio, 2019;Rehman et al., 2015;Rehman et al., 2016).
Khyber Pakhtunkhwa (KP) is one of the four Pakistani provinces situated in the state's northwestern areas. KP is positioned predominantly on the Iranian plateau and includes the intersection where the slopes of the Hindu Kush Mountains approach the Indus-watered hills in South Central Asia on the Eurasian plate. The economy is agricultural in Khyber Pakhtunkhwa; therefore, irrigation is about one-third of the farming area. The legendary Khyber Pass links the province to Afghanistan, while Circle Bakote Abbottabad Kohala Bridge is a central crossing point in the east over the Jhelum River. Kabul, Swat, Chitral, Kunhar, Siran, Panjgora, Bara, Kurram, Haroo, Gomal, and Zhob are the main basins that traverse the province. The province could be classified physically into two zones, including the northern one extending from the Hindu Kush range to the Peshawar basin border and the southern one extending from Peshawar to the Derajat basin. The north area is cold and snowy with heavy rainfall in winters and good summers with the Peshawar basin exemption, hot in summer, and cold in winter with the Peshawar basin exemption.
In this present research paper, the two practical approaches, like Kolmogorov-Smirnov and Anderson-Darling test, are used to propose statistical probability distributions for the irrigation system and water resource management. This highlights and explores the probability distribution that best fits and describes the reference evapotranspiration ET o monthly and the seasonal basis using the most appropriate estimation method FAO-56(PM) and temperature-based estimation equation Hargreaves-Samani ET o for our selected study area known as Khyber Pakhtunkhwa, Pakistan. Focusing on irrigation dimension and considering the water resources' significance and the need for fair water use in agriculture planning in the study region. The reference evapotranspiration trends analysis has done by using modified Mann-Kendall (Hamed & Rao, 1998) and the Sen Slope estimator on a monthly, seasonal, and annual time basis. The rest of this paper is a brief description of the study area, the statistical probability distributions, and the methods used to estimate evapotranspiration in 'Methods and Material'. 'Selection Criteria for Statistical Probability Distribution' describes selected locations' results and a brief description of the goodness of fit and Statistical probability distribution selection criteria. The results and discussion have described in 'Results and Discussion'. The summary and conclusions have discussed in 'Summary and Conclusion'.

Site description and meteorological data
Khyber Pakhtunkhwa has a territory of 101, 741 km 2 with about approximately 36 million population, making KP the third largest provincial economy of Pakistan. The entire region varies based on topography, from the north's green fields to the south's dry, rocky zones. With profoundly warm summers to cold winters, the climate can be extreme. Despite these weather extremes, agriculture in the study area matters much and feasible. Due to its diverse agroecological diversity, KP is especially vulnerable to the effects of climate change (Gul et al., 2018). In this research study, twelve climate stations distributed in the Khyber Pakhtunkhwa zone has appointed for our exploration. The cities concerned are Balakot, Cherat, Chitral, Kakul, Parachinar, Dir, Drosh, DI Khan, Kohat, Peshawar, Risalpur, and Saidu Sharif (Gul et al., 2020). A set of average monthly data collection of wind speed, air temperature, relative humidity and precipitation at these sites under standard grass height. We used Pakistan Meteorological Department (PMD) well-irrigated climate data to obtain the FAO-56(PM) and Hargreaves-Samani dependent ET o values from twelve centers covering the period from 2000 to 2009. The overview of the various climate stations is presented in Table 1 and illustrated graphically in Fig. 1.

Methodology for proposed analysis scheme
This analysis used a ten-year data set of reference evapotranspiration estimates using the FAO-56(PM) and Hargreaves-Samani (HS) estimation methods. Also, the estimated ET o adjusted for Twenty-one Statistical probability distribution, as follow Burr, Johnson SB, Generalized Extreme Value (GEV), Log-Pearson-3(LP3), Normal, Logistic, Burr(4P), Weibull(3P), Generalized Extreme Value(4P), Gamma, Lognormal-3P (LN(3P)), Gumbel Min, Gen. Gamma, Gamma (3P), Cauchy, Lognormal, Weibull, Rayleigh (2P), Gumbel Max, Rayleigh and Generalized Pareto. Moreover, we used the Maximum Goodness method to adjust the best fit selected probability distributions through Kolmogorov-Smirnov and Anderson-Darling goodness of fit statistic. For this purpose, the maximum likelihood (MLE) has used the Bayesian Information Criteria (BIC) and Akaike Information Criteria (AIC) to estimate the adapted probability distribution parameters. Therefore, to contribute to all studied regions and attain parameter values, the required distribution, which provided excellent adjustment to most ET o data, was the most appropriate and picked. We used Kolmogorov-Smirnov and Anderson-Darling goodness of fit statistics to confirm whether the distribution used for all the data series conformed to the data. Furthermore, we have used two data analysis tools such as MS Excel and Evapotranspiration packages, an open-source R (R.3.6.1) statistical software program for computation of FAO-56(PM) and Hargreaves-Samani methodologies. Also, for adjusting different distribution models in Non-censored or Censored data, the fit ''Distrplus'' package was used in the statistical software package (R 3.6.1) and R Studio. In this context, we used Arc GIS software for generating the spatial distribution maps, while Bayesian Kriging strategies have been used for interpolation purposes. The methodologies are discussed step by step in the following sub-sections.

Penman-Montieth equation
The Food and Agriculture Organization (FAO) of the United Nations and ASCE have advised the Penman-Montieth equation (American Society of Civil Engineers) as the single and standard technique to estimate reference evapotranspiration and to assess other equations in a situation where the appropriate information is available (Allen et al., 1998). This strategy is physically dependent and can usually be enforced without any further alteration of input parameters (Gong et al., 2006). The FAO-56(PM) standard scheme (Allen et al., 1998) has been used in the current study to measure ETo for the recent data analysis; e a = is the actual vapor pressure [kPa].
= is the slope of the saturation vapor pressure-temperature curve [kPa ( o C −1 )]. T = is the daily air temperature [ o C], U 2 = is the wind speed at 2 m hight [ms −1 ]. Allen et al. (1998) known as the Penman-Monteith approach, proposed a complete collection of equations according to the available climate data set and time phase calculation.

Hargreaves-Samani equation
The Hargreaves and Samani (1985) equation is a famous symbolize edition to estimates ET o values based on temperature (Hargreaves & Allen, 2003). The Hargreaves-Samani procedure is theoretically comparable editions, too Hargreaves & Samani (1982). The Hargreaves-Samani process intended to make the earlier editions transparent by reducing the number of measured air temperature data and using extraterrestrial radiation (Ra) as a replacement for sunlight or radiation data (Hargreaves & Allen, 2003). Therefore, the FAO adopted the Hargreaves-Samani estimation method for ETo estimation, whereas air temperature alone is merely accessible (Allen et al., 1998;Hargreaves & Allen, 2003). The mathematical form of the Hargreaves-Samani equation given by Allen et al. (1998): Where; ET HS =ET estimated by using the above HS equation (mm per day) R a =extraterrestrial radiation (mm per day) T mean = mean air temperature in degree Celsius T max = daily maximum temperature in degree Celsius T min = daily minimum temperature in degree Celsius

SELECTION CRITERIA FOR STATISTICAL PROBABILITY DISTRIBUTION
In this paper, our projected scheme for the design and analysis of statistical probability distribution for an irrigation system and water resource management tries to discover the related logic based on the appropriate test procedure's goodness. The process to determine whether a sample of n observations (x 1, x 2,... x n ) can be considered a sample from a given specified distribution or not is known as the Goodness-of-fit test.
There is a wide range of tests available in the literature to determine whether a sample could have been drawn from a specific distribution. The test are Kolmogorov-Smirnov (KS), Anderson-Darling (AD), Chi-sq (χ2 test), Shapiro-Wilk (SW), Hosmer-Lemeshow (HL), and Henze-Zirkler (HZ). However, in the current investigation, the Kolmogorov-Smirnov and Anderson-Darling are used to select probability distributions. The KS and AD tests' principal advantage is their ability to detect variations in the probability distribution's overall shape because the test can detect differences across all scales (Darling, 1957). Thus, it is also suitable for small samples (Pettitt, 1976). Strength in the area of engineering is that For that purpose, both KS and AD are appropriate in the following perspectives: (1) when disparities are apparent, but normality is maintained; (2) when the sample size is small; (3) when variances are similar but not symmetrical; (4) when shift between two distributions does not affect only the upper or lower extremities.
In essence, we calculate the ''distance'' between the sample of experience data and the distribution we are examining, such as test statistics. We then measure the distance to a specific threshold value, such as a critical value. Only if the test statistic is smaller than the critical value with the significance amount (alpha = 0.05) can the configured probability distribution be considered a good fit. Since the statistics on the goodness-of-fit test describe the distance between the data and the distributions fitted. Thus, the distribution with the lowest statistical value is the one that best matches the distribution of probability. Moreover, the top three out of twenty-one distributions are taken for every single station.
The AD test also has additional benefits over the KS test. First of all, the distribution tails are susceptible to differences. Secondly, the AD test is more effective, even in large sample sizes, in detecting minimum differences (Engmann & Cousineau, 2011). Therefore, the top two probability distributions were selected for the whole study area, according to the Anderson-Darling(AD) goodness-of-fit test.

Kolmogorov-Smirnov (KS) test
The Kolmogorov-Smirnov (KS) framework was created in the early 1930s by Kolmogorov (1933) and continued by Smirnov (1939) as a method for detecting differences in hypothetical distributions. The KS statistic for a cumulative distribution H(x) for a given data set of x data is H(x) will be rejected if KS n is larger than the critical value KS α at (alpha = 0.01, 0.05, and 0.10). If the observed distribution exceeds the theoretical expectations, the null hypothesis will be rejected. For more detail, see this article (Massey Jr, 1951).

Anderson-Darling (AD) test
This test was developed by Anderson & Darling (1952) as a new statistical test for detecting sample distribution departures from normality. Assuming that H is the observed distribution and Hn is the function of the observed (sample) function, then the quadratic function of the differences between H and Hn is given by the empirical function (Stephens, 1974). Quadratic EDF statistics is.
Which is obtained when the weight function is w(x) = [H (x) (1 -H (x))] −1 , with which AD distance places more weight on observations in the tails of the distribution. The AD test is defined as the null hypothesis (H o ): The data follow a specific distribution versus the Alternative statement H a : The data do not follow the specific distribution. We can assess whether the observed data sample x i (i = 1, 2, 3, n) comes from some specified probability distribution. The test indicates that with a given underlying distribution, and assuming that the data is drawn from this distribution is coherent with it or not.

The Modified Mann-Kendall test
The Mann-Kendall (Kendall, 1975;Mann, 1945) test is the most commonly utilized nonparametric tool for assessing dynamics in the data analysis of meteorological time series. The null hypothesis is in the MK test; the knowledge is independent and ordered at random. However, autocorrelation (positive) within the data seeks to increase the probability of trend detection when there is no trend in reality and conversely. After all, it is generally a well-known phenomenon that some researchers have addressed this issue by ignoring the autocorrelation consequence. A modified form of the MK test (Hamed & Rao, 1998) has been used in the current study. The benefit of using the modified MK test is that the evidence of autocorrelation is robust. Adjusted variance (Var(S)) is used in this method to determine Z statistics from the standard Mann-Kendall test (Hamed & Rao, 1998). The mathematical form of the MK test is first determined using the S mk statistics given below.
In the above formulas, n determines the amount of the data, and according to the test principle, the data values over time jth and kth, respectively, must be greater than 10, x i and x j , and in sgn x j − x i , is the sgn function; as Var (S mmk ), the variance value for modified MK, W reflects modified coefficient of auto-correlated data, r k present the kth auto correlated coefficient, x Indicates mean of the series, P represents the number of bound groups, and ti represents the number of degrees one ties in the Mann-Kendall test evaluation rate; Similarly, equation 3's positive value stipulates a growing drift, whereas a negative sign suggests a declining trend in the sequence. The standard normal distribution and the expected Z cal value are compared according to the confidence limits (α =5% or α =1%) (Mann, 1945;Kendall, 1975;Helsel & Hirsch, 2002;Kampata, Parida & Moalafhi, 2008). If the calculated Z cal value is more than |Z MK | < |Z (1−α/2) |, and the null hypothesis (H o ) is accepted.

Trend slope
The Sen slope is used to investigate the trend line, and its magnitude used in this research, based on research done by Theil (1950) and Sen (1968) calculated as described in the following equation: Where, 1<s<t<n, Q β Illustrates the trend line estimator, and x t , Reflects the t th data observed. The positive value of Q β Demonstrates the upward trend direction and its magnitude, and a negative expresses a declining pattern (Yue, Pilon & Phinney, 2003).

Performance analysis for statistical probability distributions
In recent decades, reference evapotranspiration, drought, precipitation, and aridity index and its concern probability theory have seen renewed popularity (Blain, 2012;Blain & Meschiatti, 2014;Huang et al., 2017;Jiang et al., 2020;Nam & Choi, 2014;Rahman et al., 2013;Wu, Xiong & Wu, 2017;Zhang et al., 2015). Hence, it is vital to understand the role of regional ET probability distributions in both dry and wet conditions. This could provide a scientific foundation for the management and allocation of regional water resources. Thus, the analysis was carried out for the mean monthly ET o values by FAO-56(PM), and Hargreaves-Samani (HS) estimates for twelve locations. A total of 120 months of the test data set has use for distribution adjustment. Tables 1 and 2 show the ET o estimates' descriptive statistics for the entire region's study period. Further, the Tables S2 and S3 show the best fit data adjustment to each weather station's probabilistic models based on the AD test (p < 0.05). According to the AD test, we have checked twenty-one probability distribution to tested stations and applied on a series of ET o based on FAO-56(PM) and Hargreaves-Samani estimates. The top three fitted distributions out of twenty-one tested probability distributions were selected for every included station in the first task. For this purpose, we used the AD test in our analysis to determine FAO-56(PM) based ETo Probability distribution models, Hargreaves-Samani (HS) based ET o Probability distribution models, and Seasonal established ET o Probability distribution models for every location in the selected region. Table 2, results revealed that the Johnson SB distribution performed the best by occupying 50% of the total stations. However, the Generalized Pareto showed a 25% inclusion in the study region. In Contrast, we applied the same type of procedure for evapotranspiration estimates based on HS methodology and, the results have shown in the Table S2. The Generalized Pareto showed a 33% presence among all the stations. Besides, Generalized Gamma(4P) and Generalized Extreme Value contributed 16% in overall stations. Moreover, the probability distributions scheme has also been applied in distinct seasons such as Winter, Spring, Dry Summer, Monsoon, and Autumn. The results designate that the Cauchy distribution occupied 25%, Burr and Logistic 16% respectively among all the stations in the winter season. However, in the Spring season, the Burr contributed 33%, while Johnson SB and Generalized Extreme Value showed 16% of the total stations. During Pre-Monsoon, the Johnson SB and log Pearson-3 showed their presence at 41.6% and 25%, respectively. Apart from this, Johnson SB, Burr, and Cauchy performed 33%, 25%, and 16% in overall stations during the red summer. While in Post monsoon, the Burr showed 41.6%, Johnson SB, and Burr(4P) with 16% results among all the study regions. In the second stage, the same procedure was checked by the Kolmogorov-Smirnov(KS) test statistic with (p < 0.05) and the top three best-fitted distribution selected for every station. According to the KS test in the Tables S4 and S5, the results demonstrate that the Johnson SB distribution performed the best by occupying 58.33% of the whole study regions. Likewise, Generalized Pareto showed a 25% enclosure in the study area on top. Hence, Generalized Pareto with 50% performance makes it the second best-fit distribution in the entire study region. In comparison, the HS ET o estimates indicate that the Johnson SB occupies 91.6% presence among all the stations. The Generalized Pareto contributed 75% in overall stations holding the second position. However, based on seasonal ET o estimates, Cauchy and Burr contributed 41.6% and 25% respectively during winter, while Generalized Extreme Value showed 41.6% and Gamma with 16% out of the overall study region in spring. In pre-monsoon, Generalized Extreme Value and Gumbel-Min showed their presence at 25% and 16%, respectively. While in red summer Generalized Extreme Value with 25% and Gumbel-Max showed its presence at 16%. Furthermore, the Burr and Generalized Extreme value distribution performed 41.6% and 16%, respectively, during post-monsoon.
In the present research analysis, according to the evidence provided by our goodness of fit tests, the top two Statistical probability distribution which offered the best adjustment to the more significant part of the study region based on mean monthly ET o values were the Generalized Pareto and Johnson SB (4 parameters lognormal distribution). However, the Burr, Johnson SB, and GEV showed the best fit and great adjustment during the covered period. Therefore, the Johnson SB and Generalized Pareto took top place for the entire Khyber Pakhtunkhwa region. On the other hand, Burr, Johnson SB, and GEV have shown stunning results in different seasons in the whole sample space. For this reason, the scale parameters lambda and alpha (λ&α) and shape parameters Gamma and Beta (γ & β) for both the probability distribution function were obtained and evaluated by the Anderson Darling test at (p < 0.05). The parameter values are shown in Tables 3 and 4, for every weather station. In the meantime, It is exciting to note that the scale parameter lambda for Johnson SB contains the minimum and maximum values (7.62 & 12) with no significant variation except Risalpur. The scale parameter value approach (63.6), the

Reference evapotranspiration trends analysis
Pakistan has ranked 18th among the most susceptible countries in the Global Climate Risk Index in 2011. The state is on the third highly critical threat in 180 nations, impacted by weather losses (Kreft, Eckstein & Melchior, 2013). Consequently, it is essential to consider ET o trends and their implications for territorial dry and wet circumstances. This could provide an objective validity for regional water resource management and allocation, as well as for scientific judgment linked with flood and drought disaster prevention (Ma et al., 2017). This research aims to conduct a temporal trend analysis of reference evapotranspiration in Khyber Pakhtunkhwa using the FAO-56 Penman-Monteith equation to assess their trends over the study period. Therefore, we applied the Mann-Kendall trend test and the Sen slope estimator to find more information about the regional reference evapotranspiration trends. The MK and Sen's slope estimator (Q β ) are determined using the mean monthly evapotranspiration results of ten years. Table 5 summarizes the findings for 12 stations in this respect. This list demonstrates two-tailed confidence patterns (5%) for the overall data series in all the locations. According to the Mann Kendall test at a 5% level of significance, the entire study region showed insignificant upward and downward mean monthly evapotranspiration trends in the concerned period. To accurately identify the trend's magnitude and direction, we used the Sen Slope estimator. Based on the Sen slope estimator, five study locations, including Drosh, Kakul, Kohat, Parachinar, and Risalpur, showed an increasing trend in the overall testing period.   In Contrast, seven study regions, such as Balakot, Cherat, Chitral, DI Khan, Dir, Peshawar, and Saidu Sharif, showed a decreasing trend. In conclusion, the Seasonal distribution of ET o trends in the study area is illustrated graphically in Fig. 4. The black triangle symbol shows the statistically significant increasing or decreasing trends, and the green triangle indicates the insignificant trends. The overall monthly, seasonal, and annual trends are discussed one by one in the following sub-sections.

Monthly evapotranspiration trend analysis
The evaluation of trend was carried out each month at every station; the results showed in Tables 6, 7, 8 and 9. According to the output, in January, 8 out of 12 (66.6%) locations showed insignificant negative trends, while Kohat, DI Khan, Risalpur, and Peshawar indicated positive drifts with an overall 33.3%. Similarly, in February, overall, 9 (75%) stations revealed a negative trend, and three (25%) showed a positive direction. Further, March's month showed 8 (66.6%) out of 12 locations showed insignificant negative trends than Cherat. Drosh, Parachinar, and Peshawar indicated positive behavior of movement, while in April total of 10 (83.3%) study regions out of 12 showed an overall negative trend. Likewise, both showed an overall (50%, 33.3%) negative and (50%, 66.6%) positive trend

Seasonal and annual evapotranspiration trends analysis
The seasons in Pakistan include winter (December to February), Spring (March to April), Pre-monsoon (Dry Summer) period (May to June), Monsoon period (July to September), and Post-Monsoon (autumn) from (October to November). The seasonal evapotranspiration trends in the Northwest region of Pakistan have checked using the same framework (MK and Q β ) shown in Tables 10 and 11. The statistically significant and insignificant trends are graphically presented in Fig. 4, respectively. The reference evapotranspiration showed negative patterns over 91.6 percent of stations overall, based on Mann Kendall and Sen slope seasonal analysis during the winter season. Only 8.33 percent of stations showed a positive trend. Balakot, Parachinar, and Saidu Sharif showed statistically significant negative trends results with Sen slope values (−0.068 mm, −0.334 mm, and −0.067 mm) per season at a 5% significance level. Overall an insignificant negative trend in evapotranspiration was observed in 75% of stations, while 25% showed positive behavior during the Spring season. Further, Balakot showed a negative statistically significant trend with Sen Slope (−0.457 mm) in the Dry Summer, while 66% of stations revealed statistically insignificant trends in the same time season. In Monsoon, Drosh showed a positive statistical significant trend  with Sen slope value (0.944 mm), while 58.3% of stations showed insignificant positive trends during the study period. However, almost 83.33% of stations revealed positive statistically insignificant trends; the rest showed negative autumn trends. Comparatively, the Annual evapotranspiration trends showed negative behavior over 75% of stations. Balakot showed a statistically significant negative trend with Sen slope value (−0.508 mm) over an annual period. In contrast, a statistically significant positive trend was present in the Drosh region with a Sen slope value (0.95 mm). Finally, a climatic analysis of temperature and reference evapotranspiration for the whole study region is shown in Fig. 5.

SUMMARY AND CONCLUSION
This paper has provided a search of the best fit distribution for the reference evapotranspiration and a thorough analysis of twelve meteorological stations of Khyber Pakhtunkhwa. We applied the twenty-one statistical probability distribution at each station and selected the most appropriate according to KS and AD fit tests approach. We also checked the ET trends through Mann-Kendall and Sen slope trend analysis. In general, we found that Johnson SB and Generalized Pareto probability distributions' performance is a better fit and adjusted to the mean monthly ETo. In seasonal ETo estimates, the Burr, Johnson SB, and GEV showed the best results.
Besides, compared to the monthly time series, the seasonal and annual time series appear to have a more significant trend variation. At 91.66% of stations, the winter season tended to decrease, while mixtures of rising and falling trends has observed in spring and Figure 5 Climograph of temperature and reference evapotranspiration. The blue color represents minimum temperature, the red color represents maximum temperature, while the bars stand for reference evapotranspiration. It is important to note that this is the minimu (A-L).
Full-size DOI: 10.7717/peerj.11597/ fig-5 autumn. Also, 83.3 percent of stations in the Post Monsoon season indicated rising drifts. The analysis found a reducing trend pattern identified on average in winter, spring, and annual cycles. The preliminary findings demonstrated that 11 out of 12 stations on average exhibited a downward trend in ET, with three statistically significant tests observed in winter and Monsoon at Balakot, Drosh, Parachinar, and Saidusharif. Analogously, the high rate of positive patterns in the warm seasons could lead to good crops and greenery in the whole area with sufficient rain. Thus, the high return of the flow distribution would be appropriate. The province's agriculture and economy rely predominantly on rainfall as a significant water source; any deficit, reduction, and decrease in ET o or decrease in precipitation may contribute to water supplies scarcity, desertification, intermittent drought, and damaged savannas. It also ensures that agriculture will need more natural water consumption, which is already constrained. Furthermore, as a recommendation, a reasonable adaptive capacity such as irrigation maintenance, soil water ecological restoration, reliable agricultural production, change in plant species, sufficient crop intensities would be required to cope with the potential impacts. Further research must commence for assessing the proposed scheme to other locations in the country. In summing up, these research findings might be helpful in water resource management, water engineering, and policymakers for future planning.