大辽河口海域环境容量计算和承载力预警技术研究

鲍思佳，汤立君，陶平，于彩虹，宋成文，邵秘华

（大连海事大学环境科学与工程学院，辽宁 大连 116026）

1 引言

环境容量是在特定时空内有效可交换水体的同化能力，以及在规定环境目标下所能容纳污染物的量。换言之，排入海域污染物浓度不超过环境管理标准的前提下所能允许的最大排放量，是海洋环境容纳污染物的最大容量[1]。河口海域水体污染物成分往往是人为因素影响在自然本底值上的加合，这种加合超过一定限度时，就会改变环境质量[2-3]。因此，环境容量是环境质量评价、环境监测和环境承载力预警的基础，也是河口海域资源开发利用的环境依据。关于环境容量的研究，目前国内外学者已进行了很多工作，海域环境容量的研究也越来越受到重视[4]。

环境容量值的确定有多种方法。20世纪80年代末一般采用浓度测定方法，21世纪随着系统网络应用技术的发展以及对海洋环境容量计算研究的不断深入，其计算方法主要有箱式模型、分担率法和最优化法[5]。特别是随着水环境承载力评估的发展，从初期的局部河流水环境范围内的污染物自净能力研究，再到湖泊、河口和海湾水域的污染物控制政策的实施，水环境承载力的研究逐渐受到关注，同时也发展了海洋环境承载力评估研究[6-8]，并建立了海域环境承载力评估量化方法，如指数评价法、生态足迹法、系统动力学方法和统计学动态模型等[9]。

河口作为海洋和陆地的交界，港口和航运等经济活动频繁，生态系统脆弱。大辽河口位于辽东湾北部海域，为浑河和太子河汇流后在辽东湾的入海口区域，是典型的三角洲河口[10]。由于大辽河口得天独厚的丰富资源，因此成为辽宁省海岸带综合开发的国家重要节点。但大辽河口沿岸企业入海污染严重，同时汇集了辽宁省6个城市的生活污水，废水废物携带量位居全省河流之首[11]，环境治理迫在眉睫。本文选取大辽河口海域为研究对象开展环境承载力研究。刘娟等[12]指出大辽河口及河口邻近海域的主要污染物是化学需氧量（Chemical Oxygen Demand，COD）、无机氮、磷酸盐和重金属等；柴宁[13]也对大辽河水系主要污染物COD和氨氮多年时空变化特征进行了研究，认为两种污染物主要来源于工业和城市的点源排放；而周丹卉[14]在对大辽河口典型污染物时空分布的研究中发现，大辽河口水体的有机物和氮污染最为严重，且季节差异明显。因此，本文选取具有代表性的3种污染要素：COD、无机氮和磷酸盐作为大辽河口海域的主要污染要素开展海域环境承载力预警技术研究。

综上所述，本文以2019年大辽河口海域水环境调查资料为基础，将响应系数改为响应浓度函数，采用反演迭代法计算该海域3种主要污染物的环境容量，再依据海洋功能区划水质标准对大辽河口海域环境承载力进行评估和预警，为该海域的海洋管理部门制定海域污染防治政策提供科学技术支持。

2 计算方法

2.1 数值模型构建

2.1.1 水质模型

鉴于本文研究的海域为近岸浅水区，表层和底层物理、化学和生物性质差别较小，水平方向和垂直方向上的变化可忽略[15]，因此采用深度平均。二维流体动力学方程和水质控制模型方程如下：

式中：ζ为平均海平面以上的瞬时水位高度，单位：m；H为水深，单位：m；Q为径流量；g为重力加速度，单位：m/s2；U和V分别为垂直平均流速在x和y轴上的分量，单位：m/s；f为柯氏参数，f=2ωsinφ，ω为地转角速度，φ为地理纬度；C为谢才系数，C=；n为曼宁系数（取0.025～0.035）；t为时间；x和y平面取在未扰动的平均海平面上，z轴垂直向上，构成右手坐标系（见图1）。式（1）—（3）的初始条件从静止水状态开始，U=V=ζ=0。边界条件分两类：沿岸闭边界，取法向流速等于零（Vn=0）；开边界各点水位为时间的已知函数，即ζ=ζ(t)。

图1 计算坐标系

2.1.2 差分方程

