川藏铁路季节性粗颗粒冻土边坡水热过程模拟

梁树，谢强，赵文，叶唐进, ，蒋良文

川藏铁路季节性粗颗粒冻土边坡水热过程模拟

梁树1，谢强1，赵文1，叶唐进1,2，蒋良文3

(1. 西南交通大学 地球科学与环境工程学院，四川 成都 610031；2. 西藏大学 工学院，西藏 拉萨 850000；3. 中铁二院工程集团有限责任公司，四川 成都 610032)

以川藏铁路东段的季节性粗颗粒冻土边坡为研究对象，通过建立一个考虑积雪、渗流和年循环气温作用下的带相变的瞬态水热耦合的饱和-非饱和的多孔多相介质数值计算模型，结合野外实测，分析季节性粗颗粒冻土边坡的温度场、水分场的分布特征、冻结深度及其影响因素。研究结果表明：积雪消融入渗改变季节性粗颗粒冻土的水分场，地表可形成最大0.8 m的暂态饱和区。水分场的动态变化提高了热传递速度，增强了冻结能力，边坡冻结深度增大60%，冻结速率增大30%，融化速率增大200%。地下水热对流作用抑制土体冻结，加速土体融化，其中坡脚地下水出露边坡的冻融深度为地下水深埋边坡的63%，冻结速率为79%，融化速率增大1倍。川藏铁路新都桥地区季节性粗颗粒冻土边坡的冻结深度为1.0 m，最大可达到1.9 m。边坡不同位置的冻结深度不同，在坡肩处最高，坡脚处最低；进入融化期后，因冻结深度小及地下水热流作用，坡脚处最先融化。

川藏铁路；季节性冻土；季节性粗颗粒冻土；冻融深度；水热耦合；数值模拟

粗颗粒土是典型的多孔介质，其孔隙大，孔隙连通率高，渗透系数大，有利于流体的流动。流体流动过程中产生的热对流，成为粗颗粒冻土热传递的重要方式。季节性粗颗粒冻土主要研究冻融循环作用下的表层非饱和−饱和粗颗粒土体，其热传递过程受到积雪、积雪融化、地下水渗流、外界温度等多方面的影响，是一个多孔多相介质带相变的渗流场、温度场复杂问题[1−3]。对于渗流场，非饱和粗颗粒土冻结后，孔隙中的水相变，渗流主要受到2个方面的影响。一是液态水变为固态冰以后，冰填充孔隙，降低土体的孔隙率，进而降低土体的渗透性。细颗粒土中的孔隙小，水变成冰体积增大，有可能充满孔隙或隔断过水通道，因此细颗粒土的渗透性受到的影响很大，往往可将其视为不透水层进行简化处理[4]。而粗颗粒土的孔隙大，水变成冰的过程中难以充满孔隙和隔断过水通道。因此，季节性粗颗粒冻土的渗透性相对复杂，受到含水率、孔隙率、温度、水的流速等多方面的综合影响，与热传递过程协调变化。二是负温对水的基质吸力产生较大的影响，而基质吸力是非饱和渗流的重要驱动力，因此，非饱和渗流受负温的影响[5]。对于温度场，由于粗颗粒冻土渗流的复杂性，以渗流为基础的热对流作为热传递的主要方式，热场(温度场)受到的影响也很大。此外，温度变化，土体的导热系数、比热容均随之变化，加上相变潜热，使得温度场受到的影响因素也较多。川藏铁路康定−林芝段途径川西新都桥、理塘等高寒地区，这些地区的坡麓、沟槽等平缓地带是线路通过的首选。初步调查表明，线路首选的这些平缓地带多为坡洪积层，以粗颗粒土为主，浅层土体季节性冻融，为粗颗粒冻土。虽然冻土地区一般遵循“宁填不挖”的基本原则，但在实际工程中，如车站、隧道进出库等，仍然不可避免开挖，产生季节性粗颗粒冻土边坡。这些地区的这类边坡地下水埋深浅、流速快，地表积雪易消融入渗，不同于其他细颗粒土或多年冻土问题，其针对性的水热耦合研究资料还不够充分。本文对川藏铁路新都桥地区某G318国道附近的季节性粗颗粒冻土边坡，进行为期1 a的野外实测。结合多孔多相介质的带相变的渗流场理论、温度场理论及其耦合理论，建立了一个考虑积雪、渗流和年循环气温作用下的含相变的瞬态水热耦合的饱和−非饱和的季节性粗颗粒冻土边坡数值模型，分析边坡的温度场、水分场的分布特征、冻结深度及其影响因素。

