Tsunami Intrusion and River Ice Movement

Abstract

A two-dimensional wave model coupled with ice dynamics is developed to evaluate ice effects on shallow water wave propagation on a beach and in a channel. The nonlinear Boussinesq equations with ice effects are derived and solved by the hybrid technique of the Godunov-type finite volume method and finite difference method with the third-order Runge–Kutta method for time integration. The shock capturing method enables the model to simulate complex flows over irregular topography. The model is capable of simulating wave propagations accurately, including non-hydrostatic water pressure and wave dispersions. The ice dynamic module utilizes a Lagrangian discrete parcel method, based on smoothed particle hydrodynamics. The Boussinesq wave model is validated with an analytical solution of water surface oscillation in a parabolic container, an analytical solitary wave propagation in a flat channel, and experimental data on tsunami wave propagations. The validated model is then applied to investigate the interaction between ice and tsunami wave propagation, in terms of ice attenuation on tsunami wave propagations over a beach, ice deposition on the beach driven by the tsunami wave, and ice jam formation and release in a coastal channel with the intrusion of the tsunami wave. The simulated results demonstrated the interactions between tsunami waves and surface ice, including the maximum run up, ice movement along the beach, and ice jamming in a channel.

1. Introduction

Nonlinear shallow water wave equations are commonly used in global tsunami wave simulations and forecasting. However, they only account for nonlinear properties of the wave and neglect dispersions, which are essential in wave transformation in shallow water. Boussinesq-type equations (BE) considering non-hydrostatic water pressure and wave dispersions are commonly used by coastal engineers.

The classical, standard BE applies for weakly nonlinear and weakly dispersive waves, which was derived in terms of integrated horizontal velocities with linear dispersions. It is restricted to shallow waters with $kH<0.75$, where $k$ is the wavenumber and $H$ is the water depth [1]. By including higher order terms and more accurate vertical velocity distributions, the BE has been improved so that the linear dispersion can be extended to deep water conditions [2,3,4,5,6,7,8,9,10]. Recently, new forms of BE have been developed with higher-order accuracy, which extends the applicability ranging from shallow to deep waters [11,12]. BE is usually modified by adding dissipation terms or the artificial viscosity in momentum equations to account for energy dissipation of wave breaking in the surf zone [13].

The extreme weather might significantly affect wave patterns, flow conditions, sediment transport, and bathymetric changes in coastal region-related climate changes [14,15,16]. The tsunami wave, as an example, was found to interact with sea ice flexure and trigger the Sulzberger Ice Shelf calving in the Honshu earthquake that occurred on 11 March 2011 [17]. The Honshu tsunami wave produced wave run-up heights of 5 m and structure damages from the floating sea ice with thickness of about 0.3 m in the Kuril Islands were observed [18]. Although a tsunami accompanied by floating ice is destructive, the interaction between tsunami wave propagation and floating ice movement has not been quantitatively studied. However, significant progress has been made in ice dynamic studies, especially in terms of mathematical modeling. Shen et al. [19] developed a one-dimensional river ice model to examine the leading-edge propagation for ice jam and potential flooding risks. Later, the two-dimensional ice dynamic model was developed to explore the interactions between unsteady flows and ice movement [20,21,22,23]. Such techniques contributed to the understanding of ice jamming formation and release, and ice cover breakup [20,21,22,23].

In the literature, the BE is rarely applied to investigate wave propagations with ice effects, while the ice dynamic models have not been coupled with the tsunami wave. The main objective of this study is to investigate the interaction between surface ice transport and tsunami wave propagation. This paper first proposes a two-dimensional Boussinesq wave model that is coupled with ice dynamics. Then, the model is validated with analytical solutions for oscillatory flows, theoretical solitary wave, and an experimental tsunami wave. The model is finally applied to simulate tsunami wave propagation with ice, including surface ice transport over a beach, ice deposition on the beach, and ice jam in a coastal channel. Model results showed that surface ice can significantly affect wave propagations along the beach and tsunami waves and vice versa.

