Exploring Urban Heat Distribution and Thermal Comfort Exposure Using Spatiotemporal Weighted Regression (STWR)

: With rapid urbanization, many cities have experienced significant changes in land use and land cover (LULC), triggered urban heat islands (UHI), and increased the health risks of citizens’ exposure to UHI. Some studies have recognized residents’ inequitable exposure to UHI intensity. However, few have discussed the spatiotemporal heterogeneity in environmental justice and countermeasures for mitigating the inequalities. This study proposed a novel framework that integrates the population-weighted exposure model for assessing adjusted thermal comfort exposure (TCEa) and the spatiotemporal weighted regression (STWR) model for analyzing countermeasures. This framework can facilitate capturing the spatiotemporal heterogeneities in the response of TCEa to three specified land-surface and built-environment parameters (i.e., enhanced vegetation index (EVI), normalized difference built-up index (NDBI), and modified normalized difference water index (MNDWI)). Using this framework, we conducted an empirical study in the urban area of Fuzhou, China. Results showed that high TCEa was mainly concentrated in locations with dense populations and industrial regions. Although the TCEa’s responses to various land-surface and built-environment parameters differed at locations over time, the TCEa illustrated overall negative correlations with EVI and MNDWI while positive correlations with NDBI. Many exciting spatial details can be detected from the generated coefficient surfaces: (1) The influences of NDBI on TCEa may be magnified, especially in rapidly urbanizing areas. Still, they diminish to some extent, which may be related to the reduction in building construction activities caused by the COVID-19 epidemic and the gradual improvement of urbanization. (2) The influences of EVI on TCEa decline, which may be correlated with the population increase. (3) Compared with NDBI, the MNDWI had more continuous and stable significant cooling effects on TCEa. Several mitigation strategies based on the spatiotemporal heterogeneous relationships also emanated. The effectiveness of the presented framework was verified. It can help analysts effectively evaluate local thermal comfort exposure inequality and prompt timely mitigation efforts.


