CN111368460B - 一种模拟潜流带温度场时空分布规律的方法 - Google Patents

Abstract

本发明公开了一种模拟潜流带温度场时空分布规律的方法，包括：运用软件选择模型方程和函数，模型方程包括潜流带的渗流场模型方程、潜流带的温度场模型方程，函数包括土壤水力函数；先将河道断面的几何模型导入到软件当中，再对河道断面进行区域划分，然后定义模型方程，并设置河道断面模型的边界条件；对河道断面进行网格剖分并插入域点探针，经过计算得到潜流带温度场时空分布规律；采用均方根误差、决定系数和相对误差对潜流带温度场时空分布规律的模拟结果进行精度评估。

Description

一种模拟潜流带温度场时空分布规律的方法

技术领域

本发明属于潜流带温度模拟方法技术领域，涉及一种模拟潜流带温度场时空分布规律的方法。

背景技术

潜流带是河流底部及其周边连接地表水体和地下水的饱和沉积物层，在此区域会发生地表水与地下水的交换，地表水与地下水的交换作用伴随着热量的传递，对潜流带温度场的时空分布规律进行分析是研究潜流带动植物的生存环境、微生物的新陈代谢、土壤的呼吸等水环境问题的基础。

由于太阳辐射和大气温度的影响，潜流带中不同点的温度时序曲线具有每日波动、随季节增减、随土壤含水量变化而变化等特征。随着自动化温度观测设备改进、信号分析与处理方法的发展以及计算机数值计算与仿真分析能力的提升，潜流带温度的监控与分析得到了硬件上的支持。在研究潜流带温度分布时，需要进行的土壤导热系数的确定，现有的很多模型无法将土壤的饱和与非饱和性质表现出来，进而造成潜流带的模拟温度不精确。

相关研究表明，潜流带内发生着强烈的生物地球化学反应，而温度在大多反应中起到了催化的作用，这些反应的发生对于改善水环境质量有着重要的影响，同时潜流带也是很多鱼类生物生存以及繁衍后代的主要场所。目前，中国正大力推广的海绵城市建设中的雨水花园、人工湿地等措施也涉及到一些和潜流带生物环生存境息息相关的问题，但是并没有对潜流带温度场时空分布规律进行模拟的方法。

发明内容

本发明的目的是提供一种模拟潜流带温度场时空分布规律的方法，能得到精确的潜流带温度场时空分布规律。

本发明所采用的技术方案是，一种模拟潜流带温度场时空分布规律的方法，具体按照以下步骤实施：

步骤1、运用软件选择河道断面的模型方程和函数，模型方程包括潜流带的渗流场模型方程、潜流带的温度场模型方程，函数包括土壤水力函数；

步骤2、先将河道断面的几何模型导入软件中，再对河道断面进行区域划分，然后定义模型方程，并设置河道断面模型的边界条件；

步骤3，对河道断面进行网格剖分并插入监测点，计算得到潜流带温度场时空分布规律；

步骤4、采用均方根误差、决定系数和相对误差对潜流带温度场时空分布规律的模拟结果进行精度评估。

本发明的特点还在于：

步骤1中运用COMSOL软件选择模型方程和函数，潜流带的渗流场方程为：

式中：θ为含水量，K_v为水力传导率；h为压力水头，H为总水头，C为土体容水度，

σ为多孔介质的孔隙率，S_s为弹性贮水率，Q_s为渗流源汇项，

为拉普拉斯方程，D_T为水动力弥散系数，t为时间；

潜流带的温度场方程为：

上式中，c为土体比热容，ρ为土体等效密度，

为拉普拉斯方程，λ为土体等效导热系数，c_w为水的比热容，ρ_w为水的密度，v为水的平均流速，T为初始温度场，Q_h为温度场源汇项；

土壤水力函数为：

上式中，θ_s为饱和含水率，θ_r为残余含水率，n和α为VG模型参数，其中m＝1-1/n，l为经验拟合参数，S_e为相对饱和度，S_e＝(θ-θ_r)/(θ_s-θ_r)，K_s为饱和水力传导率。

