不同冻土条件下渠道地基热状况模拟研究

韩洪武，穆彦虎，虞 洪，丁泽琨，陈 领

(1. 青海省水利水电工程局有限责任公司，青海 西宁 810001； 2. 中国科学院西北生态环境资源研究院 冻土工程国家重点实验室，甘肃 兰州 730000； 3. 中国科学院大学，北京 100049)

多年冻土指地表以下一定深度内地温在0 ℃及0 ℃以下至少连续存在两年的岩土层。我国为世界第三冻土大国，多年冻土分布面积约占国土面积的22%，主要分布在青藏高原、大小兴安岭及西部高山地区[1]。作为一种客观特殊地质体，多年冻土是地质历史和气候变迁背景下受区域地理环境、地质构造、岩性、水文和植被特征等因素共同影响下，通过地气间物质与能量交换发育而成的。从物质组成来看，冻土是由矿物颗粒、冰、未冻水和气体组成的一种四相体材料[2]。与常规岩土不同，冻土的物理、力学、水力性质及伴随冻结和融化过程中的热、力、变形现象均与其温度变化过程密切相关[3]。

在多年冻土区进行工程建设，工程活动对包括局地地形、地貌、植被、浅层岩土体类型、地表水体、径流等局地因素的改变将不可避免地打破原有的地气能量平衡，进而影响到下伏多年冻土水热状况，引发多年冻土的退化过程，具体可表现为活动层厚度增加、地温升高、融化夹层发育[4-7]。对于渠道工程，涉及的开挖工程和过水等将对活动层季节冻融过程和下部多年冻土水热状况产生长期且显著的影响，而这一过程反过来将大幅降低冻土地基承载力，引发显著的沉降变形过程，由此威胁到水工结构物的长期稳定性，甚至结构功能失效[8]。目前，受工程实践需求所限，国内外冻土区渠道相关研究主要集中在季节冻土区，包括渠基土的水热监测与模拟、衬砌结构的冻胀机理与过程、适应冻胀的衬砌结构优化与措施研发、循环冻融作用对渠基土、混凝土等材料工程性能的影响等方面[9-14]。而在多年冻土区，相关方面的研究十分有限，主要集中在俄罗斯、加拿大等高纬度多年冻土区的一些水利大坝的设计、建设及运营过程中坝体的热状况、冻融界限及包括冻胀、融沉等坝体变形的监测和模拟[15-19]。

以高海拔多年冻土区某渠道工程为研究对象，开展渠道工程和多年冻土地基相互作用研究。针对长期过水及气候变暖影响，在考虑冻融土体内水分迁移、冰水相变及土体未冻水含量与温度非线性关系基础上，构建冻融土体水-热耦合数学模型，利用数值模拟手段开展渠道下伏冻土地基长期热状况演化过程模拟研究，并通过已有现场监测数据对模型进行验证。

1 水热耦合计算模型

为研究渠道下部多年冻土地基长期热状况演化规律，需要建立冻融土体水热耦合数学模型。模型建立过程中的基本假设包括：土体均匀、连续且各向同性；冻结、融化土体内的水分迁移均满足达西定律；考虑水势梯度作为冻融土体中液态水迁移的驱动力，忽略温度梯度的影响。

1.1 质量守恒方程

根据质量守恒定律，在不考虑水汽迁移的情况下，冻土中总等效体积含水量θ可以写成：

式中：θ为总等效体积含水量（m3/m3）；θu为液态水体积含量（m3/m3）；θi为体积冰含量（m3/m3）；qlh为土壤水势梯度下液态水通量密度（m/s）；ρi和ρw分别为冰密度和水密度（kg/m3）；t为时间。

