流域水文模型在面源污染模擬與管控中的應用研究進展

謝 暉,邱嘉麗,董建瑋,高田田, 4,賴錫軍

1 中國科學院南京地理與湖泊研究所,中國科學院流域地理學重點實驗室,南京 210008 2 中國科學院地理科學與資源研究所,中國科學院陸地水循環及地表過程重點實驗室,北京 100101 3 南京師范大學海洋科學與工程學院,南京 210023 4 南京信息工程大學應用氣象學院,南京 210044

面源,或稱之為非點源,是以面狀分布并排放污染物進而造成水體污染的發生源。隨著點源污染控制水平的不斷提高,面源污染成為全球水體污染的重要污染來源[1—2]。按照污染類型,面源污染主要分為城市面源和農業面源污染。城市面源污染伴隨著城市化的快速發展,不透水下墊面比例增加,暴雨徑流及攜帶的面源污染物超出城市排水系統的消納和處理能力,對城市水環境帶來沖擊。近年來雨污分流改造和海綿城市建設不斷推進,城市面源污染逐漸得到緩減和控制。我國作為粗放型農業生產大國,對化肥農藥的依賴度極高,不合理的水肥管理導致的農業面源污染問題尤為突出,威脅著水環境健康、飲用水和糧食安全[3]。第二次全國污染源普查數據顯示,農業源排放的化學需氧量、總氮和總磷分別占水污染物排放總量的49.8%、46.5%和67.2%,表明農業面源污染仍然是造成我國水污染的重要原因。我國高度重視農業面源污染防治工作,自2015年以來,國家相關部委相繼出臺了《關于打好農業面源污染防治攻堅戰的實施意見》、《重點流域農業面源污染綜合治理示范工程建設規劃(2016—2020年)》等一系列政策性措施,習近平總書記也強調要“以釘釘子精神推進農業面源污染防治”。2021年3月生態環境部和農業農村部聯合印發了《農業面源污染治理與監督指導實施方案(試行)》,明確了“十四五”至2035年重點區域農業面源污染防治的工作目標和主要任務,對面源污染的監測和評估提出了更高的要求。

面源污染的發生和遷移包含若干動態過程,如降雨徑流、土壤侵蝕、污染物累積與沖刷、土壤溶質淋溶等,具有隨機性、復雜性、不確定性及時滯性等特點,獲取一定精度和密度的監測數據需要高強度且復雜的野外監測,現有監測網絡在點位分布和技術手段等方面還不足以準確評估面源污染對水體污染的貢獻,分布式和半分布式流域水文模型逐漸盛行。目前國內外應用于面源污染模擬的此類模型主要有SWAT模型(Soil and Water Assessment Tool)[4]、SWMM模型(Storm Water Management Model)[5]、AnnAGNPS模型(Annualized Agricultural Non-Point Source Pollution Model)[6]、GWLF模型(Generalized Watershed Loading Function)[7]、SPARROW模型(SPAtially Referenced Regressions On Watershed attributes)[8]、MIKE SHE模型[9]和HSPF模型(Hydrological Simulation Program-FORTRAN)[10—11]等。這些模型既有相通之處,又具備各自的獨特性,用戶可以根據研究區特色、模擬尺度、數據完備性、研究目的等多方面進行模型的適用性選擇。本文針對HSPF模型進行綜述。HSPF模型是美國環保署重點推薦的用于模擬流域水文及水質的綜合模型。它起源于20世紀60年代的斯坦福水文模型(Stanford Watershed Model),隨后水質模塊被逐漸引入,并最終由Robert Carl Johanson于1980年正式提出[11]。發展至今,HSPF模型的開發和維護主要在美國環保署的資助下進行,現已更新至12.5 版本,并且被集成至綜合模擬平臺BASINS(Better Assessment Science Integrating Point and Non-point sources)。該平臺集成多種模型、數據處理及分析工具,實現了HSPF模型的可視化配置、運行及結果分析。HSPF模型現階段的主要開發方RESPEC公司強調了該模型的兩個特點：“comprehensive(綜合)”和“flexible(靈活)”,其綜合性體現在模型在連續過程、動態事件以及穩態條件下均可模擬流域水文、河道水力以及水質的復合過程,模擬范圍涵蓋了面源污染的產生排放、遷移轉化及水體輸運等主要過程；靈活性體現在模型可根據實際情況調整模擬步長、計算單元的劃分方式、陸域與水體的承接關系、污染物遷移的模擬方法等,以應對面源污染在不同尺度需求下的模擬。這兩個重要的特點突出了HSPF模型在面源污染模擬方面的優勢,使其在全球范圍內得以廣泛運用。

近年來HSPF模型逐漸被用于我國流域水環境過程模擬,應用過程中也面臨著模型機理和參數本地化、模型構建精細化、模型結果不確定性較大、落地應用困難等方面的挑戰。本文圍繞HSPF模型在面源污染模擬與管控中的研究進展進行了綜述。首先簡介該模型的結構功能及基本原理,綜述其在變化環境下的模擬方法和成果;其次,對應對參數識別、不確定性分析、措施效果評估和總量控制的思路和方法進行了總結;最后,分析了新形勢下該模型的延伸拓展,提出了未來HSPF模型可進一步深化的發展動向。

1 HSPF模型結構及功能、基本原理

1.1 結構及功能

HSPF模型的運行主要依賴兩個文件,即用戶控制輸入(User Control Input,UCI)文件和流域數據管理(Watershed Data Management,WDM)文件。前者是ASCII文本文件,以模塊化形式記載流域屬性和模型配置。WDM文件是一種非格式化、32位架構的二進制文件,它存儲HSPF模型所必須的氣象驅動數據及模擬結果產生的時間序列數據。在最新版的BASINS 4.5下,這兩個文件的生成和編輯可通過SARA Time series Utility工具和winHSPF軟件實現,HSPF模型的搭建環境及流程如圖1所示。