步骤2具体包括：

步骤2.1、将河道断面的几何模型导入COMSOL软件中；

步骤2.2、根据土质情况的的不同，将河道断面分为两个区域，并设置两个区域的材料属性参数和介质参数；

步骤2.3、将土体等效导热系数λ修改为Ren模型中导热系数随含水率及土壤参数变化的表达式：

λ＝(λ_sat-λ_dry)·Ke+λ_dry (5)；

λ_sat＝0.53C_sand+0.1γ_d (6)；

λ_dry＝0.087C_sand+0.019γ_d (7)；

Ke＝exp(ω-θ^-μ) (8)；

ω＝0.493C_sand+0.86C_silt+0.014C_om+0.778 (9)；

μ＝0.736C_clay+0.006C_om+0.222 (10)；

当砂砾质量分数C_sand的取值为：0<C_sand<1时，干土容重11<γ_d<20，γ_d和ρ_b之间的关系式为：

γ_d＝g·ρ_b (11)；

上式中，λ_dry为干土导热系数，λ_sat为饱和土导热系数，Ke为克斯特差值系数，γ_d为干土容重，ω和μ为曲线形状因子，θ为土壤含水率，C_sand为砂土质量分数，C_silt为壤土质量分数，C_clay为粘土质量分数，C_om为有机质质量比,ρ_b为体积密度；

对潜流带的温度场进行定义：选择达西速度场作为速度场，选择水作为流体材料，并输入干土容重γ_d、体积密度ρ_b、砂粒质量分数C_sand、黏粒质量分数C_silt、粉粒质量分数C_cla、有机质质量比C_om；并设置横向和纵向的热弥散性D_L、D_T；

步骤2.4、对潜流带的渗流场方程进行定义：选择水作为流体材料，并输入饱和液体体积分数θ_s和残余液体体积分数θ_r；选择水力传导率作为渗透率模型，并输入饱和水力传导率K_s；

步骤2.5、设置河道断面模型的边界条件，边界条件包括渗流场边界条件和温度场边界条件，渗流场边界条件为：河床最高水位至最低水位之间为水头边界条件，最高水位之上以及断面左右两边的边界为无流通边界条件，河道断面底部为透水边界条件；

温度场边界条件为：河床最高水位至最低水位之间为水温边界条件，断面底部边界以及左右两边界为绝热边界，最高水位之上与大气接触的部位为大气边界条件；

先将实测环境温度、水头和水温随时间变化的序列导入河道断面模型中，再将气温序列施加在大气边界上，同时将水温、水头序列分别导入水温边界、水头边界。

步骤3具体包括：

步骤3.1、采用自由剖分三角形网格中的流体动力学网格方式对河道断面进行多次剖分，直至网格质量接近于1；

步骤3.2、设定监测点并插入域点探针，将监测点的位置坐标导入河道断面模型中，通过计算得到潜流带温度场时空分布规律。

步骤4中均方根误差、决定系数和相对误差计算公式分别为：

上式中，O_i为试验值，S_i为模拟值，n为样本个数，

为试验平均值。

本发明的有益效果是：

一种模拟潜流带温度场时空分布规律的方法，运用软件建立河道断面模型，经过精确模拟得到潜流带温度场的时空分布规律数据，有利于潜流带生物生存环境研究的进行，对潜流带动植物的生存环境、微生物的新陈代谢、土壤的呼吸等生态环境方面的研究有重要意义。

附图说明

图1是本发明一种模拟潜流带温度场时空分布规律的方法的河道断面图；

图2是本发明一种模拟潜流带温度场时空分布规律的方法中监测点T6-0.10m模拟温度与实测温度对比图；

图3是本发明一种模拟潜流带温度场时空分布规律的方法中监测点T6-0.20m模拟温度与实测温度对比图；

图4是本发明一种模拟潜流带温度场时空分布规律的方法中监测点T6-0.50m模拟温度与实测温度对比图；

图5是本发明一种模拟潜流带温度场时空分布规律的方法中监测点T6-0.75m模拟温度与实测温度对比图；

