. Multi-hazard loss estimation for shaking and tsunami using stochastic rupture sources. International Journal of Disaster Risk Reduction

This study develops a probabilistic multi-hazard loss estimation methodology for coastal communities which are exposed to cascading shaking-tsunami hazards caused by o ﬀ shore mega-thrust subduction earthquakes. The method captures the common source e ﬀ ects for simulating shaking and tsunami hazards and risks. It facilitates the quantitative evaluation of multi-hazard loss for coastal communities by accounting for uncertainties associated with loss estimation based on novel stochastic earthquake source modeling and state-of-the-art shaking-tsunami fragility modeling. A case study is set up to illustrate an application of the developed method to Miyagi Prefecture by focusing on possible large seismic events in the Tohoku region of Japan. The quantitative multi-hazard impact assessment results serve as e ﬀ ective means to make decisions regarding shaking-tsunami disaster risk reduction.


Introduction
A mega-thrust subduction earthquake, originated at major plate boundaries, can cause tremendous damage and loss to society, infrastructure, and assets, as revealed during the 2004 Indian Ocean, 2010 Maule Chile, and 2011 Tohoku Japan events.Remarkable characteristics of the catastrophic impact are that: (i) multiple hazards are triggered by a single event, affecting the population and built environment in cascades, and (ii) spatial scale of the impact can be very large, evolving over time.
The main hazard processes due to off-shore large subduction earthquakes are ground shaking and tsunami [10,16], while other secondary hazards, such as aftershock, liquefaction, landslide, and fire, can be significant [6].Moreover, earthquake-triggered hazards and their impact scale differently in terms of earthquake magnitude and are highly dependent on earthquake rupture processes [21].For instance, ground shaking will cause damage to structures that are distributed widely in a region.On the other hand, tsunami damage is localized along coastal areas and will increase drastically with earthquake magnitude [19].These features result in spatiotemporally correlated financial losses, which are challenging to accommodate and diversify through conventional risk sharing mechanisms, such as insurance and reinsurance [17].For seismic risk management, it is important to capture the multi-hazard loss generation mechanism and its uncertainty in assessing regional earthquake impact because the damage and loss patterns for shaking and tsunami are different.Thus emergency responses may need to be organized and operated differently.
An integrated multi-risk assessment due to cascading hazards is a major challenge [11,29,52].Several recent studies proposed multi-risk approaches by integrating existing impact assessment methodologies in different ways [31,35,39,40,53].Three essential components of the multi-risk framework are: (i) joint probability of hazards, (ii) timevariant multi-hazard dependent vulnerability, and (iii) combination of losses from different hazards in a coherent manner.
Multi-hazard modeling and impact assessment for ground shaking and tsunami have been attempted in literature (e.g.[7,33,48]).In the context of the present study, it is important to define the term 'multihazard' clearly.It is attributed to common source features of strong shaking and tsunami that are captured through physical parameters; simulations of such physical hazard processes are conducted sequentially by maintaining the source dependency.From this perspective, studies that perform two disjoint analyses for shaking and tsunami are not focused upon.Maeda et al. [33] developed a physics-based simulation method of generating seismic and tsunami waves to model cascading shaking-tsunami hazards.A limitation of the proposed method was its high computational cost, making it practically unfeasible to conduct numerous runs of inland tsunami inundation simulations at high resolutions.De Risi and Goda [7] developed a multi-hazard methodology for conducting probabilistic hazard analyses for shaking and tsunami by taking into account stochastic earthquake ruptures as trigger mechanism of strong ground motions and tsunami waves.Their approach was innovative in that a wide range of seismologically plausible earthquake sources was represented by adopting stochastic synthesis of slip distribution [17] and probabilistic source scaling relationships for mega-thrust subduction earthquakes [20].Recently, Park et al. [48] have extended a logic-tree-based probabilistic seismic hazard analysis into a probabilistic seismic and tsunami hazard analysis for the Cascadia subduction zone in the Pacific Northwest.In that approach, uncertainties of earthquake scenarios were captured by branches of the logic tree and were propagated through ground-motion and tsunami simulations.Although the previous investigations have expanded the current capability to carry out probabilistic multi-hazard assessments for shaking and tsunami, their combined impact on the built environment has not been assessed.
It is important to point out that multi-hazard mapping of the relevant parameters at a selected return period level, such as peak ground acceleration (PGA) map for shaking and inundation depth map for tsunami evaluated at a 500-year return period level, is useful, but these two maps will not inform users of possible combined impact in terms of damage and loss, which are more relevant for risk management purposes.Currently, a probabilistic method to evaluate the combined shaking-tsunami loss at regional scale, focusing upon a portfolio of buildings, is lacking, although such methods are available for shaking damage alone (e.g.[3,14,55]) and for tsunami damage alone (e.g.[21,47]).
This study develops a probabilistic multi-hazard loss estimation framework for shaking and tsunami by adopting the stochastic rupture approach [17,20] to capture the common source effects.The proposed method essentially builds upon the multi-hazard method of De Risi and Goda [7] for evaluating seismic and tsunami hazards concurrently.Earthquake shaking simulations are performed using ground motion prediction equations (GMPEs) together with applicable spatial correlation models to represent realistic spatial distribution of seismic effects at building locations [14], while tsunami inundation simulations are conducted by evaluating nonlinear shallow water equations with initial boundary conditions caused by the earthquake rupture.The adopted methods facilitate the realistic representation of spatially correlated hazard values.To enable quantitative evaluation of shaking-and tsunami-triggered damage and loss to structures, seismic and tsunami fragility models and building cost models are integrated.The formulation of the multi-hazard impact assessment is equivalent to that of performance-based earthquake engineering methodology (e.g.[5,24,50]).In this context, the developed multi-hazard loss estimation method in this study can be regarded as the performance-based earthquake-tsunami engineering methodology, particularly applicable to multiple buildings at regional level.The developed framework for the multi-hazard shaking-tsunami loss estimation is demonstrated through a case study for Miyagi Prefecture in Japan, which was devastated by the 2011 Tohoku tsunami [10].The case study is illustrative, rather than comprehensive evaluation of multi-hazard risks, and considers seismic sources off the Tohoku region only; it does not take into account crustal seismic sources for shaking-related loss nor distant tsunami sources (e.g.Chile and Cascadia) for tsunami-related loss.
The paper is organized as follow.Section 2.1 presents the formulation of the multi-hazard loss model.Relevant model components are introduced in Section 2.2.In Section 3, loss estimation results for four coastal cities and towns in Miyagi Prefecture are discussed by deriving multi-hazard loss exceedance curves.The main purposes of the analyses are to demonstrate the novel loss estimation methodology through a realistic case study and to highlight the multi-hazard earthquake impact with regards to conventional single-hazard approaches.The new tool is particularly useful for investigating individual contributions from different hazards to total economic loss.Finally, Section 4 draws key conclusions from this study and mentions further applications of the developed loss estimation framework in earthquake disaster risk management.

