基于WRF模型的澎溪河流域非點(diǎn)源污染預(yù)測(cè)

張萬順，卜思凡，彭 虹，夏 函，朱 磊，劉 馨

(1.武漢大學(xué)資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430072； 2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072； 3.武漢大學(xué)中國發(fā)展戰(zhàn)略與規(guī)劃研究院，湖北 武漢 430072；4.武漢大學(xué)水利水電學(xué)院，湖北 武漢 430072)

流域非點(diǎn)源污染過程受自然和人類活動(dòng)等多重影響[1-5]，非點(diǎn)源污染精細(xì)化模擬預(yù)測(cè)研究對(duì)保障流域水環(huán)境安全和提高生態(tài)環(huán)境監(jiān)管水平具有重要意義，是流域水環(huán)境管理領(lǐng)域的研究熱點(diǎn)[6-7]。降水是影響流域非點(diǎn)源污染過程的關(guān)鍵因素[8-11]，精準(zhǔn)反映局部氣象過程的空間異質(zhì)性成為制約非點(diǎn)源污染模擬預(yù)測(cè)精度的瓶頸。傳統(tǒng)的研究多將氣象站點(diǎn)的降水?dāng)?shù)據(jù)進(jìn)行空間插值，導(dǎo)致降水等氣象條件局部差異的均化，造成在反映區(qū)域氣象過程的空間異質(zhì)性方面存在一定誤差[12-15]。綜合考慮大氣和陸面過程的氣象研究與預(yù)報(bào)(weather research and forecasting, WRF)模型具有分辨率高、參數(shù)化方案豐富、時(shí)空連續(xù)性強(qiáng)等特點(diǎn)，在流域洪水預(yù)報(bào)、徑流模擬、水質(zhì)預(yù)報(bào)等領(lǐng)域得到廣泛應(yīng)用[16-18]。WRF模型是基于非靜力平衡的數(shù)值模型，可作為全球模式和區(qū)域模式進(jìn)行天氣現(xiàn)象的數(shù)值模擬。WRF模型可有效地將降水?dāng)?shù)據(jù)降尺度至精細(xì)化格點(diǎn)單元[19-20]，作為水文模型的邊界條件，精準(zhǔn)模擬流域非點(diǎn)源過程，反映局部空間單元上非點(diǎn)源污染的異質(zhì)性。Lee等[21]通過耦合WRF模型和分布式水文模型(soil and water assessment tool, SWAT)研究不同流量條件下氣象數(shù)據(jù)分辨率對(duì)流域徑流模擬的影響，顯著提高了模擬精度。Yuan等[22]構(gòu)建了集成空氣質(zhì)量(community multiscale air quality, CMAQ)模型、WRF模型和SWAT模型的流域水質(zhì)綜合評(píng)估方法，對(duì)大型流域氮的遷移轉(zhuǎn)化評(píng)估表現(xiàn)出較好的模擬效果。

三峽庫區(qū)流域特殊的降水、地形地貌和土地利用的空間異質(zhì)性和復(fù)雜性加劇了流域水文和水質(zhì)過程的不確定性，流域非點(diǎn)源污染問題突出。本文以三峽庫區(qū)典型支流澎溪河為研究區(qū)域，基于WRF模型，構(gòu)建了能反映流域3 km×3 km精細(xì)格點(diǎn)降水的流域非點(diǎn)源污染模擬預(yù)測(cè)模型，能精準(zhǔn)模擬流域非點(diǎn)源過程，準(zhǔn)確反映局部空間單元非點(diǎn)源污染的異質(zhì)性，為水資源科學(xué)管控和水環(huán)境精細(xì)化治理提供參考。

1 研究區(qū)概況

澎溪河位于北緯30°50′～31°42′、東經(jīng)107°56′～108°54′，是三峽庫區(qū)上游北岸的一級(jí)支流，干流全長(zhǎng)約182 km，面積5 172.5 km2，流域內(nèi)有東河、桃溪河、南河和普里河等支流，如圖1所示。流域海拔高度為148～2 549 m，地勢(shì)北高南低。流域處于亞熱帶季風(fēng)區(qū)，降雨相對(duì)充沛，氣候潮濕，多年平均降水量和溫度分別為800～1 500 mm和18.5 ℃。根據(jù)支流入?yún)R情況將流域劃分為東河、桃溪河-南河、普里河-澎溪河干流3個(gè)區(qū)域。

圖1 澎溪河流域

2 模型構(gòu)建

2.1 WRF模型

