深度平均法论文研读：二维泥沙输运模型

原文

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$$

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

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

侵蚀速率 (Partheniades(1965))

$$E_f = \begin{cases} M\left( \dfrac{\tau_b}{\tau_e} - 1 \right), & \text{if } \tau_b > \tau_e \text{ and } C_{sb} > 0, \\ 0, & \text{otherwise}. \end{cases}$$

$$E_p = \begin{cases} M\left( \dfrac{\tau_b}{\tau_{e,p}} - 1 \right), & \text{if } \tau_b > \tau_{e,p} \text{ and } C_{sb} = 0, \\ 0, & \text{otherwise}. \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).

连续方程的推导

引入自由表面的运动学条件和底部运动学条件

$$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$$