土壤凍融過程中動態蓄水容量產流計算方法與流程

本發明涉及一種新的流域土壤凍融過程中動態蓄水容量產流計算方法，屬于地球物理下水文分支技術領域。

背景技術：

流域產流作為水循環中最重要的一個環節，產流計算方法是流域水文模擬、水資源計算和洪水預報的重要理論基礎。上世紀三十年代horton提出了超滲產流的概念，即當雨強大于下滲能力時產生地表徑流，二十世紀六十年代河海大學趙人俊提出了蓄滿產流的概念，即當降雨量滿足土壤蓄水容量后產生徑流(包括地表徑流、壤中流和基流)，自此以后發展的水文模型都以超滲或蓄滿作為產流計算的理論方法，近年來也有人認為流域內在不同時間會有兩種產流機制共存的產流方式，提出混合產流計算方法。

多年和季節凍土廣布于高寒地區，其水熱傳輸貫穿于寒區流域的產流、入滲和蒸散發過程中，是寒區水文過程的核心環節。而常用的單一蓄滿產流、單一超滲產流以及蓄超混合產流方式都不能涵蓋高寒地區土壤凍融過程中的產流機理，新的流域土壤凍融過程中動態蓄水容量產流計算方法不僅可以提高高寒地區水文模擬和水資源計算精度，同時也大大豐富了寒區水文學的內涵。已有研究表明凍土不會阻止融雪/降雨水分下滲，融雪水/降雨在土壤凍結和不凍結情況下都會下滲到土壤中。關于凍土下滲產流問題，已有許多試驗和理論研究，但是目前還沒有能夠用于實際且滿足寒區流域徑流預報的下滲強度和凍土透水程度的試驗資料。在國內外廣泛應用的水文模型(如vic模型、swat模型、新安江模型等)和陸面模式(如clm)中，至今沒有考慮土壤凍融動態蓄水容量與產流的關系和相應的凍土產流計算方法。

技術實現要素：

本發明所要解決的技術問題是針對上述現有技術存在的問題，而提供一種提高高寒地區水文模擬和水資源計算精度的土壤凍融過程中動態蓄水容量產流計算方法。

為解決上述技術問題，本發明采用的技術方案是：

土壤凍融過程中動態蓄水容量產流計算方法，主要步驟如下：

(a)根據土壤溫度分布計算流域網格內不同時間的凍土深度；

(b)分析凍土深度的逐日空間分布狀況，得到流域不同網格內包氣帶達到的隨時間分布的田間持水量w’m；

(c)找出流域不同網格內最大的田間持水量w’mm，并計算流域最大持水量wm，流域最大持水量wm為：

(d)根據得到的不同網格內最大的田間持水量w’mm和流域最大持水量wm，計算流域土壤凍融產流：

當p-e＞0，則產流，否則不產流，

產流量計算方法為：

若p-e+a＜w’mm則局部產流，有

若p-e+a≥w’mm，則全流域產流，有

r＝p-e-(wm-w0)(13)

式中，r為產流量；p為降水量；e是蒸發量；w0為流域初始土壤蓄水量；a、b為參數。

采用土壤水熱耦合遷移模型計算步驟a中的凍土深度，其中土壤水熱耦合遷移模型為：

其中，θu、θi分別為土壤中未凍水、冰的體積含量，t、z分別為時間及空間坐標，d(θu)、k(θu)分別為非飽和凍土水分擴散率與導水率，ρi、ρw分別為冰和水的密度，t為土壤溫度，cvs、λ分別為土壤體積熱容量、熱導率，l為融化潛熱，θmax(t)為相應土壤負溫度(t)條件下可能的最大未凍水含量；

凍土深度由溫度小于0位置處的空間坐標確定。

基于土壤水熱耦合遷移模型推求的網格點凍土深度，通過方程(2)轉換為土壤蓄水量：

w0＝h*θu(2)

其中，w0為土壤蓄水量，h為網格消融層深度。

