川西北高寒区冻融交替作用后土壤水热运移模拟研究

宋 洁, 李晓宁, 赵 丽, 樊 敏

(1.西南科技大学 环境与资源学院, 四川 绵阳 621010; 2.河北邯郸市瀚林环境评估有限公司, 河北 邯郸 056000)

川西北高寒区处于长江和黄河的发源地，能够在黄河和长江上游的水源涵养与补给、生态平衡中发挥着重要作用[1]。高寒地区土壤长期经历昼夜及季节性冻融作用的影响，土壤内部结构易于发生变化：孔隙度增大，容重减小，破坏团聚体水稳性，并减小其抗剪强度，增大土壤可蚀性，造成土壤退化[2-4]。同时，土壤周期性冻融作用对动植物生存存在影响：土壤热量传输，水分相变和盐分堆积，改变了土壤的水分传输能力，从而影响化学物质的运移[5]，导致植被退化，野生动物的栖息地质量和生物多样性下降。因此，土壤水热运移规律的研究对于综合评价地表、地下水资源，有效地利用土壤水、热资源，合理解决高寒区资源的开发、保护生物多样性和土壤沙化防治与植被恢复等实际问题都具有重要意义。尚伦宇等[6]揭示了土壤水热变化对青藏高原地表能量的影响。于炜[7]分析了科尔沁地区沙坨地和草甸地两种典型地貌条件下的土壤冻融期内水热动态变化过程及联系。随着科学技术发展，许多学者在研究土壤水热运移时将模型模拟方法运用其中，Harlan[8]最早提出了基于非饱和土迁移机制的土壤冻结过程中水热耦合运移模型，郑秀清等[9]采用包括水迁移和热对流迁移的水热耦合数值模拟模型，模拟天然条件下土壤的季节性冻融过程以及其中的水热迁移规律。李瑞平等[10]利用SHAW模型研究了内蒙古河套灌区盐渍化土壤冻融期水热盐的动态变化；王宇祥等[11]利用HYDRUS软件模拟科尔沁沙地沙丘—草甸相间区土壤水分动态变化指出：流动沙丘和草甸地降雨与表层土壤水分呈极显著相关。目前，国内外关于高寒区冻融土壤研究多在融雪入渗水热运移规律[12]、冻融土壤物理结构分析[13]、气候对土壤冻融影响等问题的试验与模拟研究，且国内研究多在藏北、东北地区进行，对于川西北高寒区冻融土壤水热运移对土壤沙化进程的影响研究较少。本文以川西北高寒区经过室内冻融交替作用后沙化和天然草地土壤为研究对象，利用HYDRUS软件模拟土壤剖面水热变化，明确川西北高寒区土壤水热运移特性，探明冻融交替条件作用后沙化和天然草地土壤体积含水率和温度随土层深度变化的规律，为该地区土壤沙化治理与预防，维护生态系统平衡提供科学依据和理论支持。

1 材料与方法

1.1 研究区概况

研究区位于阿坝藏族羌族自治州红原县(31°50′—33°22′N，101°51′—103°23′E)，位于青藏高原东部，地势西北高、东南低，海拔在3 000～4 000 m范围内。研究区内季节性冻土分布较广，气候属大陆性高原寒温带季风气候，四季变化划分不明显，冬长夏短，寒冷气候占据全年大部分时间；但日照长，太阳幅射强烈，早晚温差极大，极端最高气温24.6℃，极端最低气温达-22.8℃，年平均气温2.9℃[4]。阿坝藏族羌族自治州红原县干雨季节分明，雨热同季，降水量较为丰沛，年降水量可达860.8 mm，年均积雪期可达76 d。研究区主要的土壤类型为沼泽土、亚高山草甸土和风沙土。按照土壤质地的国际分类制，研究区沙化草地土壤属于砂质土；天然草地土壤属于砂壤质土[4]。沼泽、草甸、灌丛和森林为该区主要植被类型[4]。

1.2 供试材料

1.2.1 土壤采集及物理性质测定 研究采用的土壤取自阿坝藏族羌族自治州沙化最严重的地区之一红原县瓦切乡。依据《天然草地退化、沙化盐渍化的分级指标》(GB19377—2003)选取沙化草地和植被覆盖良好的天然草地2类。沙化草地的平均盖度13%，物种丰富度低，极多裸露砂粒，干燥、人为扰动大；天然草地土层较厚，平均盖度95%以上，植物物种丰富。植被类型主要有四川蒿草(Kobresiasetchwanensis)、细叶亚(Ajaniatenuifolia)、小柴胡群丛(Bupleurumtenue)等。