Introduction
With population growth and rapid socio-economic development, most cities' urban land use and land cover (LULC) have experienced significant changes.Many natural surfaces (e.g., forests, water) with low heat storage properties transformed into urban materials (e.g., buildings, bridges, and roads) with high heat storage properties [1], accelerating urban heat islands (UHI) [2].UHI generally refers to the higher air or surface temperature in urban areas than in surrounding rural areas.UHI can usually be categorized into atmospheric UHI (AUHI) and surface UHI (SUHI) [3].Our study focuses on the latter for the relatively easier accessibility of spatial and temporal data [4].Tan et al. [5] reported that people who were exposed to an environment of high UHI intensity for a long time are more likely to suffer from insomnia, listlessness, and other symptoms, leading to an increase in related diseases (e.g., cardiovascular, cerebrovascular, and respiratory system diseases) and deaths [6].
Exposure to high UHI intensity poses increasing health risks for citizens [7], especially given the growing consumption of global energy fossil fuels and rising emissions of greenhouse gases.With the development of the social economy, urban residents' awareness of high-quality life increased.Their expectations for less exposure to high UHI-intensity environments are growing daily.Many urban green spaces, ecological corridors, and other facilities have also been put into construction [8].However, the environmental justice concerns about urban heat islands remain [9].How to optimize urban LULC and arrange urban public green space more scientifically and rationally is still an unresolved issue closely related to the well-being of the citizens.
Scientifically measuring the population's thermal comfort exposure (TCE) is a prerequisite for formulating mitigation measures.TCE refers to the population's exposure to UHI intensity, with no population equating to no exposure risk.Therefore, TCE should be measured by considering the population distribution.However, the spatial distribution of the population and urban heat islands is inconsistent [10].On the one hand, areas with high UHI intensity may exhibit low TCE due to sparse populations.Estoque et al. [11] assessed populations' exposure and vulnerability to UHI intensity in 139 cities in the Philippines.They found that Luzon exhibited shallow TCE despite having higher UHI due to low population density.On the other hand, in some areas, the density of the population may increase TCE.Ref. [12] reported that although the UHI intensity is not high in areas with dense, tall buildings in Chongqing, China, the TCE is significant due to high population density.
In addition, some dense population areas also have a large amount of production activities or energy consumption.These will enhance UHI, increasing the populations' exposure to UHI intensity in these areas [13].Therefore, a method that links population and UHI is required to assess TCE.Wu et al. [14] employed a population-weighted exposure model to determine the changes in heat exposure across 398 major cities in the United States from 2000 to 2020.The model considered the spatial interactions between heat exposure and population distribution, enabling better quantification of the exposure to urban thermal environments in both temporal and spatial dimensions.Additionally, the model can facilitate long-term monitoring and assessment of heat island exposure at fine scales.However, the original TCE was easily affected by extreme local population values and may have fewer differences when the total population was significant.Therefore, this study adopts an adjusted thermal comfort exposure (TCEa) using a mean population of denominators.
The main factor affecting UHI is the land cover type (e.g., urban green vegetation, buildings, water bodies, etc.) [15][16][17].Scholars generally believe optimizing LULC can effectively mitigate exposure to urban heat islands [18,19].Chen et al. [20] analyzed the impact of underlying surface changes on UHI and found that increasing green space areas can effectively alleviate the scope and intensity of local heat islands.Chen et al. [21] investigated the driving factors of UHI in the urban agglomeration of the Yangtze River Delta in eastern China.The normalized difference built-up index (NDBI) [22] significantly impacts UHI.Sun and Chen [23] reported that water bodies could effectively reduce UHI through evaporation and heat absorption.
Exploring the relationship between the driving factors and TCEa is the key to mitigating TCEa.Wang et al. [24] utilized the ordinary least square (OLS) model to investigate the influencing factors of UHI in the Beijing-Tianjin-Hebei megaregion in China.They found that the enhanced vegetation index, population density, and gross domestic product are the main factors leading to UHI.However, TCEa is spatially heterogeneous [25].The global OLS model assumes that the entire study area is affected by the same impact, ignoring the spatial correlation and heterogeneity between variables [26].Liu et al. [27] used the geographically weighted regression (GWR) model to analyze the spatial relationship between UHI and potential driving forces in Tokyo, Japan.They found that a negative correlation between forest proportion and UHI is shrinking, while the proportion of urban fabric areas hurts UHI.Although the GWR model can detect spatial heterogeneity of TCEa, it does not consider the temporal dimension and, therefore, cannot capture temporal heterogeneity [28].The relationships between TCEa and land-surface and built-environment parameters are complex in space and time [29].Few studies discussed spatiotemporal heterogeneity.A novel spatiotemporal weighted regression (STWR) model [30], characterized by utilizing numerical difference rates rather than time intervals for weighting, could be a solution to fill the gap by quantifying local spatiotemporal non-stationarity.
Therefore, our study aims to develop a new framework to explore TCEa's spatiotemporal heterogeneity and corresponding mitigation strategies and verify their effectiveness through empirical research.