根据Harlan模型，冻土中的液态水迁移与非饱和融土中液态水迁移相似[20]。冻土中的液态水迁移可以用Richards方程来描述，且主要受水势梯度和导水系数的影响[21]。在未冻土和冻土中，水势梯度和导水系数有差异显著，但冻结土体中水分迁移仍可认为遵循达西定律[22]。因此，只考虑水势梯度作为冻结土体和未冻土体中液态水迁移的驱动力，则液态水的通量密度[23]可以写成:

式中：Klh为水势梯度下土体液态水导水系数（m2/s）；y为纵向坐标（m）；h为压力水头（m）。考虑将总等效体积含水量作为变量，冻融土体中液态水质量守恒方程可以写为：

1.2 水力参数

冻融土体中液态水迁移与其压力水头相关，可通过土水特征曲线（SWCC）计算压力水头。本文中，冻融土体的水力参数可以用van Genuchten模型和Mualem模型[24-25]进行描述：

式中：Se为有效饱和度；Ks为饱和土体导水系数（m/s）；θl、θs、θr分别为液态水含量、饱和液态水含量和残余含水量（m3/m3）；α为土体进气值的倒数（m−1）；m= (1 – 1/n)、n、l为经验参数，文献 [25]建议l取 0.5。

1.3 能量守恒方程

已有试验研究表明，与热传导和冰水相变潜热相比，土体中对流传热量级非常小，在冻融土体的传热分析中一般可以忽略[21]。因此，仅考虑热传导和冰水相变潜热，则冻融土体的能量守恒方程[26]可以写为：

式中：T为温度；Cm、λm分别为土体等效体积比热容（J/m3/K）和等效导热系数（W/m/K）；L为冰水相变潜热（3.34×105J/kg）。根据显热容法，土体等效体积比热容和等效导热系数可以写为：

式中：Tf± ΔT为相变区间；Cu和λu分别为冻结温度以上土体比热容和导热系数；Cf和λf分别为冻结温度以下土体比热容和导热系数；Ls为含水岩土介质单位体积相变潜热。

1.4 联系方程

对于冻结区，以上质量守恒、能量守恒两个方程涉及温度、体积含冰量、液态水含量3个未知量，未知量大于方程数，需要增加联系方程。土体冻结后，未冻水的含量取决于温度、压力、水盐度、矿物学、土壤比表面积和土壤表面化学成分等因素。基于已有理论和试验研究结果[21]，采用下式来表示土体冻结过程中的最大未冻水含量：

式中：a和b为土体试验参数。体积液态水含量和含冰量可通过等效体积含水量和温度确定：

式中：Tf为土体冻结温度(℃)。已有研究表明，土体的冻结温度并不是一个定值，并且只有土体中的液态水大于最大未冻水含量时才会有冰生成。

2 数值计算模型及边界条件

2.1 数值计算模型

以高海拔多年冻土区某渠道为研究对象，工程区平均海拔4 500 m以上，区域多年平均气温为−5.1 ℃，年平均降水量313.8 mm，属于大片连续多年冻土区。根据现场钻探及相关资料，区域内活动层厚度（季节融化层深度）在2~3 m，岩性大部分为含砾砂土及砾砂，呈散体状。多年冻土含冰量特征以富冰冻土为主，局部为饱冰-含土冰层，多分布于多年冻土上限至5 m深度范围。

考虑到渠道结构对称性，以渠道中心线为对称轴，取渠道结构一半建立物理模型。根据现场工程地质勘查结果和渠道设计方案，渠基土自上而下简化为3层，即砂砾土、粉质黏土和强风化泥岩，建立数值计算模型如图1所示。图1中，渠道底面宽度为25 m，渠道边坡按1∶3放坡，坡长15 m，渠内水深依据目前实测水深设置为1 m。考虑到渠道工程沿水流方向包括结构尺寸及气温和水温差异性不大，渠道地基长期热状况的模拟可以简化为二维非稳态温度场计算问题。

图1 渠道物理模型（单位：m）Fig. 1 Physical model of the canal (unit: m)

