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

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

Info

Publication number
CN112650970B
CN112650970B CN202011512457.2A CN202011512457A CN112650970B CN 112650970 B CN112650970 B CN 112650970B CN 202011512457 A CN202011512457 A CN 202011512457A CN 112650970 B CN112650970 B CN 112650970B
Authority
CN
China
Prior art keywords
unsaturated
solute
zone
saturated
unsaturated zone
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011512457.2A
Other languages
English (en)
Other versions
CN112650970A (zh
Inventor
毛威
朱焱
杨洋
赵天兴
杨金忠
伍靖伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202011512457.2A priority Critical patent/CN112650970B/zh
Publication of CN112650970A publication Critical patent/CN112650970A/zh
Application granted granted Critical
Publication of CN112650970B publication Critical patent/CN112650970B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

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

Description

饱和-非饱和水分及溶质运移双重迭代耦合方法及装置
技术领域
本发明属于地下水分及溶质计算领域,具体涉及区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法及装置。
技术背景
区域尺度饱和-非饱和水分及溶质运移规律一直是地下水资源与环境的研究重点,关系到地下水污染、土壤盐碱化、地下水开采、农业灌区灌溉排水等问题,建立区域尺度饱和-非饱和水分及溶质运移模型是一种有效的研究手段。非饱和带土壤水分是地下水系统的重要补给来源和潜层地下水的主要消耗途径,同时也是地下水污染的主要源头与缓冲区,因此需要将二者进行统一的分析与模拟。
饱和带地下水及溶质与非饱和带地下水及溶质的时空运移尺度存在较大差异,因此研究中常采用拟三维方法予以处理。基于拟三维方法的区域尺度饱和-非饱和水分及溶质模型目前仍然较少,仅有UZF-MODFLOW-MT3D、HYDRUS-MODFLOW-MT3D、Q3D等几种。UZF-MODFLOW-MT3D模型与HYDRUS-MODFLOW-MT3D模型基于松散耦合方法。该方法较为简单且易于实现,然而存在难以避免的质量误差问题,且不能保证耦合计算的精度。Q3D模型基于完全耦合方法,该方法虽然可以保证耦合计算的精度,但是实现方法极为复杂,且稳定性与收敛性较差。目前这些方法都难以得到准确有效的计算结果,从而无法获得准确、可靠的区域尺度饱和-非饱和水分及溶质运移信息。
发明内容
本发明是为了解决上述问题而进行的,目的在于提供一种基于饱和-非饱和水分运动系统的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置,能够有效提高区域尺度饱和-非饱和水分及溶质的计算精度、计算效率与模型稳定性,得到准确、可靠的区域尺度饱和-非饱和水分及溶质运移信息。
本发明为了实现上述目的,采用了以下方案:
<方法>
本发明提供了一种区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于,包括以下步骤:
步骤1.收集研究区的水文地质气象数据资料;将研究区中非饱和带离散为L个分区,饱和带离散为m×n个有限差分网格;根据研究区数据资料确定离散后各区域网格的参数信息;将这些参数信息作为输入参数,对于各个应力期各个时间步(从第一个应力期第一个时间步开始按顺序计算),均采用以下步骤2~11计算得到对应的饱和-非饱和水分及溶质运移情况;
步骤2.设置非饱和带厚度:
设当前正处于第t个计算时间步的第p个迭代步,应力期长度为ΔT,上一时间步计算得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带厚度能够得到整个研究区非饱和带的厚度为L维数组;
步骤3.运行非饱和带水分运移模块:
运行非饱和带水分运移模块,计算时间段为从t时刻至t+ΔT时刻;
步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求:
计算非饱和剖面的所有时刻的Courant数,记为Cr,若Cr>1,则折减时间步长并返回步骤3重新计算非饱和带水分运移模型,直至满足Cr≤1;
步骤5.计算ΔT时段内的地下水平均补给速率R:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量;依次对各区域中地下水补给量进行计算,得到该应力期中整个研究区的地下水补给速率,将之写为与地下水有限差分网格一一对应的m×n维数组形式,记为Rt,p
步骤6.运行饱和带地下水运动模块:
将地下水补给速率Rt,p传递给地下水运动模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深均为m×n维数组;
步骤7.进行水分运动模块收敛性判断:
max(|Ht,p-Ht,p-1|)<εw,
式中,εw为容忍误差;若两个公式同时满足,则判断为收敛,进入步骤8,进行溶质运移模块的计算;否则返回至步骤2开始该时间步的下一个迭代过程计算非饱和带水分模块,直至满足该收敛条件;
步骤8.设置非饱和带溶质下边界条件:
设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件为:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带溶质下边界条件能够得到整个研究区非饱和带溶质的下边界条件为L维数组;
步骤9.运行非饱和带溶质运移模块:
运行非饱和带溶质运移模块,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
步骤10.运行饱和带溶质运移模块:
将地下水补给浓度传递给溶质运动模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/>为m×n维数组;
步骤11.进行溶质运移模块收敛性判断:
式中,εs为容忍误差;若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至步骤8开始该时间步的下一个迭代过程计算非饱和带溶质模块,直至满足该收敛条件。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,还可以具有以下特征:在步骤1中,研究区的水文地质气象数据资料包括:水文地质数据、气象数据、土地利用类型数据、河渠水位数据等资料;根据研究区数据资料确定的参数信息为:包含地表高程、底板高程、河渠位置及水位空间的地理信息,和地下水位、地下水浓度的初始条件信息,以及地下水饱和渗透系数、土壤水水力参数、溶质纵向弥散度等信息。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,还可以具有以下特征:在步骤2中,对于第一个应力期的第一个时间步,采用根据步骤1获取的输入参数所确定的初始条件进行计算;若为第一个迭代步,则采用上一个时间步计算得到的结果作为非饱和带厚度。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,还可以具有以下特征:在步骤6中,地下水运动模型为MODFLOW模型。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,还可以具有以下特征:在步骤8中,于第一个应力期的第一个时间步,采用根据步骤1获取的输入参数所确定的初始条件进行计算;若为第一个迭代步,则采用上一个时间步计算得到的结果作为饱和带潜水含水层地下水浓度。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,还可以具有以下特征:在步骤10中,溶质运动模型为MT3DMS模型。
<装置>
进一步,本发明还提供了一种自动实现上述<方法>的区域饱和-非饱和水分及溶质运移计算装置,其特征在于,包括:
离散部,将研究区中非饱和带离散为L个分区,并将饱和带离散为m×n个有限差分网格,其中L、m、n均为正整数;
参数信息确定部,根据研究区的水文地质气象数据确定离散后各区域网格的参数信息;
非饱和带厚度设置部,设当前正处于第t个计算时间步的第p个迭代步,应力期长度为ΔT,上一时间步计算得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带厚度得到整个研究区非饱和带的厚度为L维数组;
非饱和带水分运移模型计算部,运行非饱和带水分运移模型,计算时间段为从t时刻至t+ΔT时刻非饱和带水分运移;
非饱和带水盐运移一致性判断部,对非饱和带水分运移模型计算部的计算结果是否满足溶质运移计算的需求进行判断:计算非饱和剖面的所有时刻的Courant数,记为Cr,若Cr>1,则折减时间步长,并使非饱和带水分运移计算部重新计算非饱和带水分运移模型,直至满足Cr≤1;
地下水平均补给速率计算部,计算ΔT时段内的地下水平均补给速率R:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量;依次对各区域中地下水补给量进行计算,得到该应力期中整个研究区的地下水补给速率,将之写为与地下水有限差分网格一一对应的m×n维数组形式Rt,p
饱和带地下水运动模型计算部,将地下水补给速率Rt,p传递给地下水运动模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深均为m×n维数组;
水分运动模型收敛性判断部,采用下式对收敛情况进行判断:
max(|Ht,p-Ht,p-1|)<εw,
式中,εw为容忍误差;若两个公式同时满足,则判断为收敛,进入非饱和带溶质下边界设置部,进行溶质运移模型的计算;否则返回至非饱和带水盐运移一致性判断部开始该时间步的下一个迭代过程计算非饱和带水分运动,直至满足该收敛条件;
非饱和带溶质下边界设置部,设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为/>第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件为:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带溶质下边界条件得到整个研究区非饱和带溶质的下边界条件为L维数组;
非饱和带溶质运移模型计算部,运行非饱和带溶质运移模型,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
饱和带溶质运移模型计算部,将地下水补给浓度传递给溶质运动模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/>为m×n维数组;
溶质运移模型收敛判断部,采用下式对收敛情况进行判断:
式中,εs为容忍误差;若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至非饱和带溶质下边界设置部开始该时间步的下一个迭代过程计算非饱和带溶质模型,直至满足该收敛条件。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移计算装置,还可以包括:图像生成部,根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的结果数据,生成反映饱和-非饱和水分运移情况的图像,和反映溶质运移情况的图像。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移计算装置,还可以包括:动态信息生成部,根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的结果数据,生成反映饱和-非饱和水分运移动态过程的动态信息,和反映溶质运移动态过程的动态信息。
优选地,本发明提供的区域饱和-非饱和水分及溶质运移计算装置,还可以包括:控制部,控制离散部、参数信息确定部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部的运行。
发明的作用与效果
与现有技术相比,本发明的有益效果在于:
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,上一时间步计算得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>如果为第一个应力期的第一个时间步,则采用初始条件进行计算,如果为第一个迭代步,则采用上一个时间步计算得到的最终结果。以面积为Sl的第l个子区域为例,非饱和带厚度为,
其中,指在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度。对于整个研究区,其非饱和带厚度为/>(56维数组)。
步骤3.运行非饱和带水分运移模块:
运行非饱和带水分运移模块UBMOD,计算时间段为从t时刻至t+ΔT时刻。
步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求:
计算非饱和剖面的所有时刻的Courant数(记为Cr),假如Cr数超过1,则折减时间步长并返回步骤3重新计算非饱和带水分运移模型,直至满足Cr数判别条件(Cr≤1)。该步骤保证了非饱和带水分运动模块与非饱和带溶质运移模块计算时间步长的一致性。
步骤5.计算ΔT时段内的地下水平均补给速率R,以区域l为例:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量。同该方法类似,可以得到该应力期中整个区域的地下水补给速率,将之写为与地下水有限差分网格一一对应的数组形式,记为Rt,p(300×100维数组)。
步骤6.运行饱和带地下水运动模块;
将地下水补给速率Rt,p传递给MODFLOW模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深(均为300×100维数组)。
步骤7.进行水分运动模块收敛性判断;
根据以下公式进行收敛性判断:
max(|Ht,p-Ht,p-1|)<εw, (3)
式中,εw为容忍误差,取为0.05m。若两个公式同时满足,则判断为收敛,进入步骤8溶质运移模块的计算;否则返回至步骤2开始该时间步的下一个迭代过程计算非饱和带水分模块,直至满足该收敛条件。
步骤8.设置非饱和带溶质下边界条件:
假设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算的得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为/>如果为第一个应力期的第一个时间步,则采用初始条件进行计算,如果为第一个迭代步,则采用上一个时间步计算得到的最终结果。以面积为Sl的第l个子区域为例,非饱和带溶质下边界条件为:
其中,指在第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件。对于整个研究区,其非饱和带溶质下边界条件为/>(56维数组)。
步骤9.运行非饱和带溶质运移模块:
运行非饱和带溶质运移模块,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
步骤10.运行饱和带溶质运移模块:
将地下水补给浓度传递给MT3DMS模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/>(300×100维数组)。
步骤11.进行溶质运移模块收敛性判断:
式中,εs为容忍误差,取为0.05kg/m3。若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至步骤8开始该时间步的下一个迭代过程计算非饱和带溶质模块,直至满足该收敛条件。
根据以上方法计算得到全区域2008年的饱和-非饱和带的土壤水盐运移过程,将该计算结果与实测值进行比较:
图5展示了饱和带不同观测井地下水埋深、地下水浓度的实测值与模拟值对比结果,图6展示了按照土地利用类型划分的土壤监测点非饱和带含水率、土壤水浓度的实测值与模拟值对比结果。以平均误差(MRE)、均方根误差(RMSE)、纳什效率系数(NSE)与决定系数(R2)为评价指标。图5中区域平均地下水位的MRE为24.8%,RMSE为0.183m,NSE为0.760,R2为0.832,地下水浓度的MRE为8.8%,RMSE为0.265g/L,NSE为0.011,R2为0.512。对于不同观测井,其地下水位的MRE小于50.0%,RMSE小于0.40m,NSE大于0.40,R2大于0.60,地下水浓度的MRE为小于33.0%,RMSE小于0.7g/L,NSE大于0.00,R2大于0.01。图6中耕地不同日期的土壤水分含量MRE小于25%,RMSE小于0.07cm3/cm3,R2大于0.50,土壤盐分含量MRE小于20%,RMSE小于0.07g/L,R2大于0.20。城镇不同日期的土壤水分含量MRE小于33%,RMSE小于0.08cm3/cm3,R2大于0.50,土壤盐分含量MRE小于40%,RMSE小于0.4g/L,R2大于0.90。荒地不同日期的土壤水分含量MRE小于35%,RMSE小于0.08cm3/cm3,R2大于0.60,土壤盐分含量MRE小于50%,RMSE小于0.2g/L,R2大于0.80。
根据以上计算结果可知,本方法对区域尺度饱和-非饱和水分运动及溶质运移的拟合精度良好,计算结果可以准确反应区域尺度水分及溶质的运移转化特征。
本实施例中,水分耦合模块与溶质耦合模块的迭代时间步长如图7所示,水分耦合模块的平均迭代次数仅为1.9次,溶质耦合模块的平均迭代次数仅为2.1次,本实施例的计算时长为524s。
以上数据证实了本方法具有较好的迭代效率、模型可靠性与计算稳定性。
进一步,本实施例还提供能够自动实现上述方法的区域饱和-非饱和水分及溶质运移计算装置,该装置包括控制离散部、参数信息确定部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部、输入显示部以及控制部。
离散部将研究区中非饱和带离散为L个分区,并将饱和带离散为m×n个有限差分网格,其中L、m、n均为正整数。
参数信息确定部,根据研究区的水文地质气象数据确定离散后各区域网格的参数信息。
非饱和带厚度设置部用于设置非饱和带厚度:设当前正处于第t个计算时间步的第p个迭代步,应力期长度为ΔT,上一时间步计算的得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带厚度得到整个研究区非饱和带的厚度为L维数组。
非饱和带水分运移模型计算部运行非饱和带水分运移模型,计算时间段为从t时刻至t+ΔT时刻非饱和带水分运移。
非饱和带水盐运移一致性判断部对非饱和带水分运移模型计算部的计算结果是否满足溶质运移计算的需求进行判断:计算非饱和剖面的所有时刻的Courant数,记为Cr,若Cr>1,则折减时间步长,并使非饱和带水分运移计算部重新计算非饱和带水分运移模型,直至满足Cr≤1。
地下水平均补给速率计算部计算ΔT时段内的地下水平均补给速率R:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量;依次对各区域中地下水补给量进行计算,得到该应力期中整个研究区的地下水补给速率,将之写为与地下水有限差分网格一一对应的m×n维数组形式Rt,p
饱和带地下水运动模型计算部将地下水补给速率Rt,p传递给地下水运动模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深均为m×n维数组;
水分运动模型收敛性判断部采用下式对收敛情况进行判断:
max(|Ht,p-Ht,p-1|)<εw,
式中,εw为容忍误差;若两个公式同时满足,则判断为收敛,进入非饱和带溶质下边界设置部,进行溶质运移模型的计算;否则返回至非饱和带水分运移模型计算部开始该时间步的下一个迭代过程计算非饱和带水分部,直至满足该收敛条件;
非饱和带溶质下边界设置部用于设置非饱和带溶质下边界条件:设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为/>第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件为:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带溶质下边界条件得到整个研究区非饱和带溶质的下边界条件为L维数组。
非饱和带溶质运移模型计算部运行非饱和带溶质运移模型,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
饱和带溶质运移模型计算部将地下水补给浓度传递给溶质运动模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/>为m×n维数组。
溶质运移模型收敛判断部采用下式对收敛情况进行判断:
式中,εs为容忍误差;若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至非饱和带溶质下边界设置部开始该时间步的下一个迭代过程计算非饱和带溶质模型,直至满足该收敛条件。
图像生成部根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据,生成反映饱和-非饱和水分运移情况的图像和反映溶质运移情况的图像。
动态信息生成部根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据,生成反映饱和-非饱和水分运移动态过程的动态信息和反映溶质运移动态过程的动态信息。
输入显示部用于让操作员输入操作指令,并根据操作指令对相应部的数据进行显示,例如,根据操作指令对图像生成部生成的反映饱和-非饱和水分运移情况的图像和反映溶质运移情况的图像进行显示,根据操作指令对动态信息生成部生成的反映饱和-非饱和水分运移动态过程的动态信息和反映溶质运移动态过程的动态信息进行显示。
控制部用于控制离散部、信息输入部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部的运行。
综上,本实施例所提供的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置,将土壤水分与地下水分、土壤溶质与地下水溶质通过双重迭代的方法予以耦合,从而可以精确获得区域尺度饱和-非饱和水分及溶质运移过程的运动状态与饱和带、非饱和带交互通量。该双重迭代耦合方法保证了计算精度、迭代效率与模型鲁棒。此外,本发明中为了保证状态变化较为剧烈的非饱和带水分及溶质运移过程计算的稳定性、精度与可靠性,因而通过Courant数将非饱和带水分运动模块与非饱和带溶质运移模块的计算时间步长严格保持一致。本发明为多过程数值求解耦合方法的研究提供了新的思路,为区域水分及溶质运移过程的研究提供了新方法与新工具。
以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的饱和-非饱和水分及溶质运移双重迭代耦合方法及装置,包括:通过Courant数判断控制非饱和带水分运移模块计算时间步长、非饱和带水分运动模块与饱和带水分运动模块的迭代耦合、非饱和带溶质运移模块与饱和带溶质运移模块的迭代耦合、非饱和带水分运动模块与非饱和带溶质运移模块的计算时间步耦合等技术内容,并不仅仅限定于在以上实施例中所描述的内容,而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换,都在本发明的权利要求所要求保护的范围内。

Claims (10)

1.区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于,包括以下步骤:
步骤1.收集研究区的水文地质气象数据资料;将研究区中非饱和带离散为L个分区,饱和带离散为m×n个有限差分网格;根据研究区数据资料确定离散后各区域网格的参数信息;将这些参数信息作为输入参数,采用以下步骤2至11计算得到各应力期各时间步下的饱和-非饱和水分及溶质运移情况;
步骤2.设置非饱和带厚度:
设当前正处于第t个计算时间步的第p个迭代步,应力期长度为ΔT,上一时间步计算得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带厚度能够得到整个研究区非饱和带的厚度 为L维数组;
步骤3.运行非饱和带水分运移模块:
运行非饱和带水分运移模块,计算时间段为从t时刻至t+ΔT时刻;
步骤4.判断非饱和带水分运动模块的计算结果是否满足溶质运移计算的需求:
计算非饱和剖面的所有时刻的Courant数,记为Cr,若Cr>1,则折减时间步长并返回步骤3重新计算非饱和带水分运移模型,直至满足Cr≤1;
步骤5.计算ΔT时段内的地下水平均补给速率R:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量;依次对各区域中地下水补给量进行计算,得到该应力期中整个研究区的地下水补给速率,将之写为与地下水有限差分网格一一对应的m×n维数组形式,记为Rt,p
步骤6.运行饱和带地下水运动模块:
将地下水补给速率Rt,p传递给地下水运动模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深均为m×n维数组;
步骤7.进行水分运动模块收敛性判断:
max(|Ht,p-Ht,p-1|)<εw,
式中,εw为容忍误差;若两个公式同时满足,则判断为收敛,进入步骤8,进行溶质运移模块的计算;否则返回至步骤2开始该时间步的下一个迭代过程计算非饱和带水分模块,直至满足该收敛条件;
步骤8.设置非饱和带溶质下边界条件:
设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为/>第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件为:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带溶质下边界条件能够得到整个研究区非饱和带溶质的下边界条件 为L维数组;
步骤9.运行非饱和带溶质运移模块:
运行非饱和带溶质运移模块,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
步骤10.运行饱和带溶质运移模块:
将地下水补给浓度传递给溶质运动模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/> 为m×n维数组;
步骤11.进行溶质运移模块收敛性判断:
式中,εs为容忍误差;若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至步骤8开始该时间步的下一个迭代过程计算非饱和带溶质模块,直至满足该收敛条件。
2.根据权利要求1所述的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于:
其中,在步骤1中,研究区的水文地质气象数据资料包括:水文地质数据、气象数据、土地利用类型数据、河渠水位数据资料;根据研究区数据资料确定的参数信息为:包含地表高程、底板高程、河渠位置及水位空间的地理信息,和地下水位、地下水浓度的初始条件信息,以及地下水饱和渗透系数、土壤水水力参数、溶质纵向弥散度信息。
3.根据权利要求1所述的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于:
其中,在步骤2中,对于第一个应力期的第一个时间步,采用根据步骤1获取的输入参数所确定的初始条件进行计算;若为第一个迭代步,则采用上一个时间步计算得到的结果作为非饱和带厚度。
4.根据权利要求1所述的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于:
其中,在步骤6中,地下水运动模型为MODFLOW模型。
5.根据权利要求1所述的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于:
其中,在步骤8中,于第一个应力期的第一个时间步,采用根据步骤1获取的输入参数所确定的初始条件进行计算;若为第一个迭代步,则采用上一个时间步计算得到的结果作为饱和带潜水含水层地下水浓度。
6.根据权利要求1所述的区域饱和-非饱和水分及溶质运移模型双重迭代耦合方法,其特征在于:
其中,在步骤10中,溶质运动模型为MT3DMS模型。
7.区域饱和-非饱和水分及溶质运移计算装置,其特征在于,包括:
离散部,将研究区中非饱和带离散为L个分区,并将饱和带离散为m×n个有限差分网格,其中L、m、n均为正整数;
信息输入部,让操作员输入离散后各区域网格的参数信息;
非饱和带厚度设置部,设当前正处于第t个计算时间步的第p个迭代步,应力期长度为ΔT,上一时间步计算得到的非饱和带厚度为上一个迭代步计算得到的非饱和带厚度为/>在第t个计算时间步第p次迭代中,第l个子区域的非饱和带厚度:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带厚度得到整个研究区非饱和带的厚度 为L维数组;
非饱和带水分运移模型计算部,运行非饱和带水分运移模型,计算时间段为从t时刻至t+ΔT时刻非饱和带水分运移;
非饱和带水盐运移一致性判断部,对非饱和带水分运移模型计算部的计算结果是否满足溶质运移计算的需求进行判断:计算非饱和剖面的所有时刻的Courant数,记为Cr,若Cr>1,则折减时间步长,并使非饱和带水分运移计算部重新计算非饱和带水分运移模型,直至满足Cr≤1;
地下水平均补给速率计算部,计算ΔT时段内的地下水平均补给速率R:
式中,r代表区域l中非饱和带水分运移模型每个时刻计算得到的地下水补给量;依次对各区域中地下水补给量进行计算,得到该应力期中整个研究区的地下水补给速率,将之写为与地下水有限差分网格一一对应的m×n维数组形式Rt,p
饱和带地下水运动模型计算部,将地下水补给速率Rt,p传递给地下水运动模型,并在t时刻至t+ΔT时刻运行该模型,获得地下水位Ht,p和地下水埋深均为m×n维数组;
水分运动模型收敛性判断部,采用下式对收敛情况进行判断:
max(|Ht,p-Ht,p-1|)<εw,
式中,εw为容忍误差;若两个公式同时满足,则判断为收敛,进入非饱和带溶质下边界设置部,进行溶质运移模型的计算;否则返回至非饱和带厚度设置部开始该时间步的下一个迭代过程计算非饱和带水分运动,直至满足该收敛条件;
非饱和带溶质下边界设置部,设当前正处于第t个计算时间步的第q个迭代步,上一时间步计算得到的饱和带潜水含水层地下水浓度为上一个迭代步计算得到的饱和带潜水含水层地下水浓度为/>第t个计算时间步第q次迭代中,第l个子区域的非饱和带溶质下边界条件为:
式中,Sl为第l个子区域非饱和带的面积;根据各子区域的非饱和带溶质下边界条件得到整个研究区非饱和带溶质的下边界条件 为L维数组;
非饱和带溶质运移模型计算部,运行非饱和带溶质运移模型,计算时间段为从t时刻至t+ΔT时刻,得到地下水补给量的浓度
饱和带溶质运移模型计算部,将地下水补给浓度传递给溶质运动模型,并在t时刻至t+ΔT时刻运行该模型,获得饱和带潜水含水层地下水浓度为/> 为m×n维数组;
溶质运移模型收敛判断部,采用下式对收敛情况进行判断:
式中,εs为容忍误差;若两个公式同时满足,则判断为收敛,进入下一个时间步的计算;否则返回至非饱和带溶质下边界设置部开始该时间步的下一个迭代过程计算非饱和带溶质模型,直至满足该收敛条件。
8.根据权利要求7所述的区域饱和-非饱和水分及溶质运移计算装置,其特征在于,还包括:
图像生成部,根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据,生成反映饱和-非饱和水分运移情况的图像,和反映溶质运移情况的图像。
9.根据权利要求8所述的区域饱和-非饱和水分及溶质运移计算装置,其特征在于,还包括:
动态信息生成部,根据非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部的数据,生成反映饱和-非饱和水分运移动态过程的动态信息,和反映溶质运移动态过程的动态信息。
10.根据权利要求9所述的区域饱和-非饱和水分及溶质运移计算装置,其特征在于,还包括:
控制部,控制离散部、信息输入部、非饱和带厚度设置部、非饱和带水分运移模型计算部、非饱和带水盐运移一致性判断部、地下水平均补给速率计算部、饱和带地下水运动模型计算部、水分运动模型收敛性判断部、非饱和带溶质下边界设置部、非饱和带溶质运移模型计算部、饱和带溶质运移模型计算部、溶质运移模型收敛判断部、图像生成部、动态信息生成部的运行。
CN202011512457.2A 2020-12-19 2020-12-19 饱和-非饱和水分及溶质运移双重迭代耦合方法及装置 Active CN112650970B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011512457.2A CN112650970B (zh) 2020-12-19 2020-12-19 饱和-非饱和水分及溶质运移双重迭代耦合方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011512457.2A CN112650970B (zh) 2020-12-19 2020-12-19 饱和-非饱和水分及溶质运移双重迭代耦合方法及装置

Publications (2)

Publication Number Publication Date
CN112650970A CN112650970A (zh) 2021-04-13
CN112650970B true CN112650970B (zh) 2024-03-22

Family

ID=75358710

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011512457.2A Active CN112650970B (zh) 2020-12-19 2020-12-19 饱和-非饱和水分及溶质运移双重迭代耦合方法及装置

Country Status (1)

Country Link
CN (1) CN112650970B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114444252B (zh) * 2021-11-25 2022-08-26 中国科学院生态环境研究中心 基于环境容量及自然消减模型的土壤环境承载力计算方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
CN105022913A (zh) * 2015-06-01 2015-11-04 中国水利水电科学研究院 一种降雨入渗补给地下水临界埋深计算方法
KR20160000261A (ko) * 2014-06-24 2016-01-04 한국원자력연구원 지하수 내 오염물 반응 이동 모사 시스템 및 모사 방법
CN106777968A (zh) * 2016-12-14 2017-05-31 中国水利水电科学研究院 一种大深埋条件下地表水补给地下水的计算方法及装置
CN109376433A (zh) * 2018-10-26 2019-02-22 北京市水文总站 基于土壤非饱和水和地下水耦合的区域水流运动模拟方法
CN109522655A (zh) * 2018-11-19 2019-03-26 武汉大学 基于变饱和水分运动系统的区域地下水补给量计算方法
WO2020206908A1 (zh) * 2019-04-12 2020-10-15 中国科学院南京土壤研究所 基于PaaS平台的超低功耗土壤近地无线传感系统及使用方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
KR20160000261A (ko) * 2014-06-24 2016-01-04 한국원자력연구원 지하수 내 오염물 반응 이동 모사 시스템 및 모사 방법
CN105022913A (zh) * 2015-06-01 2015-11-04 中国水利水电科学研究院 一种降雨入渗补给地下水临界埋深计算方法
CN106777968A (zh) * 2016-12-14 2017-05-31 中国水利水电科学研究院 一种大深埋条件下地表水补给地下水的计算方法及装置
CN109376433A (zh) * 2018-10-26 2019-02-22 北京市水文总站 基于土壤非饱和水和地下水耦合的区域水流运动模拟方法
CN109522655A (zh) * 2018-11-19 2019-03-26 武汉大学 基于变饱和水分运动系统的区域地下水补给量计算方法
WO2020206908A1 (zh) * 2019-04-12 2020-10-15 中国科学院南京土壤研究所 基于PaaS平台的超低功耗土壤近地无线传感系统及使用方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
井渠结合区饱和-非饱和带水均衡模型研究;潘敏;伍靖伟;于健;杨金忠;;中国农村水利水电;20151015(第10期);全文 *

Also Published As

Publication number Publication date
CN112650970A (zh) 2021-04-13

Similar Documents

Publication Publication Date Title
Xu et al. Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas
US20120253673A1 (en) Method and apparatus for three dimensional dynamic measurements in water system
Cuthbert et al. Linking soil moisture balance and source-responsive models to estimate diffuse and preferential components of groundwater recharge
Howes et al. Modeling runoff and runon in a desert shrubland ecosystem, Jornada Basin, New Mexico
CN111766189B (zh) 一种基于水力刺激的堤防隐伏渗漏通道三维层析扫描方法
CN114861502A (zh) 基于Modflow模型的三维动态地下水污染模拟的安全饮用水区域确定方法
CN115238550A (zh) 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法
Banejad et al. Numerical simulation of groundwater flow and contamination transport in Nahavand plain aquifer, west of Iran
CN112650970B (zh) 饱和-非饱和水分及溶质运移双重迭代耦合方法及装置
Baroncini‐Turricchia et al. Integrating MRS data with hydrologic model‐Carrizal Catchment (Spain)
Rossi et al. Depression storage and infiltration effects on overland flow depth-velocity-friction at desert conditions: field plot results and model
Alaboz et al. Assessment of various pedotransfer functions for the prediction of the dry bulk density of cultivated soils in a semiarid environment
Motarjemi et al. Important factors when simulating the water and nitrogen balance in a tile-drained agricultural field under long-term monitoring
CN113484210B (zh) 一种强风化层弥散度现场尺度试验测定方法
Mazi et al. A groundwater-based, objective-heuristic parameter optimisation method for a precipitation-runoff model and its application to a semi-arid basin
CN115983158B (zh) 一种地下水模型与二维水动力模型松散耦合的方法
CN111898257A (zh) 区域暗管布局及排水排盐数值模拟方法和装置
Jhorar Estimation of effective soil hydraulic parameters for water management studies in semi-arid zones: integral use of modelling, remote sensing and parameter estimation
CN109724570B (zh) 地下跌水的跌水量、跌水宽度、坎上水层厚度的计算方法
Kisekka et al. Simulating water table response to proposed changes in surface water management in the C-111 agricultural basin of south Florida
Li et al. Study on pollutant model construction and three-dimensional spatial interpolation in soil environmental survey
Kahsay Groundwater resource assessment through distributed steady-state flow modeling, Aynalem Wellfield, Mekele, Ethiopia
Budhathoki Modelling snowmelt infiltration processes in seasonally frozen ground
Berhe Modeling groundwater resources of the eastern fringe of Botswana Kalahari: integrated numerical groundwater flow model
Liu et al. Time Series Data Inversion and Monitoring Method for Cross-Hole ERT Based on an Improved Extended Kalman Filter

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant