NUMERICAL MODELLING OF TIDES AND TSUNAMI WAVES PROPAGATION Текст научной статьи по специальности «Строительство и архитектура»

doi: 10.18720/MCE.82.21

Numerical modelling of tides and tsunami waves propagation

Численное моделирование приливов и распространения волн цунами

Д-р физ.-мат. наук, гл. науч. сотр. В.И. Климович*,

д-р геол.-мин. наук, гл. науч. сотр.-консультант О.К. Воронков, д-р техн. наук, вед. науч. сотр. В.М. Давиденко,

АО "ВНИИГ им. Б.Е. Веденеева", Санкт-Петербург, Россия канд. техн. наук, доцент С.Н. Иванов,

Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, Россия

Key words: tide; tsunami wave; long wave; Ключевые слова: прилив; волна цунами; numerical algorithm; hydraulic structure; bay; длинная волна; численный алгоритм; shallow water model гидротехническое сооружение; залив; модель

мелкой воды

Abstract. Modelling of tides and tsunami waves propagation is important problem for the hydraulic engineering practice since it allows to predict characteristics of these phenomena (load, values of run up and flooding zones, propagation times to structures and so on) and to make reasonable choice of design solution for coastal structures. Numerical modeling is effective way to simulate tides and tsunami waves propagation. Results of simulations of tides and tsunami waves propagation on the basis of shallow water model are presented and analyzed in the article. It was considered explicit and implicit numerical scheme for solution of shallow water model equations, which were verified on the number of experimental results including experimental data on distribution of dam-break wave propagation on a dry and sloping bottom. It is shown that the numerical data of tides in Persian Gulf are well co-ordinated with the data of tidal tables, both on time and on level of tides in various points. Numerical results of tsunami waves propagation in the Krasheninnikov Bay and Kozmino Bay are presented at various angles of tsunami at open boundary of computed area. Numerical results of tide currents for Krasheninnikov Bay are presented for various velocities of a wind. Developed programs may be used, particular, for numerical calculations of impact characteristics of tsunami waves on structures in coastal and offshore zones according requirements of Russian normative documents.

Аннотация. Моделирование приливно-отливных течений и распространения волн цунами является важной задачей для практики гидротехнического строительства, так как позволяет предсказать характеристики этих явлений (нагрузки, величины заплесков и зон затопления, времена добегания до сооружений и т.п.) и сделать обоснованный выбор конструктивных решений для прибрежных сооружений. Численное моделирование является эффективным инструментом для определения характеристик приливно-отливных явлений и распространения волн цунами. В настоящей работе приведены и проанализированы результаты расчетов приливно-отливных течений и распространения волн цунами на основе численного решения уравнений мелкой воды. Рассмотрены явная и неявная численные схемы решения уравнений мелкой воды, которые были верифицированы на большом количестве экспериментальных результатов, в том числе и на экспериментальных данных по распространению волн прорыва по сухому и наклонному дну. Показано, что данные расчетов по моделированию приливов в Персидском заливе удовлетворительно согласуются с данными приливных таблиц, как по времени наступления приливов, так и по величине приливов в различных точках. Представлены результаты расчетов распространения волн цунами в акваториях бухты Крашенинникова и бухты Козьмино при различных углах подхода волн цунами к открытой границе расчетной области. Приведены результаты расчетов приливно-отливных течений для бухты Крашенинникова при различных заданиях скоростей ветра. Разработанные программные средства могут, в частности, применяться для проведения численных расчетов по определению характеристик воздействия волн цунами на

V.I. Klimovich*, O.K. Voronkov, V.M. Davidenko,

"B.E. Vedeneev VNIIG", JSC, Saint Petersburg,

Russia

S.N. Ivanov,

Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russia

сооружения, расположенные в прибрежной и береговой зоне, согласно требованиям Российских нормативных документов.

1. Introduction

Numerical simulation of tidal currents and tsunami wave propagation is an important problem for the practice, allowing predicting the characteristics of these phenomena (loads, run up, propagation time to structures and so on), creating warning systems for tsunami waves and making a reasonable choice of design solutions for the design of various offshore hydraulic structures. Tidal currents and tsunami waves are traditionally referred to long-wave processes, since the lengths of tidal waves and tsunami waves are estimated at tens of kilometers.

The study of tidal currents and the modeling of tsunami wave propagation have been devoted to a significant number of publications, for example, [1-21]. In [1-4] fundamental theoretical studies on various mathematical models describing long-wave processes, shelf effects for long-wave processes, the interaction of tidal waves and tsunami waves with islands, the energy of ocean tides, and the resonance phenomena of long-wave processes are presented. At present, the model of shallow water has been most widely used for numerical modeling of tides and tsunami waves [7-13], which makes it possible to efficiently calculate marine areas with characteristic dimensions, measured in hundreds and thousands of kilometers, to take into account the movement of the flow on a dry bottom, to simulate flow in multiply connected areas, while providing an acceptable for practice accuracy of the calculation results. Large-particle models (SPH) [14, 15] and three-dimensional models under the assumption of hydrostatic pressure distribution [5, 6] are also used to simulate the tides and propagation of tsunami waves, and even one-dimensional models [16] that allow for warning systems give a rapid assessment of tsunami waves transformation. In [17-21] questions of physical modelling of tsunami waves and their interaction with hydraulic engineering constructions are considered. Different types of tsunami generators were used: pneumatic generators [17, 20, 21], landslide generators [18, 19].

Now it is represented that the most perspective direction for calculation of tidal currents and tsunami wave propagation is using procedures of parallel calculations on the computer together with using processors of graphic cards [22], thus time of calculations can be essentially reduced and consideration of 3D models of a currents for water areas of the big sizes is possible.

