Numerical modeling of submarine turbidity currents over erodible beds using unstructured grids

Accepted Manuscript

Numerical modelling of submarine turbidity currents over erodible beds using unstructured grids ´ Xin Liu, Julio Angel Infante Sedano, Abdolmajid Mohammadian PII: DOI: Reference:

S1463-5003(17)30046-X 10.1016/j.ocemod.2017.03.015 OCEMOD 1194

To appear in:

Ocean Modelling

Received date: Revised date: Accepted date:

19 December 2016 10 February 2017 30 March 2017

´ Please cite this article as: Xin Liu, Julio Angel Infante Sedano, Abdolmajid Mohammadian, Numerical modelling of submarine turbidity currents over erodible beds using unstructured grids, Ocean Modelling (2017), doi: 10.1016/j.ocemod.2017.03.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

CR IP T

7

• Time step constraints are proposed to guarantee the non-negativity of computed fluid depth and sediment concentration.

AN US

6

• A special discretization for bed-slope source term is developed in order to obtain the well-balanced property.

M

5

ED

4

• A 2-D coupled numerical model is developed to simulate turbidity current over erodible bed with wetting and drying.

PT

3

CE

2

Highlights

AC

1

1

ACCEPTED MANUSCRIPT

9

Numerical modelling of submarine turbidity currents over erodible beds using unstructured grids

10

´ Xin Liua,∗, Julio Angel Infante Sedanoa , Abdolmajid Mohammadiana

8

a

12

University of Ottawa, Ottawa, ON, K1N6N5, Canada

CR IP T

11

Abstract

ED

M

AN US

Second-order central-upwind schemes proposed by Bryson et al. [5] for the Saint-Venant system have two very attractive properties: well-balanced and positivity preserving, which are originally designed for constant fluid density and fixed beds in [5]. For the turbidity current system with variable density over erodible beds, such desired properties can be obtained by developing a well-balanced and positivity preserving central-upwind scheme following the ideas in [5]. To this end, in this paper, a coupled numerical model for twodimensional depth-averaged turbidity current system over erodible beds is developed using finite volume method on triangular grids. The proposed numerical model is second-order accurate in space using piecewise linear reconstruction and third-order accurate in time using a strong stability preserving Runge-Kutta method. Applying the central-upwind method to estimate numerical fluxes through cell interfaces, the model can successfully deal with sharp gradients in turbidity flows. The developed numerical model can preserve the well-balanced property over irregular bottom, guarantee the non-negative turbidity current depth over erodible beds, and preserving the positivity of suspended sediment. These features of the developed numerical model and its robustness and accuracy are demonstrated in several numerical tests. Keywords: Numerical modelling, turbidity current, well-balanced, positivity preserving, wetting-drying, central-upwind method, unstructured grid.

15

1. Introduction

17 18 19 20 21

CE

Turbidity current is a type of underwater density current driven by the negative buoyancy forces due to the suspended sediment. Oceanic turbidity currents can be caused by tectonic disturbances of the sea floor, and the slumping of sediment that has piled up at the top of convergent plate margins, continental slopes and submarine canyons. Their profound impacts on the continental shelves, sea floor and ocean environment have attracted intensive research interests.

AC

16

PT

14

13

∗

Corresponding author ´ Email addresses: [email protected] (Xin Liu ), [email protected] (Julio Angel Infante Sedano), [email protected] (Abdolmajid Mohammadian) Preprint submitted to Elsevier

March 31, 2017

ACCEPTED MANUSCRIPT

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64

CR IP T

28

AN US

27

M

26

ED

25

PT

24

A large number of laboratory tests have been conducted to study the turbidity currents, see [41, 12, 43, 3, 13, 42, 1, 2, 6, 16, 44, 46, 24, 37, 8, 11]. However, laboratory studies are usually constrained by their small scale and can only be roughly applied to large-scale geophysical events. Only a few field observations of turbidity currents are reported in the literature due to the unpredictability and destructive nature of such events, see [23, 17, 39]. Due to the aforementioned constrains of experimental and observational studies, numerical investigations based on theoretical and physical models have attracted significant attention because of their affordable computational cost and available details of time-dependent flow structure. Three-dimensional (3-D) models can resolve the vertical structure of turbidity current, particularly for the current front that may be considerably non-uniform in vertical. Several 3-D numerical investigation of density currents have been performed, see [21, 26, 38, 9, 40]. However, the computational costs of such 3-D numerical models are very expensive, especially when applied to large scale cases. In addition, the physics of turbulent density currents has not been well understood; this also limits the applications of the 3-D turbidity current models. An alternative way to numerically investigate the submarine turbidity current is to use the depth-averaged two dimensional (2-D) numerical models which are attractive due to the fact that such 2-D models are well-understood, of lower computational cost and easy to use. One-layer 2-D depth averaged models have been widely studied. Parker et al. [42] derived the depth-averaged governing equations for 2-D turbidity currents which provides the basic model for numerical studies. Choi [7] developed a 2-D numerical model for turbidity currents using a dissipative-Galerkin finite element method. A deforming grid generation algorithm was used to track the current front. However, such techniques may become very expensive if the solution in the far-field is of interest because a large number of elements are required to produce accurate solutions. This prevents the extension of such a model to large-scale studies. Bradford and Katopodes [4] developed a 2-D finitevolume method for turbid underflows considering non-uniform sediment. Roe scheme is used to estimate fluxes through cell interfaces. Their model didn’t include the friction of ambient water which may lead to overestimation of propagating speed. Imran and Parker [22] developed a 2-D numerical model describing the formation of a submarine fan by a spreading turbidity current using finite difference method. In those models, the bed evolves in response to the exchange of suspended sediment in the turbidity current. However, the momentum and mass exchange along with sediment exchange are not considered in their studies. To fully consider the feedback impacts of bed deformation, Hu and Cao [18] have modified the governing equations for turbidity current and extended it to 2-D cases [19, 20] based on finite volume method. Those modified equations are also adopted in the current study. One desired property of such numerical models for slender flows, especially for turbidity currents, is positivity preserving. Although applying CourantFriedrichsLewy (CFL) condition is sufficient to guarantee the numerical stability of finite volume methods, using the time step size satisfying CFL condition can not guarantee the non-negativity of current depth and volumetric concentration of sediment. None of the above-mentioned numerical models for turbidity currents incorporates the positivity preserving property. In applications, the current depth and volumetric concentration may be forced to be zero if the computed value

CE

23

AC

22

3

ACCEPTED MANUSCRIPT

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91

92

93 94 95 96 97

98

99

CR IP T

71

AN US

70

M

69

ED

68

PT

67

is negative. This leads to a loss of certain amount of mass and momentum. In this study, the positivity preserving properties for both depth and concentration of turbidity current are guaranteed using the developed numerical model. The developed numerical model in the current study is an extension of the well-balanced positivity preserving central-upwind methods developed by Bryson et al. in [5]. However, the well-balanced and positivity preserving techniques in [5] were originally designed for Saint-Venant system with constant fluid density and fixed bed. Consequently, they can not be directly applied to the studied turbidity current system. To overcome this, the well-balanced and positivity techniques in [5] are modified in our proposed schemes, and a numerical model with such desired properties is developed. In this paper, a numerical model based on the finite volume method, is developed to solve 2-D coupled depth-averaged turbidity current system consisting of sediment transport, hydrodynamic and morphodynamic models. The resulting hyperbolic system of balance laws consists of five coupled equations. The central-upwind method [31, 28, 29, 30, 27] is adopted to estimate the numerical fluxes through cell interfaces. In order to guarantee the well-balanced property, special discretization for bed-slope source term is developed in this paper. Using non-negative reconstruction for bed, desingularization of point values and special treatment at wet-dry fronts, the model can deal with wetting and drying over complex geometry while obtaining the well balanced property at wet-dry boundaries. The positivity preserving property of the proposed model is proved. It is noted that, due to the existence of ambient water, there is no true dry bed involved. The wetting-drying process mentioned in the current paper refers in particular to the wetting and drying process of turbidity current in the ambient water. This paper is organized as follows. In §2, the governing equations and closures are presented, for which in §3 the well-balanced and positivity preserving techniques for turbidity current system are developed. Several numerical tests are conducted and presented in §4. Some concluding remarks complete the study in §5. 2. Governing Equations

