水熱耦合模擬中的凍土滲透系數模型選擇

王青林， 陳 磊， 明 鋒

（1. 青海省湟源公路工程建設有限公司， 青海 西寧 810000； 2. 青島農業(yè)大學 建筑工程學院，山東 青島 266000；3. 中國科學院 西北生態(tài)環(huán)境資源研究院 凍土工程國家重點實驗室， 甘肅 蘭州 730000）

0 引言

我國是第三凍土大國，凍土區(qū)面積約占全國面積的53.5%［1］。受凍土區(qū)氣溫變化的影響，不僅土中原位水發(fā)生相變產生約9%的體積膨脹，而且未凍結土層中的孔隙水在溫度梯度作用下將向凍結鋒面遷移并凍結成冰。由此將引起地表出現(xiàn)大幅度隆起，局部變形量甚至超過40 cm，這將導致修建在寒區(qū)環(huán)境中的建筑物產生凍害［2-4］。為解決因凍脹引起的凍害問題，研究者提出了諸多水熱耦合模型來評估凍脹［5-7］。在這些凍脹模型中，凍土滲透系數是控制水分遷移速率及遷移量的關鍵參數。因此，選擇合理的滲透系數模型就成為準確預測凍脹的前提條件。

獲取凍土滲透系數的方法有兩種：試驗測試和理論推導。盡管飽和凍土滲透系數試驗存在費時、費力、離散性大等缺陷，但諸多研究人員仍得到了凍土滲透系數隨溫度的變化規(guī)律［8-12］?；谝延性囼灲Y果，國內外學者提出了一系列經驗模型和理論模型來彌補試驗方法的不足。經驗模型的特點是形式簡單，但沒有嚴格的推導過程，僅僅是根據試驗結果所給出擬合公式，通常具有一個或多個經驗參數［13］。而且，經驗模型中凍土滲透性系數是未凍水含量或者溫度的函數。在水熱耦合模型中，大多采用經驗模型來描述凍土滲透系數變化規(guī)律，如應用較多的冰阻模型［7］。該模型認為飽和凍土的滲透系數由于受到了冰的阻力作用，致使其降低至相同液態(tài)水含量條件下非飽和未凍土滲透系數的1/I（I是冰阻系數，其常為含冰量或溫度的函數）。經驗模型中的參數是通過試驗結果擬合得到，所以這些參數僅適用于特定土體。當將其用于其他土體時，模型有效性將有所降低［14］。所以，在不同水熱模型中，通常會選擇不同的經驗模型來提高預測精度。相對而言，理論模型具有嚴格的推導過程，其模型參數具有明確的物理意義。由于理論模型減少了對經驗系數的依賴，其普適性更強，更值得在水熱模型中推廣使用。因此，本文對凍土滲透系數理論模型在水熱耦合模型中的有效性進行討論。

由于水熱耦合模擬結果對凍土滲透系數極為敏感，所以合理選擇凍土滲透模型就變得十分重要。為評價理論模型在水熱耦合模型中的適用性，本文以四組飽和凍土的滲透系數理論模型為例，首先將滲透系數理論模型與經驗模型預測結果進行對比，然后選擇預測效果較好的滲透系數模型帶入水熱耦合模型，以凍脹量變化為判據評價其適用性。

1 凍土滲透系數理論模型

在凍土滲透實驗的基礎上，人們提出了諸多經驗模型及理論模型，表1 列舉了四組飽和凍土的滲透系數理論模型。根據模型特點，將四組模型分別命名為分形模型［13］、離散統(tǒng)計模型［14］、連續(xù)統(tǒng)計模型［15］和等價模型［16］，各模型的數學表達式及主要計算參數如表1所示。

表1 凍土滲透系數模型Table 1 Hydraulic conductivity models for frozen soil

