Parameterizing Anthropogenic Heat Flux with an Energy-Consumption Inventory and Multi-Source Remote Sensing Data

Anthropogenic heat (AH) generated by human activities is an important factor affecting the urban climate. Thus, refined AH parameterization of a large area can provide data support for regional meteorological research. In this study, we developed a refined anthropogenic heat flux (RAHF) parameterization scheme to estimate the gridded anthropogenic heat flux (AHF). Firstly, the annual total AH emissions and annual mean AHF of Beijing municipality in the year 2015 were estimated using a top-down, energy-consumption inventory method, which was derived based on socioeconomic statistics and energy consumption data. The heat released from industry, transportation, buildings (including both commercial and residential buildings), and human metabolism were taken into account. Then, the county-scale AHF estimation model was constructed based on multi-source remote sensing data, such as Suomi national polar-orbiting partnership (Suomi-NPP) visible infrared imaging radiometer suite (VIIRS) nighttime light (NTL) data and moderate resolution imaging spectroradiometer (MODIS) data. This model was applied to estimate the annual mean AHF of the counties in the Beijing–Tianjin–Hebei region. Finally, the gridded AHF data with 500-m resolution was obtained using a RAHF parameterization scheme. The results indicate that the annual total AH emissions of Beijing municipality in the year 2015 was approximately 1.704× 1018 J. Of this, the buildings contribute about 34.5%, followed by transportation and industry with about 30.5% and 30.1%, respectively, and human metabolism with only about 4.9%. The annual mean AHF value of the Beijing–Tianjin–Hebei region is about 6.07 W·m−2, and the AHF in urban areas is about in the range of 20 W·m−2 and 130 W·m−2. The maximum AHF value is approximately 130.84 W·m−2, mostly in airports, railway stations, central business districts, and other densely-populated areas. The error analysis of the county-scale AHF results showed that the residual between the model estimation and energy consumption statistics is less than 1%. In addition, the spatial distribution of RAHF results is generally centered on urban area and gradually decreases towards suburbs. The spatial pattern of the RAHF results within urban areas corresponds well to the distribution of population density, building density, and the industrial district. The spatial heterogeneity of AHF within urban areas is well-reflected through the RAHF results. The RAHF results can be used in meteorological and environmental modeling for the Beijing–Tianjin–Hebei region. The results of this study also highlight the superiority of Suomi-NPP VIIRS NTL data for AHF estimation.


