Reduction in human activity can enhance the urban heat island: insights from the COVID-19 lockdown

The COVID-19 lockdowns drastically reduced human activity, emulating a controlled experiment on human–land–atmosphere coupling. Here, using a fusion of satellite and reanalysis products, we examine this coupling through changes in the surface energy budget during the lockdown (1 April to 15 May 2020) in the Indo-Gangetic Basin, one of the world’s most populated and polluted regions. During the lockdown, the reduction (>10%) in columnar air pollution compared to a five year baseline, expected to increase incoming solar radiation, was counteracted by a ∼30% enhancement in cloud cover, causing little change in available energy at the surface. More importantly, the delay in winter crop harvesting during the lockdown increased surface vegetation cover, causing almost half the regional cooling via evapotranspiration. Since this cooling was higher for rural areas, the daytime surface urban heat island (SUHI) intensity increased (by 0.20–0.41 K) during a period of reduced human activity. Our study provides strong observational evidence of the influence of agricultural activity on rural climate in this region and its indirect impact on the SUHI intensity.


Introduction
Human-induced changes in the Earth's surface climate have traditionally been difficult to constrain [1,2], particularly since these changes typically occur at time scales similar to natural decadal perturbations. In contrast, the COVID-19 lockdowns-the restrictions placed by various governing bodies as a response to the COVID-19 pandemic in 2020caused unprecedented slowdown in human activity [3], were short in duration, and yet intense enough to produce measurable effects. They can thus serve as natural experiments on the anthropogenic control on surface climate.
An ideal study area to explore the results of this perturbation experiment is the Indo-Gangetic Basin (IGB), one of the most densely populated regions on the planet with high levels of air pollution [4]. Covering the majority of North India, the IGB has a subtropical monsoon climate, and is a global hot spot for land-atmosphere coupling [5]. From late March to end of May, a countrywide lockdown was imposed in India as a response to COVID-19. This lockdown strictly restricted people's movement outside their homes, suspended educational, industrial, and hospitality services, and limited all transportation systems [6].
Being both heavily cultivated and rapidly urbanizing [7,8], the IGB is strongly influenced by anthropogenic changes in land use and land cover [9,10]. Human influence on surface climate is distinct for urban and rural areas. Urbanization modifies the biophysical properties of the surface due to replacement of natural landscapes with built-up structures. Consequently, cities are usually characterized by higher temperatures compared to their surroundings-the urban heat island (UHI) effect [11]. The UHI is commonly calculated as the temperature difference between the city and a non-standard rural area around the city. Rural areas are influenced by land use and land management practices, which is primarily due to agriculture in this region.
The reduction in atmospheric aerosols during the lockdown (up to 45% reduction for some Indian states) is well-documented [12] and, all else remaining constant, would increase incoming surface radiation. Since urban and rural areas may have different levels of pollution, this radiative forcing change can impact the UHI intensity [13,14]. However, the UHI, and surface climate in general, is also modulated by non-radiative pathways [15]. The lockdown restrictions delayed crop harvesting [16], which would allow the rural surfaces to maintain high vegetation cover and can enhance latent heat flux (λE), a non-radiative pathway of surface heat dissipation. Previous studies have noted that the seasonality of the UHI in this region is influenced by the variability in surface vegetation in the rural area [17,18]. Agricultural influence on surface climate in this region is not well-captured by land-surface models due to inaccurate representation of vegetation properties and the poorly constrained influence of irrigation on the hydrological cycle [10,19]. Thus, beyond the widely studied changes in atmospheric composition [20,21], the lockdowns provide a unique opportunity to ask broader questions about human-landatmosphere interactions in the IGB. For example, the role of these interactions in modulating the ensuing South Asian monsoon circulations is critical to the livelihood of over a billion people dependent on this region [22].
Here, we examine the changes in the surface temperature in the IGB during the lockdown with a focus on the UHI, allowing us to separate human influence on urban and rural surface climate, particularly due to air pollution and land use. First, using a suite of satellite observations, we isolate changes in atmospheric and surface properties over urban areas and their rural periphery. Second, since satellite observations are restricted by the presence of clouds and do not directly provide estimates of the surface energy budget, we use a reanalysis dataset to gain a mechanistic understanding of the observed perturbations in rural surface climate.

