基于GWLF模型的于橋水庫流域氮磷負(fù)荷估算及來源變化解析

李澤利，羅娜，2，蘇德岳，高鍇*，趙興華，2，梅鵬蔚，張震

（1.天津市生態(tài)環(huán)境監(jiān)測中心，天津 300191；2.天津天濱瑞成環(huán)境技術(shù)工程有限公司，天津 300191；3.天津市引灤工程于橋水庫管理處，天津 301900）

于橋水庫位于天津北部薊州區(qū)城東4 km 處，是引灤入津輸水工程的重要調(diào)蓄水庫，同時(shí)也是天津市工農(nóng)業(yè)和生活用水的重要水源地[1]，自1983 年引灤入津以來一直為天津市人民生產(chǎn)和生活提供優(yōu)質(zhì)水源。然而自2014 年開始，于橋水庫富營養(yǎng)程度加劇[2]，2015—2017 年夏秋季庫區(qū)暴發(fā)大規(guī)模藍(lán)藻水華，嚴(yán)重威脅了天津市城市供水安全。氮、磷作為浮游植物賴以生長的重要營養(yǎng)物質(zhì)，參與光能轉(zhuǎn)化代謝過程，是導(dǎo)致湖庫富營養(yǎng)化最為重要的兩個(gè)因子[3-4]。由于氮、磷的輸入與流域內(nèi)人類活動密切相關(guān)，若要控制湖庫富營養(yǎng)化和藍(lán)藻水華，首先要明確進(jìn)入到水體中氮、磷營養(yǎng)物質(zhì)的來源及負(fù)荷[5]。于橋水庫的匯水主要來自引灤調(diào)水，沙河、黎河、淋河3 條入庫河流和庫周30 余條溝渠及大氣降水等。研究表明，引灤源頭潘家口、大黑汀水庫水質(zhì)惡化，以及流域氮磷面源污染加重，是于橋水庫水質(zhì)惡化、水生態(tài)平衡遭到破壞的重要原因[6-8]。為了保護(hù)和改善于橋水庫水質(zhì)，“水十條”實(shí)施以來，津冀兩?。ㄊ校┱扇《囗?xiàng)舉措來加強(qiáng)流域的保護(hù)力度[9]：①2016 年雙方簽訂《關(guān)于引灤入津上下游橫向生態(tài)補(bǔ)償?shù)膮f(xié)議》[7]，天津市根據(jù)跨界斷面水質(zhì)達(dá)標(biāo)情況撥付補(bǔ)償金；②2017 年河北省對潘家口、大黑汀水庫養(yǎng)魚網(wǎng)箱進(jìn)行全面拆除[10]；③天津市在入庫溝渠截污、農(nóng)村污水納管和農(nóng)業(yè)面源污染防治方面投入了大量資金，流域氮磷營養(yǎng)鹽輸入情況發(fā)生了較大變化。解析于橋水庫流域氮磷營養(yǎng)鹽來源、輸入負(fù)荷及變化趨勢，對區(qū)域水環(huán)境管理具有重要意義。

目前國內(nèi)外在面源污染負(fù)荷估算上，機(jī)理模型的使用具有主導(dǎo)地位[11]。趙磊[12]總結(jié)了國內(nèi)外應(yīng)用廣泛的GWLF、SWAT、AGNPS 等模型的主要優(yōu)缺點(diǎn)和適用范圍，表明GWLF 模型具有結(jié)構(gòu)簡單、參數(shù)易于獲取等特點(diǎn)，同時(shí)該模型的水文模塊原理與SWAT、AGNPS 相同，在我國北方灤河流域[13-14]、南方新安江流域[15]都取得了很好的應(yīng)用效果。本研究首先利用GWLF 模型的水文學(xué)模塊模擬于橋水庫各個(gè)子流域的降雨產(chǎn)流量，然后乘以子流域出口斷面的氮磷濃度，進(jìn)而估算氮磷面源污染負(fù)荷，最后針對2016—2017 年于橋水庫氮磷污染的來源、負(fù)荷和變化趨勢進(jìn)行分析，以期為流域水環(huán)境管理提供技術(shù)支持。

1 材料與方法

