Climate and forest form divergences amplify spatiotemporal heterogeneity and environmental inequity in China’s forest cooling services

Driven by the coupling of rapid urbanization and climate change, Urban Heat Island (UHI) risks have intensified across China. However, given the nation’s vast geographic expanse, the cooling services of Urban Forests (UF) exhibit significant spatial heterogeneity. Neglecting human-centric perspectives and environmental justice may exacerbate gaps in cooling services and lead to inequitable distribution of environmental benefits. This study investigates 368 Chinese cities using the GTNNWR model to quantify three core indicators: Cooling Efficiency (CE), Cooling Capacity (CC), and Cooling Benefit (CB). We systematically explore the spatiotemporal heterogeneity, driving mechanisms, and equity evolution of cooling services. The results indicate that: (1) UF demonstrate a pronounced cooling effect with an average CE of 3.61, implying that for every one-standard-unit increase in the UFCI, the LST decreases by 3.61°C; however, the average CC (0.55) and CB (0.33) remain relatively low (CC and CB are both dimensionless). (2) Driving mechanisms exhibit multidimensional differentiation: natural factors exert negative driving effects by regulating thermal baselines and vegetation growth, while socio-economic factors exert positive driving effects through anthropogenic heat emissions and land-use changes. (3) The Gini coefficient for cooling service equity is 0.844, categorized as highly unequal, though this inequality has gradually moderated in recent years. To enhance urban thermal resilience and ensure a balanced supply of cooling services, this study proposes sustainable development strategies from the dual perspective of climate zones-forest forms, providing scientific support for differentiated urban thermal environment governance.

Driven by the dual pressures of global warming and rapid urbanization in China, the continuous expansion of high-density built-up areas has exacerbated the Urban Heat Island (UHI) effect, posing a severe threat to the quality of human settlements and public health (Tuholske et al., 2021). As a nature-based solution, urban forest (UF) plays an irreplaceable role in mitigating heat stress and regulating microclimates through multiple ecological processes like canopy shading, transpiration and air circulation regulation (Li H. et al., 2024). Consequently, its cooling services have become a core focus of urban ecological planning and management.

The UF refers to tree-based urban ecosystems located within and around cities, encompassing all woodlands, trees groups, and individual trees (Konijnendijk et al., 2005). Previous studies show that UF cooling services exhibit significant spatiotemporal heterogeneity rather than a uniform distribution (Yu et al., 2025). This heterogeneity is co-constrained by macro-climates and micro-forest forms. Specifically, modulated by climatic conditions like temperature and precipitation, the interaction between cooling effect and forest coverage does not follow a simplistic linear decline, but instead demonstrates complex threshold dynamics and fluctuating patterns (Cheng et al., 2022; Zhang J. et al., 2026). On the other hand, tree morphology determines its capacity to intercept solar radiation, thereby defining its shading ability (Still et al., 2022). Research by Moser-Reischl et al. indicated that the cooling intensity of different tree species or canopy structures can vary by up to four times (Moser-Reischl et al., 2019). Ma and Wang noted that broadleaf forests and mixed forests provide better cooling effects than needleleaf forests (Ma and Wang, 2024). However, existing research has predominantly focused on local effects within a single city or specific climate zones, rarely considering the dual interactive mechanisms of macro-climate zones and micro-forest forms at a national scale, limiting the precise formulation of cross-regional climate adaptation planning and the effective allocation of woodland resources.

Building on the elucidation of spatiotemporal heterogeneity of cooling services, exploring distribution patterns of this ecological service within urban spaces is the core essence of achieving equality in urban environmental wellbeing (Li Y. et al., 2024). Traditional environmental justice research mostly relies on nighttime static registered population to depict social attributes. However, urban heat exposure and heat stress risks typically peak during the daytime active period of populations (Li et al., 2026). To bridge the supply-demand gap in cooling services caused by the current emphasis on UF construction over services, we focuses the connotation of environmental justice on universal spatial entitlement and spatial supply-demand matching. We innovatively construct a three-dimensional evaluation index system for cooling: CE is employed to measure the reduction in LST resulting from a unit increase in UF, aiming to characterize the inherent physical cooling potential of woodlands (Marando et al., 2022); CC refers to magnitude of absolute regional surface temperature reduction that existing UF can provide relative to unforested built-up areas, reflecting ecological baseline capacity of a city to resist thermal stress at a macro scale; CB quantifies actual cooling wellbeing accrued by populations under realistic heat exposure scenarios, an indicator that deeply depicts the spatial coupling characteristics between UF ecological supply and dynamic population (POP) distribution (Li Y. et al., 2024). By establishing a heat mitigation framework of ecological supply - spatial transfer - population benefit, these three indicators objectively reveal the spatial mismatch and equity characteristics of cooling services in terms of quality, quantity, and secure accessibility (Wolch et al., 2014; Xie et al., 2026).

While quantifying spatial mismatches and environmental equity provides a necessary baseline, delving into the underlying driving mechanisms of this spatiotemporal heterogeneity is fundamental to the development of targeted climate adaptation strategies (Wang et al., 2025). At present, explainable machine learning and spatial econometric methods have been widely applied to explore the driving factors of the urban thermal environment and identify critical thresholds. For instance, an analytical framework integrating the CatBoost model with SHAP has successfully decoupled the non-linear effects of urban morphology on the thermal environment (Yang et al., 2026a). However, when dealing with large-scale urban spaces, cooling services exhibit inherent characteristics such as profound spatial non-stationarity and complex scale-dependent interactions. Traditional geographically weighted regression and machine learning methods struggle to adaptively account for these spatial constraints (Zhang X. et al., 2026). Consequently, it is necessary to introduce explainable deep learning-based approaches to capture these localized dynamic features (Yang et al., 2026b). To this end, we introduce the Geographically and Temporally Neural Network Weighted Regression (GTNNWR) model, which incorporates a spatiotemporal weighting framework. By leveraging the non-linear mapping capabilities of neural networks, this model adaptively generates spatiotemporal weight matrices. It exhibits enhanced robustness and scalability when processing spatiotemporal sequential data, overcoming the “black-box” nature of traditional machine learning models and achieving both high fitting accuracy and mechanistic interpretability (Que et al., 2020).

