Intercity connectivity enhances urban mobility resilience to extreme rainfall

Data sources and preprocessingMobility data

Two Baidu location-based service datasets were used. As of 2023, these data were derived from more than 120 billion daily location-based service requests collected from over 1.1 billion mobile devices across China, covering approximately 77.5% of the national population and showing a gender composition close to official census statistics47,48. The first dataset, from Baidu Migration (2020–2024), provided daily city-level intracity mobility intensity, intercity inflow/outflow indices and origin–destination flow proportions. Intercity indices were validated against Tencent Migration data to confirm linear scalability (Supplementary Fig. 10), consistent with previous studies49. The second dataset (October 2013 to May 2014) consisted of a daily population flow matrix among 2,862 counties in China50 and was aggregated to the city level. To match Baidu Migration data, which reports only the top 100 destinations per origin, we retained for each city-day the top 100 intercity connections by flow volume. Spatially, the analysis covered 366 mainland Chinese cities, excluding Sansha, Hong Kong, Macau and Taiwan because of data incompatibility or low population.

Meteorological and covariate data

Historical daily precipitation, maximum temperature and wind data (1981–2023) were obtained from ERA5-Land51 (0.1° resolution). Future daily precipitation projections (2025–2050) under the SSP245 and SSP585 scenarios were obtained from an ensemble of Coupled Model Intercomparison Project Phase 6 (CMIP6) climate models52 (Supplementary Table 11). All data were averaged to the city level in Google Earth Engine using boundaries consistent with the mobility data in 2023. ERA5 data (1981–2020) were used to bias-correct future projections via quantile mapping53.

To isolate connectivity effects, we compiled city-level covariates spanning geographical, land use, socioeconomic and governance dimensions (Supplementary Table 5 and Supplementary Method 1). Future socioeconomic projections under SSP245 and SSP585 were obtained from public datasets, including total population54, ageing ratio55, urbanization ratio55, GDP56 and the shares of industrial and service-sector output57.

Quantifying mobility resilience

We identified city-level extreme rainfall events based on daily precipitation. An extreme rainfall day is defined as any day with precipitation exceeding 50 mm or surpassing a city-specific 3-year return level. Return levels were estimated via a peak-over-threshold generalized Pareto distribution using 1981–2020 data58. Spatially and temporally adjacent extreme rainfall days were grouped into a single event. Each event affecting a city was treated as a city-event exposure, which served as the unit of analysis.

For a given city-event, the first extreme rainfall day was defined as Day 0. Baseline intracity mobility intensity was calculated as the average mobility level on non-holiday days with minimal rainfall (≤2 mm d−1) during the 14 days preceding Day 0, with separate weekday and weekend baselines. City-event exposures preceded by another extreme rainfall within 14 days were excluded to avoid overlapping impacts.

For each city-event exposure (where i denotes the city and j denotes the event), we quantified four metrics: maximum impact (MIi,j, percentage change), recovery time (RTi,j,days to return to baseline), recovery shape (Si,j, integrated curvature) and total performance loss (TPLi,j, cumulative deviation treated as a dimensionless metric, with 1 equivalent to 100%) (Extended Data Fig. 1). We examined the process-based components of TPL using a log-linear decomposition:

$$\mathrm{ln}\left({{\rm{TPL}}}_{i,\,j}\right)={\beta }_{0}+{\beta }_{1}\mathrm{ln}\left({{\rm{MI}}}_{i,\,j}\right)+{\beta }_{2}\mathrm{ln}\left({{\rm{RT}}}_{i,\,j}\right)+{\beta }_{3}\mathrm{ln}\left({S}_{i,j}\right)+{\varepsilon }_{i,\,j}$$

(1)

Mobility change curves and descriptive resilience metrics for 2023 are shown in the main figures to represent non-pandemic conditions.

Impact analysis

We calculated 15 intercity connectivity indicators from annually averaged daily mobility networks (2022–2023), spanning interaction activity, network centrality, connectivity preference and interaction evenness (Extended Data Table 1). All indicators were relative measures (proportions or normalized indices) to avoid scale-dependent effects.

To estimate conditional associations between intercity connectivity and mobility resilience, we used a multitreatment DML framework with fixed effects in four steps.

First, we residualized all resilience metrics (Y), connectivity treatments (D) and covariates (X) by partialling out city and year fixed effects, isolating within-city temporal variation. All resilience metrics were log-transformed, allowing the coefficients to be interpreted as proportional changes in resilience associated with relative shifts in connectivity.

