# Numerical Simulation for 2D Sedimentation Model by an Upwind Discontinuous Galerkin Procedure

Authors: JinFeng Jian; HuanZhen Chen; BaoHai Shi
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We reformulate the mathematical model for the 2D sedimentation in an estuary as a coupled nonlinear differential system. Combining the mass-conservation character of the discontinuous Galerkin method and the jump-capturing property ofLesaint-Raviart upwind technique, we design an upwind discontinuous Galerkin finite element method, which obeys the local mass conservation and possesses good stability. Our theoretical analysis shows that there exists a unique solution to the numerical procedure and the discrete solution permits O
(
h
k
+
Δ
t
) convergence rate. Numerical experiments are conducted to verify our theoretical findings. This may provide a theoretical principle for better understanding of the mechanism and morphological characters of sedimentation at estuaries.


---

## Body

## 1. Introduction

This article is initiated from our project “Numerical Simulation for the Sedimentation at Yellow River’s Estuary.” It is known that the estuary of Yellow River is a weak-tidal, arenose, frequently swinging, and accumulated firth. Its sedimentation is heavily influenced by the reduced amount of water but a relatively large amount of sand from the upstream and the oceanic dynamics such as the tide, current, and wave in the littoral region. Numerical modeling and simulation for the flow and sedimentation are very important means of quantitatively predicting the formation and development of the delta and the evolution of the fluvial bed, as well as the silt deposition transport, and then providing a theoretical principle for better understanding of the mechanism and morphological characters of the sedimentation. This will be of great significance for the further development of Yellow River delta.Based upon the literature of the sedimentation studies [1–5] for estuaries, theoretically we reformulate the mathematical model for 2D sedimentation as a combination of 4 governing equations: the flow continuity equation, the flow motion equation, the sediment continuity equation, and the bed deformation equation, in which the first two equations describe the water dynamics of the littoral region and the last two are for the dynamics of sediment deposition.From the point of view of numerical and applicable issue, there are two important ingredients which should be strengthened in designing an ideal numerical procedure for this model. The first one is that the procedure should obey the local mass/momentum conservation law to reflect the physical characteristics of the dynamics of the flow and sediment, and the second is that the procedure should recognize the flow direction so that the procedure can capture the jumps or interfaces and maintain the good stability since the velocity of the fluid is much larger than that of the sediment deposition. For this purpose, we recall the mass-conservation character of the discontinuous Galerkin method and the jump-capturing property of upwind technique and attempt to use these merits so as to design a strongly stable and locally conservative numerical procedure for the mathematical model.The goals of this article are as follows: (1) in the flow continuity equation, we reformulate the mathematical model by taking the bed height Z
b varying with time so that the flow continuity equation is naturally coupled with the sediment continuity equation and the bed deformation equation. By doing so, the mathematical model can much better describe the practical sediment process compared to the simplified case of Z
b-constant [6–8] and (2) by introducing the Broken Sobolev spaces and applying the discontinuous Galerkin method to discrete those space-derivative terms to design a nonsymmetric interior penalty Galerkin (NIPG) method. Considering the fact that the flow velocity U is much larger than that of the sediment deposition, we pay our special attention to those terms involving U by insertingLesaint-Raviart upwind technique [9]. Therefore, an upwind discontinuous Galerkin finite element method is proposed for the reformulated mathematical model. (3) With the help of the induction argument, we prove that the upwind discontinuous Galerkin finite element method possesses good stability, from which the existence and uniqueness of the discrete solution follow. (4) Combining the proof that the discrete depth H
n of the fluid possesses a positive lower bound, we prove that the approximate solutions for the depth of fluid H, the velocity U, the sediment concentration S, and the bed height Z
b are of O
(
h
k
+
Δ
t
) convergence rates, respectively. (5) Numerical experiments are conducted to verify the theoretical findings, stability, local conservation property, and convergence rates.The achievement of these five goals constitutes the main content of the article, presented in Sections2–6, respectively. For simplicity of presentation, we will use W
m
,
p
(
E
) to denote the usual Sobolev spaces that provided the norm ·
m
,
p
,
E and the seminorm ·
m
,
p
,
E for any 2D domain E. When p
=
2 and E
=
Ω, the subscript p and Ω will be omitted. We also let (
·
,
·
) and ·
,
· stand for the inner products on E and on the edges ∂
E, respectively, and · denotes their L
2-norms. C is a generic constant and ε is a small positive constant, which may take different values at different places but independent of the triangulation parameters h and Δ
t. We define L
p
(
a
,
b
;
X
)
=
{
f
:
∫
a
b
f
(
·
,
t
)
X
p
d
t
<
∞
}, where if p
=
∞, the integral is replaced with the essential supreme.


## 2. Reformulation of the Mathematical Model

LetΩ be a bounded polygonal domain in R
2, and let [
0
,
T
] be a time interval. Then, the sediment transportation model at an estuary can be described by the following 4 governing equations defined over the space-time domain Ω
×
0
,
T; here we ignore the Coriolis force, since it does not affect the correctness of the model.The  flow  continuity  equation is(1)
∂
H
∂
t
+
∇
·
U
H
=
0
,the  flow  motion  equation is(2)
∂
U
∂
t
+
U
·
∇
U
+
g
∇
Z
+
τ
U
-
μ
Δ
U
=
0
,the  sediment  continuity  equation is(3)
∂
S
∂
t
+
U
·
∇
S
-
c
0
Δ
S
+
a
ω
H
S
-
S
∗
=
0
,and the  bed  deformation  equation is(4)
∂
Z
b
∂
t
-
a
ω
γ
′
S
-
S
∗
=
0
,where the unknowns H
,
Z
,
S
,
U
=
(
u
,
v
)
T and Z
b
=
Z
-
H denote the depth of water, the water level, the sediment concentration, the vertical mean velocity, and the bed height, respectively; S
∗ is the sediment carrying capacity, g is the gravitational acceleration, τ
=
g
(
u
2
+
v
2
/
c
2
H
), c
=
(
1
/
n
)
H
1
/
6 denotes the Chezy coefficient, n denotes Manning’s roughness coefficient, μ is the viscosity coefficient, c
0 is the sediment diffusion coefficient, a is the coefficient of saturation recovery, ω is the settling velocity of sediment, and γ
′ is the sediment dry bulk density.In the uncoupled sedimentation model [5] and the shallow water system (SWS) [10, 11], (1) is replaced by ∂
Z
/
∂
t
+
∇
·
(
U
Z
)
=
0. It means that the flow continuity equation is not coupled with the bed deformation equation; in other words, the bed height Z
b was assumed to be independent of time in the flow continuity equation. Therefore, the shallow water equations (1) and (2) and the sediment deposition models (3) and (4) can be solved separately. We remark that this assumption will lead to a simple calculation, but it is unreasonable in the real process of sediment transport. In this paper, we assume Z
b changes with time in shallow water equations, so that the flow continuity equation is much naturally coupled with the sediment continuity equation and the bed deformation equation.Without loss of generality, the boundary of the domain is decomposed into three parts: the inflow partΓ
1, the rigid boundary Γ
2, and the outflow part Γ
3 (see Figure 1), and the initial-boundary conditions are prescribed as(5)
H
x
,
y
,
0
=
H
0
x
,
y
,
U
x
,
y
,
0
=
U
0
x
,
y
,
S
x
,
y
,
0
=
S
0
x
,
y
,
Z
b
x
,
y
,
0
=
Z
b
0
x
,
y
,
(6)
U
=
U
^
,
U
^
·
n
<
0
,
S
=
S
^
,
H
=
H
^
,
on
Γ
1
,
U
=
0
,
∂
S
∂
n
=
0
,
on
Γ
2
,
U
=
U
^
,
U
^
·
n
>
0
,
∂
S
∂
n
=
0
,
on
Γ
3
,where the initial values H
0
,
U
0
,
S
0
,
Z
b
0 and the boundary conditions U
^
,
H
^
,
S
^ are given data, andn denotes the unit outward normal vector to ∂
Ω.Figure 1
The illustration ofΩ.


## 3. The Upwind Discontinuous Galerkin Finite Element Procedure

In this section, we shall revisit the Broken Sobolev spaces and modify those terms involvingU byLesaint-Raviart upwind technique and then apply the discontinuous Galerkin method to discrete those space-derivative terms in (1)–(4) to design an upwind nonsymmetric interior penalty Galerkin (NIPG) method.


### 3.1. Revisiting the Broken Sobolev Spaces

We begin this subsection by recalling some notations of the triangulation for the given domainΩ. We let ε
h be a quasi-uniform subdivision of Ω with element E; let Ω
h
=
⋃
E
∈
ε
h
E, h
E denote the diameter of E, and h
=
max
⁡
{
h
E
}. Further we need to describe the edges of an element and the functions defined on them. We denote, by e, any edge of an element E
∈
ε
h, and then e belongs to either its interior or boundary of Ω. We also denote, by Γ
h, the set of interior edges. If e
∈
∂
Ω, the unit normal vector n
e of E is taken to be outward to ∂
Ω, or if e is a common edge of two neighbor elements E
e
+ and E
e
-, we understand that n
e is oriented from E
e
- to E
e
+. For any function v of H
1
(
ε
h
), we denote by v
± its traces on E
e
+ and E
e
- and define its average and jump by {
v
}
=
(
v
+
+
v
-
)
/
2 and [
v
]
=
v
-
-
v
+, respectively. When e
∈
∂
Ω, we naturally define {
v
}
=
[
v
]
=
v
-
. With the help of these notions, we present the definitions of the Broken Sobolev spaces and the functionals on the whole edges. The Broken Sobolev spaces H
s
(
ε
h
) for any real number s
≥
1 are defined by (7)
H
s
ε
h
=
v
∈
L
2
Ω
:
v
∣
E
∈
H
s
E
,
∀
E
∈
ε
h
,equipped with the norm (8)
v
H
s
ε
h
=
∑
E
∈
ε
h
v
H
s
E
2
1
/
2
. Here v
H
s
E
2
=
∑
α
≤
s
D
α
v
0
,
E
2
1
/
2 with D
α
v are being the α
t
h-order generalized derivatives of v. The functionals with respect to the averages and jumps on the whole edges are defined by (9)
J
1
:
H
1
ε
h
2
×
H
1
ε
h
2
⟶
R
,
J
1
V
,
W
=
∑
e
∈
Γ
h
∪
∂
Ω
σ
e
1
e
V
,
W
e
,
J
2
:
H
1
ε
h
×
H
1
ε
h
⟶
R
,
J
2
v
,
w
=
∑
e
∈
Γ
h
∪
Γ
1
σ
e
2
e
v
,
w
e
,where e is the length of e and the selected nonnegative real numbers σ
e
1
,
σ
e
2 are called penalty parameters.


### 3.2. Variational Formulation for the Mathematical Model

Letφ be any function in H
1
(
ε
h
)
. We multiply (1) by suitable test functions φ on each element E, apply the Green formula, and then sum the resulting equations over the domain to obtain the variational form for the flow continuity equation(10)
H
t
,
φ
-
∑
E
∈
ε
h
H
U
,
∇
φ
E
+
∑
e
∈
Γ
h
H
U
·
n
e
,
φ
e
+
∑
e
∈
Γ
1
H
^
U
^
·
n
e
,
φ
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
U
^
·
n
e
,
φ
e
=
0
,
∀
φ
∈
H
1
ε
h
.Similarly, we obtain the variational forms for the rest of the three equations (2)–(4): (11)
U
t
,
Φ
+
U
·
∇
U
,
Φ
-
∑
E
g
H
+
Z
b
,
∇
Φ
+
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
H
,
Φ
·
n
e
+
∑
e
∈
Γ
h
∪
∂
Ω
g
Z
b
,
Φ
·
n
e
+
τ
U
,
Φ
+
∑
E
μ
∇
U
,
∇
Φ
E
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
U
·
n
e
,
Φ
e
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
Φ
·
n
e
,
U
e
+
J
1
U
,
Φ
=
∑
e
∈
∂
Ω
μ
∇
Φ
·
n
e
,
U
^
e
+
∑
e
∈
∂
Ω
σ
e
1
e
U
^
,
Φ
e
-
∑
e
∈
Γ
1
g
H
^
,
Φ
·
n
e
,
∀
Φ
∈
H
1
ε
h
2
,
S
t
,
ψ
-
∑
E
S
∇
·
U
,
ψ
E
-
∑
E
S
U
,
∇
ψ
E
+
∑
e
∈
Γ
h
∪
∂
Ω
S
U
·
n
e
,
ψ
e
+
∑
E
c
0
∇
S
,
∇
ψ
E
-
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
S
·
n
e
,
ψ
e
+
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
ψ
·
n
e
,
S
e
+
a
ω
H
S
-
S
∗
,
ψ
+
J
2
S
,
ψ
=
∑
e
∈
Γ
1
c
0
∇
ψ
·
n
e
,
S
^
e
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
,
ψ
e
,
∀
ψ
∈
H
1
ε
h
,
Z
b
t
,
ϕ
=
a
ω
γ
′
S
-
S
∗
,
ϕ
,
∀
ϕ
∈
H
1
ε
h
. Consequently, the variational formulation corresponding to models (1)–(4) is defined as to find the map H
,
U
,
S
,
Z
b
:
0
,
T
↦
(
H
1
ε
h
)
2
×
(
H
1
(
ε
h
)
)
2
×
H
1
(
ε
h
)
×
H
1
(
ε
h
), satisfying (10)-(11).


### 3.3. The Upwind Discontinuous Galerkin Finite Element Scheme

