Climate-Induced Perspective Variations in Irrigation Schedules and Design Water Requirements for the Rice–Wheat System

Conceptualizing the implications of climate change for crop evapotranspiration (ETc) and subsequent net irrigation water requirement (NIWR) is critical to sustaining Pakistan’s agriculture and food security. In this article, future ETc, NIWR, and design water requirements (DWR) were projected for the rice–wheat system of Punjab, Pakistan. Consistently increasing temperatures signify an impending hotter transition in the future thermal regime, accompanied by a substantial increase in monsoon rainfall. Future climate warming accelerated ETc and NIWR, which were compensated by 2–5 and 1–2 additional irrigations during the rice and wheat seasons, respectively. Future rice and wheat required 13–18 and 2–5 irrigations per season, respectively. Effective rainfall increments did not compensate for the warming-driven higher ETc and NIWR because of uneven and erratic rainfall distribution. Rainfall occurrence and the duration of peak irrigation demand were mismatched, resulting in surplus rainwater availability during the future rice season. The results suggest that DWR for 5- and 10-year return period droughts during the baseline period (965 and 1000 mm, respectively) should be revised to accommodate the additional 100–200 mm of irrigation water per season; otherwise, the study area will face an acute water shortage in the future.


Introduction
According to the recently released Sixth Assessment Report (AR 6 ) of the Intergovernmental Panel on Climate Change (IPCC), human-induced anthropogenic interventions have unequivocally changed the Earth's climate, which is evident from global surface warming, sporadic rainfall patterns, intensification of the water cycle, and an overall energy imbalance in the atmosphere [1]. Evidence of human-inflicted climate change in triggering the slew of cataclysmic weather and climate extremes and their attribution to human influences has strengthened since AR 5 [1]. The literature includes many studies addressing the assessment and mitigation of climate change impacts across various aspects of human life [2,3]. Agriculture is the primary target of adverse climate change impacts due to its indisputable higher dependence on the climate for survival [4]. Climate-driven natural calamities such as floods, droughts, and warm and dry spells, exacerbated by climate warming or rainfall changes [1], are reshaping crop growth patterns, irrigation demand and supply, and particularly yields [2,3,[5][6][7][8].
For an agriculture-dependent economy like Pakistan, these climate change trends are alarming because of its overreliance on irrigated agriculture and low adaptive capacity [5,6]. Pakistan ranked fifth on the Global Climate Risk Index among the countries frequently experiencing extreme climate events during 1999-2019. A severe drought struck the country in 2018-2019, and 2020 was the wettest recorded year [9]. Ineffective use of scanty water supplies and climate change are key drivers of low agricultural production and food shortages in the country [3,10]. Rice and wheat constitute up to 62% of the

