Assessment of the peak tsunami amplitude associated with a large earthquake occurring along the southernmost Ryukyu subduction zone in the region of Taiwan

The southernmost portion of the Ryukyu Trench near the island of Taiwan potentially generates tsunamigenic earthquakes with magnitudes from 7.5 to 8.7 through shallow rupture. The fault model for this potential region dips 10^∘ northward with a rupture length of 120 km and a width of 70 km. An earthquake magnitude of M_w 8.15 is estimated by the fault geometry with an average slip of 8.25 m as a constraint on the earthquake scenario. Heterogeneous slip distributions over the rupture surface are generated by a stochastic slip model, which represents the decaying slip spectrum according to k⁻² in the wave number domain. These synthetic slip distributions are consistent with the abovementioned identical seismic conditions. The results from tsunami simulations illustrate that the propagation of tsunami waves and the peak wave heights largely vary in response to the slip distribution. Changes in the wave phase are possible as the waves propagate, even under the same seismic conditions. The tsunami energy path not only follows the bathymetry but also depends on the slip distribution. The probabilistic distributions of the peak tsunami amplitude calculated by 100 different slip patterns from 30 recording stations reveal that the uncertainty decreases with increasing distance from the tsunami source. The highest wave amplitude for 30 recording points is 7.32 m at Hualien for 100 different slips. Compared with the stochastic-slip distributions, the uniform slip distribution will be highly underestimated, especially in the near field. In general, the uniform slip assumption only represents the average phenomenon and will consequently ignore the possibility of tsunami waves. These results indicate that considering the effects of heterogeneous slip distributions is necessary for assessing tsunami hazards to provide additional information about tsunami uncertainties and facilitate a more comprehensive estimation.

1 Introduction

Almost all destructive tsunamis are generated by shallow earthquakes that occur within subduction zones. Numerous destructive tsunami events, including the M_w 9.1 Sumatra earthquake in 2004 (Lay et al., 2005), the M_w 8.8 Chile earthquake in 2010 (Lay et al., 2010; Fritz et al., 2011) and the M_w 9.0 Tohoku earthquake in 2011 (Goda et al., 2015; Goda and Song, 2016), all of which occurred in subduction zones, have occurred recently. The island of Taiwan, which is located at the convergent boundary between the Philippine Sea Plate and the Eurasian Plate, is constantly under possible threat of a tsunami. The convergence rate in this area is approximately 80–85 mm yr⁻¹ (Seno et al., 1993; Yu et al., 1997; Sella et al., 2002; Hsu et al., 2009, 2012). Thus, earthquakes occur frequently in and around Taiwan. The shallow earthquakes that occur in the Manila Trench to the south and the Ryukyu Trench to the northeast are particularly tsunamigenic, and earthquakes occur more actively in the southernmost Ryukyu Trench than in the northern Manila Trench (Wu et al., 2013). The most well-known historic tsunami events that occurred in northeastern Taiwan are the 1867 Keelung earthquake (M_w 7.0) (Tsai, 1985; Ma and Lee, 1997; Cheng et al., 2016; Yu et al., 2016) and the 1771 Yaeyama (Japan) earthquake (M_w ∼ 8) (Nakamura, 2009). Accordingly, these historic recordings demonstrate that Taiwan is under a potential tsunami threat. Furthermore, the 2011 Tohoku earthquake induced a powerful tsunami that destroyed coastal areas and caused nuclear accidents (Mimura et al., 2011). As there are four nuclear power plants along the coast of Taiwan, it is necessary to carefully estimate the tsunami hazard in addition to the hazards of compound disasters.