圖1 基于BASINS 4.5的HSPF模型搭建環境及流程Fig.1 Environment and process to develop HSPF model based on BASINS 4.5BASINS: 點源和面源污染評估系統 Better Assessment Science Integrating Point and Non-point Sources；WDM: 流域數據管理 Watershed Data Management；HSPF: 流域水文模型HSPF Hydrological Simulation Program-FORTRAN；UCI: 用戶控制輸入 User Control Input；DEM: 數字高程模型 Digital Elevation Model

對于一個目標流域,HSPF模型會將其地理空間概化為3類模擬單元：透水陸面(PERLND)、非透水陸面(IMPLND)和水體(RCHRES),每個子流域中則是由不同數量和種類的3類模擬單元組成,子流域間再由拓撲關系相聯。HSPF模型對模擬單元劃分時可考慮將氣象驅動、土地利用和土壤類型等環境因子作為空間異質性的劃分準則,并將不同因子映射于每個模擬單元,形成水文響應單元。這種劃分方法體現了該模型半分布式的架構,有利于表征面源污染輸出的時空分異特征。

HSPF模型功能特點可很好地應對面源污染的模擬。從研究區特征來說,由于HSPF模型對陸面進行了透水和非透水劃分,因此可對農業流域、城市流域或混合復雜的用地特征流域進行面源污染模擬。當然,模型對非透水陸面的模擬過程中考慮的物理過程相比于專門針對城市面源污染模擬的SWMM模型來說較為簡單,不能體現管網匯流這一城市下墊面特色,但HSPF模型也為模型使用者提供了一種簡易的城市面源污染模擬方法。從模擬的時間尺度上看,由于HSPF模型受連續氣象數據的驅動,因此可模擬面源污染的長時間連續過程。單次降雨事件帶來的面源污染輸出同樣受到業界的關注,由于HSPF模型模擬時間步長可在小時至分鐘調整,輸出時間序列可細至每小時,因此也可用于模擬單次降雨事件下的面源污染完整過程。從模擬對象上看,HSPF模型成熟的水文模擬框架能以較高的精度模擬徑流和演算河道洪水。懸浮泥沙在陸面的分離搬遷以及在河道中的動力輸運過程也在HSPF模型中得以刻畫。常見的面源污染物如氮磷的累積沖刷、形態間轉化、各相間分配、作物和藻類吸收等基本過程均在模型中有所考慮,因此HSPF模型對面源氮磷污染模擬具有很好的適用性。此外,HSPF模型還有專門模擬農藥流失和擴散的PEST模塊。對于其他化學物質,HSPF模型依據物質的特性(如保守性、衰減性、吸附性等),可在模塊中新增污染物并設定相應的參數進行模擬,已有研究利用HSPF模型模擬過大腸桿菌、汞、有機碳的面源污染過程并取得了較好的效果[12—14]。

1.2 基本原理

徑流：徑流由PERLND模塊和IMPLND模塊分別計算,前者是HSPF模型中的核心模塊,涉及到的物理過程最多。透水陸面的剖面自上而下分層為地表蓄積、壤中流蓄積、上土壤層蓄積、下土壤層蓄積和地下水蓄積。降水在經冠層截留后落在地表。HSPF模型沿用斯坦福模型的概化方法,將落地雨的潛在去向劃分為：形成地表滯蓄增量、形成壤中流滯蓄增量或直接下滲,且入滲量隨地塊面積呈線性變化。具體的下滲容量分配概化如圖2所示。

圖2 HSPF模型下滲容量分配(改自Bicknell等[10])Fig.2 Infiltration distribution concept in HSPF model (modified from Bicknell et al. [10])

落地雨經以上分配后,各土壤蓄積進行動態的水量平衡計算。其中從地表蓄積中產生的地表徑流以Chezy-Manning方程計算；壤中流按線性水庫法計算；補給河道流量的地下水部分計為有效地下水,由地下水蓄積量、地下水蓄泄系數、地下水退水系數、地下水水力坡度構建的非線性方程計算。非透水陸面不考慮土壤分層,產流只包含地表徑流,采用的算法與透水陸面中一致。河道和湖庫中的徑流演算在RCHRES模塊中進行,但認為水體具有單向流動性,單層且完全混合。在以上的假定條件下,RCHRES采用水力函數表FTABLE結合插值法模擬水深、容積和流量的關系。

泥沙：HSPF模型考慮了粒徑的影響,將泥沙細化為沙粒、粉粒和粘粒三種組分。在透水陸面上對泥沙輸移量的模擬過程包括泥沙與土壤基質的吸附、分離,以及分離泥沙和土壤基質的沖刷,沖刷量的計算依賴于地表徑流量、沖蝕系數和指數。非透水陸面泥沙輸移量的模擬僅考慮堆積和沖刷過程。河道中泥沙首先被區分為粘性泥沙(粉粒和粘粒)和非粘性泥沙(沙粒)。對于粘性泥沙,用戶須設置懸浮泥沙的沉降臨界剪切力以及床面泥沙的起動臨界剪切力。當水流剪切力介于沉降臨界剪切力與起動臨界剪切力之間時,陸面進入水體的粘性泥沙以及河道中的懸浮粘性泥沙將通過平流過程輸運；當水流剪切力大于起動臨界剪切力時,河床上的粘性泥沙將被沖起并輸運；當水流剪切力小于沉降臨界剪切力時,粘性泥沙將沉降至河床。對于非粘性泥沙,模型將首先計算水流對非粘性泥沙的輸運能力。當該輸運能力超出陸面進入水體的與河道中的懸浮非粘性泥沙之和時,河床上的非粘性泥沙將被沖起并輸運；反之,非粘性泥沙將沉降至河床。

