第2章 分布式水循環模擬模型EasyDHM

2.1 水文模擬

2.1.1 產流理論

本書提出的EasyDHM產流模型是在Wetspa模型（Water and Energy Transfer Between Soil，Plants and Atmosphere）、SWAT模型（Soil and Water Assessment Tool）、新安江模型等產流模型的基礎上，通過集成創新而提出的一種產流模型，其在不同地區、不同水文地質條件下均具有通用性，其產流過程概化見圖2-1。EasyDHM模型在垂向上劃分了4層，即植被冠層、地表層、土壤層和地下水含水層。首先降雨進入植被冠層，發生冠層截留，穿過冠層的水分則會進入地表層，一部分超滲的水量會發生地表填洼，一部分形成地表徑流，其他水分則通過入滲進入到土壤，土壤水側向流出形成壤中流，土壤水下滲進入地下水含水層，形成地下徑流。地表徑流、壤中流和地下徑流之和即為總徑流。植被冠層、地表填洼、土壤水和地下水均會產生蒸散發。

EasyDHM產流模型除了支持植被冠層、地表過程、土壤水過程及地下水過程模擬外，還可實現寒區水文模擬，如積雪/融雪、凍土/凍融模擬等。

2.1.1.1 冠層過程模擬

一般，冠層截留損失在降水初期相對較高，隨后會逐漸趨于零。在EasyDHM產流模型中，認為冠層截留達到最大存儲量時，截留率才會下降。如果在某一時段內，總降水量比冠層截留蓄水量大，則降水先會充滿所有冠層截留容量，剩余降水則會繼續下降；否則全部降水均認為被冠層截留，剩余的冠層截留量則在之后的時段里從降水中除去，其表達式為：

式中：Ii（t）為t時段內i單元的冠層截留損失，mm；Ii，0為單元截留蓄水能力，mm；SIi（t-1）為t-1時間段的單元截留蓄水量，mm；Pi（t）為單元降水量，mm。

在某一單元中，截留蓄水量平衡的計算公式為：

式中：SIi（t-1）和SIi（t）分別為t-1時間段和t時間段i單元的截留蓄水量，mm；EIi（t）為i單元截留量中的蒸發部分，mm。

當截留量為0或是在降水過程中，EIi（t）=0；當Pi（t）=0且EP&gt；SIi（t-1）&gt；0時，EIi（t）=SIi（t-1），其中EP為潛在蒸發量，mm；其他條件下，EIi（t）=EP。

截留蓄水能力是有關葉面積指數和植被種類的函數。由于截留蓄水能力隨時間不斷變化，為方便模型計算，建立了一個簡單的正弦變化曲線?；陂L期測量值的統計分析，其經驗方程為：

式中：Ii，min為單元網格i中的最小截留蓄水能力，mm；Ii，max為最大截留蓄水能力，mm；d為日序數；指數b控制正弦變化曲線的形狀，可根據局部條件進行調整。

EasyDHM模型中假定一天中每小時截留蓄水能力恒定，因此截留蓄水能力為日函數。

2.1.1.2 地表過程模擬

扣除截留后的降雨到達地表后，部分降雨將發生地表入滲、填洼、地表徑流及填洼蒸發等過程，這些統稱為地表過程。

在EasyDHM中，計算地表入滲量和扣除入滲量后的剩余降雨量的方法有兩種：一種是與土壤含水量有關的修改系數法；另一種是SWAT模型中采用的SCS曲線法。修改系數法的思想是采用超滲產流，即在特定降雨過程中，當降雨強度超過了土壤下滲能力時產生的地面徑流，用該方法計算剩余降雨量和地表入滲量公式為：

式中：PEi（t）為一定時間間隔中i網格上超滲降水量，mm；Pi（t）為t時段的降水量，mm；Ii（t）為植被冠層截留損失，mm；θi（t）為網格內t時間的土壤含水量，m3/m3；θs，i為土壤孔隙度，即飽和含水量；a為與降水強度有關的指數；Ci為網格潛在降水超滲系數或潛在徑流系數；Fi（t）為單元網格的入滲量，mm。

SCS曲線法又稱曲線數值法（curve number method），是美國農業部水土保持局（Soil Conservation Service，SCS）于20世紀50年代年開發研制的，其幾乎可以在所有的土地利用類型和土壤類型的區域使用。該方法計算剩余降雨量和地表入滲量的公式為：

式中：St（t）為滯留參數。

St（t）的計算公式為：

式中：CNi（t）為SCS曲線參數。