Second, we assessed residual confounding using high-dimensional diagnostics. Using residualized data, Elastic Net and Random Forest models showed that controls explained substantial variation in TPL (cross-validated R2 = 0.38–0.41) but limited connectivity variation (Elastic Net R2 < 0.06, Random Forest R2 < 0.43; Supplementary Table 12), underscoring the necessity of rigorous confounding adjustment. When using only hazard and governance indicators, explanatory power for TPL remained robust (R2 = 0.39–0.40) whereas predictability of connectivity dropped precipitously (R2 < 0.05 for Elastic Net and R2 < 0.13 for Random Forest), mitigating reverse causality concerns.

Third, we applied DML orthogonalization to obtain de-biased estimates of connectivity effects under high-dimensional controls34. The formula is:

$$\widetilde{{Y}_{{ij}k}}=g\left(\widetilde{{X}_{{ij}k}}\right)+{u}_{{ij}k}$$

(2)

$$\widetilde{{D}_{ike}}={m}_{e}\left(\widetilde{{X}_{ikt}}\right)+{v}_{ike}$$

(3)

where \(\widetilde{{Y}_{{ij}k}}\), \(\widetilde{{X}_{{ij}k}}\) and \(\widetilde{{D}_{ike}}\) denote variables in year k after removing city and year fixed effects. The nuisance functions \(g(\bullet )\) and \({m}_{e}(\bullet )\) were estimated using Elastic Net regularization to stabilize estimation under multicollinearity. The residuals \({u}_{{ij}k}\) and \({v}_{ike}\) were subsequently used to identify the orthogonalized association between intercity connectivity and mobility resilience.

Fourth, because connectivity indicators were correlated, we grouped indicators by correlation structure (Supplementary Fig. 5) and estimated multitreatment DML models. From each correlation group, we selected one statistically significant indicator in the single-treatment DML analysis and jointly included the selected indicators from different groups in the multitreatment DML framework:

$${u}_{{ij}k}=\sum _{e=1}^{E}{\theta }_{e}{v}_{{ike}}+{\varepsilon }_{{ij}k}$$

(4)

where \({\theta }_{e}\) represents the estimated conditional association of connectivity indicator e, holding other indicators constant.

To address sample-splitting randomness, all DML estimations were repeated across 100 random splits. Following stability-selection logic59, we classified indicators as high-stability significant if >90% of iterations yielded P < 0.1, low-stability significant if 35–90% and otherwise insignificant. The 35% lower bound acts as a heuristic noise filter, conceptually analogous to tail probabilities outside 1 s.d. Sensitivity analyses with alternative thresholds (85%/30%) did not alter conclusions. Multisplit inference with aggregated P values confirmed this classification34 (Supplementary Tables 4 and 7).

Robustness was evaluated through: (1) single-treatment and multitreatment DML comparisons; (2) alternative indicator combinations; (3) replacing Elastic Net learners with random forest learners; and (4) placebo tests.

Sensitivity of resilience to intercity connectivity changes

We conducted counterfactual simulations to evaluate how deviations from the 2023 connectivity baseline altered TPL. Three primary counterfactual configurations were defined: (1) socioeconomic rollback (2013–2014) representing early-stage connectivity; (2) policy-constrained mobility (2021) representing a network under policy restrictions; and (3) socioeconomic evolution (2024) representing the most recent configuration under ongoing socioeconomic development. A strict-restriction configuration (21–29 February 2020 in the early phase of the COVID-19 outbreak) was used for sensitivity analyses only.

Connectivity indicators derived from these diverse historical and policy-driven configurations \({D}_{{ie}}^{\rm{cf}}\) were combined with coefficients estimated from the multitreatment DML models. These simulations were intended as what-if analyses of the sensitivity and directional response of resilience to collective changes in intercity connectivity, rather than as forecasts.

For each city i during a rainfall event j, the predicted relative TPL change (\(\Delta {\rm{TPL}}_{{ij}}\)) was:

$$\Delta {{\rm{TPL}}}_{{ij}}=\exp \left(\mathop{\sum }\limits_{e=1}^{E}{\theta }_{e}\left({D}_{{ie}}^{\rm{cf}}-{D}_{{ie}}^{2023}\right)\right)-1$$

(5)

where \({\theta }_{e}\) represents the coefficient of the e connectivity indicator derived from 100 DML iterations and \({D}_{{ie}}^{\rm{cf}}\) and \({D}_{{ie}}^{2023}\) denote the counterfactual and 2023 baseline configurations, respectively. Using observed event-specific TPL (\({{\rm{TPL}}}_{{ij}}^{{\rm{obs}}}\)) as baseline, predicted counterfactual TPL was:

$${\mathrm{TPL}}_{{ij}}^{\mathrm{cf}}=\,{\mathrm{TPL}}_{{ij}}^{\mathrm{obs}}\times \,(1\,+\,\Delta {\mathrm{TPL}}_{{ij}})$$

(6)

where \({{\rm{TPL}}}_{ij}^{{\rm{obs}}}\) is the observed TPL for city i during event j under the 2023 connectivity configuration.