污染物：污染物在陸面的遷移過程可通過兩種方式進行模擬,一種是基于PQUAL和IQUAL子模塊分別在透水單元和非透水單元中進行模擬,模擬方法主要基于經驗方程。吸附態污染物的輸移通量計算由泥沙沖刷量與沖刷系數的乘積獲得,沖刷能力系數是需要率定的關鍵參數。溶解態污染物的傳輸路徑有地表徑流、壤中流和基流,地表徑流中的污染物輸移通量由地表污染物蓄積量、地表徑流量以及地表徑流對污染物的沖刷速率共同決定。壤中流和基流攜帶的污染物質量計算受壤中流和基流中污染物濃度影響,且污染物濃度值可設置為每月變化。另一種是通過NITR、PHOS和PEST子模塊,以質量守恒為基礎分別考慮氮、磷和農藥在土壤蓄積中的轉化、淋溶以及遷移等復雜的遷移轉化過程。前者允許用戶設置多達10種污染物,糞大腸菌群主要通過這種方式模擬。降水發生之前,糞大腸菌群在地表的累積量取決于細菌的積累速率,而最大累積量取決于土地利用類型；隨著降水的發生,地表徑流的沖刷效應使累積的細菌進入河道,糞大腸菌群負荷與徑流量密切相關,而壤中流或地下水的貢獻主要根據這些傳輸介質中菌群的恒定濃度計算。后者是更為復雜的模擬方式,僅能夠模擬透水陸面中的氮、磷和農藥三種污染物的遷移轉化過程。這三種污染物對透水陸面的來源輸入主要來自大氣干濕沉降和施肥。NITR子模塊中考慮的氮的形態包括氨態氮、硝態氮、活性有機氮和非活性有機氮,模擬主要過程包括氮與泥沙間的吸附與解吸、硝化與反硝化、礦化與固定、固氮與揮發、作物吸收與氮回收等。PHOS子模塊用于模擬磷酸鹽和有機磷在土壤蓄積中儲存量的變化,主要模擬過程有磷酸鹽與泥沙的吸附與解吸、礦化與固化、作物吸收。PEST子模塊用于模擬溶解態、吸附態和晶體態農藥在土壤中的吸附、解吸以及降解過程。以上過程采用一階反應動力學方程、FREUNDLICH方程和溫度修正的ARRHENIUS方程等進行表達。

水體的污染物輸移過程也可以采用簡單或復雜兩種模擬方式。簡單的模擬方式主要基于通用的GQUAL子模塊進行模擬,其中不涉及生物化學過程、溶解態污染物考慮其降解過程、吸附態污染物考慮隨泥沙的沉降、懸浮、降解以及吸附和解吸過程等。復雜的模擬方式是為了滿足需要考慮氮磷在水體中的生物化學轉化作用的用戶需求,主要基于RQUAL子模塊,該子模塊首先模擬溶解氧和生化需氧量的平衡過程,進而模擬無機氮磷的平衡過程,以及氮磷與浮游生物相關的過程,最終計算污染物經隨流輸移從RCHRES的出口匯出量。

2 HSPF模型在變化環境下的面源污染模擬

2.1 氣候變化

氣候變化在改變流域水文循環過程的同時,也影響著流域污染物的物理化學過程、遷移轉化過程及水體對污染物的稀釋能力。氣候因子,如降雨、溫度是流域面源污染的主要驅動力,因此預測未來氣候情景下流域面源污染的發展趨勢對未來流域水環境的管理和規劃具有重要意義。對氣候變化情景模擬的主要思路是將氣候預測系統同流域模型相耦合。一般包括以下步驟：流域模型率定與驗證；構建氣候變化情景；將氣候數據輸入到模型中,模擬氣候情景下的流域水文過程和面源污染輸移過程。HSPF模型在評估氣候變化對流域水文的影響方面具有很好的適用性。Gizaw等[15]采用HSPF模型評估RCP4.5(Representative Concentration Pathway 4.5)和RCP8.5(Representative Concentration Pathway 8.5)情景下的氣候變化對埃塞俄比亞四個主要流域徑流量的影響,結果表明氣溫和降水的升高會導致明顯的徑流增加。在面源污染方面,張鵬飛[16]通過對密云水庫流域內四個氣象站1961—2000年40年的氣象資料分析,采用統計分析和線性回歸的方法預測流域氣候變化趨勢,并分別從任意情景設置和歷史資料兩方面確定了 25種氣候情景和3種水文情景,之后將不同情景下的氣象數據輸入至建立的HSPF模型, 模擬密云水庫流域不同氣候情景下水文與面源污染的變化趨勢。Yi等[17]采用HSPF模型評估了RCP4.5和RCP8.5情景下未來2020s、2050s和2080s內單位面積面源總氮和總磷排放強度的變化。Kim等[18]通過全球氣候模式(GCMs)預測未來氣候情景,結合美國地質調查局的土地利用預測數據,采用HSPF模型分析了氣候變化和城市化對博伊西河流域水文和水環境的影響,結果表明氣候變化與城市化的綜合影響可導致面源污染輸出呈現顯著的季節變異性。

2007年,美國環保署在BASINS系統中集成了氣候評估工具CAT(Climate Assessment Tool),允許用戶采用CAT創建的氣候變化情景作為HSPF模型的輸入,為氣候變化對面源污染的影響研究以及水資源管理提供了便利。CAT的主要功能是依據任意基準期的溫度和降水變化規律預測未來氣候的可能變化情景,基準期氣候數據包括歷史實測數據與氣候模式輸出數據[19]。Zhou等[20]采用HSPF-CAT耦合模型評估了未來降水變化對尖峰嶺雨林流域水文過程的影響。Taner等[21]在四個GCM的降尺度輸出數據的基礎上,運用HSPF-CAT耦合模型評估了氣候變化對徑流量和面源污染的影響,結果表明極端降水事件發生頻率的增加提高了洪峰量和營養鹽流失量。總之,CAT為氣候變化情景的創建與修正提供了便利,通過BASINS平臺可較為方便地與HSPF模型耦合,使HSPF-CAT成為評價氣候變化對面源污染影響的有效工具。