1.1 研究區(qū)概況

于橋水庫正常蓄水面積為86.8 km2，控制流域面積2 060 km2，其中上游河北省境內(nèi)面積占79%，下游天津市境內(nèi)面積占21%。流域?qū)儆跍貛О霛駶櫞箨懶约撅L(fēng)氣候，四季分明，多年平均氣溫12.5 ℃，平均降水量670 mm，主要集中在6—9 月。流域內(nèi)主要入庫河流為沙河、黎河（引灤調(diào)水通道）和淋河，此外庫周還有30余條入庫溝渠，多為季節(jié)性河流。

研究區(qū)內(nèi)設(shè)有4 個(gè)水文站，分別是大黑汀分水閘、前毛莊、水平口和淋河橋站，其中大黑汀分水閘站觀測引灤調(diào)水流量，前毛莊站觀測黎河自產(chǎn)流量與引灤調(diào)水流量，淋河橋站觀測淋河子流域自產(chǎn)流量，但受上游龍門口水庫蓄水影響長期斷流，水平口站觀測沙河流域自產(chǎn)流量，常年有水流動。研究區(qū)內(nèi)設(shè)有4個(gè)水質(zhì)監(jiān)測斷面，分別是隧洞出口、黎河橋、沙河橋和淋河橋，其中隧洞出口為大黑汀水庫經(jīng)引灤暗渠進(jìn)入于橋水庫流域的第一個(gè)監(jiān)測斷面，黎河橋、沙河橋和淋河橋分別為黎河、沙河、淋河的跨省界監(jiān)測斷面（圖1）。此外，庫周還設(shè)置38 個(gè)溝渠水質(zhì)觀測點(diǎn)[16]，逐月開展水質(zhì)監(jiān)測工作。

1.2 研究方法

于橋水庫流域氮磷污染來源從空間上可分為引灤調(diào)水、沙河、黎河、淋河、庫周地表水、庫周地下水和庫區(qū)大氣降水7 部分，從管理角度可匯總為引灤調(diào)水源、入境河流源（沙河、黎河和淋河）和境內(nèi)源（庫周地表水、地下水和大氣降水）3 大部分。用流量乘以氮磷濃度是計(jì)算負(fù)荷的最直接方法，但由于流域內(nèi)各條河流的水文站和水質(zhì)斷面的位置在空間上并不全部重合，且沙河、黎河下游和庫周均無流量監(jiān)測數(shù)據(jù)，無法直接計(jì)算。為此，本研究選取氣象、水文、下墊面數(shù)據(jù)比較齊全且受人類干擾較少的沙河水平口子流域，使用GWLF模型模擬產(chǎn)流過程并率定水文參數(shù)，然后將水文參數(shù)推及到整個(gè)流域[17]，模擬每個(gè)水質(zhì)斷面子流域產(chǎn)流量，再乘以相應(yīng)斷面的氮磷濃度，進(jìn)而估算整個(gè)流域的氮磷污染負(fù)荷。

圖1 于橋水庫流域河網(wǎng)及子流域示意圖Figure 1 Sketch map of river network and subbasins in Yuqiao reservoir watershed

于橋水庫流域水量和氮磷污染負(fù)荷估算所需的數(shù)據(jù)見表1。

1.2.1 GWLF模型原理

GWLF 模型是1987 年由美國賓夕法尼亞大學(xué)開發(fā)的一個(gè)半分布式、半經(jīng)驗(yàn)式的流域氮磷污染負(fù)荷模型[18]，該模型對地表徑流量的計(jì)算基于日氣象數(shù)據(jù)，使用徑流曲線數(shù)模型（SCS-CN）計(jì)算目標(biāo)流域內(nèi)不同土地利用類型的徑流深，作為表征地表徑流量的模擬結(jié)果。其中最重要的參數(shù)是表征地表產(chǎn)流能力的徑流曲線數(shù)（Curve number，CN）值，該值與土地利用類型、土壤類型密切相關(guān)。模型對地下潛流的計(jì)算基于日水量平衡，以淺層飽和區(qū)非重力自由水的遷移轉(zhuǎn)化為核心，將淺層飽和區(qū)看作一個(gè)簡單的線型水庫進(jìn)行模擬，通過兩個(gè)關(guān)鍵參數(shù)退水系數(shù)r（指淺層飽和區(qū)水體橫向轉(zhuǎn)移成為地下潛流水的比例）和滲漏系數(shù)s（指淺層飽和區(qū)水體縱向向下轉(zhuǎn)移進(jìn)入深層飽和區(qū)的比例）調(diào)整水量平衡。將逐日的地表徑流量與地下水潛流量相加，可得到逐日的河川徑流量，即到達(dá)流域出口的產(chǎn)流量[15]。

