劉 超， 胡智華， 姚天啟， 李 港

(1.重慶交通大學(xué)河海學(xué)院， 重慶 400074； 2.南京水利科學(xué)研究院生態(tài)環(huán)境研究所， 江蘇 南京 210029；3.貴州民族大學(xué)生態(tài)環(huán)境工程學(xué)院， 貴州 貴陽 550025)

洪澤湖流域地處經(jīng)濟(jì)高速發(fā)展的蘇北平原，流域內(nèi)工農(nóng)業(yè)以及水產(chǎn)養(yǎng)殖業(yè)發(fā)展迅速，為區(qū)域經(jīng)濟(jì)發(fā)展做出了強(qiáng)有力的貢獻(xiàn)，由此衍生出較為嚴(yán)重的環(huán)境問題，如地表水環(huán)境受到強(qiáng)烈的人為活動(dòng)干擾，氮磷質(zhì)量濃度遠(yuǎn)高于目標(biāo)水質(zhì)標(biāo)準(zhǔn)。洪澤湖屬于典型過水型湖泊，入湖河流是湖區(qū)水資源的主要補(bǔ)給方式，因此改善湖區(qū)水質(zhì)最關(guān)鍵手段是控制外源輸入。從1990—2018年的時(shí)間過程中，洪澤湖水質(zhì)整體呈現(xiàn)出V類到劣V類再到V類的變化趨勢(shì)，主要水質(zhì)超標(biāo)指標(biāo)為TN和TP[1]。研究表明2000年以來洪澤湖實(shí)施的各種治理和管控措施對(duì)水質(zhì)改善具有一定的成效，但湖區(qū)整體水質(zhì)仍然未達(dá)到III類目標(biāo)水質(zhì)標(biāo)準(zhǔn)，且空間差異仍然十分明顯。因此，進(jìn)一步開展降低洪澤湖水質(zhì)空間差異和改善洪澤湖整體水質(zhì)的研究具有重要意義。

洪澤湖流域涉及淮安市和宿遷市的6個(gè)區(qū)縣，各個(gè)區(qū)縣經(jīng)濟(jì)發(fā)展水平以及產(chǎn)業(yè)側(cè)重點(diǎn)各不相同，導(dǎo)致污染物的來源存在較大差異。同時(shí)受不同入湖河流水量以及湖區(qū)自身地形復(fù)雜致使湖泊各個(gè)區(qū)域水動(dòng)力條件有著較大差異，各入湖河流對(duì)洪澤湖不同區(qū)域水質(zhì)的影響程度不一。目前關(guān)于洪澤湖水質(zhì)改善的研究相對(duì)比較宏觀，如徐嘉興等[2]研究了洪澤湖流域土地利用組成，發(fā)現(xiàn)加大集中生活廢水處理力度以及減少水產(chǎn)養(yǎng)殖面積，可降低入湖河流污染負(fù)荷。閻敏等[3]研究表明農(nóng)業(yè)面源是入湖河流污染負(fù)荷主要來源，控制耕地面積有利于湖區(qū)水質(zhì)的改善。已有研究都相對(duì)缺乏以洪澤湖水質(zhì)作為目標(biāo)進(jìn)行精細(xì)化的探索，只有識(shí)別入湖源與湖區(qū)水質(zhì)在時(shí)空上的定量響應(yīng)關(guān)系，才能更加準(zhǔn)確實(shí)施控源措施以及提升水質(zhì)的改善效率。因此，準(zhǔn)確識(shí)別各入湖河流與湖區(qū)水質(zhì)在時(shí)間和空間上的定量響應(yīng)關(guān)系亟待解決。

本研究通過環(huán)境流體動(dòng)力學(xué)模型(Environmental fluid dynamics code, EFDC)建立洪澤湖水動(dòng)力水質(zhì)模型，將概化的7條主要入湖河流作為邊界，選取洪澤湖主要的水質(zhì)超標(biāo)指標(biāo)TN和TP作為模擬輸出目標(biāo)，探究不同入湖源對(duì)湖區(qū)水質(zhì)的定量響應(yīng)關(guān)系，優(yōu)化洪澤湖外源輸入控制措施，為提高洪澤湖水質(zhì)改善效率提供科學(xué)的理論依據(jù)[4]。