2. Model Formulation

A two-dimensional Boussinesq wave model is coupled with the ice dynamic model to study the interactions between ice movement and tsunami wave propagation. The wave model is solved first with given surface ice thickness, areal ice concentration, and ice velocity, based on the extended Boussinesq equations considering ice mass and ice frictions on the flow dynamic. At each coupling time (every 10 min, for example), hydrodynamic parameters, including flow velocity fields, and water surface elevations solved from the wave model are passed to the ice dynamic model. Then, the ice model calculates the ice movement in the coupling time with the given flow velocity and water surface level. Afterward, the calculated ice thickness, concentration, and ice velocity are fed back to the wave model to consider ice influence on mass and momentum conservations. Such processes are repeated until the simulation is completed. Since the smoothed particle hydrodynamic is used to solve the ice dynamics, the ice properties in the Lagrangian field is interpolated to the Eulerian flow field, using the inverse distance weighted method. The coupling processes are the same as the two-dimensional river ice model [20].

2.1. Hydrodynamics Equations

Based on the two-dimensional extended Boussinesq equations from Madsen and Sørensen [24] and the consideration of ice effects on continuity and momentum equations proposed by Shen et al. [20], the conservative form of Boussinesq equations can be expressed as:

$$\frac{\partial \mathit{U}}{\partial t}+\frac{\partial \mathit{F}}{\partial x}+\frac{\partial \mathit{G}}{\partial y}=\mathit{S}$$

(1)

With $\mathit{U}$ denoting vector variable, $\mathit{F}$ and $\mathit{G}$ components of fluxes, and $\mathit{S}$ the source vector, as detailed in the following:

$$\mathit{U}=\left(\begin{array}{c}H\\ P\\ Q\end{array}\right);\mathit{F}=\left(\begin{array}{c}{q}_{tx}\\ \frac{q_{tx}^2}{H_t}+\frac{g}{2}\left(z_s^2-2z_s{z}_b\right)\\ q_{tx}{q}_{ty}/H_t\end{array}\right);\mathit{G}=\left(\begin{array}{c}{q}_{ty}\\ q_{tx}{q}_{ty}/H_t\\ \frac{q_{ty}^2}{H_t}+\frac{g}{2}\left(z_s^2-2z_s{z}_b\right)\end{array}\right)$$

(2)

In Equation (2), $P, Q$ are the reconstructed flow discharge components, following the technique from Tonelli and Petti [25]:

$$P=q_{tx}-\left(\frac{1}{3}+B\right)H^2\left(\frac{{\partial }^2{q}_{tx}}{\partial x^2}+\frac{{\partial }^2{q}_{ty}}{\partial x\partial y}\right)-H\frac{\partial H}{\partial x}\left(\frac{1}{3}\frac{\partial q_{tx}}{\partial x}+\frac{1}{6}\frac{\partial q_{ty}}{\partial y}\right)-\frac{1}{6}H\frac{\partial H}{\partial y}\frac{\partial q_{tx}}{\partial x}$$

(3)

$$Q=q_{ty}-\left(\frac{1}{3}+B\right)H^2\left(\frac{{\partial }^2{q}_{ty}}{\partial y^2}+\frac{{\partial }^2{q}_{tx}}{\partial x\partial y}\right)-H\frac{\partial H}{\partial y}\left(\frac{1}{3}\frac{\partial q_{ty}}{\partial y}+\frac{1}{6}\frac{\partial q_{tx}}{\partial x}\right)-\frac{1}{6}H\frac{\partial H}{\partial x}\frac{\partial q_{ty}}{\partial y}$$

(4)

The source term $\mathit{S}$ is given as:

$$\mathit{S}=\left(\begin{array}{c}\frac{\partial }{\partial t}\left({{Nt}_i}^{\prime }\right)\\ \frac{1}{\rho }\left[\frac{\partial T_{xx}}{\partial x}+\frac{\partial T_{yx}}{\partial y}+{\tau }_{sx}-{\tau }_{bx}\right]-g{z}_s\frac{\partial z_b}{\partial x}+Bg{H}^3\left(\frac{{\partial }^3\eta }{\partial x^3}+\frac{{\partial }^3\eta }{\partial x\partial y^2}\right)\\ +H\frac{\partial H}{\partial x}\left(2BgH\frac{{\partial }^2\eta }{\partial x\partial x}+BgH\frac{{\partial }^2\eta }{\partial y\partial y}\right)+Bg{H}^2\frac{\partial H}{\partial y} \frac{{\partial }^2\eta }{\partial x\partial y}\\ \frac{1}{\rho }\left[\frac{\partial T_{yy}}{\partial y}+\frac{\partial T_{xy}}{\partial x}+{\tau }_{sy}-{\tau }_{by}\right]-g{z}_s\frac{\partial z_b}{\partial y}+Bg{H}^3\left(\frac{{\partial }^3\eta }{\partial y^3}+\frac{{\partial }^3\eta }{\partial y\partial x^2}\right)\\ +H\frac{\partial H}{\partial y}\left(2BgH\frac{{\partial }^2\eta }{\partial y\partial y}+BgH\frac{{\partial }^2\eta }{\partial x\partial x}\right)+Bg{H}^2\frac{\partial H}{\partial x} \frac{{\partial }^2\eta }{\partial x\partial y}\end{array}\right)$$

(5)

Figure 1 shows a sketch of the wave domain with the floating ice and ice thickness as $t_i$. As indicated by the figure and the above equations, $H=h+\eta$, total water depth; $\eta =$ water depth above the reference level; $h=$ water depth below the reference level; $z_b$ = bed elevation; $z_s$ = water surface elevation; $q_{tx}, q_{ty}$ = components of total unit width water discharge; $t, x, y$ = temporal and spatial coordinates; $N=$ areal ice concentration; $t_i^{\prime }=$ equivalent submerged ice thickness; $H_t=$ an equivalent water depth for the total unit water discharge, $q_t$; $\rho =$ water density; $B$ = dispersion coefficient and 1/15 is suggested by Madsen and Sørensen [24]; $T_{ij}={\epsilon }_{ij}\left(\frac{\partial q_i}{\partial x_j}+\frac{\partial q_j}{\partial x_i}\right)$ is the viscous stress; ${\epsilon }_{ij}$ = eddy viscosity coefficients; ${\tau }_{sx}, {\tau }_{sy}$ = components of surface shear stress, including wind stress and ice shear stress; ${\tau }_{bx},{\tau }_{by}$ = components of bed shear stress; and $g$ = gravitational acceleration. The bed shear stress components are calculated through the empirical equations related to Manning’s coefficient.

Equations (1) to (5) return to the shallow water equations with ice effects when wave dispersion terms are ignored, and they become the BE proposed by Tonelli et al. [25] by ignoring the ice effects. The flux gradient term and source terms are reformulated in a well-balanced form to avoid difficulties in applying to steep topography [26]. In comparison to the version of Boussinesq equations proposed by Tonelli and Petti [25], the main modifications presented in Equation (1) lie in the source terms that consider flow viscosity $T_{ij}$, ice effects on mass conservation ${{Nt}_i}^{\prime }$, and ice friction effects on momentum changes ${\tau }_{sx}, {\tau }_{sy}$.

2.2. Numerical Method