氣候變化對面源污染的影響不僅表現在平均態氣候變化的影響上,更高強度的擾動事件對面源污染的產輸量有著更為嚴重的影響,這種極端氣候事件包括洪澇、干旱及極端氣溫等。目前運用HSPF模型評估極端氣候事件對面源污染的影響研究相對較少。Ouyang等[22—23]采用HSPF-CAT耦合模型評估了未來可能的氣溫升高強度(增加1°C和2°C)與極端降水事件(雨強增加10%和20%)對美國亞祖河流域泥沙、硝酸鹽和磷酸鹽負荷的協同影響。Qiu等[24]結合土壤含水量的變化特征,采用HSPF模型評價了密云水庫上游流域干旱與極端降水條件對面源污染的協同作用,結果表明,長期干旱后的極端降水事件可能導致累積的沉積物、營養物、糞大腸菌群與其他污染物大量流失,污染物濃度嚴重超出水質標準,對下游水體的水環境健康造成威脅。這項研究突出了HSPF模型在評估極端氣候對面源污染研究中的適用性,也指出了適應極端氣候的流域管理體系的重要性。

2.2 土地利用/覆被變化

土地利用/覆被變化能夠帶來水文條件、下墊面特征、人類活動強度的空間差異,進而影響污染負荷產生和傳輸的時空分異。HSPF 模型可被用于模擬歷史和未來土地利用/覆被變化情景的影響,也能用于模擬土地利用/覆被與氣候變化雙重驅動下的流域水文與面源污染輸出特征,并表現出較好的模擬精度[25—28]。Bello等[29]采用HSPF模型評估了植被覆蓋變化對穆爾河流域水質的影響,結果發現BOD、硝態氮和磷酸鹽對植被覆蓋變化非常敏感。Bai等[27]通過HSPF模型評估了東江流域土地利用變化下流域面源污染的時空變化特征,發現城市化趨勢增加了BOD與總磷的負荷量,城市生活污水與城市雨洪是這兩種污染物的主要來源；而流域內總氮主要來源于氮肥的施用,流域耕地面積減少使得氮肥施用量減少,進而導致總氮負荷減少。Liu和Tong[30]采用HSPF模型量化了不同寬度(60 m,90 m和120 m)河岸帶情景下徑流、硝酸鹽及總磷負荷的變化,研究結果指導了河岸帶森林和濕地緩沖區的規劃方案設計。Mohamad等[31]采用HSPF模型評估了農業用地擴張對面源污染的影響,結果表明農業用地增加60%使得年均總磷負荷和總氮負荷分別增加3.8%和5.3%,農業用地增加2倍使得總氮和總磷負荷量增加約2倍。在太湖西苕溪流域,Luo等[32]首先利用CA-Markov模型生成2020s和2030s的土地利用情景,同時結合GCM的降尺度未來氣候情景,協同驅動HSPF模型來模擬流域徑流和營養鹽輸出。

綜上,氣候變化與土地利用/覆被變化對流域水環境過程有著顯著的影響,HSPF模型不僅可有效地模擬單一變化情景下的水文與面源污染過程,還可模擬多種氣候變化與未來土地利用變化組合情景對面源污染的影響。

3 模型參數識別與不確定性分析

3.1 參數識別

面源污染過程涉及到水循環、土壤侵蝕、污染物遷移轉化等多種物理化學過程,HSPF模型在概化這些過程時所用的參數眾多,包括具有明確物理機制的確定型參數(如河道底泥氨氮在厭氧環境下的釋放速率等),也包括集總式的概念型參數(如可帶走土壤中90%污染物的地表徑流速率等),針對前者的取值可通過實測獲取,而后者的取值則需要進行參數識別。參數識別是面源污染模型建模過程中的關鍵步驟,通常包含參數敏感性分析及參數率定。

參數敏感性分析有助于識別對模擬結果具有顯著影響的模型參數,在后續參數率定中可重點針對敏感性參數進行調整,可以一定程度上減少模型校驗的工作量。此外,由于不同類型參數代表著流域系統內不同子系統過程,被識別的高敏感度參數也有助于甄別關鍵過程并指導高效削減的管控方案設計。常用于HSPF模型的參數敏感性分析方法包括Morris、LH-OAT、RSA、Sobol′s等[33—36]。在參數調整過程中,需要注意的是模型參數的敏感度可能在不同的水雨情條件和不同時間尺度下發生變化,這就需要考慮參數的動態識別分析[37]。Xie等[38]將HSPF模型、傅立葉振幅敏感性分析、時間變化與多尺度結合的方法結合應用,識別了水沙參數在不同水文情勢和不同時間尺度下的敏感度,對關鍵參數的監測獲取及污染削減措施實施提出了指導建議。