1 渗流及热传导理论

1.1 非饱和渗流理论

根据渗流连续性方程，有：

式中：k(k)为()方向的渗透系数；为源汇项；为体积含水率；为时间；为总水头；=μ/γ+, μ为基质吸力；γ为水的重度；为位置水头。

当考虑水的相变时，基质吸力随负温的变化而变化。Black等[6]建立了基质吸力与负温的线性 关系：

式中：Δμ为基质吸力的变化量；为土的温度(低于0 ℃)；系数1 110 kPa/℃包含了冻融潜热、对流换热及热传导。

1.2 含相变的热传导方程

基本热传导方程为：

式中：为温度；k(k)为()方向的导热系数；q为热源汇项；为储热系数，应考虑相变，引入未冻结含水率w，并有关系：

式中：为土的体积比热容，冻结前后不同；为水的潜热；w为未冻结含水率。

1.3 含相变的热对流理论

当考虑渗流过程中的对流传热时：

式中：c为水的体积比热容；q(q)为()方向的水流量。

1.4 计算参数

1.4.1 土水特征曲线

粗颗粒土的土水特征曲线可采用grain size- Modified Kovacs模型[7]。

式中：为孔隙率；∞为毛细水上升高度，与吸力有关；，和μ为与材料有关的系数；C为矫正系数。

1.4.2 水力传导曲线

水力传导曲线表示土的渗透系数与吸力之间的关系，采用适用于各类土的经验公式的Fredlund模型[8]。

式中：k为非饱和渗透系数；k为饱和渗透系数；为最低负吸力；为最高负吸力；为到的间隔；′(e)为土水特征曲线的斜率。

1.4.3 未冻结含水率

未冻结含水率是指未冻结含水率与空隙中全部水的百分比，用于表达温度在略低于0℃时冻融过程中的热量变化情况。Anderson等[9]提出了一个经验计算公式。

式中：1和为与土的性质相关的常数。

1.4.4 导热系数

饱和状态下的导热系数[10]：

式中：sat为饱和导热系数；K为土颗粒导热系数；K为冰的导热系数；K为水的导热系数；为孔隙率；w为未冻结含水率。

非饱和状态下的导热系数为：

dry为干燥土体的导热系数，可通过经验公式确定：

K为系数，取值为：

为饱和度。

1.5 边界条件

1.5.1 积雪消融入渗边界条件

积雪出水量与环境温度、雪层厚度、辐射及雪的状态有关。根据陈凯[11]对季节性积雪在融雪期的雪盖出流机制的研究结果，积雪出水量的计算方法为：

式中：h为雪层中滞后的剩余水量的深度；W为存储的剩余水量；C3，C4为经验系数；ρ为液态水密度；d为雪层厚度；ρ为雪层中冰的密度。

受环境温度、积雪厚度等影响，不同阶段的积雪消融入渗量不同。本文将积雪消融入渗分为3个阶段：一是积雪初期，降雪量小，温度较高，雪不易积垫，易融化，雪融化后就地入渗，积雪消融入渗量等于降雪量。二是积雪期，本阶段环境温度较低，融雪速度相对较慢，雪融水以入渗土体为主，入渗量影响因素多，不易计算。三是融雪期，融雪期的温度回升，雪融化速度快，出水量大，雪融水往往形成地表径流，积雪入渗量等于最大土体最大入渗量。为了简化计算，按照工程保守原则，将积雪期的雪融水入渗量按照土体最大入渗量取值。

1.5.2 年循环地面温度边界条件

