鮑思佳，湯立君，陶平，于彩虹，宋成文，邵秘華

（大連海事大學(xué)環(huán)境科學(xué)與工程學(xué)院，遼寧 大連 116026）

1 引言

環(huán)境容量是在特定時空內(nèi)有效可交換水體的同化能力，以及在規(guī)定環(huán)境目標(biāo)下所能容納污染物的量。換言之，排入海域污染物濃度不超過環(huán)境管理標(biāo)準(zhǔn)的前提下所能允許的最大排放量，是海洋環(huán)境容納污染物的最大容量[1]。河口海域水體污染物成分往往是人為因素影響在自然本底值上的加合，這種加合超過一定限度時，就會改變環(huán)境質(zhì)量[2-3]。因此，環(huán)境容量是環(huán)境質(zhì)量評價、環(huán)境監(jiān)測和環(huán)境承載力預(yù)警的基礎(chǔ)，也是河口海域資源開發(fā)利用的環(huán)境依據(jù)。關(guān)于環(huán)境容量的研究，目前國內(nèi)外學(xué)者已進(jìn)行了很多工作，海域環(huán)境容量的研究也越來越受到重視[4]。

環(huán)境容量值的確定有多種方法。20世紀(jì)80年代末一般采用濃度測定方法，21世紀(jì)隨著系統(tǒng)網(wǎng)絡(luò)應(yīng)用技術(shù)的發(fā)展以及對海洋環(huán)境容量計算研究的不斷深入，其計算方法主要有箱式模型、分擔(dān)率法和最優(yōu)化法[5]。特別是隨著水環(huán)境承載力評估的發(fā)展，從初期的局部河流水環(huán)境范圍內(nèi)的污染物自凈能力研究，再到湖泊、河口和海灣水域的污染物控制政策的實(shí)施，水環(huán)境承載力的研究逐漸受到關(guān)注，同時也發(fā)展了海洋環(huán)境承載力評估研究[6-8]，并建立了海域環(huán)境承載力評估量化方法，如指數(shù)評價法、生態(tài)足跡法、系統(tǒng)動力學(xué)方法和統(tǒng)計學(xué)動態(tài)模型等[9]。

河口作為海洋和陸地的交界，港口和航運(yùn)等經(jīng)濟(jì)活動頻繁，生態(tài)系統(tǒng)脆弱。大遼河口位于遼東灣北部海域，為渾河和太子河匯流后在遼東灣的入?？趨^(qū)域，是典型的三角洲河口[10]。由于大遼河口得天獨(dú)厚的豐富資源，因此成為遼寧省海岸帶綜合開發(fā)的國家重要節(jié)點(diǎn)。但大遼河口沿岸企業(yè)入海污染嚴(yán)重，同時匯集了遼寧省6個城市的生活污水，廢水廢物攜帶量位居全省河流之首[11]，環(huán)境治理迫在眉睫。本文選取大遼河口海域?yàn)檠芯繉ο箝_展環(huán)境承載力研究。劉娟等[12]指出大遼河口及河口鄰近海域的主要污染物是化學(xué)需氧量（Chemical Oxygen Demand，COD）、無機(jī)氮、磷酸鹽和重金屬等；柴寧[13]也對大遼河水系主要污染物COD和氨氮多年時空變化特征進(jìn)行了研究，認(rèn)為兩種污染物主要來源于工業(yè)和城市的點(diǎn)源排放；而周丹卉[14]在對大遼河口典型污染物時空分布的研究中發(fā)現(xiàn)，大遼河口水體的有機(jī)物和氮污染最為嚴(yán)重，且季節(jié)差異明顯。因此，本文選取具有代表性的3種污染要素：COD、無機(jī)氮和磷酸鹽作為大遼河口海域的主要污染要素開展海域環(huán)境承載力預(yù)警技術(shù)研究。

