The issue of pedestrian evacuation from hazardous events has been studied with a focus on efficiently moving large groups from enclosed spaces (e.g., vehicles, buildings, underground malls) to limited exits, as reviewed by [2]. These insights inform the design of exits, evacuation routes, and signage systems.

Since the early 2000s, disasters such as urban fires, floods, typhoons, tsunamis, and hazardous material accidents have heightened attention to horizontal evacuations across road networks. [3] conducted a comprehensive global review of evacuation planning over the past 30 years, while [4] provided a systematic review of simulation-based optimization studies.

When road network conditions are static and link travel times constant, choosing routes with the shortest travel time is efficient. However, as traffic density increases, movement becomes constrained, leading to slower speeds and longer travel times. Thus, link travel times are modeled as increasing functions of traffic volume, diverging near critical volume thresholds to reflect congestion. By integrating this with conservation laws, continuity conditions, and cost functions, dynamic traffic assignments can be computed; for example, the model by [5].

Traffic flow, defined as density times spatial average speed, increases under free-flow conditions but decreases in congested regions once critical density is surpassed. On roads with uniform width and slope, the dynamics of congestion buildup and dissipation have been theoretically analyzed using kinematic wave theory and the fundamental diagram (FD). Parameter estimation enables simulation of evolving traffic flow. In traffic planning, goals typically include maximizing completed trips within a timeframe, minimizing total travel time, or minimizing journey completion time. Maintaining high flow rates requires keeping density below critical thresholds. Hence, many models focus on subcritical capacity conditions. [6] proposed a linear programming approach using a linear link performance function. Related optimization-based models include those by [7], [8], [9], and [10].

The aforementioned traffic distribution models assume that density remains below critical levels and that congestion is avoided. However, during disaster evacuations, urgency often compels users to enter downstream roads as long as any capacity exists, increasing the likelihood of congestion beyond critical density. If congestion queues can be restricted to safe zones, a rapid reduction in vehicle volume within upstream hazard areas may mitigate damage. Moreover, if congestion in hazardous zones clears before the hazard strikes due to timely upstream evacuation, its negative impact can be avoided. Therefore, when planning evacuations using objective functions such as expected evacuees or expected casualties, congestion is not inherently detrimental and may even yield advantages. It is thus desirable to use dynamic models capable of capturing traffic conditions that include congestion. [11] developed a Cell Transmission (CT) model for a single-lane setting, representing time-dependent changes in traffic flow per road segment (cell) using a flow–density relationship aligned with kinematic wave theory. [12] demonstrated that by modeling cell FD relationships as piecewise linear functions and defining inter-cell connections as linear functions, the overall traffic flow, including upstream congestion due to downstream saturation, can be formulated as a linear programming problem and solved efficiently.

In this study, we extend the [12] model to construct an optimization model. The formulation employs the notations presented in Table 1. Indices for the cell set A are denoted as a , b and for the node set N as r , s . The node adjacent downstream of cell a is represented as a + , and the sets of cells adjacent upstream and downstream of node s are represented as I s and O s , respectively. Discrete time is represented as t = 0 , … , T with one time unit defined as the duration for free flow to pass through one cell. The endogenous variables are y t a : queued flows in cell a at the end of time t , corresponding to p t i − v t   i in our model; v t   a ​: exit flow from cell a during time t ; u t   a : inflow rate into cell a at the end of time t or at the beginning of time t + 1 ; and z t r s : travel demand between origin-destination pair r s departing at time t . Other parameters include C a : flow capacity of cell a ; H a : holding capacity of cell a ; D a : demand originated from cell a , and δ : the ratio of the wave speed to the free flow speed. Nie model is formulated as follows: min ∑ t = 1 T ∑ a ∈ A y t a + v t a ≡ min ∑ t = 1 T ∑ a ∈ A p t a (1) Subject to  y t + 1 a = y t a + u t a − v t a , t = 0 , ⋯ T − 1 , ∀ a ϵ A (2) ∑ a ∈ O s u t a = z t r s + ∑ a ∈ I s v t a , t = 0 , ⋯ T − 1 , ∀ s ϵ N , r ≠ s (3) v 0 a + y 0 a = D a ≥ 0   g i v e n , ∀ a ϵ A (4) 0 ≤ u t a ≤ C a , t = 0 , ⋯ T − 1 , ∀ a ϵ A (5) 0 ≤ v t a ≤ C a , t = 0 , ⋯ T − 1 , ∀ a ϵ A (6) u t a ≤ δ H a − y t a − v t a , t = 0 , ⋯ T − 1 , ∀ a ϵ A (7) y t a > 0 ⇒   v t a = C a , or  u t b ≤ C b , ∀ b ∈ O a + , or  u t b = δ H b − y t b − v t b , ∀ b ∈ O a + , t = 0 , ⋯ T − 1   ∀ a ϵ A (8)

