一種新的土壤凍融特征曲線模型

付子騰， 吳青柏， Miles DYCK， 何海龍

（1.西北農(nóng)林科技大學(xué)資源環(huán)境學(xué)院，陜西 楊凌 712100； 2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院凍土工程國家重點(diǎn)實(shí)驗(yàn)室，甘肅 蘭州 730000； 3.Department of Renewable Resources，University of Alberta，AB，T6G 2H1，Canada）

0 引言

凍土是指溫度保持在或低于水的冰點(diǎn)，并含有冰的巖石和土壤，分為短時(shí)凍土（凍結(jié)時(shí)間數(shù)小時(shí)至半月）、季節(jié)凍土（凍結(jié)時(shí)間在半月至數(shù)月）和多年凍土（凍結(jié)時(shí)間在兩年以上）［1］。在凍土中，冰和液態(tài)水（未凍水）共存，其含量隨凍融過程中溫度的變化而變化。未凍水含量和含冰量的變化顯著影響凍土的物理和機(jī)械性能，從而決定寒區(qū)工程建設(shè)（如：房屋、道路、石油管線等）的成?。?-3］。同時(shí)，凍融循環(huán)過程也會(huì)帶來如土壤凍融侵蝕、洪水、污染物的遷移、地下水上升引起的土壤鹽漬化等一系列環(huán)境問題［4-7］。因此，掌握凍融過程中土壤水熱耦合的變化機(jī)理對生產(chǎn)生活有至關(guān)重要的作用。

通常，純水的凍結(jié)發(fā)生在某一確定的溫度下［8］，具有確定的凝固點(diǎn)，但在濕潤的土壤中，即使當(dāng)溫度低于自由水的冰點(diǎn)時(shí)，土壤中所含水分也不會(huì)立刻轉(zhuǎn)化為固態(tài)，而是隨著溫度的降低表現(xiàn)出相對平滑的轉(zhuǎn)變，以曲線的形式降低［9-12］。由于土顆粒表面能的作用，土體中始終會(huì)保持一定數(shù)量的液態(tài)水被稱作未凍水［1］，當(dāng)溫度繼續(xù)降低，更多的未凍水會(huì)轉(zhuǎn)化成冰，此時(shí)毛管力和吸附力會(huì)被放大，剩余的未凍水會(huì)以薄膜水的形式存在［13］。當(dāng)溫度低于一定值時(shí)，隨著溫度的進(jìn)一步降低，未凍水含量保持不變，該溫度下的未凍水含量被稱為殘余水含量θres［14-16］。凍融過程中未凍水含量θl隨溫度T（℃）的變化曲線被定義為土壤凍融特征曲線（SFTC），代表了凍土中液態(tài)水含量和能量狀態(tài)之間的關(guān)系［17］，其中未凍水含量θl是指土壤孔隙中沒有凍結(jié)的那部分含水量（用體積含水量表示，單位為cm3·cm-3）。一般來說，決定凍土中未凍水含量大小的主要因素是溫度，其次是孔隙水壓力和土質(zhì)（包括土顆粒的礦物化學(xué)成分、分散度、含水量、密度、水溶液的成分和濃度等）［18］。溫度作為主要影響因子，直接決定凍土中未凍水含量的高低，其他因素雖不是主要影響因子，但也顯著影響著凍融特征曲線的變化。在一定的溫度范圍內(nèi)，未凍水含量隨孔隙水壓力和溶質(zhì)濃度的增大而增大；相同條件下溶質(zhì)濃度越大的土壤殘余水含量也越高，此外，土體液態(tài)水中溶質(zhì)的存在還會(huì)引起水分凝固點(diǎn)的降低［9］，從而影響土壤凍融特征曲線。土壤的質(zhì)地和結(jié)構(gòu)包括比表面積、礦物組分、土壤顆粒級配等，這些性質(zhì)會(huì)影響毛管力和吸附力的大小，進(jìn)而決定凍土中未凍水含量及殘余水含量［19］。

有關(guān)凍融特征曲線數(shù)學(xué)表達(dá)的研究已有很多，許多學(xué)者也各自通過實(shí)驗(yàn)數(shù)據(jù)和理論推導(dǎo)建立了相關(guān)的模型［15，20-34］。這些模型一般包括兩類：一類是基于土壤基本物理性質(zhì)提出的半經(jīng)驗(yàn)公式［14-15，23］，一類則是通過類比水分特征曲線而提出計(jì)算未凍水含量的公式［7，22，24-25］。這些模型雖然在各自文章中能較好地?cái)M合實(shí)驗(yàn)數(shù)據(jù)，但在實(shí)際應(yīng)用中仍然有其局限性：（1）很多模型只適用于飽和砂土，而且忽視了溶質(zhì)濃度所引起的凝固點(diǎn)降低現(xiàn)象［14-15］；（2）認(rèn)為未凍水含量的變化只依賴于溫度和土壤屬性，忽視了土壤中初始含水量的影響［34］，而Tice和Anderson［35］的研究表明，在給定的溫度下，土壤未凍水含量會(huì)隨著土壤初始含水量的增加而增加；（3）部分模型形式復(fù)雜，參數(shù)繁多且只能應(yīng)用于某種單一的土壤中，如果推廣到其他土壤中，則預(yù)測結(jié)果偏離實(shí)際值［36］。

