Tsunami hazard and risk assessment for multiple buildings by considering the spatial correlation of wave height using copulas

It is necessary to evaluate aggregate damage probability to multiple buildings when performing probabilistic risk assessment for the buildings. The purpose of this study is to demonstrate a method of tsunami hazard and risk assessment for two buildings far away from each other, using copulas of tsunami hazards that consider the nonlinear spatial correlation of tsunami wave heights. First, we simulated the wave heights considering uncertainty by varying the slip amount and fault depths. The frequency distributions of the wave heights were evaluated via the response surface method. Based on the distributions and numerically simulated wave heights, we estimated the optimal copula via maximum likelihood estimation. Subsequently, we evaluated the joint distributions of the wave heights and the aggregate damage probabilities via the marginal distributions and the estimated copulas. As a result, the aggregate damage probability of the 99th percentile value was approximately 1.0 % higher and the maximum value was approximately 3.0 % higher while considering the wave height correlation. We clearly showed the usefulness of copula modeling considering the wave height correlation in evaluating the probabilistic risk of multiple buildings. We only demonstrated the risk evaluation method for two buildings, but the effect of the wave height correlation on the results is expected to increase if more points are targeted.

1 Introduction

Probabilistic hazard and risk assessment methods of disasters are developed mainly in the field of nuclear safety focused on countermeasures relative to severe accidents at nuclear power plants. Among them, a variety of probabilistic tsunami hazard assessment (PTHA) and probabilistic tsunami risk assessment (PTRA) methods for tsunami disasters have been rapidly developed since the 2000s (e.g., Geist and Parsons, 2006; Annaka et al., 2007; González et al., 2009; Thio et al., 2010; Løvholt et al., 2012, 2015; Goda et al., 2014; Fukutani et al., 2015; Park and Cox, 2016; De Risi and Goda, 2017; Grezio et al., 2017; Davies et al., 2018). The main purpose of a PTHA is to assess the likelihood of a given measure of tsunami hazard metrics (e.g., maximum tsunami wave height) being exceeded at a particular location within a given time period. The most basic outcome of such an analysis is typically expressed as a hazard curve, which shows the exceedance level of the hazard metric with the probability. This is often expressed as a rate of exceedance per year. A PTHA can be expanded to a PTRA by combining hazard assessment with loss evaluation of a target. Several studies have proposed a method of PTRA for an individual site in a local area. Detailed risk assessment is undoubtedly important in terms of grasping the risk of exposing assets located in a local area.

However, probabilistic risk evaluation methods are also utilized in cases to evaluate risks for multiple buildings (e.g., Kleindorfer and Kunreuther, 1999; Chang et al., 2000; Grossi and Kunreuther, 2005; Goda and Hong, 2008; Salgado-Gálvez et al., 2014; Scheingraber and Käser, 2019). With respect to businesses that own a building portfolio, including factories and offices over a wide area, it is extremely important in risk-based management decisions to evaluate the detailed risks posed by the building portfolio. A portfolio means a collection of assets held by an institution or a private individual. By quantitatively assessing the risks posed by the building portfolio, for example, it is possible to identify assets held that have a large impact on the overall risk and to compare the amount of risk held over time, which leads to support for decision-makers.

When evaluating physical risks for multiple buildings over a wide area, it is necessary to evaluate the aggregate risk for the buildings that are located at a distance. In these types of cases, it is necessary to evaluate the risk by considering the spatial correlation of hazards. For example, let us consider assessing the risk of two buildings located at two sites. When the positive correlation of hazards between two sites is strong, the hazard at one site tends to be large if the hazard at another site is large. In this case, the hazards at the two target sites both increase, and as a result, the aggregate risk for the two buildings considering the hazard correlation increases. Conversely, when the positive correlation of hazards is small, the hazard at one site is not necessarily large, even if the hazard at another site is large. In this case, compared to the former case, the hazards at the two target sites are smaller, and as a result, the aggregate risk for the two buildings is smaller if we assume that the vulnerability of the two buildings is equal. Therefore, analyses that do not consider the spatial correlation of hazards involve the risk of underestimating the risk over a wide area. It is clear that the difference of aggregate risk between two cases becomes more prominent as the number of target sites increases. Analyses that consider the spatial correlation of hazards are relatively advanced in the field of earthquake hazard and risk assessment (e.g., Boore et al., 2003; Wang and Takada, 2005; Park et al., 2007) albeit insufficient in the field of tsunami hazard and risk assessment. Analyses that consider the hazard correlation using copulas are used in hydrological/earthquake modeling (e.g., Goda and Ren, 2010; Goda and Tesfamariam, 2015; Salvadori et al., 2016) although there is a paucity of the same in tsunami modeling.