结合现场钻孔岩芯实测干密度和含水量，以及室内冻土热物理参数研究结果[21-22, 26-27]，表1给出了3层土体的热物理参数。

表1 土体热物理参数Tab. 1 Physical-thermal parameters of soil layers

结合现场工程地质勘查和多年冻土地温实测数据，渠道沿线多年冻土地温分布差异显著，10~15 m深度年平均温度（TMAGT）最高为−0.3 ℃，最低可达−1.5 ℃。依据《青藏铁路多年冻土勘测暂行规定》中关于冻土含冰量的划分，渠道沿线冻土体积含冰量（iv）情况则包括少冰（iv≤10%）、多冰（10%

表2 不同含冰类型土的热物理参数Tab. 2 Physical-thermal parameters of frozen soils with different ice contents

表3 不同含冰类型土的视比热容Tab. 3 Apparent specific heat of frozen soils with different ice contents 单位：J·kg−1·℃−1

2.2 边界条件

为考虑气候变暖影响，收集整理了距工程区最近气象站数据。结果表明，自1961年至2018年期间工程区气候呈现显著的暖湿化趋势，气温上升速率为0.33 ℃/(10 a)。根据附面层理论、现场观测和工程区气温上升速率，将数值模型中边界LC和CE设置为天然地表温度边界：

从实际传热过程看，渠道过水截面部分为第三类对流换热边界。然而，受限于边界对流换系数选取困难，目前相关模拟研究中，常常以实测水体温度作为边界条件，即采用第一类边界，在简化问题的同时确保计算结果的可靠性[28]。为此，本文在渠道过水截面部分采用第一类边界条件。已有研究表明，湖泊水体温度与环境气温存在良好的线性关系[29]：

式中：Tw(t)、Ta(t)为同一时间尺度上的水温和气温；A和B为常数。

根据林战举等[30]对工程附近区域多个湖泊的实际监测水温和气象站气温数据，可确定式（13）中的A、B值。

在不考虑地下水输入、人工热输入和防风设置等因素影响时，同一地区A、B值差异不大[29]。因此，根据工程区气象站气温监测数据，并考虑区域气温升温速率，可将边界AB和BL设定为渠道水温边界：

模型中KG为对称边界，EF为绝热边界。底边GF为热流边界，根据现场测温所得该深度地温梯度和土层热物理参数确定热流值为0.03 W/m2。

渠道施工过程中，浅层土体的开挖导致下伏多年冻土直接暴露在空气中，势必对多年冻土的热状况产生一定影响。然而，相较于后期长期运营过程，施工过程持续时间短，且所产生的热扰动相对较小，因此在本文模拟预测中未考虑施工开挖过程对下伏多年冻土的热影响。

3 模拟结果及分析

3.1 模型验证

图2给出了渠道开挖1年后天然地面不同深度地温的数值模拟结果与现场监测结果的对比情况。由监测结果可以看出，天然地面冻土上限保持在2.1 m左右。数值模拟结果显示，最大季节融化深度时的多年冻土上限同样在2.1 m左右，与现场监测结果接近。自天然地面以下15 m深度处，现场监测结果和数值模拟结果均显示多年冻土地温为−0.3 ℃。从两者的对比看，本文数学模型及所采用的热物理参数和边界条件能够较好地反映该地区多年冻土热状况，包括活动层厚度和下伏多年冻土地温均表现出较好的一致性。但从活动层范围内及多年冻土顶板附近的地温对比看，模拟结果和实测结果仍存在一定的差异，这些差异与数值模型中土层的简化及实测过程中测点深度的确定有关。

图2 开挖1年后天然场地不同深度地温实测与模拟对比Fig. 2 Field observed and numerical simulated temperature profiles at natural ground after the canal construction

3.2 不同年平均地温条件下渠道地基温度场