Urban-rural Delineations for region of interest
We only consider the Indian portion of the IGB (figure S1(a); ∼50 000 km 2 (available online at stacks.iop.org/ERL/16/054060/mmedia)), to avoid uncertainties arising from variations in lockdown periods in other countries. An urban cluster database was developed for this region by generating polygons covering contiguous groups of high-density urban pixels from the Global Human Settlement Index for 2015 [23]. The corresponding normalized rural reference (roughly equal to area of urban cluster) for each of these clusters was created using an iterative buffering procedure with a step size of 300 m. Our methodology generates 1420 urban-rural delineations for the region.

Satellite observations
Multiple satellite-derived products were used to extract urban-rural differentials in relevant variables (details in table 1), including land surface temperature (LST) for calculating surface UHI (SUHI) [24], reflectance data from Moderate Resolution Imaging Spectroradiometer (MODIS) to estimate surface vegetation [25], and metrics of air pollution, including aerosol optical depth (AOD) from MODIS [26] and individual pollutant estimates from the Sentinel-5P TROPOspheric Monitoring Instrument (TRO-POMI) sensor [27]. We also extract cloud fraction (CF) from Sentinel-5P since clouds strongly affect the radiative budget and can be higher over cities [28]. Finally, the black-sky albedo (BSA) and whitesky albedo (WSA), the reflectivity of the surface for direct beam and diffuse radiation, respectively, were extracted from MODIS [29]. These can be combined to derive total surface albedo (α; see section 2.3), which is known to vary between urban and rural areas [30]. Although the MODIS observations are daily, the multi-day (see table 1) composites are used to reduce cloud contamination [31].
The normalized difference vegetation index (NDVI) is a proxy for green vegetation [33] and used here to estimate the impact of the lockdown on surface vegetation cover. We calculate NDVI from the near infrared (NIR) and RED bands of the 8 day composite MODIS surface reflectance product, available for 1 km × 1 km grids (table 1), as: The lockdown in India started from midnight 24 March and continued in a limited capacity till 7 June. To remove the noise from the transition periods, we considered 1 April to 15 May 2020 to be the lockdown case. The five year (2015-2019) mean of the data from 1 April to 15 May was treated as the baseline (only 1 April to 15 May 2019 for TROPOMI due to data unavailability). Since the satellite observations are relatively coarse (see table 1), only urban clusters with an area of above 10 km 2 were considered. This threshold, along with cloud screening, leaves 382 clusters for the MODIS and Sentinel-derived CF data (figure S2) and 302 clusters for the Sentinel-derived air pollutant data.
Urban and rural means of all the variables of interest were extracted after regridding to 300 m European Space Agency Climate Change Initiative (ESA CCI) grids using the Google Earth Engine platform [34]. The urban values were calculated as the spatial means of all the urban pixels, as defined by the ESA CCI land cover data [32], within an urban Land cover Yearly 300 m [32] cluster. The corresponding rural values are the spatial means of the non-urban, non-water pixels (from the ESA data) in the rural references. The urbanrural differential in LST is the SUHI, while for the other variables, we use the subscript urb-rur . We also calculate the averages of each variable (and their differences) weighted by the urban cluster areas. Since larger urban areas are known to have higher SUHI intensity, area weighing gives us regional mean SUHI (versus the urban cluster mean SUHI).

