2.1. Description of study area

The study was conducted for a region of Pakistan that is sparse in data and is located in subtropical zone of southern Asia. The climatic conditions of Pakistan are semiarid with hot summers and mild winters (^{Ahmad et al., 2018b}). The entire region of Pakistan was divided into 10 × 10 km grids with a total of 8188 grids to take into account the soil and climatic spatial heterogeneity.

2.2. HAPPI scenarios

Climate data of maximum temperature (Tmax), minimum temperature (Tmin), precipitation (Pr) and solar radiation were downloaded from the HAPPI project (^{HAPPI, 2013}), whose generated data describe how the climate might be different from the current one in a 1.5 and 2.0 ºC warmer scenario than the pre-industrial conditions (^{Mitchell et al., 2017}). The baseline period of the HAAPI dataset is 2006-2015, while the projected climate period ranges between 2030 and 2040 (^{IPCC, 2018}). The output of the following four GCMs under HAPPI scenarios were used to assess the climate impacts: the Community Atmosphere Model version 4 (CAM4), developed by the Federal Institute of Technology (ETH) Zurich with spatial resolution of 1.87º × 2.5º (^{Neale et al., 2010}); the European Centre for Medium-Range Weather Forecasts, version 6 (ECHAM6), developed by the Max Planck Institute for Meteorology, Hamburg, with spatial resolution of 1.87º ×1.87º, including a modified version of the land component (^{Stevens et al., 2013}); the Model for Interdisciplinary Research on Climate, version 5 (MIROC5), developed by the National Institute for Environmental Studies, Tsukuba, Japan with spatial resolution of 1.4º ×1.4º (^{Shiogama et al. 2014}); and the Norwegian Earth System model, version 1 (NorESM1), developed by the Climate Center with spatial resolution of 0.94º ×1.25º. The latter is an ocean model which includes an advanced module for aerosols and aerosol-cloud radiation interactions (^{Bentsen et al., 2013}).

The GCMs biases were corrected by the methods described by ^{Hempel et al. (2013)} and ^{Frieler et al. (2017)}. Climate data were generated from the ensemble simulations under various initial conditions. Climate data for Pakistan were extracted from global climatic data and disaggregated to 10 × 10 km grids. The annual mean ambient temperature (TAV) and annual amplitude (AMP) were also calculated from the climate data. The CO₂ mixing ratio for the baseline period was 390 ppm, while it was 423 and 486 ppm for the 1.5 and 2.0 ºC additional warming scenarios, respectively (^{Mitchell et al., 2017}).

2.3. Soil data

Soil properties data were acquired from ISRIC-World soil information (^{Hengl et al., 2014}). Soil profile data of 0-5, 5-15, 15-30, 30-60 and 60-100 cm depths were taken at a 1-km spatial resolution. The soil parameter includes saturation percent (SP%), bulk density (gcm^-3), cation exchange capacity (cmol kg^-1), organic carbon (%), and soil pH. Global soil data for each of the above parameters were aggregated to a 10-km resolution. The drainage upper limit (field capacity) was calculated as 0.50 × SP and the lower limit (permanent wilting point) as 0.25 × SP (^{Miller and Kissel, 2010}).

The soil profile data of 0-5 cm depth is displayed in Figure 1, which indicates that more available water and organic carbon were found in north Pakistan as compared to central and southern Pakistan. Higher bulk density and soil pH were observed in central and southern Pakistan, as compared to northern Pakistan. Higher percentage of silt and clay in soil were recorded in the Indus basin irrigation system.

Fig. 1 Spatial distribution of soil properties in Pakistan. 

2.4. Model calibration and evaluation

The CERES-Wheat model under the shell of DSSAT was calibrated and evaluated using field data. The experiment was conducted at Ayub Agriculture Research Institute, Faisalabad during the Rabi season of 2016-2017 and 2017-2018.

The responses of three wheat cultivars (Faisalabad-2008, Lasani-2008 and Sehar-2006) were evaluated under four different levels of nitrogen (0, 50, 100 and 200 kg ha^-1). The genetic coefficients of the model were adjusted at non-stressed treatment of 100 kg N ha^-1 for three cultivars using generalized likelihood uncertainty estimation (GLUE) following ^{Ahmad et al. (2018a)}. After calibration, the model was evaluated using experimental data for the period 2017-2018.

2.5. Impact assessment of 1.5 and 2.0 ºC additional warming scenarios using a gridded approach