Assume thatk is a positive integer and P
k
(
E
) denotes the space of polynomials of degree less than or equal to k on E. Then, the finite element spaces used later are chosen to be the subspace of H
s
(
ε
h
) defined by (12)
D
k
ε
h
=
v
∈
L
2
Ω
:
v
∣
E
∈
P
k
E
,
∀
E
∈
ε
h
. For time discretization we divide the time interval [0, T] by the nodes t
n
(
n
=
0,1
,
…
,
N
) with time step Δ
t
=
T
/
N and let u
n
(
x
,
y
) be the node-function value of u
(
t
,
x
,
y
).Considering the dominance of convection in a real sedimentation process, we hope that the numerical procedure should recognize the direction of the flow, so as to design a stable algorithm. For this purpose, we introduce the upwind values forH and S on the interior edge e by H
up and S
up, respectively. First we define the upwind values of H and S on the interior edge e by H
up and S
up, respectively as follows: (13)
H
up
=
H
-
if
U
·
n
e
>
0
,
H
+
if
U
·
n
e
≤
0
,
S
up
=
S
-
if
U
·
n
e
>
0
,
S
+
if
U
·
n
e
≤
0
.In discretization, we use backward Euler’s scheme to approximate the time derivatives and use the nonsymmetric interior penalty Galerkin method to the diffusion part, in which the average values ofH
,
S on interior edge are replaced by H
up
,
S
up and the nonlinear convective term U
·
∇
U is handled by theLesaint-Raviart upwind technique. By doing so, the desired stability can be guaranteed. For easy computation we use, for example, U
n
H
n
+
1 instead of U
n
+
1
H
n
+
1 to linearize the nonlinear terms.Upon these ideas, we propose our fully upwind discontinuous Galerkin finite element procedure (UWDG) as follows: find the map{
H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
}
:
{
t
n
}
n
=
0
N
↦
D
k
(
ε
h
)
×
(
D
k
(
ε
h
)
)
2
×
D
k
(
ε
h
)
×
D
k
(
ε
h
) satisfying (14)
H
h
n
+
1
-
H
h
n
Δ
t
,
φ
-
∑
E
∈
ε
h
H
h
n
+
1
U
h
n
,
∇
φ
E
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
φ
e
+
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
,
φ
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
h
n
+
1
U
h
n
·
n
e
,
φ
e
=
0
,
∀
φ
∈
D
k
ε
h
,
(15)
U
h
n
+
1
-
U
h
n
Δ
t
,
Φ
+
∑
E
U
h
n
·
∇
U
h
n
+
1
,
Φ
E
+
∑
E
U
h
n
·
n
e
U
h
n
+
1
,
i
n
t
-
U
h
n
+
1
,
e
x
t
,
Φ
i
n
t
∂
E
-
-
∑
E
g
H
h
n
+
1
+
Z
b
h
n
+
1
,
∇
Φ
+
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
H
h
n
+
1
,
Φ
·
n
e
e
+
∑
e
∈
Γ
h
∪
∂
Ω
g
Z
b
h
n
+
1
,
Φ
·
n
e
e
+
τ
h
n
U
h
n
+
1
,
Φ
+
∑
E
μ
∇
U
h
n
+
1
,
∇
Φ
E
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
U
h
n
+
1
·
n
e
,
Φ
e
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
Φ
·
n
e
,
U
h
n
+
1
e
+
J
1
U
h
n
+
1
,
Φ
=
∑
e
∈
∂
Ω
μ
∇
Φ
·
n
e
,
U
^
n
+
1
e
+
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
,
Φ
e
-
∑
e
∈
Γ
1
g
H
^
n
+
1
,
Φ
·
n
e
e
,
∀
Φ
∈
D
k
ε
h
2
,
(16)
S
h
n
+
1
-
S
h
n
Δ
t
,
ω
-
∑
E
S
h
n
+
1
∇
·
U
h
n
,
ω
E
-
∑
E
S
h
n
+
1
U
h
n
,
∇
ω
E
+
∑
e
∈
Γ
h
S
h
n
+
1
,
up
U
h
n
·
n
e
,
ω
e
+
∑
e
∈
Γ
1
S
^
n
+
1
U
h
n
·
n
e
,
ω
e
+
∑
e
∈
Γ
2
∪
Γ
3
S
h
n
+
1
U
h
n
·
n
e
,
ω
e
+
∑
E
c
0
∇
S
h
n
+
1
,
∇
ω
E
-
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
S
h
n
+
1
·
n
e
,
ω
e
+
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
ω
·
n
e
,
S
h
n
+
1
e
+
a
ω
H
h
n
S
h
n
+
1
-
S
∗
,
ω
+
J
2
S
h
n
+
1
,
ω
=
∑
e
∈
Γ
1
c
0
∇
ω
·
n
e
,
S
^
n
+
1
e
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
,
ω
e
,
∀
ω
∈
D
k
ε
h
,
(17)
Z
b
h
n
+
1
-
Z
b
h
n
Δ
t
,
ϕ
=
a
ω
γ
′
S
h
n
+
1
-
S
∗
,
ϕ
,
∀
ϕ
∈
D
k
ε
h
.The initial values H
h
0
,
U
h
0
,
S
h
0
,
and
Z
b
h
0 also are clarified by their L
2-orthogonal projections: (18)
H
0
-
H
h
0
,
φ
=
0
,
∀
φ
∈
D
k
ε
h
,
U
0
-
U
h
0
,
Φ
=
0
,
∀
Φ
∈
D
k
ε
h
2
,
S
0
-
S
h
0
,
ω
=
0
,
∀
ω
∈
D
k
ε
h
,
Z
b
0
-
Z
b
h
0
,
ϕ
=
0
,
∀
ϕ
∈
D
k
ε
h
.Here τ
h
n
=
g
(
U
h
n
l
2
/
c
2
H
h
n
)
,
∂
E
-
=
{
x
∈
∂
E
:
{
U
h
n
}
·
n
e
<
0
}, v
i
n
t, and v
e
x
t denote the interior and exterior traces of v on ∂
E
-. If ∂
E
-
∈
∂
Ω, then v
e
x
t
=
v
^
,
v
i
n
t
=
v.


## 3.1. Revisiting the Broken Sobolev Spaces

We begin this subsection by recalling some notations of the triangulation for the given domainΩ. We let ε
h be a quasi-uniform subdivision of Ω with element E; let Ω
h
=
⋃
E
∈
ε
h
E, h
E denote the diameter of E, and h
=
max
⁡
{
h
E
}. Further we need to describe the edges of an element and the functions defined on them. We denote, by e, any edge of an element E
∈
ε
h, and then e belongs to either its interior or boundary of Ω. We also denote, by Γ
h, the set of interior edges. If e
∈
∂
Ω, the unit normal vector n
e of E is taken to be outward to ∂
Ω, or if e is a common edge of two neighbor elements E
e
+ and E
e
-, we understand that n
e is oriented from E
e
- to E
e
+. For any function v of H
1
(
ε
h
), we denote by v
± its traces on E
e
+ and E
e
- and define its average and jump by {
v
}
=
(
v
+
+
v
-
)
/
2 and [
v
]
=
v
-
-
v
+, respectively. When e
∈
∂
Ω, we naturally define {
v
}
=
[
v
]
=
v
-
. With the help of these notions, we present the definitions of the Broken Sobolev spaces and the functionals on the whole edges. The Broken Sobolev spaces H
s
(
ε
h
) for any real number s
≥
1 are defined by (7)
H
s
ε
h
=
v
∈
L
2
Ω
:
v
∣
E
∈
H
s
E
,
∀
E
∈
ε
h
,equipped with the norm (8)
v
H
s
ε
h
=
∑
E
∈
ε
h
v
H
s
E
2
1
/
2
. Here v
H
s
E
2
=
∑
α
≤
s
D
α
v
0
,
E
2
1
/
2 with D
α
v are being the α
t
h-order generalized derivatives of v. The functionals with respect to the averages and jumps on the whole edges are defined by (9)
J
1
:
H
1
ε
h
2
×
H
1
ε
h
2
⟶
R
,
J
1
V
,
W
=
∑
e
∈
Γ
h
∪
∂
Ω
σ
e
1
e
V
,
W
e
,
J
2
:
H
1
ε
h
×
H
1
ε
h
⟶
R
,
J
2
v
,
w
=
∑
e
∈
Γ
h
∪
Γ
1
σ
e
2
e
v
,
w
e
,where e is the length of e and the selected nonnegative real numbers σ
e
1
,
σ
e
2 are called penalty parameters.


## 3.2. Variational Formulation for the Mathematical Model

Letφ be any function in H
1
(
ε
h
)
. We multiply (1) by suitable test functions φ on each element E, apply the Green formula, and then sum the resulting equations over the domain to obtain the variational form for the flow continuity equation(10)
H
t
,
φ
-
∑
E
∈
ε
h
H
U
,
∇
φ
E
+
∑
e
∈
Γ
h
H
U
·
n
e
,
φ
e
+
∑
e
∈
Γ
1
H
^
U
^
·
n
e
,
φ
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
U
^
·
n
e
,
φ
e
=
0
,
∀
φ
∈
H
1
ε
h
.Similarly, we obtain the variational forms for the rest of the three equations (2)–(4): (11)
U
t
,
Φ
+
U
·
∇
U
,
Φ
-
∑
E
g
H
+
Z
b
,
∇
Φ
+
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
H
,
Φ
·
n
e
+
∑
e
∈
Γ
h
∪
∂
Ω
g
Z
b
,
Φ
·
n
e
+
τ
U
,
Φ
+
∑
E
μ
∇
U
,
∇
Φ
E
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
U
·
n
e
,
Φ
e
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
Φ
·
n
e
,
U
e
+
J
1
U
,
Φ
=
∑
e
∈
∂
Ω
μ
∇
Φ
·
n
e
,
U
^
e
+
∑
e
∈
∂
Ω
σ
e
1
e
U
^
,
Φ
e
-
∑
e
∈
Γ
1
g
H
^
,
Φ
·
n
e
,
∀
Φ
∈
H
1
ε
h
2
,
S
t
,
ψ
-
∑
E
S
∇
·
U
,
ψ
E
-
∑
E
S
U
,
∇
ψ
E
+
∑
e
∈
Γ
h
∪
∂
Ω
S
U
·
n
e
,
ψ
e
+
∑
E
c
0
∇
S
,
∇
ψ
E
-
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
S
·
n
e
,
ψ
e
+
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
ψ
·
n
e
,
S
e
+
a
ω
H
S
-
S
∗
,
ψ
+
J
2
S
,
ψ
=
∑
e
∈
Γ
1
c
0
∇
ψ
·
n
e
,
S
^
e
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
,
ψ
e
,
∀
ψ
∈
H
1
ε
h
,
Z
b
t
,
ϕ
=
a
ω
γ
′
S
-
S
∗
,
ϕ
,
∀
ϕ
∈
H
1
ε
h
. Consequently, the variational formulation corresponding to models (1)–(4) is defined as to find the map H
,
U
,
S
,
Z
b
:
0
,
T
↦
(
H
1
ε
h
)
2
×
(
H
1
(
ε
h
)
)
2
×
H
1
(
ε
h
)
×
H
1
(
ε
h
), satisfying (10)-(11).


## 3.3. The Upwind Discontinuous Galerkin Finite Element Scheme

Assume thatk is a positive integer and P
k
(
E
) denotes the space of polynomials of degree less than or equal to k on E. Then, the finite element spaces used later are chosen to be the subspace of H
s
(
ε
h
) defined by (12)
D
k
ε
h
=
v
∈
L
2
Ω
:
v
∣
E
∈
P
k
E
,
∀
E
∈
ε
h
. For time discretization we divide the time interval [0, T] by the nodes t
n
(
n
=
0,1
,
…
,
N
) with time step Δ
t
=
T
/
N and let u
n
(
x
,
y
) be the node-function value of u
(
t
,
x
,
y
).Considering the dominance of convection in a real sedimentation process, we hope that the numerical procedure should recognize the direction of the flow, so as to design a stable algorithm. For this purpose, we introduce the upwind values forH and S on the interior edge e by H
up and S
up, respectively. First we define the upwind values of H and S on the interior edge e by H
up and S
up, respectively as follows: (13)
H
up
=
H
-
if
U
·
n
e
>
0
,
H
+
if
U
·
n
e
≤
0
,
S
up
=
S
-
if
U
·
n
e
>
0
,
S
+
if
U
·
n
e
≤
0
.In discretization, we use backward Euler’s scheme to approximate the time derivatives and use the nonsymmetric interior penalty Galerkin method to the diffusion part, in which the average values ofH
,
S on interior edge are replaced by H
up
,
S
up and the nonlinear convective term U
·
∇
U is handled by theLesaint-Raviart upwind technique. By doing so, the desired stability can be guaranteed. For easy computation we use, for example, U
n
H
n
+
1 instead of U
n
+
1
H
n
+
1 to linearize the nonlinear terms.Upon these ideas, we propose our fully upwind discontinuous Galerkin finite element procedure (UWDG) as follows: find the map{
H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
}
:
{
t
n
}
n
=
0
N
↦
D
k
(
ε
h
)
×
(
D
k
(
ε
h
)
)
2
×
D
k
(
ε
h
)
×
D
k
(
ε
h
) satisfying (14)
H
h
n
+
1
-
H
h
n
Δ
t
,
φ
-
∑
E
∈
ε
h
H
h
n
+
1
U
h
n
,
∇
φ
E
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
φ
e
+
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
,
φ
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
h
n
+
1
U
h
n
·
n
e
,
φ
e
=
0
,
∀
φ
∈
D
k
ε
h
,
(15)
U
h
n
+
1
-
U
h
n
Δ
t
,
Φ
+
∑
E
U
h
n
·
∇
U
h
n
+
1
,
Φ
E
+
∑
E
U
h
n
·
n
e
U
h
n
+
1
,
i
n
t
-
U
h
n
+
1
,
e
x
t
,
Φ
i
n
t
∂
E
-
-
∑
E
g
H
h
n
+
1
+
Z
b
h
n
+
1
,
∇
Φ
+
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
H
h
n
+
1
,
Φ
·
n
e
e
+
∑
e
∈
Γ
h
∪
∂
Ω
g
Z
b
h
n
+
1
,
Φ
·
n
e
e
+
τ
h
n
U
h
n
+
1
,
Φ
+
∑
E
μ
∇
U
h
n
+
1
,
∇
Φ
E
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
U
h
n
+
1
·
n
e
,
Φ
e
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
Φ
·
n
e
,
U
h
n
+
1
e
+
J
1
U
h
n
+
1
,
Φ
=
∑
e
∈
∂
Ω
μ
∇
Φ
·
n
e
,
U
^
n
+
1
e
+
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
,
Φ
e
-
∑
e
∈
Γ
1
g
H
^
n
+
1
,
Φ
·
n
e
e
,
∀
Φ
∈
D
k
ε
h
2
,
(16)
S
h
n
+
1
-
S
h
n
Δ
t
,
ω
-
∑
E
S
h
n
+
1
∇
·
U
h
n
,
ω
E
-
∑
E
S
h
n
+
1
U
h
n
,
∇
ω
E
+
∑
e
∈
Γ
h
S
h
n
+
1
,
up
U
h
n
·
n
e
,
ω
e
+
∑
e
∈
Γ
1
S
^
n
+
1
U
h
n
·
n
e
,
ω
e
+
∑
e
∈
Γ
2
∪
Γ
3
S
h
n
+
1
U
h
n
·
n
e
,
ω
e
+
∑
E
c
0
∇
S
h
n
+
1
,
∇
ω
E
-
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
S
h
n
+
1
·
n
e
,
ω
e
+
∑
e
∈
Γ
h
∪
Γ
1
c
0
∇
ω
·
n
e
,
S
h
n
+
1
e
+
a
ω
H
h
n
S
h
n
+
1
-
S
∗
,
ω
+
J
2
S
h
n
+
1
,
ω
=
∑
e
∈
Γ
1
c
0
∇
ω
·
n
e
,
S
^
n
+
1
e
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
,
ω
e
,
∀
ω
∈
D
k
ε
h
,
(17)
Z
b
h
n
+
1
-
Z
b
h
n
Δ
t
,
ϕ
=
a
ω
γ
′
S
h
n
+
1
-
S
∗
,
ϕ
,
∀
ϕ
∈
D
k
ε
h
.The initial values H
h
0
,
U
h
0
,
S
h
0
,
and
Z
b
h
0 also are clarified by their L
2-orthogonal projections: (18)
H
0
-
H
h
0
,
φ
=
0
,
∀
φ
∈
D
k
ε
h
,
U
0
-
U
h
0
,
Φ
=
0
,
∀
Φ
∈
D
k
ε
h
2
,
S
0
-
S
h
0
,
ω
=
0
,
∀
ω
∈
D
k
ε
h
,
Z
b
0
-
Z
b
h
0
,
ϕ
=
0
,
∀
ϕ
∈
D
k
ε
h
.Here τ
h
n
=
g
(
U
h
n
l
2
/
c
2
H
h
n
)
,
∂
E
-
=
{
x
∈
∂
E
:
{
U
h
n
}
·
n
e
<
0
}, v
i
n
t, and v
e
x
t denote the interior and exterior traces of v on ∂
E
-. If ∂
E
-
∈
∂
Ω, then v
e
x
t
=
v
^
,
v
i
n
t
=
v.


