越冬期麥田根區水熱耦合模型參數原位估算與檢驗

程 強 徐 嬙 陳 超 薛緒掌 王忠義 孫宇瑞

(1.中國農業大學信息與電氣工程學院， 北京 100083； 2.北京農業智能裝備技術研究中心， 北京 100097)

0 引言

越冬期作物根區充足的蓄水量能夠減小作物對于春季降水的依賴程度，有利于春季冬小麥的生長并有可能增加產量[1]。土壤中水分平衡的精確預測不僅是決定水資源可用性的關鍵，也對其優化管理起到了重要作用[2]。

在季節性凍土區，越冬期冬小麥根區水分的補給渠道除降水和田間灌溉外，因上層土壤凍結驅動下層土壤水分向上運移也是冬小麥根區水分補給的重要來源。隨著溫度降至0℃以下，土壤上層液態水變成冰，降低了該層的水勢。在水勢梯度作用下，下層土壤的水分開始向上運移，形成一個水熱耦合運移機制。以往許多學者用復雜的數學模型對該機制進行了研究[3-8]。

在越冬期土壤的水熱耦合運移過程中，土壤未凍水含量、含冰量和土壤溫度是必不可少的3個關鍵變量[8]。土壤凍融特征曲線(Soil freezing characteristic curve, SFC)用于表示凍土中土壤未凍水含量和溫度之間的關系，在模擬凍土的水熱耦合運移過程中起到了重要作用。由于SFC和土壤水分特征曲線(Soil moisture characteristic curve, SMC)之間的相似性[9]，一些學者利用SMC來估算SFC進行凍融過程中的水熱耦合運移模擬[10]。為提高模型的預測精度，一些學者[9,11]將針式時域反射計(TDR)和溫度傳感器結合獲取SFC，但是在傳感器安裝過程中容易破壞土壤結構。而介電管式傳感器能夠原位測定土壤未凍水含量，提高了測量精度[12-13]。基于原位測量的土壤未凍水含量和溫度數據，得到的水熱參數有具體的物理意義，比參數辨識(反演)法所使用的模型結構簡單，并取代田間取樣后在實驗室測定土壤參數的傳統做法。

此外，由于冬季作物進入休眠期，基本不考慮蒸騰作用，而土壤表面蒸發量作為邊界通量也對根區的蓄水和熱量平衡起著重要作用。土壤上邊界的蒸發量不僅影響地表熱通量的變化[14]，還降低了表層土壤的未凍水含量。然而，以往的研究往往認為冬季地表蒸發量可能很小甚至可能為零[15-17]。與這些研究不同的是，FLERCHINGER等[4]在其模型中引入了一個地表蒸發方程，以完善越冬期水熱耦合運移(Simutaneous heat and water, SHAW)機制，并在之后的研究中用氣象數據對模型進行了測試[14,18-19]。 SCOTT等[1]利用HYDRUS軟件[7]來模擬越冬期的非飽和流，通過考慮由最大向上流速或彭曼潛在蒸發計算得到地表蒸發通量(以較小者為準)。SAITO等[20]和ZENG等[21-22]提出了水-熱-蒸汽耦合的數值模型，并證明地表蒸發量對于土壤中水分運移有重要影響。ZHANG等[23]提出的耦合模型中計算了凍融過程中液態水和蒸發量對于水熱耦合運移的貢獻率，結果證明水蒸氣在所有的土壤深度處貢獻率都超過了15%，驗證了蒸發量作為土壤凍融模型的上邊界條件，在模擬過程尤其是在淺層土壤中是不可忽略的。

本文基于原位測量的土壤未凍水含量和溫度(2011—2012年)，提出一種改進型的凍融過程中水熱耦合運移模型參數估算方法，并將估算參數用于模型精度的驗證過程(2012—2013年)。考慮到地表蒸發量對土壤水分運移與再分布的重要性，本文擬構建土壤凍融過程水熱耦合運移模型用于模擬越冬期小麥根區蓄水與能量交換，并探討冬季地表蒸發量對土壤水熱耦合模型預測精度的影響。

1 材料與方法