1.2.2 GWLF模型構(gòu)建

GWLF 模型產(chǎn)流模擬的輸入數(shù)據(jù)主要包括氣溫、降水、土地利用類型及面積等。首先基于ArcGIS 平臺加載研究區(qū)的數(shù)字高程圖（DEM），根據(jù)實(shí)際進(jìn)行河網(wǎng)矯正后，以各個(gè)水質(zhì)監(jiān)測斷面為出口分別生成沙河橋、黎河橋、淋河橋子流域和庫周子流域（圖1），然后與土地利用類型圖疊加，分別提取出每個(gè)子流域中各個(gè)土地利用類型的面積，作為模型的輸入數(shù)據(jù)；氣溫和降水?dāng)?shù)據(jù)以日為步長作為模型輸入數(shù)據(jù)，水平口水文站實(shí)測流量數(shù)據(jù)以月為步長作為模型校準(zhǔn)數(shù)據(jù)。與產(chǎn)流過程有關(guān)的參數(shù)，如各類土地利用類型的徑流CN值、地下水退水系數(shù)r、滲漏系數(shù)s和不同月份的蒸發(fā)覆蓋因子等，先采用模型推薦值輸入，然后進(jìn)行率定。

1.2.3 GWLF模型參數(shù)率定

利用沙河水平口水文站月觀測流量進(jìn)行水文參數(shù)率定，模擬時(shí)間跨度為2006年1月—2017年12月，共144 個(gè)月。使用美國康奈爾大學(xué)基于GWLF 模型原理開發(fā)的ReNuMa 軟件[15-16，19]進(jìn)行參數(shù)校準(zhǔn)，該軟件在Excel文件下運(yùn)行，利用了Excel強(qiáng)大的規(guī)劃求解宏，基于月實(shí)測數(shù)據(jù)校準(zhǔn)一個(gè)或多個(gè)參數(shù)，使得模擬結(jié)果與實(shí)測結(jié)果盡可能接近（誤差平方和最?。20]。將2006 年的數(shù)據(jù)作為模型“預(yù)熱期”，使模型初始條件穩(wěn)定平衡；2006—2013 年作為模型校準(zhǔn)期，校準(zhǔn)的參數(shù)主要有徑流曲線數(shù)CN 值、地下水退水系數(shù)r和滲漏系數(shù)s等；2014—2017 年作為驗(yàn)證期，驗(yàn)證模型參數(shù)的可信程度。采用Moriasi 等[21]推薦的統(tǒng)計(jì)指標(biāo)和評價(jià)標(biāo)準(zhǔn)判斷模型性能和適應(yīng)性，包括判定系數(shù)（R2）、納氏系數(shù)（NSE）、百分比偏差（PBIAS），其評價(jià)標(biāo)準(zhǔn)見表2，后兩者的計(jì)算公式如下：

式中：Qoi為觀測值；Qpi為模擬值為多年平均觀測值；n為樣本個(gè)數(shù)。

如表2 所示，水平口子流域徑流模擬除了驗(yàn)證期的NSE＜0.75 外，在全模擬期的R2、NSE 和PBIAS 分別為0.77、0.76和-3.4%，模型性能總體表現(xiàn)非常好。從圖2 的2006—2017 年的降水過程與徑流擬合效果圖來看，模擬流量和實(shí)測流量過程線吻合得較好，GWLF 模型適用于于橋水庫流域的徑流模擬。使用率定好的模型分別在沙河橋、黎河橋和庫周溝渠子流域建模，通過變換土地利用類型和面積等輸入數(shù)據(jù)模擬各個(gè)子流域出口的產(chǎn)流量，其中庫周溝渠子流域模擬出的地表徑流和地下潛流分開參與下一步的負(fù)荷量估算。

