冻土多物理场耦合作用理论与模型研究现状

王朋宾，夏胡熙，姜海强，马勤国

（华南理工大学土木与交通学院，华南理工大学亚热带建筑科学国家重点实验室，华南理工大学华南岩土工程研究院，广州 510640）

冻土和盐渍土在我国分布广泛，季节冻土覆盖北方大部分地区[1]。盐渍土与季节冻土分布具有较好一致性，主要分布在新疆、甘肃、内蒙古、黑龙江、陕西等省区。土体冻胀变形导致建筑物产生不均匀沉降、道路路面不平，温度变化引起盐分结晶可导致盐胀、建筑物被腐蚀等问题。冻土中水分、温度、盐分和应力相互影响，引起作物冻胀和盐胀等病害严重制约寒旱区经济发展。冻土性质研究是解决此问题基础，研究冻胀盐胀机制和水热盐力耦合作用机理对于解决冻土病害问题具有重要意义。

1 冻土多场耦合理论基础

1.1 水分迁移驱动力

土体冻结过程中，水分从暖端向冷端迁移，引起水分重分布。同时，溶解在水中的盐随水分发生迁移和结晶，导致土体冻胀和盐胀[2]。冻土水分迁移驱动力是冻土多场耦合研究基础问题之一。单向冻结试验是研究冻结过程土体水盐迁移特性重要方法。何菲等开展室内单向冻结试验显示，水分迁移现象随初始含水量增加和顶板负温降低变得明显，证明初始含水率和温度梯度对水分迁移有显著影响[3]。翟金榜等开展室内大尺寸单向冻结试验，结果显示，当有外界水源补给时，土体含水量增大且冻结区含水量增加较多[4]。王文华采用直接法和冻融法开展毛细水上升试验，研究碳酸盐渍土水分迁移特性，结果显示盐分含量高土样毛细现象较弱，毛细上升不明显[5]。张彧等选取青海省察尔汗至格尔木高速公路一段作现场试验，研究高氯盐含量盐渍土水盐迁移特性[6]。结果显示，热效应是影响水盐迁移主要因素，且土体热效应作用随深度增加而减弱，最终形成盐分聚集区。肖泽岸等利用单向冻结试验研究氯化钠盐渍土水盐迁移特性发现，盐分在冻结锋面处聚集抑制水分迁移和冻胀[7]。

土壤水分发生迁移是因水分和热量不平衡，包括物理、化学和力学作用。为解释土壤水分迁移现象，研究人员提出各种水分迁移驱动力理论。Everett提出毛细理论，认为冰水界面存在压力差，这种压力差引起的毛细吸力是水分迁移驱动力[8]。Taylor等模拟未冻水迁移过程，提出未冻水含量梯度是水分迁移驱动力[9]。在补给盐溶液条件下，徐敩祖对不含盐冻土进行冻结试验，发现盐溶液在水势梯度作用下发生迁移[10]。季雨坤建立水-热-力耦合冻胀模型时，将基于压力势梯度的冻吸力作为水分迁移驱动力，体现荷载与冻胀间相互作用[11]。系统中能量差被认为是物质传递驱动力，吴道勇从吉布斯自由能角度研究冻土中物质传递驱动力[12]。土水势从势能角度解释土中水分迁移驱动力，是一种在冻土和盐渍土中使用较广泛的水分迁移驱动力理论。在不同初始和边界条件下，土水势可由温度势、重力势、溶质势、基质势和压力势中某几项组成。使用土水势理论时，往往根据具体情况简化。原国红考虑非饱和土孔隙连通性及土体内温度势引起的水分迁移量较小，忽略压力势和温度势影响[13]。但这种简化也存在一些问题，例如土体内并非每一点压力均等于大气压力，压力势对土中水分迁移也有影响。土体组分、含水率和密度等性质也影响土水势，但现有研究常常忽略这一点。温度影响水黏滞力和表面张力，Lai和Zhang等均考虑温度势对水分迁移的影响[14-15]。土中水分迁移受势能、应力和温度等多种因素影响。基于不同假设和试验条件，国内外学者研究水分迁移驱动力。但由于冻土本身结构和性质复杂，现在尚无一种理论可完美解释冻土中水分迁移驱动力。

