原文

A depth-averaged two-dimensional sediment transport model for environmental studies in the Scheldt Estuary and tidal river network - ScienceDirect

附PDF: https://pan.baidu.com/s/1XmZkQXuFonfgSAhfdp6Wrw?pwd=pnkh

工作内容

针对SLIM有限元模型(Second-generation Louvainla-Neuve Ice-ocean Model)的二维深度平均分量和一维截面平均分量设计了泥沙模块, 并将其应用于斯海尔德盆地的潮汐部分. 该泥沙输运模块侧重于细粒、粘性沉积物. 它是进行环境生物地球化学研究的必要工具, 其中细粒沉积物动力学起着至关重要的作用。

流体力学和盐度模型

河口部分使用二维深度平均方程, 河流部分使用一维截面平均方程

浅水方程(Shallow water equations)

$$
\frac {\partial \eta}{\partial t} + \nabla \cdot(Hu) = 0
$$

$$
\frac{\partial u}{\partial t} + u \cdot (\nabla u) + fk \times u = -g \nabla \eta - \frac{1}{\rho}\nabla P_a + \frac{1}{H}\nabla \cdot (H\nu (\nabla u)) + \frac{\tau_s - \tau_b}{\rho H}
$$

其中 t 是时间, $H=h+\eta$ 是水深, h 是水柱参考高度; $f=2\omega\sin \phi$为科里奥利参数, 以$\omega$为地球角速度, 以$\phi$为纬度; k为单位向上矢量; G为重力加速度; $\rho$ 是水的密度, 假设水的密度是恒定的, $P_a$ 是水面的大气压力; $\nu$ 为水平涡动粘度; $\tau_s$和$\tau_b$分别为表面和底部应力矢量.

深度平均盐度运输方程

$$
\frac{\partial (HS)}{\partial t} + \nabla\cdot (HuS) = \nabla \cdot (Hk\nabla S)
$$

其中k为扩散系数

水平涡粘性系数

根据 Smagorinsky(1963)的工作,水平涡粘性系数 $\nu$ 通常按照以下形式给出:

$$
ν = C_s^2 * Δx * |S|
$$
这里,

$C_s$ 是 Smagorinsky 常数,一般取值在 0.1 到 0.2 之间;
$Δx$ 是网格间距,代表了模型分辨率;
$|S|$ 是应变率张量 S 的模,它描述了速度场的空间变化,具体为:
$$
S = (∇u + (∇u)^T) / 2
$$
其中 $∇u$ 是速度梯度张量,而上标 T 表示转置。应变率张量的模可以通过下式计算:

$$
|S| = sqrt(2 * S_{ij} * S_{ij})
$$
这里的 $S_ij$ 表示应变率张量的分量。

曼宁-斯特里克勒公式(the Chézy–Manning–Strickler formulation)

$$
\tau_b = \rho g n^2 \frac{u \begin{Vmatrix} u \end{Vmatrix} }{H^{1/3}}
$$

其中 n 是曼宁系数. 在大陆架等于$0.0235\ s\ m^{-1/3}$ , 在Antwerp附近, 从河口线性增加到$0.028\ s\ m^{-1/3}$

扩散系数 k (Okubo)

$$
k = c_k \Delta^{1.15}
$$

其中 $c_k$ 为常数, $\Delta$ 是网格的特征长度尺度(即二维网格中三角形的最长边, 或一维网格中线段的长度). 对其值进行校准以准确地表示Scheldt中的盐度变化, 从而得到 $c_k = 150\ m^{0.85}s^{-1}$.

泥沙模型

泥沙模块考虑了三层: 悬浮沉积物所在的水柱, 由底部新沉积的沉积物组成的一层, 以及后者下方的一个母层。虽然考虑了三层, 但该模块仅由两个相互作用的变量组成, 即$C_{ss}$, 悬浮沉积物的深度平均浓度[$kg\ m^{−3}$], 以及$C_{sb}$, 新鲜层中底层沉积物的浓度[$kg\ m^{−2}$]。母层是沉积物的无限来源, 只有当$C_{sb}$局部为零时才会被侵蚀。

悬浮中的沉积物是通过平流和扩散输送的, 而底部的沉积物则不是。在区域的二维部分, 悬浮泥沙浓度和新沉积泥沙浓度服从如下方程:
$$
\frac{\partial (HC_{ss})}{\partial t} + \nabla \cdot (HuC_{ss}) = \nabla \cdot (Hk\nabla C_{ss}) + E_f + E_p - D
$$

$$
\frac{\partial C_{sb}}{\partial t} = D - E_f
$$

image-20241208170632962

式中$E_f$为新鲜层沉积物的侵蚀速率, $E_p$为母层沉积物的侵蚀速率, $D$为新鲜层沉积物的沉积速率。

假设母层从未供应, 并且两个底层之间没有交换。此外, 也没有考虑床质输运。

侵蚀速率 (Partheniades(1965))

$$
E_f = \begin{cases} \begin{aligned}
&M\left( \frac{\tau_b}{\tau_e} -1\right) & & if\ \tau_b > \tau_e\ and\ C_{sb} > 0 \newline
&0 & &otherwise
\end{aligned}
\end{cases}
$$

