Comparative numerical simulation of tsunami waves during the January 1, 2024 Noto Peninsula earthquake, Japan

Full Text

Introduction

It is well known that Japan is among ten most earthquake-prone countries in the world ¹. Its location within the Pacific Ring of Fire has resulted in numerous destructive earthquakes throughout its history. One notable event occurred on 17 January 1995, when the M = 7.3 earthquake struck the western region of the Japanese island of Honshu resulting in over 6400 fatalities [1]. Another significant earthquake, with M = 9, followed by a tsunami measuring between 3 and 15 meters in height, struck the coasts of Honshu and Hokkaido islands on 11 March 2011. As a result of this natural disaster ², more than 28 thousand people died and large areas were contaminated with radioactive substances from the Fukushima I Nuclear Power Station [2, 3].

High seismic activity is stipulated by the location of the archipelago at the collision of several lithospheric plates, two of which – the Philippine and Pacific – are in motion. This leads to a significant number of intense Earth tremors usually accompanied by earthquakes and tsunamis. Several thousand earthquakes occur in Japan per year, an average of 18 earth tremors per day. However, strong earthquakes accompanied by tsunamis occur much less frequently: once every 10 years – with M = 8, once a year – with M = 6 (see 1, ³, ⁴ and [4]).

We consider a strong earthquake (M = 7.6) occurred on the Noto Peninsula in Ishikawa Prefecture, Japan, on 1 January 2024. As a result of a series of Earth tremors, massive destruction was caused in the towns of Suzu, Wajima, Noto and Anamizu. At the same time, significant damage was also recorded in the neighboring prefectures of Toyama and Niigata. As a result of this earthquake, a tsunami hazard was declared in several regions of Honshu Island. Japan Meteorological Agency issued a tsunami warning map, according to which, presumably, tsunami waves of 5 m in height approached the western coast of the island (Fig. 1).

This paper presents the results of numerical simulation of tsunami wave generation by a seismic source formed within the keyboard model of an earthquake source [5]. The greatest uncertainty in computation possible tsunami development scenarios is associated with the initial motion of the bottom during an earthquake. When using keyboard model, specifying the displacements and velocities of the bottom during an earthquake provides a solution to the problem of tsunami wave formation. In this case, the key-block motion is simulated according to the aftershock stage of the process. After the end of the aftershock stage, the blocks stop. During the strongest earthquakes with a rupture length of several hundred kilometers, their sources can cover many blocks.

In this paper, we consider a five-block earthquake source as well as generation of a tsunami source at different localizations of a seismic source and different kinematic processes of key-block motion in the source. Analysis of wave characteristics obtained during simulation of different localizations and dynamics of the seismic source leads to conclusions about the ambiguity of the selected process model.

The work aim is a comparative numerical simulation of generation and propagation of tsunami waves induced by the M = 7.6 earthquake source on the Noto Peninsula, Ishikawa Prefecture, Japan, on 1 January 2024.

 

Fig. 1. Tsunami hazard map for Honshu Island published by the Japan Meteorological Agency on 1 January 2024 ³

 

Problem statement

Aftershock stage of the process

Using the data from [6–9], we have analyzed the sequence of the aftershock stage of earthquake process. Aftershock data were obtained from the U.S. Geological Survey ⁵ and are represented in Fig. 2. After the main earthquake shock on 01.01.2024 at 07:10:09 UTC, 42 more shocks with M > 4.5 occurred within three days (Table 1).

 

Fig. 2. Aftershocks after the 01.01.2024 earthquake on the Noto Peninsula (yellow star is the earthquake epicenter; circles with numerals indicate the order of aftershocks ⁵)

 

Table 1. All aftershocks at M > 4.5 from 01.01.2024 to 03.01.2024 after the earthquake on 1 January 2024

s/n	M	Time (UTC)	s/n	M	Time (UTC)	s/n	M	Time (UTC)
0	7.5	07:10:09	15	4.9	08:22:10	30	4.6	17:33:30
1	6.2	07:18:41	16	4.7	08:29:02	31	4.5	19:27:14
2	4.6	07:27:28	17	4.6	08:42:04	32	4.6	19:42:06
3	4.6	07:28:15	18	4.5	08:48:21	33	4.5	22:13:30
4	4.7	07:29:04	19	5.5	09:03:48	34	5.4	01:17:31
5	5.2	07:39:49	20	4.6	09:06:13	35	4.8	06:57:54
6	4.8	07:42:43	21	5.6	09:08:17	36	4.6	08:13:40
7	4.7	07:45:30	22	5.2	09:30:21	37	4.8	17:21:47
8	4.8	07:48:14	23	5.0	09:39:59	38	5.3	01:54:34
9	4.7	07:54:21	24	5.1	09:54:26	39	4.9	03:54:13
10	5.6	07:56:47	25	4.6	09:49:15	40	4.7	09:48:00
11	5.0	08:02:44	26	4.6	10:06:54	41	4.6	15:36:53
12	5.1	08:07:10	27	4.7	10:50:35	42	4.8	19:38:41
13	4.5	08:15:38	28	4.6	11:35:32
14	4.8	08:17:46	29	4.5	13:19:57

 