## 4. Stability Analysis

In this section, we shall analyze the stability of the upwind DG scheme and prove the existence and uniqueness of the discrete solution by borrowing some lemmas from the literature on finite elements and making some reasonable hypotheses.


### 4.1. Some Lemmas and Hypotheses

To conduct convergence analysis in the next section we need to define theL
2-orthogonal projections for the unknowns {
H
,
U
,
S
,
Z
b
}
. These projections are maps {
H
~
,
U
~
,
S
~
,
Z
b
~
}
:
[
0
,
T
]
↦
D
k
(
ε
h
)
×
(
D
k
(
ε
h
)
)
2
×
D
k
(
ε
h
)
×
D
k
(
ε
h
) satisfying (19)
H
-
H
~
,
φ
=
0
,
∀
φ
∈
D
k
ε
h
,
U
-
U
~
,
Φ
=
0
,
∀
Φ
∈
D
k
ε
h
2
,
S
-
S
~
,
ω
=
0
,
∀
ω
∈
D
k
ε
h
,
Z
b
-
Z
b
~
,
ϕ
=
0
,
∀
ϕ
∈
D
k
ε
h
.Following the practical background of a real sedimentation process, we shall make some reasonable hypotheses on the exact solutions of the model.(H1)
There exists a unique solution{
H
,
U
,
S
,
Z
b
} to the models (1)–(4). Then, one can deduce, at least, the existence of the variational formulations (10)-(11).
(H2)
There is a positive constantH
∗ such that H
∗
≤
H
(
x
,
t
).
(H3)
H
0
(
x
)
∈
H
k
+
1
(
Ω
), U
0
(
x
)
∈
(
H
k
+
1
(
Ω
)
)
2, S
0
(
x
)
∈
H
k
+
1
(
Ω
), and Z
b
0
(
x
)
∈
H
k
+
1
(
Ω
).
(H4)
H
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
∩
W
∞
1
(
Ω
)
), U
(
x
,
t
)
∈
L
2
(
0
,
T
;
(
H
k
+
1
(
Ω
)
)
2
∩
(
W
∞
1
(
Ω
)
)
2
), S
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
∩
W
∞
1
(
Ω
)
), and Z
b
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
).
(H5)
∂
2
H
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), ∂
2
U
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), ∂
2
S
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), and ∂
2
Z
b
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
).
(H6)
∂
U
/
∂
t
∈
L
2
(
0
,
T
;
(
H
1
(
Ω
)
)
2
∩
(
L
∞
(
Ω
)
)
2
) and ∂
H
/
∂
t
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
).Next we introduce the following lemmas which are crucial to our convergence analysis.Lemma 1 (see [12–14]).
Assume that (H4)–(H6) hold andH
~
,
U
~
,
S
~ and Z
b
~ are defined in (19). Then, the following error estimates hold for 0
≤
p
≤
∞ and 0
≤
l
≤
k
+
1:(20)
H
-
H
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
H
k
+
1
,
p
,
Ω
,
U
-
U
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
U
k
+
1
,
p
,
Ω
,
S
-
S
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
S
k
+
1
,
p
,
Ω
,
Z
b
-
Z
b
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
Z
b
k
+
1
,
p
,
Ω
.Lemma 2 (trace theorem [12, 15]).
Assume thatv
∈
H
s
(
E
) and E
∈
ε
h
. Then, there exists a constant C
>
0 such that (21)
v
0
,
e
≤
C
e
1
/
2
E
-
1
/
2
v
0
,
E
+
h
∇
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
1
,
∇
v
·
n
e
0
,
e
≤
C
e
1
/
2
E
-
1
/
2
∇
v
0
,
E
+
h
∇
2
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
2
. In particular, if v
∈
P
k
(
E
), then (22)
v
0
,
e
≤
C
h
-
1
/
2
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
1
,
∇
v
·
n
e
0
,
e
≤
C
h
-
1
/
2
∇
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
2
.Lemma 3 (discrete Gronwall inequality [15, 16]).
Assume that the nonnegative sequencesa
n
,
b
n and c
n satisfy (23)
a
n
+
Δ
t
∑
i
=
0
n
b
i
≤
B
+
C
Δ
t
∑
i
=
0
n
a
i
+
Δ
t
∑
i
=
0
n
c
i
,
∀
n
≥
0
. Then, we have(24)
a
n
+
Δ
t
∑
i
=
0
n
b
i
≤
e
C
n
+
1
Δ
t
B
+
Δ
t
∑
i
=
0
n
c
i
.Lemma 4 (inverse inequalities [12, 16]).
Letχ
∈
D
k
(
ε
h
)
. Then, there is a constant K
0
>
0 such that(25)
χ
L
∞
≤
K
0
h
-
1
χ
,
∀
χ
∈
D
k
ε
h
,
χ
H
1
≤
K
0
h
-
1
χ
,
∀
χ
∈
D
k
ε
h
.Lemma 5 (boundedness).
Assume thath is sufficiently small. Then, there is a constant C
>
0 such that, for ∀
n
≥
0, (26)
(a)
H
~
n
L
∞
Ω
≥
H
∗
2
,
(b)
H
~
n
L
∞
Ω
+
U
~
n
L
∞
Ω
2
+
S
~
n
L
∞
Ω
≤
C
.Proof.
Combining the definition ofH
~, the hypotheses (H2) and (H4), and Lemma 1, we obtain (27)
H
~
n
L
∞
Ω
≥
H
n
L
∞
Ω
-
H
~
n
-
H
n
L
∞
Ω
≥
H
∗
-
K
h
H
n
1
,
∞
,
Ω
. Let h
0
=
min
⁡
{
1
,
H
∗
/
2
K
H
n
1
,
∞
,
Ω
}; then if h
<
h
0, we have (28)
H
~
n
L
∞
Ω
≥
H
∗
2
.On the other hand, by the triangle inequality, we can deduce (29)
H
~
n
L
∞
Ω
≤
H
n
L
∞
Ω
+
H
~
n
-
H
n
L
∞
Ω
≤
1
+
K
h
H
n
1
,
∞
,
Ω
,
U
~
n
L
∞
Ω
≤
U
n
L
∞
Ω
+
U
~
n
-
U
n
L
∞
Ω
≤
1
+
K
h
U
n
1
,
∞
,
Ω
,
S
~
n
L
∞
Ω
≤
S
n
L
∞
Ω
+
S
~
n
-
S
n
L
∞
Ω
≤
1
+
K
h
S
n
1
,
∞
,
Ω
. Let C
=
1
+
K
h
0
max
⁡
{
H
n
1
,
∞
,
Ω
,
U
n
1
,
∞
,
Ω
,
S
n
1
,
∞
,
Ω
}; then if h
<
h
0, we have (30)
H
~
n
L
∞
Ω
+
U
~
n
L
∞
Ω
2
+
S
~
n
L
∞
Ω
≤
C
.


### 4.2. Stability Estimate

Preceding the proof of the stability for (14)–(18), we first present a bounded estimate for H
h
,
U
h
,
S
h, and Z
b
h. Since its proof involves an induction argument with the error bounds of these unknowns we here only give its presentation in Lemma 6, whose detailed proof can be found at the last part of the next section.For the concise exposition we decompose the errors in the following manner: let(31)
H
h
-
H
=
e
H
+
θ
H
≕
H
h
-
H
~
+
H
~
-
H
,
U
h
-
U
=
e
U
+
θ
U
≕
U
h
-
U
~
+
U
~
-
U
,
S
h
-
S
=
e
S
+
θ
S
≕
S
h
-
S
~
+
S
~
-
S
,
Z
b
h
-
Z
b
=
e
Z
b
+
θ
Z
b
≕
Z
b
h
-
Z
b
~
+
Z
b
~
-
Z
b
.Lemma 6 (boundedness of approximate solution).
There exists a small positive constanth
0 such that when 0
<
h
≤
h
0 and Δ
t
=
o
(
h
), for ∀
n
≥
0, (32)
(a)
H
h
n
L
∞
Ω
≥
H
∗
4
,
(b)
e
H
n
L
∞
Ω
+
e
U
n
L
∞
Ω
2
+
e
S
n
L
∞
Ω
≤
H
∗
4
.Theorem 7 (stability estimate).
There exists a constantC independent of h and Δ
t such that for all m
(
0
<
m
≤
N
), (33)
H
h
m
+
U
h
m
+
S
h
m
+
Z
b
h
m
≤
C
H
h
0
+
U
h
0
+
S
h
0
+
Z
b
h
0
+
Δ
t
1
/
2
∑
E
∈
ε
h
∇
U
h
0
+
Δ
t
1
/
2
J
1
U
h
0
,
U
h
0
+
C
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
0
,
e
2
+
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
0
,
e
2
+
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
Γ
1
H
^
n
+
1
0
,
e
2
1
/
2
.Proof.
Settingφ
=
H
h
n
+
1 in (14), using Green’s formula for the second term on the left-hand side, we have (34)
-
∑
E
∈
ε
h
H
h
n
+
1
U
h
n
,
∇
H
h
n
+
1
E
=
-
1
2
∑
E
∈
ε
h
U
h
n
,
∇
H
h
n
+
1
2
E
=
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
H
h
n
+
1,2
E
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
∂
Ω
U
h
n
·
n
e
,
H
h
n
+
1,2
e
.Using [
a
b
]
=
[
a
]
{
b
}
+
{
a
}
[
b
], we see that (35)
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
U
h
n
·
n
e
,
1
e
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
H
h
n
+
1
e
=
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
H
h
n
+
1
e
=
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
H
h
n
+
1
2
e
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
.Combining (34) and (35), we rewrite (14) as the following form: (36)
H
h
n
+
1
-
H
h
n
Δ
t
,
H
h
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
H
h
n
+
1
2
e
=
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
H
h
n
+
1,2
E
+
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
+
1
2
∑
e
∈
Γ
1
U
h
n
·
n
e
,
H
h
n
+
1,2
e
-
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
,
H
h
n
+
1
e
-
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
·
n
e
,
H
h
n
+
1,2
=
∑
i
=
1
5
I
i
. Next we estimate the values of the five terms I
i. By Lemma 6 and Young’s inequality, we obtain (37)
I
1
≤
ε
∑
E
∈
ε
h
∇
U
h
n
2
+
C
H
h
n
+
1
2
.Using the trace theorem, we have (38)
I
2
≤
C
∑
e
∈
Γ
h
σ
e
1
e
1
/
2
U
h
n
L
2
e
·
σ
e
1
e
-
1
/
2
H
h
n
+
1
L
2
e
≤
ε
J
1
U
h
n
,
U
h
n
+
C
H
h
n
+
1
2
,where J
1
(
U
h
n
,
U
h
n
) is the bilinear form defined in Section 3.1.
Similarly,(39)
I
3
+
I
4
+
I
5
≤
ε
J
1
U
h
n
,
U
h
n
+
C
H
h
n
+
1
2
. Combining (
1
/
2
)
(
a
2
-
b
2
)
≤
(
a
-
b
)
a and (37)–(39), (36) becomes (40)
1
2
Δ
t
H
h
n
+
1
2
-
H
h
n
2
≤
ε
∑
E
∈
ε
h
∇
U
h
n
2
+
C
H
h
n
+
1
2
+
ε
J
1
U
h
n
,
U
h
n
.Similarly, setting Φ
=
U
h
n
+
1 in (15), ω
=
S
h
n
+
1 in (16), and ϕ
=
Z
b
h
n
+
1 in (17) and using the same techniques, we have (41)
1
2
Δ
t
U
h
n
+
1
2
-
U
h
n
2
+
τ
1
/
2
U
h
n
+
1
2
+
∑
E
μ
1
/
2
∇
U
h
n
+
1
2
+
J
1
U
h
n
+
1
,
U
h
n
+
1
≤
C
∑
e
∈
Γ
1
H
^
n
+
1
L
2
e
2
+
H
h
n
+
1
2
+
Z
b
h
n
+
1
2
+
U
h
n
+
1
2
+
U
h
n
2
+
ε
∑
E
∈
ε
h
∇
U
h
n
+
1
2
+
C
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
L
2
e
,
1
2
Δ
t
S
h
n
+
1
2
-
S
h
n
2
+
∑
E
c
0
1
/
2
∇
S
h
n
+
1
2
+
J
2
e
S
n
+
1
,
e
S
n
+
1
≤
ε
J
1
U
h
n
,
U
h
n
+
∑
E
∈
ε
h
∇
U
h
n
2
+
∑
E
∈
ε
h
∇
S
h
n
+
1
2
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
L
2
e
+
C
S
h
n
+
1
2
+
U
h
n
2
+
H
h
n
2
,
1
2
Δ
t
Z
b
h
n
+
1
2
-
Z
b
h
n
2
≤
C
S
h
n
+
1
2
+
ε
Z
b
h
n
2
. Multiplying (40)-(41) by 2
Δ
t, summing both sides with respect to n over n
=
0,1
,
2
,
…
,
m
-
1, and using the discrete Gronwall inequality (Lemma 3), we complete the proof of Theorem 7.From Theorem7, we see that if the initial-boundary conditions (5)-(6) are zero, then H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
(
∀
n
>
0
) are all equal to zero. This implies that the solution of the UWDG scheme (14)–(18) is unique and thus the existence.Theorem 8 (existence and uniqueness).
For∀
n
>
0, there exists a unique solution {
H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
} to the UWDG scheme (14)–(18).


