The role of diffraction effects in extreme run-up inundation at Okushiri Island due to 1993 tsunami

The tsunami generated on 12 July 1993 by the Hokkaido–Nansei–Oki earthquake (Mw= 7.8) brought about a maximum wave run-up of 31.7 m, the highest recorded in Japan during the 20th century, near the Monai Valley on the west coast of Okushiri Island (Hokkaido Tsunami Survey Group, 1993). To reproduce the extreme run-up height, the three-dimensional non-hydrostatic model (Flow Science, 2012), referred to here as the NH-model, has been locally applied with open boundary conditions supplied in an offline manner by the three-dimensional hydrostatic model (Ribeiro et al., 2011), referred to here as the H-model. The area of the H-model is sufficiently large to cover the entire fault region with one-way nested multiple domains. For the initial water deformation, Okada’s fault model (1985) using the sub-fault parameters is applied. Three NH-model experiments have been performed, namely without islands, with one island and with two islands. The experiments with one island and with two islands give rise to values close to the observation with maximum runup heights of about 32.3 and 30.8 m, respectively, while the experiment without islands gives rise to about 25.2 m. The diffraction of the tsunami wave primarily by Muen Island, located in the south, and the southward topographic guiding of the tsunami run-up at the coast are, as in the laboratory simulation (Yoneyama et al., 2002), found to result in the extreme run-up height near Monai Valley. The presence of Hira Island enhances the diffraction of tsunami waves but its contribution to the extreme run-up height is marginal.


Introduction
The tsunami generated on 12 July 1993 by the Hokkaido-Nansei-Oki earthquake (M w = 7.8) produced the worst local tsunami-related death toll in Japan in 50 years, bringing about the maximum wave run-up of 31.7 m near the Monai Valley on the west coast of Okushiri Island (Hokkaido Tsunami Survey Group, 1993;Shuto and Matsutomi, 1995).The location of the study area is shown in Fig. 1.The runup value was the largest recorded in 20th century Japan and was among the highest among tsunamis of non-landslide origin in the world before the occurrence of the 51 m run-up height at Leupung Beach twin peaks in Lhoknga by the 2004 Sumatra-Andaman earthquake.The 1993 tsunami was also recorded in Korea (Choi et al., 1994(Choi et al., , 2003)).
Several works were done in the past on run-up heights following the 1993 earthquake.Abe (1995) estimated the localmean run-up height in segment intervals of about 40 km using an empirical method for the near-field tsunami warning developed by Abe (1989Abe ( , 1993)).Takahashi et al. (1995) applied an inversion method as well as a model based on shallow water equations along with the comparison of various source models.Myers and Baptista (1995) applied the depthintegrated finite element model originally developed by Luettich et al. (1991) with use of Okada's fault model (1985).All these works, however, fail to reproduce the extreme runup height at Okushiri Island even though overall agreements were somehow reasonable.A maximum run-up height rea- sonably close to the observation was obtained by Titov and Synolakis (1997) on the basis of the depth-averaged shallow water model using a variable grid system about 50 m from the shoreline.It is, however, noted that their model was frictionless.
To obtain a better understanding of the mechanism of local occurrence of extreme run-up at Okushiri Island, a laboratory model of 1/400 scale closely resembling the actual bathymetry and topography of Monai Valley was constructed in a 205 m long, 6 m deep and 3.5 m wide tank at the Central Research Institute for Electric Power Industry (CRIEPI), Abiko Research Laboratory, Japan.Furthermore, a free surface analysis code developed in CRIEPI was applied to the run-up in the Monai district (Yoneyama et al., 2002).Takahashi (1996) commented that the code application was intended to simulate the laboratory measurements rather than the run-up of the real event.Kakinuma (2005) (Flow Science, 2012) with high resolution to reproduce the extreme run-up height near Monai Valley.However, the application of such a sophisticated high-resolution model with use of a radiation condition for sea regions significantly larger than the whole source region is hardly possible.We have therefore constructed a modeling system which is composed of two models (one is a local model, the other is a regional model) nested in a one-way offline manner.It is addressed such that the hydrostatic model in a series of regional domains is used for the propagation of tsunami waves and the non-hydrostatic model is used for the reproduction of extreme run-up heights on the local domain with high resolution.The local model covering the area near Monai Valley is based on the non-hydrostatic z coordinate model that solves fully nonlinear dispersive Navier-Stokes equations by the finite difference method (Flow Science, 2012).Specifically, a true volume-of-fluid method was used to compute free surface motion and the fractional area/volume obstacle representation technique to model complex geometric regions.To describe turbulence in a viscous fluid, the updated k model with the help of renormalization group analysis was used.The regional model covering the entire fault region is based on the three-dimensional hydrostatic σ coordinate model using finite volumes for a spatial integration and a semi-implicit algorithm for integration in time (Ribeiro et al., 2011).For computing inundation in land, the moving boundary scheme described by Martins et al. (2001) is used.In essence, a specific cell is considered a dry cell when the instantaneous total depth at the computing point of sea surface elevation and the value of neighboring sea surface elevation plus the local depth of the considered cell are smaller than a prescribed minimum depth.In the study the minimum depth was set to 0.1 m.Three one-way domains of different grid sizes dynamically nested were used in the regional model to provide the local model with open boundary conditions in an offline manner.The three-dimensional non-hydrostatic and hydrostatic models are hereafter referred to as the NH-model and H-model, respectively.The initial water deformation is computed by Okada's fault model (1985) using the fault parameters of DCRC-17a (Takahashi et al., 1995).The approach taken in this study is advantageous in that tsunami waves of both the interior and outer domain origin can be incorporated and the excessive computational load can be relaxed.The NH-model has been early used by authors to reproduce extreme wave characteristics of the 2004 and 2011 tsunamis (Kim et al., 2013a, b).
Simulations using the NH-model have been carried out to investigate the effect of the islands of Muen (south) and Hira (north) in the near-shore region off the coast of Monai Valley on tsunami wave characteristics in the near-shore region.Cases with two islands (representing the real topography), with one island (as in the physical model) and without islands have been compared to examine the effects on wave focusing and extreme run-up heights.The set-down of the sea surface elevation near the islands is examined in association with eddy generation.A description on the limitations of using the H-model for the extreme run-up is briefly included.