土壤采集原状土，每15 cm为一层，共取4层土壤样品，用于装填实验室水热运移模拟土柱。每层土壤分别用环刀重复取样3次，用于测定土壤容重和含水率。采用环刀法测定土壤含水率和容重，土壤物理性质见表1。

表1 供试土壤物理性质

1.3 研究方法

1.3.1 试验装置 本研究采用3个自制土柱土壤积水入渗装置进行入渗试验(图1)，将原状采集土壤各层分别风干后过2 mm筛，沙化草地和天然草地土壤各层含水率和容重按照野外测定数据设置(表1)，采用干堆法，根据公式(1)在每次装入土壤后，压实使其达到规定的土柱高度15 cm，保证试验土柱各层土壤干容重与天然土壤干容重相近。

W=V·γ(1+S)

(1)

式中：W为土壤质量(g)；V为每次装入的土体体积(cm3);γ为天然土体干容重(g/cm3);S为室内土壤含水量(%),室内土壤含水量则根据采样时的实测值,采用称重法,保持土壤重量含水率与实测值一致。沙化草地和天然草地土柱装填高度均为 60 cm,每隔12 cm开设孔洞(12 cm,24 cm,36 cm,48 cm和60 cm),用于插入土壤水分传感器探头(型号ECH2O),水分传感器探头孔洞对侧设置排气孔(φ=0.2 cm)。土柱底部设置排气孔板和排水室(h=0.2 cm,φ=10 cm),出水由烧杯直接截取,利用马氏瓶稳压供水,供水水温控制为室温25℃。

1.3.2 冻融条件设定 根据川西北高寒区实际冻融现状，设置连续反复冻融(模拟昼夜)模式，依据冻结温度设置为-11℃，每天冻结12 h，在空调室内恒温25℃进行融化，消融12 h的冻融条件，连续冻融6 d，第6天融化后进行水热运移试验。

图1 土柱试验装置示意图

1.3.3 试验步骤 经6 d连续反复冻融作用后，对土柱进行恒定水头入渗试验。土柱分别在12 cm，24 cm，36 cm，48 cm和60 cm处插入ECH2O监测探头，数据采集仪自动记录剖面土壤含水率和温度变化实时监测值，读取频率1次/ min。控制马氏瓶内水面高度，确保土柱表面水分入渗恒定水头2 cm。当排水室有水流溢出时，试验结束。由于试验持续时间较短，试验分析中不计蒸发对入渗过程的影响。

2 HYDRUS模型构建

2.1 基本方程

(2)

式中:θ为土壤体积含水量(cm3/cm3);t为时间(min);z为一维垂直向坐标(cm);K(h)为土壤的非饱和含水率(cm/s);h为压力水头(cm);z为土柱纵剖面空间坐标(cm),原点在图层的上边界,向下为正。本试验不考虑作物根系对土壤水的吸收作用,故S为0。

土壤水分特征曲线和非饱和导水率用Van Genuchten方程拟合：

(3)

(4)

(5)

式中:θ为土壤体积含水率；θs为土壤饱和含水率(cm3/cm3)；θr为土壤残余含水率(cm3/cm3)；K为土壤的非饱和导水率(cm/s)；Ks为饱和导水率(cm/s)；Se为无量纲的相对含水量；α，m(m=1-1/n)，n,l(一般取0.5)均为拟合参数。

2.1.2 热运移方程 土体热量运移采用的基本方程(仅考虑液态水运动)为：

C(θ)=Cnβn+Cwβw+Caav

(6)

式中:Cp(θ)为多孔介质比热容〔J/(cm3·℃)〕;Cn为固相比热容〔J/(cm3·℃)〕;Ca为气体比热容〔J/(cm3·℃)〕;Cw为液体比热容〔J/(cm3·℃)〕;λ(θ)为土壤导热率(cm2/s);q为水分通量(g/s)；T为土壤温度(℃)。

2.2 初始条件与边界条件

2.2.1 初始条件 在研究土壤入渗过程中的水热运移规律时，初始条件为各层土壤初始含水率，不考虑蒸发、降雨和地下水对试验的影响情况下，t=0时，土体剖面水流模型初始条件：