根据赵永峰等[12]在新都桥地区的实测及资料收集，确定年循环地面温度边界条件为：

2 数值模型的建立和验证

2.1 模型建立

本文采用TEMP/W与SEEP/W有限元软件进行季节性粗颗粒冻土边坡的水热数值计算。以川藏铁路新都桥段线路附近的G318国道2 905 km处的季节性粗颗粒冻土边坡为原型建立数值模型。原型边坡的岩土介质为砾土，边坡坡高4 m，坡度42°，地下水在坡脚溢出。按照原型破边建立的数值模型1:1的数值模型，模型尺寸长×宽=32 m×12 m，网格尺寸0.2 m。

首先建立年均地面温度作用下的稳态温度场模型和考略地下水的稳态非饱和渗流模型。在此基础上，建立一个年循环地面温度作用及地下水作用下的瞬态温度−非饱和渗流耦合模型，计算1 a后的某一时刻模型的状态作为初始条件。最后考虑积雪消融入渗，建立一个考虑积雪和渗流含相变的季节性粗颗粒冻土边坡的水热耦合模型。模型计算流程及边界条件见表1。通过实验、工程经验及反复试算，模型计算采用的基本参数见表2。

表1 模型建立过程

表2 模型基本参数

2.2 模型正确性验证

野外监测了边坡坡肩及坡后20 m处(半无限边界)的不同深度的温度。监测元件为PT100，监测2016年至2017年整个冻融循环周期。将2016年至2017年的冻融阶段的温度监测结果与数值模拟的结果进行对比分析，见图1~2。从图中可以看出，计算结果与实测结果的变化趋势基本一致；地下−4 m处数值计算与实测的温度均总体保持在7.8 ℃附近，表明已进入恒温带，季节性温度变化带小于4 m。地下−1.0 m处与−0.2 m处数值计算的温度先降低后升高，实测温度在计算的温度曲线附近波动，差值普遍小于1。计算结果与实测结果之间的误差一方面是由于简化计算模型的精度所致，另一方面是本文研究的是一个冻融循环作用下的边坡温度场问题，将地面温度边界简化为简单的正弦函数曲线，这与实际存在一定的差别。总体而言，数值计算结果能较好地反映边坡实际温度场分布及变化规律，证明了所建模型及参数取值的合理性和正确性。

图1 坡后20 m处不同深度的温度曲线

图2 坡肩处不同深度的温度曲线

3 川藏铁路季节性粗颗粒冻土边坡的水热数值计算

3.1 含水率分布特征

图3为不同深度的含水率分布特征。初始状态(2016−11−01)下，埋深0~2.5 m段处于非饱和状态，含水率随着深度的增加而增加；埋深超过3 m，土体位于地下水位以下，处于饱和状态；埋深2.5~3 m，因毛细力的作用，土体仍处于饱和状态。冻结期前期(2016−11−01~2016−12−09)，在降雪、融雪入渗的作用下，土体含水率逐渐升高。本阶段的降雪量往往较小，且白昼温度大多位于0 ℃以上，融雪速度快，雪不易积垫，往往就地消融入渗，难以形成地表径流。故本阶段的积雪消融入渗量较小，土体含水率提升有限。在冻结期中期(2016−12−09~ 2017−02−23)，积雪量大(期间最大积雪厚度可达30 cm)，加上太阳辐射和白昼气温的影响，融雪速度往往较大，积雪消融入渗量大，导致地表形成暂态饱和区，厚度约为0.8 m。在冻结期后期(2017− 02−23~2017−03−14)，积雪融化速度加快，更高海拔地区的积雪也逐渐加入融化，地表径流流量大，地表仍然保持着暂态饱和区，厚度变化不大。

图3 不同深度的含水率变化曲线

3.2 温度场分布特征

图4~6为温度场云图。初始状态下，温度随深度的增加而升高。冻结前期~冻结中期，浅层地表的温度降低，并逐渐进入冻结状态，从图中可以看出，坡肩处的冻结深度最大，坡脚处最小。坡脚处的冻结深度最小，是因为该处地下水埋深最小，温度较高的地下水热对流作用强。坡肩处的冻结深度大是因为受到温度较低的原始地面和坡面双重作用。进入冻结后期，土体逐渐融化，从图可以看出，最先融化的部位为坡脚，这也是因为该处地下水埋深最小，温度较高的地下水热对流作用强。

