Dynamic Multilevel Spatiotemporal Relationships between Heatwaves and Environmental Drivers as well as Emotional Responses in Chinese Urban Agglomerations

1. Introduction

Since the early 21st century, the combined impacts of global warming and rapid urbanization have significantly increased frequency, intensity, and spatial extent of extreme high-temperature events in urban areas, rendering heatwaves one of the most critical environmental challenges globally [1,2]. The World Health Organization (WHO) has recognized high temperatures as among the most severe environmental health threats of the 21st century, emphasizing that prolonged heat exposure can cause systemic physiological damage, including immunosuppression and metabolic disorders [3,4]. The combined effects of elevated temperatures and urbanization have intensified the duration, spatial reach, and extremity of heatwaves [5], posing significant risks to human habitats and public health [6]. The rapid advancement of urbanization has made urban agglomerations the "hotspots" of heatwave events by altering land cover, increasing heat emissions and intensifying the urban heat island effect [7,8,9]. The complex thermal environments and vulnerable socioeconomic structures within these densely populated urban agglomerations further amplify the health risks associated with heatwaves [10,11]. Therefore, investigating the interactions between heatwaves and the urban agglomeration environment, as well as their multidimensional impacts on human health, is of critical importance for developing effective strategies to mitigate the adverse effects of heatwaves.

A substantial body of research has explored the complex interplay between heatwaves and urban environmental factors. These studies encompass both natural environmental variables—such as urban heat island effects, land use change, vegetation cover, hydrological conditions, and topography [7,12,13], and socioeconomic factors, including population density, economic development, and infrastructure distribution [8,14]. The spatial distribution and intensity of heatwaves are shaped by the joint influence of multiple natural and anthropogenic factors across varying spatiotemporal scales, resulting in pronounced spatiotemporal variability [5,15]. Remote sensing technologies, with their extensive spatial coverage, long time series, and high spatiotemporal resolution, have become essential for monitoring, assessing, and elucidating the spatiotemporal dynamics of heatwaves [16,17]. For instance, land surface temperature (LST) data derived from remote sensing imagery are widely used to quantify the spatiotemporal characteristics of heatwaves and to reveal the specific impacts of urbanization on heatwave intensification[18,19]. Vegetation indices (e.g., Fractional Vegetation Cover [FVC], Normalized Difference Vegetation Index [NDVI]) and land use/land cover (LULC) data further facilitate analyses of the relationships among vegetation degradation, land use change, and heatwave exposure [20,21].

However, most existing studies are constrained to a single spatiotemporal section and lack systematic time series analyses, limiting the ability to comprehensively capture the dynamic evolution of heatwaves [11,22]. Research on the mechanisms by which environmental factors influence heatwaves often overlooks issues of multicollinearity and has limited capacity to disentangle multilevel, cross-scale drivers, thereby constraining explanations of the complex origins of heatwave exposure [9,14,23,24]. Geographical Convergent Cross Mapping (GCCM) model, grounded in causal inference theory, quantifies dynamic information flows among variables and effectively detects causal effects in multivariate systems, overcoming the confounding effects of multicollinearity [25]. This enables robust elucidation of the dynamic causality between environmental factors and heatwave evolution. The Bayesian Generalized Linear Mixed Model (Bayesian GLMMs) offers a flexible hierarchical structure to integrate cross-scale and multilevel data [26,27], allowing for multidimensional analysis of the drivers of heatwave exposure risk. The integration of these two methods provides a scientifically rigorous analytical framework for investigating the spatiotemporal dynamics of relationships between heatwave and environment.

In the domain of heatwave impacts on human health, prior research has primarily focused on physiological health outcomes, resulting in a relatively mature theoretical framework [28]. In contrast, studies on the psychological health impacts of heatwaves, particularly on emotional responses, remain limited. Due to the challenges of acquiring large-scale mental health data, most existing studies rely on small-sample surveys or clinical experiments, restricting effective research samples to narrow spatiotemporal scopes and hindering continuous population-level monitoring [29,30]. Moreover, conventional health indicators are inadequate for capturing the temporal dynamics of emotional fluctuations [31]. Consequently, the cross-scale interactive effects between mental health and environmental climate events—especially the temporal lag effects of emotional responses to heatwaves—remain insufficiently explored [32,33], and their underlying mechanisms are not fully understood [11].

Emotional responses constitute a critical dimension of individual and societal adaptation to heatwaves. During heatwave periods, the public may experience heightened negative emotions such as anxiety and depression [30,34,35], which not only affect individual psychological states but can also alter social behaviors and, in turn, influence societal resilience [23,36]. The proliferation of internet and social media platforms has created new avenues for individuals to express opinions and emotions [37]. Social media big data, with its extensive temporal coverage and broad spatial reach, has become an important resource for studying public emotional dynamics in response to both social and environmental stimuli [38]. Leveraging social media big data to investigate emotional responses to heatwave events can effectively address limitations of small sample sizes, short observation periods, and restricted spatial scales in traditional emotion research, providing a promising approach for examining the emotional impacts of heatwaves.

As one of the world’s most populous countries, China faces particularly acute challenges in coping with extreme climate events due to the rapid development and high population density of its urban agglomerations [14]. Understanding the spatiotemporal mechanisms linking heatwaves, urban environmental variables, and emotional responses in Chinese urban agglomerations is essential for optimizing urban spatial structure and developing effective urban management strategies to mitigate the impacts of heatwaves. This study focuses on China’s five major urban agglomerations, aiming to quantify the spatiotemporal heterogeneity of heatwave exposure, analyze dynamic spatiotemporal patterns of public emotional responses using social media data, and elucidate the underlying mechanisms. The main research objectives are: (i) to construct heatwave exposure indicators and environmental factors for urban agglomerations from 2014 to 2023 using remote sensing data, and to employ GCCM to analyze the causality between heatwaves and environmental factors, identifying key environmental drivers; (ii) to apply Bayesian GLMMs to elucidate the driving mechanisms of heatwave exposure risk; and (iii) to construct a dynamic sentiment index based on geotagged Weibo data, and, in combination with modeling, to reveal the association patterns between heatwave events and public emotional fluctuations. This research aims to address gaps in the emotional dimension of heatwave studies, provide new perspectives on the multilevel impacts of extreme heat events, and offer scientific evidence for public health planning and climate adaptation strategies.

2. Data

2.1. Study Area

This study focuses on the five major urban agglomerations in China: the Beijing–Tianjin–Hebei Urban Agglomeration (BTHUA), Chengdu–Chongqing Urban Agglomeration (CCUA), Triangle of Central China (TCC), Yangtze River Delta Urban Agglomerations (YRDUA), and Pearl River Delta Urban Agglomeration (PRDUA). These regions represent the most densely populated and economically developed urban clusters in China, occupying a central role in the country’s urbanization strategy [39]. Their geographic diversity and distribution across China’s north-south and east-west axes provide a comprehensive perspective for investigating the spatiotemporal evolution of heatwaves and associated emotional responses.

2.2. Data Sources

The primary datasets used in this study include temperature data, remote sensing data, population data, and social media data, as detailed in Table 1.

2.2.1. Temperature Data

This Research utilized the temperature_2m band from the ERA5-Land dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) [40] to extract heatwave-related indicators. ERA5-Land, derived from the ERA5 reanalysis, offers a globally consistent climate dataset with high spatial resolution (0.1°), providing an accurate representation of temperature variability across China.

Figure 1. (a) 2023 Estimates of population distribution in China, (b) Major urban agglomerations in China.

Figure 1. (a) 2023 Estimates of population distribution in China, (b) Major urban agglomerations in China.

2.2.2. Remote Sensing Data