The calibrated CERES-Wheat model was used for gridded climate analysis. The crop management practices and fertilizers recommended by the government of Pakistan such as nitrogen (N), phosphorus (P) and potassium (K) were used as input for the crop model. The management practices were divided into three categories: irrigated (four irrigations of 75 mm), rainfed (without irrigation) and partially irrigated (two irrigation of 75 mm) as shown in Table I. The spatial analysis of the model was carried out to estimate the impact of 1.5 and 2.0 ºC additional warming under HAPPI scenarios. The output yield of each grid was displayed via geospatial maps of Pakistan. The detail methodology is presented in Figure 2.

Fig. 2 Methodology for assessing the impact of 1.5 ºC and 2.0 ºC additional warming on wheat yield. 

The impacts of climate change on wheat yield were computed by the difference of yield with 1.5 and 2.0 ºC additional warming from the baseline yield, as described in Eq. (1):

where Y _c is the percent change in yield, Y _f is the future simulated yield, and Y _b is the baseline simulated yield.

Multiple linear yield estimation empirical models were developed by utilizing the climate variables for arid, semiarid and humid areas of Pakistan. Temperature, precipitation, and yield differences were computed for the baseline and future period of each GCM and the least absolute shrinkage and selection operator (LASSO) regression was employed to develop the empirical model. About 70% of the data was used for the development of the model, whereas the remaining 30% was used for validation following ^{Kohavi (1995)}.

3.1. Estimated genetic coefficient of the CERES-Wheat model

The estimated genetic coefficients of wheat cultivars during model calibration are displayed in Table II. The P1V and P1D coefficients are related to phenology when moisture is a non-limiting factor for germination and the rate of crop development depends upon temperature. The values of vernalization (P1V) were set to 9-0 days for spring wheat cultivars. The photoperiod requirements (P1D) were similar for Lasani-2008 and Sehar-2006, while lower for Faislabad-2008. A large thermal time of 430 ºC day^-1 at grain filling (P5) was estimated for Faislabad-2008 and Sehar-2006 as compared to Lasani-2006.

The G1 and G2 coefficients are related to yield traits and showed a compensatory effect since increasing the number of grains (G1) results in decreasing grain size (G2). G3 is related to plant height and biomass production and did not show any difference between wheat cultivars. The phyllochron interval (PHINT) reflects the appearance of leaves on stem and is temperature dependent. The wheat cultivars did not show any differences regarding this marker during the vegetative stage.

3.2. Calibration and evaluation of the CERES-Wheat model

The CERES-wheat model was calibrated using field observations on phenology, growth, and yield attributes (Table III). A close agreement was observed in simulation of anthesis and maturity days. The model simulated anthesis with no difference with field observations for Faislabad-2008 and Sehar-2006, whereas a difference of only one day was observed for Lasani-2008. The simulated maturity days were close to observed values for Lasani-2008 and Sehar-2006; however, the model showed a two-day early maturity for Faisalabad-2008 (Table III). The model over simulated the leaf area index maximum (LAI^max) for all cultivars with an error ranging from 5.76 to 9.25%.

The model also simulated wheat yield satisfactorily. A close match was found between simulated and observed values with an error ranging from 0.52 to 1.36%. A close fit was also recorded in the simulation of biological yield for Faislabad-2008 and Sehar-2006, while the model over simulated the biomass for Lasani-2008 with an error of 6.94%.

Evaluation results indicate that the model performed reasonably well for anthesis days at different levels of nitrogen for three cultivars of wheat, as shown in Table IV. Nevertheless, a four-day difference between simulations and field observations was recorded with the zero nitrogen dose for Sehar-2006. The simulation of maturity days was suitable, but a 2-3 days difference was recorded relative to the observed values at high dose of nitrogen (200 kg N ha^-1). The LAI simulation was satisfactory, but a higher error was recorded at zero and 200 kg N ha^-1. The model showed a close match between simulations and observed data at various doses of nitrogen and cultivars, with an error ranging from 0.59 to 11.60%. Regarding the biological simulation, the model overestimated the biomass at zero nitrogen, which could be due to presence of initial nitrogen.

3.3. Climate change projections under HAPPI scenarios

The projected winter temperature and precipitation in all GCMs under HAPPI scenarios displayed different patterns over the entire region of Pakistan (Figs. 3, 4 and 5). The spatial pattern of warming is similar under the 1.5 and 2.0 ºC additional warming scenarios, and the magnitude of warming is higher over the entire region of Pakistan under the 2.0 ºC additional warming scenario. The spatial pattern of Tmax and Tmin revealed a lower increase in northern Pakistan, while a higher increase was projected in central and southern Pakistan (Figs. 3 and 4). CAM4 showed a higher increase in Tmax and Tmin, whereas ECHM6 showed a lower increase in Tmax and Tmin over the entire region, especially in KP and Gilgit Baltistan. NorESM1 displayed a higher increase in Tmax in Punjab (Fig. 3), with increase in Tmin was lower (Fig. 4). The Tmax ensemble mean showed an increase over the entire region of Pakistan under the 1.5 ºC additional warming scenario, and a greater increase was observed in Punjab under the 2 ºC additional warming (Fig. 3). The Tmin ensemble mean showed a greater increase in Gilgit-Baltistan in both the 1.5 and 2.0 ºC additional warming scenarios (Fig. 4).

