考虑土体非均质传热的河岸带水热耦合模型

2022-03-02张文兵沈振中徐力群周聪聪杨金孟

水科学进展 2022年1期

张文兵，沈振中，任杰，徐力群，周聪聪，杨金孟

(1. 河海大学水文水资源与水利工程科学国家重点实验室，江苏南京 210098；2. 上海海事大学海洋科学与工程学院，上海 201306；3. 河海大学水利水电学院，江苏南京 210098；4. 西安理工大学省部共建西北旱区生态水利国家重点实验室，陕西西安 710048)

河岸带是指介于水生生态系统与陆地生态系统的过渡带[1-2]，是河流系统的重要屏障，能有效调蓄洪水、削减污染和保护水土环境，对维护河流生态系统健康具有重要作用[3-6]。河流地表水通过河岸带沉积层与地下水发生水热交换的区域被称为河岸带潜流层[7]，它是河岸带边缘效应和功能的重要体现区[8]。河岸带潜流层地下水与地表水交换过程时刻伴随着能量的传递和交换，温度作为能量的直观载体，是反映地表水与地下水交换过程时空变化的主要表征因子，具有无污染、易观测和受扰动小等特点，是理想的天然示踪剂[9]。近年来，随着温度自动化观测设备的革新及水热耦合机理研究的深入，对河岸带地表水与地下水交换速率及过程模式的分析已从传统的水文学及水动力学方法逐渐发展到温度示踪法[10]。现有的河岸带水热交换过程研究多侧重于试验手段[11-15]，对河岸带水热耦合模型的研究鲜见报道。因此，有必要构建合适的河岸带水热耦合模型以评估其内部水热动态变化过程，有助于揭示河岸带内水分运移和热量交换规律，对反演河岸带污染物的迁移过程以及揭示河岸带的环境效应具有重要意义。

土体热性质参数作为水热耦合模型的重要驱动因子，是决定模型是否成功和准确的关键[16]，其中，有效导热系数是影响和决定土体在传热过程中温度分布的重要参数[17]。在现有的水热耦合模拟中，土体有效导热系数均被视作固定值，例如：Su等[18]基于GA-VS2DH开发的水热耦合模拟方法，研究了大汶河和秦淮河2个不同介质类型的河岸带水热交换过程；Ren等[19]基于COMSOL构建了二维河岸带饱和-非饱和水热耦合模型，分析了河岸带温度场在不同季节的变化规律；Ren等[20]基于洞庭湖河岸带某典型断面的实测温度和水位资料，验证了所提水热耦合模型的合理性，这些模型均未考虑水热耦合模拟中土体的非均质传热问题。已有研究表明[16-17]，土体有效导热系数与土体类型、质地、粒径分布、孔隙度以及饱和度等因素均有关，并且孔隙度和饱和度的影响最大，这对涉及非饱和问题的河岸带水热交换研究至关重要。因此，对于河岸带水热耦合模拟研究，还需考虑因土体各部分含水率不同而导致的非均质传热问题。土体有效导热系数模型是通过建立导热系数与含水率之间的关系，进而实现对不同饱和状态土体有效导热系数的预测，若将其引入河岸带水热耦合模型中，则能实现对河岸带水热交换模拟过程中土体非均质传热问题的考虑。

本文在饱和-非饱和渗流及多孔介质传热理论基础上，引入土体有效导热系数模型，构建考虑土体非均质传热的河岸带水热耦合模型，给出其在COMSOL商业有限元软件中的实现方法，并通过野外原型观测资料验证和对比分析不同有效导热系数模型下河岸带水热耦合模型的模拟效果，以期为河岸带水热耦合模型的构建及有效导热系数模型的选取提供借鉴。

1 河岸带水热耦合数学模型

1.1 流动方程

河岸带饱和-非饱和瞬态渗流场方程采用Richards方程描述：

(1)