Introduction
The world is undergoing unprecedented urbanization.The process of urbanization is often accompanied with urban area expansion, urban population growth, and energy consumption increase.The increasing human activities release more and more anthropogenic heat (AH) into the atmosphere.AH is involved in the energy exchange between the surface and atmosphere, hence the surface energy balance is impacted by AH directly.Much research in recent years has shown that AH is one of the key contributing factors to forming urban heat islands and has a significant influence on the regional environment and climate [1][2][3][4].Therefore, refined AH parameterization is required for both regional climate change and urban heat island studies.
Anthropogenic heat flux (AHF) is a sophisticated measurement of the heat released into the environment by human activities, i.e., AH emissions generated per unit time and area [5].Currently, there are three main approaches to estimating AHF: the energy budget closure method, building energy modeling, and the energy-consumption inventory approach [6].AHF is considered as a parameter of the surface energy balance equation in the energy budget closure method [7] (p.435).Many researchers have used this method for AHF estimation at different spatial scales [8][9][10][11].The building energy modelling approach explicitly calculates the energy consumption from different types of buildings and evaluates heat rejection based on the energy consumption model.AH from buildings is estimated accurately based on this method, but it ignores the heat emissions from transportation, industry, and so on [12][13][14].The energy-consumption inventory approach is a classical method for regional AHF estimation.It is generally classifies into either bottom-up or top-down approaches [6].The bottom-up approach relies on detailed statistics, such as traffic volume, building height, and occupied areas.The city-scale or larger spatial-scale AHF is scaled-up based on the AHF of single buildings [3].However, these detailed statistics are often difficult to obtain.Furthermore, the accuracy of AHF estimation is often subject to the level of detail and reliability of the statistics.However, the top-down method is usually used to estimate AHF based on socioeconomic data and energy consumption data over a large area.Compared to the detailed statistical data required for the bottom-up method, the statistics at city-scale or larger spatial-scales are more operable.Therefore, the AHF estimation over a large area is more easily implemented based on the top-down method [10,[15][16][17].
The statistics used in the top-down method are based on administrative divisions.Administrative divisions are divided into different levels (e.g., provincial-level, prefecture-level, county-level, etc.).The lower level corresponds to the smaller spatial scale, and the corresponding statistics are more difficult to collect.Therefore, the spatial scale of AHF results depends on the scale of statistics in the top-down method.In order to obtain AHF estimations in different areas and different spatial scales, large amount of statistics need to be collected (e.g., socioeconomic data and energy consumption data) [18][19][20].An AHF estimation model needs to be established to realize large spatial-scale AHF parameterization.In addition, AHF estimation based on the top-down method is the mean AHF value within the administrative divisions.The real spatial distribution of AH cannot be reflected by the mean AHF value.Then, spatial proxy data can be used to assign the mean AHF to a finer spatial scale (e.g., a grid-scale with a spatial resolution from about 1 km to 9 km), such as population density (PD) [15], gross domestic product (GDP) [21], land use types [22], etc.In other words, the refined anthropogenic heat flux (RAHF) can be obtained by distributing AHF to every grid.RAHF is needed in regional climate change and urban climate studies.Therefore, finding appropriate spatial proxy data to obtain the RAHF is the focus of current scholars.
Nighttime light (NTL) data can detect low intensity light emitted by city residents.It has a unique superiority in the monitoring of human activities at night, and it is also one of the ideal spatial proxy data of socioeconomic statistics and energy consumption data [23][24][25][26].Currently, DMSP/OLS (defense meteorological satellite program's operational linescan system) NTL data is being used for AHF estimation.For example, Chen et al. [27] estimated the global grid-scale AHF based on DMSP/OLS NTL data using the statistical regression relation between AHF and NTL data.Yang et al. [28] further used DMSP/OLS NTL data to extract urban areas, and then the statistical regression relation between the urban areas and energy consumption was established to estimate long-term AHF.In more recent studies by Dong et al. [15], DMSP/OLS NTL data was used to correct population density data, and then global grid-scale AHF was estimated based on the corrected population density data.However, due to the saturation of DMSP/OLS NTL data, the accuracy of AHF estimation in the saturated region needs to be improved.Hence, Ma et al. [29] and Chen et al. [30] calculated the human settlement index (HSI) based on DMSP/OLS NTL data corrected by vegetation index to alleviate the saturated problem to some extent, but there is still some need for further refinement on the spatial scale.
The National Aeronautics and Space Administration (NASA) launched the Suomi National Polar-orbiting Partnership (Suomi-NPP) satellite at the end of 2011.A new generation of global NTL data, hereinafter referred to as Suomi-NPP VIIRS NTL data, is obtained by a visible infrared imaging radiometer suite (VIIRS) equipped with the Suomi-NPP satellite.Compared with DMSP/OLS NTL data, Suomi-NPP VIIRS NTL data has a higher resolution and a greater range of radiance values.The low intensity light information can be captured by Suomi-NPP VIIRS NTL data to reflect more precisely the intensity and spatial distribution of human nighttime activities [31,32].
This research aims to establish the RAHF parameterization scheme over a large area and verify the availability of Suomi-NPP VIIRS NTL data in AHF parameterization.First, the annual total AH emissions and annual mean AHF of each district and county in Beijing municipality were calculated using the energy-consumption inventory approach.Then, a county-scale AHF estimation model was built based on the parameters characterizing human activities.In addition, multi-source remote sensing data, such as Suomi-NPP VIIRS NTL data and moderate resolution imaging spectroradiometer (MODIS) data were used to generate the spatial proxy data for grid-scale AHF parameterization within the administrative division.Finally, the grid-scale RAHF with 500-m resolution was parameterized to provide basic data support for regional climate change and urban climate studies.

Estimation of County-Scale Annual Mean AHF with a Top-Down Energy-Consumption Inventory
According to the emission sources of heat, AH emissions are divided into four parts: the energy consumption from industry (labeled as Q I ), buildings (Q B , including both commercial and residential buildings), transportation (Q V ), and human metabolism (Q M ).Total AHF (Q F ) is equal to the sum of the four parts [33].
The four parts were calculated after a comprehensive analysis of the statistical indicators of the China Statistical Yearbook and the method of Grimmond et al. [16].Q I is distributed according to the proportion of the second and third industrial GDP of each district and county.Q B from commercial buildings is distributed according to the proportion of the third industry GDP of each district and county, and that from residential buildings is distributed according to the proportion of the population.Q V is calculated based on the total amount of civil automobiles from socioeconomic statistics with an assumption of a car driving 2.5 × 10 4 km per year, consuming 12.7 L of fuel for every 100 km, and discharging fuel heat of 45 kJ•g −1 .Q M is calculated based on the population from socioeconomic statistics.Previous studies have divided one day into an active state (7:00~23:00), at which the metabolic heat emission is 171 W per person, and sleeping state (23:00~7:00) with 70 W metabolic heat emission per person [17].See Appendix A for details.
Although there is a time lag between energy consumption and actual heat emissions released into the atmosphere, due to the lack of accurate energy consumption efficiency data, this study makes the same assumption as previous studies: the energy consumption eventually translates into anthropogenic sensible heat with no time lag or consideration for environmental loads, and particular locations that emit latent heat such as cooling towers are not taken into account [17,34,35].

County-Scale AHF Estimation Model from Multi-Source Remote Sensing Data
The energy-consumption inventory approach is time-consuming and labor-intensive for AHF estimation over a large area.Moreover, the statistics are difficult to obtain for a large area.Thus some researchers proposed an AHF estimation model based on HSI to achieve AHF estimation over a large area [30].HSI is calculated based on the highly negative correlation between vegetation index and urban impervious surface.It is often used for the extraction of built-up areas at the regional or global scale [36].The equation is as follows: where NTL nor is the normalized NTL data, and NDVI max is the maximum normalized difference vegetation index (NDVI).NDVI and NTL data are complementary in the reflection of human activities and the extraction of built-up areas.This supplementary relation can effectively reduce the saturation phenomenon of NTL data.High-value HSI pixels often correspond to a low vegetation cover index, high NTL values and high AHF values.
The relationship between AHF and HSI is established using statistical regression.The coefficient of determination (R square) is used for the verification of the relationship between AHF and HSI, i.e., Equation (3).It is the proportion of the variance in the dependent variable that is predictable from the independent variable, hence the equation is shown as follows: n is the total number of samples.AHF is the mean of the AHF values computed as described in Section 2.1.AHF * is the AHF modeled value (i.e., AHF computed using the relationship between AHF and HSI).A high R square value indicates the good reliability of the forecast and goodness of fit.

RAHF Parameterized Model
There is a significant correlation between NTL data and the quantity and quality of the urban factors, including social factors (e.g., population size, population density, urbanization, etc.), economic factors (e.g., GDP, energy consumption, urban expansion, etc.), and ecological factors (e.g., carbon emissions, surface coverings, urban disasters, etc.) [37,38].The classification of Suomi-NPP VIIRS NTL data can reflect the intensity of human activities at different levels.RAHF parameterization based on the classification of Suomi-NPP VIIRS NTL data is conducive to generating grid-scale AHF with about 500-m spatial resolution.
A K-means clustering algorithm was used to determine the category k of RAHF modeling.The goal of this clustering algorithm is to achieve the highest similarity in class and the largest difference between classes [39]: where k is the category number of the NTL data, NTL k i is the NTL radiance of the pixel i in the category k, NTL k is the mean NTL radiance in the category k, n is the total number of the pixels in the category k, and m is the total number of the categories.Category m is taken for RAHF parameterization when the sum of variances of all categories attains its smallest value (m equals 30 in this study).This is done by seeking to minimize each NTL data category's average deviation from the mean value of the NTL data category, while maximizing each NTL data category's deviation from the means of the other NTL data groups.Finally, to get the AHF value of each NTL data category, the HSI of each category is usually the independent variable in the relation model between AHF and HSI mentioned in Section 2.2.1.
where AHF k is the AHF of the category k.It is a function about the HSI of the category k.The HSI is calculated by the Equation ( 1).The relationship between AHF and HSI has some uncertainty in the statistical regression procedure.The relationship is built based on the AHF calculated by the top-down energy-consumption inventory approach, which varies according to study area.Such as the timetable, the car driving miles per year, industrial structure and so on, these factors are going to lead to the diversity of this relationship model building.Therefore, ƒ HSI k i is not a general function.
According to the related data of the study area in this paper, the detail of ƒ HSI k i is shown as the Equation ( 6) in Section 4.2.

Procedure of the RAHF Parameterization Scheme
The annual mean AHF estimation at the county-scale based on the energy-consumption inventory approach is shown in a flowchart (Figure 1).AH is divided into Q I , Q B , Q V , Q M .Using a top-down energy-consumption inventory approach, the annual total AH emissions and annual mean AHF are estimated based on the socioeconomic statistics and energy consumption data of each district and county of Beijing municipality in the year 2015.The county-scale AHF estimation modelling based on multi-source remote sensing images is presented in Figure 1.First, the HSI of the Beijing-Tianjin-Hebei region is calculated based on Suomi-NPP VIIRS NTL data and MODIS NDVI data.Then, the mean HSI in each district and county of Beijing municipality is calculated by zonal statistics.The correlation analysis between the mean HSI and AHF of Beijing municipality districts and counties is performed to obtain the statistical regression relation.The statistical regression relation between them is fitted to get the AHF estimation model.The mean HSI in the counties of the Beijing-Tianjin-Hebei region is substituted into the AHF estimation model to get the county-scale AHF.
Finally, the Suomi-NPP VIIRS NTL data is classified to get the AHF categories (i.e., k, the category number of NTL data).The grid-scale RAHF of the Beijing-Tianjin-Hebei region is parameterized, integrating the AHF estimation model and the HSI at grid scale.The HSI is calculated using Equation (1).

Study Area and Data
The Beijing-Tianjin-Hebei region, the economic development center in northern China, is located at 113 • 27 E-119 • 50 E, 36 • 05 N-42 • 40 N, and includes Beijing municipality, Tianjin municipality and 11 prefecture-level cities of Hebei province (including Shijiazhuang, Cangzhou, Tangshan, Xingtai, Qinghuangdao, Baoding, Zhangjiakou, Handan, Langfang, Hengshui, and Chengde city).In this region, the city size is large and the urban heat island effect is evident.In 2014, the resident population was about 110 million people, and the regional GDP was about 6.65 trillion yuan.The total energy consumption was about 353.4 million tons of standard coal.The AH is mainly from the waste heat of automobile exhaust emissions, the energy consumption of industrial production and the various energy consumptions of buildings (e.g., winter heating and summer air-conditioning/refrigeration) [4].
The basic data for the energy-consumption inventory approach was derived from the statistical yearbook of social economy and energy consumption in Beijing municipality, Tianjin municipality and Hebei province in 2015, including China Energy Statistics Yearbook, China Population and Employment Statistics Yearbook, and the statistical yearbook of energy and environment in Beijing municipality, Tianjin municipality and Hebei province [40][41][42][43].
Multi-temporal Terra MODIS NDVI data (MOD13A1 product) and Suomi-NPP VIIRS NTL data were used in this research.All the selected data sets were acquired in 2015.Terra MODIS NDVI images (MOD13A1 product) with 500 m spatial resolution were downloaded from the United States Geographic Survey (USGS) [44].The product is a 16-day composite product [45].Multi-temporal NDVI mosaic images between April and October 2015 were collected to get high vegetation coverage data.MODIS NDVI mosaic images of h26v04, h26v05, h27v04, and h27v05 scenes were obtained to contain the full extent of the Beijing-Tianjin-Hebei region.The MODIS NDVI images with sinusoidal projection were reprojected to transverse mercator projection, and the quality control was performed based on the QC subset of the MOD13A1 product to obtain reliable NDVI images.
The monthly composites of Suomi-NPP VIIRS NTL data were downloaded from the National Oceanic and Atmospheric Administration/National Geophysical Data Center (NOAA/NGDC) [46].The monthly composite records the radiance value (in W•cm −2 •sr −1 ) and excludes any data impacted by stray light [47].The NTL data with geographic (Lat/Lon) projection were reprojected to transverse mercator projection, and the nearest neighbor resampling algorithm was used during the reprojection procedure.The method of Ma et al. [48] was used to remove the outliers of NTL data.The annual mean NTL data were calculated using a simple average value of the monthly mean NTL after the quality control.

Annual Total AH and Annual Mean AHF Results of Beijing Municipality
The annual total AH emissions and annual mean AHF results in the districts and counties of Beijing municipality are shown in Figure 2. The bar graph in Figure 2 shows the annual total AH emissions from different sources of Beijing in 2015, reflecting the characteristics of the annual total AH emissions and proportions of different emission sources.On the one hand, the annual total AH emissions in each district and county are quite different, and the annual total AH emissions of urban areas are significantly higher than those of suburbs.On the other hand, the emission source of AH is mainly from vehicular traffic, industry and residential buildings.The annual mean AHF in the districts and counties of Beijing municipality is in the range of 0.24 W•m −2 to 111.66 W•m −2 , as shown in the line graph in Figure 2.
According to the statistical analysis of the results, the annual total AH emissions of Beijing municipality in 2015 is about 1.704 × 10 18 J.Of this, the building sector contributes about 34.5%, followed by the transportation and industrial sectors with 30.5% and 30.1% respectively, and human metabolism with only 4.9%.The largest annual total AH emission is from the Chaoyang and Haidian districts, accounting for 39.5% of the annual total AH emissions of Beijing municipality.The smallest annual total AH emission is from Mentougou and the Yanqing counties, accounting for only 2% of the annual total AH emission of Beijing municipality.The results are in general agreement with the annual total AH emissions of Beijing municipality calculated by Wang et al. [19]. The

County-Scale Annual Mean AHF Estimation of the Beijing-Tianjin-Hebei Region
After the statistic regression analysis, it was found that there is a strong correlation between the annual mean AHF and mean HSI (R 2 = 0.989) (Figure 3).The AHF in Figure 3 The AHF estimation model was established based on the above statistical regression relation, and the AHF values of the counties in the Beijing-Tianjin-Hebei region were calculated using this AHF estimation model (Figure 4).The lowest AHF value in the 206 counties of the Beijing-Tianjin-Hebei region is 0.92 W•m −2 , and the highest AHF value is 125.93 W•m −2 .The average AHF is about 6.07 W•m −2 , and the AHF value of the high-value cluster area is between 80 W•m −2 and 130 W•m −2 .In order to quantitatively analyze the difference between the AHF value in urban areas and that in suburbs, a prefecture-level AHF was obtained from the county-scale AHF.Then the average AHF values of the districts (i.e., urban area) and counties (i.e., suburbs) in the cities of each municipality were calculated respectively (Table 1).The districts are the areas with a high degree of urbanization that are geographically located within the central city.Conversely, the counties are geographically far away from the city's center.The AHF value of an urban area is related to the city scale.The highest AHF value is in the urban area of Tianjin municipality with 39.21 W•m −2 , and the smallest is in Chengde city with 1.02 W•m −2 .According to statistics, the sum of the populations of Shijiazhuang and Tangshan cities accounts for 25.1% of the total population of Hebei province.The sum of these two cities' GDP accounts for 37.6% of that in Hebei province.At the same time, the sum of these two cities' AHF accounts for 44.15% of the total AHF in the 11 prefecture-level cities of Hebei province.This is indicative that AHF is closely related to population density, and to some extent to industrial and economic development levels.

RAHF Results and Validation of the Beijing-Tianjin-Hebei Region
The AHF estimation results shown in Figure 4 only represent the mean AHF within the administrative division and it cannot reflect the spatial heterogeneity of AHF.In order to reflect the spatial distribution characteristics of AHF within the administrative divisions, a grid-scale AHF with higher spatial resolution was obtained using the RAHF parameterization scheme.The results shown in Figure 5 conclusively demonstrate that the grid-scale RAHF is more reasonable than the annual mean AHF at the county-scale (Figure 4).The county-scale AHF results only reflect the average AHF value of the county area.The former more clearly reflects the AHF spatial distribution in the Beijing-Tianjin-Hebei region, and intuitively demonstrates the intensity of AH emissions.
It also can be found that the distribution of AHF in the Beijing-Tianjin-Hebei region is not uniform and has obvious aggregation characteristics.The distribution of AHF is consistent with that of human settlements.The magnitude of AHF is closely related to the level of economic development and population density.
The contribution of different AHF categories to the total AHF parameterization result is presented in Table 2.The AHF category of 0-10 W•m −2 has the highest contribution to the overall AHF of the Beijing-Tianjin-Hebei region, reaching 39.03%.It is mainly in the large non-urban area.Secondly, a high AHF value, which is more than 100 W•m −2 , is mainly in the city center with small area, and its contribution to the overall AHF is approximately 21.66%.As per the remaining AHF categories, it was shown that the contribution decreases gradually with the increase of the AHF value.
Considering the local characteristics of the RAHF results in the cities of Beijing, Tianjin and Shijiazhuang (the capital of Hebei province) as an example (Figure 6), it is noted that the spatially-distributed AHF in the urban area is well represented, especially for the major roadways around the urban area.As shown in Figure 6d-f, compared with the same range of high-resolution remote sensing images in Google Earth, most of the high AHF values are in airports, railway stations, industrial and central business districts.The typical regional energy consumption and population density is high and often accompanied with great heat emissions.This phenomenon further illustrates the validity of the RAHF results in the representation of local detail information.
Currently, the sensible heat and latent heat generated by human activities are difficult to clearly distinguish, hence it is hard to verify AHF results based on the situ measured flux data [5,6,50].The validation is often conducted using a comparative study with the results of other studies.
In addition, there are some AHF results at the global-scale or for the entire China.They can be used as validated reference data.However, there is no absolute comparability because of the inconsistency in research scale.Due to the lack of research on county-scale AH emissions for the Beijing-Tianjin-Hebei region, the credibility of the RAHF parameterization scheme is demonstrated by comparing with the AHF results for Beijing municipality.
Using the energy-consumption inventory approach, Tong et al.In this study, the residual between official energy consumption statistics and the AHF model results was calculated for the verification of the county-scale AHF estimation model [27].The county-scale statistics in the Beijing-Tianjin-Hebei region cannot be collected completely to verify the county-scale AHF estimation results; hence the AHF result calculated from the inventory method at the prefecture-scale was regarded as the actual AHF value to verify the model estimation results computed using Equation (6).Based on the data provided by the Bureau of Statistics, the real annual total AHF at the prefecture-scale were obtained using the energy-consumption inventory method.Afterwards the model estimation results were scaled up to get the prefecture-scale AHF for AHF verification on the same scale.Finally, on the prefecture-level, the residual between the AHF computed based on official statistics and the AHF result scaled up from the county-scale AHF estimation result computed using Equation ( 6) was calculated to verify the accuracy of the county-scale AHF estimation model.The verification results are presented in Table 3.The mean residual of all municipalities is 0.65%.Therefore, it is possible to consider using the correlation between AHF and reliable spatial proxy data to estimate AHF when there is lack of socioeconomic statistics and energy consumption data.

Comparison of Three County-Scale AHF Estimation Models
NTL data is confirmed to be one of the ideal spatial proxy data in the spatialization of socioeconomic data [23][24][25][26].Grid-scale population density is often used as spatial proxy data for AHF parameterization because of its close relation with AH [15,33,52,53].Thus, here we used normalized Suomi-NPP VIIRS NTL data (VIIRS nor ), population density (PD), and the human settlement index (HSI) as spatial proxy data to estimate county-scale AHF respectively.The feasibility of applying the three kinds of spatial proxy data in the county-scale AHF estimation was made by comparing the regression models between the three kinds of spatial proxy data and AHF (or the annual total AH emissions), respectively.
The mean and cumulative values of the three kinds of spatial proxy data were counted correspondingly as parameters for regression analysis.
First, the averages of the three kinds of spatial proxy data for each district and county in Beijing municipality were calculated to make a statistical regression with the annual mean AHF respectively.The R 2 values are 0.963, 0.995 and 0.989, respectively.It is proven that the averages of the three kinds of spatial proxy data have a significant correlation with annual mean AHF respectively.The specific fitting equation is shown in Table 4, and the fitted curves are the red lines shown in Figure 7a-c, correspondingly.The county-scale AHF's in the Beijing-Tianjin-Hebei region were estimated respectively by the three AHF estimation models mentioned above.Then, the prefecture-level AHFs were obtained from the county-scale AHFs (Figure 7d).It can be seen that the difference between the three results is small in the areas with a low AHF value and greater in the areas with a high AHF value (such as Tianjin municipality and Tangshan city).Compared with the AHF computed based on socioeconomic statistics and energy consumption data, which is regarded as the actual value, the estimation result based on HSI is closer to the actual value.The residuals of the three estimation models were calculated by using the accuracy verification method mentioned in Section 4.3.The residuals are 0.9% (AHF-VIIRS nor ), 0.7% (AHF-PD) and 0.6% (AHF-HSI), respectively.It is shown that the accuracy of the estimation model based on HSI is slightly higher.
Secondly, the cumulative values of the three kinds of spatial proxy data for each district and county in Beijing municipality were calculated to produce a statistical regression with the annual total AH emissions, respectively (Figure 7a-c).The results show that there is no significant correlation between them.As such, the cumulative value used for the model is not feasible in this study.
Comparison of the results for the three county-scale AHF estimation models suggest that VIIRS nor , grid-scale PD, and HSI can all be used to estimate county-scale AHF to some extent.According to the estimation of county-scale AHF in the Beijing-Tianjin-Hebei region as spatial proxy data, HSI can achieve slightly better performance than PD and VIIRS nor .

Accuracy Comparison between RAHF Results and Other AHF Products
This paper compares the spatial morphology of the grid-scale AHF results estimated by the RAHF parameterization scheme, the AHF results of Flanner et al. [53] and the AHF results estimated by the large scale urban consumption of the energy model (LUCY model) [52].The difference between the RAHF estimation and other AHF products can be revealed in this comparison.The LUCY model estimated the heat emissions from buildings, transportation and human metabolism using the energy-consumption inventory method.National energy consumption and traffic data, high-resolution population density data, and temperature data were used in this model [57].The AHF results with a 1-km resolution for the Beijing-Tianjin-Hebei region in 2015 were estimated based on various energy economic statistics and high-resolution population density data.The AHF results in typical areas are shown in Figure 8d-f Comparing the AHF results of these three AHF parameterization schemes, it can be seen that the AHF spatial variability of the RAHF parameterization scheme is the largest.Furthermore, the high values from the typical areas shown in Figure 6 are unapparent in the Flanner and LUCY model AHF results.The reason is that AHF is the anthropogenic heat released per unit area.For the AHF results, the value of AHF represents the mean value in one grid.Airports are often located in suburbs.The non-artificial land around airports had no anthropogenic heat emission, hence if the size of the grid is large, the AHF value of this grid will be attenuated by the non-artificial area.Therefore, the characteristics of AHF within an urban area can be clearly expressed by the RAHF result shown in Figure 6.In addition, land cover types and impervious surface percent (ISP) data were also used to further quantitatively analyze the accuracy of the RAHF results and other AHF products in this study.
First, the mean and maximum values of the three AHF results in different land cover types were calculated.The land cover data is the 300-m resolution global land cover classification system (LCCS) from the European Space Agency (ESA) [58].The land cover data was resampled to 500-m resolution and reprojected to transverse mercator projection.The statistical results are shown in Table 5.The mean AHF value in the settlement type is higher than that in the other land cover types.
The mean AHF in the settlement type of the RAHF parameterization scheme is 11.1 W•m −2 , and those of the LUCY model and Flanner are 1.24 W•m −2 and 3.00 W•m −2 , respectively.For the RAHF results, the gap of the mean AHF between the settlement type and other land cover types is the highest, up to 25 times, and those of the results of LUCY model and Flanner et al. [53] are the second and the lowest, respectively.Furthermore, Flanner's model results for agriculture and other land cover types are slightly higher.The mean AHF values of the three results in each ISP category were calculated based on ISP data.ISP is the proportion of total coverage by impervious surfaces in a total land area.The ISP in the Beijing-Tianjin-Hebei region was obtained by referring to the method of Xian et al. [59].The method is also adopted by the USGS to produce ISP products for America [60].
The ISP results have been divided into four categories: high-density impervious cover (70-100%), medium-density impervious cover (40-70%), low-density impervious cover (10-40%), and natural cover (0-10%).The statistical results are shown in Table 6.It can be noted that the highest mean AHF values for the three results are all for high-density impervious cover.It is revealed that the AHF value for larger ISP areas is greater, which also means more AH emissions.Xiao et al.'s [61] results showed that the mean AHF value in the high-density impervious cover area is 41.3 W•m −2 .The mean AHF value of RAHF results in the high-density impervious cover area is 33.57W•m −2 .The accuracy of RAHF results is significantly higher than that of the other two results.

Novelty and Opportunities of Suomi-NPP VIIRS Nighttime Light Data in RAHF Parameterization
The single-valued phenomenon within the urban area is caused by the saturation problem of DMSP/OLS NTL data (Figure 9a).However, differences of human activity intensity within urban areas can be performed by the new generation of NTL data, i.e., Suomi-NPP VIIRS NTL data (Figure 9b).In this study, the regression relation between the NTL data and AHF was indirectly constructed by the HSI, calculated based on NTL data and the NDVI.The negative correlation between the vegetation index and human activity reflected by the NTL data was used to effectively eliminate the saturation of the NTL data.Therefore, the residual between the estimated energy consumption and actual energy statistics is low when AHF estimation model is applied to an AHF estimation in a large area (Table 3).The spatial distribution of AHF within an urban area based on Suomi-NPP VIIRS NTL data is more in line with the distribution of the anthropogenic heat emission sources than that of the previous study [29,62], as shown in Figures 6 and 8.

Conclusions
In this paper, the 500-m resolution grid-scale AHF over the Beijing-Tianjin-Hebei region in the year 2015 was estimated based on a RAHF parameterization scheme.The RAHF parameterization scheme was established to provide AHF data for regional climate change and urban climate studies.The conclusions are as follows: 1.
For the county-scale AHF estimation model, the mean residual between the AHF estimation result and AHF computed based on the top-down energy-consumption inventory of all municipalities is less than 1%, indicating that the model can be used to estimate the county-scale AHF in the Beijing-Tianjin-Hebei region.

2.
According to the statistical regression analysis, Suomi-NPP VIIRS NTL data, grid-scale PD, and HSI can all be used to estimate county-scale AHF to some extent.For the estimation of county-scale AHF in the Beijing-Tianjin-Hebei region, the residuals were 0.9%, 0.7% and 0.6%, respectively.HSI can achieve slightly better performance than Suomi-NPP VIIRS NTL data and grid-scale PD.

3.
The spatial proxy data used to get grid-scale AHF within the districts and counties was established from multi-source remote sensing images.The 500-m resolution grid-scale RAHF was ultimately generated.According to a comparison of RAHF results and other AHF products, the RAHF parameterization scheme can get more refined AHF, and it has obvious advantages in the representation of spatial detail.Furthermore, the RAHF results of different underlying surfaces are also more reasonable.4.
After population density, GDP, land use data and DMSP/OLS NTL data, the validity of the Suomi-NPP VIIRS NTL data used in RAHF parameterization has been confirmed.It can produce more accurate and higher spatial resolution AHF results than DMSP/OLS NTL data.The results of this study indicate that Suomi-NPP VIIRS NTL data has a good application prospect for AHF estimation.
The RAHF parameterization scheme proposed in this paper used multi-source remote sensing data to simplify the labor-intensive county-scale inventory approach to some extent and has an ideal accuracy.However, the amount of training samples used to build the model was relative small, due to the deficiency of county-scale statistics.The model was preliminarily applied in the Beijing-Tianjin-Hebei region.The uncertainty still cannot be ignored in the inventory method.Further validation is needed for AHF parameterization in other areas.This is also a research direction for the future.
annual mean AHF of Beijing municipality in 2015 is 18.33 W•m −2 based on the area of each district and county, which is in line with the rule that the annual mean AHF of mid-latitude cities is approximately in the range of 15 W•m −2 to 50 W•m −2 [49].The largest annual mean AHF is in Xicheng district with 111.66 W•m −2 , and the smallest is in Yanqing county with 0.24 W•m −2 .The sum of the annual total AH emissions from the Dongcheng and Xicheng districts is smaller than that from the Chaoyang and Haidian districts.

Figure 2 .
Figure 2. The annual total AH emissions and the annual mean AHF in the districts and counties of Beijing municipality.
is the average AHF in each district and county of Beijing municipality, and the mean HSI was calculated based on the zonal statistics of the grid-scale HSI in each district and county of Beijing municipality.The statistical regression relation between them is AHF = 48.287(HSI) 2 − 17.716(HSI) + 2.541

Figure 3 .
Figure 3.The statistical regression relation between the annual mean anthropogenic heat flux (AHF) and the mean human settlement index (HSI).

Figure 4 .
Figure 4.The annual mean AHF of each county in the Beijing-Tianjin-Hebei region in the year 2015.

[ 4 ]
results showed that the AHF in urban areas of Beijing municipality was about 75-105 W•m −2 , and the AHF near Chongwenmen in Beijing municipality was between 130 W•m −2 and 170 W•m −2 .Nie et al. [51] estimated the AHF in Beijing municipality as having a peak value of 20 W•m −2 for low-intensity residential areas, 50 W•m −2 for high-intensity residential area and 90 W•m −2 for commercial areas.GDP and population density were used to obtain the street-scale AHF in Beijing municipality in Wang et al.'s [19] research.The results showed that the high AHF value was about 60-130 W•m −2 .AHF results based on a RAHF parameterization scheme is in the order of 80-130 W•m −2 in the urban area of Beijing municipality, and 130.84 W•m −2 in Chongwenmen in the same municipality.The results are in agreement with the historical research.

Figure 5 .
Figure 5.The RAHF results for the region in the year 2015.

Figure 6 .
Figure 6.RAHF results for typical cities in the Beijing-Tianjin-Hebei region and the validation of the high value area distribution ((a) AHF results for Beijing city; (b) AHF results for Tianjin city; (c) AHF results for Shijiazhuang city; (d-f) The corresponding high-resolution remote sensing images from Google Earth).Note: 1 Beijing-Capital International Airport; 2 Beijing Nanyuan Airport; 3 Beijing West Railway Station; 4 Beijing South Railway Station; 5 Beijing Railway Station; 6 Beijing East Railway Station; 7,19,20 Central Business District; 8 Tianjin Binhai International Airport; 9 Tianjin Railway Station; 10 Tanggu Railway Station; 11 Tianjin South Railway Station; 12-14,18 Industrial District; 15 Shijiazhuang Zhengding International Airport; 16 Shijiazhuang Railway Station; 17 Shijiazhuang North Railway Station.

Figure 7 .
Figure 7.The statistical regression model between annual mean AHF and the average and cumulative values of the three kinds of proxy data respectively, and the accuracy comparison of each model.((a) The correlation between anthropogenic heat flux (AHF) and anthropogenic heat emission (AHE) with Suomi-NPP VIIRS nighttime light data normalized value (VIIRS nor ); (b) The correlation between AHF and AHE with population density (PD); (c) The correlation between AHF and AHE with human settlement index (HSI); (d) Estimation and residual between official statistics and estimation at the prefecture-level in the Beijing-Tianjin-Hebei region).Note: The meanings of AHF-MVIIRS nor , AHF-MPD and AHF-MHSI are shown in Table 4. AHE-CVIIRS nor represents AHE and cumulative Suomi-NPP VIIRS nighttime light data normalized value (CVIIRS nor ) at the county level in Beijing municipality; AHE-CPD represents AHE and cumulative population density (CPD) at the county level in Beijing municipality; AHE-CHSI represents AHE and the cumulative human settlement index (CHSI) at the county level in Beijing municipality.
Flanner et al. [53] used the energy consumption data of the US Bureau of Energy Information Statistics [55] for the period of 2005-2040 and the population density data of the Columbia University Socioeconomic Data Application Center [56] to estimate the global grid-scale AHF with 3-km resolution from 2005 to 2040.The AHF results in typical areas in 2015 shown in Figure 8a-c correspond to the cities of Beijing, Tianjin and Shijiazhuang city and their surrounding areas, respectively.The range of AHF results in these areas is from 0 W•m −2 to 49.58 W•m −2 , and the high AHF value is about 38.70-49.58W•m −2 .
. The range of AHF results in these areas is from 0 W•m −2 to 79.89 W•m −2 , and the high AHF value is within 42.31-79.89W•m −2 .AHF results with a 500-m resolution estimated by the RAHF parameterization scheme in typical areas are shown in Figure6a-c.The range of the AHF results in these areas is from 0 W•m −2 to 130.84 W•m −2 , and the high AHF value is about 80-130.84W•m −2 .

Figure 8 .
Figure 8.The AHF results of Flanner et al. and the LUCY model.(a) AHF results for Beijing city; (b) AHF results for Tianjin city; (c) AHF results for Shijiazhuang city.

Figure 9 .
Figure 9.The profile of the nighttime light data in Beijing city.((a) DMSP/OLS (defense meteorological satellite program's operational linescan system) nighttime light data; (b) Suomi-NPP VIIRS (Suomi national polar-orbiting partnership visible infrared imaging radiometer suite) nighttime light radiance).

Table 1 .
The annual mean AHF of municipalities, districts and counties in the Beijing-Tianjin-Hebei region.

Table 2 .
The contribution of different AHF categories to the total AHF in the Beijing-Tianjin-Hebei region.

Table 3 .
Comparison between the AHF estimation results and the AHF computed based on official statistics.

Table 4 .
Comparison of three statistical regression models for county-scale AHF estimation.VIIRS nor represents the anthropogenic heat flux (AHF) and the mean Suomi-NPP VIIRS nighttime light data normalized value (VIIRS nor ) at the county level in Beijing municipality; PD represents the mean population density at the county level in Beijing municipality; his represents the mean human settlement index at the county level in Beijing municipality.

Table 5 .
The mean and maximum values of different AHF results for different land cover types.

Table 6 .
The mean AHF of different results in different impervious surface percent (ISP) categories.