Fig. 3 Changes in Tmax during winter from baseline under the 1.5 ºC (top row panels) and 2.0 ºC (bottom row panels) additional warming scenarios in Pakistan. 

Fig. 4 Changes in Tmin during winter from baseline under the 1.5 ºC (top row panels) and 2.0 ºC (bottom row panels) additional warming scenarios in Pakistan. 

Fig. 5 Changes in precipitation during winter from baseline under the 1.5 ºC (top row panels) and 2.0 ºC (bottom row panels) additional warming scenarios in Pakistan. 

The precipitation trend showed an increase under the 1.5 and 2.0 ºC additional warming scenarios. Changes in the spatial pattern of precipitation was also similar in both scenarios in Pakistan. Higher precipitation was recorded in KP and Gilgit-Baltistan, whereas Punjab, Sindh and Baluchistan had lower precipitation (Fig. 5). Wheat is a winter (rabi) crop sown in November-December and harvested in March-April, therefore this crop fully relies on irrigation because of low precipitation during winter. Increased rainfall during critical stages (crown root initiation, flowering, and physiological maturity) enhances the yield (^{Mudasser et al. 2001}; ^{Ahmad et al. 2019b}). Thus, high precipitation in KP and Gilgit-Baltistan, as shown in Figure 5, enhanced the yield, while less precipitation in lower in Punjab, Sindh and Baluchistan caused a reduction.

Furthermore, HAPPI scenarios indicated that Tmax (Tmin) will increase by 2.60 ºC (0.66 ºC) under the 1.5 ºC warming scenario, whereas Tmax (Tmin) will increase by 1.45ºC (1.44ºC) under the 2.0 ºC warming scenario (Table V). The ensemble means showed that Tmax will increase by 0.85 ºC, Tmin will increase by 1.05 ºC and precipitation by 204 mm under 1.5 and 2.0 ºC additional warming during 2030-2040 (Table V). Thus, higher temperatures, coupled with less precipitation caused reduction in wheat yield, because higher temperatures increase levels of water stress in plant cells, crop water requirement, and respiration.

3.4 Impact of 1.5 and 2.0 ºC additional warming on wheat yield

The impacts of the 1.5 and 2.0 ºC additional warming scenarios include both positive and negative relative changes in wheat yield in different districts of Pakistan (Figs. 6 and 7). The spatial pattern of yield shows that higher yield was recorded in KP, Punjab and Sindh provinces, while lower yield in Baluchistan and Gilgit Baltistan in all GCMs under the 1.5 and 2.0 ºC additional warming scenarios (Fig. 6). A maximum wheat yield of 2500-3000 kg ha^-1 was recorded in KP, Punjab and Sindh. The high yield was a result of good soil quality and availability of water due to the presence of the Indus and Kabul rivers, which are near to the productive zone for wheat cultivation in current conditions and under the 1.5 and 2.0 ºC additional warming scenarios. New potential areas for wheat cultivation are emerging in Gilgit Baltistan, KP and Baluchistan under the 1.5 and 2.0 ºC additional warming scenarios (Fig. 6).

Fig. 6 Spatial patterns of mean wheat yield in different GCMs under baseline (top row panels), 1.5 ºC (middle row panels) and 2.0 ºC (bottom row panels) additional warming scenarios in Pakistan. 

The spatial pattern of changes in wheat yield (%) at district level showed a similar trend in both the 1.5 and 2.0 ºC additional warming scenarios, though a higher increase in yield was found in the 2.0 ºC warming scenario (Fig. 7). A reduction in yield was recorded in Punjab, Sindh and a few districts in Baluchistan, whereas the areas in KP, Baluchistan and Gilgit Baltistan showed an increase in yield under both the 1.5 and 2.0 ºC additional warming scenarios as compared to baseline. A high reduction in yield was recorded in CAM4 and MICROC5. CAM4 showed a decrease in yield in Baluchistan and in a few areas of KP under the 1.5 ºC warming scenario, whereas MIROC5 showed a decrease in yield in Baluchistan under both the 1.5 and 2.0 ºC additional warming scenarios (Fig. 7). ECHAM6 and NorESM1 displayed a decrease in yield in Punjab and an increase in yield in KP, Baluchistan and Gilgit Baltistan (Fig. 7).

