秦皇岛筏式养殖对水动力和污染物输运的影响

匡翠萍， 董智超， 顾 杰， 单云驰

(1.同济大学 土木工程学院， 上海 200092； 2.上海海洋大学 海洋生态与环境学院， 上海 201306； 3.华东电力设计院有限公司， 上海 200063)

随着海洋经济的快速发展，水产养殖产品和新食品供应的需求不断增加[1]，海产养殖规模不断扩大，我国海水养殖规模从1978年至2014年30余年间增长了22倍[2]，总产量占全球海水养殖的61.2%[3]，近岸海域水体富营养化程度不断加剧[4]，近海养殖已经成为海洋水环境污染的主要原因之一[5].随着近海筏式扇贝养殖产量增加和养殖面积扩大，大面积海水养殖网笼的布置，对海域水动力条件产生了一定的阻碍作用，改变了养殖区及邻近海域的水动力特性，由于浅海养殖阻碍了海水的正常运动，减弱了水动力强度，水体交换能力下降，交换周期增长，极大改变了区域的污染物扩散和水环境质量[6]，2011年5月桑沟湾爆发了大规模的赤潮[7]，2016年7月和8月秦皇岛发生了两次赤潮现象[8]，可见筏式养殖带来的生物沉积和对水体交换能力的减弱在一定程度成为了赤潮爆发的诱因[9].因此，筏式养殖对区域水动力和污染物输运影响关系的研究显得极其重要.

通过对养殖区流场的观测，可以发现海水养殖活动对养殖海域动力环境有显著的影响.Boyd和Heasman[10]测量了南非Saldanha湾贝类养殖区内外的流速，养殖区外的流速是养殖区内的6倍.李铁军等[11]对浙江三门湾筏式养殖区前后进行观测，并使用潮汐观测资料进行调和分析，养殖区前后减少的流速变化率约为6%～22%.项福亭等[12]对蓬莱湾扇贝养殖区流速调查发现1975年至1990年15年间养殖中心区的最大流速减少了近75%，最小流速减少了近90%.陈清满[13]发现枸杞岛海域紫贻贝养殖设施对养殖海域水交换产生严重阻滞作用.赵东波等[14]研究了罗源湾养殖区内潮动力结构的特征及其对养殖活动的响应，高密度的养殖设施和生物的存在，阻碍了海湾水体的流动，产生了双边界层结构并减弱了潮汐不对称效应.张泽华等[15]以黑泥湾为研究对象，分析了浅海海带筏式养殖对周边海域悬浮体动力特征的影响作用，表明筏式养殖改变了潮流流速和悬浮体质量浓度的垂向分布特征，抑制了海底泥沙的再悬浮作用.季如宝等[16]在对贝类养殖的海湾生态系统研究中指出，贝类密集区生物的沉降作用非常明显.Hatcher等[17]在加拿大Upper South Cove贻贝养殖区进行了测定，养殖区的沉降量往往是非养殖区的2倍以上.Bouchet和Sauriau[18]对太平洋贝类养殖区使用生态质量评估发现，悬吊养殖系统比底播养殖有机物富集量更高.

随着计算机技术的发展，数学模型在海水养殖数值模拟研究上得到了广泛应用.Panchang等[19]认为水动力和粒子追踪数学模型，是研究养殖区污染物的扩散和运移规律及其物质归宿是较为有效的手段.Fan等[20]运用水动力模型研究养殖区对水流的垂向结构的分布的影响，流速垂向结构分布主要由养殖区拖曳力和底摩阻共同控制.Durate等[21]基于三维模型计算加利西亚河口湾的水动力特性，通过对余流场进行分析发现筏式养殖可减少近40%的余流流速，这将助于有害藻华的发生.樊星等[22-23]基于水动力模型和现场观测，研究养殖密度对潮流垂直结构的影响机制，提出潮流上边界层的概念，通过建立双阻力模型得到养殖区潮流上边界层平均阻力系数，边界层的厚度受到养殖密度影响.史洁和魏浩[24]使用优化的POM模型对桑沟湾养殖区进行模拟，发现桑沟湾高密度的养殖使流速平均减小40%，平均半交换时间增长71%.