WRF模型是完全可壓縮的非靜力中尺度模型，模型水平方向采用Arakawa-C網(wǎng)格，垂直方向采用地形跟隨坐標(biāo)系；采用Runge-Kutta時(shí)間積分方案求解非靜力歐拉方程，在中尺度天氣模擬預(yù)報(bào)方面表現(xiàn)出突出的優(yōu)勢(shì)[23]。本文設(shè)置雙層嵌套區(qū)域，格點(diǎn)水平分辨率分別為9 km和3 km，WRF格點(diǎn)設(shè)置和澎溪河流域6個(gè)氣象站點(diǎn)分布如圖2所示，其中，澎溪河流域內(nèi)共設(shè)置568個(gè)3 km×3 km格點(diǎn)。將ERA-interim再分析數(shù)據(jù)集作為WRF模型的初始場(chǎng)和邊界條件，采用國家氣象觀測(cè)站數(shù)據(jù)進(jìn)行數(shù)據(jù)同化校正，模型運(yùn)行時(shí)間為2009—2012年，輸出時(shí)間間隔為24 h。

圖2 WRF模型格點(diǎn)

2.2 非點(diǎn)源模型

SWAT模型是基于過程的分布式水文模型，通過將流域劃分為子流域和水文響應(yīng)單元(hydrological response unit, HRU)模擬預(yù)測(cè)復(fù)雜流域的水文過程和產(chǎn)流產(chǎn)污規(guī)律[7]。模型以HRU為基本計(jì)算單元，再匯集到子流域中，將多個(gè)子流域匯集成整個(gè)流域進(jìn)行計(jì)算，模擬流域徑流過程和污染物遷移轉(zhuǎn)化規(guī)律。

本文分別將氣象站降水和WRF模型輸出降水作為SWAT模型的輸入，采用單向耦合方式將WRF模型內(nèi)層嵌套的3 km×3 km精細(xì)格點(diǎn)上模擬的降水過程通過數(shù)據(jù)處理和尺度轉(zhuǎn)換輸入SWAT模型，進(jìn)行流域徑流和非點(diǎn)源污染模擬。輸入SWAT模型的污染源包括點(diǎn)源和附加非點(diǎn)源，其中點(diǎn)源主要包括工業(yè)點(diǎn)源和城鎮(zhèn)污水處理廠，非點(diǎn)源主要包括農(nóng)村生活污水、畜禽養(yǎng)殖和農(nóng)業(yè)化肥等。

2.3 模型率定與驗(yàn)證

本文采用Pearson相關(guān)系數(shù)(r)和均方根誤差(RMSE)2個(gè)指標(biāo)對(duì)WRF模型降水模擬效果進(jìn)行評(píng)價(jià)[23]，采用決定系數(shù)(R2)評(píng)價(jià)SWAT模型對(duì)徑流量、總磷(TP)和總氮(TN)的模擬效果[8]。

選取澎溪河流域周圍6個(gè)氣象站點(diǎn)2009—2012年逐日連續(xù)降水觀測(cè)數(shù)據(jù)對(duì)WRF模型降水模擬效果進(jìn)行評(píng)價(jià)。如表1所示，日尺度降水模擬的r為0.27～0.44，RMSE為7.43～12.83 mm/d；月尺度降水模擬的r為0.58～0.89，RMSE為65.76～175.39 mm/月。模擬效果與現(xiàn)有研究結(jié)果接近[23-24]，模擬效果較為可靠，可用于研究降水對(duì)流域徑流和污染負(fù)荷時(shí)空分布的影響。

表1 WRF模型模擬降水驗(yàn)證結(jié)果

采用澎溪河流域溫泉水文站的逐日徑流監(jiān)測(cè)數(shù)據(jù)對(duì)SWAT模型徑流過程模擬效果進(jìn)行評(píng)價(jià)，以2011年為模型率定期，2012年為驗(yàn)證期，模擬效果如表2和圖3所示。采用WRF模型輸出降水的徑流模擬精度均高于采用氣象站降水，驗(yàn)證期模型徑流模擬精度提高了27%，徑流模擬結(jié)果滿足模型模擬精度要求。因此，基于WRF模型輸出降水構(gòu)建的SWAT模型可用于模擬澎溪河流域徑流過程。

(a) WRF降水模擬徑流

表2 SWAT模型徑流和水質(zhì)率定驗(yàn)證結(jié)果

