徐婉珍

摘要 基于EFDC建立潘家口水庫水動力與富營養(yǎng)化模型，模擬了2006～2007年潘家口水庫內的水動力、營養(yǎng)物質循環(huán)以及藻類生長代謝，并對模型進行了校準與驗證。結果表明，TN的模擬效果最好，最大相對誤差為4.97%；TP的模擬效果次之，最大相對誤差為-15.77%；而Chla的模擬效果最差，最大相對誤差18.29%。模擬與實測結果對比表明，基于EFDC建立的潘家口水庫富營養(yǎng)化模型，可以比較真實地反映庫區(qū)內營養(yǎng)物質的循環(huán)以及藻類的生長代謝過程。

關鍵詞水動力；水質；富營養(yǎng)化；潘家口水庫；EFDC

中圖分類號S181文獻標識碼A文章編號0517-6611（2015）34-115-04

潘家口水庫位于海河流域灤河水系的灤河干流上，地跨河北省的承德、唐山兩市，是天津市、唐山市重要的水源地，承擔著提供兩地工業(yè)、農業(yè)、城市用水的重任。截至2013年，潘家口水庫已經向天津市、唐山市以及灤下灌區(qū)提供優(yōu)質水資源近366億m3，其中向天津市供水167億m3，向唐山以及灤下灌區(qū)供水199億m3，成為天津市和唐山市可靠的供水水源，為天津市和唐山市的可持續(xù)發(fā)展做出了重要貢獻[1]。

隨著區(qū)域經濟的發(fā)展，人民生活水平的提高，潘家口水庫所在流域的點源污染與面源污染的輸入越來越被關注與重視。上游地區(qū)城市生活污水對水源地造成嚴重污染。農作物使用的大量農藥、化肥，隨著上游的來水直接進入庫區(qū)，造成庫區(qū)氮、磷含量明顯增加，富營養(yǎng)化程度升高。水庫周邊大量的工廠也導致了大量污染物的輸入。同時，網箱養(yǎng)魚成為重要的內源，造成的威脅不容小覷，在經濟利益的驅使下，潘家口水網箱養(yǎng)魚從2003年至今呈現(xiàn)快速擴大的趨勢，養(yǎng)殖密度幾乎超過了水庫的承載能力。投放的餌料以及產生的魚類糞便中含有大量的氮、磷，加速了水庫的富營養(yǎng)化進程[2]。

水庫作為調節(jié)、蓄水和控制水資源的人造系統(tǒng)，是實現(xiàn)水資源合理開發(fā)利用的有效途徑，而水庫相對封閉，水體停留時間較長，營養(yǎng)鹽循環(huán)受到干擾，大量的氮磷元素排入水體，引起水體中藻類、浮游生物的迅速繁殖，從而導致水體富營養(yǎng)化，嚴重影響居民的飲水安全[3]。為明晰潘家口水庫內的富營養(yǎng)化過程，筆者利用EFDC模型建立潘家口水庫富營養(yǎng)化模型，對潘家口水庫內的水動力、營養(yǎng)物質的循環(huán)以及藻類的生長代謝做了較系統(tǒng)的模擬。

1模型簡介

EFDC（The Environmental Fluid Dynamics Code）模型是美國國家環(huán)保署（USEPA）推薦的三維地表水水動力模型[4-5]，可實現(xiàn)河流、湖泊、水庫、濕地系統(tǒng)、河口和海洋等水體的水動力學和水質模擬，是一個多參數(shù)有限差分模型。EFDC模型采用MellorYamada 2.5階紊流閉合方程，根據需要可以分別進行一維、二維和三維計算。模型包括水動力、水質、有毒物質、底質、風浪和泥沙模塊，用于模擬水系統(tǒng)一維、二維和三維流場、物質輸運（包括水溫、鹽分、粘性和非粘性泥沙的輸運）、生態(tài)過程及淡水入流，可以通過控制輸入文件進行不同模塊的模擬。模型在水平方向采用直角坐標或正交曲線坐標，垂直方向采用σ坐標變換，可以較好地擬合固定岸邊界和底部地形。在水動力計算方面，動力學方程采用有限差分法求解，水平方向采用交錯網格離散，時間積分采用二階精度的有限差分法，以及內外模式分裂技術，即采用剪切應力或斜壓力的內部模塊和自由表面重力波或正壓力的外模塊分開計算。外模塊采用半隱式三層時間格式計算，因傳播速度快，所以允許較小的時間步長。內模塊采用考慮了垂直擴散的隱式格式，傳播速度慢，允許較大的時間步長，其在干濕交替帶區(qū)域采用干濕網格技術。該模型提供源程序，可根據需要對源程序進行修改，從而達到最佳的模擬效果。

