灤河潘家口水庫(kù)控制流域水量與污染物負(fù)荷模型構(gòu)建及率定

劉 超，李建柱，徐彬鑫

(1.河北省邢臺(tái)水文勘測(cè)研究中心，河北 邢臺(tái) 054000；2.天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350；3.水利部海河水利委員會(huì)，天津 300170)

水資源是人類(lèi)生活、生產(chǎn)和生態(tài)環(huán)境的重要保障，中國(guó)的水資源僅為世界平均水平的1/4，且降水的多年變化和年內(nèi)集中程度均很大[1]，徑流多以洪水的形式發(fā)生，較大一部分徑流不能被利用。隨著社會(huì)經(jīng)濟(jì)的發(fā)展，點(diǎn)源污染和非點(diǎn)源污染導(dǎo)致的水環(huán)境問(wèn)題也逐漸凸顯，進(jìn)一步加劇了我國(guó)水資源的供需矛盾。由于治理能力的不斷提升與管理措施的不斷嚴(yán)格，點(diǎn)源污染基本得到控制，非點(diǎn)源污染已成為造成水體污染的主要原因。

流域作為自然形成的水資源動(dòng)態(tài)系統(tǒng)，是水循環(huán)過(guò)程及水行政管理中的重要單元。流域問(wèn)題的研究，是水資源開(kāi)發(fā)利用，合理配置優(yōu)化，綜合決策管理等宏觀問(wèn)題研究的落腳點(diǎn)[2]。首要任務(wù)是獲取流域詳實(shí)的水文、氣象信息，但受限于地理?xiàng)l件與歷史因素，同時(shí)考慮到經(jīng)濟(jì)成本，觀測(cè)數(shù)據(jù)往往集中在流域的中下游區(qū)域[3]，且在時(shí)間上也不具有連續(xù)性，客觀上限制了流域決策管理的精細(xì)化。因此，流域水文、水動(dòng)力、水質(zhì)模型成為補(bǔ)充缺測(cè)數(shù)據(jù)的重要手段，以及流域水資源綜合管理的有效工具。

潘家口水庫(kù)是天津市、唐山市重要的水源地，承擔(dān)著天津、唐山工業(yè)、農(nóng)業(yè)、城市生活和環(huán)境用水的重任[4]。截至2018年底，累計(jì)向天津市供水187億 m3，為天津市經(jīng)濟(jì)可持續(xù)發(fā)展做出了重要貢獻(xiàn)。2000年之前，潘家口水庫(kù)各斷面水質(zhì)均為III類(lèi)水標(biāo)準(zhǔn)，近年來(lái)，隨著社會(huì)經(jīng)濟(jì)的快速發(fā)展，潘家口水庫(kù)上游工業(yè)和城市生活污水、農(nóng)業(yè)生產(chǎn)中大量農(nóng)藥化肥的施用、水庫(kù)網(wǎng)箱養(yǎng)魚(yú)及水庫(kù)周邊旅游設(shè)施造成的污染日益嚴(yán)重。2016-2017年，潘家口水庫(kù)由于水質(zhì)較差，天津中斷了引灤供水。因此，構(gòu)建潘家口水庫(kù)控制流域污染物負(fù)荷模型，模擬分析潘家口水庫(kù)的污染物負(fù)荷，對(duì)于精準(zhǔn)分析污染物來(lái)源，并提出針對(duì)性措施和對(duì)策可提供科學(xué)依據(jù)。

1 研究區(qū)概況

灤河發(fā)源于河北省承德市境內(nèi)的巴彥古爾圖山北麓，流經(jīng)河北省、內(nèi)蒙古自治區(qū)及遼寧省等地區(qū)，最終于河北省樂(lè)亭縣匯入渤海[5]。灤河流域位于39°10′N(xiāo)～42°30′N(xiāo)，115°30′E～119°15′E，控制面積為44 750 km2。潘家口水庫(kù)是灤河上的大型水庫(kù)，是引灤入津的水源地，目前承擔(dān)著保障天津市河道生態(tài)需水的重要任務(wù)，流域位置及水系情況如圖1。

根據(jù)2018年國(guó)控?cái)嗝嫦嚓P(guān)河湖水質(zhì)監(jiān)測(cè)資料，灤河流域重要斷面水質(zhì)為Ⅲ～劣Ⅴ類(lèi)。主要是一些月份的水質(zhì)情況不滿(mǎn)足各自水環(huán)境控制單元與水功能區(qū)的控制目標(biāo)。其中，瀑河大桑園控制斷面的水質(zhì)最差時(shí)為劣Ⅴ類(lèi)，化學(xué)需氧量(COD)、五日生化需氧量(BOD)等指標(biāo)明顯超標(biāo)；灤河干流的上板城大橋控制斷面水質(zhì)為Ⅴ類(lèi)，高錳酸鹽指數(shù)、化學(xué)需氧量等指標(biāo)超標(biāo)。研究流域內(nèi)水質(zhì)控制斷面2018年水質(zhì)情況見(jiàn)表1。

圖1 研究區(qū)域地理位置及水系圖

2 流域污染物負(fù)荷模型構(gòu)建