Tsunami source formation within the keyboard model of an earthquake seismic source

Since the mechanism of seismic source formation is unknown, then it is possible to calculate approximate dimensions of a seismic source and its dynamics by specifying the source mechanism from tectonic considerations based on known magnitude values of the earthquake and aftershocks and also using the distribution of M > 4.5 aftershocks. To calculate the source dimensions, the following formulas are applied [10]

 $\left\{\begin{array}{l}\lg L=0.59 M-\mathrm{2.44,} \\ \lg W=0.32 M-\mathrm{1.01,}\end{array}\right.$                      (1)

where M is earthquake magnitude; L is length of the rupture at the source, km; W is width of the rupture plane, km.

The maximum vertical displacement of the wave surface above the earthquake source is found by the formula ⁶

lg(H) = 0.8М – 5.6,                        (2)

where Н is maximum height of wave surface vertical displacement above the earthquake source, m. Estimated displacements of these characteristics applied for simulating the tsunami source 6 [10] are presented below:

 

Mw	L, km	W, km	S, km2	Н, m
7.5	78 ± 18	43 ± 6	3354 ± 440	4.8

 

For simulating, we took source length to be 78 km and the width to be 43 km. Since the ocean surface will rise by the same value as the block on the bottom has displaced due to fluid incompressibility and pressure hydrostaticity, and from the Iida formula the displacement was obtained to be 4.8 m, then we took the maximum vertical shift of the block in the earthquake keyboard source to be 4.8 m.

Mathematical statement of the problem

The process of tsunami wave generation is considered in the shallow-water theory approximation. For simulating, we used the equations that describe a nonlinear system of shallow-water equations in a 2D formulation (see, for example, ⁷ and [11]):

 $\left\{\begin{array}{c}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial v}{\partial y}+g\frac{\partial \eta }{\partial x}=\mathrm{0,}\\ \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+g\frac{\partial \eta }{\partial y}=\mathrm{0,}\\ \frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x}\left[\left(\eta +H-B\right)u\right]+\frac{\partial }{\partial y}\left[\left(\eta +H-B\right)v\right]=\frac{\partial B}{\partial t}.\end{array}\right.$         (3)

We assume that axis z is directed vertically upwards; x, y are spatial coordinates; u(x, t), v(у, t) are horizontal velocity components, η(x, у, t) is free surface disturbance relative to the undisturbed level; H is maximum depth of the basin; function В(х, у, t) determines the basin bottom change (taking into account the characteristics of the dynamic seismic source). At the initial moment of time, the parameters of dynamic seismic source (coordinates, displacement rate of the key-blocks) are specified in the solution domain. We assume that the water area with the initial bottom shape is at rest before the generation start, i.e., the velocity and disturbance of free surface were absent:

η(x, y, 0) = 0; u(x, 0) = 0; v(x, 0) = 0.

At the last seaward point at 5 m depth, a condition of total reflection (vertical wall) was set, which makes it possible to record the maximum and minimum wave level shift at this depth. Computational domain applied for these calculations was 125.01°– 147.00°E, 30.01°– 55.00°N with a grid of 30' ≈ 759 m step. In the numerical solution, we used a scheme constructed by analogy with the Sielecki one [3].

Numerical simulation of tsunami waves within the keyboard model of a seismic source

Figure 3 represents the computed water area with localization points of tide-gauge stations [7].

Four scenarios with different localizations of earthquake source (Fig. 4) were considered in this work. Figure 4, a shows scenario 1 with block 2 in the Toyama Bay area southeastwards of the Noto Peninsula. In Fig. 4, b, the same block is located in the Sea of Japan northwestwards of the Noto Peninsula. The location of the remaining blocks in the earthquake source is similar for these scenarios. Figure 4, c shows the source localization for scenario 3. It can be seen that the block located northwestwards of the peninsula took a triangular shape stretching along the entire peninsula. The shape of the remaining blocks also changed. In scenario 4, the shape of blocks 1–3 did not change, but direction of the location of blocks 4 and 5 changed. This change in the localization, size and orientation of the blocks is due to a comparison of simulation results with in situ data and the data from other authors. A total of 14 simulation options were carried out, with 4 of them being presented in this paper.

