西伯利亚布鲁斯地下水溢流积冰的边坡水热过程模拟

戴长雷，于 淼，张晓红，Nadezhda Pavlova，尉意茹

（1. 黑龙江大学 寒区地下水研究所，黑龙江 哈尔滨 150080；2. 黑龙江大学 水利电力学院，黑龙江 哈尔滨 150080；3. 黑龙江大学 中俄寒区水文和水利工程联合实验室，黑龙江 哈尔滨 150080；4. 俄罗斯东北联邦大学 自然科学学院，雅库茨克 677000；5. 俄罗斯科学院 西伯利亚分院麦尔尼科夫冻土研究所，雅库茨克 677000）

1 研究背景

地下水溢流积冰也被称为冰湖、冰丘、涎流冰等，是指地下水从地表（或河冰裂缝）溢流出后，受低温影响冻结而成的积冰体，是寒区特有的一种水文地质现象［1］。受温度、降水、局部地形和区域地质等因素影响，溢流积冰呈周期性生消，但并非每年一次，每年的发育程度也不尽相同［2］。同样它也是一种典型的水文地质灾害，在部分地区，溢流积冰漫延开来，会侵害道路或建筑物，以及对基础设施的性能产生负面影响［3-5］。这种现象多发生在北半球典型寒区，如俄罗斯西伯利亚、加拿大、美国阿拉斯加和中国北方等。

Tolstikhin［6］是早期研究地下水溢流积冰的俄罗斯科学家之一，通过多年监测指出，多年冻土区的地下水直接由河床底部沿着河冰裂缝溢出，或在河岸边坡溢出。俄罗斯科学院西伯利亚分院麦尔尼科夫冻土研究所（中国常称之为“西伯利亚冻土所”）Шепелёв［7］通过分析萨哈共和国（雅库特）多处溢流积冰的形成特征及作用特点，认为地下水溢流积冰现象主要受气候因素、冻土因素和水源补给通道等因素的影响。美国阿拉斯加大学Yoshikawa 等［8］、Kane［9］和美国国家航空航天局Hall等［10］通过对积冰区域多年观测以及遥感监测，发现积冰体的发育过程受冬季气温、冻结前降水量以及冻结后的降雪量影响较大。加拿大地质调查局Vellette等［11］通过对加拿大地盾积冰区的钻探和地面温度测量，明确了溢流积冰导水通道的存在及其规模。加拿大麦克马斯特大学Ming-ko Woo［12］对溢流积冰的发生，控制和预防的现有资料进行了总结。戴长雷等［1］、刘月［13］采用实地踏勘和水文地质物探等方法对地下水溢流积冰形成特征进行研究，指出流域下垫面中的浅薄含水层和冻土层是溢流积冰现象形成的关键因素。张浩［14］通过室内外实验结合数值模拟，分析涎流冰的发育机理，并提出相应的防治措施。

国内外相关研究较一致认为，水源、低温和过水通道是地下水溢流积冰发育的主要影响因素。但溢流积冰发生地多处于高寒且具有一定危险性的缺资料地区，相关监测数据与资料获取难度较大。而通过地下水溢流积冰发育过程中研究区温度场与渗流场的耦合模型，可以较为直观地反映出各影响因素在积冰体发育过程中所起到的作用［15］。因此，在实地勘测的基础上，采用Geo-studio软件结合气象数据对不同时期边坡的温度场、渗流场的分布特征进行耦合模拟并分析其对地下水溢流积冰发育的影响，以期为中国地下水溢流积冰相关研究提供参考。

2 材料与数据预处理

2.1 研究区域概况 布鲁斯（61°20'N，129°04'E）位于西伯利亚东部，隶属于俄罗斯萨哈共和国汉加拉斯基区，是西伯利亚地下水溢流积冰代表性区域之一（图1）。该地属于温带大陆性气候，多年平均气温为-9.3 ℃，多年冻土层厚度为200～420 m，活动层厚度为2.0～5.8 m，多年平均降水量约为300 mm［16］。溢流积冰发生在布鲁斯的细长山谷内，山谷呈东西走向，长约1.5 km，两侧边坡植被覆盖密集，溢流积冰出水点位于南侧边坡底部，海拔高度约为150 m。南坡垂高25 m，坡角75°，边坡上层岩性以中、细砂为主；中层岩性主要为砂质粉土，稳定性较好；底部岩层主要为寒武纪碳酸盐岩［17］。距山谷约8.2 km处有一湖泊，在暖季对边坡进行地下水补给。为研究边坡内部地下水渗流规律，在研究区内钻探一口探测井，该探测井距溢流泉（出水口）1.75 km，井深27.4 m，水文地质剖面图如图2所示。