地表下滲后，剩余降雨量部分進行填洼，部分成為地表徑流參與匯流。某一時刻地表填洼量與地表填洼能力、地表填洼水儲量及累積剩余降雨量有關，其計算經驗公式為：

式中：SDi（t）為地表填洼水儲量，mm；SDi，0為地表填洼能力，mm；PCi（t）為累積剩余降雨量，mm，即t時刻時計算單元i的累積剩余降雨量。EasyDHM產流模型僅根據式（2-8）模擬初始時刻的地表填洼水儲量，其后按更新填洼需水量公式（2-9）計算。t時刻的地表填洼量為：

進行地表填洼后的累積剩余降雨量則為：

地表徑流RSi（t）為計算單元i在計算時刻t的剩余降雨量減去地表填洼量，其計算公式為：

根據EasyDHM產流模型，認為填洼水量均通過蒸發逐漸消逝。填洼水的蒸發量計算式為：

式中：PE為潛在蒸發量，mm；EIi（t）為截留蒸發量，mm。

當發生填洼蒸發后，即需要更新填洼水儲量，其計算公式為：

2.1.1.3 土壤水過程模擬

EasyDHM產流模型中，土壤水的模擬采用分層土壤模型，逐層計算土壤層的下滲和壤中流。分層土壤模型將土壤分為若干層（層數可自行設定），從上到下，假設土壤水在每一層均勻分布，當某一層的土壤含水量θi達到田間持水量θfc，i時形成重力水，而其下層土壤未飽和，該層土壤水往下層移動，當最下層土壤水蓄滿后即向土壤層下的包氣帶滲漏，繼而對地下含水層進行補給。

分層土壤模型將每個計算單元進行了垂向分層，對于任一單元的第i層土壤，其每天可下滲水量SW0，i為：

式中：θi為第i層土壤的水含量，mm；θfc，i為第i層土壤的田間持水率，mm。

另外，若第i層土壤溫度較低，低于某一臨界溫度時，則視為凍土，凍土層土壤水不發生下滲。實際上，并不是所有可下滲水量均進入下一層土壤，第i層土壤的實際下滲量Wp，i為：

式中：Δt為計算時間步長；TTi為下滲持續時間。

TTi計算式為：

式中：θs，i為第i層土壤的飽和含水率，mm；Ki為該層土壤滲透系數，mm/s。

由于對土壤進行了分層處理，則可對土壤垂向土層分布多樣化的情況進行分析，如某層土壤為不透水層，則水體下滲到此層即停止。

當某層水蓄滿且該層存在一定的坡度時，則會發生橫向流動，即形成壤中流，否則其上層土壤水發生橫向流動，如此逐層上推，直至最上層土壤發生蓄滿產流。土壤水橫向流動Qlat，i計算式為：

式中：H0，i為土壤層等效水頭，mm；φd，i為土壤層的排水孔隙度，即土壤層孔隙度θs，i與田間持水量θfc，i的差值；Li為土壤層橫向等效長度，m；vlat，i為橫向流動速度，是滲透系數的分量；αi為土壤層坡度，其值一般較小，可認為sinαi=tanαi，即sinαi為等效水頭與橫向距離的比值。若土壤層的一段裸露于地表，則橫向流動量將貢獻地表徑流。

計算土壤水蒸發量時，首先應區分出不同深度土壤層需要的蒸發量。土壤深度層次的劃分決定土壤允許的最大蒸發量，其計算公式為：

式中：Esoil，z為z深度處蒸發需要的水量，mm；z為地表以下土壤層中心的深度，mm；Es為最大可能土壤水分蒸發量，mm；式中的系數項是為了保證50%的蒸發所需水分來自土壤表層0～10mm，95%的蒸發所需水分來自0～100mm土壤深度范圍內。

2.1.1.4 地下水徑流過程模擬

儲存在地下飽和區的水量即為地下水，其按含水層是否具有自由表面分為潛水含水層和承壓含水層。由于承壓含水層一般不發生水平方向的水量交換，并且沒有蒸發，不會對地表徑流產生直接影響。因此，模型僅對潛水含水層進行模擬。

當對河床狀況了解較少時，流域的地下水模擬可使用線性水庫或非線性水庫的概念，即認為整個地下含水層是一個地下水庫。地下水出流計算公式為：

式中：QGi為計算單元i的出口平均地下水出流量，m3/s；SGi為t時刻計算單元i的地下水儲水量，mm；m為指數，m=1為線性水庫，m=2為非線性水庫；cg為考慮計算單元面積后的地下水回歸系數。對每一個計算單元i，其地下水水量平衡表達式如下：

