Cotton is a high-value commercial crop farmed worldwide as a sustainable natural fiber source for the textile industry (Amanet et al., 2019; Aslam et al., 2020). It supplies raw materials to the expanding textile industry, cottonseed oil for culinary use, and protein-rich oil cake leftovers to livestock (Munir et al., 2020; Xun et al., 2021). Cotton cultivation is intricately associated with climatic conditions and irrigation, and the overall nitrogen released during the growing season has a significant influence (Wang et al., 2018). As a result, timely and accurate information on cotton’s spatial and temporal distribution pattern in the face of climate change is essential for crop management, regional crop response, and agricultural policy making (Li et al., 2011). Anthropogenic activities are thought to have contributed to a 1 °C increase in the global temperature compared to pre-industrial levels (Tarabon et al., 2020). Many studies based on the carbon emission scenarios estimated that the earth’s temperature could increase by 2–5 °C by the end of this century (Forster et al., 2007). Cotton growth is also affected by relative humidity in addition to temperature and precipitation, all known to impact cotton growth significantly (Li et al., 2019). Cotton grows successfully at an average temperature of 28 °C in China and 41 °C in Sudan during the crop season (Zahid et al., 2016). During the cotton growing season, the average temperature in the Indus River Basin region of Pakistan (IRBP) is around 37 °C. On the other hand, heat stress has historically been a serious barrier to cotton-producing countries, including Pakistan, Syria, and India. (Chen et al., 2020). Furthermore, temperature increases have resulted in increased cotton cultivation budgets in areas where cotton is already grown at temperatures close to 40 °C (Jans et al., 2021).

Climate change is one of the most significant issues affecting agricultural production systems (Javid et al., 2019). The impact of climate change on crop yields has been widely acknowledged at regional and global levels and found that climate factors influencing crop growth at regional and global levels (Schlenker and Roberts et al., 2009; Wheeler et al., 2013; Thornton et al., 2014; Asseng et al., 2015; Ray et al., 2015), as it is known to have a complex impact on crop growth and production (Rosenzweig et al., 2014). Dependence on precipitation and temperature directly affects crop growth (Schmidhuber et al., 2007). Both temperatures and the frequency of extreme weather events are increasing in most parts of the world, particularly in Pakistan (Naveed et al., 2021). The predicted climate changes in Pakistan are expected to impact crops yield in arid and semi-arid regions negatively (Ahmad et al., 2015; Qin et al., 2015; Rasul et al., 2016; Abbas et al., 2017). Drought has increased in recent years (Mazdiyasni and AghaKouchak, 2015), resulting in decreased crop yields (Lesk et al., 2016; Zipper et al., 2016). Historically, low seasonal rainfall and high temperatures with shifting trends at regional levels have been observed in the Indus River basin cotton belt. Simultaneously, adaptation measures such as shifting planting dates, developing new crop varieties, and changing crop growth patterns can reduce the harshness of climatic effects on agricultural production, including the negative effects on crop yields (Kumar et al., 2013; Osborne and Wheeler et al., 2013; Challinor et al., 2014). Farmers grow crops in watered, nutrient-rich soil, which causes changes in crop distribution patterns. Thus, it is very important to understand how the spatial distribution pattern of the cotton yield has been modified as a consequence of historical climate changes.

The historical relationships between cotton yield and climate are well understood using statistical models based on previous crop yield and climate observations (Lobell and Burke, 2010). Many studies have investigated the long-term impact of climate change on crop yields (Basche et al., 2016; Leng and Huang, 2017; Li et al., 2019), but relatively few studies have investigated crop yield responses to inter-annual climate change (Ray et al., 2015). Temporary changes in agricultural production, in particular, may have a more significant impact on food prices, farmer income, and food security than long-term changes in crop yields (Godfray et al., 2010; Hertel et al., 2010; Iizumi et al., 2013). However, due to changes in technology and spatial distribution patterns enabled by competitive land use, likely, long-term changes in crop yields will not be realized for some crop regions (Leng and Huang, 2017). As a result, to study the effect of the changes in cotton spatial distribution patterns can help not only in investigating the implications of the corresponding role in mitigating the effect of historical climate on production but also in comprehending the uncertainty in predicting future cotton yields. We have applied the moving average method to investigate the long-term stability trend in cotton yield variation. Simple least square regression and piecewise regression models were used to detect the trend-turning point of the corresponding climate variables. We performed a partial correlation analysis that excluded covariation between the studied climate factors. The long-term variability relationship between these variables was exposed by applying a 15-yr moving window analysis on standardized values of cotton yield and climate factors. Finally, a multiple regression model was used to identify the climate factors that significantly impact cotton yield at the tehsil level (administrative unit) in the IRBP.

