薛聯(lián)青，沈海岑，張 敏，汪 露

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098; 2.皖江工學(xué)院水利工程學(xué)院，安徽 馬鞍山 243000；3.江蘇省洪澤湖管理委員會辦公室，江蘇 淮安 223100)

洪澤湖是淮河流域重要的調(diào)蓄湖泊，淮河下游地區(qū)重要的水源地，同時也是南水北調(diào)東線工程的主要調(diào)節(jié)湖泊[1-2]。隨著區(qū)域人類活動加劇，洪澤湖面臨著自由水面面積縮減、局部富營養(yǎng)化、植被覆蓋度降低等多重環(huán)境問題[3-8]。湖泊的換水周期決定著輸移容量和水環(huán)境容量，是判斷湖泊污染物輸移能力的重要指標(biāo)[9]。風(fēng)場和吞吐流對淺水湖泊的水體交換過程有重要影響[10]，洪澤湖的湖面寬闊，風(fēng)場和吞吐流共同塑造其流場形態(tài)[11]，從而對洪澤湖污染物的擴散、遷移產(chǎn)生影響。因此探究不同風(fēng)場和不同季節(jié)來水條件下洪澤湖換水周期的時空分布，有利于識別水動力因素對污染物運移的影響[12]。

國內(nèi)外眾多學(xué)者對湖庫、海灣等水體換水能力開展了研究。Cucco等[13]構(gòu)建了威尼斯?jié)暫S水動力模型，通過潟湖內(nèi)釋放的被動示蹤劑的殘差函數(shù)來確定水體停留時間；Bocaniov等[14]利用三維水動力-富營養(yǎng)化模型計算了3種水體傳輸時間，研究發(fā)現(xiàn)營養(yǎng)鹽損失率與所有運輸時間尺度指標(biāo)均呈統(tǒng)計學(xué)顯著相關(guān)；李云良等[15]基于對流擴散理論計算不同季節(jié)鄱陽湖的換水周期和傳輸時間，研究發(fā)現(xiàn)季節(jié)水情對換水周期有較大影響；胥瑞晨等[16]研究得出水齡和水體停留時間存在63%的比例關(guān)系，并計算出太湖在不同風(fēng)場條件下不同方位湖區(qū)的水體交換率和半交換周期；陶磊等[17]基于歐拉觀點建立了渤海灣水齡計算模型，分析了潮汐和季風(fēng)兩個因素對渤海灣水體交換的影響。目前，對洪澤湖水體交換能力的研究僅停留在定量計算階段[18]，缺少對水體交換驅(qū)動機理的研究。本文基于MIKE21建立洪澤湖二維水動力模型，模擬計算洪澤湖不同湖區(qū)的換水周期，考慮洪澤湖的流場呈現(xiàn)風(fēng)生流和吞吐流雙重特點，分析不同季節(jié)和不同風(fēng)場條件下?lián)Q水周期的空間變化特征，從流場的角度分析水體交換的驅(qū)動機理，可為洪澤湖的水環(huán)境治理提供科學(xué)依據(jù)。

1 研究區(qū)概況

洪澤湖的多年平均水位(蔣壩站)為12.55 m，平均水深1.9 m，最大水深4.5 m。多年平均入湖水量為342億m3，入湖河流主要分布在湖泊西岸和南岸，包括淮河干流、懷洪新河、徐洪河、濉河、新汴河、老濉河和老汴河(圖1)，其中，淮河干流每年入湖水量超過70%。多年平均出湖水量為313億m3，出湖河道為二河、入江水道和蘇北灌溉總渠，其中，入江水道為主要出湖河道，控制工程為三河閘，出湖水量占總量的60%～70%。洪澤湖的不同湖區(qū)呈現(xiàn)不同的流場形態(tài)，流場形態(tài)對大型過水性湖泊的污染物運移方式有重要影響[19]，不同湖區(qū)的水質(zhì)狀況差別較大，結(jié)合湖泊特點將洪澤湖劃分為成子湖區(qū)、二河區(qū)、臨淮區(qū)、老子山區(qū)、三河區(qū)、中心湖區(qū)和溧河洼區(qū)7個湖區(qū)，如圖1所示。

圖1 洪澤湖分區(qū)及水系分布

2 研究方法

2.1 水動力模型原理