式中：SGi（t）和SGi（t-1）分別為t時刻和t-1時刻的子流域地下水儲水量，mm；Fi為計算單元的面積，m2；EGi（t）為計算單元i的地下水儲水量的計算平均蒸散發量，mm；QGi（t）為計算單元i的地下水出流量，m3/s；Qi，rch（t）為單元i的補給量，mm。Qi，rch（t）由上覆土壤層滲漏量Qi，srch（t）和上一時段的補給量Qi，rch（t-1）加權平均求得：

式中：δgw為有效排水時間，由地質條件決定。

EasyDHM模型中，地下水蒸發量公式為：

式中：EGi（t）為計算單元i的地下水儲水量的平均蒸散發量，mm；EP為潛在蒸散發量，mm；cv為作物系數；EIi（t）為計算單元i的截留水蒸發量，mm；EDi（t）為計算單元i的填洼水蒸發量，mm；ESi（t）為計算單元i的土壤水蒸發量，mm；cd是一個可變參數。cd的計算公式為：

式中：SGi（t）為t時刻計算單元i的地下水儲水量，mm；SGs，0為計算單元i的地下水儲水量的最大值，mm。

2.1.1.5 潛在蒸散發計算

EasyDHM模型提供了Penman—Monteith、Priestley—Taylor和Hargreaves三種計算潛在蒸散發能力的方法，另外還可以讀入實測的水面蒸發能力等直接得到的逐時段的潛在蒸散發數據［145］。本研究采用世界氣象組織推薦的無論在干旱還是濕潤地區計算精度均較高的Penman—Montieth公式法計算潛在蒸發量。

Penman-Montieth公式需要輸入的資料為輻射、氣溫、風速和相對濕度，其公式為：

式中：ET0為潛在蒸散發能力，mm；Δ為飽和水汽壓斜率，kPa/℃；Rn為凈輻射量，MJ/m2；G為土壤熱通量，MJ/m2；ρ為空氣密度，g/m3；D為飽和水汽壓差，kPa；ra為邊界層阻力，s/m；L為汽化潛熱，MJ/kg；γ為濕度計常數。

空氣密度ρ計算公式為：

式中：T為氣溫，℃；PB為大氣壓力，kPa。PB的計算公式為：

式中：Elev為計算點的高程，m。

邊界層阻力ra的公式：

式中：zd為零平面位移高度，m；zd=0.702H0.979，其中H為冠層高度，m；zo為蒸散面粗糙長度，m，zo=0.131H0.997；V為日均風速，m/s。

2.1.2 匯流理論

本次研究采用馬斯京根法對河道匯流進行演算，該法模擬了沿渠道長度柱蓄（Prism Storage）和楔蓄（Wedge Storage）組成的蓄水容量。河段槽柱蓄與楔蓄如圖2-2所示。

洪水波行進到某個河段槽，當入流量大于出流量便形成了楔形蓄水體；當洪水波退去，在河段槽便出現了出流量大于入流量的負楔蓄。另外，對于楔蓄水體，河段槽內始終包含一個體積為流域長度上橫截面不變的柱蓄水體。

總的蓄水容量為：

式中：Vstored為蓄水容量；qin為入流量；qout為出流量；K為穩定流情況下的河段傳播時間；X為流量比重因素。該公式可重新整理為：

流量比重因素X的下限為0.0，上限為0.5。這個因子是楔蓄量的函數。對于水庫式蓄水，沒有楔蓄，X=0.0；而對于一個完全的楔蓄，X=0.5；對于河流，X落在0.0～0.3之間，其平均值接近0.2。

對于蓄水容量的定義可以加入連續公式并簡化為：

式中：qin，1為該時間段開始時的入流量；qin，2為該時間段結束時的入流量；qout，1為該時間段開始時的出流量；qout，2為該時間段結束時的出流量。

其中

C1+C2+C3=1

為用體積單位表示所有值，式兩邊都要乘以該時間段：

為了保持數值穩定和避免負出流量的計算，必須滿足以下條件：

流量比重因素X的值由使用者輸入，蓄水時間常數的值估計如下：

式中：K為穩定流情況下的河段傳播時間，s；coef1和coef2為由使用者輸入的權重系數；Kbnkfull為穩定流情況下渠道蓄滿水的河段傳播時間；K0.1bnkfull為渠道蓄滿1/10水量時河段傳播時間。Cunge（1969）提出了一個計算Kbnkfull和K0.1bnkfull的公式：