In summary, taking Chinese cities as study areas and adopting a dual perspective of climate zones and forest forms, we utilizes the GTNNWR model to address following three core questions: What are the spatiotemporal heterogeneity characteristics of cooling services in Chinese cities?

How is this pattern co-regulated by climate zones and forest forms?

How are the resulting environmental equity issues of urban cooling services manifested and evaluated?

These findings provide theoretical support for tailored urban climate adaptation planning at a global scale, offering scientific evidence for optimizing green space layouts in high-density cities and enhance overall environmental health wellbeing.

Study area

We focus on Chinese cities using the Nature Cities dataset (Cattaneo et al., 2024). We defined boundaries for 368 sample cities based on having “a population exceeding 20,000 within a 1-h commuting circle,” to reflect actual human activity patterns more accurately than traditional administrative boundaries (Figure 1).

Given China’s significant geographical diversity, we classified the cities into four climate zones based on Aridity Index (AI) and mapped them into four major spatial units using key boundaries like the Qinling-Huaihe Line (Peng, 2023). Regarding the ecological baseline, the average UF coverage in the study area is 27.13%, exhibiting a distinct North-South disparity. Based on forest cover (FC) characteristics, forest forms were further subdivided into six categories (Cui et al., 2024), providing a solid foundation to analyze cooling service patterns across various climates and forest forms.

Datasets and processing

We utilized LST to characterize UHI effect, with data derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) annual average daytime LST product (MOD11A1). This product utilizes dual thermal infrared channels to achieve atmospheric correction and quantitative inversion of surface temperature. It is an authoritative and universal data source for large-scale and long-term studies on the surface thermal environment. To mitigate interannual temperature fluctuations, we averaged pixel values to obtain summer LST data (June to August) for each year (Zhang et al., 2025).

Additionally, two categories of auxiliary data—natural and anthropogenic—were simultaneously collected (Table 1). Natural environmental variables included precipitation (PREC), wind speed (WS), relative humidity (RHU), Normalized Difference Water Index (NDWI), digital elevation model (DEM), and topographic slope (Slope), serveing as proxy variables for characterizing urban background climate and ecological conditions (Fu et al., 2025). Socio-economic variables encompassed nighttime light (NTL), population distribution (POP), and Normalized Difference Built-up Index (NDBI), utilized to quantify human activities intensity and urbanization impacts (Shi et al., 2025).

All datasets underwent systematic preprocessing and were uniformly organized into six cross-sectional time periods: 2001, 2005, 2010, 2015, 2020, and 2022. Data were projected to the WGS84 UTM Zone 48N coordinate system with a 1 km × 1 km spatial resolution. Subsequently, prior to formally proceeding with spatial analyses and machine learning modeling, we uniformly applied Z-score standardization to all explanatory variables. This step was taken to eliminate inherent dimensional discrepancies and magnitude differences among the driving factors, and to ensure the stable convergence of the subsequent core models, thereby transforming them into dimensionless values with horizontal comparability (Pramoedyo et al., 2026).

Meanwhile, Spearman correlation analysis (p < 0.001) and multicollinearity diagnostics (VIF <7.5) confirmed significant correlated without severe redundancy among candidate variables (Supplementary Table S1), satisfying regression modeling prerequisites (Supplementary Figure S1).

Research methods

Research framework

This study uses a systematic methodological framework (Figure 2), combining remote sensing data with advanced statistical analysis to comprehensively analyze the coupling mechanism between UF and LST. The framework consists of three key steps. The first step is data integration and indicator construction, including remote sensing image data collection, preprocessing, and the use of Principal Component Analysis (PCA) to extract the Urban Forest Comprehensive Index (UFCI), standardizing the response, feature, and control variables. The second step involves spatiotemporal pattern identification, using landscape pattern indices and the mean−standard deviation method to quantify the spatiotemporal evolution of UFCI and the temperature grading characteristics of LST, thereby clarifying their distribution differences at multiple spatial and temporal scales. The third step is the quantification and evaluation of cooling services. The GTNNWR model is selected for model performance evaluation, utilizing its ability to handle spatiotemporal non-stationarity to derive and quantify Cooling Efficiency (CE), Cooling Capacity (CC), and Cooling Benefit (CB), revealing the inequality of cooling services. This framework provides theoretical support and decision-making reference for mitigating urban heat island effects and optimizing green space layout.

Construction of the urban forest comprehensive index (UFCI)

The Normalized Difference Vegetation Index (NDVI), maximum Leaf Area Index (LAI), and Fractional Vegetation Cover (FVC) characterize vegetation features through greenness, three-dimensional canopy structure, and planar coverage, respectively, offering highly complementary information (Supplementary Table S2). To overcome single-indicator limitations, we utilized orthogonal decomposition and adaptive weighting mechanism of Principal Component Analysis to integrate these three variables into the Urban Forest Comprehensive Index (UFCI) (Greenacre et al., 2022).

Results indicate that cumulative variance contribution rate of PC1 and PC2 reached 88.53%, with PC1 yielding an eigenvalue of 2.000, demonstrating exceptional explanatory power. Consequently, PC1 was selected as UFCI to comprehensively represent three-dimensional total vegetation volume in quantity, thickness, and quality (Table 2).