本文回顧了現(xiàn)有的凍融特征曲線模型，基于水分特征曲線與凍融特征曲線的相似性和平行關(guān)系，提出一種新的凍融特征曲線表達(dá)式。新模型與van Genuchten［37］水分特征曲線（SWCC）相似，考慮到了初始含水量對未凍水含量的影響，包含兩個(gè)參數(shù)α、β，其中參數(shù)α和土壤自身的物理性質(zhì)有關(guān)，參數(shù)β則受初始含水量的影響。此外，本文還根據(jù)已有數(shù)據(jù)對模型作了修正，考慮到溶質(zhì)濃度對未凍水含量的影響，加入?yún)?shù)Tf，使其能適用于不同溶質(zhì)濃度的土壤中。下面將對模型的確立過程、參數(shù)擬合、以及模型的驗(yàn)證作詳細(xì)的說明，同時(shí)作模型參數(shù)的敏感性分析，以研究模型參數(shù)α和β的具體物理意義。

1 模型回顧

凍融特征曲線既能為研究土壤凍融過程和高寒地區(qū)水文模型提供理論支持，同時(shí)也是進(jìn)行寒區(qū)工程設(shè)計(jì)及施工建設(shè)的理論基礎(chǔ)，能否精確模擬某種土壤中的水熱耦合運(yùn)移對于寒區(qū)生產(chǎn)建設(shè)至關(guān)重要。許多學(xué)者已基于土壤的各種基本參數(shù)建立了相關(guān)的經(jīng)驗(yàn)?zāi)Ｐ?，包括線性、分段、指數(shù)和冪函數(shù)模型等［36，38］。但由于早期對于凍結(jié)過程中恰當(dāng)?shù)臄?shù)學(xué)表示形式存在分歧，到現(xiàn)在仍未解決，以至于回顧這些模型和經(jīng)驗(yàn)公式時(shí)，現(xiàn)有模型結(jié)構(gòu)之間由于適用條件和范圍的不同而表現(xiàn)出非常明顯的差異。

1.1 Anderson and Tice（1972）模型

土壤物理特性（如土壤干密度、土壤顆粒大小、土壤孔隙度、比表面積等）都會(huì)影響凍土中未凍水含量，進(jìn)而影響土壤凍融特征曲線的變化［37］。因此，根據(jù)土壤的物理特性，尤其是比表面積，可以建立一個(gè)簡單的模型來預(yù)測土壤中未凍水含量［28］：

式中：θl為體積未凍水含量（cm3·cm-3）；ρd和ρw分別為土壤干密度和標(biāo)準(zhǔn)狀態(tài)下未凍水的密度（g·cm-3）；S為土壤比表面積(m2·g-1)；Ac為塑性指數(shù)與土壤顆粒小于2 μm 含量所占百分?jǐn)?shù)的比值；T表示溫度，本文中全部使用攝氏溫度（℃）來表示。

Anderson［20，39］通過等溫量熱計(jì)和單位表面積歸一化來確定土壤未凍水含量，根據(jù)實(shí)際測定的未凍水含量數(shù)據(jù)總結(jié)了不同含水量與比表面積之間的關(guān)系，提出了一種更為簡單的冪函數(shù)模型：

式中：a和b為土壤特征參數(shù)與土壤比表面積有關(guān)。Kozlowski［14］給出了參數(shù)和土壤表面積之間的關(guān)系：

該式成立的隱含條件是假定在無窮低的溫度下殘余水含量為零［40］。聯(lián)立式（2）和（3）：

式（4）表明，只需測定土壤比表面積和干密度即可描述未凍水含量隨溫度的變化特征。因此，盡管冪函數(shù)模型提出已有幾十年歷史，但由于其形式簡單便于應(yīng)用，各種文獻(xiàn)及資料中仍普遍采用該方法來估算凍土中的未凍水含量［1］。

1.2 Mckenzie（2007）模型

Mckenzie 等［15］提出了一個(gè)線性分段函數(shù)來近似地描述土壤凍融特征曲線：

式中：m為線性參數(shù)，Tres為未凍水含量減少到殘余水含量的溫度，假定為-12 ℃［15］。此外，Mckenzie等［15］還給出了另一個(gè)指數(shù)函數(shù)來預(yù)測未凍水含量：

