生態(tài)清淤對(duì)于橋水庫水質(zhì)影響的數(shù)值模擬

程 揚(yáng)， 賴錫軍

(1.中國科學(xué)院 南京地理與湖泊研究所 中國科學(xué)院流域地理學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210008； 2.中國科學(xué)院大學(xué)，北京 100049)

1 研究背景

湖泊水庫流域是人類社會(huì)經(jīng)濟(jì)發(fā)展的重要區(qū)域。流域內(nèi)頻繁的人類活動(dòng)對(duì)湖泊水庫的生態(tài)環(huán)境形成了很大的壓力，水體日趨嚴(yán)重的富營養(yǎng)化現(xiàn)象就是其表現(xiàn)出來的突出問題之一[1]。水體富營養(yǎng)化威脅流域內(nèi)的供水安全和生態(tài)系統(tǒng)健康，嚴(yán)重制約了流域經(jīng)濟(jì)社會(huì)的可持續(xù)發(fā)展[2]。內(nèi)源污染是研究湖庫富營養(yǎng)化問題不可忽略的因素。即使外源污染逐漸得到有效控制，湖泊底泥的氮磷釋放仍有可能使水體保持在富營養(yǎng)化狀態(tài)[3]。生態(tài)清淤將污染底泥從湖底清除，從根本上控制內(nèi)源性污染，被廣泛應(yīng)用于國內(nèi)外的湖泊水庫治理中[4]。然而，不同湖泊實(shí)施清淤后對(duì)營養(yǎng)鹽釋放的控制效果不同，對(duì)生態(tài)清淤是否能從根本上改善富營養(yǎng)化狀態(tài)，國內(nèi)外仍存在爭(zhēng)議[5-6]。準(zhǔn)確認(rèn)識(shí)動(dòng)態(tài)變化環(huán)境下沉積物-水界面營養(yǎng)鹽交換過程及其影響下水體水質(zhì)的變化是評(píng)估生態(tài)清淤工程能否實(shí)施的關(guān)鍵。目前，國內(nèi)對(duì)清淤效果的評(píng)價(jià)多為清淤后評(píng)估驗(yàn)證[7-8]，而在清淤工程開展前缺乏必要的研究和論證，不僅為湖泊治理工作帶來隱患，同時(shí)也造成了資源資金的浪費(fèi)。在生態(tài)清淤工程實(shí)施前利用模型模擬和預(yù)測(cè)清淤效果，可為生態(tài)清淤工程的決策提供重要支撐。

于橋水庫是天津市主要的飲用水源地。作為引灤入津工程的重點(diǎn)調(diào)蓄水庫，其在天津市供水系統(tǒng)中發(fā)揮著極其重要的作用。目前，于橋水庫整體水質(zhì)呈中營養(yǎng)-輕富營養(yǎng)狀態(tài)，且富營養(yǎng)化有加重的趨勢(shì)[9]。水庫水生態(tài)系統(tǒng)已遭到破壞，底泥對(duì)污染物的截留(自凈)能力已趨于飽和，底泥污染初顯。在不斷推進(jìn)水庫上游污染治理且取得一定效果的形勢(shì)下，對(duì)于橋水庫進(jìn)行清淤是進(jìn)一步改善水質(zhì)的重要備選手段。

本研究以于橋水庫生態(tài)清淤為例，基于二維淺水方程和對(duì)流-擴(kuò)散-反應(yīng)方程構(gòu)建了考慮底泥釋放的水動(dòng)力-營養(yǎng)物質(zhì)輸移模型，分別模擬了不同工況、不同清淤范圍條件下清淤工程開展前后地表水環(huán)境的時(shí)間和空間變化情況，以期對(duì)清淤工程的設(shè)計(jì)和實(shí)施提供數(shù)據(jù)支撐和合理建議。

2 考慮沉積物-水界面物質(zhì)交換的二維水動(dòng)力和營養(yǎng)物質(zhì)輸移模型

2.1 二維對(duì)流-擴(kuò)散-反應(yīng)方程