In this study, we assume the occurrence of a large earthquake in the Sagami Trough in Japan that significantly affects the metropolitan area and evaluate the tsunami risk of two buildings located at distant locations by considering the spatial correlation of the tsunami wave height between the two sites. The objective of this study involves evaluating the frequency distribution of the tsunami height via the response surface method and evaluating the spatial correlation of the tsunami heights and damages by using various copulas. Specifically, we analyze the frequency distribution (marginal distribution) of tsunami height via the response surface method and target two steel buildings located at Oiso and Miura along the Sagami Bay, Kanagawa Prefecture, in Japan. Subsequently, we derive the joint distribution of tsunami wave heights between two sites by using various copulas and the marginal distributions, convert it to the joint distribution of damage by applying a damage function, and evaluate the expected value of the aggregate damage probability for the target buildings. Finally, we confirm the extent to which the expected value of the aggregate damage probability fluctuates in a case where the spatial correlation of tsunami wave height is considered and a case where it is not considered.

Section 2 provides an outline of the response surface method and tsunami hazard and risk assessment method for multiple buildings using copulas. Section 3 describes a case where the proposed method is applied to the Sagami Trough area. The final conclusions are discussed in Sect. 4.

2 Methodology

Figure 1 shows a flowchart of tsunami hazard and risk assessment considering the correlation of tsunami wave heights in this study. Herein, the risk assessment target points only correspond to two points: Oiso and Miura, Kanagawa Prefecture, in Japan. Figure 2 shows the location of these points. First, we simulate the tsunami wave heights considering the uncertainty at the target sites by numerical tsunami simulations via nonlinear long-wave equations. Based on this, we construct a response surface and apply probability distributions to obtain a frequency distribution of tsunami wave heights. This distribution becomes a marginal distribution for a joint distribution of tsunami wave heights of two target points. Separately, we estimate appropriate copula via maximum likelihood estimation from the simulation results of the tsunami wave height considering uncertainty. Subsequently, we obtain a joint distribution of tsunami wave heights from the estimated copula and the marginal distributions of tsunami wave height. Furthermore, we obtain a joint distribution of damage probabilities by applying the tsunami damage function.

Figure 1 Flowchart of probabilistic tsunami hazard and risk assessment considering the spatial correlation of tsunami wave height. Numbers in the parentheses indicate the section numbers escribed.

The outline of the response surface method and copula modeling used in this study is explained below. The response surface method is a statistical combination method to determine an optimum solution using the lowest number of measurement data possible. The basic idea is based on a reliability-based design scheme developed in the research field of geomechanics (e.g., Honjo, 2011). Generally, the response surface model is given by Eq. (1) as follows:

$$\begin{array}{}\text{(1)}& y=f\left(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n\right)+\mathit{\epsilon },\end{array}$$

where explanatory variables correspond to x_i (i=1, 2, 3, …, n), response (object variable) corresponds to y, and error corresponds to ε. It should be noted here that a response surface is generated for a certain point. Therefore, it is necessary to generate a large number of response surfaces with spatial meshes in order to evaluate the spatial inundation height and flow depth variability, but such an analysis is outside the scope of this study. Tsunami hazard assessment has many uncertainties in each process of tsunami generation, propagation, and run-up. Even considering only the earthquake source parameters that are the basis for calculating the initial displaced water level of the tsunami, there are fault length, fault width, fault depth, slip amount, rake, strike, and dip. The temporal and spatial changes of all these parameters more or less affect the tsunami hazard assessment. Numerous studies on the effect of earthquake source parameters on the initial displaced water level of tsunamis have been conducted (e.g., Hwang and Divoky, 1970; Ward, 1982; Ng et al., 1991; Pelayo and Wiens, 1992; Whitmore, 1993; Geist and Yoshioka, 1996; Geist, 1999, 2002; Song et al., 2005). These studies reported that fault slip was an important factor governing tsunami intensity. In addition, the Sagami Trough, which is the target earthquake of this study, has a complex crustal structure in the area where the Pacific Plate, the Philippine Sea Plate, and the North American Plate meet. Therefore, the depth where the Sagami Trough earthquake occurs is considered uncertain. Therefore, in this study, we decided to consider only the tsunami hazard uncertainty caused by the changes of slip amount and fault depth as an example. The heterogeneity of fault slip is an equally important factor, but we did not consider nonuniform slip distribution for purposes of simplicity. It is an important issue in the future to evaluate the heterogeneity of fault slip using response surface methodology. This is true for both slip heterogeneity and other fault parameters. For the above reasons, we model maximum tsunami wave height considering tsunami wave uncertainty with Eq. (2) after conducting a tsunami numerical simulation with a nonlinear long-wave equation. This formula is following the tsunami hazard evaluation method proposed by Kotani et al. (2016) that applied a reliability analysis framework using the response surface method proposed in Honjo (2011). The expression is as follows:

$$\begin{array}{}\text{(2)}& h(S,D)=aS+bD+cSD+d{S}^{\mathrm{2}}+e,\end{array}$$

where h(S, D) denotes the tsunami wave height; S denotes the slip; D denotes the fault depth; and a, b, c, d, and e denote the undetermined coefficients. It should be noted that an error term is not included in Eq. (2). An example of the error term is to consider an error due to modeling. For example, Kotani et al. (2016) quantified the modeling error as the difference between the observed tsunami height and the numerically simulated tsunami height. The modeling error of the numerical analysis was also considered as one of the tsunami hazard uncertainties. However, the main purpose of this study is to propose a tsunami damage assessment method for multiple buildings using a copula considering wave height correlation. Therefore, the modeling error is also ignored for simplification in this study.

This response surface method has an advantage that the probability distribution of the objective variable can be easily evaluated by applying an appropriate probability distribution to the explanatory variable and performing a Monte Carlo simulation. Although the tsunami numerical simulation considering uncertainty usually has a high calculation cost to conduct vast numbers of simulation cases, it is possible to significantly reduce the simulation cost by using the response surface method.

Figure 2 (a) Major subduction-zone earthquakes around the Japanese islands including the Sagami Trough earthquake, the Nankai Trough earthquake, and the Tohoku-type earthquake (yellow area); (b) two targets points, Oiso and Miura, Kanagawa Prefecture, for tsunami hazard and risk assessment.

Figure 3 A simple synthetic example of a copula in a bivariate case. (a) Joint distribution; (b, c) are distribution functions of each variable (marginal distribution) and (d) is a copula distributed over [0, 1].

The foundation of the copula theory corresponds to the Sklar theorem (Sklar, 1959). A copula is a multivariate distribution whose marginals are all uniform over [0, 1]. Given this in combination with the fact that any continuous random variable can be transformed to be uniform over [0, 1] by its probability integral transformation, copulas are used to separately provide multivariate dependence structure from the marginal distributions. Let F be a n-dimensional distribution function with marginals F₁, …, F_n and H be a joint distribution function. There exists a n-dimensional copula C such that for all x in the domain of F, the following expression holds (Sklar, 1959):

$$\begin{array}{}\text{(3)}& H\left(x_{\mathrm{1}},\mathrm{\dots },x_n\right)=C\left\{F_{\mathrm{1}}\left(x_{\mathrm{1}}\right),\mathrm{\dots },F_n\left(x_n\right)\right\}=C\left(u_{\mathrm{1}},\mathrm{\dots },u_n\right),\end{array}$$

where $u_i=F_i\left(x_i\right)\in [\mathrm{0}$, 1], i=1, …, n. Figure 3 shows a simple synthetic example of a copula in a bivariate case. Figure 3a is a joint distribution function, Fig. 3b and c are distribution functions of each variable (marginal distributions), and Fig. 3c is a copula distributed over [0, 1]. Joe (1997) and Nelsen (1999) proposed the two comprehensive treatments on the topic. The two most common elliptical copulas correspond to the Gaussian copula and the t copula whose copula functions in the bivariate case correspond to Eqs. (4) and (5).

$$\begin{array}{}\text{(4)}& C\left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)={\mathrm{\Phi }}_{\mathrm{\Sigma }}\left({\mathrm{\Phi }}^{-\mathrm{1}}\left(u_{\mathrm{1}}\right),{\mathrm{\Phi }}^{-\mathrm{1}}\left(u_{\mathrm{2}}\right)\right)\end{array}$$

$$\begin{array}{}\text{(5)}& C\left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)=t_{\mathrm{\Sigma },\mathit{\nu }}\left(t_{\mathit{\nu }}^{-\mathrm{1}}\left(u_{\mathrm{1}}\right),t_{\mathit{\nu }}^{-\mathrm{1}}\left(u_{\mathrm{2}}\right)\right)\end{array}$$