式中：n為擬合參數(shù)；T0為標(biāo)準(zhǔn)狀態(tài)下自由水的冰點(diǎn)。相比式（5）、式（6）對未凍水含量與溫度之間關(guān)系的描述更符合實(shí)際的實(shí)驗(yàn)結(jié)果［15］。

1.3 Kozlowski（2007）和Zhang（2008）模型

Kozlowski［14］對6 種 黏 土 進(jìn) 行 差 示 掃 描 量 熱 法實(shí)驗(yàn)后提出了土壤凍融特征曲線的半經(jīng)驗(yàn)?zāi)Ｐ停撃Ｐ陀扇糠纸M成：

式中：wl、wtotal、wres分別為用重量含水量表示的未凍水含量、土壤總含水量和殘余水含量；Tf代表凝固點(diǎn)的降低量。用體積含水量替換重量含水量則式（7）可變形為：

在這里Kozlowski 通過土壤塑限和含水總量來計(jì)算Tf：

式中：θp為土壤塑限（%）。Kozlowski認(rèn)為殘余水含量與土壤質(zhì)地和結(jié)構(gòu)等有關(guān)，并通過實(shí)驗(yàn)證明，殘余水含量θres和比表面積S之間存在以下關(guān)系：

實(shí)驗(yàn)結(jié)果表明，當(dāng)土壤溫度在0~-1 ℃之間時(shí)，土壤中的未凍水含量急劇下降，近似呈線性趨勢。因此，當(dāng)負(fù)溫接近于0 ℃時(shí)，未凍水含量又可以表示為土壤溫度的線性函數(shù)［41］：

1.4 Bai and Lai（2018）模型

Bai 等［16］考慮到初始含水量和殘余水含量，提出一個(gè)分段函數(shù)模型來預(yù)測未凍水含量：

式中：T0為標(biāo)準(zhǔn)狀態(tài)下自由水的冰點(diǎn)，σ為模型的擬合參數(shù)。該模型可用于不同初始含水量的土壤。

2 研究方法和材料

2.1 新的凍融特征曲線模型

在融土中，隨著土壤含水量的變化，土壤水吸力也會(huì)發(fā)生變化［42］。與凍融特征曲線的定義類似，土壤含水量與土壤水能量狀態(tài)即土壤水的基質(zhì)勢之間的曲線關(guān)系被稱為土壤水分特征曲線（SWCC）［36］，是土壤持水特性的表征。很多學(xué)者已比較了SFTC和SWCC之間的相似性［14-15，23，43-45］。早在20世紀(jì)60年代，Koopmans和Miller［17］就證實(shí)當(dāng)土壤水分完全由毛細(xì)管力或吸附力控制時(shí)，SWCC 和SFTC 之間存在平行關(guān)系；Black 和Tice［40］的研究表明，在飽和土體下，相同密度和體積的土壤的SWCC和SFTC可進(jìn)行定量比較，有學(xué)者結(jié)合Clapeyron與SWCC 導(dǎo)出了非飽和土壤中的SFTC［36］?；跉w一化的無量綱土壤含水量與土壤基質(zhì)勢之間的關(guān)系可表示為［18］：

式中：a、n、m為未確定的參數(shù)；m= 1-1/n；h為土壤基質(zhì)勢；用壓力水頭表示；Θ可由下式計(jì)算：

式中：θ、θr、θs分別為實(shí)際含水量、飽和度和殘余水含量。聯(lián)立式（13）和式（14），即可得到土壤水分特征曲線的Van Genuchten［37］表達(dá)式：

注意到土壤凍融特征曲線和水分特征曲線之間的相似性：（1）都是描述土壤液態(tài)含水量隨土壤某個(gè)物理量變化的函數(shù)；（2）長期的實(shí)驗(yàn)數(shù)據(jù)表明，土壤SFTC 的變化特征同SWCC 相似（圖1）；（3）凍結(jié)-融化過程的SFTC 和脫濕-吸濕過程中的SWCC都會(huì)產(chǎn)生滯后現(xiàn)象［42-43］；（4）SWCC 中土壤含水量居于殘余水和飽和水含量之間，類似地SFTC 中未凍水含量居于殘余水和凍土中初始含水量（或總含水量）之間。因此，歸一化理論同樣可以用于描述SFTC：

式中：θinit、θres分別為初始含水量和殘余水含量（cm3·cm-3）；θl為未凍水含量。將式（15）中的自變量h替換為溫度T，即得：

由式（16）和（17）即可導(dǎo)出形式類似于Van Genuchten［18］的歸一化公式：