图4 2016−11−10(冻结前期)温度场

图5 2017−01−10(冻结中期)温度场

图6 2017−03−10(冻结后期)温度场

图7为不同深度的温度分布特征。初始状态下，土体温度随深度的增大而升高；但当埋深大于3 m时，由于土体进入地下水位以下，土体的温度变化较小。在冻结前期，表层土体温度迅速降低，至2016−12−09已形成冻结层，厚度约为0.2 m。出现冻结层以后，在非冻结−冻结分界线处表现出温度随深度变化曲线的拐点。冻结层内的土体温度随着深度的增大呈近线性增大，非冻结层的温度随深度增大呈对数形式升高。2017−02−23，土体表层温度开始回升，至2017−03−14以后，地面温度达到0 ℃以上，土体开始解冻。至2017−04−02，土体解冻基本结束，地面温度上升至平均气温以上，至2017−04−21，温度回升已导致地下3 m以上所有土体上升到平均气温以上，土体温度随着温度的升高而降低，本次冻融过程结束。

图7 不同深度的温度变化曲线

3.3 冻融深度特征

冻融深度变化特征见图8，图中纵坐标为冻结−未冻结分界线的埋深，横坐标为时间。某一时刻有一个埋深时，这个埋深为冻结深度；某一时刻有2个埋深时，较大的埋深代表冻结深度，较小埋深代表进入融化期后边坡的融化深度。

从图中可以看出，2016年11月底，土体开始冻结，随后冻结深度逐渐增大，到2017−03−11时达到最大，深度约为1.2 m。此时地表已于2017− 03−06开始解冻，即最大冻结深度的时间要稍晚于开始解冻时间约5 d。土体进入解冻期后，地下水及地面温度均高于冻结土体温度，冻结层的底部和顶部同时解冻，至2017−03−31土体全部解冻。土体解冻从2017−03−06~2017−03−31，共25 d；土体冻结从2016−11−26~2017−03−11(2017−03−06开始解冻，但底部的冻结线仍在向下移动)，共约100 d，解冻速度约为冻结速度的4倍。

图8 有无积雪入渗的冻结-未冻结分界线深度曲线

3.4 积雪消融入渗对土体冻融的影响

积雪消融入渗对土体冻融的影响(见图8)主要考虑：1)在冻结初期，雪融水的温度低于土体温度，积雪消融入渗带走土体热量，有助于土体冻结。 2) 在冻结后期，土体开始融化，冻结的土体温度低于雪融水的温度，积雪消融入渗给土体带来热量，有助于土体融化。对比有积雪消融入渗和无积雪消融入渗的季节性粗颗粒冻土的冻融深度，可以看出：1) 有积雪消融入渗的最大冻结深度为1.2 m，无积雪消融入渗为0.75 m，有积雪消融入渗的冻结深度是无积雪消融的1.6倍。2) 有无积雪消融入渗的边坡土体最大冻结深度的时间不同，无积雪消融入渗时，在外界温度最低点时，土体冻融深度达到最大，冻结层底部边界线保持稳定，直至地面温度回升到0 ℃上以后约25 d才快速上升；有积雪消融入渗时，在外界温度回升到0 ℃之后5 d，土体冻结深度达到最大，随后冻结层底部边界线快速上升。这反映了外界温度达到最低温度以后，虽然地面温度仍低于0 ℃，但土体冻结层的底部变化难以向下移动；但如果此时有积雪消融入渗的作用，热对流将大大提高冻结能力。3) 冻结深度随时间的曲线斜率反映了土体的冻结速率和融化速率。有无积雪消融入渗的冻结速率分别为1.2/100=0.012 m/d和0.75/80=0.009 4 m/d，有积雪消融入渗的冻结速率是无积雪消融入渗的1.3倍。有无积雪消融入渗的融化速率分别为1.2/20=0.06 m/d和0.75/35=0.02 m/d，有积雪消融入渗的冻结速率是无积雪消融入渗的3倍。

