闞光遠，李紀人，喻海軍，丁留謙，何曉燕，梁 珂

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

1 研究背景

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

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

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

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

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

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

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

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

圖1 陸氣界面水熱平衡動力學(xué)過程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)流模型計算獲得；WE可由包氣帶土壤水動力學(xué)模型、飽水帶地下水動力學(xué)模型結(jié)合達西定律計算獲得；ET可由定量遙感技術(shù)結(jié)合熱動力學(xué)方法計算獲得。

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

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

(2)

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

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

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

(3)

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

(4)

式中RETc為攝氏溫標(biāo)下的RET值。由以上參數(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

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

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

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

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

(6)

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

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

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

(7)

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

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

(8)

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

2.2.4 地表水運動數(shù)值模擬

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

(9)

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

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

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

(10)

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

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

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

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

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

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

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

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

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

3.2 徑流模擬

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

圖5 實測與模擬的徑流過程Fig.5 Observed and simulated hydrographs

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

圖6 實測與模擬的徑流散點圖Fig.6 Scatter plot of observed and simulated hydrographs

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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