Impact of land use land cover on variation of urban heat island characteristics and surface energy fluxes using WRF and urban canopy model over metropolitan cities of India

The present study focuses on the assessment of urban heat island (UHI) and urban energy fluxes over Indian metropolitan cities. For this purpose, a high-resolution numerical Weather Research and Forecasting (WRF) model and a coupled single-layer urban canopy model (WRF-UCM) has been implemented over Delhi, Kolkata and Hyderabad to evaluate the performance of these models in capturing the impact of urban areas in simulating the UHI over different land use land cover (LULC) classes. The model initial conditions are obtained from NCEP-FNL reanalysis data. The model simulations reveal that the air temperature at 2 m height (T2) and UHI intensity are overestimated by WRF whereas the WRF-UCM simulations are matching well with ERA5-Land reanalysis data considered as observations. For the WRF-UCM simulation, the statistical errors in T2 and UHI intensity are within the desirable limits and the index of agreement (IoA) is more than 0.8 over the built-up areas (BAs). By comparing with the satellite observations, the relative surface UHI is noticed to be better captured over the BA by WRF-UCM simulation. The sensible heat flux simulated by WRF is highly overestimated whereas WRF-UCM is in good agreement with the satellite observation over the BA. Interestingly, the latent heat flux is reasonably well simulated over all the LULC classes. The results derived from the present study have shown the performance of WRF-UCM in simulating the spatial variation of UHI and fluxes of sensible and latent heat-related varied LULC classifications over metropolitan cities of India.