Model configuration
Figure 2 shows the model domains for the regional spherical coordinate H-model as well as the local Cartesian coordinate NH-model and the source fault region of tsunami generation, superimposed on the Google map.Three domains denoted by D1, D2 and D3 are subdomains of the regional H-model with different grid sizes, while D4 is the model domain for the local NH-model.The subdomain D1 covers the sea region of 137 • 00 E to 140 • 45 E and 40 • 00 N to 44 • 30 N, the subdomain D2 the sea region of 138 • 30 E to 139 • 45 E and 41 • 45 N to 43 • 15 N and the domain D3 the sea region of 139 • 22.5 E to 139 • 27.5 E and 42 • 02.5 N to 42 • 10.5 N. It is noted that the subdomain D1 is sufficiently large enough to cover the whole source region which lies approximately in the north-south direction, extending from 41 • 45 N to 43 • 15 N.The domain D4 covers the sea region of 4 km (east-west direction) by 4 km (north-south direction).
For the largest subdomain, D1, 450 by 540 horizontal meshes have been allocated with a resolution of 30 s (approximately 900 m).The second and third domains (D2 and D3) have the 450 by 540 and 150 by 240 horizontal mesh systems with resolutions of 10 and 2 s, respectively.The regional model uses a total of 10σ layers in the vertical direction.A one-way dynamic nesting has been used between subdomains of the regional H-model; the subdomains D1 and D2 are nested using 1 : 3 grid interpolation, while the subdomains D2 and D3 using 1 : 5 grid interpolation.Values from the subdomain D3 are interpolated in an offline manner to define the open boundary values for the local NH-model using approximately 1 : 5 grid interpolation.The local NH-model employs a z grid in the vertical direction (a total of 40 grid cells above mean sea level and 43 grid cells below the mean sea level).Table 1 summarizes the grid-related information used.Time steps of 1, 1 and 0.2 s have been used for numerical integration on D1, D2 and D3, respectively.The local NH-model calculates the best time step on D4 during simulation.Drying and wetting scheme are used in all subdomains except for the subdomain D1.
The bathymetry composed by the data set of GEBCO (Jones, 2003) and the Japanese Central Disaster Management Council (CDMC, 2003) was used to define the model water depth.The data of GEBCO cover the global ocean with 30 arcsec resolution, while the data of CDMC cover the coastal area of Japan with four different resolutions ranging from 1.3 km to 50 m.These data were composed and thereafter interpolated to the center points of each model.The fine topography data additionally provided from the Geospa- tial Information Authority of Japan in 10 m resolution were collected for the computation of the tsunami wave run-up using a wet-and-dry scheme.The masking work was followed using Google Earth to define wet or dry areas.In detail, the coast line was digitized using the satellite image using Google Earth and the location of each grid center point was then used for the determination of wet or dry areas.Comparison of the high-resolution seafloor bathymetry information existed before the event with bathymetric surveys after the event allowed a meaningful characterization of the seafloor deformation triggered by the tsunami.
For the determination of the initial water deformation we used Okada's fault model (1985), which was also used by Myers and Baptista (1995).The water deformation ( U ) is composed of the vertical displacement (U z) and the verti-  cal movement (U h) converted by horizontal movement (U x, Uy) at the bottom (H ).
Note that the initial elevations and velocities are assumed to be 0. The parameters for the calculation of water deformation taken from the DCRC-17a source model (Takahashi et 3 Numerical simulation

Model simulation using the regional model
To obtain the open boundary conditions for the local NHmodel used for the reproduction of the extreme run-up heights near the Monai Valley, the regional H-model has been initially applied with use of Okada's fault model (1985) for Table 2. Fault parameters of the 1993 Hokkaido-Nansei-Oki earthquake tsunami (Takahashi et al., 1995).Figure 5 shows the maximum run-up heights, computed from the application of the regional H-model, compared with the observed maximum run-up heights (Shuto and Matsutomi, 1995)   subdomain D3, we can see that the numerical simulation gives considerable disagreements in regions north of Monai Valley and south of Okushiri Island.However, the computed maximum inundation heights are underestimated in Monai Valley where extreme heights of 23 and 31.7 m were observed.The reason for this underestimation might partly be the hydrostatic model used and partly the poor model resolution.There is, as expected, the increase in the maximum run-up heights as the model resolution increases.
As a way to examine the local topographic effects, we have introduced a local amplification ratio that is estimated using the highest maximum run-up height by dividing the mean value and averaging the maximum run-up heights in the peripheral region.The peripheral region is defined in the neighboring region of 1 km from the highest maximum runup point to the north and south, excluding the nearby region of 200 m (see Fig. 3 for the neighboring region where the average of the run-up heights is calculated).The observed highest maximum and mean values are 31.7 and 15.9 m, respectively, and the local amplification ratio of observation is thus 1.96.The regional three-dimensional model produced the highest maximum of 23.5 m (blue cross within the black circle with second highest value) and neighboring mean of 15.6 m, giving the local amplification ratio of 1.51, while the model application by Titov and Synolakis (1997) produced the highest value of 29.7 m and mean value of about 17 m in the neighboring region, giving the ratio 1.75.That is, the present model underestimates the extreme-run up at the Monai Valley compared to the observations and results of Titov and Synolakis (1997).The three-dimensional model on D3 has a resolution of about 60 m, while the model by Titov and Synolakis (1997) has a variable grid system but with about 50 m from the shoreline.That is, the model resolutions of both models are comparable.It is then speculated that the presence of bottom friction causes the difference in the computed maximum run-up height near the Monai Valley.The difference in the moving boundary schemes is of course one of the causes.
The computed time series of tsunami wave heights and associated velocities are stored for the use in the local model.The vertical variation of the velocity was found to be insignificant.

Model simulation using the local model
From Fig. 1 we can guess intuitively that the presence of two small islands (Muen Island in the south and Hira Island in the north) in the near-shore region off the Monai Valley plays an important role in producing unusual extreme runup heights.To investigate the effects of islands, simulations have been carried out with three cases, namely cases without island, with one island (Muen Island only as in the laboratory model) and with two islands (which represents real topography) using the local NH-model.A no-slip condition is imposed at the sea bottom and k-turbulence closure is used.The present study has a clear distinction from those by Yoneyama et al. (2002) and Kakinuma (2005) in that the model has been applied with real topography and with more realistic open boundary conditions.
Figure 6 shows the topographic configuration of the central part of the local model domain with two islands, one island and without islands.The detailed bathymetry of 1 m intervals is added for reference purpose.The two islands are located about 500 m away from the coast near the Monai Valley region, while the distance between the two island centers is about 300 m.The western open boundary of the subdomain D4 is located 1.5 km from the two islands.The land heights of Muen and Hira islands are about 10 and 50 m, respectively.The slope of coastal land height at Monai Valley is estimated to be approximately 30 • .
Figure 7 displays the distribution of the maximum sea surface elevation computed using the local model with the three different topographic configurations shown in Fig. 6.It is evident that calculation in the presence of two islands gives rise to the highest value of the extreme run-up heights at Monai Valley.Examination of results gives a maximum of 32.3 m and the amplification ratio of about 1.90, showing a good agreement with the observed values (31.7 m and 1.96, respectively).
Calculations with one island (Muen Island) produces the extreme height of about 30.8 m, which is slightly smaller than the observed value.Calculations by Kakinuma (2005) of the same island produced a height of 30.6 m.The slight difference between cases with one island and two islands may be because the land height of Hira Island is relatively low enough to be mostly inundated by tsunami waves.Calculation without islands is found to give a maximum runup height of about 25.2 m, indicating that the presence of is-  lands plays an important role in generating the extreme runup height near the Monai Valley.
Some interesting features are additionally shown to the offshore direction.For the case without islands there is a tendency for the maximum sea surface elevation to gradually increase from Okushiri Island in the offshore direction.The presence of islands, however, alters the tendency.Specifically, the tendency of the maximum sea surface elevation to increase on the west sea region of islands appears for the cases with one island and with two islands.This may be due to the reflection effect by islands.Furthermore, there is a hint of set-down of the sea surface elevation near the Muen Island to the southeast direction, approximately parallel to the isobaric contours over the steep bottom slope.The set-down phenomena will be discussed below in more detail.
For instructive purposes, we examine the snapshots of time-varying sea surface elevation fields computed at four different times in the presence of two islands (Fig. 8).It is evident that at t = 205 s the incoming tsunami waves approach the coastal area near the Monai Valley region, and at t = 220 s the coast region is entirely flooded and run-up starts from the north of Monai Valley.The two snapshots of t = 220 s and t = 240 s show the southward transition and focusing of the tsunami run-up to the southern end of the concave coast, resulting in the extreme run-up height near the Monai Valley.
In addition, the occurrence of set-down of the sea surface elevation is evident from Fig. 8.The set-down may be divided into two: one deep set-down right behind Muen Island and another, relatively shallow, elongated set-down between Muen and Hira islands.The reason why the intense set-down and the rotational flows do not appear on the lee-side of the island might be the bottom topographic effects.The horizontal distribution of the set-down can be shown in more detail from streaklines at t = 240 s superimposed on the sea surface elevation fields (Fig. 9).Streaklines around the islands computed for cases with two islands, with one island and without islands have been compared by releasing a sequence of particles.It is evident that the calculation with two islands produces a counterclockwise rotational flow associated with the deep set-down southwest of Muen Island and the shallow elongated set-down between Muen and Hira islands.It is noted that the formation of rotational flows is not accompanied by the elongated set-down, probably because the Hira Island with relatively low topographic height is almost flooded.Calculation with one island only shows the setdown southwest of Muen Island with the size and intensity less than that of the two island case, implying the presence of some interaction between the two set-down processes.The streakline pattern at the coast shows focusing on the southern end of the concave coast line.This may be attributed to the slanted layout of the coast.We note that none of the intense set-down of the sea surface elevation and formation of rotational flows appear in the study of tsunami run-up around an idealized conical island by Choi et al. (2007).The fact that the ratio of the incoming wave amplitude to water depth was about 0.1, as well as the consideration of a small island and use of a flat bottom, might be the cause.Calculation without islands shows the refraction rather than diffraction.
Figure 10 shows the snapshots at t = 220 and 240 s of the sea surface elevation and the superimposed instantaneous flow velocities near the Monai Valley coast computed with two islands, with one island and without islands.It is seen from the comparison of upper diagrams that calculation without islands gives rise to the earliest start of the tsunami runup, followed by the calculation with one island and finally, slightly later, the calculation with two islands.The maximum run-up heights at t = 240 s are in the reverse order; that is, the late arrival near the Monai Valley coast leads to the occurrence of higher run-up height than the early arrival.The late arrival in fact implies more diffraction, resulting in more focusing of tsunami waves.
We finally examine the vertical velocity distribution.Figure 11 displays the instantaneous vertical fields at t = 240 s for three different cases.It is evident that the cases with one island and with two islands produce the intense vertical velocities significantly stronger than the case without islands.The model gives maximum values of about 5 ms −1 in the presence of islands and 3 ms −1 without islands near the coast of the Monai Valley.Calculation with two islands produces slightly stronger vertical velocities than calculation with one island but, again, the difference is marginal.Comparison with two-dimensional model results on the same grids (not shown) indicates that the use of a non-hydrostatic model produces considerably stronger horizontal and vertical advections, resulting in considerably higher values of extreme run-up heights.

Conclusions
The occurrence of extreme run-up height (31.7 m) at Monai Valley west of Okushiri Island has been investigated in this study using the local NH-model with the open boundary conditions supplied in a one-way offline manner from the multi-domain regional H-model.The well-known Okada fault model (1985) has been used for the determination of the initial water deformation using the fault parameters of DCRC-17a (Takahashi et al., 1995).The modeling procedure used in this study is found to be efficient, in that accurate results at the target area can be obtained with reduced computational load compared with the use of the local model over the large area covering the entire fault region.Errors in nesting the H-and NH-models might have a marginal influence on the accuracy of results because there is no significant vertical structure in flow velocity.
Three experiments using the NH-model, namely experiments without islands, with one island and with two islands, have been compared.We could see that the experiment with Muen (south) and Hira (north) islands west of Okushiri Island coast line give rise to the extreme run-up height of about 32.3 m.Model experiments indicate that the shape and layout of the coast line as well as the diffraction of tsunami waves by the two islands has led to focusing of tsunami waves to the direction of Monai Valley, giving the extreme run-up height there.It has been noted that use of the H-model produces an extreme run-up height significantly smaller than the observed 31.7 m at Monai Valley, regardless of the presence of islands.The underestimation of the local vertical velocity and thereby the vertical advection at the Monai Valley might be the main cause.The model experiments with the islands and their diffraction properties demonstrated once again that the islands do not protect the coastline behind them from the incoming tsunami front.For a better and more efficient reproduction of the extreme run-up height in the future, application of three-dimensional σ coordinate NH-models, which are recently in use among the oceanographic community, obviously needs to be considered without any nesting.It may be worth comparing performance between hydrostatic and nonhydrostatic models with different turbulence schemes.Strict tests of the drying and wetting scheme are, however, required prior to the application.

Figure 1 .
Figure 1.Location map of Monai Valley showing Hira and Muen islands off the western coast of Okushiri Island, Hokkaido, Japan (Google Earth).
reproduced the laboratory experiment using two three-dimensional models (STOC-IC and STOC-VOF) contained in the modeling system STOC (Storm surge and Tsunami simulator in Oceans and Coastal areas).The two models were based on Reynoldsaveraged Navier-Stokes equations, and STOC-VOF reproduced a maximum run-up height of 30.6 m.In this study, numerical simulations with real topography have been made to investigate the tsunami wave characteristics in the near-shore region of Okushiri Island and the extreme run-up height (31.7 m) near Monai Valley.Results of previous studies led to the use in this study of a Reynoldsaveraged Navier-Stokes equation-based three-dimensional model

Figure 2 .
Figure 2. Domains for the tsunami model experiments following the 1993 Hokkaido-Nansei-Oki earthquake and location of tsunami initial conditions.Subdomains D1, D2 and D3 are for the H-model, while the subdomain D4 is for the NH-model.
Figure 3 shows the bathymetry of the domain D3, including land height contours up to a maximum of 50 m.Negative values are used for water depth, while positive values are used for land heights.It is clear that the water depth in D3 is mostly less than 100 m except for the open sea region along the western open boundary.The presence of two small islands (Muen Island in the south and Hira Island in the north) and a slightly concave form of the coast line characterize the Monai Valley region.

Figure 3 .
Figure 3.The bathymetry of the third domain (D3).The square box presents the fourth domain (D4) for the NH-model.The black solid lines of 20 m intervals represent water depth, and the blue solid lines of 10m intervals represent the land heights up to 50 m.The blue arrows indicate the regions where the mean run-up height value is calculated for the local amplification ratio.

Figure 4 .
Figure 4.The total vertical component of water deformation computed using the Okada fault model (1985).The water deformation is computed by use of the three sub-fault parameters (DCRC-17a) summarized in Table2.

Figure 5 .
Figure 5. Left: the observed and computed maximum inundation heights; right: the locations of inundation observations (black open circles) and the locations of computation points in the subdomains D1 (red marks), D2 (green marks) and D3 (blue marks) on the west coast of Okushiri Island.
on the west coast of Okushiri Island.To examine the effects of model resolution on the maximum run-up heights, the results computed at D1, D2 and D3 grid points in the vicinity of the west coast of Okushiri Island are presented together.On the right diagram of Fig. 5, the black open circles represent survey locations, while the small red, green and blue crosses represent the computation nodes of the subdomains D1, D2 and D3, respectively.The black open circles on the left represent the observed maximum run-up heights at the survey points on the right diagram, while the red, green and blue crosses on the left diagram show the computed results of the maximum water heights at the corresponding nodes shown on the right.Considering the results from the

Figure 6 .
Figure 6.Bathymetry and topography used (a) with two islands (b is a detailed contour map with two islands), (c) with one island and (d) without islands.

Figure 7 .
Figure 7.The maximum sea surface elevation computed (a) with two islands, (b) with one island and (c) without islands.Red contour lines indicate original coastlines.

Figure 8 .
Figure 8.The snapshots of time-varying sea surface elevation distribution computed with two islands.

Figure 9 .
Figure 9. Streaklines in the vicinity of islands computed at t = 240 s (a) with two islands, (b) with one island and (c) without islands.

Figure 10 .
Figure 10.Snapshots of sea surface elevation and instantaneous total flow velocities at t = 240 and 240 s computed with two islands (left panels), with one island (middle panels) and without islands (right panels).

Figure 11 .
Figure 11.Instantaneous vertical velocities computed at t = 240 s (a) with two islands, (b) with one island and (c) without islands.

Table 1 .
Information on grids and time steps.