The objective function Equation 1 represents the total number of vehicles on the road and corresponds to minimizing total travel time. Equation 2 expresses temporal conservation of vehicles, while Equation 3 ensures conservation of flow at nodes. Equation 4 captures initial demand, and Equations 5, 6 constrain exit and inflow rates by flow capacity. Specifically, Equation 7 limits inflow based on available holding capacity. Equation 8, known as the ordered solution property (OSP), ensures that in congested cells, outflow is maximized. This prevents the “flow holding back” issue, where vehicles remain in upstream cells despite available downstream capacity.

The model aligns with a trapezoidal FD. Since y t a cannot be negative, the following condition must hold: v t a ≤ p t a , t = 0 , ⋯ T − 1 , ∀ a ϵ A (9)

Equation 9 implies that the outflow rate during a unit time cannot exceed the number of vehicles present at the beginning of the period, corresponding to the increasing portion (left edge) of the trapezoidal FD. Equation 6 sets the upper limit of outflow (top edge), while Equation 5 restricts inflow based on remaining capacity (right congested edge).

The FD models congestion in a single cell, considering vehicles as forming a point queue above it. The Nie model Equation 1 captures how rising cell density limits inflow (via Equation 7), which in turn constrains outflow from upstream cells through conservation at adjacent cells (Equation 3), allowing for upstream propagation of congestion.

Although the Nie model comprises linear expressions, Equation 8 introduces a conditional constraint that deviates from standard linear programming. [12] excluded Equation 8 from the initial formulation and applied an iterative algorithm to accelerate outflows from congested cells until the OSP condition was satisfied.

While these studies focus on single road sections, scaling to larger areas presents challenges in handling detailed road networks. To address this, MFD models were developed to express the relationship between the number of vehicles in an area (spatial density) and traffic volume (trip departures and arrivals). Studies by [13], [14], [15], [16], and [17] have explored the fundamental characteristics of such models. [18] created an agent-based simulation incorporating a trapezoidal MFD to estimate dynamic traffic distributions at the urban scale. These MFD-based efforts are reviewed by [19]. As an application, [20] proposed an evacuation optimization model using a nonlinear MFD. [21] applied a triangular MFD, incorporating departure time choices to analyze optimal transitions in road usage as urban flooding evolves over time.

Compared to the above research on vehicle evacuations, pedestrian evacuation studies in road networks remain relatively limited. [22] and [23] adapted dynamic traffic assignment models for vehicular traffic to pedestrian contexts. [24] empirically studied crowd dynamics during the stoning rituals, contributing to flow–density modeling for pedestrians. [25] investigated maximum pedestrian flow under free-flow conditions, below critical capacity. [26] used a CT model to simulate pedestrian movement in confined areas such as public squares. [27] enhanced the applicability of trapezoidal pedestrian MFDs by introducing local density classifications around individuals. [19] reviewed pedestrian MFD developments and confirmed the feasibility of pedestrian evacuation modeling.

As described above, the development of zonal MFD models, evolved from CT models originally applied to single road segments, has advanced in pedestrian traffic modeling. The application of trapezoidal MFDs to evacuation contexts has also progressed. This study focuses on the trapezoidal FD-based model by [12], notable for its computational efficiency and operability as a linear programming formulation. By adapting this approach into a zonal MFD model for pedestrian evacuation, we aim to build an optimal, efficient, and tractable evacuation model.

Unlike [12], this study’s model does not incorporate the OSP condition. As discussed later, our objective function minimizes the expected number of tsunami victims during evacuation. In scenarios where a time buffer exists before tsunami arrival, it may be optimal to move toward zones closer to the coastline, despite higher eventual risk, due to the temporary equivalence in disaster probability across zones. As the hazard evolves, the rational direction of evacuation changes dynamically, and the OSP condition may not always be appropriate. Therefore, excluding the OSP condition allows the model to retain a standard linear programming structure, simplifying computation of the optimal solution.