MIKE21水動力模型基于Navier-Stokes方程，并遵從Boussinesq假定。二維非恒定流淺水方程為

(1)

(2)

(3)

其中

2.2 水動力模型建立及驗證

采用洪澤湖2019年90 m×90 m的DEM數(shù)據(jù)，結(jié)合實測水深資料，在MIKE模型中進行修正，得到湖底高程數(shù)據(jù)，利用Mesh Generator工具將湖區(qū)劃分為19 475個非結(jié)構(gòu)三角形網(wǎng)格，共計10 888個節(jié)點。將洪澤湖的主要入湖口和出湖口作為開邊界，綜合考慮洪澤湖水系特征，入湖河道選擇淮河干流(包括池河)、懷洪新河、新汴河、濉河和汴河，出湖河道選擇三河、二河、蘇北灌溉總渠和徐洪河。上邊界(入湖邊界)輸入實測流量數(shù)據(jù)，下邊界(出湖邊界)輸入實測水位數(shù)據(jù)。采用2016年的實測水位數(shù)據(jù)對模型進行驗證，模擬時段為2016年1月1日8：00至2016年12月31日8：00。初始水位設(shè)為2015年12月31日蔣壩站水位13.12 m，初始流速設(shè)為0 m3/s，最小時間步長設(shè)為1 s。根據(jù)洪澤湖實際情況，設(shè)定干水深為0.006 m，濕水深為0.1 m，淹沒水深為0.05 m。主要參數(shù)設(shè)置為：曼寧系數(shù)45 m1/3/s、CFL數(shù)0.8、渦黏系數(shù)0.28。驗證結(jié)果見圖2，計算求得尚咀、老子山、高良澗閘和蔣壩4個驗證點的Nash系數(shù)分別為0.988 6、0.990 2、0.995 3和0.988 4，滿足模擬的精度要求。

圖2 模型驗證結(jié)果Fig.2 Results of model validation

2.3 換水周期計算原理

很多學(xué)者從不同角度定義湖庫的水體交換能力，比如水齡、示蹤劑傳輸時間、換水周期等[20-22]。湖泊換水周期是評價水體中由對流擴散導(dǎo)致的物質(zhì)交換速率的常用概念，本文根據(jù)湖泊特點，選擇湖泊換水周期來定量研究洪澤湖的水體交換能力。在湖泊中投入不可降解的示蹤劑，假設(shè)示蹤劑濃度變化呈現(xiàn)指數(shù)規(guī)律衰減[23]，則湖泊平均換水周期計算公式為

(4)

其中

式中：τ為湖體平均換水周期；rt為t時刻示蹤劑質(zhì)量濃度與初始時刻質(zhì)量濃度之比；ρt為t時刻示蹤劑質(zhì)量濃度；ρ′0為初始時刻示蹤劑質(zhì)量濃度。

由于現(xiàn)實中物質(zhì)濃度不可能降為0，當(dāng)t=τ時，示蹤劑濃度降為初始濃度的e-1或37%[24]，定義示蹤劑濃度衰減至初始濃度的37%所需的時間為換水周期。計算換水周期時，將初始示蹤劑空間質(zhì)量濃度場設(shè)為1 mg/L，入湖口質(zhì)量濃度設(shè)為0 mg/L，出湖口設(shè)為自由出流，通過編寫Python3.9代碼讀取每個非結(jié)構(gòu)網(wǎng)格點質(zhì)量濃度降解為0.37 mg/L的時刻，用空間連接功能導(dǎo)入到ArcGIS10.5中，生成湖泊換水周期的空間分布，并通過分區(qū)進行區(qū)域空間統(tǒng)計，計算出各個湖區(qū)的平均換水周期。

2.4 模擬情景設(shè)置

水文和風(fēng)場條件是研究淺水湖泊水體交換過程的基礎(chǔ)[25]。為研究不同季節(jié)、風(fēng)向?qū)闈珊w交換能力以及污染物運移的影響，設(shè)置幾種數(shù)值模擬情景(表1)，模擬時間均為365 d，即1個水文年。情景1模擬2013—2019年平均來水、氣象條件下的換水周期，將示蹤劑(初始濃度場)的投放時間設(shè)定為1月1日、4月1日、7月1日和10月1日，表征洪澤湖春、夏、秋和冬季的換水能力[17]。情景2設(shè)定實際風(fēng)場、無風(fēng)和典型主導(dǎo)風(fēng)場為1.88 m/s的東南風(fēng)與東北風(fēng)，模擬不同風(fēng)場條件對湖泊各區(qū)域換水周期的影響。