采用二維對(duì)流-擴(kuò)散-反應(yīng)方程模擬氮、磷營養(yǎng)物質(zhì)在湖泊水庫內(nèi)部的輸移降解過程[10-11]。其方程為：

(1)

式中：h為水深，m；C為水體中氮、磷的濃度，mg/L；u、v分別為x、y方向的流速分量，m/s；Dx、Dy分別為氮、磷在水體中x,y方向的擴(kuò)散系數(shù)，m2/s；K為氮、磷在水體中的綜合降解系數(shù)，s-1；S為氮、磷的源匯項(xiàng)。

為了考慮底泥釋放和吸收等內(nèi)部源匯項(xiàng)，源項(xiàng)S需增加考慮沉積物底泥與表層水之間氮磷的交換。令Fex為底沉積物-水界面氮、磷交換通量，其紊動(dòng)擴(kuò)散理論[12]為：

(2)

考慮沉積物-水界面的氮磷交換通量，則由公式(1)可得垂向平均平面二維氮、磷營養(yǎng)物質(zhì)輸移方程，其表達(dá)式為：

(3)

2.2 水動(dòng)力-水質(zhì)耦合的守恒型方程

將二維淺水方程與二維氮、磷營養(yǎng)物質(zhì)輸運(yùn)方程耦合，形成的方程組可以寫成以下守恒形式[13]：

(4)

其中：

(5)

(6)

(7)

(8)

式中：U為守恒物理向量；Fx，F(xiàn)y分別為x，y方向的通量向量；S為源匯項(xiàng)向量；S0x，S0y分別為x，y方向的床底坡降；Sfx，Sfy分別為x，y方向的摩阻坡降；g為重力加速度，m/s2；cω為風(fēng)的阻力系數(shù)；ρa(bǔ)為空氣的密度，kg/m3；ω為風(fēng)速，m/s；α為風(fēng)速與y軸的夾角。

2.3 數(shù)值方法

采用非結(jié)構(gòu)網(wǎng)格有限體積法求解該二維水動(dòng)力-氮磷營養(yǎng)鹽輸移的控制方程組[13-14]。首先對(duì)方程在任意控制體積做體積分，后根據(jù)高斯定理將體積分化成面積分，再進(jìn)行空間離散，時(shí)間導(dǎo)數(shù)項(xiàng)采用一階顯格式離散，離散后方程的法向通量采用HLLC算法[13]計(jì)算。

2.4 邊界條件

以上數(shù)值方法可用于計(jì)算域內(nèi)部單元界面的求解，而對(duì)于計(jì)算域的邊界單元，其左狀態(tài)為已知狀態(tài)，而右狀態(tài)為未知狀態(tài)，此時(shí)的數(shù)值通量求解就變成了邊界黎曼問題[15]。一般可以通過根據(jù)局部流態(tài)適當(dāng)選定的輸出特征的相容關(guān)系和指定邊界條件來確定未知狀態(tài)，常見的邊界條件可以為水位、單寬流量和水位-流量關(guān)系。

對(duì)于污染物的擴(kuò)散輸移，存在三類邊界：水流入口處可給定邊界濃度過程或污染負(fù)荷過程；水流出口處可給定自由輸出的Neumann邊界條件；面域上，可給定污染負(fù)荷，加入源項(xiàng)中計(jì)算。

3 于橋水庫模型構(gòu)建

3.1 于橋水庫庫區(qū)及計(jì)算網(wǎng)格

于橋水庫位于天津市薊州區(qū)北部的州河上，是國家重點(diǎn)大型水庫(圖1)。正常蓄水水位21.16 m，平均水深 4.3 m，東西長(zhǎng)約30 km，南北寬約8 km，總庫容15.59×108m3，控制流域面積2 060 km2。于橋水庫流域四面環(huán)山，流域內(nèi)小河流眾多，匯集成較大的3個(gè)水系——黎河、沙河與淋河。庫區(qū)水下地形東高西低，在低水位時(shí)有灘地出露。為了保障于橋水庫的飲用水供水安全，天津市在東側(cè)河口建設(shè)了于橋水庫入庫河口濕地工程。該工程匯集了黎河、沙河和引灤入津的來水，經(jīng)過凈化處理后進(jìn)入于橋水庫。因此，本次計(jì)算中僅考慮水庫的主庫區(qū)，不納入濕地區(qū)塊。

