謝中宇，劉涵心，陳界仁

（1.國家海洋局第二海洋研究所工程海洋學(xué)重點(diǎn)實(shí)驗(yàn)室，杭州國海海洋工程勘測設(shè)計(jì)研究院，浙江杭州 310012；2.河海大學(xué)水利水電學(xué)院，江蘇南京 210098）

基于數(shù)學(xué)模型的樂清灣環(huán)境容量算法應(yīng)用

謝中宇1，劉涵心2，陳界仁2

（1.國家海洋局第二海洋研究所工程海洋學(xué)重點(diǎn)實(shí)驗(yàn)室，杭州國海海洋工程勘測設(shè)計(jì)研究院，浙江杭州 310012；2.河海大學(xué)水利水電學(xué)院，江蘇南京 210098）

在實(shí)測資料分析的基礎(chǔ)上建立了樂清灣潮流水環(huán)境數(shù)學(xué)模型,對樂清灣區(qū)域水動(dòng)力環(huán)境及污染物擴(kuò)散輸移進(jìn)行數(shù)值模擬，以COD為代表，計(jì)算分析了樂清灣在大麥嶼港區(qū)規(guī)劃中滿足海洋功能區(qū)劃水質(zhì)標(biāo)準(zhǔn)下的COD環(huán)境容量。結(jié)果表明：樂清灣完成一次水交換時(shí)間為90 d，排污附近海區(qū)完成一次水交換為30 d。大潮期72 h排污口排放污染物COD影響范圍大致為北至樂清灣底，西南至甌江北口，南至洞頭峽，東至鹿西島以東15 km的海域，灣內(nèi)COD濃度為1.231～1.232 mg/ L。允許混合區(qū)的COD剩余環(huán)境容量大潮期為8 148.99～8 164.98 t/a；小潮期為8 144.64～8 162.75 t/a。大麥嶼港區(qū)規(guī)劃實(shí)施后滿足樂清灣COD環(huán)境容量要求。

環(huán)境容量；污染物擴(kuò)散；水體交換；數(shù)學(xué)模型；樂清灣

樂清灣位于浙江省東南部沿海，甌江口北側(cè)，三面環(huán)陸屬典型半封閉性海灣，岸線總長約184.7 km，海灣形似葫蘆，口門寬約21 km，縱深達(dá)40 km，按地形地貌形態(tài)分為內(nèi)灣、中灣和外灣三大部分[1]，岸線以下海灣總面積463.6 km2。樂清灣海區(qū)屬半封閉性強(qiáng)潮海灣，控制本海區(qū)潮波運(yùn)動(dòng)的是以M2分潮為主的東海前進(jìn)潮波系統(tǒng)，潮波由洞頭洋經(jīng)灣口諸海峽傳入灣內(nèi)。根據(jù)海域?qū)崪y資料統(tǒng)計(jì)，灣口海域最高、最低潮位分別為4.30 m和-3.81 m；平均高、低潮位分別為2.77 m和-2.25 m，量值與坎門站比較接近。

由于土地資源開發(fā)、圍塘養(yǎng)殖、航運(yùn)等需要,臺(tái)州港大麥嶼港區(qū)規(guī)劃于玉環(huán)半島西側(cè)岸段建設(shè)包括普竹、連嶼、大巖頭在內(nèi)的三個(gè)圍填海作業(yè)區(qū)及附屬碼頭，填海面積375.1 hm2，新建碼頭泊位28個(gè)。為研究該工程后樂清灣區(qū)域特定污染水的環(huán)境容量，利用計(jì)算分析結(jié)果科學(xué)分配和控制排放總量，本文基于MIKE21數(shù)學(xué)模型，建立了包括樂清灣、隘頑灣、溫州灣在內(nèi)的大范圍水動(dòng)力模型，利用實(shí)測水文觀測資料對模型進(jìn)行驗(yàn)證，在此基礎(chǔ)上利用資料分析和數(shù)值模擬，了解樂清灣的水動(dòng)力環(huán)境。其次，在水動(dòng)力模型基礎(chǔ)上建立區(qū)域水交換模式，并計(jì)算碼頭工程區(qū)污染物擴(kuò)散過程及分析影響，從而估算區(qū)域內(nèi)環(huán)境容量。