表1 模擬情景設(shè)置

3 結(jié)果與分析

3.1 換水周期時空分布規(guī)律

圖3為不同季節(jié)換水周期空間分布，圖4為各湖區(qū)不同季節(jié)的換水周期。由圖3、圖4可見，洪澤湖換水周期空間異質(zhì)性明顯，整體呈現(xiàn)出南部湖區(qū)水體更新時間快、北部湖區(qū)水體更新時間慢的特點，換水周期從南到北的梯度變化，主要與洪澤湖的水文情勢密切相關(guān)，淮河干流入流超過入湖流量的70%，且從南部老子山處入湖，故淮河上游的入流增強了南部湖區(qū)的水體交換能力。

圖3 不同季節(jié)換水周期的空間分布

圖4 各湖區(qū)不同季節(jié)的換水周期

冬、春兩季全湖平均換水周期較長，分別為75 d和60 d，秋、夏兩季的全湖平均換水周期較短，分別為49 d和31 d。冬、春兩季的換水周期空間分布相似，淮河入湖口(老子山區(qū))湖泊換水能力較強，換水周期約為10 d，呈現(xiàn)從淮河入湖區(qū)到二河、三河出湖區(qū)遞增的特點，三河區(qū)和中心湖區(qū)的換水周期均值為45 d左右，二河區(qū)的均值為83 d左右，出湖口局部水流交換較快，在10 d以內(nèi)。成子湖區(qū)的換水能力最弱，湖水流動緩慢，均值為157 d左右，成子湖區(qū)最北部的換水周期最長，達250 d。西部湖區(qū)溧河洼區(qū)和臨淮區(qū)在冬季的換水周期均值分別為80 d和46 d，春季較冬季換水能力增強，換水周期分別縮短了31.25%和37.5%。

夏季屬于豐水期，淮河干流入湖流量占全年的40%左右，上游大量來水導(dǎo)致全湖水體更新能力顯著增強，換水周期均值為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以內(nèi)，三河區(qū)、臨淮區(qū)和中心湖區(qū)總體為10～20 d，二河區(qū)為46 d，較夏季增加了16 d。成子湖區(qū)換水周期分布的空間異質(zhì)性明顯，南部徐洪河入湖口換水周期在10 d以內(nèi)，該區(qū)域面積較夏季增大，最北部的換水周期較其他季節(jié)明顯增加，最長達298 d，超過200 d的湖泊面積較其他季節(jié)顯著增加，這可能是由于南水北調(diào)東線工程運行之后，徐洪河作為南水北調(diào)東線的復(fù)線河道，承擔(dān)抽調(diào)洪澤湖水向北送的重要任務(wù)，上半年在抽調(diào)水的過程中，徐洪河口都是水流出湖狀態(tài)，根據(jù)多年平均流量顯示，7月開始徐洪河入湖口水流由出湖轉(zhuǎn)變?yōu)槿牒?，?0月后入湖流量增加，達100 m3/s，由于徐洪河湖區(qū)整體流速較緩，且風(fēng)生流在成子湖區(qū)有顯著作用，常出現(xiàn)風(fēng)生環(huán)流，因此夏季、秋季在徐洪河入湖口形成換水周期小于10 d的條帶，成子湖北部區(qū)域受環(huán)流影響較大，換水周期仍較長，超過200 d。

3.2 換水周期對風(fēng)場的響應(yīng)

風(fēng)場對換水周期的空間分布會產(chǎn)生一定的影響[26-27]，圖5為不同風(fēng)場作用下的換水周期，圖6為各湖區(qū)不同風(fēng)場作用下的換水周期。由圖5、圖6可見，換水周期總體呈現(xiàn)從東北到西南遞減的特征，對比實際風(fēng)場和無風(fēng)時洪澤湖換水周期的空間變化，可驗證風(fēng)場是影響洪澤湖污染物運移的重要驅(qū)動力。對比1.88 m2/s的東南風(fēng)和東北風(fēng)作用下的換水周期分布，可以看出不同的風(fēng)向使得換水周期發(fā)生相應(yīng)變化。

圖5 不同風(fēng)場作用下的換水周期

