鄱陽湖湖泊流域系統(tǒng)水文水動力聯(lián)合模擬*

李云良，張 奇，姚 靜，李相虎

(1:中國科學(xué)院南京地理與湖泊研究所湖泊與環(huán)境國家重點實驗室，南京 210008)

(2:中國科學(xué)院大學(xué)，北京 100049)

鄱陽湖位于長江中下游南岸，在優(yōu)質(zhì)淡水供給、洪水調(diào)控和生物多樣性保護等方面起著重要作用[1].鄱陽湖流域主要由贛江、撫河、信江、饒河和修水五大子流域和湖區(qū)平原區(qū)構(gòu)成.鄱陽湖接納流域來水，調(diào)節(jié)后由北端湖口注入長江(圖1a).流域內(nèi)降雨和蒸發(fā)具有顯著的時空變化，流域人類活動強烈，極端氣候事件頻發(fā)，多種因素的相互疊加，致使流域內(nèi)洪澇災(zāi)害頻發(fā)[2-4].特別是極端氣候事件的發(fā)生，影響了流域水資源的時空分布和湖泊水量變化.已有研究表明，流域五大水系的入湖徑流量對湖泊水位起著主控作用，其次，鄱陽湖與長江季節(jié)性的水量交換也對湖泊水位有著不同程度的影響作用[5-6].湖泊水位的變化，導(dǎo)致鄱陽湖水面面積呈現(xiàn)明顯的萎縮與擴張，嚴(yán)重威脅湖區(qū)平原區(qū)農(nóng)業(yè)生產(chǎn)和管理及濕地生態(tài)系統(tǒng)健康等方面[5].鄱陽湖水情和水位變化，主要受季節(jié)性降雨和流域人類活動及長江來水等多因素影響，這些因素的疊加影響了流域降雨-徑流過程和湖泊泄水條件，從而影響湖泊的水量平衡[6-7].由此可知，鄱陽湖湖泊流域之間具有顯著的水文水動力聯(lián)系，主要表征為流域來水對鄱陽湖水文的驅(qū)動作用.在目前極端氣候事件頻發(fā)和人類活動加速背景下，研究鄱陽湖水位和水量對流域徑流過程的響應(yīng)成為一個緊迫的任務(wù).為此，本文提出一個聯(lián)合模擬方法，將作為一有效工具模擬分析鄱陽湖水文水動力對流域徑流的響應(yīng).

鄱陽湖湖泊流域系統(tǒng)自然屬性復(fù)雜，如流域尺度較大(16.2×104km2)，多河入湖，湖泊岸線及湖盆地形復(fù)雜，湖泊水動力過程受流域與長江的共同作用，這些都為水文水動力過程的模擬帶來挑戰(zhàn).水文或水動力模型往往在單一的流域或地表水體(水庫、湖泊等)具有較強的模擬能力，但多數(shù)情況下與周圍水系割裂，無法切實模擬復(fù)雜水文水動力系統(tǒng)[8].本文所提出的基于不同模擬功能的子模型的有效聯(lián)合擬解決鄱陽湖湖泊流域的水文水動力整體模擬.針對鄱陽湖及其流域，已有眾多研究成果且強調(diào)不同的目的[6-7，9-11].大部分研究主要采用水文模型側(cè)重于流域降雨-徑流過程的模擬計算[7，10]及其對氣候變化和人類活動的響應(yīng)[6]，但均未涉及到湖泊的相關(guān)研究以及湖泊對流域來水的響應(yīng)，由于流域與湖泊作用關(guān)系密切，單方面的流域模擬研究顯然是不足的.葉許春采用分布式水文模型WATLAC 探討未來不同的氣候變化情景對鄱陽湖流域徑流的可能影響，繼而采用湖泊水量平衡原理，計算流域徑流變化對湖泊水位的影響.該方法沒有考慮湖泊的水動力特性，只能在平均意義上描述湖泊水位的變化，存在一定的局限性[11].國內(nèi)外關(guān)于鄱陽湖湖泊流域系統(tǒng)的水文水動力聯(lián)合模擬研究幾乎沒有.本文的主要研究目的是，通過不同子模型的聯(lián)合來實現(xiàn)大尺度鄱陽湖湖泊流域系統(tǒng)的水文水動力整體模擬;闡述該聯(lián)合模擬模型的驅(qū)動方法、率定及其潛在的用途.本文提出的方法也可用于其它類似湖泊流域的水文水動力過程研究，為定量模擬湖泊水情對氣候變化和人類活動的響應(yīng)提供有效工具.

