寒区公路路堑边坡开挖水热耦合模拟

2022-11-01石振武李博梁宸晨

科学技术与工程 2022年26期

石振武，李博，梁宸晨

(东北林业大学土木工程学院，哈尔滨 150036)

为满足寒区高等级公路对线形要求，对部分沿线土体进行开挖，挖方取土破坏了地表植被，切断了地下水径流，导致路堑边坡坡体应力重分布及浅层热平衡变化[1-2]，坡体内水分向冻结区迁移，在新形成的路堑边坡坡面出现地下水露头，水分迁移导致土体水分场重分布，使土体的强度发生劣化，造成边坡侵蚀和滑塌灾害现象[3]。

边坡土体的水热传递过程受到降水、温度、风速、融雪、辐射、地下水渗流等多方面的影响[4]，是一个带相变的渗流场、温度场复杂问题[5-6]。在各项不同气候条件的作用下，水的基质吸力、土体的导热系数、比热容均发生改变，加上相变潜热，在边坡开挖的条件下，土体水热变化更加复杂。

目前，中外学者对土体水热变化研究已经进行了大量水热耦合模拟[7-9]。刘宏超等[10]基于K最近邻算法(K-nearest neighbor, KNN)的机器学习方法进行土壤水热状况的模拟，表明KNN模型能够较为精确地对青藏高原多年冻土区唐古拉地区表层土壤水热状况进行表征。刘扬李等[11]构建分布式水循环模型进行水循环过程模拟，结果表明模型模拟接近实测。国外学者也对寒区的融雪、冻土等水文过程进行模拟[12-13]，模拟结果较为理想。

大部分学者进行水热状况模拟多针对青藏高原冻土区，对于寒区绥化地区的研究相对较少，部分学者对黑土地区水热变化进行现场监测[14]，但在路堑边坡开挖前后，土体的水分场和温度场不易被监测，对于开挖前后的监测数据及相关地质资料不能够很好地描述边坡土体的水热状况。通过实测的气候条件构建土体边坡，进行开挖前后的水分场与温度场的耦合，可以较为直观地反映出各阶段土体边坡的水热发育状况[15]。因此，在实地勘测的基础上，现采用Geo-studio软件结合气象数据对开挖前后整个时期边坡的水分场、温度场的分布特征进行耦合分析，以此为正确判断边坡破坏因素、进行稳定性分析及制定防护措施方案提供参考。

1 渗流及对流热耦合理论

1.1 饱和/非饱和渗流理论

有限元渗流控制方程为

(1)

式(1)中：kx、ky分别为x、y方向的渗透系数；Q为施加的边界流量；t为时间；H为总水头；H=uw/γw+y；uw为孔隙水压力；γw为水的容重；y为位置水头；mw为储水曲线斜率。

土的冻结曲线采用土水特征曲线和Clapeyron方程来估计，Black等[16]建立土冻结时基质吸力与负温的线性关系，即

Δμw=-1 110ΔT

(2)

式(2)中：Δμw为基质吸力的变化量；T为低于0 ℃土的温度。

1.2 对流热耦合理论

二维对流热耦合控制方程为

(3)

式(3)中：ktx、kty分别为x、y方向的导热率；Qt为施加的边界通量；T为温度；ρwcpw、ρscps分别为水、土的体积热容；qwx、qwy为x、y方向符合达西定律的水流量；L为水的潜热；w为土的体积含水量；Wu为未冻结含水率。

2 模型建立及验证

2.1 工程概况

对绥大高速绥化至青冈段(K26+840-K27+400)部分土体进行初步开挖，路堑边坡出现季节性水体流出，导致坡面失稳滑塌，如图1所示。

图1 绥大高速边坡失稳Fig.1 The slope of the Suida expressway is unstable

2.2 模型建立

经实地勘测，该标段边坡为二级边坡，每级坡高5 m，坡比1.2，坡中平台宽1 m，上部土体为黑色淤泥质粉质黏土(黑土)，厚度约2 m，下部土体为黄土状亚砂土(黄土)；坡面失稳滑塌多出现在下级边坡，并伴有水体流出。根据实际工程建立长×高=28 m×15 m，网格尺寸0.5 m的数值分析模型，采用Geo-studio的SEEP/W、TEMP/W模块对模型的水分场、温度场变化规律进行耦合分析，模型及相关监测位置如图2所示。

图2 路堑边坡数值模型Fig.2 Numerical model of the slope of the graben

为分析不同土体边坡的水热状况，构建如表1所示的4种工况，进行差异性研究。

表1 所模拟的工况Table 1 The simulated operating condition