Landscape pattern indices

Research has shown that landscape pattern indices are widely used to describe and quantify landscape structure and composition. We employed Fragstats 4.2 to calculate six indices: percentage of landscape type area (PLAND), Largest Patch Index (LPI), Patch Density (PD), Patch Cohesion Index (Cohesion), Mean Fractal Dimension (FRAC_MN), and Split Index (SPLIT). These indices comprehensively capture spatial configuration, patch dominance, density, connectivity, morphological complexity, and fragmentation. Given the focus on forest categories, we prioritized class-level indices to systematically describe the overall spatial distribution of forests.

Mean-standard deviation method

To mitigate systematic errors induced by long-term interannual meteorological fluctuations and prevent the direct comparison of absolute LST values from obscuring the authentic spatial patterns of the UHI, we employed the mean−standard deviation method to categorize LST into relative thermal levels (Shi et al., 2025). The classification thresholds of this approach are derived autonomously, based entirely on inherent statistical characteristics of images. This effectively circumvents biases associated with subjective parameterization, ensuring objectivity and rigor of the classification results, and establishing a robust foundation for identifying urban thermal trends (Supplementary Table S3).

Spatial autocorrelation analysis

Spatial autocorrelation characterizes the numerical dependency of geographic variables in spatial distribution, meaning the statistical association between the observation value of one spatial unit and the attribute values of its neighboring units. We used binary Moran’s I to examine the spatial autocorrelation between LST and UFCI, quantifying their spatial clustering patterns and the strength of their association, statistically testing whether spatial relationships and their statistical significance existed (Luo et al., 2023).

Regression analysis

GTNNWR model

In the quantitative research of the urban heat island effect and ecosystem services, classical machine learning models such as XGBoost and Random Forest (RF) have been widely applied due to their powerful non-linear fitting capabilities and stable generalization performance. However, the performance of these algorithms often varies depending on the dataset and problem background. Furthermore, their capture of spatial effects is an implicit fitting process, which to some extent limits the reliability of model interpretation and the directness of mechanism analysis (Pramoedyo et al., 2026).

The core objective of this study is not merely to maximize global prediction accuracy, but rather to reveal the spatiotemporal heterogeneity of urban forest cooling services and their driving mechanisms, thereby providing a mechanistic basis for the differentiated layout of urban green spaces. Therefore, we selected classical spatial econometric models with explicit spatiotemporal weight structures to construct a progressive comparison framework of “OLS-GWR-GTWR-GTNNWR.”. This framework verifies the performance gains achieved by spatiotemporal structural improvements and directly outputs local regression coefficients for each spatial unit to support the quantification of spatial differentiation.

Specifically, OLS is a global linear model that assumes spatiotemporal constancy; GWR introduces a spatial kernel function to characterize spatial non-stationarity but ignores temporal dynamics; GTWR further constructs a spatiotemporal joint kernel function but remains constrained by traditional predefined kernel forms. As the core model of this paper, GTNNWR breaks through this limitation by deeply integrating a two-layer cascaded neural network—consisting of a Spatiotemporal Proximity Neural Network (STPNN) and a Spatiotemporal Weight Neural Network (STWNN)—with a global linear regression mechanism (Wu et al., 2021) (Figure 3). This integration allows for an accurate characterization of the complex spatiotemporal dynamic relationships underlying the cooling services of the UFCI. The core mathematical expressions and operational processes are detailed as follows:

First, the STPNN maps the independent spatial distance (d s) and temporal distance (d t) into a unified generalized spatiotemporal distance, effectively resolving the challenges of dimensional inconsistency across heterogeneous dimensions and non-linear interactions. The formula is:

Subsequently, utilizing the aforementioned generalized spatiotemporal distance, the STWNN adaptively learns and generates the spatiotemporal weight matrix (W i):

Finally, the generated spatiotemporal weights modulate the global average coefficients to obtain the varying coefficient β i (u i, v i, t i) for the i th spatial unit, which reflects localized features. This establishes the complete geographically and temporally weighted regression equation: where (u i, v i, t i) represents the spatiotemporal coordinates of the i th sample, x ik is the k th explanatory variable, and ε i is the residual term. To enhance the model’s generalization ability and convergence efficiency, this architecture deeply integrates optimization techniques including Dropout, PReLU, and Batch Normalization.

For large-scale long-series dataset, a 0.1% random sampling strategy was adopted to balance computational efficiency and model reliability (Miniandi et al., 2025). The data were partitioned into training, testing and validation datasets at a ratio of 75%:15%:10%. Through multi-dimensional comparisons with OLS, GWR, and GTWR models across indicators including the coefficient of determination (R 2), root mean square error (RMSE), and mean absolute error (MAE), the optimal performance of this model in characterizing complex spatiotemporal heterogeneity was screened and validated.

Cooling indicators

Cooling Efficiency (CE)

CE is defined as the marginal response rate of LST to UFCI (Hu et al., 2026), mathematically characterized as the absolute value of the slope derived from their spatial fitting (Marando et al., 2022). A higher CE indicates a greater magnitude of LST reduction achieved per unit increase in the forest index. To eliminate interference of confounding factors such as water bodies and meteorological conditions, we constructed a regression model encompassing multi-dimensional covariates. The model expression is specified as follows: LST ∼ UFCI + PREC + RHU + WS + DEM + SLOPE + NDWI + NTL + POP + NDBI. Leveraging the GTNNWR model to capture spatiotemporal non-stationarity, CE is quantified as the absolute UFCI regression coefficient, reflecting cooling efficacy variations across climate zones and forest forms. 2. Cooling Capacity (CC)