A structured, rectangular grid was implemented for the hydrodynamics, using the hybrid technique combining finite volume method and finite difference method [27]. As suggested by Toro [28] for the Godunov-type scheme, two main steps were involved in estimating advection fluxes: (1) Implementation of a high-order reconstruction method to evaluate variables at cell interfaces through interpolations, and (2) application of approximate Riemann solutions to estimate numerical fluxes through interfaces. The fluxes or convective terms were evaluated by the fourth order MUSCL-TVD scheme [29]. The remaining terms in Equation (2) were discretized using central differences or the fourth-order difference method [12,29]. Unlike the classical advection term and pressure gradient term, Equations (1) to (5) are revised to automatically satisfy flux balance in complicated topography. Based on the sketch in Figure 1, there is $z_s^2-2z_s{z}_b+z_b^2-z_b^2=H^2-z_b^2$ and the splitting surface gradient is balanced in the advection term and source term proven in Equations (6) and (7).

$$\frac{g}{2}\frac{\partial }{\partial x}\left(z_s^2-2z_s{z}_b\right)+g{z}_s\frac{\partial z_b}{\partial x}=\frac{\partial }{\partial x}\left(\frac{g}{2}{H}^2\right)+gH\frac{\partial z_b}{\partial x}$$

(6)

$$\frac{g}{2}\frac{\partial }{\partial y}\left(z_s^2-2z_s{z}_b\right)+g{z}_s\frac{\partial z_b}{\partial y}=\frac{\partial }{\partial y}\left(\frac{g}{2}{H}^2\right)+gH\frac{\partial z_b}{\partial y}$$

(7)

Following Liang et al. [26] wet and dry front treatment, the free water surface at the dry nodes equals that of the nearby wetted nodes, while the numerical solver is well-balanced, even at the wet-dry fronts. A critical flow depth of 0.001 m is used to distinguish wet and dry mesh nodes. Combining the finite volume method for the advection term and the finite difference method for other derivatives, Equation (1) was integrated over a small time-step, utilizing the third-order Runge–Kutta scheme [12]. The Runge–Kutta is an explicit scheme, in which the time step $\mathsf{\Delta }t$ should satisfy the CFL conditions for numerical stability.

$$\mathsf{\Delta }t=C\min \left(\frac{\mathsf{\Delta }x}{\left|u\right|+\sqrt{g{H}_t}}, \frac{\mathsf{\Delta }y}{\left|v\right|+\sqrt{g{H}_t}}\right)$$

(8)

where $\mathsf{\Delta }x, \mathsf{\Delta }y$ = spatial steps in a longitudinal and transverse direction, respectively; $u, v$ = velocity components in the x and y direction, respectively; and $C$ = the courant number assigned as 0.45 for the following simulations. That aside, there is no special treatment for the reflective boundaries and this study mainly focuses on the dispersion effects of a tsunami wave without considering breaking waves. Since a tsunami wave is a long wave, most simulations are completed before the wave is reflected back.

2.3. Ice Dynamic Model

The ice thickness, ice concentration, and ice frictions are first calculated from the ice dynamic model and then applied to the hydrodynamic model in the previous section. The governing equations of the ice model include the ice mass conservation, ice area conservation, and the momentum equations. The ice thickness and surface ice concentration are solved by the ice continuum equations. Detailed equations and numerical solutions are documented in Shen et al. [20] and Shen [30]. The momentum equation for ice dynamics is [20]:

$$M\frac{d\overrightarrow{V}}{dt}=\overrightarrow{R}+\overrightarrow{F_a}+\overrightarrow{F_w}+\overrightarrow{G}$$

(9)