式中：ρw为水的密度，kg/m3；Cm为容水度，m-1；g为重力加速度，m/s2；Se为土体的相对饱和度；Ss为弹性贮水率，Pa-1；p为压力，Pa；t为时间，s；∇为拉普拉斯算子；θ为体积含水率，m3/m3；ks为饱和渗透系数，m/s；kr(θ)为非饱和带相对渗透系数，m/s，是体积含水率的函数；μ(T)为水的动力黏度，Pa·s，μ(T)=0.000 024 24×10247.8/(T+133.16)，T为温度， ℃；z为计算点高程，m；Qm为渗流场源汇项，kg/(m3·s)。

土壤水分特征采用van Genuchten模型描述：

θ=θr+Se(θs-θr)

(2)

(3)

(4)

(5)

式中：θr为残余体积含水率，m3/m3；θs为饱和体积含水率，m3/m3；hc为压力水头，m；α为水分特征曲线进气值的倒数，m-1；β为水分特征曲线坡度的指示参数，通过拟合土壤水分特征曲线得到，m=1-1/β。

1.2 热量交换方程

河岸带中水体与土体之间的热量交换可采用下式表示：

(6)

式中：Cw为水的体积热容，J/(m3·℃)；n为孔隙度；Cs为土体的体积热容，J/(m3·℃)；Keff为土体有效导热系数，W/(m·℃)；DHi,j为水动力弥散系数，m2/s；u为平均流速，m/s，u=v/θ，v为Darcy渗流速度；Qs为温度场源汇项，W/m3。

式(6)左边表示温度在变饱和条件下随时间的变化，即为非稳态项；等式右边第1项表示热传导项，第2项表示热弥散项，第3项表示热对流项。

近两年外部环境不好，但对于黑龙江销售来说，最大的影响来自当地非法倒油屡禁不止。“东北地区高峰时一天能进来50台油槽车、合计1500吨的量，算下来对黑龙江销售整体销量影响能达到10%。”洪松涛说，虽然今年年初各个片区联合当地执法部门对其实施打压，但风声一过又出来了。

其中，水动力弥散系数可表示为

(7)

式中：αT为横向弥散度，m；|v|为流速矢量的大小，m/s；δij为克里格常量，当i=j时为1，否则为0；αL为纵向弥散度，m；vi、vj分别为速度矢量的第i个和第j个分量，m/s。

1.3 有效导热系数模型

目前，有关预测土体有效导热系数的模型众多，这些模型可分为经验模型和理论模型。相较于经验模型，理论模型往往预测精度不高、公式复杂且部分参数不易确定，难以直接应用[16-17]。在经验模型方面，苏李君等[21]总结了部分具有代表性的土体有效导热系数模型，并提出了优化改进模型。表1总结和归纳了文献[21]中所涉及的土体有效导热系数模型，以期应用于河岸带水热耦合模型中，并比较在河岸带水热交换模拟中的表现。

表1 土体有效导热系数经验模型总结Table 1 Summary of thermal conductivity empirical model of soils

2 河岸带水热耦合模型实现方法

针对构建的河岸带水热耦合模型，采用适合多场耦合计算的COMSOL软件实现模型的开发及应用[26-27]。对于河岸带饱和-非饱和渗流问题，选用COMSOL内置地下水流模块中的Richards方程求解计算，而对于河岸带传热问题，COMSOL内置的多孔介质传热模块未能考虑土体有效导热系数模型，因此需借助相关二次开发接口来实现对土体非均质传热的考虑。这里采用COMSOL提供的自定义偏微分方程(PDE)模块用于等效代替多孔介质传热模块求解河岸带温度变化，COMSOL提供的系数型PDE如下：

(8)

式中：M为因变量；ea、da、c、α、ϒ、β、a和f均为自定义系数。