The remainder of this paper will specifically address the three scientific questions: 1) How has cotton yield variability change in the IRBP over the past three decades? 2) To what extent can the combination of these three climate factors account for variation in cotton yield at the tehsil level? 3) How did the spatial distribution pattern of cotton yield in the IRBP affect the crop response to historical climate changes at the tehsil level?

The Indus River Basin cotton belt is in a fertile area of Pakistan that is home to almost 85% of Pakistan’s cotton production (Naveed et al., 2021), and has significantly benefited the country’s economy (Fig. 1). The region encompasses 142 075 km² and accounts for 18% of Pakistan’s land area. Geographically, it stretches from 24°56′48″N to 31°57′31″ N and 69°13′05″E to 72°55′41″E. The soil changes from clay loam to sandy, with clay dominating in the south (Naveed et al., 2021). From north to south, the elevation decreases in the IRBP. Cotton is grown in regions where the climate area is hot and low annual precipitation (Javid et al., 2019). Climate data for the monthly cotton growing season (May–September) for 84 IRBP’s tehsils were obtained (https://power.larc.nasa.gov/data-access), while the cotton yield and cultivated area data were obtained from regional offices and the Pakistan Bureau of Statistics (http://www.pbs.gov.pk) from 1990 to 2019. The crop reporting service of the Bureau of Statistics collects the crop’s yield and area information from their field assistants and at selected sample points before making the final estimates at the tehsil and zonal levels. Temperature, relative humidity, and precipitation were among the climate variables studied, as their spatiotemporal variations can affect cotton yield, justifying the importance of studying the impact of climate on cotton production (Lobell and Field, 2007). Cotton is a summer crop planted in April and May in the Indus River Basin and harvested once a year. The cotton crop cycle has a duration of 150–170 d. The cotton plant appears 6 d after sowing and starts flowering and budding at the end of June; its harvest begins in September.

Figure 1.  Study area indicating the district headquarters in the Indus River Basin, Pakistan. The abbreviated district names are defined as Jhang (JHG), Toba Tak Singh (TTS), Sahiwal (SHL), Pakpatan (PAK), Khanewal (KWL), Vehari (VHR), Multan (MUL), Muzafarghar (MGH), Lodhran (LDR), Bahawalnagar (BWN), Bahawalpur (BWP), Rajanpur (RJP), Rahimyar Khan (RYK), Ghotki (GHT), Nawabshah (NWS), Khairpur (KRP), Mirpur khas (MPK), Sanghar (SNG), and Umar kot (UKT)

To indicate the cotton crop area anomalies at the district level, we calculated the standardized cotton crop area for the study period given by Eq. (1).

where $ {Y}_{i} $ is the cotton crop area in the ith year of the given district, $\overline{Y}$ is the average area, and S.D. is the standard deviation of the cotton crop cultivated area of the respective district during 1990–2019.

Generally, standard deviations outside the range of standardized values from –1.6 to 1.6 are considered statistically significant at the 0.05 level (α = 0.05). Positive (negative) values indicate an increase (a decrease) in cotton crop area. To obtain further background information about the spatial distribution changes patterns of cotton crop areas within the district boundary, we applied the stability index (SI) proposed by Leng and Huang (2017) as Eq. (2):

where D represents the district, $ {A}_{t,y} $ the cotton crop area for the tehsil t during y year, k represents the number of tehsils in respective district, and $ {y}_{o}{}_{} $ is the reference year 2019. To measure the temporal stability, the coefficient of variation (CV) is a commonly used index (Challinor et al., 2014; Ray et al., 2015), while the SI calculates the stability of cotton crop area in both a spatial and temporal sense.

To show the overall spatial distribution variability, we calculated the percentage of relative variability of cotton yields. Further, for estimating the interannual cotton yield variability trend, we applied a 15-yr moving time window for each tehsil from 1990 to 2019 and obtained 16 samples of variability. We subtracted the mean yield from these samples and fit the result with a linear trend to calculate the percentage of annual increase or decrease in cotton yield. Cotton yield variations had both positive and negative values, so we applied a two-tailed t-test to test the significance of the linear trend. We analyzed the 30-yr cotton yield and, as per the requirement for robust statistics, we selected a 15-yr moving time window. To examine the sensitivity of the results of the selected time window, we performed additional analyses using smaller (10-yr) and larger (20-yr) moving time windows.

Previous studies have investigated the degree to which climate variables shift their trends using their time series data. Precipitation is the main climate factor that influences to change trend of other climate factors (Fu et al., 2021). To detect the trend turning point of precipitation for the cotton growing season during the study period, we applied two regression models. First, we used a least square regression model to detect the trend of the climate variables over the entire study period during the cotton growing season. The least-squares simple linear regression model is a frequently used method for calculating a linear temporal trend (Running and Nemani, 1988; Tabari and Hosseinzadeh Talaee, 2011). We were able to obtain the trend rate of climate variables and also test the significance of that trend using the least square regression equation given by Eq. (3):

where $ {y}_{i} $ are climate variables, $ {x}_{i} $ represents the time (the cotton growing seasons in year), a and b are the regression equation parameters (a is the intercept and b is the rate of change in trend), and $ {\varepsilon }_{i} $ is the random error.

Second, we fitted a piecewise linear regression model to explore the shifting trend in precipitation. The piecewise linear regression model is easily repeatable to generate a series of the sum of the square residuals that feature the best fit with one breakpoint, a technique commonly used in climate trend studies (Piao et al., 2011; Zhang et al., 2013). The trend-turning point (TTP) of precipitation was obtained by minimizing the residuals of piecewise linear fits, with the TTP significance being assessed by a t-test against the null hypothesis given in Eq. (4):

where $ {\;\beta }_{0},{\;\beta }_{1} $, and $ {\;\beta }_{2} $ are the fitted regression parameters, $ {\;\beta }_{0} $ is the intercept, $ {\;\beta }_{1} $ defines the trend before the turning point, and $ {\;\beta }_{1} $ $+ {\;\beta }_{2}$ is the trend after the turning point, while $ {y}_{i} $ represents the respective climate variable, $ {x}_{i} $ are cotton growing seasons (in year), $ \alpha $ is the estimated trend turning point, and $ {\varepsilon }_{i} $ is the error term.

We calculated the piecewise regression trend for temperature and relative humidity based on the precipitation estimated trend turning point. Then we applied the following multiple regression model to the impact of these climate variables on cotton yield in the IRBP for the entire study period, before and after the TTP.

To analyze the strength of the cotton yield variability within the cotton growing season (June–September), we analyzed the climate variables of temperature (T), humidity (H), and precipitation (P) in terms of their variability. We applied the correlation analysis to the standardized 15-yr moving observed values for all variables. Because climatic factors and their effects on crop yields can co-vary on an interannual time scale, knowledge of the actual effects of individual climate factors might assist in building more effective methods to adapt to the expected changes in climate (Leng et al., 2016). To remove the climate covariation, we applied partial correlation analysis. We also detrended both the cotton yield variability and climate variability time series and performed correlation analysis to check the strength of cotton yield relationship with climate factors. Should the strength of the relationship become weaker after eliminating the trends, this would indicate that the variability in cotton yield is strongly linked with the climate variability trends.

Next, we used the time series of the tehsil-level cotton yield as a dependent variable and climate factors as independent variables to construct the multiple regressions as given in Eq. (5):

where β₁–β₇ are the model parameters, t is the tehsil, and $ {\varepsilon }_{t} $ is the error term.

It is not necessary for a linear relationship to exist between crop growth response and climate factors in all conditions because climate variables feature both linear and nonlinear trends (Shi et al., 2013). We calculated the nonlinear effects by including quandary terms for each selected climate variable. To analyze the source of variation in cotton yield, we tested the assumption of autocorrelation at a 95% confidence interval (α = 0.05). An autoregressive model was fitted to those tehsils where autocorrelation existed. To obtain the optimal balance between the yield of the previous year and the predicted yield, we applied the Akaike Information Criteria (AIC) approach explained in (Lobell and Field, 2007). The statistical model was validated using P-values, and the model’s explanatory power was measured by a coefficient of determination, R². A low R² value indicates that the model did not accurately capture the observed crop yield response to climate factors (Lobell and Burke, 2010). We applied a two-tailed t-test to test the statistical significance of the regression model. A climate factor that demonstrated a significant relationship (α = 0.05) with the cotton yield was considered a dominant climate factor. For example, a P indicates that only precipitation had a significant impact on cotton yield, while PH indicates that both precipitation and humility significantly impact the cotton yield. The effects of changes in the spatial distribution pattern of cotton crops between tehsils on the response of district-level cotton yield can be estimated using the results of the regression models.

In most districts, there has been a significant change in the cotton area, with more anomalies exceeding one standard deviation (Fig. 2a). Specifically, the annual trend in the cotton-cultivated area increased before 2010, followed by a decrease in most districts, particularly those with high producers. During 2000–2005, most districts showed an increase in their cotton-cultivated area with an anomaly greater than one standard deviation. The pattern of cotton-cultivated areas between the tehsils shows different variations at the district level. Significant changes have been observed in all districts under investigation (Fig. 2b). Such changes in the pattern of cotton-cultivated areas correspond to higher SI values in the central districts and relatively lower SI values in the southern districts.

Figure 2.  Changes in the annual anomaly of cotton area (divided by standard deviation) compared to the mean cotton-cultivated area of the relative district (a) and stability index (SI) indicating the spatial distribution pattern of cotton-cultivated area at the district level (b) in Indus River Basin of Pakistan. The abbreviated district names are defined as Jhang (JHG), Toba Tak Singh (TTS), Sahiwal (SHL), Pakpatan (PAK), Khanewal (KWL), Vehari (VHR), Multan (MUL), Muzafarghar (MGH), Lodhran (LDR), Bahawalnagar (BWN), Bahawalpur (BWP), Rajanpur (RJP), Rahimyar Khan (RYK), Ghotki (GHT), Nawabshah (NWS), Khairpur (KRP), Mirpur khas (MPK), Sanghar (SNG), and Umar kot (UKT)

The spatial pattern of cotton yield variability and change trends from 1990 to 2019 were studied. Notably, the yearly variation in cotton yield at the tehsil level had experienced changes, significant at the 0.05 level over the last three decades (Fig. 3a). High variability of cotton yield was found in Jhang, Sahiwal, Nawabshah, and Umerkot, while low variability was observed in Bahawalpur, Bahawalnagar, and Rahimyar khan districts (Fig. 3b). The cotton yield variability had increased significantly in 43% of tehsils, mainly in the central part of the Indus River basin cotton belt, which accounted for 62% of total cotton production in the IRBP (Table 1). The results were consistent when smaller and larger time windows were used. Cotton yield variability showed a distinct spatial pattern, with one-third of the tehsils across the IRBP showing a decreasing trend.

Figure 3.  Cotton yield variability trends relative to the mean yield during 1990–2019 (a, b, c) (% / yr) and the spatial distribution pattern of percent cotton yield variability at tehsil-level (d) in the Indus River Basin of Pakistan. The star indicates those tehsils where the trend was statistically significant (α = 0.05)

The variability of the cotton yield at the tehsil level increased by up to 7.0% during the first 15 yr, while it increased by 9.5% during the last 15 yr. An increase in cotton yield variability was strongly correlated with the increasing trend in temperature, relative air humidity, and precipitation variability in the cotton-growing season, with correlation coefficients of up to 0.88, 0.79, and 0.84, respectively (Fig. 4). The time series of cotton yields and climate variability are also detrended, and the correlation analysis is repeated accordingly. We discovered that when their trend is removed, the strength of the relationship weakens. This means that recent trends in cotton yield variability are strongly linked to climate variability trends.

Figure 4.  Temporal variability of cotton yield and growing season climate variables for all the cotton-cultivated areas of the Indus River Basin of Pakistan. The x-axis is the 15-yr moving window’s central year. Correlation analysis is conducted between cotton yield variability and climate variables variability using the observed and detrended time series. Statistical significance is indicated at the 0.01 and 0.05 significance levels by ^** and ^*, respectively

Through piecewise regression analysis, we obtained the trend-turning year of the studied climate variables in the IRBP from 1990 to 2019. Our estimated trend turning point year was 2002, based on the precipitation. During the study period, the average precipitation decreased at the rate of –1.49 mm/yr, the rate increased after the TTP, and the observed average estimated precipitation trend was –0.64 mm/yr and 1.52 mm/yr before and after TTP, respectively (Fig. 5a). The temperature has increased at an average rate of 0.0256°C/yr in the IRBP during the cotton growing season from 1990 to 2019. The temperature trend became sharper after the TTP, as evidenced by the estimated trend before and after the TTP, of 0.0195 °C/yr and 0.0377°C/yr, respectively (Fig. 5b). Our results indicate the sharp increasing trend was most notable in the northern and central part of the study area. Relative air humidity also displayed an opposing trend before and after the TTP over the whole study period; the relative air humidity increased by 0.27 %/yr on average, and it showed a decreasing trend of –0.21 %/yr before 2002 and an increasing trend of 0.82 %/yr from 2002 to 2019 (Fig. 5c).

Figure 5.  Spatial distribution of the trend pattern in climate variables during 1990–2019 at tehsil level in the Indus River Basin of Pakistan. a. precipitation; b. temperature; c. relative humidity. 1990–2002, before the trend turning year; 2002–2019, after the trend turning year

Using a statistical model, we estimated the contribution of climate change to the observed trends in cotton yield variability. The variability trend in cotton yield due to studied climate variables is depicted in (Fig. 6a), which is derived from the overall trend from 1990 to 2019 to examine the importance of how climate compares to all other factors. Before the turning point (1990–2002), H (19.1%), TH (23.5%), and HP (11.8%) were the leading climate factors impacting the cotton yield at the tehsil level. While after the turning point, temperature and precipitation became important factors, as evidenced by the percentage by which T (10.6%), P (18.9%), TH (13.7%), and HP (14.5%) impacted the cotton yield. From 2002 to 2019, TP and TH became the main climate factors that controlled the cotton yield in the main cotton-cultivated districts, Khanewal, Vehari, and Lodhran. Autocorrelation was found in 13.5% of the districts for 1990–2019, 10.8% for 1990–2002, and 9.6% for 2002–2019 of the total number of tehsils, as indicated with black dots (Fig. 6a). Climate change had a higher impact on cotton yield variability in the central parts of the IRBP, where there had been an increasing variability trend in cotton yield. Climate variability influenced 10%–45% of the variability trends in cotton yield in the districts of the northern and the lower central parts of the Ghotki and Rahimyar khan districts, indicating that climate variability was not the main component controlling the variability in cotton yield in these areas. Cotton yield variability was also influenced by more than climate variables in the Multan, Bahawalpur, and Khairpur regions (Fig. 6b). Overall, temperature, relative humidity, and precipitation as a single climate factor, controlled cotton yield variability in 11% (15%), 6% (10%), and 7% (15%) of cotton growing tehsils (production), respectively, though climate variability was not the main factor in 16% (19%) of cotton growing tehsils (production). Similarly, the results from before TP (1990–2002) and after TP (2002–2019) are shown in Table 2.

Figure 6.  Spatial pattern of cotton yield variability controlled by climate variables at tehsil level (a) and the percentage of cotton yield variability explained by growing season climate variability as predicted by a statistical model (b) in the Indus River Basin of Pakistan. Statistically significant α = 0.05. The legends are explained in Table 2. Dots indicate the tehsils with autocorrelation in their cotton yield

To examine the effect of climate change on crop production, process-based simulation crop models and statistical models are commonly used to estimate the strength of the relationship between historical climate observations and yields. Crop model estimates frequently exhibited significant differences due to a lack of comprehensive spatially determined datasets on crop types, rotations, climate, land surface, and management (Van Ittersum et al., 2013; Asseng et al., 2015). We used statistical models to investigate the relationship between the trends of the main climate factors and cotton yield variability. The statistical model has limitations because of the independent factors selected for the regression analysis and the covariation between the selected climate factors (Leng et al., 2016). We removed the climate covariation from the yield-climate variability relationship and repeated the analysis. Similar results, albeit with a lower correlation, were obtained. We repeated the correlation analysis after detrending the cotton yield and climate variability time series. When the trends were removed, the strength of their relationships weakened, and some relationships lost their significance. This means that more recent variations in cotton yield are inextricably linked with variations in the climate. The results show that cotton yield variability is increasing at a rate of up to 0.17%/yr, which is much closer to the observations.

Multiple regression models with quandary terms were used to quantify the influence of climate variability on reported cotton yield trends. Overall, temperature and humidity have slightly increased during the study period while precipitation decreased in the Multan, Khanewal, Vehari, and Lodhran districts. These districts share almost 30% of the total cotton production in the study area. After the TTP (2002–2019), all three climate factors showed a sharp increasing trend, especially in the northern part of the in the Indus River Basin of Pakistan (Fig. 5). This shift in trend of climate factors significantly impacted the cotton yield, and the combined climate factors TP, TH, and HP became the main climate factors in controling the cotton yield in the tehsils that belong to the leading cotton producing areas. Before the TTP (1990–2002), climate factors showed an insignificant relationship with cotton yield in main cotton-producing areas (Fig. 6a). Our results indicate that, overall, 84% of the tehsils and 81% of the cotton production in the Indus River Basin of Pakistan were controlled by the main studied climate factors from 1990 to 2019. The proportion of the impact of climate variability in cotton yields for the study period is shown in (Fig. 6b). This indicates the importance of climate change in crop production, as climate is the main factor influencing cotton yields. The recent increase in climate variability is clearly a major cause of the observed trend in cotton yield variability across most cotton-cultivated areas. The fact that climate variability remains unexplained in 16% (1990–2019), 23% (1990–2002), and 17% (2002–2019) of the tehsils demonstrates the significance of variables left out of our research. Non-climatic factors, such as soil moisture, irrigation, conservation tillage (Karlen et al., 2013; Leng et al., 2014; Qin et al., 2015), fertilization, and multiple cropping (Seifert and Lobell, 2015; Leng et al., 2016) are not considered in this study. Indeed there is a need to better understand the non-climatic factors influencing the changing trend in cotton yield.

According to numerous studies, the likelihood of heatwave and drought occurrences has increased due to global warming (Hao et al., 2013; Leonard et al., 2014). Despite these findings, the effects of climate covariance on crop yields have received little attention. It is not easy to deduce the impact of a single climate factor without considering the effects of climate covariations (Leng et al., 2016). We conducted statistical analysis and removed such covariations to determine a single climatic factor’s variability relationship with cotton-yield variability and found significant results. Our model provides impartially accurate classifications of historical cotton yield variations, explaining more than 80% of the overall cotton yield variations for the IRBP. The remaining variability in cotton yields indicates that the suggested model cannot explain variability regarding processes such as crop infection, pollination, and dormancy, which could be key aspects of climate factors impacting crop yields. Further, we presented the dominant climate factors (T, H, and P) at the tehsil level across the Indus River Basin of Pakistan, and this pattern was found to be robust over the past three decades. Our findings could assist in developing effective crop management approaches to mitigate the negative impacts of climate change on cotton production. The computed sensitivity of crop yields of individual climate elements helps to assess crop model robustness in modeling the crop growth responses to a changing climate. This information is becoming more important as countries and regions are attempting to make fiber protection more resilient in the face of climate change.

We investigated the impact of temperature, relative air humidity, and precipitation on the variability of cotton yields at the tehsil level in the Indus River Basin of Pakistan from 1990 to 2019. We found that climate factors are the key components driving cotton yield variability during the cotton growing season. The statistically fitted models provided information about the change trend rate in cotton yield variability and the sensitivity of yield variability to the recently changing climate conditions. In addition to gaining a better understanding of the historical relationship between climate change and cotton yield variability, we could use these results to increase the flexibility of the Indus River Basin of Pakistan cotton production system. Climate change is likely to become more prevalent in areas where climate variability has historically accounted for most of the variability in cotton yields. As a result, some strategies for stabilizing cotton yields need to be devised to keep cotton yields stable in the future, as well as in light of such changes. So, the maps presented in this study can help management make better decisions about how to adapt to changes in the climate and how that affects cotton yield. This will help them become more flexible when dealing with climate change and determine the physical mechanism by which climate affects cotton yields. Further analysis of non-climatic factors, such as cropping methods, seed and soil quality, and information about farmers’ educational and financial conditions at the regional scale, could be more helpful for understanding the factors affecting the yield and variability of crops.