1 材料與方法

1.1 基礎(chǔ)數(shù)據(jù)與準(zhǔn)備

模型需要的主要數(shù)據(jù)包括:流域水文站點河道日徑流量數(shù)據(jù)、湖泊站點日水位數(shù)據(jù)、日氣象資料(降雨和蒸發(fā)皿數(shù)據(jù))、土地利用數(shù)據(jù)、土壤物理屬性數(shù)據(jù)(田間持水量、飽和滲透系數(shù)及總孔隙度)、葉面指數(shù)(來自MCD15A2)、地面高程數(shù)據(jù)(DEM)、分辨率為30 m×30 m 湖泊地形高程數(shù)據(jù)以及流域與湖泊邊界矢量圖等資料.降雨和蒸發(fā)站點(圖1a)數(shù)據(jù)按最近鄰法插值到流域和平原區(qū)計算單元，作為水文和平原區(qū)產(chǎn)流模擬的驅(qū)動因素;水文模型中土地利用(2005年)劃分6 種類型:耕地、林地、草地、水體、居民用地和其它非建筑用地;土壤總孔隙度、田間持水量、飽和滲透系數(shù)等均為1 km×1 km 空間柵格數(shù)據(jù)，可直接輸入模型，作為土壤物理屬性參與計算;基于遙感反演的葉面指數(shù)數(shù)據(jù)，因其在日尺度上變化相對較小，水文模型按月尺度輸入計算.湖泊水動力模型所需的主要數(shù)據(jù)為湖泊地形數(shù)據(jù)和岸線邊界(圖1b)，水動力模型同樣按最近鄰法獲取站點的降雨和蒸發(fā)資料.

1.2 模型描述

圖1 鄱陽湖湖泊流域系統(tǒng)結(jié)構(gòu)(a:湖泊流域平原區(qū)結(jié)構(gòu)簡圖及水文氣象站點空間分布;b:鄱陽湖岸線邊界及湖泊地形)Fig.1 Lake Poyang-catchment system (a:schematic diagram for lake-catchment-plain area and spatial distribution of hydro-climate stations;b:shorelines and bathymetry for Lake Poyang)

WATLAC 模型[12-13]是基于格網(wǎng)的以日為步長、聯(lián)合地表-地下徑流的分布式流域模擬模型，能針對湖泊集水域的特點，將相對獨立的多個子流域、多條入湖河流在同一個模型中加以模擬，無需對每個子流域分別建立獨立的模型，并便于對整個流域進行全局模型參數(shù)率定.其模擬的主要水文循環(huán)過程為:冠層截留與蒸散發(fā)、坡面流、河道匯流演算、土壤壤中流、土壤蒸發(fā)與滲漏補給、地下水蒸發(fā)、河流與地下水的相互作用等.有限差分方法MODFLOW-2005[14]的植入使得該模型能夠刻畫較為復(fù)雜的地下水流環(huán)境，是一個物理機制較為明確的地表水-地下水流實時耦合的流域模擬模型.詳細(xì)計算方法這里不再闡述，具體的模型耦合機制及數(shù)學(xué)方程已于 2009年發(fā)表[12-13]且成功應(yīng)用于不同的流域[7，9-13].

平原區(qū)產(chǎn)流的計算主要基于傳統(tǒng)的水文學(xué)方法，采用基于概念性的平原區(qū)產(chǎn)流估算方法，即降雨乘以徑流系數(shù)(流域多年統(tǒng)計結(jié)果)來獲取產(chǎn)流量.平原區(qū)這種近似估算方法雖難以得到有效驗證，但該方法已成功應(yīng)用于其它流域[15]，在缺乏資料的情況下，可作為平原區(qū)入湖徑流的近似估算.本文主要目的并不是關(guān)注于該方法計算的平原區(qū)入流過程的精度，至少提出的概念性評估方法彌補了以往平原區(qū)產(chǎn)流計算的空白，確保聯(lián)合模型的完整性及結(jié)構(gòu)的合理性.

