CN108763797A - 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 - Google Patents
一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 Download PDFInfo
- Publication number
- CN108763797A CN108763797A CN201810561011.5A CN201810561011A CN108763797A CN 108763797 A CN108763797 A CN 108763797A CN 201810561011 A CN201810561011 A CN 201810561011A CN 108763797 A CN108763797 A CN 108763797A
- Authority
- CN
- China
- Prior art keywords
- lake
- water
- grid cell
- level
- lakebed
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Aeration Devices For Treatment Of Activated Polluted Sludge (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于地下水模型的湖泊与地下水稳定流作用模拟方法,从而为湖泊和地下水稳定流相互作用的模拟提供了一种全新的技术手段。本发明克服了现有方法中地下水网格系统的剖分与湖泊刻画相互依存的问题,使地下水网格系统的剖分与湖泊刻画相互独立,一方面大大增加了方法适用性,另一方面也有效避免了模拟过程中由于垂向剖分层数较多导致的数值震荡问题。
Description
技术领域
本发明属于地下水数值模拟技术领域,具体涉及一种基于地下水模型的湖泊与地下水稳定流作用模拟方法的设计。
背景技术
湖泊是地表相对封闭可蓄水的天然或人工洼池,大都位于地下水的浅埋深带,与地下水特别是潜水具有显著的水力联系。湖泊-地下水相互作用耦合模拟在地下水数值模型中是一个重要研究方面,如当前广泛使用MODFLOW已针对性的开发出了LAK3软件包提供给用户使用。然而湖泊-地下水相互作用耦合模拟具有其难度和复杂性,一是湖泊水位受地下水的补给/排泄的影响是动态变化的,不能直接给定,需要以未知量的形式进行模拟试算;二是由于湖底的非规则性,湖泊积水区和未积水区面积随着湖泊水位的变化可以相互转化,模拟过程中需要进行动态调整。虽然LAK3针对上述问题进行了算法设计,但由于其算法机制的原因应用起来仍然具有一定的局限性。如LAK3不适用于起伏含水层系统的情况,地下水网格系统剖分时含水层的顶板和底板必须完全水平;湖泊的空间位置及其范围是由三维地下水数值网格系统中的无效单元格来描述的,当需要细化湖底高程分布时,将导致含水层的垂向剖分层数和网格单元数量的大幅增加,相应增加了计算成本。更重要的是由于MODFLOW在处理单元干-湿转化过程中采用经验处理方法,含水层垂向剖分太细时有时会由于单元干-湿转化计算的不稳定导致模拟失败。
发明内容
本发明的目的是提出一种基于地下水模型的湖泊与地下水稳定流作用模拟方法,解决现有技术中模拟湖泊和地下水的相互作用存在的如下问题:
(1)现有的湖泊地下水相互作用方法对湖泊水量平衡的刻画还不够全面;
(2)现有的湖泊地下水相互作用方法地下水网格系统的剖分与湖泊刻画相互依存,若要对湖泊进行相对精细的刻画,就需要对地下水网格系统也进行比较细致的剖分,不仅大大增加了计算量,同时容易出现含水层垂向剖分过细导致的地下水单元干-湿转换计算不稳定的问题;
(3)现有的湖泊地下水相互作用方法未考虑边界条件的连续性。
本发明的技术方案为:一种基于地下水模型的湖泊与地下水稳定流作用模拟方法,包括以下步骤:
S1、对湖泊进行离散化处理及边界连续性处理,得到各湖泊网格单元。
S2、采集并获取湖泊参数数据。
S3、根据上一次迭代下的湖泊水位判断各湖泊网格单元所处的积水状态。
S4、建立地下水数值计算矩阵方程。
S5、求解地下水数值计算矩阵方程,得到当前迭代下的地下水位。
S6、根据各湖泊网格单元所处的积水状态,结合当前迭代下的地下水位计算得到各湖泊网格单元的水量平衡项。
S7、根据各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量。
S8、根据湖泊与地下水的水分交换统计量和湖泊参数数据更新当前迭代下的湖泊水位。
S9、判断当前迭代下湖泊水位是否收敛,若是则模拟结束,否则返回步骤S3进入下一次迭代。
本发明的有益效果是:
(1)系统全面的进行湖泊水量平衡的刻画是提高湖泊地下水相互作用模拟精度的首要条件,本发明同时全面考虑了包括地下水作用在内的各项湖泊水量的源汇项,为提高模拟精度打下了坚实的基础。
(2)本发明使地下水网格系统的剖分与湖泊刻画相互独立,一方面大大增加了方法适用性,另一方面也有效避免了模拟过程中由于垂向剖分层数较多导致的数值震荡问题。
(3)边界条件的连续性一直是地下水数值模拟计算能否收敛的关键,本发明在湖底计算单元中,将其概化为倾斜的,根据湖泊水位和湖底高程之间的相对关系,将湖泊计算单元分为完全积水,部分积水和完全未积水三种状态模拟湖泊-地下水相互作用,不同状态之间的切换过程中能够保证边界条件的连续性,从而有利于计算过程的收敛。
附图说明
图1所示为一种基于地下水模型的湖泊与地下水稳定流作用模拟方法流程图。
图2所示为湖底数值离散示意图。
图3所示为湖泊网格单元内湖底高程变化示意图。
图4所示为完全积水单元地下水位高于单元内湖底高程最高值示意图。
图5所示为完全积水单元地下水位低于单元内湖底高程最低值示意图。
图6所示为完全积水单元地下水位低于单元内湖底高程最低、最高值之间示意图。
图7所示为完全未积水单元地下水位高于单元内湖底高程最高值示意图。
图8所示为完全未积水单元地下水位低于单元内湖底高程最低值示意图。
图9所示为完全未积水单元地下水位低于单元内湖底高程最低、最高值之间示意图。
图10所示为部分积水单元地下水位高于湖泊平均水位示意图。
图11所示为部分积水单元地下水位低于单元内湖底高程最低值示意图。
图12所示为部分积水单元地下水位位于单元内湖底高程最低值和湖泊平均水位之间示意图。
图13所示为部分积水单元地下水位高于单元内湖底高程最高值示意图。
图14所示为部分积水单元地下水位低于湖泊平均水位示意图。
图15所示为部分积水单元地下水位位于湖泊平均水位和单元内湖底高程最高值之间示意图。
图16所示为湖泊垂向示意图。
图17所示为湖泊水位-水面面积关系曲线示意图。
图18所示为模拟区平面图。
图19所示为本发明所述方法湖泊网格单元湖底高程刻画图。
图20所示为传统湖泊地下水相互作用方法研究区含水层剖面图。
图21所示为湖泊为半蓄水状态时两种方法湖泊水位模拟结果对比示意图。
图22所示为湖泊为蓄满状态时两种方法湖泊水位模拟结果对比示意图。
图23所示为湖泊为干涸状态时两种方法湖泊水位模拟结果对比示意图。
具体实施方式
现在将参考附图来详细描述本发明的示例性实施方式。应当理解,附图中示出和描述的实施方式仅仅是示例性的,意在阐释本发明的原理和精神,而并非限制本发明的范围。
本发明实施例提供了一种基于地下水模型的湖泊与地下水稳定流作用模拟方法,如图1所示,包括以下步骤S1-S9:
步骤S1具体包括:
S1-1、对湖泊进行离散化处理。
进行湖泊的空间离散化处理是湖泊-地下水作用耦合模拟的前提和基础,很大程度上决定了计算方法的选择,本发明实施例中处理过程为:将湖底所在的含水层网格单元标识为湖泊网格单元,并将湖底高程数据以离散的方式赋予每个湖泊网格单元,使得每个湖泊网格单元都具有其单元面积范围内的平均湖底高程值。在垂向方向上含水层可能剖分为多层,但只有湖底所在层位的网格单元为湖泊网格单元。如果某个湖泊网格单元不位于第一层,则该湖泊网格单元上方的所有网格单元都将被定义为无效单元,如图2所示。
S1-2、对湖泊进行边界连续性处理.
离散化处理过程中,各湖泊网格单元的湖底平均高程按照从下至上的顺序进行排序,各湖泊网格单元也被分为下级单元和上级单元。离散化时可能有若干个湖泊网格单元具有相同的湖底平均高程,它们属于同一级单元。以图3的湖底平均高程分布为示例,共有7级单元,其中最下级的湖泊网格单元的湖底平均高程为Lb,1,最上级湖泊网格单元的湖底平均高程为Lb,7。
本发明实施例中,假设湖底高程在湖泊网格单元内是倾斜的,具有单元内的最低值和最高值;对于非最上级和非最下级的湖泊网格单元,其湖底高程最低值为本单元湖底平均高程与下一级单元湖底平均高程的中间位置,最高值为本单元湖底平均高程与上一级单元湖底平均高程的中间位置;对于最下级的湖泊网格单元,其湖底高程最低值为其平均湖底高程减去其与上一级单元平均湖底高程差的一半;对于最上级湖泊网格单元,其湖底高程最高值为其平均湖底高程加上其与下一级单元平均湖底高程差的一半。
图3演示了不同湖泊网格单元内湖底高程的最高值和最低值的情况。例如,对于最下级单元,其单元内湖底高程的最低值为Lb,1-(Lb,2-Lb,1)/2,最高值为(Lb,2+Lb,1)/2。对于具有湖底平均高程Lb,2的湖泊网格单元,其单元内湖底高程的最低值为(Lb,2+Lb,1)/2,最高值为(Lb,2+Lb,3)/2。对于最上级单元,其单元内湖底高程的最低值为(Lb,6+Lb,7)/2,最高值为Lb,7+(Lb,7-Lb,6)/2。
根据湖泊水位的情况,湖泊网格单元被分为三种积水状态:(1)完全积水状态:湖泊水位高于其单元内湖底高程最高值时。以图3的湖泊水位为示例,具有湖底平均高程Lb,1,Lb,2,Lb,3,Lb,4,Lb,5的湖泊网格单元为完全积水状态。(2)完全未积水状态:湖泊水位低于其单元内湖底高程最低值时。以图3的湖泊水位为示例,具有湖底平均高程Lb,7的湖泊网格单元为完全未积水状态;(3)部分积水状态:为湖泊水位位于其单元内湖底高程最低值和最高值之间时。以图3的湖泊水位为示例,具有湖底平均高程Lb,6的湖泊网格单元为部分积水状态。不同积水状态的湖泊网格单元之间的相互转化,保证了边界条件的连续性。
S2、采集并获取湖泊参数数据,包括时段内的降水强度p、时段内的降水产流系数γ、降水入渗补给系数k、时段内的水面蒸发强度e0、湖泊上游河道汇入量Qsi、湖泊人工取水量W以及湖泊下泄量Qso等。
S3、根据上一次迭代下的湖泊水位判断各湖泊网格单元所处的积水状态。
若则该湖泊网格单元处于完全积水状态;若则该湖泊网格单元处于完全未积水状态;若则该湖泊网格单元处于部分积水状态;其中为上一次迭代下的湖泊水位,为单元内湖底高程的最高值,为单元内湖底高程的最低值。
S4、建立地下水数值计算矩阵方程,本发明实施例中,建立地下水数值计算矩阵方程时不考虑源汇项,具体方程为:
[A]{h}={q} (1)
其中[A]为系数矩阵,{h}为地下水数值计算矩阵,{q}为所有常数项和已知项集合。
S5、求解地下水数值计算矩阵方程,得到当前迭代下的地下水位。
将系数-Cm加入到系数矩阵[A]主对角线系数中,将加入到矩阵方程右端项{q}中,得到当前迭代下的地下水位其中Cm为该湖泊网格单元处湖底与含水层之间的综合水力传导系数,为当前湖泊平均水位,其计算公式为:
其中η∈[0,1]为隐式加权因子,为时段初,即上一次迭代时段末的湖泊水位,为本次迭代时段末的湖泊水位。
S6、根据各湖泊网格单元所处的积水状态,结合当前迭代下的地下水位计算得到各湖泊网格单元的水量平衡项。
本发明实施例中,将湖泊网格单元分为完全积水状态、完全未积水状态和部分积水状态来分别计算与地下水有关的湖泊水量平衡项,根据地下水位和湖底高程的相对关系,一共分为十二种情况进行计算,其中完全积水状态和完全未积水状态分别分为三种情况进行计算,部分积水状态分为六种情况进行计算,具体如下:
(1)若湖泊网格单元处于完全积水状态,当地下水位高于单元内湖底高程的最高值时,该湖泊网格单元上湖泊水体渗漏到地下水或地下水渗出排泄到湖泊。此时无论地下水位高于(图4a)或低于湖泊平均水位(图4b),该单元处湖泊水位与地下水位间都具有完全的水力联系,可直接由水位差和达西公式原理计算湖泊与地下水之间的渗流量。当该单元处的地下水位高于湖泊平均水位时,地下水渗出排泄到湖泊;当该单元处的地下水位低于湖泊平均水位时,该单元处湖泊水体渗漏到地下水,则该湖泊网格单元的水量平衡项为:
其中和分别为对应湖泊平均水位第m个湖泊网格单元积水面积部分地下水渗出排泄到湖泊的流量和湖泊水体渗漏到地下水的流量。
(2)若湖泊网格单元处于完全积水状态,当地下水位低于单元内湖底高程的最低值时,如图5所示,此时在该单元处湖泊为稳定渗漏状态,计算时假设渗漏流量与地下水位无关,而与湖底高程有关,该湖泊网格单元的水量平衡项为:
(3)若湖泊网格单元处于完全积水状态,当地下水位位于单元内湖底高程的最高值和单元内湖底高程的最低值之间时,如图6所示,此时该单元处的湖泊-地下水作用关系为湖泊的渗漏,但地下水位之下面积部分的渗漏与地下水位有关,地下水位之上面积部分的渗漏与地下水位无关,单元处湖泊的渗漏流量为两部分之和,该湖泊网格单元的水量平衡项为:
其中Ra为湖泊网格单元处地下水位之下面积部分所占的面积比例,为地下水位之下面积部分的湖泊渗漏流量,为地下水位之上面积部分的湖泊渗漏流量。
(4)若湖泊网格单元处于完全未积水状态,当地下水位高于单元内湖底高程的最高值时,如图7所示,此时湖泊网格单元处地下水一是将通过湖底渗出排泄,假设此时地下水渗出量在时段内完全流入湖泊地表水体;二是因湖底全部湿润,因此湖底处作用有最大潜水蒸发强度。该种情况下地下水渗出排泄流量由单元处的地下水位以及湖底高程计算,而单元处的潜水蒸发量按潜水蒸发深度0计算,该湖泊网格单元的水量平衡项为:
其中和分别为对应湖泊平均水位第m个湖泊网格单元未积水面积部分地下水渗出排泄到湖泊的流量和降水入渗量,E0为单元上作用的最大潜水蒸发强度,为该单元的面积。
(5)若湖泊网格单元处于完全未积水状态,当地下水位低于单元内湖底高程的最低值时,如图8所示,此时湖泊网格单元处地下水接受降水入渗补给和潜水蒸发,该湖泊网格单元的水量平衡项为:
其中为对应湖泊平均水位第m个湖泊网格单元未积水面积部分的潜水蒸发量,p为时段内的降水强度,k为降水入渗补给系数,Ep为单元上作用的潜水蒸发强度,计算公式为:
其中DM为潜水蒸发极限埋深,D为实际水位埋深,此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
(6)若湖泊网格单元处于完全未积水状态,当地下水位位于单元内湖底高程的最高值和单元内湖底高程的最低值之间时,如图9所示,此时湖泊网格单元湖底位于地下水位以下部分地下水渗出排泄到湖泊,位于地下水位以上部分作用降水入渗补给和潜水蒸发,该湖泊网格单元的水量平衡项为:
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
(7)若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位高于湖泊平均水位时,如图10所示,此时湖泊网格单元积水面积部分地下水渗出到湖泊,该湖泊网格单元的水量平衡项为:
其中Ra,p为湖泊网格单元积水面积部分占整个网格单元的面积比例。
(8)若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位低于单元内湖底高程的最低值时,如图11所示,此时湖泊网格单元积水面积部分湖泊水体渗漏到地下水,且与地下水位无关,该湖泊网格单元的水量平衡项为:
(9)若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位位于湖泊平均水位和单元内湖底高程的最低值之间时,如图12所示,此时湖泊水体在该单元处渗漏到地下水,但地下水位以下的积水面积部分按与地下水位有关的公式计算湖泊渗漏流量,地下水位以上的积水面积部分按与地下水位无关的公式计算湖泊渗漏流量,单元处湖泊的总渗漏流量为两者之和,该湖泊网格单元的水量平衡项为:
其中Ra,p1为湖泊网格单元中地下水位以下的积水面积部分占整个网格单元的面积比例,Ra,p2为湖泊网格单元中地下水位以上的积水面积部分占整个网格单元的面积比例。
(10)若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位高于单元内湖底高程的最高值时,如图13所示,此时未积水面积部分地下水排泄到湖泊,且同时作用有最大潜水蒸发强度,该湖泊网格单元的水量平衡项为:
其中Ra,n为湖泊网格单元未积水面积部分占整个网格单元的面积比例。
(11)若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位低于湖泊平均水位时,如图14所示,此时未积水面积部分作用降水入渗补给和潜水蒸发,该湖泊网格单元的水量平衡项为:
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
(12)若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位位于湖泊平均水位和单元内湖底高程的最高值之间时,如图15所示,此时低于地下水位的未积水面积部分作用地下水排泄和最大潜水蒸发,高于地下水位的未积水面积作用降水入渗补给和潜水蒸发,该湖泊网格单元的水量平衡项为:
其中Ra,n1为湖泊网格单元中低于地下水位的未积水面积部分占整个网格单元的面积比例,Ra,n2为湖泊网格单元中高于地下水位的未积水面积部分占整个网格单元的面积比例;
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
S7、根据各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量,计算公式为:
其中M为湖泊网格单元总数,为时段内湖泊积水区含水层向湖泊的渗出流量,为时段内湖泊积水区的渗漏流量,为时段内湖泊未积水区含水层的渗出流量,为时段内湖泊未积水区的降水入渗量,为时段内湖泊未积水区的潜水蒸发量。
S8、根据湖泊与地下水的水分交换统计量和湖泊参数数据更新当前迭代下的湖泊水位。
如图16所示,将湖泊各项源汇项表达为湖泊水位的函数,可得待求解的湖泊稳定流水量平衡控制方程:
针对公式(21),其中Qsi、W以及Qso三项为用户输入项,需要直接给定,其他各项采用如下公式计算:
其中为时段内湖泊水位为时,湖泊积水区的水面降水通量,p为时段内的降水强度,AP为时段内湖泊的平均水面面积,计算公式为:
为湖泊水位-水面面积曲线上对应的水面面积值。
本发明实施例中,假设湖泊水位-水面面积关系曲线是线性连续的,如图17所示,不同湖泊水位对应的湖泊水面面积可通过关系曲线上相邻的两个离散点的线性插值来确定。以图3的湖底平均高程分布为示例,曲线上各离散点的值确定如下:第1个离散点的湖泊水位值为Lb,1-(Lb,2-Lb,1)/2,定义为湖泊水位的最低值,对应的水面面积为0;第2个离散点的湖泊水位值为(Lb,2+Lb,1)/2,对应的水面面积为具有平均湖底高程Lb,1的湖泊网格单元面积之和;第3个离散点的湖泊水位值为(Lb,2+Lb,3)/2,对应的水面面积为具有平均湖底高程Lb,1、Lb,2的湖泊网格单元面积之和。第4~7个离散点的湖泊水位值与水面面积的关系按以上类推。最后一个离散点(第8个)的湖泊水位值为Lb,7+(Lb,7-Lb,6)/2,对应的水面面积为所有湖泊网格单元的总面积。当湖泊水位超过离散点中的最高水位时,认为湖泊水面面积保持最大面积不变。这里所说的最大面积,为所有湖泊网格单元的总面积。
其中为时段内湖泊水位为时,湖泊未积水区上的降水产流量,γ为时段内的降水产流系数,为时段内未积水区的平均面积,AT为湖泊网格单元总面积。
其中为时段内湖泊水位为时,湖泊积水区的水面蒸发通量,e0为时段内的水面蒸发强度。
结合公式(21)~(25)得到湖泊水量平衡控制方程的具体公式为:
以上湖泊稳定流水量平衡控制方程可采用牛顿迭代法对湖泊水位进行求解:
其中为上一次迭代的湖泊水位,为本次迭代的湖泊水位,F′(·)为F(·)的导数,则公式(27)只需确定F′(·)的函数形式即可进行迭代求解。
对公式(26)求导得:
公式(28)中AP′的确定方法如下:
其中表示湖泊水位-积水面积线性相关曲线上对应时的折线斜率。
的计算公式为:
其中,A为完全积水状态的湖泊网格单元中,单元处地下水位高于湖泊水位的湖泊网格单元的总数量;Ca为上述A个湖泊网格单元中第a个处湖底与含水层之间的综合水力传导系数;B为部分积水状态的湖泊网格单元中,单元处地下水位高于湖泊水位的湖泊网格单元的总数量;Cb为上述B个湖泊网格单元中第b个处湖底与含水层之间的综合水力传导系数;和分别为第b个湖泊网格单元的单元内湖底最高和最低高程。
的计算公式为:
其中,D为完全积水状态的湖泊网格单元中,单元处地下水位低于湖泊水位的湖泊网格单元的总数量;Cd为上述D个湖泊网格单元中第d个处湖底与含水层之间的综合水力传导系数;E为部分积水状态的湖泊网格单元中,单元处地下水位低于湖泊水位的湖泊网格单元的总数量;Ce为上述E个湖泊网格单元中第e个处湖底与含水层之间的综合水力传导系数;和分别为第e个湖泊网格单元的单元内湖底最高和最低高程。
的计算公式为:
综合公式(28)~(32)得:
由于A与D之和为完全积水状态的湖泊网格单元的总数量,B与E之和为部分积水状态的湖泊网格单元的总数量,则公式(33)可进一步化简为:
其中Cs为第s个完全积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数,S为完全积水状态的湖泊网格单元的总数量,Cf为第f个部分积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数,F为部分积水状态的湖泊网格单元的总数量,和分别为第f个部分积水状态的湖泊网格单元的单元内湖底最高和最低高程。
S9、判断当前迭代下湖泊水位是否收敛,收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值,若是则模拟结束,否则返回步骤S3进入下一次迭代。
通过对湖泊刻画精度的增加,传统湖泊地下水相互作用方法需要通过增加含水层的垂向剖分层数来实现,而本发明则不需要;并且随着垂向剖分层数的增加,传统湖泊地下水相互作用方法的计算出现了显著的不稳定性。由此可以看出,与传统方法相比本发明具有显著的优势。
本发明与传统湖泊地下水相互作用方法根据各自方法属性,对湖泊网格的剖分如图18~20所示,图18中A和B为模拟过程中的观测点,阴影部分表示湖泊单元,图20中阴影部分表示湖泊单元。其他所有参数和模型驱动数据均保持一致,考虑到发明中湖泊水面蒸发是影响湖泊水位主要因素,继续尝试分别增加和减少湖泊的水面蒸发速率,以期对湖泊的不同蓄水状态进行全面的模拟和验证,对比两种方法的模拟计算结果(如图21~23所示),可以看出两种方法模拟结果差别不大,从而证明了本发明方法的合理性。但在模拟过程中发现,LAK3并不是全局收敛,只有给其合适的初始值很靠近真实值时,迭代才会收敛,这一缺陷将很大程度上限制LAK3在实际工作的应用,而本发明则可以全局收敛。
本领域的普通技术人员将会意识到,这里所述的实施例是为了帮助读者理解本发明的原理,应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合,这些变形和组合仍然在本发明的保护范围内。
Claims (10)
1.一种基于地下水模型的湖泊与地下水稳定流作用模拟方法,其特征在于,包括以下步骤:
S1、对湖泊进行离散化处理及边界连续性处理,得到各湖泊网格单元;
S2、采集并获取湖泊参数数据;
S3、根据上一次迭代下的湖泊水位判断各湖泊网格单元所处的积水状态;
S4、建立地下水数值计算矩阵方程;
S5、求解地下水数值计算矩阵方程,得到当前迭代下的地下水位;
S6、根据各湖泊网格单元所处的积水状态,结合当前迭代下的地下水位计算得到各湖泊网格单元的水量平衡项;
S7、根据各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量;
S8、根据湖泊与地下水的水分交换统计量和湖泊参数数据更新当前迭代下的湖泊水位;
S9、判断当前迭代下湖泊水位是否收敛,若是则模拟结束,否则返回步骤S3进入下一次迭代。
2.根据权利要求1所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S1具体包括:
S1-1、对湖泊进行离散化处理:将湖底所在的含水层网格单元标识为湖泊网格单元,并将湖底高程数据以离散的方式赋予每个湖泊网格单元,使得每个湖泊网格单元都具有其单元面积范围内的平均湖底高程值;如果某个湖泊网格单元不位于第一层,则该湖泊网格单元上方的所有网格单元都将被定义为无效单元;
S1-2、对湖泊进行边界连续性处理:假设湖底高程在湖泊网格单元内是倾斜的,具有单元内的最低值和最高值;对于非最上级和非最下级的湖泊网格单元,其湖底高程最低值为本单元湖底平均高程与下一级单元湖底平均高程的中间位置,最高值为本单元湖底平均高程与上一级单元湖底平均高程的中间位置;对于最下级的湖泊网格单元,其湖底高程最低值为其平均湖底高程减去其与上一级单元平均湖底高程差的一半;对于最上级湖泊网格单元,其湖底高程最高值为其平均湖底高程加上其与下一级单元平均湖底高程差的一半。
3.根据权利要求1所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S2中的湖泊参数数据包括:时段内的降水强度p、时段内的降水产流系数γ、降水入渗补给系数k、时段内的水面蒸发强度e0、湖泊上游河道汇入量Qsi、湖泊人工取水量W以及湖泊下泄量Qso。
4.根据权利要求3所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S3具体为:
若则该湖泊网格单元处于完全积水状态;若则该湖泊网格单元处于完全未积水状态;若则该湖泊网格单元处于部分积水状态;其中为上一次迭代下的湖泊水位,为单元内湖底高程的最高值,为单元内湖底高程的最低值。
5.根据权利要求4所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S4中地下水数值计算矩阵方程为:
[A]{h}={q}
其中[A]为系数矩阵,{h}为地下水数值计算矩阵,{q}为所有常数项和已知项集合。
6.根据权利要求5所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S5具体为:
将系数-Cm加入到系数矩阵[A]主对角线系数中,将加入到矩阵方程右端项{q}中,得到当前迭代下的地下水位其中Cm为该湖泊网格单元处湖底与含水层之间的综合水力传导系数,为当前湖泊平均水位,其计算公式为:
其中η∈[0,1]为隐式加权因子,为时段初,即上一次迭代时段末的湖泊水位,为本次迭代时段末的湖泊水位。
7.根据权利要求6所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S6具体为:
若湖泊网格单元处于完全积水状态,当地下水位高于单元内湖底高程的最高值时,该湖泊网格单元的水量平衡项为:
其中和分别为对应湖泊平均水位第m个湖泊网格单元积水面积部分地下水渗出排泄到湖泊的流量和湖泊水体渗漏到地下水的流量;
若湖泊网格单元处于完全积水状态,当地下水位低于单元内湖底高程的最低值时,该湖泊网格单元的水量平衡项为:
若湖泊网格单元处于完全积水状态,当地下水位位于单元内湖底高程的最高值和单元内湖底高程的最低值之间时,该湖泊网格单元的水量平衡项为:
其中Ra为湖泊网格单元处地下水位之下面积部分所占的面积比例,为地下水位之下面积部分的湖泊渗漏流量,为地下水位之上面积部分的湖泊渗漏流量;
若湖泊网格单元处于完全未积水状态,当地下水位高于单元内湖底高程的最高值时,该湖泊网格单元的水量平衡项为:
其中和分别为对应湖泊平均水位第m个湖泊网格单元未积水面积部分地下水渗出排泄到湖泊的流量和降水入渗量,E0为单元上作用的最大潜水蒸发强度,为该单元的面积;
若湖泊网格单元处于完全未积水状态,当地下水位低于单元内湖底高程的最低值时,该湖泊网格单元的水量平衡项为:
其中为对应湖泊平均水位第m个湖泊网格单元未积水面积部分的潜水蒸发量,p为时段内的降水强度,k为降水入渗补给系数,Ep为单元上作用的潜水蒸发强度,计算公式为:
其中DM为潜水蒸发极限埋深,D为实际水位埋深,此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
若湖泊网格单元处于完全未积水状态,当地下水位位于单元内湖底高程的最高值和单元内湖底高程的最低值之间时,该湖泊网格单元的水量平衡项为:
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位高于湖泊平均水位时,该湖泊网格单元的水量平衡项为:
其中Ra,p为湖泊网格单元积水面积部分占整个网格单元的面积比例;
若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位低于单元内湖底高程的最低值时,该湖泊网格单元的水量平衡项为:
若湖泊网格单元处于部分积水状态,在积水面积部分,当地下水位位于湖泊平均水位和单元内湖底高程的最低值之间时,该湖泊网格单元的水量平衡项为:
其中Ra,p1为湖泊网格单元中地下水位以下的积水面积部分占整个网格单元的面积比例,Ra,p2为湖泊网格单元中地下水位以上的积水面积部分占整个网格单元的面积比例;
若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位高于单元内湖底高程的最高值时,该湖泊网格单元的水量平衡项为:
其中Ra,n为湖泊网格单元未积水面积部分占整个网格单元的面积比例;
若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位低于湖泊平均水位时,该湖泊网格单元的水量平衡项为:
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
若湖泊网格单元处于部分积水状态,在未积水面积部分,当地下水位位于湖泊平均水位和单元内湖底高程的最高值之间时,该湖泊网格单元的水量平衡项为:
其中Ra,n1为湖泊网格单元中低于地下水位的未积水面积部分占整个网格单元的面积比例,Ra,n2为湖泊网格单元中高于地下水位的未积水面积部分占整个网格单元的面积比例;
此时用于计算潜水蒸发强度Ep时的实际水位埋深D取值为:
8.根据权利要求7所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S7中湖泊与地下水的水分交换统计量的计算公式为:
其中M为湖泊网格单元总数,为时段内湖泊积水区含水层向湖泊的渗出流量,为时段内湖泊积水区的渗漏流量,为时段内湖泊未积水区含水层的渗出流量,为时段内湖泊未积水区的降水入渗量,为时段内湖泊未积水区的潜水蒸发量。
9.根据权利要求8所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S8中更新当前迭代下的湖泊水位的公式为:
其中为上一次迭代的湖泊水位,为本次迭代的湖泊水位,F(·)表示湖泊水量平衡控制方程,具体公式为:
其中,p为时段内的降水强度,γ为时段内的降水产流系数,e0为时段内的水面蒸发强度,AP为时段内湖泊的平均水面面积,AT为湖泊网格单元总面积,Qsi为湖泊上游河道汇入量,W为湖泊人工取水量,Qso为湖泊下泄量;
F′(·)为F(·)的导数,具体公式为:
其中表示湖泊水位-积水面积线性相关曲线上对应时的折线斜率,Cs为第s个完全积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数,S为完全积水状态的湖泊网格单元的总数量,Cf为第f个部分积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数,F为部分积水状态的湖泊网格单元的总数量,和分别为第f个部分积水状态的湖泊网格单元的单元内湖底最高和最低高程。
10.根据权利要求9所述的湖泊与地下水稳定流作用模拟方法,其特征在于,所述步骤S9中收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810561011.5A CN108763797B (zh) | 2018-06-04 | 2018-06-04 | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810561011.5A CN108763797B (zh) | 2018-06-04 | 2018-06-04 | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108763797A true CN108763797A (zh) | 2018-11-06 |
CN108763797B CN108763797B (zh) | 2020-03-17 |
Family
ID=64002198
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810561011.5A Active CN108763797B (zh) | 2018-06-04 | 2018-06-04 | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108763797B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110532641A (zh) * | 2019-08-06 | 2019-12-03 | 中国水利水电科学研究院 | 一种地表网格分层建模方法及系统 |
CN112529723A (zh) * | 2020-04-30 | 2021-03-19 | 中国科学院地球化学研究所 | 一种基于像元尺度的地下水补给量估算方法及系统 |
CN112819315A (zh) * | 2021-01-29 | 2021-05-18 | 国家基础地理信息中心 | 一种稳定水系的水系稳定度计算方法 |
CN116992783A (zh) * | 2023-05-23 | 2023-11-03 | 中国水利水电科学研究院 | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102567634A (zh) * | 2011-12-23 | 2012-07-11 | 中国水利水电科学研究院 | 一种基于水循环的地下水数值仿真方法 |
CN107315847A (zh) * | 2017-01-12 | 2017-11-03 | 中国水利水电科学研究院 | 一种河流与地下水耦合模拟参数的生成方法及装置 |
-
2018
- 2018-06-04 CN CN201810561011.5A patent/CN108763797B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102567634A (zh) * | 2011-12-23 | 2012-07-11 | 中国水利水电科学研究院 | 一种基于水循环的地下水数值仿真方法 |
CN107315847A (zh) * | 2017-01-12 | 2017-11-03 | 中国水利水电科学研究院 | 一种河流与地下水耦合模拟参数的生成方法及装置 |
Non-Patent Citations (4)
Title |
---|
RANDALL J.HUNT等: "《Simulating Ground Water-Lake Interactions: Approaches and Insights》", 《GROUND WATER》 * |
张奇: "《湖泊集水域地表—地下径流联合模拟》", 《地理科学进展》 * |
李慧 等: "《分布式"河道-沉陷区-地下水"水循环耦合模型——Ⅰ.模型原理与开发》", 《水科学进展》 * |
陆垂裕 等: "《面向对象模块化的分布式水文模型MODCYCLE I:模型原理与开发篇》", 《水利学报》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110532641A (zh) * | 2019-08-06 | 2019-12-03 | 中国水利水电科学研究院 | 一种地表网格分层建模方法及系统 |
CN110532641B (zh) * | 2019-08-06 | 2020-08-25 | 中国水利水电科学研究院 | 一种地表网格分层建模方法及系统 |
CN112529723A (zh) * | 2020-04-30 | 2021-03-19 | 中国科学院地球化学研究所 | 一种基于像元尺度的地下水补给量估算方法及系统 |
CN112819315A (zh) * | 2021-01-29 | 2021-05-18 | 国家基础地理信息中心 | 一种稳定水系的水系稳定度计算方法 |
CN112819315B (zh) * | 2021-01-29 | 2024-06-04 | 国家基础地理信息中心 | 一种稳定水系的水系稳定度计算方法 |
CN116992783A (zh) * | 2023-05-23 | 2023-11-03 | 中国水利水电科学研究院 | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 |
Also Published As
Publication number | Publication date |
---|---|
CN108763797B (zh) | 2020-03-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108763797A (zh) | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 | |
Mao et al. | Loosely coupled SaltMod for simulating groundwater and salt dynamics under well-canal conjunctive irrigation in semi-arid areas | |
CN102567634B (zh) | 一种基于水循环的地下水数值仿真方法 | |
Yi et al. | Two-way coupling of unsaturated-saturated flow by integrating the SWAT and MODFLOW models with application in an irrigation district in arid region of West China | |
CN110174506A (zh) | 一种喀斯特地区土壤有机碳估算方法 | |
CN104809516B (zh) | 引黄灌区水资源多目标优化配置模型及其求解方法 | |
Janssen et al. | Water losses through paddy bunds: Methods, experimental data, and simulation studies | |
CN108763798A (zh) | 一种湖泊与地下水非稳定流作用模拟方法 | |
CN104698507B (zh) | 一种采煤沉陷区水资源效应的定量方法 | |
CN117131652A (zh) | 一种预测不同气候变化影响浅层地下水海水入侵的方法及应用 | |
CN107229811B (zh) | 一种饱和软土自重固结过程的预测方法 | |
CN106320391A (zh) | 碾压式土石坝碾压次数的实验确定方法 | |
Ding et al. | The assessment of ecological water replenishment scheme based on the two-dimensional lattice-Boltzmann water age theory | |
CN106759072A (zh) | 兼顾床面稳定和快速排水的河工模型排水系统及设计方法 | |
Liu et al. | Dynamics of the soil water and solute in the sodic saline soil in the Songnen Plain, China | |
CN109147047B (zh) | 一种渗渠型傍河水源地数值模型构建方法 | |
CN110532500A (zh) | 一种区域大气水文模型中水库群参数化方案的构建方法 | |
Xu et al. | Modeling of surface runoff in Xitiaoxi catchment, China | |
CN206616513U (zh) | 兼顾床面稳定和快速排水的河工模型排水系统 | |
Abdul Jabbar Jamel | Analysis and estimation of downward seepage from lining and unlining triangular open channel | |
CN206448251U (zh) | 高低闸门联合运用的灌溉分水闸 | |
Zhang | Application of information entropy and TOPSIS coupling model in impervious design of hydraulic engineering | |
CN107038291A (zh) | 高填方涵洞涵顶垂直土压力计算方法 | |
Cui et al. | Application of AutoBank in The Dike Reinforcement Calculation of Lower Qinhe River | |
Maanjhu et al. | Groundwater simulation studies of parts of Western Yamuna Canal Command area, Haryana for planning sustainable development |
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 |