We apply ordinary least squares (OLS) and three spatial regression models, including spatial lag model (SLM), spatial error model (SEM), and spatial Durbin model (SDM), to explore the linear relationships between urban form and CO₂ emissions. The spatial models account for spatial dependence, which OLS overlooks. And below is a summary of the model formulations.

The SLM posits that a city’s CO₂ emissions are influenced not only by its own morphological characteristics but also by the emissions of neighboring cities, thereby capturing the spatial autocorrelation of the dependent variable (Gao et al., 2020). The model is expressed as follows (Equation 1): y = ρ W y + X β + ϵ , ϵ ∼ N 0 , σ 2 I (1) where y is the vector of ln (CO₂) emissions, W y represents the spatial lag of the dependent variable, ρ is the spatial autoregressive parameter, X is the matrix of morphological predictors with log transformations, β is the vector of regression coefficients, and I is the identity matrix, and ϵ is the error term. The spatial weights matrix W is defined based on queen contiguity, meaning that cities sharing a boundary are considered neighbors.

The SEM accounts for autocorrelation in error terms by assuming that unobserved factors influencing emissions are spatially correlated (Golgher and Voss, 2016). Its formulation is as follows (Equation 2): y = X β + u , u = λ W u + ϵ , ϵ ∼ N 0 , σ 2 I (2) where u is the spatially dependent error term, λ is the spatial error coefficient, and W u represents the spatially lagged errors. This model is particularly beneficial for addressing spatial dependence due to omitted variables or measurement errors with spatial patterns.

The SDM extends the SLM by incorporating spatially lagged variables, acknowledging that emissions and morphological features of neighboring cities can influence the emissions of the target city (Li and Li, 2020). The model is expressed in Equation 3: y = ρ W y + X β + W X γ + ϵ , ϵ ∼ N 0 , σ 2 I (3) where W X is the matrix of spatially lagged covariates, and γ is the vector of coefficients for these covariates. The SDM accommodates both endogenous and exogenous spatial effects via W y and W X , respectively, thereby providing a more comprehensive analysis at the expense of increased model complexity.

Model estimation is conducted using maximum likelihood, with the weight matrix W row-standardized for appropriate interpretation of spatial effects. Models are evaluated based on R², adjusted R², and Akaike information criterion (AIC). Additionally, Moran’s I test is performed to assess residual spatial autocorrelation. Variance inflation factors (VIFs) are calculated iteratively until correlations among predictors are below a threshold of 5, ensuring a stable final set of predictors.

Geographical Gaussian process regression (GGPR) represents an advanced extension of Gaussian process regression specifically designed to address fundamental challenges in geospatial data analysis, including spatial autocorrelation, small sample prediction, and model interpretability (Jiao and Tao, 2025). Traditional machine learning methods such as random forests operate under the assumption of independence and identical distribution (i.i.d.), which fundamentally contradicts the spatial autocorrelation principles described by Tobler’s First Law, consequently leading to biased predictions (Dale and Fortin, 2014). GGPR addresses this limitation through a non-parametric Gaussian process framework that explicitly captures spatial dependencies using specialized covariance functions.

The mathematical foundation of GGPR rests on the Gaussian process definition, characterized by a mean function m x and a covariance function k x , x ′ . For a given set of input features X and corresponding target values y , the GGPR model provides predictions for the posterior mean and variance at an unknown point x * through Equation 4: μ x * = K x * , X K X , X − 1 y , σ 2 x * = K x * , x * − K x * , X K X , X − 1 K X , x * (4) where K denotes the covariance matrix defined by the kernel functions. The key innovation of GGPR lies in its composite kernel design that effectively captures the multi-scale nature of urban carbon emissions. This approach combines a spatial similarity kernel (SSK) with a Matérn kernel (MK) through the formulation S S K + w × M K , where w represents a weighting parameter. The S S K parameterizes spatial similarity through geographical configuration vectors e = e 1 , e 2 , … , e n , measuring the similarity between spatial units i and j using Equation 5: S i , j = P · ∏ v = 1 n   exp ⁡ − e v i − e v j 2 2 σ v 2 (5) where P · is a weighted averaging function and σ v parameters are optimized using limited-memory Broyden-Fletcher-Goldfarb-Shanno bound or Bayesian optimization methods (L-BFGS-B) (Morita et al., 2022). This kernel captures broad geographical patterns reflecting regional energy infrastructure and policy effects. Complementing the SSK, the MK model’s local spatial autocorrelation as a function of Euclidean distance r through Equation 6: k r = 2 1 − ν Γ ν 2 ν r l ν K ν 2 ν r l (6) where hyperparameters l (length scale) and ν (smoothness) are learned during training. This component effectively captures fine-scale neighborhood effects arising from local urban morphological variations. The composite kernel approach addresses the inherent challenge that single-kernel methods face: SSK alone would miss local spatial structure, while MK alone would fail to capture broader geographical patterns. Empirical validation confirms the effectiveness of this composite kernel strategy. Across the years 2005–2020, it achieves an average reduction of 43.9% in cross-validation root mean square error (RMSE) compared to the SSK alone and an average reduction of 2.4% compared to the MK alone, demonstrating substantial improvement over the SSK and modest gains or comparable performance relative to the MK depending on the year.