The Gaussian copula is simply derived from a multivariate Gaussian distribution function Φ_Σ with mean zero and correlation matrix Σ by transforming the marginals by the inverse of the standard normal distribution function Φ. Given a multivariate centered t-distribution function t_Σ,ν with correlation matrix Σ, ν degrees of freedom, and with marginal distribution function t_ν, the t copula is derived in the same way as the Gaussian copula. The Archimedean copula is a widely used copula family. The Archimedean copulas include the Gumbel, Frank, and Clayton copulas whose copula functions in the bivariate case correspond to Eqs. (6)–(8), respectively, as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{C}_{\mathit{\theta }}\left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)=\exp \left\{-{\left[{\left(-l{\mathit{\nu }}_{\mathrm{1}}\right)}^{\mathit{\theta }}+{\left(-l{\mathit{\nu }}_{\mathrm{2}}\right)}^{\mathit{\theta }}+\right]}^{\mathrm{1}/\mathit{\theta }}\right\},\mathit{\theta }\ge \mathrm{1},\text{(7)}& {\displaystyle }\begin{array}{rl}{C}_{\mathit{\theta }}\left(u_{\mathrm{1}},u_{\mathrm{2}}\right)=& -{\displaystyle \frac{\mathrm{1}}{\mathit{\theta }}}\ln \left\{\mathrm{1}+{\displaystyle \frac{\left(\exp \left(-\mathit{\theta }{u}_{\mathrm{1}}\right)-\mathrm{1}\right)\left(\exp \left(-\mathit{\theta }{u}_{\mathrm{2}}\right)-\mathrm{1}\right)}{\exp (-\mathit{\theta })-\mathrm{1}}}\right\},\\ & -\mathrm{\infty }<\mathit{\theta }<\mathrm{\infty },\end{array}\text{(8)}& {\displaystyle }{C}_{\mathit{\theta }}\left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)={\left(u_{\mathrm{1}}^{-\mathit{\theta }}+u_{\mathrm{2}}^{-\mathit{\theta }}-\mathrm{1}\right)}^{-\mathrm{1}/\mathit{\theta }},\mathit{\theta }\ge \mathrm{1}.\end{array}$$

The Gumbel and Clayton copulas capture upper tail dependence and lower tail dependence, respectively, while the Frank copula does not exhibit tail dependence. Specifically, θ is estimated based on the maximum log-likelihood method. The copulas denote the symmetrical property with respect to diagonal lines of a unit square. To handle asymmetrical data in transformed space, we used an asymmetrical extreme-value copula (Tawn, 1988; Genest and Favre, 2007; Genest and Segers, 2009). Extreme-value copulas are characterized by the dependence function A as given in Eq. (9):

$$\begin{array}{}\text{(9)}& C\left(u_{\mathrm{1}},u_{\mathrm{2}}\right)=\exp \left[\log \left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)A\left\{{\displaystyle \frac{\log \left(u_{\mathrm{1}}\right)}{\log \left(u_{\mathrm{1}}{u}_{\mathrm{2}}\right)}}\right\}\right].\end{array}$$

An asymmetric model using the copula with three parameters as mentioned by Tawn (1988) is given by

$$\begin{array}{}\text{(10)}& A\left(t\right)={\left\{{\mathit{\theta }}^r(\mathrm{1}-t{)}^r+{\mathit{\phi }}^r{t}^r\right\}}^{\mathrm{1}/r}+(\mathit{\theta }-\mathit{\phi })t+\mathrm{1}-\mathit{\theta },\end{array}$$

where r, θ, and φ are estimated based on the maximum log-likelihood method. The special case θ=1 and φ=1 corresponds to the symmetric model proposed by Gumbel (1960), and thus this is termed as the asymmetric Gumbel copula. We use this copula for modeling asymmetrical data dependence.