有限差分方法是计算机数值将水质控制模型方程中的导数用网格节点上的函数值的差商代替进行离散，从而求解以网格节点上的值为未知数的代数方程组。

本文采用显隐方向交替（Alternating Direction Implicit，ADI）计算网格（见图2），将式（1）—（3）离散成相应的差分方程。在nΔt→（n+1/2）Δt时间半步长上，将式（1）和式（2）按行进行隐式差分格式离散，V按显示差分格式离散。对后半时间步长（n+1/2）Δt→（n+1）Δt，对式（1）和式（3）按列进行隐式差分格式离散，U按显示差分格式离散，进而建立大辽河口海域水质数值计算模型。

图2 差分网格图

2.1.3 大辽河口海域水动力模型验证

依据上述水动力模型，在2019年5月对大辽河口海域3#污染源进行潮位和流速流向的现场观测（见图3），并与模型计算值对比，得到潮位和流速流向实测值与计算值对比结果（见图4和图5）。由图可见，在一定的误差范围内，该水动力模型在辽东湾近岸浅水海区的潮位和流速流向计算值与实际观测值基本吻合，可作为海区环境容量计算的基础。

图3 模拟污染源位置图

图4 大辽河口海域潮位验证图

图5 大辽河口海域流速流向验证图

2.1.4 保守物质扩散方程

保守物质是在输移扩散过程中没有物质转化和形态改变的物质，扩散方程中无源项和汇项，方程仅考虑对流和扩散作用，只需要确定边界和物质扩散系数[16]。由海区功能区划图可看出（见图6），计算海域位于近岸浅水区，水浅垂直混合充分，可将垂直方向上的水体视作均匀不变，故本文采用保守物质扩散模型，即二维垂直平均扩散方程进行计算：

图6 大辽河口海域功能规划网格图

式中：P为扩散物质的深度平均浓度；p=；Sm为污染物源强；Kx和Ky为扩散系数，由Elder公式确定：

2.2 污染源输入负荷量-输出响应浓度线性系统验证

数学物理方程主要有3种类型：椭圆型、抛物型和双曲型。有学者认为污染物在海水中的扩散方程是椭圆型的线性方程，且扩散方程在第一边界条件下是线性的，符合叠加原理[17]。为了验证在环境容量的计算中，扩散方程是椭圆线性的数学物理方程，本文将通过以下扩散数值模拟计算出污染源输入负荷量-输出响应浓度场是线性系统。

以大辽河口海域为例，从北至南依次设置5个模拟污染源，分别为1#、2#、3#、4#和5#（见图3）。利用保守物质扩散方程的数值模型对1#—5#污染源各自做输入负荷量-输出响应浓度场模拟计算，使用的各项参数如下：

（1）计算其中任意一个污染源的输入负荷量-输出响应浓度时，要计算的污染源输入负荷量分别取“1 t/d”、“2 t/d”和“5 t/d”，其他污染源输入负荷量均为“0 t/d”；

（2）开边界浓度为“0 mg/L”（增量计算）；

（3）计算30个潮周期（浓度场已稳定）结束。

计算结果见表1。

表1 输入负荷量为1～5 t/d时5个污染源排污口处输出响应浓度（单位：mg/L）

表1中行数据表示各污染源单独输入1 t/d、2 t/d或5 t/d（其他污染源输入均为0 t/d）时，各污染源排污口处输出响应浓度和受影响浓度值；列数据表示各污染源单独输入1 t/d、2 t/d或5 t/d（其他污染源输入0 t/d）时产生的输出响应浓度，以及受其他污染源输入1 t/d、2 t/d或5 t/d时产生的影响浓度值。

用表1中数据对1#—5#污染源输入负荷量和输出响应浓度作相关分析，得到回归方程：

式中：Si为第i个污染源输入负荷量，单位：t/d；Ci为第i个污染源排污口处输出响应浓度，单位：mg/L；KCi为回归系数；bCi为回归常数。

式中：Ci为输出响应浓度，单位：mg/L；Si为污染源输入负荷量，单位：t/d；KSi为回归系数；bSi为回归常数。

式（6）和式（7）是由坐标原点引出的斜率为KCi和KSi的直线方程，bCi和bSi值很小，接近“0”，则式（6）和式（7）可简化为：

