基于水熱平衡的分布式水文模型研究與應(yīng)用

闞光遠(yuǎn)，李紀(jì)人，喻海軍，丁留謙，何曉燕，梁 珂

(1.流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100038；2.中國(guó)水利水電科學(xué)研究院，北京 100038；3.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038；4.水利部京津冀水安全保障重點(diǎn)實(shí)驗(yàn)室，北京 100038；5.北京中水科工程集團(tuán)有限公司，北京 100048)

1 研究背景

分布式水文模型在水利部“四預(yù)”建設(shè)中發(fā)揮著至關(guān)重要的作用[1-3]，“四預(yù)”涵蓋了預(yù)報(bào)、預(yù)警、預(yù)演、預(yù)案四大任務(wù)，急需分布式水文模型這一核心技術(shù)的有力支撐。實(shí)用分布式水文模型對(duì)研究區(qū)進(jìn)行概化，利用概念性元素描述產(chǎn)匯流過(guò)程，建立蓄滿/超滲/混合產(chǎn)流、分布式單位線、線性/非線性水庫(kù)等計(jì)算方法。在有資料區(qū)，這類方法具有精度高、資料種類要求低及建模簡(jiǎn)單的優(yōu)勢(shì)。但這類模型的物理基礎(chǔ)有待加強(qiáng)，預(yù)報(bào)結(jié)果的好壞很大程度上取決于參數(shù)率定過(guò)程中采用的徑流資料，在無(wú)徑流資料區(qū)(以下簡(jiǎn)稱無(wú)資料區(qū))[4-6]，模型參數(shù)無(wú)法通過(guò)徑流資料率定獲得[7-10]。此外，在人類活動(dòng)強(qiáng)度較大的地區(qū)，基于歷史資料的模擬和預(yù)報(bào)無(wú)法反映人類活動(dòng)的影響，形成了新的無(wú)資料區(qū)[11-16]，給模型應(yīng)用帶來(lái)很大困難[17-22]。當(dāng)前解決無(wú)資料區(qū)洪水預(yù)報(bào)[23-27](Predictions in Ungauged Basins，PUB)問(wèn)題的主要手段是相似流域和參數(shù)移用，這類方法尚存在不足有待改進(jìn)，如相似準(zhǔn)則不客觀且復(fù)雜多樣、依賴有資料流域的觀測(cè)數(shù)據(jù)以及統(tǒng)計(jì)一致性問(wèn)題。此外，分布式水文模型還面臨著數(shù)值求解精度和數(shù)值格式穩(wěn)定性等問(wèn)題。

PUB早已被國(guó)際水文科學(xué)協(xié)會(huì)確定為十年研究計(jì)劃，但由于問(wèn)題的復(fù)雜性和研究方法的局限性，它的徹底解決仍然任重而道遠(yuǎn)。許多水文學(xué)家認(rèn)為，解決該問(wèn)題的根本出路在于研發(fā)具有物理機(jī)制的分布式水文模型，在無(wú)資料區(qū)利用定量遙感反演流域特征與一部分模型參數(shù)，同時(shí)反演重點(diǎn)關(guān)注的水文變量，將反演的水文變量作為目標(biāo)，率定另一部分模型參數(shù)。這一思路在實(shí)際應(yīng)用中仍面臨挑戰(zhàn)，目前常用的遙感反演水文變量包括土壤濕度、蒸散發(fā)、總水儲(chǔ)量、徑流量等，這些數(shù)據(jù)存在尺度轉(zhuǎn)換和不確定性問(wèn)題，一旦脫離徑流資料，僅使用這類數(shù)據(jù)開(kāi)展參數(shù)率定，不能達(dá)到滿意的預(yù)報(bào)精度。

我們?cè)趯ふ腋玫臄?shù)據(jù)源時(shí)發(fā)現(xiàn)，以MODIS衛(wèi)星為代表的熱紅外遙感地表溫度產(chǎn)品日臻成熟，經(jīng)過(guò)質(zhì)量控制和校正，精度能夠滿足參數(shù)率定需求[28]。另一方面，在無(wú)資料區(qū)，水位資料的獲取相對(duì)容易，可由自計(jì)水位計(jì)等途徑獲得。對(duì)于入庫(kù)洪水預(yù)報(bào)等涉及大面積水體的應(yīng)用，還可以利用測(cè)高衛(wèi)星的庫(kù)區(qū)水面高程產(chǎn)品作為補(bǔ)充。隨著科技的持續(xù)進(jìn)步，SWOT等測(cè)高衛(wèi)星有望在未來(lái)提供空間分辨率百米級(jí)的水位產(chǎn)品用于監(jiān)測(cè)河道水面高程，能夠進(jìn)一步增強(qiáng)遙感水位產(chǎn)品的可用性。因此，地表溫度產(chǎn)品和水位資料在無(wú)資料區(qū)參數(shù)率定領(lǐng)域具有良好的應(yīng)用前景。但由于傳統(tǒng)分布式模型的一些制約，目前基于這兩類數(shù)據(jù)開(kāi)展模型參數(shù)率定還面臨一些技術(shù)難題，相關(guān)研究尚不多見(jiàn)。

