Computable Resilience Model for Urban Safety

Abstract

The governance of urban safety is undergoing a critical paradigm shift from static threshold defense to dynamic system resilience. However, a general full-process mathematical framework integrating shock entry, degradation, recovery, and post-shock adaptation remains limited. To bridge this theory-practice gap, this study proposes a methodological framework for a Computable Resilience Model (CRM). Moving beyond static indicator aggregation, the CRM is characterized by four core computable features: (1) macroscopic representation, (2) dynamic time-dependence, (3) scenario simulatability, and (4) multidimensional generalizability. By translating abstract system capabilities into independent dynamical control parameters, the framework decouples the shock-induced performance drop from the post-shock recovery, thereby quantifying the entire non-equilibrium trajectory of complex urban systems. Furthermore, by establishing a common temporal benchmark, the CRM converts cumulative resilience loss into a unified time-integrated metric for relative comparison across assessment objects. Ultimately, this study advances the computable representation of urban resilience, providing a practical mathematical foundation for complex risk governance.

Graphical abstract

Keywords

Computable resilience / Urban safety / Dynamical modeling / Extreme shocks / Scenario simulatability

Highlight

	● CRM maps the 4R attributes to independently interpretable dynamical parameters.
	● State-transition equations capture shock entry, degradation, recovery and adaptation.
	● Shock inputs are decoupled from hazard type to support multi-hazard applications.
	● A common time benchmark and signed loss enable consistent resilience comparison.
	● Resilience declines with rainfall intensity but varies markedly across cities.

1 Introduction

Resilience is an interdisciplinary concept that fundamentally describes a system’s dynamic process of absorbing shocks, maintaining core functions, and self-organizing to recover following a disturbance [1]. When applied to urban safety, the concept of resilience shifts from a natural ecological metaphor to an engineering concern. For cities, as complex giant systems with highly coupled physical, social, and economic networks, resilience to external shocks is not merely an intangible conceptual property. Rather, it manifests as an objective physical phenomenon through measurable parameters like infrastructure capacity, resource allocation efficiency, and socio-economic recovery trajectories. However, a significant methodological gap persists in quantifying resilience for specific urban contexts. The core challenge lies in transforming multidimensional, qualitative resilience attributes into computable and horizontally comparable mathematical variables that can directly inform engineering practice. This translation from metaphor to quantification remains a fundamental bottleneck in urban resilience research [2,3].

Driven by globalization and intensive urbanization, modern cities have evolved into complex giant systems characterized by deep interdependence among their physical spaces, infrastructure, and socioeconomic networks. While this interdependence enhances efficiency, it also fundamentally reshapes how urban risks are generated and propagated [4]. Today, external shocks to urban systems are increasingly non-stationary and cross-domain in their propagation. Their impacts are shifting from single-dimensional damage to multi-faceted systemic crises [5]. Consequently, localized damage to critical nodes can easily cascade beyond individual facilities, triggering unpredictable failures or even widespread functional collapse [6,7]. In contrast, traditional urban safety governance often relies on a static, defense-oriented paradigm based on historical data. This rigid approach faces an adaptability crisis when confronting highly uncertain, compounding systemic shocks [8]. Therefore, it is imperative to move beyond static, threshold-based defense and establish a dynamic resilience paradigm. Such a paradigm must be capable of quantitatively analyzing the entire process in complex urban systems under stress, from performance attenuation and shock absorption to systemic adaptation. This represents a core theoretical and practical priority for modern disaster risk reduction and urban safety.

As an alternative to traditional static disaster prevention models, resilience theory provides a conceptual foundation for evaluating the dynamic capacity of complex urban systems to cope with extreme shocks. However, at the macro scale, this dynamic concept has not been effectively translated into computable mathematical models. Mainstream evaluation paradigms still rely heavily on static indicator systems, such as the widely applied Disaster Resilience of Place (DROP) model [9]. As highlighted in the review by Hosseini et al., such approaches often rely on methods like the analytic hierarchy process, employing dimensionality reduction and linear weighted aggregation of multidimensional proxy variables [10]. While achieving broad coverage of urban subsystems, these methods exhibit systemic deficiencies in characterizing underlying physical mechanisms. Linkov et al. further point out that this linear superposition logic, based on cross-sectional data, conflates static preparedness resources with dynamic system performance [11]. It eliminates the temporal dimension, offering only a static snapshot and failing to capture the non-linear, coupled dynamics of shock and recovery in complex systems.

To address the limitations of static assessments, another branch of research has turned to micro-evolutionary simulations, such as those using agent-based modeling (ABM) or complex network theory, to simulate the temporal trajectories of system collapse and recovery [12–14]. Although these bottom-up, physics-driven models offer high temporal resolution, they face significant bottlenecks in large-scale engineering applications. As noted by Zio, generating and executing such micro-simulation models requires fine-grained data, leading to exponentially increasing computational complexity for real-world, large-scale urban networks [15]. Moreover, for low-frequency, high-impact extreme shocks, the numerous micro-behavioral parameters in non-stationary environments are often associated with high uncertainty and are difficult to calibrate [16]. Consequently, high computational costs and data requirements severely limit the use of micro-dynamic simulations for cross-regional comparisons and national-level policy optimization.