Formulation
A general formulation of the multi-hazard loss estimation for shaking and tsunami follows a so-called Pacific Earthquake Engineering Research (PEER) equation [5]: In Eq. ( 1), ν L (L ≥ l) is the mean annual occurrence rate that the total loss L for a portfolio of buildings caused by shaking and tsunami exceeds a certain loss threshold l.The random variables M, RS, IM, and DS correspond to earthquake magnitude, rupture source parameters (e.g.fault geometry and slip distribution), shaking-tsunami intensity measures (e.g.peak ground velocity, PGV, for shaking and inundation depth for tsunami), and building damage states (e.g.collapse and complete damage), respectively.The integration in Eq. ( 1) should be performed over all random variables.λ Mmin is the mean annual occurrence rate of earthquakes with M ≥ M min that cause major economic loss to the building portfolio, while f M is the conditional probability

|
is the shaking-tsunami fragility function, which predicts the probability of incurring a particular DS for given IM.P(L ≥ l|ds) is the damage-loss function given DS.It is noted that Eq. ( 1) describes the shaking-tsunami loss due to a single near-source region (e.g.off-shore Tohoku region for interface subduction earthquakes).When multiple sources are of interest, the formulation can be extended by including loss contributions from different source regions.Such an extended formulation can be simplified by neglecting shaking-related contributions due to distant sources (e.g.shaking damage due to earthquakes at > 300 km distances) and tsunami-related contributions due to non-tsunamigenic sources [8].
To evaluate the risk equation, numerous stochastic source scenarios are generated for several earthquake magnitudes, and synthesized rupture sources for a given magnitude range are used to evaluate the conditional loss exceedance function.Following the stochastic rupture approach combined with discrete representation of f M [7], Eq. ( 1) can be expressed as: where p Mk is the probability mass for a given magnitude range represented by m k (k = 1,…,n M ), and the conditional loss exceedance function P(L ≥ l|m k ) for the magnitude range m k is given by: Assuming that hazard intensities for shaking and tsunami are conditionally independent, f IM RS | can be expressed as: where the subscripts S and T for IM are used for shaking and tsunami, respectively.Practically, this means that for a given rupture source model, for instance, PGV and maximum inundation depth at a site of interest can be evaluated independently.Then, Eq. ( 5) can be obtained by substituting Eq. ( 4) into Eq.(3): It is noted that in Eq. ( 5), the structural damage state is dependent on both shaking and tsunami intensity values.More details on the fragility models is explained in Section 2.2.6., where the practical implementation is shown.
In the proposed method, Eq. ( 5) is evaluated numerically through Monte Carlo simulations (MCS).In such a case, P(L ≥ l|m k ) can be obtained as follows: where n MCS is the number of Monte Carlo sampling for stochastic realizations of RS, IM S , IM T , and DS, and I(•) is the indicator function and equals 1 when the loss for the j-th realization is greater than or equal to l and 0 otherwise.
It is important to highlight that the loss variable L is defined for a building portfolio which is distributed spatially along coast and is influenced by both shaking and tsunami intensity parameters IM S and IM T (which are dependent on the common earthquake rupture process RS and represent spatially-correlated random fields).For a given building portfolio consisting of n bldg structures, the total loss due to the j-th stochastic realization can be determined as follows: Fig. 1.Multi-hazard loss estimation procedure for shaking and tsunami.
where L i,j , L Si,j , and L Ti,j are the losses due to the combined effects of shaking-tsunami damage, shaking-related damage, and tsunami-related damage, respectively, for the i-th building and the j-th stochastic realization.In Eq. (7), it is assumed that the building loss is determined by a dominant contribution from shaking and tsunami damage.This is applicable to the situation of this study where the tsunami fragility model considered [9] was developed based on actual tsunami damage data that include shaking effects implicitly.By introducing the replacement cost of a building C R and damage ratios for shaking and tsunami DR S and DR T , Eq. ( 7) can be expressed as: The damage ratios for shaking and tsunami are determined based on applicable fragility functions and damage-loss ratios.Fragility functions map values of hazard intensity measures onto damage state probabilities, which are typically defined in a discrete manner.More details of the model implementation are given in Section 2.2; in particular, an overall procedure of multi-hazard loss estimation for building portfolios (i.e.Eqs.(2) and ( 5)) is explained and illustrated in Section 2.2.7.

Models
This section presents the main components of the proposed multihazard loss estimation framework for shaking and tsunami.The models discussed are developed for the building stock in Miyagi Prefecture considering off-shore tsunamigenic sources in the Tohoku region of Japan.It is noted that most of the model components have been calibrated and compared against actual observations from the 2011 Tohoku earthquake [17,[19][20][21].Therefore, the developed multi-hazard loss model is considered to produce realistic results with respect to the 2011 Tohoku shaking and tsunami damage.2) and ( 5)).More details of the model components for the earthquake occurrence, stochastic earthquake source, earthquake shaking simulation, tsunami inundation simulation, shaking-tsunami damage assessment, and multi-hazard loss estimation are given in the following subsections.The procedure shown in Fig. 1 is versatile and therefore, the model components described below can be changed and refined, depending on the specific requirements and constraints of the loss estimation.In the developed multi-hazard loss model, epistemic uncertainties associated with the key model components are not fully characterized, which have major influence on the final hazard and risk assessments (e.g.[36,54]).This should be investigated in the future work to further refine the developed multihazard impact assessment methodology.