Reanalysis data
We used two reanalysis products-the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis 5 (ERA5) [35] For each case, the centroid of the urban and rural polygons was located and k d was extracted for the ERA5 grid containing it, similar to [38]. The MERRA-2 reanalysis assimilates biascorrected satellite observations of aerosols and provides estimates of the aerosol direct radiative effect, making it ideal for studying the impacts of the COVID-19 lockdowns. Though MERRA-2 also includes estimates of k d , inter-reanalysis evaluations show that it significantly underestimates k d compared to other products [39]. Thus, MERRA-2 data were not used to calculate α. Instead, MERRA-2 was used to get all-sky estimates of perturbations in the surface and atmospheric variables (since satellites only provide clear-sky estimates) and to diagnose reasons for perturbations in the rural LST. The MERRA-2 variables were based on hourly data from 12:30 pm to 2:30 pm local time, corresponding to the 1:30 pm MODIS Aqua overpass.

Statistical analysis
To examine reasons for any potential SUHI change, we considered the temporal changes (∆, variable value during the lockdown minus value of the reference period) in AOD urb-rur , NDVI urb-rur , α urb-rur , and CF urb-rur as the predictors. For robustness, we used two methods-linear regressions and random forest (RF) regression. For the linear regressions, we considered each of these four predictors (∆AOD urb-rur , ∆NDVI urb-rur , ∆α urb-rur , and ∆CF urb-rur ) and all their combinations and subsets. Since the relationships between the predictors and ∆SUHI are not necessarily linear, we also checked the consistency of our results using RF regressions. RF regressions use ensembles of decision trees to detect non-linear relationships and are less sensitive to outliers than parametric linear models [40]. For the RF regression, to prevent overfitting, we trained the models using 70% of the data and checked the model accuracy using the remaining data. The training and accuracy assessment were repeated 50 times with different random splits of training and validation data [41]. The signal can be hard to separate from the noise in the satellite observations when examining small perturbations. Thus, we also correlate the variables after binning the data into five percentile bins, assuming that this noise is random and contributes to the unexplained variance within bins.

Intrinsic biophysical mechanism
Since the statistical analyses using satellite observations are primarily for hypothesis generation and do not necessarily reveal the physical mechanisms for the SUHI increase, we employed the theory of Intrinsic Biophysical Mechanism (IBPM) [14,15,42,43], implemented using the MERRA-2 variables [43], to diagnose and quantify the reasons for the change in the LST rur in the region. Conceptually, the total LST change (∆LST) due to a forcing agent is the sum of the changes in the blending height temperature (∆BHT), where the blending height is the height at which surface heterogeneity has negligible impact on atmospheric variables [44], and the local temperature response (∆T): ∆BHT is the result of atmospheric factors while ∆T is the surface response to atmospheric forcing. According to the IBPM theory, the local temperature response is: The terms on the right-hand side of equation (4), from left to right, are the contributions to ∆T from shortwave radiative forcing, longwave radiative forcing, energy redistribution through evaporation, energy redistribution through convection, and change in ground heat flux. Here, ∆K ↓ is the change in incoming shortwave radiation, ∆L is the change in absorbed longwave radiation, ∆G is the change in ground heat flux, and R * n is the apparent net radiation. The intrinsic climate sensitivity is given by λ 0 , f is the energy redistribution factor and ∆f 1 and ∆f 2 are the changes in f due to evaporation and convection, respectively. Check supplementary material for more details.