1.1 模型介紹

1.1.1水流通量傳輸過程

FLERCHINGER等[4]在1989年提出使用一階偏微分方程來模擬凍土和非凍土中的水熱耦合運移，其中水流通量的傳輸方程[24](向下為正)為

(1)

(2)

式中Ks——飽和導水率，m/h

θs——飽和含水量，m3/m3

b——土壤孔隙分布參數

θi——土壤冰含量，m3/m3

θl——未凍水含量，m3/m3

ρi——冰密度，kg/m3

ρl——液態水密度，kg/m3

qv——蒸汽通量，m/h

K——非飽和導水率，m/h

ψ——土壤基質勢，m

U——水通量的匯項，m3/(m3·h)

z——土壤深度，mt——時間，h

式(1)各項依次代表：體積含水率的變化量；體積含冰量的變化量；流入土壤層的凈液態水通量；進入土壤層的水蒸氣含量；由于根系吸收而形成的一個匯項。

土壤基質勢ψ可通過Campbell模型[25]計算，即

(3)

式中ψe——空氣進入勢，kPa

由Richards方程可以得到

(4)

式中D——土壤中水分的擴散率，m/h

地表蒸發量作為土壤表層的水流通量，影響了凍融過程中的水熱耦合運移，則由式(1)、(4)可得

(5)

由達西定律可得

(6)

式中ql——水流通量，m/h

在上邊界處ql為地表蒸發量。該模型以上邊界上的水流通量(即蒸發量)作為驅動力，計算下一個迭代時間上的液態水含量；另一方面，凍土中水分的流動同時會產生熱對流，進而影響土壤中的熱量傳輸。

1.1.2熱量傳輸過程

土壤基質中熱量傳輸(包括液態水產生的熱對流和蒸發吸收的潛熱)方程為

(7)

式中Cs——土壤體積熱容，J/(m3·℃)

T——土壤溫度，℃

ρv——水蒸氣密度，kg/m3

Lf——融化潛熱，kJ/kg

Lv——蒸發潛熱，kJ/kg

λs——土壤熱傳導率，W/(m·℃)

S——土壤層的源項，m3/(m3·h)

式(7)中各項依次代表：因溫度升高儲存的比熱；水凍結成冰所需要的潛熱；土壤層中的凈蒸發潛熱；進入土壤層的凈熱傳導；由于液態水通量產生的熱對流；土壤層的源項(可能包括太陽輻射和長波輻射)。當不考慮地表蒸發的影響時，忽略凈蒸發潛熱項。

土壤體積熱容Cs是土壤各組分的體積熱容之和，其計算公式為

Cs=∑ρjcjθj

(8)

式中ρj、cj、θj——土壤的第j個組成成分的密度、比熱容和體積分數

DE VRIES[26]提出了用于計算土壤的熱傳導率λs的理論，該理論認為非常濕潤的土壤可以看作是混合了液態水、土壤顆粒、冰晶以及分散的空氣的連續介質。這種理想模型下的土壤熱導率的計算公式為

(9)

式中mj、λj為土壤的第j個組成成分的權重因子、熱傳導率，例如砂粒含量、粉粒含量、黏粒含量、有機質含量、冰、液態水、空氣含量等。DE VRIES[26]討論了關于權重因子mj的計算方法。

1.2 參數估算

土壤凍融特征線(SFC)是指凍土中的溫度(0℃之下)和未凍水含量之間的關系，在構建凍融過程中土壤中水熱耦合運移，尤其是優化數值模型的參數時發揮了重要作用[12]。當土壤中出現冰時，液態水和固態冰的壓強與溫度有關，它們之間的關系為

(10)

式中pl——液相的壓強，Pa

pi——冰的壓強，Pa

π——土壤溶質的滲透壓，Pa

TK——絕對溫度，K

假設土壤在低鹽度條件下，滲透壓為零；且冰中為零壓力，則方程(10)可簡化為[27]

(11)

式中g——重力加速度，m/s2

假設未凍土中水流通量和土壤水分特征線在凍土中同樣適用，則由方程(3)、(11)可得

(12)

