土中水膜厚度變化規(guī)律及未凍水含量預(yù)測方法

萬旭升， 顏夢宇， 路建國， 晏忠瑞

（西南石油大學(xué)土木工程與測繪學(xué)院，四川 成都 610500）

我國凍土面積分布廣泛，多年凍土和季節(jié)性凍土面積之和超過國土面積的70%［1］。寒區(qū)工程建設(shè)面臨因水熱相互作用引發(fā)的凍脹融沉現(xiàn)象，嚴(yán)重制約著公路、鐵路以及構(gòu)筑物等基礎(chǔ)設(shè)施建設(shè)。水分是土體發(fā)生凍融的物質(zhì)基礎(chǔ)，未凍水含量的改變直接影響凍土的熱力學(xué)及變形特性。同時，未凍水含量隨溫度變化規(guī)律也是水熱數(shù)值計算的必要條件。

已有試驗(yàn)研究表明凍土中的未凍水主要受土質(zhì)、外界條件（包括溫度和壓力）、鹽分以及凍融過程等因素的影響［2-5］。水分凍結(jié)隨著土體顆粒粒徑減小、含鹽量增加、液體飽和度減小、土顆粒活化比表面積增大而難度增加［1］。隨著測試技術(shù)的不斷發(fā)展，核磁共振法、時域反射法和頻域反射法等［6-9］廣泛運(yùn)用在未凍水含量測試中，通過信號值與溫度的對應(yīng)關(guān)系回歸不同溫度下的未凍水含量值，建立未凍水含量-溫度關(guān)系。Anderson 和Tice［10］通過大量室內(nèi)試驗(yàn)，依據(jù)試驗(yàn)數(shù)據(jù)首次發(fā)現(xiàn)未凍水含量隨溫度降低呈現(xiàn)指數(shù)函數(shù)減小，并給出了經(jīng)驗(yàn)公式。隨后，此種關(guān)系得到了進(jìn)一步改進(jìn)及發(fā)展［11-12］。Clapeyron 方程建立了冰、水壓力與溫度之間的聯(lián)系，被廣泛運(yùn)用到水熱模型計算中［13-15］，該方程闡明了水分凍結(jié)主要取決于冰水壓力的變化。Shoop 和Bigl［16］通過土水特征曲線獲取了凍融過程中的水分凍結(jié)曲線并預(yù)測了未凍水含量隨溫度變化規(guī)律。Zhang 等［17］考慮了土孔隙中冰對水的吸附作用，改進(jìn)了一般的土體凍結(jié)曲線?；诿?xì)管壓力理論，Dall’Amico 等［18］通過將Clapeyron 方程和源于van Genuchten 模型的土-水特征曲線相結(jié)合，建立了水分凍結(jié)曲線與土水勢的關(guān)系，隨后就有學(xué)者證明了此方法的合理性［19］。Sheshukov和Niber［20］結(jié)合簡化后的Clapeyron 方程，拓展了未凍水含量預(yù)測公式，在土水特征曲線中引入了殘余含水量。也有學(xué)者根據(jù)孔隙大小分布函數(shù)，基于孔隙水從大孔隙到小孔隙的凍結(jié)順序，建立了理想成冰條件下的土中水凍結(jié)模型［21-22］。Hitchcock 等［23］基于均勻成核理論、孔隙堵塞及單一孔隙理論闡釋了液體在介孔材料凍融滯后機(jī)理。

多孔介質(zhì)中冰晶預(yù)融現(xiàn)象為未凍水含量的求解提供了新的思路。Cash 等［24］基于預(yù)融理論建立了負(fù)溫下等粒徑球狀顆粒堆積體系的液態(tài)水比例組成，并將未凍水含量簡化為過冷度的函數(shù)。Wettlaufer［25］揭示了冰-液接觸面以及冰晶表面內(nèi)在聯(lián)系，預(yù)融行為明顯受介質(zhì)中雜質(zhì)的影響。Hendrik和Wettlaufer［26］考慮了多孔介質(zhì)中雜質(zhì)濃度以及固體顆粒表面電荷密度的影響，建立了理想等顆粒體系預(yù)融理論。然而，由于土中土顆粒粒徑范圍分布較廣泛，預(yù)融理論在土中并未得到廣泛的應(yīng)用。本文在已有多孔介質(zhì)預(yù)融理論的基礎(chǔ)上，提出針對土體的體積未凍水含量計算方法，并分析了雜質(zhì)濃度以及顆粒粒徑大小對未凍水含量的影響。

