GWLF模型的改進及其在新安江流域的應用研究

齊作達，亢戈霖，王玉秋

(1.交通運輸部天津水運工程科學研究所，港口水工建筑技術國家工程實驗室，工程泥沙交通行業重點實驗室，天津 300456；2.南開大學 環境科學與工程學院 天津市跨介質復合污染環境治理技術重點實驗室，天津 300350)

1 研究背景

當今社會，水污染問題越來越受到人們的重視。近些年來，世界范圍內許多國家和政府逐步增加了對于水資源保護的投資。考慮到地表水監測人力、財力和物力的花費，非點源流域模型作為一種關鍵而有效的決策工具應運而生[1-3]。

非點源污染已經成為當前流域污染的重要組成部分，對于這類污染來源的控制與削減迫在眉睫[4-5]。與點源排放不同的是，非點源污染具有間接性的特征，并且與水文過程密切相關，因此識別起來更加困難和復雜[6-7]。為了能夠更加準確而簡潔地描述流域中水文循環和污染物歸趨，恰當的模型框架和適用性高的方程顯得尤其關鍵[8]。到目前為止，人們已經提出了許多能夠概化流域水文和污染物循環的方程，如SCS-CN徑流曲線方程、Green-Ampt方程、Muskingum洪水演進方程等[9-11]。

通用流域污染負荷方程GWLF(generalized watershed loading function)是一個集總參數流域負荷模型，被美國TMDL計劃所推薦和采納，并且在全世界各地得到了廣泛應用。模型在考慮一定的水文循環和污染物歸趨的基礎上，對于數據的需求量較小[12]，同時可以在月尺度上提供較為可靠的模擬結果，能夠為管理決策提供極大的便利和支持。近些年來，GWLF不僅在美國、歐洲等地區取得了較好的模擬結果[13-16]，在國內新安江流域、天津市于橋水庫上游流域、灤河流域等地區也有了諸多成功的應用案例[17-20]。

本次研究在水文框架改進的基礎上，進一步為GWLF模型添加了河道傳輸過程中泥沙的沉積、懸浮過程以及污染物的降解過程。在新安江流域的應用結果表明，改進版GWLF模型與原始GWLF模型相比，能夠更加準確地模擬流域內流量、泥沙量、溶解性氮負荷量和溶解性磷負荷量[21]。

2 GWLF原始模型及其改進方法

2.1 原始GWLF模型

GWLF模型是一個集總參數流域模型，模型需要輸入逐日的氣溫和降水數據，然后在日尺度上對流域的水文平衡進行連續模擬，最終累計到月尺度并輸出水量結果，泥沙和營養鹽則基于月尺度模擬并逐月輸出模擬結果。模型不考慮河道傳輸過程，并假設當月流域內產生的污染物全部傳輸到流域出口。水文產流計算采用了SCS-CN徑流曲線方程，泥沙計算公式采用了USLE通用土壤流失方程，營養鹽計算使用了經驗系數方程。更詳細的模型理論描述可以參考GWLF理論參考手冊。

2.2 GWLF模型的改進

2.2.1 GWLF模型的改進方法 GWLF模型的改進包括以下幾個方面：(1)將整個流域劃分為若干子流域，每個子流域內地表和地下的水文和營養鹽獨立計算，并匯入子流域對應的河道內；(2)模型的水文、泥沙、營養鹽基于日尺度計算；(3)添加河道傳輸過程，河道來水包括上游匯入、本子流域地表和地下水匯入以及點源的排放，河道的水文循環使用曼寧方程和變量存儲法進行描述，采用巴格諾德簡化方程表征泥沙在河道中的懸浮和沉積過程，營養鹽衰減使用與滯留時間相關的一階指數方程進行描述。改進后的模型框架示意圖見圖1。

圖1 改進后的模型框架示意圖

2.2.2 GWLF改進模型的水文、泥沙、營養鹽相關計算

(1)河道日平均流量計算

河道的日平均流量計算公式如下[22]：

(1)

將河道假設為河岸坡度為2∶1的梯形(見圖1)，河道中的水深可以用如下方程計算：

(2)

式中：H為河道水深，m，其下標j為迭代次數，H0= 0；Bi為河段i的河底寬度，m；ni為河段i的曼寧系數;Slopei為河段i的河岸坡度。

河段內的平均流速和傳輸時間可以使用如下方程表示：

(3)

TravelTimei,t=lengthi·1000/(Vi,t·3600)

(4)

式中：Vi,t為第t天河段i的平均流速，m/s；TravelTimei,t為第t天河段i的傳輸時間，h；lengthi為河段i的長度，km。

存儲系數SCi,t使用如下方程計算：

SCi,t=min(1,48/(2·TravelTimei,t+24))

(5)