Climate Data
The Pakistan Meteorological Department (PMD) operates the Sialkot weather station, located at latitude 32.51 N and longitude 74.53 E near the study area with reliable longterm data records. Data of the daily maximum temperature (Tmax), minimum temperature (Tmin), number of sunshine hours, solar radiation (Rn), rainfall (P), and monthly averaged relative humidity (RH) and wind speed (u2) were collected from the PMD during the baseline period of 1980-2010. Daily RH and u2 data were also retrieved for a grid containing the Sialkot station from the National Centers for Environment Prediction Climate Forecast System Reanalysis (CFSR) global climate dataset (https://globalweather.tamu.edu/#pubs (accessed on 1 February 2021)) [25][26][27]. Possible biases in the CFSR daily datasets were removed using a linear scale change factor bias correction procedure [17]. The recorded Tmax, Tmin, P, and Rn were combined with the bias-corrected RH and u2 data to devise the daily baseline climatology of the study area [7].
We selected eight GCMs (Supplementary Table S1) from the CMIP5 archives to project the future climate under the IPCC representative concentration pathway (RCPs) scenarios-4.5 and 8.5-during two time slices-the 2030s (2021-2050) and 2060s . The coarse-resolution daily GCM outputs were statistically bias-corrected using quantile mapping, which tailored the distribution functions of the climate variables in raw GCM time series. The cumulative distribution functions (CDFs) of the climate variables in the observed time series were mapped over raw GCM historic time series [28,29]. The same empirical CDFs were applied to the daily time series of the GCM-projected future

Climate Data
The Pakistan Meteorological Department (PMD) operates the Sialkot weather station, located at latitude 32.51 N and longitude 74.53 E near the study area with reliable long-term data records. Data of the daily maximum temperature (T max ), minimum temperature (T min ), number of sunshine hours, solar radiation (R n ), rainfall (P), and monthly averaged relative humidity (RH) and wind speed (u 2 ) were collected from the PMD during the baseline period of 1980-2010. Daily RH and u 2 data were also retrieved for a grid containing the Sialkot station from the National Centers for Environment Prediction Climate Forecast System Reanalysis (CFSR) global climate dataset (https://globalweather.tamu.edu/#pubs (accessed on 1 February 2021)) [25][26][27]. Possible biases in the CFSR daily datasets were removed using a linear scale change factor bias correction procedure [17]. The recorded T max , T min , P, and R n were combined with the bias-corrected RH and u 2 data to devise the daily baseline climatology of the study area [7].
We selected eight GCMs (Supplementary Table S1) from the CMIP5 archives to project the future climate under the IPCC representative concentration pathway (RCPs) scenarios-4.5 and 8.5-during two time slices-the 2030s (2021-2050) and 2060s (2051-2080). The coarse-resolution daily GCM outputs were statistically bias-corrected using quantile mapping, which tailored the distribution functions of the climate variables in raw GCM time series. The cumulative distribution functions (CDFs) of the climate variables in the observed time series were mapped over raw GCM historic time series [28,29]. The same empirical CDFs were applied to the daily time series of the GCM-projected future T max , T min , RH, u 2 , P, and R n . The detailed bias correction outcomes for each selected GCM were provided by Ahmad et al. [7].

Estimation of Crop Water Requirements and Irrigation Scheduling
According to the FAO-56 water balance model, the crop water requirement is equivalent to the water depth lost through crop evapotranspiration (ET c ), which is estimated from the product of reference evapotranspiration (ET o ) and crop coefficients (K c ) [30,31]. We used the FAO-developed CROPWAT model v 8.0 for ET c estimation of the rice-wheat system during the baseline and future time slices. The CROPWAT used the FAO Penman-Monteith equation for daily ET o estimation, and the potential irrigation water requirement (PIWR) was approximated as the difference between the ET c and effective rainfall (ER) at a 10-day time step [30,31]. The 10 daily K c values for the UCC rice-wheat system were referenced from Ullah et al. [32], and the ER was calculated using the United States Department of Agriculture's method [33].
For user-defined soil and irrigation management conditions, CROPWAT generates irrigation schedules based on the daily water balance of the root zone [34]. Seasonal cumulative water demand determines the net irrigation water requirement (NIWR) based on predefined irrigation scheduling criteria and the daily soil water balance. We derived the irrigation schedules and NIWR for the rice-wheat system according to the guidelines from the Directorate of Agriculture Information of the government of Punjab (http://dai. agripunjab.gov.pk/croptechnologies (accessed on 1 March 2021)) as well as the earlier field studies [11,23,35].
During any particular wheat season, three irrigations, each 75 mm in depth, are recommended. In Punjab, wheat is irrigated 1-4 times, depending on the seasonal rainfall and soil type, with an average of 2.38 irrigations applied per season [11]. In this study, the representative soil properties for a medium loam-textured soil were taken from Khaliq et al. [11] (Supplementary Table S2). Wheat irrigation was triggered when the soil moisture in the extractable root zone dropped below 80 mm and each time 75 mm water was applied through flooding. The wheat sowing and harvesting dates were 20 November and 10 April, respectively.
In the UCC command area, a ponding water depth of 50-75 mm is maintained during most of the rice season, and 250-900 mm of puddling water is used for land preparation [36,37]. The soil texture determines the rice field percolation rates, which range from 2 mm/day for heavy clay soils to 6 mm/day for sandy soils [38]. A percolation rate of 1.5 mm/day had been used for rice consumptive water demand estimation in various canal commands of the Indus Basin in Pakistan [32]. CROPWAT defines the after-puddling percolation rate as the maximum percolation rate of non-puddled soil raised to the power of 0.33. We set the puddling criterion as irrigation at a 0-mm water depth and a refill depth of 60 mm, and the scheduling criterion for rice NIWR estimation was irrigation at a 5-mm water depth and a refill depth of 60 mm. The durations of nursery raising and land preparation were 30 days and 25 days, respectively. The rice transplanting and harvesting dates were 1 July and 10 October, respectively.

Trend Analysis
Parametric and non-parametric tests can be used for trend detection in time series. Parametric tests are mostly preferred for trend detection in normally distributed time series, and non-normally distributed time series are treated using non-parametric tests [39]. The rank-based non-parametric Mann-Kendall (MK) test is widely used to detect the direction and statistical significance of trends in hydro-meteorological data. The World Meteorological Organization recommends the MK test because of its applicability to nonnormal data without being influenced by outliers or missing values [40][41][42].
The MK test indicates increasing or decreasing trends as signs of Z-scores, and the presence of a statistically significant trend is detected if |Z| > Z 1−α/2 , where α is the significance level [41]. We used the MK test for trend detection in the seasonal time series of the T max , T min , R n , P, ET o , ET c , PIWR, NIWR, and number of rainy days (N) during the baseline and future time slices. The trend magnitude was estimated using Sen's slope estimator [41,43].

Estimation of the Design Water Requirement
The frequency of droughts in Pakistan has been steadily increasing, as extreme droughts that covered two-thirds of the country were witnessed during 1952, 1969, 1971, 2000, 2001, and 2002 [44]. During the rice season, return periods of moderate, severe, and extreme droughts are in the range of 7-9, 13-19, and 30−45 years, respectively [5]. Since the UCC is a non-perennial canal that supplies irrigation water only during the rice season, the DWR were considered equivalent to the rice NIWR for 5-and 10-year return period droughts.
We conducted a frequency analysis using probability distribution functions (PDFs) from the Generalized Extreme Value (GEV), Gamma (GA), log-gamma (LGA), log-Pearson type 3 (LP3), Normal (NOR), log-Normal (LNOR), Weibull (WEI), Logistic (LOG), and Gumbel (GUM) distributions. The parameters of the distributions were estimated using the maximum likelihood method, and the goodness of fit was evaluated using the Anderson-Darling (AD) and Kolmogorov-Simirnov (KS) tests [45]. Based on the selected PDFs, the rice NIWR for the 5-and 10-year return period droughts were estimated using the Chow frequency factor method [46][47][48][49] and were subsequently considered the DWR during the baseline and future time slices. The MK test indicates increasing or decreasing trends as signs of Z-scores, and the presence of a statistically significant trend is detected if |Z| > Z1−α/2, where α is the significance level [41]. We used the MK test for trend detection in the seasonal time series of the Tmax, Tmin, Rn, P, ETo, ETc, PIWR, NIWR, and number of rainy days (N) during the baseline and future time slices. The trend magnitude was estimated using Sen's slope estimator [41,43].

Estimation of the Design Water Requirement
The frequency of droughts in Pakistan has been steadily increasing, as extreme droughts that covered two-thirds of the country were witnessed during 1952, 1969, 1971, 2000, 2001, and 2002 [44]. During the rice season, return periods of moderate, severe, and extreme droughts are in the range of 7-9, 13-19, and 30−45 years, respectively [5]. Since the UCC is a non-perennial canal that supplies irrigation water only during the rice season, the DWR were considered equivalent to the rice NIWR for 5-and 10-year return period droughts.
We conducted a frequency analysis using probability distribution functions (PDFs) from the Generalized Extreme Value (GEV), Gamma (GA), log-gamma (LGA), log-Pearson type 3 (LP3), Normal (NOR), log-Normal (LNOR), Weibull (WEI), Logistic (LOG), and Gumbel (GUM) distributions. The parameters of the distributions were estimated using the maximum likelihood method, and the goodness of fit was evaluated using the Anderson-Darling (AD) and Kolmogorov-Simirnov (KS) tests [45]. Based on the selected PDFs, the rice NIWR for the 5-and 10-year return period droughts were estimated using the Chow frequency factor method [46][47][48][49] and were subsequently considered the DWR during the baseline and future time slices.    Consistent increments in the projected T max and T min , regardless of the crop season, RCP, or time slice, supported the impending hotter transitions in the study area's future temperature regimes. Noteworthy features of the climate warming magnitude included the following: T min > T max , wheat season > rice season, RCP 8.5 > RCP 4.5, and 2060s > 2030s (Figure 2a,b and Figure 3a,b). Rice season receives monsoon rainfall, which makes up a major proportion of the annual rainfall. The rice season median P featured moderate and intense increments during the 2030s (RCP 4.5) and 2060s (RCP 8.5), respectively ( Figure 3c). The wheat season median P was almost unchanged during the 2030s, but it slightly dropped during the 2060s, particularly under RCP 4.5 ( Figure 2c). Contrary to projected climate warming, the future R n declined notably in the rice season. The R n decline magnitude was as follows: rice season > wheat season and RCP 8.5 > RCP 4.5. For the wheat season, the 2030s R n decline was greater, and for the rice season, the 2060s R n decline was greater compared with their corresponding counterparts (Figures 2d and 3d).

Projected Climate Change
Agronomy 2021, 11, x FOR PEER REVIEW 6 of 18 and hereafter represents the interquartile range (IQR). The whiskers extend up to 1.5 times the IQR, and the hollow dot outside the box represents the outliers.
Consistent increments in the projected Tmax and Tmin, regardless of the crop season, RCP, or time slice, supported the impending hotter transitions in the study area's future temperature regimes. Noteworthy features of the climate warming magnitude included the following: Tmin > Tmax, wheat season > rice season, RCP 8.5 > RCP 4.5, and 2060s > 2030s (Figures 2a,b and 3a,b). Rice season receives monsoon rainfall, which makes up a major proportion of the annual rainfall. The rice season median P featured moderate and intense increments during the 2030s (RCP 4.5) and 2060s (RCP 8.5), respectively ( Figure 3c). The wheat season median P was almost unchanged during the 2030s, but it slightly dropped during the 2060s, particularly under RCP 4.5 ( Figure 2c). Contrary to projected climate warming, the future Rn declined notably in the rice season. The Rn decline magnitude was as follows: rice season > wheat season and RCP 8.5 > RCP 4.5. For the wheat season, the 2030s Rn decline was greater, and for the rice season, the 2060s Rn decline was greater compared with their corresponding counterparts (Figures 2d and 3d).   Table 1 show the MK test Z-scores and Sen's slope for interseasonal trends of the Tmax, Tmin, Rn, and P during the baseline and future time slices. Except for an unchanged Tmax trend during the baseline rice season, the Tmax and Tmin trends significantly (α = 1%) increased regardless of the time slice, RCP, or cropping season. Sen's slopes of the wheat and rice season Tmax and Tmin trends were almost similar for a certain time slice or RCP, with an overall change range of 0.02-0.06 °C/year (Table 1). The future P trends for both crops were mainly non-significant, regardless of their direction. The wheat season P trends were downward throughout the baseline and future time slices, and these were significant only during the baseline (α = 5%) and 2030s under RCP 4.5 (α = 10%). Conversely, the rice season P trends were upward across the baseline and future time slices and only significant during the 2030s under RCP 8.5 (α = 1%). The wheat season Rn trends were highly significant (α = 1%) during the baseline, 2030s (RCP 8.5), and 2060s (RCP 4.5), but the rice seasonal Rn trends did not follow a particular pattern across the time slices of the RCPs. During the 2030s and 2060s, the rice seasonal Rn trends were significant (α = 1%) only under RCPs 8.5 and 4.5, respectively ( Figure 4 and Table 1).   Table 1 show the MK test Z-scores and Sen's slope for interseasonal trends of the T max , T min , R n , and P during the baseline and future time slices. Except for an unchanged T max trend during the baseline rice season, the T max and T min trends significantly (α = 1%) increased regardless of the time slice, RCP, or cropping season. Sen's slopes of the wheat and rice season T max and T min trends were almost similar for a certain time slice or RCP, with an overall change range of 0.02-0.06 • C/year (Table 1). The future P trends for both crops were mainly non-significant, regardless of their direction. The wheat season P trends were downward throughout the baseline and future time slices, and these were significant only during the baseline (α = 5%) and 2030s under RCP 4.5 (α = 10%). Conversely, the rice season P trends were upward across the baseline and future time slices and only significant during the 2030s under RCP 8.5 (α = 1%). The wheat season R n trends were highly significant (α = 1%) during the baseline, 2030s (RCP 8.5), and 2060s (RCP 4.5), but the rice seasonal R n trends did not follow a particular pattern across the time slices of the RCPs. During the 2030s and 2060s, the rice seasonal R n trends were significant (α = 1%) only under RCPs 8.5 and 4.5, respectively ( Figure 4 and Table 1).

Comparison of CROPWAT Results with the Literature
Before projecting future water demands, the CROPWAT simulation outcomes during the baseline were evaluated by comparing the ETc and NIWR of both crops with the published literature. In this study, the baseline rice seasonal cumulative ETc and NIWR varied in the ranges of 580-675 mm and 730-1080 mm, respectively, and for wheat, these were in the ranges of 204-260 mm and 75-300 mm, respectively.
In the Indus Basin of Pakistan, the seasonal cumulative ETc of rice and wheat could vary in the range of 540-1156 mm and 290-520 mm, respectively [13]. According to the field studies conducted across Punjab, the rice and wheat ETc ranges are 520-640 mm and 250-400 mm, respectively (Supplementary Table S3). Shakir et al. [24] reported the average ETc of rice and wheat in the UCC canal command to be 520 mm and 250 mm, respectively. Similarly, our estimated NIWR and number of irrigations for both crops was also supported by the findings of previously conducted surveys or experimental studies [11,15,23,35].

Comparison of CROPWAT Results with the Literature
Before projecting future water demands, the CROPWAT simulation outcomes during the baseline were evaluated by comparing the ET c and NIWR of both crops with the published literature. In this study, the baseline rice seasonal cumulative ET c and NIWR varied in the ranges of 580-675 mm and 730-1080 mm, respectively, and for wheat, these were in the ranges of 204-260 mm and 75-300 mm, respectively.
In the Indus Basin of Pakistan, the seasonal cumulative ET c of rice and wheat could vary in the range of 540-1156 mm and 290-520 mm, respectively [13]. According to the field studies conducted across Punjab, the rice and wheat ET c ranges are 520-640 mm and 250-400 mm, respectively (Supplementary Table S3). Shakir et al. [24] reported the average ET c of rice and wheat in the UCC canal command to be 520 mm and 250 mm, respectively. Similarly, our estimated NIWR and number of irrigations for both crops was also supported by the findings of previously conducted surveys or experimental studies [11,15,23,35].
In Punjab, three irrigations, each 75 mm in depth, are recommended for wheat cultivation, and occasionally, fourth irrigation can also be applied in the case of less seasonal P. In this study, the wheat was irrigated 1-4 times, which was in close agreement with the irrigation scheduling mostly adopted by previous studies to represent the farmers' practices (Supplementary Table S3).
In the case of rice, the recommended ponding depth is 50-75 mm, and 1000-1600 mm water is applied as 15-25 irrigations during a particular season. In this study, the ponding depth was 60 mm, and 11-16 irrigations were applied during the baseline. Khaliq et al. [11] and Gaydon et al. [23] also applied 16 irrigations, each 60 mm in depth, to represent the current rice farming practices in Punjab. A detailed comparison of CROPWAT outcomes in this study with the literature is also provided in Supplementary Table S3 and Figure S1. The comparison suggested that the simulation results fairly represent the current field-scale irrigation practices in Punjab, Pakistan.

Projected Wheat Irrigation Schedules
We projected various parameters that govern the future irrigation water management of wheat, including the ET o , ET c , ER, PIWR, NIWR, N, and number of irrigations. Figure 5 presents the mean and median variations of the above parameters during the baseline and future time slices, and the right-hand side panels show the variations relative to the respective baseline average values. The MK test Z-scores and Sen's slope of the interseasonal trends of the parameters are provided in Figure 6 and In Punjab, three irrigations, each 75 mm in depth, are recommended for wheat cultivation, and occasionally, fourth irrigation can also be applied in the case of less seasonal P. In this study, the wheat was irrigated 1-4 times, which was in close agreement with the irrigation scheduling mostly adopted by previous studies to represent the farmers' practices (Supplementary Table S3).
In the case of rice, the recommended ponding depth is 50-75 mm, and 1000-1600 mm water is applied as 15-25 irrigations during a particular season. In this study, the ponding depth was 60 mm, and 11-16 irrigations were applied during the baseline. Khaliq et al. [11] and Gaydon et al. [23] also applied 16 irrigations, each 60 mm in depth, to represent the current rice farming practices in Punjab. A detailed comparison of CROPWAT outcomes in this study with the literature is also provided in Supplementary Table S3 and Figure S1. The comparison suggested that the simulation results fairly represent the current field-scale irrigation practices in Punjab, Pakistan.

Projected Wheat Irrigation Schedules
We projected various parameters that govern the future irrigation water management of wheat, including the ETo, ETc, ER, PIWR, NIWR, N, and number of irrigations. Figure 5 presents the mean and median variations of the above parameters during the baseline and future time slices, and the right-hand side panels show the variations relative to the respective baseline average values. The MK test Z-scores and Sen's slope of the interseasonal trends of the parameters are provided in Figure 6 and Table 2, respectively. During the baseline, the average wheat seasonal ETo, ETc, ER, PIWR, NIWR, and N were 400 mm, 235 mm, 130 mm, 150 mm, 180 mm, and 16 days, respectively. On average, 2.4 irrigations were applied during the baseline wheat season.    ETo is the primary determinant of ETc rates, and it indicates the overall climate change influence on atmospheric evaporative demand [2,22]. Future climate warming steadily raised the wheat seasonal cumulative ETo and ETc up to 40 mm, with the highest positive shifts during the 2060s under RCP 8.5 relative to the baseline (Figure 5a,b). The difference between ETc and ER determines the PIWR and NIWR for a specific crop. In the 2030s, the future ER and N increased somewhat from the baseline, but in the 2060s, these were nearly identical to the baseline (Figure 5c,g). As a result, the baseline averaged PIWR and NIWR dropped in the 2030s and rose in the 2060s. Despite projected increases in ETo and ETc, positive ER contributions restricted the expected increase in irrigation frequency, with wheat requiring 2.5 and 3.3 irrigations per season in the 2030s and 2060s, respectively. The IQR of the number of irrigations was unclear, since wheat requires a maximum of four irrigations per season (Figure 5f). The results indicate that 2-3 irrigations per season will be the most common irrigation schedule for future wheat, especially in the 2060s.  ET o is the primary determinant of ET c rates, and it indicates the overall climate change influence on atmospheric evaporative demand [2,22]. Future climate warming steadily raised the wheat seasonal cumulative ET o and ET c up to 40 mm, with the highest positive shifts during the 2060s under RCP 8.5 relative to the baseline (Figure 5a,b). The difference between ET c and ER determines the PIWR and NIWR for a specific crop. In the 2030s, the future ER and N increased somewhat from the baseline, but in the 2060s, these were nearly identical to the baseline (Figure 5c,g). As a result, the baseline averaged PIWR and NIWR dropped in the 2030s and rose in the 2060s. Despite projected increases in ET o and ET c , positive ER contributions restricted the expected increase in irrigation frequency, with wheat requiring 2.5 and 3.3 irrigations per season in the 2030s and 2060s, respectively. The IQR of the number of irrigations was unclear, since wheat requires a maximum of four irrigations per season (Figure 5f). The results indicate that 2-3 irrigations per season will be the most common irrigation schedule for future wheat, especially in the 2060s. Figure 6 and Table 2 present the trends of the parameters during a specific time slice. The future ET o and ET c trends were identical to their respective absolute changes relative to the baseline. The ET o increased significantly (α = 5% and 1%) in the range of 0.61-0.79 mm/year, but positive ET c trends were mostly non-significant, except for the 2030s under RCP 4.5 (α = 1%) ( Table 3). The absolute changes in the wheat seasonal cumulative ER and N were positive compared with the baseline, but their trends highlighted inevitable drying of the future wheat seasons. Except for the 2030s under RCP 8.5, the ER and N trends significantly (α = 1%) dropped during the baseline and non-significantly declined during the future time slices. Due to the higher variability range, positive trends were detected in future PIWR, but the NIWR trends were not apparent due to a limited variation range ( Figure 6 and Table 2). Absolute changes occurred in the future PIWR and NIWR compared with the baseline, as well as their respective change trends during the baseline and future time slices behaving in opposite directions. This implies that the warming-induced increase in the ET c will eventually outweigh the positive contribution of the ER, resulting in higher PIWR and NIWR in the future.

Projected Rice Irrigation Schedules
The means, medians, and variability ranges of the rice seasonal cumulative ET o , ET c , ER, PIWR, NIWR, N, and number of irrigations during the baseline, 2030s, and 2060s are shown in Figures 7 and 8, and Table 3 presents the MK test Z-scores and Sen's slope of trends across the time slices and RCPs for the above parameters, respectively. During the baseline, the average rice seasonal ET o , ET c , ER, PIWR, NIWR, and N were 830 mm, 630 mm, 360 mm, 750 mm, 890 mm, and 35 days, respectively. On average, rice was irrigated 13 times per season during the baseline. The baseline ET o and ET c increased up to 40 mm under the influence of projected rice seasonal warming, particularly in the 2060s (RCP 8.5). The baseline ER increased gradually during the 2030s and dramatically increased throughout the 2060s, implying that heavy monsoon rainfalls could improve rainwater availability. Despite overall ER increments of up to 200 mm from the baseline, the future PIWR and NIWR rose to 300 mm by the end of the 2060s (Figure 7d,e). This can be explained by an N decline of up to 10 days, suggesting limited ER contributions to ET c -driven PIWR and NIWR increments because of erratic and unequal distribution of future monsoon rainfall. As a result, we estimated that the future rice production would require 2-5 extra irrigations per season (Figure 7f). The aforementioned calculations showed the absolute change in the respective variables when compared with the baseline average values.
The rice seasonal ET o and ET c trends did not follow a consistent pattern across the time slices, as these were negative in the baseline and 2030s (RCP 8.5) and positive in the remainder of the cases. Nonetheless, these were highly significant in the 2060s for both RCPs (α = 5% and 1%) ( Table 3). Absolute changes in the future ER relative to the baseline and trend directions were identical, showing that future rice seasons will receive more rainfall. The positive ER trends were significant (α = 1% and 5%) under RCP 8.5 and 4.5 during the 2030s and 2060s, respectively. The absolute change in the future N was negative compared with the baseline, but the Z-scores of the MK test were positive (Figure 8), and the trends were not significant, except during the 2030s under RCP 8.5 (α = 1%) ( Table 3). Future PIWR and NIWR trends were downward and mainly non-significant, notably for NIWR, but the PIWR-trend magnitudes were comparably higher than those of NIWR.

Projected Change in Design Water Requirements
The goodness of fit for nine PDFs was evaluated by the KS and AD statistical tests (α = 5%) to determine the best fitting PDF. The estimated parameters of the PDFs based on the maximum likelihood method and goodness-of-fit results during the baseline and future time slices are provided in Supplementary Tables S4-S6. All PDFs were suitable for rice NIWR, but according to both goodness-of-fit tests, the GEV and WEI were most suitable. The GEV is a continuous probability distribution described by three parameters-location, scale, and shape-and these parameters explain the distribution's shift on the horizontal axis, spread, and shape, respectively [50,51]. Among the three, the shape parameter is skewness-dependent, and it represents where most of the data lie. For the baseline and future NIWR time series, we calculated a negative shape factor, which reduced the GEV to the WEI distribution [50,51] (Supplementary  Tables S4-S6). Hence, WEI was chosen as the most suitable PDF for DWR estimation of the 5-and 10-year return periods.    Figure 9 shows the DWR for the 5-and 10-year (DWR 5 and DWR 10 , respectively) return periods during the baseline and future time slices. Despite the expected rise in seasonal cumulative ER, the baseline DWR 5 and DWR 10 (965 and 1000 mm, respectively) showed a future increase of 90 mm. The DWR 5 and DWR 10 thresholds were crossed during 8 baseline years (Figure 9a), and future DWR thresholds were crossed approximately 5 times for a single time slice, except during the 2030s under RCP 8.5 (Figure 9b-e) when the ER increments were limited (Figure 7c). These estimates were based on the assumption that future DWRs would be improved according to the projected NIWR; otherwise, the study area will experience acute water shortages in the future. The DWR 5 and DWR 10 were nearly identical across the time slices or RCPs; therefore, future DWRs should be revised according to NIWR for the 10-year return period drought. Agronomy 2021, 11, x FOR PEER REVIEW 13 of 18

Discussion
During the last 30 years, the seasonal and annual average temperatures across Punjab have changed at a rate of 0.020-0.60 °C/year [52]. The magnitudes and change rates of Tmax and Tmin projected in our study were in good agreement with the published literature [41,[52][53][54]. Figures 2 and 3 indicate that the future wheat season will experience stronger and warmer shifts than the rice season. Most weather stations in Punjab have detected substantial wheat seasonal warming during the past 30 years [52][53][54]. Ahmad and Choi [41] linked steady warming of the UCC command area since 1980 to a temperature rise phenomenon, notably witnessed in the winter months. We predicted a dramatic increase in the monsoon P magnitude or intensity in the rice season, particularly in the 2060s, while the wheat seasonal P remained stable in the 2030s and declined in the 2060s. Past trends and future projections of rainfall change patterns in Punjab also support the plausibility of receiving excessive monsoon rainfall in the rice season [7,17,41,52].
Climate change could have serious repercussions for the rice-wheat system of Punjab, Pakistan, impacting yields, growing season lengths, warming-driven ETc, ER, and NIWR behaviors, and eventually influencing the irrigation scheduling of the crops [3,6,7,17]. The IQRs in Figures 2 and 3 show that future Tmin change would be higher than the Tmax change for both seasons. Tmin is often associated with the night, and a conspicuously high nighttime temperature can drastically limit the spikelet fertility, grains per spike, and grain size in wheat [3,8]. A surge in nighttime warming also accentuates detrimental soil water stresses, and each increase of 1 °C can impart a 10% drop in rice yield [52].
Future seasonal P variation would directly impact the irrigation demands and available water supplies. Unexpectedly longer warmer spells with scarce rainfall during future wheat seasons would rapidly deplete the soil moisture, resulting in frequent irrigations

Discussion
During the last 30 years, the seasonal and annual average temperatures across Punjab have changed at a rate of 0.020-0.60 • C/year [52]. The magnitudes and change rates of T max and T min projected in our study were in good agreement with the published literature [41,[52][53][54]. Figures 2 and 3 indicate that the future wheat season will experience stronger and warmer shifts than the rice season. Most weather stations in Punjab have detected substantial wheat seasonal warming during the past 30 years [52][53][54]. Ahmad and Choi [41] linked steady warming of the UCC command area since 1980 to a temperature rise phenomenon, notably witnessed in the winter months. We predicted a dramatic increase in the monsoon P magnitude or intensity in the rice season, particularly in the 2060s, while the wheat seasonal P remained stable in the 2030s and declined in the 2060s. Past trends and future projections of rainfall change patterns in Punjab also support the plausibility of receiving excessive monsoon rainfall in the rice season [7,17,41,52].
Climate change could have serious repercussions for the rice-wheat system of Punjab, Pakistan, impacting yields, growing season lengths, warming-driven ET c , ER, and NIWR behaviors, and eventually influencing the irrigation scheduling of the crops [3,6,7,17]. The IQRs in Figures 2 and 3 show that future T min change would be higher than the T max change for both seasons. T min is often associated with the night, and a conspicuously high nighttime temperature can drastically limit the spikelet fertility, grains per spike, and grain size in wheat [3,8]. A surge in nighttime warming also accentuates detrimental soil water stresses, and each increase of 1 • C can impart a 10% drop in rice yield [52].
Future seasonal P variation would directly impact the irrigation demands and available water supplies. Unexpectedly longer warmer spells with scarce rainfall during future wheat seasons would rapidly deplete the soil moisture, resulting in frequent irrigations with smaller water depths. Since wheat production in the study area is entirely groundwaterdependent, warming-driven ET c would instigate rigorous groundwater exploitation for irrigation, straining the dwindling groundwater resources [6,7]. A higher monsoon P might be beneficial for future rice, but it would also accentuate the waterlogging and drainage problems. Ahmad et al. [7] argued that due to uneven distribution, the anticipated rise in the monsoon P would not compensate for monthly increases in ET c driven by warming, necessitating intensive irrigation applications.
We predicted a significant increase in the ET o , ET c , and subsequent NIWR of both crops as a result of future climate warming, which contradicts recent findings [3,6,7,17]. The ambient temperature primarily drives the crop phenological process, and its increase often leads to a shorter growth span [8]. Furthermore, the ET o is highly sensitive to the temperature and R n fluctuations in the study area [22], and the R n reduction could offset the warming-driven ET o and ET c rates [2,37,38]. All of the aforementioned studies attributed the ET o , ET c , and NIWR declines to a combination of growth span shrinkage and R n reduction. This study mainly focused on the worst-case scenario, which may serve as the foundation for future mitigation and adaptation strategies. Hence, climate warming impacts on crop growth span shrinkage were not evaluated in this study. We observed that seasonal R n declines primarily occurred during the rice season, but these were insufficient to offset warming-driven ET c and NIWR increases.
The study area was predicted to receive more rainfall in the future, but ER contributions to the PIWR and NIWR were limited due to intense and erratic rainfall events. This was especially noticeable during the rice season when, despite receiving heavy monsoon rainfalls, both the PIWR and NIWRE showed substantial positive shifts, and future wheat and rice required 1-2 and 2-5 more irrigations per season, respectively.
Since 1982, canal water supplies have been dwindling across Punjab, and crops are increasingly relying on groundwater and ER for their survival [10]. Our results suggested the availability of surplus rainwater because of the mismatch between rainfall occurrences and the duration of peak irrigation demands. Farmers should be encouraged to adopt rainwater harvesting to harness surplus rainwater and excess drainage from fields into public and community-operated wells. Pakistani farmers are usually unaware of actual crop water requirements and over-irrigate their crops for yield maximization [13,55]. This malpractice should be discouraged, and the adaptation of water and resource conservation strategies such as zero tillage, precision land leveling, and the furrow bed method of irrigation could play a vital role in sustaining rice and wheat production under climate change threats.
From the policy point of view, future crop water requirements could be managed by improving conveyance and water use efficiencies, irrigation technology adaptation, and shifting toward low-water intensive crops. The experimental findings of Khan et al. [20] support our suggestion of frequent irrigations with smaller water depths, as they achieved the maximum wheat yield by applying four irrigations, each 60 mm in depth, at a potential soil moisture deficit of 45 mm. Soomro et al. [37] concluded that rice cultivation on beds is a viable and efficient irrigation practice to save up to 23% of water without compromising yields. Improving the conveyance efficiency of irrigation canals can save an additional 14.8 billion cubic meters of water [13], which can be used to meet future crop water requirements. The adoption of high-efficiency irrigation systems, like drip and sprinkler irrigation systems, can deliver on multiple fronts such as water-saving, fertilizer use reduction, crop diversification, and enhancing yields and water productivity [21].
Our study included some limitations and assumptions that need to be addressed before implementing future water management measures and irrigation schedules. Due to the lack of detailed data, we did not consider the spatial variability of the soil, climate, and groundwater depth, which can produce extremely heterogeneous soil water balances at farm levels, ultimately influencing the NIWR and irrigation scheduling. A generalized version of the field conditions was presented in our study, but actual farm-level situations may vary dramatically according to crop variety, irrigation method, water quality, and soil management practices. Future studies should focus on tackling the climate change issue for maintaining and maximizing crop yields and water productivity, given the limited availability of irrigation water. Pakistan should seriously improve and standardize its protocols regarding data collection and the reporting of distribution, consumption, and production of land and water at the canal command scale [13]. The incorporation of such information would play a critical role in devising mitigation strategies to offset the adverse climate change impacts.
Despite the shortcomings, our research can contribute to developing an effective irrigation scheduling strategy for sustaining Punjab's rice-wheat system in the face of climate change threats. We estimated the future DWR, which may be used to develop an integrated water resource policy based on the future NIWR, DWR, and ER for better distribution of limited water resources based on agricultural water demand and canal water supply regulation.

Conclusions
In this study, irrigation schedules were generated based on future water requirements of the wheat-rice system of Punjab, Pakistan. Future wheat and rice seasons were distinguished by distinctly high warming rates, with the T min remaining significantly higher than the T max . A sharp rise in the erratic monsoon P was envisaged during future rice seasons, while the wheat season P remained unchanged or dropped slightly. Unexpectedly longer warmer and drier spells during future wheat seasons would rapidly deplete the soil moisture, while future rice seasons would experience a mismatch between rainfall occurrences and the duration of peak irrigation demands because of erratic and uneven monsoon rainfall distribution. Climate warming instigated a substantial increase in the ET o and ET c rates, and to compensate for consequently higher NIWR, the wheat and rice required 1-2 and 2-5 additional irrigations per season, respectively. Future DWR should be revised according to NIWR for the 10-year return period drought to accommodate the additional 100-200 mm of irrigation water per season; otherwise, the study area could experience water shortages in the future. Surplus rainwater was available during rice season that could be harnessed by adopting rainwater harvesting techniques and be used later during peak demand periods. Frequent irrigations with smaller water depths or deficit irrigation could be a viable irrigation strategy for wheat, whereas rice cultivation on beds could be an efficient irrigation practice for sustaining yields with limited water resources. Devising irrigation schedules to maximize the rainwater consumption could alleviate the anticipated stresses on surface water. Farmers should be encouraged to adopt resource conservation strategies and field management practices aimed at conserving and improving soil organic matter and water holding capacities. Developing suitable crop cultivars with better heat and water stress tolerance could be a long-term and reliable solution for maintaining yields with a limited water supply under future climate warming.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11102006/s1, Figure S1: Comparison of CROPWAT simulated irrigation day for wheat during the baseline (1980-2010) and the recommended irrigation period which is represented by the horizontal dashed lines, Table S1: Summary of GCMs used in the study, Table S2: Representative soil characteristics used in the study, Table S3: Comparison of CROPWAT results in this study with the literature, Table S4: Estimated parameters and goodness of fit results based on Kolmogorov-Smirnov (KS) and Anderson-Darling (AD) tests for the nine probability distribution functions (PDFs) during the baseline (1980-2010), Table S5