The governing equations for modeling the two-dimensional spreading of turbidity currents over erodible bed are conservation equations for fluid, sediment and bed material, respectively. A detailed derivation of the governing equations for 2-D turbidity currents can be found in [42]. In this study, the layer-averaged 2-D governing equations with nonequilibrium sediment transport adopted in [20] are used:

CE

66

AC

65

ht + (hu)x + (hv)y = ew |V | + (E − D)/(1 − p),

(1)

  u(E − D)(ρ0 − ρ) ρw − ρ g0 2 2 +(huv)y = −g 0 hZx −(1+rw )u2∗ − − uew |V |, (hu)t + hu + h 2 ρ(1 − p) ρ x (2)   0 g v(E − D)(ρ − ρ) ρ − ρ 0 w (hv)t + (huv)x + hv 2 + h2 = −g 0 hZy − (1 + rw )v∗2 − − vew |V |, 2 ρ(1 − p) ρ y (3) 4

ACCEPTED MANUSCRIPT

100

(hc)t + (huc)x + (hvc)y = E − D.

103 104 105 106 107 108 109 110 111 112 113 114

The bed evolution is computed by using the Exner equation for the conservation of bed sediment, which takes the form: E−D . (5) Zt = − 1−p

In equations (1)–(5), t is time, x and y are horizontal coordinates, h = h(x, y, t) is the thickness of the turbidity current layer, u = u(x, y, t) and v = v(x, y, t) are the x- and ycomponents of the layer-averaged velocities, respectively, c = c(x, y, t) is the layer-averaged √ volumetric sediment concentration, Z = Z(x, y, t) is the bed elevation, |V | = u2 + v 2 is the total velocity, g is the gravitational acceleration, g 0 = Rgc is the submerged gravitational acceleration, R = (ρs − ρw )/ρ is submerged specific gravity, ρw and ρs are the densities of water and sediment, respectively, ρ = ρw (1 − c) + ρs c is the density of the water-sediment mixture, ρ0 = ρw p + ρs (1 − p) is the density of the saturated bed, p p p is the bed porosity; u∗ = CD u|V | and v∗ = CD v|V | are the components of the friction source in the x- and y-directions, respectively, CD is the bed drag coefficient with a typical range 0.002-0.06 [43], rw is the ratio of upper-interface resistance to bed resistance with rw = 0.43 [42], ew is a fluid entrainment coefficient which is estimated by [42]:

CR IP T

102

115

0.00153 , 0.0204 + Ri

where Ri is Richardson number defined by:

119 120 121

122 123

PT

CE

118

(7)

Furthermore, in (1)–(5), E = ωEs is the sediment entrainment, D = ωcb is the sediment deposition, ω is the settling velocity of sediment calculated as [47] p (8) ω = (13.95ν/d)2 + 1.09(ρs /ρw − 1)gd − 13.95ν/d,

in which ν is the kinematic viscosity of water and d is the mean diameter of sediment particles, cb = rb c is the local near-bed sediment concentration, and rb is a coefficient greater than 1. Es is the near-bed concentration at capacity condition and can be estimated by the empirical formula proposed by Parker et al. [43] which is modified by Hu et al. [20]:

AC

117

ED

Rghc Ri = √ . u2 + v 2

116

(6)

M

ew =

AN US

101

(4)

5 1.3 × 10−7 km , (9) 5 1 + 4.3 × 10−7 km p = Rgd3 /ν is the particle Reynolds number and ψp

Es = ψp

p where km = CD Rep |V |/ω, and Rep is a correction coefficient.

5

ACCEPTED MANUSCRIPT

124

125 126

3. Numerical scheme Following the coupling approach developed by Fagherazzi and Sun [10], see also [32], the 2-D bed evolution equation (5) is rewritten as: (β)t + (huc)x + (hvc)y = 0.

129 130 131 132 133 134

CR IP T

128

with introducing a new variable β = (1−p)Z +hc, which is the sediment volume in a vertical column above a reference level that can be separated in the sediment volume deposited on the bed (1 − p)Z; and the sediment volume suspended in the fluid. In order to design a well-balanced numerical scheme for the governing system, the water surface level denoted by η := h + Z is introduced as a primary variable instead of h in the current study. The “well-balanced” property of the current model is defined in section §3.4. A coupling strategy is implemented for solving the turbidity currents over erodible bed in the present study; the coupled system (1)–(10) can be rewritten in the following vector form:

AN US

127

(10)

Ut + Fx + Gy = S + K.

where the variables U , and the fluxes F and G are:      hu hv η 2    (hu) g   (hu)(hv)        η − Z + 2 R(η − Z)(hc)   hu  η−Z         (hv)2 g   (hu)(hv)        η − Z + 2 R(η − Z)(hc)    η−Z U =  hv  , F =  , G =       (hc)(hv)    (hu)(hc)       hc     η−Z η−Z           (hc)(hv) (hu)(hc) β η−Z η−Z



       ,      

(12)

CE

and the source terms S and K are:

AC

136

PT

ED

M

135

(11)



 0    −gR(hc)Zx     S=  −gR(hc)Zy  ,     0 0

6

(13)

ACCEPTED MANUSCRIPT

139 140 141 142 143

√ e w u2 + v 2

   −(1 + r )C u√u2 + v 2 − u(E − D)(ρ0 − ρ) − ρw − ρ e u√u2 + v 2 w w D  ρ(1 − p) ρ   √ √ K=  −(1 + rw )CD v u2 + v 2 − v(E − D)(ρ0 − ρ) − ρw − ρ ew v u2 + v 2  ρ(1 − p) ρ    E−D  0



      .      

(14)

CR IP T

138



In this S study, the computational domain is discretized by an unstructured triangulation T := i Ti , consisting of triangular cells Ti of size |Ti |. The outer unit normal to the corresponding sides of Ti of length `ik , k = 1, 2, 3, is denote by nik := (cos(θik ), sin(θik))T . The coordinates of the center are denote by (xi , yi ), and the midpoint of the k-th side of the triangle Ti , k = 1, 2, 3 is denoted by Mik = (xik , yik ). The vertices of Ti are denoted by Viκ = (xiκ , yiκ ), κ = 12, 23, 31. Ti1 , Ti2 and Ti3 are the neighboring triangles that share a common side with Ti , see Figure 1.

i i

Vi31

Vi23

CE

PT

ED

M

Vi12

AN US

137

144

145 146 147 148 149 150

AC

Figure 1: A typical triangular cell with three neighbors.

3.1. Evaluation of numerical fluxes using the central-upwind scheme In order to solve the system (11), the flux terms: Fx +Gy , have to be computed by a stable and accurate solver, for which the central-upwind scheme is chosen in this study. The central-upwind schemes, which are Riemann problem free shock capturing schemes, enjoy the main advantages of the Godunov-type central schemes, i.e., simplicity, universality and robustness and can be applied to problems with complicated geometries. The readers are 7

ACCEPTED MANUSCRIPT

151 152 153 154

referred to [31, 28, 29, 30, 27, 5, 34, 36, 33, 35], where the central-upwind schemes are developed and extended. Integrating both sides of system (11) over the control volume Ti and applying the Green’s formula: Z Z div G dxdy = G · n ds, (15) Ti

one can get a semi-discrete scheme which reads:

CR IP T

155

∂Ti

3

1 X dU i =− `ik Hik + S i + K i . dt |Ti | k=1 156 157