最終河段的流出水量為FlowOuti,t(m3)：

(6)

(2)河道泥沙濃度計算

河道中的泥沙濃度計算公式如下：

(7)

式中：SedCi,t為第t天河段i的泥沙濃度，t/m3；SedSti,t為第t天河段i進水后的泥沙總量，t。

河道中的峰值流速計算公式如下：

velocityPki,t=pvaf·Vi,t

(8)

式中：velocityPki,t為第t天河段i的峰值流速，m/s；pvaf為峰值徑流調整系數。

河段最大可傳輸泥沙濃度計算公式如下：

SedCmaxi,t=lineari·velocityPki,texponenti

(9)

式中：SedCmaxi,t為第t天河段i的最大可傳輸泥沙濃度，t/m3；lineari和exponenti分別為線性系數和指數系數。

河段的泥沙沉積或懸浮量計算公式如下：

(10)

式中：SedDepi,t為第t天河段i的泥沙沉積量或泥沙懸浮量，t。

(3)河道中的營養鹽計算。

河道中的營養鹽衰減過程描述方程如下：

(11)

3 研究區域概況及數據來源

3.1 研究區域概況

本次研究選取新安江屯溪監測站上游流域作為研究區域，研究區域概況如圖2所示。所選研究區面積約為2 683 km2，是中國第一個跨省生態補償示范區，包括率水和橫江兩條主要支流(圖1(a))，長度分別為97 km和55 km，主河道總長385 km，下游匯入千島湖，為杭州市主要飲用水源地。流域被分為13個子流域，共有14個雨量監測站(圖1(b))。流域海拔在96～1 622 m之間(圖1(c))，氣候以亞熱帶濕潤季風氣候為主，1997-2013年期間年平均氣溫為17.2℃，年平均降水量為1 741 mm，大部分降水集中在3-8月。研究區域內的主要土地利用類型及面積占比為林地81%、耕地14.9%、草地3%(圖1(d))。

圖2 研究區域概況圖

3.2 數據來源

本次研究的14個降水站逐日數據和屯溪監測站的逐日流量、逐月泥沙數據來自《中華人民共和國水文年鑒》，逐日的氣溫數據來自中國氣象數據網，數字高程數據(DEM：30 m×30 m)來自地理空間數據云，土地利用數據來自國家第二次土地利用調查數據，逐月的溶解性氮和溶解性磷監測數據由黃山市環境監測站提供。模擬時間范圍為1997-2013年，其中1997-1999年的數據用于模型預熱。具體數據的格式和詳細內容參見表1。

表1 數據類型及來源

4 結果與分析

4.1 模型校準方法與結果討論

模擬值與實測值之間的比對評價，使用了3個國際公認的統計學指標，分別為納什系數(NSE)、決定系數(R2)和觀測標準偏差率(RSR)，其計算公式如下：

(12)

(13)

(14)

模型的校準方法采用了基于貝葉斯理論的廣義似然不確定性估計法(GLUE)，按照流量、泥沙、溶解性氮、溶解性磷的順序依次校準相關參數，綜合3項統計學指標選出最優參數組，校準后得到的參數數值如表2所示。

由表2可以看出，原始GWLF模型與改進模型校準后的各參數數值存在差異，但是差別較小，主要參數如徑流曲線系數、通用土壤流失方程系數、地表徑流氮、磷濃度和地下水氮、磷濃度差值保持在5%范圍以內，兩個模型的參數結果均在合理范圍之內，校準后的參數組可以用于模擬屯溪流域的流量、泥沙和營養鹽負荷量。

表2 原始GWLF與改進后的模型參數校準結果

原始GWLF模型與改進后模型的各項模擬結果的統計學指標如表3所示。

表3 模擬結果的各項統計學指標

由表3可知，改進模型的流量模擬統計學指標優于原始GWLF模型，NES、R2值均大于原始GWLF模型，盡管模擬逐日尺度流量結果不是GWLF模型原有的功能，但是改進模型的NES依然在0.9左右，表明改進模型在屯溪流域能夠提供較為精確的逐日流量模擬結果；原始模型泥沙的NES和R2均小于0.8，而改進后模型的NES和R2均大于0.9，同時改進模型的RSR值小于原始GWLF模型，表明改進模型在計算泥沙方面有所優化；溶解性氮的模擬結果統計學指標差異與泥沙相似，改進模型要明顯優于原始GWLF模型，但是溶解性磷的模擬結果統計學指標方面，改進模型與原始GWLF模型的差別較小，除RSR以外，另外兩項指標改進模型僅高出不到5%。

4.2 流量模擬結果