1 潮流數(shù)學(xué)模型的建立及驗(yàn)證

采用丹麥水力研究所MIKE 21（2007）模型建立目標(biāo)海域的潮流數(shù)學(xué)模型，運(yùn)用實(shí)測潮位、流速資料進(jìn)行模型驗(yàn)證。

1.1 潮流數(shù)學(xué)模型

數(shù)學(xué)模型控制方程為：

式中：ζ為潮位，h為水深，p、q分別為x、y方向上的垂線平均單寬流量，C為謝才系數(shù)為糙率系數(shù)，H=h+ζ，n=0.022，g為重力加速度g=9.81 m/s2，f為風(fēng)應(yīng)力系數(shù)，V、Vx、Vy分別為風(fēng)速及其在x、y方向的分量，柯氏力Ω=2 wsinφ，w為地轉(zhuǎn)角速度，φ為緯度取28.1°(計(jì)算區(qū)域平均維度)，pa為大氣壓力，ρw為海水密度，x、y為直角坐標(biāo)，t為時(shí)間，τxx、τxy、τyy分別為剪切應(yīng)力分量，u、v為水平流速分量，紊動(dòng)粘性系數(shù)采用Smagorinsky公式計(jì)算

為求解上述控制方程，必須要有合適的定解條件，所需的定解條件包括初始條件和邊界條件。

在海域中的開邊界上，海面的水位根據(jù)沿岸潮位站來給定。MIKE 21 HD采用交替方向隱格式(ADI)求解二維淺水潮波方程，方程矩陣采用雙消除法（Double Sweep）算法求解，該格式具有二階精度[2]。

模型對碼頭群樁的處理采用MIKE21 Pier結(jié)構(gòu)，其原理為：

式中：τp為等價(jià)剪切應(yīng)力，F(xiàn)為單個(gè)橋樁的拖曳力，n為橋樁的密度，Δx、Δy為網(wǎng)格間距。

總體拖曳力為群樁阻力和底摩擦力之和：

式中：CD為拖曳系數(shù)，ρ為海水密度，Be為橋樁的等效寬度，He為橋樁在水中的高度，V為流速。

1.2 模型配置

模擬排污口位于樂清灣口門，潮動(dòng)力環(huán)境較強(qiáng)，為了能更好地模擬樂清灣大范圍及局部潮流場情況，采用網(wǎng)格實(shí)時(shí)耦合的方法：耦合計(jì)算區(qū)分為大、中、小三層同步計(jì)算的區(qū)域，具體情況如圖1，網(wǎng)格的Y向均為正北向向東偏轉(zhuǎn)18°，網(wǎng)格精度分別為270 m、90 m、30 m，時(shí)間步長10 s，共計(jì)矩形網(wǎng)格數(shù)745 800個(gè)。

計(jì)算區(qū)域的大區(qū)北起椒江南岸-大陳島一線、南至鰲江南岸-南麂列島，計(jì)算大區(qū)開邊界上的水位，由浪璣山、下大陳、南麂和石坪4個(gè)潮位站的水位得到，潮位站具體位置如圖1。開邊界處的水位同時(shí)包含11個(gè)主要分潮（K1、O1、P1、Q1、M2、S2、N2、K2、M4、MS4和M6）。溫州三大河口的內(nèi)邊界由石仙婦、瑞安、鰲江三個(gè)岸站的水位控制（含11個(gè)分潮）；另外考慮到灣內(nèi)徑流對模型精度的影響，在模型中清江、楚門河等樂清灣主要河流入灣口處加強(qiáng)徑流源，流量0.032 6 m3/s[7]，模擬期內(nèi)不間斷釋放。樂海碼頭、坎門和洞頭為3個(gè)潮位驗(yàn)證站，S1-8為8個(gè)潮流驗(yàn)證站（國家海洋局第二海洋研究所為甬臺(tái)溫復(fù)線樂清灣大橋工程實(shí)測水文站位，觀測內(nèi)容包括流速、流向、水深和懸沙濃度，觀測時(shí)間為26～28 h[3]）。模型涉及地形主要采用本海區(qū)歷史海圖（13640、13710、13715、13741、13770、13771、13781）及業(yè)主單位提供的工程前沿測深圖。

