饱和-非饱和水分及溶质运移双重迭代耦合方法及装置

摘要

本发明提供饱和‑非饱和水分及溶质运移双重迭代耦合方法及装置，方法包括：步骤1.收集研究区的水文地质气象数据资料；将研究区中非饱和带和饱和带进行离散；采用以下步骤计算得到各应力期各时间步下的水分及溶质运移情况；步骤2.设置非饱和带厚度；步骤3.运行非饱和带水分运移模块；步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求；步骤5.计算ΔT时段内的地下水平均补给速率R；步骤6.运行饱和带地下水运动模块；步骤7.进行水分运动模块收敛性判断；步骤8.设置非饱和带溶质下边界条件；步骤9.运行非饱和带溶质运移模块；步骤10.运行饱和带溶质运移模块；步骤11.进行溶质运移模块收敛性判断。

说明书

技术领域

本发明属于地下水分及溶质计算领域，具体涉及区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法及装置。

技术背景

区域尺度饱和-非饱和水分及溶质运移规律一直是地下水资源与环境的研究重点，关系到地下水污染、土壤盐碱化、地下水开采、农业灌区灌溉排水等问题，建立区域尺度饱和-非饱和水分及溶质运移模型是一种有效的研究手段。非饱和带土壤水分是地下水系统的重要补给来源和潜层地下水的主要消耗途径，同时也是地下水污染的主要源头与缓冲区，因此需要将二者进行统一的分析与模拟。

饱和带地下水及溶质与非饱和带地下水及溶质的时空运移尺度存在较大差异，因此研究中常采用拟三维方法予以处理。基于拟三维方法的区域尺度饱和-非饱和水分及溶质模型目前仍然较少，仅有UZF-MODFLOW-MT3D、HYDRUS-MODFLOW-MT3D、Q3D等几种。UZF-MODFLOW-MT3D模型与HYDRUS-MODFLOW-MT3D模型基于松散耦合方法。该方法较为简单且易于实现，然而存在难以避免的质量误差问题，且不能保证耦合计算的精度。Q3D模型基于完全耦合方法，该方法虽然可以保证耦合计算的精度，但是实现方法极为复杂，且稳定性与收敛性较差。目前这些方法都难以得到准确有效的计算结果，从而无法获得准确、可靠的区域尺度饱和-非饱和水分及溶质运移信息。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种基于饱和-非饱和水分运动系统的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置，能够有效提高区域尺度饱和-非饱和水分及溶质的计算精度、计算效率与模型稳定性，得到准确、可靠的区域尺度饱和-非饱和水分及溶质运移信息。

本发明为了实现上述目的，采用了以下方案：

<方法>

本发明提供了一种区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，其特征在于，包括以下步骤：

步骤1.收集研究区的水文地质气象数据资料；将研究区中非饱和带离散为L个分区，饱和带离散为m×n个有限差分网格；根据研究区数据资料确定离散后各区域网格的参数信息；将这些参数信息作为输入参数，对于各个应力期各个时间步(从第一个应力期第一个时间步开始按顺序计算)，均采用以下步骤2～11计算得到对应的饱和-非饱和水分及溶质运移情况；

步骤2.设置非饱和带厚度：

设当前正处于第t个计算时间步的第p个迭代步，应力期长度为ΔT，上一时间步计算得到的非饱和带厚度为

式中，S

步骤3.运行非饱和带水分运移模块：

运行非饱和带水分运移模块，计算时间段为从t时刻至t+ΔT时刻；

步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求：

计算非饱和剖面的所有时刻的Courant数，记为Cr，若Cr＞1，则折减时间步长并返回步骤3重新计算非饱和带水分运移模型，直至满足Cr≤1；

步骤5.计算ΔT时段内的地下水平均补给速率R：

式中，r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量；依次对各区域中地下水补给量进行计算，得到该应力期中整个研究区的地下水补给速率，将之写为与地下水有限差分网格一一对应的m×n维数组形式，记为R

步骤6.运行饱和带地下水运动模块：

将地下水补给速率R

步骤7.进行水分运动模块收敛性判断：

