Seasonal Cropland Trends and Their Nexus with Agrometeorological Parameters in the Indus River Plain

: The ﬁne-scale insights of existing cropland trends and their nexus with agrometeorological parameters are of paramount importance in assessing future food security risks and analyzing adaptation options under climate change. This study has analyzed the seasonal cropland trends in the Indus River Plain (IRP), using multi-year remote sensing data. A combination of Sen’s slope estimator and Mann–Kendall test was used to quantify the existing cropland trends. A correlation analysis between enhanced vegetation index (EVI) and 9 agrometeorological parameters, derived from reanalysis and remote sensing data, was conducted to study the region’s cropland-climate nexus. The seasonal trend analysis revealed that more than 50% of cropland in IRP improved signiﬁcantly from the year 2003 to 2018. The lower reaches of the IRP had the highest fraction of cropland, showing a signiﬁcant decreasing trend during the study period. The nexus analysis showed a strong correlation of EVI with the evaporative stress index (ESI) during the water-stressed crop season. Simultaneously, it exhibited substantial nexus of EVI with actual evapotranspiration (AET) during high soil moisture crop season. Temperature and solar radiation had a negative linkage with EVI response. In contrast, a positive correlation of rainfall with EVI trends was spatially limited to the IRP’s upstream areas. The relative humidity had a spatially broad positive correlation with EVI compare to other direct climatic parameters. The study concluded that positive and sustainable growth in IRP croplands could be achieved through effective agriculture policies to address spatiotemporal AET anomalies.


Introduction
The cropland is a complex ecosystem that is subject to many biotic (diseases, pathogens, weeds, etc.) and abiotic (drought, salinity, waterlogging, etc.) stresses, as well as humaninduced land-use changes and environmental pollution. Moreover, climate change is directly (by temperature, precipitation, and carbon dioxide changes) or indirectly (by resource degradation and shifts in crop suitability) affecting cropland's productivity and is becoming an inevitable threat to its sustainability. The negative impacts of climate change on crop productivity can be eased by adopting efficient strategies and making the cropping system more resilient towards extreme events [1]. However, past studies have reported variable impacts of climate parameters on crop productivity subjected to land use patterns and geographic location [2][3][4]. Therefore, local-scale existing trends in cropland and its nexus with climate parameters are of extreme importance to layout efficient adaptation policies to limit climate change adversities on the cropping system. Regional cropland surveillance has been revolutionized since the first Earth observation satellite launch more than five decades ago. The space-borne remote sensing allowed the collection of accurate data to explore the different Earth's ecosystems on a large scale. The free access to the remote sensing data with moderate spatial resolutions provided the developing countries a way to integrate local-scale ground information into their regional policymaking for better sustainability of the resources.
In recent years, several large-scale studies have utilized remote sensing data to analyze the region's spatiotemporal vegetation dynamics and its nexus with different climatic parameters. For example, Zhang et al. [5] studied the vegetation dynamics and its driving forces in the upstream region of China's three main rivers from 1982 to 2012. The net primary productivity (NPP) and its trends were estimated and analyzed using the normalized difference vegetation index (NDVI). The study concluded that NPP in a large part of the area increased significantly, and solar radiation was a key climatic parameter contributing to this change. Similarly, Zhang et al. [6] observed a significant greening trend in global vegetation from 2001 to 2015 and attributed it to enhanced carbon sinkage in the land. Whereas Pan et al. [7], using multi-source remote sensing data, reported an increasing trend of vegetation browning hidden under overall vegetation greening in the northern hemisphere due to drought events. Two studies conducted in the Yangtze River basin concluded that vegetation trends are predominately influenced by human activities [8,9].
The cropland trend analysis may include various components of anthropogenic activities based on the region's characteristics to better understand the existing patterns and their implication in the future. For example, Wang et al. [10] recently studied the trends of croplands in China's loess plateau by analyzing the spatiotemporal patterns of field observed evapotranspiration and found a significant upward trend. The factorial analysis of variance between ET and leaf area index (LAI) pointed out that the enhanced cropping intensity was the main driver for the observed cropland trend. Many studies have also analyzed the effect of urban expansion on croplands around the big cities [11][12][13][14]. Most land-use change studies that analyzed the urban expansion on agricultural land have adopted an image classification strategy using multi-temporal images during the study period.
Apart from land-use change studies, the methods used for trend analysis of croplands are primarily based on different parametric and nonparametric statistical tests. The Mann-Kendall trend test is a nonparametric test that has been used in various regional and global studies to detect monotonic trends in vegetative areas [8,15]. The magnitude of trends can be estimated using simple linear regression slope or more robust Theil-Sen's slope estimator [16,17]. The Pearson correlation coefficient (r) with a significance level can be used to illustrate the relationship between variables and to identify the driving factors of the observed pattern [18,19]. However, the parameters required by the correlation analysis are often limited by their availability from the sources of statistics.
This study aims (1) to identify the existing seasonal trends in cropland of Indus River Plain (IRP) using vegetation index (VI) data from 2003 to 2018; (2) to correlate the temporal response of VI with various agrometeorological parameters to establish existing croplandclimate nexus in the region; (3) to analyze the dominant driving factor and underline causes. The results provide support to agricultural policymaking for improving cropland productivity and regional food security and better understand the potential impact of future climate change.