图6是本发明一种模拟潜流带温度场时空分布规律的方法中监测点T6-1.0m模拟温度与实测温度对比图；

图7是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP2-1.2m模拟温度与实测温度对比图；

图8是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP2-1.7m模拟温度与实测温度对比图；

图9是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP2-2.3m模拟温度与实测温度对比图；

图10是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP2-2.75m模拟温度与实测温度对比图；

图11是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP1-0.5m模拟温度与实测温度对比图；

图12是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP1-0.7m模拟温度与实测温度对比图；

图13是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP1-1.5m模拟温度与实测温度对比图；

图14是本发明一种模拟潜流带温度场时空分布规律的方法中监测点BP1-2.15m模拟温度与实测温度对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

一种模拟潜流带温度场时空分布规律的方法，如图1所示，具体按照以下步骤实施：

步骤1、运用软件选择河道断面的模型方程和函数，模型方程包括潜流带的渗流场模型方程、潜流带的温度场模型方程，函数包括土壤水力函数；

运用COMSOL选择模型方程和函数，潜流带的渗流场方程为：

式中：θ为土壤含水率，K_v为水力传导率；h为压力水头，H为总水头，C为土体容水度，

σ为多孔介质的孔隙率，S_s为弹性贮水率，Q_s为渗流源汇项，

为拉普拉斯方程，D_T为水动力弥散系数，t为时间；

潜流带的温度场方程为：

上式中，c为土体比热容，ρ为土体等效密度，

为拉普拉斯方程，λ为土体等效导热系数，c_w为水的比热容，ρ_w为水的密度，v为水的平均流速，T为初始温度场，Q_h为温度场源汇项；

土壤水力函数为：

上式中，θ_s为饱和含水率，θ_r为残余含水率，α和n为VG模型参数，其中m＝1-1/n，l为经验拟合参数，通常取平均值0.5，S_e为相对饱和度，S_e＝(θ-θ_r)/(θ_s-θ_r)，K_s为饱和水力传导率。

步骤2、先将河道断面的几何模型导入软件中，再对河道断面进行区域划分，然后定义模型方程，并设置河道断面模型的边界条件；

步骤2.1、将河道断面的CAD几何模型导入COMSOL软件中，如图2所示；

步骤2.2、根据土质情况的的不同，将河道断面分为zone 1和zone 2两个区域，并设置两个区域的材料属性参数和介质参数，材料类型为土壤，土壤参数包括：渗透率K，土壤动力粘度η，比热率f；介质为水，其材料属性为默认的参数；

步骤2.3、将土体等效导热系数λ为Ren模型中导热系数随含水率及土壤参数变化的表达式：

λ＝(λ_sat-λ_dry)·Ke+λ_dry (5)；

λ_sat＝0.53C_sand+0.1γ_d (6)；

λ_dry＝0.087C_sand+0.019γ_d (7)；

Ke＝exp(ω-θ^-μ) (8)；

ω＝0.493C_sand+0.86C_silt+0.014C_om+0.778 (9)；

μ＝0.736C_clay+0.006C_om+0.222 (10)；

当砂砾质量分数C_sand的取值为：0<C_sand<1时，干土容重11<γ_d<20，γ_d和ρ_b之间的关系式为：

γ_d＝g·ρ_b (11)；

上式中，λ_dry为干土导热系数，λ_sat为饱和土导热系数，Ke为克斯特差值系数，γ_d为干土容重，ω和μ为曲线形状因子，θ为土壤含水率，C_sand为砂土质量分数，C_silt为壤土质量分数，C_clay为粘土质量分数，C_om为有机质质量比，ρ_b为体积密度；

对潜流带的温度场进行定义：为了形成流热耦合选择达西速度场作为速度场，选择水作为流体材料，并输入干土容重γ_d、体积密度ρ_b、砂粒质量分数C_sand、黏粒质量分数C_silt、粉粒质量分数C_cla、有机质质量比C_om；并设置横向和纵向的热分散率D_L、D_T；