所述步驟(b)包括流域單點動態蓄水容量計算和流域空間動態蓄水容量曲線獲取：

流域單點動態蓄水容量計算：依據土壤凍融深度計算結果和土壤活動層上層土壤凍結和消融狀態，獲取流域單點逐日蓄水容量過程；

流域空間動態蓄水容量曲線：利用克里金插值法，對土壤特征參數比水容量cw(θu)、導水率k(θu)及擴散率d(θu)和土壤熱特征參數體積比熱容cvs和熱導率λ進行空間插值分析，計算流域內每個網格內的土壤凍融深度和土壤活動層上層土壤凍結和消融狀態，統計分析土壤凍融深度的逐日空間分布狀況，繪制流域空間動態逐日蓄水容量曲線。

步驟(c)中，基于動態蓄水容量曲線的凍土產流計算方法，利用蓄滿產流原理，計算得到流域入滲到土壤中的水量和產流兩部分。

本發明基于多點凍土觀測資料以及凍土下滲實驗，分析土壤凍融過程中蓄水容量變化和下滲產流機理，進而推求流域土壤凍融動態蓄水容量曲線，提出一種新的流域土壤凍融過程中動態蓄水容量產流計算方法。

高寒地區是氣候變化的敏感區，主要是因為隨著氣溫升高，冰川、融雪和凍土對水文過程影響機理更加復雜，考慮到目前用于流域水文過程模擬的模型都沒有考慮到土壤凍融過程動態蓄水容量對產流的影響，不能詳盡刻畫凍土產流過程、春汛洪水的模擬以及水資源演變計算，本發明把土壤水熱耦合遷移運動與下滲產流有機結合，提出一種新的流域土壤凍融過程中動態蓄水容量產流計算方法，既可以計算土壤非凍結期產流，也可以計算土壤凍融過程中的產流。

本發明可以依據觀測氣溫模擬逐日土壤凍結/消融、凍土深度和土壤溫度，依據降雨觀測計算融雪/降雨徑流的逐日過程，提高了春季土壤消融期的徑流模擬精度，為春汛防洪決策提供科學依據，同時也填補了現行國內外水文模型中凍土區產流計算的空白。

本發明采用以上技術方案與現有技術相比，具有以下技術效果：本發明利用土壤水熱耦合遷移數值模擬土壤凍結和消融過程、土壤不同深度溫度變化，研究凍土變化與包氣帶蓄水容量關系，提出土壤凍融動態土壤蓄水容量產流計算方法，發展動態土壤蓄水容量的產流模塊，為高寒地區土壤凍融過程中的產流計算提供了一種新的方法，該流域土壤凍融過程中動態蓄水容量產流計算方法可以提高高寒地區水文模擬和水資源計算精度，推進寒區水文學的發展。

附圖說明

圖1為本發明流域土壤凍融過程中動態蓄水容量曲線計算方法流程圖；

圖2為本發明產流計算技術路線圖；

圖3為流域蓄水容量曲線概化圖；

圖4為基于流域蓄水容量曲線的產流量計算示意圖；

圖5為多組變動蓄水容量曲線示意圖；

圖6為本發明實施例中的凍土觀測點凍土深度逐日變化過程圖；

圖7為本發明實施例中的凍土觀測點逐日氣溫變化過程圖；

圖8為本發明實施例中土壤凍融深度模擬結果與實測結果對比圖；

圖9為本發明實施例中不同深度下(5cm)的土壤溫度模擬結果與實測結果對比圖；

圖10為本發明實施例中不同深度下(20cm)的土壤溫度模擬結果與實測結果對比圖；

圖11為本發明實施例中流域逐日流量模擬結果與實測結果對比圖。

具體實施方式

下面結合附圖和具體實施例，進一步闡明本發明。

本文以黃河源區某地區為例，采用本發明的方法對該地區土壤凍融過程中的產流進行計算。

具體包括如下步驟：