CC is defined as total potential of woodland ecological patches to mitigate environmental thermal loads, with its numerical value represented as the product of physical intensity (CE) and spatial extent (UFCI) (Equation 1). Its physical essence lies in the total amount of energy regulation collectively achieved by all UF patches within a region, serving as the physical characterization of the supply side of ecosystem services. CC comprehensively reflects dual influence of quality and quantity; specifically, higher CE and larger UFCI proportions jointly determine a more robust ecological cooling supply capacity, leading to enhanced stability in the city’s overall regulation of thermal environment.

CC represents the average cooling capacity of UF nationwide. CC i is cooling capacity of the i th grid unit, UFCI i is UFCI for that grid, UFCI min is the minimum UFCI across China for that year, and CE i is the cooling efficiency of the i th grid unit, with n being the total number of grid units nationwide. 3. Cooling Benefit (CB)

CB is a human-centric indicator grounded in concept of supply-demand matching. Substantively, it represents the spatial transmission of UF cooling services, which can be used to quantify cooling benefits actually accrued by urban residents from the perspective of social wellbeing (Equations 2 and 3). Beyond merely reflecting service coverage, this indicator sensitively identifies spatial mismatches between high-density populations and high-efficiency cooling sources, thus providing a critical criterion for assessing urban thermal environmental equity and the efficacy of climate-adaptive planning.

CB i represents cooling benefit for the i th grid unit, POP i is the population of that grid unit, POP ¯ represents average population size of nationwide, and CB is the average cooling benefit of urban forests nationwide.

Gini coefficient

Building upon the quantification of actual CB obtained by residents, we employed Gini coefficient to measure spatial distribution equilibrium of CB within each city (Equation 4), aiming to accurately evaluate social distributional equity of cooling ecological services (Li Y. et al., 2024). Concurrently, referring to the United Nations’ Gini coefficient classification benchmarks, the equity of urban cooling services was classified into five tiers.

Where n is the number of CB data points, S is the total sum of all CB values, k is the k th data point in the sorted list, and x (k) is the k th cooling benefit value in ascending order.

Spatiotemporal distribution patterns of UF and UHI

Morphological evolution characteristics of UFCI

From 2001 to 2022, UF in China underwent a transformation from quantitative to qualitative changes, which can be divided into three periods: extensive expansion, ecological restoration, and high-quality development (Table 3).

From 2001 to 2005, the PLAND showed a slight decrease, indicating that the expansion of built-up land outpaced the growth of forest land. Simultaneously, the LPI dropped sharply from 6.909 to 2.780, while the SPLIT surged to a historic peak (410) and Cohesion dropped to its lowest (98.17). These abrupt changes indicated that large forest patches within cities were fragmented, and landscape integrity and connectivity were severely damaged during the extensive expansion phase of urbanization.

From 2005 to 2010, PLAND and LPI showed signs of recovery, indicating that the trend of patch destruction was curbed. The PD increased sharply, signaling a significant growth of forest patches within the city; however, the FRAC_MN remained extremely stable, suggesting that the new patches were mainly artificial and had smooth edges. Overall, this period saw the beginning of green space recovery, but the model was still predominantly artificial green-insertion with small and regular patches, marking the ecological restoration phase.

From 2010 to 2022, PD stabilized at a high level of 0.015, while PLAND and LPI continued to rise, indicating that ecological forests no longer simply pursued an increase in patch numbers. The continued decline of SPLIT confirmed that forest landscapes were gradually recombining. This marked the shift of UF construction from simple quantity accumulation to high-quality development phase.

Spatiotemporal intensity distribution characteristics of LST

Temporally, the distribution of LST exhibited significant stability. Although there were minor interannual fluctuations, the medium-temperature zone consistently dominated (Figure 4g).

Spatially, the southern regions are primarily covered by blue areas representing high and sub-high temperature zones, whereas the northern regions are predominantly covered by red areas indicating low and sub-low temperature zones (Figures 4a–f). This distribution not only reflects the north-south gradient difference in LST intensity but also shows only slight dynamic changes over time, further confirming the long-term stability of China’s urban heat environment macro-pattern. It is noteworthy that in some areas, there was a grade jump: the northern semi-arid zone evolved from a medium-temperature zone to a low-temperature zone, while the central and southern humid and semi-humid zones evolved from a medium-temperature zone to a high-temperature zone.

Spatial correlation characteristics between UF and UHI

Results indicate that Moran’s I index from 2001 to 2022 ranged from −0.506 to −0.340, indicating a significant negative spatial correlation between UFCI and LST (Figure 5). With all p-values at 0.001 and extremely high absolute Z-scores, this spatial relationship demonstrates 99% statistical confidence. However, the declining absolute values of Moran’s I suggest a gradually weakening negative spatial correlation.

The Local Indicators of Spatial Association cluster map reveals significant North-South spatial differentiation characteristics. The northern region is dominated by High-High (HH) agglomerations, represented by red areas—where a high UFCI is accompanied by a high LST. Conversely, the southern region is primarily characterized by Low-Low (LL) agglomerations, denoted by blue areas—where a low UFCI is accompanied by a low LST. Meanwhile, High-Low (HL) and Low-High (LH) clusters are sporadically interspersed between them. This local spatial agglomeration phenomenon confirms the existence of spatial dependence among the variables, which is highly consistent with their univariate clustering effects. Consequently, it emphasizes the necessity of incorporating other key factors when investigating the underlying driving mechanisms, providing a solid basis for subsequent regression modeling.

Spatiotemporal heterogeneity of cooling services

Model evaluation