To resolve the dual bottlenecks of static indicators (lacking dynamic processes) and micro-simulations (lacking macroscopic generalizability), this study proposes a methodological framework for computable resilience. This framework reconstructs the 4R attributes (Robustness, Redundancy, Rapidity, and Resourcefulness) of resilience theory [1], mapping them onto dynamical control parameters within the CRM [17]. On this basis, we develop a general computational model that encompasses standardized shock input, dynamic performance simulation, and spatiotemporal loss quantification. The framework can accommodate different types of urban shocks and provide a computable basis for resilience-oriented engineering planning.

The principal methodological innovation of the CRM is the establishment of a general computable framework that represents the full resilience process from external shock, risk entry, and performance degradation to recovery and post-shock adaptation. The framework integrates these stages within a continuous macroscopic dynamical process and translates the 4R attributes into independently interpretable control mechanisms. Through this architecture, shock buffering, recovery progression, and long-term adaptive capacity can be represented separately while remaining dynamically connected within the same system trajectory. This provides a transferable basis for simulating heterogeneous shocks, comparing resilience responses, and evaluating intervention scenarios across different urban and engineering systems.

2 The General Methodological Framework of the CRM

The purpose of the CRM is to translate the abstract concept of urban resilience into an executable dynamical system. As conceptually illustrated in Fig. 1, the CRM represents an urban system as a macroscopic performance trajectory under external shock, rather than as a static composite score. The 4R attributes—robustness, redundancy, rapidity, and resourcefulness—are translated into distinct dynamical characteristics that respectively govern the maximum degradation amplitude, shock buffering during risk entry, recovery rate during risk exit, and the terminal post-shock steady-state level. Because the generalized shock input is decoupled from specific physical hazard forms, the same computational architecture can be applied to multiple types of urban shocks. The CRM therefore provides a unified framework for describing risk entry, performance degradation, risk exit, and transition toward a post-shock steady state.

2.1 Cross-domain shock measurement and parameter calibration

To support the multi-hazard applicability illustrated in Fig. 1, the input layer of the CRM is decoupled from specific physical hazard forms. For assessment object i under scenario s, the effective external shock X_{i,s} is characterized by hazard intensity H_{i,s} and the exposure of affected assets or system components E_{i,s}:

(1)X_{i,s}=G\left(H_{i,s},E_{i,s}\right)

Here, G(\cdot ) denotes a domain-specific shock characterization function. This design allows different types of urban shocks to be incorporated within the same computational architecture. For probabilistic hazards, occurrence probability or return period is used to define or weight the corresponding shock scenario. For example, extreme rainfall can be represented by rainfall intensity under a specified return period and urban exposure, while blackouts or cyber-physical attacks can be characterized by load gaps or compromised nodes and the corresponding infrastructure exposure.

Once the shock input is established, the CRM calibrates the dynamical parameters governing system evolution. The external shock primarily enters the model through the scenario-specific vulnerability amplitude a_{i,s}:

(2)a_{i,s}=F\left(X_{i,s},C_i^{\text{rob}}\right)

Here, C_i^{\text{rob}} denotes the robustness-related resistance capacity and F(\cdot ) represents a domain-specific shock–response relationship. Thus, a_{i,s} reflects the combined effect of external shock severity and system robustness. Departing from traditional static assessments that aggregate economic, social, and engineering indicators into a single score, the CRM constructs a proxy indicator system aligned with the 4R attributes. The corresponding 4R attributes are subsequently translated into the dynamical parameters a, b, c, and d through an application-specific calibration procedure. The CRM does not prescribe a universal set of proxy indicators or a single weighting algorithm, because the observable variables and their physical meanings depend on the hazard type, system boundary, and available data. Instead, each application should explicitly define the indicator system, normalization and weighting rules, and parameter-mapping procedure while preserving the dynamical meanings summarized in Table 1. We illustrate the CRM using an extreme-rainfall case study, with its detailed operationalization provided in the Supplementary Material.

2.2 Core equations and dynamical construction of the CRM

Existing quantitative resilience models commonly use piecewise functions or single-phase degradation assumptions to characterize system evolution. Classical recovery models, such as the exponential recovery model proposed by Cimellaro et al. [18], and the linear or step models proposed by Ayyub [19,20], typically divide system responses into two sequential phases: disruption and recovery. This formulation often assumes that degradation occurs instantaneously at the moment of disturbance [21]. However, in real urban systems, physical damage, cascading failure, emergency response, and resource allocation are frequently intertwined over time. As Eisenberg et al. emphasized, simple rebound curves may not adequately capture the overlapping and non-equilibrium nature of resilience processes in complex systems [22].

Continuous dynamical models provide a more suitable basis for representing such processes. For example, Todman et al. developed a critical-damping-based dynamical model [23], while Tang and Shen proposed a recovery framework based on the gamma function family [24]. These studies effectively overcame the non-differentiable mathematical singularities present in step models but left room for further expansion regarding the physical decoupling of dynamical parameters. In gamma-family formulations, the temporal parameters b and c jointly enter the same time-dependent term, so their effects on the temporal trajectory are mathematically coupled. This parameter coupling limits independent attribution of shock buffering and recovery progression. Consequently, this limits the model’s ability to separately project and attribute engineering interventions for extreme-shock prevention. The gamma-family formulation was subsequently applied to large-scale urban extreme-rainfall resilience assessment by linking the 4R attributes to the recovery-curve parameters [25]. In that application, the temporal parameters jointly shaped the recovery trajectory and d was fixed at 1, whereas the present CRM separates risk-entry and recovery-transition dynamics and allows resourcefulness to determine the terminal post-shock state.