Changes in urban-Rural differentials during lockdown
The daytime SUHI increased (non-weighted mean from 0.23 K to 0.43 K; area-weighted mean from 0.56 K to 0.97 K) during the COVID-19 lockdown period compared to the five year baseline (figures 1 and S3), with 67% or 257 of the 382 urban clusters showing an increase (figure S4). In contrast, the nighttime ∆SUHI was statistically insignificant (p-value > 0.01). Daytime LST decreased during the lockdown (compared to baseline), though the mean decrease over rural locations was greater (−1.96 K non-weighted; −1.55 K weighted) compared to that over urban regions (−1.76 K non-weighted; −1.14 K weighted).
To explain this differential perturbation in daytime LST urb and LST rur , we consider major atmospheric (CF and AOD) and surface (NDVI and α) factors that impact the SUHI. The air quality over both urban and rural locations improved during the lockdown. This reduced air pollution is evident from both MODIS-derived AOD (∼5%-6% decrease; figures 2(a) and S5(a)) and individual air pollutants measured by TROPOMI (∼14% and ∼19% decrease in urban NO 2 and SO 2 , respectively; figure S6). AOD decreases in roughly 67% of the urban clusters (increasing in others) and for most non-urban grids in the IGB ( figure S7(a)). The small magnitude of MODIS-derived AOD urb-rur could be due to the large contributions from natural dust and biomass burning aerosols, as well as occurrence of turbulently mixed deep boundary layers, in this region during this period [45]. Moreover, there was a surprisingly large (>36% for non-weighted; >43% for weighted) increase in CF during the lockdown (figures 2(d) and S5(d)), though some regions show a decrease ( figure S7(b)).
We find large (>12%) increases in NDVI during the lockdown period. These increases are generally higher in rural references (13.3% or 0.043) than in urban areas (12.8% or 0.038; figure 2(b)), with ∆NDVI urb-rur being negative in ∼59% of the urban clusters (figures S4(b), S7(a), and S8(a)). The positive ∆NDVI rur (figure S8(b)) is consistent with the impact of the lockdown on agricultural activity, which is the predominant land use for ∼88% of this area according to the ESA CCI data ( figure S2(b)). The lockdown overlapped with the harvesting season for rabi (winter) crops, which, together with the drier conditions during this period, reduces surface vegetation during regular years, as can be seen from the moving average of MODIS-derived 8 day NDVI rur in figure 3(c). The delay in harvesting activity [16] may have temporally shifted this normal drop in NDVI rur , contributing to surface greening compared to the baseline ( figure 3(c)). The ∆NDVI rur is higher than the inter-annual standard deviation of NDVI rur during this period, also demonstrated by the largest standardized anomalies of the year during the lockdown. Most of IGB shows this surface greening, with ∆NDVI rur , ranging from −0.01 (5th percentile) to 0.11 (95th percentile), as well as a reduction in LST rur (figures 3(a) and (b)) from −4.3 K (5th percentile) to 0.4 K (95th percentile), with the SUHI increasing compared to baseline ( figure 3(d)). Simultaneously, α decreased in both urban (−2.6%) and rural (−2.9%) areas (figures 2(c) and S5(c)). Paired two-sample t-tests confirm that all the satellite-observed changes in the variables, other than that for nighttime SUHI, are statistically significant (p-value < 0.01). The 95% confidence intervals for the mean of each variable are in table S1.

Attribution of daytime SUHI enhancement
Previous studies have shown relationships between AOD urb-rur and SUHI [13,14]. Here, we find almost no associations between the perturbations in the two during the lockdown period (figure 4(a); r 2 = 0.02 for cluster; not statistically significant for binned). Similarly, AOD urb-rur and daytime SUHI are not wellcorrelated during the baseline and lockdown periods ( figure S9(b); r 2 ∼ 0). However, ∆NDVI urb-rur shows a relatively strong inverse relationship with ∆SUHI ( figure 4(b); r 2 = 0.16 for cluster; 0.76 for binned). Given the generally higher ∆NDVI rur , the pause in human activity in the rural area may have contributed to the enhanced daytime SUHI. This impact of urban-rural vegetation differentials on SUHI is consistent with previous studies [46,47] and is corroborated by the relatively strong associations between NDVI urb-rur and daytime SUHI for the two periods ( figure S9(a)). Although ∆α urb-rur shows a positive relationship with ∆SUHI (r 2 = 0.05 for cluster; 0.52 for binned), this relationship is not physically possible in isolation, since a higher α implies lower absorption of solar insolation by the surface, and thus, lower ∆SUHI. Since α and NDVI are correlated (figure S10), the positive relationship between ∆α urb-rur and ∆SUHI may be a statistical artifact of the relatively higher NDVI rur . Neither CF urb-rur and daytime SUHI intensity, nor their perturbations from the baseline to the lockdown, are correlated (figures 4(d) and S9(c)).
We use multiple linear and RF regressions to provide further statistical robustness to our findings (table S2). In general, ∆NDVI urb-rur explains the largest portion of the variability in ∆SUHI (adjusted R 2 = 0.15), followed by ∆α urb-rur , ∆AOD urb-rur , and CF urb-rur . The permutation importance scores from the RF models also support the primary control of ∆NDVI urb-rur on ∆SUHI (table S2).
The percentage changes in the atmospheric variables (AOD and CF) are higher over urban areas compared to the rural areas during the lockdown period and the surface properties (NDVI and α) change more for rural areas (figures 2 and S5). Although 2020 was already wetter than regular years [18], a further reduction in AOD over the study area may  have perturbed the regional circulation and thermodynamics [48,49], creating an environment conducive to more cloud formation. It may also be possible that the lower aerosol loading reduced the cloud-burning effect [50], thus increasing CF. Finally, the surface greening could also enhance CF [51]. We expect a greater decrease in AOD over urban clusters to correspond to an increase in K ↓ , which can enhance SUHI intensity. However, our statistical analysis does not support this hypothesis. Instead, the observed positive ∆SUHI is associated with the higher vegetation cover over rural areas. This control of vegetation cover on LST is further corroborated by the negative correlations between LST and NDVI for the urban and rural units (figure S11). The relatively weaker correlations for the cluster-level data in figure 4(b) compared to figures S9(a) and S11 (and previous studies [46,47,52]) is because we are dealing with differences of differences in figure 4(b) during a time of the year with low expected SUHI intensities ( figure 3(d)), making the signal hard to isolate from the noise.