To assess the systemic impact, we calculated the national exposure-weighted relative changes in TPL as:

$$\Delta {\mathrm{TPL}}_{\mathrm{national}}=\frac{{{\sum }}_{i,j}\,\left({\mathrm{TPL}}_{ij}^{\rm{cf}}-{\mathrm{TPL}}_{ij}^{\mathrm{obs}}\right)}{{{\sum }}_{i,j}{\mathrm{TPL}}_{ij}^{\mathrm{obs}}}$$

(7)

Fewer than 9% of counterfactual values fell outside the empirical range (Supplementary Table 13); winsorization-based sensitivity checks confirmed that results were not driven by extrapolation bias (Supplementary Table 14).

Future climate change scenario analysis

On the basis of projected changes in extreme rainfall under future SSP–RCP (representative concentration pathways) scenarios, we examined trends in (1) extreme rainfall city-days, (2) extreme rainfall events and (3) the share of city-days occurring within multicity events. All projections used CMIP6 multimodel ensemble means.

Future evolution of intercity connectivity

To quantify how future socioeconomic development translates into structural changes in intercity connectivity, we used a multivariate Random Forest regression model60 to model the joint evolution of the identified indicators, trained on city-level panel data (2013 and 2023) with fivefold cross-validation. Permutation-based variable importance identified railway connectivity as a dominant driver (Supplementary Fig. 7).

On the basis of these results, we constructed three future connectivity pathways under common SSP-based socioeconomic and population projections, differing only in railway development assumptions: (1) SSP-based pathway, with railway connectivity fixed at 2023 levels; (2) railway expansion pathway, assuming annual railway growth of 2% from 2023 and nationwide connectivity by 2035; and (3) railway equity pathway, prioritizing cities with weaker initial railway infrastructure, with annual railway growth of 5% until reaching the 2023 national average and 2% thereafter. Pathways (1) and (3) were constructed on the basis of the existing national railway development plan of China37, which targets nationwide railway connectivity by 2035.

Projected impacts on mobility resilience

Using projected connectivity trajectories, we quantified city-level relative TPL changes under extreme rainfall events. To account for uncertainty, we combined 1,000 bootstrap resamples with Monte Carlo simulations, propagating uncertainty from DML estimates and projected connectivity evolution. Predicted connectivity indicators were recalibrated using linear correction functions derived from observed and predicted values in 2023 (Supplementary Fig. 11), ensuring zero baseline deviation. Results were reported as median estimates across all cities with 95% confidence intervals (CIs).

Economic impact assessment

To translate projected TPL changes into mobility-related indirect economic losses, we first estimated a future baseline TPL assuming 2023 connectivity. We used the best-performing machine learning model (Random Forest, XGBoost or Elastic Net) trained on 2022–2023 data. The training set covered cities in nearly all provinces and includes historical events with >100-year return periods, supporting spatial-temporal generalizability under a space-for-time substitution assumption. To minimize systematic bias across different urban scales, projected baseline TPL values were calibrated using population-group-specific correction coefficients derived from the relationship between predicted and observed values in 2023 (Supplementary Fig. 12).

Projected TPL changes were translated into indirect economic outcomes by accounting for heterogeneous mobility dependence across economic activities. We disaggregated the secondary and tertiary sectors into 18 subsectors, weighting the sectoral GDP share by employment and national average wages for each city:

$${{{S}}}_{{{i}},{{j}},{{k}}}=\frac{{{{L}}}_{{{i}},{{k}}} {{{W}}}_{{{k}}}}{{\sum }_{{{k}}\in {{j}}}{{{L}}}_{{{i}},{{k}}} {{{W}}}_{{{k}}}}$$

(8)

where \({S}_{i,j,k}\) is the GDP share subsector k of sector j in city i, Li,k is the number of employed persons in subsector k in city i; Wk is the national average wage of subsector k. Sectors were further classified into four categories of mobility dependence, with sector-specific elasticity ranges linking mobility loss to output loss. Two elasticity regimes were considered: a normal case and a technology-enhanced case reflecting greater penetration of remote work and digital services (Supplementary Information and Supplementary Table 10). To address uncertainty, the main analysis randomly sampled between regimes. Two sensitivity scenarios were further examined: (1) static technology level (elasticities fixed at 2023 normal levels) and (2) accelerated technology growth (quadratic increase from normal to maximum technology-enhanced regime by 2050).

To quantify avoided economic losses and account for cascading uncertainties, 1,000 Monte Carlo iterations propagated uncertainties from climate projections, connectivity pathways, TPL predictions and elasticities. All reported economic impacts were expressed in constant 2023 GDP at purchasing power parity (PPP). Given the assumption of stable future relationships and several sources of uncertainty, estimates should be interpreted as order-of-magnitude and directional indications of mobility-related economic impacts, rather than as precise actuarial valuations.