In this study, we use the bivariate case as the tsunami wave height at two target points and model the correlation using a copula. The linear correlation coefficient (Pearson's correlation coefficient) is an index that captures the linear relation between variables and essentially cannot express the dependency between variables that are not in linear relation. Conversely, the copula is a function that expresses the correlation based on the order of the data of each variable rather than the data themselves. The order of the data is expressed by Kendall's τ (Kendall, 1938). Therefore, it is possible to quantify the nonlinear correlation between the variables. Table 1 shows theoretical value of Kendall's τ corresponding to the bivariate copulas and their parameter vectors. In this study, we show a simple evaluation method for two target points, although correlation between more points can be considered by using copulas.

Table 1 Bivariate copula, parameter vectors, and Kendall's τ.

ρ: Pearson's correlation coefficient; $D_{\mathrm{1}}\left(\mathit{\theta }\right)=\underset{\mathrm{0}}{\overset{\mathit{\theta }}{\int }}\frac{\frac{x}{\mathit{\theta }}}{\left(e^x-\mathrm{1}\right)}\mathrm{d}x$: the first Debye function.

3 Application to the Sagami Trough area

In this chapter, we demonstrate a case study where the hazard and risk assessment method described in the previous chapter is applied for two buildings located on the coast of Sagami Bay, Kanagawa Prefecture, in Japan. Section 3.1 shows the assessment target points, Sect. 3.2 shows the tsunami numerical simulation considering uncertainties, Sect. 3.3 constructs the response surface, Sect. 3.4 shows the modeling of tsunami wave height correlation using copulas, and Sect. 3.5 shows the results of the evaluation and discussion.

3.1 Risk assessment targets

Figure 2a shows major subduction-zone earthquakes around the Japanese islands, namely the Sagami Trough earthquake, the Nankai Trough earthquake, and the Tohoku-type earthquake announced by NIED (2017). Figure 2b shows the located points of tsunami hazard and risk assessment targets, namely Oiso and Miura, Kanagawa Prefecture, in Japan. The Sagami Trough earthquake covers most of the Kanto region, including the target points. Oiso is located at the approximate center of Sagami Bay coast, and Miura is located at the tip of the Miura Peninsula, which is located between Tokyo Bay and Sagami Bay. We assume a steel-framed building located at these two points and evaluate the tsunami damage probability for the two buildings.

3.2 Tsunami numerical simulation considering uncertainties

In this section, we evaluate the tsunami wave heights by considering the uncertainty at the target points.

Figure 4 (a) The 10 sources of the Sagami Trough earthquakes (NIED, 2017) and (b) initial water levels of the tsunami calculated from the fault parameters using the Okada equation (Okada, 1985). © OpenStreetMap contributors 2019. Distributed under a Creative Commons BY-SA License.

We selected 10 earthquake occurrence sources of the moment magnitude (M_w) 8 class along the Sagami Trough, which significantly affect the metropolitan area in Japan. The Sagami Trough is a 300 km long boundary between the Philippine Sea and North American plates. The assumed earthquake sources are shown in Fig. 4a. There are 10 earthquake sources, and the M_w of the sources ranges from M_w=7.9 to M_w=8.6. Source 8 has maximum M_w=8.6. The sources are used for probabilistic ground motion prediction in Japan published by NIED (2017), and thus they exhibit a 0.7 % occurrence probability in the next 30 years, and the weights of occurrence probability are used for each earthquake source. Table 2 shows the number of small faults in each source. Each small fault corresponded to a 2.5 km square, and the slip amount of the fault was set to a uniform value based on the moment magnitude (M_w) of each earthquake by using the following scaling laws of earthquakes according to Kanamori (1977):

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{Mo}=\mathit{\mu }SA,\text{(12)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{w}}={\displaystyle \frac{{\log }_{\mathrm{10}}\mathrm{Mo}-\mathrm{9.1}}{\mathrm{1.5}}},\end{array}$$

where “Mo” denotes moment magnitude (Nm), μ denotes shear modulus (Pa), S denotes slip amount (m), and A denotes earthquake source area (m²). μ was set to 3.4×10¹⁰ (Pa). In this study, we did not consider nonuniform slip distribution for purposes of simplicity. We set other fault parameters (i.e., fault depth, dip, rake, and strike) to the sources based on information published by the Cabinet Office (2013) in Japan, which were created from the crustal structure of data of the plates.

Figure 5 Tsunami numerical simulation results (a Oiso and b Miura) in the case changing the M_w (moment magnitude) and the fault depth of source 8.

Table 2 Moment magnitude, average slip, number of faults, and area in each earthquake source of the Sagami Trough earthquake.

Figure 4b shows the calculation results of the initial water level distribution of the tsunami using the Okada (1985) equation. The initial water level of up to approximately +3.5 m is distributed off to Sagami Bay and Tokyo Bay. Using the initial water level as an input value, we performed a tsunami numerical simulation via a nonlinear long-wave equation. We use the following continuity equation (Eq. 13) and nonlinear shallow water equations (Eqs. 14 and 15) as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial M}{\partial x}}+{\displaystyle \frac{\partial N}{\partial y}}=\mathrm{0},\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial M}{\partial t}}& +{\displaystyle \frac{\partial }{\partial x}}\left[{\displaystyle \frac{M^{\mathrm{2}}}{D}}\right]+{\displaystyle \frac{\partial }{\partial y}}\left[{\displaystyle \frac{MN}{D}}\right]+gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}\\ & +{\displaystyle \frac{g{n}^{\mathrm{2}}}{D^{\mathrm{7}/\mathrm{3}}}}M\sqrt{M^{\mathrm{2}}+N^{\mathrm{2}}}=\mathrm{0},\end{array}\end{array}$$

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial N}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\partial }{\partial x}}\left[{\displaystyle \frac{MN}{D}}\right]+{\displaystyle \frac{\partial }{\partial y}}\left[{\displaystyle \frac{N^{\mathrm{2}}}{D}}\right]+gD{\displaystyle \frac{\partial \mathit{\eta }}{\partial y}}\text{(16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{g{n}^{\mathrm{2}}}{D^{\mathrm{7}/\mathrm{3}}}}N\sqrt{M^{\mathrm{2}}+N^{\mathrm{2}}}=\mathrm{0},\end{array}$$

where η denotes the water level, D denotes the total water level, g denotes the acceleration due to gravity, n denotes the Manning coefficient, and M and N denote the fluxes in the x and y directions, respectively. The governing equations were discretized via the staggered leapfrog scheme (Goto and Ogawa, 1982; UNESCO, 1997). To consider wave height uncertainty, we implemented 25 cases of tsunami numerical simulation for each earthquake source. As detailed in the second chapter, this study focused on the slip amount and the fault depth among many uncertain factors. In each source, the slip amount was varied by ±0.1 times and ±0.05 times with respect to the reference case (five cases) in terms of M_w conversion based on the scaling law, and the fault depth was changed by +2.0, +1.0, −0.5, and −1.0 km with respect to the reference case (five cases) to consider the changes of the slip and the fault depth as uncertainty.

There are a total of 10 earthquake sources; thus, we implemented a total of 250 cases of tsunami numerical simulation nested in four stages of 270, 90, 30, and 10 m in the Japanese plane rectangular coordinate system IX for each simulation and executed the simulation for 3 h from the earthquake occurrence. As an example, Fig. 5 shows the numerical simulation results of nine cases around Oiso and Miura in which the M_w of source 8 is changed to ±0.1 and the fault depth is changed to +2.0 and −1.0 km. As shown in the figure, the distributions of the maximum tsunami wave height vary locally by changing the slip amount and the fault depth, and the effect of the slip amount on the maximum tsunami wave height is more dominant than the fault depth. In addition, while there is a clear positive correlation between the maximum tsunami wave height and slip amount of the earthquake, there is no clear correlation between the maximum tsunami wave height and the fault depth. Figure 6 shows the maximum tsunami wave heights of Miura and Oiso and Pearson's correlation coefficient relative to the tsunami numerical simulation results of each earthquake source. We confirmed that the correlation coefficient corresponded to at least 0.8 in any source; thus the correlation between tsunami wave height of Miura and Oiso was relatively high. The results suggest that we should assess tsunami risk considering the spatial correlation of tsunami wave height between the target points.

Figure 6 Maximum tsunami wave heights simulated from the tsunami numerical simulation at Miura and Oiso and Pearson's correlation coefficients in each earthquake source.

3.3 Construction of response surface

In this section, we construct response surfaces, which indicate maximum wave height at target sites.

With respect to the results of the maximum wave height of the tsunami numerical simulation, we regressed the response surface (Eq. 2) using the least-squares method. The explanatory variables correspond to the fault slip and the fault depth, and the objective variable denotes the maximum wave height at the target sites. We performed the regression analysis based on all combinations of four explanatory variables (2${}^{\mathrm{4}}-\mathrm{1}=\mathrm{15}$ cases) and adopted a response surface with a high coefficient of determination and the minimum Akaike information criterion (AIC) (Akaike, 1974). AIC can compare the quality of a set of statistical models to each other. The best model is the one that has the minimum AIC among all the other models. Table 3 shows the AIC values of 15 case regression analyses for Miura and Oiso, and Table 4 shows the regression coefficients of the response surface where AIC corresponds to the minimum in each earthquake source. For example, Fig. 7a and b show the response surface for the earthquake source 8 (M_w=8.6) with the highest M_w in the Sagami Trough earthquake. The blue circle denotes the maximum wave height obtained from the tsunami numerical simulations, and the red curved surface denotes the response surface. The response surfaces accurately represented the results of the tsunami numerical simulation. The response surfaces are in accordance with Eq. (16) for Oiso and Eq. (17) for Miura as follows:

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }h(S,D)=\mathrm{0.6567}S+\mathrm{0.0459}D-\mathrm{0.5189}{S}^{\mathrm{2}}+\mathrm{0.5147},\text{(18)}& {\displaystyle }& {\displaystyle }h(S,D)=\mathrm{11.1136}S-\mathrm{4.0165}{S}^{\mathrm{2}}-\mathrm{3.1327}.\end{array}$$