式中：θinit、θres分別為初始含水量和殘余水含量(cm3·cm-3）；T為溫度（℃）；α、β為模型的擬合參數(shù)。

未凍水含量的變化是多種因素綜合作用的結(jié)果，以往的模型大都忽略了初始含水量和溶質(zhì)濃度的影響。然而研究表明，當(dāng)土壤凍結(jié)時(shí)，隨著土壤水分逐漸向冰相的轉(zhuǎn)變，土壤溶液的濃度會(huì)增加（假設(shè)冰相中不含溶質(zhì)），溶質(zhì)的存在會(huì)降低土壤水分的凝固點(diǎn)，從而影響未凍水含量［47］。此外，凍融過程中土壤的凝固點(diǎn)還會(huì)隨初始含水量的增加而升高，隨溶質(zhì)含量的增加而降低［10，21，48-50］。在式（18）中，我們考慮到土體初始含水量θinit的影響，等式中的兩個(gè)參數(shù)α、β或許能解釋溶質(zhì)濃度對未凍水含量的影響，下文將進(jìn)行參數(shù)α、β敏感性分析，以探索其具體物理意義。

2.1.1 對參數(shù)α敏感性分析

如圖2（a）顯示，若保持初始含水量不變，當(dāng)固定參數(shù)β不變時(shí)，α的變動(dòng)導(dǎo)致了整個(gè)模型曲線發(fā)生變化。保持參數(shù)β等于2.0 不變，令α分別取0.5、1.0、2.0、3.0、5.0。隨著α的增加，曲線整體下凹并向左移動(dòng)。結(jié)合圖2（b）可以發(fā)現(xiàn)，新模型曲線隨參數(shù)α的變化規(guī)律與新模型在不同溶質(zhì)濃度土壤中的曲線變化規(guī)律相似這表明參數(shù)α或許和土壤溶質(zhì)濃度有關(guān)，隨著溶質(zhì)濃度的增大，α逐漸減小。

圖2 式（18）中參數(shù)α的敏感性分析（a）；式（18）在不同溶質(zhì)濃度土壤中的曲線變化（b）Fig.2 Unfrozen water content changing with temperature for various α in the Eq.（18）（a）and for various solute concentration of the Eq.（18）（b）

2.1.2 對參數(shù)β敏感性分析

保持α固定不變，曲線隨參數(shù)β的變化如［圖3（a）］所示，保持α等于2.0 不變，令β分別取1.2、1.5、2.0、3.0、5.0，從圖3 中可以看出，隨著參數(shù)β的增加，土壤凍融特征曲線向下移動(dòng)。保持初始含水量不變，曲線的最低點(diǎn)，即殘余水含量隨β的增加而減少。結(jié)合圖3（b），在相同初始含水量下，隨著黏粒含量的減小，殘余水含量也越來越低，表明β的取值可能同土質(zhì)有關(guān)。受土壤物理性質(zhì)（如比表面積、干密度、黏粒含量等）的影響，黏粒含量越小，β越大。另一方面，參數(shù)β的增加使得凍融特征曲線在0~-2 ℃的溫度范圍內(nèi)變得越來越陡峭，當(dāng)β取到5.0 時(shí)曲線幾乎垂直。這與圖2 中曲線隨α的變化規(guī)律相似，說明參數(shù)β同樣與溶質(zhì)濃度有關(guān)。

圖3 式（18）中參數(shù)β的敏感性分析（a）；式（18）在不同質(zhì)地土壤中的曲線變化（b）Fig.3 Unfrozen water content changing with temperature for various β in the Eq.（18）（a）and unfrozen water content changing with temperature for various soils with different textures of the Eq.（18）（b）

對參數(shù)α、β的敏感性分析表明，α受溶質(zhì)濃度的影響，β則與土質(zhì)和溶質(zhì)濃度之間存在聯(lián)系，可以據(jù)此建立相關(guān)的參數(shù)轉(zhuǎn)換函數(shù)來計(jì)算不同土壤中α、β的值。

2.2 數(shù)據(jù)的來源與處理

為了驗(yàn)證新模型的可行性及研究模型中參數(shù)α、β的實(shí)際物理意義，本文引用了文獻(xiàn)中已有的試驗(yàn)數(shù)據(jù)（見文后附表1），進(jìn)行參數(shù)的擬合與模型效果的評價(jià)。在進(jìn)行數(shù)據(jù)篩選時(shí)需要考慮到：1）數(shù)據(jù)的完整性。必須要能獲取到完整凍結(jié)或融化過程的數(shù)據(jù)資料；2）數(shù)據(jù)的多樣性。所選擇的數(shù)據(jù)需要來自于盡可能多的土壤，便于模型驗(yàn)證；3）優(yōu)先選擇土壤物理性質(zhì)比較詳細(xì)的數(shù)據(jù)。如土壤比表面積、黏粒含量、容重、孔隙度等，便于模型參數(shù)物理性質(zhì)的分析?；谝陨显瓌t，本文所選土樣數(shù)據(jù)包括Smith 和Tice［51］文中的6 種粉質(zhì)土（Hot Springs 粉土、Leda 粉土、Suffield 粉質(zhì)黏土、Niagara 粉土、蘭州粉土），11 種黏土以及4 種代表性土壤礦物組分（膨潤土、高嶺石、鋰輝石、礫石）的凍結(jié)實(shí)驗(yàn)數(shù)據(jù)；Spaans 和Baker［52］文中分別處于干燥和濕潤兩種狀態(tài)下Waukegan 粉砂土的凍結(jié)試驗(yàn)數(shù)據(jù)；Suzuki［10］中濕地土壤的數(shù)據(jù)；Wen 等［53］青藏高原粉質(zhì)黏土的凍結(jié)試驗(yàn)數(shù)據(jù)；Christ 和Park［54］砂土、粉砂土、粉土的試驗(yàn)數(shù)據(jù)；Wang 等［55］中Harbin 黏土、Lanzhou 黃土以及Ma 等［56］粉土和黏土的凍結(jié)實(shí)驗(yàn)數(shù)據(jù)（數(shù)據(jù)選用的詳細(xì)信息見附表2）。

考慮到模型效果評價(jià)的合理性，將以上所有數(shù)據(jù)按照來源和土壤類型整理后分為兩組，其中一組包含70%的數(shù)據(jù)，用于模型的訓(xùn)練和參數(shù)擬合，另一組包含30%的數(shù)據(jù)，用于模型的驗(yàn)證和評價(jià)。此外，本文還選取了Patterson 和Smith［57］粉質(zhì)黏土的凍結(jié)實(shí)驗(yàn)數(shù)據(jù)以研究模型是否適用于不同溶質(zhì)濃度的土壤中。該數(shù)據(jù)包含4 個(gè)對照組，即分別用0、10、20、35 g·L-1NaCl溶液飽和，然后進(jìn)行凍結(jié)處理，得到不同溶質(zhì)濃度下的凍結(jié)試驗(yàn)數(shù)據(jù)。在本文中，同樣將其分為訓(xùn)練和驗(yàn)證數(shù)據(jù)兩部分。

2.3 參數(shù)擬合預(yù)評價(jià)指標(biāo)

本文重點(diǎn)研究新模型的可行性及適用范圍。通過對實(shí)驗(yàn)數(shù)據(jù)的分析，研究新模型能否適用于不同的飽和度和不同含水量的土壤中，檢驗(yàn)新模型在不同類型土壤中的表現(xiàn)，前文中提到的五種模型被用于模型間的比較。此外，文中還討論了模型中參數(shù)的變化對擬合結(jié)果的影響。通過MathCAD（Prime 3.0，Parametric Technology Corporation），利用非線性最小二乘算法來進(jìn)行參數(shù)α、β的擬合。

采用表1所示三種模型評價(jià)指標(biāo)對模型性能進(jìn)行評價(jià)：（1）均方根誤差（RMSE），用于衡量預(yù)測值與真實(shí)值間的偏差，反映測量的精密度。均方根誤差的值越小越好；（2）平均偏差（AD），用于來測定預(yù)測值與真實(shí)值算數(shù)平均值間的差異；（3）納什效率系數(shù)（NSE），用于模型的效率評價(jià)和對擬合結(jié)果的評判，NSE 的變化范圍從-∞~1，其值越接近1，則模型的可信度越高，預(yù)測值越靠近真實(shí)值。此外，文中還給出了各模型在部分土壤中預(yù)測值與實(shí)測值的1∶1線性圖以便更直觀地展示測試數(shù)據(jù)和模型預(yù)測值之間的差異。

表1 復(fù)合介質(zhì)混合模型的性能評價(jià)指標(biāo)Table 1 Performance evaluation indexes of the composite medium hybrid model

3 模型評價(jià)

3.1 新模型在不同類型土壤中的表現(xiàn)

使用文獻(xiàn)中已有的實(shí)驗(yàn)數(shù)據(jù)擬合參數(shù)α、β以研究新模型在不同的類型土壤中的適用情況。圖4展示了新模型在7 種代表性土壤中的擬合結(jié)果，其中粉土和粉砂土是Christ 和Park［54］的野外實(shí)驗(yàn)數(shù)據(jù)；Harbin 黏土來自Wang 等［55］，其余四種土樣均來自文獻(xiàn)Patterson［58］。從圖4可以看出，對于這7種土壤樣品，用新模型所得到的擬合曲線結(jié)果與實(shí)測值較為吻合，納什效率系數(shù)（NSE）都在0.95 以上，這表明新的模型可廣泛適用于不同類型土壤中。

圖4 不同土壤中未凍水含量與隨溫度變化的擬合曲線，其中散點(diǎn)為實(shí)測數(shù)據(jù)，實(shí)線為式（18）計(jì)算的模擬值Fig.4 Unfrozen water contents of various soils changing with temperature，where the dots are measured data and lines are calculated value by the Eq.（18）

表2 給出了相關(guān)的擬合指標(biāo)，均方根誤差（AD）表明，在大部分土樣中，預(yù)測值稍低于實(shí)測值，但在統(tǒng)計(jì)學(xué)意義上誤差均在可接受范圍內(nèi)。盡管不同土壤中模型參數(shù)有所差異，但對于給定的土壤，只要能通過實(shí)驗(yàn)數(shù)據(jù)得到參數(shù)α、β的值，新模型就能準(zhǔn)確預(yù)測相同條件下該種土的土壤凍融特征曲線。

3.2 新模型不同初始含水量下的表現(xiàn)

為了研究新模型在不同初始含水量土壤中的表現(xiàn)，本文選用了Wang等［53］的實(shí)驗(yàn)數(shù)據(jù)。土樣是來自青藏高原的粉質(zhì)黏土，比表面積為8.12 m2·g-1，干密度1.56 cm3·cm-3，保持其他條件不變，設(shè)置土壤初始 含 水 量 分 別 為0.254、0.313、0.357、0.398、0.459 cm3·cm-3，分別擬合出每個(gè)初始含水量下的參數(shù)α、β的值，然后通過該參數(shù)去計(jì)算所有初始含水量下的凍融特征曲線，觀察其差異性。從圖5 中可以看出，同一土壤分別在五種初始含水量下得到的模型參數(shù)對整體的擬合結(jié)果并顯示出沒有出現(xiàn)太大的差異，但相對來說，在初始含水量較低的條件下得到的模型參數(shù)在擬合整體時(shí)會(huì)更準(zhǔn)確一些。表3 顯示不同初始含水量下得到的參數(shù)α、β值不同但差異較小，而且隨著初始含水量的變化β和α并沒有表現(xiàn)出明顯的線性的變化規(guī)律?？傮w來看，預(yù)測結(jié)果較好地吻合了實(shí)測值。表3給出的相關(guān)擬合指標(biāo)顯示，模型的NSE都在0.96以上，也就是說，在任意土壤類型和初始含水量下得到的模型參數(shù)都是有效的。

圖5 不同初始含水量下新模型的表現(xiàn)，散點(diǎn)為實(shí)測數(shù)據(jù)，實(shí)線為新模型的模擬結(jié)果Fig.5 Unfrozen water content changing with temperature for various θinit（Note：the scattered points are observed and the solid lines are simulated of the new model）

表2 新模型在7種土壤樣本中的參數(shù)α、β 及擬合指標(biāo)Table 2 A summary of the parameters α、β for the seven soils and performance measures of the new model to fit the 7 soils samples

3.3 新模型在不同溶質(zhì)濃度下的表現(xiàn)

本文選擇Patterson 和Smith［57］粉質(zhì)黏土來進(jìn)行模型在不同溶質(zhì)濃度條件下的檢驗(yàn)，土樣干密度為1.19 g·cm-3，保持初始含水量0.523 cm3·cm-3不變，取溶質(zhì)濃度為10、20、35 g·L-1NaCl 溶液，設(shè)置0 g·L-1NaCl 溶液作為對照。圖6 為式（18）在不同濃度溶液中的表現(xiàn)，從圖中可以看出，模型在低濃度時(shí)的表現(xiàn)依舊能很好地同實(shí)測值吻合，但隨著溶質(zhì)濃度的增加，擬合曲線逐漸偏離實(shí)測值，且濃度越大，預(yù)測值與實(shí)測值的差異越大，可能是因?yàn)樵摰仁酵瑯記]有包含溶質(zhì)濃度對土壤水分凝固點(diǎn)降低的影響，考慮到這一點(diǎn)，引入一個(gè)新的參數(shù)Tf，式（18）可變?yōu)椋?/p>