Perturbations to rural Background surface climate
The lack of large-scale continuous observations of meteorological variables in urban areas makes diagnosing these patterns for all-sky conditions difficult using in-situ data. Instead, we use the MERRA-2 reanalysis, which is observationally constrained by ground-level measurements of surface meteorology and satellite measurements of columnar AOD, and physically constrained by the model components [36]. MERRA-2 primarily represents the rural background since it does not incorporate urban land cover. Since, as suggested by the satellite-derived NDVI and α (figure 2), urban surfaces changed less than rural surfaces during the lockdown, the reanalysis data can be used to generate mechanistic insights about the SUHI enhancement.
The MERRA-2 reanalysis captures the direction of the changes in the region during the lockdown compared to satellite observations (table  S3). Although the midday aerosol direct radiative effect in MERRA-2 decreases by almost 25% (from −78.4 ± 12.6 W m −2 to −58.7 ± 13.3 W m −2 ) during the lockdown, with the potential to increase K ↓ by 19.8 ± 5.1 W m −2 , we find an overall reduction in K ↓ (−10.5 ± 39.8 W m −2 ). This decrease in K ↓ is due to the compensating effect of increased cloudiness during the lockdown, as well as the higher water vapor content in the atmospheric column, as seen from the higher near-surface relative humidity (table Figure 5. Contributions of different pathways to regional land surface temperature (LST) change. Sub-figure (a) shows contributions of all pathways, namely shortwave radiative forcing, longwave radiative forcing, evaporation, convection, and ground heat storage to the total calculated local temperature change (∆T) in the Indo-Gangetic Basin during the lockdown compared to the five year baseline. The corresponding ∆T, ∆LST, and change in blending height temperature (∆BHT) in MERRA-2 are also shown. The standard errors are displayed in all cases. Sub-figure (b) shows the correlations between the MERRA-2 grid-averaged leaf are index (LAI) and the corresponding ∆T from MERRA-2, the calculated ∆T using the IBPM framework, and the contributions to ∆T through the evaporative pathway and convective pathways. Each data point corresponds to a grid cell average. The lines of best fit are shown and the corresponding equations (including confidence bounds for the slopes of the lines, sample sizes, and p-values) are annotated. Values that are outside the 1-99 percentile of the total diagnosed local temperature change are considered outliers and not shown in the scatter plot. S3). Overall, the total absorbed energy by the surface decreases slightly (−10.4 ± 32 W m −2 ) during the lockdown despite the negative ∆AOD rur .
Separating the contributions from both radiative and non-radiative pathways that can change LST reveals large evaporative cooling (−1.79 ± 0.05 K) during this period ( figure 5(a)), which is expected if vegetation cover increased. The diagnosed and MERRA-2 calculated ∆T are similar in magnitude (−1.29 K versus −1.22 K) and in spatial distribution ( figure 5(b)). MERRA-2 uses prescribed vegetation, with identical leaf area index (LAI; regional mean = ∼0.79) for the two periods. However, since it is constrained by observed surface meteorology, it captures the decrease in Bowen ratio-the ratio of sensible heat flux (H) and λE-during midday (from 2.87 ± 3.02 for baseline to 1.24 ± 0.76 during the lockdown; table S3), which is an expected impact of surface greening. The higher precipitation, latent heat, cloud cover, and relative humidity point to a more intense hydrological cycle during the lockdown compared to regular years (table S3).
The satellite-observed negative ∆NDVI urb-rur suggests that the increase in evaporative cooling during the lockdown was more for rural areas compared to urban areas. This would be true even if NDVI urb and NDVI rur had changed identically, since urban areas are generally more moisture-limited. Figure 5(b) shows the correlations of ∆T from MERRA-2, the IBPM calculations, and contributions due to evaporation and convection with LAI. The negative correlations between LAI and ∆T demonstrate a stronger cooling response during the lockdown over more densely vegetated surfaces. Here, the lowest LAI grids represent relatively urbanized areas. The increasing positive temperature response through the convective pathway with LAI suggests that this is a negative feedback to the evaporative cooling [53]. This convective feedback can be understood either in terms of energy conservation or through evaporation-induced near-surface stability. A relative increase in λE under similar (or reduced) available energy requires a corresponding decrease in H. Alternatively, the additional evaporative cooling of the surface compared to the blending height renders the lower atmosphere relatively stable compared to the baseline period, impeding the dissipation of available energy via H.
Overall, ∆T is almost a third of ∆LST, with the other two-thirds attributable to atmospheric factors, including AOD and CF ( figure 5(a)). The IBPM results show that the sum of the evaporative cooling and its convective feedback accounted for roughly 79% of the midday ∆T, while the evaporative cooling alone accounted for roughly 46% of the corresponding ∆LST (of −3.92 K). The land contribution to ∆LST found here is probably a lower bound estimate since there is strong coupling between the land and the atmosphere. For instance, enhanced surface evaporation due to the increase in vegetation cover would also increase low level cloudiness through condensation feedback [54], lowering both BHT and LST. A similar theoretical diagnosis is not possible for the urban surfaces explicitly since the MERRA-2 land cover dataset does not consider urban areas. Nonetheless, these results can explain the SUHI increase, as the implicit assumption is that the surface characteristics of the rural areas changed more than those of the urban areas during the lockdown period, which is reasonable given the time scale.