水動力學方程是基于三維不可壓縮的、變密度紊流邊界層方程組，為了便于處理由密度差而引起的浮升力項，常常采用Boussinesq假設。在水平方向上采用曲線正交坐標變換，在垂直方向上采用sigma坐標變換后，得到以下控制方程：

動量方程為：

（mHu）t+（myHuu）x+（mxHvu）y+

（mwu）z-（mf+vmyx-umxy）Hv=

-myH（gζ+p）x-my（hx-zHx）pz+z（m1HAvuz）+Qu（1）

（mHv）t+（myHuv）x+（mxHvv）y+

（mwv）z+（mf+vmyx-umxy）Hu=

-mxH（gζ+p）y-mx（hy-zHy）pz+z（m1HAuvz）+Qv（2）

pz=-gHρ-ρ0ρ0=-gHb（3）

連續(xù)方程為：

（mζ）t+（myHu）x+（mxHv）y+（mw）z=0

（4）

（mζ）t+（myH∫10udz）x+（mxH∫10vdz）y=0

（5）

ρ=ρ（P，S，T）（6）

EFDC模型中的富營養(yǎng)化模塊模擬了21種物質的狀態(tài)變量，包括3種藻類中碳、氮、磷組分、硅循環(huán)、溶解氧以及糞便大腸桿菌。水體中的污染物控制方程由溶于水中的污染物、吸附于水中可溶性物質上的污染物以及吸附于水中懸浮顆粒物上的污染物共同構成。溶于水中的污染物遷移控制方程如下：

t（mxmyHCw）+x（myHuCw）+y（mxHvCw）+z（mxmywCw）

=z（mxmyAbHzCw）+mxmyH（ikidSSiχis+jKjdDDjχjD）-

mxmyH[i（KiaSSi）（ψwCw）（iS+χiS）+j（KjaDDj）（ψwCw）（jD-χjD）+γCw]

（7）

式中，Cw為溶解于水中污染物的單位體積質量濃度；χS為吸附于i類懸浮物的污染物單位體積質量濃度；χD為吸附于j類溶解物質的單位體積質量濃度；為孔隙率；ψw為可吸附溶解于水中污染物的水所占比例；Ka為吸附速率；Kd為解吸速率；γ為凈線性衰減系數(shù)；S為水中懸浮物單位體積質量濃度；D為水中溶解物質單位體積質量濃度。

EFDC的富營養(yǎng)化模塊包含3種藻類：藍藻、綠藻和硅藻，描述這些藻類動力的過程大致相同，只是在不同藻類的方程中改變了參數(shù)值。模塊中的源匯項主要是藻類的生長、新陳代謝、捕食、沉積和外來負荷，方程如下：

tBc=（Pc-BMc-PRc）Bc+Z（WScBc）+WBcV（8）

式中，Bc為藻類生物量（g C/m3）；t為時間（d）；Pc為藻類的生產率（d-1）；BMc為藻類新陳代謝速率（d-1）；PRc為藻類捕食率（d-1）；WSc為沉積速率（m/d）；WBc為藻類的外源輸入（g C/d）；V為網格單元體積（m3）。

2模型建立與數(shù)據來源

2.1模型建立