圖1 模型配置圖Fig.1 The model configuration

1.3 潮位、流速驗(yàn)證

圖2給出了潮位驗(yàn)證曲線，由圖可見，模擬所得的三個(gè)潮位與觀測結(jié)果吻合良好，高、低潮位誤差基本小于10 cm，相位基本一致，其中的差異一方面可能是由網(wǎng)格離散點(diǎn)與潮位站之間的位置差異所造成的，另一方面，可能是由于實(shí)際水深與離散網(wǎng)格水深取值的差異所致。

圖2 潮位驗(yàn)證圖Fig.2 Tidal level validation

分析計(jì)算與實(shí)測潮流流速可知，模擬所得的流場基本與實(shí)際情況相符，計(jì)算值和實(shí)測流速相對誤差小于15%；總體看來，模擬所得的8個(gè)連續(xù)站的潮流與實(shí)測流吻合得較好，區(qū)域內(nèi)流場模擬計(jì)算結(jié)果基本反映了該海域潮流和潮波的實(shí)際變化，流場模擬驗(yàn)證基本達(dá)到了后續(xù)預(yù)測計(jì)算的要求[4]。

2 水環(huán)境數(shù)值模擬及環(huán)境容量估算

為分析排污對海域水環(huán)境的影響，在潮流數(shù)學(xué)模型基礎(chǔ)上建立污染物擴(kuò)散數(shù)學(xué)模型，計(jì)算樂清灣水體交換變化、污染物擴(kuò)散特性，以COD污染物為水質(zhì)預(yù)測因子，計(jì)算分析排放典型污染物對海域水質(zhì)、水環(huán)境容量的影響。

圖3 大、小潮潮流流速、流向驗(yàn)證圖Fig.3 Tidal current velocity and direction validation of spring-neap tide

2.1 海灣水交換及水體更新周期

水交換狀況對環(huán)境容量有著極重要的影響，既代表了外來物質(zhì)的輸入能力，也反映了灣內(nèi)物質(zhì)的稀釋能力。為模擬樂清灣的水交換狀況，先假設(shè)灣內(nèi)均勻分布著濃度為1的溶解態(tài)保守性物質(zhì)，外海水的保守物質(zhì)濃度為0，先運(yùn)行潮波運(yùn)動(dòng)模式，待潮流穩(wěn)定后開始運(yùn)行溶解態(tài)保守物質(zhì)的對流-擴(kuò)散模式，因而在某一時(shí)刻灣內(nèi)保守物質(zhì)的濃度c可代表灣內(nèi)水經(jīng)過一段時(shí)間對外擴(kuò)散后剩余的濃度，即1-c可表示灣內(nèi)水體被外海置換比率。

以溶解態(tài)的保守性物質(zhì)作為灣內(nèi)水的示蹤劑，建立樂清灣海域?qū)α?擴(kuò)散型的水交換數(shù)值模式[5]，灣內(nèi)水示蹤劑的控制方程為：

其中，(u，v)深度平均流速在笛卡爾坐標(biāo)(x，y)方向的分量；c保守性溶解態(tài)灣內(nèi)水的示蹤劑濃度：k=kh+ks+kt，k擴(kuò)散系數(shù)，kh垂向結(jié)構(gòu)的余環(huán)流引起的水平輸運(yùn)的擴(kuò)散系數(shù)，ks垂向剪切引起的水平輸運(yùn)的擴(kuò)散系數(shù)，kt紊流擴(kuò)散系數(shù)，這里均取模型默認(rèn)穩(wěn)態(tài)系數(shù)值。需要指出的是徑流對于封閉海灣水交換具有一定的促進(jìn)作用，樂清灣流域有山溪性河流30余條，主要有大荊溪、清江、楚門河等，多年平均徑流量約103 000 m3。因此與水動(dòng)力模型相同的，在模型中清江、楚門河等主要河流入灣口處加徑流源強(qiáng)，流量0.032 6 m3/s，模擬期內(nèi)不間斷釋放，徑流邊界保守物質(zhì)濃度為0。