Results showed that GTNNWR model had an R2 of 0.931 in the training dataset, significantly higher than other models. The model also achieved a RMSE of 0.261 and a MAE of only 0.191, indicating a small bias between predicted and actual values and a strong explanatory power. The difference between the training, validation, and testing dataset was less than 0.1, indicating no overfitting (Table 4).

Additionally, model’s predicted values closely aligned with the real values on the 1:1 line, demonstrating excellent robustness and the ability to capture complex nonlinear spatiotemporal interactions (Figure 6).

Cooling indicators

Cooling efficiency

The average CE in Chinese cities was 3.61, indicating a 3.61 °C LST decrease per standardized unit increase in UFCI. However, the range of CE values varied widely, with the maximum (31.21) being 322 times higher than the minimum (0.10). Temporally, reflecting early fluctuations followed by a later decline (Figure 7g). Spatially, CE distribution showed a significant west-high, east-low, north-high, south-low pattern (Figures 7a–f). The long-term evolution trend further confirmed the temporal and spatial trends, showing that 58.37% of cities had a downward trend, mostly concentrated in the densely populated cities in the east and south. Conversely, 38.06% of the cities show no significant trend. While these areas are broadly dispersed nationwide, they occupy an absolute dominant position in the arid and semi-arid regions (Figure 7j), where UF are inherently scarce. (Figure 7j).

In climate zones, CE increases strictly exponentially with drought severity. CE in arid (18.117) and semi-arid regions (5.970) exceeded the annual average. The CE in arid regions is 22 times that of the humid regions, indicating positive regulatory effect of climate drought severity on CE (Figure 7h).

Regarding forest forms, deciduous forests exhibit higher CE than evergreen forests. Specifically, CE of deciduous coniferous forests is 5.96, 2.78 for deciduous broad-leaved, 0.74 for evergreen coniferous, and 0.69 for evergreen broad-leaved forests. This indicates that coniferous forests generally outperform broad-leaved forests in CE (Figure 7i).

Cooling capacity

The average CC in Chinese cities was 0.55, ranging from 0.005 to 14.15, indicating extreme imbalance in cooling capacity of UFCI. The trend was consistent with that of CE. Temporally, CC significantly declined after 2015, and remained below the overall average from 2020 to 2022 (Figure 8g). Spatially, CC gradually increased from east to west, with 67.21% of urban areas having a CC below the national average (Figures 8a–f). The long-term evolution trend showed that 41.16% of cities had no significant trend, while 55.01% showed varying degrees of decline, with a local concentration in certain areas (Figure 8j). The CC, being the product of CE and area, showed that despite the increase in PLAND, CC still mainly declined. This suggests that newly added forest areas did not meet the increasing demand for cooling services caused by deteriorating urban heat environments.

In climate zones, the average CC increases stepwise with drought severity. The maximum CE in arid zone (14.15) is 3.7 times that of the maximum in humid zone (3.82), and the minimum (0.04) is 7.5 times that of the minimum in humid zone (0.0053). CC intensity levels across climate zone are consistent with CE (Figure 8h).

Regarding forest forms, the deciduous coniferous forest remains the highest (1.36), nearly double the second-ranked deciduous broad-leaved forest (0.77). However, unlike CE, evergreen coniferous forest is lower than the evergreen broad-leaved forest (Figure 8i).

Cooling benefit

The average CB in Chinese cities was 0.33, ranging from 0 to 299.59, showing significant spatial heterogeneity. The annual average CB fluctuated rising initially and then falling, peaking in 2015 (0.36) and dropping to its lowest point in 2020 (0.31) (Figure 9g). Spatially, the overall CB was dominated by low values, with high-value areas scattered, and 84.12% of the regions were below national average, showing significant differences between regions (Figures 9a–f). The long-term evolution trend showed that the region with no significant downward trend was widely distributed in the densely populated eastern and central cities, while the regions with significant decline were concentrated in North China and parts of Southwest China. The proportion of cities with a significant downward trend (5.67%) was higher than the other two indicators, indicating that some samples showed obvious degradation in CB, while the proportion of cities with an upward trend (4.95%) was the highest among the three indicators, indicating that the increase in CB in some samples was relatively more prominent (Figure 9j).

In climate zones, similar to CE and CC, CB increases significantly along a gradient of drought severity. The CB in arid zone are relatively high (0.76), but highly dispersed CB values, whereas humid zones show lower, stable values (0.20), reflecting climate’s dominant role in shaping CB spatial patterns (Figure 9h).

Regarding forest forms, the mean values of CB are relatively small, and there is a significant pattern reversal compared to CE and CC. Deciduous broad-leaved forests perform the best (0.05), mixed forests rank second (0.04), while the deciduous coniferous forest that performed the best in CE and CC becomes the worst in CB (Figure 9i).

Spatial mismatch mechanisms of cooling services

The evolutionary patterns of cooling service processes of UFCI are deeply intertwined with the phased transitions in the development of Chinese UF, the hydrothermal differentiation characteristics across climate zones, the physiological and functional differences among forest forms, and the spatial patterns of POP. This reveals the underlying interactive logic of natural cooling capacity–ecological space supply–human demand matching.

Global dimension: mismatch between urban expansion and ecological governance

The three stages of extensive expansion, ecological restoration, and high-quality development experienced by UF from 2001 to 2022 directly dominated the coupled evolution trends of CE, CC, and CB (Yue et al., 2026). The three indicators showed a synergistic response with peak misalignment characteristic, with the misalignment becoming more pronounced as the stages evolved.