MIKE 21 模型[16]是丹麥水力研究所(DHI)開發(fā)的系列水力學(xué)軟件之一，屬于平面二維水流數(shù)學(xué)模型，該模型可用于模擬河流、湖泊、河口、海灣、海岸等大型水體流場特征.因其具有完全的物理機制且具有靈活的三角形網(wǎng)格剖分，所以能夠很好適應(yīng)鄱陽湖的寬淺型地形特點和復(fù)雜岸線特征.該模型已成功應(yīng)于不同國家及地區(qū)[17-18]，詳細(xì)的模型理論及方法，這里不再闡述.

為了保證模型空間離散分辨率與大多數(shù)基礎(chǔ)數(shù)據(jù)初始分辨率一致且能更好反映下墊面自然屬性特征，WATLAC 模型的離散空間分辨率選擇較為精細(xì)的1 km×1 km 格網(wǎng)單元，共剖分138634 個計算單元，計算域(圖1a，DEM 覆蓋區(qū)域)占整個湖泊流域系統(tǒng)面積比重約為87%.模型中僅考慮潛水含水層的模擬且概化為一層，主要用來計算河流與含水層的水量交換.其它模型詳細(xì)的構(gòu)建過程與Ye 等方法相同[9]，主要不同之處除了當(dāng)前流域模型采用更為精細(xì)的格網(wǎng)外，最主要的是當(dāng)前模型充分考慮平原區(qū)這部分產(chǎn)流的模擬計算.

平原區(qū)產(chǎn)流模擬計算域(圖1a，灰色區(qū)域)同樣采用1 km×1 km 格網(wǎng)，共剖分17556 個單元，占整個湖泊流域系統(tǒng)面積比重約11%.由于鄱陽湖流域降雨-徑流關(guān)系比較密切，本文采用降雨乘以多年平均徑流系數(shù)(α=0.58)的方法來估算平原區(qū)入湖徑流量.徑流量的估算大致分為以下兩步:首先，將日氣象站點的降雨資料按最近鄰法插值到整個平原區(qū)格網(wǎng)單元，將格網(wǎng)單元的降雨量乘以多年平均徑流系數(shù)，便可得格網(wǎng)的產(chǎn)流量;其次，假定平原區(qū)沿著湖區(qū)邊界的分布式入流與點入流有著等同的影響效果，故將計算格網(wǎng)產(chǎn)流采用最近鄰原則分配至5 個子流域出口，逐時段計算，便可得五大子流域與平原區(qū)的日合成徑流序列過程，充分體現(xiàn)全流域模擬的完整性.

鄱陽湖水動力模型計算域范圍約為東西向100 km，南北向170 km，占湖泊流域系統(tǒng)的面積比重約2%.模型采用基于有限體積法的三角形網(wǎng)格以適應(yīng)復(fù)雜的湖泊地形及岸線邊界(圖1b).靈活地變網(wǎng)格大小為70 ～1500 m，共剖分網(wǎng)格數(shù)為20450.網(wǎng)格的合理剖分及大小不僅適應(yīng)湖區(qū)內(nèi)河道深窄且流速較快、平原寬淺且水流較慢的地形特點，也提高了水動力模型的計算效率.為了更加合理的模擬，模型將站點降雨、蒸發(fā)資料考慮到水動力模型中參與計算.以湖泊各站點水位空間插值結(jié)果作為初始水位場，初始流速設(shè)定為0.為了合理描述湖泊地形特點，變曼寧數(shù)給定為30 ～50 m1/3/s 分別適應(yīng)于湖泊底部的平原區(qū)和主河道.渦粘系數(shù)采用Smagorinsky 公式[19]計算，其中Smagorinsky 因子Cs 取為0.28.模型中通過定義臨界水深來確定干濕單元，很好地去適應(yīng)湖泊水面面積隨水位變化顯著的特點.

1.3 模型框架

為了充分體現(xiàn)大尺度復(fù)雜聯(lián)合系統(tǒng)的完整模擬，湖泊水動力模型的上游邊界選在5 大子流域(贛江、撫河、信江、饒河、修水)出口斷面處(圖1b)，以5 大子流域的日模擬徑流(約占湖泊88%貢獻量)及平原區(qū)入湖徑流(約占湖泊12% 貢獻量)的合成流量以5 個點入流的方式作為水動力模型的上邊界條件，其反映了流域-湖泊相互作用關(guān)系;模型下邊界采用實測湖口水位過程線，其反映了長江-湖泊相互作用關(guān)系;其它均按閉邊界條件給定.3 個子模型的連接采用輸入-輸出的外部耦合技術(shù)將流域WATLAC 模型、平原區(qū)產(chǎn)流模型以及湖泊水動力模型MIKE 21 聯(lián)合起來，形成一個流域-平原區(qū)-湖泊結(jié)構(gòu)的水文水動力聯(lián)合模擬框架(圖2).