采用澎溪河流域水東壩監(jiān)測(cè)斷面逐月水質(zhì)監(jiān)測(cè)數(shù)據(jù)對(duì)SWAT模型水質(zhì)模擬效果進(jìn)行評(píng)價(jià)，以2011年為模型率定期，2012年為驗(yàn)證期，模擬效果如表2和圖4所示。采用WRF模型輸出降水的水質(zhì)模擬精度均高于采用氣象站降水，驗(yàn)證期模型TP質(zhì)量濃度模擬精度提高了31%，TN質(zhì)量濃度模擬精度提高了36%，水質(zhì)模擬結(jié)果滿足模型模擬精度要求。因此，基于WRF模型輸出降水構(gòu)建的SWAT模型可用于澎溪河流域非點(diǎn)源污染研究。

(a) WRF降水模擬TP質(zhì)量濃度

3 結(jié)果與分析

3.1 降水量時(shí)空分布

WRF模型模擬降水和氣象站降水過程如圖5和圖6所示。WRF模型模擬降水量與氣象站降水量峰值對(duì)應(yīng)較好，流域內(nèi)60%～80%的降水量集中在6—9月，降水量峰值多出現(xiàn)在7—8月，降水量較少的月份出現(xiàn)在11月至次年3月，呈現(xiàn)出明顯的豐枯月差異。WRF模型模擬的日降水量最大值為 141.4 mm，氣象站日降水量最大值為122.6 mm；WRF模擬的月均降水量約為120 mm，月降水量最大值為585.7 mm；氣象站月均降水量約為95 mm，月降水量最大值為340.9 mm。

(a) 奉節(jié)站

(a) 奉節(jié)站

WRF模型模擬的2011年3月、6月、9月和12月降水量空間分布如圖7所示，降水量呈精細(xì)格點(diǎn)狀分布，空間異質(zhì)性明顯，4個(gè)典型月降水量均呈現(xiàn)出從東北向西南遞減趨勢(shì)，降水量較多的區(qū)域主要分布在流域東北部和南部邊緣等區(qū)域，這些區(qū)域多為海拔較高的林地和草地，對(duì)水汽具有阻擋抬升作用，容易形成降水[25]。

(a) 3月

從各格點(diǎn)單元月降水量看，澎溪河流域3月、6月、9月和12月格點(diǎn)單元降水量最大值分別為 209 mm、699 mm、553 mm和82 mm。從區(qū)域分布上，3月最大降水量從大到小依次分布在東河、桃溪河-南河、普里河-澎溪河干流；6月和9月最大降水量從大到小依次分布在東河、普里河-澎溪河干流、桃溪河-南河；12月最大降水量從大到小依次分布在普里河-澎溪河干流、東河、桃溪河-南河。

從各區(qū)域月平均降水量來看，月平均降水量最大的區(qū)域均集中在東河，3月、6月、9月和12月東河月平均降水量分別為92 mm、268 mm、272 mm和43 mm，各區(qū)域月平均降水量從大到小依次為：東河、桃溪河-南河、普里河-澎溪河干流。

3.2 徑流量時(shí)空分布

采用WRF模型輸出的降水?dāng)?shù)據(jù)作為非點(diǎn)源模型輸入條件，模擬的澎溪河流域2010—2012年年徑流量分別為54.04億m3、42.54億m3和35.28億m3，呈下降趨勢(shì)。澎溪河流域多年平均年徑流量約為35.80億m3[26]，基于WRF模型降水?dāng)?shù)據(jù)模擬的年徑流量范圍與Shi等[7]研究中模擬的年徑流量范圍31.90億～47.30億m3接近。

基于WRF模型輸出降水模擬的澎溪河流域地表徑流深空間分布如圖8所示，整體上呈現(xiàn)東北部大、東南部次之、中部最小的規(guī)律。2010—2012年，澎溪河流域徑流深總體呈現(xiàn)出逐年遞減的趨勢(shì)，年平均徑流深依次為1 617 mm、1 270 mm和1 021 mm，平均徑流深最大的區(qū)域均集中在東河，各區(qū)域平均徑流深由大到小依次為：東河、桃溪河-南河、普里河-澎溪河干流。徑流深大于1 200 mm的子流域單元集中分布在東河?xùn)|北部，徑流深小于500 mm的子流域單元多集中在普里河北岸以及澎溪河干流附近。流域東北部地勢(shì)高，對(duì)水汽具有阻擋抬升作用，易形成降水，增加了流域產(chǎn)流量；流域南部植被覆蓋度較高，水源涵養(yǎng)能力較強(qiáng)，產(chǎn)流相對(duì)其他區(qū)域較小。

(a) 2010年

3.3 非點(diǎn)源負(fù)荷時(shí)空分布

