秦皇島筏式養殖對水動力和污染物輸運的影響

匡翠萍， 董智超， 顧 杰， 單云馳

(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) 養殖區對污染輸運有一定影響，養殖區內部影響大，逐漸向外海影響減小.養殖區存在時，污染物對流和擴散較慢，養殖區內濃度高于無養殖區工況，外海區濃度變化不大.

致謝：感謝國家海洋局秦皇島海洋環境監測中心站為本文的研究提供了充足的實測資料.