We can obtain the frequency distribution of the tsunami wave height by giving a probability distribution function that expresses the uncertainty in the explanatory variable (slip ratio S and fault depth D) of the evaluated response surface and by performing a Monte Carlo simulation.

Figure 7 Response surfaces at (a) Oiso and (b) Miura for source 8 of the Sagami Trough earthquake. The blue circle denotes the maximum wave height obtained from the tsunami numerical simulations, and the red curved surface denotes the response surface.

Table 3 Akaike information criterion (AIC) results of the regression analyses. The regression analyses were performed based on all combinations of four explanatory variables.

Table 4 Regression coefficients of each selected response surface for each earthquake source.

Figure 8 Histograms of tsunami wave height simulated from the response surface at (a) Oiso and (b) Miura for source 8 of the Sagami Trough earthquake.

As reported by Japan Society of Civil Engineers (2002), the estimated variation of M_w of an earthquake of the same magnitude is approximately 0.1. Based on the aforementioned value, we set a normal distribution with an average value of 1.0 and a standard deviation of 0.1 for the slip rate by using the scaling law. With respect to the uncertainty of the fault depth, we also set a normal distribution. The average value was set to 0.0 m, and the standard deviation was set to a random number generated from a lognormal distribution that was obtained from the seismic observation error data from October 2016 to September 2017 (N=305 030) as published by the Japan Meteorological Agency (2017). We used the lognormal distribution with an average of 0.12 km and a standard deviation of 0.65 km. We would like to note that it is essentially necessary to apply a probability distribution that appropriately expresses all possible uncertainties to the explanatory variables of the response surface, but in this study we applied a relatively limited probability distribution as uncertainty since we did not focus on discussing the details of the tsunami wave uncertainty but on the proposed tsunami hazard and risk assessment method using response surface and copulas. Figure 8a and b show the frequency distribution of the tsunami wave height obtained by the aforementioned procedure. By using the response surface method, we can significantly reduce the simulation costs for probabilistic tsunami hazard assessment considering uncertainty.