1.2.4 氮磷負(fù)荷估算

表2 基于月尺度徑流模擬的模型性能等級評價(jià)標(biāo)準(zhǔn)及評價(jià)結(jié)果Table 2 Evaluation standard of model performance levels and results based on monthly runoff simulation

圖2 沙河水平口子流域降水過程與徑流擬合效果圖Figure 2 Precipitation process and runoff fitting diagram of Shuipingkou subbasin

引灤調(diào)水、沙河、黎河、淋河的氮磷負(fù)荷分別根據(jù)大黑汀分水閘水量、沙河橋模擬水量、黎河橋模擬水量、淋河橋站實(shí)測水量乘以隧洞出口、沙河橋、黎河橋、淋河橋斷面對應(yīng)月份的總氮（TN）、總磷（TP）濃度計(jì)算得出。2017 年庫周地表水氮磷負(fù)荷估算已另文發(fā)表[16]，2016年采用相同方法進(jìn)行估算。庫周地下水中TN 濃度取自薊州區(qū)2019 年24 個(gè)千人以上地下水水源地水質(zhì)監(jiān)測結(jié)果，平均值為14.3 mg·L-1；TP 濃度取自薊州區(qū)2016—2017 年2 個(gè)國家考核地下水井監(jiān)測結(jié)果，平均值為0.04 mg·L-1。根據(jù)地下潛流水量模擬結(jié)果乘以上述氮磷濃度值，估算出庫周地下水輸入到于橋水庫的氮磷負(fù)荷。庫區(qū)降水輸入氮磷負(fù)荷采用降水量×水面面積×降水中氮磷濃度來估算，降水中TN 濃度根據(jù)2016 年天津市降水離子成分監(jiān)測結(jié)果折算得到：雨水中銨離子（以N 計(jì)）濃度年均值為2.62 mg·L-1，硝酸根離子（以N 計(jì)）濃度年均值為1.30 mg·L-1，因此TN（以N 計(jì)）年均濃度為3.92 mg·L-1；降水中TP 濃度無實(shí)測數(shù)據(jù)，根據(jù)太湖流域研究經(jīng)驗(yàn)[22]，按照0.05 mg·L-1計(jì)算。

2 結(jié)果與討論

2.1 氮磷負(fù)荷及來源構(gòu)成

表3 為2016—2017 年于橋水庫流域產(chǎn)流量、氮磷負(fù)荷估算結(jié)果。從表3 可以看出，2016 年于橋水庫流域產(chǎn)流量約4.6 億t，TN 輸入負(fù)荷2 185.9 t，TP 輸入負(fù)荷150.2 t；2017年產(chǎn)流量約4.1億t，TN 輸入負(fù)荷2 963.1 t，TP 輸入負(fù)荷39.0 t。徐寧等[23]利用SWAT 模型模擬2010—2013 年于橋水庫入庫污染負(fù)荷，結(jié)果顯示年均TN 入庫負(fù)荷為2 900.6 t，TP 負(fù)荷為163.2 t，本研究估算的TN、TP 負(fù)荷與其在量級上一致，表明GWLF模型適用于于橋水庫流域污染負(fù)荷研究。

表3 2016—2017年于橋水庫流域產(chǎn)流量、氮磷負(fù)荷估算結(jié)果Table 3 Estimation results of water yield and nitrogen and phosphorus load of Yuqiao reservoir from 2016 to 2017

將2016—2017 年于橋水庫流域產(chǎn)流量、氮磷負(fù)荷估算結(jié)果求平均值后進(jìn)行來源解析（圖3）。從圖3可以看出，產(chǎn)流量由大到小依次表現(xiàn)為引灤調(diào)水＞沙河＞黎河＞大氣降水＞庫周地表水＞淋河＞庫周地下水，分別占37.3%、18.3%、12.2%、12.0%、9.0%、7.7%和3.5%；TN負(fù)荷貢獻(xiàn)由大到小依次表現(xiàn)為引灤調(diào)水＞沙河＞黎河＞淋河＞庫周地表水＞庫周地下水＞大氣降水，分別占27.9%、25.3%、10.8%、10.7%、8.9%、8.4%和8.0%；TP負(fù)荷貢獻(xiàn)由大到小依次表現(xiàn)為引灤調(diào)水＞庫周地表水＞沙河＞黎河＞大氣降水＞淋河＞庫周地下水，分 別 占66.3%、13.0%、9.2%、7.0%、2.8%、1.1% 和0.6%。

