土中水膜厚度變化規律及未凍水含量預測方法

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

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

我國凍土面積分布廣泛，多年凍土和季節性凍土面積之和超過國土面積的70%［1］。寒區工程建設面臨因水熱相互作用引發的凍脹融沉現象，嚴重制約著公路、鐵路以及構筑物等基礎設施建設。水分是土體發生凍融的物質基礎，未凍水含量的改變直接影響凍土的熱力學及變形特性。同時，未凍水含量隨溫度變化規律也是水熱數值計算的必要條件。

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

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

1 預融理論

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

圖1 多孔介質中液態水存在形式Fig.1 Existence form of liquid water in porous medium

基于多孔介質預融理論，Cash等［24］首次建立了兩種常見堆積，簡單立方堆積（Simple Cubic Packing，SC）和立方最密堆積（Face-Centered Cubic Packing，FCC），其顆粒排列方式如圖2 所示。堆積體系中液態水的比例與薄膜厚度、顆粒半徑以及冰-液接觸面半徑之間的關系，如公式（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表示冰液接觸面的半徑；其大小會影響凍結溫度；可通過Gibbs-Thomson關系確定［27］。

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

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

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

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

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

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

2 研究方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對于固體顆粒，Hamaker 常數通常取9.86×10-19J［24］；對于冰晶接觸面而言，表面電荷密度可取0.01 C·m-2以及Hamaker 常數為-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 結果與分析

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

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

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

從圖6 可以看出，土中水膜厚度與溫度呈現出線性關系，要使水膜厚度變化小，則要更大的過冷度或者更低的溫度。K1和K2兩個點表明要使水膜厚度變成22 nm、15 nm，則需要過冷度分別達到20 K 或30 K。若土的初始凍結溫度為0 ℃，土體溫度從-20 ℃變化至-30 ℃，水膜厚度僅變化7 nm。該規律表明在-30 ℃以下土中水膜厚度變化很小，一般凍土地區土中溫度很難使水膜厚度變為幾個納米。故計算土中未凍結含水量，通?？紤]-30 ℃以內變化。

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

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

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

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

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

3.2 水膜厚度影響因素

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

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

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

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

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

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

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

4 討論

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

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

5 結論

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

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

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

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

（4）當土顆粒表面雜質濃度小于1 μmol·m-2，土中水分凍結溫度由冰-液接觸面半徑決定，該半徑隨等效粒徑減小而減小，且半徑越小，凍結溫度越低。