Numerical Study on the Turbulent Structure of Tsunami Bottom Boundary Layer Using the 2011 Tohoku Tsunami Waveform

: In this study, the tsunami-induced bottom boundary layer was investigated based on actual waveforms obtained by the GPS buoys along the coast of the Tohoku region during the 2011 Great East Japan Earthquake tsunami. The k- ω model was utilized for the numerical analysis in this study. As a result, the tsunami boundary layer thickness was found to be extremely thin compared to the water depth. The velocity distribution was similar to that of the bottom boundary layer under wind-generated waves. The ﬂow regime is located in the transition from smooth turbulence to rough turbulence. Because of this, the gradient of the ﬂow across the layer is much greater than the gradients in the steady ﬂow direction. Therefore, the bottom friction is underestimated if the steady friction factor, such as in the Manning formula, is used. This study proposes a new simple method for calculating the bottom shear stress due to an irregular tsunami based on the wave friction law, and the k- ω model results are used to validate the proposed methods.


Introduction
One of the most distinct differences in wave-induced bottom boundary layers compared with steady open channel flow is an extremely thin boundary layer thickness with a steep velocity gradient near the sea bottom. In most wave theories, inviscid treatment ignoring the shear stress effect is sufficient, even for practical engineering applications. However, when we consider phenomena related to the frictional impacts, such as wave height attenuation and sediment movement, the inviscid theories are insufficient, and consideration of bottom friction is essential. For this reason, after the pioneering experiment by [1] using an oscillating tunnel, many investigations have been undertaken regarding the bottom boundary layer formed by wind waves.
In contrast to the wind-induced wave boundary layers, Manning's roughness coefficient has generally been utilized in numerical simulations of tsunami to plan and design disaster prevention structures in various practical situations. However, according to the recent investigation by the authors, even if the wave motion is fully classified into a long wave condition, the boundary layer thickness can exhibit similar behavior under the ordinary wave boundary [2,3]. Furthermore, during the 2010 Chilean Earthquake Tsunami, typical characteristics of the wave boundary layer were observed from the measurement of the flow velocity distribution under the tsunami [4,5]. In addition, [6,7] found that the boundary layer near the tsunami source area is surprisingly in the laminar flow regime. However, in investigating the tsunami boundary layers, the waveform is usually assumed to comprise sinusoidal waves or solitary waves, which are totally different from the actual irregular tsunami waveform.
According to the linear long wave theory, the waveform of a water surface and a flow velocity waveform is the same. Moreover, according to boundary layer approximation, the pressure gradient that acts in a boundary layer can be replaced by the acceleration obtained from the flow velocity waveform. Therefore, the pressure gradient in the boundary layer under the actual tsunami, which has an irregular waveform, generally shows a complicated time variation. For this reason, it differs from the pressure gradient based on the simple waveform described by sinusoidal waves or solitary waves. It is well known that the turbulent flow transition in a boundary layer is generally greatly dependent on a pressure gradient [1].
This study presents a numerical analysis of the turbulent structure of the tsunami bottom boundary layer using the measured 2011 Tohoku Tsunami waveform using the data at several observation points of the GPS wave gauge.

GPS Wave Gauge Data
This investigation used tsunami data [8] from six locations, as shown in Figure 1a. The result for each measuring point is illustrated using the station number shown in Table 1. It is to be noted that Figure 1b,c is explained in detail later.
J. Mar. Sci. Eng. 2022, 10, x FOR PEER REVIEW to comprise sinusoidal waves or solitary waves, which are totally different from the irregular tsunami waveform.
According to the linear long wave theory, the waveform of a water surface and velocity waveform is the same. Moreover, according to boundary layer approxim the pressure gradient that acts in a boundary layer can be replaced by the accele obtained from the flow velocity waveform. Therefore, the pressure gradient in the b ary layer under the actual tsunami, which has an irregular waveform, generally sh complicated time variation. For this reason, it differs from the pressure gradient bas the simple waveform described by sinusoidal waves or solitary waves. It is well k that the turbulent flow transition in a boundary layer is generally greatly dependen pressure gradient [1].
This study presents a numerical analysis of the turbulent structure of the tsu bottom boundary layer using the measured 2011 Tohoku Tsunami waveform usin data at several observation points of the GPS wave gauge.