1.2 冻胀、盐胀机制

1.2.1 冻胀机制

冻胀由两部分组成。一部分是孔隙水原位冻结引起的体积膨胀，这部分主要是水分体相变增大量（约9%），引起的冻胀量较小。另一部分是迁移水产生的冻胀，未冻区水分向冻结缘迁移，在冻结缘与已冻土间发生相变形成冰透镜体，这部分体积增大量约为迁移水体积1.09倍。冰透镜体是典型土体冻结伴生现象，对土体冻胀起决定性作用。Taber认为吸力是冻胀驱动力，吸力大小决定土颗粒表面水膜平衡厚度，但吸力理论未在实际中得到应用[16]。Konrad等提出分凝势概念，认为水分向冰透镜体补给导致冻胀[17]。分凝势模型描述单个冰透镜体生长过程，可用于工程中冻胀量计算，仅适用于温度梯度已知状况。毛细理论将土体连通孔隙看作毛细管，认为冰透镜体生长与多孔介质物理性质相关，冰水两相间压力差是水分迁移和冰透镜体生成的原因。基于毛细理论，Peppin将热力学平衡状态下土体划分为三层，即冰层（土体完全冻结）、饱和土层和储水层（未冻土），并认为平衡状态下未冻水未发生迁移[18]。维持平衡状态所需水压力可表示为

其中，P0为上覆层压力（Pa）；ρw为水密度（kg·m-3）；Lf为冰水相变潜热（J·mol-1）；Tm为冻结温度（K）。

当温度降低，上覆土层压力不变时，维持平衡状态所需压力减小，而储水层水压力不变，储水层中水分在相对较高压力作用下向冻结区迁移、冻结。当储水层获得持续补水时，理论上可使冰透镜体不断生长。同时，增加上覆荷载也可抑制冰透镜体生长。

但是，毛细理论无法解释不连续冰透镜体的形成。针对这一问题，Miller等提出冻结缘理论[19]。冻结缘内同时存在冰和水，且无冻胀发生，冰的存在抑制水分向冻结区迁移。Zhou等建立一种预测冻胀和冰透镜体生长速率的理论模型，指出孔隙水渗入冰透镜体底部未冻水膜以补充冰透镜体生长；这一过程导致冰透镜体与未冻土或冻结缘界面处产生水力阻力，毛细理论和冻结缘理论均忽略各向异性的冰应力，无法较好解释这一阻力影响，新模型可解决这一问题[20]。冻结缘是未冻土向冻土发展的过渡区，但Watanabe在试验中利用可探测孔隙冰的拉曼光谱技术观测冰透镜体与未冻土间区域，该区域并未探测到孔隙冰[21],由此判断冻结缘不存在。例如，当土体主要发生原位冻结时不会形成冻结缘，因冻结锋面迁移速度过快，水分来不及迁移，但这种解释无法适用于所有试验。因此，关于冻结缘是否存在的争论还在继续。

目前，多孔介质中冰透镜体生长机制并不明确。Arenson利用高分辨率摄像头观察研究粉土冻结过程中冰分凝现象，提出一种冰透镜体生成假设：土体冻结产生垂直裂隙，水分通过这些裂隙迁移形成冰透镜体[22]。这种假说认为，冻结使多孔介质产生裂隙，然后冰透镜体在裂隙中生长。但Vlahou等认为，岩石孔隙和微裂纹中水分冻结，使孔隙和裂纹扩展，导致岩石开裂和生成新冰透镜体[23]。Style等试验发现，冰的生长导致土壤产生裂隙，裂隙壁承受逐渐增大冰压力，压力促使裂隙扩展，形成新冰透镜体[24]。因此，冰透镜体生长机制仍需进一步研究。

1.2.2 盐胀机制

结晶压力理论是研究盐胀问题基础。Correns最早指出结晶压力是盐分导致土体破坏的原因，并给出结晶压力表达式[25]。Lai等认为，宏观结晶应力决定土体变形，并建立宏观结晶应力表达式[26]。针对结晶压力计算模型缺少试验验证的问题，谭仁义等利用土压力传感器直接测量盐分结晶过程中土压力[27]。结果显示，实测值与宏观结晶压力计算值吻合较好，证明Lai等[26]结晶压力计算模型的准确性。