Study Area
The Indus river basin lies between latitude 24-37 • N and longitude 66-83 • E. It encompasses the territories of four countries, including China, Afghanistan, Pakistan, and India. The Indus River Plain (IRP) starts at the foothills of the upper Indus basin. It is a vast expanse of fertile soil covering about 619,000 km 2 from Himalayan Piedmont in the north to the Arabian Sea in the south. The IRP can be further divided into upper and lower IRP primarily based on river physiography. Five primary tributaries of the Indus river traverse through the upper IRP and converge into one mainstream before entering to lower IRP. Thal and Thar are two deserts that exist in the IRP (Figure 1). Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 16 Figure 1. The geographical location of the Indus river basin, its plains, and administrative boundaries, MODIS tiles coverage, and major geographical features.
About 90% of the vegetative area in IRP consists of croplands and forest cover in the north. The fraction of cropland between Pakistan, India, and Afghanistan is 63%, 36%, and 1%, respectively. The cropland in four provinces of Pakistan (Punjab, Sindh, Balochistan, and Khyber Pakhtunkhwa) and three states of India (Punjab, Haryana, and Rajasthan) is the focus of this study. The study area has two major cropping seasons. The Kharif season starts in April-May and lasts till October-November. The Rabi season starts in November-December and ends in April-May of the next year. The study period of this study is from November 2002 to October 2018. The geographical location, administrative units, and significant geographical features of the Indus river basin are shown in Figure 1.

Remote Sensing Data
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a sensor aboard the Terra and Aqua satellites. There are multiple standard data products of MODIS that can be used for various global to regional studies. Two data products from MODIS were used in this research. One is the enhanced vegetation index (EVI). The MODIS-EVI data from Terra (MOD13A1 (V6)) and Aqua (MYD13A1 (V6)) at 500 m spatial resolution was used to examine cropland trends and to evaluate cropland-climate nexus. By temporally interlacing the 16-day MOD13A1 and MYD13A1, an EVI data set with the 8-day temporal resolution was composed.
The second MODIS-based data product is the actual evapotranspiration (AET). The data product MOD16A2 (V6) with 500 m spatial and 8-day temporal resolution was used to extract AET data. The 8-day MYD16A2 (V6) was also used to fill the data gaps in the primary data product MOD16A2. The MODIS-AET product is based on the Penman-Monteith equation that employs daily meteorological data and vegetation biophysical properties like LAI (Leaf Area Index) and albedo derived from remote sensing [20]. The spatial About 90% of the vegetative area in IRP consists of croplands and forest cover in the north. The fraction of cropland between Pakistan, India, and Afghanistan is 63%, 36%, and 1%, respectively. The cropland in four provinces of Pakistan (Punjab, Sindh, Balochistan, and Khyber Pakhtunkhwa) and three states of India (Punjab, Haryana, and Rajasthan) is the focus of this study. The study area has two major cropping seasons. The Kharif season starts in April-May and lasts till October-November. The Rabi season starts in November-December and ends in April-May of the next year. The study period of this study is from November 2002 to October 2018. The geographical location, administrative units, and significant geographical features of the Indus river basin are shown in Figure 1.