Table 2 shows kinematics of the key-block motion in the earthquake source. For all scenarios, the source consists of five blocks that perform sequential motion at certain time intervals. Analysis of time sequence of aftershock occurrence makes it possible to create an estimated version of key-block motion kinematics in the seismic source.

 

Fig. 3. Estimated water area (yellow star is the earthquake epicenter; red triangles denote the location of tide-gauge stations [7])

 

In scenario 1, the blocks make sequential motion at equal time intervals of 30 s. In scenario 2, the blocks make sequential motion at equal time intervals of 30 s but with negative displacement values. In scenarios 3 and 4, the block displacements in different directions remain. In scenario 3, the blocks make sequential motion at different time intervals in the range of 60–140 s. In scenario 4, the blocks make sequential motion at different time intervals from 60 to 130 s. Fig. 5 illustrates the process of tsunami source generation in scenarios 1 and 4 which clearly corresponds to the earthquake source locations shown in Fig. 4, a and d.

 

Fig. 4. Simulation of the scenarios of seismic source different locations: 1 (а); 2 (b); 3 (c); 4 (d) (1 – 5 are the key-block numbers)

 

Table 2. Kinematics of block motion in the earthquake source in four scenarios

Parameter	Block number
	1	2	3	4	5
Scenario 1
Displacement height, m	1.5	3	1	1	1
Start time of motion, s	30	0	90	60	120
End time of motion, s	60	30	120	90	150
Scenario 2
Displacement height, m	2	1	−1	1	−1
Start time of motion, s	0	30	90	60	120
End time of motion, s	30	60	120	90	15
Scenario 3
Displacement height, m	3	1.4	−1.3	1.5	1.2
Start time of motion, s	0	90	190	350	520
End time of motion, s	30	120	220	380	550
Scenario 4
Displacement height, m	3.2	−1.1	−0.2	0.2	1.3
Start time of motion, s	0	90	190	350	440
End time of motion, s	30	120	220	380	470

 

In Fig. 5, it is seen that the water surface displacement corresponds to block displacement presented in Table 2 for scenarios 1 and 4. The main difference is associated with the change in the localization and shape of block 1 with a displacement towards the southeast (scenario 1) or northwest (scenario 4) from the Noto Peninsula. The configuration of blocks 2–5 does not differ significantly.

Fig. 6 shows the positions of wave fronts for nine points in time. It is clearly seen that 1-meter waves reach Sado Island 10 minutes after the start of generation, and the waves move in this direction faster than towards the bay; this is associated with the Toyama Trough. After 30 minutes, the wave front reached the coastal cities of Kashiwazaki and Toyama. At the 45th minute, the 0.5-meter wave front reached the northern part of Honshu Island and went around the Noto Peninsula. 1 hour 45 minutes after the earthquake, waves of ~20 cm height reached Vladivostok and the southern part of Hokkaido Island. At 05:56:30, the wave front reached Kholmsk, located in the northern part of Sakhalin Island.

The distribution of heights in Fig. 7 shows good agreement with Fig. 1 where the most hazardous areas are the coasts located in the northeastern part of the Noto Peninsula as well as the cities of Toyama, Sado and Kashiwazaki. The computation is carried out up to the 5-meter isobath.

 

Fig. 5. Generation of a tsunami source at six time points in scenarios 1 (a) and 4 (b)

 

Fig. 6. Propagation of tsunami waves across the computed water area in scenario 4

 

Fig. 7. Distribution of maximum wave heights across the computed water area in scenarios 1 (a) and 4 (b)

 

Using the computation data for scenario 4, 3D histograms of height distribution in Toyama Bay, Chuba, Tohoku, Kansai regions and Sado Island were constructed (Fig. 8). It is clearly seen that the average height of the incoming waves is 1 m in Toyama Bay and 5 m in the western part of the Noto Peninsula.

 

Fig. 8. 3D histograms of wave heights: a – Toyama Bay; b – Chuba and Kansai regions; c – Tohoku region and Sado Island (blue and green color correspond to sea and land respectively)

 