圖1 于橋水庫流域及水系

為了保證較高的計(jì)算效率，又能較好地反映于橋水庫地形特征，采用非結(jié)構(gòu)網(wǎng)格將計(jì)算區(qū)域共劃分為2 280個(gè)計(jì)算單元(圖2)。其中以四邊形網(wǎng)格為主，三角形網(wǎng)格為輔，各單元的邊長(zhǎng)從40m到250m不等，入庫河道主槽單元較為精細(xì)，而在地形變化不大的灘地和地形變化較大卻常年處于淹水狀態(tài)的區(qū)域則采用了較大尺度的單元。

圖2 于橋水庫計(jì)算網(wǎng)格

3.2 初始條件和邊界條件

將在初始時(shí)刻的入庫流量和出庫水位驅(qū)動(dòng)下計(jì)算24 h得到準(zhǔn)恒定條件下的水流和水質(zhì)條件作為計(jì)算的初始場(chǎng)。

于橋水庫的主要入庫水量和營養(yǎng)物質(zhì)來自于果河(沙河、黎河)和淋河。計(jì)算時(shí)，以果河和淋河入庫流量和營養(yǎng)物質(zhì)濃度作為邊界條件。在下游于橋水庫的出口，給定水庫出口的水位過程作為邊界條件，水質(zhì)的邊界條件設(shè)成自由輸出的無反射邊界條件。

3.3 參數(shù)率定

選擇2015年全年的水量和水質(zhì)數(shù)據(jù)對(duì)建立的于橋水庫平面二維水動(dòng)力-水質(zhì)模型進(jìn)行率定和檢驗(yàn)。本次于橋水庫清淤計(jì)算主要關(guān)注于總氮、總磷的指標(biāo)變化。模型率定的參數(shù)主要包括水動(dòng)力參數(shù)、底泥參數(shù)和水質(zhì)參數(shù)。

水動(dòng)力參數(shù)：于橋水庫水動(dòng)力較弱，且?guī)靺^(qū)大部分區(qū)域多為常年淹沒狀態(tài)。計(jì)算時(shí)將曼寧糙率系數(shù)統(tǒng)一取為0.021，風(fēng)拖曳力系數(shù)為0.0026。

底泥參數(shù)：沉積物上表層平均孔隙度為0.82，沉積物底層平均孔隙度為0.8，總氮和總磷的理想擴(kuò)散系數(shù)均為19.8×10-10m2/s；沉積物中的氮、磷釋放速率主要由沉積物間隙水和上覆水中氮、磷營養(yǎng)物質(zhì)的濃度梯度決定，根據(jù)中國科學(xué)院南京地理與湖泊研究所開展的于橋水庫底泥釋放調(diào)查和實(shí)驗(yàn)數(shù)據(jù)[16]，設(shè)定在現(xiàn)狀條件下底泥間隙水總氮和總磷濃度分別為7 mg/L和0.8 mg/L，清淤后間隙水總氮和總磷濃度分別下降至3.5 mg/L和0.01 mg/L。

水質(zhì)參數(shù)：20℃時(shí)總氮和總磷的降解系數(shù)分別為0.04 d-1和0.02 d-1，考慮降解系數(shù)隨溫度變化，根據(jù)菲爾普斯公式[17]，不同水溫條件下的降解系數(shù)值與20℃值的關(guān)系為：

KT=K20β(T-20)

(9)

式中:β為溫度校正系數(shù)。對(duì)于總氮和總磷，分別取1.07和1.04;KT和K20分別為水溫T℃和20℃時(shí)的降解系數(shù)，d-1。

3.4 模型驗(yàn)證