图3给出了少冰冻土条件下，年平均地温分别为−0.5、−1.0和−1.5 ℃时渠道开挖后第50年10月15日地基温度场分布云图。从图3（a）可以看出，当年平均地温为−0.5 ℃时，开挖后第50年10月15日天然地表处的最大季节融化深度为3.5 m，相较于初始最大季节融化深度2.1 m，增加了约1.5 m，年均增加速率为3 cm/a。这一量值与目前青藏高原活动层增加速率监测结果接近[31]。当年平均地温为−1.0、−1.5 ℃时，两者的最大季节融化深度增速较年平均地温为−0.5 ℃条件下小，反映了气候变暖条件下不同年平均地温多年冻土退化响应的差异[32]。

图3 不同年平均地温条件下渠道开挖后第50年10月15日地基温度场Fig. 3 Temperature fields of permafrost subgrades with different MAGTs on October 15, 50 years after the canal construction

从图3可以看出，受渠道过流影响，50年后渠底处、坡脚处及天然地面靠水渠的部分下部30 m范围内多年冻土已经退化，反映了渠道强烈的热侵蚀作用，且这种热侵蚀在垂向和水平向均比较显著。多年冻土的退化势必引发渠道底部和岸坡的差异融沉变形，导致渠道结构和形态的变化。不同年平均地温条件下，渠道的侧向热侵蚀存在显著差异，当多年冻土地温为−0.5 ℃时，自岸坡向外约10 m范围内下部多年冻土已退化，而当多年冻土地温为−1.0和−1.5 ℃，岸坡下部仍有多年冻土分布。这一结果体现了不同年平均地温多年冻土热惰性的显著差异。靠近渠道坡脚部分多年冻土的退化和岸坡下显著的温度梯度，可导致岸坡的滑塌和失稳，并由此引发渠道形态的改变。

3.3 不同含冰量条件下渠道地基温度场

图4给出了年平均地温为−0.5 ℃时，3种含冰量条件下渠道开挖后第50年10月15日地基温度场分布云图。可以看出，3种含冰量条件下渠道地基温度场差异显著。当渠道下部多年冻土含冰量为少冰时，渠道下部以及自坡肩向外约10 m范围内多年冻土已经退化。随着含冰量的增加，多年冻土热惰性显著增加，因此其退化速率明显减小。与不同地温条件下的对比（图3）可以发现，相较于地温，多年冻土的含冰量对气候变暖和人为扰动背景下多年冻土退化响应影响更加显著。这也反映了含冰量作为冻土工程地质条件的关键指标，决定着冻土工程地基热力稳定性和基础与结构物长期稳定性。

图4 不同含冰量下渠道开挖后第50年10月15日地基温度场Fig. 4 Temperature fields of permafrost subgrades with different ice contents on October 15 of the 50th year after the canal construction

当多年冻土含冰量为多冰和富冰时，渠道运营50年后其下部仍有多年冻土存在。但受气候变暖和过水影响，其季节融化深度显著增加，渠道下部形成了一个“锅底状”的融化盘。当多年冻土为多冰冻土时，渠道中心锅底状融化盘深度可达渠底以下约20 m，而当多年冻土为富冰冻土时，融化盘深度为渠底以下约15 m。含冰量不同时，岸坡外多年冻土的退化过程也存在明显不同，包括季节融化深度和多年冻土上限以下地温分布，含冰量越高，其季节融化深度越小，随时间增加过程越缓慢。

3.4 渠道典型位置地温状况时间演化过程

从以上分析可以看出，在气候变暖和渠道过水影响下，渠道不同位置多年冻土的退化过程和速率明显不同。为分析渠道不同位置地温状况时间演化过程，以年平均地温为−0.5 ℃，含冰量为多冰冻土为例，图5为20年内渠道中心、坡脚、坡肩和天然地面等4个位置地温-深度曲线。为节约篇幅，图5中仅给出了渠道运营第5、10、15和20年10月15日最大融化深度时的地温-深度曲线。