## 4.1. Some Lemmas and Hypotheses

To conduct convergence analysis in the next section we need to define theL
2-orthogonal projections for the unknowns {
H
,
U
,
S
,
Z
b
}
. These projections are maps {
H
~
,
U
~
,
S
~
,
Z
b
~
}
:
[
0
,
T
]
↦
D
k
(
ε
h
)
×
(
D
k
(
ε
h
)
)
2
×
D
k
(
ε
h
)
×
D
k
(
ε
h
) satisfying (19)
H
-
H
~
,
φ
=
0
,
∀
φ
∈
D
k
ε
h
,
U
-
U
~
,
Φ
=
0
,
∀
Φ
∈
D
k
ε
h
2
,
S
-
S
~
,
ω
=
0
,
∀
ω
∈
D
k
ε
h
,
Z
b
-
Z
b
~
,
ϕ
=
0
,
∀
ϕ
∈
D
k
ε
h
.Following the practical background of a real sedimentation process, we shall make some reasonable hypotheses on the exact solutions of the model.(H1)
There exists a unique solution{
H
,
U
,
S
,
Z
b
} to the models (1)–(4). Then, one can deduce, at least, the existence of the variational formulations (10)-(11).
(H2)
There is a positive constantH
∗ such that H
∗
≤
H
(
x
,
t
).
(H3)
H
0
(
x
)
∈
H
k
+
1
(
Ω
), U
0
(
x
)
∈
(
H
k
+
1
(
Ω
)
)
2, S
0
(
x
)
∈
H
k
+
1
(
Ω
), and Z
b
0
(
x
)
∈
H
k
+
1
(
Ω
).
(H4)
H
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
∩
W
∞
1
(
Ω
)
), U
(
x
,
t
)
∈
L
2
(
0
,
T
;
(
H
k
+
1
(
Ω
)
)
2
∩
(
W
∞
1
(
Ω
)
)
2
), S
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
∩
W
∞
1
(
Ω
)
), and Z
b
(
x
,
t
)
∈
L
2
(
0
,
T
;
H
k
+
1
(
Ω
)
).
(H5)
∂
2
H
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), ∂
2
U
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), ∂
2
S
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
), and ∂
2
Z
b
/
∂
t
2
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
).
(H6)
∂
U
/
∂
t
∈
L
2
(
0
,
T
;
(
H
1
(
Ω
)
)
2
∩
(
L
∞
(
Ω
)
)
2
) and ∂
H
/
∂
t
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
).Next we introduce the following lemmas which are crucial to our convergence analysis.Lemma 1 (see [12–14]).
Assume that (H4)–(H6) hold andH
~
,
U
~
,
S
~ and Z
b
~ are defined in (19). Then, the following error estimates hold for 0
≤
p
≤
∞ and 0
≤
l
≤
k
+
1:(20)
H
-
H
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
H
k
+
1
,
p
,
Ω
,
U
-
U
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
U
k
+
1
,
p
,
Ω
,
S
-
S
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
S
k
+
1
,
p
,
Ω
,
Z
b
-
Z
b
~
l
,
p
,
Ω
≤
C
h
k
+
1
-
l
Z
b
k
+
1
,
p
,
Ω
.Lemma 2 (trace theorem [12, 15]).
Assume thatv
∈
H
s
(
E
) and E
∈
ε
h
. Then, there exists a constant C
>
0 such that (21)
v
0
,
e
≤
C
e
1
/
2
E
-
1
/
2
v
0
,
E
+
h
∇
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
1
,
∇
v
·
n
e
0
,
e
≤
C
e
1
/
2
E
-
1
/
2
∇
v
0
,
E
+
h
∇
2
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
2
. In particular, if v
∈
P
k
(
E
), then (22)
v
0
,
e
≤
C
h
-
1
/
2
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
1
,
∇
v
·
n
e
0
,
e
≤
C
h
-
1
/
2
∇
v
0
,
E
,
∀
e
∈
∂
E
,
s
≥
2
.Lemma 3 (discrete Gronwall inequality [15, 16]).
Assume that the nonnegative sequencesa
n
,
b
n and c
n satisfy (23)
a
n
+
Δ
t
∑
i
=
0
n
b
i
≤
B
+
C
Δ
t
∑
i
=
0
n
a
i
+
Δ
t
∑
i
=
0
n
c
i
,
∀
n
≥
0
. Then, we have(24)
a
n
+
Δ
t
∑
i
=
0
n
b
i
≤
e
C
n
+
1
Δ
t
B
+
Δ
t
∑
i
=
0
n
c
i
.Lemma 4 (inverse inequalities [12, 16]).
Letχ
∈
D
k
(
ε
h
)
. Then, there is a constant K
0
>
0 such that(25)
χ
L
∞
≤
K
0
h
-
1
χ
,
∀
χ
∈
D
k
ε
h
,
χ
H
1
≤
K
0
h
-
1
χ
,
∀
χ
∈
D
k
ε
h
.Lemma 5 (boundedness).
Assume thath is sufficiently small. Then, there is a constant C
>
0 such that, for ∀
n
≥
0, (26)
(a)
H
~
n
L
∞
Ω
≥
H
∗
2
,
(b)
H
~
n
L
∞
Ω
+
U
~
n
L
∞
Ω
2
+
S
~
n
L
∞
Ω
≤
C
.Proof.
Combining the definition ofH
~, the hypotheses (H2) and (H4), and Lemma 1, we obtain (27)
H
~
n
L
∞
Ω
≥
H
n
L
∞
Ω
-
H
~
n
-
H
n
L
∞
Ω
≥
H
∗
-
K
h
H
n
1
,
∞
,
Ω
. Let h
0
=
min
⁡
{
1
,
H
∗
/
2
K
H
n
1
,
∞
,
Ω
}; then if h
<
h
0, we have (28)
H
~
n
L
∞
Ω
≥
H
∗
2
.On the other hand, by the triangle inequality, we can deduce (29)
H
~
n
L
∞
Ω
≤
H
n
L
∞
Ω
+
H
~
n
-
H
n
L
∞
Ω
≤
1
+
K
h
H
n
1
,
∞
,
Ω
,
U
~
n
L
∞
Ω
≤
U
n
L
∞
Ω
+
U
~
n
-
U
n
L
∞
Ω
≤
1
+
K
h
U
n
1
,
∞
,
Ω
,
S
~
n
L
∞
Ω
≤
S
n
L
∞
Ω
+
S
~
n
-
S
n
L
∞
Ω
≤
1
+
K
h
S
n
1
,
∞
,
Ω
. Let C
=
1
+
K
h
0
max
⁡
{
H
n
1
,
∞
,
Ω
,
U
n
1
,
∞
,
Ω
,
S
n
1
,
∞
,
Ω
}; then if h
<
h
0, we have (30)
H
~
n
L
∞
Ω
+
U
~
n
L
∞
Ω
2
+
S
~
n
L
∞
Ω
≤
C
.


## 4.2. Stability Estimate

Preceding the proof of the stability for (14)–(18), we first present a bounded estimate for H
h
,
U
h
,
S
h, and Z
b
h. Since its proof involves an induction argument with the error bounds of these unknowns we here only give its presentation in Lemma 6, whose detailed proof can be found at the last part of the next section.For the concise exposition we decompose the errors in the following manner: let(31)
H
h
-
H
=
e
H
+
θ
H
≕
H
h
-
H
~
+
H
~
-
H
,
U
h
-
U
=
e
U
+
θ
U
≕
U
h
-
U
~
+
U
~
-
U
,
S
h
-
S
=
e
S
+
θ
S
≕
S
h
-
S
~
+
S
~
-
S
,
Z
b
h
-
Z
b
=
e
Z
b
+
θ
Z
b
≕
Z
b
h
-
Z
b
~
+
Z
b
~
-
Z
b
.Lemma 6 (boundedness of approximate solution).
There exists a small positive constanth
0 such that when 0
<
h
≤
h
0 and Δ
t
=
o
(
h
), for ∀
n
≥
0, (32)
(a)
H
h
n
L
∞
Ω
≥
H
∗
4
,
(b)
e
H
n
L
∞
Ω
+
e
U
n
L
∞
Ω
2
+
e
S
n
L
∞
Ω
≤
H
∗
4
.Theorem 7 (stability estimate).
There exists a constantC independent of h and Δ
t such that for all m
(
0
<
m
≤
N
), (33)
H
h
m
+
U
h
m
+
S
h
m
+
Z
b
h
m
≤
C
H
h
0
+
U
h
0
+
S
h
0
+
Z
b
h
0
+
Δ
t
1
/
2
∑
E
∈
ε
h
∇
U
h
0
+
Δ
t
1
/
2
J
1
U
h
0
,
U
h
0
+
C
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
0
,
e
2
+
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
0
,
e
2
+
Δ
t
∑
n
=
1
m
-
1
∑
e
∈
Γ
1
H
^
n
+
1
0
,
e
2
1
/
2
.Proof.
Settingφ
=
H
h
n
+
1 in (14), using Green’s formula for the second term on the left-hand side, we have (34)
-
∑
E
∈
ε
h
H
h
n
+
1
U
h
n
,
∇
H
h
n
+
1
E
=
-
1
2
∑
E
∈
ε
h
U
h
n
,
∇
H
h
n
+
1
2
E
=
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
H
h
n
+
1,2
E
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
∂
Ω
U
h
n
·
n
e
,
H
h
n
+
1,2
e
.Using [
a
b
]
=
[
a
]
{
b
}
+
{
a
}
[
b
], we see that (35)
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
U
h
n
·
n
e
,
1
e
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
H
h
n
+
1
e
=
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
+
∑
e
∈
Γ
h
H
h
n
+
1
,
up
U
h
n
·
n
e
,
H
h
n
+
1
e
=
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
H
h
n
+
1
2
e
-
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
.Combining (34) and (35), we rewrite (14) as the following form: (36)
H
h
n
+
1
-
H
h
n
Δ
t
,
H
h
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
H
h
n
+
1
2
e
=
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
H
h
n
+
1,2
E
+
1
2
∑
e
∈
Γ
h
H
h
n
+
1,2
,
U
h
n
·
n
e
e
+
1
2
∑
e
∈
Γ
1
U
h
n
·
n
e
,
H
h
n
+
1,2
e
-
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
,
H
h
n
+
1
e
-
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
·
n
e
,
H
h
n
+
1,2
=
∑
i
=
1
5
I
i
. Next we estimate the values of the five terms I
i. By Lemma 6 and Young’s inequality, we obtain (37)
I
1
≤
ε
∑
E
∈
ε
h
∇
U
h
n
2
+
C
H
h
n
+
1
2
.Using the trace theorem, we have (38)
I
2
≤
C
∑
e
∈
Γ
h
σ
e
1
e
1
/
2
U
h
n
L
2
e
·
σ
e
1
e
-
1
/
2
H
h
n
+
1
L
2
e
≤
ε
J
1
U
h
n
,
U
h
n
+
C
H
h
n
+
1
2
,where J
1
(
U
h
n
,
U
h
n
) is the bilinear form defined in Section 3.1.
Similarly,(39)
I
3
+
I
4
+
I
5
≤
ε
J
1
U
h
n
,
U
h
n
+
C
H
h
n
+
1
2
. Combining (
1
/
2
)
(
a
2
-
b
2
)
≤
(
a
-
b
)
a and (37)–(39), (36) becomes (40)
1
2
Δ
t
H
h
n
+
1
2
-
H
h
n
2
≤
ε
∑
E
∈
ε
h
∇
U
h
n
2
+
C
H
h
n
+
1
2
+
ε
J
1
U
h
n
,
U
h
n
.Similarly, setting Φ
=
U
h
n
+
1 in (15), ω
=
S
h
n
+
1 in (16), and ϕ
=
Z
b
h
n
+
1 in (17) and using the same techniques, we have (41)
1
2
Δ
t
U
h
n
+
1
2
-
U
h
n
2
+
τ
1
/
2
U
h
n
+
1
2
+
∑
E
μ
1
/
2
∇
U
h
n
+
1
2
+
J
1
U
h
n
+
1
,
U
h
n
+
1
≤
C
∑
e
∈
Γ
1
H
^
n
+
1
L
2
e
2
+
H
h
n
+
1
2
+
Z
b
h
n
+
1
2
+
U
h
n
+
1
2
+
U
h
n
2
+
ε
∑
E
∈
ε
h
∇
U
h
n
+
1
2
+
C
∑
e
∈
∂
Ω
σ
e
1
e
U
^
n
+
1
L
2
e
,
1
2
Δ
t
S
h
n
+
1
2
-
S
h
n
2
+
∑
E
c
0
1
/
2
∇
S
h
n
+
1
2
+
J
2
e
S
n
+
1
,
e
S
n
+
1
≤
ε
J
1
U
h
n
,
U
h
n
+
∑
E
∈
ε
h
∇
U
h
n
2
+
∑
E
∈
ε
h
∇
S
h
n
+
1
2
+
∑
e
∈
Γ
1
σ
e
2
e
S
^
n
+
1
L
2
e
+
C
S
h
n
+
1
2
+
U
h
n
2
+
H
h
n
2
,
1
2
Δ
t
Z
b
h
n
+
1
2
-
Z
b
h
n
2
≤
C
S
h
n
+
1
2
+
ε
Z
b
h
n
2
. Multiplying (40)-(41) by 2
Δ
t, summing both sides with respect to n over n
=
0,1
,
2
,
…
,
m
-
1, and using the discrete Gronwall inequality (Lemma 3), we complete the proof of Theorem 7.From Theorem7, we see that if the initial-boundary conditions (5)-(6) are zero, then H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
(
∀
n
>
0
) are all equal to zero. This implies that the solution of the UWDG scheme (14)–(18) is unique and thus the existence.Theorem 8 (existence and uniqueness).
For∀
n
>
0, there exists a unique solution {
H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n
} to the UWDG scheme (14)–(18).