圖6 式（18）在不同濃度溶液中的擬合情況，散點(diǎn)為實(shí)測值，實(shí)線為新模型模擬結(jié)果Fig.6 Unfrozen water content changing with temperature of different solute concentrations of the Eq.（18）（Note：the scattered points are observed，the solid lines are simulated by the new model）

Tf表示土壤水分凝固點(diǎn)的降低，在模型中有兩種含義：（1）類似于α、β，看作一種擬合參數(shù)，沒有實(shí)際的物理意義；（2）表示土壤水分凝固點(diǎn)隨溶質(zhì)濃度的降低，具有實(shí)際的物理意義，并且有相應(yīng)的計(jì)算公式。本文考慮到兩種情況，對比作為參數(shù)擬合得到的Tf和作為變量計(jì)算得到的Tf，選出最適合的模型形式。

表3 五個(gè)初始含水量下分別得到的參數(shù)α、β值及其在不同初始含水量情況下的統(tǒng)計(jì)指標(biāo)（加黑數(shù)值表示該含水量下對應(yīng)的指標(biāo)值）Table 3 The values of parameter α and β corresponding to the five initial water content θinit and their statistical indexes（Note：the corresponding index under initial moisture content θinit is indicated in bold type）

計(jì)算Tf的經(jīng)驗(yàn)公式有很多，Banin 和Anderson［49］描述了土壤中溶質(zhì)的存在會(huì)降低冰點(diǎn)溫度，隨著土壤水分的凍結(jié)，溶質(zhì)逐漸被滯留在越來越小的空間中，其中，凝固點(diǎn)的降低Tf與溶質(zhì)濃度之間存在以下關(guān)系［59］：