Study Area
Figure 1 shows our study area (25 , which is the urban area of Fuzhou, a city located in southeastern China.Fuzhou has a typical subtropical humid monsoon climate, with long summers and short winters, little frost, plenty of rain, and evergreens all year round.The seasonal variation of vegetation cover in Fuzhou is slight (https://www.fuzhou.gov.cn,accessed on 15 January 2024).The average summer temperature in 2022 is 29.5 • C, and the average number of days with temperatures above 35 • C ranges from 40 to 79.According to the "Fuzhou Statistical Yearbook" (https:// tjj.fuzhou.gov.cn,accessed on 15 January 2024), as of 2021, Fuzhou's urbanization rate reached 73%, with a population of about 8.42 million.The rapid urbanization process has significantly increased Fuzhou's economic GDP and modernization.However, it also changed the land cover pattern and intensified UHI [31].

Figure 1.
The study area is in Fuzhou City, China.

Data Sources and Preprocessing
The remote-sensing data used in this study (including Landsat8 Collection 2 Level-2, L8C2L2) was sourced from the USGS (https://earthexplorer.usgs.gov,accessed on 7 January 2024).The original 30 m spatial resolution data were cloud-free images for the study area from June to September 2017 to 2022.The observation times were 16 September 2017, 15 August 2017, 3 September 2018, 6 September 2019, 22 September 2019, 22 July 2020, 7 August 2020, 23 August 2020, 27 September 2021, and 12 July 2022.The 1 km spatial resolution population data from 2017 to 2022 was collected from Oak Ridge National Laboratory (https://landscan.ornl.gov,accessed on 9 January 2024).For years with multiple images, such as in 2017, the maximum value synthesis method was used to minimize the impact of clouds and other pollutants [32].
Various preprocessing methods, including radiometric calibration, atmospheric correction, geometric correction, cropping, and resampling, were performed to make these data consistent before fitting the models.The resampled 240 m spatial resolution referred to [33] to better investigate the relationships between TCEa and land-surface parameters.In addition, Table 1 lists the formulas of LST, land-surface, and built-environment parameters calculated from these data.

Name Formulas Descriptions
Spectral bands B3, B4, B5, B6, B10 Original bands reflectance of L8C2L2.B3 is the green band; B4 is the red band; B5 is the near-infrared band; B6 is the short-wave infrared band; B10

Data Sources and Preprocessing
The remote-sensing data used in this study (including Landsat8 Collection 2 Level-2, L8C2L2) was sourced from the USGS (https://earthexplorer.usgs.gov,accessed on 7 January 2024).The original 30 m spatial resolution data were cloud-free images for the study area from June to September 2017 to 2022.The observation times were 16 September 2017, 15 August 2017, 3 September 2018, 6 September 2019, 22 September 2019, 22 July 2020, 7 August 2020, 23 August 2020, 27 September 2021, and 12 July 2022.The 1 km spatial resolution population data from 2017 to 2022 was collected from Oak Ridge National Laboratory (https://landscan.ornl.gov,accessed on 9 January 2024).For years with multiple images, such as in 2017, the maximum value synthesis method was used to minimize the impact of clouds and other pollutants [32].
Various preprocessing methods, including radiometric calibration, atmospheric correction, geometric correction, cropping, and resampling, were performed to make these data consistent before fitting the models.The resampled 240 m spatial resolution referred to [33] to better investigate the relationships between TCEa and land-surface parameters.In addition, Table 1 lists the formulas of LST, land-surface, and built-environment parameters calculated from these data. 2.3.Methods

UHI Intensity
The relative brightness temperature (T R ) is used to characterize the UHI intensity [36] and can be defined as: where LST i denotes the land-surface temperature of location i; LST mean is the average surface temperature across the study area.

Population-Weighted Exposure Model
The inequality of TCEa can be measured using the population-weighted exposure model [14].We adjusted the denominator of the formula and replaced the total population with the regional population average, making the calculation results more suitable for spatial comparative analysis without affecting the evaluation results.Its basic formulation can be defined as: where TCE a denotes adjusted thermal comfort exposure.P i denotes the population of location i; P denotes the average population across the study area; T Ri is UHI intensity of location i.

STWR
STWR is a temporally extended model based on GWR, characterized using numerically different rate attenuation weighting strategies [37].The basic formula of the STWR model can be defined as [38]: where y t i and X t ik denote the dependent and independent variables of the i-th observation point, respectively.β t k (u i , v i ) is the regression coefficient of the i-th observation point, ε t i is a random error term.X s∆t is the matrix of independent variables observed within the time interval ∆t.W ∆t (u i , v i ) is the spatiotemporal weight matrix.Each weight in the matrix can be calculated by: where w t ijST denotes the observation j's weight at t on the regression point i. α is an adjustable parameter to weigh the strength of the local temporal impacts, optimized with the bandwidth.d Sij is the spatial distance.b ST and b T are the spatial and temporal bandwidth, respectively.y t i − y t−q j represents the numerical differences between the regression point i and observation j in period t − q, and ∆t = t − (t − q) is a time interval.

Research Framework
Figure 2 illustrates the analytical framework we designed, which consists of two parts: First, the assessment of TCEa based on remote-sensing images and population data using the population-weighted exposure model.Second, the analysis of mitigation countermeasures based on spatiotemporal heterogeneity between TCEa and significant driving factors using the STWR model.
Buildings 2024, 14, x FOR PEER REVIEW 6 of 18 where w ijST t denotes the observation j's weight at t on the regression point i. α is an adjustable parameter to weigh the strength of the local temporal impacts, optimized with the bandwidth.d Sij is the spatial distance.b ST and b T are the spatial and temporal bandwidth, respectively.y i t − y j t−q represents the numerical differences between the regression point i and observation j in period t − q , and ∆t = t − (t − q) is a time interval.

Research Framework
Figure 2 illustrates the analytical framework we designed, which consists of two parts: First, the assessment of TCEa based on remote-sensing images and population data using the population-weighted exposure model.Second, the analysis of mitigation countermeasures based on spatiotemporal heterogeneity between TCEa and significant driving factors using the STWR model.

Adjusted Thermal Comfort Exposure (TCEa)
Figure 3 shows the spatial distribution of population, UHI, and thermal comfort exposure (TCEa) in the study area.The TCEa illustrated an overall pattern of low in the west and high in the east.Most high TCEa locations were mainly concentrated in dense population areas, commercial activities areas like eastern Gulou and central Taijiang, and industrial areas like the east Cangshan district.Some interesting phenomena can be found: (I) Although the UHI intensities in eastern Gulou, central Taijiang, and Jinan were not the highest, the highest TCEa appeared due to the dense population in these areas, especially in 2017.(II) Despite the relatively high UHI intensities in western Cangshan, TCEa was not high due to the relatively low population density.(III) The TCEa range in eastern Cangshan was extensive and relatively high for a long time.In contrast, for western Cangshan, the TCEa was relatively low and increased, mainly contributing to population growth, especially since 2020.

Spatiotemporal Heterogeneity
Table 2 shows the multicollinearity diagnostic results.All the Variance Inflation Factors (VIF) [39] were less than 6, indicating no significant multicollinearity exists.Table 3 shows the comparison results of OLS, GWR, and STWR model performances.From 2018 to 2022, STWR outperformed the other models with a higher yearly R-squared (R 2 ) [40] and lower residual sum of squares (RSS) [41] and corrected Akaike's information criterion (AICc) [42], indicating best-fitting results.In 2017, the STWR and GWR results were the same for the first year because no previously observed data points can be used to improve the performance of the STWR model.The optimal temporal bandwidths of STWR were all 2, indicating that the best fit of the model is to use the data of the last 2 years since 2017.The optimal yearly STWR's spatial bandwidths were smaller than those of GWR since 1998, indicating that the STWR model was more localized than GWR.This may be related to the fact that it borrowed historical observation points and resulted in more precise details.

Impacts of EVI
Figure 4 shows the yearly box violin plot of significant EVI coefficients, indicating an overall negative and upward trend of influences on TCEa, which even turns weakly positive in 2022.This implied that the negative impacts of vegetation were weakened, and the surrounding environment or other factors may affect the cooling effect.positive in 2022.This implied that the negative impacts of vegetation were weakened, and the surrounding environment or other factors may affect the cooling effect.Figure 5 shows the spatial distributions of significant EVI coefficients.Since 2020, the TCEa has increased while the negative significant coefficient of EVI has decreased.In areas where TCEa was not particularly strong, vegetation can play a critical mitigating role.However, the cooling effect of EVI fluctuated with the changes in heat exposure.In some areas with intense TCEa, EVI turned to positive significant impacts.Figure 5 shows the spatial distributions of significant EVI coefficients.Since 2020, the TCEa has increased while the negative significant coefficient of EVI has decreased.In areas where TCEa was not particularly strong, vegetation can play a critical mitigating role.However, the cooling effect of EVI fluctuated with the changes in heat exposure.In some areas with intense TCEa, EVI turned to positive significant impacts.
positive in 2022.This implied that the negative impacts of vegetation were weakened, and the surrounding environment or other factors may affect the cooling effect.Figure 5 shows the spatial distributions of significant EVI coefficients.Since 2020, the TCEa has increased while the negative significant coefficient of EVI has decreased.In areas where TCEa was not particularly strong, vegetation can play a critical mitigating role.However, the cooling effect of EVI fluctuated with the changes in heat exposure.In some areas with intense TCEa, EVI turned to positive significant impacts.

Impacts of NDBI
Figure 6 shows the yearly box violin plot of significant NDBI coefficients, indicating an overall positive and downward trend of influences on TCEa.All the median significant coefficients were above zero, while the degree and magnitude of the impact were weaker than those of EVI.Before 2019, the dispersion of significant NDBI coefficients was small, while in 2019, the urbanization process accelerated, and infrastructure construction and new buildings increased, resulting in a rapid increase in dispersion and reaching a peak.The degree of dispersion declined in the following 2020 and 2021.

Impacts of NDBI
Figure 6 shows the yearly box violin plot of significant NDBI coefficients, indicating an overall positive and downward trend of influences on TCEa.All the median significant coefficients were above zero, while the degree and magnitude of the impact were weaker than those of EVI.Before 2019, the dispersion of significant NDBI coefficients was small, while in 2019, the urbanization process accelerated, and infrastructure construction and new buildings increased, resulting in a rapid increase in dispersion and reaching a peak.The degree of dispersion declined in the following 2020 and 2021.Figure 7 shows the spatial distributions of significant NDBI coefficients.The number of significant coefficients in the zoom-in area (Figure 7b,e,h,k) was the highest in 2019, while in the latter years, the number gradually decreased.Moreover, some significant negative effects appeared in some local locations.However, except for the few negative significant points, the most positive significant points were in the fast-developing areas of the generated coefficient surfaces of those years.Figure 7 shows the spatial distributions of significant NDBI coefficients.The number of significant coefficients in the zoom-in area (Figure 7b,e,h,k) was the highest in 2019, while in the latter years, the number gradually decreased.Moreover, some significant negative effects appeared in some local locations.However, except for the few negative significant points, the most positive significant points were in the fast-developing areas of the generated coefficient surfaces of those years.

Impacts of MNDWI
Figure 8 shows the yearly box violin plot of significant MNDWI coefficients.All the median significant coefficients were below zero except for 2022, indicating that water was primarily negatively correlated with TCEa.However, since 2019, the cooling effect of water has also appeared to decline, as has the degree of dispersion.

Impacts of MNDWI
Figure 8 shows the yearly box violin plot of significant MNDWI coefficients.All the median significant coefficients were below zero except for 2022, indicating that water was primarily negatively correlated with TCEa.However, since 2019, the cooling effect of water has also appeared to decline, as has the degree of dispersion.

Impacts of MNDWI
Figure 8 shows the yearly box violin plot of significant MNDWI coefficients.All the median significant coefficients were below zero except for 2022, indicating that water was primarily negatively correlated with TCEa.However, since 2019, the cooling effect of water has also appeared to decline, as has the degree of dispersion.Figure 9 shows the spatial distributions of significant MNDWI coefficients.Significant adverse or mitigating effects were shown near most water areas.In the highly developed downtown areas, the cooling effects of water were quite sound.For instance, in the areas surrounding West Lake Park areas (Figure 9c,f,i,l), the TCEa intensities were high, but inside the park, there was almost no TCEa.Figure 9 shows the spatial distributions of significant MNDWI coefficients.Significant adverse or mitigating effects were shown near most water areas.In the highly developed downtown areas, the cooling effects of water were quite sound.For instance, in the areas surrounding West Lake Park areas (Figure 9c,f,i,l), the TCEa intensities were high, but inside the park, there was almost no TCEa.

Discussion
This study proposed a new framework that can help measure TCEa and capture the spatiotemporal heterogeneity (the spatial heterogeneity of numerical different change rates of TCEa) correlations with land-surface and built-environment parameters in the progress of urbanization.As shown in Figure 3, although the UHI intensities of some locations in Gulou were low, their TCEa were high due to the dense populations.In contrast, the UHI intensities of some locations in western Cangshan were high, but their TCEa were relatively low due to the low population density.However, as the population of these locations grew, the corresponding TCEa also increased significantly in those years.In contrast to the GWR and geographically and temporally weighted regression (GTWR) [43] models, the STWR model integrated into the framework adopts a decay weighting strategy based on the numerical difference change rate of a certain time interval rather than the time interval itself, which enables exploring the heterogeneity among change rates at different locations [38].In terms of the spatiotemporal heterogeneous response of TCEa to

Discussion
This study proposed a new framework that can help measure TCEa and capture the spatiotemporal heterogeneity (the spatial heterogeneity of numerical different change rates of TCEa) correlations with land-surface and built-environment parameters in the progress of urbanization.As shown in Figure 3, although the UHI intensities of some locations in Gulou were low, their TCEa were high due to the dense populations.In contrast, the UHI intensities of some locations in western Cangshan were high, but their TCEa were relatively low due to the low population density.However, as the population of these locations grew, the corresponding TCEa also increased significantly in those years.In contrast to the GWR and geographically and temporally weighted regression (GTWR) [43] models, the STWR model integrated into the framework adopts a decay weighting strategy based on the numerical difference change rate of a certain time interval rather than the time interval itself, which enables exploring the heterogeneity among change rates at different locations [38].In terms of the spatiotemporal heterogeneous response of TCEa to different land-surface parameters, several exciting findings were detected through the coefficients with |Z-score| greater than 1.96 SD: (1) The negative impact of EVI on TCEa may diminish around green space as the population increases.In 2019, the number of coefficients of EVI (|z-score| > 1.96) was large in the zoom-in area but dropped significantly in subsequent years.(Figure 5b,e,h,k), although TCEa in the area increased (Figure 5c,f,i,l).This may be attributed to the growing population living in the area.On the one hand, the increased population promoted TCEa and accelerated socio-economic activities [44], increasing energy consumption.On the other hand, population growth led to increasing hardened ground surfaces, which reduced water infiltration and further limited the natural cooling process of plants through transpiration [45].(2) NDBI's impact on TCEa may be amplified, especially in rapidly urbanizing areas.Still, it may diminish to a certain extent due to the reduction in urban construction activities caused by the COVID-19 epidemic and the gradual improvement of urbanization.Many positive coefficient points of EVI (|z-score| > 1.96) appeared in the zoom-in area of Figure 7b, one of the fastest-urbanizing areas in 2019, while the number has decreased since then (Figure 7c,h,k).As population growth increased TCEa in fasturbanizing areas, more natural surfaces were replaced by hard materials with high heat capacity and thermal conductivity (such as concrete, asphalt, etc.) [46], further amplifying TCEa.In 2021, urbanization in the areas reached a certain saturation point, and the impact of NDBI may reduce as few buildings change the overall heat capacity of the ground surface [47].Therefore, in future climate mitigation and adaptation planning, more attention may be paid to the retention and increase of green spaces and water bodies, as well as the improvements in energy efficiency and measures to reduce heat emissions.However, the slight increase in points in this region in 2022 may be related to increased urban construction activities and population mobility after the COVID-19 epidemic.(3) MNDWI substantially mitigated TCEa, especially in densely populated and built-up urban centers.Some significant points appeared in the zoom-in area of Figure 6b,e,h,k (West Lake, a city park in downtown Fuzhou) in 2019, indicating a good effect on mitigating TCEa in central urban areas.Moreover, as the surrounding TCEa increased, significant points increased rather than decreased.This illustrated that water bodies such as lakes and rivers had continuous and stable cooling effects on areas with dense populations and buildings.However, it is worth noting that the inaccessibility of water bodies may exaggerate its impact, as this will reduce TCEa.
In addition to these exciting findings, some constructive measures to the inequalities of TCEa can also be derived through the mitigation strategy analysis based on the spatiotemporal heterogeneity relationships: (I) Artificial lakes could be built in rapidly urbanizing areas such as the east and southeast Cangshan.The area's increasing number of buildings and manufactured surfaces (such as sidewalks) promoted the surface's ability to absorb heat, exacerbating the UHI intensities.Water bodies (e.g., lakes and rivers) have the best effect on mitigating TCEa.Moreover, lakes can make more manufactured surface fragmentation and enhance urban landscapes.(II) Vegetation coverage could be increased in the southeastern part of the study area (the east region of Gaogaishan Park).Because the TCEa in the area was much higher and denser than others, many significant negative EVI coefficient points clustered, indicating strong mitigation effects.(III) For eastern areas of Taijiang, new buildings should be designed and constructed using higher reflectivity and more effective heat dissipation materials [48].Moreover, manufacturers and socio-economic activities in this region would be recommended to improve energy efficiency and reduce heat emissions.We recommend integrating or adding natural elements such as vertical green plants [49], green spaces, and water bodies for shopping malls in this region.
Although this study achieved the above-mentioned exciting findings, it still had the following challenges: (1) The data used in this study area would affect the accuracy and reliability of the results to some extent.First, the choice of image observation time would affect the calculated LST.Ref. [50] reported that the relationship between the land-surface temperature and land cover is strong in summer.This study uses almost cloud-free image data of Fuzhou in the summer (June-September).For years with multiple images, such as in 2017, the maximum value synthesis method was used to minimize the impact of clouds and other pollutants [32].However, in 2018, 2021, and 2022, we only collected one cloud-free image in those years, and the estimated values may be relatively lower, but these effects should be limited and are not expected to have significant impacts on the results because the image data are highly correlated without obvious outliers.Second, the TCEa values calculated based on Landsat 8 Collection 2 Level 2 may be slightly higher compared with using the relatively high-precision and widely recognized MODIS data.Although the calculated values from the two data products are highly correlated, there are still certain deviations in the results due to reasons such as spatial resolution and observation angle.The Landsat surface temperature data will be slightly higher than the MODIS data, especially for urban areas in summer [51], resulting in an overestimation of the heat island intensity and TCEa, therefore exaggerating the degree of thermal comfort exposure in the city.However, this influence is relatively limited.MODIS datasets are more suitable for thermal environment analysis in large-scale, high-frequency, and relatively uniformly covered areas.Since the scope of this study is small and the spatial changes in surface cover are frequent, the Landsat dataset, which is more suitable for small-scale and low-frequency research [52], is used.Future research can consider combining Landsat and MODIS product data to reduce the impact of Landsat data overestimation and improve the accuracy of the results.(2) Some preprocessing procedures, including radiance calibration [53], atmospheric correction [54], geometric correction [55], cropping, etc., were performed on the collected remote-sensing images.All the preprocessed images and population data were resampled to a spatial resolution of 240 m, according to [33].However, the processed remote-sensing data and downscaled population data may lose the details of the original data, resulting in reduced accuracy of local regional analysis results.For example, for the areas close to the inner boundaries of water bodies, the population data will increase due to resampling, resulting in an overestimation of TCEa.In addition, the local temperature extreme value estimate will be reduced, resulting in an underestimation of the extreme thermal comfort in the area.Moreover, we also performed Z-Score to standardize the surface parameters and TCEa to eliminate potential biases between datasets and ensure the stability and comparability of the results [56].(3) The land-surface and built-environment parameters also have some limitations.First, the current EVI can only be used to measure vegetation cover and cannot distinguish different types of vegetation.The current EVI can only be used to measure vegetation cover and cannot distinguish different types of vegetation.This is because the same EVI does not necessarily mean the same vegetation type, and the same vegetation type does not necessarily have the same EVI [57].Second, the NDBI cannot be used to distinguish between different building materials.Different building (specific heat capacity) materials have different effects on TCEa.Building materials with high specific heat capacity can absorb and store more heat, which can potentially help mitigate the urban heat island effect [58], while using low-specific heat capacity materials may result in higher LST [59].This is because high-specific-heat-capacity materials, such as concrete, brick, and stone, can gradually absorb heat during the daytime and slowly release heat during nighttime [60].Moreover, Fuzhou City is still experiencing accelerated urbanization.The building materials are being replaced rapidly under the ongoing "replace the old with new" project, making it challenging to analyze their impact in different locations.Further research could also consider the impact of regional building materials and related policies.Third, the MNDWI cannot be used to identify the water body types (such as rivers versus lakes), water quality, and seasonal water level, which may also have great impacts on TCEa.
(4) Both the choice of spatial and temporal resolution will influence the calculations of LST and TCEa.On the one hand, higher spatial resolution will result in more refined TCEa values, as well as more detailed (localized) spatial varying coefficient surfaces of different independent variables generated by the STWR model.On the other hand, the higher temporal resolution of collected original remote-sensing data would result in more accurate and robust annual land-surface and built-up parameters.If we shorten the time interval of our study, for example, based on monthly rather than yearly data, the temporal bandwidths would be more precise, and more coefficient surfaces would be generated by STWR for exploring higher dynamic geospatial processes.Moreover, some of the delta observation time was less than a year, which also had a certain impact on the results of the STWR model.Although the STWR model can utilize observation data at different time intervals, this study treated the calculated land surface and built-up parameters (i.e., EVI, NDBI, and MNDWI) as annual data, ignoring the delta-t of original remote-sensing data less than a year.This may result in slightly higher estimated coefficients by the STWR model because the increase of different change rates would lead to higher local temporal weights.(5) Exploring the multiscale relationships between the land-surface parameters and TCEa in space and time remains challenging.To the best of our knowledge, there are no GWR-based models with multiple spatial and temporal bandwidths simultaneously.The current STWR only had one single spatial and one temporal bandwidth.The advanced multiscale spatiotemporal model still needs to be further developed.

Conclusions
This research proposed a new proactive and analytical framework that integrated the STWR and population-weighted exposure models to measure TCEa and analyze mitigation strategies.Using this framework made scholars follow the study outlines "population growth → urbanization → land cover change → climate change → influence analysis (including the possible impact of COVID-19)" clearer.We also conducted an empirical study in Fuzhou, China.Some conclusions can be drawn as follows: (1) The effectiveness of this framework was verified.The overall impacts of various land-surface parameters and built-environment factors on TCEa were consistent with most studies (i.e., TCEa negatively correlated with EVI and MNDWI and positively correlated with NDBI), and delicate changes in local spatial relationships at different time stages can be detected using the STWR model.It can provide complete details of spatiotemporal heterogeneity and help identify underlying patterns in TCEa process changes.
(2) The spatiotemporal heterogeneous responses of TCEa to land-surface parameters were vital in analyzing mitigation strategies.This is because TCEa at different times and locations respond differently to various land-surface parameters and are probably the opposite.Mitigation strategies must be adjusted according to time and place to promote effective improvement of local TCEa.(3) Regarding the strategies for mitigating TCEa, the negative impacts of MNDWI were continuous and stable.In contrast, the impacts of EVI and NDBI were easily affected by the fluctuation of TCEa to a certain extent.Integrating optimization methods used in local spatial locations at different time stages is worth consideration and recommendation.
In addition, due to the complexity and dynamic characteristics of TCEa, data availability, etc., this empirical study may have some limitations, including the small scope of the study area and the mere inclusion of three land-surface parameters, which may affect the reliability of the results.Future studies could consider more land-surface parameters or other drivers for a comprehensive analysis of the urban environment based on this framework.

Figure 1 .
Figure 1.The study area is in Fuzhou City, China.

Figure 3
Figure 3 shows the spatial distribution of population, UHI, and thermal comfort exposure (TCEa) in the study area.The TCEa illustrated an overall pattern of low in the west and high in the east.Most high TCEa locations were mainly concentrated in dense population areas, commercial activities areas like eastern Gulou and central Taijiang, and industrial areas like the east Cangshan district.Some interesting phenomena can be found: (I) Although the UHI intensities in eastern Gulou, central Taijiang, and Jinan were not the highest, the highest TCEa appeared due to the dense population in these areas, especially in 2017.(II) Despite the relatively high UHI intensities in western Cangshan, TCEa was not high due to the relatively low population density.(III) The TCEa range in eastern Cangshan was extensive and relatively high for a long time.In contrast, for western Cangshan, the TCEa was relatively low and increased, mainly contributing to population growth, especially since 2020.

Figure 3 .
Figure 3. Spatial distributions of UHI intensity, population, and adjusted thermal comfort exposure from 2017 to 2022.

Table 1 .
Calculation of LST and land-surface parameters based on L8C2L2 data.

Table 1 .
Calculation of LST and land-surface parameters based on L8C2L2 data.

Table 2 .
Results of multicollinearity diagnostics.

Table 3 .
Comparisons of the model performances of OLS, GWR, and STWR.