气候条件与气温、降水量、相对湿度、风速及净辐射能量有关，原始数据来源于国家科技基础条件平台-国家地球系统科学数据中心(http://www.geodata.cn)、黑龙江海伦农田生态系统国家野外科学观测研究站，多年平均的日均气象数据如图3和图4所示。

图3 日均降水量、气温与相对湿度Fig.3 Average daily precipitation, air temperature and relative humidity

图4 日均风速与净辐射通量Fig.4 Daily average wind speed and net radiation flux

日均地面温度数据与积雪消融数据根据原始数据进行处理。土地气候的相互作用采用Penman-Monteith[17]蒸发蒸腾模型进行表达，即

(4)

式(4)中：λET为潜热通量；Δ为饱和水汽压关于温度的曲线的斜率；Rn为净辐射通量；G为土的热通量；ρa为空气密度；Cp为空气定压比热容；es为饱和水汽压；ea为实际水汽压；γ为湿度计常数；ra和rc为两个关键参数。

首先建立稳态热及稳态渗流模型，作为瞬态水热耦合的基础，再根据地面温度及气候条件建立2000年1月1日—2020年3月1日的瞬态水热耦合模型，使水热状况充分发育，获得开挖前土体水热分布特征，在此基础上，进行2020年3月2日—12月31日的开挖后边坡瞬态水热耦合模拟。模型建立流程及边界条件如表2所示，模型参数如表3所示。

表2 模型建立Table 2 Establishment of the model

表3 模型参数Table 3 Parameters of the model

2.3 模型验证

将数值模拟结果与2015—2016年的气象站冻结深度监测数据进行对比，如图5所示。

从图5可以看出，模拟结果与气象监测数据变化趋势基本一致，表层土壤都于11月初开始冻结，3月末出现融化，最大冻结深度为1.92～1.95 m。但模拟结果对气候变化反应更为敏感，一方面是由于边界条件简化，日均气象数据边界与实际气候条件有差异，另一方面是由于气象站监测土体与边坡现场状况有不同。总体而言，数值模拟结果能够较好地反映土体冻结深度的变化规律，证明所建模型及参数取值的合理性及正确性。

图5 冻结深度对比曲线Fig.5 Contrast curve of freezing depth

3 结果与分析

3.1 水分场分布特征分析

图6为边坡开挖前体积含水量分布特征。可以看出，深度3.5 m体积含水量最低为5%，深度2.5～4.5 m土体体积含水量小于10%，是由于浅层黑土的渗透系数及水力传导率较小，渗入下部黄土的水分较少，黄土渗透系数较高、导水性较好，上部渗入水分小于下部土体下渗水分，体积含水量较低。表层土体水分随气候变化性较强，冬季寒冷干燥，部分液态水转化为冰体，另一部分在水势梯度的影响下发生迁移，浅层深度0～1.5 m体积含水量最终稳定在20%～22%；深度1.5～2.5 m土体受上部水分入渗影响，体积含水量下降速率最快为8%/m；深度4.5～6 m土体处于饱和与非饱和状态的过渡段，是一种暂非饱和状态，气温升高，上层冰体转化为液态水，地下水位回升，土体为饱和状态；埋深超过6 m，土体位于地下水位以下，为饱和状态。

图6 土体开挖前体积含水量分布状态Fig.6 The distribution of volumetric moisture content before excavation

边坡开挖，水分场发生改变。起初，受零下低温影响，表层水发生相变，冻结锋面自上向下逐步推进，暖端土体水分向锋面迁移，发生重分布。绘制监测点体积含水量变化曲线如图7所示。

图7 监测点体积含水量变化曲线Fig.7 Variation curve of volume water content in monitoring points

可以看出，开挖后上部监测点1、2、3受负温影响，土体冻结，体积含水量在10～15 d内稳定；监测点4体积含水量短暂升高，是由于地下水向冻结峰面迁移，土体冻胀；监测点5处于地下水位以下，体积含水量变化滞后，当水体迁移导致地下水位下降时，其发生波动，并大幅度降低。

气温上升，冻结土体边坡融化，体积减小，孔隙率增大，体积含水量降低，随着水体从土层中涌出，土中细颗粒减少，孔隙率进一步增大，边坡地下水位快速下降，体积含水量在3月中下旬骤减，监测点4、5最为显著，在4月初，已下降至25.39%、41.15%，在5月中旬进入缓慢下降趋势，呈线性平稳。11月初土体开始冻结，各监测点体积含水量保持稳定，在非冻结开挖期内分别下降18.26%、32.79%、70.78%、81.46%、78.62%。

3.2 温度场分布特征分析

图8为边坡开挖前温度场分布特征，深度越深温度越高。在冻结期内，随着环境温度的降低，表层土体的直接表征是冻结现象的发生，随之深层土体的温度不断下降。表层土体受大气温度及太阳辐射的影响，温度变化波动较大；深层土体受上部土体热传导及热对流的影响，温度变化波动较小，并随着土层深度的不断增加，温度变化幅度逐渐缩小，变化趋势趋于稳定，深度3 m土体温度达3.2 ℃。

图8 边坡开挖前温度分布状态Fig.8 Temperature distribution before slope excavation

边坡开挖，土体温度分布形式被打破，受负温影响，坡面冻结。为描述开挖坡面温度变化特征，绘制坡上部及坡中平台同一深度(0 、1 、2 、3 m)的温度变化曲线如图9所示。

从图9可知，坡上部温度成波动变化，且随着埋深的增加，波动范围呈指数衰减，分别为：-25.73～31.58 ℃、-5.02～13.94 ℃、0.57～8.19 ℃、2.05～5.30 ℃；且波动的滞后性凸显，深度0 m与1、2、3 m间的最低温度间隔分别为56、116、166 d。开挖后坡中平台埋深减少，温度的波动变化显著，同时也存在滞后性，但弱于坡上部。

图9 温度变化曲线Fig.9 Temperature change curve

对比同一深度坡上部与坡中平台，埋深越浅，二者相关性越强。开挖后，深度0 m温度曲线经过10 d发展相吻合，表明气候条件对初始温度状态不同的浅层土壤影响显著，不同土质土层对气候条件的响应性差异较小；深度1 m温度曲线变化差异是由于坡上部受负温影响，形成一定深度冻结层，弱化了热效应及上部土体的热传导；深度3 m温度曲线变化差异是由于坡中平台受地下水及地热的影响较大；深度2 m温度曲线变化差异既有冻结层因素，又受地下水及地热因素影响。

3.3 冻结深度分析

绘制具有代表性差异的坡上部、坡中平台以及坡脚位置的冻结深度数据曲线如图10所示。

图10 冻结深度变化曲线Fig.10 Curve of freezing depth change

可以看出坡上部冻土上界面最初于3月15日出现融解，而冻土下界面仍缓慢发育，在4月14—23日维持最大冻结深度1.95 m。由于地下热能对冻结土体进行热量补给，4月26日冻土下界面出现一定融化现象，但速率较慢，只有当地表温度为正，才出现消融，在消融过程中，冻结层发生双向融解。4月29日上界面融解深度至0.85 m时，速率减慢，5月2日下界面融化速率大于上界面融化速率，是由于冻结土体下界面的热量补给大于上界面的热量传输，5月14日下界面融解深度至1.52 m时，速率减慢，最终在6月10日冻结土体完全消融。

边坡开挖，坡中平台及坡脚位置土体冻结，其最大冻结深度分别出现在3月24日为0.63 m、3月3日为0.76 m。坡中平台冻结深度受气候变化影响较大，呈波动变化状态，坡脚位置受影响较小，呈逐渐消融状态，是由于其体积含水量不同，低温下土体间存在温度梯度，部分迁移水被冻结，冻结进程受土体含水率影响。边坡开挖初始至4月2日，两处冻结土体完全消融，冻结消融期为31 d。

3.4 不同工况下水热分析

3.4.1 水分场差异性分析

绘制不同工况下监测点的体积含水量数据曲线如图11所示。

图11 体积含水量变化曲线Fig.11 Volumetric moisture content change curve

从图11可以看出，在相同气候边界条件的水热发育下，不同工况的体积含水量变化状况有差异。起初工况2、4体积含水量较高是由于黑土渗透系数较低，饱和体积含水量较高，在同一气候条件下，土体的导水性较差而保水性较好，监测点处于饱和状态；工况1、3由于黄土渗透系数较高，饱和体积含水量较低，土中水向下渗流导致地下水位降低，部分监测点处于非饱和状态，体积含水量较低。

对比工况2、4可以看出，开挖前监测点体积含水量相差不大，开挖后工况2体积含水量下降速率较快，随着水出露，两者最终保持同一体积含水量。是由于工况2上部黑土的渗透性较差，体积含水量较高，边坡开挖，高体积含水量降低了冻结锋面的推进速度，下部未冻结区域较大，未冻结水含量较多，而下部土体处于饱和状态，使得土体边坡水大量涌出，监测点体积含水量下降，同时由于黑土的渗透性、质地等影响，下部黑土的体积含水量下降幅度较工况1、3并不显著，仅为7.99%～27.14%；工况1、3的体积含水量变化曲线基本重合，是由于下部黄土的渗透系数较高、导水性较好，其体积含水量的下降趋势主要受地下水影响，对于上部土体的水分传导依赖性不大，边坡开挖后由于其高渗透性，体积含水量下降显著。

3.4.2 温度场差异性分析

绘制不同工况下坡上部及坡中平台的同一深度2 m的温度数据曲线如图12所示。

从图12可以看出，工况3整体温度较高，且随温度变化反应迅速。开挖前，工况3坡上部黄土体积含水量较低，冻结迅速，但温度下降幅度与其他工况相差不大；坡中平台处于地下水位以下，属饱和状态，受地热及地下水影响，温度变化平稳。边坡开挖，由于黄土渗透系数较高，体积含水量大幅度减少，低含水率对温度响应迅速，气候温度回暖，土体温度上升快，气候温度降低，土体温度衰减迅速。工况2整体温度相对较低，是由于黑土渗透系数较低，体积含水量较高，高体积含水量对温度变化的延迟作用明显，温度上升及下降趋势慢。从整体温度变化趋势可以看出，不同工况下温度变化差异性受体积含水量影响显著。