研究人员针对盐渍土盐胀特性开展广泛研究。万旭升等室内试验显示，Na2SO4盐渍土盐胀主要发生在正温区间，负温区间盐晶体析出较少，土体主要发生冻胀[28]。吴道勇等利用不同含盐量的土壤作模型试验，结果显示水盐共同作用导致盐渍化冻土变形[29]。且对于含盐量较低土体，冻胀和融沉是变形主导因素。对于含盐量大于1%土体，盐胀和溶陷是变形主导因素。但并非所有盐渍土均发生盐胀，肖泽岸等指出，由于NaCl冰盐共晶点较低（-22℃），当温度高于共晶点时，土体冻结过程中仅产生冰晶[7]。而溶解度受温度影响较大的盐在冻结过程中更易使土体产生盐胀破坏，如硫酸盐和碳酸盐。土体盐胀受多种因素影响，充分考虑各种因素作用有助于理解盐渍土发生冻胀盐胀的机理，也有助于改进工程中防盐胀措施。

1.3 冻融条件下多场耦合作用

土体冻胀、盐胀变形受水分、温度、盐分和应力综合作用（见图1）。冻土中水分迁移是指未冻水迁移，温度是决定冻土中未冻水含量主要因素。温度降低，孔隙水冻结成冰，冰晶堵塞孔隙，导致土体渗透系数急剧减小。Horiguchi和Miller研究0～-0.35℃冻土渗透系数随温度变化，结果发现，土壤渗透系数数量级从冻结前的10-8m·s-1下降到10-12～10-13m·s-1[30]。张虎等测量含水率为50%青藏高原粉质黏土在高负温区间的渗透系数，结果显示，在-0.6～+0.1℃范围内，粉质黏土渗透系数处于8.22×10-11～7.19×10-9cm·s-1[31]，温度升高，冰晶融化，但冰水相变引起的体积膨胀对土体孔隙结构造成破坏，因此冰晶融化后土渗透系数增大[32]。由此可见，温度影响冰水相变，间接导致渗透系数改变，最终影响土中水分迁移。

图1 冻土中水-热-盐-力耦合机制Fig.1 Sketch for hydro-thermal-salt-stress coupling mechanisms in freezing/frozen soils

梁月通过数值计算方法对比不同含水条件下埋地管道周围土壤温度场，发现无水冻土整体温度明显低于饱和含水冻土，水分迁移和冰水相变对冻土温度场分布产生较大影响[33]。张泽等采用数值模拟方法对是否考虑水分迁移的隧道温度场进行计算和对比，结果显示，当仅考虑相变而不考虑水分迁移时，隧道温度场热传导速度更快，证明相变潜热对隧道温度场分布影响远大于热对流[34]。考虑到水分在土体中流动速度较慢，且与相变释放的潜热相比对流传热影响较小，一些学者在研究中忽略对流传热影响。但张明礼等研究发现，发生在含水量较大土壤的水分迁移对土体传热影响较大，土体含水量越大，对流传热作用越显著，这种情况下，对流传热影响不可忽略[35]。粗颗粒土孔隙较大，冻结过程中孔隙不易被完全堵塞。因此，相比于细粒土，粗粒土渗流问题更复杂，以渗流为基础的热对流对粗粒土温度场的影响也不能忽略。

盐分进入土体孔隙才对多孔介质造成破坏，而这一过程必须有水分参与。Bing等通过试验发现，土壤水分迁移机理在很大程度上决定盐分迁移机理[36]。证明盐分以水为载体发生迁移。未冻水含量增加可溶解更多盐，未冻水含量减少则有利于盐晶体析出。土中盐分也影响土体渗透性。盐分结晶相变对土体渗透性的影响与冰分凝类似，Na2SO4吸水结晶生成Na2SO4·10H2O不仅引起体积膨胀，还减少未冻水含量。生成的晶体堵塞土体孔隙，阻碍水分流动。晶体溶解时又释放水分子，导致未冻水含量增大。

