本發(fā)明涉及一種新的流域土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法，屬于地球物理下水文分支技術(shù)領(lǐng)域。

背景技術(shù)：

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

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

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題是針對(duì)上述現(xiàn)有技術(shù)存在的問題，而提供一種提高高寒地區(qū)水文模擬和水資源計(jì)算精度的土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法。

為解決上述技術(shù)問題，本發(fā)明采用的技術(shù)方案是：

土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法，主要步驟如下：

(a)根據(jù)土壤溫度分布計(jì)算流域網(wǎng)格內(nèi)不同時(shí)間的凍土深度；

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

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

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

當(dāng)p-e＞0，則產(chǎn)流，否則不產(chǎn)流，

產(chǎn)流量計(jì)算方法為：

若p-e+a＜w’mm則局部產(chǎn)流，有

若p-e+a≥w’mm，則全流域產(chǎn)流，有

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

式中，r為產(chǎn)流量；p為降水量；e是蒸發(fā)量；w0為流域初始土壤蓄水量；a、b為參數(shù)。

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

其中，θu、θi分別為土壤中未凍水、冰的體積含量，t、z分別為時(shí)間及空間坐標(biāo)，d(θu)、k(θu)分別為非飽和凍土水分?jǐn)U散率與導(dǎo)水率，ρi、ρw分別為冰和水的密度，t為土壤溫度，cvs、λ分別為土壤體積熱容量、熱導(dǎo)率，l為融化潛熱，θmax(t)為相應(yīng)土壤負(fù)溫度(t)條件下可能的最大未凍水含量；

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

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

w0＝h*θu(2)

其中，w0為土壤蓄水量，h為網(wǎng)格消融層深度。

所述步驟(b)包括流域單點(diǎn)動(dòng)態(tài)蓄水容量計(jì)算和流域空間動(dòng)態(tài)蓄水容量曲線獲取：

流域單點(diǎn)動(dòng)態(tài)蓄水容量計(jì)算：依據(jù)土壤凍融深度計(jì)算結(jié)果和土壤活動(dòng)層上層土壤凍結(jié)和消融狀態(tài)，獲取流域單點(diǎn)逐日蓄水容量過程；

流域空間動(dòng)態(tài)蓄水容量曲線：利用克里金插值法，對(duì)土壤特征參數(shù)比水容量cw(θu)、導(dǎo)水率k(θu)及擴(kuò)散率d(θu)和土壤熱特征參數(shù)體積比熱容cvs和熱導(dǎo)率λ進(jìn)行空間插值分析，計(jì)算流域內(nèi)每個(gè)網(wǎng)格內(nèi)的土壤凍融深度和土壤活動(dòng)層上層土壤凍結(jié)和消融狀態(tài)，統(tǒng)計(jì)分析土壤凍融深度的逐日空間分布狀況，繪制流域空間動(dòng)態(tài)逐日蓄水容量曲線。

步驟(c)中，基于動(dòng)態(tài)蓄水容量曲線的凍土產(chǎn)流計(jì)算方法，利用蓄滿產(chǎn)流原理，計(jì)算得到流域入滲到土壤中的水量和產(chǎn)流兩部分。

本發(fā)明基于多點(diǎn)凍土觀測資料以及凍土下滲實(shí)驗(yàn)，分析土壤凍融過程中蓄水容量變化和下滲產(chǎn)流機(jī)理，進(jìn)而推求流域土壤凍融動(dòng)態(tài)蓄水容量曲線，提出一種新的流域土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法。

高寒地區(qū)是氣候變化的敏感區(qū)，主要是因?yàn)殡S著氣溫升高，冰川、融雪和凍土對(duì)水文過程影響機(jī)理更加復(fù)雜，考慮到目前用于流域水文過程模擬的模型都沒有考慮到土壤凍融過程動(dòng)態(tài)蓄水容量對(duì)產(chǎn)流的影響，不能詳盡刻畫凍土產(chǎn)流過程、春汛洪水的模擬以及水資源演變計(jì)算，本發(fā)明把土壤水熱耦合遷移運(yùn)動(dòng)與下滲產(chǎn)流有機(jī)結(jié)合，提出一種新的流域土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法，既可以計(jì)算土壤非凍結(jié)期產(chǎn)流，也可以計(jì)算土壤凍融過程中的產(chǎn)流。

本發(fā)明可以依據(jù)觀測氣溫模擬逐日土壤凍結(jié)/消融、凍土深度和土壤溫度，依據(jù)降雨觀測計(jì)算融雪/降雨徑流的逐日過程，提高了春季土壤消融期的徑流模擬精度，為春汛防洪決策提供科學(xué)依據(jù)，同時(shí)也填補(bǔ)了現(xiàn)行國內(nèi)外水文模型中凍土區(qū)產(chǎn)流計(jì)算的空白。

本發(fā)明采用以上技術(shù)方案與現(xiàn)有技術(shù)相比，具有以下技術(shù)效果：本發(fā)明利用土壤水熱耦合遷移數(shù)值模擬土壤凍結(jié)和消融過程、土壤不同深度溫度變化，研究凍土變化與包氣帶蓄水容量關(guān)系，提出土壤凍融動(dòng)態(tài)土壤蓄水容量產(chǎn)流計(jì)算方法，發(fā)展動(dòng)態(tài)土壤蓄水容量的產(chǎn)流模塊，為高寒地區(qū)土壤凍融過程中的產(chǎn)流計(jì)算提供了一種新的方法，該流域土壤凍融過程中動(dòng)態(tài)蓄水容量產(chǎn)流計(jì)算方法可以提高高寒地區(qū)水文模擬和水資源計(jì)算精度，推進(jìn)寒區(qū)水文學(xué)的發(fā)展。