This study develops a mathematical optimization model to determine the most efficient evacuation strategy while ensuring physical feasibility. The target area is divided into zones based on a regional mesh, which are rectangular units approximately 500 m in length, commonly used in urban planning for population aggregation. Each zone is divided into (a) road section, where pedestrians move between zones, (b) off-road land-use section, and (c) evacuation shelter section, which include safe facilities such as shelters, disaster prevention parks, or evacuation towers.

In the model, T represents the set of time periods from earthquake onset to the end of the first tsunami wave’s arrival, Z E denotes the set of zones containing shelter sections, and Z indicates general zones without shelters. Additionally, J i is the set of zones adjacent to zone i , and adjacent zones in the cardinal directions are denoted as N i , S i , W i , and E i . Table 1 lists the symbols for all variables and parameters used in the formulation. Figure 1 illustrates a typical zone structure and associated variable definitions.

The probability that an evacuee in mesh i at time t will encounter a tsunami (hereafter referred to as “tsunami encounter probability” R t i ) is defined as follows. While a flood depth of 0.3 m or more indicates tsunami impact, we use a continuous Logit-type function (rather than a step function) to enable optimization. This function, based on predicted maximum inundation depth by [28], is expressed by Equation 10. R t i = 1 / 1 + exp − 30 * m a x i m u m   i n u n d a t i o n   d e p t h   b e f o r e   t i m e   t − 0.3 (10)

The tsunami encounter probability is set to 0 for meshes outside the inundated area, evacuation shelter sections within inundated zones, and all zones for t < 30 , before tsunami run-up begins.

Next, the expected number of tsunami victims based on a specific population distribution when evacuation does not occur is referred to as the “static risk value.” Specifically, this is calculated by first determining the average value of the tsunami encounter probability for each zone until the future and then summing the products of these averages with the given resident population, as expressed by Equation 11. S R V = ∑ i ∈ Z ∪ ZE     1 T 1 − T 0 ∑ t = 0 T 1   R t i p t i + q t i (11)

When evacuation occurs, population distribution changes over time. The dynamic risk value (DRV) estimates the expected number of tsunami victims under such conditions. It is calculated by summing, at each time point, the product of population and tsunami encounter probability, then averaging over the inundation period, as expressed by Equation 12. Minimizing the DRV is the objective of our optimization model. D R V = ∑ t = 0 T 1   ∑ i ∈ Z ∪ ZE   1 T 1 − T 0 R t i p t i + q t i (12)

This model adopts similar constraints to those in the [12] model, including continuity conditions, capacity constraints for pedestrian flows, roadway capacity limits, maximum flow constraints, shelter capacity limits, initial population per zone, departure preparation time, walking time constraints, and non-negativity of decision variables.

Traffic continuity conditions are given by Equations 13–18. p t i + u t i − v t i + l t i − n t i = p t + 1 i , ∀ t ∈ T ≡ 0 , … , T 1 − 1 , ∀ i ∈ Z ∪ Z E (13) q t i − l t i + n t i − o t i = q t + 1 i , ∀ t ∈ T , ∀ i ∈ Z E (14) q t i − l t i + n t i = q t + 1 i , ∀ t ∈ T , ∀ i ∈ Z E (15) r t i + o t i = r t + 1 i , ∀ t ∈ T , ∀ i ∈ Z E (16) u t i = m i , t S i + m i , t N i + m i , t E i + m i , t W i , ∀ t ∈ T , ∀ i ∈ Z ∪ Z E (17) v t i = m N i , t i + m S i , t i + m W i , t i + m E i , t i , ∀ t ∈ T , ∀ i ∈ Z ∪ Z E (18)

Capacity constraints in road section are expressed by Equations 19, 20. m j , t i ≤ C j , t i , ∀ t ∈ T , ∀ i ∈ Z ∪ Z E , ∀ j ∈ J i (19) u t i ≤ δ i H i − p t i , ∀ t ∈ T , ∀ i ∈ Z ∪ Z E (20)

Maximum flow rate conditions limited by the existence of pedestrians are given by Equations 21–23. ϵ 1 v t i + n t i ≤ p t i , ∀ t ∈ T , ∀ i ∈ Z ∪ Z E (21) l t i + o t i ≤ q t i , ∀ t ∈ T , ∀ i ∈ Z E (22) l t i ≤ q t i , ∀ t ∈ T , ∀ i ∈ Z (23)