傳統(tǒng)分布式水文模型要么僅考慮水量平衡過(guò)程，無(wú)法輸出熱量平衡過(guò)程對(duì)應(yīng)的通量和狀態(tài)變量(如：地表溫度)；要么雖然考慮了水熱平衡過(guò)程，但需將地表溫度作為模型輸入。另一方面，傳統(tǒng)水文模型著重于模擬流量過(guò)程，無(wú)法獲取水位、水深等水文要素的模擬結(jié)果。若需要獲取水位預(yù)報(bào)結(jié)果，通常利用水位流量關(guān)系反推水位過(guò)程，在無(wú)資料區(qū)無(wú)法應(yīng)用。因此，傳統(tǒng)分布式水文模型無(wú)法同時(shí)模擬流域地表溫度和水位的時(shí)空分布，難以利用地表溫度和水位數(shù)據(jù)提供的信息開(kāi)展參數(shù)率定，需要進(jìn)行改造。此外，傳統(tǒng)模型基于流量過(guò)程模擬誤差開(kāi)展全部產(chǎn)匯流參數(shù)的率定，未考慮狀態(tài)變量在參數(shù)率定中應(yīng)有的約束作用，沒(méi)有采用產(chǎn)流、匯流分層次率定的技術(shù)路線，導(dǎo)致徑流模擬誤差被分?jǐn)偟疆a(chǎn)流蒸散發(fā)和匯流參數(shù)中，加劇了模型參數(shù)的不確定性和異參同效問(wèn)題。

針對(duì)以上難題，本研究從水熱循環(huán)產(chǎn)匯流理論、模型構(gòu)建與求解、參數(shù)率定三個(gè)方面入手，改進(jìn)傳統(tǒng)分布式水文模型。通過(guò)研究水熱平衡強(qiáng)約束條件下陸氣界面—包氣帶—飽水帶水熱平衡動(dòng)力學(xué)過(guò)程，建立耦合熱量平衡的流域產(chǎn)匯流理論，揭示水量熱量平衡過(guò)程耦合機(jī)理和水熱循環(huán)通量量化關(guān)系。結(jié)合具有物理機(jī)制的產(chǎn)匯流計(jì)算方法，構(gòu)建系統(tǒng)完整的物理機(jī)制分布式水文—水動(dòng)力—熱動(dòng)力模型(Physically-based distributed HYdrological-HydrodYnamic-tHermodYnamic model，PHY模型)，實(shí)現(xiàn)流域地表溫度和水位的分布式數(shù)值模擬。通過(guò)明確模型結(jié)構(gòu)和模型參數(shù)的物理意義，利用野外觀測(cè)和定量遙感實(shí)現(xiàn)大部分模型參數(shù)的先驗(yàn)估計(jì)。將參數(shù)率定的目標(biāo)由徑流量轉(zhuǎn)換為地表溫度和水位，綜合利用遙感反演和臺(tái)站觀測(cè)資料，率定剩余的模型參數(shù)。在江西贛江流域開(kāi)展模型應(yīng)用和驗(yàn)證，結(jié)果表明，模型數(shù)值求解方法精度高、穩(wěn)定性好，實(shí)現(xiàn)了無(wú)資料區(qū)參數(shù)率定和水文模擬。研究成果不僅對(duì)豐富和發(fā)展流域水文模擬技術(shù)具有重要的理論意義，而且為無(wú)資料區(qū)分布式水文模型研發(fā)與應(yīng)用提供了新思路。

2 耦合水熱平衡動(dòng)力學(xué)過(guò)程的物理機(jī)制分布式水文模型

2.1 耦合熱量平衡的流域產(chǎn)匯流理論為了實(shí)現(xiàn)不依賴徑流資料的參數(shù)率定和水文模擬，將流域水量熱量平衡過(guò)程相耦合，基于水熱平衡原理，提出了耦合熱量平衡的流域產(chǎn)匯流理論，揭示了水量熱量平衡動(dòng)力學(xué)過(guò)程耦合機(jī)理，建立了水熱循環(huán)通量量化關(guān)系和流域地表溫度、蒸散發(fā)量的分布式數(shù)值模擬方法，為基于地表溫度的參數(shù)率定奠定了理論基礎(chǔ)。

2.1.1水量熱量平衡動(dòng)力學(xué)過(guò)程耦合機(jī)理 典型的陸氣界面水熱平衡動(dòng)力學(xué)過(guò)程見(jiàn)圖1。流域水量平衡過(guò)程(圖1中左半邊的水文過(guò)程)描述了涵蓋降雨—蒸發(fā)—入滲—產(chǎn)匯流的水文過(guò)程，水量收支可由以下的水量平衡方程描述：

圖1 陸氣界面水熱平衡動(dòng)力學(xué)過(guò)程Fig.1 Dynamic process of water-heat balance for interface between land and atmosphere