In 2001-2005 period, during the extensive expansion phase, rapid industrialization and urbanization led to fragmentation of forest land, while NDBI increased, weakening the cooling performance of green spaces (Yan et al., 2025). Although CE remained at a high level, fragmented forests led to CC being at a moderate-low level. Meanwhile, the misalignment between forest land edge and urban population distribution caused CB to exhibit an extreme imbalance with high values in a few regions and low values in most regions, confirming the green space expansion lagging behind urbanization argument (Yang et al., 2025), reflecting the disconnect between ecological benefits and social welfare under the extensive development model.

In 2005-2010 period, the ecological restoration phase, forest patch numbers surged due to sand prevention and control strategies such as the Three North Shelterbelt project (Li J. et al., 2024), and the spatial extent of UFCI expanded significantly. The marginal cooling effect of forests stabilized, and CE remained above the multi-year average. However, because artificial forests have lower shading and transpiration efficiency than natural forests (Beele et al., 2024), the growth of CC did not fully match the expansion of UFCI. The stable characteristic of FRAC_MN confirmed the common issue of edge rule and monotonous form in artificial forest patches (Cheng et al., 2024).

In 2010-2022 period, the high-quality development phase, the ecological priority concept promoted the improvement of UF quality. However, the deterioration speed of the urban heat environment outpaced the enhancement of forest cooling capacity, leading to a decrease in CE. This phenomenon aligns with Leng et al.‘s findings, which concluded that when UHI intensity exceeds a critical threshold, expanding tree canopy coverage does not lead to a synchronized increase in cooling effects (Leng et al., 2024). CB declined after 2015. This phenomenon underscores the importance of the vegetation supply-demand matching principle proposed by Li (Li and Sun, 2023), which states that the mere increase of green coverage cannot solve the fairness issue of cooling services (Leng et al., 2024; Li Y. et al., 2024).

Climate zone dimension: mismatch between natural background and demand distribution

Spatially, CE and CC gradually increase from east to west. CE in western arid and semi-arid regions are extremely high, essentially amplifying the relative cooling efficiency due to extreme background conditions. Regarding relative cooling efficiency, arid regions are dominated by bare soil and deserts. Forests, as scarce cooling sources, create a strong thermal contrast with surrounding environment, exhibiting an extreme oasis effect (Li and Sun, 2023). A study across eight U.S. cities along the aridity gradient provides evidence that tree cooling efficiency increases with aridity, and cooling effect is enhanced during heatwaves (Hao et al., 2025). Regarding absolute cooling capacity, dry air in arid regions severely limits tree transpiration cooling due to water constraints, and actual cooling per unit leaf area is often lower than in humid regions. Although CE is lower in humid regions, the abundant water supply allows for higher tree transpiration efficiency, making the absolute cooling capacity per unit of green space stronger. This pattern aligns with Li et al., confirming that tree cooling effectiveness is significantly regulated by background climate (Wang et al., 2025).

The spatial pattern of CB exhibits a significant reversal characteristic compared to those of CE and CC, with high-value areas predominantly concentrated in the East. Although CE and CC within eastern humid zone are relatively low, UF is extensively distributed across areas with high POP density, such as urban built-up areas and agricultural shelterbelt networks. Consequently, these ecological cooling sources are effectively translated into tangible cooling benefits for residents’ livelihoods, yielding a higher CB.

Forest form dimension: mismatch between physiological functions and spatial distribution

Forest forms further regulate the spatial differentiation of CE, CC, and CB, sometimes reversing indicator patterns within the same region. Regarding physical supply, Deciduous Needleleaf Forest (DNF) demonstrates the highest cooling potential, with its CE and CC being 8.5 and 4.6 times those of evergreen broadleaf forest, respectively. This physical advantage is primarily attributed to its specialized canopy structure; specifically, closely arranged needles efficiently intercept short-wave radiation (Speak et al., 2020), which, coupled with the lush LAI of deciduous forests in summer, significantly amplifies the canopy shading effect. This conclusion aligns with the field observations of Rahman et al., who noted that the shading rate of needleleaf forests is generally 20%–30% higher than that of broadleaf forests (Rahman et al., 2024).

However, the spatial pattern of CB exhibits a starkly opposite trend. While DNF possesses the strongest physical CC, its CB value ranks last, whereas deciduous broadleaf forest rank first. The root cause of this high supply, low efficiency paradox lies in the severe spatial mismatch between natural forest patterns and population heat exposure trajectories. DNF is largely concentrated in the sparsely populated western mountainous and plateau regions; conversely, deciduous broadleaf forest is extensively rooted in the high-POP-density plains and hilly areas of the East. The vast base of the exposed population allows its cooling services to be transformed into substantial tangible wellbeing. This significant supply-demand inversion suggests that future urban forest planning should not blindly pursue absolute cooling peaks but must instead adhere to the principle of spatial supply-demand adaptation for precise resource allocation.

Non-linear coupled driving mechanisms of multiple factors on cooling services

The multi-scale differentiation of cooling service is rooted in natural geographic patterns, responding to urbanization and ecological policies, aligning with global research findings (Li H. et al., 2024).

Global dimension: dynamic game between ecological background and human activities

Overall, with the exception of RHU, all other factors show significant nonlinear response characteristics at 95% confidence level. Among them, natural factors exhibit varying negative correlations, while socio-economic factors all demonstrate positive driving effects (Figure 10).

The interaction between natural and socio-economic factors determines cooling services strength, with policy interventions indirectly regulating their balance. DEM is the strongest negative driver, its U-shaped influence curve reflecting how topography shapes temperature background, closely linked to ecological protection policies.