綜上所述，本文以2019年大遼河口海域水環(huán)境調(diào)查資料為基礎(chǔ)，將響應(yīng)系數(shù)改為響應(yīng)濃度函數(shù)，采用反演迭代法計算該海域3種主要污染物的環(huán)境容量，再依據(jù)海洋功能區(qū)劃水質(zhì)標(biāo)準(zhǔn)對大遼河口海域環(huán)境承載力進(jìn)行評估和預(yù)警，為該海域的海洋管理部門制定海域污染防治政策提供科學(xué)技術(shù)支持。

2 計算方法

2.1 數(shù)值模型構(gòu)建

2.1.1 水質(zhì)模型

鑒于本文研究的海域?yàn)榻稖\水區(qū)，表層和底層物理、化學(xué)和生物性質(zhì)差別較小，水平方向和垂直方向上的變化可忽略[15]，因此采用深度平均。二維流體動力學(xué)方程和水質(zhì)控制模型方程如下：

式中：ζ為平均海平面以上的瞬時水位高度，單位：m；H為水深，單位：m；Q為徑流量；g為重力加速度，單位：m/s2；U和V分別為垂直平均流速在x和y軸上的分量，單位：m/s；f為柯氏參數(shù)，f=2ωsinφ，ω為地轉(zhuǎn)角速度，φ為地理緯度；C為謝才系數(shù)，C=；n為曼寧系數(shù)（取0.025～0.035）；t為時間；x和y平面取在未擾動的平均海平面上，z軸垂直向上，構(gòu)成右手坐標(biāo)系（見圖1）。式（1）—（3）的初始條件從靜止水狀態(tài)開始，U=V=ζ=0。邊界條件分兩類：沿岸閉邊界，取法向流速等于零（Vn=0）；開邊界各點(diǎn)水位為時間的已知函數(shù)，即ζ=ζ(t)。

圖1 計算坐標(biāo)系

2.1.2 差分方程

有限差分方法是計算機(jī)數(shù)值將水質(zhì)控制模型方程中的導(dǎo)數(shù)用網(wǎng)格節(jié)點(diǎn)上的函數(shù)值的差商代替進(jìn)行離散，從而求解以網(wǎng)格節(jié)點(diǎn)上的值為未知數(shù)的代數(shù)方程組。

本文采用顯隱方向交替（Alternating Direction Implicit，ADI）計算網(wǎng)格（見圖2），將式（1）—（3）離散成相應(yīng)的差分方程。在nΔt→（n+1/2）Δt時間半步長上，將式（1）和式（2）按行進(jìn)行隱式差分格式離散，V按顯示差分格式離散。對后半時間步長（n+1/2）Δt→（n+1）Δt，對式（1）和式（3）按列進(jìn)行隱式差分格式離散，U按顯示差分格式離散，進(jìn)而建立大遼河口海域水質(zhì)數(shù)值計算模型。

圖2 差分網(wǎng)格圖

2.1.3 大遼河口海域水動力模型驗(yàn)證

依據(jù)上述水動力模型，在2019年5月對大遼河口海域3#污染源進(jìn)行潮位和流速流向的現(xiàn)場觀測（見圖3），并與模型計算值對比，得到潮位和流速流向?qū)崪y值與計算值對比結(jié)果（見圖4和圖5）。由圖可見，在一定的誤差范圍內(nèi)，該水動力模型在遼東灣近岸淺水海區(qū)的潮位和流速流向計算值與實(shí)際觀測值基本吻合，可作為海區(qū)環(huán)境容量計算的基礎(chǔ)。

圖3 模擬污染源位置圖

圖4 大遼河口海域潮位驗(yàn)證圖

圖5 大遼河口海域流速流向驗(yàn)證圖

2.1.4 保守物質(zhì)擴(kuò)散方程

保守物質(zhì)是在輸移擴(kuò)散過程中沒有物質(zhì)轉(zhuǎn)化和形態(tài)改變的物質(zhì)，擴(kuò)散方程中無源項(xiàng)和匯項(xiàng)，方程僅考慮對流和擴(kuò)散作用，只需要確定邊界和物質(zhì)擴(kuò)散系數(shù)[16]。由海區(qū)功能區(qū)劃圖可看出（見圖6），計算海域位于近岸淺水區(qū)，水淺垂直混合充分，可將垂直方向上的水體視作均勻不變，故本文采用保守物質(zhì)擴(kuò)散模型，即二維垂直平均擴(kuò)散方程進(jìn)行計算：