圖2 湖泊流域系統(tǒng)水文水動力聯(lián)合模型驅(qū)動關(guān)系及模擬框架Fig.2 Driving relationship and framework of integrated hydrological and hydrodynamic models for Lake Poyang-catchment system

2 率定策略

選用2000年1月1日-2005年12月31日作為整個聯(lián)合模型的率定期.通過6 個水文站點(圖1a)的河道日徑流量來率定WATLAC 模型的地表水流部分.將上述6 個水文站點的日徑流量進行基流分割[20]，并采用加權(quán)平均值法計算整個流域的均基流指數(shù)(基流量/河道流量)，以此來評估地下水模擬效果.總之，以6個站點的河道徑流量和地下水基流指數(shù)作為多目標(biāo)函數(shù)來率定分布式水文模型.

在當(dāng)前聯(lián)合模型參數(shù)率定方面，PEST[21]參數(shù)自動率定技術(shù)主要用于WATLAC 模型，這是因為其能夠與WATLAC 模型較為方便地建立數(shù)據(jù)及參數(shù)傳遞接口.PEST 是一個非線性參數(shù)估計工具，其原理是采用GML算法來優(yōu)化模型，通過盡量少的迭代和模型調(diào)用次數(shù)使得計算值與觀測值之間的殘差平方和達到最小，以此來尋求最優(yōu)參數(shù)值.PEST 自動率定技術(shù)的嵌入使用，能夠提高模型的運行效率和整體性能，加速模型率定進程和增強模擬結(jié)果可靠性，同時也為WATLAC 模型大尺度流域、長時間連續(xù)模擬提供保障.WATLAC模型率定的主要敏感性參數(shù)為坡面流滯后系數(shù)、馬斯京根法的時間蓄量常數(shù)、地下水補給率、降雨入滲因子及潛水含水層滲透系數(shù)和給水度[12-13].鑒于PEST 優(yōu)化技術(shù)對模型參數(shù)的率定已取得成功應(yīng)用[10]，故本文沒有列出具體的參數(shù)優(yōu)化結(jié)果.而水動力模型MIKE 21 仍然采用手工試錯法率定，模型可調(diào)整的主要參數(shù)為湖床糙率系數(shù)和渦粘系數(shù).

3 結(jié)果與討論

3.1 模型率定與評估

當(dāng)前聯(lián)合模型采用基于觀測數(shù)據(jù)來獨立率定各個子模型.水文水動力聯(lián)合模型分別以五河6 個水文站點河道徑流量、流域基流指數(shù)與湖泊4 個站點水位作為目標(biāo)變量來評價模型的整體能力.

表1 為基于多目標(biāo)變量的WATLAC 模型定量化模擬結(jié)果評價.6 個水文站點擬合的納希效率系數(shù)[22](Ens)變化范圍為0.71 ～0.84，確定性系數(shù)(R2)介于0.70 ～0.88 之間，相對誤差(RE)基本控制在±10%內(nèi)，但個別站點擬合的相對誤差稍微偏大(表1).模擬值與觀測值大部分集中在45°線附近，表明模型模擬效果較好(圖3).此外，通過圖4 選取的模擬期末2005年河道日徑流量過程線可見，總體呈現(xiàn)較好的曲線擬合宏觀效果，但個別時段洪峰流量擬合存在一定的誤差，這主要是因為流域存在眾多大小型水庫改變了徑流的天然狀態(tài).WATLAC 模擬的全流域基流指數(shù)(45%)與基流分割結(jié)果(47.6%)較為接近，表明地下水對河道徑流基流貢獻量的模擬結(jié)果具有一定的可靠性，基流指數(shù)可用來直接率定分布式水文模型及間接評估地下水流模擬效果.總體來說，流域水文模型WATLAC 能夠很好地再現(xiàn)大尺度流域上的降雨-徑流響應(yīng)過程，能為湖泊水動力模型提供可靠的流量輸入條件.

表1 流域水文模型WATLAC 率定結(jié)果*Tab.1 Calibration results of catchment hydrological model WATLAC