秦皇岛养殖功能用海分布在洋河口渔港航道以南至滦河口唐山分界以北的海域，区域内底栖生物丰富，水体透明度高，水质良好，符合渔业水质标准[25-26].秦皇岛海域海湾扇贝养殖业快速发展[27]，为了降低筏式养殖对近岸海域环境的污染，需要对筏式养殖对动力环境的影响机理进行综合分析.国内外的研究学者关于筏式养殖对水平和垂向的水流结构的影响已经有了诸多研究，但对养殖区周边海域动力和水质的综合性影响分析较少.针对秦皇岛筏式养殖区，课题组已经通过物理模型试验，研究了筏式养殖区垂向流速分布特征，探究了养殖区内外流速在不同水深、来流速度、养殖密度条件下的变化情况[28]，并利用明渠流理论，以对数函数曲线与抛物线曲线对流速垂向分布进行了拟合和分析，指出养殖密度对流速垂向分布公式的影响特性[29].本文基于物理模型试验对筏式养殖垂向流速分布的研究成果，建立秦皇岛海域水动力和水质模型，研究筏式养殖区对海域的水动力和污染物输运的影响机制，为海水养殖的规划和布局及海洋环境评价提供一定的理论依据.

1 数学模型的建立和验证

1.1 控制方程

基于物理模型试验的成果对动力垂向结构进行概化，采用二维数学模型研究筏式养殖区对秦皇岛海域的动力环境水平结构的影响.水动力模型的控制方程是建立在二维浅水方程的基础上，在笛卡尔坐标系下，通过对三维水平动量方程和连续方程沿深度进行积分，得到二维浅水方程如下：

(1)

(2)

(3)

式中：η表示水位，h=η+d，h为总水深，d为静止水深；S为源流量，uS和vS表示源流入在x和y方向上的速度；t为时间；u与v表示x和y方向上的水深平均流速；g为重力加速度；f为柯氏力系数；ρ为水体密度，τsx和τsy表示x和y方向上的风应力；τbx和τby表示x和y方向上的底部摩擦力；侧应力Tij包括黏性摩擦和紊动摩擦等.通过对涡黏方程沿水深平均流速梯度方向计算如下：

(4)

水质模型的基本控制方程为

(5)

式中：C为浓度标量；kp表示污染物的衰减率；CS为源处浓度标量；FC为水平扩散项；A为水平涡黏系数.

(6)

式中：Dh为水平扩散系数.

对于营养盐等非保守物质，在随流平流、扩散迁移的同时，由于大气挥发作用、微生物降解作用、浮游生物的富集作用以及悬浮颗粒的吸附作用等，其本身的性质和质量也在发生变化.其中的生物化学作用一般与温度、盐度等条件有关，本文以COD(化学需氧量)为例，污染物降解考虑一级衰减过程，满足一级反应方程，即

(7)

式中:K(T,S)COD为降解速率，若K(T,S)COD=0，为保守物质.

1.2 模型区域和网格

为保证模型计算的合理性和准确性，本文采取大小模型嵌套模式进行计算，使用验证过的大区域模型为小区域模型提供计算边界.大区域模型为渤海模型[30]，使用大连老虎滩至烟台两个潮位站连线作为开边界[31].小区域模型包含整个养殖区和邻近海域，小区域模型沿岸南起滦河以南20 km，北至山海关以北约25 km，向外海延伸约60 km，东北和西南侧外海开边界长约50 km，东南侧外海开边界长约115 km，区域范围为119°7′53.93″E ～ 120°29′29.83″E，39°1′32.43″N ～ 40°5′55.45″N.秦皇岛计算网格和养殖区位置如图1所示，数学模型采用非结构三角网格，网格东西跨度约110 km，南北跨度约120 km，网格共计9 442个节点，三角形单元共计17 973个.网格分辨率按研究需要对近岸区域网格进行加密处理，外海网格则相对稀疏.网格最高分辨率约为15 m，最低分辨率约为7 km.近岸区域较高分辨率有利于合理表现近岸地形的变化和曲折岸线的处理，远海区较低分辨率有利于合理概化外海地形并减少模型计算时间.筏式养殖区主要集中于洋河和滦河之间，老虎石以北有少量分布.本文根据卫星遥感资料显示的养殖区密度，将养殖区概化成高密度区(Ⅰ)和低密度区(Ⅱ).