凍土中同時測量得到的溫度和液態水含量確定了SFC，根據方程(12)可以得到T<0時，θl與T之間為冪函數關系(圖1)，通過曲線擬合可得到相應的統計學參數，進而反演出ψe和b的值，分別為-127 kPa和2.5。

圖1 土壤凍融特征曲線 Fig.1 Soil freezing characteristic curves

圖2 2011—2012年和2012—2013年兩個越冬期土壤凍融上邊界條件 Fig.2 Upper boundary conditions of air temperature and soil surface evaporation for soil freezing-thawing process over two winters of 2011—2012 and 2012—2013

1.3 邊界條件

模型中的邊界條件分為第1類邊界條件和第2類邊界條件。其中第1類邊界條件是指在上邊界給出了空氣溫度值；第2類邊界條件是指在上邊界給出了向上的水流通量(即蒸發量)的值(圖2)。從2011—2012年和2012—2013年兩個冬季的空氣溫度和蒸發量的值可知，2011—2012年冬季空氣最低溫度為-18.16℃，在11月24日左右溫度降至0℃左右，之后溫度一直下降至0℃以下，進入凍結期，直至3月中旬時溫度上升至0℃以上，且存在明顯的融化潛熱階段；2012—2013年空氣最低溫度為-22.68℃，凍融期約為2012年11月29日至2013年3月1日，空氣溫度波動較大，融化潛熱階段不明顯。

1.4 研究區概況

試驗區位于北京市昌平區小湯山精準農業基地(116°26′39″E，40°10′43″N)，屬于大陸性季風氣候，年平均降雨量為587 mm，年平均溫度為11.7℃。該基地屬于季節性凍土區，每年11—12月土壤開始凍結，2—3月開始融化，越冬期(11—3月)平均降雨量約為42 mm，土壤平均溫度約為-1.1℃。按照國際制土壤質地分級標準，試驗區土壤被劃分為粉砂質壤土，其砂粒、粉粒和粘粒所占比例分別為39.9%、46.6%和13.5%，土壤飽和含水量為0.41 m3/m3，飽和導水率為0.01 m/h。

進行田間試驗時，在每個測量點處，將6根PVC管(長2 m、直徑50 mm)以2 m的間隔垂直插入土壤中，使用介電管式傳感器(Diviner 2000型, Sentek Sensor Technologies)進行不同深度土壤水分的測量，測量間隔為10 cm。土壤未凍水含量可通過介電常數計算得到[22]，測量點的液態水含量取6個點的平均值，其水平和垂直方向上的其他位置數據通過插值計算獲得。土壤溫度使用一組數字化傳感器DS18B20型(DS18B20+， Maxim Integrated)進行測量，其精度為0.5℃(在-55～125℃之間)。在距介電管式傳感器1 m遠的位置挖洞，直徑為5 cm，將溫度傳感器以10 cm的間隔放入洞中。在放置溫度傳感器時，不斷地將該位置的原土壤填回至洞中[8]。

土壤蒸發量作為模型的邊界條件，使用基地的稱重式蒸滲儀系統得到。試驗所用蒸滲儀數據為24套中型蒸滲儀(長1 m、寬0.75 m、高2.3 m)測量的平均值。每套蒸滲儀采用杠桿式稱重系統，在利用平衡塊抵消土箱和土體質量后，使用質量傳感器測量土壤中水分質量，以反映土壤的蒸發量的變化[28]。傳感器的測量頻率為每15 min一次，靈敏度為0.05～0.1 mm。

1.5 誤差分析

(1)均方根誤差(Root mean square error, RMSE)用于反映估算值和實測值之間的總體差異，對特大或特小誤差反映敏感。當RMSE越接近于0時，表明估算誤差越小，模擬精度越高。其計算公式為

(13)

式中Mi——實測值Si——模擬值

N——實測樣本數

(2)平均偏差(Mean bias error, MBE)用于反映估算值和實測值之間的平均偏差，正值表示蒸散量被高估，負值表示被低估。MBE越接近于0時，精度越高。其計算公式為

(14)