## 5. Convergence Analysis

In this section, we shall conduct error estimates for the UWDG solution and prove theO
(
h
k
+
Δ
t
) convergence rate based upon the prior-estimate techniques and induction argument. Combining these error estimates and applying the induction argument again, we give the detailed proof for the boundedness of the approximate solution stated in Lemma 6.


### 5.1. Error Equations

First, we derive the error equation of the flow continuous equation (1). Choosing the test function φ
=
e
H
n
+
1 in (10) and (14), subtracting (10) from (14), and using the definition of θ
H in (31), we derive the error equation for the flow continuous equation as follows: (42)
e
H
n
+
1
-
e
H
n
Δ
t
,
e
H
n
+
1
-
∑
E
∈
ε
h
e
H
n
+
1
U
h
n
,
∇
e
H
n
+
1
E
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
-
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
+
∑
E
∈
ε
h
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
∇
e
H
n
+
1
-
∑
e
∈
Γ
h
H
~
n
+
1
,
up
U
h
n
·
n
e
-
H
n
+
1
U
n
+
1
·
n
e
,
e
H
n
+
1
e
-
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
-
H
^
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
-
∑
e
∈
Γ
2
∪
Γ
3
H
h
n
+
1
U
h
n
·
n
e
-
H
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
. Taking into account that there is no diffusion term in the flow continuous equation (1), we apply Green’s formula to the second term on the left-hand side (43)
-
∑
E
∈
ε
h
e
H
n
+
1
U
h
n
,
∇
e
H
n
+
1
E
=
-
1
2
∑
E
∈
ε
h
U
h
n
,
∇
e
H
n
+
1
2
E
=
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
∂
Ω
U
h
n
·
n
e
,
e
H
n
+
1,2
e
. Using [
a
b
]
=
[
a
]
{
b
}
+
{
a
}
[
b
], we obtain (44)
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
U
h
n
·
n
e
,
1
e
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
. Similarly, we use Green’s formula for the second term on the right-hand side; thus (45)
∑
E
∈
ε
h
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
∇
e
H
n
+
1
=
-
∑
E
∈
ε
h
∇
·
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
e
H
n
+
1
+
∑
e
∈
Γ
1
H
~
n
+
1
U
h
n
-
H
^
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
+
∑
e
∈
Γ
h
e
H
n
+
1
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
·
n
e
,
1
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
~
n
+
1
U
h
n
-
H
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
. Combining (43)–(45), we rewrite the error equation (42) into the following form: (46)
e
H
n
+
1
-
e
H
n
Δ
t
,
e
H
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
=
-
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
+
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
+
1
2
∑
e
∈
Γ
1
U
h
n
-
U
^
n
·
n
e
,
e
H
n
+
1,2
e
-
∑
E
∈
ε
h
∇
·
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
e
H
n
+
1
+
∑
e
∈
Γ
h
e
H
n
+
1
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
·
n
e
,
1
e
-
∑
e
∈
Γ
h
H
~
n
+
1
,
up
U
h
n
·
n
e
-
H
n
+
1
U
n
+
1
·
n
e
,
e
H
n
+
1
e
+
∑
e
∈
Γ
1
H
~
n
+
1
-
H
^
n
+
1
U
h
n
·
n
e
,
e
H
n
+
1
e
-
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
-
U
^
n
·
n
e
,
e
H
n
+
1,2
.Similarly, we obtain the other three error equations (47)
e
U
n
+
1
-
e
U
n
Δ
t
,
e
U
n
+
1
+
τ
n
+
1
e
U
n
+
1
,
e
U
n
+
1
+
∑
E
μ
∇
e
U
n
+
1
,
∇
e
U
n
+
1
E
+
J
1
e
U
n
+
1
,
e
U
n
+
1
=
-
U
n
+
1
-
U
n
Δ
t
-
U
t
n
+
1
,
e
U
n
+
1
+
∑
E
U
n
+
1
·
∇
U
n
+
1
-
U
h
n
·
∇
U
h
n
+
1
,
e
U
n
+
1
E
-
∑
E
U
h
n
·
n
e
U
h
n
+
1
,
i
n
t
-
U
h
n
+
1
,
e
x
t
,
e
U
n
+
1
,
i
n
t
∂
E
-
-
τ
n
+
1
θ
U
n
+
1
,
e
U
n
+
1
+
τ
n
+
1
-
τ
h
n
U
h
n
+
1
,
e
U
n
+
1
-
∑
E
μ
∇
θ
U
n
+
1
,
∇
e
U
n
+
1
E
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
θ
U
n
+
1
·
n
e
,
e
U
n
+
1
e
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
e
U
n
+
1
·
n
e
,
θ
U
n
+
1
e
-
J
1
θ
U
n
+
1
,
e
U
n
+
1
+
∑
E
e
H
n
+
1
,
∇
e
U
n
+
1
+
∑
E
θ
H
n
+
1
,
∇
e
U
n
+
1
+
∑
E
e
Z
b
n
+
1
,
∇
e
U
n
+
1
+
∑
E
θ
Z
b
n
+
1
,
∇
e
U
n
+
1
-
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
e
H
n
+
1
+
θ
H
n
+
1
,
e
U
n
+
1
·
n
e
e
-
∑
e
∈
Γ
h
∪
∂
Ω
g
e
Z
b
n
+
1
+
θ
Z
b
n
+
1
,
e
U
n
+
1
·
n
e
e
,
e
S
n
+
1
-
e
S
n
Δ
t
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
∑
E
c
0
∇
e
S
n
+
1
,
∇
e
S
n
+
1
E
+
J
2
e
S
n
+
1
,
e
S
n
+
1
=
-
S
n
+
1
-
S
n
Δ
t
-
S
t
n
+
1
,
e
S
n
+
1
+
∑
E
S
h
n
+
1
∇
·
U
h
n
-
S
n
+
1
∇
·
U
n
+
1
,
e
S
n
+
1
E
-
1
2
∑
E
∇
·
U
h
n
,
e
S
n
+
1,2
+
1
2
∑
e
∈
Γ
h
e
S
n
+
1,2
,
U
h
n
·
n
e
e
-
∑
e
∈
Γ
1
S
^
n
+
1
U
h
n
·
n
e
-
S
~
n
+
1
U
h
n
·
n
e
,
e
S
n
+
1
e
+
1
2
∑
e
∈
Γ
1
U
h
n
·
n
e
-
U
n
·
n
e
,
e
S
n
+
1,2
e
+
∑
E
∇
·
S
n
+
1
U
n
+
1
-
S
~
n
+
1
U
h
n
,
e
S
n
+
1
-
∑
e
∈
Γ
h
e
S
n
+
1
S
n
+
1
U
n
+
1
-
S
~
n
+
1
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
·
n
e
-
U
n
·
n
e
,
e
S
n
+
1,2
-
∑
e
∈
Γ
h
S
~
n
+
1
,
up
U
h
n
·
n
e
-
S
n
+
1
U
n
+
1
·
n
e
,
e
S
n
+
1
e
-
∑
E
c
0
∇
θ
S
n
+
1
,
∇
e
S
n
+
1
E
+
∑
e
∈
Γ
h
∪
∂
Ω
c
0
∇
θ
S
n
+
1
·
n
e
,
e
S
n
+
1
e
-
∑
e
∈
Γ
h
∪
∂
Ω
c
0
∇
e
S
n
+
1
·
n
e
,
θ
S
n
+
1
e
+
a
ω
H
h
n
-
a
ω
H
n
+
1
S
∗
,
e
S
n
+
1
-
a
ω
H
n
+
1
e
S
n
+
1
+
θ
S
n
+
1
,
e
S
n
+
1
+
a
ω
H
n
+
1
-
a
ω
H
h
n
S
h
n
+
1
,
e
S
n
+
1
-
J
2
θ
S
n
+
1
,
e
S
n
+
1
,
e
Z
b
n
+
1
-
e
Z
b
n
Δ
t
,
e
Z
b
n
+
1
=
-
Z
b
n
+
1
-
Z
b
n
Δ
t
-
Z
b
t
n
+
1
,
e
Z
b
n
+
1
+
a
ω
γ
′
e
S
n
+
1
+
θ
S
n
+
1
,
e
Z
b
n
+
1
.


### 5.2. Error Estimates