圖6 大遼河口海域功能規(guī)劃網(wǎng)格圖

式中：P為擴(kuò)散物質(zhì)的深度平均濃度；p=；Sm為污染物源強(qiáng)；Kx和Ky為擴(kuò)散系數(shù)，由Elder公式確定：

2.2 污染源輸入負(fù)荷量-輸出響應(yīng)濃度線性系統(tǒng)驗(yàn)證

數(shù)學(xué)物理方程主要有3種類型：橢圓型、拋物型和雙曲型。有學(xué)者認(rèn)為污染物在海水中的擴(kuò)散方程是橢圓型的線性方程，且擴(kuò)散方程在第一邊界條件下是線性的，符合疊加原理[17]。為了驗(yàn)證在環(huán)境容量的計算中，擴(kuò)散方程是橢圓線性的數(shù)學(xué)物理方程，本文將通過以下擴(kuò)散數(shù)值模擬計算出污染源輸入負(fù)荷量-輸出響應(yīng)濃度場是線性系統(tǒng)。

以大遼河口海域?yàn)槔?，從北至南依次設(shè)置5個模擬污染源，分別為1#、2#、3#、4#和5#（見圖3）。利用保守物質(zhì)擴(kuò)散方程的數(shù)值模型對1#—5#污染源各自做輸入負(fù)荷量-輸出響應(yīng)濃度場模擬計算，使用的各項(xiàng)參數(shù)如下：

（1）計算其中任意一個污染源的輸入負(fù)荷量-輸出響應(yīng)濃度時，要計算的污染源輸入負(fù)荷量分別取“1 t/d”、“2 t/d”和“5 t/d”，其他污染源輸入負(fù)荷量均為“0 t/d”；

（2）開邊界濃度為“0 mg/L”（增量計算）；

（3）計算30個潮周期（濃度場已穩(wěn)定）結(jié)束。

計算結(jié)果見表1。

表1 輸入負(fù)荷量為1～5 t/d時5個污染源排污口處輸出響應(yīng)濃度（單位：mg/L）

表1中行數(shù)據(jù)表示各污染源單獨(dú)輸入1 t/d、2 t/d或5 t/d（其他污染源輸入均為0 t/d）時，各污染源排污口處輸出響應(yīng)濃度和受影響濃度值；列數(shù)據(jù)表示各污染源單獨(dú)輸入1 t/d、2 t/d或5 t/d（其他污染源輸入0 t/d）時產(chǎn)生的輸出響應(yīng)濃度，以及受其他污染源輸入1 t/d、2 t/d或5 t/d時產(chǎn)生的影響濃度值。

用表1中數(shù)據(jù)對1#—5#污染源輸入負(fù)荷量和輸出響應(yīng)濃度作相關(guān)分析，得到回歸方程：

式中：Si為第i個污染源輸入負(fù)荷量，單位：t/d；Ci為第i個污染源排污口處輸出響應(yīng)濃度，單位：mg/L；KCi為回歸系數(shù)；bCi為回歸常數(shù)。

式中：Ci為輸出響應(yīng)濃度，單位：mg/L；Si為污染源輸入負(fù)荷量，單位：t/d；KSi為回歸系數(shù)；bSi為回歸常數(shù)。

式（6）和式（7）是由坐標(biāo)原點(diǎn)引出的斜率為KCi和KSi的直線方程，bCi和bSi值很小，接近“0”，則式（6）和式（7）可簡化為：

將表1中各污染源輸出形成的響應(yīng)濃度代入式（6）中，計算1#—5#污染源輸入的負(fù)荷量，同給定的負(fù)荷量值進(jìn)行比較，結(jié)果見表2。

從表2可看出，由相關(guān)方程計算的各污染源負(fù)荷量和給定的負(fù)荷量只相差±0.002，這一差值在誤差范圍內(nèi)。

表2 計算值與給定值比較表