图1 研究区区位及出水点位置

图2 西伯利亚布鲁斯南侧边坡水文地质剖面

2.2 数据来源 选取距研究区较近的伊斯特（Исит）、雅库茨克（Якутск）和捷戈尤里切（Тегюльте）水文气象站测得的气温数据作为研究区气温背景参考。气温数据来自俄罗斯水文气象信息研究所—世界数据中心（http：//meteo.ru），主要为2017—2019年俄罗斯地面日均气温与降水量数据集。碳酸盐岩的强风化带厚度由俄罗斯制造的基于电阻率多电极16 通道断层扫描仪（Многоэлектродная 16-ти канальная электротомографическая станция Скала）进行探测，该仪器可根据电阻率反馈结果，识别出强风化带所处位置和规模。模拟所用水位高度由距出水点（溢流泉）1.5km处的探测井测得。研究区积雪数据采用量雪尺进行人工量测。

2.3 数据预处理 研究区与各水文气象站间海拔高度差异不大但相距较远，因此水文气象站的数据难以直接体现研究区的气温条件。在进行分析前，需要依据各水文气象站间的气温数据对研究区数据进行修正。选择距离反比权重法（IDW）对数据进行距离修正，采用插值点和样点的间距为权重后，再加权平均［18］，公式如下：

式中：X为研究区的气温；Di为插值水文气象站点与研究区的水平距离；Xi为插值水文气象站点的气温；n为插值的样点数目；p为计算距离权重所用幂指数，本研究中选择p=2。

研究区积雪深度数据由人工测得，会出现日数据缺测的情况，采用趋势预测法对其进行补差［19］。由于近出水点是模拟的关键部分，因此得到探测井到出水点模拟结果后，截取距出水点60 m范围内的模拟结果。

3 研究区域水热数值模拟

3.1 渗流及热传导理论 根据渗流连续性方程，有：

式中：kx、ky分别为x、y 方向的渗透系数；Q 为源汇项；θ 为体积含水率；t 为渗流时间；H 为总水头；H=μw/γw，μw为基质吸力；γw为水的重度；y为位置水头。

基本热传导方程为：

式中：T为温度；ktx、kty分别为x、y方向的导热系数；Qt为热源汇项；λ为储热系数，考虑相变，引入未冻结含水率wu，存在关系：

式中：c为土的体积比热容；LW为水的潜热；wu为未冻结含水率。

利用无人机机动性好、操作方便、成本低廉、无人员伤亡等特点，在暴雨或者山洪等极端气候条件下，派出无人机对矿山尾矿库、堆渣场、边坡等进行监测，在发生次生地质灾害的前期及时进行观测，发出预警信号。

3.2 水热数值模拟软件选择 随着数值模拟技术的发展，加拿大岩土公司GEO-SLOPE 开发的Geostudio 2012软件，具有土体渗流稳定分析的SEEP/W模块和具有不同土体岩性传热分析的TEMP/W模块，SEEP/W模块可对多种不同类型的土体模型进行模拟，包括二维渗流控制方程、初始条件和边界条件，其功能具有综合性，既可以分析极为复杂的渗流问题，又具有较高的精度和实用性。TEMP/W模块可以导入测量好的气象数据对瞬态能量流动进行模拟分析。因此采用两种模块间的组合对研究区主要边坡进行水热过程模拟。

3.3 模型的构建 根据布鲁斯南侧边坡为原型建立数值模型，原型边坡的岩土介质从高到低分别为细砂、中砂、砂质粉土、寒武纪碳酸盐岩（强风化带）和寒武纪碳酸盐岩，边坡坡高25 m，基岩取5 m深的碳酸盐岩，坡度75°，地下水在坡角出溢出。基本按照原型建立1∶1的数值模型，模型呈梯形结构分布，上底、下底和高的尺寸分别为51.96 m、60 m和30 m，网格尺寸为0.2 m，如图3所示。

图3 研究区边坡断面网格剖分

首先建立2017—2018年均地面温度条件下的稳态温度场模型以及考虑地下水的边坡稳态非饱和渗流模型。然后建立2017—2018年年循环地面温度条件下的瞬态温度-非饱和渗流耦合模型，计算到1a后某一时刻模型的状态作为初始条件。在此基础上，建立2018—2019年年循环气温及地下水作用下的瞬态水热耦合模型。模型计算的基本参数见表1，模型计算流程及边界条件见表2。