針對HSPF模型參數的手動率定,美國地質調查局最早研發了其水文參數的專家決策軟件HSPExp[39],之后又研發了該軟件的改進版HSPExp+[40],并增加了泥沙和營養鹽模塊的手動率定驗證功能。 Mohamad[41]利用HSPExp+對HSPF模型在馬來西亞Johor River上游流域的營養鹽模擬進行了率定驗證,并以此評估面源污染輸出對氣候變化和化肥施用情景的響應。現如今越來越多的自動率定技術被應用于HSPF模型的水文模擬研究。早期應用最為廣泛的是PEST程序[42],該程序是基于Gauss-Marquardt-Levenberg算法,解決的是模型參數取值尋優過程中的非線性最小二乘問題,由于PEST操作性強且迭代運算次數較少,現在依然在國內流域中使用較多[27,43—45]。遺傳算法則是近十年HSPF模型應用研究中被頻繁提及的自動率定算法。隨著多站點率定、多目標率定的需求,帶精英策略的非支配排序遺傳算法(Non-dominant sort genetic algorithm-II, NSGA-II)[46]被用于HSPF模型水文參數識別研究中。相較于傳統的遺傳算法,NSGA-II基于Pareto最優前沿的理論使用非支配占優的排序方法,可求解多目標函數的協同最優解。Xie等[36]針對徑流時間序列的特點定義了三種目標函數,利用NSGA-II對連續模擬和降雨場次模擬下的HSPF模型水文參數估值進行了比對。Xie和Lian[47]定義三個水文站流量模擬和觀測值的納什效率系數作為多目標函數,利用NSGA-II對HSPF模型水文參數的最優取值進行了求解。

現階段關于HSPF模型水質參數自動率定的研究非常有限,事實上隨著高頻水質自動監測站和嵌套式多站點水質監測的普及,越來越多的河流水質數據可服務于模型的水質參數自動率定。需要注意的是HSPF模型的水質參數中有很多是將多過程集總表達的概念型參數,并不具備共識的參數取值區間,如果局限于采用自動率定算法獲得這類水質參數的“最優值”,很有可能出現擬合效果較好但與實際情況偏差較大的情況,未來的研究應在自動率定的過程中充分融合模型使用者的專家經驗和先驗信息,以保證面源污染模擬結果的可靠性。

3.2 不確定性分析

隨著面源污染模型越來越多地服務于水資源管理、污染物總量控制、生態補償等領域,模型結果的可靠性問題和風險問題備受關注。由于流域內水循環及生物地球化學過程屬于降雨-徑流、面源污染物遷移轉化和河流水質水動力等高度非線性過程的復雜組合,而流域模型對這些機理過程通常概化成一般的數學方程,本質上是對真實過程的一種近似狀態的描述,這就會帶來模擬結果的不確定性。美國最大日負荷量(Total Maximum Daily Load,TMDL)計劃的具體實施中明確規定了安全余量的計算,要求如果基于模擬結果分配負荷則須要考量模擬結果不確定性,以降低流域管理的決策風險[48]。在自然科學領域,不確定性通常被歸為隨機不確定性和認知不確定性[49],而具體到面源污染模型的不確定性則概括為三類,即參數、輸入及結構不確定性[50]。

3.2.1參數不確定性

流域模擬中的參數不確定性最早在流域水文模型中發現和探討。1992年Beven和Binley[51]認為在模擬水文過程時存在“異參同效”(equifinality)現象,即模型結果受參數組合的綜合影響,不同的參數組合可能獲得相同或相近的模型響應。為詮釋和評估這種參數不確定性的影響,普適似然不確定性評估(Generalized Likelihood Uncertainty Estimation,GLUE)方法被兩位學者首次提出,后續被廣泛用于流域水文模型中,HSPF模型也得益于該方法在參數不確定性方面得以發展。GLUE方法在應用時需設置模型參數的先驗分布和用于評價參數組模擬效果的似然函數,隨后對參數組進行抽樣并重復模擬,結合似然函數閾值以及參數組的似然權重,最終獲得參數組的后驗分布以及模擬結果的不確定性區間。HSPF模型水文參數的先驗分布通常被假設為均勻分布和三角分布,而納什效率系數則常被選為篩選模擬結果的似然函數[36,47,52]。水文參數之間復雜的相關關系(例如下層土壤額定含水量與土壤滲透率呈極顯著負相關關系,上層土壤額定含水量與地下水退水系數呈極顯著正相關關系)也通過GLUE方法的應用被探討[53]。在水質模擬方面,Roostaee和Deng[54]針對硝酸鹽和磷酸根在壤中流、基流和陸面沖刷的相關HSPF模型參數的不確定性進行了探討。Mishra等[55]則利用GLUE方法針對大腸桿菌的陸面沖刷以及一階死亡率參數的不確定性進行了評估和量化。Xie等[36]考慮了場次降雨事件下面源污染過程特征,構建了一系列基于污染物在次降雨尺度下傳輸特征的似然函數并將其融合,用以評估HSPF模型泥沙參數對陸面及河道懸浮泥沙輸移模擬的不確定性。這些前人的研究成果為HSPF模型參數不確定性理論和技術提供了寶貴的經驗。但作為一種“偽”貝葉斯(pseudoBayesian)方法,學術屆對GLUE的一些假設提出了質疑,認為包括非統計意義的廣義似然函數的選取和有效性參數組的判斷篩選等,會導致參數后驗分布的推求過程和結果不具備顯著統計特征[56—57]。基于經典貝葉斯理論的馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法[58]是用于面源污染模型參數不確定性分析的另一種方向。MCMC方法的應用同樣需要確定模型參數的先驗分布,但需定義包含模擬誤差概率分布的嚴格似然函數,構建若干條獨立并行的Markov鏈搜索參數空間,通過參數組的抽樣、模擬和對比評價,不斷更新樣本信息使Markov鏈收斂,實現具有統計意義的最大后驗估計。前有學者對比過GLUE方法和MCMC方法應用于HSPF模型評估大腸桿菌參數的不確定性,結果表明二者推求的參數后驗分布與模擬不確定性區間相似,但GLUE方法需要設定較高的閾值,而MCMC方法需要更多的計算資源[59]。類似的對比結果也出現于該方法在其他流域水文模型中的應用研究中[57,60],發現MCMC方法取得的參數后驗具有略高的識別度,且可以更準確地將參數不確定性貢獻區分于模型總不確定性。后來在高效采樣算法的引入下,Markov鏈的搜索效率得到提升,例如近年來被廣泛用于環境模擬的差分進化自適應Metropolis算法(Differential evolution adaptive metropolis,DREAM)以及其改進的DREAMZS算法[60—61],具備多鏈并行抽樣能力,自適應地調整搜索步長和方向[62],使得MCMC方法比GLUE方法可更高效地收斂至后驗分布,且收斂結果通常具有更好的參數識別度。但目前集成此類高效算法的MCMC方法尚未見報道應用于基于HSPF模型的面源污染模擬研究中。