Remote sensing data from MODIS [41], Landsat 8 [42], Copernicus DEM [43], and VNL V2 [44], along with their derived products (Table 1), were used to generate a suite of natural environmental factors, including surface Albedo, Modified Normalized Difference Water Index (MNDWI), Normalized Difference Thermal Index (NDTI), Elevation, Fractional Vegetation Cover (FVC), and Nighttime Light data (NTL). All data processing was conducted on the Google Earth Engine platform.

2.2.3. Population and Administrative Boundary Data

Population distribution data were obtained from the LandScan dataset [45], which provides high-resolution (30 arc-second) global population grids. Figure 1 (a) presents the estimated population distribution of China in 2023 in the LandScan dataset. Administrative boundary data were sourced from the Ministry of Natural Resources of the People’s Republic of China to delineate the spatial extents of the five urban agglomerations.

2.2.4. Social Media Data

Using a Python 3.8-based web crawler, we collected climate change–related posts from core cities within the five urban agglomerations on Sina Weibo (one of China’s most active social media platforms) between 2019 and 2023. After filtering out irrelevant or insufficiently detailed posts, approximately 3.2 million original Weibo entries were retained for sentiment analysis.

Table 1. Datasets used in this study.

Datasets	Environmental factors	Time span	Source
ERA5-Land	Temperature	2014-2023	https://cds.climate.copernicus.eu
Administrative region		2014-2023	https://www.mnr.gov.cn/sj/sjfw/
LandScan	Population	2014-2023	https://landscan.ornl.gov
MODIS Albedo	Albedo	2014-2023	https://lpdaac.usgs.gov/products/mcd43a3v061/
Landsat 8 Collection 2 Tier 1 Raw Scenes	Modified Normalized Difference Water Index (MNDWI)[46]	2014-2023	https://www.usgs.gov/land-resources/nli/landsat/landsat-8-data-users-handbook
	Normalized Difference Thermal Index (NDTI)[47]
Copernicus DEM	Elevation	2015	https://dataspace.copernicus.eu/collections/copernicus-digital-elevation-model
MOD13Q1 V6.1	Fractional Vegetation Cover (FVC)	2014-2023	https://lpdaac.usgs.gov/products/mod13q1v061/
VNL V2	Nighttime lights (NTL)	2014-2023	https://eogdata.mines.edu/products/vnl/

Table 1. Datasets used in this study.

Datasets	Environmental factors	Time span	Source
ERA5-Land	Temperature	2014-2023	https://cds.climate.copernicus.eu
Administrative region		2014-2023	https://www.mnr.gov.cn/sj/sjfw/
LandScan	Population	2014-2023	https://landscan.ornl.gov
MODIS Albedo	Albedo	2014-2023	https://lpdaac.usgs.gov/products/mcd43a3v061/
Landsat 8 Collection 2 Tier 1 Raw Scenes	Modified Normalized Difference Water Index (MNDWI)[46]	2014-2023	https://www.usgs.gov/land-resources/nli/landsat/landsat-8-data-users-handbook
	Normalized Difference Thermal Index (NDTI)[47]
Copernicus DEM	Elevation	2015	https://dataspace.copernicus.eu/collections/copernicus-digital-elevation-model
MOD13Q1 V6.1	Fractional Vegetation Cover (FVC)	2014-2023	https://lpdaac.usgs.gov/products/mod13q1v061/
VNL V2	Nighttime lights (NTL)	2014-2023	https://eogdata.mines.edu/products/vnl/

3. Methods

3.1. Methodological Framework

The methodological framework comprises three main steps (Figure 2):

(i) Multi-source remote sensing, meteorological, demographic, and topographic data were processed through the GEE platform to generate environmental factors including FVC and NDTI, as well as heatwave indicators such as intensity, frequency, and duration. Simultaneously, social media data were acquired using Python web crawling techniques, and sentiment indices were constructed based on Baidu's open-source NLP Senta model to quantify public emotional responses to heatwave events.

(ii) Geographic Convergent Cross Mapping (GCCM) was applied to quantitatively assess the bidirectional causality strength between heatwave indicators and environmental factors and their spatiotemporal variation characteristics, thereby identifying causality between environmental factors and heatwave indicators.

(iii) Bayesian Generalized Linear Mixed Model (Bayesian GLMMs) was constructed to integrate heatwave indicators, selected environmental factors, and sentiment indices into a unified analytical framework, exploring the complex interactions and potential moderating effects between heatwaves and environment as well as sentiment, revealing the multilevel spatiotemporal response mechanisms of public sentiment to heatwaves, and providing scientific support for urban thermal environment management and public health protection.

3.2. Heatwave Indicators

3.2.1. Excess Heat Factor (EHF) Calculation

Bulleted lists look like this: The Excess Heat Factor (EHF) [1] was adopted as the primary metric for heatwave events, with several derived indicators to characterize heatwave features. EHF is calculated as follows:

$$\begin{array}{c}{EHF}_{sig}={\displaystyle \frac{T_i+T_{i-1}+T_{i+1}}{3}}-T_{95}\end{array}$$

(1)

$$\begin{array}{c}{EHF}_{accl}={\displaystyle \frac{T_i+T_{i-1}+T_{i+1}}{3}}-{\displaystyle \frac{T_{i-2}+,...,+T_{i-29}}{30}}\end{array}$$

(2)

$$\begin{array}{c}EHF={EHF}_{sig}\times max\left(1,{EHF}_{accl}\right)\end{array}$$

(3)

where $T$ is the daily mean temperature (average of daily maximum and minimum); ${\displaystyle \frac{T_i+T_{i-1}+T_{i+1}}{3}}$ is the average temperature of three consecutive days；$i$ represents the i-th day of each year，$T_{95}$ is the 95th percentile value of the daily temperature distribution at a certain location from 1981 to 2010, which is considered as the threshold for extreme high temperatures in this area；${\displaystyle \frac{T_{i-2}+,...,+T_{i-29}}{30}}$ is the average temperature of the past 30 days.

Figure 2. Workflow of this study.

Figure 2. Workflow of this study.

3.2.2. Heatwave Event Definition

In this study, a high temperature weather event is defined as an occurrence where $EHF>0$ , and a heatwave event is identified when $EHF>0$ persists for at least three consecutive days [48]. To more comprehensively reflect the severity of heatwave events, multiple heatwave indicators recommended by the WMO Commission for Climatology [49] are adopted. These indicators are widely used to describe the frequency, duration, and intensity of heatwave events.

Furthermore, to further distinguish the severity levels of high-temperature events, ${EHF}_{85}$ is introduced, representing the 85th percentile of all EHF values greater than zero at a given location during 2014–2023, arranged in ascending order, and used as the threshold for moderate heat events. When $0<EHF<{EHF}_{85}$ , the event is defined as a mild heat event (HW_day); when ${EHF}_{85}\le EHF<2{EHF}_{85}$ , it is defined as a moderate heat event (HW_sev); and when $EHF\ge {2EHF}_{85}$ , it is defined as a severe heat event (HW_ex). By accumulating EHF values during heatwave events, the annual total EHF (HW_a) is obtained to measure the overall intensity of annual heatwave events.

All heatwave-related indicators were calculated using the Google Earth Engine platform. Figure 3. (a) shows the spatial distribution of $T_{95}$ across China in this study, while Figure 3. (b)–(j) present the distributions of nine heatwave indicators in the five urban agglomerations.

Figure 3. (a) Spatial visualization of the 95th percentile value of the daily temperature distribution at a certain location from 1981 to 2010; (b)-(j) are HWN, HWD, HWF, HWa, HWM, HWday, HWsev, HWex, and HWa respectively. The five columns from left to right in each sub-graph are BTHUA, CCUA, TCC, YRDUA and PRDUA. From top to bottom are for different years from 2023 to 2014.