Remote Sensing Data
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a sensor aboard the Terra and Aqua satellites. There are multiple standard data products of MODIS that can be used for various global to regional studies. Two data products from MODIS were used in this research. One is the enhanced vegetation index (EVI). The MODIS-EVI data from Terra (MOD13A1 (V6)) and Aqua (MYD13A1 (V6)) at 500 m spatial resolution was used to examine cropland trends and to evaluate cropland-climate nexus. By temporally interlacing the 16-day MOD13A1 and MYD13A1, an EVI data set with the 8-day temporal resolution was composed.
The second MODIS-based data product is the actual evapotranspiration (AET). The data product MOD16A2 (V6) with 500 m spatial and 8-day temporal resolution was used to extract AET data. The 8-day MYD16A2 (V6) was also used to fill the data gaps in the primary data product MOD16A2. The MODIS-AET product is based on the Penman-Monteith equation that employs daily meteorological data and vegetation biophysical properties like LAI (Leaf Area Index) and albedo derived from remote sensing [20]. The spatial extent and description of tiles of the MODIS products are shown in Figure 1. The preprocessing has Remote Sens. 2021, 13, 41 4 of 15 applied to the remote sensing data sets including format conversion, mosaic, subsetting, and reprojection.

Agrometeorological Data
Due to the sparse spread of meteorological stations in the region and limited data availability, climate parameters from a newly available reanalysis dataset (AgERA5) were used for cropland-climate nexus. The AgERA5 provides gridded climate parameters at a spatial resolution of 10 km from 1979 to 2018. The AgERA5 dataset is based on fifth-generation climate reanalysis data (ERA5) produced by independent intergovernmental organization ECMWF (European Centre for Medium-Range Weather Forecasts). The detailed characteristics and general configuration of ERA5 global reanalysis can be found in Hersbach et al. [21]. Several recent studies have indicated better performance of ERA5-derived climate parameters when compared with other climate reanalysis datasets, satellite-derived, and ground-observed climate parameters [22][23][24][25].
Five climatic parameters, including total daily precipitation (PCP), mean daily temperature (TMP), average daily relative humidity (RH), average daily wind speed (WND), and mean daily solar radiation (SLR) were derived from AgERA5. The Hargreaves method [26] was used to estimate historic potential evapotranspiration (PET) using TMP and SLR of AgERA5. Additionally, two indices were also calculated to indicate water deficiency across the study area. The Aridity Index (AI) was calculated using the ratio of PCP and PET to highlight the natural deficiency of water at a given location, as defined by Barrow [27]. The Evaporative Stress Index (ESI) was derived as ESI = AET/PET and it is structurally similar to the crop coefficient (K c ) defined by Allen et al. [28].

Overall Methodology
The stepwise procedure of this study is shown in Figure 2. The analytical procedure has two main components: (a) identifying existing croplands trends in the study area, and (b) analyzing the linkage between cropland response to agrometeorological parameters to understand cropland-climate nexus. For the first part, the seasonal stacks of EVI for the Rabi and Kharif seasons were prepared. The Mann-Kendall (MK) test was used to detect EVI's monotonic trend in both seasons, whereas Sen's slope estimator was used to quantify the seasonal trend. The results of the MK-test and Sen's slope were combined to identify different classes of EVI trends.
For the second part, a Pearson correlation analysis between EVI and 9 agrometeorological parameters, including PCP, TMP, SLR, WND, RH, AET, PET, ESI, and AI, was conducted at a monthly level for both seasons. The alternate hypothesis of correlation greater than zero at the significance level of 0.05 was used in estimating the p-value. The spatial analysis of the correlation level between EVI and the agrometeorological parameters was conducted for each season to identify the key parameters responsible for the observed seasonal EVI trend. The spatial resolution of remote sensing data, EVI and AET, was down-scaled to 10 km by the pixel aggregation method, using average value, to match the climatic data spatial resolution. For the second part, a Pearson correlation analysis between EVI and 9 agrometeorological parameters, including PCP, TMP, SLR, WND, RH, AET, PET, ESI, and AI, was conducted at a monthly level for both seasons. The alternate hypothesis of correlation greater than zero at the significance level of 0.05 was used in estimating the p-value. The spatial analysis of the correlation level between EVI and the agrometeorological parameters was conducted for each season to identify the key parameters responsible for the observed seasonal EVI trend. The spatial resolution of remote sensing data, EVI and AET, was down-scaled to 10 km by the pixel aggregation method, using average value, to match the climatic data spatial resolution.