上述四組模型的詳細推導過程，可參閱文獻［13-16］。為更好體現(xiàn)四組滲透系數模型的優(yōu)缺點，下面對四組理論模型進行簡單介紹。分形模型是將分形理論應用到凍土中，可更好地表征土中孔隙真實無序的狀態(tài)。該模型結果主要取決于凍土中的最大孔徑，孔隙率、孔徑分維數和迂曲度，其中最大孔徑對凍土的滲透系數影響最大［13］。離散統(tǒng)計模型是基于土體凍結特征曲線與土體孔隙分布函數的相似性，以及非飽和土的滲透系數與土水特征曲線之間的關系而提出的，該模型的計算結果主要取決于其中有效孔隙分布［14］。連續(xù)統(tǒng)計模型則是通過建立凍結特征曲線和孔徑分布函數的關系函數，提出了滲透系數模型，該模型只需要對不同未凍水含量條件下的土體凍結特征曲線進行積分即可得到滲透系數［15］。等價模型借鑒非飽和土的滲透系數模型，認為孔隙冰與空氣對水分遷移的阻礙作用是相當的，認為飽和凍土的滲透系數等于相同液態(tài)水含量條件下非飽和未凍土的滲透系數［16］。由于部分文獻中并沒有給出飽和未凍土的滲透系數，所以只能計算其相對滲透系數進行比較。當采用相對滲透系數表述預測結果時，只需要土體凍結特征曲線和土體顆粒級配曲線。

為量化滲透系數預測值和實測值之間的差異，采用均方根誤差（RMSE）進行分析，其計算公式如下：

式中：N是實測數據數量；Krti是第i個實測相對滲透系數（由經驗模型確定）；Krpi是第i個預測相對滲透系數（由理論模型確定）。

2 結果與分析

為評價滲透系數理論模型的適用性，此處選取文獻［17-19］中的三個水熱耦合模型為例進行分析。在這三組水熱耦合模型中，不僅都采用的是飽和凍土的滲透系數經驗模型，而且文獻中都給出了詳細的模型參數。模型驗證從兩個方面開展：滲透系數變化和凍脹變化。為方便比較，將案例一、案例二和案例三中的滲透系數經驗模型分別命名為經驗模型一，經驗模型二和經驗模型三。

2.1 滲透系數變化

2.1.1 案例一

案例一為文獻［17］中的水熱耦合模型，該模型采用經驗模型一來描述凍土滲透系數變化。經驗模型一可以表示為［20］：

式中：k為飽和凍土的滲透系數；ks為飽和未凍土的滲透系數；θu為體積未凍水含量；θ為初始飽和體積含水量；b為與土顆粒級配相關的擬合參數，可以表示為：

式中：dg為幾何平均顆粒直徑；σg為幾何標準差；dg與σg可由土顆粒級配曲線獲得。

基于文獻［17］中給出的土體凍結特征曲線、土顆粒級配曲線和飽和滲透系數，利用表1中的模型，可計算得到四組模型計算的相對滲透系數（圖1）。從圖1可以看出本文所列舉的四組理論模型可以有效反映滲透系數隨溫度的變化規(guī)律，且預測值與文獻［17］中的滲透系數值比較接近。應該指出的是，分形模型存在有效的臨界溫度，對應凍結過程中的粉質黏土、粉土和粉砂分別是-18.0 ℃、-15.0 ℃和-13.7 ℃。當溫度低于該臨界溫度時，分形模型將不再適用。從圖1 中的計算結果來看，該臨界溫度較低，可滿足凍土滲透系數計算要求。離散統(tǒng)計模型高估了粉質黏土的滲透系數，但低估了粉砂土的滲透系數。整體而言，分形模型、等價模型和連續(xù)統(tǒng)計模型的表現(xiàn)較好。

圖1 不同土體的相對滲透系數隨溫度變化Fig. 1 Variation of the relative hydraulic conductivity with the temperature