1 預(yù)融理論

在一個理想等半徑球狀固體堆積體系中，0 ℃以上純水充滿固體顆粒間隙，假定固體顆粒為均質(zhì)的、各向同性且表面光滑，且在超低溫下孔隙水全部相變?yōu)楸?。發(fā)生預(yù)融（冰晶開始融化）之前，所有間隙中充滿冰晶，當(dāng)溫度上升后，冰晶開始融化，水分的存在形式有薄膜水和間隙水，如圖1 所示。薄膜水主要存在土顆粒表面和冰晶接觸面，而間隙水主要存在于冰晶與土顆粒所包圍的間隙中，其存在溫度由冰-液彎曲界面決定。

圖1 多孔介質(zhì)中液態(tài)水存在形式Fig.1 Existence form of liquid water in porous medium

基于多孔介質(zhì)預(yù)融理論，Cash等［24］首次建立了兩種常見堆積，簡單立方堆積（Simple Cubic Packing，SC）和立方最密堆積（Face-Centered Cubic Packing，F(xiàn)CC），其顆粒排列方式如圖2 所示。堆積體系中液態(tài)水的比例與薄膜厚度、顆粒半徑以及冰-液接觸面半徑之間的關(guān)系，如公式（1）所示。

圖2 等粒徑顆粒排列方式Fig.2 Arrangement of uniform size particle

式中：fp（SC）和fp（FCC）分別表示SC 和FCC 排列下顆粒占總體積的百分比（填充率）；fp（SC）=0.5236［24］；fp（FCC）=0.7405［24］。R表示顆粒半徑；dp表示冰顆粒和固體顆粒之間的水膜厚度；dgb表示冰顆粒之間的水膜厚度；fp是球體顆粒占總體積的比例。r表示冰液接觸面的半徑；其大小會影響凍結(jié)溫度；可通過Gibbs-Thomson關(guān)系確定［27］。

若不考慮系統(tǒng)雜質(zhì)及表面帶電時，固體顆粒表面薄膜厚度d可用公式（2）計算［24,26］。

式中：△T表示過冷度，即初始凍結(jié)溫度降低值（與273.15 K的偏差值）。Lwi表示冰水相變潛熱（J·mol-1）。ρl是孔隙溶液摩爾密度（mol·m-3）。Tm是純水冰點(diǎn)溫度，273.15 K。AH是Hamaker 常數(shù)［27］，不同物質(zhì)接觸形式會影響Hamaker 常數(shù)的取值，故dgb和dp存在一定差異。

對于土而言，系統(tǒng)中含有雜質(zhì)（鹽分），況且很多類型的土顆粒表面帶有電荷，故多孔介質(zhì)體系中應(yīng)該考慮雜質(zhì)和顆粒帶電量的影響。Hendrik 和Wettlaufer［26］推導(dǎo)了含雜質(zhì)、帶電堆積等粒徑顆粒系統(tǒng)中的水膜厚度與過冷度之間的關(guān)系，如公式（3）所示：

式中：c是常數(shù)，可通過公式（4）計算。

式中：nim是土顆粒表面的雜質(zhì)濃度；ε0是真空介電常數(shù)；ε是液態(tài)水的相對介電常數(shù)；NA是阿伏伽德羅常數(shù)；Tm=273.15 K；e是電子電荷；k是玻爾茲曼常數(shù)；qs是固定顆粒表面的電荷密度，孔隙溶液摩爾密度ρl在計算中可認(rèn)為是常數(shù)；R是氣體常數(shù)。公式（3）右邊第一項表示雜質(zhì)濃度影響，第二項表示色散力影響下的界面自由能，第三項表示顆粒表面電荷相關(guān)的界面自由能。