The hyperparameter optimization procedure employs Bayesian optimization with Gaussian process priors over the hyperparameter space, utilizing expected improvement as the acquisition function. The optimization procedure runs for 100 iterations with 5-fold spatial cross-validation, searching across a carefully defined parameter space that includes S S K parameters σ v ∈ [0.1, 10] for each geographical configuration dimension, Matérn kernel parameters with length scale l ∈ [10, 500] and smoothness ν ∈ [0.5, 5], and composite weight w ∈ [0.05, 0.2] initialized at 0.1. All parameters are constrained using the L-BFGS-B algorithm to ensure both numerical stability and physical interpretability throughout the optimization process.

To enhance model interpretability, GGPR is integrated with GeoShapley, a game theory-based framework that extends the Shapley value methodology to quantify spatial effects in machine learning models (Li, 2024). GeoShapley conceptualizes location features as joint players in a model prediction game, treating coordinates such as latitude and longitude as a single joint player rather than individual features. This approach enables the quantification of both intrinsic location effects and synergies between location and other features. The GeoShapley value for location features is calculated using Equation 7: ϕ G E O = ∑ S ⊆ M ∖ G E O   s ! p − s − g ! p − g + 1 ! f S ∪ G E O − f S (7) where GEO denotes a set of location features of size g .   S represents a feature subset excluding GEO with size s .   M is the set of all features, and p is the total number of features in the model. f · is the model prediction function.

For non-spatial features, the GeoShapley value is given by Equation 8: ϕ j = ∑ S ⊆ M ∖ j   s ! p − s − g ! p − g + 1 ! f S ∪ j − f S (8)

The interaction effect between the location feature set and a specific feature j is computed using Equation 9: ϕ GEO , j = ∑ S ⊆ M ∖ GEO ∪ j   s ! p − s − g − 1 ! p − g + 1 ! · Δ GEO , j (9) where the interaction term Δ GEO , j is defined using Equation 10: Δ GEO , j = f S ∪ GEO , j − f S ∪ GEO − f S ∪ j + f S (10)

The final model prediction y ^ can then be decomposed using Equation 11: y ^ = ϕ 0 + ϕ GEO + ∑ j = 1 p   ϕ j + ∑ j = 1 p   ϕ GEO , j (11) where ϕ 0 is a constant base value; ϕ GEO   represents the intrinsic effect of location features; ϕ j captures the location-invariant effect of each non-spatial feature j ; and ϕ GEO , j quantifies the spatially varying interaction effect between the location features and feature j (Li, 2024). The GeoShapley values are estimated using an extended kernel SHAP method through a weighted least squares framework, making the approach computationally feasible for practical applications (Chen Y. et al., 2025). By combining GeoShapley with GGPR, the proposed framework simultaneously captures complex spatial dependencies and interpretable feature contributions, providing comprehensive insights into spatial heterogeneity and nonlinear relationships in geospatial phenomena, as summarized in Figure 1.

The land use data from Yang and Huang is used to compute two-dimensional urban morphological indices that describe various aspects of urban spatial structure (Yang and Huang, 2021). These indices, derived from existing studies that measure the shape, arrangement, and diversity of impervious surfaces and land use types within a city (Moudon, 1994; McGarigal et al., 2002; Dai et al., 2018; Jia et al., 2019; Olsen et al., 2019; Ge et al., 2021), can be organized into four functional categories as indicated in Table 1. Compactness and concentration indicators capture how concentrated and compact the urban form is: the perimeter-area fractal dimension (PAF) measures the shape complexity of impervious surface patches, with higher values indicating more compact, less irregular shapes. Whilst the urban aggregation index (UAI) evaluates the clustering of patches, with higher values reflecting more aggregated, compact urban areas. The largest patch index (LPI) quantifies the dominance of the largest patch relative to the total urban area, where higher values suggest a consolidated urban core, whereas and the total core area (TCA) represents the sum of stable core areas within a 10 km buffer of the city center, reflecting the extent of robust urban patches. And the built area dynamic aggregation index (ADA) emphasizes the horizontal footprints of buildings and their radial concentration pattern. Fragmentation and distribution indicators measure dispersion and spatial distribution patterns. For instance, the patch density index (PDI) assesses fragmentation by calculating the number of patches per unit area, with higher values indicating a more dispersed urban form. While the mean elongation index (MEI) quantifies the average elongation of patches, where higher values suggest linear structures (e.g., roads) and lower values indicate compact forms (e.g., impervious surfaces). The patch cohesion index (COH), as a connectivity and integration indicator, measures the connectivity of patches, with higher values indicating well-connected urban areas that may enhance infrastructure efficiency, whilst the cumulative building coverage continuity (CBC) serves as an accumulated indicator of built density decay based on distance-weighted information concerning the area occupied by buildings. Finally, Shannon’s diversity index (SDI), representing the diversity indicator, evaluates the diversity of land use types within the urban boundary, with higher values indicating greater land use variety.