Earthquake occurrence model
The earthquake occurrence model is one of the most critical elements in the risk equation and is related to λ Mmin and p M in Eq. ( 2).This is because both λ Mmin and p M have direct influence on the loss ex- ceedance function ν L (L ≥ l).A popular choice for characterizing the earthquake occurrence in a given source zone is a time-independent Poisson process (λ Mmin ) for temporal occurrence combined with a truncated Gutenberg-Richter (GR) relationship [25] for earthquake magnitude (p M ).Other variations for the earthquake occurrence model include: time-dependent renewal models for temporal occurrence and characteristic magnitude models for earthquake magnitude (see [37]).A time-dependent renewal model can be implemented by specifying additional information, i.e. time horizon of the hazard-risk assessment, probability distribution for the inter-arrival time of earthquakes, elapsed time since the last event, and magnitude-recurrence model [13].
In this study, the earthquake occurrence model is based on a Poisson temporal process together with a regional GR relationship.The model set-up is the same as the seismic hazard model for the Tohoku region by the Headquarters for Earthquake Research Promotion [26], which differs from the earthquake occurrence models for other subduction zones in Japan (e.g.Nankai and Tonankai regions) where renewal-type occurrence models have been adopted.Parameters of the earthquake occurrence model are estimated based on instrumental earthquake catalogs for the target source zone off the Tohoku coast (note: the spatial process of earthquake source is characterized through stochastic earthquake rupture modeling; see Section 2.2.2).The source region approximately corresponds to the source zones of off-shore interface subduction earthquakes considered by the HERP; it is indicated as a rectangle with broken grey lines in Fig. 2a.It is important to emphasize that estimating the long-term recurrence rate for large earthquakes (> M w 8.5) involves significant uncertainty because the recording period of modern instrumental catalogs is generally short with regards to the recurrence periods of such major earthquakes.Therefore, extrapolation of the fitted magnitude-recurrence model should be considered carefully [28].Moreover, the Poisson-GR model set-up may lead to conservative estimates of the recurrence rate for large earthquakes in the Tohoku region because the constant hazard rate for the Poisson model is greater than the hazard rates indicated by renewal-type models during the early phase of a strain accumulation process.On the other hand, the strain-relaxation process due to the mega-thrust event in 2011 has not yet returned to the long-term recurrence rate of earthquakes (i.e.aftershock activity is still relatively high).A rigorous evaluation of the earthquake occurrence of tsunamigenic subduction earthquakes is outside of the scope of this study.
The regional GR relationship is obtained by analyzing seismic data from the Harvard CMT catalog (http://www.globalcmt.org/CMTsearch.html)and the NEIC catalog (http://seisan.ird.nc/USGS/mirror/neic.usgs.gov/neis/epic/code_catalog.html).Fig. 2a shows the seismicity data in the off-shore Tohoku region from the NEIC catalog, whereas magnitude-recurrence plots of the earthquake data from the two catalogs are displayed in Fig. 2b.The lower and upper cut-off magnitudes of 6.0 and 9.1 are considered.The fitted GR models indicate that the annual rate of earthquakes with M w ≥ 7.5 can be estimated to be 0.08 per year (i.e.λ Mmin ).Note that the fitted GR models shown in Fig. 2b are similar to the magnitude-recurrence model adopted by the HERP [26], which is based on the local Japan Meteorological Agency catalog.Subsequently, the discrete probability mass function for earthquake magnitude (i.e.p M ) is evaluated by considering the magnitude range between 7.5 and 9.1 (see Fig. 2c), which is divided into 8 bins with an interval of 0.2 units (note: the central value of the magnitude bin is used to represent the magnitude interval m k in Eq. ( 2)).The lower bound of M w 7.5 is suitable for off-shore tsunamigenic interface subduction earthquakes, although off-shore events with M w < 7.5 may generate damaging ground shaking at sites in the Tohoku region (however, the chance of causing major shaking damage is relatively low because the minimum source-to-site distance is about 40-50 km along the Tohoku coast and the damage threshold in terms of PGV is relatively high; see Section 2.2.6).Note also that the maximum magnitude is capped at M w 9.1 and the possibility of having greater magnitudes is neglected.It is noteworthy that the above-mentioned estimates of earthquake occurrence for large earthquakes may be biased because major historical events, such as the 869 Jogan earthquake [56], are missing in the analyzed catalogs.

Stochastic earthquake rupture model
The stochastic source modeling captures the spatial uncertainty of earthquake rupture within the target source zone.More specifically, in the developed model, the spatial source uncertainty is represented by probabilistic prediction models of earthquake source parameters and stochastic synthesis of earthquake slip [19,20].It is highlighted that the method accounts for rupture source uncertainty not only in location and geometry of the fault plane but also in earthquake slip distribution over the rupture plane.This model component is represented by f RS M | in Eq. ( 5).
First, the regional fault source model is developed by extending the fault plane geometry for the 2011 Tohoku earthquake considered by Satake et al. [51].It covers a 650 km by 250 km area and has a constant strike of 193°and variable dip angles, gradually steepening from 8°to 16°along the down-dip direction.The eastern boundary of the fault plane model approximately coincides with the Japan Trench.To characterize heterogeneous earthquake slip over the fault plane, the source zone is discretized into sub-faults having a size of 10 km by 10 km (see the top-middle panel in Fig. 1).The surface projection of the fault plane model is shown as the grey rectangle in Fig. 2a.
Next, for a given magnitude range m k , eight source parameters, i.e. fault width (W), fault length (L), mean slip (D a ), maximum slip (D m ), Box-Cox power transformation parameter for marginal slip distribution (λ), correlation length along dip (A d ), correlation length along strike (A s ), and Hurst number (H), are generated using probabilistic prediction models, which were developed based on 226 finite-fault models of the past earthquakes [20].The geometrical parameters W and L determine the size of the fault rupture, and the position of the synthesized fault plane is determined such that it fits within the source zone.The slip parameters D a and D m specify the earthquake slip statistics over the fault plane, whereas λ determines how the slip values are marginally distributed over the fault plane and is used to capture non-normal characteristics of earthquake slip (e.g. a distribution with heavier righttail than the normal distribution can be represented by λ < 1.0; see [18]).The spatial slip distribution parameters A d , A s , and H are used to characterize the heterogeneity of earthquake slip over the fault plane, represented by von Kármán wavenumber spectra [17,34].Essentially, the wavenumber spectra specify how slip values are spatially correlated over the fault plane.The multivariate lognormal prediction models by Goda et al. [20] consider correlation of prediction errors among different source parameters and thus can generate more realistic stochastic earthquake source models.
Once the spatial slip distribution parameters (i.e.A d , A s , and H) are sampled, a random slip field is generated using a Fourier integral method [45], where amplitude spectrum is represented by von Kármán spectra and its phase is uniformly distributed between 0 and 2π.To achieve slip distribution with realistic right-heavy tail features, the synthesized slip distribution is converted via Box-Cox power transformation using the simulated value of λ.The transformed slip distribution is then adjusted to achieve the target mean slip D a and to avoid very large slip values exceeding the target maximum slip D m .Due to the uncertainty in the source parameters, random sampling of W, L, and D a may result in a seismic moment M o (= μWLD a where μ is the rock rigidity) that is very different from the target moment magnitude.To avoid such an inadequate combination of W, L, and D a , sampling of the three parameters is repeated until the calculated seismic moment falls within the magnitude range (i.e.central magnitude minus/plus 0.1 units).Another important feature that is taken into account in the implemented stochastic source method is that a so-called asperity region is used to distribute more earthquake slip in the shallower part of the fault rupture zone when earthquake magnitude becomes larger.This reflects the current empirical knowledge of earthquake slip distribution for mega-thrust subduction earthquakes (e.g.[22,26,51]) and is also based on recent results from dynamic simulations of mega-thrust earthquake rupture (e.g.[43]).Further details of the stochastic synthesis can be found in Goda et al. [17,20].
In this study, 4000 stochastic source models are generated for 8 magnitude bins (i.e.500 models per magnitude range).It is noteworthy that the number of simulated source models is sufficiently large to obtain stable tsunami hazard results at the sites of interest [7].The synthesized earthquake source models, which reflect possible variability of tsunamigenic earthquakes in terms of geometry, fault location, and slip distribution, are then used to evaluate shaking and inundation hazards at locations of interest (Sections 2.2.3 and 2.2.4).To illustrate the stochastic source modeling, simulated values of L and D a of the 4000 stochastic source models are shown in Fig. 3a and b, respectively.For the fault length (Fig. 3a), it can be observed that the upper limit of 650 km (i.e.maximum length of the target source region) is reached for the M w 9.0 scenario.Due to the trade-off between L and D a in conserving the total seismic moment, simulated values of D a tend to increase for the M w 9.0 scenario (see Fig. 3b).Moreover, Fig. 3c displays three realizations of the synthesized source models for the M w 9.0 scenario, illustrating the significant variability of the geometry, location, and slip distribution of the source models.

Earthquake shaking model
Seismic intensity measures at building locations (i.e.f im rs ( | )

S
) are evaluated by using GMPEs together with suitable spatial correlation models for prediction errors.In this study, PGV is selected as IM S because most of empirical seismic fragility models for buildings in Japan [38,59,60] adopt this parameter.Among existing GMPEs, a relationship by Morikawa and Fujiwara [42] is chosen for three reasons: (i) it is applicable to mega-thrust interface subduction earthquakes in Japan, (ii) underlying data include strong motion observations from the 2011 Tohoku earthquake, and (iii) a model for PGV is available.The latter requirement excludes a regional subduction model by Zhao et al. [62] and a global subduction model by Abrahamson et al. [1].
The Morikawa-Fujiwara model computes a median PGV value for a given scenario in terms of moment magnitude (M w ) and shortest distance to rupture plane (R rup ).It also includes shallow as well as deep site correction terms, expressed as a function of average shear-wave velocity in the uppermost 30 m (V s30 ) and depth to the shear-wave velocity of 1400 m/s (D 1400 ), respectively.These additional site corrections reduce the logarithmic standard deviation of model prediction errors (i.e.sigma) significantly (from 0.340 to 0.223 in log 10 base).For each of the stochastic source models (see Fig. 3), the moment magnitude is obtained from the simulated rupture scenario, whilst the rupture distance is calculated as the minimum value of rupture distances to all sub-faults.It is noted that the distance calculation does not take into account spatial features of earthquake slip (e.g.distances to major asperities).The consequence of using such simplified distance measures is the inflated uncertainty in the developed GMPEs [18]; however, this is consistent with how empirical ground motion models are typically developed.In future, GMPEs for mega-thrust subduction earthquakes can be developed using asperity-based source representations and corresponding source-to-site distance measures.Values of V s30 are obtained from J-SHIS (http://www.j-shis.bosai.go.jp/en/), which are available at 250-m mesh grids.The value of D 1400 is set to 250 m, i.e. default value in Morikawa and Fujiwara [42].Simulation of PGV random fields is carried out by considering the 250-m grid spacing (same as the J-SHIS data).The median values at the 250-m grids are evaluated by using the Morikawa-Fujiwara equation (with M w , R rup , V s30 , and D 1400 ), whereas random error terms of the equation are simulated by considering the average spatial correlation model by Goda and Atkinson [15] and separation distance matrix of the grid points.
Simulation of PGV random fields for a given earthquake source model is illustrated in Fig. 4 by considering two example scenarios of M w 8.4 and M w 9.0.With the increase in the fault plane size, rupture distances to the target sites (i.e.grid points in southern part of Miyagi Prefecture) tend to be shorter and thus PGV values on average tend to be greater.From the simulated PGV random fields, the effects of shallow site terms (i.e.V s30 ) and spatial correlation of the prediction error terms can be observed; at locations close to the coast, V s30 values are smaller and thus PGV values tend to be amplified.To determine the PGV values at building locations (see Fig. 5), interpolation is carried out based on the simulated PGV random fields at 250-m grids.Similar shaking simulations are performed for all stochastic sources.

Tsunami inundation model
Tsunami intensity measures at building locations (i.e.f im rs ( | )

T
) are evaluated through tsunami simulations by solving nonlinear shallow water equations for initial boundary conditions of sea surface caused by earthquake rupture.In this study, inundation depth is adopted as IM T because it is the most common parameter for empirical tsunami fragility functions [58].It is noted that the momentum flux, which is a more efficient tsunami intensity measure to capture the hydrodynamic effects of tsunami waves acting on structures [47], is not considered for two reasons.Firstly, in this study, tsunami simulations are conducted at 50-m grid resolutions, which are too coarse to evaluate the flow velocity and momentum flux accurately [9].Secondly, momentum-flux-based tsunami fragility models are difficult to validate against observations.
In this study, tsunami inundation simulations are performed as follows.For each stochastic source model, initial water surface elevation is evaluated based on formulae by Okada [44] and Tanioka and Satake [57], and then tsunami wave propagation is evaluated by solving nonlinear shallow water equations with run-up [23].Note that the Okada equations do not account for the hydrodynamic response of sea water.More rigorous approaches that account for the hydrodynamic behavior of water in response to abrupt seabed deformation can be implemented using a nonhydrostatic, dispersive model (e.g.[61]), and filtering of (nonphysical) sharp peaks of seabed deformation (e.g.[12]).The computational domains are nested following a 1/3 ratio A complete dataset of bathymetry/elevation, coastal/riverside structures, and surface roughness at 50-m grid resolution is obtained from the Miyagi Prefectural Government.The bottom friction is evaluated using Manning's formula following a guideline suggested by the Japan Society of Civil Engineers [27].All bathymetry, elevation, and structural height data are defined with respect to the standard mean sea level in Japan.In the tsunami simulation, the coastal/riverside structures are represented by a vertical wall at one or two sides of the computational cells.To evaluate the volume of water that overpasses these walls, Honma's weir formulae are employed [27].The fault rupture is assumed to occur instantaneously.It is noted that for M w 9-class events, the assumption of instantaneous earthquake rupture may not be valid [51], and ideally kinematic earthquake rupture models should be considered.Although such kinematic rupture models can be implemented for stochastic earthquake sources [22], currently, probabilistic models for rupture starting points are not yet calibrated using recent large tsunami events, such as the 2004 Indian Ocean tsunami and the 2011 Tohoku tsunami.For this reason, kinematic rupture models are not implemented in this study.
Fig. 4 illustrates the simulated inundation results in southern part of Miyagi Prefecture for the two stochastic sources of M w 8.4 and M w 9.0 scenarios.It can be observed that significant increases in rupture area and earthquake slip amplitude tend to increase the extent of tsunami inundation significantly.The simulated tsunami wave heights at grid points are interpolated to obtain the wave heights at building locations; subsequently, land elevations at the building locations are subtracted from the interpolated height values to determine the inundation depths.

Building exposure data
The exposure model characterizes the assets at risk within a region of interest.In this study, a building dataset that was compiled by the Ministry of Land Infrastructure and Transportation [41] for the post-2011-Tohoku tsunami damage assessment is used.The building data span over all prefectures along the Tohoku coast, and contain formation on locations, experienced damage levels (minor damage, moderate damage, major damage, complete damage, collapse, and washed away, as defined by the MLIT), structural material (reinforced concrete, steel, wood, others, and unknown), and the number of stories (and some other details).
For the multi-hazard loss estimation, four areas in Miyagi Prefecture are focused upon: Iwanuma, Ishinomaki, Onagawa, and Shizugawa.During the 2011 Tohoku tsunami, these areas were inundated completely and most of buildings near the coast were destroyed [10].The locations of the four areas are shown in Fig. 5. Iwanuma and Ishinomaki are on low-lying coastal plain, whereas Onagawa and Shizugawa are on ria coast with complex topography.Building types that are considered in this study are low-rise wooden structures (up to 4-story buildings), for which well-calibrated shaking-tsunami fragility models are available.Zoomed maps in Fig. 5 show the spatial distribution of buildings located in Iwanuma, Ishinomaki, Onagawa, and Shizugawa.The number of wooden structures is 6096, 26,146, 1488, and 2488 for Iwanuma, Ishinomaki, Onagawa, and Shizugawa, respectively.
For the loss estimation, building cost information for the wooden buildings is needed (i.e.C R ).In this study, two sources of information, i.e. unit building cost statistics and floor area statistics, are utilized to evaluate the replacement costs of the buildings probabilistically.Using the Japanese building cost information handbook published by the Construction Research Institute [4], the mean and coefficient of variation (CoV) of the unit replacement cost for wooden buildings are obtained as 1600 USD/m 2 and 0.33, respectively (note: 1 USD is assumed as 100 yen).The adopted cost statistics are the regional average for Tohoku.Moreover, the mean and CoV of typical floor areas of wooden houses are determined as 130 m 2 and 0.33, respectively, by averaging the construction data for Miyagi and Iwate Prefectures during the period between 2012 and 2014, maintained by the MLIT (http://www.mlit.go.jp/toukeijouhou/chojou/stat-e.htm).Both unit cost and floor area are modeled by the lognormal distribution.It is noted that availability and use of detailed building cost information (unit cost and floor area) are novel aspects of the developed loss model.Based on the above building cost information, the expected total cost of all buildings in Iwanuma, Ishinomaki, Onagawa, and Shizugawa is 1268, 5438, 324, and 518 million USD, respectively.

Fragility models and damage ratios for shaking and tsunami
Fragility functions relate hazard intensity measures to probabilities of attaining damage states (i.e.f ds im im ( , )

S T
).In this study, for given values of PGV and inundation depth at building locations caused by a stochastic source model, damage ratios for shaking and tsunami, i.e.DR S and DR T in Eq. ( 7), are estimated separately for shaking and tsunami by applying seismic and tsunami fragility functions (see below).Subsequently, for each building, a greater of the estimated values of DR S and DR T , is adopted as the final damage state of the building, as indicated in Eq. ( 8).The rationale for adopting the larger value of the damage ratios is that physically shaking precedes tsunami; the tsunami fragility model by De Risi et al. [9], which is adopted in this study, was developed based on extensive tsunami damage data from the 2011 Tohoku earthquake, where the observed tsunami damage was affected by both shaking and tsunami and their effects were not distinguishable.It is noted that a possible alternative to determine the combined as well as individual effects for shaking and tsunami damage is to adopt a numerical model of a building [2,46,49] and to subject it to a suitable shaking-tsunami sequence [21].This is a future topic that is worthy of further investigations.
In this study, three empirical shaking fragility models for low-rise wooden buildings are implemented.The model by Yamaguchi and Yamazaki [60] was developed based on extensive damage data from the 1995 Kobe earthquake in Japan.Midorikawa et al. [38] developed a seismic fragility model using actual damage data from 8 major earthquakes in Japan (e.g.2004 Chuetsu earthquake) that occurred after the 1995 Kobe earthquake and prior to the 2011 Tohoku earthquake.Recently, Wu et al. [59] developed a seismic fragility model using the data from the 2011 Tohoku earthquake (note: in their analyses, near-coast locations that were inundated by the tsunami were excluded).All three models adopt PGV as seismic intensity measure.The damage states for shaking were defined as: partial damage, half collapse, and complete collapse.These damage states are closely related to the standard postearthquake damage assessment procedure in Japan, and the corresponding damage ratios can be assigned as 0.03-0.2,0.2-0.5, and 0.5-1.0,respectively [30].Fig. 6a to c show shaking fragility curves by Yamaguchi and Yamazaki [60], Midorikawa et al. [38], and Wu et al. [59] for partial-damage, half-collapse, and complete-collapse damage states, respectively (note: [38] did not develop a fragility curve for partial damage).The comparison of the three fragility models indicates that the curves for complete collapse based on the pre-2011-Tohoku and the 2011-Tohoku data differ significantly.For the base case, the three models are with an equal weight, whereas sensitivity of using different shaking fragility models is investigated in Section 3.2.
For tsunami, the empirical model by De Risi et al. [9] is adopted in this study.The model was developed based on damage observations after the 2011 Tohoku tsunami complied by the MLIT [41].The model considers five damage states for tsunami: minor, moderate, extensive, complete, and collapse (note: in [9], the washed-away damage state is integrated into the collapse damage state).According to the MLIT [41], damage ratios for the minor, moderate, extensive, complete, and collapse damage states can be assigned as: 0.03-0.1,0.1-0.3,0.3-0.5, 0.5-1.0, and 1.0 (deterministic), respectively.Fig. 6d shows tsunami fragility curves for low-rise wooden buildings for five damage states by De Risi et al. [9].Although the De Risi et al. model allows the use of bivariate tsunami intensity measures in terms of water depth and flow velocity, which improves the accuracy of predicting tsunami damage, this option is not implemented in this study because the minimum grid size of the tsunami inundation simulation (i.e.50 m) is too coarse to estimate the flow velocity at building locations accurately.Moreover, the De Risi et al. model was developed using tsunami damage data from the 2011 Tohoku earthquake with complete material and story-number information; exclusion of missing data may result in biased estimates of tsunami fragility [32].
In the MCS, the final damage state of a structure subjected to shaking and tsunami intensity values is determined as follows, noting that definitions of damage states for shaking and tsunami, as described above, are not identical.For a given PGV value, one of the three shaking fragility models is selected randomly according to assigned weights and probabilities for three damage states are calculated.By sampling a uniform random variable ranging between 0 and 1 and by comparing this simulated value with the damage state (cumulative) probabilities, the corresponding shaking damage state can be determined.Moreover, another uniform random sampling is performed to determine the shaking damage ratio DR S for the selected damage state.Similarly, for a given value of inundation depth, probabilities for five damage states are calculated, then the tsunami damage state is determined, and finally the tsunami damage ratio DR T is evaluated.By comparing the simulated values of DR S and DR T , the final damage ratio for the building and the stochastic scenario is determined based on Eq. ( 8).Subsequently, a multi-hazard loss value is calculated by sampling a value of the total replacement cost for the building (i.e.C R ) from the lognormal distributions and then by multiplying it by the final damage ratio.To obtain the portfolio multi-hazard loss for the scenario, the abovementioned MCS procedure is repeated for all buildings.

Multi-hazard loss estimation: integration of all model components
A numerical procedure of integrating the hazard and risk model components mentioned in Sections 2.2.1-2.2.6 (i.e.Eqs. ( 2) and ( 5)) is illustrated.MCS is performed by generating 500 stochastic source models for 8 magnitude bins (Section 2.2.2), and multi-hazard simulations of shaking and tsunami are carried out (Sections 2.2.3 and 2.2.4).For a building portfolio (Section 2.2.5), damage assessment and loss estimation are conducted by using relevant shaking and tsunami fragility models and building cost models (Section 2.2.6).At the end of the MCS, loss samples for all buildings are available for the 4000 stochastic source models.As indicated in Eq. ( 5), these loss samples can be used to construct the conditional probability distribution function of the total portfolio loss (i.e.P(L ≥ l|m k )) for a given magnitude bin.Subsequently, the conditional loss curves for different magnitude bins are integrated by considering λ Mmin and p M (Section 2.2.1) to obtain the multi-hazard loss exceedance curve for the building portfolio (i.e.ν L (L ≥ l); see Eq. ( 2)).Fig. 7 shows the conditional loss distribution curves for buildings in Iwanuma for two magnitude bins, i.e.M w 8.4 and M w 9.0.In the figure, three curves are included: combined shaking-tsunami loss, and shaking loss, and tsunami loss.The calculation procedures for the single-hazard (shaking or tsunami) loss is almost identical to that for the combined multi-hazard loss except that only one type of the hazards is considered in the loss estimation.Fig. 7a shows that for the M w 8.4 scenario, shaking-related loss dominates tsunami-related loss, whereas Fig. 7b shows that for the M w 9.0 scenario, the trend is reversed and the large loss is predominantly caused by tsunami damage.Differences between the multi-hazard loss curve and the single-hazard loss curve (either shaking or tsunami) provide useful information regarding the impact of ignoring multi-hazard loss.It is noteworthy that simple sum of singlehazard loss values at a given probability level does not always equate the combined multi-hazard loss value at the same level because these loss cases usually correspond to different stochastic source models.
Fig. 8 illustrates the calculation step of integrating different conditional loss distributions into the final loss exceedance curve.To evaluate the values of ν L (L ≥ l) as a function of loss threshold l, P(L ≥ l|m k ) can be obtained from the conditional probability distribution functions an example of this operation for l = 500 million USD is shown in Fig. 8. Once the conditional exceedance probabilities of the loss are evaluated for all magnitudes (Fig. 8a), they are weighted by their occurrence probabilities (i.e.λ Mmin and p M ) and are summed to obtain the unconditional loss estimate (Fig. 8b), noting that values of P(L ≥ l|m k ) for M w < 8.6 are negligible for l = 500 million USD and thus are omitted.By repeating the above procedures for different threshold values, the loss exceedance curve, which reflects a range of magnitude scenarios from M w 7.5 to M w 9.1, can be obtained.

Application
A loss exceedance curve, obtained from the developed multi-hazard loss model, quantifies potential loss in terms of frequency and severity, and provides valuable insight on the financial impact due to cascading shaking-tsunami hazards.This section presents an application of the developed multi-hazard loss estimation framework to wooden buildings in Miyagi Prefecture.In Section 3.1, loss estimation results for the four areas shown in Fig. 5 are discussed by focusing upon key features of the multi-hazard loss exceedance curves in comparison with conventional single-hazard loss exceedance curves.In particular, loss contributions from shaking and tsunami damage are focused upon.In Section 3.2, sensitivity of loss exceedance curves to shaking fragility models is investigated, whereas in Section 3.3, critical shaking-tsunami concurrent  hazard maps are derived at several return period levels T R .

Multi-hazard versus single-hazard loss exceedance curves
Fig. 9 compares loss exceedance curves for multi-hazard and singlehazard cases for Iwanuma, Ishinomaki, Onagawa, and Shizugawa.In each plot, three curves, i.e. combined shaking-tsunami loss, shaking loss, and tsunami loss, are included, and the total asset values of buildings within the areas are indicated by vertical broken lines.The difference between the vertical broken line and the lower-end-point of the loss curve at annual exceedance probability of 10 −4 can be used as an approximate indicator for the severity of a hazard of interest for a worst situation with respect to the total asset value at risk.For instance, in Fig. 9a, the maximum loss values attained by the combined loss, shaking loss, and tsunami loss are 1194, 165, and 1184 million USD, respectively, indicating that even for the extreme situations, the shaking-related loss for the building portfolio in Iwanuma may be at most 13% of the total asset value, whilst the combined shaking-tsunami and tsunami-related losses could destroy almost all buildings in Iwanuma (about 93-94%).The approximate upper limit of the shaking damage for the building portfolio is valid because the magnitude scaling of seismic intensity parameters for large earthquakes (> M w 8.2) reaches a plateau (i.e.saturation in magnitude scaling; [42]) and realistic spatial correlation of prediction errors is taken into account in the shaking simulations [15].
A common observation from all four plots shown in Fig. 9 is that shaking loss curves intersect with tsunami loss curves.At shorter return periods (e.g.T R < 200 years), shaking loss curves are generally positioned on the right-hand side with respect to tsunami loss curves, whereas at longer return periods (e.g.T R > 1000 years), their relative positions are reversed and loss contributions due to tsunami damage become dominant (i.e. the combined loss curves and the tsunami loss curves nearly coincide).At intermediate return periods, crossing of the two single-hazard curves occurs.In other words, for buildings in the four areas, shaking damage tends to occur more frequently but overall impact is somewhat limited at the local level.On the other hand, tsunami damage is relatively rare but could be devastating.It is important to emphasize that relative loss contributions the shaking and tsunami damage to total loss depend on various factors.For the shaking loss curves, the proximity of the area to the fault rupture plane and the surface soil types at the building sites are crucial because these two factors control the severity of strong shaking in the area.For the tsunami as well as combined loss curves, proportions of buildings that are located within low-elevation areas are the decisive factor because buildings at relatively high elevations or far from the coast will not be inundated by the tsunami.In this regard, buildings in Iwanuma, in extreme situations, could be washed away completely (as they are on a flat plain), while not all buildings in Onagawa will be inundated (as they are on a slope near valleys; noting that they could be still damaged by strong shaking).
Relative loss contributions of different hazards to total loss can be better understood by examining the dominant sources of loss generation as a function of earthquake magnitude.Fig. 10 shows proportions of three loss sources, i.e. no loss, shaking loss, and tsunami loss, as a function of earthquake magnitude.It is noted that the proportions shown in Fig. 10 are defined based on the number of incidences of loss generation, and they are not based on the size of the losses caused by shaking or tsunami damage.For each magnitude bin, the numbers of no loss cases, shaking-dominated cases, and tsunami-dominated cases for all stochastic source models (within the magnitude bin) and all buildings (within the selected area) are counted, and the normalized numbers are plotted (i.e. the sum of the three proportions equals 1.0).General trends of the loss contribution plots indicate that with the increase in moment magnitude, the loss proportion curves for shaking and tsunami increase, while the curve for no loss decreases (as expected).The trends for the shaking loss are concave and become flat, while those for the tsunami loss increase rapidly.These features essentially corroborate the trends observed for the loss exceedance curves for shaking and tsunami (Fig. 9).Among the four areas that are considered in this study, for the M w 9.0 scenario, tsunami-related loss contributions become particularly dominant for Iwanuma and Shizugawa, where the majority of wooden buildings are located within the low-lying flat areas near the coast (i.e.prone to complete inundation in extreme situations, as observed during the 2011 Tohoku tsunami; [10]).On the other hand, the numbers of loss incidences due to shaking and tsunami for Ishinomaki and Onagawa are similar for the M w 9.0 scenario (i.e. more or less equal chance of loss generation due to shaking and tsunami).The latter is because some proportions of the buildings in Ishinomaki and Onagawa are located in unflooded areas.

Sensitivity of loss exceedance curves to shaking fragility models
In Section 3.1, an equal weighting is adopted to capture epistemic uncertainty of seismic fragility models in the probabilistic multi-hazard loss estimation.However, as mentioned and demonstrated in Section 2.2.6 (Fig. 6a-c), variations of applicable seismic fragility models for   11, there is one loss exceedance curve for tsunami-related damage only, while there are three loss curves for the combined damage and shaking-related damage.For the shaking loss curves (in red), the effects of using different combinations of seismic fragility models can be significant; as expected, when the Wu et al. model alone is considered, the shaking loss exceedance curve is shifted to left with respect to the base case (i.e.equal weighting of the three models).For the case of the combined shakingtsunami loss (in blue), the effects due to the shaking fragility models are noticeable, especially for relatively short return periods (e.g.T R < 500 years) because at these probability levels, shaking-related damage contributes significantly to the total loss.Therefore, careful selection of suitable shaking fragility models is an important consideration in multihazard loss estimation for coastal communities (note: its impact will be greater when buildings far from the coastal line are of interest, which essentially becomes a single shaking hazard case).

Critical shaking-tsunami concurrent hazard maps
The developed multi-hazard loss estimation framework can offer new insight regarding the concurrent mapping of shaking and tsunami hazards.More specifically, shaking-tsunami concurrent hazard maps (see Fig. 4), caused by the same earthquake rupture scenario, can be produced by referring to a return period level of the multi-hazard loss exceedance curve (see Fig. 9).It is highlighted that concurrent hazard maps produced in this way differ from hazard maps that are developed separately from probabilistic seismic hazard-risk analysis and probabilistic tsunami hazard-risk analysis.From earthquake disaster risk management perspectives, concurrent shaking-tsunami hazard maps are more useful than disjoint single-hazard maps, because for the latter, the hazard/risk scenarios are not defined for the realistic multi-hazard situations and the two hazard maps do not correspond to the same earthquake rupture scenario.
To demonstrate concurrent shaking-tsunami hazard mapping at several return periods that are of practical interest in disaster risk management, three sets of an earthquake source model, PGV distribution, and inundation height distribution at return periods of 100, 500, and 1000 years are developed for Iwanuma and Onagawa.The results are shown in Figs. 12 and 13, respectively.The buildings within the selected areas are shown with grey dots, plotted over the PGV and inundation distributions.In the earthquake rupture source panel (left), the total multi-hazard loss value at the selected return period level is indicated (note: the corresponding scenario magnitude of the source model is also included in the panel), whilst in the PGV and inundation height distribution panels (middle and right), the single-hazard loss values at the return period level are indicated.It is important to note that for a given return period level, multiple sets of such concurrent shaking-tsunami hazard maps can be found; they differ in various aspects, i.e. different rupture characteristics and shaking-tsunami hazard distributions, but result in similar extent of the building portfolio loss.The results shown in Figs. 12 and 13 are the cases that are closest to the selected return period levels in terms of multi-hazard loss.
Figs. 12 and 13 clearly show that with the increase in return period levels, earthquake scenarios become more critical (i.e.larger earthquake magnitude with larger fault plane size and higher slip).Consequently, values shown in the concurrent hazard maps become severer.For the PGV distribution, hazard values become greater (i.e. more redish colors) and all buildings are subjected to non-negligible shaking.On the other hand, the inundation heights along the coast become greater and the inundation distribution is expanded to cover wider areas, thus causing more and more buildings to be flooded by relatively high tsunami waves.By comparing the results for Iwanuma and Onagawa, the effects of different topographical features (i.e.coastal plain versus ria coast), which have major influence on tsunami loss estimation, can be observed.Visualizing the multi-hazard earthquake impact for shaking and tsunami is particularly effective in discussing and communicating the hazard and risk prediction assessment results with local stakeholders.

Conclusions
This study developed a novel probabilistic multi-hazard loss estimation methodology that is applicable to multiple structures (i.e.building portfolio) located in coastal areas subjected to cascading shaking-tsunami hazards.The method is innovative in two aspects: (i) common source effects for shaking damage and tsunami damage are captured, and (ii) uncertainties associated with earthquake source modeling are fully taken into account by integrating new prediction models of earthquake source parameters and stochastic synthesis of heterogeneous earthquake slip.Subsequently, uncertainties in earthquake occurrence, rupture characteristics, multi-hazard simulations, and shaking-tsunami damage assessment, were accounted for and were propagated to estimate the total loss to a building portfolio.The method was demonstrated by applying to a realistic case study for wooden buildings located in Miyagi Prefecture, Japan, where the catastrophic tsunami struck in 2011.As part of the example, loss contributions from shaking and tsunami damage were discussed in detail; sensitivity of loss exceedance curves to shaking fragility models was investigated; and critical shaking-tsunami concurrent hazard maps were derived at several return period levels.
The developed method is particularly useful for evaluating the earthquake impact at regional and local levels.It also promotes a performance-based earthquake-tsunami engineering methodology, which is compatible with current single-hazard methods, such as performancebased earthquake engineering and performance-based tsunami engineering.The computational framework can be extended to include other major secondary geo-hazards, such as landslide and liquefaction.
Many applications can be envisaged in future.The proposed framework can be used for quantitative cost-benefit analysis of multi-hazard risk mitigation measures and will promote risk-informed management as well as financial decisions related to shaking-tsunami disaster risk reduction.In the catastrophe modeling field, fully probabilistic multi-hazard loss estimation tools will meet demands from clients in promoting more accurate assessments of multi-hazard insurance and reinsurance coverages and in designing multi-hazard financial risk transfer products.

Fig. 1
Fig.1illustrates a MCS-based computational procedure for evaluating the risk equation (i.e.Eqs.(2) and (5)).More details of the model components for the earthquake occurrence, stochastic earthquake source, earthquake shaking simulation, tsunami inundation simulation, shaking-tsunami damage assessment, and multi-hazard loss estimation are given in the following subsections.The procedure shown in Fig.1is versatile and therefore, the model components described below can be changed and refined, depending on the specific requirements and constraints of the loss estimation.In the developed multi-hazard loss model, epistemic uncertainties associated with the key model components are not fully characterized, which have major influence on the final hazard and risk assessments (e.g.[36,54]).This should be investigated in the future work to further refine the developed multihazard impact assessment methodology.

Fig. 2 .
Fig. 2. (a) Regional seismicity in the Tohoku region based on the NEIC catalog, (b) Gutenberg-Richter models for the offshore Tohoku region based on the Harvard-CMT (HCMT) and the NEIC catalogs, and (c) conditional distribution of earthquake magnitude for the off-shore Tohoku region.

Fig. 3 .
Fig. 3. (a,b) Scaling relationships for fault length and mean slip by Goda et al. [21], in comparison with the simulated fault length and mean slip of the 4000 stochastic source models, and (c) three realizations of the stochastic source models for the M w 9.0 scenario.

Fig. 4 .
Fig. 4. Stochastic earthquake source models and results for earthquake shaking and tsunami inundation simulations: (a) M w 8.4 scenario and (b) M w 9.0 scenario.

Fig. 5 .
Fig. 5. Elevation data of Miyagi Prefecture and locations of four coastal areas for the multi-hazard loss estimation.Zoomed maps display the spatial distributions of low-rise wooden buildings in the four areas (the numbers in the brackets are those for buildings for different stories in the area).

Fig. 7 .
Fig. 7. Conditional probability functions of multi-hazard and single-hazard losses for buildings in Iwanuma: (a) M w 8.4 and (b) M w 9.0.

Fig. 8 .Fig. 9 .
Fig. 8. (a) Conditional probability distribution functions the multi-hazard loss for buildings in Iwanuma for 8 magnitude bins and (b) multi-hazard loss exceedance curve for buildings in Iwanuma.