Capacity constraints of evacuation shelter and its entry is given by Equations 24, 25. r T 1 i ≤ F i , ∀ i ∈ Z E (24) o t i ≤ E t i , ∀ t ∈ T , ∀ i ∈ Z E (25)

Initial population in each zone is given exogenously by Equation 26. q 0 i = D i , ∀ i ∈ Z ∪ Z E (26)

Preparation time constraints are expressed by Equations 27, 28. l t i = 0 , ∀ t ∈ 0 , … , τ 1 − 1 , ∀ i ∈ Z ∪ Z E (27) o t i = 0 , ∀ t ∈ 0 , … , τ 1 + τ 2 − 1 , ∀ i ∈ Z E (28)

Walking time constraint is as given by Equation 29. p t i ≥ ∑ s = 1 ϵ 1   u t − s i , ∀ t ∈ ϵ 1 , … , T 1 , ∀ i ∈ Z ∪ Z E (29)

Lastly, nonnegative conditions for variables are given by Equation 30. p t i , q t i , u t i , v t i , m j , t i , l t i , n t i , o t h ≥ 0 , ∀ t ∈ 0 , … , T 1 , ∀ i ∈ Z ∪ Z E , ∀ j ∈ J i , ∀ h ∈ Z E (30)

The model also incorporates the congestion effect, where increased pedestrian density reduces inflow to a zone. This is represented using a trapezoidal MFD structure, with key constraints including Equation 21, which indicates the limits on outflow based on existing population; Equation 19, which shows the flow capacity constraints at zone boundaries; and Equation 20, which demonstrates the congestion effects restricting inflow. These collectively describe the trapezoidal-shaped MFD governing outflow volume ( u t i ).

The optimization model developed above seeks to minimize total dynamic risk across the entire region. As a result, individual evacuees may be directed along routes that diverge from those leading to the nearest shelter or that initially approach higher-risk areas. In such cases, evacuees may reject these instructions and instead follow routes they personally perceive as safer, potentially compromising evacuation effectiveness.

To address this issue, it is important to limit movement directions to those more acceptable to evacuees, such as toward the nearest shelter or higher ground. By choosing from among such acceptable directions, evacuation plans can be made more realistic and acceptable. This constraint can be implemented by Equation 31, below. m i , j t = 0 , ∀ t ∈ 1 , 2 , … , T , ∀ i , j ∉ B (31) where B denotes the set of acceptable movement directions. In the empirical analysis, we define and compare three constrained movement scenarios against the unconstrained optimal case (denoted O):

E) Movement allowed only toward the nearest evacuation shelter; if higher ground is closer, only that direction is allowed.

H) Movement allowed toward both the nearest shelter and nearest higher ground.

S) Movement allowed toward all adjacent zones with lower static risk, as well as toward the nearest shelter.

Typical movement options per mesh are: E) one direction; H) one or two directions; and S) two or three directions.

As a result of incorporating constraint (Equation 31), all model variables remain continuous and nonnegative, and both the objective function and constraints are linear. This ensures that the model remains a standard linear programming problem, solvable by general-purpose optimization solvers. In this study, we used Gurobi Optimizer v11.0.3.

The central urban area of Hachinohe City was selected as the target. The tsunami scenario used is based on a Japan Trench earthquake model, with inundation from the first wave expected within 30–60 min of the earthquake. In line with Hachinohe City’s evacuation policy, which excludes bridge crossings due to potential seismic damage, our study focuses on zones between the Mabechi and Niida Rivers (see Figure 2). The analysis included 284 regional meshes: 77 within the inundation area and 107 outside it. These meshes were treated as computational zones.

According to 2020 Census data, the inundation area has a nighttime population of 44,050. This “baseline population distribution” is used in all simulations. We aggregated the maximum capacities of designated shelters by Hachinohe City, and determined the capacity F i for each zone. Figure 2 shows the nighttime population along with the shelter capacities, displayed using concentric circles.

The tsunami disaster probability R t i for each mesh at given time intervals was determined based on the numerical simulation results of inundation depths conducted by [28] using the tsunami analysis code from Nippon Koei, Co., Ltd. Figure 3 illustrates the static risk value (SRV) calculated from the tsunami disaster probabilities R t i for each zone, representing the expected value of tsunami disaster for a single individual remaining in each mesh until the end.