上述自定义系数可根据用户的需求定义为常数或者不同类型的函数，并且函数可以是不同阶次的导数，既可以时间相关，也可以空间相关，其自变量可以是独立变量，也可以是求解变量本身。为利用系数型PDE等效热量交换方程，需将式(6)按照式(8)的形式进行整理，然后在相应的编辑框中定义各系数(即ea=0,M=T,da=θCw+(1-n)Cs,c=Keff+θCwDHi,j,α=0,γ=0,β=θCwu,a=0,f=0)，则修改后的式(8)可等效热量交换方程。需要注意的是，Keff主要通过θ(θ=nSr)建立与孔隙度和饱和度(Sr)的关系，反映非均质参数的空间分布。因此，需增加预置储存模块用于计算和储存每个时间步内不同空间位置的有效导热系数，该过程可以通过在模型树的“组件”选项下定义局部变量实现。

河岸带水热耦合模型的有限元求解方法按照渗流场和温度场的直接耦合和间接耦合，可分为直接求解算法和分步求解算法。直接求解算法是通过Newton迭代直接求解式(1)和式(8)，即在每个时间步内同时求解单元节点的压力和温度，但由于水力特性对饱和度和压力的依赖性，导致渗流场方程高度非线性，并且有效导热系数模型的引入使得这种非线性程度增强，因此在利用直接求解算法时，结果往往不易收敛。为避免模型计算结果收敛性差的问题，可采用分步式求解算法对压力(p)和温度(T)进行解耦计算，即将渗流场与温度场划分为2个独立的场进行求解，并通过变量交换实现耦合作用。分步求解过程具体可描述为：首先，初始化变量，以获得初始压力(p0)、温度(T0)和有效导热系数(Keff0)，然后在tn+1

3 模型验证及对比分析

3.1 试验数据及模型设置

为研究美国沃克湖流域相关河流的渗漏损失，美国垦务局以及弗吉尼亚州雷斯顿地质调查局的研究人员在沃克湖流域的部分河床及河岸带埋设了温度和水位监测装置，对该区域地表水与地下水温度和水位进行长期动态监测[28](图1)。首先将带有滤网包裹的PVC管打入河岸带沉积层，然后将3～4个独立的防水温度传感器分别悬挂在距离河道中心约2.8 m和5.8 m处的PVC管中，分别测量距离地表0.95 m、1.40 m、2.00 m和1.50 m、1.95 m、2.30 m处沉积层的温度变化。此外，在河道中还布设有水温、水位监测仪用于实时监测河道水位及水温变化。地表温度则是通过埋设在河岸带浅层土体中的温度传感器进行记录。上述温度与水位数据在收集过程中由数据记录仪控制、记录和存储，均为每1 h记录一次。

图1 水位、温度监测仪器布设示意Fig.1 Layout diagram of water level and temperature monitoring instruments

图2 2012年福克斯1号河道水位、水温和地表温度时序资料Fig.2 Time-series data of water level,water temperature and land surface temperature of the Fox 1 River

基于构建的河岸带水热耦合模型及其在COMSOL软件中的实现方法，结合河岸带水位和温度原型观测资料，可以实现对考虑土体非均质传热的河岸带水热耦合模型有限元求解。图3为河岸带二维模型计算区域示意图。根据渗透性的不同，计算区域大致可分为图示的3个区域。对于河岸带饱和-非饱和渗流场，模型左边界(AF)、右边界(BC)和空气与土壤接触界面(CD、EF)为无流动边界；河水与土壤接触界面(DE)为变水位边界条件，是通过将实测水位以插值函数的形式定义在模型边界上；底部边界(AB)为透水层边界；渗流场的初始条件是通过将模拟前一时段的实测压力水头进行线性插值施加于计算域。对于河岸带温度场，模型左边界(AF)、右边界(BC)和底部边界(AB)均假定为绝热边界；空气与土壤接触界面(CD、EF)及河水与土壤接触界面(DE)分别为地表温度边界和变水温边界，同样是通过实测地表温度和河道水温以插值函数的形式施加于模型边界上；温度场的初始条件为模拟前一时段各区域平均温度。在网格划分方面，采用COMSOL预定义的三角形网格单元，并将网格单元大小校准为较细化的流体动力学尺寸，共产生7 877个网格节点和15 371个网格单元。