式中：K為穩定流情況下的河段傳播時間，s；Lch為渠道長度，km；ck為指定深度處的波速。波速ck計算公式為：

式中：流速qch由曼寧公式定義。

2.1.3 水庫調度

模型將分布式水文模型與水庫調度模型相結合，實現上游來水預報與水庫調度演算的耦合。因此，研究將分布式水文模型EasyDHM與水庫調度模型充分結合，開發水庫群聯合調度模型的同時也對分布式水文模型EasyDHM進行了改進?？紤]水庫的興利、防洪調度作用后，水庫節點的水文過程能得到更合理的重現，并能對下游河道的徑流過程產生影響。此外，由于分布式水文模型存有水系、水文站、水庫的空間拓撲關系，能很方便地描述水庫群上下游間的影響關系，這也是基于分布式水文模型的水庫群聯合調度模型的優點所在。分布式水文模型與水庫群調度模型耦合過程如圖2-3所示，共有三種耦合方式，依據水庫群調度模型中是否包含產流計算和匯流計算分為三類，即產流—匯流—水庫調度的緊密耦合，匯流—水庫調度的簡化，以及兩模型分開模擬的離線耦合方式。

分布式水文模型緊密耦合水庫群調度模塊。這種方式是直接在分布式水文模型中開發一個水庫調度模塊，實現了水文模擬、水庫調度的緊密耦合。使用時，各子流域首先進行子流域內部計算單元的產流計算；然后進行子流域內部的匯流計算，匯流至子流域出口點后判斷子流域出口點是否存在水庫，如果存在水庫則進行水庫調度計算，得到水庫出流作為子流域的出流匯入到下游子流域，否則直接匯入下游子流域，水庫調度模型與水文模型所用時間步長是完全一致的。其優點是使用方便，在計算水文模擬的同時計算完水庫調度，同步輸出各模型結果，方便進行比較。缺點是模型計算時占據電腦內存較大，計算速度相對慢；同時，水文模型輸出的中間結果數據量較大，從水庫調度角度來看并不需要，增加了額外存儲空間。

基于分布式水文模型區間產流的匯流/水庫調度簡化耦合模型。和上述方式不同，它是一個耦合河道匯流模型和水庫調度模型的水庫群調度模型，需事先輸入水文模型的產流模擬結果。這種方式中，匯流模型和水庫調度模型時間步長必須一致，但可以與水文模型（即模擬產流過程）的時間步長不一致。如果時間步長不一致，則需在使用水文模型的產流結果時預先進行時間轉化。這種方式既能節省模型計算時間也能準確考慮河道匯流過程。因此，這種方式的水庫群調度模型可以直接利用分布式水文模型率定好的匯流參數，避免單獨率定水庫群之間河道參數的困難，可以方便地應用到水庫群實時調度中去。

分布式水文模型與水庫群調度模型的離線耦合方式。這種方式的水庫群調度模型和傳統的水庫群調度模型一樣，但區間來水（區間入流）是從分布式水文模型中得出的。首先，借助分布式水文模型EasyDHM計算出各水庫的天然徑流過程，依據各方案水庫間的拓撲關系以及各水庫的天然徑流過程即可得出各水庫的區間入流量，即區間入流量為區間產流量減去河道槽蓄量及損失量。由于水庫群調度模型不包含產流模塊和匯流模塊，在研究中進行相應簡化，給出以下假定：盡管考慮上游水庫的調蓄作用后，上游分區匯入的某水庫分區的流量會發生改變，但認為這種流量改變對該分區的河道槽蓄作用（包含蒸發、滲漏過程）影響不大，即河道的槽蓄過程不隨河道入口點輸入流量的改變而發生變化。這種方式是一種水庫調度模型與分布式水文模型間離線耦合的方式，水庫調度模型與分布式模型的時間步長同樣可以不同。

總的來說，以上三種方式中，第一種方式計算最為耗時，但使用最為方便，對整個流域水循環模擬過程考慮得更為詳盡，計算結果也更為精確；第二種方式是為節省第一種方式計算時間而改進的，因為在同等的氣象及參數條件下，區間產流結果是不會發生變化的，改變水庫調度規則或者設置不同水庫情景改變的只是河道匯流過程；第三種方式則為純粹的水庫群調度模型，通過輸入各水庫區間入流量進行自上而下的水庫群調度模擬，水庫間的河道匯流過程暗含在區間入流中。