步骤2.4、对潜流带的渗流场方程进行定义：选择水作为流体材料，并输入饱和液体体积分数θ_s(m³·m^-3)和残余液体含水率θ_r(m³·m^-3)；选择水力传导率作为渗透率模型，并输入饱和水力传导率K_s(m/s)；

步骤2.5、设置河道断面模型的边界条件，边界条件包括渗流场边界条件和温度场边界条件，渗流场边界条件为：河床最高水位至最低水位之间为水头边界条件，最高水位之上以及断面左右两边的边界为无流通边界条件，河道断面底部为透水边界条件；

温度场边界条件为：河床最高水位至最低水位之间为水温边界条件，断面底部边界以及左右两边界为绝热边界，最高水位之上与大气接触的部位为大气边界条件；

先将实测环境温度、水头和水温随时间变化的序列导入河道断面模型中，再将气温序列施加在大气边界上，同时将水温、水头序列分别导入水温边界、水头边界。

步骤3，对河道断面进行网格剖分并插入监测点，计算得到潜流带温度场时空分布规律；

步骤3.1、采用自由剖分三角形网格中的流体动力学网格方式对河道断面进行多次剖分，直至网格质量接近于1；

步骤3.2、选择3个监测井，共计13个监测点并插入域点探针，将监测点的位置坐标导入河道断面模型中，通过计算得到潜流带温度场时空分布规律。

步骤4、采用均方根误差、决定系数和相对误差对潜流带温度场时空分布规律的模拟结果进行精度评估，均方根误差RMSE、决定系数R2和相对误差Re计算公式分别为：

上式中，O_i为试验值，S_i为模拟值，n为样本个数，

为试验平均值。

将运用HYDRUS软件建立Chung and Horton模型的方法与本发明的方法构建模型的方法进行比较，HYDRUS软件描述土壤水力函数的方式如下，与本发明的土壤水力函数表示方法相同：

上式中，θ_s为饱和含水率(m³·m^-3)；θ_r为残余含水率(m³·m^-3)；n和α为VG模型参数，其中m＝1-1/n；l为经验拟合参数，通常取平均值0.5；S_e为相对饱和度，S_e＝(θ-θ_r)/(θ_s-θ_r)；K_s为饱和水力传导率；

潜流带的渗流场根据土壤水分运动方程，采用修正的二维Richards方程表示为：

上式中，θ为土壤体积含水率(m³·m^-3)，h为压力水头(m)，t为时间(h)，K(θ)为非饱和土壤导水率，(m·h^-1)，z为垂直坐标，x为水平坐标；

潜流带的温度场根据热量运移方程，忽略气态水扩散而仅考虑液态水的运动对土壤热量传输的影响，二维土壤热流运动的基本方程可表示为：

上式中，C(θ)为多孔介质比热容(J·m^-3·℃^-1)，C_w为液体比热容(J·m^-3·℃^-1)，λ_ij为土壤导热率(W·m^-1·℃^-1)，q_i为水分通量(m³/h)，T为土壤温度(℃)，x为水平方向坐标(m)，z为垂直方向坐标(m)；

其土壤导热率根据Chung and Horton的经验公式确定：

λ_ij＝b₀+b₁θ+b₂θ^0.5 (19)；

式中：λ_ij为土壤导热率,W/(m·℃)；b₀、b₁、b₂为热导率函数中的系数，对于不同质地的土壤，系数有所差别。

通过以上方式，一种模拟潜流带温度场时空分布规律的方法，运用软件建立河道断面模型，经过精确模拟得到潜流带温度场的时空分布规律数据，有利于潜流带生物生存环境研究的进行，对潜流带动植物的生存环境、微生物的新陈代谢、土壤的呼吸等生态环境方面的研究有重要意义。

实施例

采用如表1所示的土壤参数，运用本发明的模拟潜流带温度场时空分布规律的方法和HYDRUS软件建立Chung and Horton模型的方法分别进行仿真。

表1土壤参数

模型验证：

一.进行野外试验：