1.3 边界条件和参数设定

模型外海开边界N、E、S使用潮位过程和流速过程的复合边界条件来驱动.其中，潮位和流速过程由渤海潮流模型计算提供.模型陆向开边界为11个河口，河流边界通过河流实测月平均流量控制.COD是表征有机污染的一个综合因子，也是海域定量描述受污染程度的重要指标之一.很多研究将COD与赤潮的爆发直接联系起来，甚至将其作为海域富营养化的重要指标之一，而且COD受生物活动的影响相对来说比营养盐小，它的生化降解作用也较容易确定[31].因而本文选取COD作为计算边界.污染物COD浓度根据外海污染物平均浓度取值，N、S边界COD浓度取1.2 mg·L-1，E边界COD浓度取1.3 mg·L-1.河流入海条件，见表1,由7、8月实测月平均流量及入海污染物浓度提供.

a 研究区域网格

b 筏式养殖区范围

模型计算参数主要包括时间步长、CFL数、干湿边界判定参数、水平涡黏系数、底床摩擦力、科氏力等.模型的时间步长设置为0.01～1 s的变化步长，使CFL数小于1，可以有效地解决模型不稳定的问题；干湿边界判定参数中，干水深hdry取值0.005 m，淹没水深hflood取值0.05 m，湿水深hwet取值0.1 m；水平涡黏系数设定采用Smagorinsky亚网格尺度模型计算，Smagorinsky系数cs取常数0.28；科氏力由模型所在区域的纬度确定.

表1 秦皇岛入海河流2013年7、8月平均流量(m3·s-1)及COD浓度(mg·L-1)

底床摩擦力是二维水动力模拟中的一个重要参数，本文采用曼宁数作为衡量底床摩擦力的指标.养殖区以外海域曼宁数M根据以往该区域的研究[32]取均值为74 m1/3·s-1.秦皇岛海域筏式养殖区由多个子养殖区构成，各子养殖区间隔大于2 km定义为低密度养殖区，间隔小于2 km定义为高密度养殖区.基于物理模型试验养殖区曼宁数参考物模试验的率定结果取值[33]，相应高密度养殖区曼宁数取值为20 m1/3·s-1，低密度养殖区曼宁数取值取值为30 m1/3·s-1.

水温根据秦皇岛近岸海域实测资料，因7、8月份水温较接近，取实测平均值为25.8 ℃；海区初始COD浓度设为1.3 mg·L-1；COD的阿伦尼乌斯温度系数θCOD取1.02.根据秦皇岛近岸海域污染物分布规律，本文取COD在20 ℃时的降解系数K20=0.001～0.032 5[33].

1.4 模型验证

采用养殖区海域实测潮位、潮流和COD浓度资料对水动力模型和水质模型进行验证，验证点位如图2a所示.潮位验证(WL1点)采用2013年5月10日8:00至2013年5月12日8:00秦皇岛潮位站点的实测潮位资料，流速和流向验证(SDL07及SDL09点)采用2013年5月11日8:00至2013年5月12日8:00秦皇岛海域潮流站点实测流速和流向资料，污染物COD的验证使用2013年秦皇岛8月期间28个实测站(Q01-Q28)所测COD浓度进行验证.潮位、流速、流向及COD浓度的验证结果如图2所示.潮位、流速和流向的模型计算值与实测值相位一致、数值基本吻合，拟合较好.COD计算值取测量日期COD浓度日平均值，28个验证点位计算结果与实测值误差较小，除Q21点外，绝对误差均在±0.45mg·L-1以内，相对误差均在30%以内.

