GRMS一维河网水质模型

• 0. 前言
• 1. 河网水质控制方程
• 2. 数值求解
• 3. 汊点
• 4. 富营养化模型
  • 4.1. 模型结构
  • 4.2. 富营养化变量方程
  • 4.3. 动力转换项的具体组成形式
    • 4.3.1. 浮游植物模块
    • 4.3.2. 磷循环模块
    • 4.3.3. 氮循环模块
    • 4.3.4. CBOD-DO 子模块

0. 前言

​GRMS 一维河网水质模型基于采用有限体积方法求解 Saint-Venant（圣维南）方程组的一维河网模型，包括基础的对流扩散计算和考虑生物化学反应的富营养化计算。

1. 河网水质控制方程

​描述单一河道污染物输移扩散的控制方程为一维对流扩散方程，具体形式如下：

$$\frac{\partial (AC)}{\partial t} + \frac{\partial (QC)}{\partial x} = \frac{\partial}{\partial x}(A E \frac{\partial C}{\partial x})-S_c+S \tag{1}$$ 其中$Q$为流量（$\text{m}^3/\text{s}$）；$A$为河道过水断面面积（$\text{m}^2$）；$C$为水流输送的污染物质的浓度（$\text{g}/\text{m}^{3}$）；$E$为纵向扩散系数；$t$为时刻（$\text{s}$）；$x$为距离（$\text{m}$）；$S_c$为与污染物浓度有关的衰减源项；$S$为单位时间内、单位河长上的污染物汇流源项（$\text{g}/(\text{m}\text{s})$）。

2. 数值求解

​时间项离散为： $$\frac{\partial(AC)}{\partial t} = \frac{A_{i}^{n+1} C_{i}^{n+1}-A_{i}^{n} C_{i}^{n}}{\Delta t} \tag{2}$$ 式中：上标$n$表示当前时间步，$n+1$表示下一时间步；$A_{i}$表示单元$i$的过水断面面积（$\text{m}^2$）；$C_{i}$表示单元$i$的污染物浓度（$\text{g}/\text{m}^{3}$）；$\Delta t$表示计算时间步长（$\text{s}$）。