h(z,t)=h0(z)t=0, 0≤z≤Z,Z=66 cm

T(z,t)=t0(z)

(7)

式中:h0(z)为土体初始负压水头(cm);t0(z)为土壤初始温度(℃)。

2.2.2 边界条件 水流上边界为恒定2 cm水头边界。下边界为自由排水边界：

h(z,t)=h0(t)z= 0

(8)

式中:h0(t)为上边界边界压力水头值(cm)；k为土壤非饱和导水率(cm/min)；q0(t)为土壤下边界水通量(cm/s)。

T(z,t)=T0(t)z= 0

(9)

式中:T0(t)为进水流温度(℃);T和T0分别为土壤和下边界温度;λ为为土壤导热率〔J/(cm·K)〕；Cw为液体比热容〔J/(cm·K)〕；q0为下边界土壤水通量(cm/s)。

2.2.3 土体剖面信息 模拟土体剖面空间步长设置为1cm，60 cm土柱则剖分为60层，共60个节点，距地表不同距离处设定土壤水分、温度及水势变化监测点，含水量和温度数据获取点为4个，分别在12 cm，24 cm，36 cm和48 cm处。设置初始时间步长0.001 d，可变步长设置为0.5～0.001 d。

2.3 模型参数

2.3.1 土体水力学参数 在HYDRUS-1 D模型的水分运动模型中，输入不同土层深度的粒级分析(砂粒、粉砂粒和黏粒的百分含量)和容重值，采用RETC软件对已有试验数据进行拟合，并利用软件获取Van Genuchten模型中的土壤含水率参数(θr，θs，α，n和l)以及饱和水力传导系数(Ks)等相关参数，土体相关水力参数见表2。