3.4.1 水庫水質(zhì) 基于模型參數(shù)和2015年的邊界條件，計(jì)算得到了于橋水庫全年的水質(zhì)動(dòng)態(tài)變化過程，采用庫中心和出庫兩監(jiān)測(cè)點(diǎn)的總氮和總磷實(shí)測(cè)數(shù)據(jù)對(duì)模型計(jì)算結(jié)果進(jìn)行驗(yàn)證(圖3)。結(jié)果表明，模型較好地反演了2015年總氮、總磷濃度的上下波動(dòng)變化過程?？偟趲熘行暮统鰩靸烧镜南鄬?duì)誤差分別為17%和24%；總磷的相對(duì)誤差分別為26%和32%。

圖3 2015年于橋水庫總氮和總磷濃度實(shí)測(cè)值與模擬值對(duì)比

3.4.2 底泥釋放通量 由于實(shí)驗(yàn)數(shù)據(jù)的限制，選擇氨氮進(jìn)行底泥釋放通量的驗(yàn)證。對(duì)于模型參數(shù)，水溫為20℃時(shí)氨氮的降解系數(shù)取0.1 d-1，理想擴(kuò)散系數(shù)和溫度校正系數(shù)與總氮取值相同，其他水動(dòng)力參數(shù)、底泥參數(shù)及水質(zhì)參數(shù)與3.3節(jié)一致，沉積物間隙水中氨氮的濃度為6.0 mg/L。利用模型求得于橋水庫2015年8月庫區(qū)底泥的氨氮釋放通量155.60 mg/(m2·d)，與文帥龍等[16]對(duì)于橋水庫沉積物-水界面交換通量的實(shí)驗(yàn)研究結(jié)果129.55 mg/(m2·d)基本一致，相對(duì)誤差20.11%，能夠反映底泥向水體釋放氨氮的水平。

3.4.3 實(shí)驗(yàn)設(shè)計(jì) 水庫的水動(dòng)力、水質(zhì)背景狀況和清淤方案的選擇都會(huì)影響分析結(jié)果。綜合考慮于橋水庫實(shí)施清淤工程的開展條件和環(huán)境背景，分別設(shè)計(jì)了兩種工況和兩種方案。

根據(jù)天津市水利科學(xué)研究院與中國科學(xué)院南京地理與湖泊研究所承擔(dān)的“于橋水庫環(huán)境現(xiàn)狀診斷及其對(duì)上游調(diào)水的響應(yīng)研究”成果，于橋水庫底泥總磷分布較為平均，但總氮含量在壩前、北岸峰山至宋家營大堤段、劉相營段兩側(cè)以及南岸東部等區(qū)域較高，因此在綜合考慮施工便利條件、資源資金的合理利用和清淤工程效果的前提下，設(shè)計(jì)了庫周清淤和全庫清淤兩個(gè)方案：(1)庫周清淤：清淤范圍以常水位以下、高程17.8 m以上的庫周區(qū)域(圖4)，具體包括北岸、東岸和南岸東部，水庫北岸富營養(yǎng)底泥清除厚度為30 cm，水庫東部及南岸為40 cm。富營養(yǎng)底泥清除總面積13.00 km2，共計(jì)清除底泥488×104m3。(2)庫區(qū)完全清淤：全庫清除表層底泥，清除深度為35 cm。

圖4 庫周清淤范圍

2014年于橋水庫入庫河口濕地工程開工建設(shè)，建成后上游來水排入水庫前需先經(jīng)過庫前區(qū)域的濕地生態(tài)系統(tǒng)的凈化，入庫水體所含營養(yǎng)鹽等污染物得到一定程度的削減。基于此背景，本研究針對(duì)河口濕地啟用前后兩種工況分別進(jìn)行模擬和分析：(1)2015年實(shí)際水流和水質(zhì)(下文簡(jiǎn)稱“現(xiàn)狀”)：采用2015年實(shí)際觀測(cè)數(shù)據(jù)驗(yàn)證計(jì)算后的水動(dòng)力-物質(zhì)輸移參數(shù)，給定清淤后清淤位置的沉積物間隙水氮、磷濃度；(2)未來河口濕地正常運(yùn)行(下文簡(jiǎn)稱“未來”)：采用濕地設(shè)計(jì)的引灤入津水量54 m3/s作為入庫水量，庫水位保持在20.6m；入庫水質(zhì)總氮和總磷濃度分別為3.14 mg/L和0.11 mg/L。