将表1中各污染源输出形成的响应浓度代入式（6）中，计算1#—5#污染源输入的负荷量，同给定的负荷量值进行比较，结果见表2。

从表2可看出，由相关方程计算的各污染源负荷量和给定的负荷量只相差±0.002，这一差值在误差范围内。

表2 计算值与给定值比较表

综上所述，通过扩散方程的数值模拟计算结果和相关分析，以及式（8）和式（9）计算得到的负荷量值与给定值在误差允许的范围内相等，证明了污染源输入负荷量-输出响应浓度是线性系统。

2.3 建立污染源排放量-输出响应浓度函数的两个相关方程

根据已建立的水质模型计算大辽河口海域各污染源不同单位排放量时的输出响应浓度函数值（一个潮周期各正点浓度），再建立各污染源排放量-输出响应浓度和输出响应浓度-污染源排放量的两个相关方程，即：

式中：Si为第i个污染源污染物的负荷量，单位：t/d；KCi为第i个方程回归系数；Ci为第i个控制点处污染物响应浓度，单位：mg/L。

式中：Ci为第i个控制点处的污染物响应浓度，单位：mg/L；KSi为第i个方程回归系数；Si为第i个污染源污染物的负荷量，单位：t/d。

因此，水质模型是线性系统方程，在初始浓度和边界浓度均为0的条件下，污染源负荷量和输出响应浓度成线性关系，而且是由原点引出的一条斜率为KCi（KSi）的直线，符合叠加原理。

若预警海域有n个污染源，则式（10）化为：

式中：Si为第i个污染源污染物的负荷量，单位：t/d；KCj为第j个污染源单位输出负荷量对第j个污染源控制点产生的响应浓度系数；cj为第j个污染源输出负荷量产生的响应浓度；，是（n-1）个污染源对第n个污染源产生的影响浓度总和。

因此，可以根据各污染源污染要素的响应浓度推算出该海区污染物的负荷量，这是建立计算海区负荷量方程的基础。

2.4 反演迭代法

目前，常见的环境容量计算方法是基于线性规划的排海通量最优化法（简称线性规划法）[18-19]。该方法能够应用于污染情况复杂的海域，但计算过程繁琐。为了减少计算工作量和避免污染源间负荷量调配的过程，基于上文构建的两个可逆方程，本文采用“反演迭代法”。

反演方法在海洋学中的首次应用是在水文潮流场上，它是在已知或假定的动力学约束条件下进行统计推断的一般方法[20]。本文中，反演法是指已知方程Si=KCiCi中浓度Ci的值，倒推算相应污染源的负荷量Si的值；若有多个污染源，则各污染源相应控制点处的输出响应浓度互受影响，是一个动态趋于稳定的过程，可以通过线性系统的叠加原理来完成计算。上述两种计算过程合称为“反演迭代法”。

本文根据大辽河口海域国家功能区划水质标准，通过反演迭代法即可算得该功能区污染源的标准环境容量；再将该功能区控制点处现场水质（COD、无机氮和磷酸盐）观测浓度资料代入方程Si=KCiCi中，即可算得该功能区污染源的水质现状环境容量。各污染源的标准环境容量和现状环境容量，是评估海域环境承载力等级的关键数据。

2.5 有效浓度

在环境容量计算中，根据《水质达标方案技术指南》，一个潮周期内响应浓度的取值有两种方法：一种是取一个潮周期内最高浓度和最低浓度的平均值；另一种是取一个潮周期各潮时的响应浓度的平均值。但污染源输入-输出响应浓度场在一个潮周期内随潮时而变化，即不同潮时的响应浓度也不同。潮汐规律不同的海区，响应浓度随潮时的变化规律也不同。根据实测资料可知，大辽河口海域是不规则的半日潮海区[21]，接近标准的半日潮流变化，响应浓度随潮时成正弦函数曲线变化。本文借鉴物理电工学中正弦交流电“有效值”的概念，提出响应浓度的“有效浓度”值概念，即取一个潮周期内响应浓度最大值的0.707倍作为响应浓度进行计算。因此，响应浓度用下式计算：

式中：Ci为有效浓度；Cimax为一个潮周内响应浓度最大值；K为有效浓度系数，取值0.707。式（13）适用于标准半日潮和全日潮（或接近两种潮汐）类型海区。