ΔSM=P-ET-R-WE

(1)

式中：SM為土壤濕度；P為降水量；ET為蒸散發(fā)量；R為徑流量；WE為包氣帶與飽水帶間的水量交換量。其中R可由產(chǎn)流模型計(jì)算獲得；WE可由包氣帶土壤水動(dòng)力學(xué)模型、飽水帶地下水動(dòng)力學(xué)模型結(jié)合達(dá)西定律計(jì)算獲得；ET可由定量遙感技術(shù)結(jié)合熱動(dòng)力學(xué)方法計(jì)算獲得。

流域熱量平衡過(guò)程(圖1中右半邊的熱動(dòng)力學(xué)過(guò)程)描述了涵蓋太陽(yáng)輻射—陸面反射—土壤吸放熱的熱動(dòng)力學(xué)過(guò)程。在計(jì)算網(wǎng)格內(nèi)利用定量遙感技術(shù)反演出坡地網(wǎng)格的植被覆蓋率，將計(jì)算網(wǎng)格劃分為裸土與植被兩個(gè)部分。太陽(yáng)凈輻射量到達(dá)地面后，帶來(lái)的熱量分別被裸土和植被吸收，熱量平衡過(guò)程驅(qū)動(dòng)蒸發(fā)的產(chǎn)生，蒸發(fā)影響水量平衡過(guò)程，水量增減影響供水條件反作用于蒸發(fā)，形成互饋耦合關(guān)系，最終達(dá)到水量熱量平衡狀態(tài)，這時(shí)的地表溫度稱為典型平衡態(tài)溫度(Representative Equilibrium Temperature，RET)。描述熱量收支的熱量平衡方程式可表達(dá)為：

ΔW=Rn-G-(Hs+Hc)-(LEs+LEc)

(2)

式中：W為熱量蓄量；Rn為太陽(yáng)凈輻射量；G為土壤熱通量；H和LE分別為感熱和潛熱通量，下標(biāo)s和c分別表示裸土和植被冠層。針對(duì)植被與裸土分別建立熱量平衡方程描述熱量收支過(guò)程。利用定量遙感方法計(jì)算土壤熱通量、感熱通量及潛熱通量，蒸散發(fā)過(guò)程還受到氣象條件和下墊面供水條件的制約，這些條件的調(diào)控作用通過(guò)熱量平衡方程中的空氣動(dòng)力學(xué)阻抗、裸土阻抗及植被冠層阻抗等變量來(lái)反映。

水熱平衡過(guò)程的耦合機(jī)理分析如下：針對(duì)每個(gè)計(jì)算時(shí)刻、流域內(nèi)每個(gè)計(jì)算網(wǎng)格，可列出水量和熱量平衡方程。水量平衡方程中的未知量為土壤濕度和蒸散發(fā)量，熱量平衡方程中的未知量為感熱通量和潛熱通量，而二者均能表達(dá)為土壤濕度和蒸散發(fā)量的函數(shù)。蒸散發(fā)量是溝通水量和熱量平衡方程的關(guān)鍵變量，借助于定量遙感技術(shù)，蒸散發(fā)量能夠與潛熱通量、地表溫度建立起緊密的物理和數(shù)值聯(lián)系。因此，可以利用蒸散發(fā)量作為水量和熱量平衡方程的耦合紐帶，聯(lián)立水量和熱量平衡方程，形成水熱平衡方程組，而待求解的未知量則轉(zhuǎn)化為土壤濕度和地表溫度。通過(guò)以上分析，揭示了水熱平衡動(dòng)力學(xué)過(guò)程耦合機(jī)理，為在陸氣界面上構(gòu)建與土壤含水量相耦合的分布式熱動(dòng)力學(xué)方法計(jì)算地表溫度和蒸散發(fā)量奠定了理論基礎(chǔ)。

2.1.2 水熱循環(huán)通量量化關(guān)系 蒸散發(fā)量作為耦合水量、熱量平衡方程的紐帶變量，其參數(shù)化方案至關(guān)重要?；跓釀?dòng)力學(xué)和定量遙感方法，建立了蒸散發(fā)量與潛熱通量、地表溫度間的數(shù)值聯(lián)系，構(gòu)建的計(jì)算公式如下：

(3)

式中：LE為潛熱通量；λ為汽化潛熱；ρw為水的密度；ρa(bǔ)為空氣密度；cp為濕空氣比熱；γ為干濕表常數(shù)；fv為植被覆蓋率；ra、rabs、rc、rs分別為植被空氣動(dòng)力學(xué)阻抗、裸土空氣動(dòng)力學(xué)阻抗、冠層阻抗、土壤阻抗，阻抗的參數(shù)化方案均與土壤濕度有關(guān)；ea為蒸汽壓力，與溫度有關(guān)；e*為飽和蒸汽壓力，可由RET計(jì)算得到：

(4)