1 材料與方法

1.1 研究區(qū)域

洪澤湖是我國第四大淡水湖，屬于典型過水型湖泊，是淮河流域最大的防洪供水調(diào)蓄湖泊，也是南水北調(diào)東線工程的重要水源地，在調(diào)節(jié)氣候和保護(hù)生物多樣性等方向發(fā)揮著極其重要的作用。洪澤湖氣候四季分明，多年平均降水量為925.5 mm，每年6—9月降水量占全年總降水量的65.5%，冬季降水量僅占全年7%～8%。洪澤湖流域保護(hù)面積達(dá)到3 704 km2，湖泊保護(hù)范圍線長1 006 km。主要入湖河流有淮河、成子河、濉河、老濉河、新汴河、懷洪新河、徐洪河，淮河是洪澤湖主要的供水河道，出湖河流主要有淮河入江河道和蘇北灌溉總渠等。

1.2 水質(zhì)數(shù)據(jù)

江蘇省水文水資源勘測(cè)局在洪澤湖水質(zhì)監(jiān)測(cè)站點(diǎn)共15個(gè)，#1、#2、#3、#5監(jiān)測(cè)站點(diǎn)位于成子湖區(qū)，#4、#6、#7、#8監(jiān)測(cè)站點(diǎn)位于湖心區(qū)，#12、#13、#14、#15監(jiān)測(cè)站點(diǎn)位于溧河洼區(qū)，#9、#10、#11監(jiān)測(cè)站點(diǎn)位于入江過水通道區(qū)，各個(gè)監(jiān)測(cè)站點(diǎn)的監(jiān)測(cè)頻次為12次/a(每月1次)。本研究數(shù)據(jù)時(shí)間序列為2016—2018年，監(jiān)測(cè)指標(biāo)為DO、水位、TN和TP，DO采用水質(zhì)監(jiān)測(cè)儀直接測(cè)得，水位為當(dāng)?shù)厮徽緦?shí)時(shí)監(jiān)測(cè)數(shù)據(jù)。TN和 TP為現(xiàn)場(chǎng)分層(表、中、底)取混合水樣10L冷藏運(yùn)回實(shí)驗(yàn)室檢測(cè)，TN采用堿性過硫酸鉀消解紫外分光光度法(CB 11894—1989)，TP采用鉬酸銨分光光度法(GB 11893—1989)。

1.3 情景模擬方案設(shè)計(jì)

目前關(guān)于洪澤湖水質(zhì)變化特征的研究相對(duì)較多，已有研究表明洪澤湖水質(zhì)在時(shí)間和空間上都有著顯著的差異。利用2010—2019年洪澤湖長序列水質(zhì)監(jiān)測(cè)數(shù)據(jù)分析得到，TN和TP年平均質(zhì)量濃度變化范圍為1.39～1.86mg/L和0.08～0.171mg/L，變化幅度分別為33.81%和26.25%。為準(zhǔn)確探究入湖源在時(shí)間和空間上對(duì)湖區(qū)水質(zhì)的定量響應(yīng)關(guān)系，避免水質(zhì)在時(shí)間上的波動(dòng)以及水質(zhì)模型自身存在的誤差導(dǎo)致模擬結(jié)果變化不顯著的影響，對(duì)入湖源TN和TP質(zhì)量濃度以50%削減率作為標(biāo)準(zhǔn)情景工況，采用一次一個(gè)變量的單一控制變量法對(duì)7個(gè)情景方案進(jìn)行情景模擬，分析計(jì)算不同入湖源對(duì)湖區(qū)不同監(jiān)測(cè)點(diǎn)位TN和TP質(zhì)量濃度的改善能力。

2 水動(dòng)力水質(zhì)模型

2.1 模型原理

EFDC模型水質(zhì)模塊[5]中涉及氮的變量共有5個(gè)，分別為難溶性顆粒態(tài)有機(jī)氮(RPON)、活性顆粒態(tài)有機(jī)氮(LPON)、溶解態(tài)有機(jī)氮(DON)、氨氮(NH4)和硝態(tài)氮(包括硝酸鹽氮：NO2；亞硝酸鹽氮：NO3)。涉及磷的變量有4個(gè)，分別為難溶性顆粒態(tài)有機(jī)磷(RPOP)、活性顆粒態(tài)有機(jī)磷(LPOP)、溶解態(tài)有機(jī)磷(DOP)和TP酸鹽(PO4)。模型中各水質(zhì)變量彼此之間具有如下關(guān)系結(jié)構(gòu)。