Fig. 8, b shows that the wave height was 3 m in the southern part of the Noto Peninsula and the average wave height was 1 m further in the Chuba region. Fig. 8, c demonstrates clearly that the average height of incoming waves in the northern part of the Tohoku region was 1 m, but the wave height in the northern part of Sado Island was 2.5 m on average.

 

Analysis of results of tsunami numerical simulation during the earthquake on 01.01.2024 in Japan

Fig. 9 shows 2D histograms of tsunami wave heights for different sections of the computed coasts in four scenarios under consideration.

 

Fig. 9. 2D histograms of wave heights: a – Toyama Bay; b – Chuba and Kansai regions; c – Tohoku region and Sado Island

 

Fig. 9, a demonstrates clearly that wave height is 3 m in scenarios 3 and 4 in the area of Wadzima town; the wave height in scenario 4 is less than a meter and it is 1–1.5 m in scenarios 1, 2 and 3 in Itoigawa and Joetsu. Fig. 9, c shows that the wave heights in scenario 1 on the coasts of the Chuba and Kansai regions are on average 0.5 m less than in scenarios 3 and 4. This is due to the fact that block 1 is located inside Toyama Bay in scenario 1 and in the opposite part of the Noto Peninsula in the other scenarios. Fig. 9, c illustrates that the wave height in scenario 1 is greater than in the other scenarios near Toyama Bay (up to 38°N) and northwards of the bay the wave heights of scenario 4 predominate. This is also due to the location of the blocks, since the greatest wave heights of scenario 1 propagate from Toyama Bay and the waves of scenarios 2, 3, 4 go around Sado Island and reach the northern part of Honshu Island with large amplitudes. Table 3 provides data for all performed scenarios (scenarios 1–4) and in situ data from tide-gauge stations 5.

 

Table 3. Sea level maximum rise (cm) near the settlements located on the estimated water area coast

Settlement	In situ data 5	Scenario 1	Scenario 2	Scenario 3	Scenario 4
Saigo	27.1	34.2	24.0	20.1	21.9
Mikuni	31.5	68.7	45.1	28.6	30.2
Toyama	79.4	124.8	81.8	78.6	85.0
Kashiwazaki	40.2	121.3	71.2	53.5	65.4
Sado	32.0	37.2	20.1	31.4	24.6
Oga	27.4	62.3	26.0	31.2	31.4
Fukaura	33.9	25.2	10.0	23.4	28.5
Vladivostok	28.0	34.0	18.4	29.5	27.8
Hakodate	9.9	2.6	1.5	3.9	3.3
Wakkanai	8.5	2.8	1.8	5.7	7.5
Sosunovo	11.4	7.1	3.5	7.6	10.7
Kholmsk	12.0	1.2	1.0	3.5	4.1

 

Discussion of results

The purpose of this work was to obtain wave characteristics in the Sea of Japan waters based on simulation for which 14 scenarios were considered; 4 of them with the best wave characteristics of the process are presented in this work. Available publications on this earthquake were analyzed during the study [6–9]. A comparison of the obtained results with the data from the above publications (in most detail with the data from [6, 7] (Table 4)) performed for simulation earthquake sources, the localization of which is close to that proposed in our work, shows that the data of our numerical simulation are mainly close to the tide-gauge data at the corresponding stations with an accuracy of 5–6 cm. Where a significant difference between the tsunami amplitude and the calculated height takes place, it can be explained by complex physical processes in Toyama Bay, such as re-reflection, edge waves and bay resonance phenomena [6–9].

 

Table 4. Sea level maximum rise (cm) near the settlements located on the estimated water area coast (data from different sources)

Settlement	In situ data [7]	Scenario 4	Data from [6]
Saigo	27.1	21.9	22
Mikuni	31.5	30.2	40
Toyama	79.4	85.0	75
Kashiwazaki	40.2	65.4	60
Sado	32.0	24.6	10
Oga	27.4	31.4	18
Fukaura	33.9	28.5	19
Vladivostok	28.0	27.8	5
Hakodate	9.9	3.3	–
Wakkanai	8.5	7.5	6
Sosunovo	11.4	10.7	3
Kholmsk	12.0	4.1	3

 

Thus, at Kashiwazaki station, all models, both our and those used in [6, 7], overestimate the tsunami amplitude (from 89.4 to 227.1 cm). According to our computation, the data spread at this point was from 53.5 to 121.3 cm with a change in the localization and dynamics of the earthquake source. However, scenario 4, chosen as the most adequate, provided a calculated value of 65.4 cm, which also exceeded the observed maximum amplitude at this point.

