**Research article**
22 Nov 2019

**Research article** | 22 Nov 2019

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

Yo Fukutani Shuji Moriguchi Kenjiro Terada Takuma Kotani Yu Otake and Toshikazu Kitano

^{1},

^{2},

^{2},

^{3},

^{4},

^{5}

**Yo Fukutani et al.**Yo Fukutani Shuji Moriguchi Kenjiro Terada Takuma Kotani Yu Otake and Toshikazu Kitano

^{1},

^{2},

^{2},

^{3},

^{4},

^{5}

^{1}College of Science and Engineering, Kanto Gakuin University, Yokohama, 236-8501, Japan^{2}International Research Institute of Disaster Science, Tohoku University, Sendai, 980-8572, Japan^{3}Research and Development Center, Nippon Koei Co., Ltd., Ibaraki, 300-1259, Japan^{4}Faculty of Engineering, Niigata University, Niigata, 950-2181, Japan^{5}Civil Engineering, Nagoya Institute of Technology, Nagoya, 466-8555, Japan

^{1}College of Science and Engineering, Kanto Gakuin University, Yokohama, 236-8501, Japan^{2}International Research Institute of Disaster Science, Tohoku University, Sendai, 980-8572, Japan^{3}Research and Development Center, Nippon Koei Co., Ltd., Ibaraki, 300-1259, Japan^{4}Faculty of Engineering, Niigata University, Niigata, 950-2181, Japan^{5}Civil Engineering, Nagoya Institute of Technology, Nagoya, 466-8555, Japan

**Correspondence**: Yo Fukutani (fukutani@kanto-gakuin.ac.jp)

**Correspondence**: Yo Fukutani (fukutani@kanto-gakuin.ac.jp)

Received: 27 Apr 2019 – Discussion started: 27 May 2019 – Revised: 17 Oct 2019 – Accepted: 19 Oct 2019 – Published: 22 Nov 2019

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.

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.

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.

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:

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:

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.

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*_{1}, …, *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):

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).

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:

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):

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

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.

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.

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):