總體而言,現階段HSPF模型的參數不確定性研究多集中于水文模擬,對水質相關參數的研究比較缺乏。由于HSPF模型的水質參數種類眾多,涉及陸面、河道、水土中遷移、形態轉換等多界面多過程,應借鑒其他面源污染模型應對這種多維度和多模態參數不確定性的先進方法,該方向的研究也將為HSPF模型更好地服務面源污染管控提供理論和技術支撐。

3.2.2輸入不確定性

降雨和DEM是面源污染模擬的關鍵輸入數據。針對HSPF模型來說,Wang等[63]認為二者的不確定性會傳遞至徑流、泥沙及氮磷污染模擬并具有一定的擴大效應,該結論也與其他面源污染模型的相關研究結果較為一致[64]。目前降雨數據的獲取仍以地面雨量站網觀測為主,由于站網布設不合理以及測量誤差,會對面源污染模型的模擬結果帶來較大的不確定性[65]。增加雨量站點的數量和密度可提高降雨空間異質性的描述,從而提高徑流和面源污染模擬的精度,但雨量站點稀疏是很多流域建模時遇到的普遍現象。HSPF模型在BASINS平臺中建模時,默認選用的是臨近法或泰森多邊形法對子流域和雨量站進行匹配,但對于復雜地形流域,可能需要考慮高程、坡度、坡向等要素對降雨插值的影響。Huo等[66]采用信息熵和模糊數學方法,并在考慮了高程數據之后為HSPF模型生成了降雨輸入數據,相比常用的泰森多邊形法增強了對降雨空間異質性的描述從而提高了徑流模擬的可靠性。此外,為分析降雨數據本身的測量隨機誤差的影響,劉興坡等[67]以青龍河流域為研究區,區分降雨時間序列的隨機成分并建立隨機模擬序列驅動了HSPF模型的徑流模擬。衛星降水產品由于具有高時空采樣頻率和覆蓋率,也被逐漸用于面源污染模擬的輸入。Solakian等[68]基于地面高密度定點雨量站觀測數據評價了三種衛星降水產品(TRMM、CRMORPH和PERSIANN)的差異,并通過HSPF模型傳遞衛星降雨輸入不確定性對流量、懸浮泥沙、水溫和溶解氧模擬的影響,分析了這種影響在不同模擬空間尺度和季節間的變化。

DEM數據分辨率影響著坡度、坡長、河道坡降等流域關鍵特征參數的提取,進而改變了相關參數在模型中的輸入。DEM數據分辨率帶來的不確定性對流域水文及面源污染模擬的重要影響也已在前人的研究中得到共識[69—70]。針對HSPF模型,Roostaee和Deng[71]的研究發現當DEM數據分辨率在3.5 m至100 m區間變化時,地勢平緩地區的徑流模擬較為敏感,月均流量的相對誤差可高達50%；而在坡度陡峭的山區,這種輸入不確定性對徑流模擬影響非常小。Wang等[63]通過對DEM柵格數據進行重采樣獲取30 m至480 m的16種分辨率,評估不同面源污染因子對DEM分辨率變化的敏感性。他們的研究結果發現模型對硝酸鹽和氨氮的模擬受其不確定性影響最大,且高于對徑流模擬的不確定性影響。這表明DEM數據輸入的不確定性在徑流-面源污染模擬過程中具有傳遞和擴大效應。

3.2.3結構不確定性

針對模型結構的不確定性分析常用到多模型集合模擬的方法,即加權平均多模型結果以減少單模型由于其對物理過程概化不足造成的不確定性。權重的選取在早期的研究中較為主觀且缺乏物理意義,近年來貝葉斯模型平均(Bayesian model averaging, BMA)方法[72]多被應用于面源污染模擬的結構不確定性評價研究,該方法通過似然值代替了權重,具有更高的魯棒性。Wang等[73]在我國灤河流域評價了HSPF模型模擬氮素污染的結構不確定性,通過集合SWAT模型和GWLF模型并采用BMA方法對多模型結果進行了貝葉斯加權平均。該方法也被用于減少HSPF模型及SWAT模型單獨模擬氮磷污染的結構不確定性,通過置信區間的寬度及對觀測值的覆蓋率等不確定性量化方法提高了模型的解釋力[37]。此外,加強模型對機理過程的數學描述方式和方法也是減小模型結構不確定性的一種方向[74],但針對HSPF模型鮮有相關的研究。雖然HSPF模型已經在全球范圍被廣泛認可,它的模型結構在某些特定現實情景下還有提升的空間。例如在干旱區域非汛期,地表水和地下水的交互作用突出,而HSPF模型尚不能準確表達這種耦合情景。針對水稻田這種特殊的流域地理單元,已有研究對HSPF模型的結構進行了改進,在模型的源代碼中引入了圩埂高度和Dirac delta函數以表征水稻田內排水、灌溉和施肥過程,構建了針對水稻灌區的HSPF-paddy模型并模擬了徑流、水溫、五日生化需氧量、總氮及總磷。該改進模型的模擬結果在田塊尺度和流域尺度均得到了驗證[75—76]。

4 HSPF模型在面源污染管控中的應用

4.1 BMPs和LIDs效果評估