In this study, 3D morphological indices are computed using an innovative concentric zone approach. The analysis focuses on 286 Chinese cities and their government buildings as of 2005. We plot 30 concentric circles at 1 km intervals, radiating outward from the centers of these cities (Figure 2). The selection of government office locations as city center is significant, as urban development typically extends from historic city centers where municipal offices are located, thereby reflecting the spatial expansion patterns of Chinese cities. A 30 km concentric circular buffer is methodologically appropriate for delineating urban extent because it aligns with the spatial zone of most rapid land expansion. In Beijing, land development intensity peaked within the 10–30 km buffer from the city center, marking it as the core zone of outward urban growth and suburban transition (Huang et al., 2020). Similarly, gradient-based urban studies have operationalized buffers up to 30 km, categorizing 0–14 km as urban core, 14–24 km as suburban, and 24–30 km as fringe areas, effectively distinguishing morphological zones within the urban gradient (Ma et al., 2024). Moreover, the 30 km scale bridges neighborhood-focused analysis and city-region planning, offering both conceptual and empirical relevance.

These circles define 30 circular geographic buffer zones surrounding each city’s core. We employ building footprints and height data from 2005 to 2020 in 5-year intervals, examining building coverage area, volume, and average height statistics across various buffer zones. Each city is represented by 30 measurements per year. Figure 2 illustrates the building height distribution across the concentric zones of Shanghai. Building height data were sourced from Chen et al. (Chen P. et al., 2025), who utilized the ALOS World 3D 30 m digital surface model (DSM) alongside NASA’s digital elevation model (DEM), global artificial impervious area, and global roads inventory project data. A normalized DSM was generated by calculating the difference between the digital terrain model (DTM) and DSM on Google Earth Engine, resulting in a 30 m national map. Validation against field measurements demonstrated good accuracy for flat/mixed roofs (R² = 0.81, RMSE = 4.26–3.51 m) across 6 Chinese cities (R² = 0.64–0.79, RMSE = 4.65–5.35 m), yielding an overall accuracy of RMSE = 4.98 m and R² = 0.75 (Huang et al., 2022). Data for Chinese cities were collected for 2005, 2010, 2015, and 2020 for this analysis.

The 3D morphological indices employed in this study are determined by referring to similar 3D urban form measurement concepts from previous studies (Li and Zeng, 2024; Liu et al., 2021; Guo et al., 2022), and can be categorized into 3 groups (Table 2). Compactness and concentration indicators in three dimensions focus more on the radial concentration patterns. The volumetric dynamic aggregation index (VDAI) quantifies the radial decrease in volumetric mass of buildings, measuring how building volume is concentrated toward the urban center, and the edge-core volume ratio (ECVR) compares peripheral and core volumes to evaluate the extent of suburban sprawl relative to central concentration. Fragmentation and distribution indicators capture three-dimensional dispersion patterns: the volumetric gradient kurtosis (VGKI) assesses the peakedness of volume distribution along the urban gradient, facilitating differentiation between monocentric and polycentric structures. Whilst the volumetric gradient decay rate (VGDR) measures the average decrease in volume intensity from the city center to the surrounding areas. Connectivity and integration indicators measure synchronized development: the urban gradient synergy index (UGSI) captures the interaction between changes in building heights and coverage, revealing whether vertical and horizontal densities decline concurrently with distance.

The open data inventory for anthropogenic CO₂ (ODIAC) provides one of the most detailed global estimates of fossil fuel CO₂ emissions, with a spatial resolution of 1 km² (Oda and Maksyutov, 2011). Developed by Japan’s national institute for environmental studies, this inventory uniquely integrates satellite-based nighttime light data with an extensive database of individual power plant emissions and their locations. National totals, derived from energy statistics, are allocated to non-point sources based on calibrated nighttime light radiance, while emissions from more than 17,000 individual power plants are precisely located at their geographic coordinates. This hybrid approach transcends the spatial limitations of population-based or coarse-resolution inventories, accommodating both diffuse urban footprints and concentrated point sources. First published in 2011, ODIAC has undergone numerous refinements and updates since its inception. It now includes NASA’s Black Marble nighttime light product, enhancing estimates in response to recent developments in energy use and urbanization. The ODIAC data can support a variety of applications, including atmospheric CO₂ flux inversions and city-scale emission mapping, and serves as a benchmark for fine spatial resolution, globally consistent anthropogenic CO₂ modeling. For this study, we utilized a 30 km buffer around the city center for extracting CO₂ emissions data from 2005 to 2020, to align with the spatial and temporal scopes of the urban morphological indices.

The global Moran’s I values for CO₂ emissions in 2005, 2010, 2015, and 2020 ranged from 0.255 to 0.257, which indicated significant positive autocorrelation with z-scores exceeding 12.890 (Table 3). Local Moran’s I statistics revealed distinct spatial patterns across China. High-High clusters were concentrated in eastern and central provinces, with Shandong, Jiangsu, and Henan identified as high-emission areas. In contrast, Low-Low clusters dominated the western regions, particularly in Tibet and Qinghai, which were characterized as low-emission areas. High-Low and Low-High outliers were observed in transitional zones between these clusters, suggesting that the spatial aggregation of high CO₂ emissions was significantly influenced by economic and industrial development patterns in eastern China, as illustrated in Figure 3.