清河口斷面至潘家口水庫大壩區(qū)域水面開闊、水深較大，污染物濃度在水平及垂直方向上均有較大的差別，為了能充分揭示污染物在潘家口水庫的遷移擴散規(guī)律，需要進行三維模擬。其中，水平方向上劃分為637個單元格，垂直方向上平均劃分為4層，共計2 916個單元格（圖1）。模型共有5處入流，分別為清河口、瀑河口、石泉浩、溝門以及菜子溝門，瀑河口與菜子溝門有水質點源輸入，除此外的入流根據水質監(jiān)測數(shù)據輸入水質組分的濃度時間序列。校核點為燕子峪與潘家口兩處。水庫出流為潘家口水庫大壩，采用上層出流。

圖1潘家口水庫EFDC模型邊界條件與校核點示意圖

2.2數(shù)據來源

除2006～2007年潘家口水庫庫區(qū)的氣象資料來源于氣象資料分享網站和流域年鑒外，

潘家口水庫庫底地形圖、水文測站基本信息、2006～2007年潘家口水庫水文監(jiān)測資料（水位、流量）、水質監(jiān)測數(shù)據和點源入流的數(shù)據均來源于水利部海河水利委員會。

3模型的校核與驗證

水動力的模擬是潘家口富營養(yǎng)化模型的基礎。水動力的模擬決定了物質運輸?shù)穆窂揭约拔镔|混合的方式，這是模型首先要評價的初步內容。因為缺乏流速資料，僅以潘家口水庫大壩壩上水位進行水動力的校驗。結果表明，模擬值與實測值吻合良好（圖2），平均誤差為0.056 m。

EFDC的富營養(yǎng)化模擬了藍藻、綠藻和硅藻在太陽輻射、水溫、營養(yǎng)鹽水平影響下的生長、代謝、死亡。每種藻類都有其自己的生長、演替規(guī)律，藻動力模型通過改變與生產、代謝、布施、沉積速率等相關的關鍵參數(shù)來模擬不同的過程。細胞個數(shù)和葉綠素濃度這兩個指標可以描述水體中的藻類生物量。細胞計數(shù)、染色劑濃度都是計量藻類生物量的通用方法。直接的細胞數(shù)計量與體積計量十分費力，需要專門的生物分類以及樣品保存。單獨某種藻類的葉綠素濃度測算也十分復雜。但是，水體中的葉綠素濃度可以直接指示生態(tài)系統(tǒng)中的生物量與富營養(yǎng)化水平，因此模型采用葉綠素濃度代表總藻類生物量。

富營養(yǎng)化模型的校準關鍵在于藻類初級生產力、有機物分解、營養(yǎng)鹽濃度的模擬。模型校準的目標是選擇與藻類生長、營養(yǎng)鹽轉化、呼吸過程相關的最合適的參數(shù)組合，減小模擬值與實際值的誤差。對浮游生物體以及營養(yǎng)鹽輸移轉化過程的理解，有助于提高富營養(yǎng)化過程模擬的準確性與精確性。模型初步運行應用相關文獻推薦的經驗參數(shù)，逐漸調整個別參數(shù)，讓模擬值與實測值間接近吻合。模型的驗證是用另外一套獨立的數(shù)據，在不同的外界條件下檢驗模型反映水環(huán)境系統(tǒng)變化的能力。藻類動力過程的關鍵參數(shù)包括生長速率、新陳代謝速率以及捕食率，這3個參數(shù)描述了藻類生物量的增加與減少。極端的高溫與低溫會嚴重影響藻類濃度。此外，營養(yǎng)鹽濃度也會影響藻類的生長速率。富營養(yǎng)化模型中的各個狀態(tài)變量也取決于硝化過程、反硝化過程以及有機物的水解過程。

在模型校驗之后進行敏感性測試，可以加深對模型表現(xiàn)

的理解，并檢驗模型設置的關鍵參數(shù)。結果表明，藻類生產率是最敏感的內部參數(shù)，磷的半飽和濃度次之。敏感性分析顯示，10%藻類生產率的變化可以導致大于43.5%的葉綠素濃度變化。藻類生產濃度作為內部參數(shù)，取決于藻的種類及其特征，而磷的半飽和濃度反映了外界營養(yǎng)鹽負荷對藻生物量的影響。在模擬中，TP負荷發(fā)生10%的擾動，會導致葉綠素濃度至少35.7%的變幅；TN發(fā)生10%的擾動，導致葉綠素濃度大約13.2%的變幅，所以可以認為在模型中磷負荷的輸入是影響藻類生長最敏感的外因。