Trend Analysis
The Mann-Kendall (MK) is a nonparametric test used to indicate upward or downward monotonic trends shown by a given data set during a period of time. The null hypothesis in the MK test assumes no monotonic trend [29,30]. At a specified significance level (α), the Z-score value rejects the null hypothesis and accepts the alternative hypothesis of the monotonic trend presence. In general, α = 0.05 indicates a statistical significance. The MK test calculates the statistic S as:

Trend Analysis
The Mann-Kendall (MK) is a nonparametric test used to indicate upward or downward monotonic trends shown by a given data set during a period of time. The null hypothesis in the MK test assumes no monotonic trend [29,30]. At a specified significance level (α), the Z-score value rejects the null hypothesis and accepts the alternative hypothesis of the monotonic trend presence. In general, α = 0.05 indicates a statistical significance. The MK test calculates the statistic S as: where sgn x j − x i is a sign function as Remote Sens. 2021, 13, 41 The variance is computed as where n is the extent of the time series data, x i and x j are the data values in the time series i and j(j > i), respectively, g is the number of tied groups and t p is a tie group number that has equal values in pth group. If there are no tie groups, tie groups' summation process can be ignored [31]. The Z statistic of the MK test is defined as The value of Z ranges from (−∞, +∞). Positive and negative values of Z value indicate increasing and decreasing trends, respectively. At significance level α = 0.05 value of |Z| = ±1.960. For |Z| > 1.960, the null hypothesis is rejected, and a significant trend exists in the time series data.
The Theil-Sen estimator is another robust nonparametric technique used to quantify the slope of the observed trend. Sen's slope is a robust parameter because the outliers in the data least affect the slope value. Sen's slope estimator makes a pair-wise combination of total n(n − 1)/2 data pairs and calculates slopes for each of the pairs [32]. The median of all the slopes reflects the change in data per unit time. Mathematically Sen's slope is expressed as: where x i and x j are the data values at times i and j, respectively. The sign of β represents the data trend, while its value indicates the change per unit time. In a given period, if β > 0, the trend of given data is increasing, while β < 0, a decreasing trend is shown. The combination of Sen's slope values (β) and MK test, Z-score values at the confidence interval of 95% were used to classify the cropland trends in the region. The cropland area with β > 0 and Z > 1.96 (p < 0.05) and area β < 0 and Z < −1.96 (p < 0.05) was classified as the significantly improved and significantly decreased trend, respectively. The area with, β > 0 or β < 0 with Z values between −1.96 to 1.96 (p ≥ 0.05) was termed as having insignificantly improved or decreased cropland trends.

Rabi Season
Based on the classification scheme above, in the study area, 58% of the total cropland during the Rabi season has significantly improved over the past 16 years. In comparison, only 2% of the cropland suffered from significant degradation. A large portion (40%) of the cropland was under insignificant trends during Rabi season. The spatial variation of Z-score, value of Theil-Sen slope, and cropland trend classification during Rabi season is shown in Figure 3 trend. The fraction of the cropland area in each category of trend class within a province is shown in Table 1.
croplands in Sindh province have also demonstrated an insignificant decreasing trend. The croplands at the outskirts of IRP have improved significantly, especially in Thal and Thar regions' drylands. In Punjab's arid regions (Pak) and KP, the croplands have shown mostly positive trends with some large portion showing statistically insignificant improvement. More than 70% of the cropland in Indian Punjab and Haryana has shown a significant positive trend. Whereas, 25% of the cropland in Rajasthan showed an insignificant positive trend. The fraction of the cropland area in each category of trend class within a province is shown in Table 1.  The cropland of IRP during the Kharif season has shown spatially different trends in various parts of regions as compared to the Rabi season. Nearly 50% of the total cropland has shown a significantly increasing trend, and 2% of cropland shows a negative (p < 0.05) trend. About 47% of the croplands are under insignificant trends during this season. The spatial variation of z-score, value of Theil-Sen slope, and cropland trend classes during