圖2 進一步在雙對數坐標中對比了滲透系數模型計算值和文獻［17］中所用的滲透系數值。在《土工試驗方法標準》（GB/T 50123—2019）中，土體滲透系數的最大試驗允許差值為±2×10-n（n為自然數）。參考此誤差標準，本文確定相對滲透系數預測值的可接受試驗誤差值為±2×10-n。如圖2 所示，大多數的數據點都落在可接受試驗誤差區(qū)域內，僅少部分數據點落在可接受實驗誤差區(qū)域外。需要注意的是，有三組滲透系數數據點甚至落在一個數量級誤差區(qū)域外，如粉砂的分形模型，粉土和粉砂的離散統(tǒng)計模型?？傮w而言，本文所列四組滲透系數模型預測值和文獻［17］水熱耦合模型中所用的滲透系數值較接近。這表明本文所列四組滲透系數理論模型可代替文獻［17］中的經驗模型。

圖2 文獻［17］所用相對滲透系數與四組理論模型計算值對比Fig. 2 Comparison of the relative hydraulic conductivity used in Reference ［17］ and the values calculated by the four hydraulic conductivity models

表2匯總了這四組模型的均方根誤差值。可以看出，四組模型的均方根誤差值均較小，說明了四組模型具有良好的穩(wěn)定性。四種模型均對粉砂的預測效果較差，對粉土和粉質黏土的效果較為準確。等價模型的均方根誤差最小，其次是連續(xù)統(tǒng)計模型、分形模型和離散統(tǒng)計模型。對于文獻［17］中的土樣，等價模型的預測能力最好。

表2 模型的均方根Table 2 RMSE of the model

2.1.2 案例二

案例二以文獻［18］中的水熱耦合模型為例進行分析，并以粉質黏土為研究對象。文獻［18］中使用的滲透系數為經驗模型二，其表達式為：

式中：模型所需參數可查詢參考文獻［18］。

圖3對比了四組理論模型預測值與文獻［18］所用的相對滲透系數。如圖所示，相對滲透系數經驗模型值與分形模型和離散統(tǒng)計模型的預測值較為接近，但是與連續(xù)統(tǒng)計模型和等價模型的預測值稍有偏差。鑒于該經驗模型預測結果位于本文四組模型預測范圍內，且被O’Neill 等［21］、周國慶等［22］、周揚等［23］多次采用，表明該模型具有一定的合理性。需要指出的是，該模型的指數對不同土并不是一個常數。對于等價模型，根據顆粒級配曲線計算得到其指數n=12.08，與文獻［18］所用的n=9 有所差別，所以導致相對滲透系數預測值偏低。對分形模型，該粉質黏土分形模型有效的臨界溫度是-22 ℃，低于該溫度則分形模型不再適用?？梢钥闯觯馏w顆粒級配曲線不僅嚴重影響等價模型中的指數n的取值，而且對土體凍結特征曲線的數值和形態(tài)也有較大影響，進而影響分形模型、離散統(tǒng)計模型和連續(xù)統(tǒng)計模型的預測精度。

圖4 在雙對數坐標中對比了四組理論模型與經驗模型二的滲透系數預測結果。從圖4 中可以看出，只有分形模型和離散統(tǒng)計模型的數據點落在了可接受試驗誤差區(qū)域內，等價模型的數據點則大多落在了可接受實驗誤差區(qū)域的邊緣。連續(xù)統(tǒng)計模型的數據點雖然大多落在了可接受實驗誤差區(qū)域外，但也基本分布在一個量級的誤差區(qū)域內。因此，可以采用分形模型和離散統(tǒng)計模型來代替經驗模型二。

圖4 文獻［18］中所用的相對滲透系數與四組理論模型計算值對比Fig. 4 Comparison of the relative hydraulic conductivity used in Reference ［18］ and the values calculated by the four hydraulic conductivity models

表3 匯總了這四組模型的均方根誤差值，可以看到四組模型的均方根誤差值均較小，說明四組模型的穩(wěn)定性較好。相對而言，分形模型和離散統(tǒng)計模型的均方根誤差較小，表明這兩組模型可以代替經驗模型二。

