Temporal and Spatial Characteristics of Precipitation and Temperature in Punjab, Pakistan

Identifying the changes in precipitation and temperature at a regional scale is of great importance for the quantification of climate change. This research investigates the changes in precipitation and surface air temperature indices in the seven irrigation zones of Punjab Province during the last 50 years; this province is a very important region in Pakistan in terms of agriculture and irrigated farming. The reliability of the data was examined using double mass curve and autocorrelation analysis. The magnitude and significance of the precipitation and temperature were visualized by various statistical methods. The stations’ trends were spatially distributed to better understand climatic variability across the elevation gradient of the study region. The results showed a significant warming trend in annual Tmin (minimum temperature) and Tmean (mean temperature) in different irrigation zones. However, Tmax (maximum temperature) had insignificant variations except in the high elevation Thal zone. Moreover, the rate of Tmin increased faster than that of Tmax, resulting in a reduction in the diurnal temperature range (DTR). On a seasonal scale, warming was more pronounced during spring, followed by that in winter and autumn. However, the summer season exhibited insignificant negative trends in most of the zones and gauges, except in the higher-altitude Thal zone. Overall, Bahawalpur and Faisalabad are the zones most vulnerable to warming annually and in the spring, respectively. Furthermore, the elevation-dependent trend (EDT) indicated larger increments in Tmax for higher-elevation (above 500 m a.s.l.) stations, compared to the lower-elevation ones, on both annual and seasonal scales. In contrast, the Tmin showed opposite trends at higher- and lower-elevation stations, while a moderate increase was witnessed in Tmean trends from lower to higher altitude over the study region. An increasing trend in DTR was observed at higher elevation, while a decreasing trend was noticed at the lower-elevation stations. The analysis of precipitation data indicated wide variability over the entire region during the study period. Most previous studies reported no change or a decreasing trend in precipitation in this region. Conversely, our findings indicated the cumulative increase in annual and autumn precipitation amounts at zonal and regional level. However, EDT analysis identified the decrease in precipitation amounts at higher elevation (above 1000 m a.s.l.) and increase at the lower-elevation stations. Overall, our findings revealed unprecedented evidence of regional climate change from the perspectives of seasonal warming and variations in precipitation and temperature extremes (Tmax and Tmin) particularly at higher-elevation sites, resulting in a variability of the DTR, which could have a significant influence on water resources and on the phenology of vegetation and crops at zonal and station level in Punjab.


Introduction
Global metrics of future climate changes are directly associated with regional climatic variations.There is mounting evidence that the climate is changing in response to human activity, and the variations in temperature and precipitation recorded during the twentieth century cannot be entirely attributed to natural causes [1].It has been established that the global mean surface air temperature (hereafter mentioned as T mean ) has increased by approximately 0.72 • C during 1951-2012, with a slower increase during 1998-2012 [2].Similarly, the amount of global land precipitation has increased by 2% during the last century [3], and changes in its pattern would directly influence the water resources and agriculture of the concerned regions [4].However, the variations in temperature and precipitation are dynamic in nature and exhibit continuous change over time and space.There are many measures of climate change, and changes in annual and seasonal temperature and precipitation are significant indicators that enable the assessment of climate variability [5,6].A recent study reported that the increasing tendency of T mean could be attributed to the noteworthy variations in maximum and minimum temperatures (hereafter mentioned as T max and T min ), which could ultimately lead to surface warming and changes in climate and precipitation patterns [7][8][9].In addition to variations in the above-mentioned indicators, the diurnal temperature range (DTR) (i.e., the difference between T max and T min ) is an important meteorological indicator and index of observed climate change [10].Across the globe, there have been different variations in daily T max and T min , resulting in a smaller DTR and an increase in T mean [11].The amplification in variation of a) temperature indices and b) precipitation amounts is also highly dependent on the regional elevation gradient.Increasing evidence suggests that rate of warming trends in mountainous regions is faster as compared with the average global warming trends [12].Nonetheless, the variation in climate with respect to elevation does not follow a general trend but it, rather, exhibits different patterns from region to region.
It has been widely accepted that long-term changes in climate variables can alter ecological and hydrological processes [13], and have significant effects on the growth and development of vegetation that lead to a decline in the net agricultural productivity [14].The deleterious effects of climate change would mostly be endured by under-developed nations.Moreover, the natives of developing countries are mostly engaged in agriculture, which is a sector that is highly susceptible to climate change, and have limited resources to adapt to changing circumstances.In this context, the detection of historical variations in climate indicators is highly important for the countries where agriculture is the backbone of the economy, such as Pakistan, which is situated in one of the zones with a rapidly increasing temperature [15].In recent studies, significant climatic changes have been documented in Pakistan, indicating that the T max and T min fluctuated at rates of 0.12-0.29 and 0.10-0.37• C/decade, respectively [16,17].Previous reports have indicated variations in temperature extremes in different parts of the country; for example, positive trends were observed in T max and T min in the upper, middle, and lower Indus basin [18], while decreasing (increasing) trends in T max (T min ) were observed in the northern part of the country [19].Similarly, studies on trend assessments of precipitation have indicated increasing tendencies in annual precipitation for the northern, northeastern and northwestern regions [20,21].Conversely, [22] reported an increase in temperature and a decrease in precipitation and corresponding river flows, starting from the late 1990s, for different sub-basins in Upper Indus Basin (UIB).They estimated the trends by using six atmospheric reanalysis products due to the lack of availability of in-situ observation at UIB.Furthermore, central and southern parts of the country are mostly reported to have experienced a decreasing trend in annual precipitation [23,24].However, no comprehensive research has been conducted regarding the spatiotemporal trends of the temperature indices (e.g., T max , T min , T mean , DTR) and precipitation in various zones of Punjab, which are highly significant from the perspectives of agriculture and irrigated farming, because the major agricultural Water 2019, 11,1916 3 of 23 commodities of the country are produced in these zones, and studies have shown that these zones are highly susceptible to climate change [25].While most previous studies have not focused much on the important segments of climatic variations at the zonal and sub-basin levels, this study intends to bridge the knowledge gap through the use of more detailed and comprehensive research on the spatiotemporal variations in air temperature indices and precipitation in different irrigation zones of Punjab Province, Pakistan.The primary goal of this work is to assess the spatial and temporal trends of the air temperature indices (e.g., T max , T min , T mean , DTR) and precipitation in different irrigation zones of Punjab.The secondary objective is to analyze the trends and heterogeneity in the air temperature indices (e.g., T max , T min , T mean , DTR) and precipitation across the elevation gradient of the whole Punjab region.In addition, the spatial distribution of the annual and seasonal station trends, especially an analysis during the crop growing season (Rabi and Kharif) in different irrigation zones, is presented.
The next sections of the paper are organized as follows.Section 2 provides a brief description of the study area.The data and methods section includes the different statistical methods and techniques, while the results and discussion section presents the outcomes of the study.The conclusions of the study are presented in the last section.

Study Area
Pakistan is a sovereign country in southwest Asia situated between 24 • -37 • N and 60 • -75 • E. The country has an area of 8 × 10 6 km 2 .The country has a diverse landscape ranging from the Karakoram and Himalayan mountains in the north and northwest, with the agricultural plains of the Indus River basin in the center and the Arabian Sea along the southern coast [26].Punjab Province is Pakistan's second largest province, with geographical coordinates of 31.17 • N and 72.70 • E, as shown in Figure 1.Additionally, the province has the largest population and is the agricultural hub of the country, producing more than 50% of the country's agricultural commodities [27].The province accommodates five major rivers, namely, Jhelum, Chenab, Ravi, Bias, and Sutlej.There are seven irrigation zones in Punjab, including Thal, Sargodha, Lahore, Multan, Faisalabad, Dera Ghazi (D.G.) Khan, and Bahawalpur.These zones are distributed on the basis of flow canals and distributaries and are managed by the provincial irrigation department of Punjab.
Water 2018, 10, x FOR PEER REVIEW 3 of 23 most previous studies have not focused much on the important segments of climatic variations at the zonal and sub-basin levels, this study intends to bridge the knowledge gap through the use of more detailed and comprehensive research on the spatiotemporal variations in air temperature indices and precipitation in different irrigation zones of Punjab Province, Pakistan.The primary goal of this work is to assess the spatial and temporal trends of the air temperature indices (e.g., Tmax, Tmin, Tmean, DTR) and precipitation in different irrigation zones of Punjab.The secondary objective is to analyze the trends and heterogeneity in the air temperature indices (e.g., Tmax, Tmin, Tmean, DTR) and precipitation across the elevation gradient of the whole Punjab region.In addition, the spatial distribution of the annual and seasonal station trends, especially an analysis during the crop growing season (Rabi and Kharif) in different irrigation zones, is presented.
The next sections of the paper are organized as follows.Section 2 provides a brief description of the study area.The data and methods section includes the different statistical methods and techniques, while the results and discussion section presents the outcomes of the study.The conclusions of the study are presented in the last section.