3 大辽河口海域污染物负荷等级预警技术研究

3.1 大辽河口海域污染物负荷等级预警方程和评价体系

引用《环渤海污染压力和海上响应的统筹调控研究》对入海负荷超标率计算及等级评估方法[22]，即式（14）计算入海负荷超标率，表3海域环境承载力等级评价体系可作为海洋环境承载力评估的依据。

表3 入海负荷等级评价体系

式中：Kj为某种污染物第j个功能区入海负荷超标率；为某种污染物第j个功能区的水质现状环境容量；为某种污染物第j个功能区的水质标准环境容量；排放量单位均为t/a。

3.2 大辽河口海域环境承载力预警系统

依据上文所述计算方法和评价体系，运用Visual Basic 6.0计算机语言编写了大辽河口海域污染物承载力预警系统。使用者根据需要在界面输入污染物的各项浓度，可以计算污染物的环境容量，再由此对海域环境承载力作出判断。该系统操作简单，可以在Windows系统下直接运行。系统主要由3部分组成：现状环境容量计算，标准环境容量计算，海区承载力等级评价和预警。该系统进行承载力预警的基本流程如下：

（1）打开运行文件，进入系统主界面，选择海区，以“大辽河口”为例，点击计算（见图7）。

图7 选择海区界面

（2）通过反演迭代法，分别计算大辽河口海域3个污染源处各污染要素的现状环境容量：在程序界面中分别输入3种污染要素的现状负荷量以及本底浓度，然后得到各污染源现状环境容量计算结果（见图8a）。

（3）通过反演迭代法，分别计算大辽河口海域3个污染源处各污染要素的标准环境容量：在程序界面中分别输入3种污染要素的标准浓度以及本底浓度，得到各污染源的标准环境容量（见图8b）。

图8 大辽河口海域环境容量计算结果

3.3 大辽河口海域环境承载力预警结果

3.3.1 海域和站位布设

大辽河口海域面积736 km2。海域各污染源坐标见表4，海区功能规划网格图见图6。

表4 大辽河口海域污染源坐标

3.3.2 大辽河口海域水质调查结果

2019年5月9日对营口海域10个站位的3种化学要素（COD、无机氮和磷酸盐）进行了调查，按《海洋监测规范》方法执行，水质调查结果见表5。

表5 营口海区水质调查结果

3.4 大辽河口海域入海污染物承载力预警结果

依据反演迭代法计算得到的现状环境容量和标准环境容量，以及3.1节的预警方程对大辽河口海域环境承载力进行评价，结果见图9。该海域COD和磷酸盐是1级可载状态，无机氮标准环境容量小于0，是6级严重超载状态。

图9 大辽河口污染物环境承载力预警等级结果

4 结论

基于二维流体动力学方程和水质控制模型，运用保守物质扩散方程和反演迭代法计算出大辽河口海域3种主要污染要素的环境容量，以海洋功能区划水质标准为基准进行承载力评估，具体结论如下：

（1）以大辽河口不同类型功能区环境管理要求为基础，运用反演迭代法得到该海域3种主要污染物COD、无机氮和磷酸盐的标准环境容量和现状环境容量，并进行海域环境承载力状态评估。结果如下：COD和磷酸盐处于1级可载状态，无机氮处于6级严重超载状态。

（2）对于辽东湾海区，在第一边界条件下，使用保守物质水质扩散方程计算的各项因子的浓度与污染负荷量是线性的、齐次的和可叠加的，在此基础上构建的反演迭代法为海洋环境容量计算提供了新思路。

（3）在标准的半日潮或全日潮海区，由于响应浓度随潮时呈正弦函数曲线变化，故一个潮周期响应浓度的取值，可借鉴物理电工学理论有关正弦交流电的概念，方法为：取0.707的有效系数乘以一个潮周期内响应浓度最大值表示“有效浓度”，这一取值方法在标准半日潮和全日潮（或接近两种潮汐）类型海区具有参考价值。

（4）本文编制的海域环境承载力预警系统软件，实现了海域环境承载力等级评价的计算机程序化，基于水质调查结果，对辽东湾海域环境承载力进行分级评价和预警，以此为海洋管理部门合理规划海域内污染物的排放提供科学的调控依据。