In connection with intensive development of the Russian shelf and construction in coastal tsunami dangerous areas (see, for example, [23]) the purpose of this paper was to develop the effective algorithms of numerical simulation of tides and tsunami wave propagation including flow movement on a dry bottom for solving shallow water equations constructed on a conservative explicit and implicit numerical Roe scheme of second-order accuracy in space and time. In addition, the possibility of numerical modeling of tides and the propagation of tsunami waves at different angles of long waves to the open boundary of the calculated region was considered.

2. Methods

2.1. Governing equations

The equations describing unsteady two-dimensional shallow water flow with account of bottom friction, wind action, Coriolis forces and viscosity may be represented in the following conservative form (see, for example, [24, 25]):

dU

~дГ

Here

dF dx

д G dy

dFv dG

v

dx dy

Qs .

(1)

" h ' "Qx " "Qy "

U = Qx , F = Qx2/ h , G = QxQy / h

Qy _ QxQy / h _ Qy/ h _

(2)

Fy =

0 0 q

hTxx Pw G = hTxy Pw , Qs = Qcq т- gh ^ Pw Pw dx

hTyx Pw hTyy Pw ^ Tby Twy , dH ncQx y + y Sh Pw Pw dy

Qx hux, Qy huy ' Tbx PwCcurQx, Tby PwCcurQy, Ccur S

Q2

w cur y

C 2 h 2

C =

h

1/6

n

Twx = Tw COSl°w >

Twy = Tw sin >

Tw = Pa** ln 2

W1

z 0 = 0.0144-

V z0

PaS

The depth-averaged shear stresses Txx, Txy, Tyx, Tyy may be written as [24] (here the indices i, j correspond to x, y):

Pw

= V

d( hu{ ) + d(huj )

dx,

dx,-

(3)

In the equations (1)-(3), h - the water depth; H=h+Zb - the water surface level; Zb - the bottom surface level; Ux, Uy - velocity projections in x, y directions respectively; Vt - turbulent kinematic viscosity; Oc = 2rn sin cp - the Coriolis parameter; g> - Earth angular velocity; cp - latitude of the locality; g - free

fall acceleration; pw - water density; pa - air density; Tbx, Tby - projections of bottom friction force in x, y directions respectively, C - the Chezy coefficient; n - bottom roughness coefficient; Ccur - coefficient of bottom friction; Twx, Twy - projections of the wind tangent stresses in x, y directions respectively, k = 0.4 - the Karman constant; W - wind velocity at height zw relative to the surface elevation; Ow - angle of inclination of the wind velocity vector in x direction; zoo - water surface roughness parameter; q(x,y) - mass source intensity in the point with coordinates (x, y). The wind shear stresses (see (2)) were determined according to [25].

2.2. Numerical algorithms

2.2.1. Explicit scheme

For solving the set of equations (1), select the variables h, Qx, Qy, Has the unknowns, and consider the extended set of equations [26] by formally adding the equation for H to the set (1):

д U

dF dx

д G dy

dFv д G

y

dx dy

Qs .

(4)

Here

U =

~h ' ~QX ' Qy '

Qx , F = Qx2/ h , G = QxQy / h

Qy QxQy / h Qy/ h

H Qx _ Qy _

Fy =

" 0 " " 0 "

hTxx hTxy

Pw hTyx = Pw hTyy

Pw Pw

0 0

, Qs =

QJQ -hbL + 1™ - gh 5H Pw Pw

-acQx -TP +TP Pw Pw

- gh

dx dH

dy

q + -

dZh

dt

It is assumed in (4) that the temporal variation in the bottom geometry in each point of the computation domain is known and described by the function

dZb = f (t, x, y).

dt

In the further analysis the subdivision of the computational domain into the finite volumes (FV) will be considered as an arbitrary quadrangular node grid. The computational curvilinear node grid is formed by assigning the coordinates of the apexes of FV; the coordinates of the computation points are determined in the gravity center of FV.

Using the finite volume method, obtain for the explicit scheme:

dU 4 - - -Scv— + Z F - Fvnk )i k -QsScv = 0 . dt k=i

(5)

Here Fk = Fknx + Gkny, F"k = Fv knx + Gv kny - the convective and the diffusive fluxes

respectively, normal to k-th edge of FV; nx, ny - projections of the external normal to the edge of FV in the x,y directions respectively ; Scv - the finite volume area, ik - the length of k-th edge of FV.

The computation of the convective fluxes Fn at the edge of FV is made using the second order accuracy Roe's scheme over the spatial variables x, y:

Fn = I

2

Fn(UR) + Fn(UL)- A + AJ(Ur-Ul)

Ah =

0 0 0 0

0 0 0 c2nx

0 0 0 c2n

0 0 0 0

c =

4gh.

(6)

Here Ul and Ur are reconstructed values of the variables on the inner and outer side of the edge

of FV respectively; A = dFn/dU - the Jacobian; Ah - additional matrix accounting of the components

h dH/dx, h dH/dy in (4).

The matrix A + Ah may be represented as (see [26, 27]):

^H

"1 0 0 0" "2 0 0 0 " 1 0 1 c2 "

T -1 A*T, A* = RAR -1, T = 0 nx ny 0 , A = 0 2 0 0 , R = un - c 0 un + c 0 . (7)

0 - ny nx 0 0 0 2 0 ul -1 ul ulc 2

0 0 0 1 0 0 0 24 1 0 1 u„2 _

Here T- the matrix of transition from the coordinate system x,y to the local coordinate system n,l at the edge of FV; Xi=Un-c, A,2=Un, A3=Un+c, 2.4=0 - the eigenvalues of the matrix A*; R - the matrix of the

right-side eigenvectors A ; un = uxnx + uyny,

ul = -Uxny +uynx .

The matrix

A + A

h