Morphological characteristics varied significantly among cities, as detailed in Table 4. PAF exhibited an average fusiform shape complexity with a mean of 1.33. MEI averaged 866.65, demonstrating drastic heterogeneity across cities. TCA was notably heterogeneous with a mean of 23,313.93, while PDI displayed a wide range, reflecting urban density variations. Aggregation and connectivity were strong across the studied cities. UAI showed a mean of 85.52, and COH averaged 98.03, indicating well-connected urban areas. LPI presented the greatest variability, reaching up to 50.04, whereas SDI showed a mean of 0.262, suggesting low to moderate diversity. Both VDAI and ADA were low and consistent, while the ECVR was highly variable.

CO₂ emissions displayed significant positive correlations with several morphological indices. The strongest correlations were observed with TCA at 0.808, LPI at 0.805, and CBC at 0.739, suggesting that emissions increased with larger and more aggregated patches. Moderate correlations were found with UAI at 0.444 and SDI at 0.696, while weak-to-moderate correlations were observed with PDI at 0.121 and COH at 0.313. Weak negative correlations were noted with PAF and MEI, though the correlation with PAF was not significant. UGSI also showed a weak negative correlation at −0.199. The indices VDAI, ADA, VGDR, and ECVR showed very weak correlations with emissions. Notably, TCA emerged as a strong predictor of both LPI with a correlation of 0.954 and SDI with a correlation of 0.930, as shown in Figure 4.

The spatial regression analysis showed consistent patterns across the study period, with the SEM emerged as the superior model. These findings illuminated the fundamental spatial dependence structure underlying urban carbon emissions and showed how morphological factors operated through spatially mediated mechanisms. In 2005, whilst the OLS model yielded an R² of 0.642, the SEM outperformed other models with an R² of 0.701. The improvement over OLS (ΔR² = 0.059) indicated that approximately 6% of emission variance was explained by spatially correlated unobserved factors, suggesting that urban carbon dynamics were fundamentally influenced by spillover effects from neighboring cities through shared infrastructure networks, economic linkages, and regional industrial policies. The SEM included a spatial error term λ of 0.487 and demonstrated an AIC of 485.56, thereby effectively eliminating the residual autocorrelation. Although the SDM performed slightly better with an R² of 0.702, its higher AIC of 497.99 and the fact that only the lagged UAI proved significant reinforced the preference for the SEM (Table 5).

The temporal evolution of model performance displayed critical insights into China’s changing urban carbon landscape. For 2010, the OLS model had an R² of 0.646. The SEM produced an R² of 0.686 with a λ of 0.414 and an AIC of 486.90, showing no residual autocorrelation. The decline in spatial error coefficient from 0.487 to 0.414 suggested that inter-city spatial dependence in carbon emissions had weakened over time, which potentially reflected the increasing importance of city-specific policies and technological innovations that had differentiated emission patterns across neighboring urban areas. In 2015, the OLS model recorded an R² of 0.624. The SEM achieved an R² of 0.662 with λ of 0.395 and AIC of 506.11. The continued decline in spatial dependence coupled with slightly reduced overall explanatory power suggested the increased urban carbon heterogeneity, where local morphological characteristics had begun to dominate regional patterns. By 2020, the OLS model had an R² of 0.616. The SEM fitted with an R² of 0.648, a λ of 0.367, and an AIC of 515.27. The further reduction in spatial dependence (λ = 0.367) reflected the maturation of China’s urban system, where cities had developed increasingly distinct carbon profiles based on their specific morphological characteristics and local policy environments.

Across all years studied, the SEM consistently provided the best fit. This superiority was evidenced by the lowest AIC values, ranging from 485.56 in 2005 to 515.27 in 2020, and a lack of residual spatial autocorrelation with Moran’s I p-values exceeding 0.5, compared to OLS, SLM, and SDM models. While the SDM offered marginal R² improvements, these were accompanied by higher AIC values and generally insignificant lagged covariates. Likelihood ratio tests confirmed the significance of spatial effects. The SEM’s ability to account for spatially dependent errors while avoiding the empirical need for lagged covariates rendered it the most parsimonious model. Although the SLM outperformed the OLS model, it still exhibited residual autocorrelation, which further supported the choice of SEM. The insignificance of spatially lagged covariates ( W X ) in the SDM, despite theoretical expectations of spillover effects, suggested that carbon emissions were more influenced by spatially correlated unobserved factors (captured by the spatial error term λ) rather than direct spillovers from neighboring cities’ morphological characteristics. This pattern was consistent with emissions having been driven by regional energy infrastructure and economic structures that created spatial clustering in the error terms rather than direct morphological contagion effects.

Elasticity, defined as the percentage change in CO₂ emissions in response to a 1% change in the predictor variable, could be interpreted through the SEM coefficients as reflecting the relationship between urban compactness and CO₂ emissions. Among compactness indicators, UAI emerged as the most influential predictor of CO₂ emissions. The elasticity of emissions with respect to UAI increased from 0.754 in 2005 to 0.781 in 2020, indicating an intensifying relationship where urban clustering became increasingly carbon-intensive. This strengthening effect suggested that while compact development initially provided efficiency gains, Chinese cities may have been approaching or exceeding optimal density thresholds where agglomeration diseconomies, including congestion, concentrated industrial activities, and infrastructure stress, began to outweigh the carbon benefits of compactness.