According to the observations, the maximum amplitude in the Toyama station area was 79.4 cm; according to scenario 4, we obtained a maximum wave displacement of 85.0 cm, i.e., a difference of 5.6 cm. In [7], the maximum amplitude for this point was 79.4 cm according to simulation data, and in [6] – 75 cm. In other words, the difference between the calculated and in natural data in [6, 7] was 4.4 and 4.2 cm, respectively, i.e., 1.2 and 1.4 cm less than our corresponding data. Although the differences in both computed from [6, 7] and our are, in our opinion, within the framework of the source simulation accuracy, work [6] suggests that such a difference is possible in the presence of an additional landslide source to the seismic one, which is confirmed by the results of simulation performed by the authors of the work. The data in the remaining points in our computation have values that are quite close to the observed amplitudes, except for Wakkanai and Kholmsk. Thus, the maximum wave amplitude in Vladivostok is 28 cm and the computed values are 27.8 cm. These values are 11.4 and 10.7 cm, respectively, and their difference is 0.7 cm in Sosunovo village.

Conclusion

The computations using the block-keyboard earthquake model presented in this paper showed that it was especially important to take into account bottom geomorphology and geodynamics of tectonic processes for the earthquake that occurred in the Sea of Japan on 01.01.2024. This model permits to consider the initial stress distribution in the earthquake preparation zone as well as the dynamic transient process of forming the distribution of bottom displacements. Earthquake sources of various localizations with different shapes of the blocks that comprise it were considered. The most appropriate shape for simulating an intraplate active fault was the source with a long triangular block in the northwest of Nota Island. Numerical simulation results showed that the keyboard earthquake model made it possible to simulate adequately even such complex earthquake sources as the one that occurred in the northwest of Nota Island. A comparison of the obtained computation data with the amplitude of the maximum tsunami wave height from tide-gauge stations provided an value of less than 6 cm, except for three points.

 

¹ Regan, H., Akbarzai, S., Kobayashi, C. and Maruyama, M., 2024. ‘Battle against Time’ to Find Quake Survivors as Japan Lifts Tsunami Warnings and Death Toll Rises. 2024. [online] Available at: https://edition.cnn.com/2024/01/02/asia/japan-earthquake-tsunami-warnings-tuesday-intl-hnk/index.html [Accessed: 28 February 2024].

² TASS. Chronology of Earthquakes with Fatalities in Japan. 2024. [online] Available at: https://tass.ru/info/19667253 [Accessed: 28 September 2024].

³ Japan Meteorological Agency. Weather Map. 2024. [online] Available at: https://web.archive.org/web/20240101100413/ https://www.jmago.jp/bosai/map.html#5/38.891/141.24/&elem=warn&contents=tsunami [Accessed: 10 May 2024].

⁴ USGS. Latest Earthquakes. 2024. [online] Available at: https://earthquake.usgs.gov/earthquakes/map [Accessed: 10 May 2024].

⁵ IOC. Sea Level Station Monitoring Facility. 2024. [online] Available at: https://www.iocsealevelmonitoring.org/map.php [Accessed: 10 May 2024].

⁶ Pelinovsky, E.N., 1982. [Nonlinear Dynamics of Tsunami Waves]. Gorky: IPF AN SSSR, 226 p. (in Russian).

⁷ Voltsinger, N.E., Klevanny, K.A. and Pelinovsky, E.N., 1989. Long-Wave Dynamics of the Coastal Zone. Leningrad: Gidrometeoizdat, 272 p. (in Russian).

About the authors

Raisa Kh. Mazova

Nizhny Novgorod State Technical University n. a. R. E. Alekseev

Author for correspondence.
Email: raissamazova@yandex.ru
ORCID iD: 0000-0003-2443-149X
Scopus Author ID: 6506297372

Professor of Applied Mathematics Department, Professor

Russian Federation, Nizhny Novgorod

Aleksandr A. Martynenko

Autonomous Non-Commercial Organization of Higher Education “Central University”

Email: martynenko.busy@gmail.com

Master’s Student

Russian Federation, Moscow

Andrey А. Kurkin

Nizhny Novgorod State Technical University n. a. R. E. Alekseev; V. I. Il'ichev Pacific Oceanological Institute, Far Eastern Branch of RAS

Email: aakurkin@gmail.com
ORCID iD: 0000-0003-3828-6406
Scopus Author ID: 7003446660
ResearcherId: A-1972-2014

Vice-Rector for Research, Leading Research Associate, DSc. (Phys.-Math.), RAS Professor

Russian Federation, Nizhny Novgorod; Vladivostok

References

Supplementary files