式中：Tref為參考溫度，與溶質(zhì)種類有關(guān)，NaCl 可取62 ℃［31］。Sn為溶質(zhì)濃度（g·L-1NaCl）。Kozlowski［9］修改了式（9），給出了另一種計(jì)算Tf的方法：

式中：w為重量含水量（%）；wL為土壤液限（%）。此后，Kozlowski［60］又提出可根據(jù)土壤干密度和初始含水量來求得Tf：

式中：α和β為給定土壤的系數(shù)，是一個(gè)常數(shù)。Bodnar［61］根據(jù)試驗(yàn)數(shù)據(jù)，回歸出了一個(gè)三次多項(xiàng)式來描述Tf與溶質(zhì)濃度的關(guān)系：

用上述四種方法計(jì)算Tf，并分別與新模型耦合。表4顯示四種經(jīng)驗(yàn)公式的擬合結(jié)果在統(tǒng)計(jì)學(xué)意義上相差不大，考慮到經(jīng)驗(yàn)公式的不確定性，新模型采用第一種方法［式（21）］來計(jì)算Tf值。圖7顯示在加入Tf后，低濃度條件下，擬合曲線并未有太大變動(dòng)，但在溶質(zhì)濃度較高時(shí)，式（19）的擬合結(jié)果明顯要優(yōu)于式（18），而且相比于用經(jīng)驗(yàn)公式計(jì)算的Tf值，參數(shù)擬合得到的結(jié)果要更加精確（表5），而濃度為35 g·L-1NaCl溶液的樣本在三種方法的計(jì)算下NSE都沒有超過0.9，這可能是由于在高濃度條件下，實(shí)測數(shù)據(jù)的不準(zhǔn)確所導(dǎo)致。