$$
E_p = \begin{cases} \begin{aligned}
&M\left( \frac{\tau_b}{\tau_{e,p}} -1\right) & & if\ \tau_b > \tau_{e,p}\ and\ C_{sb} = 0 \newline
&0 & &otherwise
\end{aligned}
\end{cases}
$$

当$\tau_b$(底部应力向量的范数$\tau_b$) 高于一个阈值$\tau_e$时, 沉积物从新鲜层被侵蚀, 或者, 如果新鲜层局部是空的, 当$\tau_b$高于另一个阈值$\tau_{e,p}$; M被称为侵蚀速率参数。

沉积速率(Einstein和Krone (1962))

$$
D=w_sC_{ss}
$$

其中 $w_s$ 为沉降速度.
$$
w_S = w_{s,0} \left( \frac{C_{ss}}{C_{ss,0}}\right)^m
$$
其中 $C_{ss,0}=0.1\ kg m^{-3}$ 是 Scheldt SSC的参考值, $w_{s,0}$ 是相关的参考沉降速度, 其值是经验确定的, m 是0.5到3.5之间的系数(van Leussen, 1999).

连续方程的推导

image-20241211171404424

引入自由表面的运动学条件和底部运动学条件
$$
u_z|_{z=z_s} = \frac{dz_s}{dt} = \frac{\partial z_s}{\partial t} + \frac{\partial z_s}{\partial x} u_x + \frac{\partial z_s}{\partial y} u_y
$$

$$
u_z|_{z=z_b} = \frac{dz_b}{dt} = \frac{\partial z_b}{\partial t} + \frac{\partial z_b}{\partial x} u_x + \frac{\partial z_b}{\partial y} u_y
$$

对连续性方程沿水深方向进行积分
$$
\int_{z_b}^{z_s}\frac{\partial u_j}{\partial x_j} dz = 0 \ (j=1,2,3)
$$

$$
\int_{z_b}^{z_s}\frac{\partial u_j}{\partial x_j}dz + u_z|{z_s} - u_z|{u_b} = 0 \ (j=1,2)
$$

由牛顿-莱布尼兹公式
$$
\frac{\partial}{\partial x_j} \int_a^b f dz = \int_a^b \frac{\partial f}{\partial x_j} dz + f|b \frac{\partial b}{\partial x_j} - f|a \frac{\partial a}{\partial x_j}
$$
将上式第一项展开,可得
$$
\frac{\partial }{\partial x_j}\int
{z_b}^{z_s}u_jdz - u_j|
{z_s}\frac{\partial z_s}{\partial x_j} +u_j|{z_b}\frac{\partial z_b}{\partial x_j} + u_z|{z_s} - u_z|{z_b} = 0 \ (j=1,2)
$$
代入两个运动学条件,可得
$$
\frac{\partial }{\partial x_j}\int
{z_b}^{z_s}u_jdz - u_j|{z_s}\frac{\partial z_s}{\partial x_j} +u_j|{z_b}\frac{\partial z_b}{\partial x_j} + \frac{\partial z_s}{\partial t} + \frac{\partial z_s}{\partial x_j} u_j|{z_s} - \frac{\partial z_b}{\partial t} - \frac{\partial z_b}{\partial x_j} u_j|{z_b} = 0 \ (j=1,2)
$$
化简,得
$$
\frac{\partial }{\partial x_j} \int_{z_b}^{z_s} u_j dz + \frac{\partial z_s}{\partial x_j} - \frac{\partial z_b}{\partial x_j} = 0 \ (j=1,2)
$$
由 $H=z_s-z_b$, 进一步化简为
$$
\frac{\partial H}{\partial t} + \frac{\partial }{\partial x_j} \int_{z_b}^{z_s} u_j dz = 0 \ (j=1,2)
$$
又有
$$
\int_{z_b}^{z_s} u_j dz = HU_j \ (j=1,2)
$$
代入上式,得
$$
\frac{\partial H}{\partial t} + \frac{\partial HU_j}{\partial x_j} = 0 \ (j=1,2)
$$

$$
\frac{\partial H}{\partial t} + \frac{\partial HU_x}{\partial x} + \frac{\partial HU_y}{\partial y} = 0
$$
本篇文章中 $H = h + \eta$, $H$ 为水深,$h$ 为参考水柱高度(河床到基准面的高度),$\eta$ 为水面相对于基准面的高度,是一个瞬时变化量。

因此,第一项水深随时间的变化率写为 $\frac{\partial \eta}{\partial t}$,第二三项的为水平方向的水流通量,其通过面积可表示为 $H$, 因此此项应写为 $\frac{\partial Hu_j}{\partial x_j}$.

综上
$$
\frac{\partial \eta}{\partial t} + \frac{\partial Hu_j}{\partial x_j} = 0 \ (j=1,2)
$$
使用哈密顿算子可以表示为
$$
\frac{\partial \eta}{\partial t} + \nabla \cdot (Hu) = 0
$$