圖3 不同來源產(chǎn)流量及氮磷負(fù)荷貢獻(xiàn)比例Figure 3 Contribution ratio of water yield and nitrogen and phosphorus load from different sources

從管理角度將氮磷負(fù)荷估算結(jié)果歸并為引灤調(diào)水源、入境河流源和境內(nèi)源3 部分，并與天津市水務(wù)局2010 年發(fā)布的《于橋水庫周邊水污染源近期治理工程實(shí)施方案》有關(guān)研究結(jié)果[24]進(jìn)行比較。從表4 可以得出，與2010年相比，2017年入境河流源TN、TP貢獻(xiàn)比例分別下降10.2 個(gè)百分點(diǎn)和9.2 個(gè)百分點(diǎn)，但總體貢獻(xiàn)均在30%以上；引灤調(diào)水源TN、TP 貢獻(xiàn)比例分別上升2.9 個(gè)百分點(diǎn)和10.6 個(gè)百分點(diǎn)，總體貢獻(xiàn)均在30%以上；境內(nèi)源TN 貢獻(xiàn)比例上升7.3 個(gè)百分點(diǎn)，TP 貢獻(xiàn)比例下降1.5 個(gè)百分點(diǎn)。值得注意的是，2016年引灤調(diào)水TP 貢獻(xiàn)比例高達(dá)75.4%，還需進(jìn)一步開展研究分析。

表4 基于管理角度的于橋水庫氮磷負(fù)荷來源比較Table 4 Comparison of nitrogen and phosphorus load based on management perspective of Yuqiao reservoir

2.2 氮磷來源變化分析

氮磷負(fù)荷受水質(zhì)和水量共同影響。由于在估算庫周地表水、地下水和大氣降水等境內(nèi)源氮磷負(fù)荷時(shí)，2016年和2017年氮磷濃度數(shù)據(jù)相同，其負(fù)荷量變化主要受年降水量影響，2016 年流域平均降水量為675.6 mm，2017 年為637.7 mm，降幅為5.6%，2017 年境內(nèi)源氮磷污染負(fù)荷均低于2016年。對2016—2017年引灤調(diào)水、沙河、黎河、淋河輸入的氮磷負(fù)荷變化進(jìn)行分析如下。

2.2.1 引灤調(diào)水污染負(fù)荷變化

引灤調(diào)水水質(zhì)和水量共同決定調(diào)水污染負(fù)荷。從圖4 可以看出，引灤調(diào)水TP 濃度由2016 年初最高的0.74 mg·L-1明顯下降到2017 年的0.09 mg·L-1，降幅87.8%。同時(shí)，調(diào)水水量由2016 年（1 月、5 月、6 月之和）的1.91 億m3下降到2017 年（6 月、8 月、9 月之和）的1.34 億m3，下降了29.8%。水量和濃度的共同下降使得TP 負(fù)荷由2016 年的113.3 t 下降到2017 年的12.1 t，降幅高達(dá)89.3%，這也是表4 中引灤調(diào)水貢獻(xiàn)比例由75.4%下降到31.1%的主要原因。資料顯示2016 年11 月—2017 年5 月，國家、河北省、天津市財(cái)政共同出資清理潘家口和大黑汀水庫的養(yǎng)魚網(wǎng)箱，共清理拆解網(wǎng)箱79 687個(gè)，拆除后大黑汀水庫總磷濃度明顯下降[10]。

圖4 2016—2017年引灤調(diào)水量和污染物濃度變化Figure 4 Comparison of water diversion and pollutant concentration in Luanhe River from 2016 to 2017

