薛聯青,沈海岑,張 敏,汪 露
(1.河海大學水文水資源學院,江蘇 南京 210098; 2.皖江工學院水利工程學院,安徽 馬鞍山 243000;3.江蘇省洪澤湖管理委員會辦公室,江蘇 淮安 223100)
洪澤湖是淮河流域重要的調蓄湖泊,淮河下游地區(qū)重要的水源地,同時也是南水北調東線工程的主要調節(jié)湖泊[1-2]。隨著區(qū)域人類活動加劇,洪澤湖面臨著自由水面面積縮減、局部富營養(yǎng)化、植被覆蓋度降低等多重環(huán)境問題[3-8]。湖泊的換水周期決定著輸移容量和水環(huán)境容量,是判斷湖泊污染物輸移能力的重要指標[9]。風場和吞吐流對淺水湖泊的水體交換過程有重要影響[10],洪澤湖的湖面寬闊,風場和吞吐流共同塑造其流場形態(tài)[11],從而對洪澤湖污染物的擴散、遷移產生影響。因此探究不同風場和不同季節(jié)來水條件下洪澤湖換水周期的時空分布,有利于識別水動力因素對污染物運移的影響[12]。
國內外眾多學者對湖庫、海灣等水體換水能力開展了研究。Cucco等[13]構建了威尼斯?jié)暫S水動力模型,通過潟湖內釋放的被動示蹤劑的殘差函數來確定水體停留時間;Bocaniov等[14]利用三維水動力-富營養(yǎng)化模型計算了3種水體傳輸時間,研究發(fā)現營養(yǎng)鹽損失率與所有運輸時間尺度指標均呈統(tǒng)計學顯著相關;李云良等[15]基于對流擴散理論計算不同季節(jié)鄱陽湖的換水周期和傳輸時間,研究發(fā)現季節(jié)水情對換水周期有較大影響;胥瑞晨等[16]研究得出水齡和水體停留時間存在63%的比例關系,并計算出太湖在不同風場條件下不同方位湖區(qū)的水體交換率和半交換周期;陶磊等[17]基于歐拉觀點建立了渤海灣水齡計算模型,分析了潮汐和季風兩個因素對渤海灣水體交換的影響。目前,對洪澤湖水體交換能力的研究僅停留在定量計算階段[18],缺少對水體交換驅動機理的研究。本文基于MIKE21建立洪澤湖二維水動力模型,模擬計算洪澤湖不同湖區(qū)的換水周期,考慮洪澤湖的流場呈現風生流和吞吐流雙重特點,分析不同季節(jié)和不同風場條件下換水周期的空間變化特征,從流場的角度分析水體交換的驅動機理,可為洪澤湖的水環(huán)境治理提供科學依據。
洪澤湖的多年平均水位(蔣壩站)為12.55 m,平均水深1.9 m,最大水深4.5 m。多年平均入湖水量為342億m3,入湖河流主要分布在湖泊西岸和南岸,包括淮河干流、懷洪新河、徐洪河、濉河、新汴河、老濉河和老汴河(圖1),其中,淮河干流每年入湖水量超過70%。多年平均出湖水量為313億m3,出湖河道為二河、入江水道和蘇北灌溉總渠,其中,入江水道為主要出湖河道,控制工程為三河閘,出湖水量占總量的60%~70%。洪澤湖的不同湖區(qū)呈現不同的流場形態(tài),流場形態(tài)對大型過水性湖泊的污染物運移方式有重要影響[19],不同湖區(qū)的水質狀況差別較大,結合湖泊特點將洪澤湖劃分為成子湖區(qū)、二河區(qū)、臨淮區(qū)、老子山區(qū)、三河區(qū)、中心湖區(qū)和溧河洼區(qū)7個湖區(qū),如圖1所示。
圖1 洪澤湖分區(qū)及水系分布
MIKE21水動力模型基于Navier-Stokes方程,并遵從Boussinesq假定。二維非恒定流淺水方程為
(1)
(2)
(3)
其中
采用洪澤湖2019年90 m×90 m的DEM數據,結合實測水深資料,在MIKE模型中進行修正,得到湖底高程數據,利用Mesh Generator工具將湖區(qū)劃分為19 475個非結構三角形網格,共計10 888個節(jié)點。將洪澤湖的主要入湖口和出湖口作為開邊界,綜合考慮洪澤湖水系特征,入湖河道選擇淮河干流(包括池河)、懷洪新河、新汴河、濉河和汴河,出湖河道選擇三河、二河、蘇北灌溉總渠和徐洪河。上邊界(入湖邊界)輸入實測流量數據,下邊界(出湖邊界)輸入實測水位數據。采用2016年的實測水位數據對模型進行驗證,模擬時段為2016年1月1日8:00至2016年12月31日8:00。初始水位設為2015年12月31日蔣壩站水位13.12 m,初始流速設為0 m3/s,最小時間步長設為1 s。根據洪澤湖實際情況,設定干水深為0.006 m,濕水深為0.1 m,淹沒水深為0.05 m。主要參數設置為:曼寧系數45 m1/3/s、CFL數0.8、渦黏系數0.28。驗證結果見圖2,計算求得尚咀、老子山、高良澗閘和蔣壩4個驗證點的Nash系數分別為0.988 6、0.990 2、0.995 3和0.988 4,滿足模擬的精度要求。
圖2 模型驗證結果Fig.2 Results of model validation
很多學者從不同角度定義湖庫的水體交換能力,比如水齡、示蹤劑傳輸時間、換水周期等[20-22]。湖泊換水周期是評價水體中由對流擴散導致的物質交換速率的常用概念,本文根據湖泊特點,選擇湖泊換水周期來定量研究洪澤湖的水體交換能力。在湖泊中投入不可降解的示蹤劑,假設示蹤劑濃度變化呈現指數規(guī)律衰減[23],則湖泊平均換水周期計算公式為
(4)
其中
式中:τ為湖體平均換水周期;rt為t時刻示蹤劑質量濃度與初始時刻質量濃度之比;ρt為t時刻示蹤劑質量濃度;ρ′0為初始時刻示蹤劑質量濃度。
由于現實中物質濃度不可能降為0,當t=τ時,示蹤劑濃度降為初始濃度的e-1或37%[24],定義示蹤劑濃度衰減至初始濃度的37%所需的時間為換水周期。計算換水周期時,將初始示蹤劑空間質量濃度場設為1 mg/L,入湖口質量濃度設為0 mg/L,出湖口設為自由出流,通過編寫Python3.9代碼讀取每個非結構網格點質量濃度降解為0.37 mg/L的時刻,用空間連接功能導入到ArcGIS10.5中,生成湖泊換水周期的空間分布,并通過分區(qū)進行區(qū)域空間統(tǒng)計,計算出各個湖區(qū)的平均換水周期。
水文和風場條件是研究淺水湖泊水體交換過程的基礎[25]。為研究不同季節(jié)、風向對洪澤湖水體交換能力以及污染物運移的影響,設置幾種數值模擬情景(表1),模擬時間均為365 d,即1個水文年。情景1模擬2013—2019年平均來水、氣象條件下的換水周期,將示蹤劑(初始濃度場)的投放時間設定為1月1日、4月1日、7月1日和10月1日,表征洪澤湖春、夏、秋和冬季的換水能力[17]。情景2設定實際風場、無風和典型主導風場為1.88 m/s的東南風與東北風,模擬不同風場條件對湖泊各區(qū)域換水周期的影響。
表1 模擬情景設置
圖3為不同季節(jié)換水周期空間分布,圖4為各湖區(qū)不同季節(jié)的換水周期。由圖3、圖4可見,洪澤湖換水周期空間異質性明顯,整體呈現出南部湖區(qū)水體更新時間快、北部湖區(qū)水體更新時間慢的特點,換水周期從南到北的梯度變化,主要與洪澤湖的水文情勢密切相關,淮河干流入流超過入湖流量的70%,且從南部老子山處入湖,故淮河上游的入流增強了南部湖區(qū)的水體交換能力。
圖3 不同季節(jié)換水周期的空間分布
圖4 各湖區(qū)不同季節(jié)的換水周期
冬、春兩季全湖平均換水周期較長,分別為75 d和60 d,秋、夏兩季的全湖平均換水周期較短,分別為49 d和31 d。冬、春兩季的換水周期空間分布相似,淮河入湖口(老子山區(qū))湖泊換水能力較強,換水周期約為10 d,呈現從淮河入湖區(qū)到二河、三河出湖區(qū)遞增的特點,三河區(qū)和中心湖區(qū)的換水周期均值為45 d左右,二河區(qū)的均值為83 d左右,出湖口局部水流交換較快,在10 d以內。成子湖區(qū)的換水能力最弱,湖水流動緩慢,均值為157 d左右,成子湖區(qū)最北部的換水周期最長,達250 d。西部湖區(qū)溧河洼區(qū)和臨淮區(qū)在冬季的換水周期均值分別為80 d和46 d,春季較冬季換水能力增強,換水周期分別縮短了31.25%和37.5%。
夏季屬于豐水期,淮河干流入湖流量占全年的40%左右,上游大量來水導致全湖水體更新能力顯著增強,換水周期均值為31 d,比冬季和春季減少了近50%,南部湖區(qū)均值小于10 d。其中溧河洼區(qū)、中心湖區(qū)、臨淮區(qū)和三河區(qū)水體更新能力較春季改善顯著,換水周期均從30~60 d減少為10 d左右,東部沿岸為10~20 d;二河區(qū)換水能力較弱,為30 d左右,但較春季和冬季減少了63%;成子湖區(qū)的換水周期均值為103 d,湖區(qū)最北端的換水周期長達239 d,但總體來說,均值較其他季節(jié)減少了34.8%,其中,徐洪河入湖口的換水周期小于10 d,成子湖區(qū)南部沿岸小于90 d,中心區(qū)域為90~150 d。
秋季為退水期,全湖的換水周期均值為49 d,較夏季增加了18 d。老子山和溧河洼湖區(qū)的換水周期總體仍在10 d以內,三河區(qū)、臨淮區(qū)和中心湖區(qū)總體為10~20 d,二河區(qū)為46 d,較夏季增加了16 d。成子湖區(qū)換水周期分布的空間異質性明顯,南部徐洪河入湖口換水周期在10 d以內,該區(qū)域面積較夏季增大,最北部的換水周期較其他季節(jié)明顯增加,最長達298 d,超過200 d的湖泊面積較其他季節(jié)顯著增加,這可能是由于南水北調東線工程運行之后,徐洪河作為南水北調東線的復線河道,承擔抽調洪澤湖水向北送的重要任務,上半年在抽調水的過程中,徐洪河口都是水流出湖狀態(tài),根據多年平均流量顯示,7月開始徐洪河入湖口水流由出湖轉變?yōu)槿牒?,?0月后入湖流量增加,達100 m3/s,由于徐洪河湖區(qū)整體流速較緩,且風生流在成子湖區(qū)有顯著作用,常出現風生環(huán)流,因此夏季、秋季在徐洪河入湖口形成換水周期小于10 d的條帶,成子湖北部區(qū)域受環(huán)流影響較大,換水周期仍較長,超過200 d。
風場對換水周期的空間分布會產生一定的影響[26-27],圖5為不同風場作用下的換水周期,圖6為各湖區(qū)不同風場作用下的換水周期。由圖5、圖6可見,換水周期總體呈現從東北到西南遞減的特征,對比實際風場和無風時洪澤湖換水周期的空間變化,可驗證風場是影響洪澤湖污染物運移的重要驅動力。對比1.88 m2/s的東南風和東北風作用下的換水周期分布,可以看出不同的風向使得換水周期發(fā)生相應變化。
圖5 不同風場作用下的換水周期
圖6 各湖區(qū)不同風場作用下的換水周期Fig.6 Water exchange cycle under different wind fields in each lake area
總體來看,實際風場作用下,湖泊換水周期的均值為72 d,比無風時減少了16 d,說明實際風場對洪澤湖的換水能力起到正向作用,風生流可以促進污染物的運移擴散。風場對洪澤湖區(qū)換水周期空間分布影響最大的是溧河洼區(qū),在實際風向作用下,其換水周期均值為92 d,比無風時減少了89 d,說明實際風場加速了溧河洼區(qū)的水體流動,增強了該區(qū)域污染物的稀釋能力。臨淮區(qū)的平均換水周期為55 d,減少了21 d,臨淮區(qū)和溧河洼區(qū)均位于湖西,可以看出在無風條件即僅受吞吐流作用時,湖泊水體整體向東南、東北方向流動。中心湖區(qū)、二河區(qū)和成子湖區(qū)的平均換水周期在實際風場的作用下分別為36 d、67 d和204 d,僅分別減少了5 d、5 d和7 d,而三河區(qū)和老子山區(qū)為40 d和10 d,增加了14 d和2 d,說明東北部湖區(qū)的換水能力受風場的改善效果較小,主要進湖、出湖區(qū)域的換水能力略微減弱,其中,老子山區(qū)受風場的影響最小,說明該區(qū)域的水體更新主要受吞吐流的影響。
在東南風的作用下,湖泊換水周期的均值為73 d,比無風時減少了15 d,說明東南風對洪澤湖的換水能力起到正向作用,東南風可以促進污染物的運移擴散,從而加強了污染物的稀釋能力。受東南風影響最大的區(qū)域是溧河洼區(qū),其換水周期均值為58 d,比無風時減少了123 d,說明東南風加速了溧河洼區(qū)的水體流動,增強了該區(qū)域污染物的稀釋能力,同屬西部湖區(qū)的臨淮區(qū)為48 d,減少了28 d,換水能力也明顯增強。東南風也顯著改善了成子湖區(qū)的平均換水能力,其換水周期為158 d,較無風時減少了53 d,最高值為193 d,較無風狀態(tài)減少了47.1%。成子湖區(qū)流速緩慢,污染物容易在此聚集,受吞吐流影響較弱,所以湖流主要受風場影響,東南風加速了成子湖區(qū)的水體更新。東南風對其他湖區(qū)的換水能力產生了減弱的作用,其中,受影響最大的是二河區(qū),其換水周期均值為121 d,較無風時增加了49 d。其次是三河區(qū),換水周期均值從26 d增加到52 d,換水能力減弱了50%,中心湖區(qū)的換水周期也增加了17 d,由41 d增加到58 d。同時,老子山區(qū)的均值也由8 d增加到16 d,老子山區(qū)的淮河入湖區(qū)域的換水周期仍然為1~10 d,而其東部靠近三河區(qū)的區(qū)域換水周期最高值為58 d,導致均值增加了8 d。東南部湖區(qū)換水能力普遍減弱,可能是由于東南風的影響,淮河干流由老子山進入湖區(qū)后向三河區(qū)流動的水流受到阻礙,水體向東北方向流動的趨勢增強,加強了東北部水體污染物濃度的稀釋能力。
在東北風的作用下,湖泊換水周期的均值為106 d,比無風時增加了18 d,說明東北風對洪澤湖的換水能力起到抑制作用,東北風減緩了污染物的運移擴散速率,從而減弱了污染物的稀釋能力,同時,各個湖區(qū)換水周期的整體梯度性分布特征明顯,呈現東北高值,西南低值的局面。受東北風影響最大的是臨淮區(qū),其換水周期均值為145 d,比無風時增加了69 d,湖泊在無風狀態(tài)下,受吞吐流的影響,淮河干流入湖流向臨淮區(qū),東北風起到抑制作用,減弱了臨淮區(qū)對污染物的稀釋能力。二河區(qū)的換水能力也顯著減弱,其平均換水周期由76 d增加到145 d;中心湖區(qū)從41 d增加到69 d;三河區(qū)從26 d增加到43 d。東北風僅對成子湖區(qū)的水體更新起到正向作用,但對成子湖區(qū)平均換水能力的改善效果較小,成子湖區(qū)的換水周期均值為205 d,較無風時僅減少了2.8%;對換水周期最高值所在點位的水體改善效果較好,無風時,成子湖區(qū)最北邊的水體交換周期接近365 d,東北風作用下最高值為246 d,減少了32.6%。東北風也降低了溧河洼區(qū)的換水周期,溧河洼區(qū)的換水周期為152 d,較無風時減少了16%。東北風對老子山區(qū)的影響較小,其平均換水周期由8 d增加到11 d,進一步驗證了老子山湖區(qū)的換水能力主要受吞吐流影響。
風場、吞吐流以及湖泊形態(tài)影響洪澤湖流場的結構,流場結構對大型淺水湖泊的污染物運移過程和水體更新能力都會產生重要影響[28-29]。選取換水周期最短的老子山區(qū)以及換水周期最長的成子湖區(qū)為典型湖區(qū),從流場的角度分析洪澤湖換水周期時空變化的機理。
老子山區(qū)的流場主要受淮河干流的入流影響,由圖7可見,較枯月份淮河干流流量約為400 m3/s,老子山區(qū)的流速為0.024~0.048 m/s,較豐月份淮河干流流量約為1 500 m3/s,老子山區(qū)的流速為0.09~0.195 m/s。在枯水期,來水較少,流速較慢,老子山區(qū)水質變化受外界輸入的影響減弱,受水體自凈的影響增強,而在豐水期,流速提高了2倍,外界輸入對老子山區(qū)水質變化的貢獻度提高至90%。在來水較豐時期,淮河干流入湖流量大,導致老子山區(qū)水體流速加快,從而加快污染物的稀釋,提高湖泊的換水能力?;春痈闪鞯乃|狀況是影響老子山區(qū)水質的重要因素,但由于水體流速快,水質變化受底泥釋放的影響較小,營養(yǎng)鹽不易聚集,局部富營養(yǎng)化的風險較低。
(a) 2018年1月17日 (b) 2018年7月30日
由前文計算結果可知,成子湖區(qū)的換水周期大部分區(qū)域均超過150 d,模擬得到豐、枯來水條件下成子湖區(qū)的流場如圖8所示。由圖8可見,來水較少時,成子湖區(qū)大部分區(qū)域的流速為0.005~0.03 m/s,流速緩慢,受吞吐流影響較弱,主要受風場影響,局部呈現逆向環(huán)流;來水較豐時期,成子湖區(qū)流速略微加快,大部分區(qū)域流速為0.012~0.04 m/s,仍然存在局部逆向環(huán)流。說明無論是豐水期還是枯水期,成子湖區(qū)受吞吐流的影響較小,流速緩慢,造成成子湖區(qū)的換水周期超過100 d,最北部湖區(qū)超過200 d。由于水流交換過緩,可以推測出外界輸入污染負荷對成子湖區(qū)水質變化的影響較??;同時,流速緩慢的水動力條件下,水體中的氮磷營養(yǎng)鹽主要呈現溶解態(tài),濃度高、水流緩慢且氣候適宜的條件下,藻類易繁殖和停留,湖泊有較高的水華風險,高濃度營養(yǎng)鹽也加劇底泥釋放的風險。
(a) 2018年1月17日
a.洪澤湖換水周期空間異質性明顯,整體呈現出南部湖區(qū)水體更新快、北部湖區(qū)水體更新慢的特點,換水周期從南到北的梯度變化,主要與洪澤湖的水文情勢密切相關。老子山區(qū)流場主要受吞吐流影響,換水能力較強,換水周期約為10 d左右;成子湖區(qū)流速緩慢,流場主要受風場影響,局部呈現逆向環(huán)流,換水能力最弱,換水周期超過100 d。
b. 冬季和春季全湖平均換水周期較長,分別為75 d和60 d,秋季和夏季的全湖平均換水周期較短,分別為49 d和31 d。溧河洼區(qū)和臨淮區(qū)春季較冬季換水能力增強,換水周期分別縮短了31.25%和37.5%。夏季南部湖區(qū)的均值小于10 d,受來水量影響,溧河洼區(qū)、中心湖區(qū)、臨淮區(qū)和三河區(qū)水體更新能力改善顯著。
c.實際風場對洪澤湖的換水能力起到正向作用,無風時湖泊水體整體向東南、東北方向流動,換水周期均值為72 d,比無風時減少16 d。風場影響最大的是溧河洼區(qū),比無風時減少89 d,老子山區(qū)的水體更新受風場的影響最小。東南風對洪澤湖的換水能力起正向作用,比無風時減少15 d;主要影響溧河洼區(qū),比無風時減少123 d;東北風對洪澤湖的換水能力起到抑制作用,比無風時增加18 d;主要影響臨淮區(qū),比無風時增加69 d。