Fig. 7 Mean percent change of wheat yield relative to the baseline in different GCMs under the 1.5 ºC (top row panels) and 2.0 ºC (bottom row panels) additional warming scenarios in Pakistan at district level. 

The yield will decrease by 3.25% and 4.75% in Punjab, while in Sindh it will diminish by 17.8% and 13.8% under the 1.5 and 2.0 ºC warming scenarios, respectively. However, the yield will increase by 13% in KP, while an increase of 9.4% and 15.3% in Baluchistan under the 1.5 and 2.0 ºC additional warming scenarios is projected. Gilgit Baltistan showed an increased in yield of 10.6% under the 1.5 ºC warming scenario and 98% in 2.0 ºC scenario. The higher increase in yield at Gilgit Baltistan was due to the suitability of new areas for wheat production owing to additional warming of the 1.5 and 2.0 ºC scenarios as shown in Table VI.

Climate change impacts showed a reduction in yield of 3.2 and 4.7% in Punjab, and 17.8 and 13.8% in Sindh under the 1.5 and 2.0 ºC additional warming scenarios, respectively. However, the yield increases by 4.7 and 13% in KP, and 9.4 and 15.3% in Baluchistan under the 1.5 and 2.0 ºC additional warming scenarios, respectively. The reduction in yield was calculated using Eq. (1). Gilgit Baltistan showed greater increases in yield of 10.6 and 98% in both the 1.5 and 2.0 ºC additional warming (Table VI). The greater increased in yield under the 2.0 ºC additional warming could be due to the emergence of new areas for wheat production.

The GCMs exhibited an uncertainty in the evaluation of climate change impacts. The country average showed a higher reduction in yield of 11.7% in CAM4 and 9.8% in MIROC5 recorded under 1.5 ºC additional warming, whereas yield increases in all GCMs were regsitered under 2.0 ºC additional warming (Table VI and Fig. 8). A greater increase in yield was recorded in NorESM1 under both additional warming scenarios. CAM4 and MIROC5 have more interquartile ranges as compared to ECHAM6 and NorESM1 (Fig. 8). Considering values lower than the 25th percentile or higher than 75th percentile as outliers, ECHAM6 under the 1.5 ºC scenario (Fig. 8a) and MIROC5 under the 2.0 ºC scenario displayed more outliers (Fig. 8b).

Fig. 8 Changes in wheat yield (%) in different GCMs under (a) 1.5 ºC and (b) 2.0 ºC additional warming scenarios in Pakistan. 

3.5. Development of an empirical model for yield estimation

Empirical models for yield estimation were developed for arid (annual precipitation < 300 mm), semi-arid (300-1000) and humid (>1000) environments using LASSO regression (Table VII). The differences of Tmax, Tmin, precipitation and yields were calculated using the future GCMs relative to baseline. Seventy percent of the gridded data of all GCMs was used to develop the model and the remaining 30% to test the model. The results indicated that all climate variables significantly contributed to the development of the model, except for Tmin for arid environment and precipitation for humid environment (Table VII).

The empirical models for estimating wheat yield for arid (Y _a ), semi-arid (Y _sa ) and humid (Y _h ) environments are given in Eqs. (2), (3), and (4), respectively:

where Y is in kg ha^-1, Tmax and Tmin are in ºC and P is in mm. The results indicate that a rise of 1 ºC in Tmax will diminish the yield by 22.9 kg ha^-1 in an arid environment, 119 kg ha^-1 in a semi-arid environment and 264 kg ha^-1 in a humid environment, while a 1 ºC rise in Tmin will increase the yield by 8 kg ha^-1 in an arid environment, 73 kg ha^-1 in a semi-arid environment and 288 kg ha^-1 in a humid environment

The performance of the developed empirical model for arid environment was adequate, since the root mean square error (RMSE) value was 102.9 kg ha^-1, with an R² of 0.71. The results for semi-arid environments were also good with RMSE values of 128.2 kg ha^-1 and R² of 0.74. The developed model for humid environment showed an RMSE value of 147.5 kg ha^-1 and R² of 0.68. The developed models will be useful for stakeholders such as policy makers, researchers, farmers and academics for the estimation of yield in changing climate conditions. The yield change (kg ha^-1) can be estimated by knowing the change in temperatures and precipitation. The impact of unit change in Tmax, Tmin and precipitation on wheat yield can be calculated for arid, semi-arid and humid environments by using the developed empirical models shown in Eqs. (2), (3), and (4), respectively.