根據(jù)流域水文氣象、土地利用及污染物來(lái)源等數(shù)據(jù)，構(gòu)建SWAT水文水質(zhì)模型對(duì)流域徑流過(guò)程和污染物負(fù)荷進(jìn)行模擬分析。

2.1 SWAT模型原理

SWAT模型(Soil and Water Assessment Tool)是美國(guó)農(nóng)業(yè)部開(kāi)發(fā)的流域尺度分布式水文模型，模型可用于模擬單一流域或多個(gè)具有水文聯(lián)系流域的水文過(guò)程[6]。SWAT模型內(nèi)包含有土地利用、地下水、流域管理、河道水質(zhì)情況等多個(gè)模塊，能夠預(yù)測(cè)長(zhǎng)時(shí)間序列下流域的水文、泥沙與營(yíng)養(yǎng)物質(zhì)的變化情況[7]。

表1 研究流域重要水域控制斷面水質(zhì)

SWAT模擬的水文循環(huán)基于水量平衡方程

(1)

式中：SWt為最終土壤水含量，mm；SW0為第i天初始土壤含水量，mm；t為天數(shù)，d；Rday為第i天降水量，mm；Qsurf為第i天地表徑流，mm；Ea為第i天蒸發(fā)量，mm；Wseep為第i天土壤底層滲透量，mm；Qgw為第i天地下水出流，mm。

SWAT在降雨徑流計(jì)算中則采用SCS曲線(xiàn)法，該方程自上世紀(jì)50年代起開(kāi)始得到廣泛應(yīng)用，其方程式為：

(2)

式中：Qsurf為地表徑流量，mm；Rday為降雨量，mm；Ia為初損量，mm；S為滯留系數(shù)，mm；表達(dá)式如下：

(3)

式中：CN為模擬日的曲線(xiàn)數(shù)，初損Ia通常估計(jì)為0.2S。

表2 SWAT模型建模所需基礎(chǔ)數(shù)據(jù)

2.2 SWAT模型構(gòu)建

SWAT模型在構(gòu)建過(guò)程中需要輸入數(shù)字高程數(shù)據(jù)、土壤數(shù)據(jù)、土地利用數(shù)據(jù)與氣象數(shù)據(jù)(降水、氣溫、太陽(yáng)輻射、風(fēng)速及相對(duì)濕度等)，具體輸入數(shù)據(jù)見(jiàn)表2。在模型建立后，還需要利用實(shí)測(cè)徑流數(shù)據(jù)對(duì)模型進(jìn)行率定。具體建模過(guò)程如下：首先通過(guò)八流向法根據(jù)DEM與設(shè)定的集水面積閾值生成流域內(nèi)河網(wǎng)，并劃定出若干子流域；其次，定義流域內(nèi)土地利用情況、土壤類(lèi)型與坡度分級(jí)，通過(guò)設(shè)定各自對(duì)應(yīng)的閾值，進(jìn)而形成閾值組合以在各子流域內(nèi)生成水文響應(yīng)單元(HRUs)；之后在模型中加載氣象數(shù)據(jù)并運(yùn)行，SWAT以HRU為最小計(jì)算單元各自獨(dú)立進(jìn)行水文循環(huán)過(guò)程模擬，之后通過(guò)SCS曲線(xiàn)法計(jì)算得到地表徑流，河道水流演算則采用馬斯京根法。

2.2.1 數(shù)字高程(DEM)數(shù)據(jù)處理

潘家口水庫(kù)控制流域DEM由NASA提供的SRTM的DEM數(shù)據(jù)裁剪得到，分辨率為30 m×30 m(如圖2)，將研究區(qū)域劃分了34個(gè)子流域。

圖2 研究區(qū)DEM及子流域劃分情況

2.2.2 土地利用及土壤數(shù)據(jù)

選用2005年土地利用數(shù)據(jù)，并將土地利用數(shù)據(jù)進(jìn)行重分類(lèi)，分為林地、草地、耕地、水域、建設(shè)用地及未利用土地六大類(lèi)，

土壤數(shù)據(jù)資料來(lái)源于HWSD土壤數(shù)據(jù)庫(kù)(http://www.fao.org)，該數(shù)據(jù)庫(kù)由聯(lián)合國(guó)糧農(nóng)組織(FAO)、中國(guó)科學(xué)院南京土壤研究所(ISSCAS)、歐盟委員會(huì)聯(lián)合研究中心(JRC)等國(guó)內(nèi)外多家機(jī)構(gòu)聯(lián)合構(gòu)建。研究區(qū)域土地利用情況與土壤類(lèi)型見(jiàn)圖3。

圖3 研究區(qū)域土壤類(lèi)型及土地利用情況

2.2.3 水文氣象數(shù)據(jù)

降水?dāng)?shù)據(jù)來(lái)自三道河子、寬城、下河南等10個(gè)雨量站1998-2018年逐日數(shù)據(jù)，氣象數(shù)據(jù)選用多倫、圍場(chǎng)、承德及青龍4個(gè)氣象站同時(shí)段逐日數(shù)據(jù)，數(shù)據(jù)來(lái)源為中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn)。