本研究運(yùn)用的水動(dòng)力水質(zhì)模型是在Sigma坐標(biāo)下基于垂向靜水壓力的假設(shè)，并采用布西內(nèi)斯克方程(Boussinesq equation)近似得到如下的動(dòng)量方程和連續(xù)性方程。

x方向的動(dòng)量方程為：

(1)

y方向的動(dòng)量方程：

(2)

z方向的動(dòng)量方程:

(3)

連續(xù)性方程：

(4)

(5)

密度方程：

ρ=ρ(p,S,T,C)

(6)

(7)

P=myHuQ=mxHv

(8)

(9)

式中：(x，y)為水平方向的曲線-正交坐標(biāo)；z為垂向σ坐標(biāo)；(u，v)為(x，y)方向的水平速度分量，m/s；H為總水深，m；mx，my分別為x，y坐標(biāo)變化系數(shù)，在笛卡爾坐標(biāo)下，變化系數(shù)等于1；Patm為大氣壓強(qiáng)，Pa；p為參考密度ρ0下的附加靜水壓；b為浮力；f為科里奧利力系數(shù)，涵蓋網(wǎng)格曲率加速度；AH為水平動(dòng)量擴(kuò)散系數(shù)，m2/s；Av為垂向紊動(dòng)粘性系數(shù)，m2/s；cp為植被阻力系數(shù)；Dp為與每單位水平面積的流量相交的投影植被區(qū)域；Su和Sv為(x，y)方向動(dòng)量方程的源/匯項(xiàng)，m2/s2；Sh為質(zhì)量守恒方程的源/匯項(xiàng)，m3/s；S為鹽度，ng/L；T為溫度，℃；C為總懸浮無機(jī)顆粒質(zhì)量濃度，g/m3；U和V為(x，y)方向的質(zhì)量濃度平均速度分量，m/s；P和Q為(x，y)方向的質(zhì)量通量分量，m2/s，w為σ坐標(biāo)下的垂向速率，m/s；w*為z坐標(biāo)下的垂向速率，m/s。

2.2 模型構(gòu)建

本研究所建立的EFDC水動(dòng)力水質(zhì)模型是耦合了水動(dòng)力和水質(zhì)兩個(gè)模塊。洪澤湖岸線提取，利用正交網(wǎng)格將整個(gè)洪澤湖劃分成9 749個(gè)矩形網(wǎng)格單元，每個(gè)網(wǎng)格單元邊長為400 m×400 m，垂向采用Sigma坐標(biāo)系，平均分為3層，各層占水深比例為0.333，為了適應(yīng)水位波動(dòng)帶來的影響，設(shè)置臨界干水深為0.07 m。邊界條件有氣象邊界條件、水動(dòng)力邊界條件以及水質(zhì)邊界條件。氣象邊界條件主要包括降水、氣壓、氣溫、太陽輻射、濕度、蒸發(fā)量以及云覆蓋系數(shù)，氣象資料均來自洪澤湖周圍氣象站，獲取方式為中國氣象數(shù)據(jù)網(wǎng)。洪澤湖周邊河道復(fù)雜多樣，省略一些較小且數(shù)據(jù)不完善的河流，最后將主要入湖河流概化為7條，流量邊界采用2018年1月1日到2018年12月31日的逐日流量數(shù)據(jù)，數(shù)據(jù)來源為水文統(tǒng)計(jì)年鑒。選用洪澤湖水質(zhì)主要超標(biāo)項(xiàng)目TN和TP作為模型輸出目標(biāo)，水質(zhì)數(shù)據(jù)包含入湖河流和湖區(qū)15個(gè)監(jiān)測(cè)站點(diǎn)2016—2018年月尺度數(shù)據(jù)。初始水位和水溫均用2018年1月1日實(shí)測(cè)值。模型初始計(jì)算日期為2018年1月1日，運(yùn)行周期為364 d，時(shí)間步長為20 s。