2000-2013年逐月流量實測值與模型模擬結果對比如圖3所示。從圖3整體上看，兩個模型的模擬計算結果令人滿意，幾乎再現了實測流量的變化趨勢，且模擬值與實測值的峰值保持一致。此外，改進模型在枯水期的模擬效果要明顯好于原始GWLF模型，原始GWLF模型傾向于低估枯水期的流量，這種差異在月尺度的模擬結果中并不是十分顯著。但是在日尺度上，二者的流量模擬結果差異較大，改進后的模型極大地提高了模擬的準確性。

圖3 2000-2013年逐月流量實測值與模型模擬結果對比

由于2000-2013年14 a的日尺度模擬結果時間序列圖清晰度較差，所以選擇以2005年為例，展示兩個模型之間逐日流量模擬結果的詳細差異，如圖4所示。由圖4可見，無降雨時期，兩個模型均無法完全重現流量的變化趨勢，但是原始GWLF模型過分地低估了流量，幾乎接近于零，并且過于平緩，對于較少降水的響應較差；而在大多數情況下，改進后的模型模擬結果與實測數據非常接近，并且基本捕捉了流量的變化趨勢。這種差異的重要原因是改進模型增加了河道過程，并引入了子流域概念。原始GWLF模型還存在高估暴雨事件流量峰值的現象，同時也會高估隨后的流量，這主要是因為沒有河道過程，缺少河道存水，當日的子流域產生的徑流全部轉化為流域出口的排水量。

圖4 2005年逐日流量實測值與模型模擬結果對比

4.3 泥沙量模擬結果

2000-2013年逐月泥沙量實測值與模擬結果對比如圖5所示。由圖5可看出，兩個模型對泥沙量的模擬效果較好，基本捕捉到了泥沙產量的波動性。但是原始GWLF模型對峰值模擬準確性較低，往往極大地高估或低估了極值，改進后的模型模擬值更加接近于實測值，其對峰值模擬結果有了明顯的改善。結合表3中的統計學指標可知，改進后GWLF模型對于泥沙的模擬效果要明顯優于原始GWLF模型。這種差距說明模型水文循環結構改進優化了泥沙模擬結果，增加河道泥沙懸浮沉積過程，使得改進后模型能夠對河道泥沙的攜帶能力做出相應調整，模擬效果更好，適應性更強。

圖5 2000-2013年逐月泥沙量實測值與模擬結果對比

4.4 溶解性氮模擬結果

2000-2013年逐月溶解性氮負荷量實測值與模擬結果對比如圖6所示。

圖6表明，與實測數據相比，原始GWLF模型和改進模型的模擬結果重現性較好，均能很好地再現實測數據變化趨勢。但原始GWLF模型對極值的模擬效果不甚理想，傾向于高估峰值而低估極小值，在枯水期模型響應較差，對低值的波動性捕捉較差。改進后的模型不僅優化了水文模擬過程，還新增了污染物河道衰減指數，對峰值和低值的模擬更加準確，波動極端性得到改善。同時改進后的模型模擬結果統計指標NES和R2均大于0.9，RSR小于0.3，明顯優于原始GWLF模型模擬結果。修正后的模型在一定程度上解決了原始GWLF模型對極值模擬的缺陷問題，模型模擬效果得到進一步優化。

4.5 溶解性磷模擬結果

2000-2013年逐月溶解性磷負荷量實測值與模擬結果對比如圖7所示。由圖7可見，雖然對溶解性磷的模擬結果與溶解性氮相比不甚理想，但兩個模型基本上可以反映溶解性磷的歷史變化趨勢。原始GWLF模型對峰值模擬往往偏低，平水期模擬偏高，枯水期時間序列過于平緩，無法準確捕捉到觀測數據的每一次波動。改進后模型對極值的模擬在一定程度上彌補了原模型的不足，模擬效果得到了一定的改善。這種差異可能來自于模型機理的改進：在水文過程優化的基礎上，將原始GWLF模型中城市營養鹽負荷為固態的假設改為城市總磷營養鹽負荷的1/4為溶解性磷[23]，并添加了污染物河道衰減過程。

圖7 2000-2013年逐月溶解性磷負荷量實測值與模擬結果對比

5 結 論

(1)在原始GWLF模型中加入以變量存儲法和污染物衰減為主要核心的污染物河道傳輸模塊，使GWLF模型可以進行多子流域水質模擬，從而拓寬模型使用流域范圍。

(2)在屯溪流域對原始GWLF模型和改進后的模型進行模擬性能比較，主要校準指標為流量、泥沙、溶解性氮和溶解性磷。通過對模型模擬結果進行詳細對比發現，改進后的模型明顯克服了原始模型對極值模擬存在較大偏差的缺點，顯著提高了模型模擬精度。

(3)改進后的模型具有模擬結果可靠、用戶友好、數據需求較少等優點，可作為我國流域科學管理的支持工具之一，為流域科學和精細化的水環境管理提供技術支持。