Figure 3. (a) Spatial visualization of the 95th percentile value of the daily temperature distribution at a certain location from 1981 to 2010; (b)-(j) are HWN, HWD, HWF, HWa, HWM, HWday, HWsev, HWex, and HWa respectively. The five columns from left to right in each sub-graph are BTHUA, CCUA, TCC, YRDUA and PRDUA. From top to bottom are for different years from 2023 to 2014.

3.2.3. Population Exposure

To assess the impact of heatwaves on emotional responses, we estimated total population exposure using a grid-based approach[50]:

$$\begin{array}{c}PE=HW\times POP\end{array}$$

(4)

where $HW$ is the heatwave indicator and $POP$ is the total population in a given area, $PE$ is the total amount of population exposure to heatwaves in this area.

3.3. Geographical Convergent Cross Mapping

When exploring the influencing factors of a group of variables, the direction of causality and mirroring effect should be considered. To uncover the driving mechanisms of heatwave exposure within the five major Chinese city clusters, we employ the Geographical Convergent Cross Mapping (GCCM) model[25]. GCCM is a state-space reconstruction technique designed for causal inference from spatial cross-sectional data, which is particularly well suited for Earth system analyses where time series data may be limited or exhibit only subtle variations. In addition, GCCM performs cross-mapping prediction based on the spatial section data in the reconstructed state space, which overcomes the mirroring effect well and can fully mine the causality of the spatio-temporal section data in the Earth system (Figure 4).

3.3.1. Spatial Embedding and State-Space Reconstruction

GCCM builds on the principles of Takens’ embedding theorem and the generalized embedding theorem. GCCM uses spatial lags—values observed at a focal spatial unit and its neighbors—to construct a “shadow manifold” for each variable. Formally, for a spatial variable $X$ observed over discrete spatial units, the embedding is defined as:

$$\begin{array}{c}\psi (x,s)=\{h_s(x),h_{s\left(1\right)}(x),\dots ,h_{s(L-1)}(x\left)\right\}\end{array}$$

(5)

where $h_s\left(x\right)$ is the observation at the focal unit and $h_{\mathrm{s}\left(\mathrm{i}\right)}\left(x\right)$ represents the $i$th-order spatial lag.

Figure 4. (a) Lorenz system: Nonlinear governing equations and attractor with a butterfly-like structure. Relationships between original and shadow attractors are illustrated. (b) Cross-mapping principle: Nearest neighbors are mapped to determine causal information. If nearest neighbors remain after mapping, causal information exists (top); otherwise, they do not remain nearest (bottom). Mapping direction and causal information are inverse. (c) Phase Space Reconstruction (PSR): Reconstructs a diffeomorphic shadow attractor using a single time series from the original attractor. Spatial lags: Focal unit and spatial lags of different orders in raster data (first-order: yellow, second-order: orange, third-order: red). (d) Causal inference via convergence: GCCM results (red line) indicate causality if convergent; non-convergence (blue line) indicates no causality. Mutual causality occurs if both lines converge, while unidirectional causality occurs if only one line converges. Two non-convergent lines suggest variables belong to separate systems.

Figure 4. (a) Lorenz system: Nonlinear governing equations and attractor with a butterfly-like structure. Relationships between original and shadow attractors are illustrated. (b) Cross-mapping principle: Nearest neighbors are mapped to determine causal information. If nearest neighbors remain after mapping, causal information exists (top); otherwise, they do not remain nearest (bottom). Mapping direction and causal information are inverse. (c) Phase Space Reconstruction (PSR): Reconstructs a diffeomorphic shadow attractor using a single time series from the original attractor. Spatial lags: Focal unit and spatial lags of different orders in raster data (first-order: yellow, second-order: orange, third-order: red). (d) Causal inference via convergence: GCCM results (red line) indicate causality if convergent; non-convergence (blue line) indicates no causality. Mutual causality occurs if both lines converge, while unidirectional causality occurs if only one line converges. Two non-convergent lines suggest variables belong to separate systems.

3.3.2. Cross-Mapping Prediction

Once the state spaces (manifolds) for two spatially distributed variables (e.g., a heatwave indicator and a candidate environmental factor) are reconstructed, GCCM employs cross mapping to quantify the causality. In this framework, if two variables $X$ and $Y$ are dynamically linked, the manifold reconstructed from $X$ (denoted $M_X$) can be used to predict the corresponding states in the manifold of $Y$ ($M_Y$). Specifically, for a given spatial unit $s$, the value of $Y$ is estimated as:

$$\begin{array}{c}{\widehat{Y}}_s|M_X={\displaystyle \sum _{i=1}^{L+1}}{\omega }_{s_i}{Y}_{s_i}|M_X\end{array}$$

(6)

where the weights ${\omega }_{s_i}$ are computed based on the distances between the focal state $\psi (x,s)$ and its $L+1$ nearest neighbors in $M_X$ :

$$\begin{array}{c}{\omega }_{s_i}|M_X={\displaystyle \frac{\exp \left[-\frac{dis\left(\psi \right(x,s_i),\psi (x,s\left)\right)}{dis\left(\psi \right(x,s_1),\psi (x,s\left)\right)}\right]}{{\sum }_{i=1}^{L+1}\exp \left[-\frac{dis\left(\psi \right(x,s_i),\psi (x,s\left)\right)}{dis\left(\psi \right(x,s_1),\psi (x,s\left)\right)}\right]}}\end{array}$$

(7)

The distance between two points in the shadow manifold is defined by:

$$\begin{array}{c}dis(\psi (x,s_i),\psi (x,s))={\displaystyle \frac{1}{L}}\left(\left|h_{si}\left(x\right)-h_s\left(x\right)\right|+{\displaystyle \sum _{k=1}^{L-1}}abs\left[h_{si\left(k\right)}\left(x\right),h_{s\left(k\right)}\left(x\right)\right]\right)\end{array}$$

(8)

where the term $abs\left[h_{si\left(k\right)}\left(x\right),h_{s\left(k\right)}\left(x\right)\right]$represents the absolute difference between the $k$th-order spatial lags of the focal unit and its neighbor. In the data based on this paper, the absolute difference function is defined as follows:

$$\begin{array}{c}abs\left[h_{si\left(k\right)}\left(x\right),h_{s\left(k\right)}\left(x\right)\right]={\displaystyle \frac{1}{D}}\sum _d^D\left|u_{si(k,d)}\left(x\right)-u_{s(k,d)}\left(x\right)\right|\end{array}$$

(9)

where $u_{s(k,d)}\left(x\right)$ is the observed value in the $d$th direction of the $k$th-order lag, and $D$ is the total number of directions considered.

3.3.3. Evaluating Prediction Skill and Inferring Causality

The skill of the cross-mapping prediction is evaluated using the Pearson correlation coefficient between the observed and predicted values:

$$\begin{array}{c}\rho ={\displaystyle \frac{Cov(Y,\widehat{Y})}{\sqrt{Var\left(Y\right)Var\left(\widehat{Y}\right)}}}\end{array}$$

(10)

By increasing the library size (i.e., the number of spatial units used in the reconstruction), convergence of $\rho$ is examined. A statistically significant convergence indicates a robust causal association between $X$ and $Y$. In particular, if the cross mapping from $X$ to $Y$ (denoted $X\to Y$) yields a significantly higher $\rho$ than the inverse mapping ($Y\to X$), the dominant causal direction is inferred as $X\to Y$.