圖7 加入Tf后的新模型在不同濃度溶液中的擬合情況，Tf用經(jīng)驗(yàn)公式和參數(shù)擬合兩種方法獲取。其中實(shí)線未考慮Tf影響，短劃線Tf通過經(jīng)驗(yàn)公式計(jì)算，虛線的Tf通過擬合獲得Fig.7 Unfrozen water content changing with temperature while considering effects of Tf under different solute concentrations（Note：solid line is not considered the influence of Tf，short line Tf is calculated by the empirical formula and the dashed line Tf is obtained by fitting）

綜合來看，使用式（18）和式（19），均能準(zhǔn)確地描述不同溶質(zhì)濃度下的未凍水含量隨溫度變化的關(guān)系，但式（19）能更好地反映出溶質(zhì)濃度的變化對土壤凝固點(diǎn)降低的影響。因此，在實(shí)際應(yīng)用過程中，當(dāng)土壤溶質(zhì)濃度較低時(shí)，可直接通過式（18）得到準(zhǔn)確的凍融特征曲線，而當(dāng)土壤溶質(zhì)濃度較高時(shí)，通過給模型加入Tf［式（19）］可有效提高其擬合精度。

表4 四種經(jīng)驗(yàn)公式計(jì)算的Tf值在不同土壤溶質(zhì)濃度下的表現(xiàn)，其中參數(shù)α、β通過公式（19）擬合得到Table 4 The statistical indexes of the four empirical formulas calculating the Tf at different soil solute concentrations，of which the parameters α，β are fitted by the Eq.（19）

表5 式（18）和式（19）在不同溶質(zhì)濃度土壤中的擬合指標(biāo)，式（19）中Tf的值分別用經(jīng)驗(yàn)公式（20）計(jì)算和參數(shù)擬合Table 5 The statistical indexes of Eq.（18）and Eq.（19）in different soil solute concentrations，the values of Tf in Eq.（19）are calculated by the Eq.（20）and by parameter fitting，respectively

3.4 模型間的對比