where $M={\rho }_{i}N{t}_i=$ ice mass per unit area; ${\rho }_i=$ ice density; $\overrightarrow{V}=u_i\overrightarrow{i}+v_i\overrightarrow{j}$ ice velocity; $D/Dt=$ material derivative; $\overrightarrow{F_a}=$ wind drag at the air–ice interface; $\overrightarrow{F_w}=$ water drag at the ice–water interface; $\overrightarrow{G}=$ gravitational force due to the water surface slope; and $\overrightarrow{R}=$ internal ice resistance = $\overrightarrow{i}\left[\frac{\partial }{\partial x}\left({\sigma }_{xx}N{t}_i\right)+\frac{\partial }{\partial y}\left({\sigma }_{xy}N{t}_i\right)\right]+\overrightarrow{j}\left[\frac{\partial }{\partial x}\left({\sigma }_{yx}N{t}_i\right)+\frac{\partial }{\partial y}\left({\sigma }_{yy}N{t}_i\right)\right]$. with ${\sigma }_{xx}, {\sigma }_{xy}, {\sigma }_{yx}$ and ${\sigma }_{yy}$ for ice internal stresses at different directions. A viscoelastic-plastic (VEP) model [31] is used to calculate the internal stresses. The ice dynamics are simulated with the Lagrangian discrete-parcel method, based on smoothed particle hydrodynamics (SPH) [22,23,30], capable of modeling surface ice runs, ice jam accumulation, and release processes. The ice dynamic model has been validated with field observations and applied to many rivers worldwide [20,30].

3. Model Validations

The two-dimensional wave model is verified for three test cases: (1) The sloshing of water in a parabolic container is simulated and compared with an analytical solution to test the model capability in handling complex wet and dry boundary changes; (2) the analytical solitary wave is used to validate the model accuracy in mass and wave energy conservation; and (3) experimental tsunami wave propagation cases are used to validate the model accuracy in simulating a tsunami wave in a flume channel.

3.1. Case I: Oscillation in a Parabolic Container

The Boussinesq model is first validated with analytical solutions for water sloshing in a frictional parabolic container to validate the model’s mass conservations. The oscillating water surface in the container is damped by the bed friction as a function of the flow depth and local velocity. The parabolic bed elevation and analytical water surface is the same as in Liang and Borthwick’s study [26]. The numerical domain is 10,000 m long and 1000 m wide. The spatial step is 20 m and the total simulation time was 6000 s. The simulated water surface compared with the analytical solution is shown in Figure 2. The initial water surface difference is more than 10 m between the left and right walls. This difference and oscillatory flow is damped gradually at 2000 s, 4000 s, and 6000 s. The strong agreement between simulated and analytical solutions verifies the capability of the model in dealing with complex topography and wet and dry boundaries. The conservation of water mass during the whole simulation is also confirmed.

3.2. Case II: Conservation of Solitary Wave Propagation

The stability and conservative properties of the Boussinesq wave model is tested by the simulation of solitary wave propagation in a flat, frictionless channel, as in Orszaghova et al. [32]. The channel length is 500 m and width is 80 m. The uniform grid resolution is 5 m. The solitary wave with an amplitude of 0.01 m is assumed at $x=100$ m in the otherwise still water channel with a depth of 1 m. The total simulation time is 80 s. The simulation stopped before the wave could reach the downstream boundary, and there is no special treatment for the reflective boundary. Since the bed is frictionless, the initial solitary wave is expected to propagate downstream with the same profile. Simulated results in Figure 3 show that the wave form is preserved during its propagations at different times of 0, 20, 40, 60, and 80 s. This confirms the accuracy of the model.

3.3. Case III: Tsunami Wave in a Laboratory Flume

Yasuda [33] conducted a set of experiments on tsunami wave propagation in a laboratory flume. The channel was 34 m long and 5 m wide. The bed slope of the last 28 m was 1/250 and the first 6 m was flat, where a target wave height of 0.05 m and period of 5 s was generated by the wave paddle. The still water level was 0.4 m and the wave elevation was recorded at $x=10 \mathrm{m}$, 14 m, 18 m, and 22 m. Yasuda’s experimental run 9 is used to validate the model capability in simulating wave transformations during the propagation. Based on the experimental setup, the downstream boundary is a constant discharge of 0.0 with zero input velocities there. The measured elevation at 10 m is used as the upstream boundary condition and a reflective boundary condition is used at the downstream. The longitudinal grid resolution ranges from 0.02 to 0.1 m and the transverse grid resolution is 0.1 m. The total simulation time is 20 s. Figure 4 shows simulated wave forms under varying channel roughness in comparison with measured data at different gauging points, where the dimensionless wave height $\overline{\eta }$ is the wave height divided by the still water depth. Both the simulated and measured wave propagations show that the wave height increases as it climbs up to the shallower water channel. That is because of wave shoaling. The simulation reproduces the experimental wave propagation over a sloping channel. The single-peak wave splits into two peaks when it climbs up the sloping channel due to wave dispersion. The simulated results agree well with the measured tsunami wave at different times and different locations. It further validates the capability and accuracy of the Boussinesq wave model. That aside, sensitivity analysis is conducted with the model on the Manning’s roughness coefficient $n$ ranging from 0.011, 0.021, to 0.03. The simulated results show that the model is not sensitive to the roughness coefficient in simulating tsunami wave propagations in the flume.