2.3.2 土体热力学参数 热力学参数利用HYDURS软件的内部程序Rosetta Lite.v.1.1获得，土壤固相占总体积的比率为0.57；土壤纵向热扩散率DL(cm2/s)为5.00；土壤横向热扩散率Dr(cm2/s)为1.00；热导率函数系数(b1)为1.13×108；热导率函数系数(b2)为1.53×108；热导率函数系数(b3)为1.16×108；土壤固相热容Cn(J/(g·℃)为6.91×108；土壤有机质热容Co(J/g·℃)为1.04×108；土壤液相热容Cw(J/g·℃)为1.54×108。

2.4 模型验证与评价

2.4.1 模型验证 对实测的土体数据进行拟合，分析土柱体积含水率和温度模拟数据的模拟效果。本文采用均方根误差(RMSE，公式10)及决定系数(R2，公式11)2个指标来评价模型的模拟效果。

(10)

(11)

2.4.2 土壤体积含水率和温度动态变化的验证结果为了验证该模型模拟的效果，利用HYDRUS-1D模拟经过连续冻融作用后沙化草地和天然草地土壤含水率和温度随时间变化过程，结合实测数据，绘制土壤体积含水率和温度模拟值与实测值对比曲线图(图2和图3)。为了观察模拟数据的模拟效果，对实测的土体数据进行拟合，采用均方根误差及决定系数2个评价指标对模拟值和实测值之间的偏差进行评价。

表2 率定后土壤水分特征参数

从图2和图3可见，沙化草地和天然草地各土层体积含水率和温度模拟值与测定值接近，且变化趋势基本一致。天然草地15 cm处土壤体积含水率模拟值略低于测定值，这种差异的存在说明该层土体接近地表，容易受到外界环境因素的影响。60 cm处土壤体积含水率模拟中期数据高于测定数据，说明土体初始含水量越高，土柱在融化时间内未能完全融化，导致实测时体积含水率偏低，随着深度的增加模拟精度有所提高，这可能与接近地表的热通量受到外界的影响较多。

图2 不同土层土壤水分随时间变化的模拟值与测定值对比

2.4.3 模型模拟结果评价 经过参数的率定和验证，沙化草地和天然草地土壤体积含水率和温度在剖面各节点上的相对差值都较小，展现了比较好的拟合效果和拟合精度。模拟值与实测值的均方根误差RMSE在合理范围内，决定系数R2大于0.9。HYDRUS-1D模型能够用于模拟土壤含水率和温度的变化，具有一定的可靠性与稳定性(表3)。

3 模型应用与分析

3.1 HYDRUS-1D模拟预测土壤水分变化过程

选取2010年实测的辐射、气温、湿度等气象数据输入模型，模拟并分析土壤经过冻融循环作用后，土壤体积含水率变化特性。

图3 不同土层土壤温度随时间变化的模拟值与测定值对比

表3 模拟结果的评价

图4中显示冻融循环作用前后各层土壤体积含水率模拟变化情况。冻融前(图4A和图4B)，沙化草地0—15 cm土壤体积含水率最低，45—60 cm土壤体积含水率最高，15—30 cm和30—45 cm土壤体积含水率次之，且差异较小；天然草地各层土壤呈现随着剖面加深，土壤体积含水率逐渐降低的趋势；沙化草地表层土壤体积含水率与天然草地表层土壤体积含水率相差约10%。反复冻融循环后(图4C和图4D)，两种土壤含水率随季节变化均呈现上下波动的趋势，但天然草地各土层含水率范围在0.08%～0.25%，总体均高于沙化草地含水率的变化范围为0.06%～0.16%；天然草地各层含水率有明显差异，0—15 cm土壤含水率最低，15—30 cm土壤含水率最高；沙化草地各层土壤呈现随着剖面加深，且各土层含水率差异较小。

3.2 模拟预测冻融循环条件下土壤温度变化过程

土壤冻融状况受土壤温度影响，在温度梯度作用下，土壤温度会随着太阳辐射变化和季节更替而出现昼夜变化和季节波动。

从图5中可看出冻融循环前后，土壤温度模拟值呈现先波动上升后波动下降的趋势，60～300 d两种不同程度植被覆盖下各层模拟土壤温度变化趋势基本一致，均呈现随着剖面加深，土壤温度逐渐降低的趋势，且沙化草地土壤各层最高温度平均低于天然草地各层土壤温度1.8℃。60 d前和300 d后两种不同程度植被覆盖下各层土壤温度呈现与60～300 d相反的变化趋势，即底层温度最高，表层温度最低。沙化草地和天然草地表层土壤温度最大值均出现在6月份下旬和7月份，气温最高分别为12.4℃和15.6℃，土壤温度最低为1月份下旬。

4 讨 论

4.1 冻融作用下土壤水分运移特性模拟

通过对沙化草地和天然草地冻融前后土壤中水分运移模拟表明，冻融方式和植被覆盖度对水分运移的影响极大。郭志强等[14]研究表明，夏季昼夜温差大，土壤体积含水率在强烈蒸发、强辐射作用下出现显著变化，导致土壤昼夜冻融反复循环的发生；冉洪伍等[12]对藏北高寒草地土壤冻融过程水分变化的研究结果显示，由于积雪覆盖于土壤表层，土壤处于长期冻结状态，各类草地土壤温度结果均表现为无明显变化；另有研究表明草甸植被覆盖度的下降可以导致土壤含水量下降[15]。本研究综合考虑以上因素对土壤水分运移的影响，结合HYDRUS-1D模型分析水分运移规律及其在草地退化及沙化过程中的作用。

HYDRUS模拟预测冻融前后沙化草地和天然草地的土壤水分运移特性中，不同深度土壤的含水量分别在75 d和50 d之后才明显增大，分析原因在于一维模型忽略了土壤中空气压力的变化[16]和水流的滞后作用。沙化草地土壤体积含水量在120 d前后的变化趋势相反，考虑到因前120 d为该区的1—4月份，气温未回升，无水分入渗和降水补给，且蒸发不强烈；120～300 d为该区的5—10月份，气温回升，积雪融化，雨季来临降水量逐渐增加，气温逐渐达到一年中的最高温(15.6℃)，蒸发作用增强，导致表层的土壤体积含水率低于深层的土壤体积含水率。天然草地在100～300 d内土壤体积含水率沿着剖面自上而下呈降低趋势，分析原因为此时间段积雪融化水分入渗、降水量逐渐增加，红原县牧草生长的茂盛，土壤表面植被覆盖率高，有效减缓了土壤水分的蒸发，加之植物根系发挥了涵养水源的作用。

图4 冻融循环后各土层土壤体积含水率模拟

图5 不同冻融循环后各土层土壤温度模拟

反复冻融循环后，由于高寒区季节变化、降水差异、蒸发强度和土壤理化性质的差异，土壤水分的垂直分布和时间变化具有明显的规律。3—5月份为冻土消融前期，气温回升雪水融化入渗，各层土壤体积含水率逐渐增加，表层土壤中变化更为明显。此时，沙化草地和天然草地各层土壤体积含水率沿着土壤剖面向下呈现增加趋势。6—10月份气温升高加速深层冻土消融，夏季降水量增加，沙化草地和天然草地土壤体积含水量均呈不同程度的上升趋势；7月份左右气温达峰值，土壤体积含水率在强烈蒸发作用下出现显著变化，加之研究区海拔高、太阳辐射强，在冻融期间表层土壤和大气之间时刻发生着水分和能量的流动，导致了冻融期间昼夜冻融循环的发生，土壤体积含水率波动性明显。同时土壤质地对于冻土水力传导度的影响更为显著，植被覆盖率较高的天然草地土壤质地黏性大，冻结后各土层水力传导度的变化幅度要显著小于砂性较重的沙化草地土壤。10月份至次年3月份为土壤冻结期，积雪覆盖于土壤表层，隔绝大气与土层的直接接触，土壤处于冻结状态，水分无法入渗，沙化草地和天然草地土壤水分模拟结果均表现为无明显变化。

研究表明，植被覆盖率与土壤有机质含量、容重和孔隙度均有显著相关关系，连续反复冻融过程是土体的密实过程，随着土壤密实度增加，干容重减小，土体的孔隙率变大，在补水条件下，连续反复冻融后，冰变成水体积减小承载力下降，土体内部发生沉降。连续反复冻融下对天然草地表层土壤容重和孔隙较其他土层影响大，0—10 cm土层含水率明显低于其他土层，主要是因为冻融循环使表层土壤的毛管孔隙度下降，持水力降低。植物根系及其生物量的大小在蓄水保肥、防风固沙方面至关重要，根系生物量及其分布特征与容重、有机质的物化性质变化相关性较高，两者的变化均易于引起土壤持水性的变化，草甸植被覆盖度的下降可以导致土壤含水量下降。

4.2冻融作用下土壤温度变化特性模拟

沙化草地和天然草地土壤表层和深层温度在时间序列上呈现截然相反的变化趋势，这可能与接近地表的热通量受到外界的影响较多，另外说明土壤热参数的取值有一定的误差[17]。沙化草地和天然草地土壤温度在季节上呈现一定的周期性，3月份以前和9月份以后，土壤温度呈现上层低于下层的趋势，3—9月份土壤温度从表层到底层逐渐降低，距表层越近温度变化幅度越大[4]，加之热量损耗导致土体内部随深度增加，波动幅度呈现减小的趋势。在整个模拟周期内，土壤各层温度呈现不断波动的现象，由于冻融土壤受降水、蒸发、太阳辐射影响较大，温度变化剧烈[17]。土壤冻融循环过程中，温度随之上下波动，土壤中的颗粒、水分和气体的占比及分布随之发生变化，当土壤内的水发生固液相变化时，水或冰会对周围的土体产生不同的挤压力，这一作用力会破坏土壤颗粒之间的胶结作用，引起土壤的大孔隙、中孔隙、小孔隙和毛细孔隙位置及形态发生变化，从而也改变了土壤孔隙形态结构。太阳照射能够引起土壤的温度升高，在试验过程中土壤物理性质的变化没有引起土壤温度的明显变化。

5 结 论

基于土壤水分运动的动力学方程和土壤热流基本方程，建立水热运移数学模型，利用HYDRUS-1D软件对模型进行求解，结合实测数据对模型模拟进行评价，主要结论如下：(1) 入渗水头为2 cm时，经冻融交替作用后60cm土柱剖面各土层土壤温度和水分的测定值和模拟值吻合较好，率定后的土壤水热运移数学模型能够较好地反映出土壤水热的变化及分布状态。(2) 利用率定后的饱和含水率(θs)、残留含水率(θr)等设定HYDRUS-1 D软件中的模型运行参数，模拟预测长时间(365 d)的土壤水热运移动态。此结果有助于了解川西北高寒区冻融交替作用下土壤垂直剖面水分温度的动态变化特征，对进一步研究该地区土壤结构，沙化治理方法和维护生态平衡具有重要的理论意义和应用价值。