Numerical assessment of coastal multihazard vulnerability in Tokyo Bay

Many bays worldwide are susceptible to coastal hazards such as storm surges, river floods, and tsunamis. Because most previous studies have focused on one or two of the above-mentioned hazards, in this study, we assess coastal vulnerability based on all three hazards. To accommodate the increase in the number of cases in multihazard analysis, an efficient method based on an estimated overflow volume without computing for inundation is proposed. Subsequently, the method is validated via a comparison with inundation simulation. It is shown that when the free overflow is dominant, the result yielded by the method is consistent with that of the inundation simulation. Using Tokyo Bay as the study area, an efficient method is applied to multihazard vulnerability assessment. By comparing the overflow volume maps and maximum anomaly distribution along the coast for four types of hazards, we investigate the characteristics of different types of hazards and identify the differences between single and multiple hazards. Furthermore, we compare the differences between superposing and concurrent computation methods for multiple hazards. It is discovered that the linear superposing method tends to overestimate the total water elevation in coastal regions; however, in the coast, where the superposing method underestimates multihazard anomalies, dike upgrades must be considered.


Introduction
The Tokyo Bay area is a major political and economic center in Japan with a large population; however, it has been exposed to several types of coastal hazards, including typhoons, storm surges, river floods, and earthquake tsunamis. Among the historical storm surge events that occurred in Tokyo Bay, the typhoon in the sixth year of the Taisho era (1917) resulted in one of the largest storm surges, which exceeded the coastal defense system and resulted in extensive damage to human lives and properties (Miyazaki 2003). In recent years, Typhoon Hagibis struck Japan on October 12, 2019, passed through the eastern part 1 3 of Honshu Island, and caused serious damage, including 93 deaths and more than 87,000 damaged houses (NHK WORLD-JAPAN News 2019). Typhoon Faxai passed through Tokyo Bay on September 9, 2019, causing at least three deaths and approximately 40 injuries (The Japan Times 2019). Although the major coastal hazard in the head of Tokyo Bay is storm surges, tsunamis have also occurred in the bay. In the 2011 Tohoku earthquake tsunami, part of the coast of Tokyo Bay was damaged severely, affecting seaweed farming in Futtsu Cape and the northern part of the Banzu tidal flat; furthermore, inundation was observed in some coastal areas (Sasaki et al. 2012). In addition to storm surges and tsunamis, river floods have become major natural disasters that occur during the summer in Japan. According to the Disaster Management Report (Cabinet Office of the Japanese Government, 2015), the occurrence of heavy downpours is predicted to increase in the future. Although the Tokyo Bay area is generally well protected by coastal structures against storm surges, tsunamis, and river floods, the combination of increased sea level and increased typhoon intensity is likely to render this region more vulnerable to climate change in the future (Hoshino et al. 2012(Hoshino et al. , 2016. A hazard refers to a source or situation that may result in human injury or ill-health, property damage, environmental damage, or a combination of these (Wisner et al. 2004). Disasters are often described as a consequence of the combination of exposure to a hazard and insufficient capacity to manage impending negative consequences (Blaikie et al. 2014;Hossain and Paul 2017). Vulnerability is typically portrayed as susceptibility to harm, and the key parameters of vulnerability are the stress to which a system is exposed, its sensitivity, and its adaptive capacity (Adger 2006). Numerous studies have emphasized the technical aspects of coastal vulnerability, such as the numerical modeling of cyclone effects (e.g., Islam and Peterson 2008;Tasnim et al. 2015), as well as physical and socioeconomic aspect assessments of vulnerability (e.g., Cutter et al. 2008;Hossain 2015). Only a few studies have been performed to identify the vulnerability of Tokyo Bay to coastal hazards. Matsuda (2013) verified Tokyo's vulnerability to natural disasters. They discovered that artificial changes in natural conditions, including river courses and withdrawal of groundwater, affect natural disasters and countermeasures in the lowlands of Tokyo. Takabatake and Shibayama (2012) investigated the risk of storm surges and tsunamis in Tokyo Port via numerical simulations of storm surge, tsunami, and overland inundation. Hirano et al. (2014) investigated the maximum possible typhoon conditions based on the Intergovernmental Panel on Climate Change Special Report on Emission Scenarios A1B using an atmospheric-ocean-wave-coupled model and showed that sea-land interactions and river flows may significantly affect the depth and extent of inland inundation. The Tokyo Metropolitan Government (2018) and Chiba Prefecture Government (2018) attempted to predict the maximum inundation induced by the worst expected typhoon storm surge. However, all these studies focused on single coastal hazards, such as storm surges or tsunamis. Recently, on October 12, 2019, a typhoon (No. 19) passed through the Greater Tokyo area, which caused a M w 5.7 earthquake to occur near Tokyo Bay (https:// www. eagle news. ph/5-7magni tude-quake-hits-japan/). Although no tsunami was observed, it would be meaningful to identify the multihazard vulnerability to achieve better disaster mitigation management, considering that the probability of earthquake occurrence in Tokai, Tokai-Tonankai, or Tokyo Inland has been increasing as time progresses since its previous occurrence.
Motivated by the accuracy and extensive application of numerical models for coastal hazard management, and that the joint effect of three hazards (storm surge, tsunami, and river flood) has not been investigated, a numerical assessment of inundation vulnerability caused by these hazards was conducted in this study. Inundation simulation is a widely used and straightforward method to assess coastal vulnerability. However, it is computationally expensive, and considering the increase in the number of cases in multiple hazard analyses, an efficient method must be devised to identify the overall vulnerability and to filter out representative scenarios for a detailed analysis. Hence, an efficient method using an estimated overflow volume without computing for inundation is proposed herein. The method is validated based on a comparison with an inundation simulation. Subsequently, based on Tokyo Bay as the study area ( Fig. 1), an efficient method is applied to a multihazard vulnerability assessment. The worst multihazard case and resultant vulnerability are identified. The difference between single-hazard and multihazard vulnerabilities, including moderate multiple hazards and worst single hazards, is discussed.  (Chen et al. 2003) was adopted for the numerical simulation. FVCOM is a prognostic, unstructured grid, finite-volume, freesurface, and three-dimensional ocean model. The model was used to solve momentum, continuity, temperature, and salinity equations. In the present computation, the ratio of the internal to external time steps was set to 10 for all applications. The initial and boundary conditions for the different simulation cases are listed in Table 1. A Manning roughness coefficient of n = 0.025 was used in the model to describe the ocean bottom friction. A high-resolution grid is necessary for the complicated port and channel configurations. To reduce the boundary effects, a wide area was selected as the computational domain (Fig. 2). An unstructured mesh system was constructed using Surface Water Modeling System 11.1 (Aquaveo). Coastline data from the Global Self-consistent, Hierarchical, Highresolution Geography Database were utilized for all areas except Tokyo Bay, whose data were extracted from Google Earth during mesh generation. The mesh size varied from 50 km at the open boundary to 50 m on the inner side of Tokyo Bay, and the typical mesh resolution in the Tokai-Tonankai area (the tsunami source area) was 1 km. The numbers of mesh nodes and triangular elements in the computational domain were 141,003 and 274,482, respectively. A Digital Elevation Model (DEM) with 5 m resolution was adopted from the Geospatial Information Authority of Japan. The bathymetry data were data based on the ETOPO-1 Global Relief model and J-EEG500 of the Japan Oceanographic Data Center (a 500 m resolution mesh dataset). Bathymetry data of 1 arc-minute resolution provided by ETOPO were used for the entire domain, except for the area around Tokyo Bay and the inner Tokyo Bay area. Bathymetry data with a resolution of 500 m were used for Fig. 2 Computational mesh the area around Tokyo Bay, whereas those with a 50 m resolution were interpolated to the entire Tokyo Bay mesh.

Storm surge simulation
For the storm surge simulation, the model was driven by wind and atmospheric pressure. The typhoon atmospheric pressure field was computed using the Myers formula (Myers 1957), which is expressed as where P(r) is the pressure at a radial distance r from the typhoon center, P c (hPa) the typhoon central pressure, P 0 (= 1013.25 hPa) the ambient or environmental pressure, r the distance from the computational mesh node to the typhoon center, and r max (km) the maximum wind speed radius.
Based on a previous study (Liu and Sasaki 2019), the Mitsuta-Fujii model was used to compute the wind field, as presented in Eq. (2), and the pressure estimated using Eq. (1) was applied to the Mitsuta-Fujii model.
where U w is the total wind vector, U w1 the wind vector induced by the rotating component, U w2 the wind vector induced by the moving component, U w1 (r) the magnitude of U w1 , P(r) the pressure field calculated using the Myers formula, C 1 and C 2 dimensionless coefficients ranging from 0.6 to 0.75, f the Coriolis parameter, r the distance from the typhoon center, a the atmospheric density, and V the moving velocity vector of the cyclone center.
The earthquake moment magnitude ( M w ) and seismic moment ( M 0 ) are correlated as follows (Hanks and Kanamori 1979): The length L F (km) and width W F (km) of the earthquake fault were determined using Eqs. (4) and (5), respectively, for a specified earthquake (Leonard 2010;Samarasekara et al. 2017).
The average net slip (Wells and Coppersmith 1994;Anderson et al. 1996) was calculated as follows: where = 3.4 × 10 10 (N∕m 2 ) is the elastic modulus of the crust.

River flood simulation
On the inner side of Tokyo Bay, two major rivers (Arakawa River and Edo River) are connected. Owing to availability, only the daily discharge data of the Edo River from 1955 to 2015 were acquired from the Edo River Office of the Ministry of Land, Infrastructure, Transport, and Tourism. Figure 4 shows the annual maximum discharge (a) and return period (b). The discharge of Arakawa River was estimated using the mean discharges of Arakawa River (87.5 m 3 /s) and Edo River (82.3 m 3 /s).
For the Arakawa, Nakagawa, and Edo rivers, the flow discharge was 2800, 1000, and 2600 m 3 ∕s , respectively, which are approximately 50-year return period discharges, as shown in Fig. 4b.

Vulnerability measuring and land overflow volume
Measuring vulnerability is difficult, particularly in quantitative hazard studies, and depends primarily on the selection of indicators based on the data. In addition, current vulnerability (4) log 10 M 0 = 2.5 log 10 L F + 7.96 (5) log 10 W F = 0.667 log 10 L F + 1.24 Okada (1985) formulation assessment studies typically do not consider all the possible pressures that are associated with the diminished capacity and decreased resilience of communities (David and Colin 2000).

Fig. 3 Seismic parameters in
In most coastal vulnerability studies (Mani Murali et al. 2013;Benchekroun et al. 2015), the problem of tsunami vulnerability is addressed using empirical fragility functions and damage curves. Tsunami fragility is defined as the structural damage probability or fatality ratio with regard to the hydrodynamic features of the tsunami inundation flow, such as the inundation depth, current velocity, and hydrodynamic force. Tsunami fragility curves developed using satellite remote sensing and numerical modeling represent the overall damage in an area that exhibits many types of tsunami features (Suppasri et al. 2012). However, inundation simulation is computationally expensive, considering the increase in the number of cases in multihazard assessments. Hence, an efficient approach using overflow volume to measure coastal vulnerability is introduced. In the present study, we measured multihazard vulnerability by combining the overflow volume and maximum anomaly distribution, and four types of hazard scenarios were prepared for the vulnerability assessment (see Table 2). Figure 5 shows dike overflow process. Between Sections I-I and II-II, the Bernoulli's equation can be written as follows: where P 1 and P 2 represent atmospheric pressures (hPa), the water density ( kg∕m 3 ), g the gravitational acceleration ( m∕s 2 ), v 1 the vertically averaged horizontal velocity in Section I-I (m/s), H 1 the water surface elevation at section I-I (m), v c the overflow velocity in . 4 a Annual maximum discharge, b return period for Edo River Section II-II (m/s), H 2 the elevation of the sea dike or river bank (m), and h c the overflow depth above the sea dike or river bank (m). Rearranging Eq. (7) and then combining it with the critical flow condition, the overflow depth and overflow velocity in Section II-II can be obtained as follows: then, the overflow volume at this node is estimated by: where V nodal is the nodal overflow volume at a coastline node, m 2 ; i is the time step number, N out is the number of output steps, Δt out is the output time interval (s), and the other parameters have been introduced.
To measure the vulnerability, the overflow volume is calculated as follows: where V dike (m 2 ) is the total overflow volume of the coastline, j the number of coastline nodes, and m c the total number of coastline nodes. If H 2 > H 1 , then no overflow occurs at the node, and V nodal equals 0.

Numerical experiment setup
The overflow estimation was conducted using only the propagation method. To validate this method, the overflow volume estimated by the propagation method was compared with the results of inundation computation. The computational meshes for the propagation method and inundation computations are shown in Fig. 6. A uniform inflow was imposed at the left-side boundary. The computational domain for the inundation computation (see Fig. 6) was partitioned into two regions: propagation and inundation areas. When the experiment was started, the water first entered and propagated in the propagation area. As the water level increased, the water began to overflow the coastal dike (dashed red line). Subsequently, the inundation depth in the inundation area increased gradually. Finally, the inundation depth exceeded the dike height. During this process, the estimated dike overflow volume was compared with the inundation volume in the inundation area. As the dike height determines the time when overflow occurs (the inflow is constant), a series of dike height values was set. The dike height was varied from 0.5 to 1.5 m at an interval of 0.5 m.

Verification
In Fig. 7, the estimated dike overflow volumes and inundation volumes are shown in the left panels, whereas time series of water elevations at four observation points are shown in the right panels, for cases involving dike heights of (a) 0.5, (b) 1.0, and (c) 1.5 m.
As the dike height changed, the consistency between the estimated dike overflow volume and inundation volume varied. For instance, when the dike height was 1.5 m (see the left panel of Fig. 7c), from 900 to 1250 s, the estimated dike overflow volume was consistent with the inundation volume. However, beyond 1250 s, a significant discrepancy was observed between them. As shown in Figs. 7a, b, the consistent periods were from 280 to 600 s and from 600 to 1100 s, respectively. This phenomenon can be explained based on the time-series water elevation at four observation points ( P 1 , P 2 , P 3 , and P 4 ). As shown in the right panel of Fig. 7c, the elevation of three observation points (P 1 , P 3 , and P 4 ) did not change from 0 to 900 s, i.e., it was always equal to the elevation from bathymetry data. For example, the green line ( P 1 ) representing 1.5 m corresponded to the dike height; the blue line ( P 3 ) representing 0.2 m corresponded to the bathymetry data in the inundation area. This was because during this period, water propagated only in the propagation area and it did not reach the dike crown. However, the situation of P 2 was different; water began to flow from 0 m, and immediately after a short time period (approximately 80 s), the elevation of P 2 (yellow color line) began to increase because the inflow required only a short time to reach P 2 from the left boundary. At 900 s, P 2 reached the dike crown height, and P 3 and P 4 began to elevate. At 1250 s, P 3 and P 4 reached the same elevation as the other two points, and the overflow was completely submerged. Subsequently, the overflow volume was not consistent with the inundation volume.
Therefore, this approach can be used under free-overflow conditions. In the present study, we assume that the inundation depth does not exceed the sum of the dike height and overflow depth, which implies that the free-overflow type is dominant. Fig. 6 Mesh grid for a estimation of overflow volume using only propagation method, b inundation computation (Uniform inflow of 2 m 3 ∕s per unit length is given at the left-side boundary; P 1 , P 2 , P 3 and P 4 are assumed as 4 observation floats, they can only move vertically with the water level change at its location, P 1 is on the dike, the distance of P 2 to dike, P 3 to dike, and P 3 to P 4 is 10, 5, and 10 m, respectively) and bathymetry configuration for dike height 0.5 m 1 3

Multihazard scenarios
The Tokyo Metropolitan Government (2018)   investigate the worst storm surge condition, several parallel tracks with the same intensity were assumed in each typhoon course, as shown in Fig. 8.
The intensity of the assumed typhoon cases was fixed. The central pressure was 910 hPa, the radius of the maximum wind was 75 km, and the forward translational wind speed was 30 m/s. According to the JMA best track data archive, the selected typhoon intensity was based on the 1000-5000-year return period. As shown in Table 3, the maximum storm surge anomaly in the inner side of Tokyo Bay was greater than that in the southern part of the bay. In the multihazard assessment section, typhoon tracks T1S2, T2S3, T3S3, T4S1, and T4S2 were selected as the scenario components because these five cases generated the most significant anomalies in the northern part of the bay.
Currently, tsunamis induced by the Tokai-Tonankai-Nankai-type earthquake and an inland earthquake that may occur near Tokyo metropolitan area are major concerns in the Tokyo Bay (Sasaki et al. 2012). Particularly, the Tokai-Tonankai-type earthquake might occur in the near future owing to its short return period of approximately 120 years. To encompass the earthquake cases as much as possible, the earthquake fault area was calculated using the fixed rupture length and width of rectangular subfaults, i.e., 50 and 25 km, respectively; hence, Zone 1 (the potential maximum fault area of Tokai-Tonankai-type earthquake indicated by red mesh) and Zone 2 (the potential maximum fault area of Tokyo-inland-type earthquake indicated by blue mesh) were partitioned into 50 km × 25 km cells, as shown in Fig. 9a.
The fault rupture length and width ( L F and W F , respectively) calculated using Eq. (4) and (5), respectively, were rounded off to be multiples of 50 and 25 km, respectively (e.g., Samarasekara et al. 2017). Figure 9b shows an example of the construction of fault areas for a M w 8.0 earthquake. The total number N of possible rectangular tsunami faults for an earthquake of a certain magnitude in Zone 1 can be estimated as follows:  According to the Disaster Management Report (Cabinet Office of the Japanese Government 2015), the probability of a large earthquake (i.e., M w 8-9) occurring within 30 years in the Nankai trough and a M w 7 earthquake occurring in the southern Kanto area within 30 years is approximately 70%. By reviewing past earthquake records, it was discovered that several major earthquake tsunamis have occurred in the Tokyo Bay, including the 1703 Genroku Kanto earthquake ( M w 8.2 ) tsunami, the 1854 Ansei Tokai earthquake ( M w 8.4 ) tsunami, the 1923 Taisho Kanto earthquake ( M w 7.9 ) tsunami, the 1960 Chilean earthquake tsunami, and the 2011 Tohoku earthquake tsunami (Sasaki et al. 2012). The magnitude of devastating tsunamigenic earthquakes in the Tokai-Tonankai area and Tokyo Bay was around 8.4 and 8.0, respectively. Therefore, M w 8.4 and M w 8.0 earthquakes were selected as the scenarios for Zones 1 and 2, respectively. The number of possible tsunami fault scenarios for each earthquake magnitude was determined using Eq. (12) and are summarized in Table 4.
(15) L F = 179.47 km, W F = 55.29 km . 9 a Potential maximum fault areas for Tokai-Tonankai type earthquake (red mesh), and Tokyo inland type earthquake (blue mesh), b Fault areas for a M w 8.0 earthquake in Zone 1 ( N L F and N W F are the total numbers of the meshes covering the length and width of the potential maximum fault area, N L and N W are the total numbers of sub-faults covering the length L F and width W F for a specific earthquake fault. Green rectangle covering four mesh cells in length and two mesh cells in width represents one possible sub-fault for a M W 8.0 earthquake, and by shifting the green rectangle every time one mesh cell, from green rectangle to yellow rectangle, and then from yellow rectangle to blue rectangle, etc. Then, all possible sub-faults can be obtained.)

3
The actual earthquake slips may not be uniform throughout the fault area. To simplify the preparation procedure in the present study, a uniform slip throughout the entire rupture area was assumed. Hence, the main parameter to be modified was the magnitude of the earthquake. For tsunamigenic earthquakes that occurred in Zone 1, the fault parameters of the Keicho earthquake (Kanagawa Prefecture Government Report 2012; Takabatake and Shibayama 2012) were selected, as shown in Table 5, whereas for Zone 2, the fault parameters of the Tokyo inland earthquake were used based on the materials provided by the Tokyo Metropolitan Government Disaster Mitigation Report (2012). The values of longitude, latitude, fault length, fault width, and slip extent were determined based on the earthquake moment magnitude and the scenarios.
Considering the effect of time when tsunami occurs during typhoon movement, the multihazard scenarios were classified into two groups: Groups 1 and 2 (see Table 6). Group 1 (1-T1S2-14 to 1-T4S2-20) and Group 2 (2-T1S2-20 to 2-T4S2-20) comprise tsunami scenarios in Zones 1 and 2, respectively. In Group 1, because the fault area is distant from Tokyo Bay, it was assumed that tsunamis occurred when the typhoon center propagated to different locations during its approach to Tokyo Bay. By contrast, for Group 2, we only considered the situation where tsunami occurred when the typhoon center was located in Tokyo Bay; this is because under this situation, the multihazard anomaly was the greatest. Figure 10 shows the typhoon cases and initial condition for a M w 8.4 earthquake tsunami in Zone 1 and a M w 8.0 earthquake tsunami in Zone 2. Detailed information regarding all combinations and the typhoon center locations when earthquake occurred are provided in Table 6. Figure 11 shows the spatial variations in the sea surface elevation for the worst cases of storm surge, river flood, and the Tokai-Tonankai earthquake tsunami. To present the similarities and differences among the four types of hazard scenarios more clearly, the maximum anomaly distributions and overflow maps of the four types of hazards are shown simultaneously in Fig. 12, and the overflow volume was categorized into three levels, as described in Table 7.
Comparing the maximum anomaly distributions shown in Fig. 12a and b, it can be concluded that the river discharge enlarged the anomaly in the river channels. After combining the Tokai-Tonankai earthquake tsunami, more significant anomalies occurred at the head of the bay. However, the multihazard case combined with the Tokyo inland earthquake tsunami did not result in such a phenomenon.
Comparing the overflow volume map and the maximum anomaly distribution in each subfigure (e.g., Fig. 12a), it was discovered that the overflow volume along the coasts was similar to that of the maximum anomaly distribution, (particularly for the storm  Fig. 12a), the multihazard case combining storm surge and river flood in Fig. 12b, and the multihazard case combining storm surge, Tokyo inland earthquake tsunami, and river flood in Fig. 12d. This is because, for the former two cases, both significant anomalies and overflow volumes occurred in the head of the bay and rivers. However, for the multihazard case combining storm surge, the Tokai-Tonankai earthquake tsunami, and the river flood shown in Fig. 12c, the actual overflow volume condition cannot revealed if only the maximum anomaly distribution is considered. For

Discussion
Because investigations pertaining to the concurrence of storm surges and tsunamis are scarce, we provide a discussion pertaining to it in this section. To reveal the difference between single hazard (tsunami or storm surge) and multihazard cases, a comparison between a multihazard case (hereinafter, concurrent case or concurrent multihazard case) and a case (hereinafter, superposing case) linearly superposing each single hazard case was conducted. The anomaly difference along the coast between the concurrent multihazard case and the superposing case is shown in Fig. 13. To present the difference caused by the typhoon center locations when the tsunami was initiated, we analyzed multihazard cases that were based on the same typhoon. The difference value was obtained by subtracting the anomaly of the concurrent multihazard case from that of the superposing case. It was discovered that the difference varied spatially and significantly and that it was more significant in ports and river channels than in other coasts. In most coastal areas, the anomaly of the superposing case was greater than that of the concurrent case; however, in a few locations, the concurrent case resulted in greater anomalies than in the superposing case. Compared with the concurrent computation, the superposing method overestimated the total water level for most coastlines, particularly the Chiba Port, Yokohama Port, and water channels in Koto Ward, where the difference exceeded 0.5 m (shown in red). Although studies regarding the interaction between storm surges and tsunamis have not been reported, similar phenomena have been observed in many studies (e.g., Prandle and Wolf 1978;Johns et al. 1985;Sinha et al. 2008) pertaining to tide and storm surge interactions or tide and tsunami interactions, all of which  1 3 indicate that such a linear superposing method is not strictly correct and typically tends to overestimate the total coastal sea level elevation. From the perspective of security, the superposing method is acceptable; however, considering the economy of the design and Fig. 13 Anomaly difference of the concurrent multihazard case and superposing case: a to g Typhoon location at 14-20 when tsunami initiated construction of coastal defense infrastructure, it may not be appropriate. The merits and demerits of superposing and concurrent computation methods are listed in Table 8. Compared with studies pertaining to the tsunami-tide interaction (e.g., Lee and Kaneko 2015;Kowalik et al. 2006), tides in the open sea do not significantly affect tsunami wave propagation; as such, a linear superposing method can be applied when varying sea levels  are involved. In our study, storm surges elevated the water level in the bay, including that of coastal areas when typhoon was passing, which may weaken tsunami wave shoaling compared with the only-tsunami case (beginning from still water). Subsequently, a comparison between the multihazard case (moderate storm surge and moderate tsunami) and the single-hazard case (the worst storm surge) was conducted. As shown in Fig. 14, in the worst storm surge case, the anomalies were greater than those in the moderate multihazard case (see Fig. 14d), particularly in the north of the bay; meanwhile, for the multihazard case (see Fig. 14c), the difference in anomalies in the north and south was insignificant.
Although the moderate multihazard event may cause less damage than the worst single storm surge, the damage risk of single and multihazard events on local coastal areas may vary owing to the different spatial patterns. In Tokyo Bay, port areas such as Chiba Port and Yokohama Port, low-lying areas such as narrow channels in Koto Ward and coasts such as those in Funabashi, and the upstream areas of the investigated rivers in the northern part of the bay must be prioritized because they are more susceptible to damage compared with other places. To prevent inundation and reduce losses caused by multiple hazards, the dike heights in areas where the superposing method underestimates the multihazard anomalies must be upgraded. Early-warning and hazard prediction systems are useful for disaster mitigation. However, the current system used in the study region may not be able to directly provide the required flood information when multihazard events occur. Hence, the capability of multihazard computation may need to be incorporated to obtain more accurate hazard information. For Tokyo Bay, the current dike system can defend against concurrent moderate multihazards comprising moderate storm surges, moderate tsunamis, and river floods; however, for worst concurrent multihazards, policymakers may need to reconsider the design code or construction standard of the defending system.

Conclusions
Each year, many bays worldwide pose the risks of significant economic damage or loss of human lives because of coastal hazards, such as storm surges, river floods, and tsunamis. Significant efforts have been expended by coastal researchers to identify the risks posed by coastal hazards and the underlying mechanisms. As no studies have focused on the concurrent hazards of storm surges, river floods, and tsunamis, although they nearly occurred during the recent typhoon (No. 19) on October 12, 2019 in Tokyo Bay, the vulnerability of multihazards must be investigated owing to the importance of the bay area. Herein, an efficient method was proposed using an estimated overflow volume without computing inundation, which was validated via a comparison with an inundation simulation. It was indicated that when the free dike overflow was dominant, the proposed method yielded results that were consistent with results of the inundation simulation. Using Tokyo Bay as the study area, an efficient method was applied to the multihazard vulnerability assessment, and a multihazard vulnerability map was generated. By comparing the vulnerable coasts under different types of hazards, it was discovered that river flooding increased the overflow extent in the riversides, and the entire bay exhibited more overflow when storm surges, tsunamis, and river floods were combined. Furthermore, it was discovered that tsunamis deteriorated the overflow extent of the riverside, including upstream. The anomaly difference between the superposing and concurrent cases (storm surge and tsunami occurring simultaneously) was more significant in ports and river channels than in other coastal areas, which demonstrated the significant interaction between storm surges and tsunamis in those areas. Worst storm surges would cause greater anomalies than moderate (1 m storm surge anomaly or tsunami) multiple hazards, in which the return periods of these two types of hazards were identical. However, the damage risk of multiple and single hazards in local coastal areas might vary owing to their different spatial patterns. The superposition method may not be appropriate to identify multihazard risks compared with concurrent computations of coastal regions because it tends to overestimate the water elevation of multihazards in coasts.
his team for providing the FVCOM model. The computation was carried out using the computer resource offered under the category of General Projects by Research Institute for Information Technology, Kyushu University.

Author contributions
The design of the present study was originally developed by FL under the supervision of JS. Numerical simulation and analysis were performed by FL. The first draft of the manuscript was prepared by FL with the revision by JS. JC and YW contributed to the mesh generation as well as the design of the multihazard scenarios. JC conducted part of the analysis in river flood simulation. All the authors made comments, which were accommodated by FL. All authors confirmed and approved the final manuscript.