In our study, causality determination is based on the convergence behavior of the two mapping curves as the library size (i.e., the number of spatial units used for state space reconstruction) increases. If one mapping curve converges while the other remains low or does not converge, the direction corresponding to the convergent curve is considered the dominant causal direction, indicating a unidirectional causality. If both mapping curves converge as the library size increases, the relationship is interpreted as bidirectional. If neither mapping curve converges or both remain at low levels, no causality is inferred. A convergent curve reaching a high correlation coefficient indicates strong causality, while convergence within a moderate range suggests weak causality.

3.4. Senta Model

he open-source Senta model provided by Baidu NLP National Artificial Intelligence Open Innovation Platform is based on deep learning technology, utilizing a pre-trained BERT model as its foundation to efficiently and accurately perform sentiment analysis on text. First, unstructured Chinese text is preprocessed through tokenization and extraction of key elements (such as target, opinion, holder, and time), which are crucial for capturing contextual sentiment. Subsequently, the Baidu NLP platform’s API is used to call the Chinese sentiment lexicon, mapping lexical cues to sentiment polarity. In typical scenarios involving subjective descriptive text, a high-precision deep neural network (DNN) semantic model, integrated with an extensive sentiment lexicon (Supplementary S1), is employed to capture the hedonic state of individual emotions in Weibo updates. The DNN model computes rich semantic representations of the input text and outputs a set of indicators to facilitate sentiment determination: the probability of the text being classified as positive, the probability of being classified as negative, the overall sentiment polarity label, and the associated confidence scores for the classification. This multi-output result not only quantifies the degree of emotion but also ensures robust handling of emotional cues in the dynamic social media environment.

3.5. Bayesian Generalized Linear Mixed Models

This study employs Bayesian Generalized Linear Mixed Models (Bayesian GLMMs) to analyze the driving mechanisms of heatwave exposure risk. GLMMs can simultaneously handle fixed and random effects, and by appropriately selecting the probability distribution and link function for the response variable, they can address the challenges of non-normal data. Within the Bayesian framework, prior distributions are specified and Markov Chain Monte Carlo (MCMC) methods are used to sample from the posterior distribution, providing uncertainty estimates for model parameters and a solid basis for subsequent inference.

3.5.1. Model Structure

All continuous covariates (heatwave exposure indicators $H_{it}^{*}$ , environmental factors $E_{it}^{*}$ , and sentiment indices $S_{it}^{*}$ ) are first centered and standardized by their overall sample standard deviation:

$$\begin{array}{c}{X}_{it}={\displaystyle \frac{X_{it}^{*}-\overline{X}}{SD\left(X\right)}}\end{array}$$

(11)

where $\overline{X}$ and $SD\left(X\right)$ denote the mean and standard deviation of the variable across all observations, respectively.

Let the response variable $y_i$ follow an appropriate distribution, with its mean linked to a linear combination of covariates and random effects via a link function. The general form of the model is:

$$\begin{array}{c}g({\mu }_i)=x_i^{\top }\beta +z_i^{\top }ib\end{array}$$

(12)

where ${\mu }_i=E\left(y_i\right)$; $g(\cdot )$ is a known link function; $x_i$ is the vector of fixed-effect covariates for the i-th observation, $\beta$ is the corresponding vector of regression coefficients; $z_i$ is the covariate vector associated with random effects, and $b$ follows a multivariate normal distribution $b\sim \mathcal{N}\left(0,G\right)$ , where $G$ is the covariance matrix of the random effects.

3.5.2. Priors and Posterior Distributions

Weakly informative priors are assigned to the fixed-effect parameters $\beta$ , such as normal distributions ${\beta }_j=\mathcal{N}\left(0,{\sigma }_{\beta }^2\right)$ with large variance to reflect prior uncertainty. For the covariance matrix $G$ of the random effects, an inverse Wishart prior is used to ensure model flexibility and numerical stability. Based on the above priors and the likelihood function, the posterior distribution is:

$$\begin{array}{c}p(\beta ,b,G\mid y)\propto {\displaystyle \prod _{i=1}^n}p\left(y_i\left|x_i^{\top }\beta +z_i^{\top }ib\right.\right)p\left(\beta \right)p\left(b\left|G\right.\right)p\left(G\right)\end{array}$$

(13)

where $y_i$denotes the response for the i-th observation; $x_i$ is the fixed-effect covariate vector with corresponding parameter $\beta$ ; $b$ is the covariate vector for random effects, with random effects $b\sim \mathcal{N}(0,G)$； $G$ is the covariance matrix of the random effects, with its prior$p\left(G\right)$ typically specified as an inverse Wishart distribution.

3.5.3. Model Validation

Markov Chain Monte Carlo (MCMC) methods are used to sample from the posterior. Appropriate numbers of iterations, burn-in periods, and thinning intervals are set, and diagnostic metrics (such as R², effective sample size (MAE), etc.) are used to assess chain convergence and mixing. After sampling, model convergence is diagnosed using posterior density plots, trace plots, and autocorrelation plots. Additionally, posterior predictive checks are performed by comparing simulated and observed heatwave indicators for 2024 to evaluate model fit. For each parameter, the 95% credible interval is calculated; if the interval does not include zero, the corresponding effect is considered statistically significant.

To facilitate the interpretation and direct comparison of the relative importance of different environmental predictors, we report standardized effect sizes for all fixed effects. The standardized effect size for each predictor is defined as the estimated regression coefficient (${\beta }_k$) when all predictors and the response variable have been standardized to have zero mean and unit variance. This approach allows the effect size to be interpreted as the expected change in the response variable (in standard deviation units) associated with a one standard deviation increase in the predictor, holding other variables constant.

3.5.4. Integrated Analysis

The sentiment scores of the five urban agglomerations from 2019 to 2023 were weighted by the population distribution of each Urban Agglomeration to generate annual sentiment indices for each region. These sentiment indices, together with selected environmental factors, were incorporated into the model as interaction terms with heatwave variables. The model is specified as follows:

$$\begin{array}{c}g\left({\mu }_{it}\right)={\beta }_0+{\beta }_1{H}_{it}+{\beta }_2{S}_{it}+{\beta }_3{E}_{it}+{\beta }_4\left(H_{it}\times S_{it}\right)+{\beta }_5\left(H_{it}\times E_{it}\right)+{\beta }_6\left(H_{it}\times E_{it}\times S_{it}\right)+b_i,\end{array}$$

(14)

where $H_{it}$ denotes the heatwave exposure indicator for Urban Agglomeration $i$ in year $t$ , $S_{it}$ is the sentiment index calculated and weighted based on geo-tagged Weibo data, $E_{it}$ represents the selected key environmental factors, and $g(\cdot )$ is the appropriate link function. The random intercept $b\sim \mathcal{N}(0,{\sigma }^2)$ captures the heterogeneity among urban agglomerations. ${\beta }_0$ is the fixed intercept, and ${\beta }_1$ to ${\beta }_6$ are the regression coefficients for heatwave indicators and environmental factors in the equation. This model structure allows for the simultaneous examination of the individual and interactive effects among heatwaves, sentiment, and environmental factors, providing a statistical basis for elucidating how these factors jointly influence heatwave exposure risk.

4. Results

4.1. Causality between Heatwave and Environment

GCCM analysis revealed that most relationships exhibited asymmetric bidirectional causality: environmental factors generally showed strong causality toward heatwave indicators, whereas heatwave indicators exhibited weaker or negligible causality toward environmental factors. Intensity indicators (HWA, HWM, HWA) showed the strongest response to environmental factors, frequency indicators (HWF, HW_day) demonstrated moderate to strong responses, and the sensitivity of duration indicators (HWD) was slightly lower than that of intensity indicators (Figure 5., Figure 6.). Albedo exhibited the strongest overall causality with heatwave indicators (mean ρ ≈ 0.772); FVC and NDTI showed lower causality with most heatwave indicators compared to Albedo, but were still relatively strong among all factors (mean ρ ≈ 0.673 and ρ ≈ 0.691, respectively). MNDWI and NTL exhibited moderate causality (mean ρ ≈ 0.475 and ρ ≈ 0.533, respectively), with a greater impact on intensity indicators (HWA, HWM) than on frequency indicators (HWN, HWF, HWday). Elevation showed the weakest causality with heatwave indicators (mean ρ ≈ 0.304).