以上理論表明等顆粒多孔介質(zhì)中水分凍結(jié)問題得到了較完美的解決，然而對于不同粒徑組成的土中水分凍結(jié)依然缺乏有效的方法。

2 研究方法

由于土體顆粒分布范圍較為廣泛，例如粉質(zhì)黏土通常從幾個微米到幾個毫米，故現(xiàn)有預(yù)融理論在土中應(yīng)用受到限制。假定土顆粒為球狀顆粒，忽略表面粗糙度，視表面為理想光滑球面。基于土體顆粒級配曲線，取足夠小的顆粒粒徑范圍，通常試驗(yàn)測定的也是某一粒徑范圍Ri土顆粒質(zhì)量（體積）占總顆粒質(zhì)量（體積）比例ωi，顆粒分布及粒徑劃分如圖3所示。根據(jù)參考文獻(xiàn)［28-29］，黃土和砂土的顆粒分布曲線如圖4所示。

圖3 粉質(zhì)黏土粒徑分布Fig.3 Soil particle distribution of silty clay

圖4 黃土和砂土粒徑分布Fig.4 Soil particle distribution of loess and sand

在第i個粒徑范圍內(nèi)，認(rèn)為粒徑均勻分布，土顆粒半徑均為Ri，若不考慮土的結(jié)構(gòu)性差異，理想認(rèn)為土體是以上兩種顆粒堆積，則可利用聯(lián)合公式（1）和公式（3）求解該粒徑下土中未凍水含量。則總的未凍水含量可用公式（5）求解。

式中：ri表示土顆粒半徑為Ri情況下的冰液接觸面的半徑。然而，公式（5）計算較為繁瑣，利用了相同方法重復(fù)求解不同粒徑下的未凍水含量。針對此問題，本文根據(jù)不同粒徑土顆粒的體積占比，提出等效粒徑的方法，具體計算如公式（6）所示。

式中：Re表示土顆粒等效粒徑。

土中土顆粒排列處于SC和FCC之間，為了進(jìn)一步簡化計算，本文采用平均堆積顆粒排列，即取fp（SC）和fp（FCC）的平均值。結(jié)合公式（5）和公式（6），土中液態(tài)含水量可用公式（7）計算。

其中冰-液接觸面彎曲半徑r可通過公式（8）計算［23］。

式中：ρs是冰晶的摩爾密度；γsl是冰液表面自由能。對于冰水相變潛熱，在溫度較小范圍內(nèi)變化時，忽略其變化，取Lwi≈6000 J·mol-1。

隨著溫度升高無限接近冰點(diǎn)時，假定冰晶處于孔隙中心處，無限小可近似認(rèn)為一點(diǎn)，如圖5 所示，此時冰晶具有最大水膜厚度。則可依據(jù)顆粒排列幾何關(guān)系建立不同土顆粒堆積下的孔隙初始體積濃度與水膜表面雜質(zhì)濃度之間的關(guān)聯(lián)，如公式（9）所示。

圖5 最大水膜厚度示意圖Fig.5 Schematic diagram of maximum thickness of water film

式中：N為單位分子雜質(zhì)電解離子數(shù)，對氯化鈉而言，N=2。

為了方便計算土中初始濃度，把孔隙水的初始濃度c0考慮成SC和FCC模型的平均濃度，則雜質(zhì)濃度可用公式（10）計算：

假定初始濃度c0(SC)=c0(FCC)=c0，將公式（9）變形后帶入公式（10），可得到土中平均雜質(zhì)濃度計算公式，見公式（11）。

對于固體顆粒，Hamaker 常數(shù)通常取9.86×10-19J［24］；對于冰晶接觸面而言，表面電荷密度可取0.01 C·m-2以及Hamaker 常數(shù)為-3.3×10-22J［26］；R=8.314 J·(mol·K)-1［27］；ε0=8.85419×10-12F·m-1［27］；ε=80.4 F·m-1［27］；e=1.602×10-19C［27］；k=1.308×10-23J·K-1［27］；NA=6.022×1023mol-1［27］；ρl=5.555×104mol·m-3［24］；ρs=5.094×104mol·m-3［24］；γsl=0.03 J·m-2［24］。