​根据采用有限体积法离散的守恒格式圣维南方程的连续性方程，对流项离散为： $$\frac{\partial (QC)}{\partial x} = \frac{(f_1)_{i+1/2}^{*} C_{i+1/2}^{*}-(f_1)_{i-1/2}^{*} C_{i-1/2}^{*}}{\Delta x_i} \tag{3}$$ 式中：$\Delta x_i$表示单元长度（m）；$(f_1)_{i+1/2}^{*}$和$(f_1)_{i-1/2}^{*}$分别表示界面$i+1/2$和界面$i-1/2$处的水流质量通量（$\text{m}^3/\text{s}$）；$C_{i+1/2}^{*}$和$C_{i-1/2}^{*}$分别表示界面$i+1/2$和界面$i-1/2$处的污染物浓度（$\text{g}/\text{m}^{3}$），可根据下式确定： $$C_{i+1/2}^* = \left\{\begin{matrix} C_{i}^n & \text{ if } (f_1)_{i+1/2}^*\geq0 \\ C_{i+1}^n & \text{ if } (f_1)_{i+1/2}^* < 0 \end{matrix}\right. \\ C_{i-1/2}^* = \left\{\begin{matrix} C_{i-1}^n & \text{ if } (f_1)_{i-1/2}^*\geq0 \\ C_{i}^n & \text{ if } (f_1)_{i-1/2}^* < 0 \end{matrix}\right. \tag{4}$$

​扩散项离散为： $$\frac{\partial}{\partial x}(A E \frac{\partial C}{\partial x}) = \frac{A_{i}^{n+1} E_{i}^{n+1}}{\Delta x_i}(\frac{C_{i+1}^{n}-C_{i}^{n}}{\Delta x_{i+1/2}}-\frac{C_{i}^{n}-C_{i-1/2}^{n}}{\Delta x_{i-1/2}}) \tag{5}$$ 式中，$\Delta x_{i+1/2}$表示单元$i$与单元$i+1$之间的距离（m）；$\Delta x_{i-1/2}$表示单元$i$与单元$i-1$之间的距离（m）；纵向扩散系数$E_{i}^{n+1}=0.07 h_{i}^{n+1} u^{*}$，$h_{i}^{n+1}$表示单元$i$的水深（m），$u^*$表示摩阻流速（$\text{m}/\text{s}$），$u^{*}=\sqrt{gR}$，$g$表示重力加速度（$\text{g}/\text{s}^{2}$），$R$表示水力半径（m）。

3. 汊点

​汊点可作为内边界条件为各河段提供污染物浓度边界，汊点处污染物质量守恒关系如下： $$C_{cd}^{n+1}=\frac{\sum_{j=1}^{N_{in}} Q_{j}^{n+1} C_{j}^{n}}{\sum_{j=1}^{N_{in}} Q_{j}^{n+1}} \tag{6}$$ 式中，$C_{cd}^{n+1}$表示汊点内边界的污染物浓度；$Q_{j}^{n+1}$表示流入汊点的第$j$条河段的入流流量；$C_{j}^{n}$表示流入汊点的第$j$条河段的入流污染物浓度；$N_{in}$表示流入汊点的河段条数。

​若汊点可蓄水，则汊点处污染物质量守恒关系如下： $$V_{cd}^{n+1} C_{cd}^{n+1}-V_{cd}^{n} C_{cd}^{n}=\Delta t(\sum_{j=1}^{N_{in}} Q_{j}^{n+1} C_{j}^{n}-\sum_{k=1}^{N_{out}} Q_{k}^{n+1} C_{cd}^{n}) \tag{7}$$ 式中，$V_{cd}$表示汊点水量，$\text{m}^3$；$C_{cd}$表示汊点内边界污染物浓度；$Q_{j}$表示流入汊点的第$j$条河段的入流流量；$Q_{k}$表示流出汊点的第$k$条河段的出流流量；$C_{j}$表示流入汊点的第$j$条河段的入流污染物浓度；$N_{in}$表示流入汊点的河段条数；$N_{out}$表示流出汊点的河段条数。

4. 富营养化模型

4.1. 模型结构

​富营养化模型在美国环境保护局提出的水质模型程序 WASP（the Water Quality Analysis Simulation Program）中的 EUTRO 四个模块基础上，即 CBOD-DO 模块、浮游植物模块、氮循环模块、磷循环模块，构建溶解氧（DO）、碳生化需氧量（CBOD）、浮游植物（PHYT）、氨氮（$\text{N}\text{H}_4$-$\text{N}$）、硝酸-亚硝酸盐氮（$\text{N}\text{O}_3/\text{N}\text{O}_2$-$\text{N}$）、有机氮（DON）、无机磷（DIP）、有机磷（DOP）耦合系统关系。以叶绿素a（Chl-a）质量浓度代表浮游植物的量，DO、CBOD、PHYT、$\text{N}\text{H}_4$-$\text{N}$、$\text{N}\text{O}_3/\text{N}\text{O}_2$-$\text{N}$、DON、DIP、DOP水质指标的质量浓度分别以$C_1、C_2、C_3、C_4、C_5、C_6、C_7、C_8$这8个变量表示。四个模块之间的相互转换关系，见图2，模块中各变量的生成和消失项见表1。

4.2. 富营养化变量方程

描述富营养化变量的方程形式如下：

$$\frac{\partial (AC)}{\partial t} + \frac{\partial (QC)}{\partial x} = \frac{\partial}{\partial x}(A E \frac{\partial C}{\partial x})+S+AS_k \tag{8}$$ 其中$Q$为流量（$\text{m}^3\text{s}^{-1}$）；$A$为河道断面面积（$\text{m}^2$）；$C$为富营养化各水质变量浓度（$\text{g}\text{m}^{-3}$），分别为溶解氧、CBOD、浮游植物、氨氮、硝酸-亚硝酸盐氮、有机氮、无机磷、有机磷；$E$为纵向扩散系数；$t$为时刻（$\text{s}$）；$x$为距离（$\text{m}$）；$S$为单位时间内、单位河长上的污染物排放量（$\text{g}\text{m}^{-1}\text{s}^{-1}$）；$S_k$为与污染物浓度有关的动力转换项，主要由化学和生物过程影响组成。方程中除$AS_k$外的其余各项求解均与水质对流扩散方程相同。

4.3. 动力转换项的具体组成形式

4.3.1. 浮游植物模块

浮游植物在系统动力循环过程中起着重要作用，它影响到水环境中氮、磷、溶解氧等状态变量。其主要反应动力学方程如公式（2）~（4）所示。

$$S_{k3}=C_3(G_p-D_p-\frac{V_{s1}}{D}) \tag{9}$$ 式中第一项为浮游植物生长项，第二项为浮游植物内源性呼吸与死亡项，第三项为浮游植物沉降项；$G_p$为浮游植物生长速率（$\text{d}^{-1}$）；$D_p$为浮游植物呼吸与死亡速率（$\text{d}^{-1}$）；$V_{s1}$为浮游植物沉降速率（$\text{d}^{-1}$）；$D$为水深（$\text{m}$）。

浮游植物生长率方程： $$G_p=K_{g\text{max}}X_{rt}X_{ri}X_{rv} \tag{10}$$ 式中，$K_{g\text{max}}$为 20 摄氏度时浮游植物最大生长速率（$\text{d}^{-1}$）；$X_{rt}$为温度调节因子；$X_{ri}$为光照限制因子；$X_{rv}$为营养盐限制因子。

浮游植物呼吸与死亡率方程： $$D_p=K_{r}(t)+K_{p}+K_{gz} \tag{11}$$ 式中，$K_{r}(t)$为浮游植物内呼吸速率（$\text{d}^{-1}$）；$K_{p}$为浮游植物代谢死亡速率（$\text{d}^{-1}$）；$K_{gz}$为浮游植物被浮游动物摄食死亡率（$\text{d}^{-1}$）。

4.3.2. 磷循环模块

可溶解的或可利用的 DIP 通过吸附-解吸机理与颗粒无机磷相互作用。浮游植物由于生长而吸收 DIP，因此 DIP 合成了浮游植物生物量。通过内源呼吸和非吞食性死亡，磷又从浮游植物生物体中返回到溶解和颗粒有机磷以及溶解无机磷。有机磷通过矿化能转化成溶解无机磷。其主要反应动力学方程如公式（5）~（6）所示：

有机磷： $$S_{k8}=D_{p}C_{3}A_{pc}F_{op}-K_{23}\theta_{23}^{T-20}(\frac{C_3}{K_{mpc}+C_3})C_8-\frac{V_{sop}(1-F_{d2})}{D}C_8 \tag{12}$$

无机磷： $$S_{k7}=D_{p}C_{3}A_{pc}(1-F_{op})+K_{23}\theta_{23}^{T-20}(\frac{C_3}{K_{mpc}+C_3})C_8-G_{p}C_{3}A_{pc}-\frac{V_{sip}(1-F_{d3})}{D}C_7 \tag{13}$$ 有机磷方程式中第一项为浮游植物内源性呼吸与死亡项，第二项为有机磷矿化作用项，第三项为有机磷沉降项；无机磷方程式中第一项为浮游植物内源性呼吸与死亡项，第二项为有机磷矿化作用项，第三项为浮游植物摄取项，第四项为无机磷沉降项；$A_{pc}$为浮游植物中磷碳比；$F_{op}$为浮游植物释放出磷中所含有机磷比例；$K_{23}$为 20 摄氏度时有机磷的矿化速率（$\text{d}^{-1}$），$\theta_{23}$为其温度调整系数；$T$为温度（摄氏度）；$K_{mpc}$为矿化作用半饱和质量浓度（$\text{mg/L}$）；$V_{sop}$为有机颗粒沉降速率（$\text{m/d}$）；$F_{d2}$为有机磷中溶解态比例；$V_{sip}$为无机颗粒沉降速率（$\text{m/d}$）；$F_{d3}$为无机磷中溶解态比例。

4.3.3. 氮循环模块

浮游植物生长吸收氨氮和硝酸-亚硝酸盐，并将其合成浮游植物生物量。吸收氮的速率是氮浓度的函数，而其浓度又与总的可利用无机氮有关。通过内源呼吸和非吞食性死亡，氮又从浮游植物生物量转化为溶解和颗粒有机氮以及氨。有机氮矿化为氨，其矿化速率又依赖于温度，而氨也可以转化成硝酸盐，其硝化速率也依赖于温度和氧气。硝酸盐在缺氧状况下，也可以转化成氮气，其反硝化速率是温度和氧气的函数。其反应动力学方程如公式（7）~（10）所示：

有机氮： $$S_{k6}=D_{p}C_{3}A_{nc}F_{on}-K_{45}\theta_{45}^{T-20}(\frac{C_3}{K_{mpc}+C_3})C_6-\frac{V_{sop}(1-F_{d4})}{D}C_6 \tag{14}$$

氨氮： $$S_{k4}=D_{p}C_{3}A_{nc}(1-F_{on})+K_{45}\theta_{45}^{T-20}(\frac{C_3}{K_{mpc}+C_3})C_6-G_{p}C_{3}A_{nc}P_{NH_3}-K_{56}\theta_{56}^{T-20}(\frac{C_1}{K_{nit}+C_1})C_4 \tag{15}$$

$$P_{NH_3}=C_4[\frac{C_5}{(K_{mn}+C_4)(K_{mn}+C_5)}]+C_4[\frac{K_{mn}}{(C_4+C_5)(K_{mn}+C_5)}] \tag{16}$$

硝酸-亚硝酸盐氮：

$$S_{k5}=K_{56}\theta_{56}^{T-20}(\frac{C_1}{K_{nit}+C_1})C_4-G_{p}C_{3}A_{nc}(1-P_{NH_3})-K_{d}\theta_{D}^{T-20}(\frac{K_{NO_3}}{K_{NO_3}+C_1})C_5 \tag{17}$$ 有机氮方程式中第一项为浮游植物内源性呼吸与死亡项，第二项为有机氮矿化作用项，第三项为有机氮沉降项；氨氮方程式中第一项为浮游植物内源性呼吸与死亡项，第二项为有机氮矿化作用项，第三项为浮游植物摄取项，第四项为硝化作用项；硝酸-亚硝酸盐氮方程式中第一项为硝化作用项，第二项为浮游植物摄取项，第三项为反硝化作用项；$A_{nc}$为浮游植物中氮碳比；$F_{on}$为浮游植物释放出氮中所含有机氮比例；$K_{45}$为 20 摄氏度时有机氮的矿化速率（$\text{d}^{-1}$），$\theta_{45}$为其温度调整系数；$F_{d4}$为有机氮中溶解态比例；$P_{NH_3}$为浮游植物摄取氮营养盐中偏好氨氮程度；$K_{nit}$为硝化作用饱和质量浓度（$\text{mg/L}$）；$K_d$为 20 摄氏度时反硝化速率（$\text{d}^{-1}$），$\theta_D$为其温度调整系数，$K_{NO_3}$为反硝化作用的半饱和质量浓度（$\text{mg/L}$）。

4.3.4. CBOD-DO 子模块

溶解氧含量与其他状态变量相结合。溶解氧的来源有大气复氧和浮游植物的光合作用。溶解氧的消耗主要有浮游植动物的呼吸作用、水体中碳质物质的氧化、硝化作用。其反应动力学方程如公式（11）~（12）所示：

CBOD： $$S_{k2}=A_{oc}D_{p}C_{3}-K_{ox}\theta_{ox}^{T-20}(\frac{C_1}{K_{BOD}+C_1})C_2-\frac{5}{4}\frac{32}{14}K_{d}\theta_{D}^{T-20}(\frac{K_{NO_3}}{K_{NO_3}+C_1})C_5-\frac{V_{sop}(1-F_{d7})}{D}C_2 \tag{18}$$

DO: $$S_{k1}=K_{2}\theta_{2}^{T-20}(C_s-C_1)+G_{p}[\frac{32}{12}+\frac{48}{14}\frac{14}{12}(1-P_{NH_3})]C{3} \\-K_{ox}\theta_{ox}^{T-20}(\frac{C_1}{K_{BOD}+C_1})C_2-\frac{64}{14}K_{56}\theta_{56}^{T-20}(\frac{C_1}{K_{nit}+C_1})C_{4}-\frac{32}{12}K_{rc}\theta_{RC}^{T-20}C_{3}-\frac{O}{D}\theta_{O}^{T-20} \tag{19}$$

CBOD 方程式中第一项为浮游植物死亡项，第二项为 CBOD 氧化项，第三项为反硝化作用项，第四项为含碳物沉降项；DO 方程式中第一项为大气复氧项，第二项为浮游植物生长项，第三项为 CBOD 氧化项，第四项为硝化作用项，第五项为浮游植、动物呼吸作用项，第六项为底泥耗氧项；$A_{oc}$为浮游植物中氧碳比；$K_{ox}$为 20 摄氏度时$\text{CBOD}$氧化速率（$\text{d}^{-1}$），$\theta_{ox}$为其温度调整系数；$K_{BOD}$为 CBOD 半饱和质量浓度（$\text{mg/L}$）；$F_{d7}$为 CBOD 中溶解比例；$K_2$为 20 摄氏度再曝气速率（$\text{d}^{-1}$），$\theta_2$为其温度校正系数；$C_s$为饱和溶解氧量（$\text{mg/L}$）；$O$为 20 摄氏度时底泥耗氧量（$\text{g}/\text{m}^2$），$\theta_O$为其温度调整系数。