Existing resilience formulations have advanced the dynamic representation of recovery, but their emphasis differs across model families. Piecewise and rebound-type curves often simplify the interaction between degradation and recovery stages, while some continuous recovery functions describe the overall trajectory through parameters that are jointly embedded in the same temporal term. These formulations can limit the independent interpretation of shock buffering, recovery progression, and long-term adaptation. The CRM extends this line of research by constructing a full-process state-transition framework in which risk entry, damage evolution, recovery transition, and post-shock adaptation are represented as connected but distinguishable dynamical processes. This structure enables the principal resilience attributes to regulate different stages of system evolution independently, thereby improving parameter interpretability, scenario controllability, and transferability across different shock and system contexts. The core equation is defined as:

(3)P(t)=P_{\text{damage}}(t)[1-R_c(t)]+d{R}_c(t)

or equivalently,

(4)P\left(t\right)=\left[1-a{S}_b\left(t\right)\right]\left[1-R_c\left(t\right)\right]+d{R}_c\left(t\right)

Here, P\left(t\right) denotes the normalized macroscopic system performance at time t. The first term, P_{\text{damage}}(t)[1-R_c(t)], captures the residual effect of the damage-only trajectory before recovery is completed. The second term, d{R}_c\left(t\right), drives the transition toward the terminal steady state. Therefore, the CRM treats resilience as a continuous state-transition process governed by the risk-entry-and-exit logic, rather than as a sequential “collapse–rebound” curve.

The damage-only trajectory is defined as:

(5)P_{\text{damage}}(t)=1-a{S}_b(t)

where a is the scenario-specific vulnerability amplitude and S_b(t) is the risk-entry function. More explicitly, a represents a_{i,s}, which is determined by the interaction between the external shock and the robustness of assessment object i under scenario s. For notational simplicity, the object and scenario subscripts are omitted in the core equations. A larger a indicates a deeper potential performance degradation under the specified shock scenario.

The risk-entry function is expressed as:

(6)S_b(t)=1-\exp (-t/b)

Here, S_b\left(t\right) represents the progress of shock-induced damage entering the system, and b is the shock-buffering time constant associated with redundancy. The parameter b has the same temporal dimension as t. A larger b corresponds to a longer characteristic timescale of risk entry, representing stronger buffering capacity that delays the translation of external shocks into functional degradation. Accordingly, increasing b slows the downward evolution of system performance without altering the eventual magnitude of the damage-only trajectory.

The risk-exit process is represented by a normalized recovery progress function:

(7)R_c(t)={\mathrm{\Phi }}_c{(t)}^2

Here, {\mathrm{\Phi }}_c(t) is a normalized recovery kernel satisfying {\mathrm{\Phi }}_c(0)=0 and \underset{t\to \mathrm{\infty }}{lim}{\mathrm{\Phi }}_c(t)=1. The parameter c is the recovery rate associated with rapidity and therefore has the dimension of inverse time. It controls the speed at which the system exits the degraded state and moves toward the terminal steady state. A larger c indicates faster recovery progress and a steeper upward curvature of the performance trajectory.

The square term in R_c(t) suppresses the initial recovery rate and prevents the recovery process from dominating the trajectory at the moment of shock. In other words, the system must first experience risk-entry and shock buffering before substantive recovery unfolds. This design preserves the chronological logic of risk-entry-and-exit: external shock first induces performance degradation through S_b\left(t\right), and the system then gradually transitions from the damage-only trajectory to the terminal steady state through R_c\left(t\right).

The terminal steady-state level is controlled by d, which is associated with resourcefulness. This parameter allows the CRM to relax the strict assumption that every system must return exactly to its pre-shock baseline. The equation satisfies two boundary conditions. At the initial moment, P\left(0\right)=1, representing the normalized pre-shock performance baseline. As time approaches infinity, \underset{t\to \mathrm{\infty }}{lim}P(t)=d.

This boundary condition is important because the final performance level is controlled only by resourcefulness parameter d, rather than being coupled with the vulnerability amplitude a. Thus, the CRM preserves the physical independence of robustness and resourcefulness. Accordingly, the CRM encompasses three possible terminal states:

1) When d=1, the system achieves elastic recovery [P(\mathrm{\infty })=P(0)];

2) When d>1, the system achieves a transformative leap and builds back better [P(\mathrm{\infty })>P(0)];

3) When , the system converges to a degraded steady state, indicating persistent or irreversible functional loss [].

The normalized recovery kernel can adopt different functional forms to accommodate different shock and recovery mechanisms, with three typical forms summarized in Table 2.