采用WRF模型輸出的降水?dāng)?shù)據(jù)作為非點(diǎn)源模型輸入條件，模擬的澎溪河流域2010—2012年TP負(fù)荷分別為1 402.50 t/a、1 140.06 t/a和911.50 t/a，均值為1 151.35 t/a；TN負(fù)荷量分別為14 503.39 t/a、11 487.39 t/a和9 288.37 t/a，均值為11 759.72 t/a，總體呈下降趨勢(shì)，與徑流量變化趨勢(shì)吻合度較高。高銀超等[26]研究中，澎溪河流域TN負(fù)荷年均值為8 335.23 t/a，石熒原等[7]研究中澎溪河流域TP負(fù)荷范圍為802～1 432 t/a，TN負(fù)荷范圍為7 674～13 775 t/a，與本文研究結(jié)果接近。

基于WRF模型輸出降水模擬的澎溪河流域氮磷負(fù)荷空間分布如圖9和圖10所示。TN、TP負(fù)荷空間分布總體呈現(xiàn)出“局部集中、靠近水體”的特點(diǎn)[27]。2010年，澎溪河流域各區(qū)域最大單位面積TP負(fù)荷為8.00 kg/hm2，各區(qū)域單位面積TP負(fù)荷由大到小依次為：東河、桃溪河-南河、普里河-澎溪河干流；2011年，各區(qū)域最大單位面積TP負(fù)荷為11.99 kg/hm2，各區(qū)域單位面積TP負(fù)荷由大到小依次為：桃溪河-南河、普里河-澎溪河干流、東河；2012年，各區(qū)域最大單位面積TP負(fù)荷為9.27 kg/hm2，各區(qū)域單位面積TP負(fù)荷由大到小依次為：桃溪河-南河、普里河-澎溪河干流、東河。2010—2012年，澎溪河流域單位面積TN負(fù)荷最大的區(qū)域均為東河，單位面積TN負(fù)荷量依次為99.58 kg/hm2、81.37 kg/hm2和62.23 kg/hm2，各區(qū)域單位面積TN負(fù)荷由大到小依次為：東河、桃溪河-南河、普里河-澎溪河干流。單位面積TP負(fù)荷高于10 kg/hm2的子流域單元主要集中于東河?xùn)|北部、桃溪河西岸、南河上游兩岸和普里河南岸；單位面積TN負(fù)荷高于80 kg/hm2的子流域單元主要集中于東河?xùn)|北部、普里河南岸。流域東北部區(qū)域海拔較高、坡度較大，降雨沖刷作用強(qiáng)，易造成水土和養(yǎng)分流失；流域南部耕地和草地面積比例較大，農(nóng)藥化肥使用量大，導(dǎo)致該區(qū)域非點(diǎn)源污染負(fù)荷增大[28]。

(a) 2010年

(a) 2010年

4 結(jié) 論

a.本文構(gòu)建的WRF模型模擬的降水呈精細(xì)格點(diǎn)狀分布，對(duì)日尺度降水模擬的相關(guān)系數(shù)最大可達(dá)到0.44，RMSE為7.43～12.83 mm/d；月尺度降水模擬的相關(guān)系數(shù)最大可達(dá)到0.89，RMSE為65.76～175.39 mm/月，模擬效果較好。

b.基于WRF模型，構(gòu)建了能反映流域3 km×3 km 精細(xì)格點(diǎn)降水的流域非點(diǎn)源污染預(yù)測(cè)模型。與由澎溪河流域6個(gè)氣象站點(diǎn)2009—2012年觀測(cè)降水資料得到的模擬結(jié)果相比，基于WRF模型輸出降水的徑流量模擬精度提高了27%，TP和TN質(zhì)量濃度模擬精度分別提高了31%和36%。

c.采用WRF模型輸出降水模擬的澎溪河流域年徑流量范圍為35.28億～54.04億m3，徑流深大于1 200 mm的子流域單元分布在東河?xùn)|北部，年徑流量表現(xiàn)出較明顯的空間異質(zhì)性。

d.基于WRF模型輸出降水模擬的TP負(fù)荷范圍為911.50～1 402.50 t/a，TN負(fù)荷范圍為9 288.37～14 503.39 t/a。單位面積TP負(fù)荷高于10 kg/hm2的子流域單元主要集中于東河?xùn)|北部、桃溪河西岸、南河上游兩岸和普里河南岸；單位面積TN負(fù)荷高于 80 kg/hm2的子流域單元多集中在東河?xùn)|北部、普里河南岸，空間分布差異顯著。