Kharif Season
The cropland of IRP during the Kharif season has shown spatially different trends in various parts of regions as compared to the Rabi season. Nearly 50% of the total cropland has shown a significantly increasing trend, and 2% of cropland shows a negative (p < 0.05) trend. About 47% of the croplands are under insignificant trends during this season. The spatial variation of z-score, value of Theil-Sen slope, and cropland trend classes during the Kharif season is shown in Figure 4. The analysis within each province revealed that Balochistan and Sindh have the lowest fraction of the cropland showing a significantly increasing trend. Moreover, both provinces have the highest fraction of croplands with a significantly decreasing trend compared to other provinces/states. The fraction of the cropland area in each category of Kharif season's trend class at the province/state level is presented in Table 2. the Kharif season is shown in Figure 4. The analysis within each province revealed that Balochistan and Sindh have the lowest fraction of the cropland showing a significantly increasing trend. Moreover, both provinces have the highest fraction of croplands with a significantly decreasing trend compared to other provinces/states. The fraction of the cropland area in each category of Kharif season's trend class at the province/state level is presented in Table 2.

Rabi Season
The pair-wise correlation matrix between EVI and 9 agrometeorological parameters is shown in Figure 5. The red color in Figure 5 indicates that Pearson's correlation (r) value is positive and significant (p ≤ 0.05). The positive and significant nexus of EVI with AET and ESI was among the highest (r = 0.78 and r = 0.70, respectively). From direct climate parameters, EVI correlated with RH and PCP (r = 0.45 and r = 0.18, respectively, at p ≤ 0.05). The nexus of EVI with other agrometeorological parameters and their mutual correlation can be seen in Figure 5. To further analyze the EVI and climate interactions, pixel-

Rabi Season
The pair-wise correlation matrix between EVI and 9 agrometeorological parameters is shown in Figure 5. The red color in Figure 5 indicates that Pearson's correlation (r) value is positive and significant (p ≤ 0.05). The positive and significant nexus of EVI with AET and ESI was among the highest (r = 0.78 and r = 0.70, respectively). From direct climate parameters, EVI correlated with RH and PCP (r = 0.45 and r = 0.18, respectively, at p ≤ 0.05). The nexus of EVI with other agrometeorological parameters and their mutual correlation can be seen in Figure 5. To further analyze the EVI and climate interactions, pixel-level correlations were divided into five categories using r and p. The description of each category is presented in Table 3. Each category's spatial analysis was done at the basin and sub-basin (state/province) level to identify the primary driver of seasonal cropland trend. level correlations were divided into five categories using r and p. The description of each category is presented in Table 3. Each category's spatial analysis was done at the basin and sub-basin (state/province) level to identify the primary driver of seasonal cropland trend.  The spatial analysis of the Rabi season's correlation at basin level revealed that 63% of the IRP croplands have a strong correlation between EVI and ESI. The strong nexus of EVI with AET is at 53% of the total IRP's cropland. No cropland area demonstrated strong nexus of EVI with PCP and RH. The fractions of cropland area at the basin and sub-basin level with defined nexus classes of EVI and agrometeorological parameters are shown in Figure 6. The nexus class analysis at the province/state level found that the highest cropland fraction of Punjab (Ind, Pak), Haryana, Sindh, and Rajasthan had strong nexus between EVI and ESI compared to AET. In comparison, Balochistan and KP's cropland had the highest fraction of non-significant nexus of EVI with all agrometeorological parameters. Further details of cropland fractions at defined nexus levels between EVI and various climatic parameters within each province/state can be seen in Figure 6. The pixellevel spatial variation of Rabi season's cropland-climate nexus is shown in Figure 7. Figure  7 depicted that low altitude croplands of Balochistan and KP had better nexus with ESI and AET. However, these two provinces' high-altitude croplands were better correlated  The spatial analysis of the Rabi season's correlation at basin level revealed that 63% of the IRP croplands have a strong correlation between EVI and ESI. The strong nexus of EVI with AET is at 53% of the total IRP's cropland. No cropland area demonstrated strong nexus of EVI with PCP and RH. The fractions of cropland area at the basin and sub-basin level with defined nexus classes of EVI and agrometeorological parameters are shown in Figure 6. The nexus class analysis at the province/state level found that the highest cropland fraction of Punjab (Ind, Pak), Haryana, Sindh, and Rajasthan had strong nexus between EVI and ESI compared to AET. In comparison, Balochistan and KP's cropland had the highest fraction of non-significant nexus of EVI with all agrometeorological parameters. Further details of cropland fractions at defined nexus levels between EVI and various climatic parameters within each province/state can be seen in Figure 6. The pixel-level spatial variation of Rabi season's cropland-climate nexus is shown in Figure 7. Figure 7 depicted that low altitude croplands of Balochistan and KP had better nexus with ESI and AET. However, these two provinces' high-altitude croplands were better correlated with TMP, PET, and SLR. The correlation between EVI and TMP was not shown in Figure 7 due to the nexus between TMP and PET. with TMP, PET, and SLR. The correlation between EVI and TMP was not shown in Figure  7 due to the nexus between TMP and PET.   with TMP, PET, and SLR. The correlation between EVI and TMP was not shown in Figure  7 due to the nexus between TMP and PET.