圖3 模型率定期(2000-2005年)模擬與觀測的河道日徑流量散點圖Fig.3 Scatter diagrams of simulated stream flow against the observed stream flow for calibration period 2000-2005

湖泊MIKE 21 模型的模擬與觀測水位的定量化評估結(jié)果表明，湖泊各站點擬合的確定性系數(shù)R2介于0.96 ～0.98，納希效率系數(shù)Ens變化范圍為0.88 ～0.98，相對誤差RE 均小于±5%，取得令人滿意的率定結(jié)果(表2).鄱陽湖4 個水文站點實測與計算水位過程線擬合良好，均很好地呈現(xiàn)了湖泊水位的季節(jié)性變化和年際變化趨勢(圖5).總體來說，湖泊水動力模型MIKE 21 不但能夠取得理想的模擬效果，而且具有較強的模擬能力去再現(xiàn)長時間的水位變化及捕捉水位的峰值變化，充分體現(xiàn)該模型能夠很好的模擬湖泊水位對流域和平原區(qū)徑流量的共同響應(yīng)過程.盡管如此，在枯水位的模擬上有著相對較大的誤差(圖5)，誤差來源最可能是湖盆地形的誤差或者說地形誤差對低水位的模擬較為敏感.

圖4 2005年鄱陽湖流域6 個水文站點模擬與觀測的河道日徑流量擬合Fig.4 Comparison of simulated and observed daily stream flows at six hydrological stations in Lake Poyang catchment in 2005

圖5 2000-2005年鄱陽湖4 個水文站點水位驗證Fig.5 Validation of lake levels at four hydrological stations in Lake Poyang during 2000-2005

圖6 2005 年4 個典型時段模擬鄱陽湖水位的空間變化Fig.6 Simulated spatial distribution of Lake Poyang water levels at four time slices in 2005

表2 湖泊水動力模型MIKE 21 率定結(jié)果*Tab.2 Calibration results of lake hydrodynamic model MIKE 21

圖7 2005年鄱陽湖淹水天數(shù)統(tǒng)計結(jié)果Fig.7 Statistical result of submerged days of the entire Lake Poyang area in 2005

3.2 鄱陽湖空間水位模擬

鄱陽湖水位在年內(nèi)的空間上存在著明顯的水位梯度，即水位空間差異性較大，特別是低水位月份，比如1月份(圖6a)，只有湖泊主河道有少量的水分布，大部分區(qū)域水深較淺，呈現(xiàn)出露狀態(tài).而水位相對較高的7月份(圖6c)，整個湖泊保持著較高的水位，且整個湖區(qū)近似呈水平面.此外，本文將2005年的逐日空間水位分布結(jié)果進行統(tǒng)計，將水深大于零的格網(wǎng)單元定義為湖區(qū)淹水狀態(tài)，反之，定義為出露狀態(tài)，空間結(jié)果見圖7.鄱陽湖整個湖區(qū)的淹水情況同樣呈現(xiàn)較大的差異性，不同的湖區(qū)有著不同的淹水情況(圖7).常年被水淹沒的區(qū)域主要分布在主河道地區(qū)，該淹水分布與低水位的空間分布相似(圖6a)，表明這些主河道區(qū)域常年有水的流動.而湖泊岸線附近的上游區(qū)域則呈現(xiàn)1 ～6 個月不等的淹水時間，表明這些地區(qū)由于湖盆地形較高，一年中的絕大部分時間則處于出露狀態(tài).圖6 和圖7 充分體現(xiàn)了鄱陽湖在不同季節(jié)河湖相交替的獨特水情動態(tài)，也表明當(dāng)前所構(gòu)建的聯(lián)合模擬模型能夠合理呈現(xiàn)鄱陽湖水位在年內(nèi)的空間變化，而這種水位空間動態(tài)變化及淹水信息對評估湖泊濕地系統(tǒng)、水生植被的生長和發(fā)展等有著重要的現(xiàn)實意義.

4 結(jié)語