TCA’s moderate significance primarily in 2010 reflected the critical period when China’s urbanization accelerated after 2008. During this phase, core area expansion became a key determinant of emission patterns, suggesting that the spatial organization of urban centers played a crucial role in carbon dynamics during rapid growth periods. The temporal dynamics of ADA demonstrated a fundamental transition in how directional urban expansion influenced emissions. ADA depicted significant negative effects in the early study years, with a 1% increase associated with 0.513% emission reduction in 2005 and 0.329% reduction in 2010. Nevertheless, this effect became insignificant by 2015–2020. This evolution suggested that during the initial phase (2005–2010), directional compact growth effectively reduced emissions through infrastructure efficiency and reduced sprawl. However, as cities matured and growth patterns stabilized, the carbon benefits of directional development diminished.

UGSI’s evolution represented perhaps the most significant finding for urban carbon management. Its negative impact intensified dramatically over time, with elasticity increasing from −0.096 in 2005 to −0.297 in 2020. This tripling effect revealed that synchronized vertical-horizontal development had become increasingly critical as Chinese cities transitioned from extensive to intensive growth models. The intensifying effectiveness suggested that as urban land became scarce and development became more three-dimensional, the coordination between building height and coverage became a primary determinant of carbon efficiency.

The GGPR filtered morphological variables for each year of the study period by computing mutual information (MI) scores between CO₂ emissions and all predictors. Morphological variables with MI scores above a threshold of 0.1 were selected, to ensure only those with significant predictive power were retained for modeling, while less informative variables were excluded if their MI scores fell below the threshold (Paninski, 2003). The GGPR model demonstrated exceptionally strong performance across all years with consistently low mean absolute percentage error (MAPE) values ranging from 1.6% to 2.4%, indicating excellent relative accuracy in predicting CO₂ emissions. While absolute error metrics increased over time, with mean absolute error (MAE) rising from 43,634 in 2005 to 77,497 in 2020 (a 78% increase) and RMSE increasing from 76,432 to 122,961 (a 61% increase), the R² values remained remarkably stable at 0.999 across all years from 2005 to 2020 (Table 6).

The GeoShapley analysis revealed complex spatial patterns in how urban morphological characteristics influenced carbon emissions across Chinese cities between 2005 and 2020, as seen in Figures 5, 6, Supplementary Figures S1, S2. Contrary to conventional expectations, several key urban form indicators demonstrated predominantly negative contributions to emissions, suggesting that certain aspects of urbanization actually mitigated carbon output. LPI exhibited a consistently negative impact on emissions, with mean GeoShapley values of −1.924 × 10⁶ (±1.256 × 10⁶) in 2005 evolving to −1.300 × 10⁶ (±1.506 × 10⁶) in 2020. This finding indicated that larger continuous urban patches facilitated more efficient infrastructure and reduced per capita emissions. Only 12.9% of cities showed positive LPI values by 2020, concentrated primarily in rapidly expanding second-tier cities. The spatial distribution showed a clear east-west gradient, with eastern coastal megacities showing less negative values while western cities demonstrated strongly negative values.

COH presented a nuanced picture, transitioning from slightly negative mean values of −0.106 × 10⁶ (±0.337 × 10⁶) in 2005 to −0.103 × 10⁶ (±0.387 × 10⁶) in 2020. Approximately 48.6% of cities showed positive COH values by 2020, indicating substantial spatial heterogeneity. Eastern metropolitan areas demonstrated positive contributions, while central and western cities predominantly showed negative values, suggesting that urban connectivity’s impact on emissions was highly context-dependent. TCA showed the strongest negative association with emissions among compactness indicators, with mean values of −1.053 × 10⁶ (±0.610 × 10⁶) in 2005 deepening to −1.226 × 10⁶ (±1.503 × 10⁶) by 2020. This increased variance indicated growing spatial differentiation, with only 21% of cities showing positive TCA values in 2020, primarily provincial capitals and major economic centers where larger urban cores concentrated high-emission activities. Fragmentation indicators revealed important regional variations. PDI demonstrated predominantly negative but weak effects with mean values of −0.069 × 10⁶ (±0.441 × 10⁶) in 2020. The quartile distribution (Q25 = −0.292 × 10⁶, Q50 = −0.059 × 10⁶, Q75 = 0.121 × 10⁶) revealed significant regional variation, with central Chinese cities, particularly in Henan and Hubei provinces showing positive PDI values, suggesting that fragmented development in these rapidly urbanizing regions increased emissions. SDI maintained consistently negative impacts, with mean values of −0.768 × 10⁶ (±0.512 × 10⁶) in 2005 shifting to −0.706 × 10⁶ (±1.059 × 10⁶) in 2020.

Meanwhile, 3D urban form indicators demonstrated evolving patterns. VDAI showed a notable reversal from negative values of −0.107 × 10⁶ (±0.335 × 10⁶) in 2005 to slightly positive values of 0.025 × 10⁶ (±0.326 × 10⁶) by 2020, representing the only indicator to transition from negative to positive mean impact. This suggested that asymmetric vertical development patterns were increasingly associated with higher emissions, particularly in rapidly verticalizing cities. ECVR showed mean values of −0.055 × 10⁶ (±0.711 × 10⁶) by 2020, with substantial spatial variation. Eastern coastal cities tended toward positive values while inland cities showed negative contributions, reflecting different peripheral development patterns.