Discussion and conclusions
The UHI effect is traditionally viewed as an outcome of the replacement of the natural landscape by built-up structures. The consequences of this land cover change are simpler to define in the abstract than to measure in practice. While cities modify their local climate as they expand, the UHI intensity is usually quantified through space-for-time substitution using snapshot measurements in the urban area and for some rural reference. For the SUHI, how to define this rural reference remains a contentious issue [46,[55][56][57]. Generally, the urban-rural delineations are more clearly constructed for less sprawling cities with very little land management surrounding the city centers. Here, we show an example of a region of the world with high human intervention in both urban and rural areas, the interruption of which leads to the seemingly counter-intuitive enhancement of the SUHI during a period of low human activity. Thus, the COVID-19 lockdown period illustrates the importance of the rural reference on the SUHI intensity, as estimated using a traditional buffer-based method, during a perturbation scenario, the relevance of which has previously only been examined for the mean climate state in this region [17,18,58].
Our results can help contextualize a larger current discussion in the urban climate community on the utility of the UHI as a metric to examine urban public health [59][60][61]. The UHI intensity is the impact of urbanization on local temperature [62]. However, urban heat stress is dependent on the absolute temperature, or more accurately, a combination of temperature, humidity, and other factors [63]. As such, the relevance of UHI for urban public health can be misleading during certain times because enhancement in UHI intensity does not necessarily imply similar enhancement in heat stress (or even temperature) in urban areas. In agreement [64], argues that mitigating the UHI should not be the goal when addressing the public health consequences of urbanization. In theory, one can reduce the UHI intensity by increasing the rural temperature, which does not change the potential heat stress in the urban area. Here, we see something similar occurring, with the SUHI increasing due to rural areas cooling down more than the urban core between the baseline and the lockdown periods, rather than due to an increase in urban temperature. However, this is only the temporal perspective. From the spatial perspective, it is also true that residents moving from rural to urban regions in the IGB were exposed to higher temperatures than they would have had they remained in rural areas. This separation of the temporal and spatial perspectives is critical to reconciling the debate in the community. The criticisms of UHI as a metric primarily pertain to the total impact of temperature on human health in urban areas. In contrast, since the UHI is an abstract isolation of the contribution to that temperature from urbanization, it remains theoretically important, assuming we establish a more consistent definition of the rural reference to facilitate accurate inter-urban and inter-study comparisons.
Lastly, our finding demonstrates the importance of human-land-atmosphere coupling on the regional climate over South Asia as a whole. Agricultural practices in this region strongly control the vegetation phenology of the croplands, modulating how energy is partitioned and dissipated from the surface through non-radiative means (figure 3(c)). As seen here, the importance of these non-radiative components on the LST is apparent, even when input energy to the system is relatively stable, providing large-scale observational evidence of previously modeled results [54]. Over 88% of the landmass in North India is agricultural ( figure S1(b)). Therefore, an increase in surface vegetation due to agriculture can lead to large regional cooling, modify cloud formation, and lower tropospheric stability. Moreover, enhanced evaporative cooling over the Indian landmass during the monsoon onset period (as seen here) can also perturb the land-sea thermal gradient-a major driver of monsoonal wind circulation [65]. Our study puts forth important observational evidence of human-induced control on surface climate, which strengthens the need for the ongoing efforts to explicitly include these in earth system models to better predict long-term climate change [66,67]. For the UHI, the inclusion of human dynamics can help constrain its future estimates, since urban and rural areas are expected to change differently in future scenarios [68].
Although our two-pronged approach using satellite observations and reanalysis product demonstrates the consistency of these perturbations, a few uncertainties remain. First, since 2020 was wetter than the baseline years (table S3), the NDVI perturbations seen during lockdown may not have been solely due to the well-documented delay in harvesting in this region [16]. The lockdown-induced anthropogenic pause could influence the natural variability, cloud cover, and rainfall [69,70], and in turn also affect the NDVI/LST. Regardless, the large, standardized anomaly in NDVI rur during the lockdown, seen in figure 3(c), strongly suggests that the lockdown played a role. The higher ∆NDVI rur (compared to ∆NDVI urb ) and ∆SUHI is also see when using 18 year (2003-2019) baseline from MODIS Aqua measurements (figures S12(a) and (b)). Moreover, NDVI differences are seen at urban-to-urban periphery scale (5-30 km), which is much smaller than the inherent spatial scale of the anticipated natural variability. Second, since the perturbations are small in magnitude, sensor noise could account for some of the variability. Our results are qualitatively replicated when we calculate the relevant variables from the MODIS Terra satellite (figures S12(c) and (d)), which has a different sensor and equatorial overpass time (∼10:30 am), indicating that the perturbations are robust and cannot be just random noise from one sensor. Third, since our study deals with regional changes using coarse satellite observations, we neither fully examine the perturbations for individual urban clusters, which can vary from the mean changes (figures S4(a) and (b)), nor characterize intra-urban distributions. Some of these limitations can be addressed with the development of better parameterized models for this region with explicit irrigation schemes, which can clearly isolate the impact of the agricultural cycle on regional climate.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.We acknowledge the Yale Institute for Biospheric Studies for a grant on studying aerosol-UHI interactions.