式中RETc為攝氏溫標(biāo)下的RET值。由以上參數(shù)化方案可實(shí)現(xiàn)蒸散發(fā)量、潛熱通量、地表溫度三者間的相互轉(zhuǎn)化，為水熱平衡方程組的聯(lián)立求解奠定了量化基礎(chǔ)。

2.1.3 地表溫度和蒸散發(fā)量分布式數(shù)值模擬 將熱量平衡方程與水量平衡方程聯(lián)立形成如下的水熱平衡方程組：

(5)

圖2 水熱平衡方程組求解流程圖Fig.2 Solution procedure for water-heat balance equations

具體求解流程為：利用時(shí)段初(t-1時(shí)刻)土壤濕度作為迭代計(jì)算的初值，代入熱量平衡方程，進(jìn)行熱量平衡牛頓—拉夫森迭代計(jì)算，獲取熱量平衡狀態(tài)下的地表溫度和蒸散發(fā)量，此為內(nèi)層熱量平衡迭代；基于水量平衡原理，以上時(shí)刻土壤濕度作為初值，扣除熱量平衡方程求得的蒸散發(fā)量，利用Richards方程和熱傳導(dǎo)方程求解土壤的產(chǎn)流、分水源及土壤濕度、溫度和地下水位的再分布過(guò)程，從而獲取時(shí)段末(t時(shí)刻)的土壤濕度，此為外層水量平衡迭代，如果尚未達(dá)到收斂條件，則將更新后的時(shí)段末土壤濕度回代入熱量平衡方程再次進(jìn)行內(nèi)層迭代，以此類推。隨著內(nèi)外層迭代的交替執(zhí)行，土壤濕度逐漸趨于穩(wěn)定收斂，同時(shí)可獲得平衡態(tài)下的地表溫度等關(guān)鍵狀態(tài)變量。算法通過(guò)水量平衡和熱量平衡過(guò)程的耦合模擬，計(jì)算各時(shí)段土柱的蒸散發(fā)量、產(chǎn)流量、土壤濕度、地下水位、土柱溫度分布和地表溫度。

2.2 水熱平衡分布式模型構(gòu)建及數(shù)值求解方法基于耦合熱量平衡的流域產(chǎn)匯流理論，構(gòu)建了PHY模型，模型包括地表溫度計(jì)算模塊、土壤水熱運(yùn)移計(jì)算模塊、地表水運(yùn)動(dòng)計(jì)算模塊、地下水運(yùn)動(dòng)計(jì)算模塊，實(shí)現(xiàn)了流域水熱通量和狀態(tài)變量的分布式數(shù)值模擬。模型模擬的水文過(guò)程包括：冠層截留、地表填洼、水面蒸發(fā)、土壤蒸發(fā)、植被散發(fā)、潛水蒸發(fā)、降雨入滲、包氣帶補(bǔ)給潛水、河道滲漏、垂向一維包氣帶土壤水運(yùn)動(dòng)、二維飽水帶潛水運(yùn)動(dòng)、一二維耦合地表水運(yùn)動(dòng)、陸氣界面及垂向一維包氣帶熱量運(yùn)移等。模擬的水文過(guò)程概化圖見(jiàn)圖3。

圖3 水循環(huán)過(guò)程概化圖Fig.3 Sketch of the water cycle

2.2.1 土壤水運(yùn)動(dòng)數(shù)值模擬 利用垂向一維Richards方程構(gòu)建能夠考慮土壤各向異性、側(cè)向壤中水徑流、變動(dòng)地下水位的土壤水運(yùn)動(dòng)數(shù)值模型，用于描述水量平衡約束條件下的包氣帶—飽水帶水分運(yùn)動(dòng)過(guò)程，實(shí)現(xiàn)給定蒸散發(fā)量條件下的產(chǎn)流、分水源、地下水位求解。為了使模型能夠同時(shí)適用于飽和、非飽和區(qū)域的數(shù)值模擬，采用水頭型控制方程并寫(xiě)為如下形式：

(6)

式中：C為容水度；h為壓力水頭；K為水力傳導(dǎo)度；Sc為源匯項(xiàng)；t為時(shí)間；z為位置坐標(biāo)。模型求解涉及的土壤水力特性由van Genuchten-Mualem模型描述。采用有限體積法對(duì)控制方程進(jìn)行離散，并采用改進(jìn)的Picard法和追趕法求解形成的非線性方程組。