本文選擇了前文中提到五種較為典型的模型與新模型對比，包括兩種線性模型［15-41］，兩種分段函數(shù)模型［14］以及早先提出的指數(shù)模型［20］。四種土壤（包括砂土、粉砂土、粉土、哈爾濱黏土）被用于對比，砂土、粉砂土、粉土以及一種哈爾濱黏土。其中粉土、粉砂土及砂土凍結(jié)數(shù)據(jù)均來Christ 和Park［54］，凍結(jié)過程未凍水含量的測定方法為時(shí)域反射儀（TDR）法；哈爾濱黏土數(shù)據(jù)來自Wang等［55］，測量方法為核磁共振法（NMR）。

圖8 和表6 顯示了各模型在四種土壤中的擬合結(jié)果。在砂土中，除了Zhang等［41］的線性模型外，其余五種包括新模型都能較為準(zhǔn)確地描述未凍水含量和溫度之間的關(guān)系。在其他三種土壤中，新模型的擬合精度最高，Bai 等［16］的模型次之，Anderson等［20］的預(yù)測結(jié)果與實(shí)測值偏離最大，預(yù)測值往往低于測值，且隨著土壤黏粒成分的增加和土壤顆粒半徑的減小，偏離值增大，可能是由于該模型忽視了初始含水量，只考慮到了土壤比表面積的影響所導(dǎo)致的。Zhang 的線性模型考慮到了初始含水量，其曲線的斜率與初始含水量和殘余水含量有關(guān)，因此相對于觀測值，該模型在未凍水含量接近殘余水含量之前總是顯示出過高地估計(jì)?？傮w上看，除了新模型和Bai 的模型外，其余模型在應(yīng)用到砂土以外的土壤中時(shí)均會(huì)較大程度地偏離真實(shí)值，而且隨著土壤顆粒級配的變化這種差異也會(huì)越發(fā)顯著。此外，部分土壤樣品被用來擬合模型的1∶1散點(diǎn)圖（圖9），研究預(yù)測值和實(shí)測值之間的差異。從圖9 可以看出，在所有的土壤樣品中，新模型的預(yù)測值都非常接近真實(shí)數(shù)據(jù)。其他模型中，Bai 模型的表現(xiàn)良好，預(yù)測值與真實(shí)值也較為接近；Kozlowski 模型由于可用數(shù)據(jù)太少，無法得出準(zhǔn)確結(jié)論，但從圖中擬合結(jié)果來看，在砂土和粉砂土中表現(xiàn)較好。Anderson 模型在未凍水含量較低時(shí)表現(xiàn)較好，在較高的未凍水含量下則出現(xiàn)明顯的偏離。

圖8 新模型與另外五種模型分別在砂土、粉砂土、粉土及黏土中的表現(xiàn)Fig.8 Unfrozen water content changing with temperature for sand（a），silt-sand（b），silt（c）and Harbin clay（d），measured and simulated by the new model and other five quoted models

圖9 6種模型在部分土壤數(shù)據(jù)的預(yù)測值和實(shí)測值1∶1散點(diǎn)圖Fig.9 The 1∶1 scatter plots of unfrozen water content simulated by the new model and the other five quoted models against the measured

4 結(jié)論

本文考慮到初始含水量和溶質(zhì)濃度對凍融特征曲線的影響，類比水分特征曲線模型，提出一種包含兩個(gè)參數(shù)的新凍融特征曲線模型。利用7種典型土壤樣本對模型適用性進(jìn)行了驗(yàn)證，結(jié)果表明新模型可適用于任意類型和飽和度的土壤中。通過加入新的參數(shù)Tf，可在統(tǒng)計(jì)學(xué)上顯著提高模型在高溶質(zhì)濃度條件下的預(yù)測精度。此外，明確模型參數(shù)α和β在模型中的物理意義有助于模型的完善和應(yīng)用，敏感性分析表明，參數(shù)α與土壤溶質(zhì)濃度有關(guān)；參數(shù)β與土壤溶質(zhì)濃度和土壤自身性質(zhì)如比表面積、干密度和黏粒含量有關(guān)，可根據(jù)這些影響因子可建立參數(shù)轉(zhuǎn)換函數(shù)，方便模型在實(shí)際中的應(yīng)用。由于收集數(shù)據(jù)的不完整，本文未能建立有效的參數(shù)轉(zhuǎn)換函數(shù)，但給出了用于測試模型的38種土壤樣品的參數(shù)值（附表1）以供參考。

表6 新模型和其他五種模型的擬合指標(biāo)Table 6 The statistical indexes of the new model and the other five quoted models

附表1 38種土壤樣品的參數(shù)α和β值及模型擬合指標(biāo)Attached Table 1 The parameters α and β of the 38 soil samples，together with determination method and index

附表2 土樣數(shù)據(jù)的來源、測定方法和基本性質(zhì)Attached Table 2 Sources，determination methods and basic properties of the soil samples