表1 各岩层土层土样的物性参数

表2 模型建立过程

4 结果与讨论

4.1 温度场分布特征 研究区温度场变化主要受气温调控，2018—2019年极低温可达到-43 ℃，自2018-10-14起研究区气温稳定达到负温，至2019-4-23结束，长达192 d的负温环境，也为溢流积冰的发育提供了环境［20］。且研究区位于西伯利亚中部，属于典型寒区，其温度场变化不仅受气温调控，还与积雪深度有关。在溢流积冰发育过程中，积雪深度在2019-2-28日达到42 cm，其保温隔热作用有效地减缓过水通道内地下水的冻结过程，或使其不被冻结。根据地下水溢流积冰所处的状态可将其发育周期分为4个不同的阶段。主要包括加速冻结期、衰弱冻结期、稳定结冰期和融化期，如图4所示。

图4 2018—2019年气温和积雪深度变化曲线

在地下水溢流积冰的发育过程中，受气温变化影响，坡体表层温度变化较为剧烈，尤其是坡肩处，因为坡肩受到坡顶和坡面温度的双重影响（图5）。在加速冻结期（2018-10-11—2018-12-1），受低温环境影响，浅层地温温度降低，并逐渐进入冻结状态。在衰弱冻结期（2018-12-1—2019-2-28），坡体内温度明显降低，坡顶和破面的冻结层厚度逐渐加深。在此期间，碳酸盐岩的强风化带温度在0 ℃附近徘徊，出水点最低温度可达-41.6 ℃。在稳定结冰期（2019-2-28—2019-4-28），受土体内部温度场控制，坡体内部温度进一步降低，一般不高于2 ℃，但由于此时气温在逐步回暖，表层温度在渐渐升高，出水点附近温度已高于-16 ℃。2019-4-28 以后为融化期，坡体表层温度进一步升高，在此期间，出水点温度已达到正温。

图5 2017—2018地下水溢流积冰发育过程中边坡温度场

采用土体初始冻结温度等温线代替冻结锋面，通过温度场数据得到冻土层厚度，进而分析坡体冻融过程［21］。由于坡体上层温度变化较为剧烈，因此对6 m内不同深度的温度分布特征进行分析，如图6所示。在冻结初期，表层土体温度迅速降低，2018-11-11已经形成冻结层，厚度约为0.1 m。至2019-5-12，地面温度达到0 ℃以上，土体开始解冻，冻土上限深度开始降低。2019-6-8冻结深度达到最大值，冻土下限深度为4.97 m，冻结速率为2.36 cm/d，此后冻土下限开始融化，即冻土下限解冻时间晚于冻土上限27 d。

图6 不同深度温度变化曲线

图6 2018—2019年地下水溢流积冰发育过程中边坡渗流场

4.2 渗流场分布特征 通过模拟计算，可以得出坡体渗水部位主要发生在碳酸盐岩的强风化带，其中坡脚处出水点的渗流速度最大，坡体内部其他部位水体流动速率较小，如图6 所示。在加速冻结期，受到同年大气降水的入渗补给，以及远处湖泊的渗流补给，坡体内水位较高。在水头压力的作用下，该岩层内水流速度较快，出水点最大水流速率可达1.23 m/d，该岩层内平均孔隙水压力最高可达44.03 kPa。此阶段降雪量较小，且平均气温低于0 ℃，地表已经开始冻结，使得融雪水入渗量较小，因此对该岩层内水流速度影响不大。在衰弱冻结期，随着冻深的逐渐增加，远处湖泊的渗流补给能力减弱，地下水水位开始下降，导致渗流速度逐渐减缓，出水点最大渗流速度和强风化带岩层内平均孔隙水压力分别在0.78 ～1.23 m/d 和4.94 ～28.82 kPa 范围内变化。在此期间，积雪量开始增加。在稳定结冰期，42 cm的积雪完全融化，但由于季节性冻结层未完全融化，融雪水仍然难以补给到破碎岩层。且地下水水量进一步消耗，水位继续下降，破碎岩层内平均孔隙水压力继续降低，最低点仅为3.98 kPa，出水点最大渗流速度维持在0.58 ～0.78 m/d。到了融化期，此时气温已经上升至0 ℃以上，远处的湖泊重新进行开始地下水补给，破碎岩层内平均孔隙水压力增加，已不低于6.31 kPa，渗流速率开始加快，一般高于0.78 m/d。