图3 模型计算区域示意Fig.3 Schematic diagram of model calculation area

参考相关文献[28-30]，给出河岸带水热耦合计模型算参数，如表2所示。其中，参数ks、θr、α、β、n、Com、Cs和Cw参照文献[28]给出；Cclay、Csand、Csilt和θs根据文献[29]中的砂壤土取值；ρb取自文献[30]中土壤质地、孔隙率和粒径分布相近的土体；横向、纵向热弥散度均取0.01 m；水的导热系数取0.598 W/(m·℃)；水的密度取1 000 kg/m3。

表2 河岸带水热耦合模型计算参数Table 2 Calculation parameters for hydrothermal coupling model of riparian zone

3.2 非均质传热建模有效性验证

为验证利用PDE模块代替多孔介质传热模块计算河岸带土体非均质传热过程的有效性，将前述6种土体有效导热系数模型分别考虑至河岸带水热耦合模型中，计算得出不同模型土体有效导热系数与体积含水率变化关系的数值解，并与相应的模型解析解对比，如图4所示。由图4可见，不同模型下土体有效导热系数与含水率变化关系的数值解与解析解吻合度高，并且具有较好的一致性。

此外，从图4显示的土体有效导热系数与体积含水率的关系还可以看出，不同饱和状态土体所对应的有效导热系数明显不同，并且当土体处于非饱和状态时(即θ<0.41)，不同模型预测的有效导热系数亦存在显著差异，表明所构建的河岸带水热耦合模型能够考虑因土体处于不同饱和状态而导致的非均质传热问题，并且能够实现对不同土体有效导热系数模型的准确计算。

图4 土体有效导热系数模型的解析解与数值解对比Fig.4 Comparison between analytical solution and numerical solution of soil effective thermal conductivity models

3.3 河岸带水热耦合模型验证及对比分析

为验证河岸带水热耦合模型模拟效果，拟将模拟的渗流场和温度场结果与现场监测结果进行对比。由于试验过程中河岸带压力传感器发生故障，未能有效获得计算时段内河岸带监测井的水头数据，给直接验证河岸带水头变化带来困难。这里借鉴Ren等[19]及Naranjo和Smith[28]对渗流场的验证方法，即在缺失实测压水头变化数据的情况下，通过比较模型计算的河岸带地下水渗流速度与实测河流水位变化规律的一致性来定性反映模拟的渗流场合理与否。此外，还可以将模拟的河岸带地下水渗流速度与水动力学法计算结果进行比较，进而分析渗流场计算结果的可靠性。图5给出了河流水位及河岸带地下水渗流速度时序变化曲线。由图5可见，基于水热耦合模型模拟得到的河岸带地下水渗流速度与河流水位变化规律基本一致，并且模拟的地下水渗流速度与水动力学法计算结果较为吻合。因而，所构建的水热耦合模型能合理反映河岸带内渗流场变化。

图5 2012年福克斯1号河流水位及河岸带地下水渗流速度时序变化曲线Fig.5 Time-series variation curves of river level and riparian groundwater velocity of the Fox 1 River

为对河岸带温度场作进一步验证，并比较不同有效导热系数模型下河岸带水热耦合模型的模拟效果，图6给出了不同模型下河岸带各测点温度模拟值与实测值对比。由图6可见，不同有效导热系数模型下河岸带水热耦合模型模拟的各测点温度变化效果存在差异，其中，Campbell有效导热系数模型下的河岸带水热耦合模型模拟的各测点温度值明显低于实测结果，表明该模型对河岸带土体传热效果预测偏低。由图4显示的土体有效导热系数与体积含水率的关系可知，在相同体积含水率下，Campbell模型预测的有效导热系数明显低于其他几种模型，这是导致该模型对河岸带温度变化预测不准的最主要原因。与其他土体有效导热系数模型相比，Campbell模型相对简单，所涉及的模型参数仅包括土体堆积密度、含水率和黏土质量分数，而土体有效导热系数主要与孔隙度和饱和度有关，该模型并未考虑孔隙度和饱和度的影响，使得土体有效导热系数预测精度不高。此外，当土体有效导热系数按饱和状态取为固定值时(Keff=2.01 W/(m·℃))，模拟得到的各测点温度值与实际偏差较大。相比之下，Johansen、Cté、Lu、改进Cté和改进Lu模型均表现出了较好的模拟效果，这些模型均是在考虑孔隙度和饱和度基础上做出的改进，并且考虑的影响因素较Campell模型更加丰富，因此在预测土体有效导热系数时具有较高的精度，进而使得模拟的河岸带温度值与实测值较为接近。由此可以看出，在利用水热耦合模型模拟河岸带水热交换过程时，选择合适的土体有效导热系数模型至关重要。