GPS Wave Gauge Data
This investigation used tsunami data [8]

from six locations, as shown in Figu
The result for each measuring point is illustrated using the station number shown in 1. It is to be noted that Figure 1b,c is explained in detail later.  In order to investigate the tsunami-induced bottom boundary layer character the estimated tsunami-induced velocity was based directly on the measured water s elevations at the GPS stations, as similar to [6]. The equation can be expressed as fo  In order to investigate the tsunami-induced bottom boundary layer characteristics, the estimated tsunami-induced velocity was based directly on the measured water surface elevations at the GPS stations, as similar to [6]. The equation can be expressed as follows: where U(t) is the flow velocity, η(t) is the water surface elevation, h is the still water depth, and g is the gravitational acceleration. Figure 2 shows the flow velocity results from Equation (1) for six stations, which indicate a highly irregular waveform.
Ut is the flow velocity, () t  is the water surface elevation, ℎ is the depth, and is the gravitational acceleration. Figure 2 shows the flow velocity results from Equation (1) for six stations dicate a highly irregular waveform.

Numerical Analysis Using the k-ω Turbulence Model
The present study used the same k-ω model as in the previous study by [2 1-D incompressible unsteady flow, the equation of motion within the boundar be expressed as follows: where is the horizontal velocity in the boundary layer, is the time, p is th  is the fluid density,  is the shear stress, x is the horizontal coordinate, a vertical coordinate taken as positive upward from the bottom surface.
At the location outside of the boundary layer, the pressure gradient term the temporal variation free-stream velocity, () Ut , obtained from Equation (1 tion, the shear stress for turbulence flow is calculated by introducing the edd model. Hence, Equation (2) can be rewritten as follows:  is the molecular kinematic viscosity, and t  is the eddy viscosity.
Although the nonlinear advection term is not considered in Equation (3), ized boundary layer equation has been successfully employed and validated ev linear wave-induced boundary layers by imposing a pressure gradient under ear waves, such as hyperbolic wave [10], solitary wave [11][12][13], saw-tooth wav cnoidal wave [15]. In addition, the same Fortran code that was developed and by our group [12,15,16] was utilized in this study to ensure the reliability of results.
The system is closed using the standard two-equation k-ω model propos and further developed by [18]. The transport equations for turbulent kinetic ene

Numerical Analysis Using the k-ω Turbulence Model
The present study used the same k-ω model as in the previous study by [2,9]. For the 1-D incompressible unsteady flow, the equation of motion within the boundary layer can be expressed as follows: where u is the horizontal velocity in the boundary layer, t is the time, p is the pressure, ρ is the fluid density, τ is the shear stress, x is the horizontal coordinate, and z is the vertical coordinate taken as positive upward from the bottom surface. At the location outside of the boundary layer, the pressure gradient term is equal to the temporal variation free-stream velocity, U(t), obtained from Equation (1). In addition, the shear stress for turbulence flow is calculated by introducing the eddy viscosity model. Hence, Equation (2) can be rewritten as follows: where ν is the molecular kinematic viscosity, and ν t is the eddy viscosity. Although the nonlinear advection term is not considered in Equation (3), this linearized boundary layer equation has been successfully employed and validated even for nonlinear wave-induced boundary layers by imposing a pressure gradient under the nonlinear waves, such as hyperbolic wave [10], solitary wave [11][12][13], saw-tooth wave [14], and cnoidal wave [15]. In addition, the same Fortran code that was developed and validated by our group [12,15,16] was utilized in this study to ensure the reliability of the model results.
The system is closed using the standard two-equation k-ω model proposed by [17] and further developed by [18]. The transport equations for turbulent kinetic energy k and the specific dissipation rate of the turbulent kinetic energy ω can be expressed as follows: Turbulent kinetic energy (k): ∂k ∂t Specific dissipation rate (ω): From k and ω, the eddy viscosity can be calculated as where the values of the model constants are given as σ kω = 0.5 and β * = 0.09, which describes the nonlinear transition of a flow regime from a laminar flow to a turbulent flow in Menter's model, σ ω = 0.5, γ = 0.553, and β = 0.075 respectively, and F 1 is a blending function, given as: where arg 1 = min max( in which d 1 is the distance of the cell to the nearest surface, and CD kω is the positive portion of the cross-diffusion term of Equation (5) defined as: Thus, Equations (3)-(5) were solved simultaneously using the free stream velocity, U, angular frequency, σ, kinematics viscosity, ν and water depth, h. An implicit finitedifference scheme with exponentially increasing grid spacing is used to solve the governing equations, allowing for more precise numerical analysis near the bottom [2]. It is determined for a detailed velocity distribution near the bottom by making a mesh interval exponentially increase from the surface of a wall to a numerical computation using an implicit scheme finite difference method. One period is divided into 150,000 as a temporal decomposition, and the magnitude of the mesh is made to increase in the vertical direction in a geometric series towards the upper part so that about ten calculating points may be arranged in a viscous sublayer. The total number of vertical layers is 101. Because the details of the calculation methodology are explained by [16], the detailed explanations are not described here. The uniform median bed sediment is set to d = 0.3 mm, and k s = 2d is the Nikuradse equivalent roughness estimation [19].
In the present study, the bottom shear stress results from the k-ω model are assumed as the true values, and other proposed methods will be compared with the k-ω model results to discuss the accuracy.

Transition from Laminar to Turbulence in the Tsunami Boundary Layer
In most of the current conventional tsunami numerical models, the tsunami wave is assumed as a longwave, and the shallow water equation can be applied. In this theory, the flow regime is classified as a turbulent flow regime, and the boundary layer thickness is developed up to the water surface. However, [6,7] found that the boundary layer near the tsunami source area is surprisingly in the laminar flow regime. Therefore, the use of the steady flow theory may not be applicable for the tsunami case. In addition, the development of a flow regime under a tsunami wave has not been studied before. In this study, the k-ω model was applied to provide detailed time variations of vertical velocity profiles.
As mentioned in [20]'s paper, Kondo [21] developed the laminar solution of the oscillating flow boundary layer under asymmetric flow velocity waveform using Laplace transformation. Although the details are omitted here, the theoretical solution is written as follows: )U(ξ)dξ (10) where u(t, z): mainstream flow velocity, U(t) is free-stream velocity, ν: kinematic viscosity coefficient, t: time; ξ is dummy integration variable. It is possible to obtain an arbitrary flow velocity waveform by numerically integrating the above equation. Figure 3 shows the comparison results between the k-ω model and Kondo's analytical laminar solution under the No.4 GPS tsunami wave at selected phases. In order to quantitatively analyze the difference between the k-ω model and the analytical solution, the periodic averaged deviation, ε, which was introduced by [22], is utilized and is written as follows: where the subscripts (lam) and (kω) denote the analytical and k-ω model, respectively; i is the time index in a wave; U m is the maximum velocity at the wave crest; M = 101 is the number of vertical layers.
as follows: (10) where ( , ) u t z : mainstream flow velocity, () Ut is free-stream velocity,  : kinematic viscosity coefficient, t: time;  is dummy integration variable. It is possible to obtain an arbitrary flow velocity waveform by numerically integrating the above equation. Figure 3 shows the comparison results between the k-ω model and Kondo's analytical laminar solution under the No.4 GPS tsunami wave at selected phases. In order to quantitatively analyze the difference between the k-ω model and the analytical solution, the periodic averaged deviation, ε, which was introduced by [22], is utilized and is written as follows: where the subscripts (lam) and (kω) denote the analytical and k-ω model, respectively; is the time index in a wave; is the maximum velocity at the wave crest; = 101 is the number of vertical layers. Figure 4 shows the time variation in the deviation difference between the k-ω model and the analytical laminar solution. The k-ω model result perfectly agrees with the laminar flow regime condition solution within 100 s of the initial stage of tsunami arrival, indicating the flow regime is laminar flow. After that, the k-ω model velocity profile tends to deviate compared to the laminar solution in the boundary layer. The flow transforms to a turbulent flow regime in an actual tsunami condition.    For determining the transition to turbulence, the spatial and temporal kinetic energy, can be computed by the k-ω model, and is shown in Figure 5. The average value of turbulent kinetic energy is normalized by the quadratic of the amplitude velocity under the crest, . Figure 6 shows the time variation in the total molecular viscosity, ν, and eddy viscosity, ν , at several elevations near the bottom. The computed values of and total viscosity remained nearly constant until the suddenly increased transition point to turbulence around 100 s.   For determining the transition to turbulence, the spatial and temporal kinetic energy, k can be computed by the k-ω model, and is shown in Figure 5. The average value of turbulent kinetic energy is normalized by the quadratic of the amplitude velocity under the crest, U m . Figure 6 shows the time variation in the total molecular viscosity, ν, and eddy viscosity, ν t , at several elevations near the bottom. The computed values of k and total viscosity remained nearly constant until the suddenly increased transition point to turbulence around 100 s. For determining the transition to turbulence, the spatial and temporal kinetic ene can be computed by the k-ω model, and is shown in Figure 5. The average valu turbulent kinetic energy is normalized by the quadratic of the amplitude velocity un the crest, . Figure 6 shows the time variation in the total molecular viscosity, ν, eddy viscosity, ν , at several elevations near the bottom. The computed values of total viscosity remained nearly constant until the suddenly increased transition poin turbulence around 100 s.  The transition from laminar to turbulence conditions can also be observed by analyzing the bottom boundary layer thickness development. Using the dimensional analysis for the partial derivative of Equation (2)     The transition from laminar to turbulence conditions can also be observed by analyzing the bottom boundary layer thickness development. Using the dimensional analysis for the partial derivative of Equation (2) for a steady flow, many studies have shown that the boundary layer under a laminar condition is proportional to the squared root of kinematic viscosity, ν, and time as δ = √ νt. Under a turbulent flow, the thickness of a boundary layer development can be estimated as δ = √ ν t t. Figure 7 shows the log-plot analysis results of the bottom boundary layer thickness development from both the k-ω model and the analytical solution. Because of the absence of overshooting velocity during the accelerating phase, the bottom boundary layer thickness is consistently defined as the distance away from the bed surface where the velocity reaches 99% of the positive free-stream velocity and 101% of the negative free-stream velocity. During the first 100 s, the boundary layer thickness develops in an exponential of 1/2 with time, and thus confirms the flow is in laminar flow. The result of the boundary layer thickness by the k-ω model when transitioning to turbulent flow shows that it is very similar to the squared function of the time. Based on the above formula for turbulent boundary layer thickness, in order to obtain this quadratic relationship, the eddy viscosity ν t should be proportional to the cubed function of the time in the turbulent flow. The transition from laminar to turbulence conditions can also be observed by analyzing the bottom boundary layer thickness development. Using the dimensional analysis for the partial derivative of Equation (2)      As discussed above, the velocity profile from the k-ω model result excellently agrees with the analytical laminar solution within 100 s of the initial stage of tsunami arrival, indicating the flow regime is in a laminar flow. Hence, the k-ω model velocity profile tends to deviate from the laminar solution in the boundary layer to transform to a turbulent flow regime in an actual tsunami condition. Although the transition from laminar flow to turbulence in the tsunami boundary layer exists, its duration is negligible.

Inspection of the Depth-Limitation Condition
The thickness of the wave boundary layer is calculated using the full range equation [23,24]. However, to apply this equation to non-sinusoidal waves such as the surveyed tsunami waveform, it is necessary to assume a certain equivalence.
The flow velocity waveform calculated by Equation (1) is shown in Figure 2. The maximum flow velocity at the outer edge of the wave boundary layer, U m , is estimated from the first wave; and the representing length in a wave boundary layer is not the water depth, h, but the maximum water particle excursion length, a m . In the case of sinusoidal waves, the following equation can be applied: where σ: the angular frequency. However, Equation (12) cannot be used in this investigation. Therefore, a method developed by [25] was utilized. The a m is calculated by integrating the positive flow velocity of the first wave, and a m is computed by dividing the integration by 2. The formula is written as.
where t 0 is the time of the previous zero crossing of the free stream velocity. Using Equation (13), the dimensionless boundary layer thickness (δ/h) can be computed, and the result is shown in Figure 1b. According to this, the ratio δ/h ranges from 0.011 to about 0.034, which means the boundary layer does not span the entire water surface for all six cases. Therefore, applying the steady friction law, such as the Manning formula, is invalid.

Inspection of the Transition to Turbulence
In this section of the study, we examined the turbulent flow transition in a wave boundary layer. According to the investigation of the hypothetical tsunami-induced flow regime of the bottom boundary layer by [7], it has been reported that a boundary layer is in a laminar flow regime near the tsunami source area. For this reason, it is crucial to clarify the location where the transition flow happened.
As previously mentioned, the relevant length in a wave boundary layer is the maximum water particle excursion length a m . Therefore, the Reynolds number is calculated by the following equation: The result is shown in Figure 1c. According to [26,27], it has been reported that the critical Reynolds number that changes to a turbulent flow is Re = 2.5 × 10 5 . It was found that all cases are larger than this critical value, and the Reynolds number shown in Figure 1c is in a turbulent flow region. Figure 8 shows the flow regime for each measuring GPS station. Here, the flow regime classification of [23] was used. According to this, all the cases belong in the transition region from a smooth turbulent flow to a rough turbulent flow.  Figure 8 shows the flow regime for each measuring GPS station. Here, the f gime classification of [23] was used. According to this, all the cases belong in the tra region from a smooth turbulent flow to a rough turbulent flow.

Inspection of the Transitional Characteristics of the Boundary Layer Thickne Friction Factor
Using the full range equation of [24], the dimensionless boundary layer th (δ/h) and wave friction factor can be estimated and are illustrated in Figures 9 respectively. It is again confirmed that all six cases are located in the turbulent tran region from a smooth turbulent flow to a rough turbulent flow. Therefore, the jud of the flow regime is important in the boundary layer under a tsunami. In addition the actual tsunamis, the friction coefficient of a transition region has a small valu pared with the friction coefficient of a rough turbulent flow.

Inspection of the Transitional Characteristics of the Boundary Layer Thickness and Friction Factor
Using the full range equation of [24], the dimensionless boundary layer thickness (δ/h) and wave friction factor can be estimated and are illustrated in Figures 9 and 10, respectively. It is again confirmed that all six cases are located in the turbulent transitional region from a smooth turbulent flow to a rough turbulent flow. Therefore, the judgment of the flow regime is important in the boundary layer under a tsunami. In addition, under the actual tsunamis, the friction coefficient of a transition region has a small value compared with the friction coefficient of a rough turbulent flow. Figure 8 shows the flow regime for each measuring GPS station. Here, the f gime classification of [23] was used. According to this, all the cases belong in the tra region from a smooth turbulent flow to a rough turbulent flow.

Inspection of the Transitional Characteristics of the Boundary Layer Thickn Friction Factor
Using the full range equation of [24], the dimensionless boundary layer th (δ/h) and wave friction factor can be estimated and are illustrated in Figures 9 respectively. It is again confirmed that all six cases are located in the turbulent tran region from a smooth turbulent flow to a rough turbulent flow. Therefore, the ju of the flow regime is important in the boundary layer under a tsunami. In addition the actual tsunamis, the friction coefficient of a transition region has a small valu pared with the friction coefficient of a rough turbulent flow.

Tsunami-Induced Velocity
The flow velocity distribution obtained by numerical computation is shown in 11. The dimensionless vertical coordinates and the flow velocity are normalized water depth, ℎ, and the maximum flow velocity . The acceleration phase in Fig consists of a uniform flow above the boundary layer and logarithmic distributio the boundary layer. The development characteristic of the wave boundary layer is to that of an infinite flat plane. However, the velocity distribution of the acceleratio differs significantly from the well-known periodic wave boundary layer (e.g., [28, stead, it is similar to the velocity distribution of a monotonous boundary layer.
Moreover, δ/h = 0.05, and the log-law distribution is not developed for the wa face like in an open channel and a large velocity gradient near the bottom. It is cl the friction law of steady flows, such as the Manning formula, cannot be used.
The overshooting of velocity is remarkable at t = 966 s during the decelerating as shown in Figure 11b. As mentioned above, the characteristic of the velocity dist ranging from the acceleration phase to the decelerating phase is closely similar to perimental result of the velocity distribution in a boundary layer under a solitar [30,31].

Tsunami-Induced Velocity
The flow velocity distribution obtained by numerical computation is shown in Figure 11. The dimensionless vertical coordinates and the flow velocity are normalized by the water depth, h, and the maximum flow velocity U m . The acceleration phase in Figure 11a consists of a uniform flow above the boundary layer and logarithmic distribution inside the boundary layer. The development characteristic of the wave boundary layer is similar to that of an infinite flat plane. However, the velocity distribution of the acceleration phase differs significantly from the well-known periodic wave boundary layer (e.g., [28,29]). Instead, it is similar to the velocity distribution of a monotonous boundary layer.

Tsunami-Induced Velocity
The flow velocity distribution obtained by numerical computation is shown in Figure  11. The dimensionless vertical coordinates and the flow velocity are normalized by the water depth, ℎ, and the maximum flow velocity . The acceleration phase in Figure 11a consists of a uniform flow above the boundary layer and logarithmic distribution inside the boundary layer. The development characteristic of the wave boundary layer is similar to that of an infinite flat plane. However, the velocity distribution of the acceleration phase differs significantly from the well-known periodic wave boundary layer (e.g., [28,29]). Instead, it is similar to the velocity distribution of a monotonous boundary layer.
Moreover, δ/h = 0.05, and the log-law distribution is not developed for the water surface like in an open channel and a large velocity gradient near the bottom. It is clear that the friction law of steady flows, such as the Manning formula, cannot be used.
The overshooting of velocity is remarkable at t = 966 s during the decelerating phase, as shown in Figure 11b. As mentioned above, the characteristic of the velocity distribution ranging from the acceleration phase to the decelerating phase is closely similar to the experimental result of the velocity distribution in a boundary layer under a solitary wave [30,31].  Moreover, δ/h = 0.05, and the log-law distribution is not developed for the water surface like in an open channel and a large velocity gradient near the bottom. It is clear that the friction law of steady flows, such as the Manning formula, cannot be used.
The overshooting of velocity is remarkable at t = 966 s during the decelerating phase, as shown in Figure 11b. As mentioned above, the characteristic of the velocity distribution ranging from the acceleration phase to the decelerating phase is closely similar to the experimental result of the velocity distribution in a boundary layer under a solitary wave [30,31]. Figure 12 shows the time variation in the bottom shear stress obtained using the k-ω model at three measuring stations (No.2,No.4,and No.5) among the six stations in Table 1. The methods used in the figures are discussed in a later section. From t = 0 s to t = 520 s, the increase in shear stress is more rapid than the flow velocity. 1. The methods used in the figures are discussed in a later section. From t = 0 s to t = 520 s, the increase in shear stress is more rapid than the flow velocity. Figure 13 shows the relationship between the instantaneous bottom shear stress 0 ( ) and the boundary layer outer edge flow velocity ( ). The result confirms that the shear stress is proportional to the squared flow velocity. Therefore, a loop-like relationship is observed, as shown in Figure 13. However, the conventional expression for accurately estimating the bottom shear stress is required for practical application instead of the k-ω model. Therefore, we propose several simple methods in the next section.    Figure 13 shows the relationship between the instantaneous bottom shear stress τ 0 (t) and the boundary layer outer edge flow velocity U(t). The result confirms that the shear stress is proportional to the squared flow velocity. Therefore, a loop-like relationship is observed, as shown in Figure 13. However, the conventional expression for accurately estimating the bottom shear stress is required for practical application instead of the k-ω model. Therefore, we propose several simple methods in the next section. 1. The methods used in the figures are discussed in a later section. From t = 0 s to t = 520 s, the increase in shear stress is more rapid than the flow velocity. Figure 13 shows the relationship between the instantaneous bottom shear stress 0 ( ) and the boundary layer outer edge flow velocity ( ). The result confirms that the shear stress is proportional to the squared flow velocity. Therefore, a loop-like relationship is observed, as shown in Figure 13. However, the conventional expression for accurately estimating the bottom shear stress is required for practical application instead of the k-ω model. Therefore, we propose several simple methods in the next section.

Simple Methods of Calculating Bottom Shear Stress
As discussed in previous sections, the tsunami-induced bottom shear stress needs to be calculated using the wave friction theory instead of the steady friction theory. However, all wave friction theories [11][12][13][14][15][16] are developed based on simple waveforms such as sinusoidal, solitary, and cnoidal waves, and a few are based on an irregular wave [23,32,33]. Therefore, we intended to develop simple methods to apply the wave friction factor, such as the full range equation by Tanaka and Thu (1994), in practice without using a complex model such as the k-ω model. One of the difficulties is how to determine the value of a m . In this study, three methods are proposed to calculate a m , as shown in Table 2. Method 1 uses Equation (13) for each wave crest phase and wave trough phase, then uses a m1 , a m2 , . . . to calculate f w1 , f w2 , . . . for each phase. Method 2 applies Equation (13) for the first maximum positive velocity phase to obtain a m1 and f w1 , then assigns these values to the subsequent phases. Method 3 calculates a m1 and f w1 based on the sinusoidal wave formula for the first wave, and these values are applied to the subsequent phases. Method 4 uses the conventional steady flow friction f c theory. Table 2. Calculation method of bottom shear stress.

Method
Explanation k-ω model (truth) Numerical solution using the turbulence model.

Method 1
Using a m1 , a m2 , . . . for each crest phase and trough phase from numerical integration of the velocity, calculate f w1 , f w2 , . . . for each.

Method 2
Using a m1 from the first positive wave, calculate f w1 , and use this single value for subsequent waves.

Method 3
Using the U m and T 1 (duration of the first positive wave), calculate a m1 = U m T 1 /2π to obtain f w1 , and use this single value for subsequent waves.

Method 4
Steady flow friction coefficient f c (from log law).
Finally, the bottom shear stress results from the k-ω model are assumed to be the true values in the present study, and other proposed methods are compared with the k-ω model results to discuss the accuracy. Figures 12 and 14 compare these four methods' tsunami-induced bottom shear stress results with the k-ω model. As can be seen, Methods 1, 2, and 3 give good results in the first positive wave phase. However, they are underestimated compared to the k-ω model results in the subsequent negative wave phase. Nonetheless, Method 1 still provides better comparison results. Method 4 always gives remarkably underestimated results compared to the k-ω model in all phases.
For a more quantitative assessment of the methods' accuracy, the root mean squared error (RMSE) was calculated for the trough and crest flow velocity phases, as shown in Figure 15. The underestimation in the subsequent trough flow velocity phase is more remarkable even with the wave friction coefficient methods. The significant difference is due to the limitation of the computation method that uses the knowledge of the sinusoidal waves. Therefore, further investigation is needed to consider the wave-like characteristics of such non-sinusoidal waveforms.

Conclusions
The characteristics of the bottom boundary layer under a tsunami were based on the surveyed tsunami waveform. Due to the irregularity of the tsun form, the bottom shear stress could not be estimated accurately using the kn sinusoidal waves. In this study, we proposed simple methods for calculating bo stress for a tsunami with irregular waveforms. The new methods are comparab k-ω model results. In the future, the effect of the history of the reflecting wave should be considered in order to increase the accuracy of the tsunami-induc shear stress estimation.

Conclusions
The characteristics of the bottom boundary layer under a tsunam based on the surveyed tsunami waveform. Due to the irregularity of th form, the bottom shear stress could not be estimated accurately using t sinusoidal waves. In this study, we proposed simple methods for calculat stress for a tsunami with irregular waveforms. The new methods are com k-ω model results. In the future, the effect of the history of the reflecting should be considered in order to increase the accuracy of the tsunami shear stress estimation.

Conclusions
The characteristics of the bottom boundary layer under a tsunami were examined based on the surveyed tsunami waveform. Due to the irregularity of the tsunami waveform, the bottom shear stress could not be estimated accurately using the knowledge of sinusoidal waves. In this study, we proposed simple methods for calculating bottom shear stress for a tsunami with irregular waveforms. The new methods are comparable with the k-ω model results. In the future, the effect of the history of the reflecting waveform shape should be considered in order to increase the accuracy of the tsunami-induced bottom shear stress estimation.
Author Contributions: N.X.T. and H.T. wrote the original draft of the manuscript, and undertook conceptualization, methodology, and interpretation of the results. X.Y. and G.L. performed the data analysis and numerical simulation as well as preparation of the figures. All authors have read and agreed to the published version of the manuscript.