土中孔隙水在过冷度驱动下冻结成冰，盐分在过饱和度驱动下结晶[12]。水盐结晶过程影响土压力。谭仁义等试验结果显示，随水分冻结和盐分结晶，土体压力先迅速增大，再逐渐减小至稳定；随土壤含盐量增加，盐结晶压力增大，冰晶压力减小[27]。土体冻结区结晶作用产生正孔隙压力，在未冻区产生负孔隙压力，正孔隙压力抑制水分迁移，破坏土体结构，使土体孔隙增大，为晶体生成提供空间；负孔隙压力则促进水分向冻结区迁移。吴道勇研究发现，土中含盐量增加会降低负孔隙压力梯度，从而抑制土中溶液迁移[12]。

温度变化引起孔隙中固-液相相互转化，包括孔隙水冻融和水合物合成与分解。盐溶解度也受温度影响。在低于32.4℃范围内，硫酸钠溶解度随温度降低而下降，溶解在水中的硫酸钠和10个水分子结合生成芒硝（Na2SO4·10H2O）。当温度升高时，析出晶体重新溶解。盐分结晶相变过程中释放的潜热对土体温度场产生影响。盐分也影响土热学性质，其他条件不变情况下，含盐量增加导致土体冻结温度降低。温度变化也对应力场产生影响。高娟等试验结果显示，温度越低，冻结砂土强度越大，原因是孔隙水冻结增加颗粒间黏结力，孔隙水压力变成有效应力[37]。压力对温度场的影响主要包括两个方面：一是压力影响土体冻结温度，土体冻结温度随压力增大而降低；二是压力影响土体温度，Wu等通过单向冻结试验发现，当施加较大外部荷载时，土柱同一位置处温度降低[38]。现有研究一般认为，外部荷载可通过孔隙溶液迁移影响热量传递，从而影响温度场。此外，在外部荷载作用下，冻土内部出现压融现象。黄永庭等通过室内试验发现，在-8～24℃正弦波动周期温度边界条件下，荷载对冻土融化速率有显著影响；无荷载作用时，冻融循环4次后，随冻融循环次数增加，融化锋面位置不变；在上覆荷载作用下，试样在第2次冻融过程中全部融化[39]。压力对土体温度场的影响及冻土压融问题仍需进一步研究。

2 未冻水含量经验公式

2.1 未冻水含量与温度关系

即使在冻结点以下，土中仍存在未冻结的液态水。土中未冻水含量显著影响冻土热力学性质和工程性质，是寒区工程稳定性重要参考指标。研究人员常用各种仪器直接测量未冻水含量，但时间和经济成本较大。为解决此问题，研究人员利用实测数据建立预测未冻水含量经验模型。未冻水含量预测模型可减少未知变量，是数值模拟研究中常用联系方程之一。

Anderson等认为未冻水含量可直接表示为温度的函数[40]：

其中，Wu、W0分别为质量未冻水含量和初始质量未冻水含量（%）；α、β为拟合参数；T、Tf分别为任意时刻土壤温度和土的初始冻结温度。

徐敩祖等认为，未冻水含量与温度保持动态平衡，并结合试验数据给出未冻水含量与温度关系[1]：

其中，Wu为质量未冻水含量（%）；Wp为塑限含水量（%）；Wl为液限含水量（%）；a、b为与土的性质有关常数；T、Tl、Tp分别为土体温度、液限含水量时冻结温度和塑限含水量时冻结温度。

Anderson等[40]和徐敩祖等[1]提出的未冻水含量模型形式简单，在建立多物理场耦合模型时被广泛使用，但质量未冻水含量在试验中难以直接测量。Zhang等假设土体温度降至-273.15℃以下时土中液态水全部冻结[41]，在此基础上，将体积未冻水含量与温度的关系（式6）应用到多年冻土地区阴阳坡路堤稳定性分析中，验证预测模型准确性。

其中，θu、θ0分别为体积未冻水含量和初始体积未冻水含量；T、Tf、Tk0分别为土体温度、孔隙水初始冻结温度和将摄氏温度转为开尔文温度的常数，取273.15；ω为与土体性质相关的拟合参数。

Lu等结合现场观测数据对比评价Anderson[40]和Zhang[41]的未冻水含量预测模型，指出在一些土样中Zhang等[41]模型得到的未冻水含量预测值比实际值偏高，但整体上预测模型仍可靠且更加方便[42]。