第一步：從中國氣象數據網(http：//data.cma.cn)上下載研究區域氣象站點逐日降水(雪)、日平均氣溫、日最高溫度、日最低溫度和0cm地溫資料，利用克里金插值法，對下載資料進行空間插值分析，生成流域內每個網格的逐日資料系列。高寒地區多以融雪與降水為主要補給源的山區性流域，其中融雪量計算采用度日因子模型：：

m＝cm×(ti-tb)+ceer(3)

式中，m為日平均積雪消融量，cm為融雪的度日因子，ti為第i個網格融雪的日平均溫度(℃)；tb為雪消融的臨界溫度，ce為雪的輻射系數，er為太陽短波輻射或凈輻射。

第二步：利用土壤水熱耦合遷移方程構建數值模擬模型，基于土壤特征參數、土壤熱特征參數和黃河源區凍土深度和溫度實驗觀測值，計算流域網格內不同深度的土壤溫度和凍融深度隨時間變化過程。根據計算土壤溫度分布變化，識別不同時間低于0℃的土壤剖面分布，可以得到不同時間的土壤凍結厚度、凍結位置和凍結鋒面，從而得到土壤活動層上層土壤凍結和消融狀態。

1)在凍融過程中非飽和土水熱耦合遷移過程的研究中，認為凍土中水分遷移規律與非飽和土壤水運動規律類似，可用含相變項的richards方程表示，以自變量為θ的richards方程：

式中，θu、θi分別為土壤中未凍水、冰的體積含量，t、z分別為時間及空間坐標(垂直向下為正)，d(θu)、k(θu)非飽和凍土水分擴散率與導水率，ρi、ρw為冰和水的密度。

該方程其特點是便于用數值模擬方法進行求解，適用于均質非飽和水分運動。

將相變潛熱作為內熱源的傳導方程為：

式中，t為土壤溫度，cvs、λ為土壤體積熱容量、熱導率，l為融化潛熱。

上述(4)和(5)為土壤凍融過程中水熱耦合遷移的兩組基本方程，但是需要求解的是三個未知函數，即θu(z，t)、θi(z，t)和t(z，t)。因此還必須補充一個聯系方程，即土壤中未凍水含水率θu與溫度t的關系方程。在一定的負溫度下，凍土中總含有部分未凍結水θu，并與負溫度、壓力等條件下處于動力平衡狀態，在凍土研究中，當外界壓力一定時，未凍水含量是溫度的函數，可表示為凍土中水與熱運動之間的聯系：

θu≤θmax(t)(6)

式中，θmax(t)為相應土壤負溫度(t)條件下可能的最大未凍水含量。

2)設定初始和邊界條件。初始條件中的含水量分布θ0(z)和溫度分布t(z)是已知的，土壤水分遷移的上邊界條件是融雪/降雨入滲或者土壤蒸散發，下邊界條件可為定水位和無限深度。土壤熱流的邊界條件在第一邊界條件時，已知地表(z＝0)處溫度隨時間的變化過程t(t)及下邊界處溫度維持不變，為t(l)＝c。

3)土壤水分特征參數計算。與土壤水分遷移相關特征參數包括土壤水分特征曲線(土水勢ψ或吸力s與土壤含水量關系)、比水容量cw(θu)、導水率k(θu)、擴散率d(θu)，各參數存在以下關系：

利用實驗或者理論方法得到其中兩個參數，其它參數即可計算得到。土壤水分特征曲線可以在野外或實驗室內測定，也可以通過vg模型計算得到。

4)土壤熱特征參數計算。土壤熱特性參數包括體積比熱容cvs和熱導率λ，可以由試驗測定，也可用半經驗半理論公式計算。

5)凍土水熱耦合遷移方程離散化要求。利用有限差分法求解，將計算區域離散化。由于土壤凍結鋒面處水的相變引起釋放大量潛熱和凍土消融過程中吸收潛熱，離散化時應取合適的距離步長和時間步長，在凍結鋒面處距離步長取小一些。