图6 不同模型下各测点温度模拟值与实测值对比Fig.6 Comparison of simulated and measured temperature values at each point under different models

为定量评价不同有效导热系数模型下河岸带水热耦合模型的模拟效果，引入均方根误差(ERMS)、决定系数(R2)和平均相对误差(EMR)作为模型性能评价指标[31]。表3给出了不同模型下河岸带各测点温度模拟值与实测值的对比统计结果。由表3可见，Campbell模型温度模拟结果的ERMS值为0.57～2.72 ℃，R2值为0.55～0.96，EMR值为3.91%～14.39%，各项性能评价指标均处于相对较低水平，该模型下的水热耦合模型低估了河岸带水热交换作用。当不考虑土壤热性质发生变化时(即土体有效导热系数取固定值)，其温度模拟结果的ERMS值为0.74～2.05 ℃，R2值为0.43～0.96，EMR值为3.54%～19.44%，各评价指标结果表现较差，同样不能较好地反映河岸带水热交换过程。相比之下，其他几种土体有效导热系数模型模拟结果的各项性能指标表现较好，其中，Johansen模型温度模拟结果的ERMS值为0.59～1.16 ℃，R2值为0.91～0.95，EMR值为3.35%～7.72%，各性能评价指标表现最好。

表3 不同模型下各测点温度模拟结果的评价指标统计Table 3 Statistics of evaluation indexes of temperature simulation results of each point under different models

由于不同模型模拟的河岸带温度变化在不同测点效果表现不一，难以综合评价各有效导热系数模型的性能表现。因此，将河岸带6个监测点处的温度实测值与模拟值进行整体比较(图7)，得到了反映模型整体性能表现的评价指标统计结果(表4)。由图7和表4可见，在整体分析情况下，Johansen、Cté、Lu、改进Cté和改进Lu模型的模拟结果偏差总体较小，这些模型均能较好地反映河岸带水热动态变化过程。相比之下，当不考虑土体非均质传热或采用Campbell有效导热系数模型时，河岸带水热耦合模型模拟效果较差，进一步表明河岸带水热耦合模型的模拟效果与有效导热系数的预测准确与否密切相关。总的来看，基于Johansen模型的河岸带水热耦合模型在模拟河岸带水热交换过程中效果表现最佳，该模型能较为真实地反映河岸带内部水热动态变化过程。

表4 整体分析下不同模型温度模拟结果的评价指标统计Table 4 Statistics of evaluation indexes of temperature simulation results of different models under global analysis

图7 不同有效导热系数模型下河岸带温度模拟值与实测值的整体比较Fig.7 Global comparison of measured and simulated temperature values of riparian zone under different models

图8给出了基于Johansen模型模拟得到的河岸带体积含水率及温度分布图。由图8可见，河岸带地下水与河流地表水进行侧向潜流交换作用并发生热量交换，其温度场在太阳辐射和地表水入侵双重影响下重新分布。河岸带温度在水平方向上大致可以分为高温区、中温区和低温区，并且在垂直方向上呈现“上暖下凉”的现象，这与李英玉等[7]和刘东升等[13]此前野外试验所揭示规律一致。