4 生態(tài)清淤的影響分析

4.1 總體影響特征

在水庫外源不變的情況下，生態(tài)清淤可降低水庫氮、磷營養(yǎng)鹽濃度(表1)，改善庫區(qū)水環(huán)境富營養(yǎng)化狀況。

表1 清淤前后水體TN、TP年平均濃度

對(duì)于TN，若入庫水質(zhì)維持現(xiàn)狀，庫周清淤和全庫清淤兩種方案可使庫中心年平均濃度從1.273降至1.270和0.908 mg/L，分別下降0.23%和28.70%；若未來河口濕地正常運(yùn)行，入庫水質(zhì)改善，兩種方案可分別使TN年平均濃度從1.144降至1.139和0.775 mg/L，分別下降0.48%和32.29%。無論是在現(xiàn)狀或未來的水質(zhì)背景條件下，庫周或全庫清淤方案下，庫區(qū)水體中TN濃度均有所下降。但是，庫周清淤方案效果有限，即使清除了庫岸區(qū)域的高濃度底泥，庫區(qū)水體TN濃度減少并不明顯，而全庫清淤治理效果顯著，現(xiàn)狀和未來背景下，TN濃度均大幅下降。

對(duì)于TP，現(xiàn)狀工況下，清淤前庫中心年平均濃度為0.066 mg/L，庫周和全庫清淤實(shí)施后，年平均濃度降至0.065和0.049 mg/L，分別下降0.52%和26.06%；未來工況下，TN年平均濃度為0.036 mg/L，庫周和全庫清淤實(shí)施后，年平均濃度降至0.035和0.019mg/L，分別下降1.64%和48.05%。總體變化規(guī)律與TN相似，除現(xiàn)狀背景下全庫清淤的效果比TN略低外，其他工況和清淤方案下的治理效果均優(yōu)于TN。

雖然清淤后于橋水庫水體中TN、TP均有所減少，但是值得注意的是，這種改善并不能簡(jiǎn)單地認(rèn)為是長(zhǎng)效的，沉積物底泥和水體中的營養(yǎng)鹽存在一個(gè)再平衡的過程，上覆水濃度變化后，底泥中的營養(yǎng)鹽含量也會(huì)隨之變化。底泥清除后，水庫來水營養(yǎng)鹽含量和本底若仍維持較高的水平，營養(yǎng)鹽仍然會(huì)在底泥累積使得底泥再次成為內(nèi)部釋放源[18]。如我國杭州西湖幾十年來經(jīng)歷了3次大規(guī)模的庫周和全湖清淤，清淤后富營養(yǎng)化狀態(tài)有所改善，但維持的時(shí)間不長(zhǎng)[19-20]。因此，水庫的清淤工程僅是富營養(yǎng)化治理的一個(gè)方面，以改善水庫富營養(yǎng)化狀況為目的水污染防治應(yīng)該綜合外部營養(yǎng)輸入的控制、內(nèi)部生物的營養(yǎng)利用等措施共同開展。

從TN、TP濃度的年內(nèi)變化情況(圖5)看，生態(tài)清淤工程實(shí)施前后變化趨勢(shì)基本不變，清淤并不能改變水庫營養(yǎng)鹽的變化過程，庫區(qū)水體營養(yǎng)鹽的變化主要取決于外部的輸入。

4.2 水質(zhì)分布空間變化

根據(jù)上文分析，以清淤對(duì)水環(huán)境改善效果較好的2015年5月1日為例，對(duì)比分析不同工況、不同范圍的清淤工程實(shí)施前后庫區(qū)總氮和總磷的空間分布情況。

現(xiàn)狀水質(zhì)背景條件下(圖6)，清淤前TN空間分布總體表現(xiàn)為庫周較高而中心范圍較低，東部入庫口向西部出庫口逐漸降低，這與徐媛等[8]對(duì)于橋水庫富營養(yǎng)化空間分布的調(diào)查研究結(jié)果一致，表明庫區(qū)水體TN濃度受水流影響較大，上游來水導(dǎo)致入庫口TN濃度居高不下[21]。庫周清淤方案實(shí)施后，可以明顯看到庫周淺水區(qū)的TN濃度下降明顯，但庫中心區(qū)域底泥向庫區(qū)水體中釋放的營養(yǎng)鹽含量依然很高，庫區(qū)水體TN濃度下降不明顯。全庫清淤后，整個(gè)庫區(qū)水體中TN含量都明顯降低，其中，東部水域TN下降較西部地區(qū)更加明顯，可能的原因是庫區(qū)西部水域本底值較低，TN濃度下降值要比東部低。