在模型的校準和驗證中，將模擬值與實測值相比較，計算平均誤差（ME）與平均相對誤差（MRE），從而評價模擬的效果。具體計算公式為：

ME=1nni=1（Oi-Pi）（9）

MRE=1nni=1Oi-PiOi×100%（10）

式中，Oi為實測數(shù)據；Pi為模擬數(shù)據；n為校準與驗證數(shù)據的數(shù)目。校核與驗證結果表明，TN的模擬效果最好，TP次之，Chla的模擬偏差最大。對于Chla，在校準時期燕子峪的最大相對誤差達18.29%，在驗證時期潘家口的最大相對誤差達18.13%；TP模擬的最大相對誤差在校準時期與驗證時期分別達-15.77%、15.52%；TN模擬效果最好，最大相對誤差出現(xiàn)在驗證時期的燕子峪，達4.97%（表1）。TN、TP、Chla模擬的相對誤差都控制在20%以內。模擬結果顯示，該模型可以描述潘家口水庫表層的藻類生長。在校準、驗證中，營養(yǎng)鹽循環(huán)的模擬有著更好的吻合度。

由圖3～6可知，2006和2007年Chla模擬值與實測值呈現(xiàn)較好的一致性，在120 d左右，低溫環(huán)境逐漸改變，藻類生命活動逐漸活躍起來，隨著氣溫回暖，生長速率、代謝速率、繁殖速率都逐漸上升，180～240 d Chla濃度達到峰值。Chla濃度的變化過程主要依賴于水溫的變化，180～240 d是北方的夏季，溫度為25～30 ℃，正好藍藻、綠藻、硅藻的最適溫度分別為27.5、20、25 ℃左右；同時，充足的太陽輻射為藻類的光合作用提供了足夠的能量。該段時間內，Chla濃度呈現(xiàn)波動變化，原因是7、8月份正是北方雨水充沛的時節(jié)，流域上游的大量來水一方面會在短時間內稀釋Chla，但又因為大量營養(yǎng)物的輸入，最終又促進了藻類的生長，同時大壩放水，也對庫內的Chla濃度造成了影響。此外，庫內該段時期內浮游動物與魚類的大量捕食，也是造成Chla濃度波動的原因之一。在240 d之后，溫度下降是造成Chla濃度下降的主要原因。

營養(yǎng)物質的濃度也受到藻類生長的影響，TN、TP濃度會在120 d左右出現(xiàn)明顯下降，因為此時溫度升高，藻類生長，開始吸收營養(yǎng)物質，導致在該時間段會出現(xiàn)氮、磷濃度的明顯下降。但隨著溫度逐漸升高，藻類生長逐漸旺盛，TN、TP濃度卻并沒有一直下降，這是因為上游來水也逐漸增大，同時帶來流域內大量的營養(yǎng)物質，隨著雨季的到來，氮、磷的濃度在7、8月份達到最高峰，因為此段時間常伴隨著水庫階段性大量的放水，導致了庫內TN、TP濃度的波動。240 d后雨季結束，氮、磷營養(yǎng)物質的輸入大幅度下降，導致240 d左右TN、TP的濃度又一次明顯的下降，而300 d后又普遍呈上升的趨勢，這是由浮游植物的死亡與分解引起的。而在0～60與330～360 d大約為期3月的時間內，庫區(qū)冰封，既沒有明顯影響庫內營養(yǎng)物質濃度的水量、營養(yǎng)物質輸入，也沒有浮游植物的生長代謝，也沒有水庫的開閘放水，處于非常平穩(wěn)的狀態(tài)，變化微小，呈緩慢下降的趨勢，可能是因為吸附在顆粒物上的營養(yǎng)物質隨著顆粒的沉淀而進入了水底沉積物中，脫離了水體。