图2 验证点位，潮位、流速、流向及COD浓度验证

为了评价模型计算结果的优劣，比较计算值和实测值之间的差异程度，本文采用数学模型中常使用的模型评价方法进行模型计算评价.百分比模型偏差法[34]用来计算模型的模型计算值M和相应实测值D之间的百分比偏差率，计算方法如下：

(8)

当P≤10，拟合结果极好；当1040，拟合结果差.使用百分比模型偏差法对秦皇岛水动力和水质模型计算效率进行评价，计算得到潮位百分比偏差为29、流速、流向百分比偏差均小于20，COD浓度28个站位总评价百分比偏差为12.5，所有站点评价均为好或非常好，该模型可以用于计算分析养殖区对海域水动力和污染物输运的影响.

2 筏式养殖对水动力的影响

为了与养殖区存在情况下的近岸海域水动力进行对比，模拟无养殖区情况作为模型的对照组.采用有养殖区流场与无养殖区流场流速相减，得到流速差值，如图3所示.I为高密度养殖区，Ⅱ为低密度养殖区.高密度养殖区流速有明显减小；低密度养殖区流速也有所减小，但减小程度比低密度养殖区小；滦河口外由于岸线走向改变沿岸流的方向，流速减小；其他区域减小幅度自近岸向外海递减.有养殖区情况下养殖区Ⅰ内流速明显减小，Ⅰ区近岸处流速减少幅度约为20%～50%，Ⅰ区向海侧流速减少幅度为10%～20%；Ⅱ区流速变化较小，减少幅度小于10%；因涨落潮流方向为NE-SW向，非养殖区的沿岸区域也有一定流速减小，减小幅度小于10%；受养殖区布置及走向影响，养殖区外侧向外海延伸区域，流速增加，增加幅度小于10%.

图3 平均流速变化场(有养殖区无养殖区)

Fig.3 Distribution of mean flow velocity change caused by the scenarios with and without culture zone

随着人们对海产品的需求不断增加，养殖面积和养殖密度都在逐步加大，养殖密度增加使养殖区对水流的阻碍更加显著.为探究养殖区密度的变化对秦皇岛近岸海域水动力影响程度和影响范围，通过改变模型区域的曼宁数场来实现养殖密度的增加.由于高密度养殖区阻水能力更强，养殖区密度越高，曼宁数越小.在模型地形区域、边界条件不变的条件下，将高密度养殖区曼宁数较初始值由20 m1/3·s-1更改为10 m1/3·s-1，低密度养殖区曼宁数较初始值由30 m1/3·s-1更改为20 m1/3·s-1[33]形成更高密度养殖区.

更高密度养殖区与无养殖区流速场进行差值的结果如图4所示.更高密度养殖区条件下养殖区(Ⅰ、Ⅱ)内流速明显减小，对流速的削减量较正常养殖情况下更大，影响范围增加，养殖区外侧出现大范围流速增大的现象.经计算，Ⅰ区近岸处流速减少幅度约为40%～80%，Ⅰ区向海侧流速减少幅度为20%～40%；Ⅱ区流速变化较小，减少幅度为10%～30%；因涨落潮流方向为NE-SW向，沿岸区域也有一定流速减小，减小幅度小于20%；向外海延伸区域，流速增加，增加幅度小于10%.当养殖密度增加时，高密度养殖区流速减少更加明显，而其外部部分海域流速增加也更明显.

图4 平均流速变化场(更高密度养殖区无养殖区)

Fig.4 Distribution of mean flow velocity change caused by the scenario with higher density culture zone and without culture zone

a P1

b P2

c P3

3 筏式养殖对污染物输运的影响