圖6 各湖區(qū)不同風(fēng)場作用下的換水周期Fig.6 Water exchange cycle under different wind fields in each lake area

總體來看，實際風(fēng)場作用下，湖泊換水周期的均值為72 d，比無風(fēng)時減少了16 d，說明實際風(fēng)場對洪澤湖的換水能力起到正向作用，風(fēng)生流可以促進污染物的運移擴散。風(fēng)場對洪澤湖區(qū)換水周期空間分布影響最大的是溧河洼區(qū)，在實際風(fēng)向作用下，其換水周期均值為92 d，比無風(fēng)時減少了89 d，說明實際風(fēng)場加速了溧河洼區(qū)的水體流動，增強了該區(qū)域污染物的稀釋能力。臨淮區(qū)的平均換水周期為55 d，減少了21 d，臨淮區(qū)和溧河洼區(qū)均位于湖西，可以看出在無風(fēng)條件即僅受吞吐流作用時，湖泊水體整體向東南、東北方向流動。中心湖區(qū)、二河區(qū)和成子湖區(qū)的平均換水周期在實際風(fēng)場的作用下分別為36 d、67 d和204 d，僅分別減少了5 d、5 d和7 d，而三河區(qū)和老子山區(qū)為40 d和10 d，增加了14 d和2 d，說明東北部湖區(qū)的換水能力受風(fēng)場的改善效果較小，主要進湖、出湖區(qū)域的換水能力略微減弱，其中，老子山區(qū)受風(fēng)場的影響最小，說明該區(qū)域的水體更新主要受吞吐流的影響。

在東南風(fēng)的作用下，湖泊換水周期的均值為73 d，比無風(fēng)時減少了15 d，說明東南風(fēng)對洪澤湖的換水能力起到正向作用，東南風(fēng)可以促進污染物的運移擴散，從而加強了污染物的稀釋能力。受東南風(fēng)影響最大的區(qū)域是溧河洼區(qū)，其換水周期均值為58 d，比無風(fēng)時減少了123 d，說明東南風(fēng)加速了溧河洼區(qū)的水體流動，增強了該區(qū)域污染物的稀釋能力，同屬西部湖區(qū)的臨淮區(qū)為48 d，減少了28 d，換水能力也明顯增強。東南風(fēng)也顯著改善了成子湖區(qū)的平均換水能力，其換水周期為158 d，較無風(fēng)時減少了53 d，最高值為193 d，較無風(fēng)狀態(tài)減少了47.1%。成子湖區(qū)流速緩慢，污染物容易在此聚集，受吞吐流影響較弱，所以湖流主要受風(fēng)場影響，東南風(fēng)加速了成子湖區(qū)的水體更新。東南風(fēng)對其他湖區(qū)的換水能力產(chǎn)生了減弱的作用，其中，受影響最大的是二河區(qū)，其換水周期均值為121 d，較無風(fēng)時增加了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，導(dǎo)致均值增加了8 d。東南部湖區(qū)換水能力普遍減弱，可能是由于東南風(fēng)的影響，淮河干流由老子山進入湖區(qū)后向三河區(qū)流動的水流受到阻礙，水體向東北方向流動的趨勢增強，加強了東北部水體污染物濃度的稀釋能力。

在東北風(fēng)的作用下，湖泊換水周期的均值為106 d，比無風(fēng)時增加了18 d，說明東北風(fēng)對洪澤湖的換水能力起到抑制作用，東北風(fēng)減緩了污染物的運移擴散速率，從而減弱了污染物的稀釋能力，同時，各個湖區(qū)換水周期的整體梯度性分布特征明顯，呈現(xiàn)東北高值，西南低值的局面。受東北風(fēng)影響最大的是臨淮區(qū)，其換水周期均值為145 d，比無風(fēng)時增加了69 d，湖泊在無風(fēng)狀態(tài)下，受吞吐流的影響，淮河干流入湖流向臨淮區(qū)，東北風(fēng)起到抑制作用，減弱了臨淮區(qū)對污染物的稀釋能力。二河區(qū)的換水能力也顯著減弱，其平均換水周期由76 d增加到145 d；中心湖區(qū)從41 d增加到69 d；三河區(qū)從26 d增加到43 d。東北風(fēng)僅對成子湖區(qū)的水體更新起到正向作用，但對成子湖區(qū)平均換水能力的改善效果較小，成子湖區(qū)的換水周期均值為205 d，較無風(fēng)時僅減少了2.8%；對換水周期最高值所在點位的水體改善效果較好，無風(fēng)時，成子湖區(qū)最北邊的水體交換周期接近365 d，東北風(fēng)作用下最高值為246 d，減少了32.6%。東北風(fēng)也降低了溧河洼區(qū)的換水周期，溧河洼區(qū)的換水周期為152 d，較無風(fēng)時減少了16%。東北風(fēng)對老子山區(qū)的影響較小，其平均換水周期由8 d增加到11 d，進一步驗證了老子山湖區(qū)的換水能力主要受吞吐流影響。

