多年冻土区风积沙与块石混合层对降雨的水热响应分析

韩风雷，喻文兵，陈 琳，胡 达

(1.重庆交通大学 省部共建山区桥梁及隧道工程国家重点实验室，重庆 400074；2.重庆交通大学 土木工程学院，重庆 400074；3.中国科学院西北生态环境资源研究院 冻土工程国家重点实验室，甘肃 兰州 730000

块石路基作为解决冻土融沉的关键技术，在青藏铁路工程得到广泛运用[1-2]。但风积沙填堵块石层后，其换热方式由强迫对流变为自然对流，从而影响青藏铁路块石路基降温性能，对冻土的长期热稳定性不利。针对风积沙环境下青藏铁路冻土块石路基的研究，主要从理论计算、室内模型试验和现场监测开展了传热机理和降温效果的分析。赖远明等[3]通过室内试验在试验路基边坡铺设防水土工膜并设置填充细砂的厚碎石层，模拟大气环境下块石路基风积沙填堵的情况，结果表明封闭块碎石在一定厚度条件下仍具有良好的降温作用，表现为可变等效导热系数的特性。吴青柏等[4]通过对青藏铁路现场工程开放和封闭边界块石路基下部冻土温度的监测，发现块石对路堤下部土体的降温作用开放边界远比封闭条件效果好，但块石路基在被风积沙或积雪填堵后，降温性能大大减弱。陈琳等[5]考虑气候升温条件下，通过数值模拟研究了青藏铁路风积沙填充块碎石层和沙层覆盖块石路基后降温效果变化特征。以上文献主要针对风积沙环境下块石路基长期热稳定性，而对降雨作用引起的块石路基下部冻土水热变化的研究较为匮乏[6-7]。

本文利用青藏高原降雨现场监测数据[8]，结合青藏铁路风积沙填堵块石路基调查，建立水热耦合模型，通过青藏高原红梁河风积沙填堵块石层探坑的现场试验进行数值计算分析，研究降雨条件下青藏铁路风积沙与块石混合层的水热响应特征以及下部冻土的水热状态，为多年冻土块石路基的长期稳定性评价、病害防治和块石路基设计提供理论基础和科学依据。

图1 青藏铁路红梁河风积沙填堵块石路基现场调查

1 青藏铁路块石路基风积沙填堵特征

青藏铁路块石路基所用块石粒径建议值为20～30 cm。采用不规则堆放，室内测试块石层的孔隙率为17%～30%。风积沙迁移填堵块石层空隙，从而改变其多孔介质结构，为研究冻土块石路基风积沙填堵堆积情况及表面分布特征，选取青藏高原红梁河地区青藏铁路断面进行现场调查分析。

风积沙填堵块石现场调查见图1，由图1可知青藏铁路块石路基风积沙填堵堆积状况，除顶部3.5 m范围无明显风积沙覆盖，其他块石层表面均存在沙层覆盖。在块石路基护坡不同位置风积沙填堵厚度差异较大。风积沙堆积范围从护栏至坡脚长度达6.7 m，基底块石层基本被完全覆盖；8.7 m范围内护坡碎石层风积沙填堵厚度大于50 cm，3.5 m范围内碎石护坡层大于40 cm深度处仍有风积沙存在。同时风积沙只能填充一定厚度的块石层，在一定深度处块石层空隙未被完全填充。说明风积沙已阻断了块石层与外界空气的流动，其边界条件发生改变。天然地表风积沙堆积厚度见图2。每隔40 cm测量风积沙覆盖厚度，距离坡脚越近覆盖厚度越大，最大、最小覆盖厚度分别为41.2、16.5 cm，基本呈指数递增，主要因为该段铁路西侧约1.5 km的高大沙丘，在西风作用下，成为此处沙害的主要来源，并在路基坡脚形成大量堆积，填堵覆盖块石护坡[9]。因此，在风沙灾害区，风积沙迁移填堵块石层结构比较严重，改变了块石路基边界特征和换热方式，影响了块石层的降温性能[10-11]。

根据现场及室内风积沙测试结果，风积沙颗粒以0.075～0.5 mm为主，其中0.075～0.25 mm占比57.3%，比表面积小，亲水性弱，由于青藏高原降水量小蒸发量大，其天然含水率较低。风积沙渗透性较好，孔隙比在0.482～0.641之间，其渗透系数为0.010～0.023 m/s[12-13]。夏季降雨通过风积沙较易渗入块石路基下部，改变了下部冻土水热状况。采用水热耦合模型，对降雨环境下冻土水热分布和变化特征进行计算分析。

图2 路基坡脚外风积沙堆积厚度

2 风积沙填堵块石层水热耦合控制方程

风积沙与块石混合层包括风积沙、块石、液态水、固态冰和空气组成的混合材料。假设其分布均匀，各向同性。对于多孔介质中水分运动问题包括液态水的运动和冰水相变，忽略气态水的迁移。根据质量守恒定律，对于变饱和多孔介质中水分运动问题可采用修正的Richards方程来描述[14-15]为

(1)

式(1)中等号左边表示多孔介质水分含量的时间变化；右边第一项分别表示压力水头、重力和温度梯度对水分入渗的影响，第二项S表示地表植被蒸发和根系吸收的源(汇)项。对于青藏高原地表裸露，S的影响可忽略不计。

为了描述非饱和土冻融循环中水力特性采用Van-Genuchten模型和Mualem模型[16]，引入独立参数m、n。

(2)

(3)

当取m=1-1/n时，式(3)为

(4)

温度梯度下非饱和土的渗透系数可表示土水势作用下液态水导水系数的温度效应[14]为

(5)

γ=75.6-0.142 5T-2.38×10-4T2

(6)

式中：GwT为无量纲的增益系数；γ为土水表面张力；γ0=71.89 g/s2为25 ℃条件下土壤吸附水分的表面张力。

非饱和土在冻结过程孔隙水凝结成冰会堵塞水流通道，使冻土的导水能力大大下降，从而影响整个水分场的分布状况。因此，为反映冰颗粒对水分流动的阻碍作用，引入阻抗系数的概念，因此冻结区的渗透系数表达式为[17]

Kfh=Kh10-ΩQ

(7)

式中：Ω为阻抗系数；Q为含冰量增大对水分流动通道的阻滞作用，可取Q=θi/(θ-θr)。根据文献[15]，本文风积沙块石层水热计算时，阻滞系数取为10，以此表示冰颗粒的阻水效果。

非饱和-饱和土壤的冻融作用的热传递过程，只考虑固体颗粒骨架间的热传递，液态水迁移和固-液相变，忽略液-气相变过程中的热量变化和气态水扩散，计入液态水迁移带来的热量,则温度场控制方程表示为[18-19]

(8)

式中：Cp为多孔介质的等效体积热容，J/m3·℃；λp为多孔介质的等效导热系数，W/m·℃；Lf为液态水冻结或冰融化潜热，约为3.34×105J/kg；qw为向下运动的液态水的通量，m/s；Cw为液态水的体积热容，J/m3·℃。

式(8)左边第一项表示单位时间混合基质材料的能量变化，第二项表示冰水相变潜热；右边第一项表示二元基质材料的热传导，第二项表示液态水流动的对流换热。

由于空气和水汽对冻土体积热容影响较小，可忽略不计，则等效体积热容可由不同组分的体积加权表示为

Cp=Csθs+Cwθw+Ciθi

(9)

式中：Cs、Cw、Ci分别为固体颗粒、液态水、冰的体积热容；θs、θw分别为固体颗粒、液态水的体积分数，θs+θw+θi=1。

同理，多孔介质的等效导热系数可表示为

λp=λsθsλwθwλiθi

(10)

式中：λs、λw、λi分别为固体颗粒、液态水、冰的导热系数。

土体冻结之后其含有的液态水并不能全部转变成固态的冰，由于颗粒表面能的影响其中一定数量水会始终保持为液态，冻土中的未冻水含量与温度之间存在着动态平衡关系，即随温度降低，未冻水含量减少，反之亦然。采用徐斅祖等[20]的经验表达式确定冻结过程中风积沙中最大未冻水含量为

θu=a|T|-b

(11)

式中：a、b均为与土的性质有关的试验常数。

3 计算模型和水热参数

为研究青藏铁路块石路基在风积沙填堵后，混合层及其下部冻土在降雨环境下的水热变化，结合青藏高原红梁河野外试验现场温度监测数据，计算模型简化见图3。计算区域为天然地表向下30 m，块石层边界向外延伸5～6倍块石层宽度。计算区域土层0～2.5 m为砂砾土，2.5～30 m为强风化泥岩。计算模型各层介质热物理参数见表1。

根据附面层理论[22]和现场观测[23]，温度边界根据红梁河不同下垫面观测数据拟合而成，天然地表AB和CD边的温度边界条件为

(12)

风积沙填堵块石表面BC边的温度边界条件为

(13)

30 m深度GH边的热流密度为0.06 W/m2[21]，边界AEG和DFH视为绝热。考虑水分入渗时水流通量携带的热流量按下式计算为

Q=qwρwCwΔT

(14)

式中：ΔT为入渗边界的温差。模拟计算中降雨的温度近似可取10 m高度处的气温或露头温度，现场实测7月份的雨水温度平均值约为5 ℃[24]。

青藏高原中部地区降雨时间一般从每年4月至10月，其中7月份降雨量比较显著，年累计降雨量为300～500 mm之间，陆面蒸发过程集中在每年5—10月份，以7月份蒸发最为强烈[25]。青藏高原中部地区降雨量不大，多以短期、高频次降雨为主。因此地表按照无积水、无径流考虑。各介质层水分运动参数见表2，保温层按隔水材料计算。土层表面蒸发强度按均一值考虑，蒸发量以7月平均日蒸发量计算，取为3.07 mm/d[8]。

图3 计算模型(单位：m)

表1 各介质热物理参数[21]

表2 各介质水分运动参数

4 降雨对风积沙块石层水热状况的影响

风积沙与块石混合层计算模型按照中心对称进行分析。根据计算区域的边界条件和水热参数，计算无降雨条件下天然土层的初始温度场和水分场，并将获得的计算值作为模型计算的初始值。模型计算初始时间从冬季最低温度开始，降雨采用青藏高原中部地区气象观测站2013年每0.5 h降雨量观测数据，年降雨量为315 mm，见图4。模拟计算自2013年1月起预测分析未来20 a内降雨环境下风积沙与块石混合层下部冻土及天然土层水热分布特征和变化规律。

图4 计算模型初始值和降雨量

图5 风积沙填堵块石层温度场分布 (单位：℃)

4.1 风积沙填堵块石层温度场特征

风积沙填堵块石层后2023年即第10年4月和10月的等温线分布见图5，红色区域表示风积沙块石层和保温板。由图5(a)可知，风积沙填堵块石层后温度较天然土层高，融化时间提前，冻土上限为-0.426 m，风积沙块石层下部冻土温度抬升，低温区向两侧退化。由图5(b)可知，天然地表进入冻结状态，风积沙块石层下部形成扇形高温区，最大融化深度为-3.141 m，天然地表的冻土上限为3.032 m，两者差距不大。风积沙填堵块石层后增强了夏季对环境温度的响应，在底部形成高温区，改变下部冻土热状况。

冻土上限及温度变化见图6、表3。由图6和表3可知，天然土层和风积沙块石层在2013年至2023年的20年内，两者冻土上限分别增加0.128、0.120 m，两者变化速率差距较小，相对而言天然土层受降雨影响大。冻土温度取天然地表以下15 m深度处，受降雨影响冻土温度呈周期函数波动并逐渐降低，天然土层和风积沙块石层下部冻土温度差异不大，但其温度最大值和最小值存在降低趋势。天然土层和风积沙块石层，20 a内两者冻土温度分别降低0.224、0.305 ℃。与冻土上限相比，降雨对土层冻土温度的影响大。风积沙的填堵改变块石层热传导系数，导致块石路基降温性能降低，但降雨作用会消耗一部分外界传递的热量，减少向下传递的热量。因此，降雨作用下，随时间冻土上限抬升，冻土温度下降。说明降雨增加了液态水对流传热，土层热传导通量减小，有利于缓解对冻土的热扰动。

不同深度处土层温度时程变化曲线见图7，土层温度呈周期函数变化且温度幅值逐渐减小，随深度增加土层温度对外界大气环境的响应逐渐减弱。表4和图8分别为不同深度处土层温度幅值和变化量。与天然土层相比，风积沙填堵块石层后冻土温度高，受降雨的影响温度降低速率大，土层温度最小值相对于最大值对降雨环境更为敏感。同时在冻土活动层深度范围内，夏季降雨能够降低土层温度。

图6 冻土上限和冻土温度时程曲线

表3 不同土层冻土上限和冻土温度变化

图7 不同深度处温度时程曲线

表4 不同深度土层温度幅值变化 ℃

图8 不同深度处第1年与第20年温度变化量

4.2 风积沙填堵块石层后水分场分布

风积沙填堵块石层后第10年4月和10月水分场分布见图9。由图9(a)、图9(b)可知，天然土层处于负温，液态水体积含量表示冻土中未冻水含量。对于风积沙块石层下部土层，深度-1.5 m和-2.5 m附近体积含冰量分别为12%和16%，其体积含冰量相对较高存在集中分布。对于天然土层，最大和最小体积含冰量分别为12%和2%，主要集中在冻土上限附近，与风积沙填堵块石层相比小4%。最大融化深度时液态水体积含量见图9(c)，由图9(c)可知，风积沙块石层下部土层液态水含量为21%，相比于4月15日增大了5%，主要集中在深度-3.0～-2.0 m冻土上限附近，冻土上限下部多年冻土区液态水体积含量基本不变。对于天然土层，活动层区域液态水含量增大，最大值为17%，主要集中在冻土上限附近。因此，风积沙填堵块石层后，受降雨影响其下部土层液态水体积含量和体积含冰量都相对增大，冰水相变释放的潜热和体积膨胀会对下部土层变形和热响应产生不利影响。

图10和表5为风积沙填堵块石层后不同深度水分含量。由图10(a)可知，总等效体积含水量深度-4.0 m处为0.15基本保持不变；深度-2.5 m处第1年与第20年相比，总等效体积含水量增加了0.052；深度-1.5 m处则减小了0.124。由于深度-4.0、-2.5、-1.5 m处土层初始总等效体积含水量分别为0.15、0.17、0.14，在风积沙填堵块石层后，其下部土层总等效体积含水量逐年增加，存在水分累积现象，深度-1.5 m处受降雨和高渗透性风积沙影响，其含水量变化幅度较大，多年冻土区基本不受降雨影响。由图10(b)可知，液态水体积含量，前5年深度-1.5 m和-2.5 m处液态水体积含量呈周期性波动逐年分别降低和增加，速率变化大。对于液态水体积含量，深度-1.5 m处第1年和第20年相比，最大值和最小值分别降低0.125、0.003；深度-2.5 m处其最大值和最小值分别增加0.041和降低0.005；深度-4.0 m处液态水体积含量呈周期函数变化，其最大值和最小值分别降低了0.004、0.006，该深度处于冻结状态液态水体积含量变化不大。由图10(c)可知，体积含冰量与液态水体积含水量变化规律基本一致。土层体积含冰量呈周期性波动，深度-1.5 m处第1年、20年体积含冰量逐渐降低；深度-2.5 m处则逐渐增加；初始阶段集中在-1.5 m处，之后逐渐迁移到-2.5 m附近，主要是夏季降雨入渗随时间迁移到下部冻土层。深度-4.0 m处等效体积含冰量变化不大。

图9 土层水分场分布(单位：m3/m3)

表5 不同深度风积沙块石层水分含量 m3/m3

图10 风积沙块石层不同深度水分变化

图11和表6为天然土层不同深度水分含量变化。由图11(a)可知，总等效体积含水量深度-4.0 m处为0.15基本保持不变；总等效体积含水量深度-2.5 m处第1年与第20年相比增加0.007；深度-1.5 m处则增加0.003。由图11(b)可知，液态水体积含量第1年和第20年相比，深度-1.5、-2.5、-4.0 m处，其随时间呈周期性波动且最大变化幅值分别为0.003、0.006、0.004。因此，天然土层液态水土层不存在水分累积现象。由图11(c)可知，等效体积含冰量与液态水体积含水量变化规律基本一致，等效体积含冰量最大变化幅值为0.008。与风积沙块石层下部土层水分变化相比，天然土层水分含量小，受降雨影响不大，冻土活动层内水分主要集中在冻土上限附近且随时间增加，冻土区不受影响。由于青藏高原蒸发量大，每半小时降雨量小，不存在水分累积现象。

图11 天然土层不同深度水分变化

表6 不同深度天然土层水分含量变化 m3/m3

5 结论

本文在调查青藏铁路风积沙填堵块石层实际情况基础上，通过建立水热耦合模型对风积沙和块石混合层在降雨条件下其下部冻土温度场和水分场变化进行分析，得到以下结论：

(1)风积沙填堵块石层后，下部土层易形成扇形高温区，增大对外界环境温度的响应。降雨作用下，风积沙块石层和天然土层下部冻土上限前5 a抬升速率大，之后平缓增加，与初始冻土上限相比20 a内分别抬升了0.128、0.120 m。

(2)天然土层和风积沙块石层下部冻土温度差异不大，随时间呈周期波动，温度幅值存在降低趋势其20a内分别降低了0.224、0.305 ℃。与冻土上限相比，降雨对冻土温度影响大。降雨增加了液态水对流传热，土层热传导通量减小，有利于缓解对冻土的热扰动。

(3)风积沙填堵块石层后，总等效体积含水量浅表层深度降低，在冻土上限附近即深度-3.0～-2.0 m处增加0.052，夏季液态水体积含量和冬季体积含冰量分别由0.17、0.12增大至0.21、0.16，且在风积沙填堵块石层下部存在水分累积。天然土层下部水分含量变化小，受降雨影响不大，由于青藏高原蒸发量大，不存在水分累积现象。

(4)在年平均气温较高区域，降雨虽然能够抬升冻土上限和降低冻土温度，但风积沙填堵块石层后其下部土层存在水分累积，会引起冷暖季节的冻胀融沉病害，应加强防排水设计措施，防止路基坡脚积水造成侧向入渗，从而保证路基稳定性。