To ascertain the normality of the frequency distributions, we performed the Kolmogorov–Smirnov test. Table 5 shows the results of p values for each source. In several cases the p values were less than 0.05, thereby indicating that the distribution of the tsunami heights does not necessarily follow a normal distribution.

Table 5 Kolmogorov–Smirnov test results.

3.4 Dependence modeling using copulas

In this section, we estimate appropriate copulas from the results of the tsunami numerical simulation considering uncertainties and evaluate the spatial correlation structure of tsunami wave height between two sites.

As confirmed in the previous section, despite the high linear correlation of the frequency distribution of the tsunami wave height in Miura and Oiso, it is observed that the normality of tsunami wave height for several sources was not secured by the normality test. The Pearson correlation coefficient did not accurately grasp the spatial correlation structure of tsunami wave height, and thus we attempt modeling using a copula. Hereafter, we only illustrate the analysis results of the earthquake source 8 (M_w=8.6) with the largest M_w as an example.

Table 6 shows the results of estimating copulas by maximum likelihood estimation for the distribution obtained by converting the numerical simulation results over [0, 1]. We considered a copula associated with the minimum AIC and Bayesian information criterion (BIC) (Schwarz, 1978) as the best-fit copula. The BIC is more useful in selecting a correct model, while the AIC is more appropriate in finding the best model for predicting future observations. In source 8, the copula with the minimum AIC and BIC corresponded to the Frank copula. We derived the joint distribution of the tsunami wave heights considering the wave height correlation using the Frank copula and the empirical cumulative distributions obtained from the histogram of the tsunami wave height evaluated in the previous section. Figure 9 shows the Frank copula over [0, 1] with 10 000 trials, Fig. 10a and b show the empirical cumulative distributions of tsunami wave height for Oiso and Miura, and Fig. 11a shows the results considering the wave height correlation. The black points denote the results of the Monte Carlo simulation. The number of simulations is 10 000. The red points denote the results of the tsunami numerical simulation using the nonlinear long-wave equation. To compare with this result, Fig. 11b shows the results without considering the wave height correlation. We independently generated the tsunami wave height by using a uniform random number and the cumulative frequency distribution of the tsunami wave height at each site without using a copula. By considering the spatial correlation of the tsunami wave heights using copula, we performed a Monte Carlo simulation that appropriately captures the nonlinear spatial correlation of the tsunami wave height. We clearly showed the usefulness of copula modeling considering the wave height correlation.

Figure 9 Selected Frank copula for source 8.

Figure 10 Empirical cumulative distributions of tsunami wave height (a Oiso and b Miura) for source 8.

Figure 11 Monte Carlo simulation results for source 8. The black points denote the results with 10 000 trials (a) considering and (b) not considering the spatial correlation of tsunami wave heights using the Frank copulas. The red points denote the results calculated from 25 cases of tsunami numerical simulation.

Table 6 Maximum likelihood estimation results of each copula for source 8.

Table 7 shows the result of estimating copulas under the same procedure for other earthquake sources. In the earthquake sources targeted in this study, four types of copula were estimated, namely the rotated Gumbel copula, asymmetric Gumbel copula, Frank copula, and Gumbel copula. The rotated Gumbel copula corresponds to a copula that rotates the ordinary Gumbel copula by 180^∘. For reference purposes, the copulas for all earthquake sources are illustrated in Fig. 12. From the characteristics of the copula mentioned before, there is a tail dependency in the wave heights due to source 1, 2, 3, 5, 7, and 9, but there is no tail dependency in the wave heights due to source 4, 6, 8, and 10. The tail dependency of the wave height could change in various ways under the effects from the relative position of the earthquake sources and the target points, the bottom and land topography.