3.3 典型湖區(qū)流場特征

風(fēng)場、吞吐流以及湖泊形態(tài)影響洪澤湖流場的結(jié)構(gòu)，流場結(jié)構(gòu)對大型淺水湖泊的污染物運移過程和水體更新能力都會產(chǎn)生重要影響[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ū)水質(zhì)變化受外界輸入的影響減弱，受水體自凈的影響增強，而在豐水期，流速提高了2倍，外界輸入對老子山區(qū)水質(zhì)變化的貢獻度提高至90%。在來水較豐時期，淮河干流入湖流量大，導(dǎo)致老子山區(qū)水體流速加快，從而加快污染物的稀釋，提高湖泊的換水能力?；春痈闪鞯乃|(zhì)狀況是影響老子山區(qū)水質(zhì)的重要因素，但由于水體流速快，水質(zhì)變化受底泥釋放的影響較小，營養(yǎng)鹽不易聚集，局部富營養(yǎng)化的風(fēng)險較低。

(a) 2018年1月17日 (b) 2018年7月30日

由前文計算結(jié)果可知，成子湖區(qū)的換水周期大部分區(qū)域均超過150 d，模擬得到豐、枯來水條件下成子湖區(qū)的流場如圖8所示。由圖8可見，來水較少時，成子湖區(qū)大部分區(qū)域的流速為0.005～0.03 m/s，流速緩慢，受吞吐流影響較弱，主要受風(fēng)場影響，局部呈現(xiàn)逆向環(huán)流；來水較豐時期，成子湖區(qū)流速略微加快，大部分區(qū)域流速為0.012～0.04 m/s，仍然存在局部逆向環(huán)流。說明無論是豐水期還是枯水期，成子湖區(qū)受吞吐流的影響較小，流速緩慢，造成成子湖區(qū)的換水周期超過100 d，最北部湖區(qū)超過200 d。由于水流交換過緩，可以推測出外界輸入污染負荷對成子湖區(qū)水質(zhì)變化的影響較?。煌瑫r，流速緩慢的水動力條件下，水體中的氮磷營養(yǎng)鹽主要呈現(xiàn)溶解態(tài)，濃度高、水流緩慢且氣候適宜的條件下，藻類易繁殖和停留，湖泊有較高的水華風(fēng)險，高濃度營養(yǎng)鹽也加劇底泥釋放的風(fēng)險。

(a) 2018年1月17日

4 結(jié) 論

a.洪澤湖換水周期空間異質(zhì)性明顯，整體呈現(xiàn)出南部湖區(qū)水體更新快、北部湖區(qū)水體更新慢的特點，換水周期從南到北的梯度變化，主要與洪澤湖的水文情勢密切相關(guān)。老子山區(qū)流場主要受吞吐流影響，換水能力較強，換水周期約為10 d左右；成子湖區(qū)流速緩慢，流場主要受風(fēng)場影響，局部呈現(xiàn)逆向環(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.實際風(fēng)場對洪澤湖的換水能力起到正向作用，無風(fēng)時湖泊水體整體向東南、東北方向流動，換水周期均值為72 d，比無風(fēng)時減少16 d。風(fēng)場影響最大的是溧河洼區(qū)，比無風(fēng)時減少89 d，老子山區(qū)的水體更新受風(fēng)場的影響最小。東南風(fēng)對洪澤湖的換水能力起正向作用，比無風(fēng)時減少15 d；主要影響溧河洼區(qū)，比無風(fēng)時減少123 d；東北風(fēng)對洪澤湖的換水能力起到抑制作用，比無風(fēng)時增加18 d；主要影響臨淮區(qū)，比無風(fēng)時增加69 d。