The simulation has been conducted using Weather Research Forecast–Advanced Weather Forecast (WRF–ARW) model version 3.9.1 (https://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/user_guide_V3.9/contents.html ). The model uses terrain-following, hydrostatic-pressure vertical coordinate with the top of the model at constant pressure surface. The horizontal grid is the Arakawa-C grid (Skamarock et al., 2005). The WRF-ARW model consists of two parts, namely, the WRF pre-processing system (WPS) and WRF sections (Figure 1). The three parts related to geogrid, ungrib, and metgrib were combined in WPS in which geogrid was used to define the model domain and grid size, and interpolates the land surface boundary condition, i.e., LULC, into the grids (Figure 1). The model has to itself define the LULC categories and external LULC are reclassified as per defined classes as per USGS 24-category Land Use. Ungrib was used to extract information from meteorological variables, and it interpolates horizontal grids into the model grid, which is processed from static data in geogrid and accomplished in metgrib. The level-4 MODIS LAI are developed every 8 days at a 1-km resolution on a sinusoidal grid, which is used directly by the “Noah” land surface model within WRF (Kumar et al., 2014). The Noah land surface scheme gets modified with land use data that adjust the parameters of the vegetation fraction, vaporization, surface roughness (albedo cover), and so on for each grid and each domain. These parameters work for downscaling the data as per new scenarios of LULC. There are various tools like NCL, GrADS, or CDO that can be used for post-processing visualizations (Samie et al., 2020).

In the study, the nested domain with ratio 1:3 was used, in which D-01 referred to as a parent domain, D-02 as an intermediate domain, and D-03 as the main domain (Figure 2A). The parent domain (first nested domain) D-01 was fixed between 40°E–117°E and 5°S–45°N with 81 km grid spacing, comprising 8,560 grids (56,162,160 km²). The intermediate domain (second nested domain) D-02 was fixed between 58°E–110°E and 1°S–43.5°N with 27 km grid spacing, comprising 44,037 grids. The main domain (third nested domain) D-03, which covered the entire Indian sub-continent and fixed between 75°E–111°E and 15°N–40°N with 9 km grid spacing, comprises a total of 210,202 grids. These three domains were used at 30 vertical pressure levels with 50 hPa as top vertical level.

The WRF model uses a large-scale forcing atmosphere and land for initialization of initial and lateral boundary conditions, and land surface boundary conditions. For WRF simulation, National Centres for Environmental Prediction (NCEP) global forecast system (GFS) reanalysis data at 1° horizontal resolution (Saha et al., 2010a; Decker et al., 2012), NCEP Climate Forecast System Reanalysis at 0.50° horizontal resolution (Saha et al., 2010b), and European Centre for Medium-Range Weather Forecast (ECMWF) reanalysis datasets at 0.75° horizontal resolution, which are available at 6-h intervals (i.e., at 00:00, 06:00, 12:00, and 18:00 h UTC), were used as boundary conditions. The details about various types and availability of boundary conditions can be found on the WRF user page WRF NCEP-NCAR GRIB DATA. In the study, initial and lateral boundary conditions were taken from ECMWF datasets (Dee et al., 2011) at a 6-hourly time scale for 7 years from 2009 to 2015 and 2009 was considered as a spin-up time.

This boundary condition at 6-h intervals in WRF was interpolated and updated as per grid spacing and model domain. The land surface boundary condition is already present in the WRF model as “Noah-modified 21-category IGBP-MODIS land use” with a horizontal grid spacing of 10 min, but it was replaced by MODIS—IGBP LULC of years 2002 and 2015. Due to the difference in the classification of WRF IGBP LULC and external LULC, it has been modified to make it conformal for analysis as shown in Table 1. The WRF model has varieties of physical parameterization schemes, which are used for different types of sensitivity analysis or simulations. Parameterization based on the various previous studies (Li et al., 2012; Paul et al., 2016; Paul et al., 2018) was considered for simulations. The model configuration and parameterization used are shown in Table 2. Finally, we conducted simulations for two different time periods of LULC, i.e., of 2002 and 2015 with an ECMWF reanalysis climatic data. The first simulation was performed with LULC of 2002, while the second simulation was performed with LULC of 2015 using the ECMWF climatic reanalysis data to quantify anthropogenic influences. Any differences in climatic conditions between 2002 and 2015 are automatically taken into consideration by the WRF model, and the resulting mean difference between the WRF simulated data for precipitation and temperature change of 2002 and 2015 was considered as changes due to LULC as a footprint of anthropogenic activities.

To test the correlation between WRF simulated output and observed data from IMD for both temperature and precipitation, the spatial correlation method has been used using Eq. 1. r x y = ∑ ( x i − x ¯ ) ( y i − y ¯ ) ∑ ( x i − x ¯ ) 2 ( y i − y ¯ ) 2 . (1) r _xy is the correlation coefficient; x _i is the sample x-variable; x̅ is the mean value of x-variable; y _i is the sample y-variable; ȳ is the mean value of y-variable. The LAI changes were estimated in terms of percentage change (relative to 2002) using yearly mean of LAI for 2002 and 2015. The significance of the mean and median, temperature, and precipitation change across different LULC types between 2002 and 2015 was checked using a two-sided t-test (for mean) and Wilcoxon signed-rank test (for median).

The reanalysis datasets as acquired from NCEP-CFSR and ECMWF were used as boundary conditions to simulate in WRF models. The time-series daily datasets of precipitation (GPM) and temperature (IMD) for the year 2015 were used for validation of these simulated datasets. Since 2009 was considered as a spin-up time, it was not included in the analysis. The high correlation in precipitation was observed in the eastern, central, southeastern, and western parts of India in contrast to the low correlation in the extreme northern parts (Figure 3). In case of temperature, a very high correlation was evident in the major parts of India that may be attributed to high anthropogenic influence in the form of LULC alterations. A similar observation of precipitation pattern was observed in all observed and simulated datasets for the year 2015. The low precipitation (∼>200 mm) was observed in the extreme northern parts, as observed in both GPM and IMD datasets, while simulated, and ERA-5 exhibited conformity of moderate precipitation (∼600 mm) in the major parts of the country (Figures 4D–G). Simulated precipitation was overestimated compared to observed datasets primarily in northeast (NE) India. While simulated precipitation was higher in the upper parts of western Ghats compared to the other three datasets (GPM, IMD, and ERA5) due to the LULC alterations (dense forest cover), IMD, ERA-5, and simulated datasets exhibited similar patterns in the remaining parts of the country, in contrast to GPM datasets. The observed (IMD), reanalysis (ERA-5), and simulated temperature data exhibited a similar pattern in major parts of India, except for a purely different pattern in extreme northern India. A comparatively less number of observations of IMD stations over the Himalayan region may lead to an error in interpolation in IMD datasets (Kishore et al., 2016; Li et al., 2018). Furthermore, the changes may be simulated more accurately as, in 2015, there was a larger snow cover area that increases the albedo and reduces the temperature over the region.

MODIS-based LULC datasets for the year 2002 and year 2015 were used separately in WRF modeling to simulate temperature and precipitation for the period 2009–2015 to deduce the changes. The study exhibited the maximum transition in the lower parts of central India, where grassland was converted into agriculture, and in western India, where barren land was converted into open shrubland (Figure 5A; Table 3). In the upper western ghats, the agricultural cover was replaced by forest cover. A recent study reported increasing greening over India despite a higher LULCC (Lal et al., 2019; Zhu et al., 2016) due to an increase in LAI (Figure 5B). The increased agricultural activity in the country led by use of fertilizers resulted in cropland LAI increase (Roy et al., 2007). The study also indicated major land transformation from shrubland to cropland (120,090 km²), followed by savannas to forest (34,832.78 km²), cropland to grassland (32,629.61 km²), and barren to grassland (30,978.41 km²) in India between 2002 and 2015 (Figure 5; Table 3). The large-scale land transformation with conclusive decrease was noticed in barren land, which was transformed to grassland (30,978.41 km²; primarily in the western India, as desert area transforming into grassland), followed by snow cover [5,451.79 km²; in Jammu and Kashmir (north India) snow cover increased], wetland (2,605.08 km²), and savannas (423.90 km²). Furthermore, the major increase in urban areas was evident at the cost of cropland (835.42 km²), following grassland (375.87 km²), and savannas (228.81 km²). The conversion of forest cover to savannas (18,247.11 km²) was noted during 2002–2015, while the large portions of savannas were converted to forest (∼34,832.78 km²) and cropland (23,595.59 km²) during the period. The detailed zone wise description of LULCC is discussed in the Supplementary Table.

The study observed a high reduction in LAI (<−5%) in south India, northern parts of north India, and lower parts of northeast India (Figure 5B). Despite browning in some regions, widespread greening was observed in major parts of India, attributed to intense agriculture practices and afforestation measures to support the recent policy, while the areas under browning conditions were largely deforested, and dominantly transformed from savannas to cropland due to agriculture intensification. Greening was observed in the regions having dominant transformation from savannas to forest (evergreen broadleaf forest), primarily tea plantations. Savannas have a higher LAI, which contains more water in their roots and leaves and thus supports high convective precipitation and transportation of moisture to the atmosphere.

The LULCC modifies various parameters of the earth surface, among which temperature is one of the important parameters. Over the Indian subcontinent, we observe a divergent pattern of temperature change, with much of western and north India showing cooling whereas eastern and east-central part of India showing warming due to LULCC (Figure 6). The temperature difference between 2002 and 2015 due to LULCC was significant across all LULC types, except urban and water bodies (Figure 7), with grasslands, open shrublands and barren lands showing decrease in temperature whereas other LULC types showing increase in temperature. Changes in vegetation cover affect the flow of energy between the ground and the atmosphere, resulting in changes in temperature. Greening primarily promotes ET (and thereby increase in latent heat flux), but it does not necessarily lower the temperature as greener canopies have lower albedo than bare soil, and absorb incident sunlight, resulting in the more sensible heat flux. LULCC does not directly alter or regulate temperature, while it is caused and triggered through a series of concomitant changes. Other than greening and browning conditions, the entrapped long-wave radiation through the transpired water by plants rather reduces the amount of shortwave solar radiation through increased cloud formation. Previous studies reported a net cooling trend due to rising greening in the past 4 decades simulated from coupled ESMs (Bathiany et al., 2014; Zeng et al., 2016; Zeng et al., 2018a; Piao et al., 2020; Zhao et al., 2020). Also, in the tropical and subtropical regions, evaporative cooling is generally larger than albedo warming, a similar agreement found with the present study with the greenness augmentation. The temperature changes influenced by LULC change was analyzed using simulated (2015, 2002) and ECMWF reanalysis datasets. India lies under tropical and subtropical regions and has a maximum temperature in southeastern India and northern India (∼>25°C), in contrast to the eastern and western Himalayan region having the lowest temperature (∼4°C).

The LULCC influences temperatures through changes in albedo that can act toward warming when it decreases, but toward cooling when it increases. The differences between the simulated datasets provide the potential impact of different LULCC on temperature (Figures 4A–C). It is observed that the central to eastern parts of India were getting warmer, while the central to western parts were getting cooler. The major country-wide transformation occurred as grassland to cropland, followed by savannas to forest cover during 2002–2015. Overall, the trend of grassland to cropland transformation affects the seasonal phenomenon, while savannas-to-forest transformation tends to have a warming effect (as observed in northeast India). The study showed cooling effects, which were dominant in the western, northern, and south India and lower parts of central India as well. A decrease in temperature (∼0.5°C) is observed in the northeastern Himalayas region primarily due to the increase in snow cover replacing barren land, and parts of barren land converted into shrubland/grassland, which affects the climatic parameters. A significant decrease in temperature (∼0.8°C) was observed in the eastern Jammu and Kashmir (J&K) bordering Tibet, whereas an increase in temperature (∼0.2°C) was noted in the northern parts of J&K. Eastern parts of J&K have had large-scale change in LULC due to anthropogenic activities and natural disasters. In general, there was an increase in savannas and forest by replacing barren land and grassland, while new grasslands were transformed from barren land. A small portion of snow cover in the northern parts of J&K was transformed to barren land, enabling higher surface roughness to trap more heat compared to the snow cover.

Another important parameter is that anthropogenic influence, viz., construction of transportation networks in snow-covered regions, traps more sunlight and contributes significantly to temperature rise. Due to extensive urbanization and mostly agrarian activities throughout the decades and decrease in LAI, the temperature primarily in northern India and the western parts of northeast India rose by ∼0.4°C. Due to the intense irrigation, the vegetation vigor increased during 2015, while LAI decreased by ∼5% in the region. The reduced LAI, counteracting the water content, induces warming in the area, other than urbanization. The increase in urban area modifies surface roughness and thus traps more albedo, which increases the temperature in certain regions. The northern parts of south India and lower parts of central India observed cooling primarily due to the significant transformation of grasslands to croplands, with an increase in greening due to increased LAI (by ∼15%–20%). Comparatively lesser water content in grassland than cropland induce lesser greening in grassland, and thus, such land modification leads to a decrease in the temperature. The conversion of shrubland to croplands was observed in western India that experienced increased precipitation by ∼100 mm (Figure 8C), with an extreme increase in greening (LAI ∼>25%; Figure 6B), leading to a significant cooling. It is notable that the cropland acts cooler due to high availability of water content in the crop field and high leaf water content in contrast to the leaves of grassland and shrubland. The surface roughness in a grassland leads to higher trapped albedo, making the surface warmer. In contrast, the higher leaf water content and less surface roughness in cropland make the surface comparatively cooler during pre-harvesting period. In India, demand and supply for irrigation have tremendously increased (Rodell et al., 2004; Roy et al., 2007) due to the advancement of machinery and tools in recent decades, and with the increased support of local governments to farmers. The remaining parts of India observed anthropogenic induced LULC change with very minimum influence on temperature. The use of simulated precipitation and temperature-based land use/land cover change can be integrated with other factors including change in emission factor, micro-climate change assessment, and pollution using WRF for improved results.

The simulated precipitation (based on the mean of 2010–2015) with reference to the LULC 2002 and LULC 2015 was used to understand the influence of LULC change on precipitation patterns by differencing the two simulated precipitation datasets (Figures 4D–G). The study exhibited an increase in precipitation amount (>300 mm) in the parts of eastern (easternmost region of northeastern Himalayas) and western Himalayas (lower part of J&K), western Ghats (coastal region of central India and south India), and northwestern parts of central India. The precipitation difference between 2002 and 2015 due to LULCC was significantly higher in 2015 (difference >0) across all LULC types, except wetlands and water bodies (Figure 9). In arid regions of western India, the shrubland to cropland, and barren land to shrubland/cropland, were the dominant land transformation with increased simulation precipitation (by 100–150 mm) during 2015. This region also observed the highest increase (>25%) in LAI and attributed to an increase in the cropland area due to intensified irrigation. It influences high leaf water content compared to barren land and shrubland, leading to high ET and more regional precipitation. The decrease in precipitation in northern India could be attributed to intensive farming, and minimum LULCC coupled with minor reduction in LAI (∼8%). Due to increasing urbanization and decreasing surface porosity, the low-level air circulations get reduced, causing instability in the boundary layer, and reduction in precipitation. The negligible changes in the precipitation (<50 mm) were observed in southern and northeastern parts of India, while very high precipitation (>300 mm) was noted in western ghats, attributed to increasing forest cover, despite decrease in the LAI. These regions have primarily been modified from dense forest cover to commercial plantation leading to less greening. Previously, major parts of western ghats were primarily covered with the evergreen needleleaf forest (ENF), which was converted into evergreen broadleaf forest (EBF) having more leaf water causing an increase in the precipitation despite a decrease in the greenness. Negligible changes in the precipitation patterns (>50 mm) were observed in the lower part of central India. This region observed a transition from grassland to cropland that does not much affect the hydrological cycle, while an increase in precipitation (∼300 mm) is evident in the western parts of central India due to the increase of forest cover coupled with increased LAI (by 20%). This may be attributed to the government forest conservation policies and its effective implementation that led to high ET.

The easternmost parts of northeast Himalayas have an exceptional increase in precipitation (∼100–150 mm), which shows similar agreement with the increase of LAI (by ∼15%), due to the transformation of snow cover and barren land to ENF cover, while the remaining parts of northeast Himalayas have no significant changes in precipitation influenced by LULC change. The major parts of northeast India observed a decrease in precipitation primarily attributed to shifting cultivation and decrease in greenness (decreased by ∼10%–15%). The reduced vegetation cover tends to reduce the latent heat flux and evaporation. The decrease in LAI coupled with reduction in albedo would reduce evaporation and transpiration, which induce a decrease in precipitation in NE India. The dense forest in lower parts of J&K observed an increase in precipitation (150–200 mm), with a slight increase in LAI (∼15%) during 2015 and reported an increasing greening pattern.

Precipitation, an important meteorological factor, is a key link to the global and regional hydrological cycle. Convective precipitation is the main source to the major parts of forest in India. The surface water losses to the atmosphere occur through ET, which include transpiration (50%–80%, varies with different plant type) and evaporation (Kemp et al., 1997; Hatfield and Dold, 2019; Wang et al., 2020; Gimeno et al., 2021). Greening expands the loss of water through an extended leaf area, and the study found similar observations in that regions with an increased LAI observed high precipitation barring a few exceptions. With the rise in LAI, a larger foliage of leaves lowers the radiation contact with bare ground surfaces and restricts soil evaporation, but increases the re-evaporation precipitation intercepted by leaves; thus, high greening can alter net evaporation. Due to intense greening, a significant increase in terrestrial evapotranspiration was observed, suggesting an intensified water exchange between land and atmosphere as feedback proportional to greening and browning as supported by various previous studies (Zhang et al., 2015; Zeng et al., 2017; Zeng et al., 2018b; Xu et al., 2020; Abbas et al., 2021; Kuhn and Butman, 2021). The overall phenomenon leads to a greening-induced ET enhancement that increases atmospheric water vapor content and promotes downwind precipitation.

The integration of WRF with the land model is widely used for accessing climate change due to anthropogenic effects; a few studies tried to address the effect of LULC on climate change. Paul et al. (2016) and Paul et al. (2018) have used the WRF model for simulating and downscaling summer monsoon precipitation in the Indian region and southern India for addressing the impact of LULCC on moisture supply especially during 1983–2011 and found that savannas and forest cover play an important role in moisture supply. Hua et al. (2015) reported the direct association of large regional responses to climate change with LULCC under RCP2.6 and RCP8.5 scenarios using the CanESM2 model and highlighted a cooling effect in the western United States and Eurasia due to the decrease in the needleleaf evergreen trees, while warming (by ∼0.1°C) in the tropics due to substantial LULC change. Cao et al. (2015) reported a widespread reduction of summer temperature (∼1°C) and extensive increase in winter temperature (∼0.8°) due to land use change during 2001–2010 in the agro-pastoral transitional zone of North China. Li X. et al. (2018) used the WRF model and reported a dramatic increase in temperature due to a decrease in the latent heat flux and an increase in sensible heat flux caused by urbanization, while an increase in precipitation primarily in the eastern part of Kolkata due to irrigation. Zhang et al. (2010) concluded that rapid conversion of croplands to urban leads to increased temperature and precipitation with decreased surface humidity, with a stronger influence in summer than in winter in China using WRF. Through incorporating Noah Urban Canopy with the WRF model, Grossman-Clarke et al. (2010) explored contributions of urban expansion to near-surface air temperature and concluded the rise in maximum temperature while the conversion of irrigated agricultural lands to developed areas. Patel et al. (2020) used the WRF model integrated with World Urban Data Analysis and Portal Tool (WUDAPT)-based Local Climate zones to simulate the heavy rain events and reported the influence of building height, building distribution, open surfaces, and other thermal parameters on precipitation at microscale in the city region. Similarly, the conversion of grassland to cropland and urban resulted in an increased air temperature (by ∼0.23°C) during summer, while afforestation leads to an increase in the air temperature (by ∼0.1°C) in the Lhasa River basin (Tibetan Plateau) during two runs of the model (1980 and 2015) because latent heat flux and sensible heat flux are highly sensitive to land cover changes, and other fluxes seem to be weaker (Li et al., 2020). A similar observation in Beijing as transformation from cropland to urban leads to increase in air temperature of ∼0.36°C during 1900–2010 as modeled using WRF (Li J. et al., 2018).