降雨入滲、產(chǎn)流、分水源、土壤水運(yùn)動(dòng)、地下水位的計(jì)算流程如下：首先進(jìn)行冠層截留和地表填洼的計(jì)算。然后由熱量平衡方程提供的時(shí)段實(shí)際蒸散發(fā)量從土柱頂端逐節(jié)點(diǎn)向下按照一定比例添加到每個(gè)節(jié)點(diǎn)的源匯項(xiàng)中，依次扣減水量，如果整個(gè)土柱蒸干，則剩余的蒸散發(fā)量由潛水蒸發(fā)提供。土層頂端的產(chǎn)流計(jì)算由下滲能力進(jìn)行控制，下滲能力基于達(dá)西定律求得。雨強(qiáng)小于下滲能力時(shí)，下滲率等于雨強(qiáng)，土柱上邊界條件采用定通量。當(dāng)雨強(qiáng)大于下滲能力時(shí)將在地表產(chǎn)生積水，上邊界條件切換為定水頭。進(jìn)行壤中水徑流計(jì)算時(shí)，土層內(nèi)部設(shè)置一定厚度的壤中水徑流敏感土層，壤中水徑流將由該敏感土層包含的節(jié)點(diǎn)產(chǎn)出，壤中水徑流量取決于水力坡度、敏感土層含水量、敏感土層厚度、側(cè)向?qū)实纫蛩?，基于達(dá)西定律計(jì)算。土柱下邊界為潛水位，根據(jù)潛水是否枯竭，在零通量和零負(fù)壓水頭邊界條件間自適應(yīng)切換。當(dāng)包氣帶補(bǔ)給潛水時(shí)，產(chǎn)生的通量即為地下徑流補(bǔ)給量，反之，為潛水蒸發(fā)。壤中水徑流敏感土層厚度、各向異性系數(shù)通過(guò)參數(shù)率定獲得。

2.2.2 土壤熱量運(yùn)移數(shù)值模擬 利用與Richards方程耦合的對(duì)流擴(kuò)散方程描述土壤熱量的垂向運(yùn)移過(guò)程，控制方程的形式如下：

(7)

式中：Ts為土壤溫度；v為土壤水滲流速度；D為擴(kuò)散系數(shù)。土柱頂端的邊界條件采用熱量平衡方程求解的地表典型平衡態(tài)溫度，對(duì)流速度采用Richards方程求解的滲流速度。由于地下水溫度通常較為穩(wěn)定，故下邊界條件采用固定的地下水恒溫。對(duì)流擴(kuò)散方程利用有限體積法進(jìn)行全隱式離散，采用追趕法求解。熱量平衡方程求解地表典型平衡態(tài)溫度時(shí)需要表層土壤濕度這一狀態(tài)變量，表層土壤濕度對(duì)應(yīng)的土層稱為熱敏感土層，熱量平衡計(jì)算中所需的土壤濕度采用熱敏感土層所含節(jié)點(diǎn)土壤濕度的平均值。熱敏感土層的厚度通過(guò)參數(shù)率定獲得。

2.2.3 地下水運(yùn)動(dòng)數(shù)值模擬 由于地下水運(yùn)動(dòng)速度通常較為緩慢，對(duì)于場(chǎng)次洪水模擬，考慮潛水層的地下水運(yùn)動(dòng)即可滿足精度要求。PHY模型的地下水模塊考慮了潛水運(yùn)動(dòng)及其與包氣帶間的交互作用。采用的二維潛水運(yùn)動(dòng)控制方程為：

(8)

式中：h為水頭；q為源匯項(xiàng)；St為貯水系數(shù)；T為導(dǎo)水系數(shù)。采用有限體積法對(duì)控制方程進(jìn)行離散，時(shí)間格式采用全隱式格式，利用共軛梯度法求解時(shí)空離散后獲得的大型稀疏線性方程組。為提升迭代求解的計(jì)算效率，采用高斯賽德?tīng)栴A(yù)處理技術(shù)降低系數(shù)矩陣的條件數(shù)以加速收斂。數(shù)值求解的上邊界與包氣帶聯(lián)通，包氣帶與飽水帶的水量交換基于潛水位和達(dá)西定律計(jì)算獲得。求解區(qū)域的下邊界為地下水第一個(gè)不透水層，采用零通量邊界條件。

2.2.4 地表水運(yùn)動(dòng)數(shù)值模擬

(1)一維河網(wǎng)模型 PHY模型將天然河道概化為矩形斷面，河道演進(jìn)采用一維水動(dòng)力模型，求解完整形式的動(dòng)力波，控制方程采用圣維南方程組，其向量形式為：

(9)

式中：U為待求解未知量；F為通量；S為源項(xiàng)。采用有限體積法對(duì)控制方程進(jìn)行二階精度MUSCL空間離散，時(shí)間離散采用二階TVD龍格—庫(kù)塔顯格式，網(wǎng)格間數(shù)值通量由HLL近似黎曼求解器計(jì)算獲得，源項(xiàng)中的底坡項(xiàng)基于靜水重構(gòu)法和中心差分離散計(jì)算獲得，摩阻項(xiàng)采用全隱式處理，利用牛頓—拉夫森迭代計(jì)算，以獲取更好的數(shù)值穩(wěn)定性。模型計(jì)算步長(zhǎng)采用基于CFL條件的自適應(yīng)時(shí)間步長(zhǎng)，保證計(jì)算的穩(wěn)定性。