The spatial distribution maps from 2005 to 2020 revealed evolving regional patterns. Eastern coastal regions consistently showed less negative or slightly positive GeoShapley values across most indicators, suggesting that mature urbanization in these areas had reached a point where further densification might increase emissions. Central regions exhibited high variance across all indicators, reflecting cities in rapid transition with diverse development trajectories. Western regions consistently showed strongly negative values across most indicators, showing that urbanization patterns in these cities helped mitigate emissions, possibly due to lower development intensity. Northeast regions demonstrated unique patterns with increasingly negative values over time, potentially reflecting economic restructuring and population decline. One of the most interesting findings was the paradigm shift in understanding urban compactness impacts. Moreover, compactness indicators (LPI, TCA, CBC) predominantly showed negative associations with emissions, with CBC values of −1.149 × 10⁶ (±0.681 × 10⁶) in 2005 evolving to −0.851 × 10⁶ (±1.191 × 10⁶) in 2020. The coefficient of variation for most indicators increased between 2005 and 2020, indicating growing spatial differentiation in how urban form influenced emissions. The impact of urban morphology on emissions proved highly dependent on regional development stage, with mature eastern cities showing different patterns than developing central and western cities. While horizontal compactness reduced emissions, asymmetric vertical development showed increasing positive associations with emissions, which clearly suggested the need for balanced three-dimensional urban planning.

Feature importance analysis revealed a fundamental reorganization of urban carbon drivers (Figure 7; Supplementary Figures S3–S5). The LPI consistently dominated as the most important predictor over the study period, with stable importance values ranging from 1.820 to 2.152 × 10⁶ between 2005 and 2020 (Figure 8). This consistency suggested that patch dominance served as a core organizing principle in urban carbon systems, persisting despite temporal variations in technology, policy, and economic conditions. The observed stability of LPI’s influence indicated that the spatial configuration of the largest urban patches produced path-dependent carbon emission patterns that remained robust in the face of shorter-term interventions. A dramatic rise of TCA was also detected: it grew from a secondary importance of 1.146 × 10⁶ in 2005 to nearly matching LPI (1.695 × 10⁶) by 2020 (Figure 8). This approximately 48% increase indicated a critical transition in urban carbon dynamics. It suggested that as Chinese cities matured, the cumulative effect of multiple urban cores became more influential than that of a single large patch. This evolution mirrored China’s shift from monocentric to polycentric urban structures, wherein the aggregate carbon impact of distributed urban centers increasingly outperformed traditional core–periphery patterns. The categorical analysis found a fundamental shift in emission drivers, with compactness and concentration indicators increasing their collective contribution from 58.13% in 2005 to a peak of 69.54% in 2015, before moderating to 65.73% in 2020 (Figure 7; Supplementary Figures S3–S5). This trajectory indicates that urban compactness has become the overwhelmingly dominant morphological characteristic influencing emissions, while connectivity and integration indicators showed a corresponding decline from 20.47% in 2005 to 16.00% in 2015, suggesting that as cities densified, the role of inter-patch connectivity diminished relative to concentration effects.

GeoShapley spatial interaction analysis indicated the emergence of location-dependent morphological effects that fundamentally challenged universal urban planning approaches. GEO exhibited the widest value distribution (−10.0 to 7.5 × 10⁶), confirming that geographical context acted as a meta-variable moderating all morphological relationships. This result implied that identical urban forms yielded markedly different carbon outcomes depending on their geographical location. Interaction terms further demonstrated the evolution toward spatially contingent urban carbon systems. The coefficient for LPI × GEO was predominantly positive and increased from approximately 2.5 × 10⁶ in 2005 to 5.0 × 10⁶ by 2020 (Figure 9; Supplementary Figure S6). This doubling effect suggested that large urban patches increasingly incurred location-specific carbon penalties, likely reflecting regional variations in energy systems, industrial structures, and transportation networks as Chinese cities expanded. Finally, the emergence of 3D morphological interactions during 2015–2020 marked a paradigm shift in urban carbon determinants. Both VDAI × GEO and ECVR × GEO acquired measurable impacts (around 1.0 × 10⁶) that had been absent in earlier periods. This finding indicated that vertical urban development began producing regionally differentiated carbon impacts. It suggested that future urban carbon management must integrate vertical development patterns and geographical context—moving beyond traditional 2D planning frameworks.

The threshold analysis discovered distinct patterns in the relationship between urban morphological indices and CO₂ emissions over the study period (Figures 10, 11; Supplementary Figures S7, S8). Among the key findings, several morphological indices demonstrated evolving relationships with carbon emissions across the examined years. LPI exhibited positive GeoShapley values when exceeding thresholds of approximately 15. And TCA consistently showed positive contributions when exceeding thresholds ranging from 100,000 in 2005–2010 to 50,000 in 2015–2020, suggesting that extensive core urban areas were associated with higher emission levels. SDI showed monotonic increases with positive contributions above thresholds of 0.6 in 2005, progressively lowering to around 0.45 in 2015 and 2020. COH exhibited positive contributions with thresholds of around 97. UAI similarly maintained positive relationships when exceeding thresholds around 88–92 across all years. The temporal dynamics of 3D morphological indicators revealed particularly complex evolutionary patterns. VGKI underwent a notable transformation, showing positive contributions when values fell below 4 in 2005 and below 6 in 2010. VDAI demonstrated non-linear relationships that shifted dramatically over time, with positive contributions observed in the range of −0.03 to 0.14 in 2005, expanding to 0.0 to 0.2 in 2010, then shifting to −0.02 to 0.23 in 2015, before transitioning to a monotonic decrease pattern by 2020.