2.2.4 污染物數(shù)據(jù)概化

農(nóng)村生活污染與畜禽污染年內(nèi)污染負(fù)荷產(chǎn)生量較為固定，因此在SWAT模型中將兩種污染分別根據(jù)其各自TN、TP、NH3-N的比例，定義新的化肥類(lèi)型，以連續(xù)施肥的形式進(jìn)行概化。施肥持續(xù)時(shí)間設(shè)置為1年，施肥頻次設(shè)置為每日一次，同時(shí)按照各子流域人口與面積加權(quán)，確定子流域污染排放量。

對(duì)于化肥污染，承德市主要種植作物為玉米、大豆、馬鈴薯，生長(zhǎng)季節(jié)約在每年3-10月份，因此設(shè)置施肥時(shí)間為每年3-9月，施肥頻次設(shè)置為每月一次。

點(diǎn)源污染數(shù)據(jù)采用各排污口排放數(shù)據(jù)輸入模型，由于收集到資料為年排放總量，因此將其概化為每月輸出量相等，即年內(nèi)12個(gè)月均勻排放至流域。

2.3 SWAT模型適用性評(píng)價(jià)

采用相關(guān)系數(shù)R2及納什效率系數(shù)(NSE)兩種指標(biāo)來(lái)評(píng)價(jià)模擬結(jié)果，其計(jì)算公式分別見(jiàn)式(4)、(5)。

(4)

(5)

2.3.1 徑流量率定驗(yàn)證結(jié)果分析

采用SWAT-CUP對(duì)模型參數(shù)進(jìn)行率定，經(jīng)過(guò)2 000次迭代得到最終模擬結(jié)果，計(jì)算各水文站率定期與驗(yàn)證期的相關(guān)系數(shù)與納什效率系數(shù)，見(jiàn)表3，潘家口水庫(kù)入庫(kù)徑流量模擬結(jié)果見(jiàn)圖4?？梢?jiàn)流域內(nèi)各控制站R2、NSE均高于0.50，說(shuō)明模型能較好地反映流域內(nèi)水文變化過(guò)程。尤其是在潘家口水庫(kù)入庫(kù)徑流模擬時(shí)，率定期和驗(yàn)證期NSE分別為0.80和0.79，表明構(gòu)建的SWAT模型對(duì)潘家口水庫(kù)水量模擬結(jié)果很好。

表3 各水文站模擬效果評(píng)價(jià)

圖4 潘家口水庫(kù)徑流模擬與實(shí)測(cè)對(duì)比

2.3.2 水質(zhì)率定驗(yàn)證結(jié)果

選擇2018年對(duì)模型進(jìn)行率定，2017年對(duì)模型進(jìn)行驗(yàn)證。對(duì)各站TN、TP模擬值與實(shí)測(cè)值的相關(guān)系數(shù)進(jìn)行計(jì)算，結(jié)果如表4、表5所示，R2均在0.57以上。從圖5可以看出，率定后的模型能較好地反映研究流域營(yíng)養(yǎng)物遷移轉(zhuǎn)化情況，但對(duì)于汛期的極值點(diǎn)模擬效果并不理想，主要原因?yàn)檠雌诘棕?fù)荷主要受泥沙影響，而本研究中由于缺乏泥沙監(jiān)測(cè)數(shù)據(jù)，未對(duì)其進(jìn)行率定驗(yàn)證，因而造成了誤差。

表4 各監(jiān)測(cè)斷面TN模擬效果評(píng)價(jià)

表5 各監(jiān)測(cè)斷面TP模擬效果評(píng)價(jià)

圖5 各控制站TN、TP模擬值與實(shí)測(cè)值對(duì)比

3 結(jié)語(yǔ)

基于DEM、土地利用、水文氣象與污染物等數(shù)據(jù)構(gòu)建了潘家口水庫(kù)控制流域的SWAT水量水質(zhì)模型，并針對(duì)水量、水質(zhì)等指標(biāo)對(duì)模型進(jìn)行了率定與驗(yàn)證，采用R2與NSE評(píng)價(jià)了模型在該地區(qū)的適用性，主要結(jié)論如下：

(1)基于水文氣象與水質(zhì)資料，選取潘家口水庫(kù)控制流域，劃分了34個(gè)子流域，1 473個(gè)水文響應(yīng)單元，構(gòu)建了水量與污染物負(fù)荷的SWAT模型，

(2)采用2000-2010年實(shí)測(cè)徑流數(shù)據(jù)對(duì)水量模型進(jìn)行率定，2011-2015年數(shù)據(jù)進(jìn)行驗(yàn)證，各水文站徑流模擬的相關(guān)系數(shù)與NSE均較高，說(shuō)明模型能較好地反映流域內(nèi)的水量變化過(guò)程。

(3)選取總氮與總磷兩個(gè)水質(zhì)指標(biāo)，以2017、2018年實(shí)測(cè)斷面數(shù)據(jù)進(jìn)行了污染物負(fù)荷模型的率定和驗(yàn)證，結(jié)果表明模型能較好地反映污染物的遷移轉(zhuǎn)化過(guò)程。