Future Wheat Yield Variabilities and Water Footprints Based on the Yield Sensitivity to Past Climate Conditions

Abstract

This article analyzed the wheat yield variabilities and water footprints under projected future climate based on wheat yield sensitivity to past (1980–2017) trends of maximum temperature (T_max), minimum temperature (T_min), solar radiation (R_n) and rainfall (P) for the semi-arid condition of Punjab, Pakistan. The past and projected future trends of the climate variables were identical featuring prominent T_min rise than T_max accompanied by R_n and P declines. Based on the past influences and the projected future trends of the climate variables, T_min was the principal driver of wheat yield decline followed by the T_max, P and R_n. Wheat evapotranspiration and net irrigation water requirement declined due to the influence of both decrease of growing season length and R_n for the first half of 21st-century, and this trend reversed during the second half of 21st-century. These phenomena were caused by different degree of climate warming between the two-time slices. Due to future wheat yield decline and a gradual rise in green water scarcity the blue and grey water footprints rise consistently particularly during 2nd-time slices. CO₂ enrichment showed limited mitigation potential of nullifying the warming-induced threats on future wheat yield and water footprints.

1. Introduction

Multi-sectoral impact assessment and mitigation of climate change had been an inspiration for numerous scientific explorations in the near past and will continue in the future. Among others, the agriculture sector is often regarded as the prime target of adverse climate change impacts due to its undeniably higher climate-reliance for its sustenance [1,2]. A rise is expected in frequencies and intensities of climate-driven natural calamities like floods and droughts, governed by climate warming and rainfall variations, all around the globe [3].

From the standpoint of an agriculture-dependent economy like Pakistan such climate change trends are eminently alarming due to the country’s over-reliance on irrigated agriculture and low adaptive capacity. In 2012, Pakistan ranked 8th among 180 other nations at Germanwatch Global Climate Risk Index; and was placed at 16th position among 170 other nation of the world at Maplecroft Food security Index [4,5]. Ineffective use of scanty water supplies and climate change are the driving factors responsible for low agricultural production and food shortage in the country [6,7]. Wheat is the staple crop; accounting for more than 50% of daily calorific intake and its availability and accessibility dictates the country’s food security conditions [8]. According to the World Food Program, despite having a good wheat harvest in 2014, up to 47% of the country’s population was food insecure led by widespread malnutrition, uneven food distribution and water shortage [9,10]. Thus, understanding the future wheat yield variabilities and water consumption patterns is critical for food security of Pakistan [10,11,12].

Climate change can affect both crop yields and water consumption patterns simultaneously. For instance, climate warming could enhance the crop water requirements/consumptions, while inducing yield declines at the same time [2]. During the last 30 years, Pakistan has witnessed a prominent surge in the annual and wheat seasonal temperatures across most weather stations [13,14]. Regarding the rainwater availability, different climate models showed contrasting tendencies of both reductions and increases in the annual rainfalls across the country [9]. Laghari, et al. [15] suggested that accelerated glacier melting might result in a short-term increase in surface water supplies in Indus river of Pakistan, the water availability would significantly decline by the 2045–2065. Agricultural productivity across various agro-ecological zones of Pakistan in terms of yields of staple crops such as wheat and rice is also projected to decline in the future [9,11,12].

For Pakistan, the future wheat yield decline and water-stress risks are most likely to be accentuated because of already existing wide mismatches between crop water requirements and rain and irrigation water availabilities [9,16,17]. The country relies on irrigated agricultural lands lying under arid to semi-arid conditions for wheat production and climate-related yield decline threats are very severe for such lands [11,12]. Food security of the growing population of Pakistan is counting on sustainable and/or improved wheat yield production from these arid and semi-arid lands facing extreme climate variations [7,9].

Water footprint (WF) concept offers a calculation framework to quantify the volumetric water consumption per unit of a crop yield production as well as the volume of water polluted during the crop production process [18,19]. In the agriculture sector, the total volumetric water consumption per unit of crop yield production (m³ton^–1) is categorized as the green water footprint (WF_green), blue water footprint (WF_blue) and grey water footprint (WF_grey). The WF_green and WF_blue represent the rain and irrigation water consumption, respectively; whereas, the WF_grey accounts for the water usage for dissipating the pollutant loads generated from various material inputs such as fertilizers, pesticides and herbicides during the crop production process [20,21,22]. Approximately 98% of Pakistan’s wheat production comes from the Indus basin which has one of the highest WF_blue in the world [22,23].

Understanding the source-wise (rain and irrigation) water consumption in crop production plays a vital role in agricultural water management [24,25]. Crop WF could be used as a multidimensional index to capture the climate-driven variations in green-, blue- and grey-water consumptions for future wheat yield production under climate change. Therefore, this study mainly aimed at predicting the future wheat variability and WFs under semi-arid conditions of Punjab which is the largest agricultural province in Pakistan. Specifically, the study objectives include: (1) to identify the most influential climate variable based of wheat yield sensitivity to past climate conditions; (2) to detect the detrimental climate change impacts on wheat yield; and (3) to understand the water consumption patterns associated with future wheat yield production for better management of green and blue water resources under semi-arid conditions of Punjab, Pakistan.

2. Materials and Methods

2.1. Study Area