where “Mo” denotes moment magnitude (Nm), *μ* denotes shear modulus (Pa), *S* denotes slip amount (m), and *A* denotes earthquake source area (m^{2}). *μ* was set to 3.4×10^{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 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:

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.

## 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:

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.

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.

## 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.

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.

## 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.

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.

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.

The earthquake source parameters of the Sagami Trough model used in this study are freely available at http://www.j-shis.bosai.go.jp/map/JSHIS2/download.html?lang=en (last access: 21 November 2019; National Research Institute for Earth Science and Disaster, 2019).

YF conceived and designed the experiments, analyzed the data, and wrote the paper with assistance and input from SM, KT, TK, YO, and TK.

The authors declare that they have no conflict of interest.

We thank two reviewers who provided us valuable comments and helped improve the manuscript. This research was partially supported by funding from the International Research Institute of Disaster Science (IRIDeS) at Tohoku University.

This research has been supported by the International Research Institute of Disaster Science (IRIDeS) at Tohoku University (Tsunami mitigation research 2).

This paper was edited by Ira Didenkulova and reviewed by Elena Suleimani and one anonymous referee.

Akaike, H.: A new look at the statistical model identification, IEEE T. Automat. Control, 19, 716–723, https://doi.org/10.1109/TAC.1974.1100705, 1974.

Annaka, T., Satake, K., Sakakiyama, T., Yanagisawa, K., and Shuto, N.: Logic-tree Approach for Probabilistic Tsunami Hazard Analysis and its Applications to the Japanese Coasts, Pure Appl. Geophys., 164, 577–592, https://doi.org/10.1007/s00024-006-0174-3, 2007.

Boore, D. M., Gibbs, J. F., Joyner, W. B., Tinsley, J. C., and Ponti, D. J.: Estimated ground motion from the 1994 Northridge, California, earthquake at the site of the interstate 10 and La Cienega Boulevard bridge collapse, West Los Angeles, California, Bull. Seismol. Soc. Am., 93, 2737–2751, https://doi.org/10.1785/0120020197, 2003.

Cabinet Office: The study meeting for the Tokyo Inland Earthquakes, available at: http://www.bousai.go.jp/kaigirep/chuobou/senmon/shutochokkajishinmodel/ (last access: 11 May 2018), 2013.

Chang, S. E., Shinozuka, M., and Moore, J. E.: Probabilistic earthquake scenarios: Extending risk analysis methodologies to spatially distributed systems, Earthq. Spect., 16, 557–572, https://doi.org/10.1193/1.1586127, 2000.

Davies, G., Griffin, J., Løvholt, F., Glimsdal, S., Harbitz, C., Thio, H. K., Lorito, S., Basili, R., Selva, J., Geist, E., and Baptista, M. A.: A global probabilistic tsunami hazard assessment from earthquake sources, Geol. Soc. Lond. Spec. Publ., 456, 219–244, https://doi.org/10.1144/SP456.5, 2018.

De Risi, R. and Goda, K.: Probabilistic Earthquaketsunami Hazard Assessment: The First Step Towards Resilient Coastal Communities, Proced. Eng., 198, 1058–1069, https://doi.org/10.1016/j.proeng.2017.07.150, 2017.

Fukutani, Y., Suppasri, A., and Imamura, F.: Stochastic analysis and uncertainty assessment of tsunami wave height using a random source parameter model that targets a Tohoku-type earthquake fault, Stoch. Environ. Res. Risk. A., 29, 1763–1779, https://doi.org/10.1007/s00477-014-0966-4, 2015.

Geist, E. L.: Local tsunamis and earthquake source parameters, Adv. Geophys., 39, 117–209, https://doi.org/10.1016/S0065-2687(08)60276-9, 1999.

Geist, E. L.: Complex earthquake rupture and local tsunamis, J. Geophys. Res., 107, ESE2-1–ESE2-15, https://doi.org/10.1029/2000JB000139, 2002.

Geist, E. L. and Parsons, T.: Probabilistic analysis of tsunami hazards, Nat. Hazards., 37, 277–314, https://doi.org/10.1007/s11069-005-4646-z, 2006.

Geist, E. L. and Yoshioka, S.: Source parameters controlling the generation and propagation of potential local tsunamis along the Cascadia margin, Nat. Hazards, 13, 151–177, https://doi.org/10.1007/BF00138481, 1996.

Genest, C. and Favre, A. C.: Everything You Always Wanted to Know about Copula Modeling but Were Afraid to Ask, J. Hydrol. Eng., 12, 347–368, https://doi.org/10.1061/(ASCE)1084-0699(2007)12:4(347), 2007.

Genest, C. and Segers, J.: Rank-based inference for bivariate extreme-value Copulas, Ann. Stat., 37, 2990–3022, https://doi.org/10.1214/08-AOS672, 2009.

Goda, K. and Hong, H. P.: Estimation of seismic loss for spatially distributed buildings, Earthq. Spect., 24, 889–910, https://doi.org/10.1193/1.2983654, 2008.

Goda, K. and Ren, J.: Assessment of Seismic Loss Dependence Using Copula, Risk Anal., 30, 1076–1091, https://doi.org/10.1111/j.1539-6924.2010.01408.x, 2010.

Goda, K. and Tesfamariam, S.: Multi-variate seismic demand modelling using copulas: Application to non-ductile reinforced concrete frame in Victoria, Canada, Struct. Safety, 56, 39–51, https://doi.org/10.1016/j.strusafe.2015.05.004, 2015.

Goda, K., Mai, P. M., Yasuda, T., and Mori, N.: Sensitivity of tsunami wave profiles and inundation simulations to earthquake slip and fault geometry for the 2011 Tohoku earthquake, Earth Planets Space, 66, 105, https://doi.org/10.1186/1880-5981-66-105, 2014.

González, F. I., Geist, E. L., Jaffe, B., Kânoğlu, U., Mofjeld, H., Synolakis, C. E., Titov, V. V., Arcas, D., Bellomo, D., and Carlton, D.: Probabilistic tsunami hazard assessment at seaside, Oregon, for near and far field seismic sources, J. Geophys. Res.-Oceans, 114, C11023, https://doi.org/10.1029/2008JC005132, 2009.

Goto, C. and Ogawa, Y.: Tsunami numerical simulation with Leap-frog scheme, Tohoku University, Tohoku, p. 52, 1982.

Grezio, A., Babeyko, A., Baptista, M. A., Behrens, J., Costa, A., Davies, G., Geist, E. L., Glimsdal, S., González, F. I., Griffin, J., Harbitz, C. B., LeVeque, R. J., Lorito, S., Løvholt, F., Omira, R., Mueller, C., Paris, R., Parsons, T., Polet, J., Power, W., Selva, J., Sørensen, M. B., and Thio, H. K.: Probabilistic Tsunami Hazard Analysis: Multiple Sources and Global Applications, Rev. Geophys., 55, 1158–1198, https://doi.org/10.1002/2017RG000579, 2017.

Grossi, P. and Kunreuther, H. (Eds): Catastrophe Modeling: A New Approach to Managing Risk, Springer, New York, 2005.

Gumbel, E. J.: Distributions des valeurs extrèmes en plusieurs dimensions, Publ. Inst. Statist. Univ. Paris, 9, 171–173, 1960.

Honjo, Y.: Challenges in geotechnical reliability based design, in: Proceedings of 3rd International Symposium on Geotechnical Safety and Risk, 2–3 June 2011, Munich, Germany, 11–27, 2011.

Hwang, L. S. and Divoky, D.: Tsunami generation, J. Geophys. Res., 75, 6802–6817, https://doi.org/10.1029/JC075i033p06802, 1970.

Japan Meteorological Agency: Centralized processing earthquake source lists, available at: https://hinetwww11.bosai.go.jp/auth/?LANG=ja (last access: 11 May 2018), 2017.

Japan Society of Civil Engineers: Tsunami Assessment Method for Nuclear Power Plants in Japan, available at: http://committees.jsce.or.jp/ceofnp/node/5 (last access: 30 August 2015), 2002.

Joe, H.: Multivariate Models and Dependence Concepts, Chapman & Hall Ltd, Boca Raton, FL, p. 424, 1997.

Kanamori, H.: The energy release in great earthquakes, J. Geophys. Res., 82, 2981–2987, https://doi.org/10.1029/JB082i020p02981, 1977.

Kendall, M.: A New Measure of Rank Correlation, Biometrika, 30, 81–93, https://doi.org/10.1093/biomet/30.1-2.81, 1938.

Kleindorfer, P. R. and Kunreuther, H. C.: Challenges Facing the Insurance Industry in Managing Catastrophic Risks, edited by: Froot, K. A., University of Chicago Press, Chicago, IL, USA, 149–194, https://doi.org/10.7208/chicago/9780226266251.001.0001, 1999.

Kotani, T., Takase, S., Moriguchi, S., Terada, K., Fukutani, Y., Otake, Yu., Nojima, K., and Sakuraba, M.: Numerical-analysis-aided probablistic tsunami hazard evaluation using response surface, J. Japan Soc. Civ. Eng. Ser. A2, 72, 58–69, https://doi.org/10.2208/jscejam.72.58, 2016.

Løvholt, F., Pedersen, G., Bazin, S., Kuhn, D., Bredesen, R. E., and Harbitz, C.: Stochastic analysis of tsunami runup due to heterogeneous coseismic slip and dispersion, J. Geophys. Res., 117, C03047, https://doi.org/10.1029/2011JC007616, 2012.

Løvholt, F., Griffin, J., Salgado-Gálvez, M.: Tsunami Hazard and Risk Assessment on the Global Scale, in: Encyclopedia of Complexity and Systems Science, edited by: Meyers, R., Springer, Berlin, Heidelberg, 1–34, https://doi.org/10.1007/978-3-642-27737-5_642-1, 2015.

Nelsen, R. B.: An Introduction to Copulas, Springer-Verlag, New York, p. 218, https://doi.org/10.1007/978-1-4757-3076-0, 1999.

Ng, M. K., Leblond, P. H., and Murty, T. S.: Simulation of tsunamis from great earthquakes on the Cascadia subduction zone, Science, 250, 1248–1251, https://doi.org/10.1126/science.250.4985.1248, 1991.

NIED – National Research Institute for Earth Science and Disaster Resilience: Japan Seismic Hazard Information Station, available at: http://www.j-shis.bosai.go.jp/map/ (last access: 11 May 2018), 2017.

Okada, Y.: Surface deformation due to shear and tensile faults in a half-space, Bull. Seismol. Soc. Am., 75, 1135–1154, 1985.

Park, H. and Cox, D. T.: Probabilistic assessment of near-field tsunami hazards: Inundation depth, velocity, momentum flux, arrival time, and duration applied to Seaside, Oregon, Coast. Eng., 117, 79–96, https://doi.org/10.1016/j.coastaleng.2016.07.011, 2016.

Park, J., Bazzurro, P., and Baker, J. W.: Modeling spatial correlation of ground motion intensity measures for regional seismic hazard and portfolio loss estimation, in: Tenth International Conference on Application of Statistic and Probability in Civil Engineering (ICASP10), Tokyo, Japan, 2007.

Pelayo, A. M. and Wiens, D. A.: Tsunami earthquakes: slow thrust-faulting events in the accretionary wedge, J. Geophys. Res., 97, 15321–15337, https://doi.org/10.1029/92JB01305, 1992.

Salgado-Gálvez, M. A., Zuloaga-Romero, D., Bernal, G. A., Mora, M. G., and Cardona, O. D.: Fully probabilistic seismic risk assessment considering local site effects for the portfolio of buildings in Medellín, Colombia, Bull. Earth Eng., 12, 671–695, https://doi.org/10.1007/s10518-013-9550-4, 2014.

Salvadori, G., Durante, F., Michele, C. D., Bernardi, M., and Petrella, L.: A multivariate copula-based framework for dealing with hazard scenarios and failure probabilities, Water Resour. Res., 52, 3701–3721, https://doi.org/10.1002/2015WR017225, 2016.

Scheingraber, C. and Käser, M.: Spatial Seismic Hazard Variation and Adaptive Sampling of Portfolio Location Uncertainty in Probabilistic Seismic Risk Analysis, Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2019-110, in review, 2019.

Schwarz, G. E.: Estimating the dimension of a model, Ann. Stat., 6, 461–464, 1978.

Sklar A. W.: Fonctions de répartition à *n* dimension et leurs marges, Publications de l'Institut de Statistique de l'Université de Paris, 8, 229–231, 1959.

Song, Y. T., Ji, C., Fu, L. L., Zlotnicki, V., Shum, C. K., Yi, Y., and Hjorleifsdottir, V.: The 26 December 2004 tsunami source estimated from satellite radar altimetry and seismic waves, Geophys. Res. Lett., 32, L20601, https://doi.org/10.1029/2005GL023683, 2005.

Suppasri, A., Mas, E., Charvet, I., Gunasekera, R., Imai, K., Fukutani, Y., Abe, Y., and Imamura, F.: Building damage characteristics based on surveyed data and fragility curves of the 2011 Great East Japan tsunami, Nat. Hazards, 66, 319–341, https://doi.org/10.1007/s11069-012-0487-8, 2013.

Tawn, J. A.: Bivariate extreme value theory: Models and estimation, Biometrika, 75, 397–415, https://doi.org/10.2307/2336591, 1988.

Thio, H. K., Somerville, P. G., and Polet, J.: Probabilistic Tsunami Hazard in California, College of Engineering, University of California, Los Angeles, CA, USA, 2010.

UNESCO: IUGG/IOC Time Project: Numerical method of tsunami simulation with the leap-flog scheme, IOC Manuals and Guides No. 35, Paris, France, available at: http://www.vliz.be/imisdocs/publications/ocrd/269372.pdf (last access: 21 November 2019), 1997.

Wang, M. and Takada, T.: Macro-spatial correlation model of seismic ground motions, in: Proceedings of ICOSSAR'05, Millpress, Rotterdam, 353–360, 2005.

Ward, S. N.: On tsunami nucleation: II. An instantaneous modulated line source, Phys. Earth Planet. Int., 27, 273–285, https://doi.org/10.1016/0031-9201(82)90057-7, 1982.

Whitmore, P. M.: Expected tsunami amplitudes and currents along the North American coast for Cascadia subduction zone earthquakes, Nat. Hazards, 8, 59–73, https://doi.org/10.1007/BF00596235, 1993.