Irrigated cropland expansion exacerbates the urban moist heat stress in northern India

Agricultural irrigation has significantly reshaped the land surface energy and water balance. Previous studies have well investigated its cooling effect on air temperature (T air). However, its effect on increasing air humidity which can intensify the humid heat was often overlooked, particularly over urban areas with high population density, high T air, and consequently greater exposure to moist heat stress. In this study, using state-of-the-art reanalysis data at a high spatial resolution (∼9 km), we evaluated how changes in area equipped for irrigation (AEI) around a city affect urban moist heat stress (UMHS) in more than 1000 cities in China and India. In addition to T air, wet-bulb temperature (T WB) and wet-bulb globe temperature (T WBG), which consider humidity and are closer to the perceived temperature, were assessed. We found that although AEI expansion lowers urban T air, irrigation increases T WB and T WBG due to increased air humidity, thereby exacerbating the UMHS. This ‘warming’ effect of irrigation is more evident in northern India where AEI has expanded significantly in recent decades, and is prominent in the pre-monsoon and post-monsoon seasons, when precipitation and air humidity are low. However, this effect is not evident in China due to the lower intensity of AEI expansion and differing climatic conditions. Overall, this study highlights the side-effect of irrigation on regional climate, providing crucial information for better understanding urban heat stress and for future city planning.