We begin this subsection by deriving the error estimate for the water depthe
H stated in (46). For this we denote the right-hand terms of (46) by I
1
+
I
2
+
⋯
+
I
9 and estimate them term by term.ForI
1, we use Cauchy-Schwarz’s inequality and Young’s inequality to give (48)
I
1
=
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
≤
ε
e
H
n
+
1
2
+
C
Δ
t
2
.By noting that U
h
n
=
e
U
n
+
θ
U
n
+
U
n and Lemma 6 (boundedness of approximate solution), we split I
2 as the sum of I
2
1
+
I
2
2
+
I
2
3 and estimate them separately; (49)
I
2
=
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
=
-
1
2
∑
E
∈
ε
h
∇
·
e
U
n
+
θ
U
n
+
U
n
,
e
H
n
+
1,2
E
=
∑
i
=
1
3
I
2
i
, in which, (50)
I
2
1
=
-
1
2
∑
E
∈
ε
h
∇
·
e
U
n
,
e
H
n
+
1,2
E
≤
C
∑
E
∈
ε
h
∇
·
e
U
n
·
e
H
n
+
1
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
C
e
H
n
+
1
2
,
I
2
2
=
-
1
2
∑
E
∈
ε
h
∇
·
θ
U
n
,
e
H
n
+
1,2
E
≤
C
∑
E
∈
ε
h
∇
·
θ
U
n
·
e
H
n
+
1
≤
C
e
H
n
+
1
2
+
C
∑
E
∈
ε
h
∇
θ
U
n
2
,
I
2
3
=
-
1
2
∑
E
∈
ε
h
∇
·
U
n
,
e
H
n
+
1,2
E
≤
C
e
H
n
+
1
2
. Collecting the estimates for these 3 terms, we obtain the bound of I
2: (51)
I
2
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
C
e
H
n
+
1
2
+
C
∑
E
∈
ε
h
∇
θ
U
n
2
.Similarly, for e on the interior edge, [
U
h
n
]
=
[
e
U
n
]
+
[
θ
U
n
]
+
[
U
n
]
=
[
e
U
n
]
+
[
θ
U
n
], so we split I
3 as the sum of I
3
1
+
I
3
2 and estimate them separately: (52)
I
3
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
e
U
n
+
θ
U
n
·
n
e
e
=
∑
i
=
1
2
I
3
i
, in which (53)
I
3
1
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
e
U
n
·
n
e
e
≤
C
∑
e
∈
Γ
h
e
U
n
0
,
e
·
e
H
n
+
1
0
,
e
≤
C
∑
e
∈
Γ
h
σ
e
1
e
1
/
2
e
U
n
0
,
e
·
σ
e
1
e
-
1
/
2
e
H
n
+
1
0
,
e
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
,
I
3
2
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
θ
U
n
·
n
e
e
≤
C
∑
e
∈
Γ
h
θ
U
n
0
,
e
·
e
H
n
+
1
0
,
e
≤
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
.Collecting the estimates for these 2 terms, we obtain the bound of I
3
(54)
I
3
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
. Similarly, we obtain the bounds for I
4
,
I
8, and I
9
(55)
I
4
+
I
8
+
I
9
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
+
C
h
-
2
∑
E
θ
H
n
+
1
0
,
E
+
h
E
∇
θ
H
n
+
1
0
,
E
2
. For I
5, we use the identity (56)
H
n
+
1
U
n
+
1
-
H
~
n
+
1
U
h
n
=
-
θ
H
n
+
1
U
n
+
1
+
H
~
n
+
1
U
n
+
1
-
U
n
+
H
~
n
+
1
U
n
-
U
h
n
, and the triangle inequality to give(57)
I
5
≤
C
e
H
n
+
1
2
+
∑
E
∈
ε
h
θ
H
n
+
1
H
1
2
+
Δ
t
2
+
∑
E
∈
ε
h
θ
U
n
H
1
2
+
e
U
n
2
+
ε
∑
E
∈
ε
h
∇
e
U
n
2
. For I
6 and I
7, we apply the equality (58)
a
b
=
a
b
+
a
b
,
a
b
=
a
b
+
1
4
a
b
, to have (59)
I
6
+
I
7
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
H
n
+
1
0
,
E
+
h
E
∇
θ
H
n
+
1
0
,
E
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
. Combining the bounds (48)–(59) for I
i
(
i
=
1,2
,
…
,
9
) to obtain the estimate for the right-hand terms then applying (
1
/
2
)
(
a
2
-
b
2
)
≤
(
a
-
b
)
a to the first term on the left-hand side of (46), we finally have the estimate for the water depth e
H
(60)
1
2
Δ
t
e
H
n
+
1
2
-
e
H
n
2
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
J
1
e
U
n
,
e
U
n
+
C
Δ
t
2
+
C
e
H
n
+
1
2
+
∑
E
∈
ε
h
θ
H
n
+
1
H
1
2
+
∑
E
∈
ε
h
θ
U
n
H
1
2
+
e
U
n
2
+
C
h
-
2
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
.Combining the hypotheses (H2)–(H6) and Lemmas1–6 and using the same method and technique as these used in (48)–(59), we can obtain the error estimates of e
U
,
e
S
,
e
Z
b for the rest of the error equations (47). Here we only give the results: (61)
1
2
Δ
t
e
U
n
+
1
2
-
e
U
n
2
+
τ
1
/
2
e
U
n
+
1
2
+
∑
E
μ
1
/
2
∇
e
U
n
+
1
2
+
J
1
e
U
n
+
1
,
e
U
n
+
1
≤
ε
∑
E
∈
ε
h
∇
e
U
n
+
1
2
+
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
e
H
n
2
+
e
U
n
+
1
2
+
e
U
n
2
+
e
Z
b
n
+
1
2
+
C
∑
E
θ
U
n
+
1
H
1
2
+
θ
H
n
+
1
2
+
θ
H
n
2
+
θ
U
n
2
+
θ
Z
b
n
+
1
2
+
C
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
∑
E
θ
Z
b
n
+
1
L
2
E
+
h
E
∇
θ
Z
b
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
U
n
+
1
L
2
E
+
h
E
Δ
θ
U
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
+
1
L
2
E
+
h
E
∇
θ
U
n
+
1
L
2
E
2
+
C
Δ
t
2
,
1
2
Δ
t
e
S
n
+
1
2
-
e
S
n
2
+
∑
E
c
0
1
2
∇
e
S
n
+
1
2
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
≤
ε
∑
E
∈
ε
h
∇
e
S
n
+
1
2
+
∑
E
∈
ε
h
∇
e
U
n
2
+
e
S
n
+
1
2
+
J
1
e
U
n
,
e
U
n
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
C
e
H
n
2
+
e
U
n
2
+
e
S
n
+
1
2
+
∑
E
θ
U
n
H
1
2
+
∑
E
θ
S
n
+
1
H
1
2
+
θ
H
n
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
+
C
h
-
2
∑
E
θ
S
n
+
1
L
2
E
+
h
E
∇
θ
S
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
S
n
+
1
L
2
E
+
h
E
Δ
θ
S
n
+
1
L
2
E
2
+
C
Δ
t
2
,
1
2
Δ
t
e
Z
b
n
+
1
2
-
e
Z
b
n
2
≤
C
e
S
n
+
1
2
+
e
Z
b
n
+
1
2
+
θ
S
n
+
1
2
+
θ
Z
b
n
+
1
2
+
C
Δ
t
2
.Considering (1)–(4) is a coupled sedimentation model, we need to combine (60) and (61) to get the error estimates of e
H
,
e
U
,
e
S
,
e
Z
b. For sufficiently small ε, we have (62)
1
2
Δ
t
e
H
n
+
1
2
-
e
H
n
2
+
1
2
Δ
t
e
U
n
+
1
2
-
e
U
n
2
+
1
2
Δ
t
e
S
n
+
1
2
-
e
S
n
2
+
1
2
Δ
t
e
Z
b
n
+
1
2
-
e
Z
b
n
2
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
τ
1
/
2
e
U
n
+
1
2
+
∑
E
μ
1
/
2
∇
e
U
n
+
1
2
+
∑
E
c
0
1
/
2
∇
e
S
n
+
1
2
+
J
1
e
U
n
+
1
,
e
U
n
+
1
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
≤
C
e
H
n
+
1
2
+
e
H
n
2
+
e
U
n
+
1
2
+
e
U
n
2
+
e
S
n
+
1
2
+
e
Z
b
n
+
1
2
+
C
∑
E
θ
H
n
+
1
H
1
2
+
θ
U
n
+
1
H
1
2
+
θ
U
n
H
1
2
+
θ
S
n
+
1
H
1
2
+
θ
Z
b
n
+
1
H
1
2
+
C
h
-
2
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
+
1
L
2
E
+
h
E
∇
θ
U
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
+
C
h
-
2
∑
E
θ
S
n
+
1
L
2
E
+
h
E
∇
θ
S
n
+
1
L
2
E
2
+
C
∑
E
θ
Z
b
n
+
1
L
2
E
+
h
E
∇
θ
Z
b
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
U
n
+
1
L
2
E
+
h
E
Δ
θ
U
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
S
n
+
1
L
2
E
+
h
E
Δ
θ
S
n
+
1
L
2
E
2
+
C
Δ
t
2
. Multiplying both sides of (62) by 2
Δ
t, summing from n
=
0 to n
=
m
-
1
(
0
<
m
≤
N
), and using Lemmas 1 and 3, we have (63)
e
H
m
+
e
U
m
+
∑
n
=
0
m
-
1
Δ
t
∇
e
U
n
+
1
2
1
/
2
+
e
S
m
+
∑
n
=
0
m
-
1
Δ
t
∇
e
S
n
+
1
2
1
/
2
+
e
Z
b
m
≤
C
^
h
k
+
Δ
t
.Here, we use C
^ to replace C to denote the generic constant in order to prove Lemma 6.By using the triangle inequality and Lemma1, we can get the main conclusion for the L
2 norm error estimate of the UWDG scheme. But the process of (63) needs Lemma 6 (boundedness of approximate solution); now we use an induction argument to prove Lemma 6 for k
>
1.Proof.
We shall use the induction argument to Lemma6. For n
=
0, considering the definitions of H
h
0
,
U
h
0
,
S
h
0
,
Z
b
h
0 and Lemma 1, it is easy to get that the lemma holds valid, that is to say, (64)
a
H
h
0
L
∞
Ω
≥
H
∗
4
,
(b)
e
H
0
L
∞
Ω
+
e
U
0
L
∞
Ω
2
+
e
S
0
L
∞
Ω
≤
H
∗
4
.Assuming Lemma 6 holds for n
=
m
(
m
=
0,1
,
2
,
…
,
N
-
1
), then repeat the derivation process of (63) above; for n
=
m
+
1, we obtain (65)
e
H
m
+
1
+
e
U
m
+
1
+
e
S
m
+
1
+
e
Z
b
m
+
1
≤
C
^
h
k
+
Δ
t
.By the converse inequality (Lemma 4), we have (66)
e
H
m
+
1
L
∞
Ω
+
e
U
m
+
1
L
∞
Ω
2
+
e
S
m
+
1
L
∞
Ω
≤
C
^
K
0
h
-
1
h
k
+
Δ
t
. Here C
^ is the generic constant in (63) and K
0 is the constant in Lemma 4.
Leth
0
=
min
⁡
{
1
,
H
∗
/
8
C
^
K
0
}; then when h
,
t satisfy 0
<
h
≤
h
0 and Δ
t
=
o
(
h
), we obtain (67)
e
H
m
+
1
L
∞
Ω
+
e
U
m
+
1
L
∞
Ω
2
+
e
S
m
+
1
L
∞
Ω
≤
H
∗
4
. On the other hand, combining Lemma 5 and the triangle inequality, we find (68)
H
h
m
+
1
L
∞
Ω
≥
H
~
m
+
1
L
∞
Ω
-
e
H
m
+
1
L
∞
Ω
≥
H
∗
2
-
C
^
K
0
h
-
1
h
k
+
Δ
t
≥
H
∗
4
.
Therefore, Lemma6 holds and (63) is true for ∀
n
(
0
≤
n
≤
N
). Considering the above analysis result and combining the triangle inequality v
-
v
h
≤
e
v
+
θ
v and Lemma 1, we can prove that the approximate solutions for the depth of fluid H, the velocity U, the sediment concentration S, and the bed height Z
b are of O
(
h
k
+
Δ
t
) convergence rates, respectively.Theorem 9.
Assume that the hypotheses (H1)–(H6) hold,k
>
1, h, and Δ
t are sufficiently small, there is constant C independent of h and Δ
t, and Δ
t
=
o
(
h
), such that the following suboptimal a priori error estimate holds: (69)
H
n
-
H
h
n
+
U
n
-
U
h
n
+
S
n
-
S
h
n
+
Z
b
n
-
Z
b
h
n
+
e
S
m
+
∑
j
=
1
n
Δ
t
∇
U
j
-
U
h
j
2
1
/
2
+
∑
j
=
1
n
Δ
t
∇
S
j
-
S
h
j
2
1
/
2
≤
C
h
k
+
Δ
t
.


## 5.1. Error Equations

First, we derive the error equation of the flow continuous equation (1). Choosing the test function φ
=
e
H
n
+
1 in (10) and (14), subtracting (10) from (14), and using the definition of θ
H in (31), we derive the error equation for the flow continuous equation as follows: (42)
e
H
n
+
1
-
e
H
n
Δ
t
,
e
H
n
+
1
-
∑
E
∈
ε
h
e
H
n
+
1
U
h
n
,
∇
e
H
n
+
1
E
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
-
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
+
∑
E
∈
ε
h
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
∇
e
H
n
+
1
-
∑
e
∈
Γ
h
H
~
n
+
1
,
up
U
h
n
·
n
e
-
H
n
+
1
U
n
+
1
·
n
e
,
e
H
n
+
1
e
-
∑
e
∈
Γ
1
H
^
n
+
1
U
h
n
·
n
e
-
H
^
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
-
∑
e
∈
Γ
2
∪
Γ
3
H
h
n
+
1
U
h
n
·
n
e
-
H
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
. Taking into account that there is no diffusion term in the flow continuous equation (1), we apply Green’s formula to the second term on the left-hand side (43)
-
∑
E
∈
ε
h
e
H
n
+
1
U
h
n
,
∇
e
H
n
+
1
E
=
-
1
2
∑
E
∈
ε
h
U
h
n
,
∇
e
H
n
+
1
2
E
=
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
∂
Ω
U
h
n
·
n
e
,
e
H
n
+
1,2
e
. Using [
a
b
]
=
[
a
]
{
b
}
+
{
a
}
[
b
], we obtain (44)
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
U
h
n
·
n
e
,
1
e
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
+
∑
e
∈
Γ
h
e
H
n
+
1
,
up
U
h
n
·
n
e
,
e
H
n
+
1
e
=
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
-
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
. Similarly, we use Green’s formula for the second term on the right-hand side; thus (45)
∑
E
∈
ε
h
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
∇
e
H
n
+
1
=
-
∑
E
∈
ε
h
∇
·
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
e
H
n
+
1
+
∑
e
∈
Γ
1
H
~
n
+
1
U
h
n
-
H
^
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
+
∑
e
∈
Γ
h
e
H
n
+
1
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
·
n
e
,
1
e
+
∑
e
∈
Γ
2
∪
Γ
3
H
~
n
+
1
U
h
n
-
H
n
+
1
U
^
n
+
1
·
n
e
,
e
H
n
+
1
e
. Combining (43)–(45), we rewrite the error equation (42) into the following form: (46)
e
H
n
+
1
-
e
H
n
Δ
t
,
e
H
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
=
-
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
+
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
+
1
2
∑
e
∈
Γ
1
U
h
n
-
U
^
n
·
n
e
,
e
H
n
+
1,2
e
-
∑
E
∈
ε
h
∇
·
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
,
e
H
n
+
1
+
∑
e
∈
Γ
h
e
H
n
+
1
H
~
n
+
1
U
h
n
-
H
n
+
1
U
n
+
1
·
n
e
,
1
e
-
∑
e
∈
Γ
h
H
~
n
+
1
,
up
U
h
n
·
n
e
-
H
n
+
1
U
n
+
1
·
n
e
,
e
H
n
+
1
e
+
∑
e
∈
Γ
1
H
~
n
+
1
-
H
^
n
+
1
U
h
n
·
n
e
,
e
H
n
+
1
e
-
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
-
U
^
n
·
n
e
,
e
H
n
+
1,2
.Similarly, we obtain the other three error equations (47)
e
U
n
+
1
-
e
U
n
Δ
t
,
e
U
n
+
1
+
τ
n
+
1
e
U
n
+
1
,
e
U
n
+
1
+
∑
E
μ
∇
e
U
n
+
1
,
∇
e
U
n
+
1
E
+
J
1
e
U
n
+
1
,
e
U
n
+
1
=
-
U
n
+
1
-
U
n
Δ
t
-
U
t
n
+
1
,
e
U
n
+
1
+
∑
E
U
n
+
1
·
∇
U
n
+
1
-
U
h
n
·
∇
U
h
n
+
1
,
e
U
n
+
1
E
-
∑
E
U
h
n
·
n
e
U
h
n
+
1
,
i
n
t
-
U
h
n
+
1
,
e
x
t
,
e
U
n
+
1
,
i
n
t
∂
E
-
-
τ
n
+
1
θ
U
n
+
1
,
e
U
n
+
1
+
τ
n
+
1
-
τ
h
n
U
h
n
+
1
,
e
U
n
+
1
-
∑
E
μ
∇
θ
U
n
+
1
,
∇
e
U
n
+
1
E
+
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
θ
U
n
+
1
·
n
e
,
e
U
n
+
1
e
-
∑
e
∈
Γ
h
∪
∂
Ω
μ
∇
e
U
n
+
1
·
n
e
,
θ
U
n
+
1
e
-
J
1
θ
U
n
+
1
,
e
U
n
+
1
+
∑
E
e
H
n
+
1
,
∇
e
U
n
+
1
+
∑
E
θ
H
n
+
1
,
∇
e
U
n
+
1
+
∑
E
e
Z
b
n
+
1
,
∇
e
U
n
+
1
+
∑
E
θ
Z
b
n
+
1
,
∇
e
U
n
+
1
-
∑
e
∈
Γ
h
∪
Γ
2
∪
Γ
3
g
e
H
n
+
1
+
θ
H
n
+
1
,
e
U
n
+
1
·
n
e
e
-
∑
e
∈
Γ
h
∪
∂
Ω
g
e
Z
b
n
+
1
+
θ
Z
b
n
+
1
,
e
U
n
+
1
·
n
e
e
,
e
S
n
+
1
-
e
S
n
Δ
t
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
∑
E
c
0
∇
e
S
n
+
1
,
∇
e
S
n
+
1
E
+
J
2
e
S
n
+
1
,
e
S
n
+
1
=
-
S
n
+
1
-
S
n
Δ
t
-
S
t
n
+
1
,
e
S
n
+
1
+
∑
E
S
h
n
+
1
∇
·
U
h
n
-
S
n
+
1
∇
·
U
n
+
1
,
e
S
n
+
1
E
-
1
2
∑
E
∇
·
U
h
n
,
e
S
n
+
1,2
+
1
2
∑
e
∈
Γ
h
e
S
n
+
1,2
,
U
h
n
·
n
e
e
-
∑
e
∈
Γ
1
S
^
n
+
1
U
h
n
·
n
e
-
S
~
n
+
1
U
h
n
·
n
e
,
e
S
n
+
1
e
+
1
2
∑
e
∈
Γ
1
U
h
n
·
n
e
-
U
n
·
n
e
,
e
S
n
+
1,2
e
+
∑
E
∇
·
S
n
+
1
U
n
+
1
-
S
~
n
+
1
U
h
n
,
e
S
n
+
1
-
∑
e
∈
Γ
h
e
S
n
+
1
S
n
+
1
U
n
+
1
-
S
~
n
+
1
U
h
n
·
n
e
,
1
e
-
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
h
n
·
n
e
-
U
n
·
n
e
,
e
S
n
+
1,2
-
∑
e
∈
Γ
h
S
~
n
+
1
,
up
U
h
n
·
n
e
-
S
n
+
1
U
n
+
1
·
n
e
,
e
S
n
+
1
e
-
∑
E
c
0
∇
θ
S
n
+
1
,
∇
e
S
n
+
1
E
+
∑
e
∈
Γ
h
∪
∂
Ω
c
0
∇
θ
S
n
+
1
·
n
e
,
e
S
n
+
1
e
-
∑
e
∈
Γ
h
∪
∂
Ω
c
0
∇
e
S
n
+
1
·
n
e
,
θ
S
n
+
1
e
+
a
ω
H
h
n
-
a
ω
H
n
+
1
S
∗
,
e
S
n
+
1
-
a
ω
H
n
+
1
e
S
n
+
1
+
θ
S
n
+
1
,
e
S
n
+
1
+
a
ω
H
n
+
1
-
a
ω
H
h
n
S
h
n
+
1
,
e
S
n
+
1
-
J
2
θ
S
n
+
1
,
e
S
n
+
1
,
e
Z
b
n
+
1
-
e
Z
b
n
Δ
t
,
e
Z
b
n
+
1
=
-
Z
b
n
+
1
-
Z
b
n
Δ
t
-
Z
b
t
n
+
1
,
e
Z
b
n
+
1
+
a
ω
γ
′
e
S
n
+
1
+
θ
S
n
+
1
,
e
Z
b
n
+
1
.