Figure 5. Causality between heatwave indicators and environmental factors，The red arrow represents strong causality, the purple arrow represents medium causality, the blue arrow represents weak causality, and the gray arrow represents feedback under strong coupling relationship.

Figure 5. Causality between heatwave indicators and environmental factors，The red arrow represents strong causality, the purple arrow represents medium causality, the blue arrow represents weak causality, and the gray arrow represents feedback under strong coupling relationship.

Temporally, the mean maximum library prediction skill ρ for the causality between heatwave indicators and environmental factors in the five urban agglomerations remained generally stable across years, despite some fluctuations. The strength of the causality between the same environmental factor and different heatwave indicators exhibited varying trends over time (Figure. 7). The causality between Albedo and heatwave indicators remained highly stable and strong throughout the study period (ρ values for different heatwave indicators in all years were above 0.6). Only the HWA indicator showed a steadily increasing and then stabilizing trend with Albedo from 2014 to 2023 (mean ρ increased from 0.7023 to 0.8500, peaking in 2020), while the causality between other heatwave indicators and Albedo did not show regular temporal trends. The causal strength of MNDWI, NDTI, Elevation, FVC, and NTL with heatwave indicators remained generally stable over time, with only minor interannual differences and no clear temporal trends.

Spatially, the strength of the causality between environmental factors and heatwave indicators varied significantly among urban agglomerations. The causal strength between Albedo and heatwave indicators was significantly higher in BTHUA and CCUA (mean ρ values of 0.792 and 0.790, respectively) than in PRDUA and YRDUA, especially for HWN and HWF, where the mean ρ values of Albedo in BTHUA and CCUA both exceeded 0.8. The impact of MNDWI on heatwaves was much higher in YRDUA and PRDUA (mean ρ values of 0.6398 and 0.6598, respectively) than in BTHUA (mean ρ value 0.4636). NDTI showed a much stronger causality with heatwaves in PRDUA and CCUA (mean ρ values of 0.7488 and 0.7203, respectively) than in TCC (0.6288). The effect of Elevation on heatwaves was most pronounced in CCUA (mean ρ value 0.4591), much higher than in BTHUA (mean ρ value 0.3012), especially for heatwave intensity (HWA), where the prediction skill ρ value for Elevation in CCUA reached 0.5865. The causal strength of FVC on heatwaves was significantly higher in PRDUA and YRDUA (mean ρ values of 0.7787 and 0.7599, respectively) than in BTHUA (mean ρ value 0.6393), particularly for heatwave intensity (HWA) and heatwave days (HWday), where the ρ values for FVC in PRDUA and YRDUA both exceeded 0.8. The impact of NTL on heatwaves was most significant in BTHUA (mean ρ value 0.6052), higher than in CCUA (0.5051).

Figure 6. Causal inference for heatwave indicators and environmental factors (three different causal intensity relationships).

Figure 6. Causal inference for heatwave indicators and environmental factors (three different causal intensity relationships).

4.2. Bayesian GLMMs Results and Reliability

4.2.1. Fixed Effects of Environmental Factors

The results of the Bayesian GLMMs showed that the influence strength of selected environmental factors on heatwave indicators was highly consistent with the findings from the causality analysis (Figure 8, Figure 9(b)). Strong causality factors—Albedo (Standardized Effect: 0.58–0.82), FVC (0.55–0.78), and NDTI (0.53–0.76)—exhibited the strongest effects, especially on heatwave intensity indicators (HWA, HWM). Moderate causality factors—NTL (0.42–0.58) and MNDWI (0.38–0.52)—showed moderate effects, with relatively stronger impacts on heatwave frequency indicators (HWN, HWF, HW_day). Weak causality factor—Elevation (0.20–0.32)—had relatively weak effects on all heatwave indicators, with the smallest impact on heatwave severity indicators (HW_a, HW_sev, HW_ex).

The effect sizes of different environmental factors on heatwave indicators also varied. The order of effect size from strongest to weakest was: intensity indicators > frequency indicators > duration indicators > severity indicators (Figure 8). The standardized effect sizes of environmental factors exhibited a consistent hierarchical pattern across the five urban agglomerations (Figure 9). Among different urban agglomerations, the average effect strength of environmental factors on heatwave occurrence ranked as follows: BTHUA, PRDUA, YRDUA, TCC, and CCUA (Figure 9(c)). Notably, unlike other urban agglomerations where Albedo was the strongest influencing factor, FVC had the strongest effect on heatwaves in CCUA, and Elevation had a higher impact than in other urban agglomerations (Figure 9(a)). In PRDUA and YRDUA, NDTI had a stronger effect on heatwaves than FVC, second only to Albedo, whereas the opposite was true in other urban agglomerations. The water index (MNDWI) had the least impact on urban heatwaves in BTHUA but was significant in TCC, while nighttime lights (NTL) had the least impact in CCUA but was significant in BTHUA and PRDUA, showing moderate effect sizes. The effect of Elevation was relatively weak (0.359, 95% CI: 0.250–0.468) and did not reach high significance in any urban agglomeration.

Figure 7. Under the premise of convergence, when the largest library size is reached, the prediction skill ρ (keep four decimals) of heatwave indicators and environmental factors is obtained.

Figure 7. Under the premise of convergence, when the largest library size is reached, the prediction skill ρ (keep four decimals) of heatwave indicators and environmental factors is obtained.

Figure 8. Fixed effect analysis: The impact of environmental factors on heatwaves; (a) Standardized Effect Size of Environmental Factors on Heatwave Indicators, (b) Effect of Environmental Factors on Heatwave Indicators.

Figure 8. Fixed effect analysis: The impact of environmental factors on heatwaves; (a) Standardized Effect Size of Environmental Factors on Heatwave Indicators, (b) Effect of Environmental Factors on Heatwave Indicators.

The effect sizes of different environmental factors on heatwave indicators also varied. The order of effect size from strongest to weakest was: intensity indicators > frequency indicators > duration indicators > severity indicators (Figure 8). The standardized effect sizes of environmental factors exhibited a consistent hierarchical pattern across the five urban agglomerations (Figure 9). Among different urban agglomerations, the average effect strength of environmental factors on heatwave occurrence ranked as follows: BTHUA, PRDUA, YRDUA, TCC, and CCUA (Figure 9(c)). Notably, unlike other urban agglomerations where Albedo was the strongest influencing factor, FVC had the strongest effect on heatwaves in CCUA, and Elevation had a higher impact than in other urban agglomerations (Figure 9(a)). In PRDUA and YRDUA, NDTI had a stronger effect on heatwaves than FVC, second only to Albedo, whereas the opposite was true in other urban agglomerations. The water index (MNDWI) had the least impact on urban heatwaves in BTHUA but was significant in TCC, while nighttime lights (NTL) had the least impact in CCUA but was significant in BTHUA and PRDUA, showing moderate effect sizes. The effect of Elevation was relatively weak (0.359, 95% CI: 0.250–0.468) and did not reach high significance in any urban agglomeration.

Figure 9. Standardized Effect Size of Environmental Factors in Urban Agglomerations, (a) Standardized Effect Size of Environmental Factors in Urban Agglomerations, (b) Mean Effect Sizes of Environmental Factors, (c) Mean Effect Size by Urban Agglomerations.

