Resilience patterns of human mobility in response to extreme urban floods

ABSTRACT Large-scale disasters can disproportionately impact different population groups, causing prominent disparity and inequality, especially for the vulnerable and marginalized. Here, we investigate the resilience of human mobility under the disturbance of the unprecedented ‘720’ Zhengzhou flood in China in 2021 using records of 1.32 billion mobile phone signaling generated by 4.35 million people. We find that although pluvial floods can trigger mobility reductions, the overall structural dynamics of mobility networks remain relatively stable. We also find that the low levels of mobility resilience in female, adolescent and older adult groups are mainly due to their insufficient capabilities to maintain business-as-usual travel frequency during the flood. Most importantly, we reveal three types of counter-intuitive, yet widely existing, resilience patterns of human mobility (namely, ‘reverse bathtub’, ‘ever-increasing’ and ‘ever-decreasing’ patterns), and demonstrate a universal mechanism of disaster-avoidance response by further corroborating that those abnormal resilience patterns are not associated with people’s gender or age. In view of the common association between travel behaviors and travelers’ socio-demographic characteristics, our findings provide a caveat for scholars when disclosing disparities in human travel behaviors during flood-induced emergencies.


Study area
The study area, the Zhengzhou region ( Figure S1), is located in Henan province, the north-central area of China (112°42′E-114°14′E, 34°16′N-34°58′N), and covers approximately 7,446 km 2 . In 2021, the region's resident population was about 12.7 million, with an urbanization rate of 79.1%. This area has a semiarid subtropical climate, an annual mean temperature of 14.35°C, and an annual mean rainfall of 640.9 mm, with the Yellow River running through it. The region is the political, economic, technological, and educational center of the province and also the most critical transportation hub in central mainland China. Since 1992, China's central government has designated Zhengzhou as one of 11 "Hinterland Open Cities", allowing access to the same preferential policies as the eastern coastal open cities. Currently, with the promotion of the "One Belt, One Road" initiative, the Zhengzhou region has become an important node on the "New Silk Road" between Europe and Asia.

Event description and timeline
The information on the event timeline was collected from the Disaster Investigation Team of the State Council. As depicted in Figure S2, the rainstorm center moved southward to the Zhengzhou region on July 19, 2021, which led to a long-lasting heavy rainfall event. The severe pluvial flood in Zhengzhou began on July 20, gradually lessened, and ended on July 23 (the event is soon reported with the name '720' Zhengzhou flood). From 8:00 am on July 17 to 8:00 am on July 23, the accumulated precipitation exceeded more than 600 mm and reached 5,590 km 2 , and Erqi, Zhongyuan, and Jinshui districts in the central urban area had an accumulated rainfall of nearly 700 mm. The heaviest rainfall period happened between the afternoon of July 19 and the morning of July 21, with a maximum daily rainfall of 624.1 mm being recorded at Zhengzhou National Weather Station on the 20 th , which is almost equal to the average annual rainfall in the Zhengzhou region (640.8 mm). This event broke the historical record of extreme meteorological observations in mainland China.
The event timeline and emergency responses are recorded in Figure S3. A quick summary of this extreme weather event is as follows: (1) This heavy rainfall caused a wide-range large-scale urban flood and was extremely intense during a relatively short period. (2) Major rivers in the region repeatedly faced grim inundation situations. (3) The inner city drainage system was unable to deal with this sudden burden. (4) Even though there were several rounds of emergency alerts before the heavy rainfall, the event still caused massive losses to society. Figure S2. Time series of the rainfall process and collection data of rescue information from July 13 to 30, 2021. (a) daily weather conditions of precipitation and visibility; and (b) the rescue information was collected from social media platforms. Different colors were used to distinguish the type of search and rescue activities. Figure S3. Chain of events displayed at the national (top), provincial (middle), and regional (bottom) level from July 19 to 29, 2021 in the Zhengzhou region. Red blocks represent early warning forecasts, blue blocks represent emergency response, and green blocks relate to disaster relief information.