野外实地观测是研究的基础，它可以为研究人员提供最为直观的科学认识，所获得数据资料不仅为室内试验参数的设置提供了参考依据，也是现场尺度参数校正与模型验证不可缺少的部分。本发明选取的数据资料来自美国弗吉尼亚州雷斯顿地质调查局(U.S.Geological Survey,Reston,Virginia)于2012年3月至2012年11月在沃克河流域(Walker River Basin)进行的现场试验。雷斯顿地质调查局研究人员在沃克河流域相关河岸带埋设了温度及水位传感器，以实现对水位和温度的动态监测。测试过程中，数据由数据记录仪控制、记录和储存，每1h记录一次。

二、模型精度评估

对本发明所布置的3个监测井中的13个监测点进行编号1～13，把所有监测点的模拟数据与实测数据进行归类整理后进行模型的精确度评估。本发明采用均方根误差(RMSE)、决定系数(R²)和相对误差(Re)对模型模拟的精度进行评价，结果见表2。RMSE的值在0到+∞范围内变化，模拟效果越好，其值越小；R²为试验值与模拟值的y＝x线性回归方程的决定系数，模拟越好，其值越接大越接近于1；Re为模拟值与试验值之间的相对误差，拟合效果越好，其值越小。

表2观测点RMSE-R²-Re结果

从表2可以看出：

(1)对于三种模型精度评估方法：Ren模型的RMSE在0.80～3.01℃之间，平均值为1.60，Chung and Horton模型的RMSE在0.70～3.92℃之间，平均值为2.01；Ren模型的R²在0.48～0.97之间，平均值为0.84，Chung and Horton模型的R²在0.39～0.98之间，平均值为0.76；Ren模型的Re在0.04～0.17之间，平均值为0.09，Chung and Horton模型的Re在0.04～0.21之间，平均值为0.11。综合三种误差分析可以得到，Ren模型的误差波动范围均比Chung and Horton模型的波动范围小，说明Ren模型模拟各个点的结果更加稳定，变化幅度不大；三个平均值也相都能看出Ren模型更为精确。需要说明的是，在进行模型构建的时候，因为没有考虑水分的蒸发与地表水的补给，所以对于每个监测井最上面与最下面的两个监测点相对于中间监测点的误差偏大。通过均方根误差(RMSE)、决定系数(R²)、相对误差(Re)对模型的评价，确定Chung and Horton模型和Ren模型模拟的结果满足精度要求。

(2)对于三个监测井位置分布：在河床的T6监测井部分，如图2-6所示，相比于经典的Chung and Horton模型，基于Ren模型模拟的潜流带温度场时空分布规律误差更小，模拟更加接近于实测值，能更精确的刻画潜流带温度动态变化过程；在河岸带的BP1与BP2监测井部分，如图7-14所示，Chung and Horton模型随位置的改变模拟效果变动幅度较大，而Ren模型所模拟的结果更具有稳定性，误差不会大幅度波动。土壤导热系数的Ren模型相对于Chung and Horton模型在模拟潜流带温度的时空分布规律时误差均在模拟有效范围之内且波动小、拟合效果更加精确，这对于准确的模拟河床以及河岸潜流带温度场时空分布具有指导意义。

Claims (4)

1.一种模拟潜流带温度场时空分布规律的方法，其特征在于，具体按照以下步骤实施：

步骤1、运用软件选择河道断面模型方程和函数，所述模型方程包括潜流带的渗流场模型方程、潜流带的温度场模型方程，所述函数包括土壤水力函数；

步骤2、先将河道断面的几何模型导入软件中，再对河道断面进行区域划分，然后定义所述模型方程，并设置河道断面模型的边界条件；

步骤2.1、将所述河道断面的几何模型导入COOMSOL软件中；

步骤2.2、根据土质情况，将河道断面分为两个区域，并设置两个区域的材料属性参数和介质参数；

步骤2.3、将土体等效导热系数λ修改为Ren模型中导热系数随含水率及土壤参数变化的表达式：

λ＝(λ_sat-λ_dry)·Ke+λ_dry (5)；

λ_sat＝0.53C_sand+0.1γ_d (6)；

λ_dry＝0.087C_sand+0.019γ_d (7)；