Figure 9. Standardized Effect Size of Environmental Factors in Urban Agglomerations, (a) Standardized Effect Size of Environmental Factors in Urban Agglomerations, (b) Mean Effect Sizes of Environmental Factors, (c) Mean Effect Size by Urban Agglomerations.

Analysis of the average effects of all environmental factors (Figure 8(b)) indicated that surface physical properties (Albedo, FVC, NDTI) had a significantly stronger impact on heatwaves than hydrological and topographic features (MNDWI, Elevation). This result supports the urban heat island theory, highlighting that urban surface material properties and vegetation cover are key factors in regulating the urban thermal environment. There were slight differences in the responses of different urban agglomerations to environmental factors (Figure 8(c)), with BTHUA and PRDUA showing stronger average effect sizes, and TCC showing weaker effects.

Figure 10. Random effect analysis: Spatial and temporal heterogeneity; (a) Random Effects of Urban Agglomerations, (b) Temporal Random Effects of Heatwave Risk.

Figure 10. Random effect analysis: Spatial and temporal heterogeneity; (a) Random Effects of Urban Agglomerations, (b) Temporal Random Effects of Heatwave Risk.

4.2.2. Random Effects of Environmental Factors

There was significant spatial heterogeneity in the environmental effects on heatwave occurrence across the five urban agglomerations (Figure 9(a)), indicating that heatwave occurrence in different urban agglomerations was significantly influenced by random effects. The estimated random effects represent the deviation of each urban agglomeration from the national average. BTHUA (+0.2384) and YRDUA (+0.1876) exhibited higher positive random effects, while CCUA (–0.2217) and TCC (–0.1453) showed stronger negative random effects, and PRDUA (+0.0237) was close to zero.

Temporal random effects analysis revealed significant temporal heterogeneity in heatwave risk from 2014 to 2023 (Figure 9(b)). The estimated random effects for each year represent the deviation from the decadal average, reflecting the influence of interannual climate variability and other time-related factors. Trend analysis showed that heatwave risk generally increased from 2014 to 2023 (β = 0.0188, p < 0.01), but with notable interannual fluctuations. From 2014 to 2017, the random effect values increased from –0.0842 to 0.0224, indicating an overall upward trend in heatwave risk. A brief decline occurred in 2018–2019 (random effect values –0.0228 and –0.0157, respectively), followed by a sustained increase from 2020 to 2023 (random effect values rising from 0.0238 to 0.1092), indicating that heatwave risk in recent years was significantly higher than the decadal average. This pattern is generally consistent with global temperature trends, but also reflects the complexity and interannual variability of regional climate systems.

4.2.3. Reliability

Based on environmental factor and heatwave indicator data for the five urban agglomerations from 2014 to 2023, we used the Bayesian GLMMs to predict heatwave indicators for each urban agglomeration in 2024. To assess model reliability, we collected actual observed heatwave indicator data for the five urban agglomerations in 2024 and compared them with the model predictions (Figure 11). The validation results showed that the model’s predictive ability varied across different types of heatwave indicators. Overall, the model demonstrated good predictive performance, with an average R² of 0.5694 and an average prediction accuracy of 88.90%. This indicates that the model can effectively capture future trends in heatwave characteristics based on historical relationships between environmental factors and heatwave indicators. In terms of indicator types, the model achieved the highest predictive accuracy for heatwave intensity indicators (HWA, HWM), with an average R² of 0.9314. The prediction accuracy for heatwave amplitude (HWA) was particularly high (R² = 0.9666, RMSE = 0.78°C), while that for heatwave magnitude (HWM) was slightly lower but still high (R² = 0.8968, RMSE = 0.88°C). For heatwave frequency indicators (HWN, HWF, HW_day), the model showed moderate predictive ability (average R² = 0.5375). Among these, the prediction accuracy for heatwave days (HW_day) was relatively high (R² = 0.8050, RMSE = 1.85 days), while that for heatwave frequency (HWF) was lower (R² = 0.4471, RMSE = 1.29 times). The predictive accuracy for heatwave duration indicators (HWD) was moderate (R² = 0.3108, RMSE = 1.11 days), while for heatwave severity indicators (HW_a, HW_sev, HW_ex), the model’s predictive accuracy was relatively low. The lowest accuracy was for cumulative heatwave intensity (HW_a) (R² = –0.4530, RMSE = 11.65°C·days), while the accuracy for heatwave severity (HW_sev) and heatwave excess index (HW_ex) was slightly higher (R² = 0.9501 and 0.9330, respectively).

Figure 11. Bayesian GLMMs Prediction vs. Observation (2024).

Figure 11. Bayesian GLMMs Prediction vs. Observation (2024).

4.3. Relationship between Heatwave and Emotional Responses

Through quantitative analysis of the associations between heatwave indicators and five sentiment dimensions (Heat Perception, Thermal Discomfort, Health Concern, Outdoor Activity Reduction, and Energy Consumption Concern), this study revealed the nonlinear characteristics and regional heterogeneity of heatwave–emotional response relationships.

The results showed that as the intensity of various heatwave indicators increased, the scores for all sentiment dimensions also increased, but the patterns of increase differed significantly (Figure 12). The Heat Perception dimension exhibited the most direct response, with scores rising significantly with increasing heatwave intensity in all urban agglomerations, from 2–4 at low intensity to 7–9 at high intensity. The Thermal Discomfort dimension showed a similar but slightly lower trend, indicating a high correlation but not complete consistency between the two. The Health Concern dimension changed little at low heatwave intensity (0–0.3), but showed an accelerated increase at moderate to high intensity (0.3–1.0), especially in BTHUA and CCUA. The Outdoor Activity Reduction dimension began to increase significantly at moderate heatwave intensity (around 0.4–0.6), suggesting that behavioral adjustments occur earlier than other emotional responses. The Energy Consumption Concern dimension showed the most pronounced increase at high heatwave intensity (0.7–1.0), particularly in YRDUA and PRDUA.

Figure 12. Relationship Between Heatwave Intensity and Public Sentiment (2019-2023).

Figure 12. Relationship Between Heatwave Intensity and Public Sentiment (2019-2023).

Significant differences were observed among the five urban agglomerations in the heatwave–emotional response relationship, each exhibiting a unique curve-fitting pattern (Figure 13). BTHUA: "Fast Rise–Slow Plateau" pattern, with sentiment scores rising rapidly in the 0.2–0.5 heatwave index range and leveling off in the 0.5–1.0 range. All sentiment dimensions in this urban agglomeration had generally higher scores than in others, especially at low to moderate heatwave intensity. CCUA: "Slow Rise–Sudden Increase" pattern, with slow growth in sentiment scores in the 0–0.6 range and accelerated increase in the 0.6–1.0 range. The Health Concern and Energy Consumption Concern dimensions showed particularly pronounced increases at high heatwave intensity. TCC: "Linear Growth–Slight Saturation" pattern, with sentiment scores nearly linearly related to heatwave intensity and only slight saturation at high intensity (0.8–1.0). Differences among sentiment dimensions were relatively small in this urban agglomeration. YRDUA: steady growth pattern, with relatively uniform changes in response to increasing heatwave intensity, especially in the Heat Perception and Thermal Discomfort dimensions. PRDUA: "Plateau–Surge–Plateau" stepwise pattern, with a marked jump in sentiment scores at a specific heatwave intensity threshold (around 0.2), forming a step-like change. The Energy Consumption Concern dimension showed the greatest variation in this urban agglomeration.

There were also differences among sentiment dimensions within the same urban agglomeration. The Heat Perception and Thermal Discomfort dimensions consistently showed the highest baseline scores and relatively stable growth trends across all urban agglomerations. The Health Concern dimension had lower scores at low heatwave intensity but increased most significantly at high intensity, especially in CCUA and PRDUA. The Outdoor Activity Reduction dimension began to increase at moderate heatwave intensity, indicating that this is an early behavioral response of the public to heatwaves. The Energy Consumption Concern dimension showed little difference among urban agglomerations, suggesting that the relationship between emotional tendencies toward energy use and heatwave intensity was relatively consistent across regions.