The road lengths, number of road intersections within each mesh, and the number of roads between meshes were obtained from the open data source, OpenStreetMap. From this data, the road capacity H i and the flow capacity C i j for pedestrians between meshes were defined as Equations 32, 33, respectively. H i = R o a d   L e n g t h   w i t h i n   Z o n e 1 + log ⁡ 10 N u m b e r   o f   I n t e r s e c t i o n s   w i t h i n   Z o n e (32) C i j = N u m b e r   o f   R o a d s   b e t w e e n   Z o n e s   i   a n d   j × 40   p e r s o n s   p e r   m i n u t e s (33)

Pedestrian free-flow speed was set at 5 km/h. Preparation time before evacuation was fixed at 15 min, based on data from the Great East Japan Earthquake.

The four movement scenarios tested were: E) Toward the nearest shelter only; H) Toward the nearest shelter and high ground; S) Toward adjacent zones with lower static risk and the nearest shelter; O) Optimal: all physically possible directions allowed. These scenarios are illustrated in Figure 4 (panels E, H, S, and O).

Using the same baseline population, the SRV without evacuation is 2.1 × 10⁴, meaning nearly half of the 44,050 nighttime residents would encounter the tsunami without action. Table 2 summarizes results: O) Optimal: DRV = 4,883 people (11.1%); S) Safe direction: DRV = 5,167 people (11.7%); H) High ground: DRV = 6,278 people (14.3%); E) Shelter only: DRV = 9,259 people (21.0%).

Under O), shelter occupancy is highest (16,552 people), while 24,717 evacuate outside the inundation zone, fewer than in S or H. Total evacuation distance toward safer directions is 32,332 person-km, higher than S (32,068) and H (25,942). However, O) also involves 3,160 person-km of movement toward danger directions, more than the other three scenarios. This suggests that although O) minimizes risk, it may be difficult to implement due to evacuee resistance to moving toward danger.

Scenario S) results in the highest number of people evacuated outside the inundation area and limits dangerous movement to 1,436 person-km, making it a more acceptable option. Its DRV is only 1.06 times that of O), implying that, with prior explanation, this strategy could be implemented effectively.

Scenario E) assumes movement only to the nearest shelter, disregarding capacity constraints. As a result, 9,508 evacuees remain in inundated zones, producing a DRV of 9,259, which is 1.9 times that of O). The shelter occupancy rate is just 39.1%, indicating underutilization of safer, though less proximate, shelters.

Scenario H) outperforms E but is less effective than S. It achieves a DRV 1.28 times higher than O) and a casualty rate of 14.3%. However, it is implementable through basic education encouraging residents to identify and remember the nearest high ground.

Temporal changes in SRV(t) for all four scenarios are shown in Figure 5. Risk begins decreasing at t = 16, when evacuations start. At t = 30, divergence in rates becomes apparent, with E) plateauing first, followed by H) and then S). The final SRV values correspond to the DRV ranking in Table 2.

Figure 6 compares cumulative inter-zone traffic volumes and population on roads per zone. In E, movement is concentrated along fewer paths, resulting in higher congestion. In contrast, H, S, and O show broader traffic dispersal, with increasing diversity of routes from H to O, demonstrating the effectiveness of broader movement allowances in distributing traffic load.

Based on the above results, several implications for pedestrian evacuation planning emerge. In this study area, if no evacuation occurs, nearly half of the nighttime population is expected to be affected by the tsunami. To prevent congestion that could impede evacuation, it is essential to diversify evacuation traffic by encouraging route choices beyond the default direction to the nearest shelter, particularly across spatial scales greater than 500 m.

If evacuees are restricted to moving only toward the nearest evacuation shelter (E), the traffic dispersion effect is insufficient, resulting in a high expected casualty rate of 21%. In contrast, if residents are familiar with multiple safe evacuation directions from each zone (S), they can choose routes that are less prone to congestion. This behavioral flexibility supports traffic dispersion and reduces the expected casualty rate to 11.7%, which is close to the 11.1% achieved by the optimal scenario (O).

Therefore, disaster prevention education should emphasize not only the importance of prompt evacuation but also the existence of multiple safe directions from each zone. Residents should be encouraged to choose directions that avoid traffic concentration to improve safety. These insights should be integrated into public communication strategies and reinforced through regular evacuation drills.