## 5.2. Error Estimates

We begin this subsection by deriving the error estimate for the water depthe
H stated in (46). For this we denote the right-hand terms of (46) by I
1
+
I
2
+
⋯
+
I
9 and estimate them term by term.ForI
1, we use Cauchy-Schwarz’s inequality and Young’s inequality to give (48)
I
1
=
H
n
+
1
-
H
n
Δ
t
-
H
t
n
+
1
,
e
H
n
+
1
≤
ε
e
H
n
+
1
2
+
C
Δ
t
2
.By noting that U
h
n
=
e
U
n
+
θ
U
n
+
U
n and Lemma 6 (boundedness of approximate solution), we split I
2 as the sum of I
2
1
+
I
2
2
+
I
2
3 and estimate them separately; (49)
I
2
=
-
1
2
∑
E
∈
ε
h
∇
·
U
h
n
,
e
H
n
+
1,2
E
=
-
1
2
∑
E
∈
ε
h
∇
·
e
U
n
+
θ
U
n
+
U
n
,
e
H
n
+
1,2
E
=
∑
i
=
1
3
I
2
i
, in which, (50)
I
2
1
=
-
1
2
∑
E
∈
ε
h
∇
·
e
U
n
,
e
H
n
+
1,2
E
≤
C
∑
E
∈
ε
h
∇
·
e
U
n
·
e
H
n
+
1
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
C
e
H
n
+
1
2
,
I
2
2
=
-
1
2
∑
E
∈
ε
h
∇
·
θ
U
n
,
e
H
n
+
1,2
E
≤
C
∑
E
∈
ε
h
∇
·
θ
U
n
·
e
H
n
+
1
≤
C
e
H
n
+
1
2
+
C
∑
E
∈
ε
h
∇
θ
U
n
2
,
I
2
3
=
-
1
2
∑
E
∈
ε
h
∇
·
U
n
,
e
H
n
+
1,2
E
≤
C
e
H
n
+
1
2
. Collecting the estimates for these 3 terms, we obtain the bound of I
2: (51)
I
2
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
C
e
H
n
+
1
2
+
C
∑
E
∈
ε
h
∇
θ
U
n
2
.Similarly, for e on the interior edge, [
U
h
n
]
=
[
e
U
n
]
+
[
θ
U
n
]
+
[
U
n
]
=
[
e
U
n
]
+
[
θ
U
n
], so we split I
3 as the sum of I
3
1
+
I
3
2 and estimate them separately: (52)
I
3
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
U
h
n
·
n
e
e
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
e
U
n
+
θ
U
n
·
n
e
e
=
∑
i
=
1
2
I
3
i
, in which (53)
I
3
1
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
e
U
n
·
n
e
e
≤
C
∑
e
∈
Γ
h
e
U
n
0
,
e
·
e
H
n
+
1
0
,
e
≤
C
∑
e
∈
Γ
h
σ
e
1
e
1
/
2
e
U
n
0
,
e
·
σ
e
1
e
-
1
/
2
e
H
n
+
1
0
,
e
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
,
I
3
2
=
1
2
∑
e
∈
Γ
h
e
H
n
+
1,2
,
θ
U
n
·
n
e
e
≤
C
∑
e
∈
Γ
h
θ
U
n
0
,
e
·
e
H
n
+
1
0
,
e
≤
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
.Collecting the estimates for these 2 terms, we obtain the bound of I
3
(54)
I
3
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
. Similarly, we obtain the bounds for I
4
,
I
8, and I
9
(55)
I
4
+
I
8
+
I
9
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
+
C
h
-
2
∑
E
θ
H
n
+
1
0
,
E
+
h
E
∇
θ
H
n
+
1
0
,
E
2
. For I
5, we use the identity (56)
H
n
+
1
U
n
+
1
-
H
~
n
+
1
U
h
n
=
-
θ
H
n
+
1
U
n
+
1
+
H
~
n
+
1
U
n
+
1
-
U
n
+
H
~
n
+
1
U
n
-
U
h
n
, and the triangle inequality to give(57)
I
5
≤
C
e
H
n
+
1
2
+
∑
E
∈
ε
h
θ
H
n
+
1
H
1
2
+
Δ
t
2
+
∑
E
∈
ε
h
θ
U
n
H
1
2
+
e
U
n
2
+
ε
∑
E
∈
ε
h
∇
e
U
n
2
. For I
6 and I
7, we apply the equality (58)
a
b
=
a
b
+
a
b
,
a
b
=
a
b
+
1
4
a
b
, to have (59)
I
6
+
I
7
≤
ε
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
C
h
-
2
∑
E
θ
H
n
+
1
0
,
E
+
h
E
∇
θ
H
n
+
1
0
,
E
2
+
C
h
-
2
∑
E
θ
U
n
0
,
E
+
h
E
∇
θ
U
n
0
,
E
2
. Combining the bounds (48)–(59) for I
i
(
i
=
1,2
,
…
,
9
) to obtain the estimate for the right-hand terms then applying (
1
/
2
)
(
a
2
-
b
2
)
≤
(
a
-
b
)
a to the first term on the left-hand side of (46), we finally have the estimate for the water depth e
H
(60)
1
2
Δ
t
e
H
n
+
1
2
-
e
H
n
2
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
≤
ε
∑
E
∈
ε
h
∇
e
U
n
2
+
J
1
e
U
n
,
e
U
n
+
C
Δ
t
2
+
C
e
H
n
+
1
2
+
∑
E
∈
ε
h
θ
H
n
+
1
H
1
2
+
∑
E
∈
ε
h
θ
U
n
H
1
2
+
e
U
n
2
+
C
h
-
2
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
.Combining the hypotheses (H2)–(H6) and Lemmas1–6 and using the same method and technique as these used in (48)–(59), we can obtain the error estimates of e
U
,
e
S
,
e
Z
b for the rest of the error equations (47). Here we only give the results: (61)
1
2
Δ
t
e
U
n
+
1
2
-
e
U
n
2
+
τ
1
/
2
e
U
n
+
1
2
+
∑
E
μ
1
/
2
∇
e
U
n
+
1
2
+
J
1
e
U
n
+
1
,
e
U
n
+
1
≤
ε
∑
E
∈
ε
h
∇
e
U
n
+
1
2
+
J
1
e
U
n
,
e
U
n
+
C
e
H
n
+
1
2
+
e
H
n
2
+
e
U
n
+
1
2
+
e
U
n
2
+
e
Z
b
n
+
1
2
+
C
∑
E
θ
U
n
+
1
H
1
2
+
θ
H
n
+
1
2
+
θ
H
n
2
+
θ
U
n
2
+
θ
Z
b
n
+
1
2
+
C
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
∑
E
θ
Z
b
n
+
1
L
2
E
+
h
E
∇
θ
Z
b
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
U
n
+
1
L
2
E
+
h
E
Δ
θ
U
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
+
1
L
2
E
+
h
E
∇
θ
U
n
+
1
L
2
E
2
+
C
Δ
t
2
,
1
2
Δ
t
e
S
n
+
1
2
-
e
S
n
2
+
∑
E
c
0
1
2
∇
e
S
n
+
1
2
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
≤
ε
∑
E
∈
ε
h
∇
e
S
n
+
1
2
+
∑
E
∈
ε
h
∇
e
U
n
2
+
e
S
n
+
1
2
+
J
1
e
U
n
,
e
U
n
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
C
e
H
n
2
+
e
U
n
2
+
e
S
n
+
1
2
+
∑
E
θ
U
n
H
1
2
+
∑
E
θ
S
n
+
1
H
1
2
+
θ
H
n
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
+
C
h
-
2
∑
E
θ
S
n
+
1
L
2
E
+
h
E
∇
θ
S
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
S
n
+
1
L
2
E
+
h
E
Δ
θ
S
n
+
1
L
2
E
2
+
C
Δ
t
2
,
1
2
Δ
t
e
Z
b
n
+
1
2
-
e
Z
b
n
2
≤
C
e
S
n
+
1
2
+
e
Z
b
n
+
1
2
+
θ
S
n
+
1
2
+
θ
Z
b
n
+
1
2
+
C
Δ
t
2
.Considering (1)–(4) is a coupled sedimentation model, we need to combine (60) and (61) to get the error estimates of e
H
,
e
U
,
e
S
,
e
Z
b. For sufficiently small ε, we have (62)
1
2
Δ
t
e
H
n
+
1
2
-
e
H
n
2
+
1
2
Δ
t
e
U
n
+
1
2
-
e
U
n
2
+
1
2
Δ
t
e
S
n
+
1
2
-
e
S
n
2
+
1
2
Δ
t
e
Z
b
n
+
1
2
-
e
Z
b
n
2
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
H
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
H
n
+
1,2
e
+
τ
1
/
2
e
U
n
+
1
2
+
∑
E
μ
1
/
2
∇
e
U
n
+
1
2
+
∑
E
c
0
1
/
2
∇
e
S
n
+
1
2
+
J
1
e
U
n
+
1
,
e
U
n
+
1
+
J
2
e
S
n
+
1
,
e
S
n
+
1
+
1
2
∑
e
∈
Γ
h
U
h
n
·
n
e
,
e
S
n
+
1
2
e
+
1
2
∑
e
∈
Γ
1
U
^
n
·
n
e
,
e
S
n
+
1,2
e
+
1
2
∑
e
∈
Γ
2
∪
Γ
3
U
^
n
·
n
e
,
e
S
n
+
1,2
e
≤
C
e
H
n
+
1
2
+
e
H
n
2
+
e
U
n
+
1
2
+
e
U
n
2
+
e
S
n
+
1
2
+
e
Z
b
n
+
1
2
+
C
∑
E
θ
H
n
+
1
H
1
2
+
θ
U
n
+
1
H
1
2
+
θ
U
n
H
1
2
+
θ
S
n
+
1
H
1
2
+
θ
Z
b
n
+
1
H
1
2
+
C
h
-
2
∑
E
θ
H
n
+
1
L
2
E
+
h
E
∇
θ
H
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
+
1
L
2
E
+
h
E
∇
θ
U
n
+
1
L
2
E
2
+
C
h
-
2
∑
E
θ
U
n
L
2
E
+
h
E
∇
θ
U
n
L
2
E
2
+
C
h
-
2
∑
E
θ
S
n
+
1
L
2
E
+
h
E
∇
θ
S
n
+
1
L
2
E
2
+
C
∑
E
θ
Z
b
n
+
1
L
2
E
+
h
E
∇
θ
Z
b
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
U
n
+
1
L
2
E
+
h
E
Δ
θ
U
n
+
1
L
2
E
2
+
C
∑
E
∇
θ
S
n
+
1
L
2
E
+
h
E
Δ
θ
S
n
+
1
L
2
E
2
+
C
Δ
t
2
. Multiplying both sides of (62) by 2
Δ
t, summing from n
=
0 to n
=
m
-
1
(
0
<
m
≤
N
), and using Lemmas 1 and 3, we have (63)
e
H
m
+
e
U
m
+
∑
n
=
0
m
-
1
Δ
t
∇
e
U
n
+
1
2
1
/
2
+
e
S
m
+
∑
n
=
0
m
-
1
Δ
t
∇
e
S
n
+
1
2
1
/
2
+
e
Z
b
m
≤
C
^
h
k
+
Δ
t
.Here, we use C
^ to replace C to denote the generic constant in order to prove Lemma 6.By using the triangle inequality and Lemma1, we can get the main conclusion for the L
2 norm error estimate of the UWDG scheme. But the process of (63) needs Lemma 6 (boundedness of approximate solution); now we use an induction argument to prove Lemma 6 for k
>
1.Proof.
We shall use the induction argument to Lemma6. For n
=
0, considering the definitions of H
h
0
,
U
h
0
,
S
h
0
,
Z
b
h
0 and Lemma 1, it is easy to get that the lemma holds valid, that is to say, (64)
a
H
h
0
L
∞
Ω
≥
H
∗
4
,
(b)
e
H
0
L
∞
Ω
+
e
U
0
L
∞
Ω
2
+
e
S
0
L
∞
Ω
≤
H
∗
4
.Assuming Lemma 6 holds for n
=
m
(
m
=
0,1
,
2
,
…
,
N
-
1
), then repeat the derivation process of (63) above; for n
=
m
+
1, we obtain (65)
e
H
m
+
1
+
e
U
m
+
1
+
e
S
m
+
1
+
e
Z
b
m
+
1
≤
C
^
h
k
+
Δ
t
.By the converse inequality (Lemma 4), we have (66)
e
H
m
+
1
L
∞
Ω
+
e
U
m
+
1
L
∞
Ω
2
+
e
S
m
+
1
L
∞
Ω
≤
C
^
K
0
h
-
1
h
k
+
Δ
t
. Here C
^ is the generic constant in (63) and K
0 is the constant in Lemma 4.
Leth
0
=
min
⁡
{
1
,
H
∗
/
8
C
^
K
0
}; then when h
,
t satisfy 0
<
h
≤
h
0 and Δ
t
=
o
(
h
), we obtain (67)
e
H
m
+
1
L
∞
Ω
+
e
U
m
+
1
L
∞
Ω
2
+
e
S
m
+
1
L
∞
Ω
≤
H
∗
4
. On the other hand, combining Lemma 5 and the triangle inequality, we find (68)
H
h
m
+
1
L
∞
Ω
≥
H
~
m
+
1
L
∞
Ω
-
e
H
m
+
1
L
∞
Ω
≥
H
∗
2
-
C
^
K
0
h
-
1
h
k
+
Δ
t
≥
H
∗
4
.
Therefore, Lemma6 holds and (63) is true for ∀
n
(
0
≤
n
≤
N
). Considering the above analysis result and combining the triangle inequality v
-
v
h
≤
e
v
+
θ
v and Lemma 1, we can prove that the approximate solutions for the depth of fluid H, the velocity U, the sediment concentration S, and the bed height Z
b are of O
(
h
k
+
Δ
t
) convergence rates, respectively.Theorem 9.
Assume that the hypotheses (H1)–(H6) hold,k
>
1, h, and Δ
t are sufficiently small, there is constant C independent of h and Δ
t, and Δ
t
=
o
(
h
), such that the following suboptimal a priori error estimate holds: (69)
H
n
-
H
h
n
+
U
n
-
U
h
n
+
S
n
-
S
h
n
+
Z
b
n
-
Z
b
h
n
+
e
S
m
+
∑
j
=
1
n
Δ
t
∇
U
j
-
U
h
j
2
1
/
2
+
∑
j
=
1
n
Δ
t
∇
S
j
-
S
h
j
2
1
/
2
≤
C
h
k
+
Δ
t
.