6)不同深度土壤溫度、未凍水含量和冰含量的計算。利用中心差分格式對凍土水熱耦合遷移方程進行數值求解，可以計算出不同深度的土壤溫度、未凍水含量和冰含量隨時間變化過程。

7)土壤凍融深度的計算；根據計算土壤溫度分布變化，識別不同時間低于0℃的土壤剖面分布，可以得到不同時間的土壤凍結厚度、凍結位置和凍結鋒面，從而計算出流域網格內不同深度的土壤溫度隨時間變化過程，為動態蓄水容量計算提供數據。

第三步：統計分析土壤凍融深度的逐日空間分布狀況，可得到不同網格下包氣帶達到田間持水量隨時間的分布，繪制若干組不同時間下的流域動態蓄水容量曲線(圖5)。

基于土壤水熱耦合遷移模型推求的單個網格點凍土深度，可通過方程(2)轉換為土壤蓄水容量。

流域上各個網格點包氣帶厚薄及土壤特性一般不相同，當全流域處于最干旱狀態時，各處的包氣帶的缺水量不一定，即各處的包氣帶達到田間持水量不一樣，其中最大的田間持水量為w’mm。將全流域面積看做1，以包氣帶田間持水量為縱坐標，小于等于某一田間持水量所占的流域面積比重為橫坐標α，所得到的曲線(如圖3)稱為流域蓄水容量曲線：

曲線所包圍的全部面積等于流域平均蓄水容量或最大持水量wm。

式中w’m為流域某處包氣帶達到的田間持水量，α值表示流域中≤w’m的流域面積所占的比重，b為流域蓄水容量曲線的方次，一般取值0.2～0.4，表征蓄水容量分布不均勻性的參數，b越大代表流域蓄水容量分布越不均勻。

第五步：利用蓄滿產流原理，基于動態蓄水容量曲線的流域產流計算方法，計算得到流域入滲到土壤中的水量δw和產流兩部分。如圖4，若初始土壤含水

量為w。，則

當p-e＞0，則產流，否則不產流，產流量計算方法為：

若p-e+a＜w’mm則局部產流，有

δw＝p-e-r(12)

若p-e+a≥w’mm，則全流域產流，有

r＝p-e-(wm-w0)(13)

式中，w0～為流域初始土壤蓄水量(mm)；r～為產流量(mm)。

在本實施例中，選擇黃河源區某區域作為研究區域，黃河源區一般是指河源至唐乃亥之間的區域，海拔高度在3000m以上，地處青藏高原的東北部，地理位置在95°50′～103°30′，32°20′～35°50′n之間。流域內屬高原大陸性氣候，主要為濕潤半濕潤氣候區，多年平均氣溫為-4-5.2℃，年日照時數為2250-3131小時，平均風速3-4.5m/s。

為了驗證本發明方法的實施，選擇7月1日到翌年的6月30日為一周期，這個時間段能完全把黃河源區氣象站的凍融期包括在內，實測數據包括該地區1997-2007年的凍土資料、地表溫度和徑流資料，圖6和圖7分別是該地區凍土觀測點的凍土深度和氣溫的部分逐日過程，基于該實測資料對構建的土壤水熱耦合遷移模型進行參數率定和驗證，從而模擬區域單個網格點的土壤凍融深度和不同深度的土壤溫度變化(圖8-10)，從模擬結果來看，模型模擬效果較好，可為繪制流域空間動態逐日蓄水容量曲線提供數據。

利用本發明提出的土壤凍融動態土壤蓄水容量產流計算方法，發展動態土壤蓄水容量的產流模塊，將模擬的地表徑流量和實際觀測值作對比，如圖11所示，模擬的地表徑流量與實測的流量比較接近，相對誤差為4％，確定性系數為0.89，模擬精度較高，說明該發明提出的研究方法在高寒區具有較好的適用性。