Kharif Season
For the Kharif season, the pair-wise correlation matrix is shown in Figure 8. Like Rabi season, the EVI is strongly correlated with AET and ESI (r = 0.75 and r = 0.71, respectively, at p ≤ 0.05). Another positive and significant correlation of Kharif season's EVI is with RH (r = 0.53). The correlation of AI and PCP with EVI is r = 0.42 and r = 0.41, respectively ( Figure 8). The spatial analysis of nexus categories at basin level discovered that area fraction under strong nexus class was highest between EVI and RH (37%). The area fraction of total cropland under strong nexus class for AET and ESI was 36% and 28%, respectively. However, the cumulative fraction of strong and relatively strong nexus for Kharif season' cropland was between EVI and AET ( Figure 9).

Kharif Season
For the Kharif season, the pair-wise correlation matrix is shown in Figure 8. Like Rabi season, the EVI is strongly correlated with AET and ESI (r = 0.75 and r = 0.71, respectively, at p ≤ 0.05). Another positive and significant correlation of Kharif season's EVI is with RH (r = 0.53). The correlation of AI and PCP with EVI is r = 0.42 and r = 0.41, respectively ( Figure 8). The spatial analysis of nexus categories at basin level discovered that area fraction under strong nexus class was highest between EVI and RH (37%). The area fraction of total cropland under strong nexus class for AET and ESI was 36% and 28%, respectively. However, the cumulative fraction of strong and relatively strong nexus for Kharif season' cropland was between EVI and AET ( Figure 9).