表3 模型的均方根Table 3 RMSE of the model

2.1.3 案例三

案例三以文獻［19］的水熱耦合模型為例進行分析。該耦合模型采用的是經驗模型三，其表達式為：

從公式（5）可以看出，經驗模型三的表達式簡單，只需要溫度和飽和滲透系數就可以預測凍土滲透系數。圖5給出了本文所列舉的滲透系數模型計算的相對滲透系數。如圖所示，除等價模型外，其余三組理論模型與經驗模型三的滲透系數預測值較為接近。需要指出的是，該粉質黏土對應的分形模型有效臨界溫度為-1.40 ℃。對比前兩個案例，此粉質黏土的分形有效臨界溫度偏高，這是因為該土樣的凍結特征曲線斜率過大，使得大部分水分被迅速凍結。等價模型與經驗模型三的預測結果相差最大，可能是由于土體顆粒級配曲線不合理，導致計算得到的指數n偏離極大所致。盡管經驗模型三存在經驗參數，但仍適用于描述此粉質黏土樣品滲透系數變化。

圖5 相對滲透系數隨溫度的變化曲線Fig. 5 Variation of the relative hydraulic conductivity with the temperature

圖6 進一步在雙對數坐標中對比了本文所列舉的四組滲透系數模型和經驗模型三計算的相對滲透系數數值。如圖6 所示，連續(xù)統(tǒng)計模型的數據點基本都落在可接受試驗誤差范圍內；分形模型的數據點也大多落在可接受試驗誤差區(qū)域內；離散統(tǒng)計模型只有少部分數據點落在可接受試驗誤差內，但是大部分都落在了一個量級的誤差區(qū)域內；等價模型預測結果則幾乎全部不在可接受試驗誤差區(qū)域內。這說明等價模型并不適用于此土樣，而其余三組理論模型可以代替經驗模型三建立凍土滲透系數變化趨勢。

圖6 文獻［19］中所用的相對滲透系數與四組模型計算值對比Fig. 6 Comparison of the relative hydraulic conductivity used in Reference ［19］ and the values calculated by the four hydraulic conductivity models

表4匯總了這四組模型的均方根誤差值。由表4 可見，除了等價模型，其余三組模型的均方根誤差值均較小，說明了這三組模型的可靠性。其中，連續(xù)統(tǒng)計模型的均方根誤差最小，其次為分形模型和離散統(tǒng)計模型，而等價模型的均方根誤差最大。因此，可以采用分形模型和連續(xù)統(tǒng)計模型預測該凍結粉質黏土樣品的滲透系數。

表4 模型的均方根Table 4 RMSE of the model

2.2 凍脹變化規(guī)律

凍土滲透系數模型大多應用于預測凍土中的水熱分布及凍脹、融沉變形。從前一節(jié)的分析結果來看，分形模型在三個案例中的適用性更強。因此，本節(jié)采用分形模型代替文獻［18］中的經驗模型二，并利用原文中的相關參數對凍脹試驗結果進行預測（圖7）。如圖7 所示，在凍結前期（0～50 h），無論是經驗模型二還是分形模型，其預測的凍脹量與實測凍脹量相差很小。在凍結后期（50～90 h），經驗模型二的預測結果偏小，而采用分形模型的預測結果偏大。由于經驗模型二的滲透系數預測值略小于分形模型的預測值，所以采用分形模型的凍脹預測結果偏大。但從凍脹預測誤差來看，分形模型的預測結果接近實測值（最大相對誤差小于10%），可證明分形模型的有效性。

圖7 不同滲透系數模型下的凍脹量預測結果Fig. 7 Prediction of the frost heave with different hydraulic conductivity models

3 討論

3.1 理論模型的優(yōu)缺點