Spatial-temporal changes in overall mobility
We discretize the study area into 1 km-by-1 km grids and subsequently use the collective human mobility OD data at pre-, during, and post-disaster stages to analyze the frequency of average daily flux.
It is found that the frequency distribution of daily average flux is extremely polarized, with most of the grids occupied with fewer OD movements. We found two critical points in the distribution plot at 50 and 100 intensity levels. Thus, we reclassify the mobility data into three groups: (1) less than 50, (2) more than or equal to 50 but less than 100, and (3) more than or equal to 100 OD trips ( Figure S4).
To support Figure 1 in the main manuscript, we also analyze the mobility network changes together with the network centrality measures for < 50 and 50 < < 100 in Figures S5 and S6. Similar patterns as discussed in the main manuscript can be found at these two levels of intensity.  (d-f) exceedance probability with size for degree and betweenness in pre-event, during-event, and post-event stages respectively; (g-j) percentage change of degree, edge weight, betweenness, similarity from July 13 to 30. Note that the red blocks highlight the during-event process.

Resilience quantification metric
We use the time series systemic performance curve to quantitatively evaluate human mobility resilience (Hong et al., 2021;Kontokosta & Malik, 2018). As shown in Figure S7, a typical resilience curve reflects the system behavior and describes the continuous system performance level before, during, and after an external disruptive event (Poulin & Kane, 2021). The shaded area refers to the profile of performance loss, and the performance-based resilience can be calculated as follows (Ouyang & Duenas-Osorio, 2012): where is the score of the human mobility resilience defined in this paper,  We used the percentage change of mobility flux in this study as the y-axis of the performance profile, which features excellent effectiveness and simplicity for applications. We convert the hourly total flux to plot the resilience curves of human flow networks. Since the resilience metric we adopted has a value range of [0,1], it is intuitive and easy to conduct cross-comparison. Thus, the percentage change can be calculated as follows: where % is the hourly total flux of each day and # is the hourly average flux of July 16 to 18 (with the pre-event stage as the benchmark). The normalization can be performed as: The steps for preparing the performance profile (i.e., the resilience curve) and calculating the resilience scores are shown below.
Step 1: We calculate the total flux (i.e., the total edge weight of the network) per hour per day for each demographical group. For each group, there are 408 data points for a total of 24 hours per day for the full-length observation window of 17 days.
Step 2: Then, we set a pre-event mobility-performance benchmark to facilitate the conversion from mobility flux to percentage changes using Eq. (2). The benchmark is the averaged hourly flux from July 16 to 18 (i.e., three consecutive days before the disaster hit).
Step 3: Next, we use a simple moving average algorithm to smooth and denoise the time-series curves with a sliding window size of 24 hours.
Step 4: Finally, we apply Eq. (3) to normalize the data and Eq. (1) with the 'trapez' function to calculate the resilience scores for all groups. All calculations were performed in R (v 4.2.2).

One-way ANOVA test on the mobility difference and resilience scores among demographical groups
We applied a one-way ANOVA test to statistically verify the differences among the percentage changes of the hourly total flux in all demographic groups between 13th July and 29th July. By doing so, we could test the statistical differences in those calculated resilience scores as well. The null hypothesis is: "The means (of the percentage changes of hourly total flux) of all the demographic groups are equal", and the alternative hypothesis is: "The means of two or more groups are significantly different from that of the others". After the ANOVA test, we then applied Tukey's HSD test at the 95% confidence level to further indicate the differences in population mean in a pairwise manner for a closer check. The test returns results on the difference in the observed means (Diff), the Lower bound and Upper bound of the confidence interval, and p-values (Table S1). R (v 4.2.2), and Microsoft Excel (v16.16.1) were jointly used here.
From the ANOVA and Tukey's HSD results, we can see that overall there are differences among demographic groups, and this difference was most statistically significant between the pair among gender groups. Eight pairs among age groups also show a statistically significant difference, while pairs between age groups having a minor difference in resilience scores (<0.05) show a statistically insignificant difference in population means, suggesting that those age groups with similar population attributes yield a similar level of resilience.