in (6) has the form |A + AH| = T 1R |A|R 1T , in which the matrix |A| is the diagonal matrix with the modules of the eigenvalues on the diagonal.

Due to the fact that the normal specific discharges Qn can be prescribed at the boundaries of the computation domain (e.g. for the solid boundary Qn=0), the relationship (6) may be conveniently written at the boundary of FV using the variables V = TU in the local coordinate system n,l:

1

Fn = 2

Fn (VR) + Fn (VL )-T-1 AA(Vr-Vl )

(8)

In (8) the reconstruction of variables V (and not U) is making for the boundary of FV.

The reconstruction of the variables V for the edge of FV in case of quadrangular computation grid was performed in the following way [26]. For the computation point with number (ij) considered were the

four nearest adjacent points with numbers (i±1j±1) and the derivatives ax =dV/dx, Pp =dV/dy

(k=1,...,4) were determined by turns for four combinations of point ((ij) with two adjacent points (e.g.

(i,j+1) and (i+1,j). The determination of the slopes ax, Py in FV for the computation point (ij) was made

by applying one of the known limiters (minmod, superbee, van Leer, van Albada [26]) to the computed

values of the slopes a^,pp . The reconstructed value of variable Vl at k-th edge of the inward side of FV

for the computation point (ij) is determined using the relationship (xk,yk- the mid coordinates of the k-th edge):

VL = Vij +ax(xk - xi,j) + py(yk -yi,j).

The reconstructed value of variable VR at k-th edge of the outward side of FV for the computation point (ij) is determined in a similar way.

It should be noted that the computation of the slopes ax, Py of FV adjacent to the boundary of the computation domain is made using the non-reflecting boundary conditions [26].

The computation of the values un, ui, c at the boundary of FV, required for the determination of the elements of matrices R, A, R1 (see (7),(8)), was performed with the aid of Roe's averages with entropy correction.

The approximation of the component h dH/dx was made using the following scheme. From the reconstructed values of hi, hR, Hl, Hr the expressions hL dHL /dx , hR dHR /dx were computed by the central-difference way, and their half sum was assumed as h dH/dx. The derivative h dH/dy was determined in a similar way.

The selection of additional variable H and the extended set of equations (4) makes it possible to obtain a trivial solution for the Roe's scheme using arbitrary curvilinear grids for any bed geometry, i.e. at ux =uy =0 , it follows that H=const.

It should be noted that the above computation scheme of convective flows at the edges of FV also permits to make reconstruction in terms of the velocities un, ui [28, 29], and not the specific discharges Qn, Qi. As noted in [29], the choice of variables un, ui for reconstruction permits to enhance the stability of the algorithm, especially in the presence of dry bed.

The diffusion flow F"k is computed by central difference scheme.

The temporal integration of equation (5) was made with second order accuracy using the predictor-corrector scheme.

The condition of calculation stability for the considered explicit scheme of equations (1) has the form:

Cr =

Ai

min

Ax

f \ Ay

< 1.

J i, j

(9)

7i,j J

Here Cr - Courant number, At - differential step on time, Ax, Ay - differential steps on spatial coordinates x, y.

The condition (9) means that for transition to the following time level (p+1) from current time level (p) the time step At is taking on the basis of above presented relation taking into account the specified parameter Cr <1.

The following procedure was used for making consideration of the dry bed flow. With the condition h < hmin (hmin is selected in the range of 10-6-10-12 m) satisfied in some point (ij), the flow velocities Vx, Vy at that point were assumed zero. During the computations of convective flows at the edge of FV for which the adjacent FV is dry (h < hmin), there was made the analysis of probability of the water entrance into the FV through the edge. In this case, if the water level in the adjacent dry volume exceeds the one of the analysed volume, the condition of solid-wall was imposed ((un)R=-(un)L, (u)r=(u)l, hR=hL,HR=HL). Such approach to the dry bed problem made it possible to perform computations for the bed drying and submersion of the dried bed as well as for multiply connected domains.

2.2.2. Implicit scheme

The integration of equation (4) over FV area gives the following equation for the implicit scheme:

+ ey\F? - Fn di k=i

v,k

p+i

£k + (1-0)Z\Fkn -Fvnk k=i

U +

0QP+1 + (i- 0)QP Scv = o. (io)

At 0=1 we receive the Euler implicit scheme of the first order accuracy with respect to time; at 0=1/2 we receive the Crank-Nickols implicit scheme of the second order accuracy with respect to time.

Equation (10) linearized with respect to time with the second order accuracy of temporal step At may be represented as:

AU,

Ai

bJL + 0LS AUi , +0LS AU = RHS, P1 i,j S2 i'j

Ls1AUitj = Bp1 AU— j + CS1AUU j + Dp1AUi+1, j, LpAij = Bp^^ j-1 + CpW^ j + D&2AUU j+1 , (ii)

(

RHS = -

— I

S

V Scv k=1

Fk - Fv,k

^ k - Qp

Here AUi = Up,+1 - Up- difference of variables vector in a point i, j of considered FV during

i I 1,1 1,1

time Ai.

The matrices Bp , Cpi, Dp

correspond to the components from the convective and diffusion

fluxes along the opposite edges in one direction, the matrices B^, Cp , Dp^ - in another direction. Matrices Bp , Cpx, Dp, for rectangular coordinates in direction pi = x have the form:

Be =■

Ax\(A + Ah)r-IA + aH 11.-1/2,. + Bv,p,

2Ax

= ^ A + AH )L + 1A + AH Il-"2,У + A + AH )L +1A + AH I]f +1/2,j } +

1QS - Qs)

2 dJ

+ C

vS

1 = 2^[(A + AH)R "IA + AH11+1/2,j + DvSl

DSi =

Bv ,S1 =

Ax2

Ax2

1/2, j

1/2, j

Dv ,S1 =

Ax2

(V )

■ +1/2, j

Ax2

(V)

■ +1/2, j

Cv ,S1 = - Bv,S1 - Dv ,S1, QS =

0

- gh— dx

, dH - gh—

dy

0

form.

The matrices Bs , CSl, Dsin direction S2 = y for rectangular coordinates can be written in similar The right-side part of (11) is determined in the same way as in the explicit scheme.

After partial factorization (11) with the second order accuracy with respect to At, we receive (I -identity matrix):

(I + 6AtLSi ) (I + 0AtLS2) AUJU j = AtRHS.

(12)

Equation (12) can be effectively solved using the Thomas algorithm for block three-diagonal matrices in one and another direction successively [30]. Note that in (12) the values of temporal increments of variables in points located outside the calculation domain were determined with account of non-reflecting boundary conditions.

The choice of time step for the implicit scheme does not require implementation of condition (9). However, keeping in mind that the stability analysis of difference schemes is spent on a basis of linearized equations, in practice the Courant number, at which stability of the implicit scheme is observed, is limited and it is defined by selection way. It is necessary to notice that the implicit scheme has the increased stability in comparison to the explicit scheme and Courant number, at which stability of this difference scheme is observed for a specific object, can be much more 1.

Boundary conditions on the open boundary of calculating area for simulations of tides and tsunami wave propagation were taken according [31] and were considered angle of tsunami wave on open boundary.

Verification of examined numerical algorithms was made on the basis of comparison of calculation results with experimental data and with analytical solutions for examples presented in [32]. Let's notice that verification of the presented algorithms included also comparison of numerical solutions with experimental data on distribution of dam-break wave including wave propagation on a dry bottom and sloping bottom.

3. Results and Discussion

3.1. Tidal currents in the Persian Gulf

The computation of tidal currents was made for the Persian Gulf water area. The boundary condition was described as water surface level in Ormuz Strait in the area of Jask Bay (the right-side open boundary in Figure 1) according to the Tide Tables [33] per 02-04 July 1996.

The computations were made with 201x140 mesh, n = 0.012, the Coriolis parameter -Qc = 0.0000685 1/s, wind speed - W = 0 m/s. For the explicit scheme it was assumed that Cr = 0.5, for the implicit scheme - Cr = 1.5-2.5. The computation time was 10 cycles of diurnal tides which provided the periodic regime. Figure 2 presents the results of computations of dynamics of water level in a number of localities on the shore of Persian Gulf under calm conditions for one tidal cycle. The comparison between the computed data on the water level and the available data on the tidal level variations in Dubai, BenderLenge, Shat-Al-Arab [33] indicated their good agreement. In particular, the computed time of approach of the higher high level (HHL) for Shat-Al-Arab with respect to the HHL time for Bender-Lenge (« 12 hours)

practically corresponds to the data of [33]. For Shat-Al-Arab the difference between HHL and lower low level (LLL) was computed as 3.5 m; according to [33], this difference amounted to 3.5 m; the difference between the lower high level (LHL) and higher low level (HLL) was computed as 1.6 m and according to [33] - 1.4 m. Figure 1 shows the vectors of depth-averaged flow velocity under calm conditions at the time corresponding to HHL in the area of Jask Bay.

■200 -100 0 100 200 300 400 ¡00 600 700

— -25 mlit 0 jaiKBay © Sbai-Al-Arab © l>nh;n © Bender-Lcngi 1

Figure 1. Distribution of velocity vectors in Persian Gulf at the time corresponding to HHL in Jask Bay

H, m

Figure 2. Water level tidal dynamics in Persian Gulf (—•--Jask Bay, —a--Shat-Al-Arab,

—♦--Dubai, —□--Bender-Lenge)

3.2.Distribution of tsunami wave in Kozmino Bay

Calculations of the tsunami wave propagation in the Kozmino Bay were made with the direction of the wave from the west. The period of the tsunami wave was set equal to 35 minutes, the wave height at the entrance to the bay was 0.845 m. The calculations were made on a rectangular grid of dimension 343х293 (the difference step on the coordinates was 5 m). The order of accuracy of the difference approximation with respect to spatial coordinates and time was equal to 2. The Courant parameter was chosen equal to 0.4 for the explicit scheme. The roughness coefficient was chosen to be n = 0.035, the

Coriolis parameter - Qc = 0.0000987 1/s, wind speed - W = 0 m/s. The initial conditions were given as H = 0, Vx,y = 0.

In the calculations near the shore of the Kozmino Bay checkpoints were set, in which the dynamics of water surface levels during the tsunami wave propagation were recorded. The coordinates of these checkpoints (coordinated with the coordinate system in Figures 3-4) are presented in Table 1.

Table 1. Coordinates of Kozmino Bay checkpoints

No. of 1 2 3 4 5 6

checkpoints

x, m 770 675 590 820 955 1200

y, m -120 -190 -300 -445 -600 -880

The results of calculating the tsunami wave propagation in the Kozmino Bay with the direction of the wave from the west are shown in Figures 3-6. Comparison of Figure 3 and Figure 4 shows that coastal line changes at various stages of tsunami wave impact. From Figures 5-6 we can calculate the time when tsunami wave can reach checkpoints. For example, for checkpoint 1, which was initially on the shore, the time of tsunami reach was 0.08 hour after tsunami reached the open boundary (inlet of the Bay). From Figures 5-6 we can see that the maximum of shore level where tsunami wave can run up was about +1 m.

02 04 06 0.8 1 12 1.4 16 13 -m- = 0.6 m/s x, km

Figure 3. Distribution of tsunami current vectors in Kozmino Bay for the time corresponding to the minimum tsunami wave level at the entrance of bay (colour background - module of average velocity)

Figure 4. Distribution of tsunami current vectors in Kozmino Bay for the time corresponding to the maximum tsunami wave level at the entrance of bay (colour background - module of average velocity)

i/

'■>A \

*....... *........ i'------- *-----

\

V

\

0 0 36 0 12 O 18 O. 24 3 0.3 S 0. •12 O. 48 O. 54 O

1. h

Figure 5. Water level dynamics during the tsunami wave propagation in the Kozmino Bay

(checkpoints 1 - • , 2 - ■ , 3 - ▲)

w

/

A

O 0 36 0 2 IS O. O 3 O 36 042 O. 43 O 54 06

l.h

Figure 6. Water level dynamics during the tsunami wave propagation in the Kozmino Bay

(checkpoints 4 - •, 5 - ■, 6 - ▲)

3.3. Tidal currents and distribution of tsunami wave in the Krasheninnikov Bay

Calculations for the Krasheninnikov Bay were made on a rectangular grid of dimension 251x226 (the difference step in the coordinates was 40 m). The order of accuracy of the difference approximation with respect to spatial coordinates and time was equal to 2. The Courant parameter was chosen equal to 0.5 for the explicit scheme. The roughness coefficient was set to 0.035, the Coriolis parameter was Qc = 0.000133 1/s.

Calculations of tidal currents in Krasheninnikov Bay were carried out by specifying the water levels at the open border (I-I in Figures 7-9), corresponding to the daily tide. The range of the water levels at the open border was 4.85 m.

■ t шШ^^ШШк^'Ш Н11 in {г 4

■Л--:- ' ^

-■ 0.3 m/s x, km

Figure 7. Distribution of tidal current vectors in Krasheninnikov Bay at the moment of time corresponding to the maximum value of the water surface level at the open boundary I-I (W = 0 m/s)

• Q.15 m/s x, km

Figure 8. Distribution of tidal current vectors in Krasheninnikov Bay at the moment of time corresponding to the minimum value of the water surface level at the open boundary I-I (W = 0 m/s)

Figure 9. Distribution of tidal current vectors in the Krasheninnikov Bay at the moment of time corresponding to the maximum value of the water surface level at the open boundary I-I

(W = 40 m/s)

The results of calculations of tidal currents in Krasheninnikov Bay in the absence of wind (W=0 m/s) are shown in Figures 7-8. Figure 9 shows the results of calculations of tidal currents in the presense of the northeast (NE) wind with a speed W=40 m/s. Comparison of Figure 7 and Figure 8 shows that coastal line changes at various stages of tide. Comparison of Figure 7 and Figure 9 shows that impact of wind speed on current velocity can be conciderable.

Tsunami wave calculations for the Krasheninnikov Bay were carried out by specifying the water levels at the open border (I-I in Figures 10-13), corresponding to wave with the height 6.54 m and period 1800 s. The initial conditions were given as H = const, Vx,y = 0.

Figure 10. The distribution of free water surface levels in Krasheninnikov Bay for the time corresponding to the maximum tsunami wave level at the open boundary I-I

Figure 11. Distribution of tsunami current vectors in Krasheninnikov Bay for the time corresponding to the maximum tsunami wave level at the open boundary I-I

Figure 12. The distribution of free water surface levels in Krasheninnikov Bay for the time corresponding to the minimum tsunami wave level at the open boundary I-I

Figure 13. Distribution of tsunami current vectors in Krasheninnikov Bay for the time corresponding to the minimum tsunami wave level at the open boundary I-I

Figures 10-13 show the results of tsunami wave propagation calculations for different time moments in the absence of wind (W = 0 m/s). Comparison of Figure 11 and Figure 13 shows that coastal line changes at various stages of tsunami wave impact. We can see from Figures 11, 13 that currents in the Krasheninnikov Bay had very complicated character. A considerable vortex arised near the open boundary I-I, where current velocities were about 9-10 m/s. The maximum of shore level where tsunami wave can run up was about +6.5 m (see Figure 12).

4. Conclusions

1. Conservative explicit and implicit numerical schemes Roe of 2nd order of accuracy on space and time for shallow water model are developed.

2. The presented explicit and implicit numerical algorithms make it possible to model both subcritical and supercritical flows in multiply-connected domains taking into account the movement of the flow on a

dry bottom including the computation of tidal currents and propagation of tsunami wave. In a number of cases the implicit algorithm proves more effective compared to the explicit algorithm, since it permits to make computations with larger Courant numbers and requires only one step of integration with respect to time. But in some cases expenses of machine time for implicit scheme can be more than for explicit scheme, especially it concerns cases of application OpenMP algorithms of parallel programming technology.

3. Numerical results of tides and tsunami waves propagation for various angles of long waves on open boundary of computed area and for various wind velocities are presented. The obtained numerical results for tides show good agreement with the data of Tidal Tables both on time and on level of tides in various points. Developed algorithms and programs may be used, particular, for numerical calculations of impact characteristics of tsunami waves on structures in coastal and offshore zones according requirements of [34].

References

1. Yefimov, V.V., Kulikov, E.A., Rabinovich, A.B., Fine, I.V. Volni v pogranichnih oblastiah okeana [Ocean boundary waves]. Leningrad: Gidrometeoizdat, 1985. 280 p. (rus).

2. Marchuk, G.I., Kagan, B.A. Dinamika okeanskih prilivov [Dynamics of Ocean Tides]. Second edition. Leningrad: Gidrometeoizdat, 1991. 471 p. (rus).

3. Nekrasov, A.V. Prilivnie volni v okrainnih moryah [Tide waves in the edge seas]. Leningrad: Gidrometeoizdat, 1975. 247 p. (rus).

4. Shokin, Yu.I., Chubarov, L.B. Ocherk istorii issledovaniya problemi tsunami v Sibirskom otdelenii Rossiyskoy akademii nauk [Brief history of tsunami problem investigation in Siberian Department of Russian Academy of Science] Computational Technologies. 1999. Vol. 4. No. 5. Pp. 70-105. (rus).

5. Sammartino, S., Lafuente, J.G., Garrido, J.C.S., De los Santos, F.J., Fanjul, E.A., Naranjo, C., Bruno, M., Calero, C. A numerical model analysis of the tidal flows in the Bay of Algeciras, Strait of Gibraltar. J. Continental Shelf Research. 2014. Vol. 72. Pp. 34-46.

6. Androsov, A.A. Regionalnie modeli mirovogo okeana -obschiy podhod k modelirovaniyu [Regional model of word ocean - common strategy for modeling]. Abstract of doctoral phys.-math. science thesis. 2010. 46 p. (rus).

7. Borrero, J.C., LeVeque, R.J., Greer, S.D., O'Neill, S., Davis, B.N. Observations and Modelling of Tsunami Currents at the Port of Tauranga, New Zealand. Proc. Australasian Coasts & Ports Conference. 15-18 September 2015. Auckland. New Zealand.

8. Wang, L., Ying, C. Research on astronomical tide and tsunami coupled model. IOP Conf. Series: Earth and Environmental Science. 2018. Vol. 153. 022018.

9. Hill, D.F., Griffiths, S.D., Peltier, W.R, Horton, B.P., Tornqvist T.E. High-resolution numerical modeling of tides in the western Atlantic, Gulf of Mexico, and Caribbean Sea during the Holocene. J. Geophysical Research. 2011. Vol. 116. C10014.

10. Kowalik, Z., Proshutinsky, A. Tsunami-tide interactions: A Cook Inlet case study. J. Continental Shelf Research. 2010. Vol. 30. Pp. 633-642.

11. Androsov, A.A. Nelineinoe vzaimodeistvie voln tsunami i prilivov [Nonlinear interaction of tsunami waves and tides]. Bulletin of civil engineers. 2010. Vol. 23. No. 2. Pp.181-188. (rus).

12. Lavrentiev, M.M., Lysakov, K.F., Marchuk, An.G., Romanenko, A.A Bistroe chislennoe modelirovanie tsunami [Tsunami fast numerical modeling]. Bulletin of cybernetics. 2016. No. 2, Pp. 92-102. (rus).

13. Lavrentiev, M.M., Romanenko, A.A., Lysakov, K.F. Modern Computer Architecture to Speed-Up Calculation of Tsunami Wave Propagation. Proceedings of the Eleventh Pacific/Asia Offshore Mechanics Symposium. Shanghai. China. October 12-16. 2014. Pp. 186-191.

Литература

1. Ефимов В.В., Куликов Е.А., Рабинович А.Б., Файн И.В. Волны в пограничных областях океана. Л.: Гидрометеоиздат, 1985. 280 с.

2. Марчук Г.И., Каган Б.А. Динамика океанских приливов. 2-е изд. Л.: Гидрометеоиздат, 1991. 471 с.

3. Некрасов А.В. Приливные волны в окраинных морях. Л.: Гидрометеоиздат, 1975. 247 с.

4. Шокин Ю.И., Чубаров Л.Б. Очерк истории исследования проблемы цунами в Сибирском отделении Российской академии наук // Вычислительные технологии. 1999. т. 4. № 5. С. 70-105.

5. Sammartino S., Lafuente J.G., Garrido J.C.S., De los Santos F.J., Fanjul E.A., Naranjo C., Bruno M., Calero C. A numerical model analysis of the tidal flows in the Bay of Algeciras, Strait of Gibraltar // J. Continental Shelf Research. 2014. Vol. 72. Pp. 34-46.

6. Андросов А.А. Региональные модели мирового океана - общий подход к моделированию. Автореферат дисс. на соиск. ученой степени д.ф.-м.н. 2010. 46 с.

7. Borrero J.C., LeVeque R.J., Greer S.D., O'Neill S., Davis B.N. Observations and Modelling of Tsunami Currents at the Port of Tauranga, New Zealand // Proc. Australasian Coasts & Ports Conference. 15-18 September 2015. Auckland. New Zealand.

8. Wang L., Ying C. Research on astronomical tide and tsunami coupled model // IOP Conf. Series: Earth and Environmental Science. 2018. Vol. 153. 022018.

9. Hill D.F., Griffiths S.D., Peltier W.R, Horton B.P., Tornqvist T.E. High-resolution numerical modeling of tides in the western Atlantic, Gulf of Mexico, and Caribbean Sea during the Holocene // J. Geophysical Research. 2011. Vol. 116. C10014.

10. Kowalik Z., Proshutinsky A. Tsunami-tide interactions: A Cook Inlet case study // J. Continental Shelf Research. 2010. Vol. 30. Pp. 633-642.

11. Андросов А.А. Нелинейное взаимодействие волн цунами и приливов // Вестник гражданских инженеров. 2010. Т. 23. № 2. C. 181-188.

12. Лаврентьев М.М., Лысаков К.Ф., Марчук Ан.Г., Романенко А.А. Быстрое численное моделирование цунами // Вестник кибернетики. 2016. № 2. C. 92-102.

13. Lavrentiev M.M., Romanenko A.A., Lysakov K.F. Modern Computer Architecture to Speed-Up Calculation of Tsunami Wave Propagation // Proceedings of the Eleventh Pacific/Asia Offshore Mechanics Symposium. Shanghai. China. October 12-16. 2014. Pp. 186-191.

14. Dao M.H., Xu H., Chan E.S., Tkalich P. Modelling of tsunami-like wave run-up, breaking and impact on a vertical wall by SPH method // J. Nat. Hazards Earth Syst. Sci. 2013. Vol. 13. Pp. 3457-3467.

15. Шокин Ю.И., Бейзель С.А., Рычков А.Д., Чубаров Л.Б. Численное моделирование наката волн цунами на побережье с использованием метода крупных частиц //

14. Dao, M.H., Xu, H., Chan, E.S., Tkalich, P. Modelling of tsunami-like wave run-up, breaking and impact on a vertical wall by SPH method. J. Nat. Hazards Earth Syst. Sci. 2013. Vol. 13. Pp. 3457-3467.

15. Shokin, Yu.I., Beyzel, S.A., Richkov, A.D., Chubarov, L.B. Chislennoe modelirovanie nakata voln na poberezhie s ispolzovaniem metoda krupnih chastits [Numerical simulation of the tsunami runup on the cost using the method of large particles]. Mathematical modelling. 2015. Vol. 27. No. 1. Pp. 99-112. (rus).

16. Adesoji, O.A. Development of a discrete simulation model for tsunami. J. Math. Comput. Sci. 2018. Vol. 8. No. 1. Pp. 98-113.

17. Belyaev, N.D., Lebedev, V.V., Nudner, I.S., Mishina, A.V., Semenov, K.K., Shchemelinin, D.I. Eksperimentalnie issledovaniya vozdeystviya voln tipa tsunami na grunt u osnovaniy morskih gravitatsionnih platform [Experimental study of tsunami-type waves impact on soil at foundations of offshore gravity platforms]. Magazine of Civil Engineering. 2014. No. 6. Pp. 4-12. (rus).

18. Mazinani, I., Ismail, Z., Ahmad Mustafa, A.H., Saba, A. Experimental Investigation on Tsunami Acting on Bridges. Int. Journal of Civil and Environmental Engineerig. 2014. Vol. 8. No. 10. Pp. 1092-1095.

19. Harahap, I., Vo, H. Generation, propagation, run-up and impact of landslide triggered tsunami: a literature review. Applied Mechanics and Materials. 2014. Vol. 567. Pp. 724-729.

20. Nudner, I.S., Lebedev, V.V., Semenov, K.K., Belyaev, N.D., Beizel, S.A., Chubarov, L.B. Issledovaniya rasprostraneniya volni tsunami v buhte [Studies of tsunami waves propagation in a bay]. Trudi XIII Vserossiyskoy conferencii "Prikladnie tehnologii gidroakustiki I gidrofiziki" [Proc. XIII AllRussian conference "Applied technologies of hydroacoustics and hydrophysics"]. Saint-Petersburg. 2016. Pp. 159-161. (rus).

21. Allsop, W., Chandler, I., Zaccaria, M. Improvements in the physical modelling of tsunamis and their effects. Proceedings of the 5th International Conference on the application of Physical Modelling to Port and Coastal Protection, Coastlab14. Varna, Bulgaria. 2014.

22. Prokofyev, V.A. Primenenie graficheskogo processora dlya uskoreniya resheniya trehmernih prikladnih zadach gidravliki otkritih potokov [GPU utilization for speed up the solution of the three-dimensional engineering tasks of open flows hydraulics]. Mathematical modelling. 2017. Vol. 29. No. 8. Pp. 74-94. (rus).

23. Belash T.A., Yakovlev A.D. Seismic stability of a tsunami-resistant residential buildings. Magazine of Civil Engineering. 2018. 80(4). Pp. 95-103. doi: 10.18720/M CE. 80.9

24. Hsu, C.A. SEC-HY21: a numerical model for two-dimensional open channel flows. Proc. XXIX IAHR Congr. Beijing. 2001. Theme D. Pp. 821-827.

25. Bowden, K.F. Physical oceanography of coastal waters. Ellis Horwood, New York, 1983. 324 p.

26. Kulikovskii, A.G., Pogorelov, N.V., Semenov, A.Y. Matematicheskie voprosi chislennogo resheniya giperbolicheskih sistem [Mathematical problems on numerical solution of hyperboplic systems]. Moscow: Nauka, 2001. 608 p. (rus).

27. Chen, Z., Wang, G., Wang, Z. Numerical solution of the two-dimensional unsteady depth-averaged flow and solute transport. Proc. XXIX IAHR Congr. Beijing. 2001. Theme D. Pp. 725-733.

28. Sanders, B.F. High-resolution and non-oscillatory solution of the St. Venant equations in non-rectangular and non-prismatic channels. Journal of Hydraulic Research. 2001. Vol. 39. No. 3. Pp. 321-330.

Математическое моделирование. 2015. Т. 27. № 1. C. 99-112.

16. Adesoji O.A. Development of a discrete simulation model for tsunami // J. Math. Comput. Sci. 2018. Vol. 8. No. 1. Pp. 98-113.

17. Беляев Н.Д., Лебедев В.В., Нуднер И.С., Мишина А.М., Семенов К.К., Щемелинин Д.И. Экспериментальные исследования воздействия волн типа цунами на грунт у оснований морских гравитационных платформ // Инженерно-строительный журнал. 2014. № 6. С. 4-12.

18. Mazinani I., Ismail Z., Ahmad Mustafa A.H., Saba A. Experimental Investigation on Tsunami Acting on Bridges // Int. Journal of Civil and Environmental Engineerig. 2014. Vol. 8. No. 10. Pp. 1092-1095.

19. Harahap I., Vo H. Generation, propagation, run-up and impact of landslide triggered tsunami: a literature review // Applied Mechanics and Materials. 2014. Vol. 567. Pp. 724-729.

20. Нуднер И.С., Лебедев В.В., Семенов К.К., Беляев Н.Д., Бейзель С.А., Чубаров Л.Б. Исследования распространения волны цунами в бухте // Труды XIII Всероссийской конференции «Прикладные технологии гидроакустики и гидрофизики», Санкт-Петербург, 2016.

C. 159-161.

21. Allsop W., Chandler I., Zaccaria M. Improvements in the physical modelling of tsunamis and their effects // Proceedings of the 5th International Conference on the application of Physical Modelling to Port and Coastal Protection, Coastlab14, Varna, Bulgaria. 2014.

22. Прокофьев В.А. Применение графического процессора для ускорения решения трехмерных прикладных задач гидравлики открытых потоков // Математическое моделирование. 2017. Т. 29. № 8. C. 74-94.

23. Белаш Т.А., Яковлев А.Д. Сейсмостойкость цунамистойких жилых зданий // Инженерно-строительный журнал. 2018. № 4(80). С. 95-103. doi: 10.18720/MCE.80.9

24. Hsu C.A. SEC-HY21: a numerical model for two-dimensional open channel flows // Proc. XXIX IAHR Congr. Beijing. 2001. Theme D. Pp. 821-827.

25. Bowden K.F. Physical oceanography of coastal waters. Ellis Horwood, New York, 1983. 324 p.

26. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем. М: Наука. 2001. 608 с.

27. Chen Z., Wang G., Wang Z. Numerical solution of the two-dimensional unsteady depth-averaged flow and solute transport // Proc. XXIX IAHR Congr. Beijing. 2001. Theme

D. Pp. 725-733.

28. Sanders B.F. High-resolution and non-oscillatory solution of the St. Venant equations in non-rectangular and non-prismatic channels // Journal of Hydraulic Research. 2001. Vol. 39. № 3. Pp. 321-330.

29. Prokofiev V.A. Applying of different modifications of a finite volume method for simulation of opened flows // Proc. IAHR Symp. "Hydraulic and hydrological aspects of reliability and safety assessment of hydraulic structures". St.-Petersburg. 2002. paper C10.

30. Fletcher C. Computational Techniques for Fluid Dynamics. Vol.1. Berlin, Heidelberg: Springer-Verlag. 1988.

31. Климович В.И. О граничных условиях в рамках модели мелкой воды для задач, связанных с численным моделированием волновых процессов // Доклады 8-й научно-технической конференции «Гидроэнергетика. Новые технологии и разработки». Санкт-Петербург, 2014. C. 124-131.

32. Klimovich V.I. Verification of 2D numerical algorithms for Dam-Break flows simulation // Proc. IAHR Symp. "Hydraulic and hydrological aspects of reliability and safety

29. Prokofiev, V.A. Applying of different modifications of a finite volume method for simulation of opened flows. Proc. IAHR Symp. "Hydraulic and hydrological aspects of reliability and safety assessment of hydraulic structures". St.-Petersburg. 2002. paper C10.

30. Fletcher, C. Computational Techniques for Fluid Dynamics. Vol.1. Berlin, Heidelberg: Springer-Verlag. 1988.

31. Klimovich, V.I. O granichnih usloviyah v ramkah modeli melkoy vodi dlia zadach, sviazannih s chislennim modelirovaniem volnovih protsessov [About boundary conditions within shallow water model for numerical modeling of wave processes]. Docladi 8 nauchno-tehnicheskoy konferencii "Gidroenergetika. Novie tehnologii i razrabotki" [Proceedings of the 8th scientific-technical conference "Hydroenergetics. New technologies and developments"]. Saint-Petersburg. 2014. Pp. 124-131. (rus).

32. Klimovich, V.I. Verification of 2D numerical algorithms for Dam-Break flows simulation. Proc. IAHR Symp. "Hydraulic and hydrological aspects of reliability and safety assessment of hydraulic structures", St.-Petersburg. 2002. paper C30.

33. Tablitsi prilivov na 1996 god (Tom 3: Zarubezhnie vodi. Severniy ledovitiy, Atlanticheskiy i Indiyskiy okeani). [Tide tables at 1996 year (v.3 Foreign waters. Arctic, Atlantic and Indian Oceans)]. St.-Petersburg: Central Department of Navigation and Oceanography, 1995. 586 p. (rus).

34. SP 292.1325800.2017. Zdanya I sooruzhenya v tsunami opasnih rayonah. Pravila proectirovanya. [Russian Set of Rulles SP 292.1325800.2017. Buildings and structures on tsunami hazardous areas. Regulations of design] Moscow, 2017. 138 p. (rus).

assessment of hydraulic structures". St.-Petersburg. 2002. paper C30.

33. Таблицы приливов на 1996 год (т.3: Зарубежные воды. Северный ледовитый, Атлантический и Индийский океаны). СПб: Главное управление навигации и океанографии, 1995. 586 с.

34. СП 292.1325800.2017. Здания и сооружения в цунами опасных районах. Правила проектирования. М., 2017. 138 с.

Vitaly Klimovich*,

+7(812)493-93-53; klimvita@yandex.ru Oleg Voronkov,

+7(812)493-93-90; voronkovok@vniig.ru

Vyacheslav Davidenko, +7(812)493-93-33; davidenkovm@vniig.ru

Sergey Ivanov,

+7(961)804-50-78; S.Ivanov1948@gmail.ru

Виталий Иванович Климович*,

+7(812)493-93-53;

эл. почта: klimvita@yandex.ru

Олег Константинович Воронков,

+7(812)493-93-90;

эл. почта: voronkovok@vniig.ru

Вячеслав Михайлович Давиденко,

+7(812)493-93-33;

эл. почта: davidenkovm@vniig.ru

Сергей Николаевич Иванов,

+7(961)804-50-78;

эл. почта: S.Ivanov1948@gmail.ru

© Klimovich, V.I.,Voronkov, O.K.,Davidenko, V.M.,Ivanov, S.N., 2018