In Pakistan, Punjab province is the largest wheat producer in terms of cultivation area and yield (Figure 1a). This work focused on the district Faisalabad, Punjab, Pakistan. The Pakistan Meteorological Department (PMD) has established a regional agro-meteorological (RAMC) experimental site in this district to monitor the climate-induced influences on wheat yield variations (Figure 1b). The study district was selected due to the data available from the wheat field trials conducted at the RAMC site from 2004–2017 (annual crop reports available at http://namc.pmd.gov.pk/crop-reports.php). Winter and summer are the cropping seasons in the study area lasting from November–April and May–October, respectively. Wheat is the major winter crop and its sowing starts from middle of November to the end of December and harvesting lasts from mid-March to mid-April [14]. In this area, canal water is the primary irrigation water source as the poor-quality groundwater is not suitable for irrigated agriculture [26,27]. Arid to semi-arid climate dominates the agricultural lands and up to 70% of the annual rainfall occurs as monsoon rainfalls during the summer season. Average wheat seasonal maximum and minimum temperatures ranged from 23–27 °C and 3–8 °C, respectively; whereas wheat seasonal cumulative rainfall varied from 85–210 mm. Well-drained, sandy clay loam is the dominant soil texture in this area [28]. The methodology followed to achieve the study objectives is explained in Figure 2.

2.2. Sensitivity Analysis of Wheat Yield to Climate Variables

A sensitivity analysis was conducted to identify the most influential climate variable dictating the past wheat yield trends in the study area. PMD provided a daily climate data set comprising of maximum temperature (T_max), minimum temperature (T_min), solar radiation (R_n), relative humidity (RH), wind speed (u₂) and rainfall (P) from 1980–2017. For the study district, the wheat production, cultivation area and yield (Y) data from 1980–2017 were also extracted from the Crop Reporting Service database, Government of Punjab, Pakistan (http://www.crs.agripunjab.gov.pk). Wheat seasonal (Nov–Apr) averaged T_max and T_min and seasonal cumulative R_n and P were included in the wheat Y sensitivity analysis. Change trends of wheat Y and the climate variable during the 1980–2017 were examined via simple linear regression with the year as an independent variable using the Equation (1):

$${\mathrm{V}}_{\mathrm{x}}={\mathrm{T}}_{\mathrm{x}}\times \mathrm{Year}+{\mathrm{I}}_{\mathrm{x}}$$

(1)

where V_x represents the climate variable or wheat Y; T_x is the slope of the regression line representing the trend of the climate variable and I_x represents the intercept of the regression model. Statistical significance of the T_x was evaluated at 5% significance level.

Crop yields mostly exhibit increasing trends (T_Y) due to the improvements which occurred over time in the crop production and management practices and seed quality etc. [29]. The T_Y (kgha⁻¹/year) comprised of a climate component (T_{Y_Cli}) which is only influenced by the climate factors and a management component (T_{Y_Man}) which defines the collective influences of crop production and management practices and other random factors (negligible). In order to establish climate-yield relationships, the T_{Y_Cli} and T_{Y_Man} must be isolated. The T_{Y_Man} can be removed from the T_Y by adopting de-trending techniques in which the inter-annual sensitivity of crop Y to climate variables is considered as a proxy for long-term climate change influences [30]. The first-difference method is a common de-trending approach which describes the absolute difference in Y or climate variables between two consecutive years [31]. We used the first-difference approach to de-trend the wheat T_Y during the 1980–2017 for the subsequent isolation of T_{Y_Cli} and T_{Y_Man}. The wheat T_Y was divided into the two components as follows [32,33]:

$${\mathrm{T}}_{\mathrm{Y}}={\mathrm{T}}_{\mathrm{Y}\_\mathrm{Cli}}+{\mathrm{T}}_{\mathrm{Y}\_\mathrm{Man}}+\mathrm{e}$$

(2)

where T_Y is actual wheat Y trend (kgha⁻¹/year); T_{Y_Cli}, T_{Y_Man} and e represent the T_Y component influenced by the climate change, crop production and management practices and random factors (negligible), respectively. Wheat Y sensitivity to the climate variables was defined using the simple multiple linear regression. A relationship was established between the first-difference values of wheat Y and climate variables as follows [30,31,34]:

$$\Delta \mathrm{Y}={\mathrm{S}}_{\mathrm{Tmax}}\times \Delta {\mathrm{T}}_{\max }+{\mathrm{S}}_{\mathrm{Tmin}}\times \Delta {\mathrm{T}}_{\min }+{\mathrm{S}}_{\mathrm{Rn}}\times \Delta {\mathrm{R}}_{\mathrm{n}}+{\mathrm{S}}_{\mathrm{P}}\times \Delta \mathrm{P}+\mathsf{\epsilon }$$

(3)

where ΔY, ΔT_max, ΔT_min, ΔT_Rn and ΔT_P represent the first-difference values of the corresponding variables; S_Tmax, S_Tmin, S_Rn and S_P are the coefficients of regression model representing the wheat Y sensitivity to the T_max (kgha⁻¹/°C), T_min (kgha⁻¹/°C), R_n (kgha⁻¹/MJm⁻²) and P (kgha⁻¹/mm), respectively; whereas, the ε is the intercept of the regression model. Finally, the T_{Y_Cli} was estimated using the following [31]:

$${\mathrm{T}}_{\mathrm{Y}\_\mathrm{Cli}}={\mathrm{S}}_{\mathrm{Tmax}}\times {\mathrm{T}}_{\mathrm{Tmax}}+{\mathrm{S}}_{\mathrm{Tmin}}\times {\mathrm{T}}_{\mathrm{Tmin}}+{\mathrm{S}}_{\mathrm{Rn}}\times {\mathrm{T}}_{\mathrm{Rn}}+{\mathrm{S}}_{\mathrm{P}}\times {\mathrm{T}}_{\mathrm{P}}$$

(4)

where the T_Tmax, T_Tmin, T_Rn and T_P represent the trends of respective climate variable from the Equation (1) during the 1980–2017; while other parameters are same as defined in the Equation (3). The T_{Y_Man} was estimated indirectly using the Equation (2). Different climate variables could have varying relative contribution towards T_{Y_Cli} depending on their influences and change trends. For instance, the relative contribution of T_max (RC_Tmax) towards wheat T_{Y_Cli} from 1980–2017 in the study district was estimated as [31]:

$${\mathrm{RC}}_{\mathrm{Tmax}}=\frac{{\mathrm{S}}_{\mathrm{Tmax}}\times {\mathrm{T}}_{\mathrm{Tmax}}}{\left|{\mathrm{S}}_{\mathrm{Tmax}}\times {\mathrm{T}}_{\mathrm{Tmax}}\right|+\left|{\mathrm{S}}_{\mathrm{Tmin}}\times {\mathrm{T}}_{\mathrm{Tmin}}\right|+\left|{\mathrm{S}}_{\mathrm{Rn}}\times {\mathrm{T}}_{\mathrm{Rn}}\right|+\left|{\mathrm{S}}_{\mathrm{P}}\times {\mathrm{T}}_{\mathrm{P}}\right|}\times 100\%$$

(5)

where all the parameters in Equation (5) are same as defined in the Equations (1)–(4). Finally, the relative contribution of the four key climate variables was estimated using the Equation (5).

2.3. Future Climate Change Projection

In this study, future climate change projections were also included to predict their influences on wheat Y variations and water footprint. Outputs from nine global circulation models (GCM), listed in

Table 1. Summary of global circulation models (GCMs) used in the study.

GCM	Institution	Horizontal Resolution
BCC-CSM-1.1	Beijing Climate Center, China Meteorological Administration	~1.25° × 1.875°
BNU-ESM	College of Global Change and Earth System Science, Beijing National University	~2.8° × 2.8°
GFDL-ESM2M	NOAA/Geophysical Fluid Dynamic Laboratory (GFDL)	~2.0° × 2.5°
CCSM4	US National Center for Atmospheric Research	~0.9° × 1.25°
HadGEM2-ES	UK—Meteorological Office—Hadley Center	~1.25° × 1.875°
inmcm4	Russian Institute of Numerical Mathematics (INM)	~1.5° × 2.0°
MIROC5	University of Tokyo, Japanese National Institute for Environmental Studies (NIES), and Japan Agency for Marine-Earth Science and Technology (JAMSTEC)	~1.4° × 1.4°
MPI-ESM-LR	Max Plank Institute of Technology (low resolution)	~1.9° × 1.875°
MPI-ESM-MR	Max Plank Institute of Technology (mixed resolution)	~1.9° × 1.875°

, forced under two representative concentration pathway scenarios (RCP 4.5 and 8.5) were employed to produce the future climate during two time slices: 2030s (2021–2050) and 2060s (2051–2080). The biases in the GCM-simulated and station-observed time series of T_max, T_min, R_n, RH, u₂ and P during the baseline period of 1980–2010 were corrected through the quantile mapping (QM) approach.

QM approach tailors the cumulative distribution functions (CDF) of the GCM-simulated and station-observed time series of a specific climate variable during the baseline so that both the time series exhibit consistent means, variances and other statistical properties. Same empirical CDFs are then applied to the GCM-simulated future time series. More precisely, the misrepresented values in the GCM-simulated time series are replaced by the station-observed time series values given both the values have identical percentiles. Extensive discussion about the working principals, applications, and extensions of the QM and its performance comparison with other bias-correction approaches can be referenced from the literature [29,35,36,37]. In this work, a daily QM based bias-correction of the six major climate variables was conducted during the baseline and future time slices. CDFs of the daily T_max, T_min, R_n and P for the wheat season (Nov–Apr) during the baseline (1980–2010) before and after the application of QM are shown in Supplementary Figure S1. The bias-corrected outcomes of wheat seasonal T_max, T_min, R_n and P are also presented in Supplementary Figure S2.

2.4. Wheat Yield, Evapotranspiration and Net Irrigation Water Requirement

The water-driven Aquacrop model v5.0 developed by Food and Agriculture Organization (FAO) was used to simulate the soil water balance, wheat growth and Y process while maintaining robustness, simplicity and accuracy [38,39]. By incorporating the atmospheric evaporative demands via FAO Penman–Monteith reference crop evapotranspiration (ET_o) and keeping track of incoming and outgoing water fluxes, the Aquacrop executes a daily water balance to determine actual crop evapotranspiration (ET); which in turn is further subdivided into transpiration (Tr) and evaporation (E) depending on the crop canopy cover and soil moisture status as follows [40,41,42]:

$$\mathrm{ET}=\mathrm{E}+\mathrm{Tr}$$

(6)

$$\mathrm{E}={\mathrm{K}}_{\mathrm{r}}\times {\mathrm{K}}_{\mathrm{e}}\times {\mathrm{ET}}_{\mathrm{o}}$$

(7)

$$\mathrm{Tr}={\mathrm{K}}_{\mathrm{S}}\times {\mathrm{Kc}}_{\mathrm{Tr}}\times {\mathrm{ET}}_{\mathrm{o}}$$

(8)

where: K_r is the evaporation reduction coefficient, varying between 0 and 1, with the smaller value being indicative of insufficient water availability in the top 30 cm soil layer relative to the atmospheric evaporative demands; K_e is the soil evaporation coefficient, varying between 0 and 1, with smaller and larger values being proportional to and representatives of the fraction of soil surface covered by the canopy cover and exposed to the open atmosphere (no canopy cover), respectively; K_s is the water stress coefficient being equal to 1 and 0 at the field capacity and permanent wilting point, respectively; and the Kc_Tr is the crop transpiration coefficient (0–1). Aquacrop transforms the seasonal cumulative Tr from the Equation (8) into the above-ground biomass (B) by multiplying it with the water productivity (WP*). The harvestable proportion of the B, which is the final crop yield (Y), is estimated from the product of B and crop-specific harvest index (HI) as follows:

$$\mathrm{B}={ \mathrm{WP}}^{\ast }\times {{\displaystyle \sum }}^{ }\mathrm{Tr}$$

(9)

$$\mathrm{Y}=\mathrm{HI} \times \mathrm{B}$$

(10)

where: B and Y are tonha⁻¹; WP* is kg per m² per mm of water transpired and normalized for the atmospheric CO₂ concentration and ET_o, and HI is % [38,39,40,41,42].

In this study calibrated and validated Aquacrop input files, prepared using the experimental data collected from the RAMC Faisalabad site, were used to simulate the wheat ET, net irrigation water requirement (NIWR) and Y under the baseline and projected future climate. Detailed information about the experimental layout, wheat cultivars and soil type data used for Aquacrop calibration and performance evolution can be referenced from Ahmad, et al. [43]. The calibrated Aquacrop model was run in conjunction with bias-corrected climate change projections to simulate the future wheat NIWR, ET and Y. The sowing date was fixed as 20 November and the future NIWR was estimated by generating an automatic irrigation schedule in Aquacrop irrigation module with no water stress. Irrigation was set to trigger automatically when 50% of the total available water in the soil has been depleted to bring the root zone total soil moisture back to field capacity. Two simulation sets were run for modelling the wheat Y responses to atmospheric CO₂ concentration: the first accounting for only the climate variations while keeping the CO₂ concentrations at the baseline level of 390 ppm; and the second accounting for the combined effect of climate variations and elevated CO₂ concentrations according to RCP 4.5 and 8.5 storylines [43].

2.5. Wheat Water Footprint

Frequently used [25,39,44] water footprint concept [18,19] was employed in this study to examine the wheat Y and water-usage (both P and NIWR) responses to projected future climate change. The term “water footprint” refers to the water consumption per unit production of a specific product. In the agriculture sector, total water footprint (WF_total) designates the volumetric water consumption per unit of crop yield production (m³ton⁻¹); and is further categorized as the WF_green, WF_blue and WF_grey. The WF_green and WF_blue represent the P and irrigation water consumption, respectively; whereas, the WF_grey accounts for the water usage for dissipating the pollutant loads generated from various material inputs such as fertilizer, pesticide and herbicide applications during the crop production process [20,21,22]. The WF_blue and WF_green of future wheat Y production in the study area were estimated as follows [20,21,45]:

$${\mathrm{WF}}_{\mathrm{blue}}=10\times \frac{{\mathrm{ET}}_{\mathrm{blue}}}{\mathrm{Y}}=10\times \frac{\mathrm{NIWR}}{\mathrm{Y}}$$

(11)

$${\mathrm{WF}}_{\mathrm{green}}=10\times \frac{{\mathrm{ET}}_{\mathrm{green}}}{\mathrm{Y}}=10\times \frac{\left(\mathrm{ET}-\mathrm{NIWR}\right)}{\mathrm{Y}}$$

(12)

where ET_blue and ET_green represent the P and irrigation water loss via evapotranspiration, respectively; factor 10 transforms the water depth (mm) into the water volume per area (m³ha⁻¹); Y is the Aquacrop-simulated wheat yield (tonha⁻¹) from the Equation (10). The degradative water usage for the wheat production-related pollutants dilution was estimated as follows [20,21,46]:

$${\mathrm{WF}}_{\mathrm{grey}}=\frac{\mathsf{\alpha }\times \mathrm{AR}}{{\mathrm{c}}_{\max }\times {\mathrm{c}}_{\mathrm{nat}}}\times \frac{1}{\mathrm{Y}}$$

(13)

where AR is the chemical (fertilizer, pesticide and herbicide etc.) application rate (kg ha⁻¹); α is the leaching-runoff fraction (%); and c_max and c_nat are the maximum acceptable and natural concentration (mg L⁻¹) of the considered pollutant in water, respectively.

Most WF_grey-related studies primarily focus on the pollution caused by Nitrogenous fertilizer application; which in turn governs the nitrate (NO₃) leaching into the surface and sub-surface water bodies [3,22,23]. In this study, considering the Nitrogen AR for irrigated wheat production in Punjab [47]; the Equation (13) was employed to estimate only the Nitrogen related WF_grey by using the α = 10%; c_max = 10 mgL⁻¹; and c_nat = 0 [22,23,46]. Finally, the WF_total were estimated as follows

$$\mathrm{W}{\mathrm{F}}_{\mathrm{total}}=\mathrm{W}{\mathrm{F}}_{\mathrm{green}}+\mathrm{W}{\mathrm{F}}_{\mathrm{blue}}+\mathrm{W}{\mathrm{F}}_{\mathrm{grey}}$$

(14)

The projected future wheat Y, NIWR and ET from the Equations (6)–(10) transformed the WF_total into WF_blue, WF_green, and WF_grey components by using the Equations (11)–(14). For climate change impact assessment on wheat Y and WFs, baseline period was from 2004–2017 for which the field trial data were available from the RAMC Faisalabad site. For WF_grey estimations, the Nitrogen AR was assumed to remain unchanged in the future. Finally, variations in the 2030s- and 2060s-WFs from the respective baseline-averaged values, with and without the CO₂ enrichment, were examined under the moderate (RCP 4.5) and extreme (RCP 8.5) emission scenarios.

3. Results and Discussion

3.1. Wheat Yield Sensitivity to Historic Climate

Inter-seasonal T_Tmax, T_Tmin, T_Rn and T_P estimated from Equation (1) for the wheat season (Nov–Apr) from 1980–2017 are shown in Figure 3. T_Tmin was statistically significant and the T_Tmax, T_Rn and T_P were statistically non-significant during the study period. T_Tmin and T_Tmax were increasing; whereas, the T_Rn and T_P were declining over time (1980–2017). Both T_Tmin and T_Tmax were greater than anticipated global average temperature change of 0.13 °C decade⁻¹ reported in the Intergovernmental Panel on Climate Change 5th Assessment Report (Figure 3a,c) [25,48]. In terms of magnitude and consistency, the T_Tmin was distinctly higher and prominent than those of the other climate variables. The respective change trends of T_min, T_max and P from their long-term average values were prominent after 2000; whereas, the negative R_n change trends were prominent during the last 5 years of the study period (Figure 3b,d,f,h).

The wheat production, cultivation area and Y trends during the 1980–2017 are shown in Figure 4. Wheat T_Y was increasing significantly at the rate of 0.04 tonha⁻¹/year or 36 kgha⁻¹/year which seems to contradict the regional and global scale studies [11,12,49,50] constantly highlighting the warming-driven wheat Y decline risks. Generally, the long-term crop T_Y is positive because of improvements in crop production and management practices and seed quality; thus hiding the negative climate change influences [33].

In our study area, the wheat T_{Y_Clim} was estimated to be −7 kgha⁻¹/year representing the wheat Y loss rate solely attributed to climate change in the study area. Despite the clearly evident harmful shifts in the thermal and rainfall regimes over course of the study period, the wheat T_Y was positive because of a gradual rise in cultivation area (Figure 4c) and improved production technologies. However, from 2010 and onwards the wheat Y change from its long-term average value was almost stagnant while the cultivation area was constantly rising (Figure 4c–f). This implied that T_{Y_Cli} has started to surpass the T_{Y_Man} and the current wheat cultivation and management practices would lag behind in sustaining the current wheat positive T_Y under future climate change trends.

Wheat Y sensitivity to climate variables and their RC, depending on their respective influences and trends were estimated using the Equations (3)–(5) and the results are shown in Figure 5. S_Tmin was highest followed by the S_Tmax, S_P and S_Rn; and wheat Y was negatively related to T_min and T_max and positively related to R_n and P, respectively. Compared to S_Tmax and S_Tmin, the S_P was noticeably lower as wheat production in this area is primarily irrigated (Figure 5a). Contribution rate analysis showed that all the four climate variables had negative RC in an order according to their respective influences. T_min was the major contributor because of its highest S_Tmin and T_Tmin followed by the T_max, P and R_n. Considering the difference between the S_Tmin and S_Tmax, the RC_Tmin and RC_Tmax differ greatly as the T_max, despite being the second most influential climate variable, had a much lower T_Tmax than that of T_min (Figure 5b). RC_P and RC_Rn were also negative regardless of their positive relations with wheat Y as indicated by the S_P and S_Rn in Figure 5a. Regardless of their mode of influence, the past trends of the four key climate variables were changing so as to contribute negatively towards wheat T_Y in the study area (Figure 5b).

3.2. Projected Future Climate Change

The projected changes in monthly, seasonal and annual T_max, T_min, R_n and P, after bias correction and averaged across nine GCMs, forced under RCP 4.5 and 8.5 during the baseline (1980–2010), the 2030s and 2060s are shown in Figure 6. Future climate of the study area in terms of T_max, T_min, R_n and P changed according to the past (1980–2017) climate change trends. Hence, the study area would undergo an intense and hotter thermal regime shift featuring prominent T_min rise than T_max accompanied by substantial seasonal cumulative R_n and P declines.

Anticipated climate warming (both T_max and T_min) was clearly evident at monthly, seasonal and annual scales irrespective of the GCM, RCP scenario (4.5 and 8.5) or future time slice (the 2030s and 2060s) (Figure 6a–d). The monthly T_max and T_min increased marginally during the first two months of the wheat season and then increased dramatically during the rest of the growing season. Change magnitude of the projected climate warming featured various noteworthy trends: T_min > T_max, wheat season > annual, RCP 8.5 > RCP 4.5 and 2060s > 2030s. The temperature change projections were supported by the last 30-year wheat seasonal average temperature inclines which were detected across most weather stations in Pakistan [13,14,51]. Trends and magnitudes of future climate warming projected in our study were also supported by the previous studies predicting the 21st-century climate of Pakistan [12,28,52].

Future wheat was particularly exposed to climate warming during the middle of the growing season from January to March. This could lead to a shortened grain filling duration but accelerating the gran filling rates which would not compensate for the probable Y reduction induced by shortened grain filling duration. T_min is often associated with nighttime and prominent nighttime temperature increase could affect the spikelet fertility, grains per spike, grain size, and quality. Generally, wheat is particularly vulnerable to heat stress at anthesis and grain-filling stages [50]; for instance, the wheat Y in Pakistan is very sensitive to the mid-March higher temperature during grain filling stage [12].

Future monthly R_n also followed the temperature change trends but in a downward direction. Compared to the baseline, the future monthly R_n declined slightly during the first half; and greatly during the latter half of wheat season (Figure 6c) regardless of time slice or RCP. Contrary to the temperature, future R_n decline was prominent at the annual scale for both RCPs and time slices. For RCP 4.5, the seasonal- and annual-scale future R_n-decline was higher than its RCP 8.5 counterpart during both the time slices. Overall, the projected R_n decline with respect to baseline was: Annual > wheat season, RCP 4.5 > RCP 8.5 and 2030s > 2060s (Figure 6c,d).

The magnitude and direction of monthly, seasonal and annual P change in future varied greatly across GCMs. The future P declined or remained unchanged for most of the GCMs while rising P tendencies were envisaged by few GCMs. However, the GCM-ensemble wheat seasonal and annual P decreased with respect to the baseline (Figure 6g,h). Future monthly P increased from Nov–Jan; and then it declined drastically during the rest of the wheat season for both RCPs and time slices (Figure 6g,h). Hence, the future wheat would face severe P shortages along with excessive warming and R_n declines during the middle of the growing season which could decelerate the photosynthesis rates, increase the tiller infertility rate at T_max = 36 °C and reduce the HI at T_min >14 °C [12,49].

Similar to the R_n decline, future drying was also pronounced at the annual scale. Future wheat seasonal P decline was notably higher under RCP 4.5 than RCP 8.5; while the annual P decline remained almost the same under the two RCPs during both time slices. On average, 45–53% and 13–30% seasonal P decline; whereas 56%–60% and 53%–59% annual P decline was predicted by the end of the 21st century under the RCP 4.5 and 8.5, respectively.

3.3. Future Climate Change Impacts on Wheat ET and NIWR

Wheat growing season length (GSL) recorded at the study site from 2004–2017 ranged between 116–145 days with an average of 130 days; while the seasonal cumulative growing degree days during this period ranged between 1545–1841 °C with an average of 1698 °C. From here on the term “baseline” refers to 2004–2017 period for which the wheat experimental data were available. Baseline average wheat seasonal cumulative ET_o, ET and NIWR (simulated by Aquacrop) and P received during the GSL were 342, 285, 188, and 98 mm, respectively. Figure 7 presents the Aquacrop simulated GSL, ET_o, ET, P and NIWR averaged during the baseline, the 2030s and 2060s. Projected wheat seasonal ET_o, ET, and NIWR variations from the baseline during the 2030s and 2060s under the two RCPs with and without the CO₂ enrichment are shown in Figure 8. The seasonal cumulative ET_o, ET and NIWR for the wheat declined in the future which was odd as climate warming is often associated with accelerated ET rates and consequently higher NIWR [1,25]. The three parameters showed higher declining tendency during the 2030s particularly under RCP 4.5. Overall, the future wheat seasonal ET_o, ET and NIWR had the minimum, medium and maximum declining likelihood, respectively, across the time slices or RCPs. Ambient temperature primarily controls the crop phenological progress and its increase often reduces the GSL due to rapid heat unit assimilations during crop growth stages [50]. Sultana, et al. [11] predicted a 12–14-day wheat-GSL shrinkage against future climate warming of up to 5 °C for the arid and semi-arid conditions of Pakistan. Bhatti, et al. [53] identified higher declining probabilities of ET and NIWR for the early maturing wheat in Punjab, Pakistan. In our study, future climate warming reduced the wheat-GSL by 10–12 days for RCP 4.5 and 13–20 days for RCP 8.5 by the end of 2080. This could explain the unexpected decline in the wheat seasonal cumulative ET_o, ET and NIWR.

Because of moderate climate warming and the subsequent GSL shortening, the wheat seasonal cumulative ET_o and ET declined only under RCP 4.5 during the 2030s; however, for rest of the cases, they increased because of intense climate warming despite a significant GSL shortening. During the 2030s, average ET_o change was –32 mm and +95 mm; and during 2060s, it was +65 mm and +125 mm for RCP 4.5 and 8.5, respectively (Figure 8a). Future wheat ET change patterns were comparable with that of ET_o but with slightly smaller change magnitudes. Contrary to conspicuously higher ET_o change, mildly positive changes in future wheat ET were projected regardless of the RCP or time slice. During the 2030s, average ET change was −20 mm and +1 mm; and during 2060s, it was +3 mm and +12 mm for RCP 4.5 and 8.5, respectively (Figure 8b).

Difference between the crop ET and seasonal P determines the NIRW; which, in our study, showed the maximum declining likelihood compared to the ET_o or ET. An average NIWR change of −32 mm and −45 mm during the 2030s; and −17 mm and +5 mm during the 2060s was shown for RCP 4.5 and 8.5, respectively (Figure 8c). The unexpectedly larger decline in 2030s-NIWR under RCP 8.5 was because of a positive seasonal-P change of +28 mm during this duration. The 2030s-P increased by 4% and 30%; and the 2060s-P decreased by 7% and 20% under RCP 4.5 and 8.5, respectively. It should be noted here that the future P-change trends discussed in Section 3.3 and shown as Figure 6h correspond to the general winter cropping season from Nov–Apr (180 days); whereas, the abovementioned future P-change trends correspond specifically to the wheat GSL (baseline average value = 130 days) as simulated by the Aquacrop.

Future wheat seasonal ET_o, ET and NIWR were most likely to increase during the 2060s particularly under RCP 8.5. Furthermore, the difference between the 2030s- and 2060s-IQRs of the ET_o, ET and NIWR change was larger for RCP 4.5 and smaller for RCP 8.5, respectively. Positive changes in wheat seasonal ET_o, ET and NIWR were projected to begin at the start of the 2060s for the RCP 4.5; and during the middle of 2030s for RCP 8.5, respectively. Results suggested that frequent irrigations would be necessary to compensate for the warming-driven ET rates during a shortened growth span in the second half of the 21st century. CO₂ enrichment did not induce any noticeable difference in the aforementioned ET_o, ET and NIWR change patterns during the 2030s; however, during the 2060s, the CO₂ enrichment partially contributed in limiting the warming-driven ET_o, ET and NIWR increase tendencies regardless of RCP scenario (Figure 8a–c).

In addition to the probable contraction in wheat-GSL, the R_n decline could also be a leading cause of the drops in future ET_o and ET identified in this study. Cloud cover, air pollutants and manmade aerosols affect the sunshine hours and R_n. An increased concentration of manmade aerosols and other air pollutants in the atmosphere can greatly reduce the incident R_n on a specific area [54]. This could explain the R_n decline detected in our study during the baseline and future time slices. Energy availability at the crop surface in the form of R_n plays a critical role in determining the water removal rate as ET_o or ET [55]. ET_o is highly sensitive to R_n and is directly related to R_n change patterns across various agro-ecological zones in Pakistan [14,56]. Based on a seasonal scale sensitivity analysis, we found that both the wheat ET_o and ET were highly sensitive and positively related to R_n followed by T_max, T_min and P, respectively (results not shown here). This implied that future R_n decline could also limit the warming-driven ET_o and ET increments. Our results showed that a combination of GSL and R_n decline would counteract the warming-driven ET_o and ET increments during the 2030s; but it would not withhold these ET_o and ET increments during the 2060s, ultimately resulting in higher NIWR during this duration. Similar conclusions were drawn by the recent studies from various parts of the world, in which the observed declines in ET and NIWR of various crops, despite a clear climate warming, were attributed either to GSL shortening and/or R_n declines [1,3,25,57].

3.4. Future Climate Change Impacts on Wheat Y and WFs

The average baseline wheat Y was 4.6 tonha⁻¹. Comparison between average wheat Y simulated by the Aquacrop during the baseline, the 2030s and 2060s under the two RCPs with and without CO₂ enrichment is shown in Figure 9. Moderate to severe wheat Y-declines at the baseline CO₂ concentration of 390 ppm were detected during the 2030s and 2060s, respectively, regardless of RCPs. A potential rise in CO₂ concentrations from 415–530 ppm (RCP 4.5) and from 420–750 ppm (RCP 8.5) partially and fully mitigated the wheat Y-declines prompted by future climate warming by end of the 2060s. Under RCP 8.5, despite the severe warming, 8–17% increments in future wheat Y were detected during the two time slices. Hence, by doubling the CO₂ concentration up to 750 ppm under RCP 8.5, beneficial CO₂ enrichment influences surpassed the negative climate warming impacts (Figure 9). These deductions were supported by findings of global-scale studies aimed at examining a composite climate change impact on future wheat Y resulting from the interplay between CO₂ enrichment and climate warming [58,59].

Wheat WFs for each season of the baseline period estimated using the Equations (11)–(14) are presented in Figure 10. The average WF_green, WF_blue, WF_grey and WF_total during the baseline were 210, 407, 347, and 964 m³ton⁻¹, respectively. Comparison between the average WFs during the baseline, the 2030s and 2060s are shown in Figure 11. Figure 12 shows the wheat WFs change from the baseline, averaged across nine GCMs, with and without the CO₂ enrichment during the time slices and RCP scenarios.

Without the CO₂ enrichment, the future WFs increased because of the wheat Y drop prompted by climate warming. Regarding the time slices, the 2060s-WFs increments were higher than those of the 2030s for both RCPs; whereas, for RCPs, the 8.5-WFs increments were higher than those of the 4.5 for both time slices except WF_green which showed an opposite changing tendency. Wheat WF_blue increased gradually to compensate for the green water scarcity originating from the projected drop in future P by the end of 2080. Under the enhanced CO₂ concentrations, the wheat-WFs reduced negligibly during the 2030s and notably during the 2060s particularly under RCP 8.5 because of wheat-Y increments stimulated by CO₂ enrichment (Figure 11).

Despite having similar change trends, change magnitudes of the three WFs were not the same. Without CO₂ enrichment, the 2030s-WF_green increments were larger than the 2060s-WF_green for both RCPs; indicating more green water availability per unit wheat Y or more Y production per unit of green water input. The 2030s duration featured moderate P increments which could explain the more green-water availability in this case. Because of higher wheat Y decline and a gradual rise in green water scarcity, the 2060s-WF_green constantly declined irrespective of the RCPs or CO₂ enrichment.

Contrary to the WF_green, the rises in WF_blue and WF_grey during the 2060s were higher than those of the 2030s because of the higher blue and grey water consumption per unit of the wheat Y production ultimately leading to higher WF_total. WF_blue and WF_grey increments during the 2030s were marginal because of a combination of moderate risks of climate warming, drying and wheat Y decline. During the 2060s, the warming-driven ET increments combined with a heightened risk of Y and P decline resulted in substantially higher blue and grey water consumption per unit of wheat Y production (Figure 12a–d). Among the three WFs, only the WF_grey exhibited a constantly increasing pattern irrespective of the time slice, RCP or atmospheric CO₂ concentration (Figure 12c). Although the WF_grey declined because of wheat Y improvement caused by CO₂ enrichment, its future changes remained consistently positive compared to the baseline. Results suggested that 50–150 m³ton⁻¹ of extra greywater would be required to dissipate the pollutant load generated from Nitrogen application during future wheat production (Figure 12c).

The WF_green, WF_blue and WF_total estimates produced in this study were slightly lower than the benchmark wheat WFs estimates at basin scale in Pakistan [22,23,60]. We calibrated the Aquacrop model using the field-scale wheat-trials data which were conducted under the optimum conditions such as without moisture, fertilizer or any other input restriction. Wheat Y from such trials could be very high compared to the farmer’s field and regional and basin-scale Y estimates. In our study, slightly higher wheat Y estimates were simulated by the Aquacrop, thus; producing slightly lower WFs values. Moreover, we only consider the optimum NIWR for WF_blue estimation without considering the conveyance losses or excessive field-scale irrigation water applications by the farmers which are an integral part of WF-estimation framework at a larger spatial scale. However, replacing the Aquacrop simulated wheat Y with regional statistical Y would produce similar WFs outcomes according to the benchmark values provided in the aforementioned studies. This work mainly examines the wheat-WFs change-patterns in place of their absolute values to approximate the future climate change impacts on wheat water management in the study area. Based on various climate change probabilities examined in this study, the most unlikely and extreme scenario of RCP 8.5 exhibited the prospects of CO₂ enrichment convincingly nullifying the future warming during the 2nd half of the 21st century. For the rest of the cases, the mitigation potential of CO₂ enrichment against future climate warming remained limited. This implied that CO₂ enrichment could partially compensate for the detrimental warming impacts on the future wheat Y and WFs in semi-arid regions of Punjab, Pakistan.

4. Conclusions

Wheat yield sensitivity to past climate change trends and its variability under projected future climate change was quantified for semi-arid conditions of Punjab, Pakistan. Identical past and projected future climate change trends were shown with wheat yield being inversely related to growing seasonal averaged T_min and T_max, and directly related to growing seasonal cumulative R_n and P. T_min was the principal climate driver contributing to the past and future wheat yield trends followed by the T_max, P and R_n. Warming-driven accelerated crop development process enhanced the yield failure risks in the future early-maturing wheat crop; however, the CO₂ enrichment partially compensated the detrimental warming impacts on wheat yield. Wheat ET and NIWR declined due to the influence of both decrease of growing season length and R_n for the first half of 21st-century, and this trend reversed during the second half of 21st-century. Rampant green water scarcity exerted extra pressure of blue water resources, particularly during 2nd-time slices while extra greywater was required to dissipate the pollutant load generated from Nitrogen application. Devising irrigation schedules to maximize the green water consumption could alleviate the anticipated stresses on blue water resources. Adaptation of field management practices to maintain and improve soil organic matter and water holding capacities can also reduce the grey and blue water consumptions by facilitating the nutrient and blue water availability to the crop. Development of suitable crop cultivars having better heat tolerance and R_n interception capacity could be a long-term and reliable solution for maintaining wheat yields under future climate warming.