Taking the exponential recovery kernel as an example, Fig. 2 illustrates the sensitivity of the four parameters. The parameters b and c represent different temporal characteristics and therefore should not be directly compared numerically. Specifically, b is a shock-buffering time constant with the same temporal dimension as t, whereas c is a recovery rate with the dimension of inverse time. The chronological ordering between risk entry and recovery is instead embedded in the mathematical structure of the CRM. For the exponential kernel, as t\to 0^{+}, S_b(t)=t/b+O(t^2), whereas R_c(t)={[1-\exp (-ct)]}^2=c^2{t}^2+O(t^3). Thus, shock-induced degradation develops as a first-order process near the initial moment, while the recovery transition begins as a second-order process with R_c^{\mathrm{\prime }}(0)=0. This structure prevents recovery from dominating immediately after shock onset without imposing an additional numerical constraint between b and c. Generally, the proxy indicators for a, b, and c can be directionally normalized using min–max transformation and then mapped to their physical parameter ranges using rank-percentile mapping. For d, a centered Z-score transformation is recommended after aggregating the resourcefulness indicators, so that the reference mean corresponds to the neutral terminal state d = 1. The detailed calibration procedure adopted for the illustrative extreme-rainfall application is provided in the Supplementary Material.

2.3 Spatiotemporal quantitative output under a common evaluation horizon

After establishing the core equations and parameters, the CRM maps the simulated performance trajectory into spatiotemporal resilience indicators. A sufficiently long numerical simulation window T is first specified to capture the candidate shock–recovery trajectories. Within this computational window, the CRM identifies the recovery time of each assessment object and subsequently establishes a common integration benchmark for horizontal comparison. This procedure ensures that resilience is evaluated over a complete modeled shock–recovery process while all objects within the same comparison set share an identical temporal domain.

2.3.1 Dynamic recovery time calibration

When calculating recovery time (T_{\mathrm{r}\mathrm{e}\mathrm{c}}), traditional non-equilibrium dynamical models often use a fixed, extremely small rate-of-change threshold to determine steady-state truncation (e.g., {10}^{-6}). However, this rigid threshold may lead to a numerical truncation paradox in engineering calculations. Specifically, when encountering high-intensity shocks, the recovery curve is indefinitely prolonged due to its asymptotic characteristics, leading to distorted recovery times. Conversely, under minor perturbations, the system is prone to being prematurely judged as stable.

Therefore, this study adopts a dynamic relative change threshold method based on degradation depth. Assuming that the system reaches a minimum performance P_{\mathrm{m}\mathrm{i}\mathrm{n}} after being subjected to a shock force, the degradation depth is defined as \mathrm{\Delta }P=1-P_{\mathrm{m}\mathrm{i}\mathrm{n}}. The dynamic stabilization threshold is then defined as:

(8){ϵ}_{\mathrm{d}\mathrm{y}\mathrm{n}}=\gamma \cdot \mathrm{\Delta }P

where \gamma is a dimensionless relative convergence coefficient. It is recommended to set \gamma ={10}^{-3}. Recovery is considered stable when the absolute performance change between adjacent numerical time nodes falls below {ϵ}_{\mathrm{d}\mathrm{y}\mathrm{n}} for N consecutive time steps after the performance nadir. The dynamic recovery time T_{\mathrm{r}\mathrm{e}\mathrm{c}} is defined as the first time at which this stability criterion is satisfied, where \mathrm{\Delta }t denotes the numerical time step:

(9)

2.3.2 Data-adaptive common time benchmark for horizontal comparison

Because recovery duration can differ substantially among assessment objects, using individual T_{\mathrm{r}\mathrm{e}\mathrm{c}} values as integration boundaries would produce resilience measures over inconsistent temporal domains. Conversely, prescribing a fixed horizon before the recovery timescale is known may truncate slowly recovering systems before their shock–recovery process has stabilized. This is particularly relevant to the general CRM framework, for which the characteristic recovery timescale is not necessarily known before application.

Therefore, when no externally prescribed engineering or planning horizon is available, the CRM adopts a data-adaptive common benchmark defined by the longest recovery time within the comparison set:

(10)T_{max}=\underset{(i,s)\in S}{max}\{T_{\mathrm{r}\mathrm{e}\mathrm{c},i,s}^{}\}(s=1,2,\dots ,S)

Here, S denotes the complete set of all combined scenarios. This definition ensures that all assessed trajectories are integrated over the same temporal domain while allowing every case to reach its modeled recovery state. Importantly, T_{max} should be interpreted as a comparison-set-specific temporal benchmark rather than a universal physical constant. Changing the composition of the comparison set may change its numerical value and hence the absolute magnitude of the resulting resilience scores. Therefore, direct comparison of resilience values across independently defined samples requires the use of an identical temporal benchmark. Within a consistently defined comparison set, however, T_{max} provides a common basis for relative cross-case comparison. A leave-one-city-out sensitivity analysis further confirms that although T_{max} may vary slightly with sample composition, the relative resilience rankings remain stable across the three rainfall scenarios (Fig. S1).

Where an externally specified engineering horizon is available, such as a regulatory recovery target or an operational planning cycle, this external horizon can instead be adopted as the common integration boundary. Likewise, for infrastructure systems with sufficiently complete observed shock–recovery trajectories, the CRM parameters can be calibrated through inverse fitting to observed system performance, following the general trajectory-fitting strategy used in recovery-function studies [17,24].

2.3.3 Resilience loss integral and time-integrated resilience metric

Within the unified common time interval [0,T_{\mathrm{m}\mathrm{a}\mathrm{x}}], the resilience loss of the system is defined as the signed cumulative deviation between the ideal pre-shock performance baseline and the simulated performance trajectory [P(t)]. It is calculated by numerical integration:

(11)RL={\int }_0^{T_{max}}\left[1-P\left(t\right)\right]\mathrm{d}t

Here, RL is a net resilience loss rather than a purely positive geometric area. When , the integrand is positive, representing accumulated performance loss. When P\left(t\right)>1, the integrand becomes negative, indicating that adaptive gains above the pre-shock baseline are automatically deducted from the total loss. This definition is consistent with the state-transition CRM, in which the terminal steady-state parameter d allows the system to converge to degraded recovery, elastic recovery, or adaptive improvement.

In the absence of a shock, the reference cumulative performance over T_{max} is P(0)\cdot T_{max}. The time-integrated resilience metric is therefore defined as the reference cumulative performance minus the signed resilience loss:

(12)R=P(0)\cdot T_{max}-RL

3 Theoretical Implications and Computable Features of the CRM

After the mathematical structure of the CRM has been established, its theoretical implications can be further clarified. As conceptually illustrated in Fig. 1, the CRM translates the classic 4R attributes from qualitative descriptors into dynamical characteristics that shape the system performance trajectory.

3.1 From 4R attributes to dynamical characteristics of the performance trajectory

Guided by the CRM formulation established above, this study reinterprets the classic 4R resilience framework proposed by Bruneau et al. [1] from a dynamical perspective. Through the risk-entry-and-exit logic [26], robustness, redundancy, rapidity, and resourcefulness are associated with different stages of the system trajectory from shock exposure to post-shock adaptation. It should be noted that the model’s pre-shock phase is represented mathematically as an ideal baseline. In reality, however, a system maintains a state of dynamic equilibrium prior to a shock, with its performance fluctuating within a defined threshold around a safety benchmark [2,25].

3.1.1 Robustness determines the initial degradation amplitude

Robustness represents the physical resistance of a system against extreme shocks and reflects its capacity to constrain functional degradation after shock energy breaches the system’s defenses [27]. Its effect depends on the interaction between hazard intensity and system defense capability. For example, during the 2021 Texas winter storm in the United States, the widespread failure of power grids lacking freeze protection instantaneously caused a 30 GW supply gap and cascading failures [28]. Similarly, during the extreme rainstorm in Zhengzhou, China, on 20 July 2021, the two-hour rainfall intensity exceeded the local 1000-year design level, flooding 18 stations across five metro lines, suspending train operations, and trapping passengers [29]. Thus, stronger robustness helps constrain the depth of initial performance degradation under extreme shocks.

3.1.2 Redundancy governs the shock buffering capacity

Redundancy, derived from standby nodes or alternative paths within a network, constitutes a system’s capacity to buffer against extreme shocks. From the perspective of the risk-entry-and-exit concept, it characterizes the risk-entry process through which external shocks propagate into the system. Networks with high redundancy can absorb and redistribute disruptive forces through topological reconfiguration, thereby suppressing cascading failures. Empirical evidence indicates that complex power grids or transportation systems meeting N-1 or N-2 security standards can rapidly redistribute network flows through backup circuits when critical nodes are damaged [30,31]. Such redundancy can delay the propagation of disruption and reduce cumulative functional loss.

3.1.3 Rapidity governs the recovery rate

Rapidity reflects the efficiency with which a system mobilizes resources and coordinates recovery after performance degradation. It is jointly driven by macro-governance and resource coordination capability, and it characterizes the efficiency of the risk-exit process. For instance, recovery after the 2008 Wenchuan earthquake demonstrates high rapidity, as a nationwide response system enabled rapid resource mobilization and reconstruction [32,33]. Conversely, recovery after Hurricane Katrina in 2005, relying more on decentralized markets and insurance, exhibited delayed, long-tail trajectories [34,35]. These contrasting cases illustrate the importance of response and resource-mobilization efficiency in shaping the speed of post-shock recovery.

3.1.4 Resourcefulness shapes the asymptotic steady-state

Resourcefulness represents the long-term adaptability of a system following extreme shocks and reflects its capacity for reorganization, learning, and post-disaster improvement [36]. This perspective is closely related to the “Build Back Better” (BBB) paradigm advocated by the United Nations Sendai Framework for Disaster Risk Reduction (2015–2030) [37]. For example, after the 2011 Great East Japan Earthquake, Japan abandoned the in-situ reconstruction of single-layer seawalls and shifted toward multi-layer defense system upgrades, achieving a new steady state with positive gains [38].

Together, these four attributes correspond to distinct intervention points along the shock–recovery–adaptation process, linking resilience theory with practical engineering measures for reducing degradation, buffering disruption, accelerating recovery, and enhancing long-term adaptation.

3.2 Core computable features of the CRM

Based on the above model construction, the CRM exhibits four core features required for genuine computability: macroscopic representation, dynamic time-dependence, scenario simulatability, and multidimensional generalizability. These features respectively concern the scale of system representation, the continuity of temporal evolution, the ability to project intervention scenarios, and the transferability of the model across hazard types and spatial scales.

3.2.1 Macroscopic representation

The framework adopts a phenomenological dynamical paradigm. This avoids over-extrapolating microscopic hazard mechanisms and instead directly characterizes the macroscopic evolution of complex urban systems under shock. This transforms resilience from a qualitative concept into a measurable performance state with engineering dimensions, enabling the quantitative representation of the overall urban resilience level.