where Hik is the normal flux through the side `ik . Using the central-upwind scheme, Hik is estimated by: out out ain ain ik F (Uik (Mik )) + aik F (Ui (Mik )) ik G(Uik (Mik )) + aik G(Ui (Mik )) + sin(θ ) ik out out ain ain ik + aik ik + aik  ain aout  − inik ikout Uik (Mik ) − Ui (Mjk ) . aik + aik (17) where Ui (Mik ) and Uik (Mik ) are the reconstructed values at Mik using a piecewise linear reconstruction:

160

Ui (Mik ) :=

164 165

166 167

Uik (Mik ) :=

min(U i ,U ik ) ≤ Ui (Mik ) ≤ max(U i ,U ik ),

˚ (x, y), U

(19)

k = 1, 2, 3

(20)

b x )i = (U b y )i = 0 is used in the reconstruction (18). Furthermore, is not satisfied for cell Ti , (U Z is also reconstructed over Ti using (18), and β is reconstructed by: ˚ y) := (1 − p)Z(x, ˚ ˚ y) + hc(x, β(x, y),

(x, y) ∈ Ti .

(21)

For the dry cells in which η i = Z i , the reconstructions of η, hu, hv and hc are recalculated

168 169

lim

(x,y)→Mik ;(x,y)∈Tik

(18)

PT

163

˚ (x, y), U

CE

162

lim

(x,y)→Mik ;(x,y)∈Ti

(x, y) ∈ Ti

b x )i and (U b y )i are component-wise approximations of the in which the limited gradients (U gradients Ux (xi , yi , t) and Uy (xi , yi , t), respectively, computed using a nonlinear slope limiter to minimize the oscillations of reconstruction (18), see §3.2 for the computations of Ux , Uy , b x )i and (U b y )i . In order to maintain numerical stability when reconstruction (18) is applied, (U the appearance of local extrema is checked. If the condition:

AC

161

M

of U , that is:

b x )i (x − xi ) + (U b y )i (y − yi ), ˚ (x, y) := U i + (U U

ED

159

AN US

Hik = cos(θik )

158

(16)

as: ˚ y), ˚ η (x, y) = Z(x,

˚ hu(x, y) = 0,

˚ hv(x, y) = 0, 8

˚ hc(x, y) = 0,

(x, y) ∈ Ti

(22)

ACCEPTED MANUSCRIPT

171 172 173 174 175 176

˚ y) is computed by the piecewise linear reconstruction (18), and (21) is recomwhere Z(x, puted based on (22). Cell averages of the source terms S i and K i in (16) need to be discretized in an appropriate form. In this study, special discretization of the source terms are proposed in order to guarantee the well-balanced property of the proposed scheme, as described in §3.4. out In (17), ain ik and aik are one-sided local speeds. In order to avoid division by zero, the in definitions of aik and aout ik in [5] are modified in this study. They are defined by:

CR IP T

170

ain ik = − min{λ1 [Jik (Ui (Mik ))], λ1 [Jik (Uik (Mik )], −β}, aout ik = max{λ5 [Jik (Ui (Mik ))], λ5 [Jik (Uik (Mik )], β},

178

The analytical expressions of the eigenvalues of (24) are:

λ1 =

181 182 183 184

are the eigenvalues of

0

(24)



 0     0 ,   0    0  0  0     0 .   0    0

(25)

(hu) cos(θik ) + (hv) sin(θik ) p (hu) cos(θik ) + (hv) sin(θik ) − gR(hc), λ2 = λ3 = , η−Z η−Z (hu) cos(θik ) + (hv) sin(θik ) p λ4 = 0, λ5 = + gR(hc). η−Z (26) Moreover, in order to correctly predict “lake at rest” steady states (see the descriptions in §3.4), the fourth component of (17) is computed using a centered scheme at interface `ik when ηi (Mik ) = ηik (Mik ), ci = cik and ui (Mik ) = uik (Mik ) = vi (Mik ) = vik (Mik ) = 0, where the velocities u and v at Mik are given by (33).

AC

180

CE

PT

ED

M

179

in which β=10−8 is a small positive number, λ1 [Jik ] ≤ · · · ≤ λ5 [Jik ] the Jacobian matrix: ∂G ∂F + sin(θik ) , Jik = cos(θik ) ∂U ∂U where  0 1 0 0  2 g  − (hu) + g R(hc) 2(hu) 0 R(η − Z)  (η−Z)2 2 η−Z 2   ∂F hv hu − (hu)(hv) 0 = (η−Z)2 η−Z η−Z  ∂U  hc hu  − (hu)(hc) 0  (η−Z)2 η−Z η−Z  hc hu − (hu)(hc) 0 (η−Z)2 η−Z η−Z  0 0 1 0  hv hu  − (hu)(hv) 0  (η−Z)2 η−Z η−Z  (hv)2 2(hv) ∂G  g g − (η−Z) 0 R(η − Z) 2 + 2 R(hc) = η−Z 2  ∂U  hc hv  − (hv)(hc) 0  (η−Z)2 η−Z η−Z  hc hv − (hv)(hc) 0 (η−Z)2 η−Z η−Z

AN US

177

(23)

9

ACCEPTED MANUSCRIPT

187 188 189

190

191 192 193 194

A fully discrete scheme can be obtained from the semi-discretization (16) by integrating the time-dependent terms using a stable and sufficiently accurate ODE solver. In the current study, the third-order strong stability preserving (SSP) Runge-Kutta solver, see, e.g., [14, 15], is used. To maintain the numerical stability, the time step size should satisfy the CFL condition, see [29], which can be expressed as:   rik 1 , (27) ∆t < min out 3 i,k max(ain ik , aik )

CR IP T

186

where rik , k = 1, 2, 3 are the three corresponding altitudes of triangle Ti ∈ T .

3.2. Estimations of the Gradients   The gradients (∇U )i = (Ux )i , (Uy )i of cell Ti are estimated by taking the x- and yderivatives of the planes passing through the points (xi1 , yi1 , Ui1 ), (xi2 , yi2 , Ui2 ) and (xi3 , yi3 , Ui3 ) (the gray area in Figure 2): (Ux )i =

(yi3 − yi1 )(Ui2 − Ui1 ) − (yi2 − yi1 )(Ui3 − Ui1 ) , (yi3 − yi1 )(xi2 − xi1 ) − (yi2 − yi1 )(xi3 − xi1 )

AN US

185

M

(xi2 − xi1 )(Ui3 − Ui1 ) − (xi3 − yi1 )(Ui2 − Ui1 ) . (Uy )i = (xi2 − xi1 )(yi3 − yi1 ) − (xi3 − xi1 )(yi2 − yi1 )

(28)

(xi2, yi2,Ui2)

(xi1, yi1,Ui1)

195 196 197 198 199 200

AC

CE

PT

ED

(xi3, yi3,Ui3)

Figure 2: Stencil for calculating nonlimited derivatives.