3 結(jié)果與分析

3.1 土中水膜厚度與未凍水含量變化

一般土顆粒表面帶有電荷，以青藏高原粉質(zhì)黏土為研究對象，取0.1 C·m-2［30］，該類土的初始凍結(jié)溫度為-0.2 ℃［31］。已有試驗(yàn)研究表明該冰點(diǎn)溫度歸咎于土中離子濃度［2,32］，基于溶液冰點(diǎn)計算公式［32］可反算出初始離子濃度約0.07 mol·L-1。假定土中鹽分類型為典型氯化鈉鹽，N=2，根據(jù)公式（9a）估算nim為4.5×102μmol·m-2。以圖3粉質(zhì)黏土顆粒級配曲線為例，土顆粒等效半徑為8 μm。依據(jù)公式（3）可建立土顆粒水膜厚度與凍結(jié)溫度降低之間的關(guān)系，計算結(jié)果如圖6所示。圖中斜線表示過冷度隨水膜厚度的變化規(guī)律曲線，線上點(diǎn)K1、K2表示不同水膜厚度時過冷度值。

圖6 土中過冷度與水膜厚度關(guān)系Fig.6 Curve of relationship between degree of supercooling and thickness of water film in soils

從圖6 可以看出，土中水膜厚度與溫度呈現(xiàn)出線性關(guān)系，要使水膜厚度變化小，則要更大的過冷度或者更低的溫度。K1和K2兩個點(diǎn)表明要使水膜厚度變成22 nm、15 nm，則需要過冷度分別達(dá)到20 K 或30 K。若土的初始凍結(jié)溫度為0 ℃，土體溫度從-20 ℃變化至-30 ℃，水膜厚度僅變化7 nm。該規(guī)律表明在-30 ℃以下土中水膜厚度變化很小，一般凍土地區(qū)土中溫度很難使水膜厚度變?yōu)閹讉€納米。故計算土中未凍結(jié)含水量，通?？紤]-30 ℃以內(nèi)變化。

類似地，冰晶間的水膜厚度與過冷度可依據(jù)公式（3）建立相互關(guān)聯(lián)，彎曲冰-液界面的半徑可通過公式（8）求解。將水膜厚度與彎曲液面半徑的值代入公式（7）即可求解不同過冷度下的土中未凍水含量，當(dāng)取等效半徑的0.28 倍計算時，計算值與吻合數(shù)據(jù)較好，計算結(jié)果如圖7 所示。圖中試驗(yàn)數(shù)據(jù)均來自青藏高原粉質(zhì)黏土［15,28,31-33］，可以發(fā)現(xiàn)試驗(yàn)值與計算值吻合較好，尤其在劇烈冰水相變階段（-2～0 ℃）。圖7表明改進(jìn)的公式可以求解負(fù)溫下粉質(zhì)黏土中未凍水含量。

圖7 粉質(zhì)黏土未凍水含量計算值與試驗(yàn)值對比Fig.7 Comparison between the calculated and experimental values of unfrozen water content in silty clay

為了進(jìn)一步驗(yàn)證模型的適用性，選取不同類型黃土、砂土及黏土試驗(yàn)值進(jìn)行對比。不同類型土的顆粒表面電荷量差異較?。?1］，故在計算中均取0.1 C·m-2。通過公式（6）可計算出黃土和砂土的等效粒徑分別為10 μm和56 μm。黏土的等效粒徑大約為4 μm［2,34］。