## 6. Algorithm Flow and Computational Cost

In dealing with practical issues, there are three steps: collecting the river terrain data, selecting the hydrologic parameter, and program calculation. These steps are carefully displayed in the following text and Figure2.Figure 2
The flow chart.Step 1.
Collect the river terrain data and determine the computational domain, the initial-boundary conditions, and the time range.Step 2.
Obtain the hydrologic parameters by properly selecting the sediment carrying capacity, Manning’s roughness coefficient, the viscosity coefficient, the sediment diffusion coefficient, the coefficient of saturation recovery, the settling velocity of sediment, and the sediment dry bulk density.Step 3.
Program calculation step includes(1)
determining the time and space step to get triangular subdivision;
(2)
calculatingH
h
n
+
1
,
U
h
n
+
1
,
S
h
n
+
1
,
Z
b
h
n
+
1, and the unknowns on the n
+
1 time level, by inserting H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n, and the knowns on the n time level into the fully upwind discontinuous Galerkin finite element procedure (UWDG);
(3)
judging whether the scheduled time can be reached, if it is reached, output the results and then the program ends.Now we analyze the computational cost of the UWDG scheme. Because of the lack of continuity constraints between mesh elements for the test functions, we must define the basic functions for each triangle element. We usually use the monomial functions as the basic functions; in 2D, the basic functions areϕ
i
(
x
,
y
)
=
x
I
y
J
,
I
+
J
=
i
,
and
0
≤
i
≤
k
. Let N denote the total number of elements E; if k
=
1,2 and 3, then the total degrees of freedom is 3
N
,
6
N, and 10
N, respectively. N can be understood as the scale of the computation, and N
=
O
(
h
-
2
) in 2D; here h denote the maximum diameter of elements.From Step3 above, in order to calculate the unknowns on the n
+
1 time level of H
h
n
+
1
,
U
h
n
+
1
,
S
h
n
+
1
,
Z
b
h
n
+
1, we need the knowns on the n time level of H
h
n
,
U
h
n
,
S
h
n
,
Z
b
h
n. Thus for n
=
0,1
,
…
,
N, the computational complexity T
(
n
+
1
) has (70)
T
n
=
0
=
O
1
,
T
n
+
1
=
T
n
+
O
1
.There are three steps to calculating then time level: generate meshes, form the total coefficient matrix, and solve linear equations. Let T
1
(
n
)
,
T
2
(
n
), and T
3
(
n
) denote the computational complexity of three steps above; if we use Delaunay triangulation algorithm to generate meshes and use iterative methods to solve linear equations, the computational complexity T
(
n
) has (71)
T
n
=
T
1
n
+
T
2
n
+
T
3
n
=
O
N
2
+
O
N
+
O
N
2
=
O
N
2
.From Theorem9, the choice of the time step must satisfy Δ
t
=
o
(
h
) to ensure the convergence of the UWDG scheme. Let Δ
t
=
O
(
h
1
+
α
)
(
α
>
0
); then N
=
T
/
Δ
t
=
O
(
h
-
1
-
α
), and we have (72)
N
=
O
N
1
+
α
/
2
. Combining (70)–(72), the computational cost is O
(
N
(
5
+
α
)
/
2
), for example, α
=
1, that is to say, Δ
t
=
O
(
h
2
), the computational cost is O
(
N
3
); α
=
2, Δ
t
=
O
(
h
3
), the computational cost is O
(
N
3.5
).


## 7. Numerical Experiments

In this section, we perform two numerical experiments, one is for unsteady state case and the other is for stationary-state case, to verity the theoretically proven convergence results.Example 1 (unsteady-state case).
Let the domainΩ be the unit square and the time interval be [
0,1
]. The analytic and smooth solution is prescribed to be(73)
H
=
t
+
1
sin
⁡
π
x
sin
⁡
π
y
+
1
,
U
=
10
x
sin
⁡
2
t
sin
⁡
π
x
sin
⁡
π
y
,
10
y
cos
⁡
2
t
sin
⁡
π
x
sin
⁡
π
y
,
S
=
t
2
+
1
sin
⁡
π
x
sin
⁡
π
y
,
Z
b
=
a
ω
γ
′
1
3
t
3
+
t
sin
⁡
π
x
sin
⁡
π
y
-
S
∗
t
.The boundary conditions, initial conditions, and the source terms f
H
,
f
U
,
f
S
,
and
f
Z
b are then completely determined. We also let μ
=
0.01, c
0
=
0.01, and the coefficient of saturation recovery a, settling velocity of sediment ω, and the sediment dry bulk density γ
′ are all equal to one.We choose the piecewise quadratic polynomials (k
=
2) and the piecewise linear polynomials (k
=
1) in the finite element space D
k
(
ε
h
), respectively. The errors and convergence orders in L
2 norm are displayed in Tables 1 and 3. The L
2
(
0,1
;
H
1
) norm and convergence orders between the approximate solution and the analytic solution of U and S are displayed in Table 2. The numbers in brackets are the theoretically predicted convergence orders.Table 1
Errors and convergence orders witht
=
1 and k
=
2.
(a)
h
Δ
t
H
-
H
h
Order
U
-
U
h
Order
1
/
2
1
/
16
1.0805e − 00
2.5258e − 01
1
/
4
1
/
128
5.7656e − 02
4.228 (2)
2.8826e − 02
3.131 (2)
1
/
8
1
/
1024
1.4262e − 02
2.015 (2)
3.6902e − 03
2.966 (2)
1
/
16
1
/
8192
3.8824e − 03
1.898 (2)
6.0486e − 04
2.789 (2)
(b)
h
Δ
t
S
-
S
h
Order
Z
b
-
Z
b
h
Order
1
/
2
1
/
16
7.0693e − 02
4.1697e − 02
1
/
4
1
/
128
8.3987e − 03
3.073 (2)
5.2262e − 03
2.996 (2)
1
/
8
1
/
1024
9.5399e − 04
3.138 (2)
6.5837e − 04
2.988 (2)
1
/
16
1
/
8192
1.2387e − 04
2.946 (2)
8.9013e − 05
2.886 (2)Table 2
L
2
(
0,1
;
H
1
) errors and convergence orders with k
=
2.
h
Δ
t
U
-
U
h
L
2
(
0,1
;
H
1
)
Order
S
-
S
h
L
2
(
0,1
;
H
1
)
Order
1
/
2
1
/
16
3.9640e − 00
9.3795e − 01
1
/
4
1
/
128
1.0055e − 00
1.979 (2)
2.5767e − 01
1.863 (2)
1
/
8
1
/
1024
2.6188e − 01
1.941 (2)
6.6384e − 02
1.956 (2)
1
/
16
1
/
8192
6.6941e − 02
1.968 (2)
1.6740e − 02
1.988 (2)Table 3
Errors and convergence orders witht
=
1 and k
=
1.
(a)
h
Δ
t
H
-
H
h
Order
U
-
U
h
Order
1
/
2
1
/
8
5.7279e − 00
1.8258e − 00
1
/
4
1
/
32
1.7836e − 00
1.683 (1)
5.6527e − 01
1.692 (1)
1
/
8
1
/
128
5.5102e − 01
1.694 (1)
1.2068e − 01
2.228 (1)
1
/
16
1
/
512
1.7294e − 01
1.671 (1)
2.6058e − 02
2.211 (1)
1
/
32
1
/
2048
6.1343e − 02
1.495 (1)
5.7693e − 03
2.175 (1)
(b)
h
Δ
t
S
-
S
h
Order
Z
b
-
Z
b
h
Order
1
/
2
1
/
8
5.1268e − 01
3.1549e − 01
1
/
4
1
/
32
1.6250e − 01
1.657 (1)
1.0736e − 01
1.555 (1)
1
/
8
1
/
128
4.1680e − 02
1.963 (1)
2.7807e − 02
1.949 (1)
1
/
16
1
/
512
1.0622e − 02
1.972 (1)
6.9953e − 03
1.991 (1)
1
/
32
1
/
2048
2.6583e − 03
1.998 (1)
1.7455e − 03
2.003 (1)By the data in Tables1–3, we see that the UWDG scheme works well and the convergence order is greater than the results theoretically predicted by (69), even for k
=
1 which is not included in Theorem 9.For simplicity of presentation, we use Figures3 and 4 to describe the contours of the sediment concentration S and the velocity field U at t
=
1, respectively, and omit those of the remaining unknowns.Figure 3
Contours ofS: (a) is the numerical solution for k
=
2 and (b) is the exact solution, t
=
1.
(a)
(b)Figure 4
Velocity field: (a) is the numerical solution fork
=
2 and (b) is the exact solution, t
=
1.
(a)
(b)Example 2 (stationary-state case).
We further consider nonconstant bathymetry but stationary case. We take the domainΩ to be the unit square and the time interval to be [
0,1
]. The analytic and smooth solution is prescribed to be(74)
H
=
sin
⁡
π
x
sin
⁡
π
y
+
1
,
U
=
10
x
sin
⁡
π
x
sin
⁡
π
y
,
10
y
sin
⁡
π
x
sin
⁡
π
y
,
S
=
sin
⁡
π
x
sin
⁡
π
y
,
Z
b
=
a
ω
γ
′
t
sin
⁡
π
x
sin
⁡
π
y
-
S
∗
t
.The boundary conditions, initial conditions, and the source terms f
H
,
f
U
,
f
S
,
and
f
Z
b are then completely determined by the analytic and smooth solution. We let μ
=
0.01 and c
0
=
0.01 and also let the coefficient of saturation recovery a, settling velocity of sediment ω, and the sediment dry bulk density γ
′ be all equal to one. Further we assume that Manning’s roughness coefficient is equal to zero; this forces the ted shear stress τ to be equal to zero and (2) turns to be a vortex without friction. The numerical results are displayed in Tables 4–6. The numbers in brackets are the theoretically predicted convergence orders.Table 4
Errors and convergence orders witht
=
1 and k
=
2.
(a)
h
Δ
t
H
-
H
h
Order
U
-
U
h
Order
1
/
2
1
/
16
1.9997e − 00
6.9881e − 01
1
/
4
1
/
128
7.8374e − 01
1.351 (2)
7.0298e − 02
3.313 (2)
1
/
8
1
/
1024
1.6547e − 01
2.243 (2)
6.6213e − 03
3.408 (2)
1
/
16
1
/
8192
3.4031e − 02
2.282 (2)
7.3425e − 04
3.172 (2)
(b)
h
Δ
t
S
-
S
h
Order
Z
b
-
Z
b
h
Order
1
/
2
1
/
16
5.6456e − 02
3.8771e − 02
1
/
4
1
/
128
4.1348e − 03
3.771 (2)
3.7886e − 03
3.355 (2)
1
/
8
1
/
1024
5.8254e − 04
2.827 (2)
4.6080e − 04
3.039 (2)
1
/
16
1
/
8192
6.1174e − 05
3.251 (2)
5.7319e − 05
3.007 (2)Table 5
L
2
(
0,1
;
H
1
) errors and convergence orders with k
=
2.
h
Δ
t
U
-
U
h
L
2
(
0,1
;
H
1
)
Order
S
-
S
h
L
2
(
0,1
;
H
1
)
Order
1
/
2
1
/
16
7.2990e − 00
5.0641e − 01
1
/
4
1
/
128
2.0015e − 00
1.866 (2)
1.3001e − 01
1.962 (2)
1
/
8
1
/
1024
4.5136e − 01
2.148 (2)
3.8540e − 02
1.754 (2)
1
/
16
1
/
8192
1.0421e − 01
2.115 (2)
8.4163e − 03
2.195 (2)Table 6
Errors and convergence orders witht
=
1 and k
=
1.
(a)
h
Δ
t
H
-
H
h
Order
U
-
U
h
Order
1
/
2
1
/
8
1.283e + 01
3.9782e − 00
1
/
4
1
/
32
3.3491e − 00
1.938 (1)
1.1136e − 00
1.837 (1)
1
/
8
1
/
128
1.0218e − 00
1.712 (1)
1.5430e − 01
2.852 (1)
1
/
16
1
/
512
3.2672e − 01
1.645 (1)
2.9011e − 02
2.411 (1)
(b)
h
Δ
t
S
-
S
h
Order
Z
b
-
Z
b
h
Order
1
/
2
1
/
8
2.9590e − 01
2.5971e − 01
1
/
4
1
/
32
1.1303e − 01
1.388 (1)
9.7248e − 02
1.417 (1)
1
/
8
1
/
128
2.7077e − 02
2.062 (1)
2.3966e − 02
2.021 (1)
1
/
16
1
/
512
5.9637e − 03
2.183 (1)
5.8229e − 03
2.041 (1)Our numerical tests show that the convergence accuracy of the proposed UWDG scheme in this paper is at least as sharp asO
(
h
k
+
Δ
t
); the results are theoretically predicted by (69), even for k
=
1 which is not included in Theorem 9.


## 8. Conclusions

In this paper, we design an upwind discontinuous Galerkin finite element method for the 2D sedimentation in an estuary as coupled nonlinear differential system, which obeys the local mass conservation and possesses good stability. There exists a unique solution to the numerical procedure and the discrete solution permitsO
(
h
k
+
Δ
t
) convergence rate. Numerical experiments are conducted to verify our theoretical findings.



---

*Source: 2487136-2017-02-20.xml*