4. Model Applications

Due to the limited data on tsunami waves, the Boussinesq wave model has not been applied to natural tsunami wave studies. More field observations are needed for understanding tsunami wave propagations in the coastal region. With the preceding model validations, the two-dimensional Boussinesq wave model is further applied to investigate ice attenuations on tsunami wave propagations over a beach, ice deposition on the beach driven by the tsunami wave, and ice jamming and release caused by the tsunami wave in a coastal channel. They aim to illustrate the interactions between tsunami wave propagation and ice movement. In the future, flume experiments might be conducted to investigate the solitary wave propagations along a sloping beach with the influence of the ice. Such experiments will provide useful data for understanding the interaction between wave propagation and ice movement, and calibrating or validating numerical models.

4.1. Ice Effect on Tsunami Wave Propagations over a Beach

A 1500 m long beach with a constant slope of 0.002 is used to show the interaction between the tsunami wave and ice. A 2 m still water level is assumed to be the initial sea water level without wave. A surge increasing from 2 m to 2.7 m within 4 h and remaining at 2.7 m afterward is used as the upstream boundary at $x=$ 1500 m. The flow velocities at upstream nodes are calculated from the flow equations. The bed roughness coefficient is 0.025. With the same flow boundary, four separated simulation runs are conducted to show the effect of ice concentration. Along with the incoming surge, surface ice with a thickness of 0.1 m and a concentration (N) of 0.0, 0.2, 0.4, and 0.6, respectively, is issued at the upstream boundary. The downstream boundary at $x=0 \mathrm{m}$ is supposed to be the discharge, but these boundaries nodes are dry and zero velocities are assigned in the dry beach nodes. The downstream boundary is a dry beach and is reflective. No wave-breaking is considered. The spatial grid resolution is 5 m, and the width of the domain is 100 m. The ice roughness increases linearly with its thickness from 0.02 to 0.06. The uniform grid resolution is 1 m and the total simulation time is 5 h. The simulated final stages are shown in Figure 5, which illustrates wave propagation along the beach for various incoming ice concentrations. When the ice concentration increases, the wave moves slower and its run-up distance is shorter. The ice attenuates the wave, as wave energy dissipates faster with the increasing ice concentration. In addition, ice climbs higher on the beach when the ice concentration is lower. That means the wave energy is stronger in pushing ice up on the beach, in terms of lower ice concentration. The existence of a large amount of ice causes more energy loss, due to interactions between waves and ice.

4.2. Ice Deposition on the Beach