最佳管理措施(Best management practices, BMPs)和低影響開發(Low impact developments practices,LIDs)是能夠有效削減面源污染的管理政策或技術手段[77]。為避免盲目的規劃和實施管控方案,需要科學精確地評估BMPs和LIDs的污染物削減效率。面源污染模型是常用的評估工具,研究學者通常采用模型中的措施評估模塊或調整措施評估參數的方式,在模型中反映BMPs和LIDs實施的情景并對其效果進行模擬分析。模型模擬相較于實地跟蹤監測實驗具有較大的靈活性和時效性[78]。BMPs和LIDs的實施單元通常為場地、田間、小區等,通過流域面源污染模型模擬,可將小尺度下的評估結果推廣至子流域和流域尺度,也有助于分析管控措施的尺度效應,因此面源污染模型在評估措施削減效果時具有較大的優勢[79]。

針對管理型措施,HSPF模型中的SPEC-ACTION模塊可用于評估其削減效果。用戶在該模塊中可設置和調整農業耕作時間、收割時間、化肥施用方案等農業活動的情景,進而模擬免耕、少耕、施肥管理等管理型措施的水質改善效應。針對工程型措施,HSPF模型中內置了BMPRAC模塊,用戶可自主定義單個工程BMP或LID措施(包括植草溝、過濾帶、干濕沉淀池、岸邊緩沖帶)在流域單元中的數量及污染物去除率[80]。Risal等[81]采用該模塊模擬了農田中過濾帶及沉淀池在流域尺度下對泥沙、總氮和總磷的削減效率。但BMPRAC模塊中直接設置措施效果的處理方式具有主觀性,且不能反應措施受降雨、污染源強、地形等多因素影響下削減效果的時空異質性。相比之下,調整措施評估參數則是一種過程型的模擬方法,可以充分結合面源污染模型本身的結構框架,利用模型對污染物遷移轉化的描述能力評估措施的污染物削減效果。謝暉[82]分析了四種農業BMPs(免耕、等高種植、梯田和植被過濾帶)對面源污染物的削減過程和原理,將措施在流域單元中的數學表達與HSPF模型參數關聯,以此評估了三峽庫區小流域中四種措施對面源氮磷污染的削減能力。但是這種過程型的措施評估方法與HSPF模型結合的研究較少,相比之下,SWAT模型的在該方向上具有廣泛的應用成果。這種差別主要是由于面向用戶的易操作性以及應用經驗積累的差距。Arabi等[83]首次系統提出了SWAT模型評估農業管控措施的過程型模擬方法；美國德州水資源學院在2011年編撰了基于SWAT模型的管控措施模擬手冊[84]；專門用于評估LIDs效果的SWAT-LID模塊也被研發[85],這些學者的研究經驗和成果都有助于SWAT模型在該領域的推廣應用。對于HSPF模型來說,它與SWAT模型有很多共通之處,包括流域地理單元的概化、機理過程的數學表達等等,后續的研究可借鑒SWAT等其他面源污染模型的研究經驗,在措施評估模塊的開發上有所突破,將有利于HSPF模型高效準確地評估BMPs和LIDs的水質改善效應。

4.2 流域污染物總量控制

流域尺度下的污染物總量控制,通常是指在水質標準的約束下,充分考慮社會經濟條件和水體自凈能力,控制流域允許排放的最大污染物總量,并對各單元需要削減的污染物負荷進行合理分配[86]。這種水環境管理的理念在美國TMDL計劃、歐盟水框架指令中得以發展,并且逐漸在我國多個流域內推廣實踐。

污染物總量控制須要明確并量化區域污染物排放對水體水質狀態的影響,流域面源污染模型則可建立這種驅動-響應關系,并從降雨-徑流、面源污染遷移轉化、河流水質水動力等多角度分解污染物排放過程。HSPF模型正是當前美國實施TMDL計劃的推薦模型,并在Cannon River、Chesapeake Bay、Moosy Creek等多個流域開展了總量控制的可行性研究和實踐[13,80,87]。此外,HSPF模型的不確定性對污染總量控制的影響也被研究。Mishra等[55]區分考慮了HSPF模型的隨機不確定性和認知不確定性,利用兩相蒙特卡洛隨機模擬的方法評估了大腸桿菌污染總量控制的情景方案。Ahmadisharaf和Benham[88]結合GLUE方法和HSPF模型,設計了基于水質目標風險的大腸桿菌污染削減情景評價框架,以正向模擬的思路評估了各污染源在不同削減比例下污染物總量超標的概率。另一方面,在一定條件約束下,以反向優化的思路將污染削減量及控制措施分配于不同空間單元成為總量控制領域的研究熱點[89—90]。如前所述,HSPF模型具備靈活的面源污染模擬和措施效果評估能力,能夠模擬污染削減措施在不同尺度下的環境效益,因此可將其與優化算法或工具進行耦合,實現成本投入-環境效益協同目標下的污染物總量控制。在美國San Antonio River流域內的多個城市型流域,HSPF模型在優化工具SARA Load Reduction Tool和SARA Enhanced BMP Tool的結合下服務了該區域的污染削減和雨洪管理規劃[91]。SARA Load Reduction Tool工具可優化分配不同區域的污染負荷削減量,SARA Enhanced BMP Tool可確定控制措施的優化配置(種類、數量和位置等)以確保分配的污染負荷削減量得以實現,HSPF模型則是不同方案制定和反饋迭代過程中的核心過程與紐帶。基于反向優化的措施空間配置過程同樣受到模型不確定性的影響,Xie等[38]將HSPF模型的參數不確定性引入措施優化配置過程中,在多種水質達標的置信水平下開展了“措施效果評估-不確定性分析-措施優化配置”的一體化研究,研究思路及相關成果服務于減少農村小流域污染物總量控制的決策風險。

5 新形勢下HSPF模型的延伸