從三個案例來看，本文所列四組理論模型可以有效的反映文獻中滲透系數的變化規(guī)律，而且凍脹算例也證實了理論模型的可靠性。由于選擇不合理的凍土滲透系數模型，將導致水熱分布預測結果不準確，進而影響凍脹預測精度。因此，若在水熱耦合模型中，采用滲透系數理論模型，不僅可有效避免經驗公式參數的選取，而且還能保證預測精度。但需要注意的是，滲透系數理論模型也存在一定的局限性：（1）最大孔徑對分形模型的預測結果影響較大，且該模型存在溫度下限，但該下限可滿足凍脹模型的要求；（2）離散統(tǒng)計模型雖然解決了分形模型的缺點，但其計算結果并不是一條平滑曲線，難以通過一個公式擬合，因此很難應用到數值計算中。（3）連續(xù)統(tǒng)計模型所需參數較少，但對部分土體的適用性較差。（4）等價模型不僅沒有嚴格的理論推導，而且對部分土體的預測效果較差。從四組模型對三個案例的預測結果來看，不同理論模型表現(xiàn)出不同的預測效果，而分形模型在三個案例中表現(xiàn)較為均衡，其適用性最強。因此，建議在水熱耦合模型中采用分形模型來描述凍土滲透系數變化。

3.2 凍融過程中的滲透系數滯后性

從表1 中可以看出，凍土滲透系數不僅可以采用溫度為變量，而且也可以采用未凍水含量為變量。圖8 以經驗模型一（Tarnawski 和Wagner 模型）為基礎，計算得到凍結與融化過程中的土體相對滲透系數變化規(guī)律［20］。如圖8 所示，在相同溫度條件下，計算得到的粉質黏土的相對滲透系數最大，其次是粉土，粉砂的相對滲透系數最小；而在相同體積未凍水含量下，正好相反，粉砂的相對滲透系數最大，其次是粉土，而粉質黏土的相對滲透系數最小。這是因為Tarnawski 和Wagner 模型中，土體凍結特征曲線和土體顆粒級配曲線共同決定了滲透系數的大小。土體顆粒級配曲線不僅決定了模型參數2b+ 3 的大小，而且也影響土體中的未凍水含量θu。依據土體顆粒級配曲線，計算得到粉質黏土，粉土和粉砂的模型參數2b+ 3 分別是12.22、10.83和8.52。如圖8（a）所示，計算出來的融化過程的滲透系數要小于凍結過程的滲透系數，這是因為土體在融化過程中的未凍水含量要小于其在凍結過程中的未凍水含量。但是隨著溫度的降低，滯后效應對滲透系數的影響也在逐漸減小。但是采用未凍水含量為變量時，滲透系數在凍融過程并沒有體現(xiàn)出滯后效應［圖8（b）］，這和以往的研究結果一致［24］。因此，建議以未凍水含量為變量表示凍土滲透系數模型。

圖8 凍結和融化過程中相對滲透系數隨溫度（a）和體積未凍水含量（b）的變化曲線Fig. 8 Variation of the relative hydraulic conductivity with the temperature （a） and unfrozen water content （b） during the freezing and thawing process

4 結論

本文通過對四組飽和凍土的滲透系數模型和三個水熱耦合模型中選取的經驗模型的計算結果進行對比，討論了四組飽和凍土的滲透系數模型的預測能力，得出以下主要結論：

（1）水熱耦合模型中所用的滲透系數經驗模型數值和本文所列舉的四組理論模型數值大多都非常接近，建議使用滲透系數理論模型來描述凍土滲透系數變化規(guī)律。

（2）四組理論模型均能有效描述凍土滲透系數變化，但理論模型預測結果受土體顆粒級配曲線影響明顯，推薦使用分形模型計算凍土滲透系數。

（3）用溫度表示的滲透系數表現(xiàn)出滯后效應，而用未凍水含量表示的滲透系數并沒有表現(xiàn)出滯后效應，建議在凍土滲透系數模型中以未凍水含量為變量。