Network generation modeling: Preferential attachment model
The ability to generate random networks that can accurately represent real social contact networks is essential to the study of natural or human-made systems. The Barabási-Albert model and its extensions attempt to reflect reality by generating random scale-free networks using a preferential attachment mechanism. Here, we choose the dual-Barabási-Albert (DBA) model (Moshiri, 2018) rather than Barabási-Albert (BA) model because the latter allows only two inputs, the total number of nodes and the number of edges, to attach from a new node to existing nodes. Therefore, the BA model can only describe the strength of the network generation process and cannot be applied to measure generation uncertainties of the network's growth under different growth rates (i.e., fast growth or slow growth ). The DBA model is better suited to observe such uncertainties with two controled parameters (i.e., 1 and 2).

Model construction
We simulate the growth and shrinkage processes of the mobility networks using the DBA model, a preferential attachment-based network generation model, to capture the evolutionary processes of the mobility networks at pre-, during-, and post-event stages. The parameters of the model encapsulated in the Python package 'Networkx' are set as follows: where is the number of nodes, 1 is the number of edges to link the new node to existing nodes with a probability of , 2 is the number of edges to link each new node to existing nodes with a probability of (1 − ) 1 ( Figure S8). On the contrary, the parameter 2 was set as a constant (i.e. 2 =1), which characterizes the generation of non-multi-edge networks. The reason for such a setting is that the magnitude of 2 is not the main concern here as long as it can effectively represent a slow-enough growth rate. 2 = 1 represents a situation where only one new edge is generated in the network at each time step (i.e., as a strong contrast to a much larger 1 that represents a much faster growth rate). Moreover, by controling 2 we are able to lay a benchmark for effective comparisons in later analyses, otherwise, we would not be able to compare the performance of based on bivariate ( 1, 2) pairs. The specific parameter settings of the DBA models for different demographical groups are shown in the following Table S2.

The fitting process of the parameter
For each demographic group, we perform 136 (files with average OD data in 3h) times 50 (values of between 0 and 1) times 50 (Monte Carlo simulations) = 340,000 simulation runs, which would cost 14 hours for a single-core CPU with full loading under normal conditions.
The specific fitting operation for parameter is as follows.
Step 1 -Data preparation: set the time window of the observed data to 3h and calculate the average number of trips. Re-organize the data columns into 'Sources' (i.e. the origin), 'Destinations' (i.e. the destination), and 'Weights' (i.e. the intensity of an OD connection).
Step 2 -Network modeling: construct a network model of the empirical data using 'Networkx' Python package and count the number of edges in the network.
Step 3 -Fitting: based on the parameters in Table S2, obtain the optimal parameter having the minimum discrepancy between the DBA-generated network and the ground-truth network by iterating through the range of .
Step 4 -Monte Carlo simulations: perform 50 iterations and obtained the converged values as the final outcomes.

Definition of the quadrant plot and the outliers in drawdown and drawup processes
The resilience curve uses each day minus the previous day and calculates a change ratio, generating a curve showing the change in the flow of people each day, and finally using the resilience triangle method to measure its resilience. Furthermore, for analyzing the resilience patterns (e.g. bathtub curves) in this section, we first calculate the daily average mobility flow

Community detection algorithm for identifying sub-network clusters
Here, we test how mobility patterns might change at a smaller spatial scale (the sub-cluster level of the mobility network) to see if the ratio between normal and abnormal resilience patterns could be scale-invariant. To do so, we first apply the modularity algorithm to detect the communities (i.e. the sub-network clusters) in the mobility networks under normal conditions (i.e. the pre-event condition) and classify the spatial grids into different clusters.
We apply the modularity algorithm for detection, which is often used in optimization methods for detecting community structure in networks. It can be expressed as the following algorithm (Newman, 2006): As a result, seven clusters of the pre-disaster aggregated network are detected by applying Gephi with the modularity algorithm at a resolution = 2. We found a significant spatial coupling relationship between the administrative boundaries and the structure of the clusters, which indicates that most of the human mobility activities are confined by administrative boundaries (i.e. most movement occurs within the clusters, that is, within the administrative boundaries), while the movement across cluster divisions accounts for only a small number but with relatively stable intensity and frequency.