3.2.2 Dynamic time-dependence

Moving beyond traditional static measures, the CRM is built on continuous time-dependent equations. It treats shock absorption and system response as a continuous, intertwined process, thereby simulating the full non-equilibrium dynamical trajectory. This includes the initial transient performance degradation, the subsequent non-linear recovery, and the eventual transition to a new steady state.

3.2.3 Scenario simulatability

The model incorporates a parameter-policy mapping, endowing it with forward-looking, scenario-projection capabilities. Specifically, the framework allows system trajectories under different intervention strategies to be simulated by independently adjusting the underlying dynamical control parameters. For example, infrastructure redundancy can influence shock buffering, while emergency resource investment can enhance recovery rapidity. This provides a quantifiable testbed for evaluating intervention strategies and supporting resilience-oriented planning in urban safety governance.

3.2.4 Multidimensional generalizability

By abstracting the common dynamical characteristics of system performance evolution, the framework decouples its core computational architecture from specific hazard types and spatial scales. The CRM supports applications across scales from individual facilities and infrastructure networks to urban agglomerations. It can also be adapted to different shocks, including extreme climate events, geological disasters, public health emergencies, and cyber-physical attacks. This provides a common computational basis for addressing highly uncertain systemic risks.

4 Empirical Application and Case Study

4.1 Case introduction

To illustrate the empirical implementation of the CRM, this section selects nine case cities using a two-dimensional grouping scheme based on geographic region and city size. This sample covers the three major geographic regions (east, central, and west) and three size gradients (megacity/super-large city, large city, and small city), forming a comparative case matrix (Table 3). Regarding scenario settings, the study uses the 2022 urban system as the reference state and considers 10-, 100-, and 500-year extreme-rainfall scenarios. The general CRM framework is operationalized using a rainfall-specific indicator system and calibration procedure. The scenario-dependent shock indicator enters the robustness dimension and therefore allows a to vary across return periods, whereas b, c, and d characterize city-specific redundancy, rapidity, and resourcefulness. The complete indicator definitions, weighting procedure, parameter calibration, numerical settings, and representative numerical trajectory data are provided in the Supplementary Material.

4.2 Case results

To quantify the resilience characteristics of the nine illustrative cities, the dynamic relative change threshold based on degradation depth is applied. With the relative convergence coefficient set to 10⁻³ and the stability condition required for 10 consecutive time steps, the common integration benchmark for the nine-city comparison set is determined as T_{max}=46.5 days. This benchmark is subsequently used to calculate the resilience metric R. As discussed in Section 2.3.2, the benchmark is comparison-set dependent; however, the leave-one-city-out sensitivity analysis confirms that the relative resilience rankings remain stable in the present application (Fig. S1).

Fig. 3 reveals clear differences in the performance trajectories and resilience levels of the nine cities. Within the selected case matrix, Shenzhen and Wenzhou consistently occupy the upper range of resilience across the three rainfall scenarios. Shenzhen records R values of 48.1, 46.8, and 45.3 under the 10-, 100-, and 500-year scenarios, respectively, while the corresponding values for Wenzhou are 46.3, 44.6, and 43.8. Ningde also exhibits relatively stable resilience, with R decreasing from 44.7 to 43.2 as rainfall intensity increases. Overall, the three selected eastern cities show relatively high or stable resilience within the case matrix, characterized by comparatively shallow performance degradation and limited reductions in resilience under intensified rainfall shocks.

Resilience differences across city-size categories are more heterogeneous. Among the megacity and super-large-city cases, Shenzhen maintains consistently high resilience, whereas Chengdu shows greater sensitivity to increasing shock intensity, with R declining from 44.7 under the 10-year scenario to 33.6 under the 500-year scenario. Zhengzhou remains at an intermediate level, decreasing from 44.3 to 38.7. This pattern is also broadly consistent with the vulnerability exposed during the 2021 Zhengzhou extreme rainfall event, which caused severe casualties and substantial property losses [7,39]. Among the large-city cases, Wenzhou exhibits relatively high and stable resilience, while Mianyang shows a more pronounced decline from 42.5 to 30.9. For the small-city cases, Ningde remains comparatively stable, whereas Jinchang records lower resilience levels across the three scenarios. These contrasts indicate that the modeled resilience differences reflect the combined effects of robustness, redundancy, rapidity, and resourcefulness rather than urban scale alone.

Increasing rainfall intensity produces a general decline in resilience across all nine cases, but the magnitude of the response varies substantially. From the 10-year to the 500-year scenario, Mianyang and Chengdu exhibit the largest decreases in R, by approximately 11.6 and 11.1, respectively, followed by Chuzhou with a decrease of about 6.8. By contrast, Ningde, Wenzhou, Shenzhen, and Yueyang show smaller reductions of approximately 1.5, 2.5, 2.9, and 3.1, respectively. These differences reflect heterogeneous sensitivity to intensified rainfall shocks and illustrate how the CRM translates variations in external shock intensity and city-specific dynamical attributes into differentiated performance trajectories and cumulative resilience outcomes.