max(|H

式中，ε

步骤8.设置非饱和带溶质下边界条件：

设当前正处于第t个计算时间步的第q个迭代步，上一时间步计算得到的饱和带潜水含水层地下水浓度为

式中，S

步骤9.运行非饱和带溶质运移模块：

运行非饱和带溶质运移模块，计算时间段为从t时刻至t+ΔT时刻，得到地下水补给量的浓度

步骤10.运行饱和带溶质运移模块：

将地下水补给浓度

步骤11.进行溶质运移模块收敛性判断：

式中，ε

优选地，本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，还可以具有以下特征：在步骤1中，研究区的水文地质气象数据资料包括：水文地质数据、气象数据、土地利用类型数据、河渠水位数据等资料；根据研究区数据资料确定的参数信息为：包含地表高程、底板高程、河渠位置及水位空间的地理信息，和地下水位、地下水浓度的初始条件信息，以及地下水饱和渗透系数、土壤水水力参数、溶质纵向弥散度等信息。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，还可以具有以下特征：在步骤2中，对于第一个应力期的第一个时间步，采用根据步骤1获取的输入参数所确定的初始条件进行计算；若为第一个迭代步，则采用上一个时间步计算得到的结果作为非饱和带厚度。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，还可以具有以下特征：在步骤6中，地下水运动模型为MODFLOW模型。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，还可以具有以下特征：在步骤8中，于第一个应力期的第一个时间步，采用根据步骤1获取的输入参数所确定的初始条件进行计算；若为第一个迭代步，则采用上一个时间步计算得到的结果作为饱和带潜水含水层地下水浓度。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法，还可以具有以下特征：在步骤10中，溶质运动模型为MT3DMS模型。

<装置>

进一步，本发明还提供了一种自动实现上述<方法>的区域饱和-非饱和水分及溶质运移计算装置，其特征在于，包括：

离散部，将研究区中非饱和带离散为L个分区，并将饱和带离散为m×n个有限差分网格，其中L、m、n均为正整数；

参数信息确定部，根据研究区的水文地质气象数据确定离散后各区域网格的参数信息；

非饱和带厚度设置部，设当前正处于第t个计算时间步的第p个迭代步，应力期长度为ΔT，上一时间步计算得到的非饱和带厚度为

式中，S

非饱和带水分运移模型计算部，运行非饱和带水分运移模型，计算时间段为从t时刻至t+ΔT时刻非饱和带水分运移；

非饱和带水盐运移一致性判断部，对非饱和带水分运移模型计算部的计算结果是否满足溶质运移计算的需求进行判断：计算非饱和剖面的所有时刻的Courant数，记为Cr，若Cr＞1，则折减时间步长，并使非饱和带水分运移计算部重新计算非饱和带水分运移模型，直至满足Cr≤1；

地下水平均补给速率计算部，计算ΔT时段内的地下水平均补给速率R：

式中，r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量；依次对各区域中地下水补给量进行计算，得到该应力期中整个研究区的地下水补给速率，将之写为与地下水有限差分网格一一对应的m×n维数组形式R

饱和带地下水运动模型计算部，将地下水补给速率R

水分运动模型收敛性判断部，采用下式对收敛情况进行判断：

max(|H

式中，ε

非饱和带溶质下边界设置部，设当前正处于第t个计算时间步的第q个迭代步，上一时间步计算得到的饱和带潜水含水层地下水浓度为

式中，S

非饱和带溶质运移模型计算部，运行非饱和带溶质运移模型，计算时间段为从t时刻至t+ΔT时刻，得到地下水补给量的浓度

饱和带溶质运移模型计算部，将地下水补给浓度

溶质运移模型收敛判断部，采用下式对收敛情况进行判断：

式中，ε

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以包括：图像生成部，根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的结果数据，生成反映饱和-非饱和水分运移情况的图像，和反映溶质运移情况的图像。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以包括：动态信息生成部，根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的结果数据，生成反映饱和-非饱和水分运移动态过程的动态信息，和反映溶质运移动态过程的动态信息。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以包括：控制部，控制离散部、参数信息确定部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部的运行。

发明的作用与效果

与现有技术相比，本发明的有益效果在于：

1.本发明将土壤水分与地下水分、土壤溶质与地下水溶质通过双重迭代的方法予以耦合，从而获得区域尺度土壤水及溶质、地下水及溶质精确的运动状态与交互过程。该方法保证了区域饱和-非饱和水分及溶质计算过程中的计算精度、迭代效率与模型鲁棒性，有利于解决地下水污染、土壤盐碱化、地下水开采等问题。

2.本发明在对非饱和带水分及溶质运移过程的描述中，时刻保持非饱和带水分运动模块与非饱和带溶质运移模块的计算时间步长严格一致，保证了状态变化较为剧烈的非饱和带水分及溶质运移过程的计算稳定性、计算精度与可靠性，有利于精确描述地表附近土壤中污染物运移、土壤盐分累积情况等，且为区域水分及溶质运移过程的研究提供了新方法与新工具。

3.本发明所提供的装置能够根据输入的实测参数自动计算得到非饱和带水分及溶质运移信息，并且还能够通过图像生成部生成反映饱和-非饱和水分运移情况的图像和反映溶质运移情况的图像，通过动态信息生成部生成反映饱和-非饱和水分运移动态过程的动态信息和反映溶质运移动态过程的动态信息，使得区域水分及溶质运移过程能够直观地呈现，有利于高效、精确地研究区域水分及溶质运移过程。

附图说明

图1为本发明实施例中涉及的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法的流程图；

图2为本发明实施例中涉及的土壤水分运动模块、土壤溶质运移模块、地下水分运动模块、地下水溶质运移模块关系示意图；

图3为本发明实施例中涉及的永联试验区位置、土地利用类型、DEM高程与监测点土质信息示意图；

图4为本发明实施例中涉及的永联试验区的2008年气象数据示意图；

图5为本发明实施例计算得到的地下水及盐分浓度与观测数据的对比结果示意图；

图6为本发明实施例计算得到的土壤水及盐分浓度与观测数据的对比结果示意图；

图7为本发明实施例双重迭代耦合计算过程中水分模块与溶质模块的迭代次数示意图。

具体实施方式

以下结合附图对本发明涉及的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置进行详细地说明。

<实施例>

本实施例中，以中国内蒙古河套灌区永联试验区区域尺度水盐运移过程的计算为例进行说明。本实施例目的为根据研究区的水文地质信息及2008年的气象数据，模拟全区域2008年的饱和-非饱和带的土壤水盐运移过程。如图1和2所示，本实施例所提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法包括以下步骤：

步骤1.研究区资料的收集与输入条件的设置：

研究区土地利用类型共分为56个分区，每个分区均计算其单独的垂向一维非饱和带水盐运动过程，饱和带在水平向划分为300×100个有限差分网格，计算其三维的水盐运移过程。计算时间为从2008年5月1日起至2008年10月31日止。本实施例的实测研究区域位置、实测土地利用类型分区、实测地表高程与实测监测点土壤质地情况见图3所示，2008年该区域的气象数据见图4所示。

研究区的地表高程、底板高程、河渠位置及水位等空间地理信息根据实际观测值设置。地下水水位、地下水浓度等研究对象的初始条件，根据区域10口观测井的实测数据插值得到区域数据。地下水饱和渗透系数、土壤水水力参数、溶质纵向弥散度等计算参数采用率定验证结果得到，分别为上层含水层的饱和渗透系数为1.06m/d，下层含水层的饱和渗透系数为3.5m/d，土壤水饱和渗透系数与土质有关，处于0.06到1.06m/d之间，饱和带地下水盐分纵向弥散度为10m，非饱和带土壤盐分的纵向弥散度为0.2m。之后采用步骤2至11开始第各个应力期各个时间步的具体计算工作，对于每个时间步均按照步骤2至11进行计算。

步骤2.设置非饱和带厚度：

假设当前正处于第t个计算时间步的第p个迭代步，应力期长度为ΔT，上一时间步计算得到的非饱和带厚度为

其中，

步骤3.运行非饱和带水分运移模块：

运行非饱和带水分运移模块UBMOD，计算时间段为从t时刻至t+ΔT时刻。

步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求：

计算非饱和剖面的所有时刻的Courant数(记为Cr)，假如Cr数超过1，则折减时间步长并返回步骤3重新计算非饱和带水分运移模型，直至满足Cr数判别条件(Cr≤1)。该步骤保证了非饱和带水分运动模块与非饱和带溶质运移模块计算时间步长的一致性。

步骤5.计算ΔT时段内的地下水平均补给速率R，以区域l为例：

式中，r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量。同该方法类似，可以得到该应力期中整个区域的地下水补给速率，将之写为与地下水有限差分网格一一对应的数组形式，记为R

步骤6.运行饱和带地下水运动模块；

将地下水补给速率R

步骤7.进行水分运动模块收敛性判断；

根据以下公式进行收敛性判断：

max(|H

式中，ε

步骤8.设置非饱和带溶质下边界条件：

假设当前正处于第t个计算时间步的第q个迭代步，上一时间步计算的得到的饱和带潜水含水层地下水浓度为

其中，

步骤9.运行非饱和带溶质运移模块：

运行非饱和带溶质运移模块，计算时间段为从t时刻至t+ΔT时刻，得到地下水补给量的浓度

步骤10.运行饱和带溶质运移模块：

将地下水补给浓度

步骤11.进行溶质运移模块收敛性判断：

式中，ε

根据以上方法计算得到全区域2008年的饱和-非饱和带的土壤水盐运移过程，将该计算结果与实测值进行比较：

图5展示了饱和带不同观测井地下水埋深、地下水浓度的实测值与模拟值对比结果，图6展示了按照土地利用类型划分的土壤监测点非饱和带含水率、土壤水浓度的实测值与模拟值对比结果。以平均误差(MRE)、均方根误差(RMSE)、纳什效率系数(NSE)与决定系数(R

根据以上计算结果可知，本方法对区域尺度饱和-非饱和水分运动及溶质运移的拟合精度良好，计算结果可以准确反应区域尺度水分及溶质的运移转化特征。

本实施例中，水分耦合模块与溶质耦合模块的迭代时间步长如图7所示，水分耦合模块的平均迭代次数仅为1.9次，溶质耦合模块的平均迭代次数仅为2.1次，本实施例的计算时长为524s。

以上数据证实了本方法具有较好的迭代效率、模型可靠性与计算稳定性。

进一步，本实施例还提供能够自动实现上述方法的区域饱和-非饱和水分及溶质运移计算装置，该装置包括控制离散部、参数信息确定部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部、输入显示部以及控制部。

离散部将研究区中非饱和带离散为L个分区，并将饱和带离散为m×n个有限差分网格，其中L、m、n均为正整数。

参数信息确定部，根据研究区的水文地质气象数据确定离散后各区域网格的参数信息。

非饱和带厚度设置部用于设置非饱和带厚度：设当前正处于第t个计算时间步的第p个迭代步，应力期长度为ΔT，上一时间步计算的得到的非饱和带厚度为

式中，S

非饱和带水分运移模型计算部运行非饱和带水分运移模型，计算时间段为从t时刻至t+ΔT时刻非饱和带水分运移。

非饱和带水盐运移一致性判断部对非饱和带水分运移模型计算部的计算结果是否满足溶质运移计算的需求进行判断：计算非饱和剖面的所有时刻的Courant数，记为Cr，若Cr＞1，则折减时间步长，并使非饱和带水分运移计算部重新计算非饱和带水分运移模型，直至满足Cr≤1。

地下水平均补给速率计算部计算ΔT时段内的地下水平均补给速率R：

式中，r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量；依次对各区域中地下水补给量进行计算，得到该应力期中整个研究区的地下水补给速率，将之写为与地下水有限差分网格一一对应的m×n维数组形式R

饱和带地下水运动模型计算部将地下水补给速率R

水分运动模型收敛性判断部采用下式对收敛情况进行判断：

max(|H

式中，ε

非饱和带溶质下边界设置部用于设置非饱和带溶质下边界条件：设当前正处于第t个计算时间步的第q个迭代步，上一时间步计算得到的饱和带潜水含水层地下水浓度为

式中，S

非饱和带溶质运移模型计算部运行非饱和带溶质运移模型，计算时间段为从t时刻至t+ΔT时刻，得到地下水补给量的浓度

饱和带溶质运移模型计算部将地下水补给浓度

溶质运移模型收敛判断部采用下式对收敛情况进行判断：

式中，ε

图像生成部根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据，生成反映饱和-非饱和水分运移情况的图像和反映溶质运移情况的图像。

动态信息生成部根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据，生成反映饱和-非饱和水分运移动态过程的动态信息和反映溶质运移动态过程的动态信息。

输入显示部用于让操作员输入操作指令，并根据操作指令对相应部的数据进行显示，例如，根据操作指令对图像生成部生成的反映饱和-非饱和水分运移情况的图像和反映溶质运移情况的图像进行显示，根据操作指令对动态信息生成部生成的反映饱和-非饱和水分运移动态过程的动态信息和反映溶质运移动态过程的动态信息进行显示。

控制部用于控制离散部、信息输入部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部的运行。

综上，本实施例所提供的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置，将土壤水分与地下水分、土壤溶质与地下水溶质通过双重迭代的方法予以耦合，从而可以精确获得区域尺度饱和-非饱和水分及溶质运移过程的运动状态与饱和带、非饱和带交互通量。该双重迭代耦合方法保证了计算精度、迭代效率与模型鲁棒。此外，本发明中为了保证状态变化较为剧烈的非饱和带水分及溶质运移过程计算的稳定性、精度与可靠性，因而通过Courant数将非饱和带水分运动模块与非饱和带溶质运移模块的计算时间步长严格保持一致。本发明为多过程数值求解耦合方法的研究提供了新的思路，为区域水分及溶质运移过程的研究提供了新方法与新工具。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置，包括：通过Courant数判断控制非饱和带水分运移模块计算时间步长、非饱和带水分运动模块与饱和带水分运动模块的迭代耦合、非饱和带溶质运移模块与饱和带溶质运移模块的迭代耦合、非饱和带水分运动模块与非饱和带溶质运移模块的计算时间步耦合等技术内容，并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。