運(yùn)用建立的樂清灣海域污染物對流-擴(kuò)散數(shù)值模型模擬水體交換過程，根據(jù)運(yùn)行10～90 d后樂清灣表層保守物質(zhì)濃度分布情況可知，受半封閉狹長地形影響，從灣口向?yàn)车?，樂清灣水域中的水交換由快變慢，且差別很大，具體模擬灣內(nèi)各區(qū)域的海水置換率變化情況見表1。一般認(rèn)為海域水體交換置換比率達(dá)到70%即基本完成一次水交換，根據(jù)水交換模擬結(jié)果，排污口所在的玉環(huán)半島西側(cè)海域，完成一次水交換的時(shí)間不到30 d。模擬結(jié)果如圖4所示。

圖4 水交換模擬30 d保守物質(zhì)分布Fig.4 Distribution of conservative material modelled in 30 days

表1 模擬灣內(nèi)各區(qū)域海水置換率變化（%）Tab.1 The replacement rate of the sea area(%)

2.2 海域污染物擴(kuò)散數(shù)值模擬

污水?dāng)U散模擬采用Mike21AD（2007）模型進(jìn)行計(jì)算，模擬預(yù)設(shè)排污口位于樂清灣口門，距岸邊250 m，不間斷連續(xù)排放，根據(jù)港區(qū)控制性詳規(guī)[8]擬定排污口日排放量為8 000 m3/d。污水?dāng)U散計(jì)算以COD作為水質(zhì)預(yù)測評價(jià)因子研究所排廢水對周邊環(huán)境的影響，排放源強(qiáng)為30 mg/L，在污染物總量不變的前提下，排污口源強(qiáng)統(tǒng)一按流量0.1 m3/s加入。根據(jù)實(shí)測資料得到海域COD本底值分布狀況，計(jì)算可得COD平均本底值為1.23 mg/L，在計(jì)算時(shí)將海域COD本底值分布狀況加入模型作為污水?dāng)U散計(jì)算初始濃度。同時(shí)將上述排污口源強(qiáng)加入模型連續(xù)運(yùn)行365 d，得到將一年的污染物全部排放以后的海域狀態(tài)，再分別進(jìn)行大潮和小潮期污染物擴(kuò)散過程的模擬，并分析區(qū)域內(nèi)COD污水最大增量分布情況，可知：

（1）大潮期：72 h內(nèi)COD污水影響區(qū)域大致為北至樂清灣底，西南至甌江北口，南至洞頭峽，東至鹿西島以東15km海域，整個(gè)樂清灣內(nèi)的污水濃度在1.231～1.232 mg/L之間，其中江巖山以北的區(qū)域污水濃度均大于1.231 5 mg/L。

（2）小潮期：72 h內(nèi)COD污水影響區(qū)域大致為北至樂清灣底，西南至甌江北口，南至洞頭島西岸，東至鹿西島以東7 km海域，整個(gè)樂清灣內(nèi)的污水濃度在1.231～1.232 5 mg/L之間，其中玉環(huán)電廠廠址以北的區(qū)域污水濃度均大于1.232 mg/L。

關(guān)于排污口初始稀釋度的預(yù)測，根據(jù)潮汐潮流數(shù)值模擬的結(jié)果，排污口附近流速取大潮平均流速為0.81 m/ s，取排污口流量為8 000 m3/d、單管排放、21.1 m的平均水深及海水密度1.025 t/m3，計(jì)算得出排污口初始稀釋度約為1 241.5，可見由于排污口處灣口主槽內(nèi)，水深流急，排污口附近海域污水?dāng)U散條件優(yōu)良。

圖5 大潮期COD污水最大增量分布Fig.5 Maximum increment of COD in Spring tide period

表2 污水最大濃度增量等值線包絡(luò)面積（km2）Tab.2 The envelope area of the maximum concentration increment of sewage

2.3 工程海域環(huán)境容量估算

根據(jù)《中國環(huán)境保護(hù)標(biāo)準(zhǔn)匯編-環(huán)境質(zhì)量與污染物排放2003》中有關(guān)污水海洋處置工程污染控制標(biāo)準(zhǔn)（GB18486-2001）的章節(jié)[6]，污水海洋處置工程污染物的混合區(qū)規(guī)定：若污水排往開敞海域或面積≥600 km2(以理論深度基準(zhǔn)面為準(zhǔn))的海灣及廣闊河口，允許混合區(qū)范圍：Aa≤3.0 km2。若污水排往面積＜600 km2的海灣，混合區(qū)面積必須取按以下兩種方法計(jì)算所得值（Aa）的小者：