3.5 地下水水位对冻结深度的影响

地下水作为热源热对流边界，往往抑制土体的冻结，加速土体的融化。对比地下水坡脚出露和地下水埋深大(坡脚处埋深8 m)的季节性粗颗粒冻土边坡的冻融深度，见图9，可以看出：1)坡脚地下水出露边坡的最大冻结深度为1.2 m，地下水深埋时为1.9 m，前者是后者的63%。2) 冻结深度随时间的曲线斜率反映了土体的冻结速率和融化速率。地下水坡脚出露和地下水深埋时边坡的冻结速率分别为1.2/100=0.012 m/d和1.9/125=0.015 2 m/d，前者是后者的79%。地下水坡脚出露和地下水深埋时边坡的融化速率分别为1.2/20=0.06 m/d和1.9/65 =0.03 m/d，前者比后者增大1倍。

图9 不同地下水深浅的冻结-未冻结分界线深度曲线

4 结论

1) 积雪消融入渗提高了季节性粗颗粒冻土边坡热交换的速度，增强了冻结能力。本文计算结果为季节性粗颗粒冻土边坡的冻结深度增大60%，冻结速率增大30%，融化速率增大200%。

2) 地下水作为热对流边界，会抑制土体冻结，加速土体融化。本文计算结果为坡脚地下水出露边坡的冻融深度为地下水深埋边坡的63%，冻结速率为79%，融化速率增大1倍。

3) 川藏铁路新都桥地区季节性粗颗粒冻土边坡的冻结深度为1.2 m，最大可达到1.9 m。边坡不同位置的冻结深度不同，在坡肩处最高，坡脚处最低；进入融化期后，因冻结深度小及地下水热流作用，坡脚处最先融化。

4) 季节性粗颗粒冻土边坡白昼时积雪出水量大，产生积雪消融入渗，在地表形成暂态饱和区，厚度约为0.8 m。

[1] 张玉芝, 杜彦良, 孙宝臣, 等. 季节性冻土地区高速铁路路基冻融变形规律研究[J]. 岩石力学与工程学报, 2014, 33(12): 2546−2553. ZHANG Yuzhi, DU Yanliang, SUN Baochen, et al. Roadbed deformation of high-speed railway due to freezing-thawing process in seasonally frozen regions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(12): 2546−2553.

[2] 张玉芝, 杜彦良, 孙宝臣. 季节性冻土地区高速铁路路基地温分布规律研究[J]. 岩石力学与工程学报, 2014, 33(6): 1286−1296. ZHANG Yuzhi, DU Yanliang, SUN Baochen. Temperature distribution in roadbed of high-speed railway in seasonally frozen regions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014(12): 2546− 2553.

[3] 张媛, 董建华, 董旭光, 等. 季节性冻土区土钉边坡支护结构冻融反应分析[J]. 岩土力学, 2017, 38(2): 574− 582. ZHANG Yuan, DONG Jianhua, DONG Xuguang, et al. Analysis of freezing and thawing of slope improved by soil nailing structure in seasonal frozen soil region[J]. Rock and Soil Mechanics, 2017, 38(2): 574−582.

[4] 习阿幸, 刘志辉, 卢文君. 干旱区季节性冻土冻融状况及对融雪径流的影响[J]. 水土保持研究, 2016, 23(2): 333−339. XI Axing, LIU Zhihui, LU Wenjun. Processes of seasonal frozen soil frezing-thawing and impact on snowmelt runoff in arid area[J]. Research of Soil and Water Conservation, 2016, 23(2): 333−339.

[5] 吴青柏, 蒋观利, 蒲毅彬. 青藏高原天然气水合物的形成与多年冻土的关系[J]. 地质通报, 2006, 25(2): 29−33. WU Qingbai, JIANG Guanli, PU Yibin. Relationship between permafrost and gas hydrates on Qinghai-Tibet Plateau[J]. Geologial Bulletin of China, 2006, 25(2): 29−33.