The use of the gradients estimated by (28) in reconstruction (18) can lead to oscillations and instabilities. In order to minimize the numerical oscillations, the gradients used in the piecewise linear reconstruction (18) have to be suppressed by a slope limiter. In the present study, the multidimensional slope limiter developed by Jawahar and Kab x )i and (U b y )i : math [25] is used to obtain the limited gradients (U b x )i = Λi1 (Ux )i1 + Λi2 (Ux )i2 + Λi3 (Ux )i3 , (U b y )i = Λi1 (Uy )i1 + Λi2 (Uy )i2 + Λi3 (Uy )i3 , (U 10

(29)

ACCEPTED MANUSCRIPT

where Λi1 , Λi2 and Λi3 are weights which are defined by: Λi1 = Λi2 = Λi3 =

207 208 209 210 211 212 213 214 215 216 217 218 219 220

3.3. Nonnegative Piecewise Linear Reconstruction for Z Notice that the use of the reconstruction (18) for η and Z cannot guarantee that the reconstructed values of the water surface level ηi (Mik ) are above the reconstructed bed level Zi (Mik ). This introduces unphysical negative water depths at Mik . Therefore, the correction of the piecewise linear reconstruction for Z or η is necessary to ensure that hi (Mik ) := ηi (Mik ) − Zi (Mik ) ≥ 0, k = 1, 2, 3. In [5], a 2-D correction algorithm is applied to η. However, such correction of η may fail in maintaining the “lake at rest” steady states at wetting and drying fronts for the studied system. Thus, the correction algorithm from [5] is modified in the present paper to reconstruct Z. Notice that such corrections should be taken only for those problematic cells. ˚ iκ , yiκ ) ≥ ˚ For the wet cells in which Z i < η i , it is impossible to have Z(x η (xiκ , yiκ ) for all three vertices of cell Ti . Thus, only two problematic cases should be considered: ˚i (xiκ , yiκ ) > • Case 1: There is only one vertex, for example, with index κ = 12, for which Z ˚ ηi (xiκ , yiκ ). ˚i (xiκ , yiκ ) = ˚ In this case, Z ηi (xiκ , yiκ ) is first set for the problematic vertex with κ = 12. ˚i (xiκ , yiκ ) = ˚ Then, basd on the conservation requirement, Z ηi (xiκ , yiκ ) − 32 (η i − Z i ) are set for other two vertexes with κ = 23, 31. The original piecewise linear reconstruction (18) for Z is replaced with a new one defined by: x−x ˚ y) − Z i y − yi Z(x, i = 0, (x, y) ∈ Ti . xi12 − xi yi12 − yi (31) ˚ η (x , y ) − Z i i12 i12 i xi23 − xi yi23 − yi ˚ ηi (xi23 , yi23 ) − 32 (η i − Z i ) − Z i

222 223 224 225 226

AC

CE

221

where the parameter ε is a small positive number introduced to prevent division by zero. In this paper, ε = 10−14 is used.

AN US

206

k(∇U )i1 k22 k(∇U )i2 k22 + ε , k(∇U )i1 k42 + k(∇U )i2 k42 + k(∇U )i3 k42 + 3ε

M

205

(30)

ED

204

k(∇U )i1 k22 k(∇U )i3 k22 + ε , k(∇U )i1 k42 + k(∇U )i2 k42 + k(∇U )i3 k42 + 3ε

PT

202 203

k(∇U )i2 k22 k(∇U )i3 k22 + ε , k(∇U )i1 k42 + k(∇U )i2 k42 + k(∇U )i3 k42 + 3ε

CR IP T

201

˚i (xiκ , yiκ ) > • Case 2: There are two vertices, for example, with indexes κ = 12, 23, for which Z ˚ ηi (xiκ , yiκ ). ˚i (xiκ , yiκ ) = ˚ In this case, Z ηi (xiκ , yiκ ) is first set for these two problematic vertices with ˚i (xiκ , yiκ ) = ˚ κ = 12, 23. Then, basd on the conservation requirement, Z ηi (xiκ , yiκ )−3(η i −Z i ) is set for another vertex with κ = 31. The original piecewise linear reconstruction (18) for 11

ACCEPTED MANUSCRIPT

229 230 231 232

233 234 235 236 237

= 0,

(x, y) ∈ Ti .

The above described correction yields a new linear reconstruction for Z which is con˚i (x, y) is not greater than the reconstructed servative, and guarantees that the corrected Z ˚ ηi (x, y) over the whole computational domain. In order to avoid the division by very small water depth h := η − Z, as proposed in [30], the desingularized point values of u, v and c are computed by (the i, k indexes are omitted): √ √ √ 2 h (hu) 2 h (hv) 2 h (hc) u= p , v=p , c= p , (33) 4 4 4 4 4 h + max(h , δ) h + max(h , δ) h + max(h4 , δ) where δ is a prescribed tolerance; δ = 10−16 is used in the present study. Notice that the choice of the value of δ is relevant with the precision of computer number format. In the current study, the computer program for testing the proposed model is double precision. Using the desingularized point values computed in (33), the x- and y-discharges and fluxes are recalculated by:

(34)

ED

M

(hu) := h · u, (hv) := h · v, (hc) := h · c,  T g F = (hu), (hu) · u + R(η − Z)(hc), (hv) · u, (hc) · u, (hc) · u , 2 T  g G = (hv), (hu) · v, (Dv) · v + R(η − Z)(hc), (hc) · v, (hc) · v , 2

and the one-sided local speeds (23) are computed using the following formulas: p p ⊥ ain gRi (Mik )(hc)i (Mik ), u⊥ gRik (Mik )(hc)ik (Mik ), 0}, ik = − min{ui (Mik ) − ik (Mik ) − p p ⊥ aout gRi (Mik )(hc)i (Mik ), u⊥ gRik (Mik )(hc)ik (Mik ), 0}, ik = max{ui (Mik ) + ik (Mik ) + (35) where Ri (Mik ) and Rik (Mik ) are computed using ci (Mik ) and cik (Mik ) estimated by (33), ⊥ respectively; the normal velocities denoted by u⊥ i (Mik ) and uik (Mik ) at the point Mik are given by: u⊥ i (Mik ) := ui (Mik ) cos(θik ) + vi (Mik ) sin(θik ), (36) u⊥ ik (Mik ) := uik (Mik ) cos(θik ) + vik (Mik ) sin(θik ).

241

242 243 244 245 246

AC

239 240

CE

PT

238

(32)

CR IP T

228

Z is replaced with a new one defined by: x−x ˚ y) − Z i y − yi Z(x, i xi12 − xi yi12 − yi ˚ ηi (xi12 , yi12 ) − Z i xi23 − xi yi23 − yi ˚ ηi (xi23 , yi23 ) − Z i

AN US

227

3.4. Discretization of the Source Terms A desired property of discretization of bed-slope source term S i is well-balanced property, which guarantees that the proposed discretization of S i exactly balances the numerical fluxes under the “lake at rest” steady states. The “lake at rest” steady solutions should satisfy that: η ≡ ηc = Const, c ≡ cc = Const, u ≡ v ≡ 0. (37) 12

ACCEPTED MANUSCRIPT

250

251

 3 g X 1 Rik (Mik )(hc)ik (Mik ) + Ri (Mik )(hc)i (Mik ) − ηc `ik cos(θik ) |Ti | k=1 2 2  Rik (Mik )(hc)ik (Mik )Zik (Mjk ) + Ri (Mik )(hc)i (Mik )Zi (Mjk ) (2) + Si = 0 − 2

(38)

 3 g X 1 Rik (Mik )(hc)ik (Mik ) + Ri (Mik )(hc)i (Mik ) − `ik sin(θik ) ηc |Ti | k=1 2 2  Rik (Mik )(hc)ik (Mik )Zik (Mjk ) + Ri (Mik )(hc)i (Mik )Zi (Mjk ) (3) − + Si = 0 2

(39)

CR IP T

249

and

in which

(2) Si

1 2 Rh c cos(θik ) ds = 2

Z

Ti

(∂Ti )k

g ≈− |Ti |

Z

RhcZy dxdy.

(40)

Ti

Rhchx dxdy +

Z

1 2 Rh cx dxdy + 2

Ti

Z

1 2 h cRx dxdy. 2

(41)

Ti

Based on the fact that h = η − Z and hx = ηx − Zx , (41) can be further reformulated as: −

Z

RhcZx dxdy =

CE

Ti

Z 3 X 1 k=1

−

Z

Ti

2

(Rhcη − RhcZ) cos(θik ) ds −

(∂Ti )k

1 R(η − Z)2 cx dxdy − 2

Z

Z

Rhcηx dxdy

Ti

1 (η − Z)hcRx dxdy. 2

(42)

Ti

where (∂Ti )k is the k-th side of the triangle Ti , k = 1, 2, 3. Next, one can apply the midpoint (2) rule to the integrals on the RHS of (42) and obtain the following discretization of S i :

AC

256

Ti

ED

k=1

255

RhcZx dxdy,

(3) Si

Following the procedures given in [5], the Green’s formula (15) is first applied to the vector field G = ( 12 Rh2 c, 0), and one can obtain: Z 3 X

254

Z

PT

252 253

g ≈− |Ti |

AN US

248

out In the “lake at rest” states (37), one can get ain ik = aik by (35). After substituting the “lake at rest” states into (16), one can conclude that the well-balanced discretization of S i should satisfy the following conditions:

M

247

(2) Si

 3 g X 1 Rik (Mik )(hc)ik (Mik ) + Ri (Mik )(hc)i (Mik ) = `ik cos(θik ) ηi (Mik ) |Ti | k=1 2 2  Rik (Mik )(hc)ik (Mik )Zik (Mjk ) + Ri (Mik )(hc)i (Mik )Zi (Mjk ) − 2     1 1 − gR hc (ηx )i − gR (η i − Z i )2 (cx )i − g (η i − Z i ) hc (Rx )i . 2 2 i i 13

(43)

ACCEPTED MANUSCRIPT

(3)

Similarly, one can obtain the discretization of S i : (3) Si

258 259 260 261

 3 g X 1 Rik (Mik )(hc)ik (Mik ) + Ri (Mik )(hc)i (Mik ) ηi (Mik ) = `ik sin(θik ) |Ti | k=1 2 2  Rik (Mik )(hc)ik (Mik )Zik (Mjk ) + Ri (Mik )(hc)i (Mik )Zi (Mjk ) − 2     1 1 − gR hc (ηy )i − gR (η i − Z i )2 (cy )i − g (η i − Z i ) hc (Ry )i . 2 2 i i

The gradients ηx , ηy , cx , cy , Rx and Ry are calculated using (28). Since ηx ≡ ηy ≡ 0 for η ≡ constant, cx ≡ cy ≡ 0 and Rx ≡ Ry ≡ 0 for c ≡ constant in the “lake at rest” states, the quadratures (43)–(44) satisfy (38)–(39), respectively. The source term K i is discretized straightforwardly using the midpoint rule:

268

269 270 271 272

273 274 275

276 277 278 279 280 281

M

267

3.5. Positivity preserving property of the proposed scheme Based on the fact that the maximum amount of deposition in one time step should not be greater than the amount of suspended sediment in the fluid, in the present study, the following condition is specified: (45) ∆tD ≤ σ(hc)i ,

ED

266

PT

265

where σ is a coefficient specified in the range 0 < σ ≤ 1, and σ = 0.5 is used in this study. Based on the empirical relation D = ωrb ci given in §2, using that (hc)i = hi ci , (45) can be reformulated as: ! 1 hi ∆t ≤ min . (46) i 2 ωrb

CE

264

and it does not affect the well-balanced property of the developed numerical scheme. Notice that, for the turbidity flows driven by suspended sediment over irregular beds, the “lake at rest” steady status given by (37) can not be preserved in the real-world applications due to the fact that sediment deposition will cause the density differences and horizontal advection. Such steady status is only an ideal condition defined to facilitate the developments and tests of “well-balanced” property. In order to preserve such steady status, one have to force E=D=0, see test 4.1.

AC

263

AN US

K i = K(xi , yi ), 262

(44)

CR IP T

257

It can be noted that (46) may lead to very small ∆t when water depth hi is close to zero in any one cell of the whole computational domain. To avoid such inefficiency in practical applications, a threshold hi = 1 × 10−6 is set in the current study, which implies that (46) is only used for those cells in which the water depth is no lesser than the specified threshold. Notice that this specified threshold is only used for (46) to ensure the positivity preserving property; it does not affect dealing with wetting and drying by proposed model. 14

ACCEPTED MANUSCRIPT

282 283

Using the forward Euler temporal discretization, and subtracting the fourth and fifth components from the first component of (16), one can obtain: n+1

hi

n

= hi −

3  ∆t X `ik cos(θik )  in out (hu) (M ) (hu) (M ) + a a j jk jk jk ik ik out |Ti | k=1 ain ik + aik

(47)

CR IP T

3  ∆t X `ik sin(θik )  in out − a (hv) (M ) + a (hv) (M ) jk jk j jk ik ik out |Ti | k=1 ain ik + aik

3 √  ∆t X ain aout  E−D + , `ik inik ikout hik (Mik ) − hi (Mjk ) + ∆tew u2 + v 2 + ∆t |Ti | k=1 aik + aik 1−p

285 286

where the subscripts n and n + 1 represent the time t = tn and t = tn+1 , respectively. 3 1P n hi (Mjk ), (47) can be reformulated Based on (34) and (36), and the fact that hi = 3 k=1 as: n+1 hi

AN US

284

3 √  ∆t X `ik ain  E ⊥ 2 2 = hik (Mik ) in ikout aout − u (M ) + ∆te ik w u + v + ∆t ik ik |Ti | k=1 aik + aik 1−p

n1

289

3 X

hi (Mik )

 in o ∆tD ∆t `ik aout ik ⊥ a + u (M ) − ≥ ik ik i out |Ti | ain 1−p ik + aik

n1  in  1 ci o ∆t `ik aout ik ⊥ hi (Mik ) − aik + ui (Mik ) − ≥ 0. out 3 |Ti | ain 61−p ik + aik k=1

3 X

 in  1 ci ∆t 1 ci 1 ∆t `ik aout ik aik + u⊥ ≤ `ik aout ≤ . i (Mik ) + ik + in out |Ti | aik + aik 61−p |Ti | 61−p 3

291 292

(49)

out Using u⊥ ik (Mik ) ≤ aik from (35), condition (49) gives:

AC

290

3

−

CE

k=1

n1

ED

288

⊥ One can get aout ik ≥ uik (Mik ) from (35) and therefore, the first three terms on the RHS 3 1P n hi (Mjk ) and (45), one can obtain the following of (48) are nonnegative. Using hi = 3 k=1 n+1 condition to guarantee that hi is nonnegative:

PT

287

(48)

M

 in o ∆tD ∆t `ik aout ik ⊥ + hi (Mik ) − . a + u − (M ) ik ik i in out 3 |T 1 − p i | aik + aik k=1 3 X

(50)

Based on (50) and the fact that |Ti | = 0.5 rik lik , the following limitation has to be satisfied n+1 to ensure that hi computed by (47) is nonnegative, that is: " #  1 rik 0.5ci  ∆t ≤ min 1 − . (51) out 6 i,k max(ain 1−p ik , aik ) 15

ACCEPTED MANUSCRIPT

294 295 296 297

Another desired property of the proposed numerical model is that the positivity preservn+1 ing property for sediment concentration cn+1 . Based on the fact that hi ≥ 0 guaranteed i n+1 by (51), in order to prove that cn+1 ≥ 0, one only needs to ensure that (hc)i ≥ 0. Using i the forward Euler temporal discretization, the fourth component in equation (16) can be written in a fully discrete form: n+1

(hc)i

3  ∆t X `ik cos(θik )  in aik (hu)jk (Mjk )cjk (Mjk ) + aout = (hc)i − ik (hu)j (Mjk )cj (Mjk ) out in |Ti | k=1 aik + aik n

CR IP T

293

3  ∆t X `ik sin(θik )  in aik (hv)jk (Mjk )cjk (Mjk ) + aout − ik (hv)j (Mjk )cj (Mjk ) out in |Ti | k=1 aik + aik

+

3  ∆t X ain aout  `ik inik ikout hik (Mik )cik (Mik ) − hi (Mjk )ci (Mjk ) + ∆t(E − D). |Ti | k=1 aik + aik

AN US

(52) 3 1P Based on (34) and (36), and the fact that (hc)i = (hc)i (Mjk ), (52) can be reformulated 3 k=1 as: 3  n+1 ∆t X ain hik (Mik )cik (Mik )  out `ik ik (hc)i = aik − u⊥ ik (Mjk ) + ∆tE in out |Ti | k=1 aik + aik (53) 3 i h1 X
`ik aout ⊥ − ∆tD. hi (Mjk )ci (Mjk ) − in ik out ain ik + ui (Mik ) 3 a + a ik ik k=1

299

300 301 302

304

305 306

AC

CE

PT

303

⊥ One can get aout ik ≥ uik (Mik ) from (35) and therefore, the first two terms on the RHS of (53) are nonnegative. Accordingly, the sum of the last two terms on the RHS of (53) n+1 ≥ 0. Using (34) and (45), one can get needs to be nonnegative to guarantee that (hc)i 3 n 1P (hc)i = hi (Mjk )ci (Mjk ), and this gives the following condition: 3 k=1 3 h1 X
i `ik aout ⊥ hi (Mjk )ci (Mjk ) − in ik out ain + u (M ) − ∆tD ≥ ik ik i 3 a + a ik ik k=1 (54) 3 h1 i X
`ik aout 1 ⊥ hi (Mjk )ci (Mjk ) − in ik out ain ≥ 0. ik + ui (Mik ) − 3 a + a 6 ik ik k=1

ED

298

M

n

out Using u⊥ ik (Mik ) ≤ aik from (35), condition (54) gives:   in ∆t ∆t `ik aout 1 ik ⊥ + u (M ) ≤ a `ik aout . ik ik i ik ≤ in out |Ti | aik + aik |Ti | 6

(55)

Based on (55) and the fact that |Ti | = 0.5 rik lik , the following limitation has to be satisfied n+1 to ensure that (hc)i computed by (52) is nonnegative, that is: " # 1 rik ∆t ≤ min . (56) out 12 i,k max(ain ik , aik ) 16

ACCEPTED MANUSCRIPT

307 308 309 310 311 312

ci
< 1, conditions (27), (51) and (56) can be Based on the fact that 0.5 ≤ 1 − 0.5 1−p combined into one single limitation which is (56). Notice that (56) is a stricter CFL condition than (27) and can ensure the stability of the numerical method. Finally, to guarantee the positivity preserving property of the proposed numerical scheme while maintaining the numerical stability, ∆t should satisfy the restrictions (46) and (56) at the same time for each time step.

330

4. Numerical Experiments

320 321 322 323 324 325 326 327 328

331 332

333 334 335 336 337 338

339 340

AN US

319

M

318

ED

317

In this section, we demonstrate the performance of the developed numerical model on four test problems.

PT

316

4.1. “Lake at rest” steady status over irregular beds with wet-dry interfaces This test case is designed to demonstrate the well-balanced property of the proposed scheme, which is guaranteed by the proposed discretization of the bed slope source terms (43)–(44). The well-balanced property is evaluated not only in the wet domain but also at the boundaries between wet and dry zones. To this end, two islands with different altitudes are defined as  Z(x, y) = max 0 m, 1 m − 0.8 · [(x − 4 m)2 + (y − 2 m)2 ], (57) 0.5 m − 1.2 · [(x − 2 m)2 + (y − 2 m)2 ]

CE

315

AC

314

CR IP T

329

3.6. “Lake at rest” steady solutions at wetting and drying fronts The nonnegative piecewise linear reconstruction and desingularization described in §3.3 and positivity preserving property described in §3.5 can successfully deal with the wetting and drying, and it does not require a minimum fluid depth and concentration to be set. The “well-balanced” discretization of bed slope source term described in §3.4 can successfully guarantee that the “Lake at rest” steady solutions (37) are satisfied for wet domain. In order to further guarantee the “lake at rest” steady states at wetting-drying fronts, the following treatments are used. When calculating the numerical fluxes using (17) for interfaces `ik between a quiescent wet cell and a dry cell, the gradients of variables at such interfaces are temporarily zeroed when η i of the wet cell is no higher than the Z i of the neighboring dry cell, using a zero-order extrapolation condition. This treatment is important to eliminate the unphysical waves and oscillations at wetting-drying fronts under the “lake at rest” status (37). Moreover, using the weighted slope limiter (29) can not guarantee that the “lake at rest” solution (37) is satisfied at wetting and drying fronts due to the fact that the gradient ηx and ηy of the neighboring dry cell is not zero. To avoid this, (b ηx )i = (b ηy )i = 0 is used in (18) for the quiescent wet cells which have at least one dry neighboring cell.

313

in a [0 m, 2 m]×[0 m, 1 m] computational domain, and sediment entrainment and deposition are not considered in this test. 17

ACCEPTED MANUSCRIPT

341 342 343 344 345

AN US

CR IP T

346

A quiescent lake around the two islands is defined with an initial water surface elevation of 0.6 m, so that the test can involve the wet-dry interfaces for the big island with a height of 1 m. The initial sediment concentration in the water is set to 0.01, the mean diameter of the sediment is set to 0.1 mm and the porosity of the saturated bed is 0.45. The domain is divided into 9600 triangular cells. Figure 8 shows the 3-D views of the simulated still water surface at t = 50 s using the proposed well-balanced model. An undisturbed water surface is observed throughout the entire simulation.

Figure 3: Simulated quiescent water surface over irregular bottom with wetting-drying boundaries at t = 50 s using the proposed well-balanced model.

351 352

(2) Si

354 355 356 357 358 359 360 361 362

3

g X Zi (Mik ) + Zik (Mik ) = −R(xi , yi )hc(xi , yi ) `ik sin(θik ) . |Ti | k=1 2

(58)

Using the same grid, the test runs 50 s for the well-balanced model however can only run 1.03s for the non-well-balanced model due to the large oscillations generated at the wet-dry boundaries. Figures 9 and 10 show the simulated discharges at y=2 m cross sections and simulated volumetric concentrations at x=2 m cross sections using the well-balanced and the non-well-balanced schemes, respectively. It can be clearly observed that spurious waves are generated by the non-well-balanced scheme at the plateau of the submerged island and the wet-dry boundaries as shown in Figure 9 (b), which leads to a wrong prediction of sediment concentration in these areas as illustrated in Figure 10 (b). This clearly demonstrates the importance and necessity of the well-balanced property of the current model guaranteed by the proposed scheme (43)–(44).

AC

353

3

g X Zi (Mik ) + Zik (Mik ) = −R(xi , yi )hc(xi , yi ) `ik cos(θik ) , |Ti | k=1 2

CE

(3) Si

ED

350

In order to show the necessity and importance of the proposed well-balanced discretization (43)–(44), a non-well-balanced version of the current proposed model is introduced and applied to the same numerical test for comparison purposes. The non-well-balanced version (2) (3) of the proposed model is obtained by replacing the bed-slope terms S i and S i in (38) and (39) by a straightforward midpoint rule discretization:

PT

349

M

347 348

18

CR IP T

ACCEPTED MANUSCRIPT

(a) Well-balanced model

(b) Non-well-balanced model

PT

ED

M

AN US

Figure 4: Simulated and exact discharges at y=2m cross sections

Figure 5: Simulated and exact volumetric concentrations at x=2m cross sections

365 366 367 368 369 370 371 372 373 374

CE

364

4.2. Sediment-driven turbidity current in a channel over initially dry beds This experimental test is used to verify the accuracy of the proposed numerical model. The laboratory tests, modeling the turbidity currents over submarine canyons and the formed deposition fans, were reported in [13]. The experiments were conducted in a 30 cm wide and 11.6 m long channel which consists of a 5 m long inclined bed with a slope of 0.08, followed by a 6.6 m long horizontal bed. The flume is filled by clear water as ambient fluid, and the initial bed is treated as a dry bed with h=0 m. A sediment-laden bottom current of known density was released at the top of the inclined bed at a known rate, and the current starts to flow downslope at t=0 s. The inlet current thickness h0 was set to be 3 cm with a certain layer-averaged inlet velocity u0 , see Figure 6. Notice that, since the initially dry bed is non-erodible, the updated bed level can not be less than the initial bed level. In the numerical simulations, a 11.6 m× 0.4 m computa-

AC

363

19

CR IP T

ACCEPTED MANUSCRIPT

Figure 6: Sketch of the laboratory flume.

375 376

tional domain is equally divided into 18560 triangular cells, and it is found that finer mesh leads to essentially same numerical results. The modelling parameters of several selected experimental runs in [13] are given in Table 1.

u0 (cm/s)

c0

d (mm)

CD

ω (cm/s)

p

ρs /ρw

ψp

Time (min)

GLASSA2

8.3

0.00339

0.03

0.01

0.085

0.5

2.5

0.15

30

GLASSA5

8.3

0.00394

0.03

0.01

0.085

0.5

2.5

0.15

33

GLASSA6

11

0.00449

0.03

0.01

0.085

0.5

2.5

0.15

30

GLASSA8

11

0.00579

0.03

0.01

0.085

0.5

2.5

0.15

28

382 383 384 385 386 387 388 389 390 391 392 393 394 395

ED

4.3. 2-D axisymmetric particle-driven turbidity currents Two-dimensional experiments of axisymmetric particle-driven current by Bonnecaze et al. [3] is numerically investigated using the proposed model. The experiments were performed to determine the radius of a fixed-volume gravity current at different times and its resulting deposition pattern. A plan view sketch of the laboratory tank is given in Figure 8. At t=0 s, the lock gate was suddenly lifted to release the turbid fluid stored in the upstream reservoir, i.e. the rectangular part with gray color in Figure 8. The downstream side of the lock gate was initially filled with clear water. The initial thickness of the well-mixed turbidity in the reservoir is 0.14 m. The suspension for the current was made by non-cohesive silicon carbide particles with a diameter of 37 µm, a density of ρs =3217 kg/m3 , a settling velocity of 0.17 cm/s and a porosity of 0.5. Three experimental runs with different initial sediment concentrations (0.0193, 0.00966 and 0.00506) were utilized. In the numerical simulations, the computational domain is discretized into 12174 unstructured triangular cells. The drag coefficient CD is set to 0.02, ψp is set to 0.15 and rb is set to 1 according to [4].

PT

381

CE

380

Figure 7 shows the comparisons between computed and measured results of deposition patterns against distance. Good agreements can be observed for selected runs. This shows that the proposed model can accurately simulate the sediment-driven turbidity.

AC

379

M

Runs

377 378

AN US

Table 1: Modelling parameters of selected experimental runs in [13]

20

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

ED

M

Figure 7: Computed and measured sediment deposit against distance.

397 398 399 400 401 402 403 404 405

Figure 9 (a) shows the computed bed increase over initial bed (gray area) at the end of the simulation when all the sediment deposit on the bed. As expected, it can be observed that much of the sediment settle out in the reservoir and the downstream area near the lock gate; the sediment-driven current front disappear at around x=1.77 m position because all the particles are settled out. Figure 9 (b) shows the vector of the current velocity at t=21 s; the distribution of the velocity is as expected, and no oscillation is observed at the wetting and drying front. In Figure 9 (b), it can be seen that the rectangular area on the downstream side of the gate contains the highest concentration, the concentration decreases gradually in the radial part due to the propagation and spreading. The left column of Figure 10 shows the comparisons between the computed and measured

AC

396

CE

Figure 8: Plan view sketch of the laboratory tank.

21

CR IP T

ACCEPTED MANUSCRIPT

Figure 9: Run with initial concentration c=0.00966: (a) Computed bed level increase when the test is done; (b) Simulated velocities and concentration of the turbidity current at t=21 s.

411 412 413 414 415 416 417

418 419 420 421 422 423 424

4.4. Submarine canyon-fan over irregular dry bed This numerical experiment is used to test the performance of the proposed model when solving problems involving wetting and drying over irregular topography and submarine canyon-fan transition, i.e. internal hydraulic jump. A submarine canyon followed by an initial dry plain is modeled over a 15 m×20 m domain, which consists of a 4 m long inclined bed with a slope of 0.08 in the x-direction, followed by a 11 m long plain. As shown in Figure 11, three humps are defined on this submarine plain by: n   Z(x, y) = max 0 m, 0.7 m − 0.6 · (x − 6 m)2 + (y − 10 m)2 ,    o 0.2 m − 0.1 · (x − 8 m)2 + (y − 5.5 m)2 , 0.2 m − 0.1 · (x − 8 m)2 + (y − 14.5 m)2 .

AC

425

AN US

410

M

409

ED

408

PT

407

sediment deposit as a function of radial distance for different initial concentrations. Very good agreement is observed, particularly in the radial part of the flume. The right column of Figure 10 shows the comparisons between the predicted and measured front positions as a function of time for different initial concentrations. A general good agreement between the computed and measured results can be observed. It can be seen that the model can correctly capture the trend of the front wave prorogation. Both computed and measured results show that the speed of the current front decreases as it propagates downstream due to the sediment deposition which reduces the driving force for the turbidity current; the farthest distance the current that can travel in the radial flume is reduced if the lower initial concentration is utilized. Moreover, in Figure 10, the simulated results of the proposed model are compared with the results in Bradford et al. [4], it shows that our proposed model yields more accurate results than the model developed by Bradford et al. in [4].

CE

406

426 427 428 429

The whole computational domain is initially dry and erodible. The inlet of turbidity flow is located in the center of the upstream sidewall as indicated in Figure 11. The simulated parameters of the inflow are: h=0.05 m, u=v=0 m/s, c=0.04. Other parameters are: p=0.5, d=0.05 mm, ρs =2500 kg/m3 , ρw =1000 kg/m3 , CD =0.004, ψp =1, rb =1.6. 22

AN US

CR IP T

ACCEPTED MANUSCRIPT

430 431 432 433 434 435 436

AC

CE

PT

ED

M

Figure 10: Simulated and measured sediment deposit (left column) and front position (right column) for fixedvolume releases of turbidity fluid.

Figure 11: Plan view of the computational domain.

Figure 12 shows the evolution of the turbidity current over the initially dry plain with humps at different times. It shows that the hump close to the inclined bed with higher altitude is washed over at the beginning and it then splits the upstream turbidity current. Subsequently, other two humps with lower altitude are washed over and flooded. Wetting and drying processes over irregular topography are existing during the whole simulation, and no oscillations are generated at wet-dry fronts. This can be further verified by Figure 13 which shows the vectors of current velocities at different times. No unphysical velocities are 23

ACCEPTED MANUSCRIPT

observed at the wetting and drying fronts during the whole simulation. This demonstrates the stability and robustness of the proposed numerical schemes.

AN US

CR IP T

437

PT

ED

M

Figure 12: Evolution of turbidity current over irregular dry bed (grey) at different times.

439 440 441 442 443 444 445 446 447 448 449

Richardson number (Ri) is an important parameter governing the behavior of turbidity currents [45]. A critical value 1 of Ri is generally used, so that the range Ri<1 corresponds to a high-velocity supercritical turbid flow regime, and the range Ri>1 corresponds to a low-velocity subcritical turbid flow regime. The change from supercritical flow to subcritical flow is accomplished via an internal hydraulic jump. Abrupt increase of water level near the inclined bed can be observed in Figure 12 which indicates the internal hydraulic jump. In order to further investigate the transition of the simulated turbidity flow, the Richardson number is calculated at y=9.5 m cross-section at different times, see Figure 14 (a). It shows sharp increases of Ri from less than 0.2 to more than 2 at around x=5 m position. This demonstrates the ability of simulating the super- and sub-critical turbidity flow and capturing internal hydraulic jump by proposed numerical model. Figure 14 (b) illustrates

AC

438

CE

Figure 13: Computed velocity vectors of turbidity current at different times.

24

ACCEPTED MANUSCRIPT

the bed changes at t=540 s. It can be seen that most of the suspended sediment settle out on the sloped bed near the inlet, which is in accord with expectation.

CR IP T

450

Figure 14: Richardson number (a) at y=9.5 m at different times, and simulated bed changes (b) at 540 s. 451

458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473

474

475 476

AN US

457

M

456

ED

455

A 2-D coupled numerical model is successfully developed to simulate sediment-driven turbidity current over irregular erodible bed with wetting and drying conditions. The proposed model is second order accurate in space and third order accurate in time. The model is built based on the complete mass and momentum conservation laws of the water-sediment mixture. After being verified by several numerical tests, the proposed model shows its robustness and satisfied accuracy. The developed numerical model incorporates two significant features: (i) One desired feature of the proposed models is well-balanced property which guarantees that the discretized bed-slope source terms exactly balance the numerical flux terms at steady status. In the developed numerical model, this property is successfully obtained by proposing a special discretization of bed-load source term. Furthermore, the special reconstruction of variables at wet-dry boundaries guarantee that the well-balanced property is also obtained at wetting and drying interfaces. The numerical tests clearly demonstrate the importance and necessity of the well-balanced property. (ii) Another desired feature of the proposed models is positivity preserving property which guarantees that the computed fluid depth and volumetric concentration of suspended sediment are non-negative in any time step during the simulation. In order to obtain this property, the technique of tracking the wetting and drying front is first developed using non-negative bottom reconstruction and desingularization of point values. Based on the central-upwind method and developed wetting/drying technique, the positivity preserving property of the proposed numerical model is mathematically proved under a specific CFL condition which guarantees the positivity preserving proprty while maintaining the numerical stability.

PT

454

CE

453

5. Conclusion

AC

452

6. Acknowledgement This research was supported by Natural Sciences and Engineering Research Council of Canada (NSERC).

25

ACCEPTED MANUSCRIPT

484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527

CR IP T

483

AN US

482

M

481

ED

480

[1] J. Alexander and T. Mulder. Experimental quasi-steady density currents. Mar. Geol., 186(3):195–210, 2002. [2] J. H. Baas, W. Van Kesteren, and G. Postma. Deposits of depletive high-density turbidity currents: a flume analogue of bed geometry, structure and texture. Sedimentology, 51(5):1053–1088, 2004. [3] R. T. Bonnecaze, M. A. Hallworth, H. E. Huppert, and J. R. Lister. Axisymmetric particle-driven gravity currents. J. Fluid Mech., 294:93–121, 1995. [4] S. F. Bradford and N. D. Katopodes. Hydrodynamics of turbid underflows. I: Formulation and numerical analysis. J. Hydraul. Eng., 125(10):1006–1015, 1999. [5] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova. Well-balanced positivity preserving centralupwind scheme on triangular grids for the Saint-Venant system. M2AN Math. Model. Numer. Anal., 45(3):423–446, 2011. [6] A. Cantelli, C. Pirmez, S. Johnson, and G. Parker. Morphodynamic and stratigraphic evolution of self-channelized subaqueous fans emplaced by turbidity currents. J. Sediment. Res., 81(3):233–247, 2011. [7] S. U. Choi. Layer-averaged modeling of two-dimensional turbidity currents with a dissipative-galerkin finite element method part i: Formulation and application example. J. Hydraul. Res., 36(3):339–362, 1998. [8] M. R. Chowdhury and F. Y. Testik. Laboratory testing of mathematical models for high-concentration fluid mud turbidity currents. Ocean eng., 38(1):256–270, 2011. [9] S. A. El-Gawad, C. Pirmez, A. Cantelli, D. Minisini, Z. Sylvester, and J. Imran. 3-d numerical simulation of turbidity currents in submarine canyons off the Niger Delta. Mar. Geol., 326:55–66, 2012. [10] S. Fagherazzi and T. Sun. Numerical simulations of transportational cyclic steps. Comput. & Geosci., 29(9):1143–1154, 2003. [11] J. J. Fedele and M. H. Garc´ıa. Laboratory experiments on the formation of subaqueous depositional gullies by turbidity currents. Mar. Geol., 258(1):48–59, 2009. [12] M. Garcia, W. Yu, and G. Parker. Experimental study of turbidity currents. Adv. Aerodyn., Fluid Mech., and Hydraul., pages 120–127, 1985. [13] M. H. Garc´ıa. Hydraulic jumps in sediment-driven bottom currents. J. Hydraul. Eng., 119(10):1094– 1117, 1993. [14] S. Gottlieb, D. Ketcheson, and C.-W. Shu. Strong stability preserving Runge-Kutta and multistep time discretizations. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2011. [15] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. SIAM Rev., 43:89–112, 2001. [16] M. A. Hallworth and H. E. Huppert. Abrupt transitions in high-concentration, particle-driven gravity currents. Phys. Fluids, 10(5):1083–1087, 1998. [17] A. E. Hay. Turbidity currents and submarine channel formation in rupert inlet, british columbia: 1. surge observations. J. Geophys. Res.: Oceans, 92(C3):2875–2881, 1987. [18] P. Hu and Z. Cao. Fully coupled mathematical modeling of turbidity currents over erodible bed. Adv. Water Res., 32(1):1–15, 2009. [19] P. Hu, Z. Cao, and G. Pender. Well-balanced two-dimensional coupled modelling of submarine turbidity currents. In Proc. Inst. Civ. Eng.-Maritime Eng., volume 165, pages 169–188. Thomas Telford Ltd, 2012. [20] P. Hu, Z. Cao, G. Pender, and G. Tan. Numerical modelling of turbidity currents in the Xiaolangdi reservoir, Yellow River, China. J. Hydrol., 464:41–53, 2012. [21] H. Huang, J. Imran, and C. Pirmez. Numerical model of turbidity currents with a deforming bottom boundary. J. Hydraul. Eng., 131(4):283–293, 2005. [22] J. Imran, G. Parker, and N. Katopodes. A numerical model of channel inception on submarine fans. J. Geophys. Res.: Oceans, 103(C1):1219–1238, 1998. [23] D. L. Inman, C. E. Nordstrom, and R. E. Flick. Currents in submarine canyons: An air-sea-land interaction. Annu. Rev. Fluid Mech., 8(1):275–310, 1976.

PT

479

CE

478

References

AC

477

26

ACCEPTED MANUSCRIPT

537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578

CR IP T

536

AN US

535

M

533 534

ED

531 532

PT

530

[24] M. Janocko, M. B. J. Cartigny, W. Nemec, and E. W. M. Hansen. Turbidity current hydraulics and sediment deposition in erodible sinuous channels: laboratory experiments and numerical simulations. Mar. and Petrol. Geol., 41:222–249, 2013. [25] P. Jawahar and H. Kamath. A high-resolution procedure for Euler and Navier-Stokes computations on unstructured grids. J. Comput. Phys., 164(1):165–203, 2000. [26] A. Kassem and J. Imran. Three-dimensional modeling of density current. ii. flow in sinuous confined and uncontined channels. J. Hydraul. Res., 42(6):591–602, 2004. [27] A. Kurganov and D. Levy. Central-upwind schemes for the saint-venant system. M2AN Math. Model. Numer. Anal., 36:397–425, 2002. [28] A. Kurganov, S. Noelle, and G. Petrova. Semi-discrete central-upwind scheme for hyperbolic conservation laws and Hamilton-Jacobi equations. SIAM J. Sci. Comput., 23:707–740, 2001. [29] A. Kurganov and G. Petrova. Central-upwind schemes on triangular grids for hyperbolic systems of conservation laws. Numer. Methods Partial Differential Equations, 21:536–552, 2005. [30] A. Kurganov and G. Petrova. A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system. Commun. Math. Sci., 5:133–160, 2007. [31] A. Kurganov and E. Tadmor. New high resolution central schemes for nonlinear conservation laws and convection-diffusion equations. J. Comput. Phys., 160:241–282, 2000. [32] S. Li and C. J. Duffy. Fully coupled approach to modeling shallow water flow, sediment transport, and bed evolution in rivers. Water Resour. Res., 47(3), 2011. ´ Infante Sedano, and A. Mohammadian. A coupled two-dimensional numerical model for [33] X. Liu, J. A. rapidly varying flow, sediment transport and bed morphology. J. Hydraul. Res., 53(5):609–621, 2015. ´ Infante Sedano, and A. Mohammadian. A robust coupled 2-d model for rapidly varying [34] X. Liu, J. A. flows over erodible bed using central-upwind method with wetting and drying. Can. J. Civ. Eng., 42(8):530–543, 2015. ´ Infante Sedano, and A. Kurganov. Three-dimensional shallow water [35] X. Liu, A. Mohammadian, J. A. system: A relaxation approach. J. Comput. Phys., 333:160–179, 2017. ´ Infante Sedano. Well-balanced central-upwind [36] X. Liu, A. Mohammadian, A. Kurganov, and J. A. scheme for a fully coupled shallow water system modeling flows over erodible bed. J. Comput. Phys., 300(C):202–218, November 2015. [37] K. Motanated and M. M. Tice. Hydraulic evolution of high-density turbidity currents from the brushy canyon formation, eddy county, New Mexico inferred by comparison to settling and sorting experiments. Sediment. Geol., 337:69–80, 2016. [38] F. Necker, C. H¨artel, L. Kleiser, and E. Meiburg. High-resolution simulations of particle-driven gravity currents. Int. J. Multiphase Flow, 28(2):279–300, 2002. [39] W. R. Normark. Observed parameters for turbidity-current flow in channels, reserve fan, lake superior. J. Sediment. Res., 59(3), 1989. [40] A. Ooi, N. Zgheib, and S. Balachandar. Direct numerical simulation of three-dimensional gravity current on a uniform slope. Procedia Eng., 126:372–376, 2015. [41] G. Parker. Conditions for the ignition of catastrophically erosive turbidity currents. Mar. Geol., 46(34):307–327, 1982. [42] G. Parker, Y. Fukushima, and H. M. Pantin. Self-accelerating turbidity currents. J. Fluid Mech., 171:145–181, 1986. [43] G. Parker, M. Garcia, Y. Fukushima, and W. Yu. Experiments on turbidity currents over an erodible bed. J. Hydraul. Res., 25(1):123–147, 1987. [44] O. E. Sequeiros, M. I. Cantero, and M. H. Garcia. Sediment management by jets and turbidity currents with application to a reservoir for flood and pollution control in chicago, illinois. J. Hydraul. Res., 47(3):340–348, 2009. [45] J. S. Turner. Buoyancy effects in fluids. Cambridge University Press, 1979. [46] B. Yu, A. Cantelli, J. Marr, C. Pirmez, C. O’Byrne, and G. Parker. Experiments on self-channelized subaqueous fans emplaced by turbidity currents and dilute mudflows. J. Sediment. Res., 76(6):889–902, 2006.

CE

529

AC

528

27

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

[47] R. Zhang and J. Xie. Sedimentation research in China: Systematic selections. China water and Power Press, Beijing, China, 1993.

CE

580

AC

579

28