Ice deposition on the beach is examined in this case. The upstream boundary conditions are provided at $x=1500 \mathrm{m}$ by a solitary wave increasing from 2.0 m to 2.5 m with a period of 600 s, as shown in the insert of Figure 6. Since the downstream boundary is the dry beach, the wave is reflected back there without considering wave breaking processes, the same as that in Section 4.1. The single solitary wave is used to represent the tsunami wave. Initially, the water is at rest and the velocity is zero in the whole domain. The water surface is covered by ice with an initial concentration of 0.3 and thickness of 0.1 m. Three simulation runs are conducted for ice deposition on the beach affected by different ice-bed frictions, respectively. Three ice-bed friction coefficients (defined as $f_{ib}=\tan {\varphi }_{ib}$ and ${\varphi }_{ib}$ is the ice-bed friction angle), 1.04, 2.08, and 5.20, are examined and their corresponding ice-bed friction angles are 46°, 64°, 79°. The friction force between ice and the beach increases with the growing ice-bed friction coefficient. When the ice-bed friction is larger, ice deposition on the beach is larger but the ice run-up is smaller. For a given ice-bed friction coefficient of 6.04 or ${\varphi }_{ib}$ = 81°, the water surface is initially covered by ice parcels with three ice concentrations, 0.2, 0.4, and 0.6. Using the same solitary wave, the final ice depositions on the beach are shown in Figure 7 for simulated results at 3600 s. The ice top indicates the top elevation of the floating ice, and its bottom location is represented by ice bottom with “Surface”, indicating water surface level in all following figures. When ice concentration increases, the ice deposition area decreases, but the maximum accumulated ice thickness increases. The tsunami wave driven ice deposition is sensitive to the ice-bed friction coefficient and the initial ice concentration.

4.3. Ice jamming in River Channel with Incoming Tsunami Wave

The ice jam dynamics in a rectangular channel of 15 km long and 600 m wide with the intrusion of a tsunami wave is simulated by the model. The bed slope of the first 5 km is 0.003, and the remaining 10 km is 0.0003. The bed roughness coefficient is 0.03. The grid resolutions are 30 m in the longitudinal direction and 40 m in the transverse direction. The upstream boundary conditions are a constant water discharge of 2000 m³/s with surface ice thickness of 0.4 m and concentration of 0.6. The velocities at upstream cross-section nodes are calculated based on the inflow discharge and simulated water depth. The downstream boundary condition is a solitary wave with 4 m height and a 900 s period. No special treatment is used for the reflective boundary as the wave dies out before it reflects back. The initial flow condition and longitudinal ice profiles are shown in Figure 8. The total simulation time is 3600 s.

The ice thickness distribution, velocity profiles, longitudinal water surface, and ice profiles are shown at simulation times of 7.5, 12, 13, 14.3, 18, and 22.5 min in Figure 9. The top plots show the ice thickness accumulations in color from the ice jam and water velocity vectors in the two-dimensional domain. The longitudinal water level, bed elevation, and ice thickness profiles are shown in the lower plots with enlarged plots for the area in the green box. The downstream boundary is shown in the middle right plot with a red line for the intruding wave and green points for wave height at current simulating time. When the intruding wave propagates and meeting flow and ice blocks from upstream, ice thickening appears at the wave front and gradually grows into an ice jam at 12 min. The maximum thickness of the ice jam exceeds 6.5 m and the jam bottom reaches the channel bed. After that, ice jams are partially released and then reformed upstream, due to accumulations of ice with decreasing wave height from 13 to 22.5 min. The ice jam moves further upstream as the wave propagates upstream. During the propagation, the wave attenuates and the size of the ice jam decreases. Finally, the wave dies out and the ice jam is released.

5. Conclusions

A Boussinesq-type wave model coupled with ice dynamics was developed for studying ice effects on tsunami wave propagation on beaches and in channels. The model was first validated with analytical solutions and flume data. It was then used to study the interaction between tsunami waves and surface ice. The study showed that floating ice will slow down the wave propagation and shorten the run-up distance. Ice deposition on the beach is dominated by the ice-bed friction and the availability of ice. Depositions are more likely to occur with higher ice-bed friction. Ice deposition area decreases with the increase in ice concentration, but its maximum thickness increases with ice concentration. When a tsunami wave propagates into a channel with floating ice, an ice jam may develop. A preliminary study was conducted on simulating interactions between surface ice movement and tsunami wave propagations. More quantitative comparisons with field observations and experimental data are needed in the future.