由于PHY模型應(yīng)用于天然流域的水文模擬，故所涉及的河網(wǎng)通常為樹(shù)狀結(jié)構(gòu)，即多個(gè)上游河段匯入單一下游河段，因此汊點(diǎn)的處理較為方便。模型將汊點(diǎn)處理為多個(gè)上游河段匯入的一條單一河道，數(shù)值計(jì)算時(shí)需要獲取通過(guò)汊點(diǎn)河段上斷面的總的質(zhì)量通量和動(dòng)量通量，以及汊點(diǎn)河段的平均底坡項(xiàng)。通過(guò)上斷面的總通量由多個(gè)黎曼求解器計(jì)算獲得，平均底坡項(xiàng)的計(jì)算基于上游各河段平均河底高程和平均水位、汊點(diǎn)河段下斷面河底高程和水位以及改進(jìn)的靜水重構(gòu)法和中心差分離散求得。

(2)二維坡面模型 PHY模型采用二維地表水動(dòng)力模型處理坡面匯流計(jì)算，控制方程采用完整的二維淺水方程，其向量形式為

(10)

式中：U為待求未知量；E和G為x和y方向通量；S為源項(xiàng)。采用有限體積法對(duì)控制方程進(jìn)行二階精度MUSCL空間離散，時(shí)間離散采用二階TVD龍格—庫(kù)塔顯格式，網(wǎng)格間數(shù)值通量由HLLC近似黎曼求解器計(jì)算獲得，源項(xiàng)中的底坡項(xiàng)基于靜水重構(gòu)法和中心差分離散計(jì)算獲得，摩阻項(xiàng)采用全隱式處理，利用牛頓—拉夫森迭代計(jì)算，以獲取更好的數(shù)值穩(wěn)定性。模型計(jì)算步長(zhǎng)采用基于CFL條件的自適應(yīng)時(shí)間步長(zhǎng)，保證計(jì)算的穩(wěn)定性。

(3)一二維耦合及地表水與包氣帶、飽水帶水量交互 一維河道與二維坡面間的水量交互采用側(cè)向連接，模型耦合基于河底高程、坡面高程、河道水位、坡面水位，利用HLL近似黎曼求解器計(jì)算坡面與河道間的側(cè)向質(zhì)量通量，交互的水量添加入水動(dòng)力學(xué)模型的源項(xiàng)中，保證水量平衡和計(jì)算的穩(wěn)定性。

一二維模型與包氣帶、飽水帶模型的水量交互采用側(cè)向連接，模型耦合基于河底高程、坡面高程、地下水不透水基巖高程、河道水位、坡面水位、潛水位，利用達(dá)西定律計(jì)算交互的水量，并將水量添加到動(dòng)力學(xué)模型的源項(xiàng)中，實(shí)現(xiàn)水量平衡和數(shù)值計(jì)算的穩(wěn)定。

2.3 無(wú)徑流資料區(qū)模型參數(shù)率定PHY模型由蒸散發(fā)、產(chǎn)流分水源、土壤水熱運(yùn)移及匯流等計(jì)算模塊組成，各模塊的模型參數(shù)分為流域特征參數(shù)和敏感參數(shù)兩類。流域特征參數(shù)(如地形、河寬、河槽深度等)通過(guò)遙感協(xié)同定量反演、先驗(yàn)估計(jì)結(jié)合統(tǒng)計(jì)規(guī)律綜合確定；模型的敏感參數(shù)(如敏感土層厚度、土壤各向異性系數(shù)、坡地和河道糙率等)通過(guò)參數(shù)率定獲取。將質(zhì)量控制后的MODIS地表溫度和控制斷面實(shí)測(cè)水位(或測(cè)高衛(wèi)星反演的水位)作為參數(shù)率定的目標(biāo)，分別與模型輸出的RET和水位進(jìn)行比對(duì)，基于地表溫度和水位模擬誤差最小化原則，針對(duì)產(chǎn)流蒸散發(fā)過(guò)程和匯流過(guò)程，分兩個(gè)層次先后調(diào)整相應(yīng)模塊的模型參數(shù)，獲取流域尺度的最優(yōu)等效參數(shù)，實(shí)現(xiàn)不依賴徑流資料的參數(shù)率定。

3 PHY模型的應(yīng)用

3.1 研究區(qū)域與數(shù)據(jù)

3.1.1 研究區(qū)域 本研究的實(shí)驗(yàn)流域?yàn)榻魇≮M江流域及流域內(nèi)部嵌套的4個(gè)子流域(研究流域DEM和計(jì)算單元網(wǎng)格類型見(jiàn)圖4。贛江是長(zhǎng)江的第七大支流，同時(shí)是江西省最大的河流，位于我國(guó)東南部。贛江流域控制站外洲水文站以上集水面積為81 258 km2，流域內(nèi)海拔高度11～1997 m，地形較為復(fù)雜，以山地丘陵為主，山地占流域面積的43.9%，低丘(海拔200 m以下)崗地占31.5%。贛江流域?qū)儆趤啛釒駶?rùn)季風(fēng)氣候，2003—2009年平均降水量約為1500 mm。年內(nèi)降水多集中在4—6月份，暴雨洪水產(chǎn)生的最大流量多出現(xiàn)在5—7月。四個(gè)嵌套子流域的出口水文站分別為新田、賽塘、白沙、棟背，此外，流域內(nèi)包含萬(wàn)安和峽江兩個(gè)大(一)型水庫(kù)。