(3)d值表示模擬值和實測值之間的一致性和精確度，取值范圍0≤d≤1，當d越接近1時，代表模型的模擬值精確度越接近100%。其計算公式為

(15)

2 結果與討論

2.1 2011—2013年兩年冬季凍融過程的實測和模擬結果比較

春季冬小麥返青期之前其根系分布范圍約為0～20 cm[29-30]，因此根據土壤10 cm和20 cm處未凍水含量和溫度的模擬值和實測值的對比情況來分析冬小麥根區的水熱運移(圖3、4)。在2011—2012年冬季，基于原位測量的土壤未凍水含量和土壤溫度(負溫條件時)得到的SFC參數和考慮了地表蒸發量的水熱耦合運移模型進行數值模擬，土壤溫度模擬值與實測值擬合度較高，模擬效果較好(d值均在0.95左右，表1)。由于凍土中未凍水含量主要受溫度影響，未凍水含量的模擬精度也達到較高水平，其中，10 cm處d值為0.930，RMSE為0.036 m3/m3，20 cm處d值為0.870，RMSE為0.055 m3/m3。在溫度逐漸下降至0℃，土壤20 cm處的未凍水逐漸凍結的過程中，模擬值出現一定的低估，但是土壤溫度的模擬值和實測值之間該現象并不明顯。出現這種情況的原因可能是在SFC的標定過程中(圖1)，使用整個凍融過程的凍土溫度和未凍水含量數據進行曲線擬合，但SFC曲線在凍結-融化過程中存在滯后效應，土壤凍結時速度較快，難以滿足穩態或準穩態的假設條件[10,31]，在凍結階段，相同的溫度條件下擬合曲線的未凍水含量值相對于實測值偏低。在融化階段，由于不同深度處土壤溫度的高估導致未凍水含量出現不同程度的高估，如2011—2012年冬季20 cm處融化階段的溫度高估了0.3℃，但是未凍水含量卻上升了約0.1 m3/m3。0℃附近較小的溫度變化可能引起未凍水含量的較大差異(圖1)，這主要是受融化潛熱的影響，而擬合的SFC曲線在融化階段的相對高估也是影響未凍水含量預測精度的原因之一。

將2012—2013年的土壤表層蒸發量、空氣溫度、土壤初始溫度和未凍水含量作為輸入數據，利用數值模型和已估算的SFC參數可得到未凍水含量的預測值，將該預測值和實測值進行對比分析來驗證模型的適用性。整個凍融過程中，預測值能夠較好地追蹤到土壤中溫度和未凍水含量的動態變化，土壤溫度的模擬值整體偏低，在土壤凍結后以及融化過程中精度相對提高(圖3、4)。未凍水含量在10 cm和20 cm處模擬值的d值分別為0.924和0.774，RMSE分別為0.046 m3/m3和0.071 m3/m3(表2)。在凍結階段，土壤溫度的模擬值出現了明顯的低估，從而導致未凍水含量在凍結過程中的低估(圖3)。

圖3 2011—2012和2012—2013兩個越冬期土壤10 cm和20 cm處未凍水含量模擬值與實測值比較 Fig.3 Comparisons of simulated and measured soil unfrozen water contents at 10 cm and 20 cm over two winters of 2011—2012 and 2012—2013

圖4 2011—2012年和2012—2013年兩個越冬期土壤10 cm和20 cm處溫度模擬值與實測值比較 Fig.4 Comparisons of simulated and measured soil temperatures at 10 cm and 20 cm across two winters of 2011—2012 and 2012—2013

2012—2013年土壤凍融過程水熱耦合運移模型的數值模擬結果具有較高精度(表2)，驗證了這種基于原位測量數據的改進型模型參數估算方法的可行性與適用性。傳統的田間取樣方法是指在觀測點周圍取土樣，然后在實驗室環境下測定水熱參數。然而，土壤的空間變異性會引起較大的誤差，取樣過程還會對土壤的結構造成破壞，影響測量精度。此外，從田間取樣到實驗室測定參數還耗費了大量的人力和時間成本。CHENG 等[12]比較了由原位測量的SFC推導出的土壤水分特征曲線(0.806≤R2≤0.994)和田間取樣后實驗室壓力盤測定的土壤水分特征曲線(0.980≤R2≤0.997)，結果表明該方法的估算精度與傳統方法相比在重疊區差異很小，證明了該方法的有效性。

