Impact of urbanization on surface energy balance components over metropolitan cities of India during 2000–2018 winter seasons

The present study attempts to evaluate the urban energy balance components concerning increasing urbanization and artificial surfaces over Indian metropolitan cities during the 2000–2018 winter seasons by using Landsat 7 and 8 satellite imageries. The results indicate that the estimated ranges of the energy fluxes are in the typical values reported in the earlier literature over global cities. The sensible heat flux (SHF) increased considerably, and the latent heat flux (LHF) slightly decreased during the study period. The mean SHF over the built-up areas (BA) and the dry lands (DL) of Delhi record a maximum increase of 28.2 Wm−2and 39.7 Wm−2 during the study period. The inland cities have high values of SHF over DL than the coastal cities, and the LHF is high over all the land use classes for the west coast cities. The SHF (LHF) shows a positive (negative) correlation with the land surface temperature. The SHF (LHF) is about 19–33% (1.9–15%) of the net radiation flux, and the residual heat flux is about 60 to 80% of the net radiation flux. The present study advocates that the substantial changes of the surface energy balance parameters have a profound influence on the energy exchange mechanism and could affect regional climatic change.


Introduction
The explosive population rise and fleeting industrialization in the developing countries of tropical and subtropical climate have led to brisk and abrupt expansion of the artificial surfaces and the warmer environment over the cities, the phenomena of urban heat island (UHI) (Oke 1982;Voogt and Oke 2003;Sobrino et al. 2012;Roth 2013). The impact of the expanding cities is not just limited to the city scale but also influences the regional and the global climate through radiation balance and greenhouse gas emission (Roth 2007;Grimmond 2007;Seto et al. 2012;Kuang et al. 2015). The higher temperature experienced due to the UHI effect over the cities is generally the result of the declining vegetation cover, replacement of natural surfaces by artificial materials, anthropogenic heating, and the increasing pollution in the atmosphere (Grimmond 2007;Grimmond et al. 2016). The artificial materials (such as asphalt, concrete, marbles, tiles) of high heat capacities and conductivities store more heat and increase the surface temperature (Yang et al. 2019), whereas reduction in the surface cooling through evapotranspiration is due to the loss of the vegetation cover Yamaguchi 2005, 2007;Kuang et al. 2015). The artificial impervious surfaces have significant contributions in altering the surface energy fluxes and could influence the regional climate dynamics. It is essential to analyze the surface heat balance over the urban and rural domain to understand the boundary layer dynamics over the fast developing cities (Grimmond 1992;Hanna et al. 2011;Kotthaus and Grimmond 2014) and quantitatively investigate the influence of heat fluxes on the temperature rise associated with UHI (Kato et al. 2008).
The concept of surface energy balance was utilized in studying the surface dynamics during the end of the nineteenth century, and much significant work has been spotted regarding urban climatology in terms of the heat balance involving observation data and numerical modeling during the second half of the twentieth century (Sundborg 1951;Munn 1966;Tag 1968;Myrup 1969;Probald 1971;Oke et al. 1972). These studies indicated that the observational results were coarse and incomplete, and the model results were rather simple. The first study over tropical or subtropical cities related to energy balance was conducted for Mexico City (Oke et al. 1992), and next couple of decades, a few selective studies related to urban energy balance (UEB) are noticed (Goldreich 1992;Salmond 1999Salmond , 2005Roth 2007). At the beginning of the twenty-first century, many scientific researchers focused on deducing the energy balance components over the urban regions using the ground observations and numerical modeling (Pearlmutter et al. 1999(Pearlmutter et al. , 2006(Pearlmutter et al. , 2005Crawford and Bluestein 2000;Moriwaki and Kanda 2004;Offerle et al. 2005;Hanna et al. 2011;Kotthaus and Grimmond 2014;Templeton et al. 2018).
The geometry, structure, and spread of an urban area have an important role in quantifying the changes in the fluxes over the surface and the storage. The high-resolution remote sensing techniques are highly effective in recent years for monitoring the changes in the land use land cover (LULC) patterns and map the abrupt and rapid changes in the urban areas that influence the corresponding surface heat fluxes (Kato et al. 2008). Kato and Yamaguchi (2005) separated the artificial increase in the sensible heat fluxes using the surface heat balance assumption over an urban area of Nagoya, Japan, using ASTER and Landsat data. In 2007, Kato and Yamaguchi proposed to consider the storage heat as the heat flux between the surface and the canopy layer, and the results reasonably agreed to earlier studies Spronken-Smith 2002;Grimmond et al. 2004; Moriwaki and Kanda 2004;Christen and Vogt 2004;Offerleet al. 2005). Roupioz et al. (2016) attempted to generate a time series of radiation fluxes of better spatial and temporal resolution using remote sensing data over the heterogeneous topography of the Tibetan Plateau. Yang et al. (2016) estimated the land surface temperature (LST) from remotely sensed data with and without considering the urban geometry to evaluate the turbulent heat flux. Consideration of urban geometry appears to lower the values of sensible heat. Brenner et al. (2017) used a high-resolution thermal imager mounted on an unmanned aircraft to estimate LST and turbulent energy fluxes. In 2018, Chrysoulakis et al. investigated the potential of the satellite observations data (MODIS and Sentinel 3) downscaled to 100 m × 100 m to estimate UEB fluxes at a local scale along with the suitable ground observations. In 2019 Yang et al. employed a numerical microclimate model-TUF-3D (Temperatures of Urban Facets in 3-D)to obtain radiometric and complete surface temperatures for better estimation of SHF. Hrisko et al. (2021) utilized multispectral satellite radiance for quantifying the storage heat through the gradient-boosted regression trees method. Wetherley et al. (2021) investigated the variability of urban energy fluxes for land cover and climate gradient over Los Angeles, USA, using remote sensing and energy balance model. Many remote sensing modeling methodologies are adopted in recent studies to estimate the influence of urbanization over the surface heat fluxes (Kuang et al. 2015;Eswar et al. 2017) and address the energy-balance-closure problem (Nelli et al. 2019;Mauder et al. 2020). Earth observations and ground measurements are utilized in estimating the surface energy balance over urban areas. The ground observation method requires a sufficient number of weather stations, a vast amount of data, lots of workforces, and huge cost within a small observational area for a productive study (Hanna et al. 2011;Kotthaus and Grimmond 2014;Templeton et al. 2018), whereas the earth observation data is easily accessible and covers a larger area. Recent studies indicated the potential of earth observation in accounting for the changes in the fluxes is high (Weng et al. 2014;Chen and Hu 2017;Eswar et al. 2017;Rahman and Zhang 2019;Liang et al. 2019).
The phenomena of UHI is studied widely using earth observations and remote sensing tools, but the corresponding heat fluxes studies are limited to selective studies (Chrysoulakis et al. 2018). According to the review of literature, the study corresponding to surface energy fluxes is emphasized in western tropical/subtropical cities. India is the seventh largest country by area and the second most populated country in the world and has many highly populated metropolitan cities with rapidly spreading infrastructure. The smaller cities are developing widely for the recent couple of decades, leading to modifications in the regional climate and hence summons a thorough investigation of energy fluxes over Indian metropolitan cities. A few recent studies of the surface energy fluxes related to the Asian summer monsoon over India (Mohanty et al. 2003;Mohanty et al. 2005;Bhatla et al. 2011Bhatla et al. , 2016 are the significant initial landmark in this direction. In our earlier studies, the pattern of land surface temperature (LST) and UHI intensities over the different LULC class changes are analyzed over 11 (7) Indian metropolitan cities for winter (summer) seasons during 2000-2018 (Sultana and Satyanarayana 2018. Considering the fast development of the cities in India, it is essential to quantify the surface energy balance and investigate its influence on the temperature rise associated with UHI over these cities. The present study attempts to quantify the various surface energy balance components over the different Indian metropolitan cities and investigate the influence of rapid urbanization on variation in the UEB, specifically on sensible and latent heat fluxes, using Landsat 7 and Landsat 8 data.

Data
Landsat 7 and Landsat 8 imageries (cloud-free conditions) of the representative day of the winter seasons (December-February) in 3 to 6 years' intervals during 2000-2018 are utilized for the present study over each study region. The Landsat images are obtained from the USGS Earth Explorer user interface (https:// earth explo rer. usgs. gov/). The details of the data used for the present study are given in Table 1. For the present study, the imageries are pre-processed through spatial, spectral, and radiometric corrections using ERDAS IMAGINE 2014.
The LULC classification of the study areas for the study periods is retrieved using the multispectral bands of the Landsat (Sultana and Satyanarayana 2018. The reflectance and the spectral radiance are also estimated for the multispectral bands (Bands 1-5 and 7 of Landsat 7 and Bands 2-7 of Landsat 8) are used for estimating the albedo (Allen et al. 2002) and the Normalized Difference Vegetation Index (NDVI) (Carlson and Ripley 1997;Sultana and Satyanarayana 2018). The at-satellite brightness temperature can be calculated using the infrared (IR) thermal bands (Band 6 of Landsat 7 and Bands 10 and 11 of Landsat 8). The brightness temperature along with the NDVI is employed in estimating the emissivity-corrected LST over the study regions. The albedo, LST, and the LULC classification with the meteorological data are utilized in evaluating the components of the surface energy balance, discussed in Sect. 3. The relative digital elevation model (DEM) data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is used for estimating the obstacle height (https:// mirad or. gsfc. nasa. gov/).
The meteorological data parameters such as the nearsurface air temperature (T near-surface air ) data are collected from NCEP/DOE 2 Reanalysis data products, Physical Science Division (PSD), NOAA Earth System Research Laboratory. The available data are 4 times daily averaged with a spatial coverage of 2.5 × 2.5 degree grids, where the temperature is collected at 2 m height (https:// psl. noaa. gov/ data/ gridd ed/ data. ncep. reana lysis2. html/). Precipitable water (w) data is obtained from the University of Wyoming (http:// weath er. uwyo. edu/ upper air/ sound ing. html/), the upper air soundings data. Other meteorological data (wind speed, atmospheric temperature, relative humidity, dew point temperature) are obtained from local weather stations present approximately at the center of our chosen study region from Weather Underground atmospheric soundings (https:// www. wunde rgrou nd. com/).

Methodology
Investigating the influence of rapid urbanization on the surface energy balance requires different methodologies to estimate the energy components. The expression for the surface energy balance is (Oke 1982(Oke , 1988Masson et al. 2002): where Q* is the net all-wave radiation flux (Wm -2 ), Q H and Q E are the turbulent sensible and latent heat flux (Wm -2 ), and Q F is the anthropogenic heat flux (Wm −2 ). The ΔQ S is the total change in heat storage in the surface (Wm −2 ), and ΔQ A is the total change in advected heat flux (Wm −2 ), which is considered to be small in comparison to other components and could be neglected.
The net radiation flux can be calculated using the surface energy balance algorithms for land-SEBAL algorithm (Allen et al. 2002). The turbulent heat fluxes are estimated by the use of the aerodynamic resistance method-ARM (Nishida et al. 2003;Yamaguchi 2005, 2007;Kato et al. 2008). The heat storage and anthropogenic heat flux over urban surfaces cannot be estimated separately by remote sensing data analysis as they depend on surface characteristics and anthropogenic activities. Hence, the Q F and ΔQ S are combined as the residual heat flux, Q R , and estimated from the energy balance equation as follows Yamaguchi 2005, 2007).
The net radiation also called net flux at the surface is the actual radiant energy available to earth at its surface. It can be calculated by subtracting the total outgoing fluxes from the total incoming fluxes (Allen et al. 2002). It can be expressed as: where R ISW is the incoming short wave (SW) radiation reaching the earth surface (Wm −2 ), R ILW is the incoming longwave (LW) radiation emitted by the atmosphere (Wm −2 ), R OSW is the SW radiation reflected the atmosphere (Wm −2 ), and R OLW is the outgoing LW radiation that emitted by earth surface into the atmosphere (Wm −2 ). The Eq. (3) can be rewritten as: where α is the albedo of the earth surface and ϵ 0 is the planetary surface emissivity. The additional term in Eq. (4), (1 ε 0 ) R ILW , represents the fraction of incoming LW radiation that is reflected (Allen et al. 2002). (1) R ISW is estimated using the relation given in the following (Allen et al. 2002): where G sc is the solar constant (1367 Wm −2 ), τ is the atmospheric transmissivity, D r is the Earth-Sun distance corrected for eccentricity of planet orbit, and θ is the Zenith angle. The D r and θ can be obtained using Metadata available with the Landsat 8 imageries. The atmospheric transmissivity can be estimated by the equation as follows (Allen et al. 2007): where P is the local atmospheric pressure (K Pa), w is the precipitable water (mm), and K t is the unit less air turbidity coefficient (0 < K t ≤ 1; for clear air, 0.5 for polluted air).
Surface albedo can be defined as the ratio of reflected solar radiance to the amount of irradiance received. The albedo is estimated from the following expression used by Allen et al. (2002): where α is the planetary albedo corrected for atmosphere effects, α TOA is the planetary albedo without correction for corresponding atmosphere effects, and α atm is the atmospheric albedo (a typical constant value of 0.03 is taken based on Bastiaanssen 2000).
The α TOA can be estimated using the method given by Allen et al. (2002) and Silva et al. (2016): where ρ λ is the reflectivity of the multispectral bands and ω λ is the weighting coefficient for each band, which are computed using methods mentioned in the Landsat User Handbooks (Landsat 7 Data Users Handbook 2019; Landsat 8 Data user Handbook 2019).
Both the incoming and the outgoing LW radiation are estimated using the Stefan-Boltzmann equation given as follows (Allen et al. 2002;An et al. 2017): where σ is Stefan-Boltzmann constant (5.67 × 10 −8 in W m −2 K 4 ) and ϵ air is the atmospheric emissivity computed using the relation which is a function of atmospheric transmissivity (Bastiaanssen 1995;Allen et al. 2002), The planetary surface emissivity ϵ 0 is calculated using the NDVI (Carlson and Ripley 1997), and LST is obtained using the thermal bands of Landsat data. The methodology to calculate LST and ϵ 0 from Landsat imageries is discussed in detail by Sultana and Satyanarayana (2018. Sensible heat flux is the conductive transfer of heat flux from the Earth's surface to the atmosphere. Sensible heat flux can be estimated using the equation Yamaguchi 2005, 2007;Kato et al. 2008): where ρ is the air density (1.225 kg m −3 ), C P is the specific heat of air at constant pressure (1005 kJ Kg −1 K −1 ), T atm is the Atmospheric air temperature (K), and R a is the aerodynamic resistance (sm −1 ) for heat and momentum exchange and can be estimated by using the equation as given by Brutsaert (1982): where Z ref is the reference height (m) at which wind speed and temperature are measured, U is the wind speed (m s −1 ) at the reference height, and D 0 is the displacement height (m). Z oh and Z om are the roughness length (m) for heat and momentum exchange, respectively, ψ m and ψ h are the stability correction functions for heat and momentum exchange respectively depending on Monin-Obukhov length, and K (= 0.4) is the von Karman constant. For the present study, typical values of Z oh and Z om were considered for different LULC classes (Table 2) as proposed by Yamaguchi (2005, 2007); Kato et al. (2008) as it is quite difficult to acquire specific required data to calculate the Z oh and Z om values by the methods proposed earlier (Macdonald et al. 1998;Grimmond and Oke 1999;Duijm 1999). The displacement height D 0 is estimated by the expression given by Macdonald et al. (1998) and Kato and Yamaguchi (2005): where z h is the obstacle height (m), A is a constant (= 4.43 for staggered arrays), and λP is the plan area density of obstacles, i.e., the ratio of the built-up area to the total area. The obstacle height is being estimated from ASTER DEM data. The stability correction function is calculated as a function of Monin-Obukhov length by the proposed methodology (Paulson 1970;Webb 1970;Allen et al. 2002;Prueger and Kustas 2005). Latent heat flux is the flux of energy associated with the evaporation or transpiration of water from the Earth's surface to the atmosphere and vice versa. Latent heat flux is estimated from the following equation Yamaguchi 2005, 2007;Kato et al. 2008): where γ is the psychrometric constant (0.67 hPa K −1 ), e s is the saturated water vapor pressure (hPa) at surface temperature, e a is the atmospheric water vapor pressure (hPa), and R S is the stomatal resistance. Stomatal resistance can be computed using the proposed method as follows (Nishida et al. 2003;Kato et al. 2008): where R smin is the minimum stomatal resistance and R cuticle is the canopy resistance related to the diffusion through the cuticle layer of leaves (10 5 s m −1 ). The typical values of R smin over different LULC classes are considered as proposed by Kato et al. (2008) and Yamaguchi (2005, 2007), which are summarized in Table 3. The f 1 (PAR) and f 2 (PAR) are given by (Nishida et al. 2003): where T n , T x , and T 0 are the minimum, maximum, and optimal temperatures for stomatal activity, and the typical values are considered as 2.7 °C, 45.3 °C, and 31.1 °C, respectively. A is the parameter concerning photon absorption efficiency at low light intensity (152 µmol m −2 s −1 ), and PAR is the photosynthetic active radiation (W m −2 ), which can be calculated as f * R ISW (where f is 2.05 µmol m −2 s −1 ).

Results and discussion
In this section, the energy components estimated over the study region are analyzed with respect to the changes in the LULC classes, specifically over BA and DL possessing the higher LST. Additionally, the mean of the energy balance components during the study period are analyzed over the DL, BA, and the vegetated areas. The proportional sharing of the net radiation flux in heat fluxes are estimated for each case and are compared with the earlier studies. In our earlier studies Satyanarayana 2018, 2020), the cities under consideration are classified into 5 LULC classes, water bodies (WB), built-up areas (BA), dry lands (DL), crop/grass land (CGL), and dense vegetation (DV), and the methodology for the classification, variation in the LULC classes, and the accuracy assessment of these classes are also discussed in detail. The LST pattern and the UHI intensities with the variation in LULC classes over the cities are also analyzed.

Estimated surface energy components
During the study, the net radiation flux and the latent heat flux (LHF) over the BA and the DL are comparatively lower than the vegetated areas and the WB. The WB and wetlands possess the highest values for the LHF and net radiation flux. On the other hand, the sensible heat flux (SHF) has higher values over the BA and DL than the vegetation and the water bodies. The WB has the lowest range of SHF values (low positive or negative). As the present study focuses on the high LST possessing BA and DL, in this section, the analysis is conducted over these LULC classes only, and hence, Figs. 2-12 depict the energy balance components only over the DL and the BA, and the portions over the vegetated areas and the water bodies are removed.

Kolkata
Net radiation flux The spatial distributions of the net radiation flux during the study period for BA and the DL are shown in Fig. 2a Residual heat flux Fig. 2d shows the spatial distribution of the residual heat flux over BA and the DL, and the ranges over these regions are noticed to be lower in comparison to the other LULC classes (vegetation and water bodies). The values are noticed to be higher in the outer region of the BA during the latter part of the study period due to the presence of the DL.

Visakhapatnam
Net radiation flux Fig. 3a  Residual heat flux Fig. 3d depicts the spatial distribution of the residual heat flux over BA and the DL.
The ranges of residual heat over the DL are observed to be lower than that of BA, and the ranges over water and vegetation are higher than both the BA and DL ranges.

Latent heat flux
The spatial distribution of LHF over BA and DL during the study period is depicted in Fig. 4c. The LHF over the city appears to be mostly uniform except for some exceptions. The DL and the pavements are noticed to have the lowest range of values less than 60Wm −2 during the study period, whereas BA has little higher values of ranges Residual heat flux Fig. 4d shows the spatial distribution of the residual heat flux over BA and the DL, and the ranges over these regions are noticed to be lower, where the DL has the lowest range of values followed by moderate values over the BA.

Kochi
Net radiation flux The spatial distribution of net radiation flux over BA and DL for Kochi is shown in Fig. 5a Residual heat flux Fig. 5d shows the spatial distribution of the residual heat flux over BA and the DL for Kochi during the study period. The BA is noticed to possess the lowest values of residual flux followed by the DL with moderate values.

Delhi
Net radiation flux The spatial distribution of net radiation flux for Delhi over BA and DL is shown in Fig. 6a Residual heat flux The DL and BA are noticed to store lesser residual heat in comparison to vegetated areas and water bodies during the study period. The residual heat flux over the DL is noticed to be lower than the BA in the case of Delhi as shown in Fig. 6d.

Chandigarh
Net radiation flux The spatial distribution of net radiation flux for Chandigarh over BA and DL is depicted in Fig. 7a. The figure indicates that the lowest range of values of net radiation flux is mostly over the DL, pavements, and fallow lands (blue shades). The ranges over these areas are 229.   Residual heat flux The spatial distribution of the residual heat flux over the BA and the DL are depicted in Fig. 7d. The values of residual heat flux over the BA and the DL are noticed to be lower, and the DL and the industrial areas have the lowest range of values.

Jaipur
Net radiation flux Fig. 8a shows the spatial distribution of net radiation flux over BA and DL of Jaipur. The dry barren lands, sandy-rocky region, and fallow lands in the outer region are noticed to possess a lower range of values.  Latent heat flux The spatial distribution of LHF over BA and DL during the study period is depicted in Fig. 8c Residual heat flux Fig. 8d shows the spatial distribution of the residual heat flux over BA and the DL for Jaipur and the ranges over these regions are lower than the vegetation and water bodies. The DL and the fallow lands possess the lowest range of values followed by moderate to higher values over the BA.

Lucknow
Net radiation flux The spatial distribution of net radiation flux for Lucknow over BA and DL is depicted in Fig. 9a.

Nagpur
Net radiation flux The spatial distribution of net radiation flux over the BA and the DL of Nagpur are shown in Fig. 10a

Residual heat flux
The spatial distribution of the residual heat flux over BA and the DL of Nagpur are depicted in Fig. 10d. The values of residual heat flux over the BA and the DL are lower in comparison to that of the vegetation and the water bodies with DL having the lowest range of values.

Hyderabad
Net radiation flux The spatial distribution of net radiation flux over BA and DL for Hyderabad is shown in Fig. 11a Fig. 11d shows the spatial distribution of the residual heat flux over BA and the DL for Hyderabad during the study period. The DL is noticed to possess the lowest values of residual flux followed by the BA with moderate values which are much lower than that of the vegetated and the water bodies.

Bengaluru
Net radiation flux The spatial distribution of the net radiation flux over the BA and DL are shown in Fig. 12a

Variation in surface energy component over different LULC classes
During the study period, the artificial surfaces or the BA over the study regions are noticed to increase significantly Satyanarayana 2018, 2020). The surface heat fluxes are therefore noticed to vary considerably during this period. In this section, the mean net radiation flux (Q*), the mean SHF (Q H ), the mean LHF (Q E ), and the mean residual heat flux (Q R ) with respect to different LULC classes such as BA, DL, and vegetation over the cities under consideration are discussed. The variations in the mean values of the energy components, ∆Q*, ∆Q H , and ∆Q E , during the study period for these cities are tabulated in Table 3. The ∆Q*and ∆Q H are noticed as positive, whereas the ∆Q E is negative for all the study regions.

Visakhapatnam
The mean of the energy components with respect to different LULC classes (BA, DL, and vegetation) over Visakhapatnam is shown in Fig. 13b

Mumbai
From Fig. 13c

Jaipur
The mean of the energy components with respect to different LULC classes (BA, DL, and vegetation) is depicted in Fig. 14c for Jaipur for the study period. The Q* over BA

Hyderabad
For Hyderabad, Fig. 15b  The above analysis indicates that the mean SHF Q H is increasing significantly over the BA and the DL during the study period with a maximum increase of 28.2 Wm −2 (39.7 Wm −2 ) over the BA (DL) of Delhi. The BA of Kolkata (24.4 Wm −2 ) and Mumbai (23.0 Wm −2 ) also have higher values of ∆Q H . On the other hand, the mean LHF Q E is noticed to slightly decrease (< 20 Wm −2 ) over all the LULC classes during the study period. The Q* is also increased considerably during the study period.
In the present study, the cities are selected for different parts of India, such as coastal cities in east and west, northern inland cities, and southern peninsular cities, and the UEB components over these cities are expected to be influenced by their geographic conditions. In this section, we tried to analyze the UEB components based on the different regions of India. The UEB components over the coastal cities are given in Table 4, and over the inland cities are given in Table 5. The results indicate that, over northern cities, the net radiation flux and the residual heat flux are lower than that of the rest of the cities. The tables also show that the SHF over the DL of the northern and the southern inland cities are higher (125-155 Wm −2 ) than the coastal cities (90-120 Wm −2 ), whereas over the BA, the west coast cities have high SHF values (125-150 Wm −2 ) than the rest of the regions. The LHF value is high (around 30-45 Wm −2 ) over both the BA and the DL for the west coast cities and BA of the northern cities. For the east coast and the southern cities, the LHF is around 20-30 Wm −2 . Over the vegetation, the LHF is more than 100 Wm −2 for the west coast cities and about 70-90 Wm −2 for the rest of the cities. The SHF over the vegetations remains low (40-65 Wm −2 ) for all the cities under consideration with slightly higher values for the inland cities.

SHF and LHF values with respect to LST over different LULC classes
The above results indicate that the LHF values are low due to the absence of moisture over the BA and the DL of the study regions and the SHF is high over these areas. The vegetated areas possess high LHF values and low SHF values. In our earlier studies, the estimated LST values over the BA and the DL are higher than the vegetations, and the hotspots are identified mostly over the DL and dense BA for the study regions Satyanarayana 2018, 2020). In this section, the SHF and the LHF values are analyzed with respect to the LST ranges over different LULC classes (Tables 6,7,8,9,10,11). For Kolkata (Table 6), Kochi  (Table 7), and Hyderabad (Table 11), the SHF (LHF) values are higher (lower) over the BA than the DL, but the LST is noticed in the earlier studies to be higher over the DL than that of BA Satyanarayana 2018, 2020), whereas for the rest of the cites, SHF (LHF) have higher (lower) ranges over the DL than that of the BA like the LST ranges. On the other hand, the vegetated areas possess higher (lower) LHF (SHF), and thus, LST ranges are noticed to be considerably lower than the BA and DL. The ranges of the SHF and the LHF over the BA and DL are overlapping for all the cities like the LST values are noticed to overlap over these regions.

Estimation of proportion of the heat fluxes to the net radiation flux
In this section, the ratios of the mean heat fluxes to the mean net radiation flux are analyzed and compared with earlier studies. The ratio of the mean SHF, Q H , mean LHF, Q E , and mean residual heat flux, Q R , to the mean net radiation flux, Q*, is estimated for the urban regions of the cities under  of the net radiation flux and is decreasing for all the cities during the study period. The residual heat flux is about 60 to 80% of the net radiation flux. For Hyderabad, the residual heat is higher among the cities under consideration, due to the dense BA and the rocky terrain of granites over the suburban area. Delhi and Kochi also possess a higher proportion of residual heat, which could be attributed to the dense BA and the soil types of the cities. The proportion of the heat flux to net radiation flux found in the present study is compared with the earlier studies for the winter season in Table 13. It indicates that the present study agrees well with Oke et al. (1999), Moriwaki and Kanda (2004), Yamaguchi (2005, 2007), Kato et al. (2008), and Hanna et al. (2011). The residual flux as estimated by Kato and Yamaguchi (2005) appeared to be over-estimated. The SHF and LHF over the commercial area by Kato and Yamaguchi (2007) appear to be underestimated, and the residual heat flux is overestimated. The comparison of present results with the published ground observation might appear unsuitable as they are mean daytime values over different cities; the comparison could be considered

Summary and conclusion
The present study is concentrated on the variation of the UEB components over the different LULC classes, especially over the BA and the DL.   lowest over DL. The residual heat flux has higher ranges over the WB and the vegetation due to the high heat capacity of water, whereas the BA has moderate values and DL possesses the lowest values. For all the cities under consideration, the net radiation flux and the SHF are increased for all the LULC classes during the study period. The SHF over the BA and DL increased significantly during the study period, but the LHF decreased over these regions. The increase and decrease in the SHF and LHF can be attributed to the increase in artificial materials and the loss of moisture in the surface. The inland cities have high SHF over DL than the coastal cities, and the western coastal cities have high SHF over BA than the rest of the cities. The LHF is high over all the LULC classes for the western coastal cities than that of the rest of the cities. For the northern cities, the net radiation flux and the residual heat flux are lower than in the rest of the cities.
The analysis of the relation between SHF/LHF with the LST indicates that the SHF (LHF) values are higher (lower) over the BA than the DL for Kolkata, Kochi, and Hyderabad, whereas the LST ranges are higher over the DL than that of BA. For the rest of the cities, the SHF and LHF ranges are linear with the LST ranges. In the present study, the estimated proportion of heat flux to net radiation flux indicates that the SHF is about 19-33%, with higher values over Mumbai, Kochi, and Delhi. The LHF is about 1.9-15% of the net radiation flux for the study regions. In addition to that, the results also indicate that about 60-80% of net radiation flux ends up as residual heat flux. The proportion of SHF is slightly increased (by 4%) for Delhi and Kochi, and the proportion of LHF decreased during the study period. The comparison between the present results with the earlier studies appears to be in good agreement.
The surface energy balance components estimated in the present study are based on earth observation and could not be validated with ground observations due to the unavailability of high-resolution data. In addition to that, the energy balance components estimated over the water bodies contained errors for few study areas. Hence, the satellite imageries and the ground observations of higher resolution  are essential for better estimation of the energy balance components and the surface dynamics over the fast-developing cities. Inferences from the results derived from the present study have shown that substantial variations of surface energy fluxes related to different LULC classes could precisely influence the energy exchange mechanism from the BA and results in changes in the regional climate on the city scale.  Kato et al. (2008) 12 July 2004 0.40-0.70 0.60-0.80 Kato and Yamaguchi (2007) 2 January 2004 (commercial) 0.08 0.00 0.91 2 January 2004 (residential) 0.30 0.00 0.70 Kato and Yamaguchi (2005) 8