綜上所述，通過擴(kuò)散方程的數(shù)值模擬計算結(jié)果和相關(guān)分析，以及式（8）和式（9）計算得到的負(fù)荷量值與給定值在誤差允許的范圍內(nèi)相等，證明了污染源輸入負(fù)荷量-輸出響應(yīng)濃度是線性系統(tǒng)。

2.3 建立污染源排放量-輸出響應(yīng)濃度函數(shù)的兩個相關(guān)方程

根據(jù)已建立的水質(zhì)模型計算大遼河口海域各污染源不同單位排放量時的輸出響應(yīng)濃度函數(shù)值（一個潮周期各正點(diǎn)濃度），再建立各污染源排放量-輸出響應(yīng)濃度和輸出響應(yīng)濃度-污染源排放量的兩個相關(guān)方程，即：

式中：Si為第i個污染源污染物的負(fù)荷量，單位：t/d；KCi為第i個方程回歸系數(shù)；Ci為第i個控制點(diǎn)處污染物響應(yīng)濃度，單位：mg/L。

式中：Ci為第i個控制點(diǎn)處的污染物響應(yīng)濃度，單位：mg/L；KSi為第i個方程回歸系數(shù)；Si為第i個污染源污染物的負(fù)荷量，單位：t/d。

因此，水質(zhì)模型是線性系統(tǒng)方程，在初始濃度和邊界濃度均為0的條件下，污染源負(fù)荷量和輸出響應(yīng)濃度成線性關(guān)系，而且是由原點(diǎn)引出的一條斜率為KCi（KSi）的直線，符合疊加原理。

若預(yù)警海域有n個污染源，則式（10）化為：

式中：Si為第i個污染源污染物的負(fù)荷量，單位：t/d；KCj為第j個污染源單位輸出負(fù)荷量對第j個污染源控制點(diǎn)產(chǎn)生的響應(yīng)濃度系數(shù)；cj為第j個污染源輸出負(fù)荷量產(chǎn)生的響應(yīng)濃度；，是（n-1）個污染源對第n個污染源產(chǎn)生的影響濃度總和。

因此，可以根據(jù)各污染源污染要素的響應(yīng)濃度推算出該海區(qū)污染物的負(fù)荷量，這是建立計算海區(qū)負(fù)荷量方程的基礎(chǔ)。

2.4 反演迭代法

目前，常見的環(huán)境容量計算方法是基于線性規(guī)劃的排海通量最優(yōu)化法（簡稱線性規(guī)劃法）[18-19]。該方法能夠應(yīng)用于污染情況復(fù)雜的海域，但計算過程繁瑣。為了減少計算工作量和避免污染源間負(fù)荷量調(diào)配的過程，基于上文構(gòu)建的兩個可逆方程，本文采用“反演迭代法”。

反演方法在海洋學(xué)中的首次應(yīng)用是在水文潮流場上，它是在已知或假定的動力學(xué)約束條件下進(jìn)行統(tǒng)計推斷的一般方法[20]。本文中，反演法是指已知方程Si=KCiCi中濃度Ci的值，倒推算相應(yīng)污染源的負(fù)荷量Si的值；若有多個污染源，則各污染源相應(yīng)控制點(diǎn)處的輸出響應(yīng)濃度互受影響，是一個動態(tài)趨于穩(wěn)定的過程，可以通過線性系統(tǒng)的疊加原理來完成計算。上述兩種計算過程合稱為“反演迭代法”。

本文根據(jù)大遼河口海域國家功能區(qū)劃水質(zhì)標(biāo)準(zhǔn)，通過反演迭代法即可算得該功能區(qū)污染源的標(biāo)準(zhǔn)環(huán)境容量；再將該功能區(qū)控制點(diǎn)處現(xiàn)場水質(zhì)（COD、無機(jī)氮和磷酸鹽）觀測濃度資料代入方程Si=KCiCi中，即可算得該功能區(qū)污染源的水質(zhì)現(xiàn)狀環(huán)境容量。各污染源的標(biāo)準(zhǔn)環(huán)境容量和現(xiàn)狀環(huán)境容量，是評估海域環(huán)境承載力等級的關(guān)鍵數(shù)據(jù)。

2.5 有效濃度