隨著近年來流域水環境綜合管理從理論轉向實踐,國家和地方逐漸推動流域水環境綜合管理從信息化向智慧化發展,智慧流域的理念及其信息化平臺建設成為重要的方向。地方生態環境部門不僅需要水利、水環境等監測數據的在線掌控,更希望通過仿真模擬實現河湖水質的實時模擬、預報預警和污染源解析等功能,因此對面源污染模型在信息化平臺中的集成應用與服務需求愈發強烈[92—93]。HSPF模型在現階段的操作使用非常依賴BASINS平臺,雖然BASINS平臺集成了該模型的輸入數據處理、模型配置及結果展示等功能,但BASINS的大部分圖形用戶接口是由Visual Basic 6編譯,這與微軟Windows現在的主流開發平臺Visual.Net并不兼容,且并不能在Windows以外的操作系統(例如類Unix系統)下運行[94]。這些外在因素減弱了HSPF模型自身固有的靈活性以及在信息化平臺的可操作性。而對于HSPF模型本身,其最新版本的源代碼依然保留了大量Fortran 77語法,內存管理復雜且固化,程序范式嚴謹但可讀性較差,一定程度上限制了該模型的維護、升級和移植。HSPF模型運行的核心輸入文件是WDM文件,存儲空間有限,且創建和修改需要訪問HSPF模型鏈接庫中的許多子程序,用戶通常要借助BASINS平臺中的組件如WDMUtil、SARA Timeseries Utility等特定工具進行編輯。以上這些HSPF模型固化的特征與現代化環境模擬要求的可復制性和開放性相悖[95],也正是HSPF模型亟需發展的方向。

在流域管理信息化平臺建設和現代化環境模擬的背景下, HSPF模型在不同學者和技術人員的共同努力下得到了延伸和發展。PyHSPF[94]是最早的典范,PyHSPF是以Python語言為核心的開源跨平臺軟件套裝,它沿襲了HSPF version 12.2中水和水質模擬的函數庫,利用Python聯接了函數庫、子流域屬性和驅動數據,并構建了包含數據預處理、參數自動率定、批量模擬等功能的Python模塊和類,為提升HSPF模型模擬效率[96]、模塊功能自主研發[97]以及應用編程接口開發[95]提供了多種便利。此外,HSP2(Hydrologic Simulation Program-Python)模型作為HSPF模型的Python版本被研發,該模型使用Python重編譯了HSPF模型的核心功能模塊,使用跨平臺數據存儲文件HDF5存儲和讀寫模型的輸入和輸出數據,并利用交互式Web應用Jupyter Notebook和JupyterLab 維護HSP2的開發、升級和協作[98]。HSP2模型的主要優勢在于其模型結構和源代碼脫離了特定操作系統和WDM文件的束縛,增強了子流域、土地利用等流域屬性和模型參數等的動態表達,從而有利于用戶設計多種變化環境的驅動情景,且模型代碼可與CPU/GPU異構式并行加速環境兼容。雖然HSP2模型尚處于前期開發階段,已有研究成功利用該模型在東南亞源頭小流域模擬了河流大腸桿菌濃度的高頻動態變化[99]。

6 HSPF模型的發展動向及展望

縱觀HSPF模型的發展歷程,該模型自上世紀80年代被研發后經歷了蓬勃的發展和積累。現如今HSPF模型是美國環保署重點推薦的流域水質模型,是各州執行TMDL計劃的主要工具之一。雖然該模型在國內的起步較晚,但從最早的徑流模擬已擴展至集成多種污染物、多過程的水量水質綜合模擬,到如今服務于太湖、巢湖、滇池等流域的污染管控和規劃工作,可以說HSPF模型在我國國情和流域背景下的應用有了飛速的發展和實質性突破。盡管HSPF模型具備了多種優點,但現階段依然存在諸多局限性,作為面源污染模擬的重要工具,該模型的后續完善更應該落腳至解決生態環境領域的實際問題,這就更加需要廣大科研工作者的投入和堅持。對于未來HSPF模型的發展動向及展望,筆者認為可以著重關注以下幾個方面：

(1)模型本地化改進。HSPF模型雖然已在全球范圍內得到成功的應用,但是針對特定的應用場景,如果完全依賴現階段模型本身的結構功能,還不足以支撐我國河湖管理工作。例如化學需氧量是我國地表水環境質量考核的重要指標,但HSPF模型暫未具有直接模擬化學需氧量的功能模塊。此外,閘門和泵站是人類活動干擾流域自然系統的常見水工構筑物,在我國平原圩區十分常見,對此HSPF模型的模擬概化步驟并非友好。因此在未來,應該針對特定的應用場景,改進模型對實際過程的數學表達,或開發全新的功能模塊、或耦合其他環境模型,提升HSPF模型在我國流域的適用性。

(2)大尺度精細化模擬。HSPF模型已在大尺度流域(例如166, 000 km2的Chespeake Bay)[64]獲得了成功應用,但在大尺度模擬的同時還需考慮計算單元的細化精度。在我國河長制的背景下,計算單元的劃分應結合行政區劃,服務河長制的網格化精細管理,計算單元由粗糙單元向精細單元轉變是大趨勢。對于HSPF模型而言,大尺度精細化模擬需要提高計算效率,基于Python的多線程模塊和并行計算擴展包如PyCUDA和Numba可以與pyHSPF模型和HSP2模型兼容[96],未來可設計并行架構的云計算環境,將有助于HSPF模型在大尺度復雜流域的應用。

(3)人工智能和模型的互饋集合模擬。當下數字水利、數字流域的建設促進了多源觀測資料呈幾何級增長,這些基礎大數據在人工智能的催化下,可幫助解決HSPF模型所需的高精度輸入以及多維度參數自適應尋優問題。此外在少資料地區,流域模型的模擬結果可為人工智能算法提供訓練樣本[100],從而延伸數據驅動模型的研究范式。因此對于HSPF模型而言,未來應充分結合大數據統計與智能算法,打破學科間隔閡,互饋集合應用以提升面源污染的模擬與預測水平。