The analysis also revealed patterns of threshold evolution and emergence of new factors over time. PDI showed positive contributions below thresholds of 1.7 in 2010, within the range 1.1–1.9 in 2015, and below 1.4 in 2020, indicating that moderate patch density levels were associated with higher emissions. ECVR showed monotonic decrease patterns with positive contributions below 13 in 2010, dropping to below 4 in 2015 and below 2 in 2020, highlighting the growing importance but diminishing threshold of edge effects in urban carbon dynamics. ADA maintained relatively stable positive contribution ranges across all years, typically between −0.02 and 0.06, suggesting consistent impacts of directional horizonal urban development patterns on emissions. Most significantly, the analysis revealed that CBC consistently showed positive contributions above thresholds ranging from 213 in 2005 to values around 150–160 in later years, while VGDR transitioned from positive contributions above −0.3 in 2010 to showing no clear patterns in subsequent years. The general trend toward lower thresholds and narrower positive contribution ranges for many indices over time potentially indicated improved urban planning practices and more efficient urban forms that had successfully reduced fragmentation-related emissions by 2020.

Furthermore, bootstrap resampling with 1,000 iterations quantified threshold uncertainty and temporal stability of morphological metrics from 2005 to 2020. Coefficient of variation (CV) and 95% confidence intervals were computed, with stability classified based on cross-year CV and confidence interval (CI) overlap. LPI, COH, and UAI exhibited high temporal stability (CV < 0.10, CI overlap > 60%). LPI maintained a stable threshold around 15 (CV = 0.067), with narrow confidence intervals, confirming its reliability as a threshold indicator. COH and UAI showed stable thresholds (CV = 0.084 and 0.093, respectively). Whilst VDAI and ECVR displayed higher variability, indicating context-dependent thresholds. Of 15 metrics, 7 achieved high to very high stability, with CV < 0.10, 4 showed moderate stability, with CV ranges from 0.10 to 0.15), and 4 had low to very low stability, with CV > 0.15.

This study discovers a fundamental duality in how urban morphological characteristics influence carbon emissions, as observed through a triangulation of spatial regression, threshold modeling, and GeoShapley analysis. While spatial regression models consistently identify UAI as the strongest positive predictor of emissions, with elasticity ranging from 0.754 to 0.781. GeoShapley analysis presents a contrasting picture, revealing predominantly negative UAI contributions (−0.312 × 10⁶ to −0.316 × 10⁶) across spatial contexts. Rather than contradiction, this divergence reflects the differing lenses of analysis: spatial regression captures system-wide average relationships, indicating that more aggregated cities emit more carbon in total due to economic clustering and intensified activity, whereas GeoShapley highlights local efficiency patterns, showing that aggregation often reduces per capita emissions by enabling infrastructure sharing, reduced travel distances, and spatial compactness.

This apparent paradox extends to other compactness indicators such as LPI, TCA, and CBC. GeoShapley reveals negative associations between these metrics and emissions, yet threshold analysis shows that once values exceed critical limits, such as LPI > 15 or TCA > 50,000, these same traits begin to increase emissions. This nonlinear tipping behavior suggests that cities initially benefit from agglomeration economies, but beyond certain morphological thresholds, diseconomies such as congestion, longer intra-urban trips, and concentrated industrial outputs offset early gains. Over time, these thresholds have shifted, TCA’s critical value declined from 100,000 to 50,000 between 2005 and 2020, which indicates adaptive morphological recalibrations under policy and land-use constraints.

Our finding that 3D morphological indicators emerged as significant carbon emission predictors by 2015–2020 represents a critical advancement beyond traditional 2D-focused studies. While earlier research by Fang et al. analyzing 30 Chinese provincial capitals relied exclusively on 2D morphological variables achieving R² values of 0.13–0.18 (Fang et al., 2015), our GGPR approach incorporating 3D morphological variables achieved R² values consistently above 0.99. This surpasses even recent studies like Tian et al. examining Wuhan, which introduced advanced 3D indicators including sky view factor and street topology but achieved lower explanatory power using random forest models (Tian et al., 2024). Particularly critical is our discovery that UGSI’s emission reduction impact tripled from −0.096 to −0.297, demonstrating that synchronized vertical-horizontal development has become the most powerful morphological mitigation factor, a finding that extends beyond Wang et al.’s analysis of 53 Chinese cities, which found building height increased energy consumption while building volume reduced it, but did not capture the synergistic effects of synchronized development (Wang et al., 2024).