Probabilistic tsunami hazard analysis (PTHA) is a modification of probabilistic seismic hazard analysis (PSHA) (Cornell, 1968; SSHAC, 1997), and it is intended to forecast the probability of tsunami hazards for a given region as comprehensively as possible. The recurrence rates of earthquakes have typically been estimated using the Gutenberg–Richter relationship (Gutenberg and Richter, 1944) for a defined source region in consideration of tsunamis triggered by earthquakes. The assessment of the wave height is one of the primary differences between PTHA and PSHA. PSHA assesses the ground motion based on empirical attenuation relationships (Wang et al., 2016), while PTHA assesses tsunami wave heights using empirical approaches or tsunami simulations (Geist, 2002; Geist and Parsons, 2006, 2009). Geist and Parsons (2006) mentioned that the tsunami wave height follows a definable frequency–size distribution over a sufficiently long period of time within a given coastal region (Soloviev, 1969; Houston et al., 1977; Horikawa and Shuto, 1983; Burroughs and Tebbens, 2005). This method is of great use in establishing the tsunami probability for a region if there is an extensive catalog of observed tsunami wave heights. However, given the wide distribution of global tsunamigenic earthquakes within seafloor regions throughout subduction zones, the tsunami records obtained from coastal gauges or/and ocean buoys are too sparse to comprehensively assess the associated hazards and the recording time, since their deployment is too short to enable a study on the recurrence intervals of tsunamis and earthquakes. Consequently, because the existing tsunami catalog is limited, simulations represent an effective approach. Conventional tsunami simulation adopts a simple source approximation and applies elastic dislocation theory to calculate the deformation of the seafloor surface assuming a uniform slip over the entire fault surface (Okada, 1985; Okal, 1982). However, the complexities of earthquake rupture processes play a substantial role in the generation of tsunamis. Conventional approaches are therefore unable to capture various features of short-wavelength tsunamis in the near field (Geist, 2002; Geist and Parsons, 2009). The results of previous studies that simulated tsunamis originating from historical earthquakes around Taiwan (Ma and Lee, 1997; Wu et al., 2008) using uniform slip models agreed only with long-wavelength observations. For the purposes of hazard mitigation, it is critical to predict the amplitudes of tsunamis along various coastlines for a given earthquake as accurately as possible. To make such predictions, the effects of the rupture complexity must be taken into consideration. Recent developments in PTHA have included the adoption of stochastic slip distributions of earthquakes to determine the overall probability of particular tsunami heights (Geist and Parsons, 2006, 2009). The adoption of stochastic slip distributions is able to quantify the variations in reasonable evaluations of the probabilities of specified tsunami heights at individual locations resulting from a specific fault.

In this study, we assess the heights of tsunamis along the coastline of Taiwan generated by the potential tsunamigenic zone at the southernmost end of the Ryukyu subduction zone. This potential zone is located close to Taiwan, and at least 10 earthquakes (M_w > 7) have occurred over the past 100 years (Hsu et al., 2012), the largest of which was the M_w 7.7 in 1920 (Theunissen et al., 2010). For this area, the plausible magnitude of the strongest earthquake is determined within the range between 7.5 and 8.7 (M_w) (Hsu et al., 2012). The fault zone is bounded by the Longitudinal Valley Fault to the west and the Gagua Ridge to the east (Hsu et al., 2012). This fault geometry with a defined rupture length and width is employed herein, and an earthquake with a magnitude of 8.15 is used in the tsunami simulation. The stochastic slip model is invoked to describe the uncertainty in the rupture pattern over the fault plane to enable a more realistic assessment of the tsunami probability.

2 Earthquake scenario and tsunami simulation

2.1 Assessment of seismic parameters

The estimated maximum magnitude of a possible earthquake scenario is essential for establishing the fundamental seismic conditions of the tsunami simulation. The scenario of a potential rupture fault extending to a depth of 13 km proposed by Hsu et al. (2012) occurs along the southernmost Ryukyu trench with a rupture length of 120 km, a width of 70 km and a dip of 10^∘. Kanamori and Anderson (1975) investigated the relation between the rupture area and moment, and revealed that most of the average stress drops (Δσ) vary between 10 and 100 bars. The average stress drop for most interplate earthquakes is approximately 30 bars, and thus, we set an average stress drop of 30 bars. The stress drop and seismic moment (M₀) relation along a dip-slip fault is described as follows (Kanamori and Anderson, 1975):

$$\begin{array}{}\text{(1)}& M_{\mathrm{0}}={\displaystyle \frac{\mathit{\pi }\left(\mathit{\lambda }+\mathrm{2}\mathit{\mu }\right)}{\mathrm{4}\left(\mathit{\lambda }+\mathit{\mu }\right)}}\mathrm{\Delta }\mathit{\sigma }{W}^{\mathrm{2}}L,\end{array}$$

where W and L are the width and length of the rupture plane. We can obtain the moment for this scenario under an average stress drop of 30 bars with the assumed rupture geometry. In Eq. (1), μ denotes the rigidity and λ is the Lamè parameter. We assume that the crust is elastic and homogeneous. Hence, $\mathit{\mu }=\mathit{\lambda }=\mathrm{30}$ GPa (Fowler, 2004; Piombo et al., 2007). Additionally, the seismic moment can be represented by the rupture area and average slip as follows (Lay and Wallace, 1995):

$$\begin{array}{}\text{(2)}& M_{\mathrm{0}}=\mathit{\mu }A\stackrel{\mathrm{‾}}{D}.\end{array}$$