Introduction
Irrigation is one of the most important anthropogenic land management activities, with strong biogeophysical effects on the land surface and regional climate (DeAngelis et al 2010, Yang et al 2020a). Since the 1950s, to meet the food demands of the rapidly increasing population, the global area equipped for irrigation (AEI) has increased substantially, from 111 Mha in 1950 to 306 Mha in 2005 (Siebert et al 2015). In particular, China and India show the largest AEI increases among all countries due to their huge populations (figure S1 available online at stacks.iop.org/ERL/17/054013/mmedia).
Irrigation increases soil moisture and adds additional atmospheric water vapor through evapotranspiration, which reduces the Bowen ratio, lowers surface temperature, increases air humidity, and may also influence regional precipitation. The cooling effects of irrigation on near-surface air temperature (T air ) and land surface temperature have been investigated in numerous studies in recent years, using both meteorological observations (Bonfils and Lobell 2007, Yang et al 2020a, 2020b and climate model simulations (Lobell et al 2008, Cook et al 2010, Kang and Eltahir 2019, Thiery et al 2020. Some studies have noted that irrigation can reduce human exposure to heat stress by alleviating extreme T air (Lobell et al 2008, Thiery et al 2020. However, in recent years, it has been increasingly recognized that the assessment of heat stress should consider not only T air , but also humidity, as high air humidity prevents humans and livestock from dissipating heat through evaporation of sweat (Pal and Eltahir 2015, Im et al 2017, Raymond et al 2020. For example, by analyzing 783 cases of heat-related human mortality in 36 countries, Mora et al (2017) investigated the climatic conditions associated with human death and found that T air and relative humidity (RH%) are the most important two factors that determine the lethality during heat events. Thus, moist heat stress is closely associated with human health, safety, and productivity risk characterizations (Kjellstrom et al 2016, Sherwood 2018. Wet-bulb temperature (T WB ) and wet-bulb globe temperature (T WBG ) are two widely used temperature indices that quantify the humid heat in an environment. T WB measures the joint effect of T air and air humidity on the perceived temperature of the human body, while T WBG also considers the influences of wind speed and solar radiation (equations provided in text S1). The values of these two indices are directly correlated with human health and are used to measure the intensity of heat stress. For example, a T WB of 35 • C will cause the cooling mechanism of the human body out of function and result in death within a few hours (Sherwood andHuber 2010, Hanna andTait 2015). Even a much lower T WB value of 28 • C caused deaths during the 2003 European and 2010 Russian heatwaves (Raymond et al 2020). Besides, different T WBG values are also used as thresholds for limiting outdoor military, occupational, and athletic activities (Japan Sports Association 2013).
In recent years, studies have begun to reappraise the impact of irrigation on heat stress through effects on heat and humidity. Kang and Eltahir (2018) and Krakauer et al (2020) conducted regional and global climate simulations respectively, with and without irrigation, and found that irrigation tends to increase extreme T WB as well as the intensity and frequency of moist heatwaves. Mishra et al (2020) investigated the effects of irrigation on humidity and extreme moist heat stress in India using satellite data and model simulations. They demonstrated that intensive irrigation in India increased air humidity and raised the moist heat stress metrics, indicating impacts on about 37-46 million people in South Asia. However, that study did not explicitly investigate urban areas or the influence of rapid expansion of irrigated cropland.
The effects of irrigation on moist heat stress have not been fully investigated to date, and several key concerns remain unsolved. For example, the few relevant studies have relied on numerical experiments and simply compared simulations with and without irrigation. Even for the state-of-the-art climate models, large biases still exist in these simulations (Guo et al 2019, Chen et al 2021. Moreover, due to inadequate representation of basic physical processes as well as temporal and spatial discretization, the responses of different models to irrigation may include large uncertainties (Cook et al 2014, Fyfe et al 2021. Therefore, investigating and verifying the results through examination of in-situ and reanalysis data is of great value. When considering heat and humidity effects, great attention should be paid to urban areas due to the urban heat island effect, high population density, and consequently high risk of moist heat stress (Daniel et al 2018, Katzfey et al 2020. However, no studies have specifically investigated the effects of irrigation on urban moist heat stress (UMHS). Therefore, using the latest reanalysis and in-situ datasets and focusing on urban areas, this study evaluates the effects of irrigated cropland expansion on UMHS for the first time. For more than 1000 cities in China and India that experienced rapid expansions of irrigated cropland and increases in urban population, we investigated the relationship of changes of AEI around cities with UMHS from 1980 to 2005. Moreover, the spatial diversity and the seasonal patterns of irrigation's effect on UMHS are investigated.

Data
We used four geophysical, two socio-economic, and two meteorological datasets to calculate AEI expansion around cities and measure its influences on the urban climate. Information about these datasets is provided in table 1. The climate classifications of all cities are shown in figure S2, which is determined by their coordinates (i.e. latitude and longitude) in Kottek et al (2006), including 14 climate types in China and 10 climate types in India.
State-of-the-art global reanalysis meteorological data were obtained from ERA5-Land (Munoz-Sabater et al 2021), which provides the evolution of land variables with very high temporal (hourly) and spatial (0.1˚, ∼9 km) resolution for the period of 1981-2020. T air , dewpoint temperature, downward solar radiation, wind speed, air pressure, and precipitation are used in this study. From the original meteorological variables from ERA5-Land, specific humidity (Q) and relative humidity (RH%) are calculated (text S2). In addition, to validate the reliability of ERA5-Land data, in-situ observations of T air , Q, RH%, and T WB from the Hadley Centre Integrated Surface Database (HadISD) (Dunn 2019) are used (text S3 and figures S3-S6).

Calculation of change in the area equipped for irrigation around each city
To calculate the change in AEI around each city, an impact area is defined for each city. Specifically, a square area around each city (referred to as a city square) is first determined (text S4 and figures S7 and S8), then the AEI proportion (AEI/total area of city square × 100 [%]) and its change from 1980 to 2005 are calculated for each city square. Due to the  (2019) differing sizes of cities, city square size varies among cities. In general, the size of a city square is approximately ten times the urban area of the city in 2005. In addition, city squares with large water surface areas are excluded from analysis, as water surfaces can have a strong influence on regional air humidity (Condie and Webster 1997). This algorithm obtains the AEI proportion change in each city square for 1036 cities (647 in China and 389 in India).

Moist heat stress indices
Two indices, wet-bulb temperature (T WB ) and wetbulb globe temperature (T WBG ), are used to estimate UMHS. T WB assumes that the skin is completely wet and unclothed, and can be determined as a function of T air and RH% (Stull 2011). T WBG is a more comprehensive index that is the weighted sum of the T air , T WB, and the globe temperature T G . T G measures the radiation factor (e.g. solar radiation) in human heat loss and is determined by radiant heat, T air , and wind speed. T WBG applies to realistic conditions of hard exertion when some skin is wet and exposed (Sherwood 2018). The calculation of T WBG in this study is based on Liljegren et al (2008). More detailed information on T WB and T WBG and their calculation are provided in text S1. The UMHS is defined as the intensity and duration of high values of T WB and T WBG over urban areas. We investigated the monthly-average daily maximum T WB and T WBG in each grid that includes a city coordinate at the resolution of ERA5-Land (∼9 × 9 km) to quantify UMHS. The changes of monthly-average daily maximum T air , T WB, and T WBG in the target area are calculated between 1981-1990 and 2001-2010 (the latter period minus the former, hereinafter abbreviated ∆T air , ∆T WB , and ∆T WBG , respectively). In addition to ∆T WB and ∆T WBG , changes in the number of hot days are also evaluated. The number of hot days is defined as the number of days per month on which the daily maximum temperature index values exceed a given threshold: 35 • C, 25 • C, and 31 • C for T air , T WB , and T WBG , respectively. These thresholds are based on previous research (Sherwood and Huber 2010, Japan Sports Association 2013) and represent the warning level that may cause heavy heat stress to the human body.

Change in the area equipped for irrigation around each city
The results show that cities in northwestern India and the North China Plain (NCP) have the highest AEI proportions in 1980, with AEI proportions of some cities being greater than 80% in northwestern India and 60% in the NCP (figure 1(a)). The national average AEI proportions of cities in China and India were 21.5% and 25.3% in 1980, respectively, with India's average value being slightly higher than China's (figure 1(c)).
In terms of changes in AEI proportion between 1980 and 2005 (hereinafter ∆AEI), northern India (mainly the Ganges River basin, hereinafter abbreviated NI) has gone through significant AEI expansion over those 25 years, with ∆AEI in many cities greater than 30% and a regional average ∆AEI of 14.1% (red box 1 in figure 1(b)). In China, apparent AEI expansion also occurred in the NCP, Northeastern China Plain (NECP) (red boxes 2 and 3 in figure 1(b), respectively), some cities in northwestern China (Xinjiang Province), and Inner Mongolia, with increases greater than 20% in a few cities. Regional average ∆AEI in the NCP and NECP are 6.8% and 8.2%, respectively. In general, ∆AEI has been much larger in India than in China, with national average values of 7.8% and 3.1%, respectively. Moreover, apparent AEI decreases are observed for cities in northwestern India, with the highest AEI proportion occurring in 1980. The maximum AEI proportion decrease of cities in this region is 22.0%.

Change in urban moist heat stress
We investigated ∆T air , ∆T WB , ∆T WBG (figures 2(a)-(c) and (g)-(i)) and changes in the number of hot days based on these temperature indices (figures 2(d)-(f) and (j)-(l)). The results for the hottest month in China (July) and India (May) are shown. The urban maximum T air in most Chinese cities shows an apparent increase in July (figure 2(a)), and the national average value for all Chinese cities is around 0.5 • C (figure 2(g)). The urban ∆T air values are slightly lower in the NCP and NECP than in other areas (figure 2(g)). In India, however, the ∆T air in May differs among cities located in different regions.
The urban maximum T air shows an apparent increase in cities in northwestern India, where an apparent decrease in AEI is observed ( figure 1(b)), as well as in the central region near the east coast. By contrast, apparent decreases in urban maximum T air are observed for cities in NI, which have experienced significant AEI expansion, as well as in southern India. In general, the increase of maximum urban T air is not as intense in India as in China, and the national average value for India is around 0 • C. In terms of the number of hot days based on T air , similar changes are observed in regions with significant AEI changes (i.e. NI, NCP, and NECP) as well as in the other regions of each country, except that the increase in hot days in the NECP is slightly smaller than in other regions of China (figure 2(j)).
The results for T WB and T WBG differ from those for T air . Although the apparent cooling of maximum T air is observed in NI, a significant increase in maximum T WB (∼0.5 • C) occurs in this region ( figure 2(b)). This increase is much larger than the increase in other regions (∼0.1 • C) of India ( figure 2(h)). This increase of urban maximum T WB in NI is caused mainly by the significant increases in urban maximum Q and RH% (figures S9(a) and (b)). In addition to the increase of daily maximum T WB , the number of hot days based on T WB also shows an apparent increase (∼4 d month −1 ) in NI compared to other regions of India (∼1 d month −1 ). However, for cities in China, no difference in ∆T WB is evident between regions with apparent AEI increases (NCP and NECP) and other regions. Furthermore, due to the decrease of Q in the NECP ( figure S9(a)), its maximum T WB and number of hot days show slight decreases.
For T WBG , which further considers downward solar radiation and wind speed, similar results to T WB were observed in both China and India. The increases in maximum daily T WBG and the corresponding number of hot days in NI are much greater than in other regions in India (figures 2(c), (f), (i) and (l)). The time series of T air , Q, T WB, and T WBG for cities in India also indicate that cities in NI experienced apparent T air cooling along with significant increases of Q, T WB, and T WBG compared with other regions (figure S10).

Correlation between expansion of the area equipped for irrigation and changes in urban moist heat stress
To investigate the relationship between AEI expansion and UMHS changes in recent decades, we analyzed the correlations of ∆T air , ∆T WB , and ∆T WBG with ∆AEI in NI, NCP, and NECP. The corresponding Spearman correlation coefficients are denoted as the ∆T air -∆AEI, ∆T WB -∆AEI, and ∆T WBG -∆AEI correlation coefficients, respectively. Figure 3 shows a scatterplot of temperature metric changes versus ∆AEI for cities located in NI, NCP, and NECP. Due to the differing climatic conditions of these cities, the results are presented separately for each climate classification, including Csa (warm temperate, summer dry, hot summer), Cwa (warm temperate, winter dry, hot summer), and Dwa (snow, winter dry, hot summer). The AEI changes of cities under each climate classification in China and India are shown in figure S11.
For cities in NI, the ∆T air values show significant negative correlations (p-value ⩽ 0.01) with ∆AEI in the surrounding area (figures 3(a) and (b)). The ∆T air -∆AEI correlation coefficients for Csa and Cwa in NI are −0.56 and −0.31, respectively. By contrast, ∆T WB and ∆T WBG in urban areas show positive correlations (p-value ⩽ 0.01) with ∆AEI, with ∆T WB -∆AEI and ∆T WBG -∆AEI correlation coefficients for both climate classifications greater than 0.5. These results indicate that the increase of AEI in the NI region tends to reduce the daily maximum urban T air in the hottest month, but, intensified irrigation activities also affect other climate variables (e.g. increased air humidity), which increases the risks of moist heat stress in urban areas.
However, this phenomenon is not evident in the NCP and NECP regions of China, and the corresponding correlation coefficients between  (l)). This difference between India and China may result from the following four factors. First, the expansion of AEI is more significant for cities in NI. The average value of ∆AEI in NI is approximately double those in NCP and NECP. This difference makes AEI expansion a more critical factor in shaping the regional climate in NI than in the other areas. Second, the impact of artificial moisture supply from irrigation should be greater under drier conditions. The air is much drier in NI in May (with RH% around 20%-40%) than in the NCP and NECP in July (with RH% around 60%-80%) (figure S12). Third, NI may have larger irrigation water consumption compared with NCP and NECP. For example, Siebert and Döll (2010) showed that the annual irrigation consumption in NI (mainly 400-1922 mm yr −1 ) is apparently higher than NCP and NECP (mainly 250-600 mm yr −1 and 150-400 mm yr −1 , respectively). In addition, the difference in regional irrigation method/efficiency may also influence. For example, Fu et al (2022) found that by implementing water-saving techniques, irrigation water use has been substantially reduced in Northwest China in recent decades, which mitigates the influence of irrigation on regional evapotranspiration, temperature, and humidity. Similarly, such water-saving techniques are also substantially applied in NCP and apparent improvement in irrigation efficiency is found in NCP in recent decades (Zhou et al 2020). Forth, for cities of the NCP and NECP regions, other factors such as urban area expansion and global warming may have strong impacts on the regional climate that mask the influence of AEI expansion.
To further investigate the seasonal variation of the effect of AEI expansion on daily maximum urban T air and UMHS, the ∆T air -∆AEI, ∆T WB -∆AEI, and ∆T WBG -∆AEI correlation coefficients for 12 months are plotted in figure 4. Only correlation coefficients that meet the significant criterion (p-value < 0.05) are shown in color. The results indicate that the cooling effect of AEI expansion on urban T air generally occurs throughout the year in Cwa climate zones in NI, while it is apparent only during winter (November-February) and the pre-and earlymonsoon (May-June) seasons in Csa zones in NI. These results are in accordance with the findings of Mishra et al (2020), who reported that irrigationinduced cooling is more dominant during the premonsoon and post-monsoon seasons in NI based on satellite data. For T WB , the enhancement driven by AEI expansion is most apparent during the preand early-monsoon seasons (February-June) for both Cwa and Csa climate zones in NI and during the post-monsoon season (October-November) for Cwa. In monsoon months with abundant precipitation and air humidity, the effect of AEI expansion on UMHS is insignificant. Similar to T WB , the influence of AEI expansion on T WBG in NI is stronger in the pre-monsoon (April-May) and post-monsoon (October-November) seasons. In the NCP (Cwa) and NECP (Dwa) regions of China, the effect of AEI expansion on maximum urban T air and UMHS is unclear, and no seasonal variation pattern can be identified.
To directly evaluate the correlations between ∆AEI and influencing factors of moist heat metrics, heatmaps of correlation coefficients of ∆AEI with the change in daily maximum Q, solar radiation, wind speed, and daytime mean RH% in urban areas were constructed (figure S13). Changes in Q and RH% in NI show strong positive correlations with AEI expansion during the pre-and early-monsoon season (January-June) and the post-monsoon season (October-November). In addition, wind speed change shows negative correlations with AEI expansion in NI (March-June, September-December), which may be driven by the increased roughness of the land surface caused by cropland expansion around the city. The change in solar radiation has no apparent correlation with AEI expansion.

Relative contribution of the area equipped for irrigation expansion to urban moist heat stress
In addition to AEI change, changes in urbanized areas and regional precipitation may also affect urban T air and humidity. Changes in the urban area proportion of each city square and regional precipitation in recent decades are shown in figures S14 and S15, respectively. To evaluate the contributions of AEI expansion to UMHS compared with those two factors, we applied multivariate linear regression to fit urban ∆T air , ∆Q, and ∆T WB as functions of ∆AEI, urban area proportion change (∆UA), and precipitation change (∆P) (text S5). Using the fitted regression model, the contribution of each factor to urban ∆T air , ∆Q, and ∆T WB can be measured based on a change (increase) of one standard deviation for each factor. The regression model was fitted for 70 cities under the Cwa climate in NI in May, which is both the hottest month and when UMHS changes have strong correlations with ∆AEI (as shown in figures 4 and S13). Figure 5 shows the performance of the regression models (a)-(c) and the contribution of each factor to the output (d)-(f). The results indicate that the regression model can accurately simulate metrics changes from ERA5-Land data and that the variability of metric changes is explained well by three input factors, all of which have R 2 values near or greater than 0.5. Moreover, the regression models have small root mean square error (RMSE) values. In terms of the contribution of each factor to the metrics changes, an increase in urban area proportion tends to increase urban maximum T air in NI, while increases in AEI and precipitation tend to decrease urban maximum T air , with increased precipitation showing a larger cooling effect on T air than ∆AEI. For ∆Q and ∆T WB , increases in both AEI and precipitation tend to increase urban maximum Q and T WB , with AEI expansion making the larger contribution to the increases of Q and T WB . An increase in urban areas has a slight negative effect on Q, and its impact on T WB is negligible compared to the effects of AEI expansion and increased precipitation.

Verifying results with HadISD in-situ data
For validation, we investigated the changes in daily maximum T air , Q, T WB, and daily average RH% of the hottest month between the periods of 1981-1990 and 2001-2010 at 372 stations in China and 72 stations in India using the HadISD dataset (figure S16). Although climate stations have a much lower density than cities, especially in India, similar regional patterns can be observed. In addition, despite the cooling of T air in NI is less evident in the station data than in reanalysis data due to the influences of local conditions around the station and the low station density, much greater increases of air humidity and T WB are observed in NI than in other regions of India. These results indicate that NI has experienced a strong intensification of moist heat stress in recent decades. As NI has the highest population density in India and very high gross domestic product (figure S17), the enhancement of UMHS due to AEI expansion may lead to high exposure and risk for outdoor workers, impacting economic activities.

Discussion
The near-surface meteorological variables of ERA5-Land are used to conduct the main analysis, which uses similar data inputs with ERA5 that combines multi-source observations and model forecasts. Admittedly, assumptions and inadequacies of the model may induce biases to the outputs , which may also partly result in its deviation with respect to ground-based observations (figure S16). However, vast amounts of radiosonde, satellite-based, and in-situ observations are assimilated into the estimation to obtain the most plausible state of the atmosphere. For example, in terms of near-surface field, from 1979 to 2019, there are about one billion observations for each of surface air temperature and relative humidity were processed by the land data assimilation module of ERA5-Land (Hersbach et al 2020). The assimilated data with these comprehensive observations can well reflect the time-evolving and the spatial diversity of the nearsurface climate state, which provides reliable support to measure the influence of AEI change on regional climate.
The Köppen-Geiger climate classification of Kottek et al (2006) is used to classify the cities when conducting correlation analysis, which represents the climate state of 1951-2000. Admittedly, the climate classification is changing under global warming, and some updated versions of Köppen-Geiger climate classification with higher-resolution were proposed based on different input datasets (Beck et al 2018). However, the focus of our research is to detect the influence of irrigation in specific regions with intensive irrigated cropland expansion (e.g. NI), and the climate classification is only used to mitigate the potential influence of climate difference in correlation analysis. Therefore, although there are various versions of climate classifications, we believe it will not change our results and conclusions in how irrigated cropland expansion influences the region and urban climates.
Notably, changes in UMHS can be influenced by global warming and internal climate variability in addition to the regional influence of AEI expansion (Im et al 2017, Mishra et al 2020. For example, El Niño-Southern Oscillation affects monsoon systems, altering temperature and humidity in India and China (Revadekar et al 2009, Zhang et al 2018, Kiran Kumar and Singh 2021. Nonetheless, the main objective of this study is investigating the influence of AEI expansion on UMHS. In-depth quantification of the contributions of various factors influencing temperature and humidity in these regions is a critical research direction for future studies.
In addition to moisture from anthropogenic irrigation activities, urbanization can lead to additional anthropogenic moisture emissions through multiple pathways (e.g. thermal power plants), with effects on humidity. In this study, the specific humidity Q is negatively correlated with ∆UA in NI, demonstrating an urban dry island effect (Hao et al 2018). However, the influence of urbanization on humidity remains under debate and opposing results have been obtained in different cities (known as urban wet and dry island effects) (Wang and Gong 2010, Sailor 2011, Hao et al 2018. Future comprehensive analyses of UMHS should consider these effects.

Conclusions
This study presents novel findings regarding the influence of irrigated cropland expansion around cities on UMHS in more than 1000 cities in India and China. Using state-of-the-art reanalysis and in-situ climate data, we conducted the first analysis focused specifically on urban areas and found that during 1980-2005, cities in NI, NCP, and NECP have experienced significant increases in the AEI in their suburb areas of 14.1%, 6.8%, and 8.2%, respectively. In NI, AEI expansion tends to result in urban T air cooling. Conversely, as irrigation activity increases air humidity, AEI expansion shows positive correlations with the increases of T WB and T WBG , indicating enhancement of the UMHS. This effect is most apparent in the premonsoon and post-monsoon seasons, when precipitation is low. However, no such effects of AEI expansion are observed for cities in NCP and NECP in China due to the differences in AEI expansion intensity and climatic conditions. The results of multivariate linear regression analysis support AEI expansion having stronger impacts on exacerbating the UMHS than precipitation changes and urban area expansion in NI over recent decades.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.