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

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 of Lesaint - 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 (ℎ 𝑘 + Δ𝑡) 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.


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 weaktidal, 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][2][3][4][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   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   -constant [6][7][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 inserting Lesaint-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   of the fluid possesses a positive lower bound, we prove that the approximate solutions for the depth of fluid , the velocity U, the sediment concentration , and the bed height   are of O(ℎ  + Δ) 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 Sections 2-6, respectively.For simplicity of presentation, we will use  , () to denote the usual Sobolev spaces that provided the norm ‖ ⋅ ‖ ,, and the seminorm | ⋅ | ,, for any 2D domain .When  = 2 and  = Ω, the subscript  and Ω will be omitted.We also let (⋅, ⋅) and ⟨⋅, ⋅⟩ stand for the inner products on  and on the edges , respectively, and ‖ ⋅ ‖ denotes their  2 -norms. is a generic constant and  is a small positive constant, which may take different values at different places but independent of the triangulation parameters ℎ and Δ.We define   (, ; ) = { : ∫   ‖(⋅, )‖    < ∞}, where if  = ∞, the integral is replaced with the essential supreme.

Reformulation of the Mathematical Model
Let Ω be a bounded polygonal domain in  2 , and let [0, ] 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, ]; here we ignore the Coriolis force, since it does not affect the correctness of the model.
The flow continuity equation is the flow motion equation is the sediment continuity equation is and the bed deformation equation is where the unknowns ,,, = (,V)  and   =  −  denote the depth of water, the water level, the sediment concentration, the vertical mean velocity, and the bed height, respectively;  * is the sediment carrying capacity,  is the gravitational acceleration,  = ( √  2 + V 2 / 2 ),  = (1/) 1/6 denotes the Chezy coefficient,  denotes Manning's roughness coefficient,  is the viscosity coefficient,  0 is the sediment diffusion coefficient,  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 /+ ∇ ⋅ () = 0.It means that the flow continuity equation is not coupled with the bed deformation equation; in other words, the bed height   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   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 where the initial values  0 ,  0 ,  0 ,  0 and the boundary conditions Û, Ĥ, Ŝ are given data, and n denotes the unit outward normal vector to Ω.

The Upwind Discontinuous Galerkin Finite Element Procedure
In this section, we shall revisit the Broken Sobolev spaces and modify those terms involving U by Lesaint-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  ℎ be a quasi-uniform subdivision of Ω with element ; let Ω ℎ = ⋃ ∈ ℎ , ℎ  denote the diameter of , and ℎ = max{ℎ  }.Further we need to describe the edges of an element and the functions defined on them.We denote, by , any edge of an element  ∈  ℎ , and then  belongs to either its interior or boundary of Ω.
We also denote, by Γ ℎ , the set of interior edges.If  ∈ Ω, the unit normal vector n e of  is taken to be outward to Ω, or if  is a common edge of two neighbor elements  +  and  −  , we understand that n e is oriented from  −  to  +  .For any function V of  1 ( ℎ ), we denote by V ± its traces on  +  and  −  and define its average and jump by {V} = (V + + V − )/2 and [V] = V − − V + , respectively.When  ∈ Ω, 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   ( ℎ ) for any real number  ≥ 1 are defined by equipped with the norm Here ‖V‖ 2   () = (∑ ||≤ ‖  V‖ 2 0, ) 1/2 with   V are being the  th -order generalized derivatives of V.The functionals with respect to the averages and jumps on the whole edges are defined by where || is the length of  and the selected nonnegative real numbers  1  ,  2  are called penalty parameters.

Variational Formulation for the Mathematical Model.
Let  be any function in  1 ( ℎ ).We multiply (1) by suitable test functions  on each element , apply the Green formula, and then sum the resulting equations over the domain to obtain the variational form for the flow continuity equation Similarly, we obtain the variational forms for the rest of the three equations ( 2)-( 4): Consequently, the variational formulation corresponding to models (1)-( 4) is defined as to find the map {, , ,   } : - (11).

The Upwind Discontinuous Galerkin Finite Element
Scheme.Assume that  is a positive integer and   () denotes the space of polynomials of degree less than or equal to  on .Then, the finite element spaces used later are chosen to be the subspace of   ( ℎ ) defined by For time discretization we divide the time interval [0, T] by the nodes   ( = 0, 1, . . ., ) with time step Δ = / and let   (, ) be the node-function value of (, , ).
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 for  and  on the interior edge  by  up and  up , respectively.First we define the upwind values of  and  on the interior edge  by  up and  up , respectively as follows: 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 of ,  on interior edge are replaced by  up ,  up and the nonlinear convective term  ⋅ ∇ is handled by the Lesaint-Raviart upwind technique.
By doing so, the desired stability can be guaranteed.For easy computation we use, for example,    +1 instead of  +1  +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 The initial values  0 ℎ ,  0 ℎ ,  0 ℎ , and  0 ℎ also are clarified by their  2 -orthogonal projections: Here and V ext denote the interior and exterior traces of V on

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.

Some Lemmas and Hypotheses.
To conduct convergence analysis in the next section we need to define the  2orthogonal projections for the unknowns {, , ,   }.These projections are maps { H, Ũ, S, Z } : Following the practical background of a real sedimentation process, we shall make some reasonable hypotheses on the exact solutions of the model.
Next we introduce the following lemmas which are crucial to our convergence analysis.
Then, there is a constant  0 > 0 such that Lemma 5 (boundedness).Assume that ℎ is sufficiently small.Then, there is a constant  > 0 such that, for ∀ ≥ 0, Proof.Combining the definition of H, the hypotheses (H2) and (H4), and Lemma 1, we obtain Let On the other hand, by the triangle inequality, we can deduce

Stability Estimate.
Preceding the proof of the stability for ( 14)-( 18), we first present a bounded estimate for  ℎ ,  ℎ ,  ℎ , and  ℎ .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.

Convergence Analysis
In this section, we shall conduct error estimates for the UWDG solution and prove the O(ℎ  + Δ) 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  =  +1   in ( 10) and ( 14), subtracting ( 10) from ( 14), and using the definition of   in (31), we derive the error equation for the flow continuous equation as follows: 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) Similarly, we use Green's formula for the second term on the right-hand side; thus Combining ( 43)-(45), we rewrite the error equation (42) into the following form: Similarly, we obtain the other three error equations in which, in which (55) For  5 , we use the identity (61) Considering ( 1)-( 4) is a coupled sedimentation model, we need to combine (60) and (61) to get the error estimates of   ,   ,   ,    .For sufficiently small , we have Here, we use Ĉ to replace  to denote the generic constant in order to prove Lemma 6.By using the triangle inequality and Lemma 1, we can get the main conclusion for the  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  > 1.

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 Figure 2.
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) calculating  +1 ℎ ,  +1 ℎ ,  +1 ℎ ,  +1 ℎ , and the unknowns on the  + 1 time level, by inserting   ℎ ,   ℎ ,   ℎ ,   ℎ , and the knowns on the  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   (, ) =     ,  +  = , and 0 ≤  ≤ .

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.
The boundary conditions, initial conditions, and the source terms   ,   ,   , and    are then completely determined.We also let  = 0.01,  0 = 0.01, and the coefficient of saturation recovery , settling velocity of sediment , and the sediment dry bulk density   are all equal to one.
We choose the piecewise quadratic polynomials ( = 2) and the piecewise linear polynomials ( = 1) in the finite element space   ( ℎ ), respectively.The errors and convergence orders in  2 norm are displayed in Tables 1 and  3.The  2 (0, 1;  1 ) norm and convergence orders between the approximate solution and the analytic solution of  and  are displayed in Table 2.The numbers in brackets are the theoretically predicted convergence orders.
By the data in Tables 1-3, we see that the UWDG scheme works well and the convergence order is greater than the results theoretically predicted by (69), even for  = 1 which is not included in Theorem 9.
For simplicity of presentation, we use Figures 3 and 4 to describe the contours of the sediment concentration  and  the velocity field  at  = 1, respectively, and omit those of the remaining unknowns.
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 The boundary conditions, initial conditions, and the source terms   ,   ,   , and    are then completely determined by the analytic and smooth solution.We let  = 0.01 and  0 = 0.01 and also let the coefficient of saturation recovery , 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 Mathematical Problems in Engineering in Tables 4-6.The numbers in brackets are the theoretically predicted convergence orders.
Our numerical tests show that the convergence accuracy of the proposed UWDG scheme in this paper is at least as sharp as (ℎ  + Δ); the results are theoretically predicted by (69), even for  = 1 which is not included in Theorem 9.

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 permits O(ℎ  + Δ) convergence rate.Numerical experiments are conducted to verify our theoretical findings.

Figure 3 :
Figure 3: Contours of : (a) is the numerical solution for  = 2 and (b) is the exact solution,  = 1.

Figure 4 :
Figure 4: Velocity field: (a) is the numerical solution for  = 2 and (b) is the exact solution,  = 1.