Study Area
Pakistan is a sovereign country in southwest Asia situated between 24°-37° N and 60°-75° E. The country has an area of 8 × 10 6 km 2 .The country has a diverse landscape ranging from the Karakoram and Himalayan mountains in the north and northwest, with the agricultural plains of the Indus River basin in the center and the Arabian Sea along the southern coast [26].Punjab Province is Pakistan's second largest province, with geographical coordinates of 31.17°N and 72.70° E, as shown in Figure 1.Additionally, the province has the largest population and is the agricultural hub of the country, producing more than 50% of the country's agricultural commodities [27].The province accommodates five major rivers, namely, Jhelum, Chenab, Ravi, Bias, and Sutlej.There are seven irrigation zones in Punjab, including Thal, Sargodha, Lahore, Multan, Faisalabad, Dera Ghazi (D.G.) Khan, and Bahawalpur.These zones are distributed on the basis of flow canals and distributaries and are managed by the provincial irrigation department of Punjab.The annual mean temperature varies from 23 to 26 • C, with a T min of 16-19 • C and a T max of 29-33 • C. The annual mean precipitation ranges from >800 mm in the northern part to <300 mm in the southern part [28].There are two precipitation seasons in this region (i.e., the monsoon season (July-September) and winter (December-March)) [29].Overall, the northern part of the region is dominated by humid and sub-humid climates, while the central and southern parts are dominated by tropical and coastal climates, respectively.

Data
Long-term variations in the annual and seasonal air temperature indices (e.g., T max , T min , T mean , DTR) and precipitation were investigated in the present study.The historical records from 1967 to 2017, regarding 16 meteorological stations situated in different irrigation zones of Punjab, were provided by the Regional Meteorological Department of Pakistan are shown in Table 1.These stations were selected based on their historical temporal coverage, homogeneity, and completeness of records.The climatic variables were analyzed on an annual and a seasonal scale.The required seasons were defined as winter (December-February), spring (March-May), summer (June-August), and autumn (September-November) [30], and the crop growing periods, Rabi (November-April) and Kharif (May-June).The seasonal and annual climatic data records were obtained using the means of the monthly averages.