圖4 贛江流域圖Fig.4 Location map of the Ganjiang watershed

3.1.2 數(shù)據(jù)資料 實(shí)驗(yàn)流域數(shù)據(jù)的獲取主要通過(guò)以下途徑：用于模型構(gòu)建和率定的數(shù)據(jù)來(lái)自互聯(lián)網(wǎng)下載的遙感反演數(shù)據(jù)集和水位站點(diǎn)觀測(cè)數(shù)據(jù)，用于模型驗(yàn)證的數(shù)據(jù)來(lái)自流域觀測(cè)站點(diǎn)、遙感反演數(shù)據(jù)集和陸面數(shù)據(jù)同化系統(tǒng)資料集。為了降低地表溫度產(chǎn)品的不確定性，開(kāi)展數(shù)據(jù)質(zhì)量控制與校正，包括極值檢查、一致性檢查、基于實(shí)測(cè)地表溫度修正、多遙感產(chǎn)品(MODIS、AVHRR、AATSR、ETM+和ASTER)交叉比較修正、云層遮蔽像元插值與重建。模型率定時(shí)未采用徑流資料，開(kāi)展無(wú)徑流資料地區(qū)參數(shù)率定，徑流資料用于驗(yàn)證模型的模擬精度。

3.2 徑流模擬

3.2.1 徑流過(guò)程線 圖5展示的是五個(gè)徑流站點(diǎn)，包括外洲(流域出口站)、棟背、新田、賽塘、白沙的實(shí)測(cè)和模擬小時(shí)徑流過(guò)程線。由圖5中結(jié)果可見(jiàn)，降雨—徑流關(guān)系合理、良好，模擬的過(guò)程線光滑，與實(shí)測(cè)過(guò)程線擬合良好。外洲和棟背站上游分別有峽江和萬(wàn)安兩座大(一)型水庫(kù)，受水庫(kù)放水影響明顯，因此，這兩站的洪水模擬采用區(qū)間產(chǎn)匯流與上游來(lái)水演進(jìn)相結(jié)合的計(jì)算方法求得。通過(guò)過(guò)程線觀察可知，模型模擬效果良好，尤其是對(duì)于棟背站，該站上游萬(wàn)安水庫(kù)調(diào)度頻繁，閘門(mén)啟閉對(duì)下游影響明顯，徑流經(jīng)常處于波動(dòng)狀態(tài)，模型能夠很好的模擬這種波動(dòng)過(guò)程。對(duì)于新田、賽塘、白沙三站，徑流模擬效果也是令人滿意的。

圖5 實(shí)測(cè)與模擬的徑流過(guò)程Fig.5 Observed and simulated hydrographs

3.2.2 徑流模擬散點(diǎn)圖 圖6為五個(gè)徑流站點(diǎn)的實(shí)測(cè)與模擬徑流散點(diǎn)圖。由圖6中結(jié)果可知，大部分站點(diǎn)的模擬結(jié)果均較為滿意，R2值較高。僅新田站R2略低，為0.5左右。通過(guò)觀察新田站實(shí)測(cè)和模擬徑流過(guò)程線可知，預(yù)熱期和洪峰的模擬誤差較大。這幾處徑流過(guò)程均為復(fù)式洪峰，降雨也為多峰，因此導(dǎo)致較大的非線性，且前后洪水相互疊加影響，導(dǎo)致模擬難度較大。

圖6 實(shí)測(cè)與模擬的徑流散點(diǎn)圖Fig.6 Scatter plot of observed and simulated hydrographs

3.3 地表溫度模擬圖7展示的是GLDAS數(shù)據(jù)集面平均地表溫度與模型模擬的面平均地表溫度過(guò)程。從圖7結(jié)果可見(jiàn)，面平均地表溫度在整體過(guò)程上模擬結(jié)果較好，模擬的地表溫度與GLDAS再分析數(shù)據(jù)產(chǎn)品結(jié)果比較一致，地表溫度大部分?jǐn)?shù)值的擬合結(jié)果很好。但對(duì)于高溫和低溫的極值部分，地表溫度的模擬效果還有進(jìn)一步提升的空間。

圖7 GLDAS與模擬的面平均地表溫度過(guò)程Fig.7 GLDAS and simulated areal mean land surface temperature

3.4 表層土壤濕度模擬圖8展示的是SMOS衛(wèi)星反演土壤濕度與模型模擬的表層土壤濕度過(guò)程。表層土壤濕度為地面以下5 cm厚度土層的模擬結(jié)果，且為流域面平均值。模型輸出的土壤濕度為小時(shí)尺度，模型的結(jié)果經(jīng)過(guò)取日均值繪制在圖8中，SMOS衛(wèi)星反演數(shù)據(jù)為間隔數(shù)日的日尺度產(chǎn)品。由圖8結(jié)果可知，土壤濕度的模擬效果較好，從趨勢(shì)上和具體點(diǎn)據(jù)上均較為滿意。