High-altitude areas are strictly controlled, ensuring the integrity and connectivity of forests, which amplifies cooling effects. In low-altitude plains, despite high human activity, policies like urban greenway network construction improved the diffusion efficiency of cool air. In medium-altitude regions, where urban expansion is prioritized, construction land occupies forests, resulting in diluted cooling effects. This pattern strongly aligns with Liu et al.‘s terrain-policy-heat environment coupling theory (Liu et al., 2022). WS had a weak influence, as it was counteracted by the effects of urbanization policies and natural factors: on one hand, high-density urban development significantly weakened wind speed due to the wind blocking effect of buildings; on the other hand, some cities, such as Shanghai and Guangzhou, implemented wind corridor planning to mitigate wind speed decay in local areas, but this had not yet formed a dominant effect at global scale (Li et al., 2025).

Socio-economic factors’ positive driving mainly stem from the increasing dominance of human activity heat emissions. Yang et al.‘s study on Tianjin confirmed that the rapid increase in NDBI is a key driver of the decline in cooling capacity (Yang et al., 2023).

Climate zone dimension: regulation of driver strength by hydrothermal background

The intensity ranking of driving factors across climate zones generally aligns with global scale, but slight regional differences in natural constraints and human-driven effects exist, primarily driven by hydrothermal conditions. In semi-arid regions, the cooling advantage of high-altitude areas was weakened by precipitation scarcity. In contrast, RHU played a more significant cooling role, as dry air in arid regions means even slight humidity increases significantly enhanced vegetation’s transpiration CC (Li H. et al., 2024). Human activities (NTL and POP) had a more prominent warming effect in humid regions. This conclusion is consistent with Shi et al.‘s findings, confirming that cooling services in cities in humid regions are more susceptible to human disturbance (Shi and Jia, 2024).

Forest form dimension: regulation of driving mechanisms by physiological functions

Forest forms’ physiological structures, especially canopy structure and leaf characteristics, directly dictate tree responses to external drivers. Deciduous coniferous forests exhibit optimal shading and cooling due to dense leaf arrangement, and effect further amplified by high-altitude low-temperature backgrounds. RHU in deciduous coniferous forests has a greater cooling effect than in others, confirming Li et al.‘s research that coniferous forest leaf structures are more sensitive to humidity changes (Li J. et al., 2024).

Evolutionary trends in environmental equity of cooling services

The average Gini coefficient for all cities was 0.844, with high values concentrated in the middle and lower reaches of the Yangtze River. Over time, the Gini coefficient fluctuated in high value range of 0.82–0.86, but showed a decreasing trend after 2015, indicating improved distributional equity in some regions. In 2001, only Naqu City had inequality; from 2005 to 2010, all cities were highly unequal; in 2015, cities like Neijiang and Zhoukou entered a moderately unequal category, and by 2022, 8 cities had been classified, confirming the trend of alleviating cooling service inequality. Consequently, we selected Naqu, Zhoukou, and Neijiang as typical case studies to analyze changes in cooling service fairness and their driving mechanisms (Figures 11a,b).

Representing semi-arid regions (Figure 11c), Naqu City has a high plateau climate that leads to an extremely high CE (11.37), reflecting extreme ecological vulnerability (Figure 11d). The fragmented and isolated forest landscape structure amplified the oasis effect (Shi et al., 2022). The Gini coefficient increased from 0.49 in 2001 to 0.939 in 2022, with fairness continuously deteriorating. The point-source research-based greening policies, limited by ecological red lines and technical barriers, failed to achieve universal cooling resource coverage, increasing fairness gaps.

Representing semi-humid regions, Zhoukou City has a plain forest network structure that provides a natural foundation for fairness. In 2015, the Gini coefficient dropped to 0.48 and stabilized around 0.46, thanks to the Create Forest action released in 2018 (Figure 11e). The agriculture-forest interwoven, forest-city integration strategy effectively filled the cooling service coverage gap, achieving a leap from insufficient green space to wide-area coverage (Zhu et al., 2022).

Representing humid regions, Neijiang City has favorable hydrothermal hydrothermal conditions and vegetation background, maintaining stable CE values (0.89–1.26). However, due to extreme dry-heat events, CE spiked to 2.73 in 2022 (Figure 11f). Due to topographical limitations, forests are mainly distributed in mountainous areas (Li et al., 2022), resulting in a near green, no cooling phenomenon. The Gini coefficient fluctuated around 0.5 for a long time, dropping to 0.42 in 2022, the lowest value among the three typical cities. As a node in Chengdu-Chongqing economic circle, Neijiang has developed high-density urban areas on limited flat land. Policies such as green insertion and 300-m green coverage have extended forest resources from large parks to community micro-green spaces, effectively reducing the Gini coefficient.

Recommendations for sustainable management

Based on spatiotemporal heterogeneity, driving mechanisms, and environmental equity evolutionary characteristics of cooling services, targeted sustainable development strategies are proposed to enhance urban heat environment resilience and achieve balanced cooling service provision.

First, optimizing the spatial configuration of UF to achieve precise supply-demand matching remains a core priority. Target-oriented configurations of tree species characterized by high CE should be implemented within core UHI zones based on the spatial distribution of LST. At the macro scale, trans-regional ecological corridors should be constructed across administrative boundaries to connect intra-urban forest patches with peripheral ecological forests. This approach aims to enhance the COHESION and minimize the SPLIT, thereby facilitating urban-rural integration and corridor connectivity to establish a universally covered and seamlessly interconnected cooling service network (Li and Liu, 2025; Hu et al., 2026). At the micro scale, the supply of forest land should be expanded in regions with high population density and low CB values. By developing 15-min green space living circles, a population-oriented forest allocation can be realized, ultimately mitigating environmental exposure inequality.