式中L為擴(kuò)展器長度,m。

式中A0為計(jì)算至灣口位置的海灣面積，m2。

對于重點(diǎn)海域和敏感海域，劃分污水海洋處置工程污染物的混合區(qū)時(shí)還需要考慮排放點(diǎn)所在海域的水流交換條件、海洋水生生態(tài)等。

根據(jù)《中國海灣志》記載，樂清灣岸線以下海灣總面積為463.6 km2，其中灘涂面積（岸線至理論深度基準(zhǔn)面)約為220.8 km2，約占海灣總面積48%；水域面積（理論深度基準(zhǔn)面以下）約為242.8 km2，占海域總面積52%左右，平均水深10 m左右[7]。因此按照規(guī)定根據(jù)式（7）計(jì)算，得到允許混合區(qū)面積Aa為0.6 km2。

由溫州市海洋環(huán)境功能區(qū)劃，玉環(huán)半島連嶼至橫址山前沿為四類環(huán)境功能區(qū)，樂清灣內(nèi)的其他水域?yàn)槎惌h(huán)境功能區(qū)，考慮到污染物擴(kuò)散模擬的結(jié)果，污水在牛頭頸排海后將影響包括灣底養(yǎng)殖區(qū)域在內(nèi)的樂清灣內(nèi)大部分海域，因此本規(guī)劃按照二類水質(zhì)標(biāo)準(zhǔn)進(jìn)行環(huán)境容量計(jì)算。二類水質(zhì)標(biāo)準(zhǔn)為COD≤3.0 mg/L，根據(jù)評價(jià)海域的水環(huán)境質(zhì)量現(xiàn)狀監(jiān)測結(jié)果，COD的海域平均值為1.23 mg/L，將目標(biāo)濃度減去水質(zhì)現(xiàn)狀值，在二類海域內(nèi)排污區(qū)的COD污染物允許濃度增量為1.77 mg/L。

由上述允許混合區(qū)面積，將其概化成以排污口為圓心，面積0.6 km2，半徑437 m的圓型區(qū)域，討論排污口COD排放濃度30 mg/L連續(xù)排放的情況下該區(qū)域內(nèi)的環(huán)境容量。根據(jù)水體交換過程模擬的結(jié)果，允許混合區(qū)所在海域完成一次水交換的時(shí)間小于30 d；污染物擴(kuò)散模式運(yùn)算的起算點(diǎn)為水動(dòng)力模塊基礎(chǔ)上加入海域平均值和排污口源強(qiáng)連續(xù)運(yùn)行1 a后的海域狀態(tài)，因此可以認(rèn)為，起算時(shí)海域內(nèi)污染物濃度達(dá)到動(dòng)態(tài)平衡狀態(tài)。統(tǒng)計(jì)污染物擴(kuò)散模式的計(jì)算結(jié)果，排污管道建成運(yùn)營后，允許混合區(qū)內(nèi)的COD大潮期平均濃度增量為0.003 mg/L，最大濃度增量為0.006 5 mg/L；小潮期平均濃度增量為0.003 5 mg/L，最大濃度增量為0.007 4 mg/L，將污染物允許濃度增量乘以允許混合區(qū)的水體體積，可以得到區(qū)域內(nèi)的環(huán)境容量見表3。這里由排污口流量和污水濃度得排污管排海COD污水87.6 t/a，通過計(jì)算可以得到本規(guī)劃排污對允許混合區(qū)環(huán)境容量的貢獻(xiàn)值。

表3 剩余環(huán)境容量估算及污水濃度增量對其貢獻(xiàn)率Tab.3 Estimation of residual environmental capacity and the contribution rate of sewage concentration increment

3 結(jié)語

