**Bayesian Variable Selection for Pareto Regression Models with Latent Multivariate Log Gamma Process with Applications to Earthquake Magnitudes**

#### **Hou-Cheng Yang 1, Guanyu Hu 2,\* and Ming-Hui Chen <sup>2</sup>**

<sup>1</sup> Department of Statistics, Florida State University, Tallahassee, FL 32306, USA; hou-cheng.yang@stat.fsu.edu

<sup>2</sup> Department of Statistics, University of Connecticut, Storrs, CT 06269, USA; ming-hui.chen@uconn.edu

**\*** Correspondence: guanyu.hu@uconn.edu

Received: 10 January 2019; Accepted: 9 April 2019; Published: 12 April 2019

**Abstract:** Generalized linear models are routinely used in many environment statistics problems such as earthquake magnitudes prediction. Hu et al. proposed Pareto regression with spatial random effects for earthquake magnitudes. In this paper, we propose Bayesian spatial variable selection for Pareto regression based on Bradley et al. and Hu et al. to tackle variable selection issue in generalized linear regression models with spatial random effects. A Bayesian hierarchical latent multivariate log gamma model framework is applied to account for spatial random effects to capture spatial dependence. We use two Bayesian model assessment criteria for variable selection including Conditional Predictive Ordinate (CPO) and Deviance Information Criterion (DIC). Furthermore, we show that these two Bayesian criteria have analytic connections with conditional AIC under the linear mixed model setting. We examine empirical performance of the proposed method via a simulation study and further demonstrate the applicability of the proposed method in an analysis of the earthquake data obtained from the United States Geological Survey (USGS).

**Keywords:** earthquake hazard; DIC; CPO; model selection

#### **1. Introduction**

The earthquake magnitude data has become increasingly popular over the last decade. Statistical models for earthquake have been proposed since 1800s. Since large earthquakes are rare, it is difficult to fit simple linear models. Many different parametric models (Gamma model, Weibull model) have been considered to analyze earthquake magnitudes, but some earthquakes with very small magnitudes are not reported by seismic centers. The Pareto-type distribution is a popular choice for analyzing earthquake magnitudes data (e.g., [1–3]), as the Pareto distribution is a heavy-tailed distribution with a lower threshold. In statistical analysis, a regression model is used to connect dependent covariates of earthquakes to the magnitude of the earthquake. A generalized linear model strategy can be used for the Pareto regression. Existing seismology literatures pay less attention to spatially dependent structure on earthquake magnitudes. They just built simple linear regression models or generalized linear models to explore covariates effects on earthquake magnitudes [4]. Hu and Bradley [5] proposed using the Pareto regression with spatial random effects for earthquake magnitudes, but they did not consider the model selection problems. In order to have more explicit understanding of dependent covariates of earthquake magnitudes, variable selection approaches should be considered in a Pareto regression model with spatial random effects.

Variable selection and Bayesian statistics have received widespread attention and become increasingly important tools in the field of environment and ecology [6,7]. For hierarchical spatial model, it is difficult to do inference for latent variables. Bayesian approach provides a convenient way for estimating latent variables in hierarchical models. Compared with the frequentist approach, a Bayesian approach can bring some prior information on parameters of the model. It is an important part of a statistical analysis. In practice, we may want to measure how good a model is for answering a certain question or comparing different models to see which model is best suited. There are many popular variable selection criteria, including Akaike's information criterion (AIC) [8] and Bayesian information criterion (BIC) [9], Bayes factor, conditional predictive ordinate (CPO) [10,11], L measure [12], and the deviance information criterion (DIC) [13]. Chen et al. [14] provide the connections between these popular criteria for variable subset selection under generalized linear models. However, there are some difficulties for Bayesian variable selection to carry out because of the challenge in assigning prior distributions for the parameters. In order to tackle this issue, we consider the multivariate log-Gamma distribution (MLG) based on Bradley et al. [15], which is conjugate with the Pareto distribution [5]. Hence, the Bayesian approach to variable selection is straightforward for our model. Consequently, we use CPO and DIC criteria to carry out Bayesian variable selection for Pareto regression models due to the performance of the conjugate priors (see [16], for a discussion).

Both CPO and DIC are criteria-based methods and they have some advantage over other criteria. Compared with regularized estimation approach, these two criteria consider goodness of fit of the candidate models. Furthermore, compared with negative log probability density or RMSE for predictions, these two criteria consider the model complexity. Like the AIC or BIC, these two criteria compromise the tradeoff between the goodness of fit and model complexity. The CPO provides a site-specific model fit metric that can be used for exploratory analysis and can be combined at the site to generate a logarithm pseudo marginal likelihood (LPML) as an overall model fit measure. The CPO is based on leave-one-out-cross-validation. It estimates the probability of observing data on one particular location in the future if after having already observed data. The LPML is a leave-one-out cross-validation with log likelihood as the criteria which can be easily obtained from an Markov chain Monte Carlo (MCMC) output (see [17]). More details about two criteria will be discussed in Section 2.2. The major contribution of this paper is that we introduce two Bayesian model selection criteria in generalized linear model with spatial random effects. Furthermore, we exam the relationship between the two criteria with conditional AIC (cAIC) in random effects model. Other than the variable selection problem in regression model, our criteria can also be used in model selection in the presence of spatial random effects. In general, our proposed criteria can select important covariates and random effects model simultaneously.

The remaining sections of this article are organized as follows. Section 2 introduces our proposed statistical model, and review two Bayesian model assessment Criteria including LPML and DIC [13]. In Sections 3 and 4, we present MCMC scheme and a simulation study for two scenarios, and use two criteria to select true model. In Section 5, we carry out a detailed analysis of the US earthquake dataset from United States Geological Survey (USGS) and use two criteria to select the best model(s). Finally, Section 6 contains a brief summary of this paper. For ease of exposition all proofs are given Appendix A.

#### **2. Methodology**

#### *2.1. Pareto Regression with Spatial Random Effects*

In many regression problems, normality may not be always held. Generalized linear models allow a linear regression model to connect the response variable with a proper link function. For some heavy tailed data with minimum value, it is common to use the Pareto model to fit these data. From the expression of Gutenberg–Richter law, it is possible to derive a relationship for the logarithm of the probability to exceed some given magnitude. The standard distribution used for seismic moment is the Pareto distribution. The Pareto distribution has a natural threshold. In practice, people do not take more consideration on "micro" (magnitude from 1–1.9) or "minor" (magnitude from 2–2.9) earthquakes. Compared with exponential distribution, Pareto distribution is a heavy tailed distribution. Heavy tailed distributions tend to have many outliers with very high values. The heavier the tail, the larger

the probability that you will get one or more disproportionate values in a sample. In earthquake data, most recorded earthquakes have a magnitude around 3–5, but sometime there will have some significant earthquakes with large magnitude. Hu [5] used Pareto regression to model earthquake magnitudes, since the Pareto distribution is a heavy tailed distribution with a threshold. Earthquake magnitude data also has a threshold, since people consider earthquake only over a certain magnitude. Based on the generalized linear model setting, we can build Pareto regression model as

$$f(z) = \exp(\mu(\mathfrak{s})) z\_{\mathfrak{m}}^{\exp(\mu(\mathfrak{s}))} z^{-1 - \exp(\mu(\mathfrak{s}))} \quad z \ge z\_{\mathfrak{m}\_f} \tag{1}$$

where *<sup>s</sup>* ∈ *<sup>D</sup>* ⊂ R<sup>2</sup> is a spatial location, *<sup>μ</sup>*(*s*) = *<sup>β</sup>*<sup>0</sup> + *<sup>β</sup>*1*X*1(*s*) + ... + *<sup>β</sup>pXp*(*s*), *Xi*(*s*) is *<sup>i</sup>*-th covariate on location *s* and *zm* is the minimum value of the response variable. Under this model, the log shape parameter is modeled with a fixed effects term.

The model in Equation (1) does not include spatial random effects. Consequently, it is implicitly assumed that *α*(*s*) and *α*(*w*) are independent for *s* = *w*. But for many spatial data, it is not realistic to assume that *α*(*s*) and *α*(*w*) are independent. We can add the latent Gaussian process in the log-linear model so that the generalized linear model becomes a generalized linear mixed model (GLMM). Specifically, we assumed

$$\log \{ a(\mathbf{s}) \} = \mathcal{J}^{\prime} X(\mathbf{s}) + w(\mathbf{s}) \quad \mathbf{s} \in D,\tag{2}$$

$$\mathcal{W} \sim N\{\mathbf{0}, \sigma\_w^2 H(\phi)\},\tag{3}$$

where *W* is an *n*-dimensional vector of (*w*(*s*1), ... , *w*(*sn*)) , *H*(*φ*) is a *n* × *n* spatial correlation matrix, and {*s*1, ... , *sn*} ∈ *D* are the observed spatial locations. The natural strategy to consider spatial correlation is to use in light of Tobler's first law that "near things are more related than distant things" [18]. Spatial random effects allow one to leverage information from nearby locations. Latent Gaussian process models have become a standard method for modeling spatial random effects [19]. Based on Gaussian process structure, the nearby observations will have higher correlation.

For the latent Gaussian process GLMM, we can build the following hierarchical model:

$$\begin{aligned} \text{Data Model } &: \mathcal{Z}(s\_i) | \mathcal{W}, \boldsymbol{\theta}, \sigma^2, \boldsymbol{\phi}^{\text{ind}} \sim \text{Pareto}(\mathcal{Z}\_{\text{w}}, \boldsymbol{\epsilon}^{\text{ul}(s\_i)}); \ i = 1, \dots, n \\ \text{Process Model } &: \mathbb{W} | \boldsymbol{\phi}, \sigma^2\_w \sim \text{MVN} \{ \mathbf{0}, \sigma^2\_w \boldsymbol{\mu}(\boldsymbol{\phi}) \} \\ \text{Parameter Model } & 1: \boldsymbol{\theta}\_j \sim N(0, \sigma^2\_j); \ j = 1, \dots, p \\ \text{Parameter Model } & 2: \sigma^2\_j \sim \text{IG}(a\_1, b\_1); \ j = 1, \dots, p \\ \text{Parameter Model } & 3: \sigma^2\_w \sim \text{IG}(a\_2, b\_2) \\ \text{Parameter Model } & 4: \boldsymbol{\phi} \sim \text{IG}(a\_3, b\_3). \end{aligned}$$

where "IG" is a shorthand for inverse gamma, "MVN" is a shorthand for multivariate normal, and "N" is a shorthand for a univariate normal distribution. For the Pareto regression model, the normal prior is not conjugate. A proper conjugate prior for the Pareto regression will facilitate the development of an efficient computational algorithm. Chen and Ibrahim [16] proposed a novel class of conjugate priors for the family of generalized linear model. But they did not show the connection between their conjugate prior and gaussian prior. Bradley et al. [20] proposed the multivariate log-gamma distribution as a conjugate prior for Poisson spatial regression model and established a connection between a multivariate log-gamma distribution and a multivariate normal distribution. The multivariate log-gamma distribution is an attractive alternative prior for the Pareto regression model due to its conjugacy.

We now present the multivariate log-gamma distribution from Bradley et al. [20]. We define the *n*-dimensional random vector *γ* = (*γ*1, ..., *γn*) , which consists of *n* mutually independent log-gamma random variables with shape and scale parameters organized into the *n*-dimensional vectors *α* ≡ (*α*1, ..., *αn*) , and *κ* ≡ (*κ*1, ..., *κn*) , respectively. Then define the *n*-dimensional random vector *q* as follows:

$$
\mathfrak{q} = \mathfrak{p} + \mathcal{V}\gamma,\tag{5}
$$

where *<sup>V</sup>* ∈ R*<sup>n</sup>* × R*<sup>n</sup>* and *<sup>μ</sup>* ∈ R*n*. Bradley et al. [20] called *<sup>q</sup>* the multivariate log-gamma random vector. The random vector *q* has the following probability density function:

$$f(q|\mathbf{c}, \mathbf{V}, \boldsymbol{\mu}, \kappa) = \frac{1}{\det(\mathbf{V})} \left( \prod\_{i=1}^{m} \frac{1}{\Gamma(a\_i) \kappa\_i^{a\_i}} \right) \exp[\mathbf{a}^\prime \boldsymbol{V}^{-1} (\boldsymbol{q} - \boldsymbol{\mu}) - \kappa^{(-1)'} \exp\{\boldsymbol{V}^{-1} (\boldsymbol{q} - \boldsymbol{\mu})\}]; \quad q \in \mathcal{R}^n, \tag{6}$$

where "det" represents the determinant function. We use "MLG (*μ*, *V*, *α*, *κ*)" as a shorthand for the probability density function in Equation (6).

According to Bradley et al. [20], the latent Gaussian process is a special case of the latent multivariate log-gamma process. If *β* has a multivariate log-gamma distribution MLG(**0**, *<sup>α</sup>*1/2*V*, *<sup>α</sup>***1**, 1/*α***1**). When *<sup>α</sup>* → <sup>∞</sup>, *<sup>β</sup>* will converge in distribution to the multivariate normal distribution vector with mean **0** and covariance matrix *VV* . *α* = 10,000 is sufficiently large for this approximation. MLG model is a more saturated model than Gaussian process model. For the Pareto regression model, the MLG process is more computationally efficient than the Gaussian process. In following hierarchical model, we refer to *β* and *W* as following an MLG distribution with *q*, **0***<sup>p</sup>* and **0***<sup>n</sup>* being the first parameter of MLG corresponding to *μ*, and **Σ**1/2 *<sup>W</sup>* and **<sup>Σ</sup>**1/2 *<sup>β</sup>* are the second parameter of MLG like *V*.

In order to establish conjugacy, we build a spatial GLM with latent multivariate log gamma process as follows:

$$\begin{aligned} \textbf{Data Model}: & Z(s\_i) | W, \boldsymbol{\theta}, \boldsymbol{\sigma}^2, \boldsymbol{\phi}^{\text{ind}} \in \text{Pareto}(Z\_m, \boldsymbol{\epsilon}^{\mu(s\_i)}); \ i = 1, \ldots, n \\ \textbf{Process Model}: & W | \boldsymbol{\phi}, \boldsymbol{\sigma}\_w \sim \text{MLG}(\mathbf{0}\_n, \boldsymbol{\Sigma}\_W^{1/2}, a\_W \mathbf{1}\_n, \kappa\_W \mathbf{1}\_n) \\ \textbf{Parameter Model}: & \mathbf{1}: \boldsymbol{\mathcal{B}} \sim \text{MLG}(\mathbf{0}\_p, \boldsymbol{\Sigma}\_\boldsymbol{\beta}^{1/2}, a\_\beta \mathbf{1}\_p, \kappa\_\beta \mathbf{1}\_p) \\ \textbf{Parameter Model}: & \mathbf{2}: \boldsymbol{\sigma}^2 \sim \text{IG}(a\_1, b\_1); \\ \textbf{Parameter Model}: & \mathbf{3}: \boldsymbol{\sigma}\_w^2 \sim \text{IG}(a\_2, b\_2) \\ \textbf{Parameter Model}: & \mathbf{4}: \boldsymbol{\phi} \sim \text{IG}(a\_3, b\_3). \end{aligned}$$

where *Zm* defined baseline, *μ*(*si*) = *X*(*si*)*β* + *W*, **Σ***<sup>W</sup>* = *σ*<sup>2</sup> *<sup>w</sup>H*(*φ*), **Σ***<sup>β</sup>* = *σ*2diag(*p*), *α<sup>W</sup>* > 0, *αβ* > 0, *κ<sup>W</sup>* > 0, and *κβ* > 0.

#### *2.2. Bayesian Model Assessment Criteria*

In this section, we consider two Bayesian model assessment criteria, DIC and LPML. In addition, we introduce the procedure to calculate DIC and LMPL for the Pareto regression model with spatial random effects. Let *β*(*M*) denote the vector of regression coefficient under the full model *M*. Also let *β*(*m*) and *β*(−*m*) denote the corresponding vectors of regression parameters included and excluded in the subset model *m*. Then, *β*(*M*) = *β* = ((*β*(*m*)) ,(*β*(−*<sup>m</sup>*)) ) holds for all *m*, and *β*(−*M*) = ∅.

#### 2.2.1. DIC

The deviance information criterion is defined as

$$\text{DIC} = \text{Dev}(\theta) + 2p\_{\text{D}\_{\text{L}}} \tag{8}$$

where Dev(¯ *<sup>θ</sup>*) is the deviance function, *pD* <sup>=</sup> Dev(*θ*) <sup>−</sup> Dev(¯ *θ*) is the effective number of model parameters, and ¯ *θ* is the posterior mean of parameters *θ*, and Dev¯ (*θ*) is the posterior mean of Dev(*θ*). To carry out variable selection, we specify the deviance function as

$$\text{Dev}(\mathcal{J}^{(m)}) = -2 \sum\_{i=1}^{n} \log(f(\mathcal{J}^{(m)}) | D\_i), \tag{9}$$

where *Di* = (*Yi*, *Xi*, *W*ˆ *<sup>i</sup>*), *f*(.) is the likelihood function in Equation (7), *W*ˆ *<sup>i</sup>* is the posterior mean of the spatial random effects on location *si*, *β*(*m*) is the vector of regression coefficient under the *m*-th model. In this way, the DIC criterion is given by

$$\text{DIC}\_{\text{ll}} = \text{Dev}(\vec{\mathcal{S}}^{(m)}) + 2p\_{\text{D}}^{(m)},\tag{10}$$

where

$$\mathcal{D}p\_D^{(m)} = \text{Dev}(\bar{\hat{\mathcal{B}}}^{(m)}) - \text{Dev}(\bar{\hat{\mathcal{B}}}^{(m)}),\tag{11}$$

where *<sup>β</sup>*¯(*m*) <sup>=</sup> *<sup>E</sup>*[*β*(*m*)|*D*], and Dev¯ (*βm*) = *E*[Dev(*β*(*m*))].

#### 2.2.2. LPML

In order to calculate the LPML, we need to calculate CPO first [14]. Then LPML can be obtained as

$$\text{LPML} = \sum\_{i=1}^{n} \log(\text{CPO}\_i)\_\prime \tag{12}$$

where CPO*<sup>i</sup>* is the CPO for the *i*-th subject.

Let *D*(−*i*) denote the observation data with the *i*-th observation deleted. The CPO for the *i*-th subject is defined as

$$\text{CPO}\_{i} = f(\text{Y}\_{i}|\mathbf{X}\_{i}, D\_{(-i)}) = \int f(\text{Y}\_{i}|\mathbf{X}\_{i}, \mathfrak{F}) \pi(\mathfrak{F}|D\_{(-i)}) d\mathfrak{F},\tag{13}$$

where *<sup>π</sup>*(*β*|*D*(−*<sup>i</sup>*)) is the posterior distribution based on the data *<sup>D</sup>*(−*i*).

From Chapter 10 of Chen et al. [21], CPO in (13) can be rewritten as

$$\text{CPO}\_{i} = \frac{1}{\int \frac{1}{f(y\_{i}|\mathcal{G}, \mathcal{W}, \mathbf{X}\_{i})} \pi(\mathcal{G}|D)d\mathcal{B}}. \tag{14}$$

A popular Monte Carlo estimate of CPO using Gibbs samples form the posterior distribution is given as *<sup>D</sup>* instead of *<sup>D</sup>*(−*i*). Letting {*βb*, *<sup>b</sup>* = 1, ··· , *<sup>B</sup>*} denote a Gibss sample of *<sup>β</sup>* from *<sup>π</sup>*(*β*|*D*) and using (14), a Monte Carlo estimate of CPO−<sup>1</sup> *<sup>i</sup>* is given by

$$\begin{array}{l} \text{C-P-O}\_{i} \quad \text{is given by} \\\\ \left(\widehat{\text{CPO}}\_{i}\right)^{-1} = \frac{1}{B} \sum\_{b=1}^{B} \frac{1}{f(\text{Y}\_{i}|\mathfrak{F}\_{b}, \text{X}\_{i}, \widehat{\text{W}}\_{i})}. \end{array} \tag{15}$$

So the LPML defined as

$$\text{LPML}\_{\text{ll}} = \sum\_{i=1}^{n} \log(\widehat{\text{CPO}}\_{i}).\tag{16}$$

In the context of variable selection, we select a subset model, which has the largest LPML value and/or the smallest DIC value. In practice, if we have two different results based on two criteria, we will choose both models which were selected by two criteria as the best models. In addition, we can

do more diagnostics for the two candidate models. DIC compromises the goodness of fit and the complexity of the model. The CPO is based on leave-one-out-cross-validation. The LPML, the sum of the log CPO's, is an estimator for the log marginal likelihood.

#### *2.3. Analytic Connections between Bayesian Variable Selection Criteria with Conditional AIC for the Normal Linear Regression with Spatial Random Effects*

The Akaike information criterion (AIC) has been applied to choose candidate models in the mixed-effects model by integrating out the random effects. A conditional AIC was proposed to be used for the linear mixed-effects model [22] under the assumption that the variance-covariance matrix of random effects is known. Under the this assumption, we establish analytic connections of DIC and LPML we proposed in Section 2.3 with cAIC. We have the following linear regression model with spatial random effects:

$$\mathbf{y}\_{i} = \mathbf{X}\_{i}\boldsymbol{\mathfrak{B}} + \boldsymbol{w}\_{i} + \boldsymbol{\varepsilon}\_{i} \tag{17}$$

where *β* is a *p* × 1 vector of fixed effects, *wi* is spatial random effects for individual *i*. The cAIC is defined as:

$$\text{cAIC} = -2\text{log}(\hat{\beta}|\mathbf{X}, \mathbf{y}, \mathbf{\hat{w}}) + 2k,\tag{18}$$

where *X* is with full rank *k*. Having the MLE of *β*, we can have

$$\text{cAIC} = -n \log(\frac{1}{2\pi\sigma^2}) + \frac{1}{\sigma^2} \text{SSE} + 2k,\tag{19}$$

where SSE = (*y* − *y***ˆ**) (*y* − *y***ˆ**), *y***ˆ** = (*y*ˆ1, ..., *y*ˆ*n*) , *y*ˆ*<sup>i</sup>* = *Xiβ***ˆ** + *w*ˆ*i*.

From [14], we can have DIC and LPML for the linear regression model with spatial random effects as follows

$$\text{DIC}(a\_0) = -n \log(\frac{1}{2\pi\sigma^2}) + \frac{1}{\sigma^2} \text{SSE}^\* + \frac{2(1+a\_0)}{1+2a\_0} 2k,\tag{20}$$

and

$$\text{LPML}(a\_0) = -n \log(\frac{1}{2\pi\sigma^2}) + \frac{1}{\sigma^2} \text{SSE}^\* + \frac{(1+a\_0)}{1+2a\_0}k + R,\tag{21}$$

where SSE<sup>∗</sup> is calculated by posterior mean, *a*<sup>0</sup> = 0 with conjugate prior for likelihood model, *<sup>R</sup>* <sup>=</sup> <sup>−</sup>2(1+*a*0)<sup>2</sup> <sup>1</sup>+2*a*<sup>0</sup> *R*∗, *R*<sup>∗</sup> is the remainder of Taylor expansion. So in the conjugate prior condition, our proposed Bayesian variable selection criterion is similar with cAIC for the linear regression model with spatial random effects.

#### **3. MCMC Scheme**

The algorithm requires sampling the all parameters in turn from their respective full conditional distributions. We assume that *β*, *W* are independent a priori. We further assume *β* ∼ MLG(**0***p*, **Σ**1/2 *<sup>β</sup>* , *αβ***1***p*, *κβ***1***p*) and *<sup>W</sup>* <sup>∼</sup> MLG(**0***n*, **<sup>Σ</sup>**1/2 *<sup>W</sup>* , *αW***1***n*, *κW***1***n*). Thus, sampling from *p*(*β*|·) ∝ exp *α <sup>β</sup>Hββ* − *κ <sup>β</sup>* exp(*Hββ*) and *p*(**W**|·) ∝ exp *α <sup>W</sup> HWW* − *κ <sup>W</sup>* exp(*HWW*) is straightforward. For **Σ***<sup>W</sup>* = *σ*<sup>2</sup> *<sup>w</sup>H*(*φ*) and **<sup>Σ</sup>***<sup>β</sup>* = *<sup>σ</sup>*2diag(*p*), we assume that *<sup>σ</sup>*<sup>2</sup> ∼ IG(*a*1, *<sup>b</sup>*1), *<sup>σ</sup>*<sup>2</sup> *<sup>w</sup>* ∼ IG(*a*2, *b*2) and *<sup>φ</sup>* <sup>∼</sup> IG(*a*3, *<sup>b</sup>*3), that is, *<sup>p</sup>*(*σ*2|*a*1, *<sup>b</sup>*1) <sup>∝</sup> MLG(**0**, **<sup>Σ</sup>**1/2 *<sup>β</sup>* , *αβ***1***p*, *κβ***1***p*) × IG(*a*1, *<sup>b</sup>*1), *<sup>p</sup>*(*σ*<sup>2</sup> *<sup>w</sup>*|*a*2, *b*2) ∝ MLG(**0**, **Σ**1/2 *<sup>W</sup>* , *<sup>α</sup>w***1***n*, *<sup>κ</sup>w***1***n*) <sup>×</sup> IG(*a*2, *<sup>b</sup>*2), and *<sup>p</sup>*(*φ*|*a*3, *<sup>b</sup>*3) <sup>∝</sup> MLG(**0***n*, **<sup>Σ</sup>**1/2 *<sup>W</sup>* , *αw***1***n*, *κw***1***n*) × IG(*a*3, *b*3). The sampling scheme for these three parameters is not straightforward. We use a Metropolis–Hasting algorithm to sampling from three parameters. The other difficulty is how to compute the

log-determinant of a matrix. Because we are using a log-likelihood function, the formula for the log-likelihood involves the expression log(det(Σ*β*)) or log(det(Σ*W*)). To compute the logarithm of a determinant, we encourage not try to compute the determinant itself. Instead, computing the log-determinant directly. For a matrix with a large determinant, the computation of the log-determinant will usually be achieved, however, the computation of the determinant might cause a numerical error. The method is given by

$$\begin{aligned} \log(\det(\Sigma\_{\beta})) &= 2 \ast \sum (\log(\text{diag}(L\_{\beta}))),\\ \log(\det(\Sigma\_{W})) &= 2 \ast \sum (\log(\text{diag}(L\_{W}))), \end{aligned} \tag{22}$$

where the *L<sup>β</sup>* is the Cholesky root of matrix Σ*β*, *LW* is the Cholesky root of matrix Σ*W*, and "diag" denotes a column vector whose elements are the elements on the diagonal of matrix. The derivative details for the full conditional distributions given in Appendix A.

Note that *αβ*, *κβ*, *αW*, *κW*, *a*1, *b*1, *a*2, *b*2, *a*3, and *b*<sup>3</sup> are prespecified hyperparameters. In this article, we use *αβ* = 10,000, *κβ* = 0.0001, *α<sup>W</sup>* = 1, *κ<sup>W</sup>* = 1, *a*<sup>1</sup> = 1, *b*<sup>1</sup> = 1, *a*<sup>2</sup> = 1, *b*<sup>2</sup> = 1, *a*<sup>3</sup> = 1 and *b*<sup>3</sup> = 1. For more flexibility, we can also assume *αβ*, *κβ*, *α<sup>W</sup>* and *κ<sup>W</sup>* each following a Gamma distribution with suitable hyperparameters.

#### **4. Simulation Study**

The spatial domain for the two simulation studies are chosen to be D ∈ [0, 50] × [0, 50]. The locations *si* is selected uniformly over *D* ∀*i* = 1 ... 50. We present the two different simulation settings and generate 100 replicate data sets for each scenario. We assume *β* = (*β*1, *β*2, *β*3) so that we have seven candidate models. We generate *W* from a multivariate normal distribution with mean zero and covariance **Σ***<sup>W</sup>* = *H*(*φ*). We set *φ* = 1 and fix *σ*<sup>2</sup> = 1 in both Simulations 1 and 2. We generate the elements of *X*(*si*) independently from the uniform distribution U(0,1). We define the baseline threshold (scale parameter) equal to three in both simulations.

#### *4.1. Simulation for the Connection between Multivariate Log Gamma and Multivariate Normal Distribution*

In this section, we examine the connection between the multivariate log-gamma distribution and the multivariate normal distribution. First, we draw the quantile-quantile (QQ)-plot in Figure 1 to show the normality of *q* generated from MLG(**0**, *α*1/2*V*, *α***1**, 1/*α***1**), when *α* = 10,000. In addition, we use the Kolmogorov–Smirnov test to examine the connection for one dimensional data. We use a multivariate two-sample test [23] for multivariate dimensional data. We generated one data set of size 100 from the multivariate log-gamma distribution and another data set of size 100 from the multivariate normal distribution and then calculated the *p*-value from the multivariate two-sample test for comparing these two data sets. Then, we repeated this process 1000 times. We found that 992 out of these 1000 *p*-values were larger than the significance level of 0.05. That is, in 992 of 1000 times, we did not reject the null hypothesis that the two samples were drawn from the same distribution.

#### *4.2. Simulation for Estimation Performance*

In this simulation study, our goal was to examine the estimation performance of the hierarchical model. We set *β* = (1, 1, 1). We estimated the parameters in this simulation and report the bias (bias = 1 *<sup>m</sup>* <sup>∑</sup>*<sup>m</sup> j*=1(*β*(*j*) *<sup>i</sup>* − *β*<sup>∗</sup> *<sup>i</sup>* )), the standard error (SE) (SE = 1 *<sup>m</sup>* <sup>∑</sup>*<sup>m</sup> j*=1(*β*(*j*) *<sup>i</sup>* <sup>−</sup> *<sup>β</sup>*¯ *i*)2 1/2 , where *β*¯ *<sup>i</sup>* = <sup>1</sup> *<sup>m</sup>* <sup>∑</sup>*<sup>m</sup> <sup>j</sup>*=<sup>1</sup> *<sup>β</sup>*(*j*) *<sup>i</sup>* ), and the mean square error (MSE) (MSE = <sup>1</sup> *<sup>m</sup>* <sup>∑</sup>*<sup>m</sup> j*=1(*β*(*j*) *<sup>i</sup>* − *β*<sup>∗</sup> *<sup>i</sup>* )2) in Table 1, where *<sup>β</sup>*<sup>∗</sup> *<sup>i</sup>* is the true value of *βi*.


**Table 1.** Estimation performance.

We try to predict the parameters close to true mean value of our target random variable and the variance is how scattered for our predictions. From Table 1, using the MLG prior for *β*, we got a reasonable estimation result because it achieves low bias and low variance simultaneously. Besides, we calculated the coverage probability for each variable, it indicates the 94% coverage probability for each parameter.

**Figure 1.** QQ-plot.

#### *4.3. Simulation for Model Selection*

In this simulation study, our goal was to study the accuracy of our model selection criteria. We have two different simulations in this section. Simulation 1: we set true *β* = (3, 0, 0) and calculated the difference between the true model and other candidate models for both criteria. In Figure 2, a difference beyond zero means that the true model had smaller DIC than the candidate model and the difference below zero means that the true model had higher LPML than the candidate model in Figure 2. The true model had the smallest DIC and the largest LPML in 99 of 100 simulated data sets. Simulation 2: we set true *β* = (1, 0, 0) and the results are shown in Figure 3. In each simulation, we have seven candidate models and one of them is true model and denote the true model as model 5. In Figures 2 and 3, the y-axis is the difference between "candidate model *i*" with true model. The true model had the smallest DIC in 81 of 100 simulated data sets and the largest LPML in 80 out of 100 simulated data sets. For each replicate dataset, we fit our model with 5000 Markov chain Monte Carlo iterations and treated the first 2000 iterations as burn-in. From Figures 2 and 3, in both simulation studies, we find that DIC and LPML yielded relatively consistent model selection results.

**Figure 2.** Deviance information criterion (DIC) and logarithm pseudo marginal likelihood (LPML) difference between candidate models and true model (model 5) of Simulation 1 ((**left**) DIC, (**right**) LPML).

**Figure 3.** DIC and LPML difference between candidate models and true model (model 5) of Simulation 2 ((**left**) DIC, (**right**) LPML).

#### *4.4. Simulation for Model Comparison*

In this simulation study, our goal is to evaluate the accuracy of our model selection criteria for different spatial random effects model. In this section, we generate the spatial random effects from MLG, *<sup>W</sup>* <sup>∼</sup> MLG(**0***n*, **<sup>Σ</sup>**1/2 *<sup>W</sup>* , *αW***1***n*, *κW***1***n*), where *α<sup>W</sup>* = *κ<sup>W</sup>* = 1. Other settings are same with previous simulations. We generated 100 data sets in these settings. Then, we compared the model fitness based on two following priors:

$$\begin{aligned} \mathbf{Prior} \, \mathbf{1}: \mathcal{W}|\phi\_{\prime}\sigma\_{w} &\sim \text{MLP}(\mathbf{0}\_{n\prime}\Sigma^{1/2}\_{\mathsf{W}}, a\_{\mathsf{W}}\mathbf{1}\_{n\prime}\kappa\_{\mathsf{W}}\mathbf{1}\_{n}),\\ \mathbf{Prior} \, \mathbf{2}: \mathcal{W}|\phi\_{\prime}\sigma\_{w} &\sim \mathcal{N}(\mathbf{0}\_{n\prime}\Sigma\_{\mathsf{W}}).\end{aligned} \tag{23}$$

For each replicate dataset, we fit our model with 5000 Markov chain Monte Carlo iterations and treat first 2000 iterations as burn-in. Then, we calculated the difference of DICs and the difference of LPMLs between these two priors. In Figure 4, the values below zero in the left plot imply that prior 1 has smaller DIC than prior 2. Also, the values above zero in the right plot in Figure 4 indicate that prior 1 has higher LPML than prior 2. The results shown in Figure 4 that we have a better result when we use the MLG prior than the Gaussian prior.

**Figure 4.** DIC and LPML difference ((**left**) DIC, (**right**) LPML).

#### **5. A Real Data Example**

#### *5.1. Data Description*

We analyzed seven days of US earthquake data collected in 2018, which includes *n* = 228 earthquakes that have magnitudes over *Zm* = 2.44 (https://earthquake.usgs.gov/). We present the earthquake data in Figure 5. We find the data most lie in seismic belts. In Figures 6 and 7, we present the histogram of this data and the scatter plot of this data set. In this analysis we have three variables (depth, gap, rms). The depth is where the earthquake begins to rupture. The gap is the largest azimuthal gap between azimuthally adjacent stations (in degrees). RMS is the root-mean-square (RMS) travel time residual, in sec, using all weights.

**Figure 5.** Map of US earthquake data.

**Figure 6.** Histogram of US earthquake data collected in 2018.

**Figure 7.** Scatter plot of earthquakes magnitudes, depth, gap and root-mean-square (RMS).

#### *5.2. Analysis*

We consider the model in Equation (7) and specify *αβ* = 105 and *κβ* = 10<sup>−</sup>5. These choices lead to an MLG that approximates a multivariate normal distribution. This choice of hyper-parameters will give an approximately normal prior on *β*. Inverse gamma priors are chosen for variance parameters *σ*<sup>2</sup> *w* and *σ*2, which is a usual choice of the variance parameters in Bayesian analysis. The full conditionals in the Appendix A are used to run a Gibbs sampler. We have seven candidate models in total, and *β* = (*β*1, *β*2, *β*3) =(depth, gap, rms). The number of iterations of the Gibbs sampler is 15,000, and the number of burn-in iterations is 10,000. The trace plots of posterior samples are provided in the Appendix B to show the convergence of MCMC chain. We also compare to a model when *W* approximates to Normal. The "DIC*N*" and "LPML*N*" denote the DIC and LPML for a model when *W* approximates to normal respectively. Furthermore, we calculated the log probability density (LPD) for candidate models. Based on the results in Table 2, the three criteria selected the same model with *β*<sup>1</sup> and MLG spatial random effects. Our proposed criteria had consistent results with the LPD.

From Table 2, we know that the model with *β* = (*β*1, 0, 0) has the smallest DIC and largest LPML. We also report the posterior estimates under the best model in Table 3 according to both DIC and LPML.


**Table 2.** Deviance information criterion (DIC) logarithm pseudo marginal likelihood (LPML), and log probability density (LPD) of candidate models.

**Table 3.** Posterior estimation under the best model.


From these posterior estimates, the model we select just contains depth as the important covariates and 95% credible interval does not contain zero. We see that as the depth increases, the expected value of earthquakes magnitudes increases. The other two covariates, gap and RMS, have no significant effects on earthquake magnitudes. In other words, from these seven-day earthquake data, deep earthquakes will have bigger magnitudes than shallow earthquakes. From the posterior estimates of *φ* and *σ*<sup>2</sup> *<sup>w</sup>*, we can find that there exists spatial correlation of earthquake magnitudes between different locations. In addition, using MLG as spatial random effects increases the goodness of fit of regression model in this data. This result is consistent with the earthquake literature [2].

#### **6. Discussion**

In this paper, we propose a Bayesian variable selection criterion for a Bayesian spatial-temporal model for analyzing earthquake magnitudes. Our main methodological contributions are to use the multivariate log-gamma model for both the regression coefficients and spatial random effects and to do variable selection for regression covariates with spatial random effects. Both DIC and LPML have a good selection power to choose the true model. But Bayesian model assessment criteria such as DIC and LPML do not perform well in the high-dimensional case, because the number of candidate models is very large when the number of covariates increases a lot. Developing a high-dimensional variable selection procedure is one of the important future works. The other future work is to fit other earthquake magnitudes models such as the gamma model or the Weibull model. In addition, we need to propose some Bayesian model assessment criterion to select the true data model for earthquake magnitudes. For the nature hazards problem, we need to incorporate the temporal dependent structure of earthquakes. Recently, the ETAS model [24] (combining the Gutenberg–Richter law and the Omori law) has been widely studied. Modelling earthquake dynamics is an important approach for preventing economic loss caused by an earthquake. Incorporating self-exciting effects in our generalized linear model with spatial random effects is another important future work. Furthermore, we only consider earthquake information as the covariates in our model. It will increase the predictive accuracy for us to combine more geographical information such as fault line information or crustal movement in the future.

**Author Contributions:** Data curation, G.H.; Formal analysis, H.-C.Y.; Investigation, H.-C.Y.; Methodology, G.H.; Project administration, G.H.; Software, H.-C.Y.; Supervision, M.-H.C.; Visualization, H.-C.Y.; Writing—original draft, G.H.; Writing—review & editing, M.-H.C.

**Funding:** Chen's research was partially supported by NIH grants #GM70335 and #P01CA142538. Hu's research was supported by Dean's office of College of Liberal Arts and Sciences in University of Connecticut.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Full Conditionals Distributions for Pareto Data with Latent Multivariate Log-Gamma Process Models**

From the hierarchical model in Equation (7), the full conditional distribution for *β* satisfies:

$$\begin{split} &f(\boldsymbol{\varPhi}) \propto f(\boldsymbol{\varPhi}) \prod f(\boldsymbol{Z}|\cdot) \\ & \propto \exp\left[\sum\_{i} (\mathbf{X}(\mathbf{s}\_{i})^{\prime}\boldsymbol{\varvarbeta} + \mathcal{W}(\mathbf{s}\_{i})) - \sum\_{i} (\log(\boldsymbol{Z}(\mathbf{s}\_{i}) - \log(\boldsymbol{Z}\_{\text{W}})) \exp(\mathbf{X}(\mathbf{s}\_{i})^{\prime}\boldsymbol{\upvarbeta} + \mathcal{W}(\mathbf{s}\_{i})) \right] \\ & \times \exp\left\{a\_{\beta}\mathbf{1}\_{p}^{\prime}\boldsymbol{\Sigma}\_{\text{fl}}^{-1/2}\boldsymbol{\upbeta} - \kappa\_{\beta}\mathbf{1}\_{p}^{\prime}\exp(\boldsymbol{\Sigma}\_{\text{fl}}^{-1/2}\boldsymbol{\upbeta})\right\}. \end{split} \tag{A1}$$

Rearranging terms we have

$$f(\boldsymbol{\beta}|\cdot) \propto \exp\left\{\mathfrak{a}\_{\beta}^{\prime} H\_{\beta} \boldsymbol{\mathfrak{f}} - \mathfrak{k}\_{\beta}^{\prime} \exp(H\_{\beta} \boldsymbol{\mathfrak{f}})\right\},\tag{A2}$$

which implies that *f*(*β*|·) is equal to cMLG(*Hβ*, *αβ*, *κβ*), which is a shorthand for the conditional MLG distribution used in [20].

Similarly, the full conditional distribution for *W* satisfies:

$$\begin{split} & f(\mathsf{W}|\cdot) \propto f(\mathsf{W}) \prod f(\mathsf{Z}|\cdot) \\ & \propto \exp\left[\sum\_{i} (\mathbf{X}(s\_{i})^{\prime}\boldsymbol{\mathcal{B}} + \boldsymbol{\mathcal{W}}(s\_{i})) - \sum\_{i} (\log(\boldsymbol{Z}(s\_{i}) - \log(\boldsymbol{Z}\_{m})) \exp(\mathbf{X}(s\_{i})^{\prime}\boldsymbol{\mathcal{B}} + \boldsymbol{\mathcal{W}}(s\_{i})) \right] \\ & \times \exp\left\{a\_{\mathsf{W}} \mathbf{1}\_{n}^{\prime} \boldsymbol{\Sigma}\_{\mathsf{W}}^{-1/2} \boldsymbol{\mathcal{W}} - \kappa\_{\mathsf{W}} \mathbf{1}\_{n}^{\prime} \exp(\boldsymbol{\Sigma}\_{\mathsf{W}}^{-1/2} \boldsymbol{\mathcal{W}}) \right\}. \end{split} \tag{A3}$$

Rearranging terms we have

$$f(\mathcal{W}|\cdot) \propto \exp\left\{\mathfrak{a}\_{\mathcal{W}}^{\prime} H\_{\mathcal{W}} \mathcal{W} - \mathfrak{a}\_{\mathcal{W}}^{\prime} \exp(H\_{\mathcal{W}} \mathcal{W})\right\},\tag{A4}$$

which implies that *f*(*W*|·) is equal to cMLG(*HW*, *αW*, *κW*). Thus we obtain the following full-conditional distributions to be used within a Gibbs sampler:

$$\begin{aligned} &\mathcal{B} \sim \text{cMLG}(H\_{\beta}, \mathfrak{a}\_{\beta}, \mathfrak{k}\_{\beta}) \\ &W \sim \text{cMLG}(H\_{W}, \mathfrak{a}\_{W}, \mathfrak{k}\_{W}) \\ &\sigma^{2} \propto \text{MLG}(\mathbf{0}, \mathfrak{L}\_{\beta}^{1/2}, \mathfrak{a}\_{\beta} \mathbf{1}\_{P}, \kappa\_{\beta} \mathbf{1}\_{P}) \times \text{IG}(a\_{1}, b\_{1}) \\ &\sigma\_{w}^{2} \propto \text{MLG}(\mathbf{0}, \mathfrak{L}\_{W}^{1/2}, \mathfrak{a}\_{W} \mathbf{1}\_{n}, \kappa\_{W} \mathbf{1}\_{n}) \times \text{IG}(a\_{2}, b\_{2}) \\ &\phi \propto \text{MLG}(\mathbf{0}\_{n}, \mathfrak{L}\_{W}^{1/2}, \mathfrak{a}\_{w} \mathbf{1}\_{n}, \kappa\_{w} \mathbf{1}\_{n}) \times \text{IG}(a\_{3}, b\_{3}), \end{aligned} \tag{A5}$$

where "cMLG" is the conditional multivariate log gamma distribution from [20]. A motivating feature of this conjugate structure is that it is relatively straightforward to simulate from a cMLG. For *σ*2, *σ*<sup>2</sup> *w* and *φ*, we consider using a Metropolis–Hasting algorithm or slice sampling procedure [25].

The parameters of the conditional multivariate log gamma distribution are organized into in Table A1.


**Table A1.** Parameters of the full conditional distribution.

#### **Appendix B. Trace Plot in Real Data Analysis**

**Figure A1.** (**upper left**) Trace plot for *β*; (**upper right**) Trace plot for *φ*; (**lower left**) Trace plot for *σ*2; (**lower right**) Trace plot for *σ*<sup>2</sup> *w*.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Undrained Cyclic Laboratory Behavior of Sandy Soils**

#### **Francesco Castelli 1, Antonio Cavallaro 2, Salvatore Grasso 3,\* and Valentina Lentini <sup>1</sup>**


Received: 1 August 2019; Accepted: 2 December 2019; Published: 11 December 2019

**Abstract:** The complex cyclic shear stress path experienced by the soil during an earthquake, which could also induce liquefaction phenomena, can be approximated in the laboratory only by using sophisticated testing apparatuses. Cyclic triaxial tests have been widely used, especially for coarse grained soils, as in this study. In the framework of the design for the seismic retrofitting of the "Ritiro viaduct" foundations along the A20 motorway connecting Messina with Palermo (Italy), a soil liquefaction study was also carried out. With this aim, a detailed geological and geotechnical characterization of the area was performed by in situ and laboratory tests, including seismic dilatometer Marchetti tests (SDMTs), the combined resonant column (RCT) and cyclic loading torsional shear tests (CLTSTs), and undrained cyclic loading triaxial tests (CLTxTs). In particular, the paper presents the results of cyclic triaxial tests carried out on isotropically consolidated specimens of a sandy soil. The seismic retrofitting works include the reinforcement of the foundation and replacement of the decks with newly designed type and structural schemes, mixed steel, and concrete with continuous girder. During the investigation, data were acquired for the characterization of materials, for the definition of degradation phenomena with the relative identification of possible causes, and for the estimation of the residual performance characteristics of the building. The structural campaign of investigations necessary to determine all of the key parameters useful for a correct definition of the residual performance capabilities of the work was divided into two phases: One in situ and one in the laboratory.

**Keywords:** in situ tests; laboratory tests; soil liquefaction; cyclic triaxial tests

#### **1. Introduction**

The present program of in situ investigations and laboratory tests originates from the static and seismic retrofitting works of the "Ritiro viaduct". The "Ritiro viaduct" (Figure 1), of the A20 Messina-Palermo (Italy) motorway, represents a vital node for the viability of Messina, the main connection to the motorway junctions. The works include the reinforcement of the foundations and replacement of the decks with newly designed type and structural schemes, mixed steel, and concrete with continuous girder.

During the investigation program, data were acquired for the characterization of materials, for the definition of degradation phenomena with the relative identification of possible causes, and for the estimation of the residual performance characteristics of the building. The structural campaign of investigations necessary to determine all of the parameters useful for a correct definition of the residual performance capabilities of the work has been useful to design the retrofitting works.

The campaign phase, in order to obtain as much data as possible without being extremely invasive towards the structural elements, was articulated through the execution of semi-destructive tests, such as the removal of concrete carrots, and non-destructive tests by means of execution of sclerometric tests, of pull-out with post-inserted grafts and in situ micro-seismic reliefs.

The area of the "Ritiro viaduct" was influenced by liquefaction phenomena during the 28 December 1908 Messina and Reggio Calabria historical earthquake. The evaluation of the liquefaction potential during earthquakes is an important subject in seismically active regions.

After the Niigata earthquake in 1964, which caused a lot of damage due to liquefaction, several studies [1–9] were performed in order to understand the cyclic behavior of sands.

During liquefaction, granular cohesionless saturated soil (gravel, sand, and low plasticity silt) loses its strength for a short interval of time, but long enough to cause significant failures. The liquefaction effects are usually evident on the ground surface (sand boils, large deformation, or fracture of the ground, etc.). The liquefaction susceptibility of soil deposits may be estimated by comparison between resistance profiles (i.e., Standard Penetration Test (SPT) blow count or Cone Penetration Test (CPT) tip resistance with depth) and critical values or estimation of a liquefaction safety factor (as a ratio of liquefaction resistance and action) function of depth.

The liquefaction potential can be obtained using either an estimation based on the maximum acceleration at the ground surface by a semi-empirical equation, or a dynamic calculation, including the reduction of soil stiffness due to built-up of pore pressure. Alternatively, the liquefaction resistance can be assessed by cyclic undrained laboratory tests on undisturbed or reconstituted specimens [10–12] or from correlations with the resistance measured by in situ tests (i.e., SPT or CPT).

In order to study the possible amplification phenomena of the "Ritiro viaduct" site, a comprehensive laboratory and in situ investigation was carried out to obtain a soil profile, with special attention being paid to the variation of the shear modulus (G) and damping ratio (D) with depth. This paper tries to summarize this information in a comprehensive way in order to provide a case record of site characterization for seismic response analysis.

**Figure 1.** The "Ritiro Viaduct" in Messina (Italy).

#### **2. Geology and Seismicity of the Area**

In Italy, recent strong earthquakes include the earthquake in central Italy that occurred on 24 August 2016 [13], the Emilia Romagna earthquake that occurred on 29 May 2012 [14–16], the L'Aquila earthquake that occurred on 6 April 2009 [17,18], the San Giuliano di Puglia earthquake that occurred on 31 October 2002 [19], the Umbria–Marche earthquake sequence of September–October 1997 [20], and the Catania earthquake that occurred on 13 December 1990 in southeastern Sicily.

The area under study is placed in the northeastern part of Sicily (Italy), at about 5 km from Messina old town. An example of its geological features is shown in Figure 2.

The area is covered by metamorphic complex in facies, with pan-African relicts in facies of granulitic deposits. The Messina Strait is located in the middle of the Calabrian Arc, one of the most seismically active areas of the Italian region and of the entire Mediterranean basin. The Strait's area itself has been struck several times in the past, though not with the same violence as in 1908.

Among Italian regions, Sicily is one of the most seismically active areas. In the past, strong earthquakes occurred in southeastern Sicily (1169, 1693) and in northeastern Sicily (1908). The MW 7.1, 28 December 1908, Messina Straits earthquake (Figure 3) was the deadliest earthquake in recent European history, and also one of the first to be investigated with modern instrumental data [21]. The shaking was distinctly felt in Albania, Montenegro, and the Greek Ionian islands, about 400 km to the east and northeast of the Strait; and in Malta, about 250 km to the south. The earthquake was catastrophic in the epicentral area and was immediately followed by fires and by a large tsunami. Messina and Reggio Calabria were almost completely destroyed, buildings were severely damaged (Figure 3) over an area in excess of 6000 km2 [21], and liquefaction phenomena occurred in the area.

According to [22] and [23] the main events reported in the historical catalogues and measured by Mercalli- Cancani- Sieberg scale (MCS) start around 91 BC (Reggio Calabria, IX–X MCS), 17 AD (Reggio Calabria, VIII–IX MCS), 361–363 (Messina Strait, X MCS), 853 (Messina, IX–X MCS), 1172 (Messina, VIII MCS), 1494 (Messina, VII–VIII MCS), 1509 sequence (area of Reggio Calabria, VIII MCS), 1659 (Southern Calabria, X MCS), 1783 sequence (Southwestern Calabria, X–XI MCS), 1894 (Southern Calabria, IX MCS), 1905 (Southwestern Calabria, XI MCS), and 1907 (Southern Calabria, VIII–IX MCS).

As well known, the characteristics of soil shaking are strongly influenced by local geological, geomorphological, and geotechnical conditions. These can modify significantly the amplitude, frequency, and duration of the seismic motion corresponding to the seismic hazard with reference to outcropping bedrock with horizontal topographical surface.

**Figure 2.** Geological synthesis of Messina area [24].

**Figure 3.** Messina Strait region with the 1908 earthquake epicenter and shocked localities with different colors according to Mercalli-Cancani-Sieberb scale [25].

To this aim, seismic microzonation defines the local seismic hazard through the identification of zones characterized by homogeneous seismic behavior that also includes the earthquake-induced effects, such as the slope instability, liquefaction in saturated granular soil, etc.

According to the Italian Guidelines for the Seismic Microzonation of the Department of Civil Defence, it is possible to distinguish three susceptibility classes: (1) Stable zones where no local effects occur (outcropping bedrock with a low steep morphological surface–slope angle\15◦); (2) stable zones but susceptible of local seismic amplification, where seismic input motion amplifications are expected, as a result of the local litho-stratigraphic and morphological structure; (3) unstable zones, where the expected seismic effects can be attributed to irreversible soil deformations (slope instability, liquefaction, etc.).

A zoom of the Seismic Microzonation map of the area, and particularly of the "Ritiro Viaduct", is shown in Figure 4. In the map, the zones identified as homogeneous are characterized by similar parameters as lithological and litho-technical characteristics, depth of bedrock, geomorphological conditions, etc. The area of the "Ritiro Viaduct" is within the local amplification stable zones (Zone 0, Zone 1, and Zone 5) of the Messina microzonation map.

Therefore the seismicity that affected the area, but also the effects induced by earthquakes [26,27], constitute a hazard that should not be ignored in the design or in the seismic retrofitting works of infrastructure.

**Figure 4.** Microzonation map of Messina in the area of the "Ritiro viaduct".

#### **3. Site Characterization Program and Basic Geotechnical Soil Properties**

The site investigation was performed within the area of the "Ritiro viaduct" and reached a maximum depth of 35.00 m. Laboratory tests were performed on nine disturbed samples retrieved by means of a 101 mm tube sampler. To evaluate the geotechnical characteristics, the following in situ and laboratory tests were performed in the foundation soil located in the area of the "Ritiro viaduct": 28 Boreholes, 7 seismic dilatometer Marchetti tests (SDMTs) [28–31]; 9 particle size analysis, 25 direct shear tests (DST), 1 undrained triaxial test (UTxT), 8 cyclic loading torsional shear tests (CLTSTs) [32,33], 8 resonant column tests (RCTs) [34–36], and 3 cyclic loading triaxial tests (CLTxTs). The investigation program follows the approach and the methodology used in other test sites in the seismic areas of Catania and Messina [37–41]. On the basis of laboratory tests, the "Ritiro viaduct" deposits mainly consist of a normal consolidated grey/dark grey or dark brown sand with silty gravel or gravelly silt.

The results obtained by particle size analysis are shown in Figure 5 and Table 1. In Table 1, the soil samples were divided into Type A and Type B according to the depth and the presence of minor or major gravel. The uniformity coefficient (Cu), defined as the ratio of D60 to D10, varies in the range of 5.48–52.24 and it points out the considerable non-homogeneity of the particle size. The maximum dry density was evaluated using a vibrating table, available at the geotechnical laboratory of the University Kore of Enna, according to ASTM (4253-83) standard.

Physical parameters were derived from standard classification tests performed on the samples retrieved by geotechnical survey. The values of soil unit weight and the minimum (emin) and maximum (emax) void ratios are summarized in Table 2. Most of the samples are coarse-grained soils, classifiable as silty sands to gravelly sands, showing a lower percentage of clayey material.

**Figure 5.** Profile of grain size distribution of the "Ritiro viaduct" area.



**Table 1.** *Cont*.

Type of soil A = soil with minor gravel; Type of soil B = soil with major gravel; H = depth; dmax = maximum diameter; D60 = the particle-size diameter for which 60% of the sample was finer; D10 = the particle-size diameter for which 10% of the sample was finer; Cu = uniformity coefficient.

**Table 2.** Mechanical characteristics for the "Ritiro viaduct" area.


Type of soil A = soil with minor gravel; Type of soil B = soil with major gravel; H = depth; γ = total unit weight; D60 = the particle-size diameter for which 60% of the sample was finer; Cu = uniformity coefficient; Gs = specific gravity; eo = initial void ratio; emin = void ratio of soil in densest condition; emax = void ratio of soil in loosest condition; Dr = relative density; CF = clay fraction.

The value of the natural moisture content wn prevalently ranges from between 22 to 35% (Figure 6). Characteristic values of strength parameters are (cohesion) c' = 5–24 kPa and (angle of shear resistance)' = 23–40◦, obtained by direct shear test, (undrained cohesion) cu = 228 kPa obtained by a undrained triaxial test performed on a cohesive sample, with a water content of 17%, retrieved at the depth of 14.65 m; Gs (specific gravity) ranged between 2.65 and 2.75, while eo ranged between 0.29 and 0.67 (Table 2). Another stratum with a water content of 12% was found at the depth of about 19.25 m. Figure 6 shows index properties of the "Ritiro viaduct" area.

**Figure 6.** Particle size analysis of the "Ritiro viaduct" area.

#### **4. Soil Properties by Laboratory Tests**

The equivalent shear modulus Geq and damping ratio D of "Ritiro viaduct" deposits were determined in the laboratory by means of a Resonant Column/Cyclic Loading Torsional Shear apparatus [42–44]. These dynamic parameters represent the basic data for the studies of local seismic response [45–48]. This apparatus was supplied at the Geotechnical Laboratory of the University Kore of Enna (Figure 7). A resonant column test consists of exciting one end of a confined solid or hollow cylindrical soil specimen. The specimen is fixed at the bottom (fixed-free test) and it is excited in torsion or flexure at the top by means of an electromagnetic drive system. Once the fundamental resonant frequency is established from measuring the motion of the free end, the velocity of the propagating wave and the degree of material damping are derived. The shear modulus is then obtained from the derived velocity Vs (in case of torsion) and the density of the sample.

**Figure 7.** Resonant column/cyclic torsional shear apparatus.

The equivalent shear modulus Geq is the unload–reload shear modulus that is evaluated from RCT in function of velocity Vs and density ρ of the sample, while Go is the maximum value or also "plateau" value as observed in the G-log(γ) plot. G is the secant modulus. Generally, G is constant until a certain strain limit is exceeded. This limit is called the elastic threshold shear strain (γ<sup>t</sup> e) and it is believed that soils behave elastically at strains smaller than (γ<sup>t</sup> e). The elastic stiffness at γ < γ<sup>t</sup> <sup>e</sup> is thus already defined as Go. Damping ratio D is defined by Equation (1):

$$D = \frac{\Delta W}{4\pi W'} \tag{1}$$

where ΔW is the area enclosed by the unloading–reloading loop and represents the total energy loss during the cycle, where W is the elastic stored energy. The RCT apparatus used is a fixed-free resonant column apparatus [49]. It enables the specimen consolidation under both isotropic and anisotropic stresses. It is composed of a drive system, a support system, and a base plate. The solid or hallow cylinder specimen is fixed at the bottom and its constraint at the base is due to the friction existing between the specimen and the porous synthesized bronze stone [50].

Torsional forces are applied at the top of the drive system, realized in aluminum. It is an electrical motor constituted of four magnets connected with the top of the sample and eight coils placed on the inox steel annular base, which is strictly linked to the support system. The weight of the motor is counterbalanced by a spring. A programmable function generator (PGF) excites the electrical motor of [51].

The support system, in addition to permitting the placement of the drive system, may possibly put the proximity transducers in and the filling in of water for saturated specimen tests. It is realized a Plexiglas cell pressure, to permit the isotropic consolidation using an air pressure source controlled with a manual pressure regulator. The base and the top plates are connected by three vertical rods inside the cell. In the resonant column test (RCT), the function generator produces a sinusoidal electric signal which is increased from an amplifier and transformed into torsional mechanical stress by an electromagnetic motor, consisting of eight coils and four magnets, connected to the head of the specimen. The magnetic field of the coils interacts with the magnets connected to the plate that transmits the torsional oscillation to the head of the specimen. Because the frequency of excitation is varied, the dynamic response of the sample varies in terms of amplitude. The latter is recorded by a accelerometer connected to the guide plate and proximity transducers measuring the relative movement between the plate and coils. Resonance frequency values were used for the calculation of the cutting module and the deformations. The decay curve, produced by interrupting the torsion excitation in resonance, allowed to evaluate the damping of the material following the amplitude decay method during the decrement of free vibration [42].

In the cyclic loading torsional shear test (CLTST), a cyclic torque is applied to the specimen by means of a torsional motor, to which a sinusoidal signal is sent at a fixed frequency of 0.5 Hz, continuously monitoring the torsion and angular deformation. The rotation of the specimen is measured thanks to the use of two proximity transducers. The data of the signal sent (proportional to the shear stress) and the corresponding torsion (proportional to the cutting deformation) are acquired simultaneously. The specimen shear module is then determined based on the average slope of the stress–strain hysteresis, while the damping is related to the cycle area of hysteresis as the ratio between the area enclosed by the unloading–reloading loop, and represents the total energy loss during the cycle and W is the elastic stored energy.

In the present work, solid cylindrical specimens were reconstituted by using tapping [50], in order to obtain the required relative and a good uniformity during the deposition.

The mold was assembled and a little depression was applied to let the membrane adhere to the inside surfaces. The material was placed in the mold using a funnel-pouring device. The soil was placed as loosely as possible in the mold by leaving the soil from the spout in a steady stream, holding the pouring device upright and vertical, and maintaining constant the fall height. It was possible to obtain different values of relative density changing the height of deposition. In order to realize high values of relative density, it could be necessary to beat delicately the mold surface during the deposition. Each sample was reconstituted with fresh sand. Each specimen was subjected to an isotropic load achieved in a Plexiglas pressure cell, using an air pressure source. The axial strain was measured by using a high-resolution proximity transducer, which monitors the aluminum top-cap displacement. Shear strain was measured by monitoring the top rotation with a couple of high-resolution proximity transducers. During a resonant column test, the proximity transducers were not able to appraise the value of the targets displacements, because of the high frequency of the oscillations. Then rotation on the top of the specimen was measured by means of an accelerometer.

The laboratory test conditions and the obtained small strain shear modulus Go are listed in Table 3. After the saturation phase, obtained by applying an appropriate back-pressure value, the undisturbed specimens were isotropically reconsolidated to the best estimate of the in situ mean effective stress. The same specimen was first subject to RCT, then to CLTST after a rest period of 24 h with opened drainage. CLTST was performed under stress control condition by applying a torque, with triangular time history, at a frequency of 0.1 Hz. The size of solid cylindrical specimens were radius = 25 mm and height = 100 mm.


**Table 3.** Test condition for the "Ritiro viaduct" area.

Type of soil A = soil with minor gravel; Type of soil B = soil with major gravel; H = depth; U = undrained; RCT = resonant column test; CLTST = cyclic loading torsional shear tests; CLTxT = cyclic loading triaxial test; Go (1) from RCT, Go (2) from CLTST after 24 h, Go (3) from seismic dilatometer Marchetti test (SDMT).

The initial shear modulus Go of soil is mainly influenced by the state of the soil, expressed by a combination of the void index e (or by the relative density Dr) and by the soil structure that reflects the deposition and the subsequent structural transformation processes such as aging, diagenesis, and cementation. For higher strain levels, the shear modulus G depends also, and especially on the strain level, on the stress–strain history and on the strain rate [52].

The Go values (Go (1) and Go (2)), reported in Table 3, indicate moderate influence of strain rate, even at very small strain where the soil behavior is supposed to be elastic. Values of shear modulus G (MPa) and damping ratio D (%) versus from RCT and CLTST tests are reported in Figures 8–11, respectively. The initial shear modulus obtained during CLTST shows the effect of the soil degradation because of RCT; this effect tends to become negligible with the shear strain build up (Figure 12). Meanwhile, the D values obtained during RCT and CLTST follow the same trend and are thus comparable (Figure 13). The damping ratio values obtained from RCT by amplitude decay and CLTST method are quasi constant until a strain level of about 0.01%, higher values of D have been obtained from strain level higher than 0.01%. It is possible to see that the damping ratio from RCT and CLTST, at very small strains, is so equal to about 2%. Greater values of D are obtained from RCT for the strain level of about 0.1%.

It is supposed that RCT provides larger values of D than CLTST because of the rate (frequency) effect, in agreement with data shown by [53,54]. According to these researchers, the nature of soil damping in soils can be linked to the following phenomena:


Finally, higher values of the initial shear modulus (Go (3)) were obtained from SDMTs. Generally, the small strain stiffness, determined in the laboratory on high quality reconstructed sample using appropriate apparatuses and procedures, is very close to that obtained in situ from seismic tests. Probably, in the case of the "Ritiro viaduct", disturbance phenomena occurred during reconstruction operations and differences in stress conditions determined lower values of the initial shear modulo in the laboratory.

The experimental results were used to determine the empirical parameters of the equation proposed by [55] to describe the shear modulus decay with shear strain level (Figure 14a,b and Figure 15a,b):

$$\frac{\mathbf{G}(\boldsymbol{\gamma})}{\mathbf{G}\_{\mathbf{0}}} = \frac{1}{1 + \boldsymbol{\alpha}\boldsymbol{\gamma}(\%)^{\boldsymbol{\beta}}},\tag{2}$$

where G(γ) = strain dependent shear modulus; γ = shear strain; α, β = soil constants.

Equation (2) allows the complete shear modulus degradation to be considered with strain level. The values of soil constants α and β obtained from RCTs and CLTSTs for soil type A and B are listed in Table 4.

**Figure 8.** G–γ curves from RCTs.

**Figure 9.** D–γ curves from RCTs.

J **>@**

**Figure 12.** G–γ curves from RCTs and CLTSTs.

**Figure 15.** G/Go–γ curves from CLTSTs. (**a**) Type A; (**b**) Type B.


**Table 4.** Soil constants for the "Ritiro viaduct" area from RCTs and CLTSTs.

As suggested by [55], the inverse variation of damping ratio with respect to the normalized shear modulus has an exponential form as that reported in Figure 16a,b and Figure 17a,b for the "Ritiro viaduct" area:

$$\mathrm{D}(\gamma)(^{\circ}\!\!\!\!\!\/) = \mathfrak{n} \cdot \exp\left[-\lambda \cdot \frac{\mathrm{G}(\gamma)}{\mathrm{G}\_{0}}\right] \,,\tag{3}$$

where D(γ) = strain dependent damping ratio; γ = shear strain; η, λ = soil constants.

**Figure 16.** D–G/Go curves from RCTs. (**a**) Type A; (**b**) Type B.

**Figure 17.** D–G/Go curves from CLTSTs. (**a**) Type A; (**b**) Type B.

The values of soil constants η and λ obtained from RCTs and CLTSTs for soil type A and B are listed in Table 4.

Equation (3) assumes maximum value Dmax = 45% for G(γ)/Go = 0 and minimum value Dmin = 3.18% for G(γ)/Go = 1. Therefore, Equation (3) can be re-written in the following normalized form:

$$\frac{\mathbf{D}(\boldsymbol{\gamma})}{\mathbf{D}(\boldsymbol{\gamma})\_{\text{max}}} = \exp\left[-\boldsymbol{\lambda} \cdot \frac{\mathbf{G}(\boldsymbol{\gamma})}{\mathbf{G}\_{\text{o}}}\right] \tag{4}$$

Figure 18a,b and Figure 19a,b show a comparison between the enveloping curves of the experimental data obtained during the RCTs and CLTSTs. It is possible to observe how, in general, the CLTSTs, with respect to the RCTs, determine a behavior of the soil characterized by a wider elastic field with the same level of deformation reached. This phenomenon may be due to the intergranular reassembly of the sandy soil due to the effect of the RCTS. Moreover, overall, in the case of CLTSTs higher values of D are observed.

**Figure 18.** G/Go–γ curves from RCTS and CLTSTs by [55]. (**a**) Type A; (**b**) Type B.

**Figure 19.** D–G/Go curves from RCTS and CLTSTs by [55]. (**a**) Type A; (**b**) Type B.

To perform triaxial tests, triaxial cells consisting of a structure were used for stainless steel and a Plexiglas cylinder (Figure 20). The maximum isotropic operating pressure was 1 MPa. For fluid of confinement, water was used. Sample drainage was allowed, when necessary, through porous stones placed on the two vertical load distribution bases. The specimen was subsequently placed in saturation by applying a back-pressure under an effective pressure isotropic, enough to prevent swelling. The test was then saturated by performing one measurement of parameter B. The value of approximately 0.95 was taken as an indirect measure of the complete saturation of the material. In the case of a value that is too low, saturation was prolonged for a further period of time, in some cases increasing counter-pressure again until a satisfactory value of B was reached and the specimen was brought, in several steps, to the final effective consolidation tension.

All test operations during saturation (until the Skempton B parameter reached at least 0.98) and isotropic consolidation were controlled by a panel that adjusts the confinement pressure and the counter-pressure and allows measurement of pressures and pressure volume variations of the specimen by means of a pressure transducer and a volumometer. The height variations of the specimen were detected by means of a displacement transducer. The application of cyclic loads took place by means of a contrast structure equipped with an electro-pneumatic system, which allowed to apply to the specimen a constant sinusoidal load of constant amplitude. The size of solid cylindrical specimens were radius = 35 mm and height = 140 mm. The laboratory test conditions and the obtained small strain shear modulus Eo are also listed in Table 3. During the cyclic triaxial test, the load sequence was characterized by steps of 40 strain controlled load cycles.

During the cyclic loading triaxial tests (CLTxTs) (Figure 21), the soil sample showed a rapid decrease of its mechanical characteristics at strain levels of about 10<sup>−</sup>2%. It seems that it is not possible to investigate the values of the modulus of normal elasticity (Young's modulus) at very low strain levels (less than 10<sup>−</sup>3%) due to undesired deformations caused by the deformability of the mechanical structure (system compliance) of the triaxial apparatus.

**Figure 20.** Cyclic triaxial apparatus.

**Figure 21.** E/Eo–ε curves from CLTxTs.

Low values of the Young's modulus were obtained in correspondence with the test performed on the S17C7 sample, probably due to the low initial value of Dr.

The initial damping values are around 1%, while the maximum values are between 6 and 8% at a strain level of about 1% (Figure 22).

During the cyclic loading triaxial tests (CLTxT), unload–reload cycles became unstable and degradation phenomena of material occured when a certain limit strain was exceeded (Figures 23 and 24). This limit strain is defined as volumetric threshold shear strain and is rate-dependent. The degradation caused a decrease of stiffness, an increase of D and pore pressure build-up with the increase of N because of cyclic material degradation, as obtained from a CLTST and CLTxT on "Ritiro viaduct" soil [56].

Figure 25 shows the pore pressure build up during CLTxT. The pore pressure build up during CLTxT was so negligible at low strain. On the contrary, at strain level of about 0.15%, it is possible to observe an important increase of pore pressure due to the degradation phenomenon.

**Figure 25.** Δu–ε curves from CLTxTs.

Values of shear modulus G (MPa) and damping ratio D (%) versus γ (%) from CLTST and CLTxT are reported in Figures 26 and 27, respectively. In the case of CLTST, higher values of G were always obtained, compared to those of the CLTxT, also as a function of lower initial investigated strain levels. The trend of the G modulus seems to align only for strain levels higher than 0.1%, even if the results of the two types of tests are comparable only for a strain interval between 0.01 and 0.1%. This difference on G can probably be attributed to the high interstitial pressure values obtained during the CLTxT tests (Figure 25) due to the low initial value of Dr. Higher values of D were obtained from strain levels higher than 0.01%. Moreover, at the same strain level, in the case of CLTSTs, higher values of D were observed.

**Figure 26.** G–γ curves from CLTST and CLTxT.

**Figure 27.** D–γ curves from CLTST and CLTxT.

#### **5. Soil Properties by in Situ Tests**

 The use of in situ tests as complementary experimental techniques to laboratory experiments is now commonly recognized [57]. The small strain (γ ≤ 0.001%) shear modulus, Go, was thus determined from SDMT. The SDMT provides a simple means for determining the initial elastic stiffness at very small strains and in situ shear strength parameters at high strains in natural soil deposits. Moreover, it was attempted to assess Go by means of empirical correlations, based either on penetration test results or on laboratory test results [57]. The SDMT [17,58] provides a simple means for determining the initial elastic stiffness at very small strains and in situ shear strength parameters at high strains in natural soil deposits [32,59]. This apparatus was also used in offshore conditions by [60,61]. The test is conceptually similar to the seismic cone (SCPT). First introduced by [62], the SDMT was subsequently improved at Georgia Tech, Atlanta, USA [63–65]. A new SDMT system has recently been developed in Italy. The seismic modulus is a cylindrical instrumented tube, located above the DMT blade [66], housing two receivers at a distance of 0.50 m (see Figure 28). The test configuration "two receivers"/"true interval" avoids the problem connected with the possible inaccurate determination of the "first arrival" time sometimes met with the "pseudo interval" configuration (just one receiver). Moreover, the pair of seismograms recorded by the two receivers at a given test depth correspond to the same hammer blow and not to different blows in sequence, which are not necessarily identical. The adoption of the "true interval" configuration considerably enhances the repeatability in the Vs measurement (observed repeatability Vs ≈ 1–2%). VS is obtained as the ratio between the difference in distance between the source and the two receivers (S2–S1) and the delay of the arrival of the impulse from the first to the second receiver (Δt). Vs measurements are obtained every 0.5 m of depth. The shear wave source at the surface is a pendulum hammer (≈10 kg), which hits horizontally a steel rectangular base pressed vertically against the soil (by the weight of the truck) and oriented with its long axis parallel to the axis of the receivers, so that they can offer the highest sensitivity to the generated shear wave.

Source waves are generated by striking a horizontal plank at the surface that is oriented parallel to the axis of a geophone connects by a co-axial cable with an oscilloscope [63,64]. The measured arrival times at successive depths provide pseudo interval Vs profiles for horizontally polarized vertically propagating shear waves. In Figure 28, the SDMT scheme for the measure of Vs is shown, while Figure 29 shows an example of seismograms obtained by SDMT at various test depths at the site of the "Ritiro viaduct" (it is a good practice to plot side-by-side the seismograms as recorded and re-phased according to the calculated delay). Vs may be converted into the initial shear modulus Go by the theory of elasticity.

**Figure 28.** Seismic dilatometer equipment (**a**). Schematic layout of the flat dilatometer test (**b**) and of the seismic dilatometer test (**c**).

**Figure 29.** Example of seismograms obtained by SDMT at the site of the "Ritiro viaduct".

In three SDMT test verticals, only the seismic measurements were carried out, in pre-holes performed by means of a probe and filled with gravel (with grains of a diameter strictly between 5 and 15 mm), statically advancing with a penetrometer having a maximum thrust capacity equal to 20 tons. The noticeable difference between the density of the in situ material and the filling material of the pre-hole made the interpretation of the results particularly difficult. The combined knowledge of Go and of the one dimensional modulus M (from DMT) may be helpful in the construction of the G–γ modulus degradation curves [67–71].

A summary of SDMT parameters is shown in Figure 30, where:


**Figure 30.** Results of the SDMTs (SDMT 4a) in terms of geotechnical parameters.

Figure 31 shows the values of Go obtained in situ from SDMT and those measured in the laboratory from RCT performed on reconstructed solid cylindrical specimens, which were isotropically reconsolidated to the best estimate of the in situ mean effective stress. The Go values are plotted in Figure 31 against depth. In the case of laboratory tests, the Go values are determined at shear strain levels of less than 0.001%. A comparability exists between the laboratory and in situ test results. On average, the ratio of Go (Lab) to Go (Field) by RCT and SDMT was equal to about 1.80 at the depth of 19.75 m.

**Figure 31.** Go obtained from SDMT, RCT, and empirical correlations.

It was also attempted to evaluate the small strain shear modulus by means of the following empirical correlations based on penetration tests results or laboratory results available in literature.

$$\begin{array}{c} \text{(a)} \quad \text{ [Z]} \end{array}$$

$$\mathbf{G}\_{\mathbf{o}} = \frac{530}{\left(\sigma\_{\mathbf{v}}^{\prime}/\mathbf{p}\_{\mathbf{a}}\right)^{0.25}} \frac{\gamma\_{\mathbf{D}}/\gamma\_{\mathbf{w}} - 1}{2.7 - \gamma\_{\mathbf{D}}/\gamma\_{\mathbf{w}}} \mathbf{K}\_{\mathbf{o}}^{0.25} \cdot \left(\sigma\_{\mathbf{v}}^{\prime} \cdot \mathbf{p}\_{\mathbf{a}}\right)^{0.5},\tag{5}$$

where Go, σ'v and pa are expressed in the same unit; pa = 1 bar is a reference pressure; γ<sup>D</sup> and Ko are, respectively, the unit weight and the coefficient of earth pressure at rest, as inferred from SDMT results according to [66];

(b) [57]

$$\mathbf{G}\_{\mathbf{o}} = \frac{600 \cdot \sigma\_m^{0.5} \mathbf{P\_a^{0.5}}}{\mathbf{e^{1.3}}},\tag{6}$$

where σ'm = (σ'v + 2 · σ'h)/3; pa = 1 bar is a reference pressure; Go, σ'm, and pa are expressed in the same unit. The values for parameters which appear in Equation (6) are equal to the average values that result from laboratory tests performed on quaternary Italian clays and reconstituted sands. A similar equation was proposed by [73] for Holocene clay deposits.

Equation (6) incorporates a term which expresses the void ratio; the coefficient of earth pressure at rest only appears in Equation (5). However only Equation (5) tries to obtain all the input data from the SDMT results. The Go values obtained with the methods above are also plotted against depth in Figure 28. The method by [57] was applied considering a given profile of void ratio. The coefficient of earth pressure at rest was inferred from SDMT.

Since the purely dilatometric data only investigate the first meters of depth, the use of Equation (5) is limited. On the whole, Equation (6) seems to provide the most accurate trend of Go with depth, but is not able to analyze stratigraphic variation along the depth, as can be seen in Figure 28. The results obtained by SDMT are comparable with the data of the RCT tests and they are able to identify the stratigraphic variations.

#### **6. Analysis of the E**ff**ects on the Physical Environment**

The results of cyclic triaxial tests carried out can be used for the seismic retrofitting works, which include also the reinforcement of the foundation of buildings. Data include characterization of materials, definition of degradation phenomena with the relative identification of possible causes, and estimation of the residual performance characteristics of the building. To take into account the analysis of the effects on the natural environment, according to the *"Manual for Zonation on Seismic Geotechnical Hazards"*, local seismic response, slope instability, and liquefaction of the area were analyzed using the results of cyclic triaxial tests and the results of other laboratory tests, including the combined resonant column tests (RCT) and cyclic loading torsional shear tests (CLTST). Local seismic response analysis using ONDA and DEEPSOIL computer codes was performed on the Ritiro Viaduct [74]. Results of the numerical analyses are presented as comparisons in terms of maximum acceleration profiles, maximum shear strain profiles, response spectra, surface response seismograms, Fourier spectra, and amplification ratios. The local seismic response (LSR) analysis was performed by using seismograms obtained for the 28 December 1908 Messina and Reggio Calabria earthquake scenario. The results of cyclic triaxial tests performed on samples were also used for the evaluation of the liquefaction resistance of soils of the Ritiro Viaduct [12]. The complex cyclic shear stress path experienced by the soil during an earthquake can be reproduced in the laboratory only by using sophisticated testing apparatuses. Cyclic triaxial tests have been widely used to assess soil liquefaction potential, especially for coarse-grained soils, as in this study.

#### **7. Conclusions**

In the framework of the design for the seismic retrofitting of the "Ritiro Viaduct" foundations along the A20 motorway, connecting the cities of Messina and Palermo, located in one of the most hazardous Italian seismic areas, a detailed geotechnical characterization was carried out. Indeed, the seismic effects induced by earthquakes play an important role in the planning and construction, as well as in the seismic retrofitting works of important infrastructures, such as the "Ritiro viaduct".

This paper focuses on a comprehensive laboratory and in situ investigations carried out to obtain a soil profile, with special attention to the variation of the shear modulus (G) and damping ratio (D) with depth. A detailed geological and geotechnical characterization of the area was performed by in situ and laboratory tests, including seismic dilatometer Marchetti tests (SDMTs), the combined resonant column (RCT) and cyclic loading torsional shear tests (CLTSTs), and undrained cyclic loading triaxial tests (CLTxTs).

The RCT and CLTST show a moderate influence of strain rate, even at very small strain, where the soil behavior is supposed to be elastic, while the D values obtained during RCT and CLTST follow the same trend and are thus comparable. It is possible to see that the damping ratio from RCT and

CLTST, at very small strains, is equal to about 2%. Moreover, higher values of the initial shear modulus Go were obtained from SDMTs. Probably, in the case of the "Ritiro viaduct", disturbance phenomena occurred during reconstruction operations and differences in stress conditions determined lower values of the initial shear modulus in the laboratory.

The experimental results were used to design the seismic retrofitting work and, in particular, to determine the empirical parameters of the proposed equation to describe the shear modulus decay and damping ratio build-up with shear strain level.

During the cyclic loading triaxial tests (CLTxTs), the soil sample showed a rapid decrease of its mechanical characteristics at strain levels of about 10<sup>−</sup>2%. It seems that it is not possible to investigate the values of the modulus of normal elasticity (Young's modulus) at very low strain levels (less than 10<sup>−</sup>3%) due to undesired deformations caused by the deformability of the mechanical structure (system compliance) of the triaxial apparatus. Low values of the Young's modulus were obtained in correspondence with the test performed on one sample, probably due to the low initial value of Dr. The initial damping values were around 1%. During CLTxT, unload–reload cycles became unstable and degradation phenomena of material occurred when a certain limit strain was exceeded. This limit strain is called volumetric threshold shear strain and it is rate-dependent. The degradation caused a decrease of stiffness, an increase of D, and pore pressure build-up with the increase of N because of cyclic material degradation. Finally, the in situ obtained results by SDMT, though higher, are comparable with the data of the RCT tests and they are able to identify the stratigraphic variations.

**Author Contributions:** F.C., A.C., S.G. and V.L. carried out the investigation and prepared the original manuscript according to the following percentages: 25% F.C., 25% A.C., 25% S.G. and 25% V.L.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Machine Learning Methods for Seismic Hazards Forecast**

#### **Valeri G. Gitis and Alexander B. Derendyaev \***

The Institute for Information Transmission Problems, Moscow 127051, Russia **\*** Correspondence: wintsa@gmail.com

Received: 25 April 2019; Accepted: 10 July 2019; Published: 12 July 2019

**Abstract:** In this paper, we suggest two machine learning methods for seismic hazard forecast. The first method is used for spatial forecasting of maximum possible earthquake magnitudes (*Mmax*), whereas the second is used for spatio-temporal forecasting of strong earthquakes. The first method, the method of approximation of interval expert estimates, is based on a regression approach in which values of *Mmax* at the points of the training sample are estimated by experts. The method allows one to formalize the knowledge of experts, to find the dependence of *Mmax* on the properties of the geological environment, and to construct a map of the spatial forecast. The second method, the method of minimum area of alarm, uses retrospective data to identify the alarm area in which the epicenters of strong (target) earthquakes are expected at a certain time interval. This method is the basis of an automatic web-based platform that systematically forecasts target earthquakes. The results of testing the approach to earthquake prediction in the Mediterranean and Californian regions are presented. For the tests, well known parameters of earthquake catalogs were used. The method showed a satisfactory forecast quality.

**Keywords:** machine learning; expert estimate; maximum possible magnitudes of earthquakes; one class classification; seismic hazard; seismic zoning; earthquake forecasting

#### **1. Introduction**

Tectonic earthquakes are invariably preceded by a period when stresses increase in the Earth. This process forms anomalous changes in the geological environment near the source of the expected earthquake [1–3]. To describe the seismotectonic properties of the geological environment, various types of data are used: Earthquake catalogs, time series of geodetic [4], geophysical [5] and geochemical measurements [6], and aerospace observations [7]. The success of seismic hazard forecast is largely influenced by both completeness of the description of the spatial and spatio-temporal properties of the seismic process, and the possibility of their joint analysis. In our approach to joint analysis, all available data on the properties of the process are converted into grid fields [8].

Seismic zoning is prerequisite for seismic hazard assessment [9]. The most important and complex problem of seismic zoning is to map the maximum possible magnitudes of earthquakes (*Mmax*). The values of *Mmax* cannot be measured instrumentally. Two assumptions are used to construct a digital map of *Mmax*: (1) The assumption of large earthquake repetition [10] and (2) the assumption that the values of *Mmax* depend on the properties of the geological environment [11,12].

The statistical approach uses only the first assumption. This means that the *Mmax* map is calculated using only those earthquakes whose epicenters fall into a sliding spatial window. The methods of extreme statistics are used for the estimation of *Mmax* [13–16]. These methods require a sufficiently large number of observations, which may be unavailable for some zones in the region. To improve this method, [17] used the second assumption, and estimated *Mmax* by earthquake epicenters within geologically homogeneous zones identified by an expert geologist. As in the previous approach, this method does not provide for extrapolating *Mmax* values from zones in which there are many epicenters of strong earthquakes and therefore estimates of *Mmax* are fairly accurate, to zones with similar seismotectonic properties, but with a small number of strong earthquakes.

The history of seismic observations is very short in relation to the speed of tectonic processes, and earthquakes with magnitudes close to maximum occur relatively rarely. To compensate for this effect, attempts are made to extrapolate reliable estimates to areas with similar seismotectonic properties of the geological environment. Reference [18] used the *Mmax* mapping method based on a solution of a group of experts. One of the possible algorithmic approaches applies cluster analysis [19]. A cluster analysis program divides the region into zones with similar values of geological and geophysical characteristics. These zones may consist of several isolated areas. Next, the maximum magnitudes of earthquakes recorded in one area of the zone are extrapolated to all other areas. The disadvantages of this approach are related to the fact that the zoning of a region into quasi-homogeneous zones is largely determined a set of selected features, the method of measuring the similarity between clusters, the type of clustering algorithm and, finally, the criterion of stopping the clustering process.

We describe the method of approximation of interval expert estimates which is a regression approach to the construction of a forecast map of *Mmax* [20–22]. To compile the map, dependence of *Mmax* on properties of geological environments **x** = (*x*1, ... , *xI*) is used. The values at the training sample points are determined using expert knowledge. To this end, experts choose the most studied points of the region with different seismicity and geological conditions. The expert indicates the boundaries of the interval in which, in his opinion, lies the value of *Mmax*, and evaluates the values of the confidence that *Mmax* cannot exceed the lower and upper limits of the interval. In the assessment, the expert uses historical seismic data, instrumental data on the maximum magnitude of an earthquake in the vicinity of the point in question and data on the properties of the geological zone to which this point belongs. The algorithm generalizes the least squares approximation algorithm.

The task of predicting an earthquake is to determine the time, location, and magnitude of a future earthquake. Earthquake prediction studies are conducted in many directions. They include the study of the rock failure and earthquake precursor phenomena, the study of stochastic models for earthquake prediction, machine learning methods, and testing earthquake prediction algorithms [1–3,23–29]. At the same time, there are a number of works in which it is stated that earthquakes cannot be predicted [30].

Here we suggest a new method of machine learning, called the method of the minimum area of alarm, and describe a web-based platform that predicts earthquakes in automatic mode (http://distcomp.ru/geo/prognosis/). Our method solves the one-class classification problem (other methods can be found, for instance, in [31–33]). Our training sample set includes rare anomalous objects (the epicenters of target earthquakes) and grid fields of properties of the seismic process (field of features). The method allows one to detect the largest number of the target earthquakes for the training set, provided that the size of the spatio-temporal alarm area does not exceed a specified value. We present the results of testing the approach on the data of the Mediterranean and California regions.

#### **2. The Method of Approximation of Interval Expert Estimates**

Let the seismotectonic properties of the region under study be represented by a set of spatial grid fields of features **X**1, **X**2, ..., **X***I*, and the values of the maximum possible magnitudes of earthquakes (*Mmax*) be represented by a sample set of expert estimates. The task is to find from these data the function *F*(**x**), which approximates the values of *Mmax* at the sample set, where **x** = (*x*1, ..., *xI*) is the vector with the values of the fields of features. The *Mmax* map is the *F*(**x**) values calculated for all grid nodes of a region.

The type of expert evaluation should be convenient and straightforward for unambiguous understanding by all participants of the expert survey and should enable the expert to formalize his knowledge about the value of the forecast fully. These requirements correspond to interval expert estimates:

$$Q\_{qn} = (m\_{qn}^{(1)}, m\_{qn}^{(2)}, w\_{qn}^{(1)}, w\_{qn}^{(2)}),\tag{1}$$

where *m*(1) *qn* , *<sup>m</sup>*(2) *qn* are the interval boundaries within which all the values of *Mmax* at the point *n* are the most probable and equally possible, in the opinion of the *q*-th expert, *m*(1) *qn* <sup>≤</sup> *<sup>m</sup>*(2) *qn* ; *<sup>w</sup>*(1) *qn* <sup>&</sup>gt; 0, *<sup>w</sup>*(2) *qn* > 0 are the weighs on which the *q*-th expert indicates the degree of his confidence in the possibility that the value of *Mmax* may be less or greater than the corresponding interval boundary *<sup>m</sup>*(1) *qn* or *<sup>m</sup>*(2) *qn* .

We can assume that the expert estimate *Q* corresponds to some function of the subjective probability density *f*(*Y*, *Q*), which reflects the expert's opinion about the value of *Y* at a given sample point. This function takes a constant value within the interval [*m*(1), *m*(2)] and decreases with the weights *w*(1) and *w*(2) respectively to the left and right of the interval boundaries:

$$f(\boldsymbol{Y}, \boldsymbol{Q}) = \mathbb{C} \cdot \exp\{- (\boldsymbol{w}^{(1)} \frac{|\boldsymbol{m}^{(1)} - \boldsymbol{Y}| + \boldsymbol{m}^{(1)} - \boldsymbol{Y}}{2} + \boldsymbol{w}^{(2)} \frac{|\boldsymbol{m}^{(2)} - \boldsymbol{Y}| - \boldsymbol{m}^{(2)} + \boldsymbol{Y}}{2})^{p}\},\tag{2}$$

where *p* ≥ 1, and *C* is defined by the condition ∞ −∞ *f*(*Y*, *Q*)*dy* = 1.

Suppose that there is a training sample {*Qqn*, **x***n*}, where *q* and *n* represent the expert and sample number. It is required to approximate the function *Y*(*x*) in a certain class of functions *F*(**x**, *θ*) : *θ* ∈ **Θ**, where **Θ** is the domain of admissible values of the vector *θ*.

Let's replace *Y* in (2) with the value of the forecast function *F*(**x**, *θ*) and consider the function

$$r(\mathbf{x}, \boldsymbol{\theta}) = -\ln f(\mathbf{Y}, \mathbf{Q}) + \ln \mathbf{C} = (w^{(1)}\frac{|m^{(1)} - \mathbf{Y}| + m^{(1)} - \mathbf{Y}}{2} + w^{(2)}\frac{|m^{(2)} - \mathbf{Y}| - m^{(2)} + \mathbf{Y}}{2})^p. \tag{3}$$

The function *r*(**x**, *θ*) determines the penalty for the inaccuracy of the approximation of the expert judgment *Q* by the value *F* of the forecast function. To estimate *θ*, the average penalty on the set is minimized. The estimation has the form

$$\theta = \arg\min\_{\theta \in \Theta} \sum\_{n} \sum\_{q} r(F(\mathbf{x}\_n, \theta), Q\_{qn}). \tag{4}$$

It is obvious that if the forecast function *F*(**x**, *θ*) is linear in the parameters, then the functional (4) is convex. If the domain **Θ** of admissible values of the vector is also convex, then it is possible to use iterative gradient algorithms for estimation.

It is easy to see from (3) and (4) that in case of *m*(1) = *m*(2) and *w*(1) = *w*(2) for all expert estimates the estimation algorithm (4) coincides with the method of the least absolute errors for *p* = 1, and with the method of least squares for *p* = 2. It was shown in [22] that under certain assumptions, the estimate (4) is an estimate of the maximum likelihood.

The method of approximation of interval expert estimates was repeatedly used to construct the maps of *Mmax* in a number of regions, in particular, Bulgaria [34], Caribbean and Middle America Region [35], Central Europe [12,36], Costa Rica [37], the Caucasus [21], and North Caucasus [38]. In these papers, the dependences *Mmax*(**x**) were always estimated in a class of the sum of piecewise linear functions of geological and geophysical features. This estimation allows one to interpret the *Mmax* map as the sum of non-linearly transformed fields of features.

For each of the above regions, from 10 to 200 geological and geophysical fields were analyzed. 3–4 of the most informative fields were selected from this set using the stepwise regression method. Prediction functions are the sum of piecewise linear dependencies on the values of these fields. The sum of nonlinearly transformed fields defines the *Mmax* field. This is convenient for the seismotectonic interpretation of the *Mmax* map.

Interpretation of the *Mmax* map by a specialist allows a qualitative assessment of the accuracy of determining the maximum possible seismic hazard but is not an assessment of its accuracy. By definition, *Mmax* values cannot be measured instrumentally. Statistical estimates of *Mmax* estimates are possible only in areas with sufficiently high seismic activity. In the considered method, the values of *Mmax* are replaced by interval expert estimates. These estimates are approximated by a nonlinear function of geological and geophysical fields. The accuracy of the *Mmax* forecast is determined by the deviations of the *Mmax* forecast map values from expert estimates. For the above regions, from 100 to 400 expert evaluations were used. The number of parameters estimated during training ranged from eight to 14 in each region. For those regions, who did not participate in the training, the average approximation errors of expert estimates are in the range from 0.2 to 0.34.

#### **3. Method of the Minimum Area of Alarm**

Let the properties of the seismic process are described by the spatial and spatio-temporal fields of features in a single coordinate grid with a step Δ*x* × Δ*y* × Δ*t*. The values of these fields at the nodes of the grid *<sup>n</sup>* <sup>=</sup> 1, ... , *<sup>N</sup>* correspond to the vectors of the *<sup>I</sup>*-dimensional feature space **<sup>f</sup>**(*n*) <sup>=</sup> { *<sup>f</sup>* (*n*) *<sup>i</sup>* }. A spatio-temporal forecast field **Φ** is a function of the fields of features. It is trained using retrospective data: (1) A sample set of target earthquakes *q* = 1, ... , *Q* with the magnitudes *M* ≥ *M*<sup>∗</sup> and (2) a set of grid fields of features *Fi*, *i* = 1, ... , *I*, which describes spatial (quasi-stationary) and spatio-temporal (dynamic) properties of the seismotectonic process.

The method of the minimum area of alarm uses the following data model.


We will call the base vectors of the feature space the vectors for which *<sup>f</sup>* <sup>≥</sup> *<sup>f</sup>*(*q*) componentwise. The nodes of the grid of the forecast field **<sup>Φ</sup>** with the values *<sup>φ</sup>* <sup>≥</sup> *<sup>φ</sup>*(*q*) we will call the base nodes of the forecast field.

From the assumption that anomalous refers only to the largest values of the fields of features and the monotonicity condition, it follows that the earthquake forecast can be carried out using the simplest threshold decision rule. If the value of the forecast field *φ*(*n*)≥*θ*, then spatio-temporal alarm cylinders are created at all base nodes of the forecast field with the values *<sup>φ</sup>* <sup>≥</sup> *<sup>φ</sup>*(*q*). The alarm cylinder of the grid node *n* with the coordinates (*x*(*n*), *y*(*n*), *t* (*n*)) has the center of the base in the node (*x*(*n*), *y*(*n*), *t* (*n*)), the base radius *R*, and the element [(*x*(*n*), *y*(*n*), *t* (*n*)),(*x*(*n*), *y*(*n*), *t* (*n*))]. From this, it follows that for a given value of the threshold *θ* an earthquake with the epicenter coordinates and time (*x*∗, *y*∗, *t* ∗) will be detected if and only if the cylinder with the center of the base (*x*∗, *y*∗, *t* ∗), the radius *R*, and the element [(*x*∗, *y*∗, *t* <sup>∗</sup> − *T*),(*x*∗, *y*∗, *t* ∗)] contains at least one grid node with the value *<sup>φ</sup>*(*n*) <sup>≥</sup> *<sup>θ</sup>*. This cylinder will be called a precursor cylinder.

The alarm field detects an earthquake if its epicenter falls within an area consisting of a combination of alarm cylinders (alarm area). The quality of the forecast field at threshold *θ* is determined by two indicators: (1) The fraction of correctly detected events *Q*∗(*θ*) from all *Q* events *U*(*θ*) = *Q*∗(*θ*)/*Q* (probability of detection) and (2) the fraction of number of grid nodes, falling in the alarm area *L*∗(*θ*), from the number of all grid nodes *L* of the analyzed area *V*(*θ*) = *L*∗(*θ*)/*L* (alarm volume).

For training, we have a set of target events with magnitudes *M* ≥ *M*<sup>∗</sup> and a set of fields. At the first step, the algorithm should move from a set of target earthquakes to a set of target earthquake precursors. A precursor of the earthquake *q* is the vector *f*(*q*) of a feature space which has a minimum volume of alarm *v*(*q*) = *L*(*q*)/*L* among all vectors corresponding to the grid nodes of the precursor cylinder of the event *q*, where *L* is the number of all grid nodes of the analyzed area, *L*(*q*) is the number of nodes in the grid of the alarm area generated by the base points of the vector *f*(*q*) .

The algorithm for constructing the forecast field is nonparametric. There are the three most important versions of the algorithm. The first version of the algorithm is to construct the forecast field so that when the threshold *θ* decreases, the training earthquakes are detected in the sequence in which the corresponding alarm volumes increase *<sup>v</sup>*(*Q*) <sup>≤</sup> *<sup>v</sup>*(*Q*−1) <sup>≤</sup> ... <sup>≤</sup> *<sup>v</sup>*(2) <sup>≤</sup> *<sup>v</sup>*(1) (this version is selected for testing). The version consists of the following steps.


Obviously, the choice of the order of the earthquake precursors at the 2nd step of the algorithm determines the dependence *U*(*V*) obtained from the forecast field. The 2nd version of the algorithm makes it possible to optimize the forecast field so that when the next target earthquake is detected, the alarm volume increases by a minimum value. To do this, one should arrange the precursors so that, when changing from event detection *q* + 1 to event *q*, the increase in alarm volume is minimal. Here, at each transition from the previously selected event *q* + 1 to *q*, a small search through the remaining *q* events is required. The 3rd version of the algorithm allows one to optimize the forecast field so that it detects the maximum number of target earthquakes with a total alarm volume of less than or equal to the predetermined value. In this case, you need to perform a full search on the selected number of events. The 3rd version of the algorithm allows optimizing the forecast field so that it detects the maximum number of target earthquakes, provided that the total alarm volume does not exceed the specified value. In this case, you need to perform a full search for the selected number of events.

#### **4. Testing**

The purpose of testing is to verify the proposed method of the forecast. Testing is carried out in accordance with the known characteristics of the catalog of earthquakes. Exploring the possibility to improve the quality of the forecast using a wider set of characteristics of earthquake catalogs or by adding other sources of input data is beyond the scope of this work.

The method of minimum area of alarm was tested on the platform of automatic earthquake forecast (http://distcomp.ru/geo/prognosis). The system tests the data with a constant step Δ*t*. On each step (at time *t*) the raster fields of features are computed, the alarm area is trained based on data before the time *t*, and the system tests for time since *t* till *t* + Δ*t* if the alarm area covers an epicenter of the target earthquake. Then at time *t* + Δ*t*, the training time is increased by Δ*t*, the alarm zone is updated and the test is repeated.

Testing of the forecast method should provide an opportunity to compare different methods of solving the problem on the same indicators of the forecast quality. In this method, we use two quality indicators: The probability of detecting the target events from the test interval *U* = *Q*∗/*Q* and the volume of alarm *V* = *L*∗/*L*. The number of target events *Q* is determined by a set of test samples, the number of target events detected *Q*∗ is determined by the results of the forecast, the analysis area and its size *L* is selected at the beginning of the test, the size of the alarm zone *L*∗ is determined by the training data.

In the following test experiments, the area of analysis was constructed in the following way: Any point is included in the area if in a circle around this point with radius *R* = 100 km for the period 1984–1993 there are more than 300 earthquake epicenters. This condition allows one to select a seismically active area for analysis but does not ensure its seismic homogeneity. Therefore, the indicator of the volume of alarm obtained during testing should be considered only in the context of the selected area of analysis. At the same time, the choice of the field of analysis according to a formal rule makes it possible to compare the results of the forecast obtained using various methods and according to different data.

One way to assess the quality of a forecast is to compare a regular forecast obtained by the algorithm under analysis (regular forecast) with a random one. We will assume that the forecast is random if the values of the forecast field are selected from a segment in accordance with a uniform distribution. Obviously, for this probabilistic model, the alarm volume *Vr* is equal to the probability of a random prediction *Ur*. It follows from this that comparing the probability of a regular forecast *U* with the probability of a random forecast *Ur* for the same alarm volumes *V* = *Vr* is equivalent to comparing *U* with the corresponding alarm volume *V*. If, at the same time, a sample of target events were cleaned of aftershocks and foreshocks, then by proposing the independence of target events and using the binomial distribution model for them, we could build a confidence interval for estimating *U*.

In the number of articles, the results of a regular forecast are compared with the results of a forecast by a stationary field. In papers [39–41] the regular forecast is compared with the forecast by the 2D field of seismic activity (or earthquake epicenter density). The result of the comparison makes it possible to evaluate the efficiency of a regular forecast in relation to the forecast by the field *F*, which is based only on the spatial heterogeneity of the seismic process. Comparison of results can be done in two ways. In one method the probabilities of regular prediction of *Q* target earthquakes are compared with the results of predicting the same earthquakes by a stationary field *F* (for example, *F* is a 2D field of the earthquake epicenter density). Another method uses the Gutenberg–Richter model [42]. In the beginning, the catalog of earthquakes with the following conditions is constructed: (1) The epicenters of earthquakes are in the area of analysis, (2) the magnitudes exceed the representative, and (3) the depth of the epicenters does not exceed the values specified for the target earthquake. It is assumed that *b*-value is the same for the entire area of analysis and the earthquake catalog agree well with the spatial distribution of seismicity. The alarm field *V*(*θ*) is calculated by the field *F*. Then, in accordance with the alarm field, the dependence *N*∗(*V*), is calculated, where *N*∗ is the number of epicenters in the alarm zone, and *V* is the alarm volume. The dependence *N*∗(*V*) is normalized to the number of all *N* epicenters in the analyzed area. According to the Gutenberg–Richter law and the assumption *b* =const, we have *N*∗(*V*) = *C* exp(*d*<sup>∗</sup> − *bm*), *N*(*V*) = *C* exp(*d* − *bm*), and *μ*(*V*) = *N*∗/*N* = exp(*d*<sup>∗</sup> − *d*). Thus, the value *μ*(*V*) does not depend on the magnitude of earthquakes. It shows the proportion of earthquakes with a magnitude higher than a given, which fall into the alarm zone. Consequently, the value of *μ*(*V*) in the scope of our model is equal to the probability of forecasting the target earthquakes using the stationary field *F*. If the field *F* is the density field of the earthquake epicenters, then the field obtained according to the Gutenberg–Richter law is denoted by the letter *μ*.

Testing was performed for two regions: The Mediterranean and California. The Mediterranean region: 10◦–30◦ E, 34◦–47◦ N. Input data: Earthquakes for the period from 27.05.1983 till 14.02.2018 with magnitudes *M* ≥ 2.7 and depths of hypocenters *H* ≤ 160 km from the International Seismological Centre catalog (see Materials and Methods). Target earthquakes: Magnitudes *M* ≥ 6.0 and hypocenter depths *H* ≤ 60 km. California region: 126◦–114◦ W, 32◦–43◦ N. Input data: Earthquakes for the period of 01.01.1983–15.02.2018 with magnitudes *M* ≥ 2.0 and depths of hypocenters *H* ≤ 160 km from the NEIC USGS catalog (see Materials and Methods). For the forecast, the target earthquakes with magnitudes *M* ≥ 5.7 have been selected.

The following six fields of features were analyzed for forecasting:


The estimation of 3D fields of *F*<sup>1</sup> and *F*<sup>2</sup> is performed with the method of local kernel regression. The kernel function for the *n*-th earthquake has the form *Kn* = [cosh2(*rn*/*R*)<sup>2</sup> cosh2(*tn*/*T*)]−1, where *rn* < *R*, *tn* < *T* are the distance and time interval between the *n*-th epicenter of the earthquake and the node of the 3D grid of the field, = 2, *R* = 50 km, *T* = 100 days for *F*<sup>1</sup> and *R* = 100 km, *T* = 730 days for *F*2.


To estimate the field of *F*3, *F*4, *F*5, the Student's *t*-statistic was used, which is defined for each grid node as the ratio of the difference of average values of the current (196 days) and background (3650 days) intervals to the standard deviation of this difference. Positive t-statistics values correspond to higher values on the test interval.

• *F*<sup>6</sup> is the 2D field of the density of earthquake epicenters: Kernel smoothing in the interval 1988–2008 the parameter *R* = 50 km.

The grid fields for the Mediterranean were calculated in a grid step Δ*x* × Δ*y* × Δ*t* = 0.2◦ × 0.13◦ × 49 days. The forecast field was trained from 1998 until the next step of the forecast after 2008. The radius of the alarm cylinder is *R* = 20 km, and the element is *T* = 50 days. Testing is performed in 2008–2019. There are 11 target earthquakes in the analysis area. We used the method of stepwise selection to find the most informative fields of features. The algorithm selected the *F*<sup>3</sup> and *F*<sup>6</sup> fields to construct the alarm field.

We compare the earthquake prediction probabilities obtained using different fields of features in Table 1: *U*<sup>1</sup> is the forecast probability using the earthquake density field 2D (*F*6), *U*<sup>2</sup> is the probability using the 3D field of negative earthquake density anomalies (*F*3), *μ*(*V*) is the probability of forecast by 2D field of earthquake epicenters density, obtained using the Gutenberg–Richter model, and *U*<sup>3</sup> is the probability using *F*<sup>3</sup> and *F*<sup>6</sup> fields. We can see that the highest probability of a successful forecast occurs when the fields *F*<sup>3</sup> and *F*<sup>6</sup> are used together. When *V* = 0.2 (*Ur* = 0.2), the ratios for the prediction probability obtained with *F*<sup>3</sup> and *F*<sup>6</sup> fields to the prediction probabilities obtained with 2D earthquake density field (*F*6) and for the field calculated using the Gutenberg–Richter ratio, are equal respectively *U*3(0.2)/*U*1(0.2) = 0.91/0.64 = 1.49 and *U*3(0.2)/*μ*(0.2) = 0.91/0.41 = 2.2. Table 1 shows the values of two types of alarm volumes: *Vlearn* is the alarm volume received in accordance with the training data, and *Vtest* is the alarm volume corresponding to the alarm volume *Vlearn* but observed on the test data. You can see that when testing in almost all cases, except for testing the 2D-field *F*6, the volumes of *Vtest* are greater than *Vlearn*. This is explained by the fact that the number of recorded earthquakes in a region changes over time (Figure 1). The number of earthquakes is influenced by the development of a seismic network and natural changes in the seismic process. The Figure 1 shows that the number of earthquakes increases significantly in the test interval. The same anomaly appears on the plot of the time series of average values of the density of earthquake epicenters throughout the analysis area (Figure 2). An increase in the density of earthquake epicenters leads to an increase in the field values of the *F*<sup>3</sup> function, which ultimately leads to an increase in the volume of anxiety during testing.

The grid fields for California were calculated in a grid step Δ*x* × Δ*y* × Δ*t* = 0.125◦ × 0.11◦ × 49 days. The radius of the alarm cylinder is *R* = 14 km, and the element is *T* = 100 days. Testing of the earthquake forecast was performed for the interval 2009–2018. There were nine target earthquakes. The algorithm selected tree fields of features for the construction of the alarm field: *F*4, *F*5, and *F*6.

Table 2 shows the probabilities of earthquake forecast for California.

Figures 3 and 4 show the test results for both regions. They depicted polygons selected as the area of analysis, and circles are the target epicenters of earthquakes in 2009–2018 with *M* ≥ 6.0 for Mediterranean and *M* ≥ 5.7 for California.

*Vlearn***: Alarm Volumes for Learning Interval 0.01 0.05 0.1 0.15 0.2** Test indicators *Vtest Utest Vtest Utest Vtest Utest Vtest Utest Vtest Utest* Field *F***6** 0.00 0.00 0.03 0.09 0.10 0.09 0.15 0.45 0.20 0.64 Field *F***3** 0.00 0.00 0.02 0.09 0.20 0.36 0.29 0.45 0.35 0.64 *μ***(***V***)**: probability for the field *F***<sup>6</sup>** obtained by the model - 0.00 - 0.1 - 0.26 - 0.35 - 0.41 Fields *F***3** and *F***6** 0.03 0.00 0.08 0.36 0.15 0.45 0.24 0.55 0.32 0.91

**Table 1.** Comparison of the probabilities of earthquake forecast for the Mediterranean.

**Figure 1.** Histogram of the number of earthquakes with a magnitude greater than 2.7 and a depth of hypocenters less than 160 km.

1985 1990 1995 2000 2005 2010 2015 2020

**Figure 2.** Time series of average values of the density of earthquake epicenters throughout the analysis area (trend line).


**Table 2.** Probabilities of earthquake forecast for California.

**Figure 3.** Area of analysis (marked with thick black line) and tested target epicenters of earthquakes in 2009–2018 in the Mediterranean region. Shades of grey indicates the minimum volume of alarm with which the epicenter was forecasted. Darkness of grey decreases in accordance with the volume of the alarm: 0.05, 0.1, 0.15, 0.2. A white color indicates that an earthquake is not forecasted with an alarm volume of less than 0.2.

**Figure 4.** Area of analysis (marked with thick black line) and tested target epicenters of earthquakes in 2009–2018 in California. Shades of grey indicates the minimum volume of alarm with which the epicenter was forecasted. Darkness of grey decreases in accordance with the volume of the alarm: 0.05, 0.1, 0.15, 0.2. A white color indicates that an earthquake is not forecasted with an alarm volume of less than 0.2.

#### **5. Discussion**

The method of approximating the interval expert estimates compiles the *Mmax* map, assuming a repetition of strong earthquakes and the existence of a relationship between *Mmax* and the properties of the geological environment *x* = (*x*1, ... , *xI*). At first, experts independently estimate *Mmax* values in a set of the most studied points of a region. The algorithm approximates the dependence of *Mmax* on geological and geophysical features by the function *F*(*x*). The dependence *F*(*x*) is defined as the sum of the non-linear functions of each of the feature. The *Mmax* map is the *F*(*x*) values calculated for the whole region. The presence of the formal forecast rule *F*(*x*) allows the expert to study the contribution to the forecast of each feature and interpret the map as the sum of the nonlinearly transformed feature fields.

The method of the minimum area of alarm solves the problem of one-class classification. The method algorithm has two peculiar properties. The first relates to the data model. The model postulates two properties of anomalous objects: (1) Anomalous objects are unlikely, and some of their properties take values close to the maximum (or minimum) among the sample, and (2) the vectors of the space of features, which are componentwise larger (or smaller) of the vector corresponding to the anomalous object, can also be anomalous objects. Both these properties seem sufficiently natural. This model allows one to build a classification rule from a set of anomalous objects. In this case, normal objects are taken into account statistically through the probability of detecting anomalous objects by a random forecast. The second difference is that the algorithm allows constructing a forecast function that optimizes the probability of detecting anomalous objects in the training sample if the probability of a random forecast is not more than the predetermined value.

#### **6. Conclusions**

We considered two machine learning methods and their implementations to seismic hazard forecast. The method of approximation of interval expert estimates of *Mmax* demonstrated good seismic zoning for many seismically active regions. The method of the minimum area of alarm is the basis of an automatic earthquake forecast system. The considered results of testing suggest that the method and the forecast system might contribute to advance in the problems of earthquake forecasting.

#### **7. Supplement**

The method of minimum area of alarm is the basis of an automated web-based platform that systematically forecasts target earthquakes. We presented the results of testing the approach to earthquake prediction in the Mediterranean and Californian regions. The goal of the test was to analyze the approach, the machine learning method, and the earthquake prediction platform. For the tests, ordinary parameters of earthquake catalogs were used. The testing was performed on data that did not participate in the training and showed a satisfactory forecast quality for both regions. The web-based platform has been launched and automatically calculates the seismic hazard fields from February 2018 [43]. During the time from 1 February 2018 to 8 July 2019, four target earthquakes occurred in these regions. In the Mediterranean region, two epicenters were predicted and fell into an area with an alarm volume of up to 15% (Figure 5), and in the California region, two epicenters fell into an area with an alarm volume greater than 20% and were not predicted (Figure 6).

**Figure 5.** Screenshot of the web-based platform working window in the Mediterranean region: The map shows the alarm zone for earthquakes with a magnitude *M* ≥ 6.0 and the predicted epicenters of earthquakes of 25 October 2018 with a magnitude of 6.6 and 30 October 2018 with a magnitude of 6.2 calculated for training according to data up to 26 September 2017. The palette shows areas with different alarm volumes in percent: 1%, 5%, 10%, 15%, 20%.

**Figure 6.** Screenshot of the web-based platform working window: The map shows the alarm zone for earthquakes with a magnitude *M* ≥ 5.7 and the epicenters of earthquakes of 4 July 2019 with a magnitude of 6.4 and 6 July 2019 with a magnitude of 7.1 calculated during training according to data until 3 June 2019. The palette shows areas with different alarm volumes in percent: 1%, 5%, 10%, 15%, 20%.

#### **8. Materials and Methods**

International Seismological Centre [44] was searched using http://www.isc.ac.uk/iscbulleti n/search/bulletin/ (last accessed on 14 February 2018). This online catalog was selected for its robustness and universality. It combines data from a lot of catalogs, and every earthquake with a magnitude more than three is manually checked. NEIC USGS catalog [45] was searched using https://earthquake.usgs.gov/earthquakes/feed/ (last accessed on 15 February 2018). This catalog was chosen because of numerous registered small earthquakes (magnitude of completeness less than 1.5) in California. Plots were made using the GeoTime 3 (www.geo.iitp.ru/GT3; [8]).

**Author Contributions:** Conceptualization, V.G.G.; Formal analysis, V.G.G. and A.B.D.; Investigation, V.G.G.; Methodology, A.B.D.; Project administration, V.G.G.; Visualization, A.B.D.; Software, A.B.D.; Writing—original draft, V.G.G.; Writing—review & editing, V.G.G. and A.B.D.

**Funding:** The paper is supported by the Russian Foundation for Basic Research, project 17-07-00494.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Suitability Analysis for the Emergency Shelters Allocation after an Earthquake in Japan**

#### **Tetsuya Akamatsu \* and Kayoko Yamamoto**

Graduate School of Informatics of Engineering, University of Electro-Communications, Chofu, Tokyo 182-8585, Japan

**\*** Correspondence: a1930002@edu.cc.uec.ac.jp; Tel.: +81-42-443-5728

Received: 12 June 2019; Accepted: 24 July 2019; Published: 30 July 2019

**Abstract:** Japan is exposed to several natural phenomena; the damages caused by earthquakes were enormous in particular. An emergency shelter is a place for people to temporarily live when they cannot remain in their previous homes, and it is necessary for each community to, respectively, allocate such facilities in Japan. There are the districts that are short of such facilities, especially in rural and suburban areas, because emergency shelters mainly concentrate near large-scale stations and city centres in Japan. Against these backdrops, using geographic information systems (GIS), an applied statistical method and public open data related to population and emergency shelters, the present research aims to quantitatively conduct a suitability analysis for the emergency shelters allocation after an earthquake in Japan. Based on the results, the present research grasps the districts that are short of emergency shelters, and visually shows the places where such facilities should be newly established on the digital map of GIS. Additionally, the assessment method is reproducible in the spatial and temporal dimension. It is necessary to create an original data related to emergency shelters to raise the reliability of the results, as the present research has the limitation of data availability.

**Keywords:** emergency shelter; earthquake; disaster; weighted coefficient; suitability analysis; geographic information systems (GIS)

#### **1. Introduction**

In Japan, various disasters have been frequently occurred, and the damages especially caused by earthquakes were enormous. According to Yamaga (2019) [1], the earthquake catalog search conducted by the United States (U.S.) Geological Survey (USGS), 14.8% of the earthquakes larger than magnitude 5.0, which occurred worldwide from 2010 to 2015 concentrated in Japan. Therefore, it can be said that Japan is one of the few earthquake-ridden countries.

Meanwhile, due to the serious damages that are caused by disasters, victims who cannot remain in their homes stay in emergency shelters until they can settle in new places. The maximum distance elderly people and children can evacuate on foot is called as the maximum evacuation distance, and it is estimated to be around 1.5 to 2.0 km (Dictionary of Housing Terms, 2019) [2]. However, the distance that elderly people over 70 years old can walk is 879 m, which is significantly shorter when compared with those of the other age groups (Cabinet Office, 2019) [3]. The physical abilities of elderly people over 70 years old drastically decrease, which makes long-distance evacuations extremely difficult. Additionally, in Japan, as emergency shelters mainly concentrate near large-scale stations and city centres, there are districts that are short of such facilities, especially in rural and suburban areas. Additionally, focusing on the communities, Civiletti et al. (2016) [4] identified the institutional and social decisions that increase the resilience of the communities that are exposed at risk, by analyzing observations during seismic sequences that occurred in Italy in the last decades, and Cerchiello et al. (2018) [5] addressed the assessment of the social vulnerability and resilience level of the city of Nablus,

Palestine. Mojica et al. (2010) [6] and Awaotona (2012) [7] pointed out that it is necessary to take the countermeasures against disasters while considering vulnerable people (elderly and disabled people). However, in Japan, due to the Basic Act on Disaster Control Measures, it is possible to already take the measures for disabled people requiring assistance at the time of disasters in each municipality (Cabinet Office, 2013) [8]. Additionally, in Japan, the aging rate of elderly population who are over 65 years old is 26.7% in 2015, which is the highest in the world. For reference, the above rates are 22.4% in Italy, 21.2% in German, 19.9% in Sweden, and 19.1% in France (Ministry of Health, Labour and Welfare, 2016) [9]. Therefore, when considering the regional characteristics, such as the ratio of elderly people and population distribution, it is necessary to evaluate not only facility location, but also capacity, this is because elderly people have special needs. Additionally, it is also necessary to grasp the districts (made up of streets and towns, and the minimum level of Japanese communities), with a lack of emergency shelters.

Based on the social and academic background mentioned above, using Geographic Information Systems (GIS), which is an applied statistical method and public open data related to population and emergency shelters, and while considering regional characteristics, such as the ratio of elderly people and population distribution, the present research aims to conduct a suitability analysis for the emergency shelters allocation after an earthquake on a district scale in Japan. In the present research, the analysis will be conducted to examine whether emergency shelters are sufficient or not in the unit of districts, and where such facilities should be newly established. Additionally, referring the results of the preceding research studies both within and outside Japan introduced in the next section, the present research will develop a method for the suitability analysis for the emergency shelters allocation. Based on the results, the present research provides effective information that can be referred to determine the locations to establish new emergency shelters.

#### **2. Literature Review**

As mentioned in the previous section, Japan is exposed to several natural phenomena, and there is an accumulation of many research studies in the field related to emergency shelters. The present literature review is related to (1) research studies that are related to the suitability analysis of emergency shelters location, and (2) research studies related to the establishment of new emergency shelters. The following will introduce the major preceding research studies in the above two study areas, and discuss the originality of the present research in comparison with the others.

Regarding the representative research studies that are related to the suitability analysis of emergency shelters location, especially in Japan, Furihata et al. (1994) [10] evaluated the location of evacuation facilities taking the spatial patterns of distributed users into consideration and using GIS. Takeuchi et al. (2002) [11] proposed a shared area in the taking emergency shelter to have considered the approaching direction of the tsunami that is generated by a huge earthquake. Kongsomsaksakul et al. (2005) [12] proposed the optimal locations of emergency shelters for the flood evacuation planning, while using the combined distribution and assignment (CDA) model and genetic algorism (GA) in the U.S. Asou et al. (2007) [13] clarified evacuation behavior while using GA to propose the optimal arrangements of emergency shelters. Wei et al. (2008) [14] presented a diagnosis model for emergency shelter planning from the viewpoint of local people using GIS. Ng et al. (2010) [15] presented a hybrid bi-level model for the optimal shelter assignment in emergency evacuations in the U.S. Tai et al. (2010) [16] used six indicators to evaluate shelter and applied a spatial statistic model with local indicators of spatial association (LISA) to the evacuation choice of residents in the case of an earthquake in Taiwan. Kitajima (2013) [17] evaluated the location of emergency shelters using network analysis. Yu et al. (2016) [18] introduced a framework for the multi-criteria satisfaction assessment of the spatial distribution of urban emergency shelters while using a GIS-based analytic hierarchy process approach in Shanghai City, China. Vecere et al. (2016, 2017) [19,20] focused on a critical review of currently available methodologies and corresponding software packages that were specifically developed for estimating the number of displaced people and those who need public sheltering and temporary

housing at the time of the 2011 Christchurch earthquake in New Zealand. Ashish et al. (2017) [21] presented a hybrid algorithm for efficiently managing location and relocation projects, by proposing a hybrid multi-objective decision model that is based on analytic hierarchy process (AHP), fuzzy set theory, and goal programming approach, referring two real case studies of Nepal earthquake. Xu et al. (2018) [22] developed a multi-objective mathematical model with four groups of the objectives, allied with a modified particle swarm optimization algorithm to solve the location-allocation problem for earthquake shelter in Beijing City, China. Nozaki et al. (2019) [23] conducted the location analysis of tsunami emergency shelter while considering inhabitants' preparedness for the coming Nankai Trough Earthquake.

Regarding the representative research studies that were related to the establishment of new emergency shelters, especially in Japan, Yamada et al. (2004) [24] proposed the planning support system for the locations of emergency shelters from the viewpoint of residents. Notsuda et al. (2005) [25] considered the optimum location of new emergency shelters and the placement of evacuees to each facility that is based on the location and capacity of existing emergency shelters, as well as the population distribution of evacuees. Nakai et al. (2012) [26] and Ikenaga et al. (2017) [27] considered the possibility of using vacant houses as emergency shelters and temporary evacuation shelters in the time of accidents. Asano et al. (2013) [28] and Miyoshi et al. (2017) [29] considered the possibility of using public and private lands as emergency shelters and safety evacuation areas. Araki et al. (2017) [30] examined setting up the patterns of non-designated emergency shelters at the time of the Great East Japan Earthquake (2011), which were based on GIS analysis and interview surveys. Sasaki et al. (2018) [31] examined the possibility of approximately 70,000 temples nationwide, complementing emergency shelters and safety evacuation areas that are expected to be in shortage when many disaster-affected residents appear during widespread disasters. Umeki et al. (2019) [32] proposed a method to determine the location of emergency shelters to aim at the reduction of the evacuation time of all victims.

In other countries, Kar et al. (2008) [33] developed a GIS-based model to determine the site suitability of emergency shelters for hurricane events, and proposed the candidate places to arrange new ones in the U.S. Alçada-Almeida et al. (2009) [34] incorporated multi-objective model into a GIS-based decision support system to locate emergency shelters during major fires in Portugal. Park et al. (2012) [35] developed a method that applied genetic optimization to determine optimal tsunami shelter locations with the goal of reducing evacuation time, thereby maximizing the probability of survival for the population in a coastal community in the U.S. Li et al. (2012) [36] developed a scenario-based bi-level programming model to optimize the selection of shelter locations, with explicit consideration of a range of possible hurricane events and the evacuation needs under each of those events in the U.S. Anhorn et al. (2015) [37] proposed a methodology to rank the suitability of open spaces for contingency planning and the placement of emergency shelter while using GIS in the immediate aftermath of a disaster in Turkey. Bayram et al. (2018) [38] proposed a scenario-based two-stage stochastic evacuation planning model that optimally locates shelter sites and assigns evacuees to nearest shelters and to shortest paths within a tolerance degree to minimize the expected total evacuation time in Turkey. Xu et al. (2016) [39] proposed a multi-criteria constraint location model to select and analyze the candidate for emergency shelters, and to determine the location of new ones while using GIS in China.

From the above, especially in Japan, research studies that are related to the suitability analysis of emergency shelters location increased after the Great Hanshin earthquake (1995), and research studies that are related to the establishment of new emergency shelters remarkably increased after the Great East Japan Earthquake. Because, at the time of the Great East Earthquake, large earthquake, tsunami, and an accident at the nuclear power station occurred at the almost same time, and it was essential to keep sufficient emergency shelters for many victims.

The present research will reveal its effectiveness by newly proposing the use of facilities, such as temporary evacuation shelters in the time of accidents, safety evacuation areas, and large-scale retail stores that are not set up to be emergency shelters as new ones, unlike the preceding research studies that were related to the establishment of new emergency shelters, as discussed in detail in Section 3. Additionally, though the preceding research studies that are related to the suitability analysis of emergency shelters location derived the optimum locations for emergency shelters, the present research will demonstrate the originality by conducting a quantitative assessment of current location of emergency shelters, and proposing an assessment method to consider regional characteristics, such as the ratio of elderly population and population distribution. Therefore, it is possible to conduct a suitability analysis for the emergency shelters allocation that appropriately reflect current conditions of each district, and decide on the location of new emergency shelters by suggesting specific sites on a small spatial scale as the unit of districts while using GIS.

#### **3. Methods**

From the results of the research studies that were related to the suitability analysis of emergency shelters location in the previous section, for the suitability analysis for the emergency shelters allocation, it is necessary to consider regional characteristics, such as location and capacity of existing emergency shelters, as well as the population distribution of evacuees. It is evident that GIS-based analysis, spatial data analysis, and an applied statistical method are effective for the above purpose. Additionally, from the results of the research studies related to the establishment of new emergency shelters in Japan, it is necessary to set new candidate facilities in both public and private lands, especially in the districts that are short of emergency shelters. On the other hand, from the results of the above research studies in other countries, it is evident that GIS-based methods are effective in conducting a suitability analysis for the emergency shelters allocation. Referring these results, in this section, the method of the present study will be proposed.

#### *3.1. Framework and Process of Analysis*

GIS will be used for the method of suitability analysis for the emergency shelter allocation after an earthquake. The framework and process of analysis in the present research are as shown below.

(1) Calculation of weighted coefficients and creation of the distribution map of emergency shelters in each district

The weighted coefficient for each district will be calculated. A distribution map for all the emergency shelters will be created in the form of the digital map of GIS. Additionally, the facility scale (area) will be added to the data of each emergency shelter.

(2) Calculation of the linear distance between each district and each emergency shelter

The linear distance between the center of each district and each emergency shelter will be calculated while using the distribution map of emergency shelters in (1).

(3) Suitability analysis for the emergency shelters allocation using weighted coefficients

Using the weighted coefficients and distribution maps of emergency shelters in (1) as well as the linear distance between each district and each emergency shelter in (2), the assessment value for each district will be calculated. In the present research, analysis will be separately conducted using the four types of weighted coefficients related to the regional characteristics that are described in detail in the next section. When determining the locations to establish new emergency shelters, the important regional characteristics might vary with areas. Based on the results, the present research provides the effective information to assist the policy and decision makers in planning new emergency shelters.

#### *3.2. Assessment Method*

3.2.1. Calculation of Weighted Coefficients and Creation of the Distribution Map of Emergency Shelters in Each District

(1) Distribution map of different types of weighted coefficients for each district

In the present research, the degree of importance is indicated by the weighted coefficient and an assessment method will be proposed. There are four types of weighted coefficients, including (a) the specialization coefficient of elderly population, (b) the ratio of permanent population, (c) the ratio of elderly population, and (d) the specialization coefficient of population density, indicating various regional characteristics. Regarding (a) (d), the specialization coefficients are the values that reveals which area has the highest against the average rates of elderly population and population density in the city.

Regarding (a), Saino (1992) [40] set a "specialization coefficient of elderly people" which divided the rate of people who are over 65 years old in each district by the same rate in the entire area. However, the health span in Japan for 2016 was 72.14 for men and 74.79 for woman and both greatly exceed the age of 70, and it is increasing every year (the 11th Japan 21 (secondary) Health Promotion Committee, 2018) [41]. Therefore, the present research sets the specialization coefficient of elderly population who are over 75 years old, instead of those who are 65 years old as the target. The specialization coefficient of elderly population is a value that reveals which area has the highest aging rate of elderly population in the city. The aging rate of elderly population can be calculated with Equation (1). Additionally, the specialization coefficient of elderly population can be calculated while adopting this in Equation (2).

$$A\_i = \frac{p\_{75i}}{p\_i} \tag{1}$$

*Ai*: Aging rate of elderly population in district *i p*75*i*: Elderly population in district *i* (persons) *pi*: Population of district *i* (persons)

$$B\_{\bar{i}} = \frac{A\_{\bar{i}}}{A} \tag{2}$$

*Bi*: Specialization coefficient of elderly population in district *i*

*Ai*: Aging rate of elderly population in district *i* (%)

*A*: Aging rate of elderly population in the City (%)

Regarding (b), the ratio of permanent population for each district is important, as there is a greater need for emergency shelters in districts with a higher population. Equation (3) is adopted to calculate the ratio of permanent population.

$$C\_i = \frac{p\_i}{p} \tag{3}$$

*Ci*: Ratio of permanent population in district *i*

*pi*: Population of district *i* (persons)

*p*: Population of the City (persons)

Regarding (c), the population distribution of elderly group is not reflected in this weighted coefficients, as (a) the specialization coefficient of elderly population only indicates which district has the highest aging rate of elderly population. Therefore, it is weighted with the ratio of elderly population in order for the weighting to reflect the population distribution of elderly group. The ratio of elderly population can be calculated with Equation (4).

$$D\_i = \frac{p\_{75i}}{p\_{75}}\tag{4}$$

*Di*: Ratio of elderly population in district *i p*75*i*: Elderly population in district *i* (persons) *p*75: Elderly population in the City (persons)

Regarding (d), districts that are larger than others and are around large-scale stations and city centres generally have larger populations. Therefore, it is important to take the population density into consideration, as it is necessary to prioritize the establishment of new emergency shelters in districts that are short of such facilities due to the high population density. Thus, new emergency shelters can be established and victims can stay in such facilities near their homes even in the districts with high population density. The specialization coefficient of population density is a value that reveals which area has the highest aging rate of population density in the city. Population density can be calculated with Equation (5), and adopting this in Equation (6), the specialization coefficient of population density can be calculated.

$$F\_i = \frac{p\_i}{a\_i} \tag{5}$$

*Fi*: Population density in district *i* (*persons*/*km*2)

*pi*: Population of district *i* (*persons*)

*ai*: Area of district *i* (*km*2)

$$G\_{\bar{l}} = \frac{F\_{\bar{l}}}{F} \tag{6}$$

*Gi*: Specialization coefficient of population density in district *i*

*Fi*: Population density in district *i* (*persons*/*km*2)

*F*: Population density in the City (*persons*/*km*2)

(2) Creating distribution maps of emergency shelters

The distribution of emergency shelters will be displayed on the digital map of GIS. In the present research, the facilities, such as temporary evacuation shelters in the time of accidents and safety evacuation areas, are newly added to emergency shelters. These facilities were actually utilized as emergency shelters after the occurrence of disasters in the past. However, it is necessary to assume that an emergency shelter should be established as a temporary house in safety evacuation areas that are located outdoors. Based on Shigenobu et al. (2013) [42], the total size (area) of temporary houses can be calculated by multiplying the size of safety evacuation areas by the area ratio that is available for temporary houses. As there are multiple types of safety evacuation areas, the area ratio according to each type is set, as shown in Table 1. According to the Basic Act on Disaster Control Measures in Japan, it is possible to establish temporary houses in safety evacuation areas, such as open spaces and urban greens spaces that are shown in Table 1. Additionally, temporary houses were actually established in the above safety evacuation areas after the occurrence of disasters in the past. Additionally, after the Great East Japan Earthquake, over 2,000 people evacuated and stayed at large-scale retail stores in Ishinomaki City of Miyagi Prefecture. Therefore, since such large-scale retail stores are extremely wide and have a high capacity, they are also added as new emergency shelters.

**Table 1.** Area ratio available for temporary houses of safety evacuation areas.


3.2.2. Calculation of the Linear Distance between Each District and Each Emergency Shelter

In the present research, the linear distance between each district and each emergency shelter, and the scale of each emergency shelter will be used when conducting analyses. Regarding the former, the evacuations can be quickly completed in areas that are close to emergency shelters. Additionally, regarding the latter, the larger an emergency shelter is, the more people it can contain. In order to calculate the former, one of the ArcGIS Pro analysis tools, called the "Generate Near Table", is used. Using this tool, on the digital map of GIS, a proximity feature (emergency shelter) that is within 879 m from one or more features (center of each district) is selected, and *n* value is obtained. For each district, the *n* value will be the average of the number of emergency shelters within 879 m, which is the average walking distance for those in their 70's. Subsequently, using this tool again, the maximum number of features (*n* value) is set for each district, and the linear distance between each district and each emergency shelter is calculated.

#### 3.2.3. Suitability Analysis for the Emergency Shelters Allocation

The present research develops an assessment method that is based on the *p*-median problem. The *p*-median problem, which is one of the facility location problems, places facilities by minimizing the total sum of the linear distance from users to their nearest facility, and it can be modeled, as shown in Equation (7). This model derives the optimum location that lessens the load for users in all districts as much as possible by changing *Xij*. Equation (7) is changed to Equation (8) in order to respond the purpose of the present research.

$$\min\_{\mathcal{X}\_{ij}} Z = \sum\_{i} \sum\_{j} w\_{i} d\_{ij} \mathcal{X}\_{ij} \tag{7}$$

*Xij* ∈ {0, 1}: allocation to facility *j* in district *i*

*wi*: Demand in district *i*

*dij*: Linear distance from district *i* to facility *j*

$$\text{minimize } H\_{\text{i}} = w\_{\text{i}} \times \frac{1}{n} \sum\_{j=1}^{n} \frac{d\_{\text{ij}}}{\sqrt{s\_{\text{ij}}}} \tag{8}$$


Specifically, *Xij* of Equation (7) is removed in order to fix the locations of emergency shelters. Additionally, by calculating the assessment value for each district, it is possible to quantitatively grasp the emergency shelters allocation in each district. Additionally, as far as the originality of the present research, the average of the linear distance to the No. *i* closest emergency shelter in each district *dij*, divided by the square root of the facility scale of the emergency shelter *sij*, will be weighted with the coefficient *wi*, and this will be the assessment value for that district.

As shown in the previous section, the four types of weighted coefficients, the linear distance between each district and each emergency shelter, and the facility scale of emergency shelters will be applied to Equation (8), the assessment value of each district will be calculated, and the results will be displayed on the digital maps of GIS. When not using a weighted coefficient, analyses are conducted based on the linear distance between each district and each emergency shelter, and the scale of such facilities.

#### 3.2.4. Application of Assessment Method

The following two types of comparisons will be conducted in the discussion section in order to verify the validity of the assessment method proposed in the present research (Section 6).

(1) Comparison of the results between multiple areas with different regional characteristics

In the present research, multiple areas with different regional characteristics will be selected as target areas, and the suitability of emergency shelters allocation will be analyzed. Therefore, it is possible to verify the validity of the assessment method by comparing results of the selected areas with and without the four types of weighted coefficients that are related to the regional characteristics.

(2) Comparison of the results with and without weighted coefficients

The assessment method in the present research will adopt the four types of weighted coefficients as introduced in Section 3.2.1. Therefore, it is also possible to verify the validity of assessment method by comparing results with and without the four types of weighted coefficients that are related to the regional characteristics, focusing on the increase and decrease of the assessment values for both cases.

#### **4. Selection of Target Areas and the Data Processing**

#### *4.1. Selection of Target Areas*

For the analysis in the present research, it is essential to obtain the public open data that are related to population and emergency shelters, which will be introduced in the next section. However, in present Japan, very few cities make public the data related to emergency shelters. Among the cities that make public the above public open data, the four target areas selected for the present research are Chofu and Fuchu Cities in Tokyo metropolis, Toda City in Saitama Prefecture, and Shobara City in Hiroshima Prefecture. Table 2 describes the outlines of target areas. As described below, there are differences concerning the four types of weighted coefficients that are related to the regional characteristics introduced in Section 3.2.1 and the distribution of emergency shelters among these four cities.

Specifically, while the aging rate of elderly population (population aged over 75 years old) in Shobara City is significantly high, the rates in Toda City re fairly lower than the national average (11.2%). Additionally, the total populations of these two cities are 35,575 and 139,616, respectively, which makes them medium and small-scale cities. On the other hand, Chofu and Fuchu Cities, which are residential areas in Tokyo Metropolis, are large-scale cities with a population of 230,303 and 260,116, respectively. However, the distribution conditions of emergency shelters in these two cities are extremely different. Specifically, though emergency shelters are evenly distributed in the entire Chofu City, they are mostly located in the central part of Fuchu City.

According to the "Ordinance Covering Measures for Stranded Persons" in Tokyo Metropolis (2013), in the time of accidents, it is necessary for the companies in central Tokyo to keep employees at the workplaces. Therefore, central Tokyo, which has a concentration of companies as well as an extremely high daytime population, is not considered to be a target area in the present research. Additionally, according to the Basic Act on Disaster Control Measures in Japan, it is preferable for victims to stay in emergency shelters near their homes at the time of earthquakes. In the case of earthquakes, aftershocks follow for a while and additional damages might be generated, and it is dangerous for victims to still stay in their homes, even in good maintenance conditions that were constructed according to the Building Standards Act in Japan. Additionally, if victims stay together in emergency shelters near their homes, it is appropriate for rescue parties and government administrators to accurately confirm their safety and efficiently deliver relief goods to them.



#### *4.2. Data Processing*

Table 3 shows the utilized data and sources, and the utilization method of data in the present research.


**Table 3.** List of utilized data.

#### **5. Results**

The section will, respectively, describe the distinctive features of the result for each target area, referring Figures 1–20. These figures show the results with and without the four types of weighted coefficients that are related to the regional characteristics described in Section 3.2.1. These results are based on the suitability analysis for the emergency shelters allocation in the unit of districts of each target area mentioned in Section 3.2.3. Referring to these figures, the comparison of the results with and without the four types of weighted coefficients will be easily conducted. Furthermore, in these figures, all of the values have no units. As clearly shown in Table 2, the area of Shobara City is tremendously larger than those of other three cities. Therefore, the scale of results for Shobara City is different from those for other three cities on the digital maps of GIS.

#### *5.1. Results for Chofu City*

Figures 1–5 shows the results and *n* value is 6 for Chofu City. As mentioned in Section 3.2.2, for each district, the *n* value is the average of the number of emergency shelters within 879 m. From these figures, it is evident that there is a large difference between districts when not using a weighted coefficient (Figure 1), and when using the ratio of permanent population as a weighted coefficient (Figure 3). Without a weighted coefficient, there are districts with extremely high assessment values in the northern part, which are far from emergency shelters. With the ratio of permanent population as a weighted coefficient, there are districts with extremely high assessment values in the southeastern and southwestern parts, which have a large number of populations.

**Figure 1.** Result for Chofu City without a weighted coefficient.

**Figure 2.** Result for Chofu City with the specialization coefficient of elderly population as a weighted coefficient.

**Figure 3.** Result for Chofu City with the Ratio of permanent population as a weighted coefficient.

**Figure 4.** Result for Chofu City with the ratio of elderly population as a weighted coefficient.

**Figure 5.** Result for Chofu City with the specialization coefficient of population density as a weighted coefficient.

#### *5.2. Results for Fuchu City*

Figures 6–10 shows the results and *n* value is 6 for Fuchu City. From these figures, it is evident that there is a large difference between districts when not using a weighted coefficient (Figure 6), and when using the ratio of permanent population (Figure 8). Without a weighted coefficient, the emergency shelters concentrated in the central part that has a low assessment value. With the ratio of permanent population as a weighted coefficient, there are districts with extremely high assessment values in the southeastern part, which are far from emergency shelters. Additionally, there are districts with low assessment values in central part, which have a large number of populations and emergency shelters.

**Figure 6.** Result for Fuchu City without a weighted coefficient.

**Figure 7.** Result for Fuchu City with the specialization coefficient of elderly population as a weighted coefficient.

**Figure 8.** Result for Fuchu City with the ratio of permanent population as a weighted coefficient.

**Figure 9.** Result for Fuchu City with the ratio of elderly population as a weighted coefficient.

**Figure 10.** Result for Fuchu City with the specialization coefficient of population density as a weighted coefficient.

#### *5.3. Results for Toda City*

Figures 11–15 shows the results and *n* value is 5 for Toda City. From these figures, it is evident that there is a large difference between districts when not using a weighted coefficient (Figure 11), and when using the ratio of permanent population and the ratio of elderly population as the weighted coefficients (Figures 13 and 14). Without a weighted coefficient, emergency shelters concentrated in the districts that have low assessment values. With the ratio of permanent population and the ratio of elderly population as weighted coefficients, there are districts with extremely high assessment values in central part, which have a large number of populations and elderly populations.

**Figure 11.** Result for Toda City without a weighted coefficient was not used.

**Figure 12.** Result for Toda City with the specialization coefficient of elderly population as a weighted coefficient.

**Figure 13.** Result for Toda City with the ratio of permanent population as a weighted coefficient.

**Figure 14.** Result for Toda City with the ratio of elderly population as a weighted coefficient.

**Figure 15.** Result for Toda City with the specialization coefficient of population density as a weighted coefficient.

#### *5.4. Results for Shobara City*

Figures 16–20 shows the results and the *n* value is 1 for Shobara City. From these figures, it is evident that there is a large difference between districts when not using a weighted coefficient (Figure 16), and when using the specialization coefficient of elder population as the weighted coefficients (Figure 17). In both these two cases, districts with extremely high assessment values are distributed in the entire Shobara City, and they are far from emergency shelters. However, most of districts have low assessment values when using the specialization coefficient of population density as weighted coefficients (Figure 20).

**Figure 16.** Result for Shobara City without a weighted coefficient.

**Figure 17.** Result for Shobara City with the specialization coefficient of elderly population as a weighted coefficient.

**Figure 18.** Result for Shobara City with the ratio of permanent population as a weighted coefficient.

**Figure 19.** Result for Shobara City with the ratio of elderly population as a weighted coefficient.

**Figure 20.** Result for Shobara City with the specialization coefficient of population density as a weighted coefficient.

#### **6. Discussion**

As mentioned in Section 3.2.4, this section will compare the results between multiple areas (Chofu, Fuchu, Toda and Shobara Cities) as well as the results with and without the four types of weighted coefficients that were related to the regional characteristics based on Figures 1–20 in order to verify the validity of the assessment method. Specifically, it is necessary to compare the results between multiple areas with different regional characteristics in order to verify the validity of the assessment method. Additionally, it is also necessary to compare the results with and without weighted coefficients, focusing on the increase and decrease of the assessment values for both cases. Furthermore, by verifying the assessment method, it is capable of proving the reliability of the results for the above four cities. Based on these two types of comparisons, focusing on regional characteristics, it is also possible to respectively grasp the suitability of emergency shelters allocation in the above four cities.

#### *6.1. Comparison of the Results between Multiple Areas with Ddi*ff*erent Regional Characteristics*

In this section, focusing on the maximum and minimum assessment values that are shown in Figures 1–20, the comparison of the results between target areas will be conducted. These assessment values are based on the suitability analysis for the emergency shelters allocation in the unit of districts of each target area mentioned in Section 3.2.3.

#### 6.1.1. Results without a Weighted Coefficient

As mentioned in Section 3.2.3, when not using a weighted coefficient, analyses are conducted based on the linear distance between each district and each emergency shelter, and the scale of such facilities. Table 4 indicates the maximum (highest results) and minimum (lowest results) assessment values without a weighted coefficient. As shown in Table 2, though Shobara City has a large number of emergency shelters, because the area is extremely large, it is difficult to evenly place emergency shelters in the entire city. Therefore, there are some districts with extremely high assessment values in Shobara City. In contrast, because the areas of Chofu, Fuchu, and Toda Cities are relatively small, and the maximum assessment values in these three cities are not so high as that in Shobara City. Additionally, regarding the minimum assessment values, as all four cities have some districts with a concentration of emergency shelters, there is no noticeable difference between the four cities.

**Table 4.** Results without a weighted coefficient.


6.1.2. Results with the Specialization Coefficient of Elderly Population as a Weighted Coefficient

Table 5 indicates the results with the specialization coefficient of elderly population as a weighted coefficient. Districts where the specialization coefficient of elderly population is 0 are not considered. Regarding Shobara City, as the maximum assessment value without a weighted coefficient is significantly higher in comparison with the other three cities, the maximum assessment value remains high, regardless of the specialization coefficient of elderly population. Additionally, regarding Chofu, Fuchu, and Toda Cities, the maximum values of the specialization coefficient of elderly population are directly proportional to the maximum assessment values. However, regarding Fuchu City, as the minimum assessment value when not using a weighted coefficient is low, and the minimum value of the specialization coefficient of elderly population is also low, the minimum assessment value is lower in comparison with the other three cities. When comparing Chofu, Fuchu, and Toda Cities, which have almost the same areas as described in Table 2, when using the specialization coefficient of elderly population as a weighted coefficient, the values are directly proportional to the assessment values in these three cities.


**Table 5.** Results with the specialization coefficient of elderly population as a weighted coefficient.

#### 6.1.3. Results with the Ratio of Permanent Population as a Weighted Coefficient

Table 6 indicates the results with the ratio of permanent population as a weighted coefficient. Districts where the ratio is 0 are not considered. The maximum and minimum values of the ratio of permanent population are directly proportional to the maximum and minimum assessment values when using the ratio of permanent population as a weighted coefficient, as it is made clear in Table 6.


**Table 6.** Results with the ratio of permanent population as a weighted coefficient.

6.1.4. Results with the Ratio of Elderly Population as a Weighted Coefficient

Table 7 indicates the results with the ratio of elderly population as a weighted coefficient. Districts where the ratio of elderly population is 0 are not considered. In a similar manner as the previous section, the maximum and minimum values of the ratio of elderly population are directly proportional to the maximum and minimum assessment values when using the ratio of elderly population as the weighted coefficient.

**Table 7.** Results with the ratio of elderly population as a weighted coefficient.


6.1.5. Results with the Specialization Coefficient of Population Density as a Weighted Coefficient

Table 8 indicates the results with the specialization coefficient of population density as a weighted coefficient. Districts where the specialization coefficient of population density is 0 are not considered. The maximum values of the specialization coefficient of population density are directly proportional to the maximum assessment values in all four cities. Additionally, regarding Toda City, as the minimum value of the specialization coefficient of population density is extremely low, the minimum assessment value is also low. On the other hand, regarding Shobara City, as the minimum value of the specialization coefficient of population density is extremely high, the minimum assessment value is also high. When comparing Chofu and Fuchu Cities, as there is a small difference between each minimum value of the specialization coefficient of population density, there is likewise a small difference between each minimum assessment value.

**Table 8.** Results with the specialization coefficient of population density as a weighted coefficient.


#### 6.1.6. Summary

This section will summarize the comparison results between multiple areas, while focusing on the maximum and minimum assessment values.

(1) Maximum assessment value

Regarding Shobara City, the maximum assessment value when not using a weighted coefficient is extremely high, and all of the values with different types of weighted coefficients are also significantly higher in comparison with the other three cities. Additionally, when comparing Chofu, Fuchu, and Toda Cities, as there is a small difference between the maximum assessment values when not using a weighted coefficient, the maximum values of weighted coefficients tend to be directly proportional to the maximum assessment values. Therefore, it is possible to verify the validity of the assessment method in the present research.

(2) Minimum assessment value

For all cities, as there is a small difference between the minimum assessment values when not using a weighted coefficient, the minimum values of weighted coefficients tend to be directly proportional to the minimum assessment values. Therefore, similar to when focusing on the maximum assessment value, it is possible to verify the validity of the assessment method in the present research.

#### *6.2. Comparison of the Results with and without Weighted Coe*ffi*cients*

Based on the results of Section 5, this section will compare the results with and without the weighted coefficients, and verify the validity of the assessment method by focusing on the increase and decrease of assessment values in these two cases.

#### 6.2.1. Comparison Concerning the Specialization Coefficient of Elderly Population

Districts with high assessment values also have high assessment values when not using a weighted coefficient, when using the specialization coefficient of the ratio of elderly population as a weighted coefficient. The same trend is seen between districts with low assessment values. Additionally, it is difficult to see any significant difference between the specialization coefficient of elderly population for each district, because the value is derived by dividing the rate of the population who are over 75 years old in each district by the same rate in the entire target area.

#### 6.2.2. Comparison Concerning the Ratio of Permanent Population

When using the ratio of permanent population as a weighted coefficient, many districts with high assessment values also have high values, while many districts with low assessment values also have low values. When not using a weighted coefficient, many districts with high assessment values are far from any stations and city centres, while many districts with low assessment values are close to stations and the city centre. Therefore, when using the ratio of elderly population, which greatly differs between the districts as a weighted coefficient, the assessment values are reversed in the cases without a weighted coefficient and with the ratio of elderly population as the weighted coefficient.

#### 6.2.3. Comparison Concerning the Ratio of Elderly Population

The trend is similar to that of the previous section when using the ratio of elderly population as a weighted coefficient. Specifically, many districts with high assessment values also have high values, while many districts with low assessment values also have low values. When not using a weighted coefficient, many districts with high values are far from any stations and city centres, while the districts with low assessment values are close to stations and city centres. Therefore, when using the ratio of elderly population that greatly differs between the districts as a weighted coefficient, the assessment values are reversed in the cases without a weighted coefficient and with the ratio of elderly population as a weighted coefficient.

#### 6.2.4. Comparison Concerning the Specialization Coefficient of Population Density

When using the specialization coefficient of population density as the weighted coefficient, districts with high assessment values also have high values, while districts with low assessment values also have low values. The specialization coefficient of population density is high around stations and city centre, and it is lower as the distance from stations and city centers increases. When not using a weighted coefficient, the assessment values are high in districts that are far from any stations and city centre, while the assessment values are low in districts close to stations and city centre. Therefore, when using the specialization coefficient of population density, which greatly differs between the districts as a weighted coefficient, the assessment values are reversed in the cases without a weighted coefficient and with the ratio of elderly population as a weighted coefficient.

#### 6.2.5. Summary

This section will summarize the comparison results with and without weighted coefficients. There is no significant difference between the assessment values with and without this weighted coefficient because there is a small difference between the specialization coefficient of elderly population in each district. However, as there is a significant difference between the other 3 weighted coefficients in each district, the weighed coefficients are directly proportional to the assessment values. In other words, when not using a weighted coefficient, districts with high assessment values have low weighted coefficient values, while districts with low assessment values have high weighted coefficient values. Therefore, the assessment values are reversed, depending on whether the other three weighted coefficients are used or not. Additionally, based on the above comparison results, it is possible to verify the validity of the assessment method in the present research.

#### **7. Conclusions**

The conclusion of the present research can be summarized in the following four points.

(1) The method in the present research modifies the p-median model that derives the best facility location and conducts a suitability analysis for the emergency shelters allocation in each district. The weighted coefficients such as the specialization coefficient of elderly population, the ratio of permanent population, the ratio of elderly population, and the specialization coefficient of population density that is related to the regional characteristics are integrated into the suitability analysis for the emergency shelters allocation in Japan, using the linear distance between each district and each emergency shelter as well as the coverage of such emergency shelters.

(2) As the quantitative data related to the above 4 types of weighted coefficients, the linear distance between each district and each emergency shelter, and the facility scale of emergency shelters are used to conduct a suitability analysis, the results are also quantitative, making it a useful indicator to analyze the suitability of emergency shelters allocation. Additionally, while there are only four types of weighted coefficients that are adopted in the present research, other regional characteristics can be adopted as weighed coefficients to expand the assessment method. Furthermore, it is possible to compare the sufficiency levels of emergency shelters between districts and point out the specific districts that are short of emergency shelters, as the suitability of emergency shelters allocation are analyzed by each district. Additionally, it is also possible to visually understand the suitability of emergency shelters allocation on a small spatial scale as the unit of districts, as the results are displayed on the digital maps of GIS.

(3) In the present research, the above four types of weighted coefficients, the linear distance between each district and each emergency shelter as well as the facility scale of emergency shelters are calculated while using open data, such as the National Census and the National Land Numerical Information. As analyses in the present study are conducted based on public information, by obtaining population data and geospatial data that are similar to the present research, analyses can be conducted while using data in other areas, as well as for the past and future. Therefore, the assessment method in the present research has a high temporal reproducibility as well as spatial reproducibility. For example, by using the "future population estimate by region in Japan" of the National Institute of Population and Social Security Research [43] as future data, the shortage or overage of emergency shelters in the future can be evaluated.

(4) The present research has the limitation of data availability. Specifically, the data that are related to emergency shelters of contiguous cities could not be used, as the necessary data concerning the facility scale of emergency shelters were not available. Therefore, the assessment values of districts that were located near the administrative boundaries between the target areas and neighboring cities could not be considered to be accurate. For this reason, it is necessary to create original data related to emergency shelters.

The following are two issues that can be considered as future research topics.

(1) Calculation of the distance between each district and each emergency shelter

In the present research, as the linear distance between each district and each emergency shelter was adopted, such a distance differs from the actual road distance. Therefore, it is possible to increase the accuracy of the results by adopting the road distance. Manrique (2013) [44] demonstrated that the road distance of narrow roads that are only passable for those on foot is 1.271 times longer than the linear distance, and the road distance of wide roads that are drivable is 1.415 times longer than the linear distance. Referring to the results, it is possible to calculate road distance based on the linear distance.

(2) Application to other facilities in the time of accidents

While the present research focused on emergency shelters, the assessment method of the present research may also be applied to the suitability analysis of temporary evacuation shelters in the time of accidents.

**Author Contributions:** T.A. collected and processed the data used for variables, and conducted a suitability analysis for the emergency shelters allocation. K.Y. carried out background work, and developed the framework and process of analysis. She also initially drafted the paper. T.A. and K.Y. contributed to write up and review, and approved the paper manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