Moreover, the seismic moment is dependent on the rupture area (A) and average slip ($\stackrel{\mathrm{‾}}{D})$; thus, the average slip can be estimated by Eq. (2), and it is calculated to be 8.25 m. Then, the seismic moment can be transformed into the magnitude M_w by the following (Hanks and Kanamori, 1979):

$$\begin{array}{}\text{(3)}& M_{\mathrm{w}}=\left({\displaystyle \frac{\log M_{\mathrm{0}}}{\mathrm{1.5}}}\right)-\mathrm{10.73}.\end{array}$$

Therefore, the maximum possible earthquake magnitude is M_w 8.15 ($M_{\mathrm{0}}=\mathrm{2.07}\times {\mathrm{10}}^{\mathrm{28}}$ dyne-cm).

2.2 Stochastic slip model

The rupture process of an earthquake is extremely complex. Seismic inversion results reveal that the slip distribution of a rupture has a heterogeneous spatio-temporal development. Consequently, using a simplified uniform slip distribution to simulate a tsunami captures only the long-wavelength portion of the tsunami field (Geist and Dmowska, 1999). In addition, the temporal description of the seismic rupture process can be ignored because the propagation velocity of the tsunami wave is substantially slower than the seismic rupture velocity (Dean and Dalrymple, 1991; Ma et al., 1991; Wang and Liu, 2006). Andrews (1980) showed that the static slip distribution is directly related to stress changes and that the spectrum of the slip distribution is proportional to k⁻² decay in the wave number domain:

$$\begin{array}{}\text{(4)}& \left|F_{s,t}\left[D_{x,y}\right]\right|\propto k^{-\mathrm{2}}\end{array}$$

where D_x,y is the slip distribution over a 2-D lattice, F_s,t is the 2-D Fourier transform, and $k=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}$ is the radial wave number. The k⁻² power law indicates that the slip distribution has self-similar characteristics; moreover, this characteristic can also be demonstrated from a fractal perspective (Tsai, 1997). Based on self-similarity, Herrero and Bernard (1994) introduced the k-square model, which leads to the ω-square model (Aki, 1967). The slip spectrum follows k⁻² decay beyond the corner radial wave number (k_c), which is proportional to 1/L_c. The L_c depends on the characteristic rupture dimension (Geist, 2002).

The heterogeneous slip distribution is proportional to k⁻² and is similar to fractional Brownian motion as a stochastic process (Tsai, 1997). The stochastic slip distribution can be described by multiplication in the Fourier domain:

$$\begin{array}{}\text{(5)}& D_{x,y}\propto F_{x,y}^{-\mathrm{1}}\left[F_{s,t}\left[X_{x,y}\right]\times k^{-\mathrm{2}}\right],\end{array}$$

where X_x,y is a random variable for the spatial distribution that randomizes the phase, and $F_{x,y}^{-\mathrm{1}}$ is the inverse 2-D Fourier transform. The random distribution of X, which is best described by a non-Gaussian distribution, especially by a Lèvy distribution, can be calculated by reversing Eq. (5) (Lavallée and Archuleta, 2003; Lavallée et al., 2006). The Lèvy distribution can be described by four parameters, namely α, β, γ and μ_L, as follows:

$$\begin{array}{ll}{\displaystyle }\mathit{\phi }& {\displaystyle }\left(t\right)=\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\exp \left(-{\mathit{\gamma }}^{\mathit{\alpha }}{\left|t\right|}^{\mathit{\alpha }}\left[\mathrm{1}+i\mathit{\beta }\mathrm{sign}\left(t\right)\tan {\displaystyle \frac{\mathit{\pi }\mathit{\alpha }}{\mathrm{2}}}\right.\right.& \\ \left.\left.\left({\left|\mathit{\gamma }t\right|}^{\mathrm{1}-\mathit{\alpha }}-\mathrm{1}\right)\right]+i{\mathit{\mu }}_{L}t\right),& \mathit{\alpha }\ne \mathrm{1}\\ \exp \left(-\mathit{\gamma }\left|t\right|\left[\mathrm{1}+i\mathit{\beta }{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\mathrm{sign}\left(t\right)\right.\right.& \\ \left.\left.\left(\ln \left|t\right|+\ln \mathit{\gamma }\right)\right]+i{\mathit{\mu }}_{L}t\right),& \mathit{\alpha }=\mathrm{1}\end{array}.\right.\end{array}$$

The parameter α, 0 < α≤2, affects the falloff rate of the probability density function (PDF) for the tail. The parameter β, $-\mathrm{1}\le \mathit{\beta }\le \mathrm{1}$ controls the skewness of the PDF, and the parameter γ, γ > 0 controls the width of the PDF. The parameter μ_L, −∞ < μ_L < ∞, is related to the location of the PDF. The Lèvy distribution is effective at describing the distribution of a random variable, i.e., X, from real earthquake events, implying that the slip distribution without self-similarity has a heavy tail behavior (Lavallée et al., 2006). Based on experiments of generating stochastic slip distributions, this heavy tail behavior affects the intensity of an extreme value (Lavallée and Archuleta, 2003).

Figure 1 (a) The spatially random variable is a truncated Lèvy distribution. The Lèvy parameters obtained from the Northridge earthquake were taken from Lavallée et al. (2006). (b) A stochastic slip distribution is generated by filtering the spatial random variable X in Fig. 1a. This slip pattern produces the highest maximum wave amplitude at Hualien station. (c) The slip spectrum is calculated from Fig. 1b. This slip spectrum decays with an exponent of −2 according to a characteristic corner radial wave number. This verifies that the synthetic slip distribution is identical to the k-square model and the condition of the rupture dimension. (d) This stochastic slip distribution produces the lowest maximum wave amplitude at Hualien station.

The stochastic slip distribution is generated by a 2-D spatially random distribution by imposing a self-similar characteristic beyond the corner radial wave number, which is constrained by the rupture dimension, in the wave number domain. In this study, the potential rupture fault is divided into 5 × 5 km² subfaults. The grid is composed of 24×14 meshes along the strike and dip directions. The produced variable with a spatially random distribution adopts the Lèvy distribution (α=1.51, β=0.2, γ=28.3, ${\mathit{\mu }}_L=-\mathrm{0.9}$), which is the dip-slip result from Lavallée et al. (2006) shown in Fig. 1a. In Lavallée et al. (2006), the slip distribution of the Northridge earthquake was divided into the dip-slip and strike-slip directions, and they were calculated by an inverse 2-D stochastic model to obtain the values of the Lèvy PDF. The values of the Lèvy PDF, which are mentioned above, are indicative of the result of the dip-slip direction. The Northridge earthquake is a thrust earthquake (Davis, 1994), and thus, it has a faulting mechanism that is approximately similar to our scenario fault model. There are no inverted slip models of past earthquakes in the study area that can be used to conduct an analysis of the Levy PDF parameters; therefore, the Lèvy distribution in Lavallée et al. (2006) is adopted in this study. From the perspective of mathematical operations, the slip distribution in Eq. (5) represents a filtered random distribution. However, for consistency with the physical behavior over the rupture surface suggested by the results of the inverse modeling, truncation of the Lèvy distribution must be performed to constrain the extreme slip value. The synthetic slip distribution (Fig. 1b) produced by the spatially random distribution in Fig. 1a is heterogeneous, and its power spectrum obeys a k-square model at high wave numbers (Fig. 1c). The average slip of this synthetic slip distribution is 8.25 m, indicating that the earthquake energy is constant as estimated above, and the maximum slip is 31.02 m. One hundred different slip distributions are produced for the tsunami simulation representing the uncertainty in the results associated with complex rupture processes. In the 100 sets of results, the maximum slip range is between 20.17 and 37.97 m. Smooth processes are not included, nor are additional regional constraints for the slip distribution. There are two reasons for this application. The first is that we do not have information regarding where the plate interface is locked or the locations of asperities often repeat in historical events. The second is that some studies reported that the asperities extend to the boundary of the fault model (Ide et al., 2011; Lay et al., 2011; Shao et al., 2011; Yue and Lay, 2011). According to these reasons, we do not apply any additional constraints for stochastic slip distributions. Similarly, the uniform slip case constitutes a complete uniform slip distribution over the whole fault plane. Figure 1b and d demonstrate the stochastic distribution of the scenario source models causing the maximum and minimum wave heights at recording station 26 (Hualien) (Fig. 2). Both patterns affecting the propagation will be discussed in Sect. 3.1.

Figure 2 The map of Taiwan shows the fault model and recording stations used in this study. The bathymetry is divided into two layers with different resolutions. The resolution of the outer layer is 4 min, and the resolution of the inner layer is 1 min. The red grid denotes the potential fault model (5×5 km² grid size). The pins represent 30 tidal gauges of the CWB. The red and blue colors indicate stations on the eastern and western sides of Taiwan, respectively, and the yellow squares represent the sites of nuclear power plants.

2.3 Numerical tsunami simulation

Figure 2 shows the computational domain, recording stations and fault model. The potential rupture fault is divided into 5 × 5 km² subfaults, and the stochastic slip distribution model is applied to determine the amount of discrete slip on each subfault. Vertical seafloor displacements caused by slip along the rupture plane are calculated using elastic dislocation theory (Okada, 1985). The Cornell Multi-grid Coupled Tsunami model (COMCOT) is used to run the tsunami simulations. COMCOT is capable of efficiently studying the entire life span of a tsunami, including its generation, propagation, runup and inundation (Wang, 2009), and it has been widely used in studying many historical tsunami events, such as the 1960 Chilean tsunami (Liu et al., 1995), 1992 Flores Island tsunami (Liu et al., 1995), 2003 Algeria tsunami (Wang and Liu, 2005), 2004 Indian Ocean tsunami (Wang and Liu, 2006, 2007) and 2006 Pingtung tsunami, Taiwan (Wu et al., 2008; Chen et al., 2008). COMCOT solves linear or nonlinear shallow water equations for spherical or Cartesian coordinates using the finite difference method. With a flexible nested grid system, it can properly guarantee both the efficiency and the accuracy from the near-coastal region to the far-field region. Two grid layers are used to simulate the propagation of tsunamis. The Manning coefficient is 0.013 in this study, which assumes a sandy sea bottom (Wu et al., 2008). The bathymetry-adopted open data from the National Oceanic and Atmospheric Administration (NOAA) that can be download from https://maps.ngdc.noaa.gov/viewers/wcs-client/, last access: 26 March 2018 (Amante and Eakins, 2009). The resolution of the outer layer is 4 min for the solution of the linear shallow water equation, and the resolution of the inner layer is 1 min for the solution of the nonlinear form of the shallow water equation. There are 30 recording stations referring to the positions of tidal gauges maintained by the Central Weather Bureau (CWB) along the coastlines of Taiwan and the outlying islands. The CWB website presents the locations of the tide stations (http://e-service.cwb.gov.tw/HistoryDataQuery/index.jsp, last access: 26 March 2018 and https://www.cwb.gov.tw/V7e/climate/marine_stat/tide.htm, last access: 26 March 2018). These locations are shifted slightly to the grid nodes to accurately record the data. Table 1 presents the locations and water depths of the recording stations in the computational mesh.

Table 1 The maximum, minimum and average wave heights with their standard deviations for the PTA probability distributions (in meters) with the maximum wave heights from the uniform slip model. The water depths at the stations in the computational mesh are also included.

3 The effect of heterogeneous slip on tsunamis

The stochastic slip model produces different slip distributions with the same fault geometry in addition to a constant average slip and a constant seismic moment. The model is used to describe the heterogeneous slip pattern of an earthquake and to further examine its effect on the tsunamis originating from the southernmost end of the Ryukyu subduction zone adjacent to Taiwan. According to the previous sections, the maximum possible earthquake magnitude is determined to be M_w 8.15 with an average slip of 8.25 m. Furthermore, the uniform slip distribution on the rupture plane is also used to simulate a tsunami to facilitate a discussion of the difference between the effects of uniform and heterogeneous slip on tsunamis.

3.1 Initial water elevation and energy propagation

The static vertical displacement of the ocean floor is modeled using elastic dislocation theory (Okada, 1985) with a static slip distribution. The vertical seafloor displacement is modeled as the initial water level, and the horizontal component of the seabed displacement is not included in the simulation. Figure 3a shows the initial water elevations produced by a uniform slip distribution, and Fig. 3b exhibits the maximum free-surface elevation during the propagation. Figure 3c and e demonstrate the initial water elevations produced by the stochastic slip distributions (Fig. 1b and d). The initial water elevation with a uniform slip distribution is simple and smooth, but those with stochastic slip models are more complex and relatively heterogeneous. Nonuniform slip causes an apparent change in the wavelength distribution of the initial free-surface elevation (i.e., the potential energy distribution), which affects the path of energy propagation. In the uniform slip scenario, the maximum free-surface elevation pattern is straightforward and clearly controlled by the topography. However, many strong and seemingly chaotic paths of wave energy appear in the nonuniform slip scenarios, and the free-surface field exhibits additional uncertainties in terms of the flow. In Fig. 3b, the maximum free-surface elevation mainly propagates toward two places where the seafloor bathymetry becomes shallow relative to the deep areas northeast of Taiwan, as shown in Fig. 2. Although the propagation paths due to the nonuniform slip distributions (Fig. 3d and f) also have the same characteristics, it is notable that the paths followed by the wave energy differ depending on the rupture pattern. To the northeast of Taiwan in Fig. 3f, there is a strong wave path connecting the two higher-elevation areas of bathymetry. However, this behavior is not observed in Fig. 3b and d. In addition, the maximum elevation on the footwall in Fig. 3d is higher than that in Fig. 3f. In Fig. 3b, the high elevation appears only along the coast on the footwall side. These results indicate that the wave energy variation depends on the rupture pattern, thereby causing differences in the wave paths and leading to completely different tsunami amplitudes.

Figure 3 Panels (a), (c) and (e) are the initial water elevations, and the color bars represent the elevation of the initial water surface. Panels (b), (d) and (f) are the maximum free-surface elevation (i.e., the distribution of the energy path), and the color bars represent the elevation of the maximum free-surface. Panels (a) and (b) display the results with a uniform slip distribution. Panels (c) and (d) display the results from Fig. 1b. Panels (e) and (f) display the results from Fig. 1d. The seafloor elevation fundamentally dominates the tsunami propagation, but the slip distribution also has a strong influence. In (a, c and e), yellow squares represent NPPs; in (b, d and f), the NPPs are represented by open squares.

3.2 Wave characteristics

Thirty stations located along the coastlines are available for recording the amplitude of the tsunami wave height. Relative to the other stations, stations 25 (Shihti), 26 (Hualien) and 27 (Suao) are situated near the potential rupture fault, and they have high wave amplitudes and enormous variations in the tsunami simulations of 100 different slip distributions; consequently, the time series of the wave heights at these stations are shown as an example (Fig. 4). The time series of the wave heights at the other stations are shown in the Supplement. The variability in the distribution of the initial free-surface elevation results in substantial phase changes and different wave heights. It is worth noting that the average of the disordered and chaotic time series produced by the 100 different slip distributions is almost identical to the results of the time series produced by the uniform case. This implies that the uniform slip distribution simply represents an average result and that it cannot represent all of the possible situations.

According to the statistical results from 100 different slip patterns (Table 1) for 30 stations, Hualien station has a maximum wave amplitude of 7.32 m, and its maximum wave amplitude interval ranges from 1.87 to 7.32 m, which constitutes the widest interval for any recording site, and the standard deviation of this distribution is 1.024 m. These findings indicate that Hualien station has a high uncertainty in this scenario. However, the maximum wave amplitudes from the uniform slip distribution are relatively lower than those from the stochastic results. Following the above findings, we need to consider whether the estimations from the uniform slip case are appropriate for hazard analysis by focusing on the maximum wave amplitude issue.

Figure 4 The time series of the wave heights recorded at stations 25 (Shihti), 26 (Hualien) and 27 (Suao). Gray lines represent the time series of 100 different slip distributions, black lines represent the averages of the gray lines, blue lines represent the 95 % confidence intervals, and red lines are the time series produced using uniform slip distribution. Parts of the wave heights at station 27 are lower than the water depths, and thus, these curves have been truncated.

3.3 The peak tsunami amplitude probability

According to the results of our simulations, we calculate the probability of the peak tsunami amplitude (PTA) at each recording station as shown in the histogram of Fig. 5. To verify the representativeness of the PTA probability distributions, another 100 sets of different slip distributions are produced and simulated under the same seismic conditions. In Fig. 5, the shapes of the PTA distributions from another 100 sets (black lines) are similar to the shapes of the histograms from the first 100 sets. These results verify the representativeness of the PTA probability distributions produced from 100 sets of slip distributions. This test also reinforces the reproducibility of our simulations and demonstrates that the number of simulations is roughly satisfactory for statistical analysis. Of course, the more slip distributions we use, the more comprehensive and stable the range we obtain.

Figure 5 The probabilities of the PTA along the coast of Taiwan (blue: stations 1 ∼ 19, red: stations 20 ∼ 30). The histograms display the PTAs derived from 100 different slip simulations. The black lines represent the results from another 100 simulations, and the orange lines represent the PTA obtained using a uniform slip distribution. The PTA probability distribution gives a clear PTA range and its occurring probability. The map of Taiwan shows the station locations and the sites of four NPPs (yellow squares).

In Fig. 5, the PTA distributions for the stations in eastern Taiwan (red markers) have much higher values than those in western Taiwan (blue markers) due to the specified location of the tsunami source. The shapes of the PTA distributions in eastern Taiwan resemble lognormal distributions, while those in western Taiwan resemble normal distributions. We suppose that the attenuation of the wave propagation causes the lognormal distributions to degenerate into normal distributions. The PTAs produced by a uniform slip distribution are generally located in the middle of the PTA distributions. Both PTA values (i.e., the value of the PTA from the uniform slip distribution and those from the stochastic slip distribution models) decrease with distance from the potential fault due to the attenuation of the wave propagation (Fig. 5 shows the results for all stations, and Fig. 6 shows the results for stations 20 through 30 in eastern Taiwan). However, some stations, e.g., stations 17, 19 and 21, do not precisely follow this trend; this could be the result of the coastal topography and the presence of an energy channel. From Fig. 3d, in comparison with the adjacent coastline, station 21 is located exactly where the wave energy gathers. In addition, broad distributions are frequently observed at promontories along the coastline and are caused by complex propagation path effects between the source region and the recording locations (Geist, 2002). There are many compound factors that affect the tsunami propagation and maximum wave height. Figure 6 presents the relation between the distance and wave height and shows the PTA distributions following Fig. 5. The x axis presents the shortest distance between the stations and fault plane. On the footwall side, stations 20 and 22, which do not directly face the energy propagation path (Fig. 3f), are located on islands off the coast of Taiwan; consequently, their PTA distributions are lower than those of stations 21 and 23, even though the distances from the potential fault are similar. On the hanging wall, station 29 is farther from the coastline of Taiwan than other stations; however, because of the real location of the station and its numerical grid setting, its PTA distribution is lower than that of station 30 (Fig. 3b). The ranges of the PTA distributions converge with increasing distance on both sides of the fault. Moreover, the PTA distributions and their average values roughly exhibit a linear decrease with increasing distance except for stations 25 and 26. In contrast, these two stations in the near field are directly affected by initial water elevation, and thus, the PTAs caused by uniform slip are quite low.

Although the seismic parameters have already been defined as constants in our experiment, there exists an uncertainty in the PTA, which is not a constant value. Hence, the uniform case cannot provide this uncertainty, and thus, the PTA could be underestimated. The results give specific PTA ranges, which represent the wave height uncertainties for the scenario of earthquakes originating from the Ryukyu Trench. It is therefore necessary to consider the effects of a heterogeneous slip distribution to comprehensively assess the tsunami hazard.

Figure 6 The relation between the distance and wave height for stations 20 to 30 in eastern Taiwan. Panel (a) is the station on the footwall side. Stations 20 and 22 (blue color) are off the shoreline of Taiwan. Panel (b) represents the stations on the hanging wall side. Both sides roughly exhibit a linear decay and range of uncertainty converging with increasing distance for the tsunami amplitude. Red bars show the PTA of the uniform slip distribution, and yellow bars show the average of the PTAs from the stochastic slip models.

4 Discussion

4.1 Tsunami

Most coastlines threatened by near-field tsunamis, such as the coasts of Chile, Japan and Indonesia, are parallel to the trench axis of the associated subduction zones. Many tsunami events, including the M_w 8.8 Chile earthquake in 2010 (Lay et al., 2010; Fritz et al., 2011), the M_w 9.0 Tohoku earthquake in 2011 (Goda et al., 2015; Goda and Song, 2016), the M_w 9.1 Sumatra earthquake in 2004 (Lay et al., 2005), and the M_w 8.1 Mentawai earthquake in 2010 (Satake et al., 2013), have occurred along these regions. However, the potential rupture fault in this study along the southernmost Ryukyu subduction zone is perpendicular to the coast of Taiwan, which directly affects the first wave motion. The first motion on the footwall is up; conversely the first motion on the hanging wall is down. As a result, the coastline retreats from the land to the sea as the first tsunami wave approaches, allowing people additional time to leave the seafront.

The effect of a heterogeneous slip distribution is important and necessary for near-field estimations (Geist, 2002 and Ruiz et al., 2015). Figure 5 shows that the PTA distributions in the near field are broad, and they narrow with increasing distance from the potential fault. The uncertainty in the near field is higher than in the far field. At most of the eastern stations, the values of the average PTA approach uniform results, but the uniform slip results at stations 25 and 26 are close to the minimum PTA (Table 1). Geist (2002) presented the average and extreme nearshore PTA calculated for 100 different slip distributions and compared them with the uniform slip results (Fig. 6a in Geist, 2002). The range of the PTA also becomes narrower with increasing distance. The values from the uniform slip distribution and the average PTA are similar, but some of the average values are close to the minimum PTA between approximately 19 and 19.5^∘ N. Similar characteristics of the average PTA and the results from the uniform case are observed in different regions. The average PTA is equal to the uniform slip result in the nearshore region, but this could be caused by other factors (e.g., distance to the tsunami source, propagation path, initial water elevation) that shift the average PTA toward the minimum PTA.

Four nuclear power plants (NPPs) are located on the island of Taiwan. According to the numerical results, we infer that the mean PTA in the coastal area of NPP4 ranges from approximately 2 to 3 m. The distribution at this plant may be wider than those at other nuclear power plants due to its position relative to the tsunami source. Moreover, NPP4 is located on the shore of a bay with a curved shape; the magnification effect from the geometrical shape of the bay may serve to enhance the PTA therein. NPP3 also exhibits this condition insomuch that the energy is concentrated at the location of the plant (Fig. 3b, d and f). For the coastal areas around NPP1 and NPP2, the PTA distributions are between 1 and 2 m. The coastlines of these two nuclear power plants slightly face the direction of tsunami propagation, and thus, their PTAs should be higher than those along adjacent coastlines (Fig. 3b, d and f). In general, under this scenario, the coastline at NPP4 has the largest threat. Although NPP3 is far from the tsunami source, it faces a wave height of approximately 1.5 m on average with a ±0.5 m range of uncertainty. However, NPP3 is closer to the Manila subduction zone, and thus, it could be threatened by a tsunami originating from the Manila Trench. In contrast, the coastlines of NPP1 and NPP2 are relatively safe and have fewer uncertainties with regard to the PTA.

The use of heterogeneous slip patterns clearly delineates the range of possible waveforms and provides more information on latent uncertainties in the wave height. The 95 % confidence intervals for the wave height from 100 sets in each time series provide us with a specific range for the amplitude of the tsunami wave (Fig. 4). According to these time series, we are aware of the periods of tsunami runup and runoff and can prepare supporting policies to reduce associated disasters. For example, a nuclear power plant uses a trench from the ocean to intake water to cool the reactor; thus, if the sea level is too low to take in water, the temperature of the reactor will rise excessively, causing a nuclear disaster. Based on the results of simulations, we can estimate how much water should be stored for tsunami runoff. This issue requires more attention in Taiwan because four nuclear power plants are located near the coast.

4.2 Stochastic slip model

The results of the tsunami simulations illustrate that the effect of the slip distribution on the rupture plane has significant effects on the wave propagation and wave height. The correctness of this slip distribution determines whether the wave height calculations represent a useful reference. However, some parameters of the stochastic models could influence the synthetic slip distributions. For instance, the exponent of the slip spectrum is associated with the roughness of the slip distribution. Higher exponential values inhibit the powers of high wave numbers, leading to smoother slip distributions; conversely, lower values lead to rougher slip distributions. In general, the k-square model needs to be followed. Furthermore, the interpolation of the slip distribution for a given geometry will affect the exponent of k (Tsai, 1997). Interpolation will smooth the original pattern. The powers of short wave numbers will be depressed and the powers of long wave numbers will be enhanced. Moreover, the random spatial variability of the slip distribution is relatively critical. According to Lavallée and Archuleta (2003) and Lavallée et al. (2006), we adopt the truncated non-Gaussian distribution for the spatial variability. This truncation limits the non-Gaussian distribution to a particular range. However, extreme truncation will cause the heavy-tailed characteristic of this distribution to become less pronounced or even disappear, similarly to a Gaussian distribution. In mathematics, the synthetic slip distribution is a filtering process insomuch that the characteristics of a heavy-tailed distribution affect the extrema of the slip distribution. The maximum slip will be greater as the truncated range increases, and the maximum slip may exceed reasonable values if the truncated range is excessively wide. Therefore, the parameters must be chosen carefully to match the observations acquired by inversion.

5 Conclusions

The maximum possible earthquake magnitude is M_w 8.15 with an average slip of 8.25 m in the southernmost portion of the Ryukyu Trench. One hundred slip distributions of the seismic rupture surface were generated by a stochastic slip model. The maximum slip range is between 20.17 and 37.97 m, and the average slip of each model is consistent with 8.25 m. A heterogeneous slip distribution induces variability in the tsunami wave heights and the associated paths of propagation. The simulated results demonstrate that the complexity of the rupture plane has a significant influence on the near field for local tsunamis. The PTA distribution provides a specific range for the wave height and the probability of occurrence in this scenario. These distributions and their average values exhibit an approximately linear decrease with increasing distance. The coastline, which is situated very close to or even atop the tsunami source, is directly affected by the rupture slip distribution. Then, the range of the PTA distribution will converge with increasing distance from the tsunami source. In this study, Hualien station, which is located directly above the source, has the widest PTA interval (1.87–7.32 m) and the highest wave amplitude. The statistical summary reveals that this station has a standard deviation of 1.63 m, which is larger than those of the other stations, and the largest uncertainty. However, the PTA caused by the uniform slip distribution is only 1.63 m, which is much lower and even below the average (3.36 m) at this station. This finding implies that a simplified earthquake source cannot completely represent the tsunami amplitudes in reality. If we adopt a uniform slip distribution to assess tsunami hazards, those hazards will be critically underestimated. Furthermore, the variances of tsunami amplitudes are imperative for assessing tsunami hazards, and the quantitative technique employed is also important.