Kharif Season
For the Kharif season, the pair-wise correlation matrix is shown in Figure 8. Like Rabi season, the EVI is strongly correlated with AET and ESI (r = 0.75 and r = 0.71, respectively, at p ≤ 0.05). Another positive and significant correlation of Kharif season's EVI is with RH (r = 0.53). The correlation of AI and PCP with EVI is r = 0.42 and r = 0.41, respectively ( Figure 8). The spatial analysis of nexus categories at basin level discovered that area fraction under strong nexus class was highest between EVI and RH (37%). The area fraction of total cropland under strong nexus class for AET and ESI was 36% and 28%, respectively. However, the cumulative fraction of strong and relatively strong nexus for Kharif season' cropland was between EVI and AET ( Figure 9).  The fraction of cropland-climate nexus class at basin and sub-basin level for the Kharif season is shown in Figure 9. The spatial analysis at the province/state level revealed that the strong nexus class of EVI with RH was limited to Punjab (Pak and Ind) and Haryana regions. For the rest of the sub-basins, the fraction of strong and relatively strong EVI nexus existed between AET and ESI. Overall, AI compared to PCP had better nexus with Kharif season's EVI at basin and sub-basin level. The positive nexus of IRP's seasonal EVI was negligible with agrometeorological parameters such as TMP, SLR, WND, and PET. The spatial variation of the Kharif season's cropland-climate nexus is shown in Figure 10. The spatial nexus of TMP and PET with EVI was similar; therefore, only EVI interactions with PET are shown in Figure 10.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 16 Figure 9. The fraction of cropland-climate nexus class at basin and sub-basin level for Kharif season. The class fraction is the seasonal correlation between EVI and 9 agrometeorological parameters as indicated on the x-axis.
The fraction of cropland-climate nexus class at basin and sub-basin level for the Kharif season is shown in Figure 9. The spatial analysis at the province/state level revealed that the strong nexus class of EVI with RH was limited to Punjab (Pak and Ind) and Haryana regions. For the rest of the sub-basins, the fraction of strong and relatively strong EVI nexus existed between AET and ESI. Overall, AI compared to PCP had better nexus with Kharif season's EVI at basin and sub-basin level. The positive nexus of IRP's seasonal EVI was negligible with agrometeorological parameters such as TMP, SLR, WND, and PET. The spatial variation of the Kharif season's cropland-climate nexus is shown in Figure 10. The spatial nexus of TMP and PET with EVI was similar; therefore, only EVI interactions with PET are shown in Figure 10.

Discussion
The seasonal trend analysis based on EVI indicated a significantly positive trend in more than 50% of the IRP croplands since 2003. This finding of cropland greening is in line with previous regional and global studies. For example, Sarmah et al. [33] analyzed seasonal vegetation trends in South Asia from 2001 to 2013 using VI data from three remote sensing sources. The study reported a significant greening trend in croplands of IRP along with other cropland dominated semi-arid regions. Similarly, Gao et al. [34] analyzed global agricultural trends using remote sensing-derived LAI datasets. The study concluded that agricultural greening trends in developing regions were higher than developed regions, mainly due to cropland expansion and cropping season extension. In the current study, the positive trend of EVI at the peripheries of IRP indicates cropland expansion. Whereas, from the visual analysis of Figures 3 and 4, a negative trend in croplands around the major cities, primarily due to urban sprawl, can be noticed. A large fraction of IRP's cropland based on time-series MODIS-EVI demonstrated insignificant