2.2 蒸發量對模型預測結果的影響

為分析越冬期地表蒸發量對土壤水分運移的影響程度，在忽略蒸發量影響的假設下運用模型模擬出一組未凍水含量的變化情況。根據2012—2013年冬季凍融過程中土壤未凍水含量有無蒸發時的對比(圖3c、3d)，可以發現在不考慮蒸發時，未凍水含量和溫度的模擬值均高估，且隨著凍融過程的進行，地表蒸發量對未凍水含量的累計影響越來越明顯。在融化階段不考慮蒸發時，比考慮蒸發影響高估約0.024 m3/m3，約占土壤未凍水含量的16%。對比不同土壤深度處蒸發量的影響，發現蒸發量對于土壤凍融過程中未凍水含量的影響隨著土壤深度的增加而逐漸減小。在模型中加入地表蒸發量作為邊界條件后，模型的模擬精度普遍提高，尤其是在土壤表層。對有無蒸發條件下模擬值與實測值之間進行誤差分析(表2)可以發現：考慮蒸發量的影響后，10 cm處的模擬值的d值由0.892增加至0.924，RMSE由0.059 m3/m3降為0.046 m3/m3，而20 cm處的影響較小。因此，對越冬期小麥根區的未凍水含量的模擬應該考慮蒸發量的影響，尤其是接近土壤表層[17]。

表1 2011—2012年冬季土壤未凍水含量和溫度誤差分析 Tab.1 Analysis of soil unfrozen water content and temperature in winter of 2011—2012

越冬期土壤凍融過程中，2012—2013年的蒸發量較前一年較高(圖2)，對模型的預測精度影響較大。出現該現象的原因是該年冬季降水量相對較大，使土壤表層含水量增加，從而增加了地表蒸發量。因此，模型的預測精度也會受到氣候條件等的影響。

表2 2012—2013年冬季有無蒸發時土壤未凍水含量和溫度模擬值誤差對比分析 Tab.2 Analysis of soil unfrozen water content and temperature with or without soil evaporation in winter of 2012—2013

3 結論

(1)提出了一種土壤凍融過程水熱耦合運移參數原位估算方法，借助2011—2012年原位測定的試驗數據擬合SFC曲線，對土壤水熱參數進行最優估算，然后檢驗了該方法的適用性。結果表明：2012—2013年冬季土壤凍融過程中模型的模擬值與實測值擬合度較高，能夠準確預測越冬期小麥根區未凍水含量與溫度的變化規律，10 cm處土壤未凍水含量和溫度的d值分別達到0.924和0.918，20 cm處達到0.774和0.833。可見，該方法估算的水熱參數能夠保證模型的預測精度，從而取代田間取樣并在實驗室測定土壤參數的傳統做法，比參數辨識(反演)法所使用的模型結構簡單，不破壞土壤結構，節省了人力和時間成本。

(2)為檢驗地表蒸發量對模型精度的影響，對2012—2013年越冬期的模擬值和實測值進行誤差分析，10 cm處的未凍水含量預測值的d值由0.892增加為0.924，RMSE由0.059 m3/m3降為0.046 m3/m3。結果表明，考慮地表蒸發量對土壤凍融過程中水熱運移的影響后，土壤表層未凍水含量預測精度明顯提高。此外，隨著土壤深度的增加，蒸發量的影響逐漸減小。

(3)除土壤邊界條件和氣候因素外，SFC曲線在凍融過程的滯后效應可能是造成模擬過程中凍結階段低估和融化階段高估的重要原因。且土壤在經歷多次凍融循環后，紋理結構和孔隙分布等物理性質可能會產生改變，進而影響SFC。因此，SFC曲線及其參數的標定和修正是土壤凍融過程中進行水熱耦合運移數值模擬的重要步驟，也需要在未來的研究中找出更加適當和精確的模型標定與參數估算方法。