2.2 未冻水含量其他表示方法

未冻水含量与温度的经验关系式为未冻水研究提供很大便利，但这些关系式也存在不足。例如，无法体现降温速率对未冻水含量变化的影响，也不能体现土壤物理性质对未冻水含量的影响。且未冻水含量与温度的经验关系适用于冻结稳定状态，而试验和实际情况下土体冻结往往不稳定。因此，除用温度表示未冻水含量，学者们也提出其他预测未冻水含量的方法。

Wang等将冻土中未冻水表示为未冻孔隙和已冻孔隙中未冻水之和，建立与土壤孔径分布相关未冻水含量三参数模型[43]：

其中，θu、θuu、θuf、θs分别为体积未冻水含量、未冻孔隙中体积未冻水含量、冻结孔隙中体积未冻水含量和饱和体积水含量；T、T0分别为温度和多孔介质中水分的凝固点；为多孔介质最小孔隙半径中水分凝固点降低值；α为多孔介质球形空隙表面积；φm为基质势；a、n、m为参数。

Mu等通过区分毛细水和吸附水冻结过程，提出一种考虑初始孔隙率影响的未冻水含量预测模型[44]：

其中，θ、θs、θamax分别为体积未冻水含量、饱和体积水含量和最大体积吸附水含量；e为初始孔隙率；m0、m1、m2、m3、k为参数；T、Tmin分别为温度和所有孔隙水被冻结时温度；rw、Lw分别为水的密度和冰的融化潜热。

Mu等提出的模型有助于提高在较宽温度范围内（-50℃至0℃）计算体积未冻水含量精度，且考虑初始孔隙率的影响[44]。模拟结果显示，随初始孔隙率降低，土体保水能力增强。Wang和Mu提出的模型精度较高，且考虑土壤物理性质对未冻水含量的影响，但这两种模型较复杂，需确定参数较多[43-44]，在一定程度上限制模型实际应用。

3 多场耦合模型研究

为进一步探究各物理场间作用机制，研究人员先后提出水-热、水-热-力等多场耦合模型，多场耦合理论成为研究热点。水-热耦合模型发展至今仍被广泛应用于模拟计算和科学研究。Harlan考虑相变潜热和水分对流传热，建立水-热耦合模型[45]。该模型本质上是对未冻土和部分冻结土壤水分流动的模拟，研究重点是水分迁移和重分布。在Harlan模型基础上，胡和平等建立一维水-热耦合模型，对冻土自然冻结过程进行模拟计算，并通过试验验证模型准确性[46]。周祖昊等将水-热耦合方程与“积雪-土壤-砂砾石层”连续体模型相结合，构建适用于青藏高原地区地质和气候条件的计算模型[47]。通过与实测值对比，验证模型准确性，为研究青藏高原地区水文循环过程提供新方法。

水-热耦合模型反映水分和温度间耦合关系。但与土体多场耦合作用相比，两场耦合仍存在巨大局限性。盐渍化冻土区存在复杂的水-热-盐和水-热-力耦合作用，仅考虑水热耦合无法满足工程实际需要。因此，国内外学者开始在水-热耦合模型基础上研究三场耦合模型。为实现水、热、力三场耦合，武建军在建立水分迁移方程时考虑冻土应变影响[48]。许强等考虑冻结过程中土体体积应变与温度间耦合作用，实现三场直接耦合[49]。Shen等假设冻结缘内冰压与温度按线性分布，建立冻土水-热-力耦合模型[50]。该模型考虑冻胀对应力的影响，可预测冻结锋面后水和冰的积聚量，但无法预测离散冰透镜体位置。杨韬等将水-热-力耦合模型与应力-温度耦合损伤本构模型相结合，预测高温冻土融化固结过程，模型考虑大应变状态应力方程，提高预测精度[51]。以上模型均忽略土中气相影响，主要考虑饱和土冻结。实际上多数土体处于非饱和状态，水分以液态水、水汽和冰三种相态存在，水汽在土中水分迁移和相变过程中起到重要作用。Yin等考虑土中气相影响，建立水、水汽、热和应力耦合模型，并将孔隙率作为描述孔隙压力分布重要变量[52]。