Data Processing
The quality of the recorded climatic data series, such as the identification of outliers and quest for missing values, is of primary importance.Therefore, to ensure the quality of the recorded data series, outliers were fixed using the adjacent gauge records, and missing values were obtained using those from nearby stations [31].The gaps in the data series were completed using the time-based interpolation method (monthly records were determined as the mean values of the same month for a period between ± two years [32].
Many approaches have been developed to identify the outliers and non-homogeneities in the records [33].In the present study, the double-mass curve technique was used for the climatic data series for each station [34].The detection or adjustment of discrepancies in the gauge records was performed by associating its period trends with those of other comparatively stable chronicles of a station [1].Non-linearity or bends can serve as a sign of the gauges being relocated or new instruments being installed [35].If the outcome of the double-mass curve analysis showed a straight line, there were no obvious break points, confirming the high temporal consistency in the data series.
To attain the characteristics of air temperature and precipitation in different irrigation zones, Punjab was divided into seven irrigation zones according to the canal command areas, which include Lahore, Multan, Sargodha, Bahawalpur, Faisalabad, Thal, and D.G. Khan (Figure 1).The duration of the analyzed period was 51 years (i.e., from 1967 to 2017).The series of annual and seasonal climatic variables (e.g., T max , T min , T mean , and precipitation) for different irrigation zones were obtained using the Thiessen polygon method, as shown in Figure 1 and Table S1 in the Supplementary Information (SI).This method was developed by [36], and it has been widely used to estimate the areal averages from meteorological and hydrological station data [37,38].According to the position of the stations, polygons are formed by the perpendicular bisectors of the lines joining nearby stations.Thus, each polygon comprises only one station, and the weights of the stations are computed by their relative areas, which are estimated using a polygon network [39].

Methods
Before applying the Mann-Whitney (MW), Mann-Kendall (MK), and Sen's slope techniques to detect the absolute change and trends in the air temperature indices and precipitation data series, the data series were tested using a time serial autocorrelation technique.In addition, the autocorrelation was removed from the data series using the pre-whitening method.Furthermore, the absolute changes in the temperature indices and trends were analyzed using the Mann-Whitney, Mann-Kendall, and Sen's slope techniques to evaluate the patterns and their significance levels in the data series of climatic variables in different irrigation zones.Moreover, inverse distance weighting (IDW) was incorporated as a spatial interpolation technique into the station data to analyze the spatial distribution of the temperature indices and precipitation in different irrigation zones.

Time Serial Autocorrelation
The MK test requires time series data to be serially independent.However, this condition was not present in the time series data.The probability of a significant trend would be added to the data series by the increased autocorrelation and would affect the results of the MK test [40].The presence of significant autocorrelation should be checked and removed before applying the MK trend test [41].Therefore, the following procedure was adopted before applying the MK trend test [42].Compute the lag-1 serial correlation coefficient (r1) of the different time series data as follows: where r1, X n , and X are the correlation coefficient at lag-1, climatic variable data series, and mean of the climatic variable time series, respectively.The confidence interval (CI) of r1 at the 5% significance level can be calculated as follows: If r1 falls into the confidence limits (lag-1 ≥ 0.254 or lag-1 ≤ −0.254), the data series is independent, and the MK test can be applied to the original data series.Otherwise, the MK test should be applied after the removal of the significant autocorrelation, which can be computed as The analyses of autocorrelation coefficients at lag-1 for the annual data series of T max , T min , T mean , DTR, and precipitation in different zones are described in Section 1, Figure S1 in the SI.

Mann-Whitney Test (MW)
The MW test is a non-parametric test used to evaluate the significance of a difference in the median or central trend of two data series.MW analysis was performed to determine the difference in means from the data series [43].Sample data do not require a normal distribution for the MW test.The governing equation of the MW test for each dataset is as follows: Here, n 1 and n 2 are the number of observations, and R 1 and R 2 are the ranks assigned in the first and second periods, respectively.For the application of the MW test, the data series were divided into two spans (i.e., 1967-1991 and 1992-2017) to obtain the absolute change between the two periods.

Mann-Kendall Test (MK)
The non-parametric trend diagnosis MK test was used in the present study.The MK test has been widely used for the detection of trends in metrological or hydrological temporal data series [44,45].The method is less sensitive to the abrupt break points and is robust against outliers and missing values [46].However, the test is sensitive to the serial correlation that may affect the test results [45,47].During the current study, the sequential correlation approach was employed prior to applying the MK test to examine the statistical significance of trends in the climatic data series.
The detailed procedure of the MK test is reported by [48,49].The standardized MK test statistic "S" is defined below: where n is the length of the series, and x i and x j are the sequential values of data.If n ≥ 10, then the S statistic is normally distributed and the variance (Var) is given by where m is the number of tied groups, and t i is the size of the ith tied group.The test statistic Z is computed by: The Z statistics in the MK test follow the normal distribution with an average of zero and a variance of one with a null hypothesis of no trend [50].The positive and negative values of the statistical Z shows increasing and decreasing trends, respectively.

Theil and Sen's Slope (TSS)
TSS is a non-parametric test used to estimate the slope magnitude in linear trends [51], and the test is frequently used for the analysis of hydro-meteorological data series [52].TSS is a least square regression and is widely used for estimating the magnitude of linear trends in large data sets [53].The basic governing equation of Sen's slope for the N pair of data is computed by: Water 2019, 11, 1916 7 of 23 where x j and x k are the data values at time j and k, respectively.If N is an odd number, then Sen's equation is: If N is an even number, then Sens's equation is: 2.4.5.Spatial Interpolation Spatial interpolation techniques can be used to visualize data variability in space and time.The technique predicts values that are identical to the measured values from the sample [54].IDW is one of the most commonly used deterministic methods for interpolating the climate variability distribution.This method has been used by many researchers for the spatial interpolation of data in different regions [55].The method estimates the characteristic values at unsampled sites using a linear arrangement of values at sampled points, and these values are weighted by a converse function of the distance from the point of attention to the sampled location.The weights can be expressed as: where d i is the distance between x o and x i , p is a power constant, and n represents the number of sample points used for the estimation.The weights decrease as distance increases, mainly when the rate of the power parameter increases; thus, adjacent samples have a heavier weight and have more influence on the estimation [56].

Absolute Changes and Trends of Temperature Indices
Figure 2 provides the graphical demonstration of annual time series and the linear trends for climatic variables for each irrigation zone as well as for of the whole region of Punjab.The results of linear trends and coefficient of determination of T min and T mean data series indicate a clear increasing trend.Meanwhile, the trends and coefficient of T max data series depicted non-significant trends and lower r-squared values in most of the zones except for the Thal zone.On the other hand, the DTR data series showed a higher decreasing trend with higher r-squared values, which can be attributed to the higher rate of the increasing trend in T min than that of T max .The highest value of r-squared in T mean and T min data series was observed for Bahawalpur and Multan zone, respectively.Moreover, the highest r-squared value in T max and DTR data series was noticed for Thal and Sargodha zone, respectively.Overall, the linear increase in the T min , T mean data series was more conspicuous than that of the T max in most of the zones and for of the overall Punjab region, during the whole study period.The most significant rapid change in temperature's linear trends was observed during the period 1997-2002 for most of the zones and for the whole Punjab region.According to [19], the country faced a severe and long lasting drought during this period.The rapid increase in temperature trends during the aforementioned period in Punjab region could be owing to the upshots of severe drought conditions in major parts of the country.
The annual and seasonal absolute change and trends of T max , T min , T mean , and DTR for all of Punjab Province obtained by the MW, MK, and Theil-Sen approaches are presented in Tables 2 and 3.Meanwhile, the annual and seasonal absolute changes and trends of T max , T min , T mean , and DTR for the irrigation zones other than the entire Punjab basin are presented in Tables S2 and S3  The spatial distribution of the annual and seasonal trends of T max is presented in Figure 3.The MK trends of the annual T max showed conflicting signals of non-significant positive and negative trends over the whole basin.In particular, the Murree and Islamabad stations, which are located at high altitudes, exhibited a significant increasing trend in T max with magnitudes of 0.55 and 0.24 • C/decade, respectively.In contrast, the highest reduction in T max was observed at Mianwali (−0.16 • C/decade), followed by Sarghoda (−0.13 • C/decade) station.In winter, T max indicated a slight declining trend at most stations except for the higher-altitude Muree and Islamabad stations.In contrast, the spring season indicated an insignificant increasing trend of Tmax at most sites.However, a significant increasing trend in Tmax was observed for the Murree, Islamabad, and Lahore stations, with rates of 0.61, 0.28, and 0.33 °C/decade, respectively.In summer, most of the gauge sites showed an insignificant declining trend except for the Muree, Islamabad, and Bahawalnagar stations, which indicated significant increasing trends.The decrease in Tmax during the summer season is in agreement with the findings of [57], who reported a decreasing trend in some parts of the country.The highest declining trend in Tmax was observed at the Mianwali station with a magnitude of −0.42 °C/decade.Meanwhile, the autumn season revealed a decreasing trend at most sites, with the highest reduction observed at the Mianwali (−0.32 °C/decade) station.In contrast, the spring season indicated an insignificant increasing trend of T max at most sites.However, a significant increasing trend in T max was observed for the Murree, Islamabad, and Lahore stations, with rates of 0.61, 0.28, and 0.33 • C/decade, respectively.In summer, most of the gauge sites showed an insignificant declining trend except for the Muree, Islamabad, and Bahawalnagar stations, which indicated significant increasing trends.The decrease in T max during the summer season is in agreement with the findings of [57], who reported a decreasing trend in some parts of the country.The highest declining trend in T max was observed at the Mianwali station with a magnitude of −0.42 • C/decade.Meanwhile, the autumn season revealed a decreasing trend at most sites, with the highest reduction observed at the Mianwali (−0.32 • C/decade) station.
On the other hand, the seasonal (crop growing) T max showed a slight increasing trend at most stations in the Rabi season, while the Kharif season revealed a minor decrease in T max .A significant negative trend was observed at the elevated gauge sites of Muree and Islamabad for the Rabi season, while two gauge sites, Muree and Bahawalnagar, showed negative trends for the Kharif season.Overall, the annual and seasonal T max trends varied from −0.45 to 0.75 • C/decade.The increasing tendency of T max at high altitudes are consistent with the previously reported work, for the northern belt of Pakistan by [58].Additionally, the decreasing trend in T max observed in this work is in conformity with the findings of [59].Furthermore, our zonal and spatial results are in good agreement with the findings of [17], where it is reported that the trend of T max at the annual scale increased at a rate of 0.12 • C/decade in Pakistan (country level).The present study highlights a significant increase in T max in the Thal zone at a high elevation compared to that at lower elevations.The increasing trend of T max in the higher elevation zone could have a major impact on water resources and agriculture from different perspectives.Nevertheless, this change could have positive implications, such as shortening the crop growing period, as suggested by [60], and it may also lead to an accelerated snow melting process and substantial changes in the amount and pattern of precipitation [18].
Water 2018, 10, x FOR PEER REVIEW 9 of 23 On the other hand, the seasonal (crop growing) Tmax showed a slight increasing trend at most stations in the Rabi season, while the Kharif season revealed a minor decrease in Tmax.A significant negative trend was observed at the elevated gauge sites of Muree and Islamabad for the Rabi season, while two gauge sites, Muree and Bahawalnagar, showed negative trends for the Kharif season.Overall, the annual and seasonal Tmax trends varied from −0.45 to 0.75 °C/decade.The increasing tendency of Tmax at high altitudes are consistent with the previously reported work, for the northern belt of Pakistan by [58].Additionally, the decreasing trend in Tmax observed in this work is in conformity with the findings of [59].Furthermore, our zonal and spatial results are in good agreement with the findings of [17], where it is reported that the trend of Tmax at the annual scale increased at a rate of 0.12 °C/decade in Pakistan (country level).The present study highlights a significant increase in Tmax in the Thal zone at a high elevation compared to that at lower elevations.The increasing trend of Tmax in the higher elevation zone could have a major impact on water resources and agriculture from different perspectives.Nevertheless, this change could have positive implications, such as shortening the crop growing period, as suggested by [60], and it may also lead to an accelerated snow melting process and substantial changes in the amount and pattern of precipitation [18].Due to significant changes in precipitation, disproportionate snow melting can cause rivers to overflow, which may consequently cause flooding in the low-elevation zones [61].Nonetheless, the seasonal negative variations in Tmax could be a possible explanation for the westerlies and monsoon periods in the different zones [62].In the case of Tmin, the results indicated the absolute positive change and increasing trends at the annual and seasonal time scales.The maximum positive absolute change and increasing trend were observed in the Bahawalpur zone during the annual and spring seasons with magnitudes of (1.01 and 1.26 °C absolute changes) and (0.35 and 0.48 °C/decade trends), respectively.Due to significant changes in precipitation, disproportionate snow melting can cause rivers to overflow, which may consequently cause flooding in the low-elevation zones [61].Nonetheless, the seasonal negative variations in T max could be a possible explanation for the westerlies and monsoon periods in the different zones [62].In the case of T min , the results indicated the absolute positive change and increasing trends at the annual and seasonal time scales.The maximum positive absolute change and increasing trend were observed in the Bahawalpur zone during the annual and spring seasons with magnitudes of (1.01 and 1.26 • C absolute changes) and (0.35 and 0.48 • C/decade trends), respectively.
Overall, the spring season showed the highest increasing trend for the nighttime temperature in the southern zone.The increasing trend in T min at the annual and seasonal scales strengthened the results of previous studies for different areas of Pakistan [17].However, the non-significant magnitude of the increasing trend in the summer season conflicted with the results of [59], where a sharp increase in trends during the summer season was reported in different zones of the country.The outcomes of the present study are in harmony with the findings of [18] for the summer season.The spatial distribution of the annual and seasonal trends in T min is presented in Figure 4.The spatial MK trends of the annual T min revealed a clear tendency of warming over the entire Punjab basin except at the Barkhan station, which exhibited a significant decreasing trend with a magnitude of −0.12 • C/decade.The highest positive trend in T min was observed for the Sargodha station with a magnitude of 0.46 • C/decade.For the winter season, the T min revealed a clear increasing tendency, except for the insignificant decreasing trend in the T min at the Murree and Barkhan stations.The warming in the winter season was also reported by [63] at the country level.Similarly, in the spring season, T min showed a significant increasing behavior at most stations located in central and southern Punjab, except for Barkhan, which had an insignificant decreasing trend.Our results are consistent with the finding of [64].The highest increase in T min was observed at the Bahawalnagar station with a magnitude of 0.58 • C/decade.This value was the highest magnitude observed at the annual and seasonal scales.Meanwhile, the trend analysis of the climatic variables in the summer season revealed a conflicting signal of negative and positive slopes at most stations in Punjab.The inclinations in T min indicated an insignificant increasing trend for most stations except D.I. Khan and Barkhan, where significant decreases were observed.The maximum increase in T min was observed at the Khanpur (0.24 • C/decade) station.The variations in the autumn season indicated a clear tendency of a conspicuous increase in T min.The maximum increasing trend in T min was observed at the Khanpur station with a magnitude of 0.59 • C/decade.The results were in accordance with the findings of [16,17] for autumn T min over the whole region.
Furthermore, the seasonal (crop growing) trends in T min revealed a very clear tendency of warming at most stations, except for the northwestern stations, which showed slight insignificant decreasing trends.The highest trend magnitude in T min was observed at the Sarghoda (0.51 • C/decade) and Khanpur (0.42 • C/decade) stations in the Rabi and Kharif seasons, respectively.Overall, the annual and seasonal T min trends varied from −0.25 to 0.60 • C/decade.The present results indicated elevated trends in T min with a faster rate than that of T max in the different irrigation zones and stations, except for the high-elevation sites, and this difference could have a serious impact on the evapotranspiration rate and consequently on crop productivity [57].In addition to direct effects on plant physiological processes, increases in T min could affect plant growth by changing the duration of the growing season or by changing the soil water availability [65].It was also suggested that nighttime warming aggravates the adverse effects of soil water stress, as reported in [66], and could reduce the productivity of the rice crop by 10% for each 1 • C increase in the growing season [67].Overall, there is a need for comprehensive research on the interaction between elevated T min and the ecosystem for the region under study.The interactions between climatic variability and different ecosystem characteristics are not considered in the present work because that is beyond the scope of the study.Overall, the spring season showed the highest increasing trend for the nighttime temperature in the southern zone.The increasing trend in Tmin at the annual and seasonal scales strengthened the results of previous studies for different areas of Pakistan [17].However, the non-significant magnitude of the increasing trend in the summer season conflicted with the results of [59], where a sharp increase in trends during the summer season was reported in different zones of the country.The outcomes of the present study are in harmony with the findings of [18] for the summer season.The spatial distribution of the annual and seasonal trends in Tmin is presented in Figure 4.The spatial MK trends of the annual Tmin revealed a clear tendency of warming over the entire Punjab basin except at the Barkhan station, which exhibited a significant decreasing trend with a magnitude of −0.12 °C/decade.The highest positive trend in Tmin was observed for the Sargodha station with a magnitude of 0.46 °C/decade.For the winter season, the Tmin revealed a clear increasing tendency, except for the insignificant decreasing trend in the Tmin at the Murree and Barkhan stations.The analysis of T mean indicated a positive absolute change and warming trends at the annual and seasonal scales, with the exception of the summer season, which showed a decreasing trend.On an annual scale, the maximum absolute change and trend was observed in the Bahawalpur zone with magnitudes of 0.64 and 0.20 • C/decade, respectively.The present results of the annual trends also confirm the findings of [68], who reported a slope of annual T mean of approximately 0.24 • C/decade in Pakistan during 1960-2007.On a seasonal scale, the maximum change and trend were observed in the Faisalabad zone during the spring season with magnitudes of 1 and 0.34 • C/decade, respectively.Meanwhile, the summer season showed a negative trend in the entire Punjab basin.The decreasing trend of T mean in the summer season was previously reported for other areas of Pakistan [63].The spatial distribution of the annual and seasonal trends of T mean is presented in Figure 5.The MK results for the annual T mean showed a warming trend over the entire Punjab basin except for the Barkhan station, which exhibited a significant decreasing trend with a magnitude of −0.05 • C/decade.The highest positive trend in the annual T mean was observed for the Muree station with a magnitude of 0.27 • C/decade.These observations contradict the findings of [60], who reported negative trends in temperature at high-altitude regions in the country.
In winter, the T mean revealed a clear increasing tendency at most stations except for Muree and Barkhan.Similarly, the spring T mean showed a significant increasing behavior at most stations, with the highest increase observed at Bahawalnagar with a magnitude of 0.44 • C/decade; however, an insignificant decreasing trend was observed at Barkhan.In contrast, the summer season revealed a conflicting signal of negative and positive slopes at most stations in Punjab, with the highest positive trend observed at Muree and Islamabad at a rate of 0.16 • C/decade.The variations in the autumn season indicated a clear tendency for a conspicuous increase in the T min.The maximum increasing trend in the T mean was observed at the Khanpur station with a magnitude of 0.32 • C/decade.In winter, the Tmean revealed a clear increasing tendency at most stations except for Muree and Barkhan.Similarly, the spring Tmean showed a significant increasing behavior at most stations, with the highest increase observed at Bahawalnagar with a magnitude of 0.44 °C/decade; however, an insignificant decreasing trend was observed at Barkhan.In contrast, the summer season revealed a conflicting signal of negative and positive slopes at most stations in Punjab, with the highest positive trend observed at Muree and Islamabad at a rate of 0.16 °C/decade.The variations in the autumn season indicated a clear tendency for a conspicuous increase in the Tmin.The maximum increasing trend in the Tmean was observed at the Khanpur station with a magnitude of 0.32 °C/decade.
Similarly, trends in Tmean during the crop growing seasons revealed a clear tendency of warming at most stations except for those in the northwest, where slight insignificant decreasing trends were observed.The highest trend magnitude in Tmean was observed at the Bahawalnagar station in the Bahawalpur zone with magnitudes of 0.26 and 0.23 °C/decade in the Rabi and Kharif seasons, respectively.The overall annual and seasonal Tmean trends varied from −0.25 to 0.45 °C/decade.The increasing trend of Tmean in different zones and stations has been related to various aspects, such as the concentration of anthropogenic activities, aerosol effects on the climate, increased cloud cover, land-use change, and industrial growth [69].It is also suggested that there is a significant effect of urbanization and land-use change on surface warming at the local scale [70].The present findings report a higher increasing trend of Tmean in the Faisalabad zone, which is a highly populated and industrial hub of the country; these findings are associated with the results of previous studies related to warming caused by urbanization and industrial development.[27] reported similar findings in different industrial and populated cities of Pakistan.In addition, various agricultural practices, such as increased irrigation and different cropping patterns, may have potential impacts on the Tmean trends in the Indian sub-continent [71].
Furthermore, the magnitude of the absolute change and trends of the DTR demonstrated a significant decreasing shift at the annual and seasonal scales for the different zones except those located at higher elevations (e.g., Thal and D.G Khan showed a non-significant increasing trend in the winter, spring, and Rabi seasons).On an annual scale, the highest negative absolute change and trend was observed in the Faisalabad and Bahawalpur zones, with magnitudes of −0.79 and −0.36 °C/decade, respectively.Similarly, on a seasonal scale, the highest declining change and trend was observed in the Sarghoda zone during the autumn season, with magnitudes of −1.45 and −0.57Similarly, trends in T mean during the crop growing seasons revealed a clear tendency of warming at most stations except for those in the northwest, where slight insignificant decreasing trends were observed.The highest trend magnitude in T mean was observed at the Bahawalnagar station in the Bahawalpur zone with magnitudes of 0.26 and 0.23 • C/decade in the Rabi and Kharif seasons, respectively.The overall annual and seasonal T mean trends varied from −0.25 to 0.45 • C/decade.The increasing trend of T mean in different zones and stations has been related to various aspects, such as the concentration of anthropogenic activities, aerosol effects on the climate, increased cloud cover, land-use change, and industrial growth [69].It is also suggested that there is a significant effect of urbanization and land-use change on surface warming at the local scale [70].The present findings report a higher increasing trend of T mean in the Faisalabad zone, which is a highly populated and industrial hub of the country; these findings are associated with the results of previous studies related to warming caused by urbanization and industrial development.[27] reported similar findings in different industrial and populated cities of Pakistan.In addition, various agricultural practices, such as increased irrigation and different cropping patterns, may have potential impacts on the T mean trends in the Indian sub-continent [71].
Furthermore, the magnitude of the absolute change and trends of the DTR demonstrated a significant decreasing shift at the annual and seasonal scales for the different zones except those located at higher elevations (e.g., Thal and D.G Khan showed a non-significant increasing trend in the winter, spring, and Rabi seasons).On an annual scale, the highest negative absolute change and trend was observed in the Faisalabad and Bahawalpur zones, with magnitudes of −0.79 and −0.36 • C/decade, respectively.Similarly, on a seasonal scale, the highest declining change and trend was observed in the Sarghoda zone during the autumn season, with magnitudes of −1.45 and −0.57• C/decade, respectively.Overall, there was a clear decreasing tendency in the DTR data series over the entire study region.
The spatial distribution of the annual and seasonal trends of the DTR is presented in Figure 6.The spatial MK trends of the annual DTR indicated a declining behavior at most stations except Muree and Barkhan, where a significant increase in the DTR was observed.The maximum increasing and decreasing trends were observed at the Murree (0.55 • C/decade) and Sarghoda (−0.58 • C/decade) stations, respectively.In the winter season, most stations showed significant declining behavior, with the highest fluctuation observed at Muree and Sarghoda, at rates of 0.82 and −0.The spatial distribution of the annual and seasonal trends of the DTR is presented in Figure 6.The spatial MK trends of the annual DTR indicated a declining behavior at most stations except Muree and Barkhan, where a significant increase in the DTR was observed.The maximum increasing and decreasing trends were observed at the Murree (0.55 °C/decade) and Sarghoda (−0.58 °C/decade) stations, respectively.In the winter season, most stations showed significant declining behavior, with the highest fluctuation observed at Muree and Sarghoda, at rates of 0.82 and −0.71 °C/decade, respectively.Similarly, in the spring season, most of the gauge sites showed negative trends except for the high-altitude Murree and Barkhan stations.In contrast, the summer DTR values were characterized by positive and negative trends over the entire Punjab basin, with the highest increase observed at Murree (0.41 °C/decade) and D.G. Khan (0.22 °C/decade), while the greatest reduction was observed at Khanpur station (−0.52 °C/decade).
Moreover, trends in the DTR during the crop growing seasons exhibited a negative change at most stations except those in the northern and northwestern sites, which showed significant increasing trends, particularly at Murree, D.G. Khan, and Barkhan.The overall annual and seasonal DTR trend varied from −0.85 to 0.85 °C/decade in the Rabi and Kharif seasons, respectively.The magnitude of the decrease in the DTR during Rabi was greater than that of the Kharif season.The highest increase and decrease were observed at Murree (0.64 °C/decade) and Sarghoda (−0.69 °C/decade) for the Rabi season, and at Murree (0.43 °C/decade) and Khanpur (−0.56 °C/decade) for the Kharif season.These results are also supported by previous findings in the region [27].Several past reports suggest that there is a strong relationship between the cloud cover and DTR.More cloud cover is associated with a decrease in the DTR [72].Similarly, the outcomes of the present study, which indicated a decreasing trend in the DTR for the autumn, summer, and Kharif seasons, indicate a negative relationship with cloud cover and the occurrence of significant precipitation during these seasons.

Annual and Seasonal Trends of Precipitation
The annual time series and linear trends of precipitation for each irrigation zone and for the entire province of Punjab are presented in Figure 7.The annual average precipitation showed a slight increasing trend in all zones, including the whole basin, during the study period 1967-2017.The  [27].Several past reports suggest that there is a strong relationship between the cloud cover and DTR.More cloud cover is associated with a decrease in the DTR [72].Similarly, the outcomes of the present study, which indicated a decreasing trend in the DTR for the autumn, summer, and Kharif seasons, indicate a negative relationship with cloud cover and the occurrence of significant precipitation during these seasons.

Annual and Seasonal Trends of Precipitation
The annual time series and linear trends of precipitation for each irrigation zone and for the entire province of Punjab are presented in Figure 7.The annual average precipitation showed a slight increasing trend in all zones, including the whole basin, during the study period 1967-2017.The corresponding regression line along with r-squared values were also plotted for each irrigation zone and for all of Punjab.The southern zones, including Bahawalpur, Multan, and D.G. Khan, showed an annual precipitation frequency of approximately 100-300 mm.However, the central zone (Faisalabad, Lahore, and Sarghoda) and the northern zone (Thal) indicated precipitation frequencies of approximately 400-700 and ≥ 1000 mm, respectively.Similarly, the entire Punjab basin showed an annual precipitation of approximately 400-500 mm.The highest and lowest r-squared values were noticed in Bahawalpur and Lahore zones, respectively.The highest fluctuation in linear trend was observed in Thal followed by the Lahore zone.In addition, the annual and seasonal precipitation trends and magnitude for the entire Punjab basin obtained by the MK test and the Theil-Sen approach are presented in Table 3.Meanwhile, the annual and seasonal precipitation trends for irrigation zones other than the entire Punjab basin are presented in Table S3 in the SI.Finally, the spatial distributions of the annual and seasonal precipitation trend of the stations are presented in Figure 8.
annual precipitation of approximately 400-500 mm.The highest and lowest r-squared values were noticed in Bahawalpur and Lahore zones, respectively.The highest fluctuation in linear trend was observed in Thal followed by the Lahore zone.In addition, the annual and seasonal precipitation trends and magnitude for the entire Punjab basin obtained by the MK test and the Theil-Sen approach are presented in Table 3.Meanwhile, the annual and seasonal precipitation trends for irrigation zones other than the entire Punjab basin are presented in Table S3 in the SI.Finally, the spatial distributions of the annual and seasonal precipitation trend of the stations are presented in Figure 8.The results of the MK test indicated the wide variability of precipitation for the annual and seasonal time scales in the majority of the zones during the period 1967-2017.The annual precipitation revealed a significant increasing trend in all irrigation zones, with a magnitude ranging from 15.01 to 34.45 mm/decade.The maximum cumulative increase in annual precipitation was observed in the Thal zone (34.83 mm/decade), followed by Sarghoda (25.55 mm/decade) and Lahore (23.45 mm/decade).Similarly, the spatial distribution of the annual station trends showed that most sites were characterized by non-significant positive trends.The results of the MK test indicated the wide variability of precipitation for the annual and seasonal time scales in the majority of the zones during the period 1967-2017.The annual precipitation revealed a significant increasing trend in all irrigation zones, with a magnitude ranging from 15.01 to 34.45 mm/decade.The maximum cumulative increase in annual precipitation was observed in the Thal zone (34.83 mm/decade), followed by Sarghoda (25.55 mm/decade) and Lahore (23.45 mm/decade).Similarly, the spatial distribution of the annual station trends showed that most sites were characterized by non-significant positive trends.
However, considerable positive tendencies were observed at low elevation for the Bahawalnagar, Bahawalpur, Khanpur, and Mianwali stations at rates of 16.4 21.57, 27.09, and 57.73 mm/decade, respectively.In contrast, a decreasing trend was detected for the high altitude Murre station, with a value of −45 mm/decade, in the Thal zone.The increasing trend of annual precipitation in southern Punjab contradicted the results of [73], who reported a decreasing trend in annual precipitation in this region.In winter precipitation, the negative trends were more obvious over the whole region, particularly in the D.G. Khan and Bahawalpur zones, including the entire Punjab basin.The negative trend in the winter season correlates with the findings of [74], who reported a decreasing trend of precipitation in different areas of Pakistan.The spatial distribution of the MK results showed that most stations revealed insignificant signs of positive and negative slopes during the winter season.Most of the winter precipitation was due to western disturbances, and consequently, most stations in the study area indicated less precipitation.The maximum decline in winter precipitation was observed at the high-altitude Murree station (−28.5 mm/decade) in the Thal zone.
whole region, particularly in the D.G. Khan and Bahawalpur zones, including the entire Punjab basin.The negative trend in the winter season correlates with the findings of [74], who reported a decreasing trend of precipitation in different areas of Pakistan.The spatial distribution of the MK results showed that most stations revealed insignificant signs of positive and negative slopes during the winter season.Most of the winter precipitation was due to western disturbances, and consequently, most stations in the study area indicated less precipitation.The maximum decline in winter precipitation was observed at the high-altitude Murree station (−28.5 mm/decade) in the Thal zone.The zonal trend analysis of spring precipitation showed insignificant negative and positive increases over the whole basin.However, prominent increasing trends were observed in the southern zone, while significant decreasing trends were observed in the northern zones of Punjab.The spatial distribution of the MK trends indicated variable fluctuations for most stations.However, a significant increase was observed for the Bahawalnagar station, at the rate of 4.41 mm/decade.In summer precipitation, the MK trend results illustrated the insignificant positive trends in different irrigation zones and all of Punjab, except for the Bahawalpur zone, which showed a significant increase during summer.The highest increase was assessed in the Lahore and Sargodha zones, with a rate of 11.76 mm/decade.Similarly, the spatial distribution of the station trends indicated an insignificant increasing trend over the entire study area.According to the MK analysis results, 13 stations showed positive trends, while three stations revealed negative tendencies.A significant increasing trend was detected at the Mianwali and Khanpur stations, with magnitudes of 35.68 and 12.14 mm/decade, respectively.
Meanwhile, autumn precipitation revealed a significant increasing trend in all zones, with a magnitude ranging from 0.78 to 9.81 mm/decade, with the highest increase assessed in the Thal zone at a rate of 9.81 mm/decade.The spatial distribution of station trends indicated that 14 out of 16 stations showed a positive trend; Barkhan and Sialkot, which exhibited negative trends, were the two The zonal trend analysis of spring precipitation showed insignificant negative and positive increases over the whole basin.However, prominent increasing trends were observed in the southern zone, while significant decreasing trends were observed in the northern zones of Punjab.The spatial distribution of the MK trends indicated variable fluctuations for most stations.However, a significant increase was observed for the Bahawalnagar station, at the rate of 4.41 mm/decade.In summer precipitation, the MK trend results illustrated the insignificant positive trends in different irrigation zones and all of Punjab, except for the Bahawalpur zone, which showed a significant increase during summer.The highest increase was assessed in the Lahore and Sargodha zones, with a rate of 11.76 mm/decade.Similarly, the spatial distribution of the station trends indicated an insignificant increasing trend over the entire study area.According to the MK analysis results, 13 stations showed positive trends, while three stations revealed negative tendencies.A significant increasing trend was detected at the Mianwali and Khanpur stations, with magnitudes of 35.68 and 12.14 mm/decade, respectively.
Meanwhile, autumn precipitation revealed a significant increasing trend in all zones, with a magnitude ranging from 0.78 to 9.81 mm/decade, with the highest increase assessed in the Thal zone at a rate of 9.81 mm/decade.The spatial distribution of station trends indicated that 14 out of 16 stations showed a positive trend; Barkhan and Sialkot, which exhibited negative trends, were the two outliers.The highest increasing trend was observed at the Islamabad station (17.26 mm/decade) in the Thal zone.
In addition, the precipitation trends during the crop growing seasons (Rabi and Kharif) indicated non-significant increasing trends in most zones except for the Bahawalpur zone, which showed a significant increase during the Rabi and Kharif seasons, while the entire Punjab basin showed a significant increase during the Kharif season.The highest increase was observed during the crop growing period at Sarghoda in the Rabi and Thal zones in the Kharif season, with rates of 5.22 and 30.92 mm/decade, respectively.The overall zonal trend magnitudes ranged from 1.23 to 5.22 and 11.68 to 30.92 mm/decade during the Rabi and Kharif seasons, respectively.Moreover, the spatial distribution of precipitation trends during the crop growing seasons (Rabi and Kharif) indicated considerable fluctuations at most stations.Furthermore, the station precipitation indicated that the magnitude of the positive trend in Kharif was more pronounced than in the Rabi season.The overall annual and seasonal precipitation trends magnitude ranged from −50 to 60 mm/decade during the Rabi and Kharif seasons, respectively.The findings of the current investigation for summer season contradict the results of [64], who reported a considerable increase in precipitation during the summer season.
Moreover, the entire Punjab basin showed a significant increasing trend of precipitation for the annual and autumn seasons, with magnitudes of 22.24 and 5.23 mm/decade, respectively.Similar results were reported by [21] in different regions in Pakistan.In addition, the Bahawalpur zone indicated a slight significant increasing trend in precipitation at all time scales except for the winter season.Altogether, the study area exhibited a significant positive shift in annual and autumn precipitation compared with that in the other seasons.

Annual and Seasonal Elevation-Dependent Trends (EDTs)
The analysis for annual EDT in temperature indices (T max , T min , T mean , DTR) and precipitation over the Punjab region is shown in Figure 9 (black dash linear fit).The geographical locations of the available stations in the region of Punjab located between the elevation range of 55-2167 (m a.s.l.) are shown in Table 1.The EDT result of T max indicates a positive gradient from lower to higher elevation with significant positive trends observed at higher elevation (above 500 m a.s.l.), with the magnitude of 0.55 and 0.24 • C/decade at Murree and Islamabad stations (p =< 0.01).On the contrary, the EDT of T min showed a negative slope from lower to higher elevation with an overall significant positive trend for most of the low elevation stations, except for the higher-altitude Barkhan and Murree stations.Furthermore, the EDT in T mean showed a moderate increasing slope from lower to higher altitude with significant positive trends for most of the stations.However, the EDT in DTR showed a steep slope from lower to higher altitude with a significant increase at higher elevation stations above 500 m a.s.l.On the other hand, the EDT analysis for precipitation indicates a negative gradient from lower to higher altitude, with overall positive trend at most of the stations at lower elevation and significant negative or no trend for the high-altitude stations, particularly Murree and Barkhan (above 1000 m a.s.l.).Overall, the results indicate an increase in the amount of precipitation at low-elevation stations, whereas, a decrease in amount of precipitation was noticed for higher-elevation stations.A decrease in precipitation for higher elevation (above 4000 m a.s.l) was also reported by [22] for different sub-basin levels in Upper Indus Basin by assembling different reanalysis products, during the period 1997-2014.As a consequence, decline in river flows was reported for different sub-basins in the UIB.Similarly, the results of EDTs for temperature indices and precipitation, obtained in this study, concur well with the findings of [75] who reported a similar increasing trend in T max and DTR, as well as a decrease in T min and precipitation trends, with a moderate positive slope in T mean from lower to higher gradient across the longitudinal extension of Nepal.However, for the region of Punjab, most stations are located at less elevation, yet the results obtained from the EDT analysis are highly significant, particularly for the futuristic studies related to the elevation-dependence warming over the region under study.
The EDT analysis presented in Figure 9 (black dash linear fit) deals with the 16 stations located between the elevation range of 55-2167 m a.s.l.However, most of the stations in this range are located below 500 m a.s.l., in exception to three stations (Murree, Islambad, and Barkhan) that are located above the said elevation.Generally, a chance of uncertainty could exist if such a smaller number of stations (above 500 m a.s.l) are utilized to perform EDT analysis,.Under such scenario, the uncertainty of the gauge data can be minimized by assembling different atmospheric reanalysis products.For instance, [22,76] quantified the uncertainty of climatic variables by using six atmospheric reanalysis products for cross validation of river flow data for poorly gauged regions.Meanwhile, another way to deal with such scenarios, where there are lesser number of stations involved, is to extend the analysis by incorporating data from a greater number of stations (if available) and compare the outcomes with the results obtained for analysis performed for fewer number of stations.Therefore, in order to address any ambiguities related to the results of our analysis for elevation-dependent temperature indices (T max , T min , T mean , DTR) and precipitation trends, we performed the additional EDT analysis by including a greater number of higher elevation meteorological stations (Table 4).The stations that cover the continuous elevation range above 500-2167 m a.s.l. were incorporated in the EDT extension analysis, as can be seen through Figure 9 (red linear fit).The outcome of EDT analysis, which included the additional stations, indicated a similar linear trends in the precipitation and the temperature indices, as was observed for the stations under the scope of this study (the trend is presented as black dash linear fit).Moreover, the additional stations in EDT extension analysis confirmed an analogous trend pattern, and data quality with overall higher range of magnitude and better statistical significance along the extension gradient of the Punjab and surrounding region at higher elevation.The spatial and seasonal distributions of EDT of temperature indices (Tmax, Tmin, Tmean, DTR) and precipitation over the Punjab region are shown in Figures 3, 4, 5, 6, and 8. Significant positive trends in Tmax were noticed at higher-elevation (above 500 m a.s.l.) sites in Thal and D.G. Khan zones, while low-elevation sites did not show any trend or negative slopes at the seasonal scale.The highest positive trends were observed in the winter season at higher-elevation Murree, Islamabad, and Barkhan sites, with the magnitude of 0.74, 0.47, and 0.48 °C/decade, respectively.The trend in Tmax is highly sensitive with respect to the elevation and can be easily detected by spatial distribution of trends.The current findings of Tmax trends are in consistence with the results of [75], wherein higher trends in temperature in Nepal at elevations above 1000 m a.s.l. were reported and were said to be homogenous for all seasons.Similarly, [77] reported the increasing trends of Tmax at higher elevation in Nepal during the 1990s, without finding any further significant dependence of these trends on the elevation levels.
Further, in our study, EDT in Tmin revealed positive slopes for most of the sites at low elevation, while negative or no trends in Tmin were observed for all seasons at high-elevation (above 500 m a.s.l) sites.In accordance with Tmin trends, the seasonal and spatial EDT of Tmean indicates the positive significant slopes in all seasons except summer.However, the higher-altitude Islamabad and Murree stations in the Thal zone showed an increasing trend during the summer season.In case of DTR, the higher-altitude (above 500 m a.s.l.) sites (i.e., Murree, Islamabad, and Barkhan stations) showed  The spatial and seasonal distributions of EDT of temperature indices (T max , T min , T mean , DTR) and precipitation over the Punjab region are shown in Figures 3-6 and 8. Significant positive trends in T max were noticed at higher-elevation (above 500 m a.s.l.) sites in Thal and D.G. Khan zones, while low-elevation sites did not show any trend or negative slopes at the seasonal scale.The highest positive trends were observed in the winter season at higher-elevation Murree, Islamabad, and Barkhan sites, with the magnitude of 0.74, 0.47, and 0.48 • C/decade, respectively.The trend in T max is highly sensitive with respect to the elevation and can be easily detected by spatial distribution of trends.The current findings of T max trends are in consistence with the results of [75], wherein higher trends in temperature in Nepal at elevations above 1000 m a.s.l. were reported and were said to be homogenous for all seasons.Similarly, [77] reported the increasing trends of T max at higher elevation in Nepal during the 1990s, without finding any further significant dependence of these trends on the elevation levels.
Further, in our study, EDT in T min revealed positive slopes for most of the sites at low elevation, while negative or no trends in T min were observed for all seasons at high-elevation (above 500 m a.s.l) sites.In accordance with T min trends, the seasonal and spatial EDT of T mean indicates the positive significant slopes in all seasons except summer.However, the higher-altitude Islamabad and Murree stations in the Thal zone showed an increasing trend during the summer season.In case of DTR, the higher-altitude (above 500 m a.s.l.) sites (i.e., Murree, Islamabad, and Barkhan stations) showed positive slopes in all seasons, while the highest trends were observed in the winter season with the magnitude of 0.82, 0.23, and 0.56 • C/decade, respectively.In addition, the low-elevation sites (less than 500 m a.s.l.) indicated negative DTR trends at the seasonal scale.The results of increasing trends in T max and DTR are related to the decrease in the cloud cover, while the precipitation at higher elevation are consistent with the previous findings by [78] for Central Asia and [79] for India, where similar results were reported, owing to a decrease in cloud cover and precipitation, with the increase in temperature trends.The seasonal and spatial EDT of precipitation showed the positive slopes at low elevation sites in all seasons, with an exception for the winter season where negative trends were observed for most of the sites.Moreover, the higher-elevation sites (above 1000 m a.s.l.) showed a negative pattern, which resulted in a negative gradient in precipitation amounts from lower to higher elevation across the Punjab region.The outcomes of this study regarding the increase in amounts of precipitation for lower-elevation regions are consistent with the findings of [80,81], wherein the increase in precipitation amounts and its impact on the increase in T min trends were discussed.The results of the higher decreasing trend of precipitation at annual and seasonal scale at higher elevation, and particularly results of Murree station (above 2000 m a.s.l), are associated with findings of [75,82], which reported a decreasing trend in annual precipitation from lower to higher elevation.The possible explanation of decreasing precipitation at higher elevation could be linked with the decrease in cloud cover and soil moisture and, ultimately, increase in day time and decrease in nighttime temperatures (T max and T min ) [83,84].However, an in-depth investigation is required to explore, in detail, the quantitative relationship between precipitation and temperature indices and to learn about the main driving factors for such asymmetric trends of temperature and precipitation for Punjab, which is a highly worthwhile region with respect to the agriculture and economy of Pakistan.

Conclusions
The present study investigated the spatiotemporal heterogeneity of air temperature (e.g., T max , T min , T mean , DTR) and precipitation on annual and seasonal time scales at seven irrigation zones of Punjab, using 51 years (1967-2017) of data from 16 meteorological stations in the region.The major findings of this paper are as follows: On an annual scale, T min and T mean showed a noteworthy increasing trend, while T max showed conflicting signals of insignificant positive and negative trends, except in the high-elevation Thal zone.Moreover, the rate of T min increased faster than that of T max , resulting in a reduction in the DTR, except for that in the Thal zone.Overall, southern Punjab, especially Bahwalapur, was the most vulnerable zone to warming, while the high-altitude Thal zone was more susceptible to changes in T max .
On a seasonal scale, warming was more pronounced during the spring, followed by winter and autumn seasons.Moreover, on a temporal and spatial scale, the summer season exhibited negative trends at most zones and gauges, with the exception of high-altitude sites.Overall, central Punjab, especially Faisalabad, was the most vulnerable zone to warming during the spring season.
The analysis of precipitation data indicated the wide variability of positive and negative trends during 1967-2017.However, the annual and autumn precipitation showed a significant increase, while the winter precipitation showed a decreasing trend in most zones, including the entire Punjab region.
The spatial and temporal elevation-dependent trend analysis solidified our conclusion regarding the increment in warming (T max ) at higher-elevation (above 500 m a.s.l.) than low-elevation stations, on both annual and seasonal scales.In contrast, the nighttime temperature (T min ) showed opposite trends at higher and lower elevation stations, with moderate increase in T mean slope from lower to higher altitude over the longitudinal extension of the study region.Furthermore, we witnessed an increasing trend in DTR and a decrease in precipitation amounts at higher-elevation stations (above 1000 m a.s.l.) and vice versa at the lower-elevation stations.The warming at higher-elevation stations could be related to the decrease in precipitation and cloud cover and the changing soil moisture conditions during the study period.
Overall, our study progresses scientific knowledge regarding spatial and temporal heterogeneity in climate variables of different irrigation zones of Punjab by using historical station data.The results provide preliminary and significant information of climate variability depending on the station elevations across the longitudinal section of Punjab and surrounding higher-elevation regions.Additionally, this study may prove significant to predict heterogeneity of climate trends for other regions in the world with similar meteorological conditions; this information could enable better agricultural management.Altogether, we found evidence of seasonal variations and warming that could influence the water resources and agricultural sectors of the whole study region.Industrialization, urbanization, and other anthropogenic interventions could be possible reasons for warming at zonal and station levels.On the basis of our analysis, we highly recommend further investigations to establish the cause of warming and to explore the relationship between temperature indices and precipitation at zonal and station levels across the elevation gradient of Punjab.

Figure 1 .
Figure 1.Location of the Punjab irrigation zones, climatic stations, and corresponding Thiessen polygons.
Figure2provides the graphical demonstration of annual time series and the linear trends for climatic variables for each irrigation zone as well as for of the whole region of Punjab.The results of linear trends and coefficient of determination of T min and T mean data series indicate a clear increasing trend.Meanwhile, the trends and coefficient of T max data series depicted non-significant trends and lower r-squared values in most of the zones except for the Thal zone.On the other hand, the DTR data series showed a higher decreasing trend with higher r-squared values, which can be attributed to the higher rate of the increasing trend in T min than that of T max .The highest value of r-squared in T mean and T min data series was observed for Bahawalpur and Multan zone, respectively.Moreover, the highest r-squared value in T max and DTR data series was noticed for Thal and Sargodha zone, respectively.Overall, the linear increase in the T min , T mean data series was more conspicuous than that of the T max in most of the zones and for of the overall Punjab region, during the whole study period.The most significant rapid change in temperature's linear trends was observed during the period 1997-2002 for most of the zones and for the whole Punjab region.According to[19], the country faced a severe and long lasting drought during this period.The rapid increase in temperature trends during the aforementioned period in Punjab region could be owing to the upshots of severe drought conditions in major parts of the country.The annual and seasonal absolute change and trends of T max , T min , T mean , and DTR for all of Punjab Province obtained by the MW, MK, and Theil-Sen approaches are presented in Tables2 and 3.Meanwhile, the annual and seasonal absolute changes and trends of T max , T min , T mean , and DTR for the irrigation zones other than the entire Punjab basin are presented in TablesS2 and S3in the SI.The results (MW and MK) indicated that T max exhibited a positive absolute change and increasing trend

Water 2018 ,
10, x FOR PEER REVIEW 8 of 23°C/decade, respectively.In contrast, the highest reduction in Tmax was observed at Mianwali (−0.16 °C/decade), followed by Sarghoda (−0.13 °C/decade) station.In winter, Tmax indicated a slight declining trend at most stations except for the higher-altitude Muree and Islamabad stations.

Figure 2 .
Figure 2. Time series of annual maximum temperature (Tmax), minimum temperature (Tmin), mean temperature (Tmean) and diurnal temperature range (DTR) in degree celcius (°C) for the Punjab irrigation zones and their linear trends.

Figure 2 .
Figure 2. Time series of annual maximum temperature (T max ), minimum temperature (T min ), mean temperature (T mean ) and diurnal temperature range (DTR) in degree celcius ( • C) for the Punjab irrigation zones and their linear trends.

Figure 3 .
Figure 3. Change per decade and signs of trends at meteorological stations for Tmax for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 3 .
Figure 3. Change per decade and signs of trends at meteorological stations for T max for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 4 .
Figure 4. Change per decade and signs of trends at meteorological stations for Tmin for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 4 .
Figure 4. Change per decade and signs of trends at meteorological stations for T min for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 5 .
Figure 5. Change per decade and signs of trends at meteorological stations for Tmean for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 5 .
Figure 5. Change per decade and signs of trends at meteorological stations for Tmean for the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.
71 • C/decade, respectively.Similarly, in the spring season, most of the gauge sites showed negative trends except for the high-altitude Murree and Barkhan stations.In contrast, the summer DTR values were characterized by positive and negative trends over the entire Punjab basin, with the highest increase observed at Murree (0.41 • C/decade) and D.G. Khan (0.22 • C/decade), while the greatest reduction was observed at Khanpur station (−0.52 • C/decade).Water 2018, 10, x FOR PEER REVIEW 13 of 23°C/decade, respectively.Overall, there was a clear decreasing tendency in the DTR data series over the entire study region.

Figure 6 .
Figure 6.Change per decade and signs of trends at meteorological stations for the DTR of the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 6 .
Figure 6.Change per decade and signs of trends at meteorological stations for the DTR of the year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.Moreover, trends in the DTR during the crop growing seasons exhibited a negative change at most stations except those in the northern and northwestern sites, which showed significant increasing trends, particularly at Murree, D.G. Khan, and Barkhan.The overall annual and seasonal DTR trend varied from −0.85 to 0.85 • C/decade in the Rabi and Kharif seasons, respectively.The magnitude of the decrease in the DTR during Rabi was greater than that of the Kharif season.The highest increase and decrease were observed at Murree (0.64 • C/decade) and Sarghoda (−0.69 • C/decade) for the Rabi season, and at Murree (0.43 • C/decade) and Khanpur (−0.56 • C/decade) for the Kharif season.These results are also supported by previous findings in the region[27].Several past reports suggest that there is a strong relationship between the cloud cover and DTR.More cloud cover is associated with a decrease in the DTR[72].Similarly, the outcomes of the present study, which indicated a decreasing trend in the DTR for the autumn, summer, and Kharif seasons, indicate a negative relationship with cloud cover and the occurrence of significant precipitation during these seasons.

Figure 7 .
Figure 7. Time series and linear trend of annual precipitation in Punjab and the irrigation zones.

Figure 7 .
Figure 7. Time series and linear trend of annual precipitation in Punjab and the irrigation zones.

Figure 8 .
Figure 8. Change per decade and trends of precipitation at meteorological stations (mm) for the entire year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Figure 8 .
Figure 8. Change per decade and trends of precipitation at meteorological stations (mm) for the entire year, winter, spring, summer, autumn, Rabi, and Kharif in Punjab for the period 1967-2017.

Water 2018 , 23 Figure 9 .
Figure 9. Elevation-dependent trend of temperature indices and precipitation for stations under the study region located at elevations above 500 m a.s.l.(black dash linear fit), and for additional stations (red linear fit) covering continuous elevation range of 500-2167 m. a.s.l, during the period 1967-2017.

Figure 9 .
Figure 9. Elevation-dependent trend of temperature indices and precipitation for stations under the study region located at elevations above 500 m a.s.l.(black dash linear fit), and for additional stations (red linear fit) covering continuous elevation range of 500-2167 m. a.s.l, during the period 1967-2017.

Table 1 .
Information of Meteorological stations in Punjab.

Table 2 .
Absolute changes in temperature ( • C) indices in Punjab Province during the second period (1992-2017) compared to the first period (1967-1991) using the Mann-Whitney (MW) test.Bold and italic showed significant change with Mann-Whitney (MW) test at 95% confidence level.

Table 3 .
Annual and seasonal MK temperatures trend ( • C/decade) indices and precipitation (mm/decade) in Punjab Province during 1967-2017 using Mann-Kendall and Sen's slope.
Bold and italic showed a significant trend with Mann-Kendall (MK) test at 95% confidence level.

Table 2 .
Absolute changes in temperature (°C) indices in Punjab Province during the second period (1992-2017) compared to the first period (1967-1991) using the Mann-Whitney (MW) test.Bold and italic showed significant change with Mann-Whitney (MW) test at 95% confidence level.

Table 3 .
Annual and seasonal MK temperatures trend (°C/decade) indices and precipitation (mm/decade) in Punjab Province during 1967-2017 using Mann-Kendall and Sen's slope.

Table 4 .
Information of additional meteorological stations for elevation-dependent trend analysis.