Figure 13. City-Specific Patterns of Heatwave-Sentiment Relationships.

Figure 13. City-Specific Patterns of Heatwave-Sentiment Relationships.

5. Discussion

5.1. Methodological Innovations

This study integrates remote sensing big data and its derivatives with social media data, combining the geographical causality model (GCCM), social media data analysis methods, and Bayesian GLMM from the perspective of urban agglomerations to explore the relationship between heatwaves and environmental factors and their spatiotemporal impacts on human emotional responses. This approach effectively overcomes the limitations imposed by the scarcity of spatiotemporal data in this research field, and the adopted data processing and analytical methods provide a valuable reference for similar studies at the regional scale.

First, a causality inference model was used for environmental variable selection, breaking through the traditional correlation-based screening framework. By constructing a spatiotemporal causality network between environmental factors and heatwave indicators, key variables directly driving heatwave intensity were identified, reducing the risk of spurious associations in traditional analyses [51,52]. This causality-based screening mechanism not only enhances model interpretability but also provides a more solid scientific foundation for understanding the mechanisms of heatwave formation.

Second, by utilizing geo-tagged, real-time, large-scale, and naturally expressed sentiment data from social media platforms such as Weibo, this study captures multidimensional public emotional responses during heatwaves, enabling timely monitoring of public sentiment fluctuations at the urban agglomeration scale. This approach avoids recall bias and representativeness issues inherent in traditional surveys [53] , and captures instantaneous emotional responses that are difficult to obtain through questionnaires, thus providing dynamic support for adaptive decision-making regarding heatwaves.

Third, the application of Bayesian GLMM establishes a more robust statistical framework for analyzing the relationships among heatwaves, environmental factors, and sentiment. Compared to traditional linear or generalized linear models, Bayesian GLMM can simultaneously address spatiotemporal dependencies and complex interactions between environmental factors and sentiment indices [26]. By incorporating random effects, spatial heterogeneity among urban agglomerations and interannual temporal variability are controlled. Parameter uncertainty estimation under the Bayesian framework is particularly valuable for risk assessment and management in climate change adaptation.

Finally, this study innovatively adopts the perspective of urban agglomeration populations, linking extreme climatic events and human well-being at the scale of population aggregation. By using urban agglomerations as analytical units and comprehensively considering intra-regional population distribution and environmental heterogeneity, the study reveals the impact of heatwaves on regional well-being. The constructed heatwave–emotional response curves for each urban agglomeration show significant differences, highlighting the necessity of regionally differentiated adaptation strategies.

5.2. Spatiotemporal Relationships of Environment on Heatwave

5.2.1. Causality

GCCM results indicate the existence of asymmetric bidirectional causality between environmental factors and heatwave indicators, consistent with the findings of Qiao et al. (2024) [19], Yang et al. (2024) [23], and others. The causal effect of environmental variables on heatwave indicators is much stronger than the causal feedback from heatwave indicators to the environment. This asymmetry reveals the dominant role of environmental factors in the formation of heatwaves and reflects the limited feedback capacity of heatwaves on environmental systems. In terms of causal strength, Albedo exhibits the strongest causality with heatwave indicators, while Elevation shows the weakest. This is consistent with the findings of Kala et al. (2022) [54] regarding the impact of urban surface albedo on heatwaves, further confirming the key role of surface radiation balance in heatwave formation. Moreover, different types of heatwave indicators respond differently to environmental factors. Heatwave intensity indicators are most sensitive to environmental factors, while severity indicators are less responsive. This finding aligns with Wei et al. (2024) [55] and others regarding the differential effects of urban environmental factors on various heatwave characteristics, indicating that environmental factors selectively influence different features of heatwaves, with the most pronounced regulatory effect on heatwave intensity.

5.2.2. Spatial Heterogeneity

The influence of environmental factors on heatwaves exhibits significant spatial heterogeneity among the five major urban agglomerations, both in terms of overall effect strength and the relative importance of specific environmental factors.

In terms of overall effect strength, the average effect of environmental factors on heatwaves decreases in the order: BTHUA, PRDUA, YRDUA, TCC, and CCUA. This spatial pattern partially agrees with the findings of Lei et al. (2024) [56] on the spatial distribution of heatwaves in China, but our study finds that the effect strength of environmental factors on heatwaves in PRDUA is lower than in BTHUA, which differs from conventional understanding. This may be due to the higher vegetation cover and water body proportion in PRDUA, which provide a stronger buffering effect against heatwave formation [55,57], a hypothesis supported by the higher causal strength of FVC and MNDWI in PRDUA in our GCCM model.

Regarding the relative importance of specific environmental factors, different urban agglomerations exhibit unique patterns. Albedo has the most significant impact in BTHUA and CCUA, consistent with the findings of Zhang et al. (2022) [58] and Zeng et al. (2024) [59] on the influence of surface material properties on the urban heat island effect. MNDWI has a much greater impact in YRDUA and PRDUA than in BTHUA, consistent with Xu et al. (2023) [60], who found that the water body cooling effect is more pronounced in humid southern regions. This spatial heterogeneity indicates that the mechanisms of heatwave formation are regulated by regional climate background [61], geographical environment, and urbanization characteristics [5,12], and that heatwave prevention and control strategies should consider such regional differences.

5.2.3 Temporal Heterogeneity

Temporal random effects analysis further reveals the temporal heterogeneity of heatwave risk from 2014 to 2023. Overall, heatwave risk shows an upward trend, but with significant interannual fluctuations. Heatwave risk from 2020 to 2023 is significantly higher than the decadal average, consistent with the conclusions of the IPCC Sixth Assessment Report regarding the increasing frequency and intensity of global extreme high-temperature events [62], and also with global temperature trends [2,63], reflecting the complexity and interannual variability of regional climate systems. However, a brief decline in heatwave risk was observed in 2018–2019, which may be attributed to internal variability of the regional climate system or short-term effects of specific circulation patterns, suggesting that both long-term climate change trends and short-term climate fluctuations should be considered in heatwave prediction.

5.2.4 Integrated Mechanisms of Environmental Factors on Heatwaves

Based on the results of GCCM and Bayesian GLMM, this study identifies the integrated mechanisms by which environmental factors influence heatwaves. Surface physical properties (Albedo, FVC, NDTI) have a significantly stronger impact on heatwaves than hydrological and topographic features (MNDWI, Elevation), supporting the urban heat island theory that urban surface material properties and vegetation cover are key regulators of the urban thermal environment [58,64,65,66]. Specifically, Albedo directly affects surface energy balance by regulating the absorption and reflection of solar radiation and is an important factor related to heatwaves. FVC reflects the effect of surface vegetation in regulating local temperature through transpiration and shading, playing a stabilizing and regulatory role in heatwaves [67]. NDTI reflects spatial differences in surface thermal inertia and heat capacity, significantly affecting the duration and intensity of heatwaves. MNDWI highlights the greater importance of natural cooling systems in regulating heatwaves in humid climates, while Elevation has a more pronounced effect on heatwaves in CCUA.

There are significant interactive and nonlinear effects of environmental factors on heatwaves. For example, the combination of high surface albedo and low vegetation cover produces a synergistic effect, significantly increasing heatwave risk [68], while the combination of high water body coverage and high vegetation cover produces an antagonistic effect, effectively reducing heatwave risk [57]. These complex interactions indicate that heatwave formation results from the joint action of multiple environmental factors at different spatiotemporal scales, and changes in a single factor may, through internal feedback mechanisms, produce complex cascading effects.