盐渍土和冻土在我国分布广泛，且相当一部分地区存在盐渍土病害与冻土病害重叠情况，而关于盐渍土冻结过程中水盐运移、传热、盐分结晶和积聚、土体变形之间相互作用的研究较少，为盐渍土提供完整四场耦合模型十分重要[15]。Zhang等在水热力耦合模型基础上引入溶质守恒方程，建立水-热-盐-力耦合模型，但模型没有考虑盐分对土壤变形和水分场的影响[53]。冯瑞玲等在盐分运移方程中分别考虑物理和化学作用影响，将等效含水量、等效含盐量作为耦合因子，建立水-热-盐-力四场耦合方程[54]。以上研究主要考虑温度对水盐迁移驱动作用，并通过经验公式建立未冻水含量与温度关系，未真正体现各物理场耦合。针对这一问题，吴道勇在研究中改进未冻水含量与温度的经验关系式，采用试验和数值计算方法对冻结盐渍土水-热-盐-力耦合问题进行深入研究，提出饱和盐渍土水-热-盐-力四场耦合模型[12]。该模型较好体现水分迁移、传热、结晶和变形之间耦合关系。Zhang等基于结晶动力学理论提出饱和硫酸盐渍土水-热-盐-力耦合模型，将水活度和溶液过饱和度作为水盐结晶相变驱动力，建立水盐晶体生长模型，并通过理论推导比较严谨的解释各物理场之间联系[15]。

4 讨论与展望

多场耦合是一个复杂过程。耦合方程中涉及的参数随土体结构、温度、含水量改变而变化。但是，许多研究中参数来源于经验表达式或室内试验，与现场情况存在差距。例如未冻水含量与温度的经验关系未体现动态变化特点，也未考虑盐分影响。而实际情况中，盐分可能与未冻水结合以晶体形式析出，导致未冻水含量减少。未来应在加强室内试验研究同时，从参数角度出发，探讨描述参数变化的准确方法，从而使冻结盐渍土多物理场耦合关系更紧密。

在考虑模型边界条件时，为简化模型，通常设置为固定温度边界、含水量边界或给定一个函数。但实际情况中，土体往往处于复杂多变状态。比如，冬季土中盐分随水迁移至地表附近，融化时由于地表强烈蒸发作用，盐分在地表结晶析出，地表附近产生盐渍化；太阳辐射强弱和环境温度变化会影响水分蒸发；降雨影响地表土体热参数，雨水下渗影响土体含水量，积雪则直接影响土体温度和冻结深度。因此，耦合模型中如何设置更加接近实际情况的边界条件，如何将试验地区气候条件纳入考虑范围是研究重点之一。

确定水分冻结温度是耦合计算的一个关键问题。许多研究设置冻结温度时考虑过冷度影响，将水分冻结温度设置为0℃以下，但过冷阶段通常被忽略，当温度降至冰点时认为水分被冻结。实际上随冰晶生长，水分相变产生的冻胀压力使土体内未冻水压力增加，水分冻结温度下降，即冻结过程中，水分冰点不断变化。开展数值模拟计算时，可进一步研究过冷阶段和降温过程中未冻水冻结温度变化带来的影响。

尽管国内外学者对多物理场耦合问题进行广泛探索，但现有模型仍未完全模拟盐胀冻胀过程，数值模拟精度有待提升。研究人员已针对土体水-热和水-热-力耦合问题进行大量研究，建立许多耦合模型，但针对水-热-盐、水-热-盐-力耦合问题的研究仍较少。原因可能是盐分在土体冻融过程中涉及更加复杂的物理化学变化，而现有本构模型无法准确描述此变化。此外，多数模型仅考虑一维情况、短期冻融或小尺寸模型，与实际长期冻结、大尺寸冻土体情况存在差距。因此，未来应继续探索冻土本构模型，考虑复杂变温条件和力学条件，建立能准确描述冻土受力变形行为的本构模型。同时，要建立多因素影响下大尺寸、多维度和长期冻结试验模型。结合现场环境，开展1∶1模型试验，得到更接近真实情况的试验结果和模拟结果。