2.3 模型參數(shù)敏感性分析

根據(jù)輸出目標(biāo)所涉及的模塊，通過敏感性分析選擇敏感性較高的參數(shù)，對(duì)所選擇的參數(shù)進(jìn)行敏感性分析，確定其最佳的取值。通過已有文獻(xiàn)確定模型參數(shù)的最小值與最大值，通過LHS(Latin hypercube sampling)方法抽樣分布函數(shù)，基于參數(shù)分布范圍利用LHS生成20組隨機(jī)參數(shù)組合，針對(duì)LHS生成的每一組參數(shù)，運(yùn)用模型進(jìn)行模擬，輸出目標(biāo)結(jié)果，評(píng)估參數(shù)不確定性對(duì)輸出目標(biāo)結(jié)果的不確定影響，最后采用標(biāo)準(zhǔn)秩回歸法進(jìn)行參數(shù)敏感性分析，確定各個(gè)參數(shù)的最優(yōu)取值，如表1所示[6]。

表1 模型校正后主要參數(shù)取值

2.4 模型率定和驗(yàn)證

為確保水動(dòng)力水質(zhì)模型的可靠度，需進(jìn)一步對(duì)模型進(jìn)行率定和驗(yàn)證，模型驗(yàn)證的目標(biāo)變量為水位和DO，模型率定數(shù)據(jù)采用2016年1月1日至2016年12月31日實(shí)測(cè)數(shù)據(jù)，模型驗(yàn)證數(shù)據(jù)采用2017年1月1日至2017年12月31日實(shí)測(cè)數(shù)據(jù)。

根據(jù)湖區(qū)3個(gè)水位站(高良澗站、蔣壩站、老子山站)實(shí)測(cè)水位數(shù)據(jù)與模擬水位數(shù)據(jù)對(duì)比分析可以發(fā)現(xiàn)模擬值與實(shí)測(cè)值較為接近，變化趨勢(shì)比較同步，模擬與實(shí)測(cè)數(shù)據(jù)決定系數(shù)(R2)均高于0.918，這說明模型在水位模擬上能夠較精確的反映真實(shí)水位情況，這為TN和TP等水質(zhì)指標(biāo)的模擬提供和可靠地支撐。利用湖區(qū)3個(gè)水質(zhì)監(jiān)測(cè)站點(diǎn)(#5、#14、#8)的DO實(shí)測(cè)數(shù)據(jù)與模擬值對(duì)比分析發(fā)現(xiàn)，模擬值與實(shí)測(cè)值的變化趨勢(shì)以及差值大小基本在可接受范圍之內(nèi)，說明建立的湖區(qū)水質(zhì)模型能夠較好的反映出湖區(qū)水質(zhì)的變化和分布情況，其模擬值與實(shí)測(cè)值決定系數(shù)(R2)最小為0.8246，最大值為0.8722。水動(dòng)力R2略高于水質(zhì)的R2，這表明水質(zhì)在模擬過程中的輸出結(jié)果相對(duì)于實(shí)測(cè)值的誤差比水動(dòng)力模擬輸出結(jié)果相對(duì)于實(shí)測(cè)值的誤差要稍大一些，出現(xiàn)這種情況主要是因?yàn)樗w中水環(huán)境因子變化復(fù)雜、彼此之間相互影響，加上水質(zhì)的模擬結(jié)果是建立在復(fù)雜的水動(dòng)力基礎(chǔ)之上，受到水動(dòng)力的影響。

3 結(jié)果與討論

3.1 入湖源與湖泊整體水質(zhì)的響應(yīng)關(guān)系

洪澤湖各入湖源總?cè)牒俊⒖偽廴矩?fù)荷以及在時(shí)間進(jìn)入湖區(qū)的水量分配皆有著較大差異，6—8月為豐水期，進(jìn)入湖區(qū)水量最大，10月為平水期，12月為枯水期，進(jìn)入湖區(qū)水量相對(duì)較小[7]。入湖水量最大為淮河入湖河流(301.4 × 108m3)，TN年平均質(zhì)量質(zhì)量濃度最高為新汴河(2.82 mg/L)入湖河流，TP年平均質(zhì)量質(zhì)量濃度最高為老濉河(0.183 mg/L)入湖河流。