For cities with d>1, the post-shock performance trajectory may eventually exceed the normalized pre-shock baseline. For example, Shenzhen has d=1.10 under the present calibration, and its 10-year scenario yields a signed resilience loss of approximately −1.61, resulting in R=48.11, slightly higher than T_{max}=46.5. This represents the adaptive-improvement state embedded in the CRM, in which post-shock gains above the baseline offset part or all of the preceding performance loss.

5 Conclusions and Future Perspectives

The United Nations Sendai Framework for Disaster Risk Reduction (2015–2030) and the Sustainable Development Goals (SDGs) have established resilience as a core indicator of modern urban development. Achieving these goals requires quantitative tools that can diagnose system performance under extreme shocks. The CRM developed in this study contributes to the transition of urban resilience assessment from static characterization toward a dynamic computable paradigm. By quantifying the evolutionary mechanisms of the entire process—encompassing external shocks, performance degradation, functional recovery, and adaptive enhancement—this model provides a standardized computable method for resilience measurement within the context of urban safety research.

At the macroscopic scale, particularly for metropolitan areas, regional clusters, and national territories, the CRM enables relative resilience comparison under standardized shock scenarios and a common temporal benchmark [7,25]. This application underscores that defining uniform exogenous shock conditions is a fundamental prerequisite for macro-evaluations to yield meaningful physical and engineering insights.

At the microscopic scale, particularly for critical infrastructure systems such as power grids, rail transit, and supply chains, the CRM can be applied using observed system-performance data or measurable resilience attributes. The dynamical characteristics associated with the 4R attributes can be quantified from corresponding indicators and further combined through additive or multiplicative formulations to obtain an overall resilience index [17,18,24,40].

Across these macroscopic and microscopic applications, several methodological limitations define the current scope of the CRM. At the macroscopic scale, aggregating multidimensional observations into four dynamical parameters may smooth heterogeneity among infrastructures, neighborhoods, and functional subsystems, while finer-scale applications depend more strongly on the availability and quality of observed shock–recovery trajectories. More generally, indicator selection, weighting schemes, parameter-mapping ranges, and recovery-kernel specification may influence the calibrated trajectories and resilience values. The data-adaptive common benchmark T_{max} is also comparison-set dependent, although the leave-one-city-out sensitivity analysis shows that the relative resilience rankings remain highly stable in the present application. These considerations provide clear directions for further development through multi-scale calibration, uncertainty propagation, empirical trajectory fitting, and the use of externally specified engineering or planning horizons where available.

The fundamental value of the CRM lies in transforming abstract governance objectives into executable engineering planning tools. In the future, this framework can be extended to heterogeneous systems such as urban energy, transportation, and logistics under complex shocks. Moreover, future empirical validation can use complete observed shock–recovery trajectories to calibrate the CRM inversely and benchmark its fitting performance against established recovery-function models. The CRM can thereby provide a common basis for cross-scale assessment of urban safety states and offer a forward-looking platform for prevention and mitigation planning against extreme shocks. Decision makers can use the CRM to simulate in advance how different resource inputs and infrastructure upgrades affect the 4R attributes and the resulting recovery trajectory, thereby supporting optimal decision-making paths and dynamic investment optimization in urban safety governance with quantitative data.

References

[1] Bruneau M , Chang S E , Eguchi R T . et al. A framework to quantitatively assess and enhance the seismic resilience of communities. Earthquake Spectra, 2003, 19(4): 733–752

[2] Li L Z , Ding Y , Yuan J F . et al. Quantifying the resilience of emergency response networks to infrastructure interruptions through an enhanced metanetwork-based framework. Journal of Management in Engineering, 2022, 38(5): 04022047

[3] Zhou Y , Xiao Y , Xu H R . Quantifying urban disaster resilience for informed mitigation strategies. Engineering, 2025, 53: 13–17

[4] Meerow S , Newell J P , Stults M . Defining urban resilience: a review. Landscape and Urban Planning, 2016, 147: 38–49

[5] Kemp L , Xu C , Depledge J . et al. Climate Endgame: exploring catastrophic climate change scenarios. Proceedings of the National Academy of Sciences of the United States of America, 2022, 119(34): e2108146119

[6] Wang Y M , Ye Z J , Jia X R . et al. Flood disaster chain deduction based on cascading failures in urban critical infrastructure. Reliability Engineering & System Safety, 2025, 261: 111160

[7] Zhuang L M , Wang M , Zhang Z T . et al. Measurement of multidimensional reconstruction fund planning on the effectiveness of urban post-disaster recovery and reconstruction using indirect economic impact. Environmental Impact Assessment Review, 2024, 105: 107365

[8] Ouyang M , Cheng Z K , Ma J X . et al. Coupled urban risks: a complex systems perspective with a people-centric focus. Engineering, 2025, 44: 44–50

[9] Cutter S L , Barnes L , Berry M . et al. A place-based model for understanding community resilience to natural disasters. Global Environmental Change, 2008, 18(4): 598–606

[10] Hosseini S , Barker K , Ramirez-Marquez J E . A review of definitions and measures of system resilience. Reliability Engineering & System Safety, 2016, 145: 47–61

[11] Linkov I , Eisenberg D A , Bates M E . et al. Measurable resilience for actionable policy. Environmental Science & Technology, 2013, 47(18): 10108–10110

[12] Buldyrev S V , Parshani R , Paul G . et al. Catastrophic cascade of failures in interdependent networks. Nature, 2010, 464(7291): 1025–1028