不考慮土中水過冷度，表面雜質(zhì)濃度與之前粉質(zhì)黏土一樣，同樣對等效半徑進(jìn)行縮小，黏土、黃土、砂土分別以等效半徑的0.24、0.3、0.36 倍計算，未凍水試驗(yàn)數(shù)據(jù)均來源于參考文獻(xiàn)［18,28,34-35］，計算結(jié)果如圖8所示。從圖8a可以看出，計算值能較好反映未凍水變化規(guī)律，當(dāng)溫度較低時，即ΔT＞10 K，預(yù)測值比試驗(yàn)值較高。對于砂土而言，計算值與試驗(yàn)值吻合較好，未凍水體積含量迅速減小并趨于穩(wěn)定，當(dāng)ΔT＞5 K時，未凍水比例基本不變。而對于黏土，雖然計算值反映了變化趨勢，但誤差較大，尤其在開始凍結(jié)溫度段，其主要原因是黏土顆粒的形狀呈片狀，與假定的球形差距較大。

圖8 不同土類型計算值與試驗(yàn)值對比Fig.8 Comparison between the calculated and experimental values of different soils

3.2 水膜厚度影響因素

不同物理化學(xué)界面能對水膜厚度的影響如圖9所示。色散力其與Hamaker 常數(shù)相關(guān)，對水膜的影響與雜質(zhì)濃度沒有關(guān)系，隨著水膜厚度減小，過冷度呈線性增大，當(dāng)水膜減小至幾個納米級時，其逐漸發(fā)揮作用。雜質(zhì)濃度對水膜厚度的變化影響很大，隨著濃度增加，相同過冷度下水膜厚度減小越來越困難，說明水分凍結(jié)難度增大。土顆粒表面電荷密度對水膜厚度的影響極易受雜質(zhì)濃度的影響，隨著濃度增大，土顆粒表面電荷對水膜的影響越來越??；相反，當(dāng)雜質(zhì)濃度減小至0.01 μmol·m-2，土顆粒表面電荷密度對土孔隙中水膜變化的作用大于雜質(zhì)濃度本身的影響。當(dāng)雜質(zhì)濃度大于1 μmol·m-2，表面電荷密度對水膜的影響非常小，在計算中可以不考慮。此外，依據(jù)公式（3），當(dāng)顆粒表面電荷密度減弱時，其對水膜的影響也隨之減小。

圖9 不同雜質(zhì)濃度及物理化學(xué)界面能對水膜厚度影響Fig.9 Effect of different impurity concentration and physical chemical interface energy

3.3 薄膜水和間隙水變化規(guī)律及影響因素

對于粉質(zhì)黏土，顆粒表面水膜，冰晶間水膜以及間隙水所引起未凍水含量的變化如圖10 所示。可以看出，土中水分凍結(jié)主要由土顆粒表面水膜厚度變化來決定。冰晶間水膜的影響相比土顆粒表面水膜而言較小，兩者計算值之間差了一個數(shù)量級。而土中間隙水對未凍水含量的影響最小，通過公式（8）得出，較高的初始濃度會減小冰-液彎曲界面的半徑，導(dǎo)致r/Re的平方減小很多，致使液態(tài)含水量很小。

圖10 水膜及間隙水對未凍水含量計算值的貢獻(xiàn)Fig.10 Construction of water film and interstitial water to calculation of unfrozen water content

當(dāng)土中雜質(zhì)濃度變化或土顆粒粒徑變化，諸如黏土比粉質(zhì)黏土的顆粒要細(xì)，在此假定不同濃度及粒徑條件，未凍水含量計算結(jié)果如圖11所示?？梢园l(fā)現(xiàn)在濃度較高的條件下（圖11a），隨著土的等效顆粒粒徑減小，未凍水含量增大，表明水在細(xì)顆粒土中更難凍結(jié)。在濃度較低情況下（圖11b），未凍水含量變化表現(xiàn)出相似的規(guī)律，但基于間隙水計算得到的未凍水含量所發(fā)揮的作用隨著顆粒減小愈發(fā)明顯，在很大溫度范圍內(nèi)，超過了基于冰晶間水膜厚度計算的未凍結(jié)含水量值。此外，計算結(jié)果表明土顆粒等效粒徑越細(xì)，初始凍結(jié)溫度較低，低濃度細(xì)顆粒的冰-液彎曲界面在水分初始凍結(jié)中占主導(dǎo)地位。