The transition we document from 2D to 3D dominance aligns with but extends beyond recent methodological advances in the field. Hong et al.’s literature review noted that although most past research has focused on 2D metrics to describe urban layout and density, 3D measures such as volumetric density and building elevation are gaining importance in shaping policy frameworks (Hong et al., 2022b). Our temporal analysis reveals the dynamic nature of this transition: VDAI shifted from negative to positive impacts over the study period, marking a critical threshold where asymmetric vertical development transformed from an efficiency strategy to an emission liability. This temporal dimension was overlooked in Wang et al.’s commercial building analysis in Beijing, which established correlations between building shape factor and carbon emissions but remained static in approach (Wang et al., 2023).

Our GeoShapley analysis advances the field by revealing location-specific variations that reconcile conflicting findings in previous research. Our findings demonstrate that the same 3D indicators can reduce per capita emissions in developing western cities while increasing them in saturated eastern megacities. The emergence of significant 3D interaction terms (VDAI × GEO, ECVR × GEO) by 2015–2020 confirms Tian et al.’s assertion that the configuration of the built environment plays a more substantial role in shaping emissions than classic determinants such as land use or population density (Tian et al., 2024).

Our findings suggest several transformative policy interventions that could reshape urban carbon management in China. Cities are suggested to implement performance-based planning frameworks that establish carbon intensity benchmarks linked to morphological thresholds. When LPI approaches critical values (15 for eastern cities, lower for western regions), planning authorities could require mandatory mixed-use development and green corridor integration to prevent over-concentration. These thresholds would serve as performance indicators rather than rigid boundaries, allowing flexibility in how cities achieve carbon reduction goals while respecting local development contexts. Furthermore, the strong negative impact of UGSI suggests cities should adopt synchronized height-density frameworks through performance-based floor-area-ratio systems. This approach would link vertical development rights to demonstrated infrastructure capacity and carbon performance metrics. New urban districts would begin with baseline height allowances that increase based on achieved performance in transportation connectivity, utility provision, and measured carbon intensity reductions. Areas experiencing rapid horizontal expansion would face height limitations until meeting specific infrastructure and carbon performance benchmarks. This performance-oriented approach ensures that vertical development proceeds in tandem with horizontal infrastructure readiness without requiring complex real-time adjustment mechanisms.

The regional heterogeneity revealed by our analysis demands differentiated performance targets rather than uniform morphological standards. Eastern megacities approaching saturation could establish carbon intensity reduction targets achievable through polycentric restructuring and vertical optimization strategies. Central cities in transition could set moderate performance goals attainable through selective densification combined with green infrastructure integration. Western cities maintaining development potential could adopt absolute emission caps while continuing compact development, with regular performance reviews to identify when diminishing returns emerge. These insights contribute to a broader theoretical shift toward adaptive urban carbon science. The first key contribution is the recognition of scale-dependent optimization: while system-level aggregation may increase total emissions, it can simultaneously enhance local carbon efficiency, requiring policies that balance absolute and per capita targets. Second, the analysis highlights threshold-mediated nonlinearity, where seemingly linear relationships break down once morphological indicators cross critical values, prompting the need for graduated response strategies. Finally, the study establishes the principle of spatially intelligent morphology, showing that identical urban forms can yield divergent carbon outcomes depending on their geographic and infrastructural context. Collectively, these findings challenge the prevailing logic of uniform planning standards and advocate for a new paradigm in which urban form is governed not by fixed metrics or automatic mechanisms, but by performance-based targets, spatial contingencies, and development-stage-specific interventions that remain implementable within existing governance frameworks.

Our study offers significant insights into urban morphology-emission dynamics in 286 Chinese cities from 2005 to 2020, yet its scope is bounded by several limitations that shape its generalizability. The reliance on city-level data, while robust for broad trends, may mask intra-urban heterogeneity, as Luqman et al. demonstrated that urban CO₂ emissions patterns differ significantly according to development level across their global city analysis (Luqman et al., 2023). Potential biases in our data sources, such as inconsistencies in satellite-derived morphological metrics or incomplete emissions inventories, could further skew results, a concern echoed by Privitera et al. in their global review of urban climate studies (Privitera et al., 2018). Moreover, the specific thresholds identified (e.g., for LPI or TCA) are tied to China’s unique urbanization trajectory, marked by rapid growth and state-driven planning, which may not translate to regions with slower development paces or decentralized governance.

Despite these constraints, the core mechanisms we discover, namely, the non-linear relationships, the rise of 3D morphology, and evolving thresholds, offer transferable value. The tension between compactness and emission efficiency could inform urban strategies in other developing nations, as Li et al. suggest in their cross-sectional analysis of urban form and socioeconomic drivers (Li et al., 2022). Similarly, the increasing relevance of vertical development might resonate in global cities embracing high-rise typologies, a trend Frolking et al. quantified by showing a profound global shift from lateral urban expansion to vertical development across around 1,550 cities, with growth rates in building volume increasing in essentially all regions even as horizontal expansion rates declined (Frolking et al., 2024). Generalizability hinges on contextual factors such as economic development, governance structures, and cultural attitudes toward sustainability. For instance, our findings may apply less readily to cities with robust green infrastructure or stringent environmental policies, where morphology-emission links might be mediated differently. Future research could address these limitations by testing our framework across diverse regions, incorporating micro-scale data, and expanding the suite of morphological indicators to include transportation networks or building efficiency, thereby refining and broadening the applicability of our conclusions for global carbon mitigation efforts.