Limitations

First, while human mobility resilience serves as a powerful proxy for functional urban recovery, it represents a necessary but not sufficient dimension of urban resilience. Mobility dynamics often exhibit strong synergies with economic resilience, particularly in service-oriented sectors where labour supply, consumption and daily economic activity depend directly on population movement29 and they can closely track the restoration of critical transportation infrastructure4. However, mobility-based metrics do not directly capture more latent dimensions of resilience, such as social capital or emergency management capacity32,33,61, which are central to long-term adaptive governance. An exclusive focus on mobility recovery may obscure important trade-offs. For example, rapid restoration of mobility supported by external support does not necessarily strengthen the long-term climate adaptation capacity of a city, leading to a ‘return to vulnerability’ rather than transformative adaptation62. More broadly, in high-vulnerability regions, recovery to predisaster levels alone is insufficient to constitute resilience, which instead requires structural shifts that enhance future coping and adaptive capacities. Accordingly, our findings should be interpreted as evidence of how intercity connectivity affects disaster-related mobility resilience, rather than a comprehensive assessment of urban resilience or long-term adaptation capacity.

Second, although our analysis uses panel fixed effects and DML to control for high-dimensional observable confounders, causally interpretable estimates remain challenging. On the one hand, unobserved confounding cannot be fully eliminated in observational studies and residual time-varying unmeasured factors may remain. For example, complex, difficult-to-quantify interventions that evolve over time, such as targeted economic development initiatives or sponge city policies63, might influence intercity connectivity and postdisaster recovery capacity, potentially biasing the results. Therefore, despite the use of causal methods, our findings should be interpreted as robust conditional associations rather than definitive causal effects. On the other hand, our estimation may also be conservative because stringent control for confounders may partial out shared effects of connectivity. For example, elite-city linkages probably encompass not only frequent demographic and economic interactions but also dense information exchange, which helps cities to better forecast disasters and enhance their resistance capabilities. However, by strictly controlling for high-dimensional confounders, our methodology may have partialled out these shared effects, isolating only the unique, orthogonal variation of intercity connectivity. Accordingly, statistical insignificance in some estimates does not preclude underlying mechanisms, although the findings should still be interpreted with caution.

Third, projections of future intercity connectivity and resilience rely on extrapolating historically observed structural relationships, which may be altered by disruptive technological changes, including digital connectivity, remote interaction and emerging transport technologies. For example, advances in electric vehicle technologies, together with continued improvements in road infrastructure, may increasingly facilitate intercity interactions and, in some contexts, partially substitute for rail-based connectivity. These dynamics, however, extend beyond the scope of this study. We therefore limit our forward-looking analysis to midcentury and interpret the results as directional insights rather than precise forecasts.

Fourth, our future scenario analysis considers extreme rainfall in isolation and does not explicitly evaluate compound hazards or more systemic cascading disruptions across connected cities. Although compound hazards can exacerbate TPL23 (Supplementary Fig. 8) and our empirical analyses include events affecting some central hubs, the observed rainfall events did not involve widespread simultaneous disruption of major hubs that could have caused observable cascading failures. This may lead to an underestimation of future losses.

Fifth, our economic assessment is partial rather than a full cost–benefit evaluation of railway construction. It excludes direct infrastructure damage, cascading supply-chain disruptions, economic and environmental costs associated with railway construction itself and any health-related costs. Our estimates should therefore be interpreted as directional and order-of-magnitude indications of the additional economic implications associated with resilience gains under an existing railway planning pathway.

Finally, the data used in this study have several limitations. With respect to mobility data, smartphone penetration rates in China vary across years, increasing from approximately 46.9% in 2013 to 78.6% in 202448,50. Although we applied consistent filtering and used relative indices, lower coverage in 2013 may have under-represented low-income and elderly populations with lower mobility, potentially biasing baseline connectivity upward. Because connectivity change is measured as the 2023 and 2013 difference, such baseline inflation would attenuate the estimated increase, rendering our estimated resilience gains conservative. In addition, the user population does not fully represent all demographic groups, such as children, and may involve potential biases across age groups and urban–rural contexts48. To protect user privacy, mobility data are aggregated at the daily and city levels and therefore do not capture traveller type, trip purpose or trip duration47. As a result, short-term inflow surges may reflect external support-related movements but cannot be treated as a direct proxy for external aid. Regarding railway data, the OpenStreetMap-based railway network in earlier years, particularly around 2013, may be incomplete, introducing uncertainty into historical railway connectivity indicators. Potential longer-term redistribution dynamics of railway development, such as positive spillover or negative drain effects44,45, are not explicitly incorporated into our scenario analysis, as their magnitude, direction and timescale remain debated and empirically heterogeneous. Collectively, these limitations may introduce uncertainty into our estimates and should be considered when interpreting the results.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.