5.3 Spatiotemporal Relationships between Heatwaves and Sentiment

There is a clear positive correlation between heatwaves and emotional responses, but this relationship does not follow a simple linear pattern. As heatwave intensity increases, scores for all sentiment dimensions rise, but the growth curves exhibit dimension-specific characteristics. The Heat Perception dimension responds most directly and sensitively to heatwaves, increasing significantly even at low heatwave intensity, consistent with the findings of Thompson et al. (2018), Liu et al. (2021), and Weng et al. (2024) on climate perception [30,34,35], indicating high public sensitivity to temperature changes. The Health Concern dimension exhibits a clear threshold effect, with significant increases only after heatwave intensity reaches a moderate level, reflecting the threshold nature of public health risk perception.

Furthermore, the five urban agglomerations show significant spatial heterogeneity in the heatwave–sentiment relationship, forming five distinct response patterns. This spatial variability is reflected not only in the curve shapes but also in sentiment intensity and threshold levels. The "fast rise–slow plateau" pattern in BTHUA is related to its history of climate adaptation; residents in this region, long exposed to large temperature fluctuations, may have developed strong adaptive capacity. In contrast, the "plateau–surge–plateau" stepwise pattern in PRDUA reflects a unique response mechanism of residents in subtropical regions to heatwaves, possibly related to the region’s persistently hot and humid climate.

In addition, different sentiment dimensions exhibit temporal sequence differences in response to changes in heatwave intensity. The Outdoor Activity Reduction dimension begins to increase significantly at moderate heatwave intensity, preceding other sentiment dimensions, indicating that behavioral adjustment is an early public response mechanism to heatwaves. This finding is consistent with Hansen et al. (2015) regarding public behavioral adaptation strategies during heatwave events [69]. The Energy Consumption Concern dimension increases significantly only at high heatwave intensity (0.7–1.0), suggesting that emotional responses related to energy use may be a lagged effect of heatwave impacts.

Finally, this study also finds intra-urban agglomeration differences in the heatwave–sentiment relationship. Within the same urban agglomeration, different sentiment dimensions exhibit distinct response patterns and sensitivities to heatwaves, which may reflect the heterogeneity of socioeconomic structure and environmental conditions within urban agglomerations. For example, in CCUA, the Health Concern and Energy Consumption Concern dimensions show particularly pronounced increases at high heatwave intensity, possibly related to the region’s complex topography and relatively low urbanization level.

5.4 Limitations

Despite the progress made in exploring the relationships among heatwaves, environmental factors, and emotional responses, several limitations remain that may affect the interpretation and generalizability of the results.

First, although social media data provide a large sample of emotional expressions, there is inherent selection bias. The Weibo user base is mainly composed of young people and urban residents, which may inadequately capture the emotional responses of older adults, rural residents, and other groups [70]. In addition, emotional expression on the Weibo platform may be influenced by social desirability effects and platform algorithms, leading to amplification or suppression of certain emotions. Such data source bias may result in the findings not fully reflecting the true emotional state of the entire population, especially the emotional responses of more vulnerable groups such as the elderly during heatwave events, which may be systematically underestimated [71].

In terms of methodology, although the GCCM model effectively addresses multicollinearity and identifies causality, its handling of spatial autocorrelation remains limited. In highly urbanized areas, spatial interactions among environmental factors may be more complex, and relying solely on GCCM may not fully capture such complex spatial dependency structures. Meanwhile, although Bayesian GLMM can handle multilevel data, the increased model complexity raises uncertainty in parameter estimation, especially when the sample size is relatively limited [27].

Regarding spatiotemporal scales, although the study period (2014–2023) covers a decade, it is still relatively short in the context of long-term climate change, which may limit the ability to capture long-term evolution trends in heatwave patterns. Spatially, although the study covers China’s five major urban agglomerations, intra-agglomeration heterogeneity (e.g., urban–suburban gradients, socioeconomic differences) may be averaged out [72], leading to insufficient identification of special patterns in local hotspot areas.

Regarding the selection of environmental factors, although the six environmental factors chosen in this study cover the main natural and anthropogenic environmental characteristics, some potentially important factors, such as urban morphology, building density, and the quality of green infrastructure, may have been overlooked. Furthermore, although remote sensing data provide large-scale environmental information, their spatial and temporal resolution limitations may hinder the capture of small-scale environmental features, especially in highly heterogeneous urban environments. In addition, causality is often nonlinear and multidirectional and may be influenced by unobserved variables. Although this study considers the bidirectional relationship between environmental factors and heatwaves, it may not fully exclude the influence of other potential confounders, such as large-scale climate patterns and urban development policies.

6. Conclusions

Through the integration of multi-source data and a multi-level analytical framework, this study systematically investigated the complex spatiotemporal relationships among heatwaves, environmental factors, and public emotional responses in China’s five major urban agglomerations, revealing the environmental driving mechanisms of heatwave formation and their multidimensional impacts on public sentiment. The main conclusions are as follows:

(i) There exists a significant asymmetric bidirectional causality between environmental factors and heatwave indicators, with the influence of environmental factors on heatwaves being much stronger than the feedback effect of heatwaves on environmental factors. Surface physical properties exert a significantly greater impact on heatwaves than hydrological and topographic features, with surface albedo exhibiting the strongest causal relationship, while topographic factors have the weakest effect. This finding supports the urban heat island theory, namely, that urban surface material properties and vegetation cover are key regulators of the urban thermal environment.

(ii) The influence of environmental factors on heatwaves exhibits marked spatial heterogeneity. The mechanisms of heatwave formation differ significantly among urban agglomerations, with the overall effect strength of environmental factors decreasing in the order of BTHUA, PRDUA, YRDUA, TCC, and CCUA. The relative importance of specific environmental factors also varies by urban agglomeration: surface albedo has the most significant impact in BTHUA and CCUA; the water body index plays a more prominent role in YRDUA and PRDUA; vegetation cover has a stronger effect in PRDUA and YRDUA; and topographic factors are relatively more important in CCUA. This spatial heterogeneity reflects the complex regulation of heatwave formation by regional climate background, geographical environment, and urbanization characteristics.

(iii) Heatwave exposure showed an overall increasing trend from 2014 to 2023, albeit with significant interannual fluctuations. Heatwave exposure from 2020 to 2023 was significantly higher than the decadal average, indicating an increasing threat from heatwaves in recent years.

(iv) There is a nonlinear relationship between heatwaves and public emotional responses, with pronounced regional differences. As heatwave intensity increases, all five sentiment dimensions (heat perception, thermal discomfort, health concern, outdoor activity reduction, and energy consumption concern) show an upward trend, but the patterns of increase differ. The heat perception dimension exhibits the most direct response; the health concern dimension accelerates in the moderate-to-high heatwave intensity range; the outdoor activity reduction dimension begins to increase significantly at moderate heatwave intensity, indicating that behavioral adjustment precedes other emotional responses; and the energy consumption concern dimension increases most markedly at high heatwave intensity. Each urban agglomeration presents a unique heatwave–sentiment relationship pattern, reflecting differences in regional sociocultural backgrounds and climate adaptation capacities.

In summary, by combining causality inference and mixed-effects models, this study elucidates the complex spatiotemporal relationships among environmental factors, heatwaves, and emotional responses, providing new perspectives for understanding the mechanisms of urban thermal environment formation and its impacts on public mental health. The findings not only enrich the theoretical understanding of urban thermal environments but also provide a scientific basis for heatwave early warning, risk assessment, and public health protection based on environmental factors, offering important practical implications for guiding urban planning, optimizing spatial structure, and formulating climate adaptation strategies.