附圖說明

圖1為本發(fā)明流域土壤凍融過程中動(dòng)態(tài)蓄水容量曲線計(jì)算方法流程圖；

圖2為本發(fā)明產(chǎn)流計(jì)算技術(shù)路線圖；

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

圖4為基于流域蓄水容量曲線的產(chǎn)流量計(jì)算示意圖；

圖5為多組變動(dòng)蓄水容量曲線示意圖；

圖6為本發(fā)明實(shí)施例中的凍土觀測點(diǎn)凍土深度逐日變化過程圖；

圖7為本發(fā)明實(shí)施例中的凍土觀測點(diǎn)逐日氣溫變化過程圖；

圖8為本發(fā)明實(shí)施例中土壤凍融深度模擬結(jié)果與實(shí)測結(jié)果對(duì)比圖；

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

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

圖11為本發(fā)明實(shí)施例中流域逐日流量模擬結(jié)果與實(shí)測結(jié)果對(duì)比圖。

具體實(shí)施方式

下面結(jié)合附圖和具體實(shí)施例，進(jìn)一步闡明本發(fā)明。

本文以黃河源區(qū)某地區(qū)為例，采用本發(fā)明的方法對(duì)該地區(qū)土壤凍融過程中的產(chǎn)流進(jìn)行計(jì)算。

具體包括如下步驟：

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

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

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

第二步：利用土壤水熱耦合遷移方程構(gòu)建數(shù)值模擬模型，基于土壤特征參數(shù)、土壤熱特征參數(shù)和黃河源區(qū)凍土深度和溫度實(shí)驗(yàn)觀測值，計(jì)算流域網(wǎng)格內(nèi)不同深度的土壤溫度和凍融深度隨時(shí)間變化過程。根據(jù)計(jì)算土壤溫度分布變化，識(shí)別不同時(shí)間低于0℃的土壤剖面分布，可以得到不同時(shí)間的土壤凍結(jié)厚度、凍結(jié)位置和凍結(jié)鋒面，從而得到土壤活動(dòng)層上層土壤凍結(jié)和消融狀態(tài)。

1)在凍融過程中非飽和土水熱耦合遷移過程的研究中，認(rèn)為凍土中水分遷移規(guī)律與非飽和土壤水運(yùn)動(dòng)規(guī)律類似，可用含相變項(xiàng)的richards方程表示，以自變量為θ的richards方程：

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

該方程其特點(diǎn)是便于用數(shù)值模擬方法進(jìn)行求解，適用于均質(zhì)非飽和水分運(yùn)動(dòng)。

將相變潛熱作為內(nèi)熱源的傳導(dǎo)方程為：

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

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

θu≤θmax(t)(6)

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

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

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

利用實(shí)驗(yàn)或者理論方法得到其中兩個(gè)參數(shù)，其它參數(shù)即可計(jì)算得到。土壤水分特征曲線可以在野外或?qū)嶒?yàn)室內(nèi)測定，也可以通過vg模型計(jì)算得到。

4)土壤熱特征參數(shù)計(jì)算。土壤熱特性參數(shù)包括體積比熱容cvs和熱導(dǎo)率λ，可以由試驗(yàn)測定，也可用半經(jīng)驗(yàn)半理論公式計(jì)算。

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

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

7)土壤凍融深度的計(jì)算；根據(jù)計(jì)算土壤溫度分布變化，識(shí)別不同時(shí)間低于0℃的土壤剖面分布，可以得到不同時(shí)間的土壤凍結(jié)厚度、凍結(jié)位置和凍結(jié)鋒面，從而計(jì)算出流域網(wǎng)格內(nèi)不同深度的土壤溫度隨時(shí)間變化過程，為動(dòng)態(tài)蓄水容量計(jì)算提供數(shù)據(jù)。

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

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

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

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

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

第五步：利用蓄滿產(chǎn)流原理，基于動(dòng)態(tài)蓄水容量曲線的流域產(chǎn)流計(jì)算方法，計(jì)算得到流域入滲到土壤中的水量δw和產(chǎn)流兩部分。如圖4，若初始土壤含水

量為w。，則

當(dāng)p-e＞0，則產(chǎn)流，否則不產(chǎn)流，產(chǎn)流量計(jì)算方法為：

若p-e+a＜w’mm則局部產(chǎn)流，有

δw＝p-e-r(12)

若p-e+a≥w’mm，則全流域產(chǎn)流，有

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

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

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

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

利用本發(fā)明提出的土壤凍融動(dòng)態(tài)土壤蓄水容量產(chǎn)流計(jì)算方法，發(fā)展動(dòng)態(tài)土壤蓄水容量的產(chǎn)流模塊，將模擬的地表徑流量和實(shí)際觀測值作對(duì)比，如圖11所示，模擬的地表徑流量與實(shí)測的流量比較接近，相對(duì)誤差為4％，確定性系數(shù)為0.89，模擬精度較高，說明該發(fā)明提出的研究方法在高寒區(qū)具有較好的適用性。