海洋環(huán)境容量作為海洋環(huán)境管理的重要手段，對其計(jì)算能夠?qū)Ｑ笏|(zhì)環(huán)境安全保護(hù)提供參考依據(jù)。本文以樂清灣大麥嶼港區(qū)規(guī)劃污染物排放項(xiàng)目為例，通過二維水質(zhì)數(shù)學(xué)模型計(jì)算該項(xiàng)目對區(qū)域COD環(huán)境容量的貢獻(xiàn)值，為后期水質(zhì)環(huán)境管理工作及同類型海洋環(huán)境容量計(jì)算提供技術(shù)借鑒。

排污項(xiàng)目所在樂清灣為浙東南典型半封閉強(qiáng)潮流海灣，通過水體交換模擬可知灣口處的水體在潮汐作用下的交換周期不到30 d。灣口處排放污染物擴(kuò)散條件優(yōu)良，影響范圍基本包括整個(gè)樂清灣海域，濃度增量除項(xiàng)目前沿均較小。根據(jù)水質(zhì)模型的預(yù)測結(jié)果，結(jié)合項(xiàng)目所處海洋功能區(qū)劃的水質(zhì)標(biāo)準(zhǔn)，計(jì)算得到本項(xiàng)目允許混合區(qū)的COD污水剩余環(huán)境容量大潮期為8 148.99～8 164.98 t/a；小潮期為8 144.64～8 162.75 t/a。

[1]黃秀清.樂清灣海洋環(huán)境容量及污染物總量控制研究[M].北京:海洋出版社,2011:16-17.

[2]李 峰,謝中宇.臺(tái)州港大麥嶼港區(qū)規(guī)劃數(shù)值模擬計(jì)算及潮間帶生物現(xiàn)狀調(diào)查報(bào)告[R].2012.

[3]國家海洋局第二海洋研究所.浙江省甬臺(tái)溫高速公路復(fù)線樂清灣大橋及接線工程數(shù)學(xué)模型試驗(yàn)報(bào)告[R]//海洋水文及工程可行性專題研究.2011.

[4]中華人民共和國交通部.JTJ/T 233-98海岸與河口潮流泥沙模擬技術(shù)規(guī)程[S].1998.

[5]董禮先,蘇紀(jì)蘭.象山港水交換數(shù)值研究Ⅱ.模型應(yīng)用和水交換研究[J].海洋與湖沼,1999,30(5):465-470.

[6]國家環(huán)境保護(hù)總局,國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局.GB 18486-2001污水海洋處置工程污染控制標(biāo)準(zhǔn)[S].2001.

[7]謝欽春,馮應(yīng)俊,等.中國海灣志第六分冊浙江省南部海灣[M].北京:海洋出版社,1993.

[8]交通部規(guī)劃研究院.臺(tái)州港大麥嶼港區(qū)控制性詳細(xì)規(guī)劃報(bào)告[R]//工程可行性專題研究.2006.

Application of Mathematical Model in the Yueqing Bay Environmental Capacity Algorithm

SHAO Zhu-feng，LIU Han-xin2，CHEN Jie-ren2
(1.The Second Institute of Oceanography,Soa,Hangzhou 310012;2.College of Water Conservancy and Hydropower Engineering,Hehai University,Nanjing 210098,China)

Based on the analysis of the measured data,the hydrodynamic and environment model is established,according to the numerical model of Yueqing bay area hydrodynamic environment and pollutant transport,calculation and analysis of the environmental capacity of typical bays have been carried out in the marine functional zoning.The calculated results is shown that:the water exchange time is 90 days in Yueqing bay, the water exchange time is 30 days near the sewage sea.The COD concentration is 1.231-1.232 mg/L at the inner bay in the tide period.The residual environmental capacity of COD is 8 144.64-8 162.75 t/ain spring tide period.The residual environmental capacity of COD is 8 148.99-8 164.98 t/ain neap tide period.The environmental capacity of Yueqing bays meet the requirements after the introduction of Damaiyu port planning.

environmental capacity;pollutant dispersion;water exchange;mathematical model;Yueqing bay

X26

A

1008－830X(2016)05－0430－06

2016-07-19

謝中宇（1983-），男，浙江蒼南人，工程師，研究方向：海洋工程數(shù)學(xué)建模與環(huán)境風(fēng)險(xiǎn).E-mail:xiezhongyu23@163.com