-
河岸带是指介于水生生态系统与陆地生态系统的过渡带[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商业有限元软件中的实现方法，并通过野外原型观测资料验证和对比分析不同有效导热系数模型下河岸带水热耦合模型的模拟效果，以期为河岸带水热耦合模型的构建及有效导热系数模型的选取提供借鉴。
-
河岸带饱和-非饱和瞬态渗流场方程采用Richards方程描述：
$$ {\rho _{\rm{w}}}\left( {\frac{{{C_{\rm{m}}}}}{{{\rho _{\rm{w}}}g}} + {S_{\rm{e}}}{S_{\rm{s}}}} \right)\frac{{\partial p}}{{\partial t}} + \nabla \left( { - {\rho _{\rm{w}}}\frac{{{k_{\rm{s}}}{k_{\rm{r}}}(\theta )}}{{\mu (T)}}\nabla \left( {p + {\rho _{\rm{w}}}gz} \right)} \right) = {Q_{\rm{m}}} $$ (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 02424×10247.8/(T+133.16)，T为温度，℃；z为计算点高程，m；Qm为渗流场源汇项，kg/(m3·s)。
土壤水分特征采用van Genuchten模型描述：
$$ \theta = {\theta _{\rm{r}}} + {S_{\rm{e}}}\left( {{\theta _{\rm{s}}} - {\theta _{\rm{r}}}} \right) $$ (2) $$ {S_{\rm{e}}} = \frac{1}{{{{\left( {1 + {{\left( {\alpha {h_{\rm{c}}}} \right)}^\beta }} \right)}^m}}} $$ (3) $$ {C_{\rm{m}}} = \frac{{\alpha m}}{{1 - m}}\left( {{\theta _{\rm{s}}} - {\theta _{\rm{r}}}} \right)S_{\rm{e}}^{1/m}{\left( {1 - S_{\rm{e}}^{1/m}} \right)^m} $$ (4) $$ {k_{\rm{r}}}(\theta ) = S_{\rm{e}}^{1/2}{\left[ {1 - {{\left( {1 - S_{\rm{e}}^{1/m}} \right)}^m}} \right]^2} $$ (5) 式中：θr为残余体积含水率，m3/m3；θs为饱和体积含水率，m3/m3；hc为压力水头，m；α为水分特征曲线进气值的倒数，m-1；β为水分特征曲线坡度的指示参数，通过拟合土壤水分特征曲线得到，m=1-1/β。
-
河岸带中水体与土体之间的热量交换可采用下式表示：
$$ \left[ {\theta {C_{\rm{w}}} + (1 - n){C_{\rm{s}}}} \right]\frac{{\partial T}}{{\partial t}} = \nabla \left( {{K_{{\rm{eff}}}}\nabla T} \right) + \nabla \theta {C_{\rm{w}}}{D_{{\rm{H}}i, j}}\nabla T - \nabla \theta {C_{\rm{w}}}uT + {Q_{\rm{s}}} $$ (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项表示热对流项。
其中，水动力弥散系数可表示为
$$ {D_{{\rm{H}}i, j}} = {\alpha _{\rm{T}}}|\mathit{\boldsymbol{v}}|{\delta _{ij}} + \frac{{\left( {{\alpha _{\rm{L}}} - {\alpha _{\rm{T}}}} \right){v_i}{v_j}}}{{|\mathit{\boldsymbol{v}}|}} $$ (7) 式中：αT为横向弥散度，m；|v|为流速矢量的大小，m/s；δij为克里格常量，当i=j时为1，否则为0；αL为纵向弥散度，m；vi、vj分别为速度矢量的第i个和第j个分量，m/s。
-
目前，有关预测土体有效导热系数的模型众多，这些模型可分为经验模型和理论模型。相较于经验模型，理论模型往往预测精度不高、公式复杂且部分参数不易确定，难以直接应用[16-17]。在经验模型方面，苏李君等[21]总结了部分具有代表性的土体有效导热系数模型，并提出了优化改进模型。表 1总结和归纳了文献[21]中所涉及的土体有效导热系数模型，以期应用于河岸带水热耦合模型中，并比较在河岸带水热交换模拟中的表现。
表 1 土体有效导热系数经验模型总结
Table 1. Summary of thermal conductivity empirical model of soils
模型类型 模型公式及参数 参数物理意义 模型假设 适用土质 Campbell模型[22] Keff=A+Bθ-(A-D)exp[-(Cθ)E],
A=0.65-0.78ρb+0.60ρb2, B=1.06ρb,
C=1+2.6/Cclay0.5,
D=0.03+0.1ρb2, E=4ρb为堆积密度；Cclay为黏土质量分数 土体有效导热系数仅与含水率、堆积密度及黏土质量分数有关 各种土质 Johansen模型[23] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=0.7lgSr+1 (0.05＜Sr≤0.1),
Ke=lgSr+1 (Sr>0.1),
Kdry=(0.137ρb+64.7)/(2 700-0.947ρb)Ksat和Kdry分别为土体在饱和及干燥状态下的导热系数；Kw和Ks分别为水和土体颗粒的导热系数；Sr为饱和度，Sr=θ/n 各因素对土体有效导热系数的影响均反映于量纲一化导热系数Ke中 各种土质 Côté模型[24] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=κSr/[1+(κ-1)Sr],
Kdry=χ10-ηnκ为反映土体种类对Ke—Sr关系的影响系数；χ和η为和颗粒形状有关的参数 量纲一化导热系数Ke与土体饱和度及种类均相关，且干土导热系数与孔隙率呈指数关系 各种土质 Lu模型[25] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=exp[α(1-Srα-1.33)],
Kdry=-0.56n+0.51α为由土壤质地决定的参数 干土导热系数与孔隙率呈线性关系 各种土质 改进Côté模型[21] Keff=(Ksat-Kdry)Ke+Kdry, Ksat=KwnKs1-n,
Ke=ΨSr/[1+(Ψ-1)Sr], Kdry=χ10-ηn,
Ψ=a1Cclay+a2Cslit+a3Csand+a4ComCclay、Cslit、Csand和Com分别为黏土、粉土、砂土质量分数及有机质质量比；a1、a2、a3和a4为拟合系数 Côté模型中系数κ与土壤颗粒组成及有机质含量呈一定比例关系 各种土质 改进Lu模型[21] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=exp[ζ(1-Srζ-1.33)],
Kdry=-0.56n+0.51,
ζ=b1Cclay+b2Cslit+b3Csand+b4Comb1、b2、b3和b4为拟合系数 Lu模型中系数α与土壤颗粒组成及有机质含量呈一定比例关系 各种土质 -
针对构建的河岸带水热耦合模型，采用适合多场耦合计算的COMSOL软件实现模型的开发及应用[26-27]。对于河岸带饱和-非饱和渗流问题，选用COMSOL内置地下水流模块中的Richards方程求解计算，而对于河岸带传热问题，COMSOL内置的多孔介质传热模块未能考虑土体有效导热系数模型，因此需借助相关二次开发接口来实现对土体非均质传热的考虑。这里采用COMSOL提供的自定义偏微分方程(PDE)模块用于等效代替多孔介质传热模块求解河岸带温度变化，COMSOL提供的系数型PDE如下：
$$ {e_{\rm{a}}}\frac{{{\partial ^2}M}}{{\partial {t^2}}} + {d_{\rm{a}}}\frac{{\partial M}}{{\partial t}} + \nabla ( - c\nabla M - \alpha M + \gamma ) + \beta \nabla M + aM = f $$ (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 < tfinal(tfinal为目标计算时间)范围内，将p0和T0代入式(1) 中计算出第一步的压力p1，再将T0、Keff0和代入式(8)中计算第一步的温度T1，最后将T1、p1和Keff0代入至有效导热系数模型中计算第一步的Keff1。重复上述步骤，直至计算出最后一个时间步的压力pi和温度Ti。
-
为研究美国沃克湖流域相关河流的渗漏损失，美国垦务局以及弗吉尼亚州雷斯顿地质调查局的研究人员在沃克湖流域的部分河床及河岸带埋设了温度和水位监测装置，对该区域地表水与地下水温度和水位进行长期动态监测[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 水位、温度监测仪器布设示意
Figure 1. Layout diagram of water level and temperature monitoring instruments
为验证所构建的河岸带水热耦合模型及评价不同有效导热系数模型在河岸带水热交换模拟过程中的性能表现，选取沃克湖支流福克斯1号灌渠2012年3—7月的实测数据进行分析。图 2给出了实测的河道水位、水温及地表温度时序资料。
图 2 2012年福克斯1号河道水位、水温和地表温度时序资料
Figure 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 模型计算区域示意
Figure 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
区域 Ks/(m·h-1) θs/(m3·m-3) θr/(m3·m-3) ρb/(kg·m-3) α/m-1 β n Cclay/% Csand/% Csilt/% Com/(g·kg-1) Cs/(J·m-3·℃-1) Cw/(J·m-3·℃-1) 1 0.004 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000 2 0.018 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000 3 0.003 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000 -
为验证利用PDE模块代替多孔介质传热模块计算河岸带土体非均质传热过程的有效性，将前述6种土体有效导热系数模型分别考虑至河岸带水热耦合模型中，计算得出不同模型土体有效导热系数与体积含水率变化关系的数值解，并与相应的模型解析解对比，如图 4所示。由图 4可见，不同模型下土体有效导热系数与含水率变化关系的数值解与解析解吻合度高，并且具有较好的一致性。
图 4 土体有效导热系数模型的解析解与数值解对比
Figure 4. Comparison between analytical solution and numerical solution of soil effective thermal conductivity models
此外，从图 4显示的土体有效导热系数与体积含水率的关系还可以看出，不同饱和状态土体所对应的有效导热系数明显不同，并且当土体处于非饱和状态时(即θ < 0.41)，不同模型预测的有效导热系数亦存在显著差异，表明所构建的河岸带水热耦合模型能够考虑因土体处于不同饱和状态而导致的非均质传热问题，并且能够实现对不同土体有效导热系数模型的准确计算。
-
为验证河岸带水热耦合模型模拟效果，拟将模拟的渗流场和温度场结果与现场监测结果进行对比。由于试验过程中河岸带压力传感器发生故障，未能有效获得计算时段内河岸带监测井的水头数据，给直接验证河岸带水头变化带来困难。这里借鉴Ren等[19]及Naranjo和Smith[28]对渗流场的验证方法，即在缺失实测压水头变化数据的情况下，通过比较模型计算的河岸带地下水渗流速度与实测河流水位变化规律的一致性来定性反映模拟的渗流场合理与否。此外，还可以将模拟的河岸带地下水渗流速度与水动力学法计算结果进行比较，进而分析渗流场计算结果的可靠性。图 5给出了河流水位及河岸带地下水渗流速度时序变化曲线。由图 5可见，基于水热耦合模型模拟得到的河岸带地下水渗流速度与河流水位变化规律基本一致，并且模拟的地下水渗流速度与水动力学法计算结果较为吻合。因而，所构建的水热耦合模型能合理反映河岸带内渗流场变化。
图 5 2012年福克斯1号河流水位及河岸带地下水渗流速度时序变化曲线
Figure 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、Côté、Lu、改进Côté和改进Lu模型均表现出了较好的模拟效果，这些模型均是在考虑孔隙度和饱和度基础上做出的改进，并且考虑的影响因素较Campell模型更加丰富，因此在预测土体有效导热系数时具有较高的精度，进而使得模拟的河岸带温度值与实测值较为接近。由此可以看出，在利用水热耦合模型模拟河岸带水热交换过程时，选择合适的土体有效导热系数模型至关重要。
图 6 不同模型下各测点温度模拟值与实测值对比
Figure 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
模型 BP1-0.95 m BP1-1.40 m BP1-2.00 m ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% Johansen模型 0.84 0.95 4.31 0.87 0.95 4.21 1.16 0.92 7.20 Campbell模型 1.98 0.73 10.53 2.45 0.60 14.39 2.72 0.55 16.16 Côté模型 0.87 0.95 4.32 1.01 0.93 5.02 1.36 0.89 8.19 Lu模型 0.87 0.95 4.31 1.01 0.93 4.99 1.35 0.89 8.15 改进Côté模型 0.92 0.94 4.31 1.10 0.92 5.52 1.41 0.88 8.49 改进Lu模型 0.89 0.95 4.30 1.04 0.93 5.17 1.38 0.88 8.28 固定值(2.01) 0.86 0.95 4.95 0.74 0.96 3.54 0.99 0.94 6.19 模型 BP2-1.50 m BP2-1.95 m BP2-2.30 m ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% Johansen模型 0.98 0.91 7.72 0.66 0.94 6.53 0.59 0.95 3.35 Campbell模型 1.00 0.90 6.10 0.57 0.96 3.91 1.27 0.79 7.76 Côté模型 0.85 0.93 6.65 0.46 0.97 4.75 0.78 0.92 4.18 Lu模型 0.88 0.93 6.91 0.50 0.97 5.15 0.73 0.93 3.91 改进Côté模型 0.85 0.93 6.59 0.45 0.97 4.65 0.80 0.92 4.29 改进Lu模型 0.85 0.93 6.62 0.45 0.97 4.70 0.79 0.92 4.22 固定值(2.01) 1.45 0.56 17.45 2.05 0.43 19.44 1.20 0.82 10.73 由于不同模型模拟的河岸带温度变化在不同测点效果表现不一，难以综合评价各有效导热系数模型的性能表现。因此，将河岸带6个监测点处的温度实测值与模拟值进行整体比较(图 7)，得到了反映模型整体性能表现的评价指标统计结果(表 4)。由图 7和表 4可见，在整体分析情况下，Johansen、Côté、Lu、改进Côté和改进Lu模型的模拟结果偏差总体较小，这些模型均能较好地反映河岸带水热动态变化过程。相比之下，当不考虑土体非均质传热或采用Campbell有效导热系数模型时，河岸带水热耦合模型模拟效果较差，进一步表明河岸带水热耦合模型的模拟效果与有效导热系数的预测准确与否密切相关。总的来看，基于Johansen模型的河岸带水热耦合模型在模拟河岸带水热交换过程中效果表现最佳，该模型能较为真实地反映河岸带内部水热动态变化过程。
图 7 不同有效导热系数模型下河岸带温度模拟值与实测值的整体比较
Figure 7. Global comparison of measured and simulated temperature values of riparian zone under different models
表 4 整体分析下不同模型温度模拟结果的评价指标统计
Table 4. Statistics of evaluation indexes of temperature simulation results of different models under global analysis
有效导热系数模型 ERMS/℃ R2 EMR/% Johansen模型 0.85 0.95 5.51 Campbell模型 1.80 0.79 9.64 Côté模型 0.91 0.95 5.42 Lu模型 0.91 0.95 5.48 改进Côté模型 0.94 0.94 5.57 改进Lu模型 0.92 0.94 5.45 固定值(2.01) 1.50 0.85 10.63 图 8给出了基于Johansen模型模拟得到的河岸带体积含水率及温度分布图。由图 8可见，河岸带地下水与河流地表水进行侧向潜流交换作用并发生热量交换，其温度场在太阳辐射和地表水入侵双重影响下重新分布。河岸带温度在水平方向上大致可以分为高温区、中温区和低温区，并且在垂直方向上呈现“上暖下凉”的现象，这与李英玉等[7]和刘东升等[13]此前野外试验所揭示规律一致。
图 8 河岸带不同时刻的体积含水率及温度分布
Figure 8. Distribution of water content and temperature in the riparian zone at different time
-
有效导热系数是影响和决定土体在传热过程中温度分布的重要参数，对河岸带水热交换的模拟有直接影响。在饱和-非饱和渗流及多孔介质传热理论基础上，通过引入土体有效导热系数模型，构建了考虑土体非均质传热的河岸带水热耦合模型，主要结论如下：
(1) COMSOL软件在材料属性、边界条件及求解器设置上极具灵活性，在利用地下水流模块模拟河岸带水流运动的基础上，通过自定义偏微分方程来实现考虑土体非均质传热的河岸带水热耦合模型开发及应用，可避免模型开发的编程求解。
(2) 考虑土体非均质传热的河岸带水热耦合模型能够实现对流场的合理反映，准确、可靠地模拟河岸带温度时空变化，对刻画河岸带水热动态变化过程具有优越性。
(3) Johansen有效导热系数模型是反映河岸带水热交换过程中土体非均质传热的最佳模型，而对于更广泛的土体有效导热系数模型性能表现，需进一步分析。
(4) 由于缺失河岸带监测井水头变化数据，目前仅对河岸带渗流场进行了间接验证。另外，因研究流域气象站点较少，缺乏降雨资料，加上模型的土-水接触边界存在一定程度的概化，导致模拟结果与实测值仍存在一定偏差，这些都需要在未来研究中进一步完善。
Hydrothermal coupling model of a riparian zone considering heterogeneous heat transfer of soil
-
摘要: 为弥补当前河岸带水热交换模拟对土体非均质传热考虑不足的缺陷，在饱和-非饱和渗流及多孔介质传热理论基础上，引入土体有效导热系数模型，构建考虑土体非均质传热的河岸带水热耦合模型。结合COMSOL软件的特点，给出河岸带水热耦合模型的实现方法及求解流程，并通过河岸带温度和水位原型观测资料验证和对比分析不同有效导热系数模型下河岸带水热耦合模型的模拟效果。结果表明：与传统不考虑土体非均质传热方法相比，该模型能够较好地反映河岸带水热交换过程。此外，基于Johansen有效导热系数模型的河岸带水热耦合模型模拟效果表现最佳，模拟结果与前人试验结果一致。研究成果可为河岸带水热交换及污染物迁移过程的深入研究提供技术支撑。Abstract: To improve the current simulations of hydrothermal exchange in riparian zones that do not take into account heterogeneous heat transfer of soil, this paper constructs a hydrothermal coupling model of riparian zones that includes heterogeneous heat transfer of soil by introducing the soil effective thermal conductivity model on the basis of saturated-unsaturated seepage and porous medium heat transfer theories. This study uses COMSOL software to implement the new method and simulate the hydrothermal coupling model of a riparian zone for different effective thermal conductivity models; the new model is verified by comparing and analyzing prototype observation data of temperature and water level in a riparian zone. The results show that the new model can better reflect the hydrothermal exchange process in a riparian zone than the traditional method that does not consider the heterogeneous heat transfer of the soil. In addition, the simulation of the hydrothermal coupling model based on Johansen's effective thermal conductivity model performs the best, and the simulation results are consistent with previous experimental results. This study can provide technical support for in-depth study of hydrothermal exchange and pollutant migration processes in riparian zones.
-
图 1 水位、温度监测仪器布设示意
Figure 1. Layout diagram of water level and temperature monitoring instruments
图 2 2012年福克斯1号河道水位、水温和地表温度时序资料
Figure 2. Time-series data of water level, water temperature and land surface temperature of the Fox 1 River
图 4 土体有效导热系数模型的解析解与数值解对比
Figure 4. Comparison between analytical solution and numerical solution of soil effective thermal conductivity models
图 5 2012年福克斯1号河流水位及河岸带地下水渗流速度时序变化曲线
Figure 5. Time-series variation curves of river level and riparian groundwater velocity of the Fox 1 River
图 6 不同模型下各测点温度模拟值与实测值对比
Figure 6. Comparison of simulated and measured temperature values at each point under different models
图 7 不同有效导热系数模型下河岸带温度模拟值与实测值的整体比较
Figure 7. Global comparison of measured and simulated temperature values of riparian zone under different models
图 8 河岸带不同时刻的体积含水率及温度分布
Figure 8. Distribution of water content and temperature in the riparian zone at different time
表 1 土体有效导热系数经验模型总结
Table 1. Summary of thermal conductivity empirical model of soils
模型类型 模型公式及参数 参数物理意义 模型假设 适用土质 Campbell模型[22] Keff=A+Bθ-(A-D)exp[-(Cθ)E],
A=0.65-0.78ρb+0.60ρb2, B=1.06ρb,
C=1+2.6/Cclay0.5,
D=0.03+0.1ρb2, E=4ρb为堆积密度；Cclay为黏土质量分数 土体有效导热系数仅与含水率、堆积密度及黏土质量分数有关 各种土质 Johansen模型[23] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=0.7lgSr+1 (0.05＜Sr≤0.1),
Ke=lgSr+1 (Sr>0.1),
Kdry=(0.137ρb+64.7)/(2 700-0.947ρb)Ksat和Kdry分别为土体在饱和及干燥状态下的导热系数；Kw和Ks分别为水和土体颗粒的导热系数；Sr为饱和度，Sr=θ/n 各因素对土体有效导热系数的影响均反映于量纲一化导热系数Ke中 各种土质 Côté模型[24] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=κSr/[1+(κ-1)Sr],
Kdry=χ10-ηnκ为反映土体种类对Ke—Sr关系的影响系数；χ和η为和颗粒形状有关的参数 量纲一化导热系数Ke与土体饱和度及种类均相关，且干土导热系数与孔隙率呈指数关系 各种土质 Lu模型[25] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=exp[α(1-Srα-1.33)],
Kdry=-0.56n+0.51α为由土壤质地决定的参数 干土导热系数与孔隙率呈线性关系 各种土质 改进Côté模型[21] Keff=(Ksat-Kdry)Ke+Kdry, Ksat=KwnKs1-n,
Ke=ΨSr/[1+(Ψ-1)Sr], Kdry=χ10-ηn,
Ψ=a1Cclay+a2Cslit+a3Csand+a4ComCclay、Cslit、Csand和Com分别为黏土、粉土、砂土质量分数及有机质质量比；a1、a2、a3和a4为拟合系数 Côté模型中系数κ与土壤颗粒组成及有机质含量呈一定比例关系 各种土质 改进Lu模型[21] Keff=(Ksat-Kdry)Ke+Kdry,
Ksat=KwnKs1-n,
Ke=exp[ζ(1-Srζ-1.33)],
Kdry=-0.56n+0.51,
ζ=b1Cclay+b2Cslit+b3Csand+b4Comb1、b2、b3和b4为拟合系数 Lu模型中系数α与土壤颗粒组成及有机质含量呈一定比例关系 各种土质 
下载: 导出CSV
表 2 河岸带水热耦合模型计算参数
Table 2. Calculation parameters for hydrothermal coupling model of riparian zone
区域 Ks/(m·h-1) θs/(m3·m-3) θr/(m3·m-3) ρb/(kg·m-3) α/m-1 β n Cclay/% Csand/% Csilt/% Com/(g·kg-1) Cs/(J·m-3·℃-1) Cw/(J·m-3·℃-1) 1 0.004 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000 2 0.018 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000 3 0.003 0.41 0.065 1 500 7.5 1.89 0.41 39 31 30 1.5 1 100 000 4 200 000
下载: 导出CSV
表 3 不同模型下各测点温度模拟结果的评价指标统计
Table 3. Statistics of evaluation indexes of temperature simulation results of each point under different models
模型 BP1-0.95 m BP1-1.40 m BP1-2.00 m ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% Johansen模型 0.84 0.95 4.31 0.87 0.95 4.21 1.16 0.92 7.20 Campbell模型 1.98 0.73 10.53 2.45 0.60 14.39 2.72 0.55 16.16 Côté模型 0.87 0.95 4.32 1.01 0.93 5.02 1.36 0.89 8.19 Lu模型 0.87 0.95 4.31 1.01 0.93 4.99 1.35 0.89 8.15 改进Côté模型 0.92 0.94 4.31 1.10 0.92 5.52 1.41 0.88 8.49 改进Lu模型 0.89 0.95 4.30 1.04 0.93 5.17 1.38 0.88 8.28 固定值(2.01) 0.86 0.95 4.95 0.74 0.96 3.54 0.99 0.94 6.19 模型 BP2-1.50 m BP2-1.95 m BP2-2.30 m ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% ERMS/℃ R2 EMR/% Johansen模型 0.98 0.91 7.72 0.66 0.94 6.53 0.59 0.95 3.35 Campbell模型 1.00 0.90 6.10 0.57 0.96 3.91 1.27 0.79 7.76 Côté模型 0.85 0.93 6.65 0.46 0.97 4.75 0.78 0.92 4.18 Lu模型 0.88 0.93 6.91 0.50 0.97 5.15 0.73 0.93 3.91 改进Côté模型 0.85 0.93 6.59 0.45 0.97 4.65 0.80 0.92 4.29 改进Lu模型 0.85 0.93 6.62 0.45 0.97 4.70 0.79 0.92 4.22 固定值(2.01) 1.45 0.56 17.45 2.05 0.43 19.44 1.20 0.82 10.73
下载: 导出CSV
表 4 整体分析下不同模型温度模拟结果的评价指标统计
Table 4. Statistics of evaluation indexes of temperature simulation results of different models under global analysis
有效导热系数模型 ERMS/℃ R2 EMR/% Johansen模型 0.85 0.95 5.51 Campbell模型 1.80 0.79 9.64 Côté模型 0.91 0.95 5.42 Lu模型 0.91 0.95 5.48 改进Côté模型 0.94 0.94 5.57 改进Lu模型 0.92 0.94 5.45 固定值(2.01) 1.50 0.85 10.63
下载: 导出CSV
-
[1] 杨平恒, 张宇, 王建力, 等. 水位变化影响下的河水-地下水侧向交互带地球化学动态[J]. 水科学进展, 2017, 28(2): 293-301. doi: 10.14042/j.cnki.32.1309.2017.02.015 YANG P H, ZHANG Y, WANG J L, et al. Influence of water level change on the geochemical dynamics of the lateral hyporheic zone between river water and groundwater[J]. Advances in Water Science, 2017, 28(2): 293-301. (in Chinese) doi: 10.14042/j.cnki.32.1309.2017.02.015 [2] WALTON C R, ZAK D, AUDET J, et al. Wetland buffer zones for nitrogen and phosphorus retention: impacts of soil type, hydrology and vegetation[J]. Science of the Total Environment, 2020, 727: 138709. doi: 10.1016/j.scitotenv.2020.138709 [3] CONANT B, ROBINSON C E, HINTON M J, et al. A framework for conceptualizing groundwater-surface water interactions and identifying potential impacts on water quality, water quantity, and ecosystems[J]. Journal of Hydrology, 2019, 574: 609-627. doi: 10.1016/j.jhydrol.2019.04.050 [4] 夏继红, 窦传彬, 蔡旺炜, 等. 河岸带蜿蜒性与植被密度对潜流驻留时间的复合效应[J]. 水科学进展, 2020, 31(3): 433-440. doi: 10.14042/j.cnki.32.1309.2020.03.013 XIA J H, DOU C B, CAI W W, et al. Compounding effect of meandering degree and vegetation density on hyporheic residence time in riparian zone[J]. Advances in Water Science, 2020, 31(3): 433-440. (in Chinese) doi: 10.14042/j.cnki.32.1309.2020.03.013 [5] 钱进, 沈蒙蒙, 王沛芳, 等. 河岸带土壤磷素空间分布及其对水文过程响应[J]. 水科学进展, 2017, 28(1): 41-48. doi: 10.14042/j.cnki.32.1309.2017.01.005 QIAN J, SHEN M M, WANG P F, et al. Spatial distribution of riparian soil phosphorus and its response to hydrologic process[J]. Advances in Water Science, 2017, 28(1): 41-48. (in Chinese) doi: 10.14042/j.cnki.32.1309.2017.01.005 [6] 李福林, 陈华伟, 王开然, 等. 地下水支撑生态系统研究综述[J]. 水科学进展, 2018, 29(5): 750-758. doi: 10.14042/j.cnki.32.1309.2018.05.015 LI F L, CHEN H W, WANG K R, et al. Comprehensive review of groundwater-dependent ecosystems[J]. Advances in Water Science, 2018, 29(5): 750-758. (in Chinese) doi: 10.14042/j.cnki.32.1309.2018.05.015 [7] 李英玉, 赵坚, 吕辉, 等. 河岸带潜流层温度示踪流速计算方法[J]. 水科学进展, 2016, 27(3): 423-429. doi: 10.14042/j.cnki.32.1309.2016.03.010 LI Y Y, ZHAO J, LYU H, et al. Investigation on temperature tracer method calculated flow rate of hyporheic layer in riparian zone[J]. Advances in Water Science, 2016, 27(3): 423-429. (in Chinese) doi: 10.14042/j.cnki.32.1309.2016.03.010 [8] 夏继红, 陈永明, 王为木, 等. 河岸带潜流层动态过程与生态修复[J]. 水科学进展, 2013, 24(4): 589-597. http://skxjz.nhri.cn/article/id/2329 XIA J H, CHEN Y M, WANG W M, et al. Dynamic processes and ecological restoration of hyporheic layer in riparian zone[J]. Advances in Water Science, 2013, 24(4): 589-597. (in Chinese) http://skxjz.nhri.cn/article/id/2329 [9] 任杰, 程嘉强, 杨杰, 等. 潜流交换温度示踪方法研究进展[J]. 水科学进展, 2018, 29(4): 597-606. doi: 10.14042/j.cnki.32.1309.2018.04.016 REN J, CHENG J Q, YANG J, et al. Advances in temperature tracer method of hyporheic exchange[J]. Advances in Water Science, 2018, 29(4): 597-606. (in Chinese) doi: 10.14042/j.cnki.32.1309.2018.04.016 [10] REN J, CHENG J Q, YANG J, et al. A review on using heat as a tool for studying groundwater-surface water interactions[J]. Environmental Earth Sciences, 2018, 77(22): 1-13. doi: 10.1007/s12665-018-7959-4 [11] WATSON J A, CARDENAS M B, FERENCZ S B, et al. The effects of floods on the temperature of riparian groundwater[J]. Hydrological Processes, 2018, 32(9): 1267-1281. doi: 10.1002/hyp.11504 [12] REN J, WANG X P, SHEN Z Z, et al. Heat tracer test in a riparian zone: laboratory experiments and numerical modelling[J]. Journal of Hydrology, 2018, 563: 560-575. doi: 10.1016/j.jhydrol.2018.06.030 [13] 刘东升, 赵坚, 吕辉. 大坝下游河岸带冬夏季水热交换特征对比[J]. 水科学进展, 2017, 28(1): 124-132. doi: 10.14042/j.cnki.32.1309.2017.01.014 LIU D S, ZHAO J, LYU H. Comparative analysis of water-heat exchange characteristics of the riparian zone downstream of dam in different seasons[J]. Advances in Water Science, 2017, 28(1): 124-132. (in Chinese) doi: 10.14042/j.cnki.32.1309.2017.01.014 [14] 姬雨雨, 陈求稳, 施文卿, 等. 水库运行对漫湾库区洲滩水热交换影响[J]. 水科学进展, 2018, 29(1): 73-79. doi: 10.14042/j.cnki.32.1309.2018.01.009 JI Y Y, CHEN Q W, SHI W Q, et al. Influence of reservoir operation on water and heat exchange in the Manwan's island[J]. Advances in Water Science, 2018, 29(1): 73-79. (in Chinese) doi: 10.14042/j.cnki.32.1309.2018.01.009 [15] 于丹青, 陈求稳, 马宏海, 等. 漫湾水库运行下库内洲滩潜流带夏季热传输特征[J]. 水利学报, 2019, 50(4): 497-505. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201904009.htm YU D Q, CHEN Q W, MA H H, et al. Study on heat transfer across island riparian zone under Manwan Reservoir operations in summer[J]. Journal of Hydraulic Engineering, 2019, 50(4): 497-505. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201904009.htm [16] ZHANG N, WANG Z Y. Review of soil thermal conductivity and predictive models[J]. International Journal of Thermal Sciences, 2017, 117: 172-183. doi: 10.1016/j.ijthermalsci.2017.03.013 [17] 张楠, 夏胜全, 侯新宇, 等. 土热传导系数及模型的研究现状和展望[J]. 岩土力学, 2016, 37(6): 1550-1562. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201606004.htm ZHANG N, XIA S Q, HOU X Y, et al. Review on soil thermal conductivity and prediction model[J]. Rock and Soil Mechanics, 2016, 37(6): 1550-1562. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201606004.htm [18] SU X R, SHU L C, CHEN X H, et al. Interpreting the cross-sectional flow field in a river bank based on a genetic-algorithm two-dimensional heat-transport method (GA-VS2DH)[J]. Hydrogeology Journal, 2016, 24(8): 2035-2047. doi: 10.1007/s10040-016-1459-y [19] REN J, ZHANG W B, YANG J, et al. Using water temperature series and hydraulic heads to quantify hyporheic exchange in the riparian zone[J]. Hydrogeology Journal, 2019, 27(4): 1419-1437. doi: 10.1007/s10040-019-01934-z [20] REN J, CHEN B, LI Y L, et al. Dynamic variation characteristics of the riparian zone temperature distribution under fluctuating water levels: field experiments and numerical modeling[J]. Journal of Coastal Research, 2020, 36(6): 1178-1196. https://www.jstor.org/stable/26952809 [21] 苏李君, 王全九, 王铄, 等. 基于土壤物理基本参数的土壤导热率模型[J]. 农业工程学报, 2016, 32(2): 127-133. https://www.cnki.com.cn/Article/CJFDTOTAL-NYGU201602019.htm SU L J, WANG Q J, WANG S, et al. Soil thermal conductivity model based on soil physical basic parameters[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(2): 127-133. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-NYGU201602019.htm [22] CAMPBELL G S. Soil physics with BASIC-transport models for soil-plant systems[M]. Amsterdam: Elsevier, 1985: 26-39. [23] JOHANSEN O. Thermal conductivity of soils[D]. Trondheim: University of Trondheim, 1975. [24] CÔTÉ J, KONRAD J M. A generalized thermal conductivity model for soils and construction materials[J]. Canadian Geotechnical Journal, 2005, 42(2): 443-458. doi: 10.1139/t04-106 [25] LU S, REN T, GONG Y, et al. An improved model for predicting soil thermal conductivity from water content at room temperature[J]. Soil Science Society of America Journal, 2007, 71(1): 8-14. doi: 10.2136/sssaj2006.0041 [26] 费康, 刘汉龙, 孔纲强, 等. 热力耦合边界面模型在COMSOL中的开发应用[J]. 岩土力学, 2017, 38(6): 1819-1826. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201706033.htm FEI K, LIU H L, KONG G Q, et al. Implementation of a thermo-bounding surface model in COMSOL[J]. Rock and Soil Mechanics, 2017, 38(6): 1819-1826. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201706033.htm [27] 焦会青, 盛钰, 赵成义, 等. 基于COMSOL软件的绿洲盐渍化土壤中多离子耦合运移模型构建[J]. 农业工程学报, 2018, 34(15): 100-107. https://www.cnki.com.cn/Article/CJFDTOTAL-NYGU201815013.htm JIAO H Q, SHENG Y, ZHAO C Y, et al. Modeling of multiple ions coupling transport for salinized soil in oasis based on COMSOL[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(15): 100-107. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-NYGU201815013.htm [28] NARANJO R C, SMITH D W. Quantifying seepage using heat as a tracer in selected irrigation canals, Walker River basin, Nevada, 2012 and 2013[R]. Virginia US: Geological Survey Nevada Water Science Center, 2016. [29] CARSEL R F, PARRISH R S. Developing joint probability distributions of soil water retention characteristics[J]. Water Resources Research, 1988, 24(5): 755-769. doi: 10.1029/WR024i005p00755 [30] KELLER T, HÅKANSSON I. Estimation of reference bulk density from soil particle size distribution and soil organic matter content[J]. Geoderma, 2010, 154(3/4): 398-406. https://www.sciencedirect.com/science/article/pii/S0016706109003644 [31] ZHANG W B, SHEN Z Z, REN J, et al. Modeling and comparative analysis of a flow and heat coupling model of the riparian zone based on thermal conductivity empirical models[J]. Journal of Hydrology, 2020, 582: 124539. https://www.sciencedirect.com/science/article/pii/S0022169419312740 -
点击查看大图
图(8) / 表 (4)
计量
- 文章访问数: 292
- HTML全文浏览量: 68
- PDF下载量: 57
- 被引次数: 0