Significance of differences in clusters using the K-S test
We then test whether the differences among clusters are statistically significant. Table S3 presents a comparative description of all clusters in terms of the outflows in the drawdown (i.e. the pre-event level minus the during-event level) process, and it can be seen that Cluster 1 (i.e. the one in the central urban area) has the highest amount of outflow during the disaster event and the least resilience score in the whole process, which verifies that the central urban area was severely hit and damaged. Clusters 2, 3, and 6 have a relatively higher resilience level, with Cluster 6 being the highest of all (1,611 grids; 13.36% of the total grids; mean = 10.665, with RI = 0.75). This result is in line with the observations from the spatial autocorrelation analysis suggesting that these clusters have relatively sparse disaster site distributions.  We conclude that, even though the differences in flow distributions among clusters are significant from a holistic viewpoint, there could also be similar mobility patterns at a smaller spatial scale.

Mobility resilience in clusters
In terms of mobility resilience in clusters, Figure S10 and Figure S11 suggest that the differences in resilience between outflow and inflow are rather minor, which indicates a highly similar mobility loss in terms of people's arrivals and departures within a cluster during the whole observation window. Furthermore, it can be argued that none of the seven clusters fully    Figure S12 and Figure S13 show that the proportion of the four quadrants in each cluster is mostly close to the ratio of 5:3:1:1 found at the global region scale, which verifies that the ratio of normal and abnormal resilience patterns could remain stable at a smaller spatial scale.

Step-wise regression results and statistical tests
We adopt step-wise multinomial logistics regression to further examine the explanatory strength of demographic attributes in relation to the multi-class of quadrants. Note that, of the total records in each attribute mobility dataset (i.e. gender or age), 75% are used in the training process, while the remaining 25% are used in the validation step for calculating models' accuracy and AUC.
The step-wise regression results for the gender group are shown in Similar findings are also found regarding age groups. The step-wise regression results for the age group are shown in Table S7, and the models' standard error and AUC are shown in  proportions of the elderly in Q1, Q2, and Q3 are slightly higher than that of the two genders, and (2) the distribution of the elderly in the abnormal quadrant is mostly found on the periphery of Zhengzhou region when compared with male or female. In summary, these findings support the discussion in the main manuscript that the demographic attributes of the travelers are not strongly associated with the distribution of heterogeneous resilience patterns. Part IV. Additional information for the discussion section

Spatial autocorrelation analysis
We further collect rescue data and conduct a Bivariate Local Moran's I test, a measure of spatial autocorrelation analysis to investigate how mobility changes are associated with specific damage hits (Moran, 1950;Li et al., 2007). The Local Moran's I can be defined as the following algorithm: where # is the Local Moran's I, and is the number of analysis units on the map.
Here, we use the location information of rescue calls as the sites that suffered from severe disaster damage, because rescue data were mostly from travelers who had been trapped in the flood or witnessed someone being severely impacted by the flood nearby and were then reported as rescue-seeking calls, which is a rational and intuitive assumption for the subsequent analysis. By setting 1km buffer zones at all rescue points, we find that, holistically, there exists a mild positive correlation (Moran's I = 0.325) between the distribution of the OD flux changes and the distribution of the disaster sites in Figure S15. Most disaster sites spatially overlap with those blue grids (decreasing outflows and calmer mobility activities) in the central urban area.
By cross-referencing Figure 7 and Figure S14, we can infer that grids with a normal resilience pattern might also be correlated with the disaster sites, which further indicates that the particular ratio between the normal and abnormal resilience patterns might be associated with a portion of the disaster-impacted areas in the whole study area.