Abstract:- Free surface flow of water over a shallow rough bed is characteristically turbulent due to
disturbances generated by the bed resistance and diverse causes. The paper presents a derivation of the
basic equations in two dimensions and their numerical solution by an extension of the method devel-
oped earlier for flow in one dimension. Starting from the three dimensional Reynolds Averaged Navier
Stokes (RANS) equations, the equations of continuity and horizontal momenta are depth averaged to
derive three equations for the free surface elevation ζand the horizontal, depth averaged velocity com-
ponents (U, V ). Certain closure assumptions are required for derivation of the equations. Principally,
the viscous stresses are neglected, while the Reynolds stresses are assumed to depend on the vertical
coordinate zonly for the shearing flow over the x, y - plane representing the plane bed. Secondly, it
is assumed that the instantaneous horizontal components of velocity (u, v) follow the 1/pth (p= 7)
power law of variation in the z- direction, in liu of the the logarithmic law of the wall. For numerical
solution of the three nonlinear equations of continuity and momenta, the equations are reformulated in
terms of the primitive “discharge” components (Q, R) of the velocity (U, V ), showing that Qand R
can be functions of ζalone. The transformed equation of continuity is treated by the Lax-Richtmyer
method. The two momentum equations on the other hand, transform in to two coupled second degree
equations in the derivatives of Qand R, which decouple in the important case of quasilinear straight
crested waves on the water surface. The decoupled equations are numerically solved by the iterative
modified Euler method, and illustrated by application to an initial elevation of a model for tsunami
propagation.
Keywords: Two-dimensional shallow-water equations, Turbulence, RANS equations, Depth averaging,
Discharge, Lax-Richtmyer method, Modified Euler Method.
Received: May 30, 2021. Revised: November 28, 2021. Accepted: December 22, 2021. Published: January 19, 2022.
1 Introduction
Shallow water equations in two dimensions are
useful in oceanic and atmospheric flows as the as-
pect ratio of the horizontal spread to vertical height
is very large. In the development of the equations,
for oceanic flows such as tsunamis, tides, storm
surges and oceanic circulations, the fluid pressure
at a point is considered hydrostatic, even though
the free surface may have vertical motion. The gov-
erning equations of mass conservation and the hor-
izontal components of momentum are formally in
terms of the surface elevation ζand the depth aver-
aged components of the horizontal velocity (U, V ).
The equations neglecting the wind forses on the
free surface of the water are (Mader [1], p. 27), the
equation of continuity
∂ζ
∂t +∂(ζU)
∂x +∂(ζV )
∂y = 0 (1)
and the momentum equations in the x, y - direc-
tions as
∂U
∂t +U∂U
∂x +V∂U
∂y +g∂ζ
∂x +τ0x= 0 (2)
∂V
∂t +U∂V
∂x +V∂V
∂y +g∂ζ
∂y +τ0y= 0 (3)
where (τ0x, τ0y) are the components of bed resis-
tance that are usually modeled in the hydraulics
Nonlinear Turbulent Two-Dimensional Shallow-Water
Equations And Their Numerical Solution
SUJIT K. BOSE
S.N. Bose National Centre for Basic Sciences, Kolkata 700106, West Bengal, INDIA.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
34
Volume 17, 2022
literature by either the Chezy or by the Manning
empirical formulae of the following forms:
τ0x=gU√U2+V2
C2ζ, τ0y=gV √U2+V2
C2ζ,(Chezy)
(4a)
τ0x=m2gU√U2+V2
ζ1/3, τ0y=m2gV √U2+V2
ζ1/3,
(Manning) (4b)
where Cand mare respectively the Chezy and the
Manning coefficients that depend on the roughnes
of the bed, and gthe acceleration due to gravity
The nonlinear partial differential equations (1) -
(3) are also treated by Mader [1] using first order
finite differencing and applying them to a number
of applications, in particular to tsunamis.
Among the more advanced methods, Fennema
and Chowdhry [2], [3] have respectively given im-
plicit and explicit type of finite difference methods.
The former is of Beam-Warming type and the lat-
ter of MacCormack-Gobuti type, both being sec-
ond order accurate in space as well as time. On
the other hand, Wilders et. al. [4] and Casulli [5]
develop fully implicit splitting and semi-implicit fi-
nite difference methods for the equations. In recent
years, Ouyang et. al [6] present a MacCormack-
TVD finite difference scheme for solution, where
as Fu and Hodges [7] present a time-centered split-
implicit method, citing some other works. Much
earlier however, Katopodes and Strelkoff [8] treat
the equations differently by the method of char-
acteristic surfaces. Finite volume methods have
also been proposed. Alcrudo and Garcia-Navarro
[9] give a high resolution Godunov type scheme,
while Peng and Tang [10] present a method for
which a Godunov type scheme is used for solving
the Riemann problem. Garcia-Navarro et. al. [11]
have also presented a multidimensional upwind
scheme. Oceanic circulation models based on the
two dimensional shallow water equations are also
commercially available as software packages.
Here, in this paper, the shallow water equa-
tions in two-dimensions are developed taking in
to consideration the turbulence generated by bed
friction and other causes. The methodology is a
generalization of the one given by Bose and Dey
[12] for the case of one dimensional open chan-
nel flows. Beginning with the three dimensional
RANS equations, it is assumed that the Reynolds
stress components at a point depend on the verti-
cal zcoordinate of the point only in a flow, which
is essentially of shearing type parallel to the bed.
It is then argued following Schlichting [13] that
the horizontal components of the time averaged
instantaneous velocity (u, v) vary as z1/p where
p= 7. The proportionality constants appearing
in the expressions for uand vare then written in
terms of the respective depth-averaged horizontal
velocities (U, V ), which are functions of the hor-
izontal coordinates (x, y) only. The one-seventh
power law is used here instead of the better known
logarithmic law of the wall, as the former leads to
explicit depth averaged equations for the flow. An
expression for the time averaged vertical velocity w
in terms of U, V and zfollows from the equation of
continuity. Similarly, depth averaging of the verti-
cal component of the momentum equation subject
to the assumption of linear variation of the stream
surface with respect to zyields the expression for
the time averaged pressure ¯p, which contains apart
from the hydrostatic pressure, nonlinear terms on
account of vertical instantaneous and convective
accelerations. Depth averaging of the continuity
and the two horizontal momentum equations then
yield the depth averaged continuity and momen-
tum equations in terms of ζ, U and V. Numerical
treatment of the nonlinear equations is then de-
veloped generalizing that given by Bose (2018) for
one dimensional flows. Expressing the three equa-
tions in terms of “discharge” components Q=ζU
and R=ζV parallel to the xand yaxes, it is
shown that the two quantities can be functions of
ζalone, that is to say, Q=F(ζ) and R=G(ζ).
The two nonlinear momentum equations then re-
duce to two coupled second degree equations in
F0(ζ) and G0(ζ), where the prime denotes differ-
entiation with respect to ζ. In the important case
of quasi-linear nearly straight crested waves par-
allel to the y-axis, the equations decouple in to a
quadratic equation in F0(ζ) and an explicit equa-
tion for G0(ζ). The averaged continuity equation
on the other hand reduces to a standard form con-
servation equation , which is treated by the second
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
35
Volume 17, 2022
order Lax-Richtmyer method [16], where in the
first order pair of ordinary differential equations
fot F(ζ) and G(ζ) is solved by the second order
iterative modified Euler method (Bose [15]. The
method is illustrated by application to an initial
elevation form proposed in Howe [17], for a model
tsunami propagation problem.
2 The Turbulent 2D Shallow Water Equa-
tions
The two-dimensional turbulent shallow-water free-
surface flow over a plane horizontal bed is governed
by the general RANS equations. As the free surface
in general, can be undulating during the flow, the
plane bed is conveniently chosen as the x, y-plane
with the z-axis drawn vertically upwards, with a
mean water depth of hunits. Due to turbulence,
the velocity at a point (x, y, z) in the fluid consists
of the time averaged components (u, v, w) and in-
stantaneous fluctuations (u0, v0, w0). The time av-
eraged components satisfy the RANS equation of
continuity ∂u
∂x +∂v
∂y +∂w
∂z = 0 (5)
for incompressible fluids. The fluctuating veloc-
ity components satisfy a similar equation. The
RANS momentum equations (Schlichting [13], p.
562) contain contributions of (u, v, w), that of the
pressure ¯p, and the viscous and Reynolds stresses.
Moreover, as the motion over the horizontal bed
is essentially of shearing type, the viscous and
Reynolds stress components become practically in-
dependent of xand y. Hence, neglecting the contri-
bution of the normal viscous stress in the vertical
direction, the RANS momentum equations read as
∂u
∂t +u∂u
∂x+v∂u
∂y +w∂u
∂z =−1
ρ
∂¯p
∂x+1
ρ
∂τzx
∂z +µ
ρ
∂2u
∂z2
(6)
∂v
∂t +u∂v
∂x +v∂v
∂y +w∂v
∂z =−1
ρ
∂¯p
∂y +1
ρ
∂τzy
∂z +µ
ρ
∂2v
∂z2
(7)
∂w
∂t +u∂w
∂x +v∂w
∂y +w∂w
∂z =−g−1
ρ
∂¯p
∂z −∂w02
∂z (8)
where ρ= density, µ= dynamical coefficient
of viscosity, g= acceleration due to gravity, and
τzx =−ρw0u0, τzy =−ρw0v0are Reynolds stresses
in which the overbars denote time averages. The
Eqs. (5) - (8) describing the flow is under deter-
mined, and so additional conditions have to be
incorporated.
3 The 1/pth Power Law
For open channel flow in one direction only, say the
x-direction, it is known that uvaries slowly as z1/p,
where p= 7 (Schlichting [13], p. 590; Bose and Dey
[12]). Similarly, in the present generalization, one
can assume that
u=p+ 1
pU(x, y, t)z
ζ1/p,
v=p+ 1
pV(x, y, t)z
ζ1/p (9)
where the velocity (U, V ) equals the depth average
of (u, v), and ζ(x, y) is the surface elevation over
the point (x, y) on the bed. Accordingly, from Eq.
(5) by integration
w=−Zz
0∂u
∂x +∂v
∂y dz =−∂U
∂x +∂V
∂y z1/p+1
ζ1/p
(10)
where it is assumed for simplicity that the forward
flux (∂U/∂x, ∂V /∂y) of the flow is much greater
than the slope factor (U/ζ) (∂ζ/∂x, ∂ζ/∂y) of the
free surface, so that the latter terms become in-
significant compared to the term written in equa-
tion (10). Physically, the approxiation (10) implies
that as the rate of expansion ∂U/∂x +∂V/∂V in-
crases the vertical velocity wdecreases, which ap-
parently is quite plausible due to the mass conser-
vation principle. The solution forms (9) and (10)
are next used in equations (5) - (8). The convec-
tive acceleration term in Eq. (8), using Eq. (5)
becomes
u∂w
∂x +v∂w
∂y +w∂w
∂z =u∂w
∂x +v∂w
∂y −w∂u
∂x +∂v
∂y
=u2∂
∂xw
u+v2∂
∂y w
v=u2∂2z1
∂x2+v2∂2z1
∂y2(11)
wherein z=z1(x, y) represents the stream surface
passing through the point (x, y, z). The two par-
tial derivatives ∂2z1/∂x2and ∂2z1/∂y2in equation
(11) are proportional to the curvature of the sur-
face z=z1(x, y) and is assumed to linearly vary
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
36
Volume 17, 2022
with z, as in the Boussinesq approximation of open
channel flows. Then,
∂2z1
∂x2=z
ζ
∂2ζ
∂x2,∂2z1
∂y2=z
ζ
∂2ζ
∂y2(12)
Hence using Eq. (9), the convective acceleration
(11) becomes
u∂w
∂x +v∂w
∂y +w∂w
∂z
=p+ 1
p2z
ζ2/p+1 U2∂2ζ
∂x2+V2∂2ζ
∂y2(13)
Inserting the expression (13) in equation (8) and
integrating over depth, one gets
p
ρ=−Zz
ζ
∂w
∂t dz −p+ 1
2pζnz
ζ2/p+2
−1o
×U2∂2ζ
∂x2+V2∂2ζ
∂y2−w02−g(z−ζ) + ¯p0
ρ(14)
where ¯p0is the constant pressure over the free
surface z=ζ. Thus the formal expression for the
pressure is also obtained. In equation (14), the last
two terms represent the ontribution of hydrostatic
pressure and the rest represent the additional con-
tribution due to turbulence.
4 Depth Averaging
As the water depth is shallow, it is appropriate to
consider the depth averages of the continuity and
horizontal momentum equations (5) - (7). Integra-
tion of equation (5) with respect to zover (0, ζ)
yields the depth averaged continuity equation (1).
Treating the momentum equations (6) and (7) in a
similar manner, the former leads to the equation
∂
∂t(ζU) + ∂
∂xZζ
0
u2dz+∂
∂y Zζ
0
uv dz
+1
ρZζ
0
∂¯p
∂x dz +τ0x
ρ= 0 (15)
where τ0x= bed resistance parallel to the x-axis.
The representations (9) yield for the first two inte-
grals of equation (15) as
Zζ
0
u2dz =α ζU2,Zζ
0
uv dz =α ζUV (16)
where α= (p+ 1)2/[p(p+ 2)]. On the other hand,
1
ρZζ
0
∂¯p
∂x dz =∂
∂x Zζ
0
¯p
ρdz (17)
as ¯p(x, y, ζ) = 0 on the free surface. Hence using
Eq. (10) with Eqs. (5),
1
ρZζ
0
∂¯p
∂x dz =γ∂
∂xhζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i
−δ∂
∂x hζ3∂
∂t∂U
∂x +∂V
∂y
−ζ3
p
∂ζ
∂t ∂U
∂x +∂V
∂y i+g∂ζ
∂x (18)
in which, γ= (p+ 1)2/[p(3p+ 2)], δ =p/(3p+ 1),
and the Reynolds normal stress term ∂(w02)/∂x is
negligiible as assumed earlier. The contribution of
the bed shear term τ0x=−µ(∂u/∂z)z=0 in equa-
tion (15) can only be estimated . Here the Man-
ning’s formula (4b) is adopted as the coefficient m
appearing in the formula has lesser variation with
respect to variation with the roughness of the bed
surface. Equation (15) therefore becomes,
∂
∂t(ζU) + α∂
∂x(ζU2) + α∂
∂y (ζUV )
+γ∂
∂xhζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i
−δ∂
∂xhζ3∂
∂t∂U
∂x +∂V
∂y −ζ2
p
∂ζ
∂t ∂U
∂x +∂V
∂y i
+gζ ∂ζ
∂x +m2gU√U2+V2
ζ1/3= 0 (19)
In a similar manner Eq. (7) leads to the averaged
equation
∂
∂t(ζV ) + α∂
∂y (ζV 2) + α∂
∂x(ζUV )
+γ∂
∂y hζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i−
δ∂
∂y hζ3∂
∂t∂U
∂x +∂V
∂y −ζ2
p
∂ζ
∂t ∂U
∂x +∂V
∂y i
+gζ ∂ζ
∂y +m2gV √U2+V2
ζ1/3= 0 (20)
Equations (1), (19) and (20) are then the depth
averaged continuity and the x,ycomponents of
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
37
Volume 17, 2022
the momentum equation.
5 Simplification for Numerical Solution
Since p= 7 approximately, α= 64/63 ≈
1.01587 ··· ≈ 1, γ = 64/161 = 0.39751 ··· ≈ 0.4 =
2/5, and δ= 7/22. The approximation for αleads
to familiar convective acceleration terms of equa-
tions (2) and (3), as can be seen in the equations
(23) and (24) derived in the following. The ap-
proximation for γis introduced to bring some el-
egance to the basic equations. The highly nonlin-
ear fifth term in both of the equations (19) and
(20) arise due to the instantaneous vertical accel-
eratio ∂w/∂t. For developing the numerical solu-
tion method, it is first noted that the two terms are
horizontal gradients that need not be large. Hence,
one obtains by iteration of (19) and (20), the first
order approximations
∂U
∂t ≈ −g∂ζ
∂x,∂V
∂t ≈ −g∂ζ
∂y and
∂U
∂x +∂V
∂y ≈ −1
ζ
∂ζ
∂t (21)
The left hand side quantities of (21) are then in-
serted in the δ-terms of equations (19) and (20).
Thus equation (19) becomes
ζ∂U
∂t −ζU ∂U
∂x +∂V
∂y +∂
∂xhζnU2−1
22 ∂ζ
∂t 2oi
+∂
∂y (ζUV ) + 2
5
∂
∂xhζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i
+7g
22
∂
∂xhζ3∂2ζ
∂x2+∂2ζ
∂y2i+gζ ∂ζ
∂x
+m2gU√U2+V2
ζ1/3= 0 (22)
A further approximation is possible in equation
(22) if one asssumes that U2>> 1
22∂ζ
∂t 2
, as the
horizontal velocity is much larger than the instan-
taneous vertical velocity of the free surface. With
this assumption, equation (22) simplifies to
∂U
∂t +U∂U
∂x +V∂U
∂y +2
5ζ
∂
∂xhζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i
+7g
22 ζ
∂
∂xhζ3∂2ζ
∂x2+∂2ζ
∂y2i+g∂ζ
∂x
+m2gU√U2+V2
ζ/3= 0 (23)
Similarly equation (20) reduces to
∂V
∂t +U∂V
∂x +V∂V
∂y +2
5ζ
∂
∂y hζ2U2∂2ζ
∂x2+V2∂2ζ
∂y2i
+7g
22 ζ
∂
∂y hζ3∂2ζ
∂x2+∂2ζ
∂y2i+g∂ζ
∂y
+m2gV √U2+V2
ζ/3= 0 (24)
The continuity and momentum equations (1)
and (23) - (24) are in physical dimensions. In order
to nondimensionalisze the equations, let the vari-
ables x, y, t, ζ and U, V be temporarily written as
ˆx, ˆy, ˆ
t, ˆ
ζand ˆ
U, ˆ
Vrespectively; then writing
ˆx=hx, ˆy=hy, ˆ
ζ=hζ, ˆ
t=tph/g, and
ˆ
U=Upgh, ˆ
V=Vpgh (25)
in which his the mean depth of water; then the
variables x, y, ζ, t and U, V become nondimen-
sional. By the transformations (25), equation (1)
remains unchanged in form, while the factor g
drops out from the equations (23) and (24). Equa-
tion (23) with V= 0 and ∂(·)/∂y = 0 corresponds
to the one dimensional case studied by Bose [14],
but here it is slightly simpler in form.
6 Transformation and Localization Assump-
tion
In order to develop the numerical procedure
for solving the equations (1) and (23) - (24), it is
noted that the continuity equation (1) relates the
instantaneous vertical velocity of the free surface
at (x, y) to the gradient of the discharge (Q, R) in
the x, y - directions defined by
Q:= ζ U, R := ζ V (26)
so that equation (1) in compact substricted notoa-
tion for partial derivatives becomes
ζt+Qx+Ry= 0 (27)
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
38
Volume 17, 2022
In terms of Q, R one has
U=Q
ζ, Ut=Qt
ζ+Q
ζ2(Qx+Ry),
Ux=Qx
ζ−Q
ζ2ζx, Uy=Qy
ζ−Q
ζ2ζy(28)
V=R
ζ, Vt=Rt
ζ+R
ζ2(Qx+Ry),
Vx=Rx
ζ−R
ζ2ζx, Vy=Ry
ζ−R
ζ2ζy(29)
where in the second equations of (28) and (29), the
transformed form of the continuity equation (27)
is used. In terms of Qand Rin place of Uand V
respectively, equations (23) and (24) become
ζ2Qt+Q(2ζQx−Qζx+ζRy) + R(ζQy−Qζy)
+2
5ζ2(Q2ζxxx +R2ζxyy) + 4
5ζ2(QQxζxx +RRxζyy)
+7
22 ζ4[ζ(ζxxx +ζxyy)+3ζx(ζxx +ζyy)]
+ζ3ζx+m2QpQ2+R2
ζ1/3= 0 (30)
and
ζ2Rt+ζ(RQx+QRx+ 2RRy)−R(Qζx+Rζy)
+2
5ζ2(Q2ζxxy +R2ζyyy) + 4
5ζ2(QQyζxx +RRyζyy)
+7
22 ζ4[ζ(ζxxy +ζyyy)+3ζy(ζxx +ζyy )]+
ζ3ζy+m2RpQ2+R2
ζ1/3= 0 (31)
For numerically solving the partial differential
equations (29) - (31) an invertibility assumption is
now made as in Bose [14]. Since the solution of
these equations is of the form ζ=ζ(x, y, t), Q =
Q(x, y, t), R =R(x, y, t), if the inverse of
these functions exist then, x=x(ζ, Q, R), y =
y(ζ, Q, R), t =t(ζ, Q, R). A sufficient condi-
tion for the inversion to be possible is that the
Jacobian of the transformation does no vanish at
points within the domain. Under this assump-
tion, Q=Q[x(ζ, Q, R), y(ζ, Q, R), t(ζ, Q, R)]
(a function of ζ, Q, R) and similarly, R=
R[x(ζ, Q, R), y(ζ, Q, R), t(ζ, Q, R)] (also a func-
tion of ζ, Q, R). The fixed point solution of this
pair iterative pair of equations, for each ζ, by elim-
ination of Rand Qrespectively, is of the form
Q=F(ζ), R =G(ζ) (32)
provided the functios Fand Gare Lipschitz con-
tinuous. Physically as was observed in Bose [14],
equation (32) implies that the discharge compo-
nents Qand Rat a point (x, y) and at time tde-
pends solely on the local elevation ζat that point
and time, as in the theory of long tidal waves. The
numerical method developed here may therefore be
considered as an extension of that theory to nonlin-
ear turbulent waves. Substituting the forms (32) in
equations (27) yields the linear conservation equa-
tion
ζt+Fx(ζ) + Gy(ζ) = 0 (33)
The same substitution in Eqs. (30), (31) however
yields a pair of coupled second degree equations in
the derivatives F0(ζ) and G0(ζ) that are difficult to
solve algebraically.
In the important case of almost linear crests
prpagating in the direction of the x-axis however,
|ζy|<< 1 and |G|<< |F|so that the two equations
can be solved approximately. To the first order of
ζyand G, one gets
a F 02−b F F 0+c= 0 (34)
and
G0=1
ζx(ζF 0−F)hζxF0G−ζxF G
ζ+2
5ζζxxy
+4
5ζζyζxx F F 0+7
22 ζ3{ζ(ζxxy +ζyyy)+3ζxxζy}
+ζ2ζy+m2F G sgn(F)
ζ4/3((35)
where
a=ζ2, b = 2ζ1 + 2
5ζζxx,
c=1−2
5ζ2ζxxx
ζx−m2sgn(F)
ζxζ1/3F2
−7
22 ζ4hζ
ζx
(ζxxx +ζxyy) + 3 (ζxx +ζyy)i−ζ3(36)
Equations (33) - (35) are now treated numerically.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
39
Volume 17, 2022
7 The Numerical Procedure
Since the discharge Qor F(ζ) is an increas-
ing function of ζ, the appropriate solution of the
quadratic equation (34) is
F0=bF +√b2F2−4ac
2a(37)
For real solution b2F2≥4ac representing undu-
lating flow surface; otherwise F0= 0 for steady
discharge. Eq. (37) is inegrated by the trape-
zoidal rule based iterative, second order, modified
Euler method (Bose [15], p. 264). The coeffi-
cients a, b, c given in equation (36) contain par-
tial derivatives of ζup to the third order partial
derivatives. These derivatives are calculated by
second order finite differences of ζover a grid of
points x= (m−1) ∆x, y = (n−1) ∆y, (m=
1,2,3,··· ;n= 1,2,3,···), where ∆x, ∆yare
the grid lengths parallel to the xand y- axes re-
spectively. The value of Fand Gat the grid points
are then obtained by the modified Euler method by
recursion as:
F(m+1, n) = F(m, n)+∆x×[F0(m, n)ζx(m, n)]
←1
2[F(m, n) + F(m+ 1, n)] (38)
G(m, n+ 1) = G(m, n) + ∆y×[G0(m, n)ζy(m, n)]
←1
2[G(m, n) + G(m, n + 1)] (39)
Equations (38), (39) determine the values of the
functions Fand Gfor a given elevation ζat a
given time t. The evolution of ζas solution of equa-
tion (33) is obtained by application of the second
order Lax-Richtmyer scheme for increasing time
t=l∆t, (l= 1,2,3,···). According to Richt-
myer and Morton [16], the time evolution of ζfrom
given step lis obtained in two steps l+ 1 and l+ 2
:
ζl+1
m,n =1
4ζl
m+1,n +ζl
m−1,n +ζl
m,n+1 +ζl
m,n−1
−∆t
2 ∆xFl
m+1,n−Fl
m−1,n−∆t
2 ∆yGl
m,n+1−Gl
m,n−1
(40)
where Fl
m,n =F(ζl
m,n), and at the second step
ζl+2
m,n =ζl
m,n −∆t
∆xFl+1
m+1,n −Fl+1
m−1,n
−∆t
∆yGl+1
m,n+1 −Gl+1
m,n−1(41)
Equation (41) determines the elevation at the next
step. In practice one can take ∆x= ∆yand
∆t=r∆x, where the factor r < 1 is so chosen
that the computation remains stable and conver-
gent for small changes in its value.
8 Application to a model Initial Surface El-
evation
The numerical method of the preceding section
is applied to a model initial surface elevation con-
sidered by Howe [17], p.332 for tsunami propaga-
tion caused by sudden uplifting of the sea bed in
the form of a ridge caused by subduction of the
bed in front. Accordingly the initial elevation of
the sea surface is assumed to follow the dislocated
bed form with a profile given by the equation
ζ= 1 + ζ01
{1 + [(x+L)2+y2]/L2}3/2
−1
{1 + [(x−L)2+y2]/L2}3/2(42)
where ζ0= maximum elevation/depression
above/below the mean surface of water, and 2L
= distance between the highest and the lowest
ponts of the initial free surface. For presentable re-
sults it is assumed that ζ0= 0.3 and L= 10. The
profile tapers to plane level ζ= 1 as |y|tends to
infinity. The initial value of ζat t= 0 is computed
for x, y ≥0 from Eq. (42), taking ∆x= ∆y= 0.1
and ∆t=r∆x, where r= 1/2. The numerical
scheme is found to yield stable and convergent val-
ues of ζfor small changes in the value of r. The
initial value of the functions Fand Gare estimated
from the following consideration. Since the crests
are nearly parallel to the y-axis, Gis negligible. If
moreover, the crests are assumed to move with a
constant sub-critical velocity c0to avoid choking
of flow, then following Eq. (33), F0(ζ) = c0, so
that F(ζ) = c0ζ. Accoridngly a value of c0= 0.2
is assumed here to compute F(ζ) at the different
grid points. The spatial derivatives of the function
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
40
Volume 17, 2022
ζare then computed by the second order finite
difference formulae:
ζx=ζ(x+ ∆x, y)−ζ(x−∆x, y)
(∆x)2(43)
ζxx =ζ(x+ ∆x, y)−2ζ(x, y) + ζ(x−∆x, y)
(∆x)2
(44)
ζxxx(x, y) = 1
∆x3[ζ(x+ 2∆x, y)−3ζ(x+ ∆x, y)
+3 ζ(x, y)−ζ(x−∆x, y)] (45)
ζxxy(x, y) = 1
2(∆x)3[ζ(x+ ∆x, y + ∆y)
−ζ(x+∆x, y −∆y)−2ζ(x, y +∆y)+2ζ(x, y −∆y)
+ζ(x−∆x, y + ∆y)−ζ(x−∆x, y −∆y)] (46)
and similar expressions for ζy(x, y),ζyy (x, y),
ζyyy(x, y), and ζxyy(x, y), for x, y > 0 or, m, n ≥2.
For x= 0 and y= 0, linear extrapolations
ζx(0, y) = 2 ζx(∆x, y)−ζx(2∆x, y) (47)
ζx(x, 0) = 2 ζx(x, ∆y)−ζx(x, 2∆y) (48)
and similar expressions for the other partial deriva-
tives are used. The time evolution of the elevation
ζfollowing the Lax-Richtmyer scheme (40), (41)
is then carried out using Eqs. (38), (39). The
corresponding value of the discharge functions F
and Gis carried out as at the initial time, i.e.
F=c0ζand G= 0. Performing a large number
of iterations over time, the computed results for
ζat times t= 0,20,40 and 60 with respect to x
for the three representative cases y= 0, y=1
2x,
and y=xare presented in the figures 1, 2, and
3 respectively. The three cases correspond to the
crests being viewed from the angles 0o, 26.57oand
45orespectively to the crests. It is found in the
computed data that the crest height diminishes
slowly with time with corresponding rise in the
trough depth. This feature diminishes with the
increase in the angle, becoming insignificant in the
third case.
-20 0 20 40 60 80
x
0
0.5
1
1.5
Figure 1. View of propagating surface from 0oangle, at times t= 0,20,40,60...
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
41
Volume 17, 2022
-20 0 20 40 60
x
0
0.5
1
1.5
Figure 2. View of propagating surfae from 26.57oangle, at times t= 0,20,40...
-20 0 20 40
x
0
0.5
1
1.5
-20 0 20 40
Figure 3. View of propagating surface from 45oangle at times t= 0,20,...
The method can be employed to generate data for
any other cross section.
As a comparison of the computed results pre-
sented above, one may compare these with those
of Howe [17], presented on page 333. The surface
profiles presented there are based on potential sur-
face wave theory according to which the elevation
ζis given by the Eq. (5.7.7), p. 330 of the text.
The notable qualitative difference between the two
approaches is that while in the theory presented
here, while the wave form realistically leans for-
ward as it propagates, the form does not replicate
this characteristic in the theory presented in Howe
[17].
9. Results and Discussion
In this paper, the turbulent shallow water free
surface flow equations in two dimensions are de-
veloped and numerically integrated by a method
developed in Bose [14]. Following Bose and Dey
[12], the basic equations of continuity and momen-
tum conservation for the flow are developed from
the Raynolds Averaged Navier Stokes equations.
The equations are expressed in terms of the depth
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
42
Volume 17, 2022
averaged horizontal velocity (U, V ) and the surface
elevation ζ, all being functions of x, y and time t.
Based on certain assumptions necessitated by the
presence of turbulence, three nonlinear differential
equations are obtained for the determination of
ζ, U and V. For numerical treatment of the equa-
tions, a slight simplification is made in the two hor-
izontal momentum equations where in the terms
arising from instantaneous vertical acceleration,
the quantities ∂U/∂t, ∂V /∂t and ∂U/∂x +∂V/∂y
are replaced by equivalent expressions in terms of
ζobtained by linearizing the two equations . A
transformation of the equations is then made by
replacing Uand Vby discharge components Q
and Rrespectively. An invertibility argument of
the functions shows that Qand Rcan be func-
tions of ζ, say Q=F(ζ) and R=G(ζ), as in
the theory of long tidal waves. The method may
therefore be considered an extension of the lin-
ear inviscid long wave theory to nonlinear waves
taking in to account the generated turbulence in
real fluids. Substitution in the continuity equation
leads to a conservation equation of standard form,
while the two momentum equations yield coupled
second degree equations in the derivatives F0(ζ)
and G0(ζ). In the important case of propagation
of nearly linear crested waves on the water surface
parallel to y-axis, F0(ζ) is given to a first order
by a decoupled quadratic equation, and G0(ζ) by
an expression containing F(ζ) and F0(ζ). The
two equations are numerically treated to compute
F(ζ) and G(ζ), following which the standard form
mass conservation equation is numerically solved
by the Lax-Richtmyer scheme. The method is im-
plemented for an almost linear initial bed elevation
model given in Howe [17] for simulating tsunami
propagation. The computed values are shown for
three different angles of view, showing very slowly
decreasing wave height. The method developed in
the paper opens up the possibility of generating
data for creating animation of the moving free sur-
face for greater visual effect.
Acknowledgement
The author is thankful to the SN Bose National
Center for Basic Sciences, Kolkata for providing
necessary facilities for undertaking this research.
References
[1] C.L. Mader, Numerical Modeling of Water
Waves, Boca Raton: CRC Press, 2004.
[2] R.J. Fennema and M.H Chaudhry, “ Implicit
methods for two-dimensional unsteady free surface
flows.” J. Hydraul. Res., vol. 27, pp. 321-332,
1989.
[3] R.J. Fennema and M.H. Chaudhry, “ Explicit
methods for 2-D transient free surface flows.” J.
Hydraul. Engg., vol. 116, pp. 610-614, 1990.
[4] P. Wilders, Th. L.van Stijn, G.S. Stelling, and
G.A. Fokkema, “ A fully implicit splitting method
for accurate tidal computation.” Int. J. Num.
Meth. Engg., vol. 26, pp. 2707-2721, 1988.
[5] V. Casulli, “ Semi-implicit finite difference
methods for two-dimensional shallow water equa-
tions.” J. Comput. Phys., vol. 86, pp. 56-74, 1990.
[6] C. Ouyang, S. He, Q. Xu, Y. Luo and W.
Zhang, W., “ A MacCormack-TVD finite difference
method to simulate the mass flow in mountain-
ous terrain with variable computational domain.”
Computers Geosci., vol. 52, pp. 1-10, 2013.
[7] S. Fu and B.R. Hodges, Numerical solution of
the shallow water equations using a time-centered
split-implicit method.” 18th Engg. Mech. Div.
Confce. (EMD2007), pp. 1-6, 2007.
[8] N.D. Katopodes and T. Strelkoff, “ Two-
dimensional shallow water-wave models.” ASCE
J. Engg. Mech. Div., vol. 105, pp. 317-334, 1979.
[9] F. Alcrudo and P. Garcia-Navarro, “ A high-
resolution Godunove-type scheme in finite volumes
for the 2D shallow-water equations.” Int. J. Num.
Meth. in Fluids, vol. 16, pp. 489-505, 1993.
[10] S-H. Peng and C. Tang, “ Development and
application of two-dimensional numerical model on
shallow water flows using finite volume method,”
J. App. Math. Phys., vol. 3, pp. 989-996, 2015.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
43
Volume 17, 2022
[11] P. Garcia-Navarro, M.E Hubbard and A.
Priestley, “ Genuinely multidimensional upwinding
for the 2D shallow water equations.” J. Comput.
Phys., vol. 121, pp. 79-93, 1995.
[12] S.K. Bose and S. Dey, “ Curvilinear flow pro-
files based on Reynolds averaging. ” J. Hydraul.
Engg., vol. 133, pp. 1074-1079, 2007.
[13] H. Sclichting, Boundary Layer Theory, Berlim:
Springer, 1979.
[14] S.K. Bose, “ A numerical method for the so-
lution of the nonlinear turbulent one-dimensional
free surface flow equations.” Comput. Geosci., vol.
22, pp. 81-86, 2018.
[15] S.K. Bose, Numerical Methods of Mathemat-
ics Implemented in Fortran, Singapore: Springer,
2019.
[16] R.D. Richtmyer and K.W. Morton, Difference
Methods for Initial Value Problems, New York: In-
terscience, 1967.
[17] M.S. Howe,, Hydrodynamics and Sound, Cam-
bridge: Camb. Univ. Press, 2007.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the Creative
Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2022.17.5
Sujit K. Bose
E-ISSN: 2224-3461
44
Volume 17, 2022