圖8 SMOS與模擬的面平均表層土壤濕度過(guò)程Fig.8 SMOS and simulated areal mean upper layer soil moisture

3.5 實(shí)際蒸散發(fā)量模擬圖9展示的是MODIS衛(wèi)星反演的與模型模擬的實(shí)際蒸散發(fā)過(guò)程，是流域面平均蒸散發(fā)量的結(jié)果。模型輸出的實(shí)際蒸散發(fā)為小時(shí)尺度，模型的結(jié)果經(jīng)過(guò)取日累積值繪制在圖9中。從圖9結(jié)果可知，實(shí)際蒸散發(fā)的模擬結(jié)果在總體趨勢(shì)上較好，過(guò)程上較為一致。遙感蒸散產(chǎn)品的結(jié)果與模型模擬結(jié)果相比存在低估現(xiàn)象。經(jīng)過(guò)文獻(xiàn)調(diào)研和分析發(fā)現(xiàn)，MODIS蒸散發(fā)產(chǎn)品在贛江流域存在低估現(xiàn)象，這一結(jié)論與本研究的結(jié)果較為一致。遙感產(chǎn)品的蒸散發(fā)結(jié)果較為平穩(wěn)，而模型輸出的結(jié)果具有波動(dòng)性。

圖9 MODIS與模擬的面平均實(shí)際蒸散發(fā)過(guò)程Fig.9 MODIS and simulated areal mean actual evapotranspiration

4 結(jié)論

本研究以強(qiáng)化分布式水文模型物理機(jī)制、改善分布式水文模型無(wú)徑流資料區(qū)參數(shù)率定效果及提高分布式水文模型精度和穩(wěn)定性為總目標(biāo)，選擇物理機(jī)制分布式水文模型為研究對(duì)象，開(kāi)展了以下研究：深入分析了流域水量熱量循環(huán)動(dòng)力學(xué)過(guò)程耦合機(jī)理，建立了水量熱量通量量化關(guān)系，提出了耦合熱量平衡的流域產(chǎn)匯流理論和流域地表溫度、蒸散發(fā)量的分布式計(jì)算方法；基于水熱平衡原理，構(gòu)建具有強(qiáng)物理機(jī)制的分布式水文—水動(dòng)力—熱動(dòng)力模型PHY及統(tǒng)一有限體積數(shù)值求解方法；利用臺(tái)站觀測(cè)、遙感反演、再分析等途徑獲取的多源數(shù)據(jù)，研究PHY模型參數(shù)先驗(yàn)估計(jì)方法和不依賴徑流資料的模型參數(shù)率定方法；將PHY模型在典型流域進(jìn)行了應(yīng)用與驗(yàn)證。

經(jīng)過(guò)數(shù)值模擬與結(jié)果分析，得到以下結(jié)論：

(1)研究探明了水熱平衡耦合機(jī)理，提出了缺徑流資料區(qū)水文模擬和參數(shù)率定的新理論和新方法。構(gòu)建的PHY模型綜合考慮了水量熱量耦合平衡過(guò)程，以蒸散發(fā)量作為橋梁實(shí)現(xiàn)了水熱平衡的深度耦合，完成了缺徑流資料區(qū)分布式水文模擬，模型具備地表溫度和水位變量的分布式數(shù)值模擬功能，利用地表溫度和水位變量作為參數(shù)率定的目標(biāo)，完成了不依賴徑流資料的模型參數(shù)率定。

(2)針對(duì)PHY模型中的水文、水熱、水動(dòng)力計(jì)算模塊，提出了統(tǒng)一有限體積數(shù)值求解方法，顯著提升了數(shù)值格式和數(shù)值方法的穩(wěn)定性、數(shù)值精度、可靠性和計(jì)算效率。對(duì)于各類初始和邊界條件，開(kāi)展了數(shù)值模擬測(cè)試，結(jié)果表明，求解方法具有很好的數(shù)值穩(wěn)定性和收斂性，能夠可靠高效的獲取合理的模擬結(jié)果。

(3)通過(guò)將模型應(yīng)用于江西贛江流域洪水預(yù)報(bào)實(shí)踐中，對(duì)模型模擬精度進(jìn)行了驗(yàn)證。研究結(jié)果表明，通過(guò)引入水熱平衡耦合計(jì)算方法、有限體積數(shù)值求解方法和多源信息，開(kāi)展參數(shù)率定和數(shù)值模擬，能夠有效提升模型的模擬精度、可靠度和數(shù)值穩(wěn)定性，實(shí)現(xiàn)了不依賴徑流資料的參數(shù)率定和水文模擬。成果在產(chǎn)匯流理論、水文遙感、最優(yōu)化方法等方面具有理論和實(shí)用價(jià)值，同時(shí)也為無(wú)資料地區(qū)模型參數(shù)率定提供了新思路。