然而，引灤調(diào)水TN 負(fù)荷不降反升。監(jiān)測結(jié)果顯示引灤調(diào)水TN 濃度由2016 年初的2.18 mg·L-1上升到2017 年的7.39 mg·L-1，增幅239%，這可能是由于調(diào)水期8—9 月正值汛期，此時(shí)段灤河流域面源污染較重[25-26]。盡管2017 年調(diào)水水量較2016 年下降了29.8%，但調(diào)水中TN 濃度的升幅較調(diào)水水量的降幅變化更大，導(dǎo)致2017年引灤調(diào)水TN 輸入負(fù)荷較2016年上升了123.3%。

2.2.2 入境河流污染負(fù)荷變化

沙河、黎河、淋河子流域的面源污染物通過降水產(chǎn)流過程進(jìn)入于橋水庫。由于三個(gè)流域具有相似的面源污染特點(diǎn)，現(xiàn)以沙河子流域?yàn)槔治鯰N、TP 負(fù)荷變化。沙河子流域2017年產(chǎn)流量較2016年上升了11.4%。從圖5 可以看出，2017 年沙河各月TP 濃度均低于2016 年，平均濃度由0.238 mg·L-1下降到0.069 mg·L-1，下降了71.0%，由于TP 濃度的降幅較產(chǎn)流量的升幅變化更大，相乘后使得2017 年TP 負(fù)荷較2016年下降了51.9%。2017 年沙河各月TN 濃度與2016年相比互有高低，平均濃度由9.95 mg·L-1下降到9.29 mg·L-1，下降了6.6%，但由于2017 年產(chǎn)流量的升幅較TN 濃度的降幅變化更大，相乘后使得2017 年TN負(fù)荷較2016年上升了20.5%。

圖5 2016—2017年沙河子流域產(chǎn)流量及氮磷濃度變化Figure 5 Comparison of water yield and TN，TP concentration in Shahe River from 2016 to 2017

2016 年津冀兩?。ㄊ校┱炗啞蛾P(guān)于引灤入津上下游橫向生態(tài)補(bǔ)償?shù)膮f(xié)議》以來，河北省加大了沙河、黎河、淋河流域污染治理力度，河道水質(zhì)明顯改善，TP 等五項(xiàng)考核指標(biāo)年均值均滿足《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》（GB 3838—2002）Ⅲ類水質(zhì)要求[7]。但由于TN 濃度尚未納入考核體系，關(guān)注程度和治理力度不大，流域內(nèi)農(nóng)業(yè)種植[27]、畜禽和水產(chǎn)養(yǎng)殖、農(nóng)村生活污水排放的氮素很容易通過水體流動輸入到于橋水庫，為藻類大量生長提供充足的營養(yǎng)物質(zhì)。

3 結(jié)論

（1）GWLF 模型對于橋水庫流域水文模擬效果較好。2016—2017 年于橋水庫流域年均產(chǎn)流量約4.4億m3，其中引灤調(diào)水貢獻(xiàn)占比37.3%，3條入境河流貢獻(xiàn)38.2%，庫周產(chǎn)流和大氣降水共占24.5%；年均總氮負(fù)荷2 574.6 t，其中引灤調(diào)水貢獻(xiàn)占比27.9%，3 條入境河流貢獻(xiàn)46.8%，庫周產(chǎn)流和大氣降水共占25.3%；年均總磷負(fù)荷94.6 t，其中引灤調(diào)水貢獻(xiàn)占比66.3%，3 條入境河流貢獻(xiàn)17.3%，庫周產(chǎn)流和大氣降水共占16.4%。

（2）與2016年相比，2017年于橋水庫流域總磷負(fù)荷減少了74.0%，主要是由于引灤調(diào)水和入境河流總磷濃度均明顯下降（降幅分別為87.8%和71.0%）。然而，全流域總氮負(fù)荷卻增加了35.6%，這是由引灤調(diào)水總氮濃度大幅上升（增幅239%）以及入境河流產(chǎn)流量增加（增幅11.4%）導(dǎo)致。

（3）于橋水庫流域總磷污染防控取得了一定效果，但總氮污染有加重的趨勢，需要加強(qiáng)流域非點(diǎn)源氮污染防治力度，同時(shí)選擇大黑汀水庫水質(zhì)較好的月份調(diào)水，以減少氮污染負(fù)荷。