Discussion
The seasonal trend analysis based on EVI indicated a significantly positive trend in more than 50% of the IRP croplands since 2003. This finding of cropland greening is in line with previous regional and global studies. For example, Sarmah et al. [33] analyzed seasonal vegetation trends in South Asia from 2001 to 2013 using VI data from three remote sensing sources. The study reported a significant greening trend in croplands of IRP along with other cropland dominated semi-arid regions. Similarly, Gao et al. [34] analyzed global agricultural trends using remote sensing-derived LAI datasets. The study concluded that agricultural greening trends in developing regions were higher than developed regions, mainly due to cropland expansion and cropping season extension. In the current study, the positive trend of EVI at the peripheries of IRP indicates cropland expansion. Whereas, from the visual analysis of Figures 3 and 4, a negative trend in croplands around the major cities, primarily due to urban sprawl, can be noticed. A large fraction of IRP's cropland based on time-series MODIS-EVI demonstrated insignificant (increasing and decreasing) trends during both crop seasons. These insignificant trends indicate irregular signals of cropland's temporal phenology. The prime reason for the signal inconsistency is coarser pixel data of MODIS compared to single farm size and cropland diversity [35].
The correlation analysis of 9 agrometeorological parameters revealed that AET and ESI had throughout strong nexus with EVI in both cropping seasons. The cropland trend and nexus results imply that AET and ESI in IRP's cropland have also improved significantly during the study period. The AET improvement in regional cropland fall in line with past global-scale research of Javadian et al. [36]. The study identified a significant increase of AET in regional croplands compared to no significant increase in natural vegetation worldwide. Similarly, Anderson et al. [37] studied the correlation of cropland productivity with satellite-driven LAI, PCP, and ESI anomalies. The study concluded that during the flash drought events, ESI was better able to capture soil moisture deficiency and its impact on grain yield. The IRP is a water stress region, especially during the Rabi season when it receives limited rainfall and has minimum flow in rivers. The spatial analysis of Rabi season's cropland-climate nexus categories at the sub-basin level revealed that ESI strongly correlates with EVI in major crop-producing regions. This means that during the water-limited cropping season, the ESI, which reflects the relative stress in AET, is a better indicator of cropland's phenology trends.
Among the direct climate parameters, RH had the highest positive nexus with EVI temporal response. The IRP's areas with significantly positive EVI trends suggest that those regions' RH level has also increased. This rise in RH is linked to intensive irrigation practices in the Indus river basin, as reported by Krakauer et al. [38]. The spatiotemporal nexus of EVI had a relatively strong linkage with climate dryness and wetness reflected as AI compared to direct PCP rates. The findings of this research demonstrated that AI in the north-eastern side of IRP has improved during the study period, especially during the Kharif season. Kimura [39] studied the satellite-driven global AI trends from 2001 to 2013 and reported a transition towards wetness in various parts of large river basins globally, including irrigated croplands of the Indus basin.
The agrometeorological parameters that include TMP, SLR, WND, and PET had a significantly negative nexus with EVI seasonal trend in the IRP. However, their mutual correlation was significantly positive (Figures 5 and 8). A future increase in any of these agrometeorological parameters resulting from global warming would harm the regional cropland's EVI response. For example, a rise in PET would increase the crop irrigation requirement to maintain optimal AET rates and ensure stable cropland productivity [40]. The more exploitation of water resources would further degrade the basin's ecosystem services and aggravate the existing water disputes at local and regional scales. In addition, future rise in TMP would lead to a rapid depletion of snow cover in the upper Indus basin and would limit the primary source of surface irrigation in the IRP [41]. Further investigations of regional snowfall trends are also required to study the linkage of snow melting and sustainability of IRP's cropland trends, keeping in view the strong relationship between EVI, AET, and ESI established in this study.
To sum up, the EVI trends in IRP is primarily affected by seasonal AET rates and relative stress in AET. Among other farm management practices, irrigation is exceptionally crucial to improve AET rates. The IRP is a heavily irrigated cropland with spatiotemporally skewed rainfall patterns. Regional crop water demands are met by conjunctive irrigation of surface and groundwater [42]. The surface water highly depends on natural conditions such as precipitation and snow melting. Similarly, the improved infrastructure for better water resource management is linked to the region's economic conditions. The positive EVI trends in the Thal and Thar's drylands indicate an increase in soil moisture, which is more related to the accessibility to a water supply through canal networks and adaptation of rainwater conservation techniques than an increase in direct rainfall.

Conclusions
This study has demonstrated the effective use of satellite remote sensing and fine resolution climate reanalysis data to evaluate the regional cropland trends and identify its agrometeorological drivers. This study concluded that most parts of IRP's cropland have significantly improved during the study period based on MODIS-EVI trends. The nexus analysis revealed that ESI had strong nexus with EVI during water limiting crop season, whereas for water excess crop season, AET was strongly correlated with EVI response. This research reiterates the importance of evapotranspiration and its stresses in controlling the cropland trends, especially in intensively irrigated regions with an arid and semi-arid climate. Knowledge of the existing cropland-climate nexus can improve the interpretation of future climate change impacts on regional cropland. This approach can be applied to other data-scarce regions to plan robust measures against changing climate and ensure sustainable productivity of regional cropland. Institutional Review Board Statement: Not applicable for this study.
Informed Consent Statement: Not applicable for this study.

Data Availability Statement:
The input data used in this research can be accessed freely from online sources. Remote sensing data used in this study can be downloaded via various web portals of NASA, e.g., https://earthdata.nasa.gov/. The climate reanalysis is available at climate data store website of European Centre for Medium-Range Weather Forecasts (https://cds.climate.copernicus.eu/). All the derived data is included in the manuscript. The trend and correlation analysis was done in MATLAB interface and scripts can be shared upon request.

Conflicts of Interest:
The authors declare no conflict of interest.