Quasi-two-layer finite-volume scheme for modeling shallow water flows with the presence of external forces

1 Quasi-two-layer finite-volume scheme for modeling shallow water flo ws over an arbitrary bed in the presence of external force K. V. Karelsky 1 , A. S. Petrosyan 1, 2 , A. G. Slavin 3 1 Space Research Institute of the Russian Academy of Sciences, 117997, Moscow, Profsoyuznaya st. 84/32, 117997, Moscow 2 Moscow Institute of Physics and Technology (S tate University), 9, Institutskii per., Dolgoprudny 41700, Moscow Region, Russia 3 Memorial University of Newfoundl and, Newfoundland and Labrador, Canada Email addresses: kkarelsk@iki.rssi.ru (K. V. Karelsky), apetrosy@iki.rssi.ru (A. S. Petrosyan), aslavin@mun.ca (A. G. Slavin). Finite-volume numerical method for study shallow water flows over an arbitrary bed profile in the presence of ex ternal force is proposed. This method uses the quas i-two- layer model of hydrodynamic flows ov er a stepwise boundary with advanced consideration of the flow feat ures near the step. A distinctive feature of the suggested model is a separation of a studied flow into two layers in calculating flow quantities near each step, and improving by th is means approximation of depth-averaged solutions of the initial three-dimensional Euler eq uations. We are solving the shallow-water equations for one layer, introducing the fictitious lower laye r only as an auxiliary st ructure in setting up the appropriate Riemann pr oblems for the upper layer. Besides quasi-two-layer approach leads to appearance of additional terms in one-layer fini te-difference representation of balance equations. These terms provide th e mechanical work made by nonhomogeneous bed interacting with flow. A notable ad vantage of the proposed method is the consideration of the properties of the process of the waterfall, namel y the fluid flow on the step in which the fluid does not wet part of the vertical wall of the step. The presence of dry zones in the vertical part of the step indicates violation of the conditions of hydrostatic flow. The quasi-two-layer approach determines t he size of the dry z one of the vertical component of the step. Consequent ly it gives an opportunity to figure out the amount of flow kinetic energy dissipat ion on complex boundary. Numerical simulations are performed based on the proposed al gorithm of various physical phenomena, such as a breakdown of the rectangular fluid column over an inclined plane, large-scale motion of fluid in the gravity field in the presence of Coriolis force over an mounted obstacle on underlying surface. Com putations are made for two dimensional dam-break problem on slope precisely conform to laboratory experiment s. Interaction of the Tsunami wave with the shore line including an obstacle has been si mulated to demonstrate the effectiveness of the developed algorithm in domains including partly flooded and dry regions. PACS codes: 47.10.-g; 47.10.A-; 47.11.-j; 47.35.-i Keywords : numerical methods, shallow water equat ions, source term, Riemann problem, arbitrary bed. ---------------------------------------- Corresponding Author: Professor Arakel Petrosyan Space Research Institute 2 Russian Academy Sciences 117997, Moscow, Profsoyuznaya st. 84/32 Tel: +7-495-3335478 Fax: +7-495-3333011 E.mail: apetrosy@iki.rssi.ru 1. Introduction The shallow water equations are derived fr om the non-stationary three-dimensional Euler equations by averaging over a vertical coordinate and taking into consideration hydrostatic pressure distribution [49]. It is expected that soluti ons of depth-averaged equations have similar properties as depth-aver aged solutions of initial fluid equations. The obtained equations are also rather complicated for constructing general analytical solutions because of their nonlinearity and bed co mplexity [33, 34, 35]. However, they can be successfully integrated numerically [10, 50, 53]. The main difficulty in numerical simulati on of nonhomogenous shallow water equations consists in their nonlinearity and their non-dive rgence property determined by nonhomogeneity of the right-hand si de of the momentum conser vation equations due to bed complexity. The presence of a non-divergent term induces highly nonlinear effects caused by stepwise change of hydrodynamic quantit ies in the areas of its sharp change in addition to nonlinear phenomena due to hyperbolic structure of shallow watyer equations. The other problem consists in compatibility of solutions of traditional de pth-averaged equations with depth-av eraged solutions of initial Euler equat ions due to significant role of shallow flows dependences on vertical coor dinate [28]. Numerical methods have been developed and effectively used in studies of shallow water flows on c omplex boundary when an external force effe cts are insignificant. The main problem in numeric al integration of hyperbolic balance equations is in developing finite-difference schemes satisfyi ng the conditions of conservation of stable states. In the context of Sain t-Venant’s equations this means the equilibrium of resting water. The schemes which satisfy such proper ties are called well-balanced. Roe modified scheme was suggested to use for satisfying the well-balanced conditions in [7]. This idea has been further developed in papers [43, 51]. The turbulence phenomena and mass sources in problems comprising nonhomogeneous bed are considered in papers [14, 15], and the improved reconstruction in [12, 13]. It was suggested in [40, 41] to use either ordinary or modified Godunov-type schemes (Riemann-solvers). This approach is well- proved and this idea has been widely applied to solve the problems with various source terms. In [5, 37] the method in which bed effects are compensated by changing of the values of fluid depths was su ggested. In [3, 9, 10, 44] wa s proposed the method that uses the scheme of hydrostatic reconstruction of the flows on cells interfaces. The phenomena appearing because of presence of the Coriolis force were simulated using this scheme, particularly the geostrophic adaptation phenomenon. The sa me phenomena were considered in paper [38] in which LeVeque method is us ed [40, 41]. The propagation of tsunami waves and flooding processes are simula ted in [19, 20, 21]. In [54] the surface gradient method of interpreting source terms in s hallow water equati ons is proposed. It is based on the accurate reconstruction of conser vative variables on cell interfaces. The method of separating flow quantities which keeps exac t balance between the flow gradients and source terms is proposed in paper [25]. The method of high-order accuracy and its im proving reconstruction is proposed in [52]. The exact partial solutions of the pr oblem including a step- wise boundary on a bed are presented in [1]. The exact partial solutions for the problem including nonhomogeneous bed are presented in [8]. Also we have to mention the works [2, 4, 6, 3 11, 18, 21, 22, 23, 24, 50, 55] which made a significant contribution to solution of problems of simulating hydrodynamic fl ows over a nonhomogeneous bed with source terms of various origin. The main difficulty aris ing in applying such methods is the lack of a unique analytical solution of the one-dimensional nonlinear Riemann problem for step bed [27]. However, in all afor ementioned works this problem was successfully solved by introducing some additional assumptions. The presence of external forces, bed nonhomogeneities often causes vertical nonhomogeneity of horizontal shal low flows [28] that changes the values of hydrodynamic quantities averaged over the depth. It is necessary to take into consideration the vertical structure of depth-averaged fluid flows near the peculiarities of the bed to describe the indicated effects more adequately in numerical solutions . Existence of well-developed and validated numerical methods makes especially attractable the reducti on of a problem in the presence an external force to solution of the problem of shallow water flows over a composed time-dependent bed [9, 10, 38, 40, 41]. In present paper we propose to use the Riemann-solver which is ad apted to the flow parameters for calculating shallow water flows over an arbitrary bed in the presence of external force. The proposed method belongs to the family of methods based on t he solution of the dam-break problem. The method consists in reducing of the problem to successive solutions of classical shallow water equations on the flat plane using God unov method with allowance for the vertical nonhomogeneity effect in calculating the fluxes through the boundaries of cells adjoining to stepwise boundaries. The vertical non homogeneity leads to the Riemann problem solution on a step based on the quasi-two-laye r shallow water model developed in [29, 30, 31, 32, 47, 48]. We are solving the shallow- water equations for one layer, introducing the fictitious lower layer only as an auxiliary structure in setting up the appropriate Riemann problems for the upper layer. Be sides quasi-two-layer approach leads to appearance of additional terms in one-layer finite-difference representation of balance equations. These terms provide the mech anical work made by nonhomogeneous bed interacting with flow. Quasi-tw o-layer model 29, 30, 31, 32, 47, 48 is extended here for a generalized time-dependent bed which represents an external force. The term time- dependent bed means that at each time mo ment bed can have different values. The main difficulty in model ing of fluid flows over a co mplex bed consists in that both partly and complete flooded domains may take place. Partly flooded domains may appear in flows when fluid depth is rela ted to the value of bed gradient when approximated by steps. In this case the step wall is partl y wet ting and the algorithm suggested in present work cons iders exactly mechanical work done only by this part of the step. Mechanical work done by nonhomogeneous bed means that the work is done in the process of interaction of non-stationar y flow with nonhomogeneous bed, with the step and with part of the step. Partly flooded domains present a r eal challenge for most finite- difference schemes and require special efforts need to be made to capture such domains. Algorithm suggested in our work av oids this difficulty in a natur al way. In present work we extend algorithms developed in [ 29, 30, 31, 32, 47, 48] to the case of multiply-connected domains including partly flooded and dry regions. Section 2 presents the shallow water model over an arbitrary bed profile in the presence of external force. The correspondi ng shallow water equatio ns are written and the idea of the suggested Riemann-solver adapt able to flow parame ters is described. Section 3 presents the quasi-two-layer model and the developed finite-difference scheme. Section 4 reviews the well known finite-diffe rence schemes for calculating shallow water flows accounting for the external force effect in the form of fictitious bed, and compares them with our scheme. We have demonstrated results of the co mputations that show the rate of formation of an adequate flux of the mass and the impulse. We compare our results with those obtained by the well-balan cing method. Section 5 presents the results of test computations with a constant exte rnal force corresponding to an inclined bed. The 4 simulation of large-scale motion of heavy flui d (gas) is performed with allowance for the Coriolis force over a complex trigonomet ric function shape bed. Calculations and comparisons with laboratory experiment of two-dimensional dam-break problem on a slope are performed. Interaction of the Ts unami wave with the shore line including an obstacle has been simulated. The main results of the work are formul ated in section 6. 2. Model of shallow water fl ows over an arbitrary bed in the presenc e of external force We consider two-dimensional shallow wa t er equations for fluid flows over an arbitrary nonhomogeneous bed with the external force in the divergence form [49, 53]: () ( ) () () () () 22 22 0 2 2 x y hu hv h tx y hu gh hu huv db gh E h tx y d x hv gh hv hvu db gh E h ty x d y ∂∂ ⎧ ∂ ++= ⎪ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − + ⎨ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − + ⎪ ∂∂ ∂ ⎩ , (1) where ) , , ( t y x h is the depth of fluid, (, , ) ux y t is the depth-averaged horizontal velocity in the x direction, (, , ) vx y t is the depth-averaged horizontal veloci ty in the y direction, (, ) bx y is the function describing the bed shape, ( ) (, ,, , , ) , x y E hu v x yt E E = is the external force, g is the gravity acceleration. The first equation in system (1) r epresents the mass conservation law; the second and third equations represent the momentum conservati on laws for corresponding velocity vector components. Let’s introduce fluxes 22 22 () / 2 , () /2 hu hv F hu gh G huv hvu hv gh ⎛⎞ ⎛ ⎞ ⎜⎟ ⎜ ⎟ Ψ= + Ψ= ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ + ⎝⎠ ⎝ ⎠ , where h hu hv ⎛⎞ ⎜⎟ Ψ= ⎜⎟ ⎜⎟ ⎝⎠ , and the source term 0 x y db Sg h E h dx db gh E h dy ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ =− + ⎜⎟ ⎜⎟ ⎜⎟ −+ ⎜⎟ ⎝⎠ . 5 Thus we rewrite equations (1) as: () () tx y FG S ∂Ψ + ∂ Ψ + ∂ Ψ = . The corresponding one-dimens ional problem is obtain ed by reducing the system (1) over one of spatial coordinates with an accuracy of projection of the momentum equation on a non-reduced spatial direction. This corresponds to the one-dimensional spatial problem of SWE over a time-dependent bed. Neglecting in (1), for certainty, the derivatives with respect to y we get the following equations: () () () 22 0 2 x hu h tx hu gh hu db gh E h tx d x ∂ ⎧ ∂ += ⎪ ∂∂ ⎪ ⎨ ∂+ ∂ ⎪ += − + ⎪ ∂∂ ⎩ (2) Here (, ) uu x t = , () bb x = and ( ,,, ) xx EE h u x t = . The problem becomes essentially two-dimensional if t he projection of external force E onto one of the coordinate axes depends either on the full velo city vector, or on the projection of a velocity vector onto other coordinate axis. Neglecting one of the spatial variables during the solving an essentially two-dimensional pr oblem results in viol ating the momentum conservation law. This, in turn, makes nece ssary to introduce some fictitious work to compensate these violations, in spite of all non-physical character of such compensation. Such situation exists, for example in modeli ng of rotating shallow water, when Coriolis force acts as an external force. Use of quasi-two-layer model allows improving representation of transversal co mponent of a velocity v ector under the assumption of its stationarity on each time step that in particular leads t o impr ovement of conservation of the potential vorticity vector. We will use th is advantage of developed scheme in future developments. If E depends only on h, u , x and t , splitting of system (1) into a set of subsystems (2) is physically justified. Existence of well-develope d and validated numerical methods makes especially attractable the reduction of a problem in the pr esence an external forc e to solution of the problem of shallow water flows over a composed time-dependent bed. Actually, to obtain a formal coincidence of the types of equations, external force must be presented in the following form: x dk gE dx = − , (3) Where (, ) kk x t = is the fictitious bed. As it was mentioned earlier x E depends on h, u , x and t . Function k depends only on x and t because h and u also depend only on x and t . Thus system (2) can be rewritten as: () () () 22 0 2 () hu h tx hu gh hu db k gh tx d x ∂ ⎧ ∂ += ⎪ ∂∂ ⎪ ⎨ ∂+ ∂ + ⎪ += − ⎪ ∂∂ ⎩ . (4) 6 Such a presentation of external force in sh allow water equations essentially narrows the class of possible solutions, excluding from c onsideration spatially non- integrable functions specifying the external force E . This condition is generally needed to write down formally system (4), and may not be signific ant in finite difference repr esentations of all functions defined in discrete points. If we introduce the difference grid with a spatial mesh size Х in the х direction and consider two adjacent cells in х , then, according to formula (3) we have in the left cell after averaging the fictitious bed with a slope tangent l E g , and with the tangent r E g in the right cell. Therefore, the corresponding di fference of average heights of a fictitious bed on the boundary of the tw o adjacent cells is: 2 rl EX g E X g + (5) The possibility of simulating the shallow wate r in the presence of external for ce by introducing the fictitious time-dependent bed has been used for a long time for constructing numerical models [9, 10, 38, 40, 41]. The fact that the methods developed for composed bed can be useful in solving the problems with an external force, apparently, causes no doubts. In sp ite of the fact that the in fluence of fictitious bed, obviously, directly leads to accounting for bed work performed over the fluid flow, the external force, by itself, could not commit any work. The method proposed in this paper makes it possible to provide the physical inte rpretation of the ment ioned formal adaptation of the methods. This makes it possible to better understand the limit s of applicabilit y of the given approach for splitting finite-differenc e schemes, being not limited to a particular applied model. The method proposed in this work allows to visualize the features of the splitting approach to numerical simulation as a whole and, t hus, provides physical ideas for improving stability criteria in t he finite-difference implementation. The underlying problem for all Godunov-type methods is the Riemann problem that is the one-dimensional formulati on of the Cauchy problem for two semi-infinite domains, in each of which the values of all hydrodynam ic parameters are cons tant at the initia l time. It is clear that in t he presence of the external y -dependent force such a problem is physically meaningful for a rather small time intervals only, si nce its general solution is sought as a set of partial self-similar solutions. As a result, the problem of calculating of the shallow water flows in the presence of external force becomes identical to the probl em of calculating the shallow water flows over a time-dependent nonhomogeneo us bed. Papers [29, 30, 31, 32, 47, 48] have offered the quasi-two-layer method for ca lculating shallow water flows over an nonhomogeneous bed. In present work we show how quasi-two-layer method can be applied in case of a time- dependent bed. The effective bed pr ofile is appr oximated by a piecewise constant function separating it into a finite number of domains with a stepwise boundary. 3. Quasi-two-layer model and finite-differe nce scheme for the equations of shallow water over an arbitrary bed profile in the presence of external force In this secti on we briefly describe basics of quasi- two-layer model [29, 31, 47] are required for understanding of suggested comput ational algorithm andr derive finite- difference scheme for the equati ons of shallow water over an arbitrary bed profile in the presence of external force. This model is extende d here for non-stationary boundaries and for flows that include partly flooded and dr y regions on complex bed. Distinctive feature of the suggest ed model is a separation of a stud ied flow into two layers in 7 calculating flow quantities near each step, and improving by th is means approximation of depth-averaged solutions of the initial three-dimensional Eu ler equations. The uniquenes s of such a separation into two layers is provided by the uniquene ss of solution of the Dirichlet problem for definin g this boundary. One of the advant ages of the developed method is that it allows one to take into a ccount both flow velocity and the height of the step at each point of space and at any time instant. We are solv ing the shallow-water equations for one layer, introduci ng the fictitious lower layer only as an auxiliary structure in setting up the appropriate Riemann problem s for the upper layer. Besides quasi-two- layer approach leads to appearance of additional terms in one-layer finite-difference representation of balance equations. These terms provide the mechanical work made by nonhomogeneous bed interacting with flow. Suggested approach modifies properly values of hydrodynamic fluxes in one-la yer finite-difference scheme for nonhomogeneous bed as well by adjusting flow to that for upper layer in really two-layer model in the vicinity of discrete steps. Below we develop finite-d ifference representation of the external force in numerical Godunov-type models for flows of shallow water based on quasi-two- layer model [29, 31, 47]. 3.1 Quasi-two-layer model The basic idea of the quasi-two-layer m odel and model application for development of finite-difference scheme for shallow water equations over fully flo oded bed is described in [30, 32, 48]. In this sect ion we extend these id eas to the case of existence partially inundated ar eas in bed shape and to time dependent bed to take into account the external force influence on the flow. We consider a shallo w non-viscous flow ov er a step of height 0 b turned to the left with no loss of generality (Fig. 1). On fig.1 dashed line represents interface between layers. If the fluid height is lower than the st ep height, the step would act as an impermeable boundary for the lower layer. Thus, if the height of a fluid is higher than that of the step, we may consider two fluid laye rs: the lower layer wit h the flow parameters 1 h and 1 u for which the step is an impermeabl e boundary and the upper one with the flow parameters 2 h and 2 u for which there is no direct influence of the step, obviously, 12 l hh h =+ . Equations for the two-la yer system are following: 1 11 2 22 2 22 11 11 1 1 1 22 22 22 2 2 () 0 , () 0 , 1 () ( ) 0 , 2 1 () ( ) 0 , 2 h hu tx h hu tx h hu hu g h g h tx x h hu hu g h g h tx x ∂∂ ⎧ += ⎪ ∂∂ ⎪ ∂∂ ⎪ += ⎪ ∂∂ ⎨ ∂∂ ∂ ⎪ ++ + = ⎪ ∂∂ ∂ ⎪ ∂∂ ∂ ⎪ ++ + = ∂∂ ∂ ⎩ (6) with the corresponding boundary condition 1 0 u = for 0 x = , 0 t ≥ (the step affects the height and velocity of the lower layer near the step) and with the initial conditions for 0 t = : * 1 hh = , * 2 l hh h =− , 1 l uu = , 2 l uu = for 0 x ≤ , 8 2 r hh = , 2 r uu = for 0 x > . (7) To solve system (6) we split it into two so that the variables 1 h and 1 u can be calculated independently of the variables 2 h and 2 u . Assuming that the wave pattern in the lower layer is formed much faster than t hat in the upper one, the layers are separated so that the lower layer interacting with the step is at rest and forms a common horizontal plane with the step. This makes it possible to disregard the term 12 gh h x ∂∂ in the third equation of system (6). We s uppose that the upper layer does not appreciably affect the lower one, i.e., we disregard the term 21 gh h x ∂ ∂ . Then the effect of the lower layer on the upper one is taken into account by changing the initial condition * 2 ,0 l hh h x =− ≤ , which is determined by the change in the lower laye r depth due to the st ep drag. As confirmed below by numerical calculations, these assu mptions allow us to describe adequately the resulting flow. The adequacy in this case should imply t he compatibility between the numerical results and theoretical ly possible ones when analytically solving the problem of arbitrary discontinuity decay over the step for shallow water equations. Thus, under the above assumptions, we have the quasi-two-layer model described by the equations where the variables 1 h , 1 u and 2 h , 2 u are now connected only by the initial parameter * h which should be chosen so that the lower layer he ight directly at the step coincides with the step height 10 0 x hb = = for 0 ≥ t . Consequently, the problem of finding * h reduces to the solution of the Dirichlet inverse problem: 1 11 11 11 1 22 () 0 , 1 () ( ) 0 , 2 h hu tx hu hu g h tx ∂∂ ⎧ += ⎪ ⎪ ∂∂ ⎨ ∂∂ ⎪ + += ⎪ ∂∂ ⎩ (8) * 1 h h = for 0 < x , 0 = t ; 1 l uu = for 0 ≤ x , 0 = t ; 0 1 = u for 0 = x , 0 ≥ t ; 1 0 hb = for 0 = x , 0 ≥ t . The solution of system (8) depends on the flui d flow direction, that is, on velocity l u : а ) 0 l u > , the flow has a form of the shock wa ve reflected to the left, on the left side of which the fl ow parameters are * h and 0 l u > . On the right side the fluid is at rest, that is, h = 0 b and u = 0, and the closing depth, determined by the step effect, is calculated as follows : * * 0 0 * 0 () () 2 l gb h ub h bh + =− ; (9) The solution is determined directly from the algebraic formula (9) that follows from Hugoniot’s relations on a hy draulic jump [37]; b) 0 l u < , the flow has a form of rarefaction wave moving to the left, on the left side of which the flow parameters are * h and 0 < l u , and on the right side as h = 0 b and u = 0 , Therefore, the solution is determined from the corresponding Riemann’s invariant, 9 which describes the rarefaction wave: 2 * 0 11 2 l hg b u g ⎛⎞ =− ⎜⎟ ⎝⎠ . (10) Expression (10) determines, in the explicit form, the depth of the lower part of a flow that is completely stopped by the step. Thus, we have found * h and for 0 x ≤ we now have: * 1 hh = and * 2 l hh h =− . The assumption that the wave pattern in the lower layer is formed much faster than that in the upper one is not valid for the case of a ra ther small step compared to the other flow parameters. In this situation it is necessary to assume that the wave pattern in the upper layer is formed faster than that in the lower one. The relati on between the layers' heights has the form: * 2 hh = and * 1 l hh h =− . For the boundary situation between the above possibilities, correct splitting of the system becomes impossible, the two-layer model no longer adequately describes the physical process and it is necessary to use the complete Euler equations as a minimum for the fluid layer in the ε -neighbourhood of the step boundary edge. The situation is e ssentially different in the case when * h is outside interval ( ] 0, l h . Physically, it is possible in two situations. First, * h may be equal to zero. This implies absence of the lower layer, and consequently absence of the fluid and/or changes of the height of the bed. Second possibility is * l hh > . It corresponds to the full deceleration of the layer at the left side, i.e. in such si tuation the bed is a non-leaking boundary for all fluid at the left and ther efore, it is necessary to accept that 2 0 h = for 0 x < . Then value 1 h for 0 x = can be found as a solution of the problem of interaction of all fluid at the left with a non-leaking boundary. After evaluation of * h one can find the values of flow quantities as the solution of classical Riemann problem for sh allow water over a smooth bed with corresponding initial conditions. Using the quasi-two-layer met hod, one can find the flow quantities on cell faces, which determine, in the difference scheme, the values of all hydrodynamic parameters in the next time layer. Thus, depending on flow param eters near the stepwise bed, the depths of fluid’s lower and upper layers are calculated, and for calculating the flow parameters on a face the classical Riemann problem for upper laye r’s parameters is solved. The proposed method is adapted to the flow parameters and allows one to ta ke into account the fluid flow features at each point of space and at each time instant. It should be noted that these features are considered in suggested m odel only in two-layer approximation in the steps vicinity. 3.2. Derivation of finite difference scheme To obtain the difference scheme we integrat e the system of equations of two-layer shallow water with the external fo rce over a nonhomogeneous bed (1): 10 ( ) ( ) ( ) ( ) ( ) () () () () () () () () 12 1 1 2 2 1 1 2 2 22 2 2 11 1 2 2 2 11 2 2 11 1 22 2 2 1 12 12 22 11 1 11 0 22 0 2 G G x x hh h u h u h v h v dxdy dt tx x y y h u gh h u gh hu h u hu v tx y t x dxdydt db hu v h h gh h E h h g yd x x hv g h hv t ∂+ ∂ ∂ ∂ ∂ ⎛⎞ ++ + + = ⎜⎟ ∂∂ ∂ ∂ ∂ ⎝⎠ ⎛⎞ ∂+ ∂ + ∂∂ ∂ ⎜⎟ ++ + + + ∂∂ ∂ ∂ ∂ ⎜⎟ = ⎜⎟ ∂∂ ⎜⎟ ++ + − + + ⎜⎟ ∂∂ ⎝⎠ ∂+ ∂ + ∂∂ ∫∫∫ ∫∫∫ () () () () 22 22 2 22 11 1 22 2 2 1 12 12 2 0 G y y hv g h hv hv u yx t y dxdy dt db hv u hh gh h E h h g xd y y ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎛⎞ ∂+ ∂ ∂ ⎪ ⎜⎟ ++ + + ⎪ ∂∂ ∂ ⎜⎟ = ⎪ ⎜⎟ ∂∂ ⎪ ⎜⎟ ++ + − + + ⎜⎟ ⎪ ∂∂ ⎝⎠ ⎩ ∫∫∫ , (11) where (, , ) Gxy t is the arbitrary nonempt y domain, which is hom eomorphic to a spatial cube, 11 1 ,, hu v are lower layer parameters, 22 2 ,, hu v are upper layer parameters. Passing in (11) to surface integrals, we obtain the following integral form: () () () () ( ) ( ) () 12 1 1 1 1 2 2 2 2 22 2 2 11 2 2 11 1 2 2 2 1 2 1 1 12 2 2 12 12 11 2 2 0 22 0 SS S S S SS x x SG S h h dxdy h u dydt h v dxdt h u dydt h v dxdt h u h u dxdy h u gh h u gh gh h dydt db h u v h u v dxdt E h h g h h dxdydt dx h v h v dxdy ++ + + + = ++ + + + + + ⎛⎞ ++ − + − + = ⎜⎟ ⎝⎠ + ∫∫ ∫∫ ∫∫ ∫∫ ∫∫ ∫∫ ∫∫ ∫∫ ∫∫∫ ∫∫ ww w w w ww w w () ( ) () () 22 2 2 11 1 2 2 2 1 2 1 1 12 2 2 12 12 22 0 S y y SG h v gh h v gh gh h dxdt db h v u h v u dydt E h h g h h dxdydt dy ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ++ + + + + ⎪ ⎪ ⎪ ⎛⎞ ++ − + − + = ⎪ ⎜⎟ ⎝⎠ ⎩ ∫∫ ∫∫ ∫∫∫ w w , (12) where S is the boundary of domain (, , ) Gxy t . We apply the integral conservation laws (12) to each cell, choosing S to be the surface determined by cell boundary at the time step τ . Taking into account the evolution of the lower layer one gets the following equations on the cell faces: 21 1 2 1 2 1 2 2 1 12 12 1 2 2 1 12 12 1 2 21 1 2 1 2 1 2 12 12 , , 0 , ; 12 , , 0 , ; 12 , , 0 , ; 12 , , 0 , ; . xx xx yy yy x xh H h k b u u U x xh H h k b u u U y yh H h k b v v V y yh H h k b v v V and h h H inside the cell −− ++ −− ++ =− ⇒ = = + = = =+ ⇒ = = + = = =− ⇒ = = + = = =+ ⇒ = = + = = +=      11 Here 12 , 12 xy ±± denote the coordinates for the bou ndary of the cell with a number ( x , y ). Here 12 1 x xx kk k −− =− , 12 1 x xx bb b −− =− , 12 1 x xx kk k ++ = − , 12 1 x xx bb b ++ = − , 12 1 yy y kk k −− =− , 12 1 yy y bb b −− =− , 12 1 yy y kk k ++ = − , 12 1 yy y bb b ++ = − . Besides, we apply the identity () 2 22 12 1 2 1 2 22 2 gh gh gh h g h h ++ ≡ + . Taking into account that fact that the bed and external force effects ar e now described by the hei ght of the lower layer 1 h (functions b(x,y) and k(x,y,t) correspondingly) one gets 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 0 xy xy y y tt tt tt xx xx xy xy y t y t xx tt yy yy xt xt y tt xy H dxdy H dxdy HU dydt H U dydt HV dxdt H U dxdt HU dxdy ττ τ ττ τ ++ ++ + + ++ =+ = = + = − −− −− − − ++ ++ =+ =− −− + =+ −− −+ − + +− = ∫∫ ∫∫ ∫ ∫ ∫ ∫ ∫∫ ∫∫ ∫         () () 12 12 12 12 2 2 2 2 12 12 12 12 12 12 12 12 12 2 2 22 1 2 1 2 2 12 12 12 12 12 1 2 1 2 xx y y t xx tt xy y t xx y x t tt xx yy yy yt xt t xx HU dxdy h U g h k b dydt h U g h k b dydt h UV dxdt H UV dxdt τ τ ττ ++ + + + ++ = −− − =+ + + + ++ −− =+ =− − − =− ⎛⎞ −+ + + + − ⎜⎟ ⎝⎠ ⎛⎞ −+ + + + − ⎜⎟ ⎝⎠ ∫∫ ∫ ∫ ∫ ∫∫ ∫∫        () () 12 12 12 12 12 12 12 2 2 2 2 12 12 12 12 12 12 12 12 12 2 2 22 1 2 1 2 2 12 12 12 0 1 2 1 2 x x xy xy x t yy tt tt xy xy x t yy x t yy xx xt yy HV dxdy HV dxdy h V g h k b dxdt h V g h k b dxdt h VU dyd τ τ τ + − ++ ++ + + ++ =+ = −− −− − =+ + + −− =+ − =− = ⎛⎞ −+ + + + − ⎜⎟ ⎝⎠ ⎛⎞ −+ + + + ⎜⎟ ⎝⎠ ∫∫ ∫∫ ∫∫ ∫ ∫ ∫∫      12 12 12 12 12 0 yy tt xx yt yt tH V U d y d t ττ ++ ++ =− −− ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ − = ⎪ ⎩ ∫∫ ∫∫   (13) Applying in (13) the mean value theor em and supposing that at cell boundary the values of all hydrodynamic par ameters represent the soluti ons of a corresponding one- dimensional problem dur ing the whole integration time step of, we obtain the difference scheme: 12 , 12 , 12 , 12 , , 12 , 12 , 12 , 12 ,, tt tt t t t t tt x y x y x y x y xy xy xy xy xy xy HU HU HV HV HH XY τ τ + −− ++ − − + + −− − =+ , () ( ) ( ) () () ( ) () () () () 2 12 , , 1 , , 1 , 2 12 , 1 , , 1 , , ,, , , 22 12 , 12 , 12 , 12 , ,1 2 ,1 2 , 2 2 2 / 2 tt t t t xyx x y x y x y x y tt t t t tt t t xyx x y x y x y x y xy xy xy xy tt t t xy xy x y xy tt xy xy xy gH i B B K K gH i B B K K HU H U HU HU HUV ττ τ −− − ++ ++ + −− + + −− ⎛⎞ +× − + + ⎜⎟ − ⎜⎟ ⎜⎟ ⎜⎟ +× − + + − ⎜⎟ = −+ Χ + ⎜⎟ ⎜⎟ +− ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ + () 12 , 12 , 12 , 12 tt t t xy xy xy HUV Y −+ + + − , 12 (14) () ( ) ( ) () () ( ) () () () () 2 ,1 2 , ,1 , ,1 2 ,1 2 ,1 , ,1 , ,, , , 22 ,1 2 ,1 2 ,1 2 ,1 2 12 , 12 , 1 2 2 2 / 2 tt t t t xy y xy xy xy xy tt t t t tt t t xy y xy xy xy xy xy xy xy xy tt tt xy xy xy xy tt xy xy x gH i B B K K gH i B B K K HV H V Y HU HU HUV ττ τ −− − ++ ++ + −− ++ −− − ⎛⎞ +× − + + ⎜⎟ − ⎜⎟ ⎜⎟ ⎜⎟ +× − + + − ⎜⎟ =− + + ⎜⎟ ⎜⎟ +− ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ + () 2, 1 2, 1 2 , 1 2, tt t t yx y x y x y HU V X ++ + − , Here τ is the time step; X and Y are spatial mesh sizes; H is the depth of fluid; U is the velocity in the x direction; V is the velocity in the y direction. , t x y H , , t x y U , , t x y V are the space averages of H  , U  , V  respectively. 12 , 12 t xy H ±± , 12 , 12 t xy U ±± , 1, 1 t x y V ± ± are the time averages of H  , U  , V  respectively. Subscripts x, y designate the values of function related to the center of masses of a cell with number (x, y) . Half-integers 12 , 12 xy ± ± designate the values of quantities at the boundary betwe en cells wit h numbers x , 1 x ± , and y , 1 y ± , respectively. Superscript t designates the number of a step in time, , x y K is the height of a fictitious bed and , x y B is the height of a bed. Respecti vely, according to formula (5): 1, , 1, , 22 xy x y xy x y E Xg E Xg KK + + + + = , ,1 , ,1 , 22 xy x y xy x y E Xg E Xg KK − − + + = , ,1 , ,1 , 22 xy xy xy xy E Xg E Xg KK + + + + = , ,, 1 ,, 1 22 xy xy xy xy E Xg E Xg KK − − + + = . The variables , x y ii assume either the value 0 – in t he case of negative difference of corresponding heights ( ) ( ) 1, , 1, , 2 xy x y xy x y BB K K ++ −+ + , ( ) ( ) ,1 , ,1 , 2 xy xy xy xy BB K K ++ −+ + , ( ) ( ) ,1 , , 1 , 2 xy x y x y x y BB K K −− −+ + , ( ) ( ) ,, 1 , , 1 2 xy xy xy xy BB K K −− −+ + of a effective bed (is equal sum of bed and fictitious bed), or the value of 0, 1 xy ii ≤ ≤ – in the case of positive difference. Variables x i , y i takes value 1 in case if * h is determined by equations (9,10) for a corresponding face of a cell and does not exceed value of depth inside the cell. Otherwise, variables x i , y i on a corresponding face are equal to the ratio of the depth formed under the total retardation of the flow on the indicated face to the corres ponding difference of heights of the effective bed. Thus unlikely algorithms developed in [29,30,31,32,47, 48], in finite-difference scheme (14) party flooded and dr y steps of effective bed are accounted. The values 12 , 12 t xy H ±± , 12 , 1 2 t xy U ±± , 12 , 12 t xy V ±± are calculated on the faces by solving the 13 corresponding Riemann problems considering that transversal fl ow velocity is transported convectively. We use convective transportation of the velocity hypot hesis in 1D problem instead of using 3 rd equation. It is clear that develope d finite-difference scheme is well- balanced type since it is balanced in case of ze ro velocities and free surface is horizontal. If we rewrite (14) in the terms of fluxes F( Ψ ) and G( Ψ ) we get: () () 1 12 , 12 , 12 , 12 , , 12 , 12 , 12 , 12 , tt x y x y xy x y xy x y xy xy xy xy t F F FBK FBK x t GG F B K F B K y + +− + − +− + − Δ Ψ= Ψ − − + − − Δ Δ −− + − Δ where () ( ) () () () ( ) () () 12 , 1 , , 1 , , 2 12 , 1, , 1, , 0 2 2 2 0 tt t t t xy x x y x y x y x y tt t t xy xx y x y x y x y gH i B B K K FBK gi B B K K ++ + + ++ ⎛⎞ ⎜⎟ ×× − + + + ⎜⎟ ⎜⎟ = ⎜⎟ ×− + + ⎜⎟ + ⎜⎟ ⎜⎟ ⎝⎠ , () ( ) () () () ( ) () () 12 , , 1 , , 1 , 2 12 , ,1 , , 1 , 0 2 2 2 0 tt t t t xy x x y x y x y x y tt t t xy xx y x y x y x y gH i B B K K FBK gi B B K K −− − − −− ⎛⎞ ⎜⎟ ×× − + + + ⎜⎟ ⎜⎟ = ⎜⎟ ×− + + ⎜⎟ + ⎜⎟ ⎜⎟ ⎝⎠ , () ( ) () () () ( ) () () ,1 2 ,1 2 ,1 , , 1 , 2 ,1 , ,1 , 0 0 2 2 2 xy tt t t t xy y xy xy xy xy tt t t yx y x y x y x y FBK gH i B B K K gi B B K K + ++ + ++ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ×× − + + + ⎜⎟ ⎜⎟ ⎜⎟ ×− + + ⎜⎟ + ⎜⎟ ⎝⎠ , 14 () ( ) () () () ( ) () () ,1 2 ,1 2 , ,1 , ,1 2 ,, 1 , , 1 0 0 2 2 2 xy tt t t t xy y xy xy xy xy tt t t yx y x y x y x y FBK gH i B B K K gi B B K K − −− − −− ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ×× − + + + ⎜⎟ ⎜⎟ ⎜⎟ ×− + + ⎜⎟ + ⎜⎟ ⎝⎠ . The studied equations are hyperbolic type . Therefore, by the analogy with the shallow water equations, in whic h the external force effect s is disregarded [37], the developed method uses the standard Courant – Friedrichs co ndition as the stability condition. This condition, however, should ta ke into consideration the velocity of propagation of perturbations, includi ng those in the fict itious lower layer. That is, the time step cannot exceed the minimal time during wh ich the perturbations in each layer pass the half of a cell. With regard to correcti on for two-dimensional statement we have: y x y x t t t t R Δ + Δ Δ Δ = τ , (15) where x t Δ , y t Δ are the minimum times of propagation of perturbations along the abscissa and ordinate axes, respectively; 1 < R is the additional factor for increasing the reliability. In calculations presented below this factor is taken equal to 0.4. The success of selecting the stability criterion in this form is confi rmed by test calculations presented below. For improving accuracy of developed fi nite-difference scheme we use 2 nd order accuracy computations in some tests (Subsecti ons 5.2,5.3,5.4). Increase of accuracy on spatial coordinate is reached by applying t he piecewise-linear reconstruction to the distribution of functions value in a cell with the use of the minmod limiter suggested for the first time by Kolgan [36,37] for the so lutions of accuracy problems of Godunov type methods in gas-dynamics: () () 11 min mod , 1 min mod ( , ) min , 2 tt t t t xx x x x FF F F W xx a b signa signb a b α +− ⎛⎞ −− = ⎜⎟ ΔΔ ⎝⎠ =+ where t x tt x x t x H FU V ⎛⎞ ⎜⎟ ≡ ⎜⎟ ⎜⎟ ⎝⎠ , 0.72 α = . The second order of accuracy in time is reached by application of two- step-by-step algorithm predictor - corrector. At a stage pr edictor there are found auxiliary values of required sizes for the whole step on time with the help of quasi-two-la yer algorithm of the first order of the accuracy. These auxiliary values are used to find the values on an intermediate step in time by us ing arithmetic averaging with va lues of the previous tim e 15 step. On the corrector step the given sizes are reconstructed in space: 2 1 2 t t x x Fx W τ + +Δ , 2 11 1 2 t t x x Fx W τ + ++ −Δ , accordingly at the left and at the right sides of the face 12 x + . Next, the values of fluid variables on borders of cells corresponding to an inter mediate time layer are found. 4. Analysis and comparison of the different methods for modeling of shallow water flows in presence the external force 4.1. Theoretical comparisions In this section we co mpare our algorithm with the other finite-difference schemes for numerical simulations of shallow water flows with in the presence of exter nal force in the form of a fictitious bed. To simplify the presentation, all deriva tions are performed in the cases of one-dimensional finite-differenc e schemes. We consider the case when the cell adjoins the step of height 0 b on the left side. We shall write down the finite-difference scheme for the wave-propagation method [40, 41] in the expl icit form. The method is based on solution of the auxiliary Riemann probl em at the center of each cell where t he difference of fluxes is selected in such a wa y that the external force effect is exactly balanced in the stationary case. Thus, at t he predictor step, the flow parameters , x x hu − − and , x x hu ++ are calculated inside a cell. Then, at the corrector step, using the flow parameters obtained at the predictor step, the flow parameters on cell’s interfaces 12 12 , ii hu −− and 12 12 , xx hu ++ are calculated (Fig. 2). In fact, such an approach corresponds to the momentum and mass redi stribution inside the cons idered cell, in order to compensate the influence of a step. The pr esence of such a sub grid discontinuity ensures the description of the non-divergen t right-hand side of the system of equations (2) owing to the additional source of mome ntum. As a result, the standard system of shallow water equations over the smooth bed wit h modified fluxes through the interfaces of cells is considered. To obtain the finite-difference scheme, we integrate over each of cell parts obtained by fictitious partition of the whole cell at the center: () ( ) 22 22 0 11 () ( ) () ( ) 0 22 left right left right GG GG dh u dh u dh dh dxdt dxdt dt dx dt dx hu hu gh dxdt hu hu gh dxdt tx tx ⎧ ⎡⎤ ⎡⎤ ++ += ⎪ ⎢⎥ ⎢⎥ ⎣⎦ ⎣⎦ ⎪ ⎨ ∂∂ ∂∂ ⎡⎤ ⎡⎤ ⎪ ++ + ++ = ⎢⎥ ⎢⎥ ⎪ ∂∂ ∂∂ ⎣⎦ ⎣⎦ ⎩ ∫∫ ∫∫ ∫∫ ∫∫ ww ww ( 16) The application of the Ostr ogradsky-Gauss formulas transfo rms (16) to the form: 16 () () () () () () () () 12 12 11 12 12 11 1 1 12 12 12 12 11 12 12 1 22 12 0 1 2 xx xx tt tt tt tt xx x x tt t t xx x x tt t t xx xx tt tt tt tt xx x x t t x hd x hd x h d x h d x hu dt hu dt hu dt hu dt hu dx hu dx hu dx hu dx hu gh dt hu ++ =+ =+ = = −− ++ + + +− + − ++ =+ =+ = = −− + + +− − + +− + − = + −−+ ⎛⎞ ++ − ⎜⎟ ⎝⎠ ∫∫∫ ∫ ∫∫ ∫ ∫ ∫∫∫ ∫ ∫ 1 22 12 11 22 22 1 2 11 0 22 t t x tt tt xx gh dt hu gh dt hu gh dt + − ++ −+ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎛⎞ ⎪ ++ ⎜⎟ ⎪ ⎝⎠ ⎪ ⎪ ⎛⎞ ⎛⎞ ++ − + = ⎜⎟ ⎜⎟ ⎪ ⎝⎠ ⎝⎠ ⎩ ∫ ∫∫ (17) The use of the mean value theo rem and performance of integrat ion allows one to rewrite system (17) in the form: 11 12 12 12 12 11 11 22 22 12 12 12 12 12 12 2 0 1 111 0 2 222 1 111 2 222 11 () ( ) 22 1 22 ttt t t t t t xxx x x x x x tt tt t t t t xx xx x x x x tt t tt t xx x xx x tt t xx x H x H x Hx Hx H U H U HU x HU x H U x H U x H U gH H U gH b HU g H ττ ττ ++ ++ −− ++ ++ ++ + −− − + Δ+ Δ− Δ− Δ+ − = Δ+ Δ− Δ− Δ+ ++ − + + ⎛⎞ ++ + ⎜⎟ ⎝⎠ 22 2 00 0 1 0 22 2 2 tt t xx x bb b HU g H ττ + ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎛⎞ ⎛ ⎞ ⎛⎞ ⎛⎞ ⎛ ⎞ ⎪ +− − + − = ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎝⎠ ⎝⎠ ⎝ ⎠ ⎪ ⎝⎠ ⎝ ⎠ ⎩ (18) Taking into account that in t he cell center the st ationary problem is solved, such that 00 11 , 22 tt tt xx x x H Hb H Hb +− =+ =− , and in the small neighborhood ε : 0 xx UU +− = = [40,41], system (18) is transformed into the following finite-difference scheme: 12 12 12 12 1 22 12 12 12 12 2 1 12 12 12 1 0 1 2 1 / 2 tt tt xx xx tt xx tt t xx x tt tt t t t x x xx x x x t x t x HU HU HH X HU g H H U UH U g H X H H gH b τ τ −− ++ + −− − ++ ++ + + ⎧ ⎛⎞ − =+ × ⎪ ⎜⎟ ⎜⎟ ⎪ ⎝⎠ ⎪ ⎛⎞ ⎪ +− ⎪ ⎜⎟ ⎨ ⎜⎟ ⎪ ⎜⎟ =× − − − + ⎪ ⎜⎟ ⎪ ⎜⎟ − ⎪ ⎜⎟ ⎜⎟ ⎪ ⎝⎠ ⎩ (19) In particular case when one-dimensional pr oblem is considered and step height is 0 b , the second equation of the obtained finite-diffe rence scheme formally coincides with the 17 second equation of the quasi-two-layer finite -difference scheme (14) under the following condition: 0 12 2 xx b HH + = + (20) where the flux quantity 12 x H + is calculated by the quasi-two-layer method. Note that in the wave-pr opagation method the flow parame ters at the cell boundary 12 12 , xx HU ++ are calculated as the solution of the dam-break problem for the initial parameters on the left 0 , 2 x x b H U + and on the right 11 , x x H U + + sides of the discontinuity plane. In the quasi-two-layer method 12 x H + is calculated as the so lution of the dam-break problem for the initia l parameters on the left * , x x H HU − , and on the right 11 , x x H U + + sides of the discontinuity plane. The finite-difference scheme proposed in papers [3,9,10,44] also interprets the term containing the external force as some fi ctitious bed. The method proposed in papers [3,9,10,44] is based on the sch eme of hydrostatic reconstruction of the flows on cell interfaces, which uses the solutions of st eady states. In the one-dimensional case the finite-difference scheme is obtained from the classical scheme for a smooth bed by adding the source terms 22 12 12 11 , 22 tt xx gH gH +− −+ into the scheme. On the interfaces, however, the corrected flow parameters are used, which take into account the hydrostatic balance. The scheme looks as follows: 12 12 12 12 1 22 2 12 12 12 12 12 1 1 1 222 12 12 12 1 2 / 11 1 22 2 tt tt xx xx tt xx tt t t t tt xx x xx t t x x x x t tt t x xx x HU HU HH X HU g H H U H U UX H H gH gH gH τ τ −− ++ + −− − ++ + + + ++ − − + ⎧ ⎛⎞ − =+ × ⎪ ⎜⎟ ⎜⎟ ⎪ ⎝⎠ ⎪ ⎛⎞ ⎨ +− − ⎜⎟ ⎪ =× + ⎜⎟ ⎪ ⎜⎟ ⎪ −+ − ⎜⎟ ⎝⎠ ⎩ (21) Using the hydrostatic balance condition [3,9,10,44 ], we find that in the case of a step of height 0 b , on the right side of the considered cell 12 0 12 , tt tt x xx x H Hb H H +− −+ =− = . Then the finite-difference scheme will be re-written in the form of: 12 12 12 12 1 22 2 12 12 12 12 12 1 1 1 22 12 0 0 1 2 / 11 22 tt tt xx xx tt xx tt t t t tt xx x xx t t x x x x t tt x xx HU HU HH X HU g H H U H U UX H H gH gb gH b τ τ −− ++ + −− − ++ + + + + ⎧ ⎛⎞ − =+ × ⎪ ⎜⎟ ⎜⎟ ⎪ ⎝⎠ ⎪ ⎛⎞ ⎨ +− − ⎜⎟ ⎪ =× + ⎜⎟ ⎪ ⎜⎟ ⎪ −+ − ⎜⎟ ⎝⎠ ⎩ (22) 18 The second equation of scheme (22), as well as that of scheme (19), formally coincides with the second equation of sc heme (14), provided that: 12 0 xx H Hb + = + (23) where the flux value 12 x H + is calculated by the quasi-two-layer method. Note that in the method propos ed in papers [3,9,10,44] 12 12 , xx HU ++ is calculated as the solution of the dam-break problem for the initial parameters on the left 0 , x x H bU − , and on the right 11 , x x H U + + sides of the discontinuity plane, whereas in the quasi-two- layer method 12 x H + is calculated as the solution of the dam-break problem for the in itial parameters on the left * , x x H HU − , and on the right 11 , x x H U + + sides of the discontinuity plane. Generally speaking, the finite-difference scheme proposed in papers [3, 9, 10, 44] can be obtained in the same manner, as the LeVeque scheme [40,41]. Indeed, it is sufficient to consider, at the predictor step in the LeVeque method [40,41], the problem with additional source correspo nding to the well-balance pr inc iple as the stationary hydrostatic problem (Figure 2, wave-propagat ion algorithm [40,41] – the dotted line, and the well-balancing algorithm [3, 9, 10, 44] – the dash-dotted line). It is important to emphasize that in paper [5], devoted to numeric al simulation of the classical Saint-Venant system, a similar approach of mass redistribution inside a ce ll gave quite satisfactory results. Unlike the situation with a non-di fferentiable nonhomogeneity of the right-hand side, in the given method the introduction of additional discontin uit y decay inside a cell is not required. Definitely, after integration of corresponding diffe rential equations, the obtained balance equations contain the necessa ry term, which describe the work made of the adjacent bed slope in the explicit form. Taking into consideration that the finite-d ifference presentations for the momentum conservation law with regard to (20) and (23) formally coincide our quasi-two-lay er algorithm is extension both me thods up to the accuracy of the method of finding flow quantities on cell interfaces. In this model 12 x H + is calculated wit h using (for the considered problem) the flow depth * x H H − on the left side from a discontinuity and 1 x H + – on the right side. Here * H , depending on the flow conditions, is chosen as a function of current hydrodynamic values. One can assume, in conformity with the physical situation, the whole spectrum of values [ 0, 0 b ] – at flow running on the fictitious bed, and 0 [, ] x bH – in the opposite case. At the same time, in the considered models the correction of flow quantities on calculated cell in terfaces does not depend on the values of hydrodynamic quantities in adjacent cells. 4.2. Numerical comparisions Presence of the term describing the influe nce of external force at the right-hand side of the equations essentially imposes cons trains on to accuracy of computations of fluxes on mesh borders. Furthermore, in the case of fast changing exte rnal force, the rate of their formation becomes especially import ant because it is essential to reflect adequately the influence of external force at the each time step. In order to demonstrate 19 the efficiency of the proposed me thod we provide the result s of numerical experiments that show the rate of formation of adequat e fluxes of the mass and the impulse. We compare our results with those obtained by a well-balanced algo rithm from [3,9,10,44] for model computations of flow on step.. Initial parameters of the num erical experiments are the following: depth of the fluid is 1 m over the entire spatial ar ea; velocity of the fluid is 5 m/s to the left of the step and 0 m/s to the right of the step (F ig. 3). Step height is 1 m. The evolution of the flux of the mass and the flux of the impulse over the step are displa yed on Fig. 4-5. Black lines indicate the results obtained by the quasi-two-layer method. Grey lines indicate the results given by the well-balanc ed algorithm from[3,9, 10,44]. As one can see, the results obtained by the quasi-two-layer me thod approaches the appropriate solution in the area of applicability of the method. Howerver the result given by the well-balanc ed algorithm [3,9,10,44] requires more than 7 time steps before one attaines adequate value of fluxes. For such period of time external force approximated by steps can be changed more the once. Therefore, there may not be enough time to form adequate fluxes during approximated force variations in time. In the paper [45] it is shown that almo st all currently existi ng methods conserve energy while passing over the step. However, papers [8,45] take into account energy losses. We show below the match between our result s with those obtained in [8 ]. Let us analyze our method in comparison to t he well-balanced algorithm in this context [3,9,10,44]. A distinctive feat ure of the propos ed method is t he automatic consideration of the properties of the proc ess of the waterfall, i.e. the fluid flow on the step in which the fluid does not wet part of the vertical wall of t he step. In this case the dissipation is caused by the lack of interaction with the solid underly ing surface for some of the fluid. The work of the gravity force to change the relative dept h of the fluid is bala nced by the change in the horizontal component of the momentum of the fluid becaus e of this lack. Indicated work is partly used in changing of vertical mom entum of the flow that is neglected in the shallow water approximation. Thus the presence of dry zones in the vertic al part of the step indicates violation of the conditions of hydrostatic flow. The quasi-two-layer approach determines the size of the dry zone of the ve rtical component of the step. Consequently it gives an opportunity to figure out the amount of kinetic energy dissipation. The problem of shock wave propagation over t he step has been s olved analytically in paper [42]. Initial profile of the shoc k wave oncoming on t he step is showed in Fig. 6. Here δ is the height of the step, 11 1 2 2 () 2 Dg h h h h =+ is the velocity of the initial shock wave, 1 v is the initial velocity of the fluid to the left of the step and 0 v is the initia l velocity of the fluid to the right of the step. The solutions in which the total energy of the flow conserves on the step and the solutions in which the total energy is lost have been considered. The analysis of the problem of the flows which occurs while onc oming of the shock wa ves on the step has showed that under self-similar solutions of s hallow water theory this problem is always solvable. However, the problem does not have a uniquely solution. In paper [42] five qualitatively different types of stable self-similar solutions for this problem have been obtained ( Mass fluxes). The total energy of the flow conserves on the step in three of them (solutions types I, III, V). In the other tw o (solutions types II, IV) the total energy of the flow is lost on the step. The domains of existence of t hese solutions are plotted on the plane of the dimensionless ke y parameters (Fig. 7). The subd omains in which both self- similar solutions exist simultaneously have been allocated. Here * 0 / ll hh h = , * 0 / h δδ = , 0.66 α ≈ , 3.214 β ≈ . Unlike most existing methods (see [45] ) in which the energy conser ves, our approach takes into account the loss of energy while passing over the step in this particular case of shock wave propagation ov er the step. This is demonstrated in 20 comparison with the results obtai ned by the well-balancing algorithm [3,9,10,44]. Particularly in Fig. 8 the result s of solving the follow ing Riemann problem are shown. The initial flow parameters are the following: to the left from the step the fluid depth is 5.8 m and the velocity is 6 m/s, to t he right from the step the fluid depth is 3. 3 m and the velocity is 0 m/s and t he step height is 1 m. The re sult obtained by quasi-two- layer method of the second order accuracy is pl otted as the black line (solution type II is the rarefaction wave to the left from the st ep, stationary shock over the step and right shock wave to the right from the step) T he result obtained by we ll-balancing algorithm [3,9,10,44] of the second order accuracy is pl otted as grey line (solution type III is the rarefaction wave to the left fr om the step, stationar y shock over the step, left shock wave to the right of the step and right shock wave to t he right from the step). In this case there a non-uniqueness of the solution exists fo r given parameters as shown in paper [42]. Solution obtained by quasi-two-layer me thod is a solution with energy loss while passing over the step, and solu tion obtained by well-balancing algorithm [3,9,10,44] is a solution without loss of energy. Fig. 9 shows the results of solving of the following Riemann pr oblem. The initial flow parameters are as following: to the left from the step the fluid depth is 3 m and the velocity is 8.45 m/s, to the ri ght from the step the fluid depth is 1.5 m and the velocity is 0 m/s and the step height is 1 m. The result obtained by quasi-two-layer method of the second order accuracy is plotted as the bla ck line (solution type IV is stationary shock over the step and right shock wa ve to the right from the step) . The result obtained by well- balancing algorithm [3,9,10,44] of the second order accuracy is plotted as gray line (solution type V is stationary shock over the st ep, left shock wave to the right from the step and right shock wave to the right from the step). In this case there a non-uniqueness of the solution exists for given parameters as shown in paper [42]. Solution obtained by quasi-two-layer method is a solution with ener gy loss while passing over the step, and solution obtained by well-balancing algorithm [3,9,10,44] is a so lution without loss of energy. We have also solved fundam ental problem of hydrod ynamics. Riemann problem on a step has been solved usi ng the quasi-two-layer method. Obtained solutions were compared with the set of exact solutions of t he same problem, obtained in paper [8]. In [8] the solutions have been obtained by solving enlarged system t hat includes an additional equation for the bottom geometry and then the principles of conservation of mass and momentum across the step were used. In order to exclude the multiplicity of solutions, one should impose that the entropy condition is fulfilled; so, to tal energy dissipates across the stationary shock wave at t he step and transition from subcritical to supercritical flow across an upward step is excluded. The quasi-two-layer method of the second order accuracy is applied for numerical modeling. Fig. 10 shows the solution of t he Riemann problem on a step with the following initial conditions: to the left fr om the step fluid hei ght is 1.462 m, Froud number is 0, to the right from the step fluid height is 0.309 m, Froud number is 0, step height is 0.2 m. Fig. 11 shows the solution of the Riemann problem on a step with the following initial conditions: to the left fr om the step fluid height is 0. 569 m, Froud num ber is 0.9, to the right from the st ep fluid height is 0.569 м , Froud number is 0, step height is 0.2 м . Good match up between the solutions obtai ned by the quasi-two-layer method and exact solutions obtained in paper [8] are shown. 5. Results of numerical simulations In order to verify the efficiency of the proposed method a several problems corresponding to different types of nonhomogeneit ies of the bed and different types of 21 external forces has been solved by using s uggested algorithm. In particular, horizontal and sloping planes (section 5.1, 5.4), and a su rface of parabolic type were considered as beds (section 5.3). Coriolis force (sections 5.2, 5.3) and a hydraulic friction force (section 5.4) are used for modeling of external forces . The results of numerical computations are compared with available data of laboratory experiments. Inte raction of Tsunami waves with the shore line including an obstacle is considered in section 5.5. 5.1. Dam-break problem of fl uid column of parallelepiped shape on a sloping plane Let us assume that at the center of an infi nite plane, situated at some angle to the horizon, k , exist the fluid column of parallelepipe d shape. Thus, one cons iders the effect of the constant external force that is equal to: Eg k = (24) For obtaining the numerical solution, we approximate the in clined plane by the system of ledges with a constant step in space. The selected problem of the fluid column breakdown is most indicative for checking t he method efficiency, since it combines in itself all possible types of solutions simultaneously. The quasi-two-layer method of the first or der accuracy is applied for numerical modeling. Comparisons were performed for t he cases of resting fl uid and fluid flowing upwards along the inclined pl ane. The spatial step in our computation is taken 0,012 m per cell. The following initial data were taken for comparisons: the size of a discontinuity of the fluid column with height of 1.5 m and square base with side of 2 m over the incline d plane (the inclination value is 1:20), covered with fluid. In the first numerical experiment the fluid was at rest, and in th e second one it was fl owing onto the plane at velocity of 2 m/s. In both numerical exper iments the in itial depth was 1 m everywhere outside the region occupied by a breaking column of fluid. The in itial flow parameters are shown by a dash-doted line. The dashed line on plots indicates the flow depth and velocity obtained by using the proposed hydrodynamic model, the black line – so lutions of combination of two Riemann problems on the sl ope [26] (Figs. 12, 13). The presented plots (Figures 12, 13) reveal almost complete analogy between the results obtained by means of the quasi-two-layer met hod and the solution of the combination of two Riemann problems on the slope [26]. This demonstrates the model efficiency in description of such natural phenomena. Oscillations on the fr ont of hydrodynamic jump occur due to intensity of this hydrodynamic jump. They are within method error and do not result in violating stability with due time. 5.2. Rotating fluid flow. The classical Rossby problem The quasi-two-layer method of the second order accuracy is applied for numerical modeling on space and in time. T he rectangular grid of the size 200 Х 10 cells is used. The classical Rossby problem is simulated as the test one [46]. The initial disturbance considered was 0 (, 0 ) (, 0 ) 0 (, 0 ) ( ) jet hx h ux vx V v x ⎧ = ⎪ = ⎨ ⎪ = ⎩ , (25) 22 where 0 h is the initial depth, V is the characteristic scale of velocity, () jet vx is the normalized profile specified as follows: 2 ( 1 tanh( 4 / 2))( 1 tanh( 4 / 2)) () (1 t a n h ( 2 ) ) jet xL xL vx ++ −− = + . L is the characteristic scale of di sturbance. Characteristic parameters 0 ,, gh f are fixed. The characteristic scale of velocity V and the characteristic scale of disturbance L are calculated from two dimensionl ess parameters: the Rossby-Kibe l ( R o ) and Burgers ( B u ) numbers: 2 2 , d VR Ro Bu f LL == , where d R is the deformation radius: 0 d gh R f = . The characteristic time scale is specified by the following formula: 2 f T f π = . The results of evolution of depth 0 h , in the case of 1 R o = , 0.25 B u = , are presented below. Figure 14 shows the evolution obtained by using of the quasi-two-layer model. One can see good coincidence of characteristic peaks of running-away acoustic -gravitational waves and the central balanced part with results presented in paper [10]. This testifies to the efficiency of using the quasi-two-layer model in the description of large-scale geophysical phenomena. Figure 15 shows the compar ison of potential vorticity values at the initia l ( 0 f tT = ) and final ( 16 f tT = ) time instants for the cl assical Rossby problem and 1 R o = , 0.5 B u = . The potential vorticity is defined the following formula: vu f xy Q h ∂ ∂ + + ∂∂ = (26) One can see that the invariant Q – the potential vortex is conserv ed with time. Note that the real time of the process equal s to 12 days, approximately. It is seen from the presented plot that the maximum of a function is shift ed to the anticyclonic region, and the minimum of the potential vorticity in creases with time. The given results are determined by purely nonlinear effects and ma tch well with those obtained in paper [39]. 5.3. Flow of a rotating fluid over mounted parabolic profile Problem was simulated, which cont ained both the nonhomogeneous bed, and the external force – the Coriolis force. The quasi-two-laye r method of the second order accuracy is applied for numerical modeling on space and time. The shallow water equations in this case are as follows: 23 () () () () () () 22 22 0 2 2 hu hv h tx y hu gh hu huv db gh fvh tx y d x hv gh hv hvu d b gh fuh ty x d y ∂∂ ⎧ ∂ ++= ⎪ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − + ⎨ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − − ⎪ ∂∂ ∂ ⎩ (27) where f is the Coriolis parameter . According to (27), we get the external forces: x Ef v = and y E fu =− . The large-scale motion of fluid in the gravitational field in the presence of Coriolis force over the mountain-like bed is considered. Typical parameters of the problem are: the linear dimensions 6 10 by 6 10 m, the mount ain height is 3 1, 2 1 0 ⋅ m, the fluid (depth is 3 21 0 ⋅ m, and the Coriolis parameter is 0.00001452 s --1 . The initial wind parameters are 0 u = m / s, 20 v = m / s. The boundary free-slip conditions are satisfied. The computational domain contains 3600 cells (60 × 60 cells). According to (5) the external forces x Ef v = and y Ef u =− , the corresponding differences of average heights of a fictitious surface on the boundary of the two adjacing cells are: 22 rl xx rl EX g E X g f vX g f v X g + + = , 22 rl yy rl EX g E X g f uX g f u X g + −− = . As a result, we found that the characterist ic time of one revolution of a system as a whole is 24 hours, which corresponds to th e natural phenomenon (the ch aracteristic time of one revolution of a system as a whole fo r the geophysical dynamics problems equals to one day [16]). Figure 16 shows the picture of fl ow at the initial time instant; Figure 17 shows the flow evolution during 24 hours. 5.4. Two-dimensional dam- break problem taking into account a hydraulic friction The computational domain repr esents a rectangle of size 3 Х 2 m, restricted by non- leaking walls on three sides. There is a wall apart 1 m from the le ft non-leaking boundary with a hole of width 0.40 m which is symmetrical with respect to the abscissa axis. The thickness of the wall is a negligible quantity and it is not used in calc ulations. A hole in a wall is initially closed and al l the fluid rests on the left of it. The fluid can freely leak through the right boundary of the computational domain. Define the Cartesian coordinate system in the following manner: the wall be longs to the ordinate axis, and the abscissa axis divides the computational domain in two symmetrical parts. The origin of coordinates coincides with the centre of the hole in a wall. All r equirements of the performed numerical modeling are taken according to requirements of l aboratory and numerical experiments from [17]. The slope angle of a bed in the perf ormed numerical and laboratory experiments varied fr om 0 % to 10 %. Coordinates of control points and the plan of their location are shown on Fig. 18. Laboratory installa tion is made from plexiglass and supplied by the mechanism for opening a shutte r of the hole in t he wall at a velocity sufficient that dynamics of uprise of the shutte r was not reflected in a pattern of formed flow. Detailed description of laboratory inst allation and the numeric al experiments can be fined in the paper of Fraccarollo and Toro [17]. The shallow water equations in this case are follows: 24 ( ) () () () () () 22 22 0 2 1 2 2 1 2 hu hv h tx y hu gh hu huv db gh u u h tx y d x hv gh hv hvu db gh v v h ty x d y λ λ ∂∂ ⎧ ∂ ++= ⎪ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − + ⎨ ∂∂ ∂ ⎪ ⎪ ∂+ ∂ ∂ ⎪ ++ = − + ⎪ ∂∂ ∂ ⎩ (28) where 4 2 3 2 gn h λ − = , 0.007 n = . External force for system (28) have the form of hydraulic friction: 1 2 x Eu u λ = , 1 2 y Ev v λ = . According (5) the corresponding differences of average heights of a fictitious surface on the boundary of the two neighboring cells are: 11 22 22 rl rr l l xx uu X g u uX g EX g E X g λλ + + = , 11 22 22 rl rr l l yy vv X g v vX g EX g E X g λλ + + = . The quasi-two-layer method of the second order accuracy is applied for numerical modelling on space and time. Results for the ca se of horizontal and sloping beds are shown below. The slope angle is 6.3 degrees. Dept h of a fluid to the left of a wall is 0.6 м at initial time moment, the fluid was at rest for all the results given below. The computational domain is divided by regular grid of size 150 × 50. Plots of fluid depths and velocities in a point 0 corresponding to the centre of dam- break are shown on Fig 19. The left plot corresponds to fluid depths, right – fluid velocities. Thin grey line corresponds to the data obtained in laborato ry experiment [17], dashed grey line - to numerical results of WAF method [17], heavy black line – t o numerical results of propos ed quasi-two-layer method. Dynamics of fluid depth for a horizontal bed in control poin t P1 is shown on Fig. 20: obtained experimentally – thin grey line, and num erically on the basis of WAF method [17] – dashed grey line, heavy blac k line – on the basis of prop osed quasi-two-layer method, thin black line – on the basis of proposed quasi-two- layer method without bed friction. One can see that numerical resu lts obtained without bed friction are slightly different from those obtained with bed friction in control point P1. Control point P1 is located very closely to localization domain of the whole initial mass of the fluid. In this particular case bed friction is important. Indeed an integral bed friction effect of the whole domain where the fluid spread affect on flow parameters in control point P1 in spite of low bed friction effect in local control point P1. Note that this integral effect of bed friction on flow parameters is the maximal distinctive in domai n where the whole initial mass of the fluid is. This total effect decreases as the cont rol point moves away from this domain and becomes negligible. Dynamics of fluid depth for a horizontal b ed in two control points are shown on Figs. 21-22: obtained experim entally – thin grey line, and numerically on the basis of WA F method [17] – dashed grey lin e, heavy black line – on the basis of proposed quasi-two- layer method. Dynamics of fluid depth for a sloping bed in two control po ints are shown on Fig. 23. Thin lines correspond to t he data obtained in laboratory experiment: dashed – point P2, 25 solid – point 0. Numerical results of proposed quasi-two-layer method shown by heavy lines: dashed – point P2, solid – point 0. On Fig. 24 numerical results of WAF method for the same points are shown. One can see t hat proposed quasi-two- layer method is in good agreement with experimental data and improves accuracy of results in comparison with the numerical results of WA F method obtained in article [17]. 5.5. Tsunami waves simulation Interaction of the Tsunami wave with the shore line with and without an obstacle has been simulated to demonstr ate the effectiveness of t he developed algorithm in domains including partly flooded and dry regions . The computational domain is divided into regular grid of size 1200 × 1. Slope angle of the bed is 0.125. Initial height of the Tsunami wave is 12 m, velocity is 5 m/s. In teraction of the Tsunam i wave with the shore line is shown on Fig. 25. Dynamics of fluid depth in various tim e steps (indicated on the wave front in seconds) is shown on the upper figure. Dynamics of fluid velocity in various time steps (indicated on the wave front in seconds) is shown on the lower figure. The significant increase of fluid velocity and Zu nami wave propagation for the distance about 900 m is observed during interaction of the Tsunami wave wit h the shore line. Interaction of the Tsunami wave with the shore includin g an obstacle of size 12 m located 100 m away from the coast has been simulated. Dy namics of fluid depth in the various time steps (indicated on the wave front in seconds) is shown on Fig. 26. Dynamics of fluid velocity in various time steps (indicated on the wave front in seconds) is shown on Fig. 27. The considerable decrease of fluid velo city after passing an obstacle and Zunami wave propagation for the distance about 750 m is observed. The decrease of destroying effect of the Tsunam i wave on the shore line by an obstacle is shown. Scaled- up area with an obstacle on the s hore is shown on Fig. 28. On the left is the dynamics of the fluid depth in various time steps (indicated on the obstacle in seconds). On the right is the Dynamics of the fluid velocity in the va rious time steps. The reflected shock wave propagating away from the obstacl e and some of the fluid fl owing around th e obstacle are observed. Thus effectiveness of the devel oped algorithm in domains including partly flooded areas and dry regions is demonstrated. 6. Conclusion. New numerical method for modeling of shal low water flows over an arbitrary bed profile in the presence of external force is proposed. The presentation of arbitrary external force by a fictitious bed along with the bed relief features is approximated by a time- dependent stepwise bed. This allows to apply the finite volume methods based on the Riemann problem solution for determination t he flow quantities. To the Riemann problem we use an advanced algorithm which is based on t he two-layer representat ion of fluid flow near the step. Such representation allows to obtain in suggested new finite difference scheme the additional terms that provide mechanical work made by nonhomogeneous bed interacting with flow. This also modifies properly values of hydrodynamic fluxes in steps vicinity. The calculation of flows in the proposed algor ithm is performed on the basis of the quasi-two-layer shallow water model , which represents the generalization of the classical single-layer model in relation to the initial system of Euler equations. It is shown that the developed finite-difference scheme bel ongs to the well-balance class. Unlike the known schemes of well-balanced type, the proposed scheme is adapted to flow parameters. The developed approach allows one to restore the flow structure throughout the space-time domain with consideration for it s vertical peculiarities near t he step. This makes it possible to determine the transve rsal component of velocity vector more correctly for the class of ex ternal forces, which determine essentially two-dimensional 26 problems, such as the problems with the Cori olis force. We will benefit this quality in future developments. The calculation of e ssentially two-dimensional problems wit h a corrected transversal velocity component w ill be accomplished in a sep arate work. A distinctive feature of the pr oposed method is the considerati on of the properties of the process of the waterfall, i.e. the fluid flow on the step in which the fluid does not wet part of the vertical wall of the st ep. The presence of dry zones in the vertical part of the step indicates violation of the c onditions of hydrostatic flow. The quasi-two-layer approach allows to determine the size of the dry zone of the vertical component of the step. Consequently it gives an opportunity to figure out the amount of kinetic energy dissipation. In order to demonstrate the efficiency of the proposed method we have demonstrated results of the com putations that show the rate of formation of an adequate flux of the mass and the impulse. We compare our results with those obtained by the well- balancing method. The Riemann problem on a step has been solved and good agreement of the obtained solutions and exact so lutions is shown for non-waterfall case. The applicability of the method is demonstrated by the example of solution of the problem of breakdown of a moving fluid column on an in clined plane. The motion of fluid in the presence of Coriolis force over a bed of the given profile is simulated. Reliability and validity of our results is shown on the basis of the solution of two- dimensional dam-break problem taking into account a hydraulic fric tion and comparison of the obtained solutions with available laborator y experiments results. Interact ion of the Tsunami wave with the shore line including an obsta cle has been simulated to dem onstrate the effectiveness of the developed algorithm in domains with partly flooded and dry regions. Acknowledgements The authors are grateful to Rupert Klein fo r reading manuscript and useful comments. References 1. Alcrudo F, Benkhaldoun F. Exac t solutions to the Rie mann problem of the shallow water equations with a bottom step, Comput. Fluids. 30 (2001) 643–671. 2. Alcrudo F, Garcia-Navarro P. A high resolution Godunov-type scheme in finite volumes for the 2D shallow water equation, In t. J. Numer. Meth. Fluids 16 (1993) 489–505. 3. Audusse E, Bouchut F, Bristeau M.-O, Klein R. & Perthame B. A Fast and Stable Well-Balanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows. SIAM J. Sci. Comp. Volume 25 (2004), Issue 6, 2050-2065. 4. Audusse E, Bristeau, F. A well-balanced positivity preserving ‘‘s econd-order’’ scheme for shallow water flows on unstructu red meshes. Journal of Computational Physics 206 (2005) 311–333 5. Belikov V. V, Semenov A. Yu . A Godunov-type method based on an exact solution to the Riemann problem for the shal low water equations, Computational Mathematics and Mathematical Physics, Vol. 37, No. 8 (1997) 974-986. 6. Benkhaldoun F, Elmahi I, Seaid M. Well-balanced finite volume schemes for pollutant transport by shallow water equations on unstructured meshes. Journal of Computational Physics 226 (2007) 180–203. 7. Bermudez A, Vazquez M. E. Upwind me thods for hyperbolic conservation laws with source terms, Comput. Fl uids 23 (8) (1994) 1049–1071. 8. Bernetti R, Titarev V. A, Toro E. F. Exact solution of the Riemann problem for the shallow water equations with discont inuous bottom geometry, Journal of Computational Physics (2007), doi: 10.1016/j.jcp. 2007.11.033. 9. Botta N, Klein R, L angenberg, Lützenkir chen S. Well balanced finite volume 27 methods for nearly hydrostatic flows. Jo urnal of Computational Physics, 196 (2004), 539-565. 10. Bouchut F, Le Sommer J, Zeitlin V. J. Frontal geostrophic adjustment and nonlinear-wave phenomena in one dimensional rotating shallow water. Part 2: high-resolution numerical simulations . Fluid Mech . 514 (2004), 35-63. 11. Caleffi V, Valiani A, Be rnini A. Fourth-order balance d source term treatment in central WENO schemes for shallow wate r equations. Journal of Computational Physics 218 (2006) 228–245 12. Castro M, Gallardo J. M, Pares C. High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow wate r systems, Math. Comput. 75 (2006) 1103– 1134. 13. Castro M. J, Garc ı a J. A, Gonzalez-Vida J. M, Mac ı as J, Pares C, Vazquez- Cendon M. E. Numerical simulation of two-layer shallow water flows through channels with irregular geometry, Journal of Comput ational Physi cs 195 (2004) 202–235. 14. Cea L, Vaz quez-Cendon M. E, Puertas J. Depth Averaged Modelling of Turbulent Shallow Water Flow with We t-Dry Fronts. Arch Comput Methods Eng (2007) 14: 303–341. 15. Cea L, Vazquez-Cendon M. E, Puertas J, Pena L. Numerical treatment of turbulent and mass source terms in the shallow wa ter equations for channels with irregular section. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS) 2004. 16. Dolzhansky F. V. Lectures on geophysica l hydrodynamics, Moscow. Institute of Numerical Mathematics R AS, 2006 – ISBN 5-901854-08-X 17. Fraccarollo L, Toro E. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. Journal of hydraulic research, 33 (1995) No. 6, 843-863. 18. Gallouet T, Herard J. M, Seguin N. Som e approximate Godunov schemes to compute shallow-water equations with topography, Comput. Fluids 32 (2003) 479– 513. 19. George D. L. Augmented Riemann Solvers for the Shallow Water Equations over Variable Topography with Stead y States and Inundation, J ournal of Computational Physics (2007), doi: 10.1016/j.jcp. 2007.10.027 20. George D. L. Finite volume me thods and adaptive refinement for tsunami propagation and inundation. Ph D thesis, University of Washington, 2006. 21. Gosse L. A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms, Comput. Math. Appl. 39 (2000) 135–159. 22. Greenberg J. M, Leroux A.-Y. A we ll-balanced scheme for the numerical processing of source terms in hyperbol ic equations, SIAM J. Numer. Anal. 33 (1996) 1–16. 23. Guinot V. An approximat e two-dimensional Riemann solv er for hyperbolic systems of conservation laws. Journal of Computational Physics 205 (2005) 292–314. 24. Hu K, Mingham C. G, Causon D. M. Nu merical simulation of wave overtopping of coastal structures using the non-linear shallow water equations, Coastal Engineering. 41 (2000) 433– 465. 25. Hubbard M. E, Garcia-Navarro P. Flux difference splitting and the balancing of source terms and flux gradients, Journal of Computational Physi cs. 2001. 165, 89 –125. 26. Karelsky K. V, Papkov V. V, Petrosyan A. S. The initial discontinuity decay problem for shallow water equations on slopes, Ph ys. Let. A. 2000. V. 271. 349-357. 28 27. Karelsky K. V, Petrosyan A. S. Problem of steady-state flow over a step in the shallow-water approximation, Izvestiy a RAN, Mechanics of Fluid and Gas, № 1, 2006, 15-24. 28. Karelsky K. V, Petrosyan A. S. Part icular solutions and Riemann problem for modified shallow water equations. Fluid Dynamics Research, 38 (2006) 339–358. 29. Karelsky K. V, Petrosyan A. S, Slavin A. G. Riemann problem for shallow water flows on step. Geophysical Re search Abstracts, Vol. 7, 08202, Vienna, Austria, 2005. 30. Karelsky K.V, Petrosyan A.S, Slavin A.G. Quasi-two-layer model for numerical analysis of shallow water flows over a step boundary. The Fi fth International Symposium on Environmental Hy draulics. Arizona, 2007. 31. Karelsky K. V, Petrosyan A. S, Slavin A. G. Quazi-two-layer model for numerical analysis shallow water flow s on step, Russian journal of Numerical Analysis and Mathematical modeling. V. 21, № 6, 2006. 539-559. 32. Karelsky K. V, Petrosyan A. S, Slavin A. G. Numerica l simulat ion of flows of a heavy nonviscous fluid with a free surface in the gravity field over a bed surface with an arbitrary profile, Russi an journal of Numerical A nalysis and Mathematical modeling. V. 22, № 6, 2007. 543-565. 33. Kirwan A. D. Jr, and Liu Juping. The Shallow-Water Equations on an F-Plane. Proceedings of the internati onal school of physics ”Enric o Fermi”, Nonlinear Topics in Ocean Physics, 1988. 34. Kirwan A. D. Jr, Mied Ri chard P, and Lipphar dt B. L. Jr. Rotating modons ov er isolated topography in a two-layer ocean, Z. angew. Math. Phys. 48 (1997) 535- 570. 35. Kizner Z, Khvoles R, Mc Williams J.C. Rotating multi poles on the f and - planes. Physics Of Fluids, 19, 016603 (2007). 36. Kolgan V. P. Application of the principle of minimum values of the derivative to the construction of finite-difference schemes fo r calculating discontinuous solutions of gas dynamics, Scientific No tes of TsAGI, 3 (1972), 68–77. 37. Kulikovskiy A. G, Pogorelov, N. V, Se menov A. Yu. Mathem atical Aspects of Numerical Solution of Hyperbolic S ystems. Chapman & Hall/CRC, London/Boca Raton (2001), 560p. 38. Kuo, A. C, and Polvan i L. M, Time-dependent fu lly nonlinear geostrophic adjustment. J. Phys. Oceanogr ., 27 (1997) 1614–1634. 39. Le Sommer J, Medvedev S.B, Plougonven R. Zeitlin V Singularity formation during relaxation of jets and fronts toward the state of geostro phic equilibrium. Communications in Nonlinear Science a nd Numerical Simulation, Volume 8, Number 3, Septem ber 2003, 415-442 (28). 40. LeVeque R. J. Ba lancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-p ropagation algorithm, Journal of Computational Physics. 1998. 146, 346 –365. 41. LeVeque R. J. Wave-pr opagation algorithms for mult idimensional hyperbolic systems. Journal of Computat ional Physics 131 (1998) 327-353. 42. Ostapenko V.V, Shinkarenk o E.V. Flows arising after shock waves passage over the step. Mechanics of Fluid and Gas, № 1, 2009, 106-122. (in Russian). 43. Pares C, Castro M. On the well-balanced property of Roe’s method for nonconservative hyperbolic systems. App lications to shallow water sy stems, ESAIM: M2AN 38 (2004) 821–852. 44. Reznik G. M, Zeitlin V, and Ben Jelloul M. Non linear theory of geostrophic adjustment. Part 1. Rotating shallow-wate r model. J. Fluid Mech. (2001), vol. 445, 93-120. 45. Rosatti G, Begnudelli L. The Riemann Problem for the one-dimensional, free- 29 surface Shallow Water Equations with a bed step: Theoretical analysis and numerical simulations. Journal of Computational Physics. 229 (2010) 760 –787. 46. Rossby C. G. On the mut ual adjustment of pressure and velocity distributions incertain simple current systems. J. Mar. Res. (1938) 1, 239–263. 47. Slavin A. G, Karelsky K. V, and Petrosy an A. S. A quasi-two-layer model for flows of shallow water over a stepped boundary. In : Proc. XLVII scientific conference of MFTI. Part VIII. Physics and Ener getics. Moscow, (2004) 30–32. 48. Slavin A. G. Hydrodynamics of nonviscous heavy fluid with a free surface over a bed surface with a complex profile. Proceedings of t he XXVIII Conference of Young Scientists, Departament of Mec hanics and Mathemat ics, Moscow State University, (2006) 188-192. 49. Stoker J. J. Water Wave s: The Mathematical Theory with Applications. New York: Interscience, 1957. 50. Toro E. F. Riemann Solvers and Nume rical Methods for Fluid Dynamics. A Practical Introduction, 2nd ed., Springer-Verlag, Berlin, 1999. 51. Vazquez-Cendon M. E. Improv ed treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geome try. Journal of Computational Physics 148 (1999) 497–526. 52. Vignoli G, Titarev V. A, Toro E. F. ADER schemes for the shallow water equations in channel with irregular bottom elevatio n, Journal of Comp utational Physics (2007), doi: 10.1016/j.jcp.2007.11.006 53. Vreugdenhil C.B. Numerica l methods for shallow-wate r flow. Dordrecht: Kluwer; 1994. 54. Zhou J. G, Causon D. M, Mingham C. G, Ingram D. M. The surface gradient method for the treatment of source terms in the shallow water equations. Journal of Computational Physics 2001; 168:1 –25. 55. Zoppou C, Roberts S. Numerical solu tion of the two-dim ensio nal unsteady dam break, Appl. Math. M odel. 2000. 24, 457–475. Figure captions Fig. 1. Quasi-two layer repres entation of shallow flow on step 30 Fig. 2. Scheme of updating flow parameters in the Wave propagation [41,42] method (dashed line) and in the we ll-balancing [3,9,10,44] method (dash-dotted line) Fig. 3. Initial parameters of the fluid flow for different finite difference schemes comparisions 31 Fig. 4. Comparisions for the mass fluxes in time over the step. Black lines indicate the results obtained by the quasi-two-layer method. Grey lines indicate the results given by the well-balancing algorithm. Fig. 5. Time evolution of the flux of the im pulse over the step. Black lines indicate the results obtained by the quasi-two-layer method. Grey lines indicate the results given by the well-balancing algorithm. 32 Fig. 6. Initial condition of the shock wave oncoming on the step (from paper [42]) Fig. 7. Diagram of the dimensionless key parameters characterizi ng non-uniqueness of flow on step (from paper [42]) 33 Fig. 8. Fluid depths during pr opagation of the flow through t he step at time 0.4 s. Black lines indicate the results obtained by the quas i-two-layer method. Grey lines indicate the results given by the well-balancing algorithm. Fig. 9. Fluid depths during pass age of the flow throug h the step at time 0.4 s. Black lines indicate the results obtained by the quasi-two-layer method. Grey lines indicate the results given by the well-balancing algorithm. 34 Fig. 10. Comparison of quasi-two-layer solution (left) with the exact solution (right ). Configuration: the left rarefacti on wave to the left from the step, the right shock to the right from the step. 35 Fig. 11. Comparison of the quasi-two-layer so lution (left) with the exact solution (right). Configuration: the left shock to the left from t he step, the right shock to the right from the step. Fig. 12. Depth and velocity of fluid flow over an inclined bed, t = 0.7 s in the cross section by the plane of symmetry, y = 0. Dash-doted line - initia l flow parameters; dashed line - 36 flow depth and velocity obtained by using th e proposed hydrodynami c model; black line – flow depth and velocity obtained by decisi on of the Riemann pr oblems on the slope Fig. 13. Depth and velocity of fluid flow over an inclined bed with fluid leakage, t = 0.7 s in the cross section by t he plane of symmetry, y = 0. Dash-doted line - initial flow parameters; dashed line - flow depth and velo city obtained by using the proposed hydrodynamic model; black line – flow depth and ve locity obtained by decision of the Riemann problems on the slope 37 Fig. 14. Results of computati ons of propagation of acoustic- gravitational waves using the quasi-two-layer method Fig. 15. Potential vort icity computations usi ng quasi-two-layer method at the initial 0, 2 f tT = (black line) and at the final 16 f tT = ( dashed line). Fig. 16. Initial conditions of the test problem of the flow of rota ting fluid over mounted parabolic profile. 38 Fig. 17. Evolution of fluid (gas ) flow under the Coriolis force effect over a mountain. Up per plots are velocity fields, lower plots are the fr ee surface depth 39 Fig. 18. Computationa l domain and coordinates of contro l points of two-dimensional dam- break problem calculations Fig. 19. Left: fluid depths on time at point O. Ri ght: fluid velocities on time at point O. Thin grey line - data obtained in l aboratory experiment, dashed grey line - numerical results of WAF method, heavy black line – numerical re sults of proposed quasi-two-layer method 40 Fig. 20. Fluid depths on time at point P1. Thin grey line - dat a obtained in laboratory experiment, dashed grey line - numerical resu lts of WAF method, heavy black line – numerical results of propos ed quasi-two-layer method Fig. 21. Fluid depths on time at point P3. Thin grey line - dat a obtained in laboratory experiment, dashed grey line - numerical resu lts of WAF method, heavy black line – numerical results of propos ed quasi-two-layer method 41 Fig. 22. Fluid depths on time at point P4. Thin grey line - dat a obtained in laboratory experiment, dashed grey line - numerical resu lts of WAF method, heavy black line – numerical results of propos ed quasi-two-layer method 42 Fig. 23. Fluid depths on time for laborator y experiment with a slopi ng bed. Thin lines correspond to the data obtained in laboratory experiment: dashed – point P2, solid – point 0. Heavy lines correspond to numerical re sults of proposed quasi-two-layer method: dashed – point P2, solid – point 0 43 Fig. 24. Fluid depths on time for laborator y experiment with a sloping bed. Numerical results of WAF method: solid li ne – point 0, dashed line – point Р 2 44 Fig. 25. Results of numerical computations of interaction of the Tsunami wave with the shore line 45 Fig. 26. Results of numerical computations of interaction of the Tsunami wave with the shore line with obstacle.Dynam ics of the fluid depth in vari ous time steps (indicated on the wave front in seconds) 46 Fig. 27. Results of numerical computations of interaction of the Tsunami wave with the shore line.Dynamics of the fluid velocity in various time steps (indicated on the wave front in seconds) 47 Fig. 28. Results of numerical computations of interaction of the Tsunami wave with the shore line with obstacles. Scaled- up area with an obstacle on the shore