圖11 濃度及顆粒半徑對未凍水含量的影響Fig.11 Effect of concentration and particle radius on unfrozen water content

4 討論

經(jīng)典的指數(shù)或冪函數(shù)關(guān)系［8］在水分劇烈相變階段（初始凍結(jié)溫度附近）存在很大的誤差，因?yàn)橹笖?shù)函數(shù)在冰點(diǎn)0 ℃不收斂，且模型中參數(shù)隨不同初始含水量變化較大。結(jié)合土水特征曲線和凍結(jié)特征曲線的預(yù)測模型［17］很大程度上依賴于初始含水量和殘余含水量，且模型中涉及的參數(shù)較多。這些未凍水含量模型中的參數(shù)沒有明確的物理意義，而本文改進(jìn)的模型物理參數(shù)明確，且考慮了土體顆粒大小、雜質(zhì)濃度、顆粒表面帶電特性，并涉及了水分凍結(jié)的物理化學(xué)微觀機(jī)理。

Dash 等［27］通過預(yù)融理論計算了等粒徑聚苯乙烯粉末及石墨化炭黑粉中未凍水含量的分布，發(fā)現(xiàn)當(dāng)半徑分別縮小0.15和0.4時計算效果與試驗(yàn)值吻合較好。土顆粒實(shí)際排列非常復(fù)雜，且顆粒形狀并非全是球狀體，且土顆粒表面呈現(xiàn)出一定粗糙度。本文在計算過程中采用了等效粒徑，當(dāng)粉質(zhì)黏土、黃土、砂土分別以等效半徑的0.28、0.3、0.36倍計算時，計算值與試驗(yàn)值吻合較好。而對于細(xì)顆粒土黏土，計算值與試驗(yàn)值差別較大。此外，在模型中未考慮細(xì)顆粒土中特殊結(jié)構(gòu)，大部分黏土中存在如蜂窩結(jié)構(gòu)、絮狀結(jié)構(gòu)等特殊形式，與自然排列下的土顆粒分布有一定差距，會影響模型計算精度。同時，計算中縮小半徑的處理方式還未有詳細(xì)闡述。所以顆粒形狀、顆粒表面粗糙度、土顆粒的實(shí)際排列和土的結(jié)構(gòu)類型等會作為水分微觀凍結(jié)方面后期研究內(nèi)容，進(jìn)而開展進(jìn)一步研究。

5 結(jié)論

本文基于多孔介質(zhì)預(yù)融理論，開展了土中未凍水體積含量隨溫度變化研究，提出了其計算方法，并分析討論了雜質(zhì)濃度和粒徑大小對水膜厚度及未凍水含量的影響，具體結(jié)論如下：

（1）對于較大顆粒土，未凍水含量模型預(yù)測效果較好。當(dāng)粉質(zhì)黏土、黃土、砂土等效粒徑分別縮小0.28、0.3、0.36倍時，未凍水含量計算結(jié)果與試驗(yàn)值接近，模型可較好地預(yù)測土中不同過冷度下液態(tài)水含量。而對于黏土，未凍水含量預(yù)測誤差較大。

（2）土顆粒表面電荷密度對水膜厚度的影響極易受雜質(zhì)濃度的影響，隨著濃度增大，土顆粒表面電荷對水膜的影響越來越小。隨著雜質(zhì)濃度增加，相同過冷度下水膜厚度增加，水分凍結(jié)難度增大。

（3）土中未凍水含量計算可以不考慮色散力的影響，水分凍結(jié)主要由土顆粒表面水膜厚度變化來決定。當(dāng)?shù)刃ьw粒減小時，通過間隙水求解的未凍水體積比逐漸增大，并且在較低濃度時發(fā)揮更大作用。

（4）當(dāng)土顆粒表面雜質(zhì)濃度小于1 μmol·m-2，土中水分凍結(jié)溫度由冰-液接觸面半徑?jīng)Q定，該半徑隨等效粒徑減小而減小，且半徑越小，凍結(jié)溫度越低。