圖1 水位和DO率定結(jié)果

為探究7個(gè)入湖源對(duì)湖區(qū)水體的影響關(guān)系，以代表性、全局性為原則在湖區(qū)布設(shè)了15個(gè)監(jiān)測(cè)站點(diǎn)。本研究以單一變量的研究方法進(jìn)行了8種情景工況(1組未削減和7組削減情景工況)的模擬，得到不同情景工況下湖區(qū)水體TN和TP的質(zhì)量質(zhì)量濃度的變化過程，并以15個(gè)監(jiān)測(cè)站點(diǎn)平均質(zhì)量質(zhì)量濃度來表示湖區(qū)水體營養(yǎng)鹽的整體質(zhì)量質(zhì)量濃度水平，從而分析入湖源對(duì)湖區(qū)整體水質(zhì)的影響能力。

削減模擬結(jié)果與未削減的基礎(chǔ)模擬結(jié)果如圖2所示。從時(shí)間上看，模擬前20 d各削減情景工況與未削減基礎(chǔ)情景工況的模擬結(jié)果未出現(xiàn)顯著差異，這主要是因?yàn)槟Ｐ瓦\(yùn)行初期模型尚未穩(wěn)定，故而受到模型初始條件的影響。初始條件是根據(jù)模擬年份第一天的實(shí)測(cè)數(shù)據(jù)進(jìn)行設(shè)定，只有在模型運(yùn)行穩(wěn)定，初始條件的影響消失過后計(jì)算出的結(jié)果才可靠。不同入湖源進(jìn)行相同比例削減模擬結(jié)果顯示，淮河對(duì)湖區(qū)TN和TP年平均質(zhì)量濃度改善效果皆為最好，主要是因?yàn)榛春尤牒髁看螅魉倏?、水?dòng)力條件好，因此污染物遷移能力強(qiáng)。

3.2 源邊界影響特征的時(shí)空間差異性

對(duì)于相同的入湖源，在不同的時(shí)間節(jié)點(diǎn)上，入湖源對(duì)水質(zhì)的改善效果也有差異，反之，不同的入湖源在相同的時(shí)間節(jié)點(diǎn)上亦如此。為探究不同入湖源在各自對(duì)應(yīng)的不同時(shí)間上對(duì)TN和TP質(zhì)量濃度的改善差異，分別選取了位于成子湖、溧河洼、湖心區(qū)、入江過水通道的4個(gè)監(jiān)測(cè)站點(diǎn)對(duì)應(yīng)第210天(豐水期)、第300天(平水期)、第350天(枯水期)進(jìn)行分析。結(jié)果如圖3～ 5所示。

在相同時(shí)間節(jié)點(diǎn)上各入湖源對(duì)不同監(jiān)測(cè)站點(diǎn)TN何TP的改善效果不同。如圖3中R3入湖源對(duì)#2號(hào)監(jiān)測(cè)站點(diǎn)TN改善率占比為40%，對(duì)#14號(hào)監(jiān)測(cè)站點(diǎn)TN改善率占比為4%，對(duì)#6號(hào)TN改善率占比為14%，對(duì)#9號(hào)TN改善率占比為21%。R3入湖源對(duì)#2、#14、#6、#9號(hào)TP改善率占比分別為65%、3%、5%、11%。同時(shí)可以發(fā)現(xiàn)入湖源對(duì)監(jiān)測(cè)站點(diǎn)水質(zhì)改善率主要表現(xiàn)為近高遠(yuǎn)低的響應(yīng)特征。但也有不同的響應(yīng)規(guī)律，如圖3中R1入湖源對(duì)#6和#9號(hào)監(jiān)測(cè)站點(diǎn)的改善率占比分別為65%和3%，表現(xiàn)出近低遠(yuǎn)高的響應(yīng)特征，這主要是因?yàn)樨S水期淮河入湖流量劇增、污染負(fù)荷增加，流速增大，在劇烈的水動(dòng)力條件下遷移方向較為徑直，導(dǎo)致出現(xiàn)對(duì)較遠(yuǎn)徑直站點(diǎn)影響更大[8]。這意味著監(jiān)測(cè)站點(diǎn)水質(zhì)改善效果主要受到入湖源水量、流速、污染負(fù)荷量以及監(jiān)測(cè)站點(diǎn)與入湖源之間距離的綜合影響。