Introduction
Cities are characterized with advanced infrastructure, job opportunities, economical progress and improved life-style, and also offer many environmental distresses such as urban heat island (UHI), modification in precipitation and wind circulation, flash flooding, drought and frequent heat waves (Roth 2000;Arnfield 2003;Shepherd 2005;Masson 2006;Kanda 2007;Li and Bou-Zeid 2013;Ramamurthy and Bou-Zeid 2017). These alterations in the urban climate further affect the agriculture and public health and increase the risks of property loss and life as well. The urban surfaces with altered thermal, radiative and aerodynamic properties possess higher temperature in comparison to the surrounding rural areas and thus generate the phenomena of UHI (Roth 2013). The surface temperature and the turbulent heat fluxes of the urban facets (roofs and walls), roads, pavements and other urban surfaces play an important role in the urban micro-climate changes (Wang et al. 2011b).
Mesoscale modelling could be quite effective in studying the impact of the urbanisation on regional climatic changes when parameterized properly with urban conditions (Masson 2006;Champollion et al. 2009;Flagg and Taylor 2011). Cities have a complex geometry of irregular building heights, diversified street canyons with the presence of water bodies and vegetation within the built-up regions, which can significantly affect the atmospheric thermal structure and circulation pattern. The simplest way of implementing urban conditions in mesoscale model is to consider it as a plain surface with urban surface properties, the 'slab model' (Fortuniak et al. 2005). The slab model has some drawbacks as it fails to reproduce characteristic traits of wind circulation, turbulence and the energy balance, and thus, the 'urban canopy model' (UCM) is conceptualized (Kusaka et al. 2001; Kusaka and Kimura 2004;Coceal and Belcher 2004;Ryu et al. 2011). The UCMs consider the three-dimensional urban geometry and include the effects of shading, multiple reflection and trapping of radiation by the urban facets. The model also incorporates the skin temperature of urban facets, the temperature profile and heat flux of street canyon and the anthropogenic heat added to the sensible heat flux (Chen et al. 2011a). The UCMs can be sorted to single-layer UCM and multi-layer UCM. The single-layer UCM considers a single layer of urban canyon (Masson 2000;Kusaka et al. 2001;Harman et al. 2004;Kanda et al. 2005;Holt and Pullen 2007;Oleson et al. 2008), and the multi-layer UCM considers multiple layers of building heights (Brown 2000;Martilli 2002;Dupont et al. 2004;Kondo et al. 2005;Chin et al. 2005). The UCM has been extensively coupled with the Weather Research and Forecasting (WRF) model (Skamarock et al. 2005(Skamarock et al. , 2008Chen et al. 2011a, b) and examined against the ground observations (Hamdi and Schayes 2007;Lemonsu et al. 2009;Grimmond et al. 2010;Grimmond et al. 2011) which was found to be efficient for urban climate modelling (Wang et al. 2011a, b).
The UHI effect results due to the replacement of natural surface by the artificial surfaces of modified properties such as albedo, roughness length, permeability, thermal conductivity and heat capacity. For simulating the impact of UHI, the urban parameterisation and land use information is crucial in the modelling studies (Lin et al. 2016). In 2004, Kusaka and Kimura coupled single-layer UCM to 2-D mesoscale model to portray the urban impact on UHI and the results acknowledged the observational traits. Lee et al. (2016) analysed different surface parameterisations in the Noah land surface model (LSM), modified LSM and UCM for WRF simulation over the urban boundary layer and noticed the UCM simulation to be closer to the observations. Holt and Pullen (2007) utilized single-layer UCM and multi-layer UCM to analyse the atmospheric response to urban forcing. Miao et al. (2009) used single-layer UCM coupled to WRF to effectively generate the spatial distribution of UHI and the diurnal variation of the intensity along with other boundary layer characteristics. Studies also indicate that the single-layer UCM coupled to the WRF model can also yield monthly averaged UHI effectively (Miao et al. 2010;Zhang et al. 2010). In 2011, Ryu, Baik and Lee used single-layer UCM, considering two facing building walls (as sunlit or shaded walls), and effectively evaluated the surface temperature and the energy budget of individual wall. Chen et al. (2011a, b) used the WRF-UCM coupled system to successfully produce the night-time UHI intensities along with wind and sea breeze patterns. Kusaka et al. (2012) used the WRF-UCM coupled system to forecast the urban climate for August 2070 over the Japanese cities Tokyo, Osaka and Nagoya. The results predicted that the monthly averaged temperatures for August 2070 are about 2.3 higher than those of the year 2000 and are comparable to 2010 temperatures. Yang et al. (2012) used the WRF model coupled with UCM over Nanjing area and effectively reproduced the urban climate characteristics such as temperature, relative humidity, frequency of precipitation and energy balance components. Li and Bou-Zeid (2014) used an improvised UCM coupled with the WRF model that generated precise surface and potential temperature profiles and UHI intensities and also provided way to mitigation strategies. The study also indicated that surface temperature is more susceptible to roughness length parameterisation during daytime and to the boundary layer schemes during the night. Lin et al. (2016) used the WRF-UCM system, considering a 2-D urban fraction and including anthropogenic heat (WRF-UCM2D) for better simulation of temperature over urban areas than using the original WRF-UCM simulation. Few recent studies found are concentrated on the impact of vegetation in urban climate modelling and mitigating strategies for excess urban heating effect through UCM (Lee et al. 2016;Morais et al. 2016;Ryu et al. 2016;Wang et al. 2021).
Significant studies over worldwide cities established the effectiveness of the WRF and UCM modelling system in portraying the UHI phenomena, but only limited studies over Indian cities are noticed in recent decades. India being a rapid socio-economically evolving country has many developed or developing cities encountering UHI adversity; thus, it is essential to have profound studies using numerical simulations over these cities. The present study is an attempt to assess the performance of the WRF and WRF-UCM simulations in reproducing the UHI effect over the major Indian metropolitan cities.

Study area
The urban areas considered for the present study are three vastly growing Indian metropolitan cities: Delhi (28° 42′ 15″ N, 77° 06′ 09″ E), a northern inland city; Kolkata (22° 34′ 21″ N, 88° 21′ 50″ E), a coastal city in the east; and Hyderabad (17° 23′ 06″ N, 78° 29′ 12″ E), a southern peninsular city. Delhi is the national capital of India with more than 11 million population, spreads over an area of 1484 km 2 (Sultana and Satyanarayana 2018) and has humid subtropical and semi-arid climate as per the Köppen-Geiger climate classification system. The city experiences four seasonal periods in a year: summer, monsoon, post-monsoon and winter with a maximum temperature of 41-45 °C in summer (March to June) and a minimum temperature of 3-6 °C in winter season (December to February) (Das 1968;Bhati and Mohan 2018). Kolkata (formerly known as Calcutta) is one of the oldest Indian cities and serves as the capital of the state of West Bengal. The city accommodates about 4.5 million urban populations (Census of India reported in 2011) and spreads over 1257.7-km 2 suburban areas (Sultana and Satyanarayana 2018). The climate over the city is categorized as tropical wet and dry with the highest temperature of 40 °C in summer months (March to June) and the lowest temperature of 9 °C in winter months (December to February). Hyderabad is the capital of the state of Telangana, India, and is the fourth most populated city of India with a population of around 6.8 million (Census of India reported in 2011) within the area of 1136.2 km 2 (Sultana and Satyanarayana 2018). The city is characterized by multiple artificial lakes (Hussain Sagar, Osman Sagar) and water tanks or ponds. The climate is mostly tropical wet and dry and slightly semi-arid (Norman et al. 1995). The maximum temperature over the city during the summer months (March to June) often exceeds 40 °C, and the minimum temperature during the winter months (November to February) remains within 15-20 °C (https:// www. world weath eronl ine. com/ hyder abad-weath er-avera ges/ andhraprade sh/ in. aspx).

Data
In the present study, three Indian metropolitan cities are considered for analysing the performance of the WRF and WRF-UCM models in the simulation of UHI and turbulent energy fluxes. For this purpose, the National Centers for Environmental Prediction (NCEP) FNL (Final) Operational Global Analysis data of 0.25° × 0.25° resolution over Delhi during 11-13 January 2013, over Kolkata during 22-24 January 2017 and over Hyderabad during 1-3 January 2017 are used. The near-surface air temperature or the temperature at 2 m height, sensible heat flux (SHF) and latent heat flux (LHF) from an ERA5-Land reanalysis dataset of resolution 0.1° × 0.1° are used as the reference value for the analysis in the present study in the absence of high-resolution ground observation. The land surface temperature (LST) (Sultana and Satyanarayana 2018), SHF and LHF (Sultana and Satyanarayana 2022) derived from the Landsat 7 and Landsat 8 from the earlier studies are also used in the study for validating the simulated results. The nearsurface air temperature from the automatic weather station located in the airports of the study regions are also used in the present study for validating the simulated outcomes. For this purpose, the automatic weather stations located in the Indira Gandhi International (IGI) Airport (28.55° N, 77.09° E), Delhi; the Netaji Subhash Chandra Bose International Airport (22.64° N, 88.44° E), Kolkata; and the Begumpet Airport (17.45° N, 78.47° E), Hyderabad, are considered.

Model description
In the present study, the Weather Research Forecasting model (version 3.9.1) with Advanced Research WRF dynamics core developed by the National Center for Atmospheric Research (NCAR) was utilized (Skamarock et al. 2008) to assess the impact of urban canopy on UHI and the urban energy fluxes. The WRF model includes fully compressible, Eulerian and non-hydrostatic equations and a prognostic set of variables. In the present study, the model employs the hybrid vertical coordinate (HVC) and hydrostatic pressure coordinate as the vertical coordinate system and the Arakawa C-grid staggering horizontal coordinate system, respectively. The third-order Runge-Kutta scheme is used as the time integration scheme. Three interactive nested domains are configured for each city considered in the present study with the spatial resolution of 18 km for outer domain (d01), 6 km for first nested domain (d02) and 2 km for the inner most second nested domain (d03) containing the city region. Figure 1a-c shows the constructed domains for Delhi, Kolkata and Hyderabad, and the model details are specified in Table 1.
To include different physical processes, distinct physical parameterisation schemes are incorporated in the WRF modelling system and are also presented in Table 1 with the other specifications for the simulations. The model incorporates the Noah LSM to involve the surface temperature and exchange of the turbulent heat fluxes at the surface in the mesoscale model. The model uses the Lin scheme for micro-physics (Lin et al. 1983), Kain-Fritsch cumulus parameterisation (Kain 2004), rapid radiative transfer model (RRTM) scheme for long-wave radiation (Mlawer et al. 1997) and Dudhia scheme for shortwave radiation (Dudhia 1989). Along with Noah LSM, the Pleim-Xiu surface layer (Pleim 2006) and asymmetric convection model (ACM2) planetary boundary layer (Pleim 2007) are also incorporated. The options for physical parameterisations excluding the Noah LSM are chosen based on the earlier studies conducted over the Indian capital region Mohan 2016, 2018). To include the urban impact, the singlelayer UCM is coupled with the Noah LSM through the urban percentage parameter or the urban fraction, representing the proportion of the impervious surfaces in the WRF sub-grid scale. Thus, two sets of numerical simulations are conducted during the study period over the three cities under consideration: a simple WRF simulation without switching on the UCM and a coupled WRF-UCM simulation (single-layer UCM being switched on) (Kusaka et al. 2001;Kusaka and Kimura 2004;Chen et al. 2011b).
For each study region, the simulations were conducted for the duration of 72 h (3 days) for both the WRF and WRF-UCM simulations as Mohan and Sati (2016) suggested that the shorter-duration simulations (2-5 days) yield better results than those of longer duration. The initial 24 h of each of the simulation is considered as the spin-up, and the rest of the 48 h is utilized for the analysis. Over Delhi, the simulations (WRF and WRF-UCM) are conducted from 11 January 2013 0000 UTC to 14 January 2013 0000 UTC. Similarly, over Kolkata, the simulations are conducted from 22 January 2017 0000 UTC to 25 January 2017 0000 UTC and those over Hyderabad are conducted from 1 January 2017 0000 UTC to 4 January 2017 0000 UTC.

Evaluation methods
The present study attempts to assess the performance of the WRF model and the WRF-UCM coupled model in reproducing the UHI effect over the urban settlements. The near-surface air temperature and the corresponding UHI intensity simulated over the built-up areas (BAs), dry lands (DLs) and the vegetation are analysed separately to evaluate the performance of the simulations. The model-generated results are compared with the corresponding ERA5-Land reanalysis data as reference values in the absence of sufficient ground observations. The grid points over the city area and the surrounding regions are identified as the BA, the DL or the vegetation locations, and the temperature and the UHI intensity over these points from the simulations and the reanalysis dataset are examined separately. The UHI intensities are estimated over these grid points with respect to the selected rural location for each city under consideration. For Delhi, the rural location is selected to be Tolni village (28.24° N, 77.16° E) in the south of the city. For Kolkata, the rural site is the Maheswarbati village (22.9° N, 88.3°E) located in the north-west direction, and for Hyderabad, it is the Kagazmaddur village (17.7° N, 78.3° E) in the north side of the city. The hourly averaged temperature and the UHI intensity are obtained from both the simulations (WRF and WRF-UCM) and are compared with the corresponding ERA5 reanalysis data.
The simulated temperatures and the UHI intensities are statistically analysed with respect to the reanalysis data to determine the extent of error in these predicted values. The statistical errors like mean bias (MB), mean absolute error (MAE) and root mean square error (RMSE) indicate the deviation of the simulated values from the reference/ observed values, and the statistical parameter like index of agreement (IoA) shows the degree of proximity of predicted values to the reference/observed values. According to Emery et al. (2001), the maximum MAE for the simulated temperature is 2 °C and the standard IoA should be more than 0.8, and the present study follows the same benchmarks.
Two different modes of measurements, observations and model predictions, have an unequal range of values and, upon comparison, fail to capture common features in spatial distribution figures. Bhati and Mohan (2016) introduced the relative UHI, a dimensionless parameter that successfully reproduced the spatial distribution and the high UHI zones. The relative UHI at a point (i,j) can be defined as where max│UHI│ is the maximum of the absolute values of the UHIs of all the points in the domain. Relative UHI ranges between − 1 and 1; the higher positive values indicate UHI hotspots, and the higher negative values are the cold spots.
In the present study, the surface UHIs (SUHIs) obtained from the satellite observations from earlier studies (Sultana and Satyanarayana 2018) are used for comparing with the simulated results. Thus, the relative SUHIs from both the simulated outputs and satellite observations are obtained through the following formulation. The relative SUHI at the point (i,j) is given by where max│SUHI│ is the maximum of the absolute values of the SUHIs of all the points in the domain. The SUHIs from the satellite observation are obtained from the LST retrieved from the infrared thermal bands of Landsat 7/8 imageries (Sultana and Satyanarayana 2018), and for the simulations, the SUHI is estimated using the simulated surface skin temperatures. (1)

Diurnal variation
The simulated air temperatures at 2 m height (T 2 ) over different grid points are averaged for different LULC classes and are plotted separately. The averaged T 2 values are compared with the corresponding averaged ERA5-Land reanalysis datasets (ERA5 now onwards) over different LULC classes (BA, DL and vegetation) for the study areas. The UHI intensities at 2 m height estimated with respect to the rural locations over these LULC classes are also compared with that estimated from ERA5.
Delhi Figure 2a-c shows the variation in the T 2 over the BA, the DL and the vegetation of Delhi as estimated from WRF and WRF-UCM simulations and the ERA5 dataset. Figure 2a indicates that the WRF-UCM-simulated T 2 over the BA is mostly matching with the ERA5 values and the WRF-simulated T 2 is higher than these values, thus found to be overestimating. Both the simulated T 2 over the DL and the vegetation are noticed to be matching well during the study period. The simulated T 2 values over the DL are found to slightly overestimate the ERA5 temperatures ( Fig. 2b) throughout the study period, whereas over the vegetated areas, the simulated values are overestimated during the daytime and underestimated during the evening and night-time (Fig. 2c). The UHI intensities estimated over the different LULC classes of Delhi are shown in Fig. 2d-f. The UHI intensity over the BA from WRF simulation is significantly higher than the ERA5 value whereas the WRF-UCM value is noticed to be closer to the ERA5 value (Fig. 2d). The maximum UHI intensity (UHI max ) over the BA is 8.5 °C for WRF simulation, 7.4 °C for WRF-UCM simulation and 6.6 °C for ERA5 dataset. Over the DL and the vegetation, both the simulated values are noticed to be closer to each other. Over DL, the simulated UHI intensities are noticed to be comparable with the ERA5 values (Fig. 2e). The UHI max over the DL are found to be 7.7 °C, 7.4 °C and 6.4 °C for WRF, WRF-UCM and ERA5, respectively. Over the vegetated areas, the simulated UHI intensities are noticed to be lower in comparison to the reanalysis values throughout the study period (Fig. 2f). The UHI max over the vegetated areas for WRF, WRF-UCM and ERA5 is found to be 5.2 °C, 5.6 °C and 6.3 °C, respectively.
Kolkata Figure 3a-c shows the simulated T 2 and corresponding ERA5 temperature over the BA, the DL and the vegetation of Kolkata, and the estimated UHI intensities over these classes are shown in Fig. 3d-f. Over the BA, the WRF-simulated T 2 is mostly higher than the ERA5 values whereas the WRF-UCM-simulated values are noticed to be matching with the ERA5 values (Fig. 3a). Figure 3b indicates that for the DL, both the simulated T 2 values are noticed to be reproducing the ERA5 values in the morning and the afternoon hours, but during the night-time and the early morning, WRF-simulated T 2 is noticed to be slightly overestimated. Over the vegetated areas as shown in Fig. 3c, both the simulated T 2 values are closer to each other and to the ERA5 values during the morning and the afternoon hours. During the night-time and the early morning hours, the simulated T 2 values are noticed to be underestimated in comparison to the ERA5derived T 2 values.
The WRF-simulated UHI intensity over the BA is found to be overestimating with respect to the ERA5 values whereas the WRF-UCM-simulated UHI intensity is noticed to be matching with the ERA5 values during the evening and night and with the overestimated WRF-simulated values during the daytime as shown in Fig. 3d  Hyderabad The comparison of diurnal variation of simulated T 2 and the UHI intensities with the ERA5 values over the BA, the DL and the vegetation of Hyderabad is shown in Fig. 4. Over the BA as shown in Fig. 4a, the WRF-UCMsimulated T 2 is noticed to be matching with the ERA5 T 2 and the WRF-simulated T 2 is slightly overestimated during the study period. Over the DL, the simulated T 2 was noticed to be matching with each other and remaining close to the ERA5 temperatures during the daytime. During the night-time, both the simulated T 2 values are slightly underestimated with respect to the ERA5 values (Fig. 4b).
Over the vegetation, the simulated T 2 values are noticed to be completely matching to each other and to the ERA5 values during the daytime, but during the night-time, the simulated values are significantly underestimated (Fig. 4c).
The simulated UHI intensities over the BA are noticed to be matching the ERA5 values during the day but overestimating during the late night and the early morning (Fig. 4d). The UHI max over the BA of Hyderabad is 6.68 °C from WRF simulation, 6.08 °C from WRF-UCM simulation and 5.58 °C from ERA5 reanalysis data. Figure 4e indicates that over the DL, both the simulated UHI intensities are slightly overestimated with respect to ERA5 values. The UHI max over the DL for WRF, WRF-UCM and ERA5 is 4.43 °C, 4.27 °C and 3.76 °C, respectively. Over the vegetation, both the simulated UHI intensities match to each other and are close to the ERA5 values during night whereas they are underestimated during the afternoon (Fig. 4f). The UHI max over the vegetation is 2.5 °C from WRF simulation, 2.54 °C from WRF-UCM simulation and 2.45 °C from ERA5 dataset.

Statistical assessment
The statistical assessment for the simulated T 2 and UHI intensity is conducted with respect to the ERA5 reanalysis dataset for different LULC classes separately for the study regions. According to Emery et al. (2001), the maximum limit for the MAE is 2 °C or less than 10% of the estimated values and the desirable value for IoA is 0.8 or more than that. The following assessments are done focusing on the same benchmark.

MB, MAE, RMSE and IoA for T 2 and UHI intensity
The simulated T 2 and UHI intensity over Delhi, Kolkata and Hyderabad are statistically assessed with respect to the ERA5 reanalysis data and presented in Table 2, Table 3 and  Table 4, respectively. The MB in T 2 and UHI intensity over Delhi is noticed to be higher in comparison to that over Kolkata and Hyderabad. The MAE in T 2 over Delhi (BA and vegetation) is 2 °C or more for WRF simulation, whereas for WRF-UCM simulation, the MAE is lower and within the limits (Table 2). Over Kolkata (Table 3) and Hyderabad (Table 4), the MAE in T 2 for both the simulations is well within the desired value of 2 °C and the MAE for the WRF-UCM simulation is found to be lower than that for the WRF simulation for all LULC classes. The percentage of MAE in T 2 for the WRF-UCM simulation is found to be lower than that for the WRF simulation and well within the limits of 10%. The RMSE for the WRF-simulated T 2 is high, and the corresponding WRF-UCM simulation values are comparatively less for all the study regions. The IoA for T 2 for the WRF-UCM simulation over all the LULC classes is higher than 0.8 for all the study regions, but for the WRF simulation, the IoA is occasionally less than 0.8. The MAE and the RMSE in the WRF-simulated UHI intensity over the BA of Delhi are significantly higher than the errors noticed in the WRF-UCM simulation (Table 2), but over the DL and the vegetation, the errors in both the simulations are close to each other. The MAE and the RMSE in the UHI intensity for the WRF simulation are higher than those for the WRF-UCM simulation for all the LULC classes over Kolkata and Hyderabad. The MAE percentages over the study regions for the WRF-UCM simulation are noticed to be less than those for the WRF simulation and are within the limits of 10%. The IoA for UHI intensity from both the simulations is noticed to be higher than the desirable 0.8 value over all the LULC classes of the study regions.

Time series of MAE and RMSE for T 2 and UHI intensity
The MAE and the RMSE in the simulated T 2 and the UHI intensity are estimated for each hour during the study period, and the hourly averaged values are plotted to analyse the variation in these errors with respect to time. Figure 5 depicts the averaged time series of the MAE and the RMSE in the simulated T 2 over the different LULC classes of the study regions, and Fig. 6 depicts the same for the estimated errors in UHI intensity. The figures suggest that the error in the WRF simulation is generally higher than that in the WRF-UCM simulation except for few occasions. Figure 5 suggests that for Delhi, the error values for both the simulations are lower during the late night and early morning hours over the BA and the DL, but for the vegetated areas, the error values are higher during night-time than during daytime. For Kolkata, the errors over all the LULC classes are mostly higher during the early morning hours and night-time. For Hyderabad, over the DL and the vegetation, the errors are higher for both the simulations during the early morning hours and night-time than during the daytime (Fig. 5), whereas over the BA, the change in the errors with time from both the simulations is irregular. The averaged time series of the MAE and RMSE in the simulated UHI intensities over all the LULC classes of Delhi are lower during the afternoon hours and significantly higher during the night and early morning hours. For Kolkata and Hyderabad, over the BA and the DL, the errors from both the simulations are lower during the afternoon hours and higher during the night and early morning hours. Over vegetation, the errors are uniformly distributed throughout the duration for Kolkata, and for Hyderabad, the errors are higher during daytime in comparison to evening and night-time. Figure 7 shows the spatial distribution of the averaged relative UHI over Delhi estimated from WRF-UCM simulation, for different hours of the day. The upper row (Fig. 7a-d) shows the spatial distribution for the day hours, and the lower row ( Fig. 7e-h) shows the spatial distribution for the night hours. Box 1 represents a dense BA, and Box 2 represents a BA with the presence of vegetation. The figure indicates that over the dense BA (Box 1), relative UHI is about 0-0.2 during 09:30 LT, increases to 0.4 during 12:30 LT and gradually increases up to 0.8 during late afternoon (14:30-16:30 LT). The relative UHI is noticed to have exceeded 0.8 during 21:30 LT and decreased slowly to 0-0.4 and negative values during late night hours. On the other hand, over Box 2, the relative UHI remains negative until late afternoon and the value increases to 0.4 during 16:30 LT and keeps increasing to more than 0.8 at late night hours (01:30 LT). During early morning hours (04:30 LT), the relative UHI was noticed to have decreased to 0.4 and kept decreasing. Over the vegetation (Box 3), the relative UHI values are noticed to be highly negative up to − 0.8 during the daytime, but during the night-time, the values increase to a lower negative value of − 0.2 until early morning (04:30 LT), and afterwards, the value is noticed to be reducing again. Over the DL (Box 4), the relative UHI is positive throughout the day and night hours, of lower values (0-0.4) during afternoon and of high values up to 0.8 during the late night-time. The hotspots are noticed during the afternoon over the dense BA and during the late night over the DL and sparse BA with vegetation.

Kolkata
The spatial distribution of the averaged relative UHI over Kolkata for different hours of the day estimated from WRF-UCM simulation is shown in Fig. 8. The day hours are presented in the upper row (Fig. 8a-d), and the night hours are in the lower row ( Fig. 8e-h). Box 1 represents a dense BA, and Box 2 represents a lesser-dense BA. During 09:30 LT, the relative UHI is negative over the BA, (− 0.6) to 0 in Box 1 and (− 0.2) to 0 over Box 2. In the afternoon (14:30 LT), the relative UHI over the BA gradually increases to a positive value of 0.4 and further increases to more than 0.8 during 17:30 LT. During night, the relative UHI values were noticed to have remained within 0.4-0.8 in both the boxes with few regions exceeding 0.8 in Box 1, and in the early morning hours (05:30 LT), they were noticed to be slowly decreasing (0.2-0.6). Box 3 represents a dense vegetated region, and the relative UHI is mostly negative throughout the day. During the morning and afternoon hours, the relative UHI value is around − 0.2 to 0 and reduces to − 0.4 during the evening (17:30 LT). The relative UHI value reduces to a further negative value of around − 0.8 during late night hours and is slightly increased to − 0.4 during early morning (05:30 LT). Over the DL (Box 4), the relative UHI is 0.2-0.6 during 09:30 LT and slightly reduces to 0.4 during afternoon hours. During the evening hours (17:30 LT) and night-time (21:30 LT-02:30 LT), the relative UHI varies between − 0.2 and 0.2 and reduces up to − 0.4 during the early morning hours (05:30 LT).

Hyderabad
The spatial distribution of the averaged relative UHI derived from WRF-UCM simulation for different hours of day is shown in Fig. 9. The upper row (Fig. 9a-d) shows the relative UHI during day hours, and the lower row ( Fig. 9e-h) shows that during night hours. Box 1 represents a section of dense BA, and Box 2 represents a less-dense BA with the presence of vegetation. Over the BA for both the boxes during 09:30 LT, the relative UHI is highly negative (− 0.8 to − 0.4). Over Box 1, the relative UHI is gradually increased to − 0.2 to 0 during the afternoon hours (12:30-14:30 LT) and to 0.2-0.6 during 16:30 LT. On the other hand, over Box 2, the relative UHI is slowly increasing during the afternoon hours and becomes − 0.2 to 0.2 during 16:30 LT. During the night-time, the relative UHI is highly positive (0.4 to 0.8 or more) over the BA and is slightly reduced to 0.2-0.6 in the early morning (04:30 LT). Over the dense vegetation (Box 3), the relative UHI is negative during the daytime but the vegetation close to the BA has positive relative UHI values during the night. Over the DL (Box 4), the

Validation of simulated results with the satellite-derived images
In this section, the spatial distributions of the simulated UHI, SHF and LHF are compared with the results of our earlier studies (Sultana andSatyanarayana 2018, 2022) derived from the Landsat 7 or Landsat 8 imageries of the study regions. The satellite images deliver land surface temperature instead of T 2 ; thus, the relative SUHI is estimated using the simulated surface skin temperature and is compared with the relative SUHI derived from satellite imageries.

Relative surface UHI
Delhi Figure 10a-c presents the relative SUHI over Delhi estimated from Landsat 7 imageries and WRF and WRF-UCM simulations. The satellite image is available for 13 January 2013 during 10:45 LT, and hence, the simulated images are retrieved during 10:30 LT of the same date for comparison. In the figure, the sparse BA close to vegetation is represented by Box 1 and the dense BA is by Box 2. Over Box 1 and Box 2, the satellite-derived relative SUHI is 0-0.1 with pockets of higher values and the WRF-UCM-simulated relative SUHI is 0-0.2 whereas the WRF-simulated value is found to be negative (− 0.2 to 0), thus underestimated. The satellite-derived relative SUHI over the DL identified in Box 3 is 0.2 or more, and both the simulated values range between 0.2 and 0.5, matching the satellite-derived results. Over the vegetation (Box 4), the satellite-derived relative SUHI is noticed to be generally negative (− 0.3 to 0), and for both the simulations, the relative SUHI value lies between − 0.4 and − 0.1. In the satellite-derived relative SUHI, few hotspots are identified with a value more than 0.6 over the dense BA, the DL and the airport area. In the simulated results, the hotspots are only identified around the airport and few dense BA outside the city.

Kolkata
The spatial distributions of relative SUHI over Kolkata estimated from Landsat 8 imageries, WRF simulation and WRF-UCM simulation are shown in Fig. 10d-   Hyderabad Figure 10g-i shows the spatial distributions of the relative SUHI estimated from Landsat 8 imageries and from WRF and WRF-UCM simulations for Hyderabad. The Landsat 8 image on 2 January 2017 is available at 10:40 LT, and thus, the simulated results of 10:30 LT of the same date are considered for comparison. The satellite-derived relative SUHI over the dense BA (Box 1) lies between − 0.1 and 0.2, and over the less dense BA (Box 2), it lies in negative value (− 0.2 to 0) with pockets of small positive values. For WRF simulation, the value lies between − 0.3 and − 0.1 over Box 1 and between − 0.5 and − 0.3 over Box 2. For WRF-UCM simulation, the relative SUHI value is in the range − 0.4 to − 0.2 over Box 1 and − 0.6 to − 0.4 over Box 2. The simulated relative SUHI over the BA appears to be slightly underestimated with respect to the satellitederived results. Over the DL (Box 3), the satellite-derived relative SUHI is 0-0.4 and the simulated relative SUHI is 0.2-0.5 for both the simulations. Over a small vegetation section (Box 4), the satellite-derived relative SUHI is − 0.5 to 0. Over Box 4, the relative SUHI is noticed to be positive (0.1-0.4) for WRF simulation, whereas the value lies between − 0.1 and 0.2 for WRF-UCM simulation. The simulations appear to have failed to capture the relative SUHI over the vegetation.

Sensible heat flux
Delhi Figure 11a-c depicts the spatial distributions of SHF derived from satellite imageries and the WRF and WRF-UCM simulations on 13 January 2013. Over the BA (Boxes 1 and 2), the SHF is noticed to have the value 110-140 W m −2 with patches of higher values (> 140 W m −2 ) for the satellitederived results. The WRF-UCM-simulated SHF values are in the range of 100-120 W m −2 and, thus, are comparable to the satellite-derived values. On the other hand, the WRFsimulated SHF possesses values higher than 140 W m −2 , which are highly overestimated. Over the DL (Box 3), the satellite-derived SHF has the value range 120-150 W m −2 , and for both the simulations, SHF values ranged from 130 to 160 W m −2 , which are comparable to the satellite-derived results. Over the vegetation areas (Box 4), the simulated SHF values ranged from 110 to 140 W m −2 , and the satellitederived SHF values ranged from 90 to 120 W m −2 .
Kolkata The spatial distributions of the SHF over Kolkata on 24 January 2017, derived from the satellite imageries and the WRF and WRF-UCM simulations, are shown in Fig. 11d-f. The SHF over the BA in Box 1 is around 120-160 W m −2 with small pockets of higher values, and that in Box 2 is 120-140 W m −2 with smaller value patches for the satellite imageries. The WRF-UCM-simulated SHF value over the BA ranges from 110 to 140 W m −2 (Boxes 1 and 2), which is comparable to the satellite-derived values. The WRF-simulated SHF, on the other hand, is 150-200 W m −2 , which is significantly overestimated with respect to the satellite-derived SHF. Over the DL (Box 3), the WRF-simulated SHF is 130-150 W m −2 and the WRF-UCM-simulated SHF is 140-160 W m −2 , whereas the satellite-derived SHF is around 100-120 W m −2 . Thus, the simulated SHF over the DL is slightly overestimated for both the simulations. The SHF over the vegetated areas (Box 4) is estimated to be 100-120 W m −2 from satellite imageries and interestingly matching to the simulated values as well.
Hyderabad Figure 11g-i shows the spatial distributions of the SHF over Hyderabad for 2 January 2017, estimated from the satellite imageries and the WRF and WRF-UCM simulations. The SHF estimated from the satellite images over the BA (Boxes 1 and 2)

Latent heat flux
Delhi The spatial distributions of LHF derived through the satellite observation and the WRF and WRF-UCM simulations are depicted in Fig. 12a-c. The figure suggests that both the simulated LHFs have almost similar spatial distribution and are also matching the satellite-derived distribution. Over the BA (Boxes 1 and 2), the satellitederived LHF values are noticed to be ranging between 20 and 60 W m −2 , and the simulated LHF values are found to be 10-50 W m −2 over the BA. Over the DL (Box 3), for the satellite-derived results, the LHF ranged from 20 to 50 W m −2 , and the simulated LHF values ranged from 30 to 60 W m −2 for both the simulations. Over the vegetation (Box 4), the LHF values from satellite images are mostly more than 100 W m −2 and the simulated LHF values are ranging from 90 to 140 W m −2 .
Kolkata Figure 12d-f shows the spatial distributions of the LHF over Kolkata estimated from the WRF and WRF-UCM simulations and from the satellite imageries. Over the BA, the satellite-observed LHF has a value less than 40 W m −2 with small pockets of higher values of 40-60 W m −2 (Boxes 1 and 2). For WRF simulation, the LHF is less than 20 W m −2 , and for WRF-UCM, the value is 20-40 W m −2 with a pocket of values between 0 and 20 W m −2 . Thus, the WRF-simulated LHF is comparable to the satellite-derived values, whereas WRF-UCM-simulated LHF is mostly matching with the values. Over the DL (Box 3), the satellitederived LHF is 20-40 W m −2 , but the simulated LHF values are highly overestimated as the value lies between 100 and 120 W m −2 . Over the vegetation, the satellite-derived LHF is

Diurnal variations over airport
In this section, the diurnal variations of the simulated T 2 , UHI intensity and urban energy balance components over the airports of the study regions are compared with the corresponding ERA5 dataset. The automatic weather stations located at the airports provide the ground observation for the T 2 , and thus, the observational values are also included in the comparison plots for T 2 . Figure 13 shows the variations of the T 2 , UHI intensity and urban energy balance components over the Delhi airport (IGI Airport) as estimated through simulation are compared with the corresponding ERA5 dataset. In Fig. 13a, the T 2 values from ERA5 reanalysis dataset and from the WRF and WRF-UCM simulations are plotted along with the ground observations available at the airport. The figure indicates that the WRF-UCM-simulated values are closer to the ERA5 values and/or the ground observation data, whereas the WRF-simulated temperature values are noticed to be overestimating throughout the study period. The UHI values estimated from simulations with respect to the rural location (Tolni village) are compared with the corresponding ERA5 values in Fig. 13b and are noticed to be highest during the night-time and small and negative during the afternoon hours. Both the simulated UHI intensities are observed to be overestimated with respect to the ERA5 dataset during evening to early morning hours, and the WRF-UCM-simulated values are noticed to be closer to the ERA5 values. Figure 13c shows the diurnal variation in the SHF derived from the simulations and the ERA5 dataset during the study Fig. 13 Diurnal variation of the a near-surface air temperature, b near-surface UHI intensity, c sensible heat flux, d latent heat flux and e surface energy fluxes over the Delhi airport, during 12 January 2013 0000 UTC to 14 January 2013 0000 UTC period. The WRF-simulated SHF is noticed to have matched the ERA5-derived SHF during the afternoon hours with peak values (~ 200 W m −2 ), and the WRF-UCM-simulated SHF is observed to have higher values around 250 W m −2 . Over a BA or urban surface (airport) of high heat conductivity, the SHF is expected to be as high as 300-400 W m −2 (Holt and Pullen 2007), and thus, the ERA5 and WRF-derived SHF is assumed to be slightly underestimated. Figure 13d shows the LHF variation as derived from the ERA5 dataset and the simulations during the study period. The WRF-simulated LHF is noticed to have peaked during afternoon with a value around 130 W m −2 whereas the WRF-UCM simulation peaks with the value around 30 W m −2 . The ERA5-derived LHF lies in between with a peak value around 75 W m −2 . Over the urban surface (airport), the evapotranspiration is expected to be low and, hence, the lower values of LHF. The low value of WRF-UCM and ERA5-derived LHF can be considered as satisfying the expectations, and the WRFsimulated LHF is assumed to be heavily overestimated. During the night, the turbulent heat flux values are simulated to be very small and found to have matched to the ERA5 values or slightly overestimated. Figure 13e depicts the variation of the urban energy balance components as derived by the WRF and WRF-UCM simulations for the study period. The net radiation flux (Q*) for both the simulations is observed to be same. The SHF (Q H ) is noticed to be slightly underestimated, whereas the LHF (Q E ) is found to be highly overestimated for the WRF simulation. Figure 14 shows the variations of T 2 , the UHI intensity and the surface turbulent heat fluxes from the simulations and the ERA5 data along with the comparison of the urban energy balance components from two simulations during the study period. The T 2 over Kolkata estimated from the WRF and WRF-UCM simulations is compared with the ERA5 reanalysis data and the ground observations collected from the airport in Fig. 14a. The figure suggests that the WRF-UCMsimulated temperature is closer to the ground observation value and the ERA5 data, but the WRF-simulated value is found to be mostly overestimated. Figure 14b shows the comparison of the simulated UHI intensity with the ERA5 data. The WRF-UCM-simulated UHI intensity is closer to the ERA5 values, and the WRF-simulated UHI intensity is noticed to be higher than these values, hence overestimated. Figure 14c shows the comparison of SHF over the airport derived from WRF and WRF-UCM simulations with the ERA5 data during 23 January 2017 and 24 January 2017. The WRF-simulated SHF and the ERA5-derived SHF reach the peak values of ~ 250 W m −2 during 11:30-13:30 LT, whereas the WRF-UCM-simulated SHF reaches the peak value of ~ 300 W m −2 around 12:30-14:30 LT. The LHF over Kolkata airport derived from ERA5 dataset and the WRF and WRF-UCM simulations for the study period are compared in Fig. 14d. The LHF from the WRF-UCM simulation slowly increases and peaks to around 50 W m −2 during 11:30-13:30 LT. The ERA5 data and the WRF-simulated value, on the other hand, increase to peak values around 140-160 W m −2 during 11:30-13:30 LT. During the night, the WRF-simulated LHF and the WRF-UCM-simulated LHF are negligible, but the ERA5-derived LHF has a significant positive value. The diurnal variations in urban energy balance components over the Kolkata airport, simulated from WRF and WRF-UCM simulations for the study period, are compared in Fig. 14e. The net radiation flux (Q*) is observed to be similar for both the simulations. The SHF (Q H ) is little underestimated, and the LHF (Q E ) is highly overestimated for the WRF simulation.

Hyderabad
The variations in T 2 and the UHI along with the urban energy balance components over the airport during the study period are compared in Fig. 15. The T 2 estimated from the WRF and WRF-UCM simulations is compared with the ERA5 reanalysis data and the ground observations collected from the airport in Fig. 15a. The figure indicates that the WRF-UCM-simulated T 2 is significantly matching with the observation data and overestimating the ERA5 values except for morning hours. The WRF-simulated T 2 , on the other hand, is noticed to be mostly overestimating the ground observation and ERA5 data. The simulated UHI at 2 m height is compared with the ERA5 reanalysis data in Fig. 15b. The WRF-UCM-simulated UHI intensity is noticed to be closer to the ERA5 values, and the WRF-simulated UHI intensity is noticed to be overestimated throughout the study period.
The simulated SHF over the airport is compared with the ERA5 reanalysis data for 2 January 2017 and 3 January 2017 in Fig. 15c. The WRF-simulated SHF has the peak value of ~ 320 W m −2 , whereas the WRF-UCM-simulated SHF reaches the peak values of ~ 380 W m −2 during 12:30-14:30 LT. The ERA5-derived SHF, on the other hand, reaches the peak value of ~ 270 W m −2 during 12:30-14:30 LT and slowly decreases to the negative values during 17:30-18:30 LT. The simulated LHF over the Hyderabad airport and the corresponding ERA5 data are compared in Fig. 15d. The WRF-UCM-simulated LHF peaks at around 40 W m −2 during 12:30-13:30 LT, and the WRF-simulated LHF reaches the peak value of ~ 90 W m −2 during 12:30-14:30 LT. The LHF from ERA5 reanalysis data, on the other hand, have a peak value of ~ 140 W m −2 . The simulated urban energy balance components over the Hyderabad airport, for the study period, are compared in Fig. 15e. Both the simulated net radiation fluxes (Q*) are matching to each other during the study period. The WRF-simulated SHF (Q H ) is little underestimated, and the LHF (Q E ) is overestimated.

Limitations and suggested improvements of the study
The results of the present study indicate that the implementation of the single-layer urban parameterisation scheme in the mesoscale WRF model has reasonably simulated the near-surface air temperature, the corresponding UHI intensity and the surface energy fluxes over the study areas compared to WRF alone. The limitation of the WRF simulation is non-inclusion of effects of urban fraction, which has an influence on the temperature and moisture transport from the surface. Even though the single-layer UCM has given improved simulation of temperature and energy fluxes, it is essential to include multi-layer UCM parameterisation schemes for improved results. The multi-layer UCM can better define the urban geometry than the single-layer UCM The better representation of the LULC distribution, especially urban canopy, is essential for better representing the energy exchange mechanism in the model. In the present study, we have used the coarse-resolution USGS land use data in the WRF model. It is essential to use improved multi-satellite land use datasets with high spatial resolution over the study region. The urban land use maps with a more number of urban surface categories can be more effective in modelling the urban effects. The urban parameters, such as the building heights, urban fraction, roughness length and drag coefficient of the urban facades, and the radiative, thermal and moisture characteristics of the building materials are required to be appropriately defined for better presentation of the urban surface and the energy exchange mechanism in the model.

Conclusion
The present study over the Indian cities, Delhi, Kolkata and Hyderabad, suggests that the WRF-UCM coupled system has significantly reproduced the diurnal variation of air temperature (T 2 ) and the corresponding UHI intensities over the BA whereas the WRF simulation results are found to be overestimating. The MAE in T 2 and UHI intensity from the WRF-UCM simulation is found to be within the desirable limits over all the LULC classes for all the cities. The errors from the WRF simulations, on the other hand, are noticed to be higher than those from the WRF-UCM. The percentage of the MAE is less than 10% for the WRF-UCM-simulated  Fig. 13, but over Kolkata airport, during 23 January 2017 0000 UTC to 25 January 2017 0000 UTC T 2 and UHI intensities over all the LULC classes for the study regions with few exceptions for vegetation over Delhi and Kolkata. The IoA for T 2 and UHI intensity is noticed to be more than 0.8 for WRF-UCM simulations over all study regions, but for WRF simulation over Delhi and Kolkata, the IoA for T 2 is noticed to be slightly less over DL and vegetation. Over Hyderabad, both the WRF and WRF-UCM simulations have the statistical parameter values within the limits.
The spatial distribution of the relative UHI over the study regions derived from WRF-UCM simulations indicates that over the BA, the UHI intensity is high during the afternoon and the evening hours, and the intensity reduces during the late night and early morning hours. Over the vegetation areas, the UHI is mostly low and negative throughout the day and night. Over the DL, the UHI intensity is high and positive during afternoon and is reduced during the night-time. The relative SUHI spatial distributions estimated from WRF-UCM simulation over the BA are noticed to reasonably reproduce the satellite-derived distributions, but the WRF-simulated results are noticed to be slightly underestimated. Over the DL and the vegetation areas, the relative SUHI derived from both the simulations is comparable to satellite observations. The spatial distribution of SHF derived from WRF-UCM simulation over the BA is slightly underestimated, whereas from the WRF simulation, it is highly overestimated with respect to the satellitederived results for all the three cities under consideration. On the other hand, the spatial distributions of LHF derived from both the simulations are significantly comparable to the satellite-derived results for all the LULC classes of the study regions except over the DL of Kolkata and Hyderabad.
The comparison of the diurnal variations of the simulated temperature at 2 m height and the corresponding UHI intensity with the ERA5 reanalysis data over the airports as BA indicates that the WRF-UCM simulation results are closer to the observed and ERA5 values than those of the WRF simulations. The diurnal variation of the SHF derived from WRF simulation is close to the ERA5 data which is lower than the expected values, and the WRF-UCM-simulated values are noticed to be higher and close to the desired values over a BA. The diurnal variation of LHF from WRF-UCM simulation is found to be lower over the airports as expected over a BA whereas the WRF-simulated LHF and the ERA5derived LHF show higher values.
The present study implemented a single-layer UCM coupled with the WRF model and utilized the USGS land use dataset to assess the effect of urban canopy on UHI and the energy fluxes of the urban surface. The results indicate that  Fig. 13, but over Hyderabad airport, during 2 January 2017 0000 UTC to 4 January 2017 0000 UTC the WRF-UCM simulation performed well in comparison to the WRF simulations. In the present study, the results obtained are compared with the ERA5 reanalysis data in the absence of high-resolution ground observations. The overall study gave confidence that single-canopy layer UCM coupled with WRF would be the best option for understanding the spatial as well as diurnal variability of UHI and surface energy fluxes. The performance of mesoscale models simulating the urbanisation impacts can be enhanced by high-resolution LULC maps with a more number of urban categories. The urban surfaces can be better represented through advanced parameters, and the urban geometry can be better captured through multi-layer models with improved parameterisation.