在環(huán)境容量計算中，根據(jù)《水質(zhì)達(dá)標(biāo)方案技術(shù)指南》，一個潮周期內(nèi)響應(yīng)濃度的取值有兩種方法：一種是取一個潮周期內(nèi)最高濃度和最低濃度的平均值；另一種是取一個潮周期各潮時的響應(yīng)濃度的平均值。但污染源輸入-輸出響應(yīng)濃度場在一個潮周期內(nèi)隨潮時而變化，即不同潮時的響應(yīng)濃度也不同。潮汐規(guī)律不同的海區(qū)，響應(yīng)濃度隨潮時的變化規(guī)律也不同。根據(jù)實(shí)測資料可知，大遼河口海域是不規(guī)則的半日潮海區(qū)[21]，接近標(biāo)準(zhǔn)的半日潮流變化，響應(yīng)濃度隨潮時成正弦函數(shù)曲線變化。本文借鑒物理電工學(xué)中正弦交流電“有效值”的概念，提出響應(yīng)濃度的“有效濃度”值概念，即取一個潮周期內(nèi)響應(yīng)濃度最大值的0.707倍作為響應(yīng)濃度進(jìn)行計算。因此，響應(yīng)濃度用下式計算：

式中：Ci為有效濃度；Cimax為一個潮周內(nèi)響應(yīng)濃度最大值；K為有效濃度系數(shù)，取值0.707。式（13）適用于標(biāo)準(zhǔn)半日潮和全日潮（或接近兩種潮汐）類型海區(qū)。

3 大遼河口海域污染物負(fù)荷等級預(yù)警技術(shù)研究

3.1 大遼河口海域污染物負(fù)荷等級預(yù)警方程和評價體系

引用《環(huán)渤海污染壓力和海上響應(yīng)的統(tǒng)籌調(diào)控研究》對入海負(fù)荷超標(biāo)率計算及等級評估方法[22]，即式（14）計算入海負(fù)荷超標(biāo)率，表3海域環(huán)境承載力等級評價體系可作為海洋環(huán)境承載力評估的依據(jù)。

表3 入海負(fù)荷等級評價體系

式中：Kj為某種污染物第j個功能區(qū)入海負(fù)荷超標(biāo)率；為某種污染物第j個功能區(qū)的水質(zhì)現(xiàn)狀環(huán)境容量；為某種污染物第j個功能區(qū)的水質(zhì)標(biāo)準(zhǔn)環(huán)境容量；排放量單位均為t/a。

3.2 大遼河口海域環(huán)境承載力預(yù)警系統(tǒng)

依據(jù)上文所述計算方法和評價體系，運(yùn)用Visual Basic 6.0計算機(jī)語言編寫了大遼河口海域污染物承載力預(yù)警系統(tǒng)。使用者根據(jù)需要在界面輸入污染物的各項(xiàng)濃度，可以計算污染物的環(huán)境容量，再由此對海域環(huán)境承載力作出判斷。該系統(tǒng)操作簡單，可以在Windows系統(tǒng)下直接運(yùn)行。系統(tǒng)主要由3部分組成：現(xiàn)狀環(huán)境容量計算，標(biāo)準(zhǔn)環(huán)境容量計算，海區(qū)承載力等級評價和預(yù)警。該系統(tǒng)進(jìn)行承載力預(yù)警的基本流程如下：

（1）打開運(yùn)行文件，進(jìn)入系統(tǒng)主界面，選擇海區(qū)，以“大遼河口”為例，點(diǎn)擊計算（見圖7）。

圖7 選擇海區(qū)界面

（2）通過反演迭代法，分別計算大遼河口海域3個污染源處各污染要素的現(xiàn)狀環(huán)境容量：在程序界面中分別輸入3種污染要素的現(xiàn)狀負(fù)荷量以及本底濃度，然后得到各污染源現(xiàn)狀環(huán)境容量計算結(jié)果（見圖8a）。

（3）通過反演迭代法，分別計算大遼河口海域3個污染源處各污染要素的標(biāo)準(zhǔn)環(huán)境容量：在程序界面中分別輸入3種污染要素的標(biāo)準(zhǔn)濃度以及本底濃度，得到各污染源的標(biāo)準(zhǔn)環(huán)境容量（見圖8b）。