鄱陽湖湖泊流域系統(tǒng)涉及尺度大、多個獨立子流域入湖、湖泊水動力特征顯著等特點，系統(tǒng)內(nèi)水文水動力過程復(fù)雜.本文嘗試了一個實現(xiàn)該系統(tǒng)水文水動力的聯(lián)合模擬的方法.該方法基于不同功能子模型的有機聯(lián)系，通過數(shù)據(jù)在模型間的時空傳遞實現(xiàn)了湖泊和流域的聯(lián)合模擬.該聯(lián)合模擬方法在鄱陽湖湖泊流域作了驗證，能有效反映湖泊水位對流域徑流過程的響應(yīng).本文的研究雖然針對鄱陽湖湖泊流域系統(tǒng)，所提出的方法具有一定的普適性，可為相似湖泊流域的水文水動力模擬提供借鑒或被直接采用.

[1]周文斌，萬金堡，姜加虎.鄱陽湖江湖水位變化對其生態(tài)系統(tǒng)影響.北京:科學(xué)出版社，2011.

[2]Shankman D，Qiao LL.Landscape changes and increasing flood frequency in China's Lake Poyang region.The Professional Geographer，2003，55(4):434-445.

[3]Shankman D，Heim BD，Song J.Flood frequency in China's Lake Poyang region:trends and teleconnections.International Journal of Climatology，2006，26:1255-1266.

[4]Shankman D.River management，landuse change，and future flood risk in China's Lake Poyang region.River Basin Management，2009，7(4):423-431.

[5]Hu Q，F(xiàn)eng S，Guo H et al.Interactions of the Yangtze River flow and hydrologic processes of the Lake Poyang，China.Journal of Hydrology，2007，347:90-100.

[6]郭 華.氣候變化及土地覆被變化對鄱陽湖流域徑流的影響[學(xué)位論文].南京:中國科學(xué)院南京地理與湖泊研究所，2007.

[7]劉 健，張 奇，左海軍等.鄱陽湖流域徑流模型.湖泊科學(xué)，2009，21(4):570-578.

[8]Xu ZY，Adil NG，Thomas JG.The hydrological calibration and validation of a complexly-linked watershed-reservoir model for the Occoquan watershed，Virginia.Journal of Hydrology，2007，345(3/4):167-183.

[9]Ye XC，Zhang Q，Bai L et al.A modeling study of catchment discharge to Lake Poyang under future climate in China.Quaternary International，2011，244:221-229.

[10]Li XH，Zhang Q，Xu CY.Suitability of the TRMM satellite rainfalls in driving a distributed hydrological model for water balance computations in Xinjiang catchment，Lake Poyang basin.Journal of Hydrology，2012，426/427:28-38.

[11]葉許春.近50年鄱陽湖水量變化機制與未來變化趨勢預(yù)估[學(xué)位論文].南京:中國科學(xué)院南京地理與湖泊研究所，2010.

[12]Zhang Q，Li LJ.Development and application of an integrated surface runoff and groundwater flow model for a catchment of Lake Taihu catchment，China.Quaternary International，2009，208(1/2):102-108.

[13]Zhang Q，Werner AD.Integrated surface-subsurface modeling of Fuxianhu Lake catchment，Southwest China.Water Resour Manage，2009，23:2189-2204.

[14]Harbaugh AW.MODFLOW-2005:The U.S.Geological survey modular ground-water model-the ground-water flow process.Geological Survey Techniques and Methods，2005.

[15]Kebede S，Travi Y，Alemayehu T et al.Water balance of Lake Tana and its sensitivity to fluctuations in rainfall，Blue Nile basin，Ethiopia.Journal of Hydrology，2006，316:233-247.

[16]DHI.MIKE 21 ＆ MIKE 3 FLOW MODEL FM Hydrodynamic Module.Water ＆ Environment ＆ Health，2007.

[17]Niemann SL，Jensen JH，Zyserman JA et al.Morphological modeling of a danish tidal inlet.Proceedings of ICCE，2006:1-14.

[18]Martinelli L，Zanuttigh B，Corbau C.Assessment of coastal flooding hazard along the Emilia Romagna littoral，IT.Coastal Engineering，2010，57:1042-1058.

[19]Smagorinsky J.General circulation experiment with the primitive equations.Monthly Weather Review，1963，91(3):99-164.

[20]Arnold JG，Allen PM，Muttiah R et al.Automated base flow separation and recession analysis techniques.Ground Water，1995，33(6):1010-1018.

[21]Doherty J.PEST:Model-independent parameter estimation.Watermark Numerical Computing，2005.

[22]Nash JE，Sutcliffe IV.River flow forcasting through conceptual models part 1—a discusstion of principles.Journal of Hydrology，1970，10:282-290.