养殖区除影响区域水动力、水质分布外，自身也可成为污染源，尤其在风暴潮等情况下，底质污染物易随水体扰动释放，造成区域污染.假定养殖区污染物初始浓度为单位浓度1 mg·L-1，其他海域初始浓度为0，海区边界及陆地边界条件污染物浓度均设置为0.根据验证率定参数，污染物水平扩散系数取40～120 m2·s-1.模型不指定污染物种类，也不考虑生化降解过程，即降解系数为0.设置了有无养殖区两种工况进行模拟和对比.

将有养殖区浓度场与无养殖区浓度场进行差值，得到浓度变化场如图6所示，图中深色区域表示养殖区导致污染物浓度增高，浅色表示养殖区导致污染物浓度降低.从有无养殖区的平均流速变化场(图3)可以看出，有养殖区条件下，养殖区内的流速明显降低，直接阻碍了养殖区内部污染物的对流扩散能力，所以较无养殖区条件下，养殖区内污染物对流扩散缓慢，污染物浓度增高的区域主要集中在养殖区内.污染物浓度降低的区域主要位于养殖区的南北两侧，这是由于秦皇岛近岸海域潮流总体特征为顺岸往复流，涨潮流为SW向，落潮流为NE向，潮流的方向决定了污染物对流扩散的方向，而养殖区的存在直接降低了污染物向养殖区的SW和NE向扩散的能力，所以有养殖区条件下养殖区南北两侧的浓度较无养殖区条件下有所降低.同时，养殖区作为污染源，本身结构影响了污染物分配，污染物浓度的差值的体现为近海大、远海小，养殖区附近大、其他区域小的特点，导致海域污染物浓度差值由养殖区向外逐渐减小，随着扩散时间的增加，有无养殖区情况污染物浓度差逐渐降低.

根据典型点P1、P2、P3的进行对比有无养殖区污染物浓度变化过程(图7).有无养殖区对污染物对流和扩散有一定的影响，影响程度由近岸到远岸逐渐减小.P1点位因在高密度养殖区内部，有养殖区存在，污染物浓度高于无养殖区情况，有无养殖区污染物浓度差较大，最高达0.14 mg·L-1，差值先增大后减小；P2点位同样受养殖区影响，有养殖区情况污染物浓度高于无养殖区情况，二者差值最高0.05 mg·L-1，差值也是先增大后减小；P3点位在养殖区外，前期污染物浓度基本一致，后期有养殖区情况污染物浓度略高，差值最高0.005 mg·L-1，呈波动性增加.可见养殖区的存在对养殖区内部污染物对流和扩散影响较大，养殖区密度越高，污染物对流扩散越慢.

图6 污染物浓度分布变化过程(养殖区无养殖区)

图7 典型点污染物浓度及浓度差变化过程对比

4 结论

基于二维秦皇岛海域水动力和水质模型，模拟和对比了秦皇岛有无养殖区条件下近岸海域水动力和污染物输运情景，探究养殖区与海域水动力和污染物输运的响应关系，得到如下结论：

(1) 养殖区的存在对海域水动力条件有明显的阻碍作用，养殖区内流速明显减小，养殖区外向外海延伸区域流速略有增大，增加幅度小于10%.

(2) 当养殖密度增加时，高密度养殖区流速减少更加明显，而其外部部分海域流速增加也更明显.高密度养殖区近岸处流速减少幅度约为20%～50%，向海侧流速减少幅度为10%～20%；低密度养殖区流速减少幅度小于10%；向外海延伸区域，流速略增加，增加幅度小于10%.其他区域流速几乎不变.

(3) 养殖区对污染输运有一定影响，养殖区内部影响大，逐渐向外海影响减小.养殖区存在时，污染物对流和扩散较慢，养殖区内浓度高于无养殖区工况，外海区浓度变化不大.

致谢：感谢国家海洋局秦皇岛海洋环境监测中心站为本文的研究提供了充足的实测资料.