4.3 边坡水热变化对溢流积冰发育的影响 根据模拟结果可知，2018—2019年研究区典型边坡活动层厚度约为4.5 m，该层在全年呈现季节性冻融的状态，在其下部衔接多年冻土层，接近多年冻土层顶部为居间不冻层，该岩层主要由寒武纪碳酸盐岩的强风化带（渗透系数较砂质粉土和碳酸盐岩明显升高）组成，终年不冻结，是溢流泉（出水点）稳定而可靠的给水通道。

溢流积冰的发育是一个复杂的过程，在加速冻结期，坡体内地下水受同年大气降水的入渗补给和远处湖泊渗流补给，地下水水位较高，出水点最大渗流速率可达1.23 m/d，水温一般为正温，此时谷底已经积累了一定厚度的积雪，地下水渗入积雪并冻结形成较大规模的底部基础结冰层。该冰层通常厚度不超过10 cm，一般积雪的质量不超过该层积冰体质量的10%［22］。并且由于雪的稀释作用，该冰层的矿化度通常较低。

在衰弱冻结期，基础结冰层发育完全，由于边坡表层含水层开始冻结，随着冻深的逐渐增加，远处湖泊的渗流补给能力减弱，地下水水位开始下降，渗流速度逐渐减缓，出水点最大渗流速度在0.78 ～1.23 m/d范围内变化，水温一般在0 ～1 ℃。最初的积雪已被冻结在基础冻结层，虽然在此期间仍伴随着降雪，但积冰体的发育速度仍逐渐衰弱。且坡体在以0.13 m/d的速度冻结，出水口也有冻结的趋势。但在水头压力的作用下，始终无法冻结（或冻结后很快被再次冲破），此时结冰依靠不间歇（或偶发性）的泉和降雪提供水源，地下水被迫溢流到基础结冰层表面［23］。

到了稳定积冰期，坡体内的地下水位进一步降低，出水点的渗流速度在0.58 ～0.78 m/d范围内，水温在0℃左右徘徊。由于研究区地下水含有一定矿物质，因此结冰点一般低于0 ℃。积雪渐渐融化，融化的雪水一部分渗入土壤，但由于冻结层未完全融化，融雪水难以补给到破碎岩层，另一部分漫流到积冰体表面［24］。在此期间，积冰体开始大范围地横向发展，而不是仅在靠近泉的出口范围内冻结。这是由于随着结冰的生长，冰体内形成了相应的输水渠道以帮助分配溢流出地表的地下水［25］。

融化期的持续时间较长，一般从该年的5月初持续到9月中下旬。随着气温回暖，并在辐射作用下冻结层逐渐融化［26］。坡体内地下水补给量增加，水位抬升，出水量增加。积冰也开始渐渐融化溢流泉中涌出的水、积冰融化的水、融雪水混合着大气降水从山谷中排放到位于勒拿河一级阶地的河迹湖。

5 结论

在西伯利亚东部布鲁斯地区实地勘测的基础上，采用Geo-studio 结合气象数据对2018—2019 年不同时期边坡典型剖面的温度场、渗流场的分布特征进行耦合模拟。但由于研究区条件限制，难以获取完整的实测数据。模拟结果与已观测到的渗流和地温实际数值变化吻合，反映出地下水在多年冻土层内的水热变化，进而分析其对地下水溢流积冰发育的影响。

（1）研究区2018-10-14 至2019-4-23，长达192 d 稳定的负温环境，是溢流积冰发育的先决条件。在溢流积冰发育期间，坡体内塑造的温度环境维持地下水接近冻结-未冻结的临界状态，有效地提高了其溢流到地表（或积冰体表面）后的冻结效率。坡体外部最低温度可达-43℃的超低温度，也为其冻结提供了足够的冷量。

（2）充沛的降水、充足的地下水补给是地下水溢流积冰形成的物质条件，形成较高的地下水位为地下水渗流提供驱动力。自加速冻结期至稳定结冰期，坡体内地下水水位持续下降，水头压力降低，碳酸盐岩强风化带平均孔隙水压力从44.03 kPa降低至3.98 kPa。随后在融化期（至2019-5-11），水位抬升，水头压力增加，平均孔隙水压力已不低于6.31 kPa。

（3）模拟计算结果表明，碳酸盐岩的强风化带是为坡体内为溢流积冰供水的主要过水通道，其地下水渗流速率较其他岩层水流速率明显加快。在溢流积冰发育期间，破碎岩层内地下水渗流速度最高可达1.23 m/d，而其他岩层内一般不超过0.05 m/d。