Second, differentiated forest planning strategies must be executed across distinct climate zones. Within the arid zone, planning should focus on the dual drivers of water conservation or efficiency enhancement and inclusive coverage. UFCI can be elevated by utilizing drought-resistant arbor-shrub composite communities, while small, contiguous forest patches should be developed in densely populated areas to bridge the equity gap characterized by “available cooling sources but absent services.” Within the semi-arid zone, synergistic efforts should be directed toward controlling the expansion of built-up areas, increasing forest land area, and enhancing shelterbelt network coverage. In peri-urban areas, preserving contiguous forest lands is essential to elevate the PLAND and the LPI. Concurrently, within built-up areas, urban greening initiatives—such as replacing illegal structures with forests—should be promoted to establish small forest patches utilizing vacant land parcels. This approach aims to increase PD, thereby enhancing the spatial coverage of cooling services. Drawing on the “plain forest network” experience of Zhoukou City, grid-like shelterbelt belts should be established between farmlands and rural settlements to achieve wide-area coverage of cooling services and arrest the equity degradation trend exemplified by Nagqu City. Within the humid zone and semi-humid zone, quality enhancement, efficiency optimization, and balanced spatial distribution must be advanced simultaneously. Given the regional characteristics of low CE but substantial CC potential, forest stand structures should be optimized by increasing the proportions of deciduous needleleaf forests and mixed forests to enhance the cooling efficiency per unit of forest land. In regions with high Gini coefficients, such as the middle and lower reaches of the Yangtze River, supplementary forest lands should be deployed in service-deficient areas, including urban-rural fringes and low-income communities. Relying on river and transportation corridors, waterfront greenways and ventilation corridors should be constructed to unblock the connectivity pathways of cooling sources, alleviate the spatial mismatch of “high population vs. low cooling,” and improve the equity of CB (Schnabel et al., 2025).

Third, optimizing forest stand structures and strengthening functional synergy strategies are essential. The mixed cultivation mode of deciduous needleleaf forests and deciduous broadleaf forests should be widely promoted to balance the complementary advantages of CE and CB. In core urban areas and densely populated zones, priority should be given to planting deciduous broadleaf forests to leverage their distributional flexibility for improved supply-demand matching. Conversely, deciduous needleleaf forests should be cultivated in peripheral ecological corridors to exploit their superior CC performance. Through the scientifically rigorous configuration of internal and external community structures, a “core service - peripheral assurance” forest stand pattern can be constructed, thereby augmenting the alignment between CE and CB and ensuring the sustainable and inclusive delivery of cooling services.

Limitations and future study

This study still has several limitations that need to be addressed. First, regarding welfare distribution, while we confirmed the spatial mismatch of cooling services at a national scale and revealed broad environmental spatial inequity, limited access to large-scale socioeconomic gridded data prevented a deeper analysis of micro-level disparities in cooling service access among different social classes or vulnerable groups. Future studies should integrate fine-grained survey data to further explore environmental justice across internal demographic structures (Quan et al., 2025). Second, concerning mechanism exploration, key micro-level variables of internal urban morphology were not included; future research must leverage multi-source data to uncover deeper cooling mechanisms. Third, at the practical application level, this study relied solely on data-driven simulation and fitting methods. Lacking validation and support from empirical field observations, the actual explanatory power and applicability of these findings require further verification.

Taking Chinese cities as the research object and adopting a dual perspective of climate zones and forest types, this study systematically revealed the spatiotemporal heterogeneity, driving mechanisms, and environmental equity evolutionary characteristics of UF cooling services. The primary conclusions are as follows: Analysis of the spatiotemporal heterogeneity of urban cooling services. Urban forest cooling services in China exhibit significant heterogeneity and phased evolutionary patterns across spatiotemporal dimensions. Throughout the research period, although the negative correlation clustering between UFCI and LST persisted, this natural physical regulation efficacy has been progressively weakened by the intensified heat emissions from continuous high-density urban expansion. This reveals a potent dynamic trade-off between “ecological supply” and “human interference.” Regarding spatial patterns, a significant “supply-demand inversion” was observed between physical supply (CE and CC) and social wellbeing (CB), reflecting the spatial imbalance between the natural ecological baseline and population heat exposure trajectories.

Coupled regulation mechanisms of climate background and forest forms. The heterogeneity of cooling services results from the combined effects of climatic background conditions, forest form characteristics, and human activities. At the macro-climatic scale, the oasis effect in the ARZ amplifies cooling potential, whereas the high-humidity conditions in the HUZ inhibit the transpirational cooling of vegetation. This vividly illustrates the threshold constraints imposed by hydrothermal conditions on vegetative cooling functions. At the micro-forest form scale, the high cooling efficiency of needleleaf forests stems primarily from the shading advantages of their canopy structures (highest CE and CC). Conversely, broadleaf forests maximize cooling wellbeing (CB) due to their superior spatial adaptation to population distributions, reflecting the matching logic between vegetation physiological traits and service demand. GTNNWR model results further demonstrate that this regulatory effect is non-linear, with its influence intensity dynamically adjusting across spatiotemporal contexts, revealing the scale dependence of the non-linear driving forces exerted by natural and socio-economic factors.

Evolution of environmental equity and extreme spatial mismatch under dynamic heat exposure. Urban cooling services have long remained in a state of severe exposure inequality, the core crux of which lies in the extreme spatial mismatch between the static distribution of ecological cooling sources and the dynamic heat-protection demands of the population. Although equity has shown a gradual alleviating trend since 2015 driven by ecological policies, the deprivation caused by “high supply, low efficiency” mismatches remains a bottleneck restricting environmental justice.

This study constructed a research framework of “global model recognition-local mechanism analysis”, providing a systematic methodology for investigating UF cooling services.