在不同時(shí)間節(jié)點(diǎn)上，入湖源對(duì)同一個(gè)監(jiān)測(cè)站點(diǎn)改善效果不同。如圖3～5中R1入湖源對(duì)#14號(hào)監(jiān)測(cè)站點(diǎn)TN的結(jié)果顯示，在第210天時(shí)改善率占比為68%，第300天時(shí)改善率占比下降到55%，第350天時(shí)改善率占比進(jìn)一步下降到42%。圖3～5中R1入湖源對(duì)#14號(hào)監(jiān)測(cè)站點(diǎn)TP的結(jié)果顯示，在第210天時(shí)改善率占比為59%，第300天改善率占比為35%，第350天改善率占比為58%。這主要是因?yàn)樵诓煌闂l件下，入湖源的水量和攜帶的負(fù)荷總量的動(dòng)態(tài)變化過程明顯[9-10]。這表明對(duì)同一個(gè)監(jiān)測(cè)站點(diǎn)而言，充分考慮入湖源的削減時(shí)間對(duì)于提高水質(zhì)改善率有著正向效應(yīng)。對(duì)同一監(jiān)測(cè)站點(diǎn)在不同時(shí)間入湖源對(duì)其改善率不同。如圖3～5中就#9號(hào)監(jiān)測(cè)站點(diǎn)TN而言，在第210天時(shí)改善率最大為R5入湖源(26%)，在第300天時(shí)改善率最大為15%，在第350天時(shí)改善率最大為R1(44%)。這意味著，提高同一監(jiān)測(cè)站點(diǎn)的水質(zhì)改善率，在關(guān)注時(shí)間的同時(shí)也要對(duì)入湖源進(jìn)行有針對(duì)性的篩選。

圖3 TN和TP模擬210 d 4個(gè)監(jiān)測(cè)站點(diǎn)解釋結(jié)果

圖4 TN和TP模擬3 000 d 4個(gè)監(jiān)測(cè)站點(diǎn)解釋結(jié)果

圖5 TN和TP模擬350 d 4個(gè)監(jiān)測(cè)站點(diǎn)解釋結(jié)果

4 結(jié) 論

淮河入湖源對(duì)洪澤湖水質(zhì)指標(biāo)TN和TP平均質(zhì)量濃度的改善效果最好，改善率分別為16.2%，22.95%；徐洪河入湖源對(duì)最大值質(zhì)量濃度改善效果最好，改善率分別為17.34%，24.08%。當(dāng)以平均水質(zhì)作為目標(biāo)時(shí)，對(duì)淮河進(jìn)行控源效果更明顯，當(dāng)以降低水質(zhì)峰值質(zhì)量濃度為目標(biāo)時(shí)，因?qū)π旌楹雍统勺雍尤牒磳?shí)施控源措施。時(shí)間差異性，不同時(shí)間節(jié)點(diǎn)入湖源對(duì)同一個(gè)監(jiān)測(cè)站點(diǎn)改善效果不同。對(duì)于成子湖在豐平水期對(duì)徐洪河進(jìn)行控源對(duì)水質(zhì)就能有較好的改善效果，在枯水期需要對(duì)淮河與徐洪河聯(lián)合控源才能有較好的改善效果。空間差異性表明，入湖源對(duì)監(jiān)測(cè)站點(diǎn)水質(zhì)改善率主要表現(xiàn)為近高遠(yuǎn)低的響應(yīng)特征。環(huán)湖邊界水質(zhì)受就近入湖源影響較大，向湖心區(qū)移動(dòng)受淮河入湖源影響逐漸增強(qiáng)。入湖河流對(duì)洪澤湖水環(huán)境健康有著極其重要的影響，在以改善湖區(qū)水質(zhì)為目標(biāo)的管理中，應(yīng)當(dāng)綜合考慮各個(gè)入湖河流總水量、流速、水量的時(shí)間分配、污染物質(zhì)量濃度等特征，根據(jù)入湖河流與湖區(qū)水質(zhì)在時(shí)空上的定量響應(yīng)關(guān)系做出精細(xì)化的控源措施。