[13] Chen R X , Wang P P , Zheng X Q . et al. Dynamic simulation of urban agglomeration network resilience under disturbances: an integrated deep learning and agent-based approach – Case study of the Beijing−Tianjin−Hebei region. Sustainable Cities and Society, 2026, 140: 107253

[14] Xu J , Ding C , Guo X W . et al. Urban waterlogging disaster risk evaluation framework based on complex network and agent-based model: a case study of Zhengzhou, China. International Journal of Disaster Risk Reduction, 2025, 131: 105885

[15] Zio E . Challenges in the vulnerability and risk analysis of critical infrastructures. Reliability Engineering & System Safety, 2016, 152: 137–150

[16] Taberna A , Filatova T , Hadjimichael A . et al. Uncertainty in boundedly rational household adaptation to environmental shocks. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120(44): e2215675120

[17] Cassottana B , Shen L J , Tang L C . Modeling the recovery process: a key dimension of resilience. Reliability Engineering & System Safety, 2019, 190: 106528

[18] Cimellaro G P , Reinhorn A M , Bruneau M . Seismic resilience of a hospital system. Structure and Infrastructure Engineering, 2010, 6(1−2): 127–144

[19] Ayyub B M . Systems resilience for multihazard environments: definition, metrics, and valuation for decision making. Risk Analysis, 2014, 34(2): 340–355

[20] Ayyub B M . Practical resilience metrics for planning, design, and decision making. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering, 2015, 1(3): 04015008

[21] Reed D A , Powell M D , Westerman J M . Energy supply system performance for hurricane Katrina. Journal of Energy Engineering, 2010, 136(4): 95–102

[22] Eisenberg D A , Seager T P , Alderson D L . The rebound curve is a poor model of resilience. PNAS Nexus, 2025, 4(3): pgaf052

[23] Todman L C , Fraser F C , Corstanje R . et al. Defining and quantifying the resilience of responses to disturbance: a conceptual and modelling approach from soil science. Scientific Reports, 2016, 6: 28426

[24] Tang L C , Shen L J . Resilience indices from a family of recovery functions. Fundamental Research, 2024, 4(1): 13–20

[25] Zhuang L M , Wang M , Liu K . et al. Assessing and optimizing urban dynamic resilience to extreme rainfall from shock to recovery. International Journal of Disaster Risk Science, 2026, 17(1): 128–147

[26] Shi P J, Jaeger C, Ye Q. Integrated risk governance. Berlin, Heidelberg: Springer, 2013

[27] Tierney K J , Bruneau M. . Conceptualizing and measuring resilience: a key to disaster loss reduction.. TR News, 2007, (250): 14–17

[28] Busby J W , Baker K , Bazilian M D . et al. Cascading risks: understanding the 2021 winter blackout in Texas. Energy Research & Social Science, 2021, 77: 102106

[29] Lyu H M , Yin Z Y , Shen S L . et al. A novel framework incorporating machine learning into GIS for flood susceptibility prediction of urban metro systems. Science China Technological Sciences, 2025, 68(6): 1620701

[30] Mattsson L G , Jenelius E . Vulnerability and resilience of transport systems – A discussion of recent research. Transportation Research Part A: Policy and Practice, 2015, 81: 16–34

[31] Ouyang M . Review on modeling and simulation of interdependent critical infrastructure systems. Reliability Engineering & System Safety, 2014, 121: 43–60

[32] Dunford M , Li L . Earthquake reconstruction in Wenchuan: assessing the state overall plan and addressing the ‘forgotten phase’. Applied Geography, 2011, 31(3): 998–1009

[33] Zhang H B , Tang G J . Paired assistance policy and recovery from the 2008 Wenchuan earthquake: a network perspective. Disasters, 2021, 45(1): 126–157

[34] Kates R W , Colten C E , Laska S . et al. Reconstruction of New Orleans after Hurricane Katrina: a research perspective. Proceedings of the National Academy of Sciences of the United States of America, 2006, 103(40): 14653–14660

[35] Fussell E , Sastry N , VanLandingham M . Race, socioeconomic status, and return migration to New Orleans after Hurricane Katrina. Population and Environment, 2010, 31(1−3): 20–42

[36] Folke C . Resilience: the emergence of a perspective for social–ecological systems analyses. Global Environmental Change, 2006, 16(3): 253–267

[37] Hallegatte S, Rentschler J, Walsh B. Building back better: achieving resilience through stronger, faster, and more inclusive post-disaster reconstruction. World Bank Group, 2018. Available online (accessed 10 Apr 2026)

[38] Suppasri A , Shuto N , Imamura F . et al. Lessons learned from the 2011 great east japan tsunami: performance of tsunami countermeasures, coastal buildings, and tsunami evacuation in Japan. Pure and Applied Geophysics, 2013, 170(6): 993–1018

[39] Shen S V . The 2021 Henan flood increased citizen demand for government-led climate change adaptation in China. Communications Earth & Environment, 2025, 6(1): 730

[40] Argyroudis S . Resilience metrics for transport networks: a review and practical examples for bridges. Proceedings of the Institution of Civil Engineers - Bridge Engineering, 2022, 175(3): 179–192

Rights & permissions

The Author(s) 2026. This article is published by Higher Education Press.