[6] Black P, Tice A. Comparison of soil freezing curve and soil water curve data for windsor sandy loam[J]. Water Resources Research, 1989, 25(10): 2205−2210.

[7] Aubertin M, Mbonimpa M, Bussière B, et al. A model to predict the water retention curve from basic geotechnical properties[J]. Canadian Geotechnical Journal, 2003, 40(6): 1104−1122.

[8] Fredlund D G, XING Anqing. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994(31): 521−532.

[9] Johansen O. Thermal conductivity of soils[D]. Trondheim, Norway. (CRREL Draft Translation, 637, 1977), ADA044002.

[10] Anderson D M, Tice A R. The unfrozen interfacial phase in frozen soil water systems, physical aspects of soil water and salt in ecosystems[J]. In Ecological Studies Springer Verlag, Berlin, 1973(4): 107−124.

[11] 陈凯. 融雪水在雪层中的冻融过程及侧渗出流研究[D]. 乌鲁木齐: 新疆大学, 2012. CHEN Kai. Research on freeze-thaw process and side flow of snowmelt in the snow[D]. Urumuqi: Xinjiang University, 2012.

[12] 赵永峰, 谢强, 王子江, 等. 川藏铁路季节性冻土区粗颗粒土边坡温度场特征[J]. 铁道建筑, 2017(5): 95−99. ZHAO Yongfeng, XIE Qiang, WANG Zijiang, et al. Characteristics of temperature field of coarse-grained soil slope in seasonal frozen soil region along Sichuan-Tibet Railway[J]. Railway Engineering, 2017(5): 95−99.

Hydrothermal process simulation of seasonal coarse-grain frozen soil slope on Sichuan-Tibet railway

LIANG Shu1, XIE Qiang1, ZHAO Wen1, YE Tangjin1,2, JIANG Liangwen3

(1. Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 610031, China;2. Tibet University Institute of Technology, Lhasa 850000, China; 3. China Railway Eryuan Engineering Group Co., Ltd, Chengdu 610032, China)

This article based on the seasonal coarse-grained frozen soil slope in the eastern section of the Sichuan-Tibet Railway, a transient hydrothermal coupled saturated-unsaturated porous multiphase media numerical model with phase transition was established. This model considered snow cover, seepage and annual cycle temperatures. The temperature field of seasonal coarse-grained frozen soil slope, moisture distribution characteristics, freezing depth and its influencing factors were analyzed combing with field measurements. The results show that the removal and infiltration of snow may change the moisture field of seasonal coarse-grained frozen soil. The surface can form a transient saturation area of up to 0.8 m. The dynamic changes of the moisture field cause the increase of heat transfer rates, the enhancement of freezing ability. It also causes slope freeze depthincreasing by 60%, the freezing rate by 30% and the melting rate by 200%. The soil freezing speed is suppressed by the boundary of the groundwater convection. The frozen-thawing depth of the outcrop slope at the slope foot is 63% of the deep-buried slope of the groundwater, the freezing rate is 79% and the melting rate is increased by 1 time. The freezing depth of seasonal coarse-grained frozen soil slope in the Xinduqiao area of the Sichuan-Tibet railway is 1.0 m, and the maximum reachable is 1.9 m. The different depths of freezing in different positions of the slope are highest at the shoulders and lowest at the foot. After entering the melting period, due to the small depth of freezing and the effect of groundwater heat flow, the first step of melting is at the foot of the slope.

Sichuan-Tibet railway; seasonal frozen soil; the seasonal coarse-grain frozen soil; freeze-thaw depth; hydrothermal coupling; numerical simulation

U213.1

A

1672 − 7029(2019)06− 1440 − 08

10.19713/j.cnki.43−1423/u.2019.06.012

2018−05−22

国家自然科学基金资助项目(41662020)；中国铁路总公司科技研究开发计划课题(2017G008-F)

赵文(1975−)，男，重庆人，副教授，从事特种岩土工程、边坡工程治理的研究；E−mail：civil1301@swjtu.edu.cn

(编辑 涂鹏)