圖8 大遼河口海域環(huán)境容量計算結(jié)果

3.3 大遼河口海域環(huán)境承載力預(yù)警結(jié)果

3.3.1 海域和站位布設(shè)

大遼河口海域面積736 km2。海域各污染源坐標(biāo)見表4，海區(qū)功能規(guī)劃網(wǎng)格圖見圖6。

表4 大遼河口海域污染源坐標(biāo)

3.3.2 大遼河口海域水質(zhì)調(diào)查結(jié)果

2019年5月9日對營口海域10個站位的3種化學(xué)要素（COD、無機(jī)氮和磷酸鹽）進(jìn)行了調(diào)查，按《海洋監(jiān)測規(guī)范》方法執(zhí)行，水質(zhì)調(diào)查結(jié)果見表5。

表5 營口海區(qū)水質(zhì)調(diào)查結(jié)果

3.4 大遼河口海域入海污染物承載力預(yù)警結(jié)果

依據(jù)反演迭代法計算得到的現(xiàn)狀環(huán)境容量和標(biāo)準(zhǔn)環(huán)境容量，以及3.1節(jié)的預(yù)警方程對大遼河口海域環(huán)境承載力進(jìn)行評價，結(jié)果見圖9。該海域COD和磷酸鹽是1級可載狀態(tài)，無機(jī)氮標(biāo)準(zhǔn)環(huán)境容量小于0，是6級嚴(yán)重超載狀態(tài)。

圖9 大遼河口污染物環(huán)境承載力預(yù)警等級結(jié)果

4 結(jié)論

基于二維流體動力學(xué)方程和水質(zhì)控制模型，運(yùn)用保守物質(zhì)擴(kuò)散方程和反演迭代法計算出大遼河口海域3種主要污染要素的環(huán)境容量，以海洋功能區(qū)劃水質(zhì)標(biāo)準(zhǔn)為基準(zhǔn)進(jìn)行承載力評估，具體結(jié)論如下：

（1）以大遼河口不同類型功能區(qū)環(huán)境管理要求為基礎(chǔ)，運(yùn)用反演迭代法得到該海域3種主要污染物COD、無機(jī)氮和磷酸鹽的標(biāo)準(zhǔn)環(huán)境容量和現(xiàn)狀環(huán)境容量，并進(jìn)行海域環(huán)境承載力狀態(tài)評估。結(jié)果如下：COD和磷酸鹽處于1級可載狀態(tài)，無機(jī)氮處于6級嚴(yán)重超載狀態(tài)。

（2）對于遼東灣海區(qū)，在第一邊界條件下，使用保守物質(zhì)水質(zhì)擴(kuò)散方程計算的各項(xiàng)因子的濃度與污染負(fù)荷量是線性的、齊次的和可疊加的，在此基礎(chǔ)上構(gòu)建的反演迭代法為海洋環(huán)境容量計算提供了新思路。

（3）在標(biāo)準(zhǔn)的半日潮或全日潮海區(qū)，由于響應(yīng)濃度隨潮時呈正弦函數(shù)曲線變化，故一個潮周期響應(yīng)濃度的取值，可借鑒物理電工學(xué)理論有關(guān)正弦交流電的概念，方法為：取0.707的有效系數(shù)乘以一個潮周期內(nèi)響應(yīng)濃度最大值表示“有效濃度”，這一取值方法在標(biāo)準(zhǔn)半日潮和全日潮（或接近兩種潮汐）類型海區(qū)具有參考價值。

（4）本文編制的海域環(huán)境承載力預(yù)警系統(tǒng)軟件，實(shí)現(xiàn)了海域環(huán)境承載力等級評價的計算機(jī)程序化，基于水質(zhì)調(diào)查結(jié)果，對遼東灣海域環(huán)境承載力進(jìn)行分級評價和預(yù)警，以此為海洋管理部門合理規(guī)劃海域內(nèi)污染物的排放提供科學(xué)的調(diào)控依據(jù)。