圖5 清淤前后水體TN、TP濃度年變化過程

圖6 現(xiàn)狀水質(zhì)背景下清淤前后水體TN、TP空間分布的變化

對(duì)于TP，清淤前其分布狀況與TN總體一致，但濃度差距較小，整體水質(zhì)相對(duì)穩(wěn)定。庫周清淤后對(duì)沿岸高濃度地帶治理效果較好，全庫清淤后整個(gè)庫區(qū)水體東西分布均勻，較清淤前明顯降低，清淤工程實(shí)施效果顯著。

河口濕地正常運(yùn)行后的未來水質(zhì)背景下(圖7)，上游來水污染物濃度大幅降低，外源污染的減少使水庫水質(zhì)較現(xiàn)狀背景得到有效改善?？臻g分布上，上游來水TN濃度的降低導(dǎo)致入庫口TN濃度大幅降低，上下游TN污染空間差異減小。TP空間分布與TN總體一致。

上述結(jié)果表明，庫周清淤仍只能使庫周淺水區(qū)和東部范圍的N、P濃度下降，全庫清淤效果明顯，東西部N、P濃度差進(jìn)一步降低，庫周淺水區(qū)和中心水體的N、P濃度基本一致。

圖7 未來水質(zhì)背景下清淤前后水體TN、TP空間分布的變化

5 結(jié)論和討論

(1)本文基于二維淺水方程和對(duì)流-擴(kuò)散-反應(yīng)方程，建立了考慮沉積物-水界面物質(zhì)交換的二維水動(dòng)力和營養(yǎng)物質(zhì)輸移模型。結(jié)果表明，模型能夠有效地模擬水庫底泥氮、磷營養(yǎng)鹽的動(dòng)態(tài)釋放及其與外界污染輸入共同作用下水庫的氮、磷營養(yǎng)鹽的時(shí)空分布特征，可用于生態(tài)清淤工程實(shí)施效果的預(yù)測(cè)。

(2)于橋水庫底泥是重要的污染來源，底泥的清除對(duì)水庫氮、磷營養(yǎng)鹽含量有較大影響。無論是在現(xiàn)狀或未來的入庫污染負(fù)荷條件下，庫周清淤或全庫清淤方案實(shí)施后，于橋水庫水體的氮、磷營養(yǎng)物質(zhì)含量均會(huì)有所下降，但年內(nèi)變化過程基本不變。庫周清淤對(duì)沿岸帶淺水區(qū)水體氮、磷有很好的改善作用，但對(duì)庫區(qū)的整體改善效果不佳(濃度降低率小于2%)；而全庫清淤后，整個(gè)庫區(qū)水體的氮、磷濃度均明顯下降(濃度降低率達(dá)到26%～48%)，治理效果相對(duì)明顯。清淤后底泥營養(yǎng)鹽仍有可能隨著上覆水中營養(yǎng)鹽含量的上升而重復(fù)積累，使得新生底泥在較短時(shí)期內(nèi)表現(xiàn)出污染“源”的特征。因此，全面的水庫污染防治應(yīng)該綜合內(nèi)源、外源污染控制、改善生態(tài)結(jié)構(gòu)等措施共同開展。

(3)該模型的建立為生態(tài)清淤工程的效益評(píng)估提供了科學(xué)預(yù)測(cè)工具。在應(yīng)用模型分析清淤的環(huán)境生態(tài)效應(yīng)時(shí)，應(yīng)充分開展沉積物-水界面物質(zhì)交換實(shí)驗(yàn)獲取準(zhǔn)確的模型參數(shù)，提高模型預(yù)測(cè)分析的可靠性。