Ke＝exp(ω-θ^-μ) (8)；

ω＝0.493C_sand+0.86C_silt+0.014C_om+0.778 (9)；

μ＝0.736C_clay+0.006C_om+0.222 (10)；

当砂砾质量分数C_sand的取值为：0<C_sand<1时，干土容重11<γ_d<20，γ_d和ρ_b之间的关系式为：

γ_d＝g·ρ_b (11)；

上式中，λ_dry为干土导热系数，λ_sat为饱和土导热系数，Ke为克斯特差值系数，γ_d为干土容重，ω和μ为曲线形状因子，θ为土壤含水率，C_sand为砂土质量分数，C_silt为壤土质量分数，C_clay为粘土质量分数，C_om为有机质质量比,ρ_b为体积密度；

对所述潜流带的温度场进行定义：选择达西速度场作为速度场，选择水作为流体材料，并输入干土容重γ_d、体积密度ρ_b、砂粒质量分数C_sand、黏粒质量分数C_silt、粉粒质量分数C_cla、有机质质量比C_om；并设置横向和纵向的热弥散性D_L、D_T；

步骤2.4、对所述潜流带的渗流场方程进行定义：选择水作为流体材料，并输入饱和液体体积分数θ_s和残余液体体积分数θ_r；选择水力传导率作为渗透率模型，并输入饱和水力传导率K_s；

步骤2.5、设置所述河道断面模型的边界条件，所述边界条件包括渗流场边界条件和温度场边界条件，所述渗流场边界条件为：河床最高水位至最低水位之间为水头边界条件，最高水位之上以及断面左右两边的边界为无流通边界条件，所述河道断面底部为透水边界条件；

所述温度场边界条件为：河床最高水位至最低水位之间为水温边界条件，断面底部边界以及左右两边界为绝热边界，最高水位之上与大气接触的部位为大气边界条件；

先将实测环境温度、水头和水温随时间变化的序列导入河道断面模型中，再将气温序列施加在大气边界上，同时将水温、水头序列分别导入水温边界、水头边界；

步骤3，对所述河道断面进行网格剖分并插入监测点，计算得到潜流带温度场时空分布规律；

步骤4、采用均方根误差、决定系数和相对误差对所述潜流带温度场时空分布规律的模拟结果进行精度评估。

2.如权利要求1所述的一种模拟潜流带温度场时空分布规律的方法，其特征在于，步骤1中运用COMSOL选择模型方程和函数，所述潜流带的渗流场方程为：

式中：θ为含水量，K_v为水力传导率；h为压力水头，H为总水头，C为土体容水度，

σ为多孔介质的孔隙率，S_s为弹性贮水率，Q_s为渗流源汇项，

为拉普拉斯方程，D_T为水动力弥散系数，t为时间；

所述潜流带的温度场方程为：

上式中，c为土体比热容，ρ为土体等效密度，

为拉普拉斯方程，λ为土体等效导热系数，c_w为水的比热容，ρ_w为水的密度，v为水的平均流速，T为初始温度场，Q_h为温度场源汇项；

所述土壤水力函数为：

上式中，θ_s为饱和含水率，θ_r为残余含水率，α和n为VG模型参数，其中m＝1-1/n，l为经验拟合参数，S_e为相对饱和度，S_e＝(θ-θ_r)/(θ_s-θ_r)，K_s为饱和水力传导率。

3.如权利要求1所述的一种模拟潜流带温度场时空分布规律的方法，其特征在于，步骤3具体包括：

步骤3.1、采用自由剖分三角形网格中的流体动力学网格方式对所述河道断面进行多次剖分，直至网格质量接近于1；

步骤3.2、选择监测点并插入域点探针，将所述监测点的位置坐标导入河道断面模型中，通过计算得到潜流带温度场时空分布规律。

4.如权利要求1所述的一种模拟潜流带温度场时空分布规律的方法，其特征在于，步骤4中所述均方根误差、决定系数和相对误差计算公式分别为：

上式中，O_i为试验值，S_i为模拟值，n为样本个数，

为试验平均值。

Images