Figure 12 Estimated optimal copulas distributed on [0, 1]² with 10 000 trials. (a) Rotated Gumbel copula for source 1, (b) asymmetric Gumbel copula for source 2, (c) rotated Gumbel copula for source 3, (d) Frank copula for source 4, (e) rotated Gumbel copula for source 5, (f) Frank copula for source 6, (g) Gumbel copula for source 7, (h) Frank copula for source 8, (i) Gumbel copula for source 9, and (j) Frank copula for source 10.

Table 7 Estimated optimal copulas, copula parameters, and Kendall's τ for each source of the Sagami Trough earthquake.

3.5 Risk assessment results and discussion

In this section, we evaluate the joint distribution of tsunami wave heights and damage probability of target buildings for the entire area of the Sagami Trough earthquake using the occurrence probability weights of each earthquake source.

Figure 13 (a) Joint distribution of tsunami wave height considering wave height correlation and (b) not considering wave height correlation. (c) Joint damage probability for the all sources of the Sagami Trough earthquake. The black points denote the Monte Carlo simulation results with 10 000 trials, and the red points denote the results simulated via tsunami numerical simulations.

Table 8 shows the occurrence probability weights of each source of the Sagami Trough earthquake published by NIED (2017). We first determine the earthquake occurrence source via uniform random numbers using the weights and then evaluate the joint distribution of the tsunami wave heights due to the determined earthquake using the estimated copula. Figure 13 shows the results of evaluation by Monte Carlo simulation with 10 000 trials. Figure 13a shows the joint distribution of the tsunami wave heights considering the spatial correlation of the wave height, and Fig. 13b shows the results without considering the spatial correlation of the tsunami wave height. Furthermore, Fig. 13c shows the joint damage probability of two buildings that transform both axes of tsunami wave heights in Fig. 13b into the damage probability by using the damage function of the steel frame (Suppasri et al., 2013) based on the assumption that a steel building exists at the evaluation target point. Table 9 shows the average value of the aggregate damage probability of two buildings, 95th percentile value, 99th percentile value, and maximum value assuming that the two buildings exhibit the same asset value. Although the expected value of the aggregate damage probability barely changed when compared with that of the no-correlation case, the aggregate damage probability of the 99th percentile value was approximately 1.0 % higher and the maximum value was approximately 3.0 % higher when considering the hazard correlation utilizing the copulas. We clearly showed the significance of considering the spatial correlation structure of tsunami wave height in evaluating tsunami risks for a building portfolio. In this study we only demonstrated the evaluation method for two points, but the effect of the wave height correlation on the evaluation result is expected to increase if more points are targeted.

Table 8 Occurrence probability weights of each source of the Sagami Trough earthquake (NIED, 2017).

Table 9 Tsunami risk assessment results.

4 Conclusion

In this study, we evaluated the aggregate tsunami damage probability of two buildings located at two relatively remote locations based on the frequency distribution of the tsunami height via the response surface method and the spatial correlation of the tsunami height by using various copulas, assuming the occurrence of the Sagami Trough earthquake that significantly affects the metropolitan area in Japan. The 99th percentile value of the aggregate damage probability was approximately 1.0 % higher, and the maximum value was approximately 3.0 % higher in the evaluation considering the spatial correlation of the tsunami wave height when compared with the evaluation without considering the spatial correlation. The results clearly show the significance of considering the spatial correlation of the tsunami hazard in evaluating tsunami risks for a building portfolio and suggest that spatial correlation modeling by copulas is effective in the case wherein nonlinear correlation of the tsunami hazard exists. In addition, the response surface method used in this study significantly reduces the numerical simulation costs for probabilistic tsunami hazard assessment considering uncertainty. In this study, we only focused on the slip amount and fault depth among many tsunami hazard uncertainties, and we evaluated them using the response surface method. It has been reported that the heterogeneity of the slip distribution of the fault has a great influence on tsunami intensity. It is a future issue to evaluate these effects with a response surface method.

The evaluation result was shown for only two buildings, but when an entity evaluates the risk of assets it owns it is assumed that there will be more target sites. It is clear that as the number of target assets increases, the percentile value and maximum value of the aggregate damage of assets become more prominent. Risk assessment that does not consider the spatial correlation of wave heights will lead to the underestimation of the risks held. The basic method shown in this study can be applied even when the number of target assets increases. It is also important to avoid underestimating the assessed risk by considering the wave height correlation using a copula. It is expected that the tsunami risk assessment method for a building portfolio over a wide area as proposed in this study can be used for probabilistic tsunami risk assessment of real-estate portfolios or business continuity plans by parties such as large companies, insurance companies, and real-estate agencies.