图5 渠道开挖后20年内渠道不同位置地温曲线Fig. 5 Temperature profiles at different locations within 20 years after the canal construction

从图5（a）可以看出，渠道中心受渠道内流动水体影响较大，第5、10、15和20年10月15日最大融化深度分别为5.2、8.1、10.2和12.2 m，年均增加速率为0.4~0.6 m/a。同时可以看出，最大季节融化深度的增加主要发生在前5年，这与渠道开挖的直接影响密切相关。整体来看，在渠道开挖和过流影响下，渠道下部多年冻土主要表现为自上而下的退化过程，即多年冻土上限下移显著而下部多年冻土地温升温不显著，这是由于在−0.5 ℃地温条件下，多年冻土已经接近剧烈相变区。相对于渠中，坡脚的多年冻土退化速率较慢，但其退化模式基本相同（图5（b））。岸坡和天然场地多年冻土的退化主要受气候变暖影响（图5（c）和（d）），尤其在前20年内，渠道内过水侧向热侵蚀尚未发展到岸坡位置。可以看出，仅在气候变暖影响下，两个位置的多年冻土退化相对缓慢，且退化模式均表现为多年冻土上限的下移和上限以下一定深度范围内多年冻土的升温。通过4个位置地温的对比，可以很好地反映开挖过水因素和气候变暖因素作用下多年冻土退化过程的差异性。

4 结 语

在多年冻土区，与道路、房建、输油气管道等工程有所不同，水利工程涉及大量的挖方工程和流动水体，因此对下伏多年冻土的扰动更加显著。以高海拔多年冻土区某渠道工程为例，利用数值模拟手段，研究了过水和气候变暖对下部多年冻土热状况的影响，并考虑多年冻土地温、含冰量两个重要因素，模拟分析了气候变暖背景下渠道下部多年冻土地温的时间演化过程，主要结论如下：

（1）在气候变暖背景下，天然场地多年冻土的退化主要表现为活动层厚度增加和多年冻土地温升高。以区域历史气象增温率模拟气候变暖过程，当多年冻土年平均地温为−0.5 ℃，含冰量为少冰时，天然场地活动层年均增加速率为3 cm/a，这一量值与目前实测青藏高原活动层增加速率接近。

（2）受渠道开挖和过水影响，渠道底部、坡脚及岸坡位置多年冻土退化显著。渠道运营50年后，当多年冻土含冰量为少冰时，−0.5、−1.0和−1.5 ℃三种多年冻土地温条件下，渠底、坡脚及天然地面靠水渠的部分下部30m范围内多年冻土已经退化，反映了渠道在垂向和水平向强烈的热侵蚀作用。不同年平均地温条件下，多年冻土热惰性不同，因此渠道的侧向热侵蚀存在一定差异。当多年冻土地温为−0.5 ℃时，自坡肩向外约9.8 m范围内下部多年冻土已退化；而当多年冻土地温为−1.0和−1.5 ℃，岸坡下部仍有多年冻土分布。

（3）含冰量是影响多年冻土退化速率的重要因素。渠道下部多年冻土含冰量由少冰增加至多冰、富冰时，其热惰性显著增加，气候变暖和过水双重作用下多年冻土退化速率明显减小。即使在年平均地温为−0.5 ℃时，渠道运营50年后其下部仍有多年冻土存在。但是，受气候变暖和过水双重影响，渠道下部季节融化深度显著增加，形成了一个“锅底状”的融化盘，融化盘大小与含冰量有关。

（4）渠底、坡脚、岸坡及天然场地多年冻土地温-深度的时间变化过程对比表明，在过水和气候变暖双重作用下，渠底和坡脚多年冻土的退化主要为多年冻土上限的快速下移，即自上而下的退化模式；而在气候变暖单独作用下，岸坡和天然场地多年冻土的退化表现为活动层缓慢下移的同时，上限以下一定升幅范围内多年冻土的缓慢升温。