CN108763798B - 一种湖泊与地下水非稳定流作用模拟方法 - Google Patents

Abstract

本发明公开了一种湖泊与地下水非稳定流作用模拟方法，详细刻画了湖泊水量平衡过程，有效降低了湖泊的空间离散化对地下水网格系统剖分的依赖性，并采用倾斜式湖底的算法解决湖泊‑地下水作用边界条件的连续性问题，一方面在不降低湖泊空间离散化精度的情况下可减少不必要的地下水含水层垂向剖分层数，从而有效降低了模拟计算的工作量，大幅增强了应用时的便利性，另一方面也有效避免了模拟过程中由于垂向剖分层数较多导致的数值震荡问题，提高了模拟计算的收敛稳定性。与现有的湖泊‑地下水相互作用模拟方法相比，本发明在应用便利性、计算稳定性上具有显著优势，可以广泛应用于湖泊与地下水相互作用的定量解析研究。

Description

一种湖泊与地下水非稳定流作用模拟方法

技术领域

本发明属于地下水数值模拟技术领域，具体涉及一种湖泊与地下水非稳定流作用模拟方法的设计。

背景技术

当湖泊和地表含水层具有直接水力联系时，可对湖泊和地下水系统产生显著影响，因此发展一种技术去定量湖泊和地下水之间的相互作用，评估一方平衡条件发生改变对另一方的产生的影响，对区域水资源管理具有重要意义。湖泊与水库等其他地表水体不同，其水位不是一个确定值，而是与其自身水量平衡过程有关，具有较大的评估难度。目前，分析湖泊与地下水相互作用的方法主要有试验法和模型模拟法。然而，相比于模型模拟法，试验法不仅会耗费大量的人力物力，且很容易出错，故较少采用。模型模拟法根据模型对湖泊水位的处理方式的不同，可以分为指定水位法和计算水位法。其中指定水位法中湖泊水位是由用户给定，尽管有湖底渗漏或其他的应力存在，湖泊水位只能保持不变或者在用户设定的范围内线性变化，这类方法最大的不足在于，不能模拟湖泊水位的变化，只能粗略刻画湖泊对地下水系统的影响。计算水位法是指湖泊水位通过计算得到，而不需要用户从外部输入，主要包括高K法、基于GFLOW的解析元素法、MODFLOW中LAK3模块法等。其中高K法是将湖泊用具有湖泊水力特点的模型网格单元或节点来表示，作为含水层的一部分赋予这些网格单元或节点较高的水力传导系数，且将地下水流动方程的解作为湖泊水位，然而，这种方法只能进行渗漏湖的模拟；解析元素法是依据泊松-拉普拉斯方法中势函数可以叠加的事实，以建立在裘布衣假设上的井的解析元模型为基础，将湖泊用一系列函数来表示，这种方法虽然具有较高的效率和精度，但是只能进行二维稳定流的模拟。以上方法对湖泊的刻画都过于粗略，且在物理机制方面存在明显不足，随着实际应用变得越来越复杂，方法本身的局限性愈发突出，需要开发一种更加先进的方法来解决湖泊与地下水相互作用的问题。

LAK3软件包作为MODFLOW的一个功能模块，是一种目前公认的通用方法，它将湖泊用模型网格系统中一系列无效单元格来表示，不仅可以模拟湖泊水位，而且能够模拟湖泊渗漏以及其他地表水体对湖泊的影响。虽然LAK3模块相比之前的方法已取得了很大的进步，但正如开发者所说，LAK3模块自身仍然存在局限性，如地下水网格系统剖分和湖泊单元的离散相互依存，为使对湖泊的刻画体现出一定的精度，需要细化整个研究区水平和垂向的网格剖分尺度，从而造成计算量大幅增加，这对大区域模拟来说是致命的，同时垂向剖分层数的增加，会带来更多的干湿单元的转化问题，该问题的处理方法是经验性的，从而使模型具有潜在的不稳定性。

发明内容

本发明的目的是提出一种湖泊与地下水非稳定流作用模拟方法，解决现有技术中模拟湖泊和地下水的相互作用存在的如下问题：

(1)现有的湖泊地下水相互作用方法对湖泊水量平衡的刻画还不够全面；

(2)现有的湖泊地下水相互作用方法地下水网格系统的剖分与湖泊刻画相互依存，若要对湖泊进行相对精细的刻画，就需要对地下水网格系统也进行比较细致的剖分，不仅大大增加了计算量，同时容易出现含水层垂向剖分过细导致的地下水单元干-湿转换计算不稳定的问题；

(3)现有的湖泊地下水相互作用方法未考虑边界条件的连续性。

本发明的技术方案为：一种湖泊与地下水非稳定流作用模拟方法，包括以下步骤：

S1、对湖泊进行离散化处理及边界连续性处理，得到各湖泊网格单元。

S2、采集并获取湖泊参数数据。

S3、计算当前湖泊平均水位。

S4、根据当前湖泊平均水位判断各湖泊网格单元所处的积水状态。

S5、不考虑源汇项建立地下水数值计算矩阵方程。

S6、求解地下水数值计算矩阵方程，得到当前迭代下的地下水位。

S7、根据各湖泊网格单元所处的积水状态，结合当前迭代下的地下水位计算得到各湖泊网格单元的水量平衡项。

S8、根据各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量。

S9、根据湖泊与地下水的水分交换统计量和湖泊参数数据计算得到当前迭代下的湖泊蓄水量。

S10、判断当前迭代下湖泊水位是否收敛，收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值，若是则模拟结束，否则返回步骤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所示为增加湖泊的刻画精度后模拟区平面图。

图24所示为增加湖泊的刻画精度后本发明所述方法湖泊网格单元湖底高程刻画图。

图25所示为增加湖泊的刻画精度后传统湖泊地下水相互作用方法研究区含水层剖面图。

图26所示为增加湖泊的刻画精度后湖泊水位模拟结果对比图。

具体实施方式

现在将参考附图来详细描述本发明的示例性实施方式。应当理解，附图中示出和描述的实施方式仅仅是示例性的，意在阐释本发明的原理和精神，而并非限制本发明的范围。

本发明实施例提供了一种湖泊与地下水非稳定流作用模拟方法，如图1所示，包括以下步骤S1-S10：

S1、对湖泊进行离散化处理及边界连续性处理，得到各湖泊网格单元。

步骤S1具体包括：

S1-1、对湖泊进行离散化处理。

进行湖泊的空间离散化处理是湖泊-地下水作用耦合模拟的前提和基础，很大程度上决定了计算方法的选择，本发明实施例中处理过程为：将湖底所在的含水层网格单元标识为湖泊网格单元，并将湖底高程数据以离散的方式赋予每个湖泊网格单元，使得每个湖泊网格单元都具有其单元面积范围内的平均湖底高程值。在垂向方向上含水层可能剖分为多层，但只有湖底所在层位的网格单元为湖泊网格单元。如果某个湖泊网格单元不位于第一层，则该湖泊网格单元上方的所有网格单元都将被定义为无效单元，如图2所示。

S1-2、对湖泊进行边界连续性处理.

离散化处理过程中，各湖泊网格单元的湖底平均高程按照从下至上的顺序进行排序，各湖泊网格单元也被分为下级单元和上级单元。离散化时可能有若干个湖泊网格单元具有相同的湖底平均高程，它们属于同一级单元。以图3的湖底平均高程分布为示例，共有7级单元，其中最下级的湖泊网格单元的湖底平均高程为L_b,1，最上级湖泊网格单元的湖底平均高程为L_b,7。

本发明实施例中，假设湖底高程在湖泊网格单元内是倾斜的，具有单元内的最低值和最高值；对于非最上级和非最下级的湖泊网格单元，其湖底高程最低值为本单元湖底平均高程与下一级单元湖底平均高程的中间位置，最高值为本单元湖底平均高程与上一级单元湖底平均高程的中间位置；对于最下级的湖泊网格单元，其湖底高程最低值为其平均湖底高程减去其与上一级单元平均湖底高程差的一半；对于最上级湖泊网格单元，其湖底高程最高值为其平均湖底高程加上其与下一级单元平均湖底高程差的一半。

图3演示了不同湖泊网格单元内湖底高程的最高值和最低值的情况。例如，对于最下级单元，其单元内湖底高程的最低值为L_b,1-(L_b,2-L_b,1)/2，最高值为(L_b,2+L_b,1)/2。对于具有湖底平均高程L_b,2的湖泊网格单元，其单元内湖底高程的最低值为(L_b,2+L_b,1)/2，最高值为(L_b,2+L_b,3)/2。对于最上级单元，其单元内湖底高程的最低值为(L_b,6+L_b,7)/2，最高值为L_b,7+(L_b,7-L_b,6)/2。

根据湖泊水位的情况，湖泊网格单元被分为三种积水状态：(1)完全积水状态：湖泊水位高于其单元内湖底高程最高值时。以图3的湖泊水位为示例，具有湖底平均高程L_b,1,L_b,2,L_b,3,L_b,4,L_b,5的湖泊网格单元为完全积水状态。(2)完全未积水状态：湖泊水位低于其单元内湖底高程最低值时。以图3的湖泊水位为示例，具有湖底平均高程L_b,7的湖泊网格单元为完全未积水状态；(3)部分积水状态：为湖泊水位位于其单元内湖底高程最低值和最高值之间时。以图3的湖泊水位为示例，具有湖底平均高程L_b,6的湖泊网格单元为部分积水状态。不同积水状态的湖泊网格单元之间的相互转化，保证了边界条件的连续性。

S2、采集并获取湖泊参数数据，包括时段内的降水强度p、时段内的降水产流系数γ、降水入渗补给系数k、时段内的水面蒸发强度e₀、湖泊上游河道汇入量Q_si、湖泊人工取水量W以及湖泊下泄量Q_so等。

S3、计算当前湖泊平均水位，计算公式为：

其中

为当前湖泊平均水位，η∈[0,1]为隐式加权因子，

为时段初，即上一时段末的湖泊水位，

为时段末的湖泊水位。

S4、根据当前湖泊平均水位判断各湖泊网格单元所处的积水状态。

若

则该湖泊网格单元处于完全积水状态；若

则该湖泊网格单元处于完全未积水状态；若

则该湖泊网格单元处于部分积水状态；其中

为单元内湖底高程的最高值，

为单元内湖底高程的最低值。

S5、不考虑源汇项建立地下水数值计算矩阵方程：

[A]{h}＝{q} (2)

其中[A]为系数矩阵，{h}为地下水数值计算矩阵，{q}为所有常数项和已知项集合。

S6、求解地下水数值计算矩阵方程，得到当前迭代下的地下水位。

将系数-C_m加入到系数矩阵[A]主对角线系数中，将

加入到矩阵方程右端项{q}中，得到当前迭代下的地下水位

其中C_m为该湖泊网格单元处湖底与含水层之间的综合水力传导系数。

S7、根据各湖泊网格单元所处的积水状态，结合当前迭代下的地下水位计算得到各湖泊网格单元的水量平衡项。

本发明实施例中，将湖泊网格单元分为完全积水状态、完全未积水状态和部分积水状态来分别计算与地下水有关的湖泊水量平衡项，根据地下水位和湖底高程的相对关系，一共分为十二种情况进行计算，其中完全积水状态和完全未积水状态分别分为三种情况进行计算，部分积水状态分为六种情况进行计算，具体如下：

(1)若湖泊网格单元处于完全积水状态，当地下水位

高于单元内湖底高程的最高值

时，该湖泊网格单元上湖泊水体渗漏到地下水或地下水渗出排泄到湖泊。此时无论地下水位高于(图4a)或低于湖泊平均水位(图4b)，该单元处湖泊水位与地下水位间都具有完全的水力联系，可直接由水位差和达西公式原理计算湖泊与地下水之间的渗流量。当该单元处的地下水位高于湖泊平均水位时，地下水渗出排泄到湖泊；当该单元处的地下水位低于湖泊平均水位时，该单元处湖泊水体渗漏到地下水，则该湖泊网格单元的水量平衡项为：

其中

和

分别为对应湖泊平均水位

第m个湖泊网格单元积水面积部分地下水渗出排泄到湖泊的流量和湖泊水体渗漏到地下水的流量。

(2)若湖泊网格单元处于完全积水状态，当地下水位

低于单元内湖底高程的最低值

时，如图5所示，此时在该单元处湖泊为稳定渗漏状态，计算时假设渗漏流量与地下水位无关，而与湖底高程有关，该湖泊网格单元的水量平衡项为：

(3)若湖泊网格单元处于完全积水状态，当地下水位

位于单元内湖底高程的最高值

和单元内湖底高程的最低值

之间时，如图6所示，此时该单元处的湖泊-地下水作用关系为湖泊的渗漏，但地下水位之下面积部分的渗漏与地下水位有关，地下水位之上面积部分的渗漏与地下水位无关，单元处湖泊的渗漏流量为两部分之和，该湖泊网格单元的水量平衡项为：

其中R_a为湖泊网格单元处地下水位之下面积部分所占的面积比例，

为地下水位之下面积部分的湖泊渗漏流量，

为地下水位之上面积部分的湖泊渗漏流量。

(4)若湖泊网格单元处于完全未积水状态，当地下水位

高于单元内湖底高程的最高值

时，如图7所示，此时湖泊网格单元处地下水一是将通过湖底渗出排泄，假设此时地下水渗出量在时段内完全流入湖泊地表水体；二是因湖底全部湿润，因此湖底处作用有最大潜水蒸发强度。该种情况下地下水渗出排泄流量由单元处的地下水位以及湖底高程计算，而单元处的潜水蒸发量按潜水蒸发深度0计算，该湖泊网格单元的水量平衡项为：

其中

和

分别为对应湖泊平均水位

第m个湖泊网格单元未积水面积部分地下水渗出排泄到湖泊的流量和降水入渗量，E₀为单元上作用的最大潜水蒸发强度，

为该单元的面积。

(5)若湖泊网格单元处于完全未积水状态，当地下水位

低于单元内湖底高程的最低值

时，如图8所示，此时湖泊网格单元处地下水接受降水入渗补给和潜水蒸发，该湖泊网格单元的水量平衡项为：

其中

为对应湖泊平均水位

第m个湖泊网格单元未积水面积部分的潜水蒸发量，p为时段内的降水强度，k为降水入渗补给系数，E_p为单元上作用的潜水蒸发强度，计算公式为：

其中D_M为潜水蒸发极限埋深，D为实际水位埋深，此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

(6)若湖泊网格单元处于完全未积水状态，当地下水位

位于单元内湖底高程的最高值

和单元内湖底高程的最低值

之间时，如图9所示，此时湖泊网格单元湖底位于地下水位以下部分地下水渗出排泄到湖泊，位于地下水位以上部分作用降水入渗补给和潜水蒸发，该湖泊网格单元的水量平衡项为：

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

(7)若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

高于湖泊平均水位

时，如图10所示，此时湖泊网格单元积水面积部分地下水渗出到湖泊，该湖泊网格单元的水量平衡项为：

其中R_a,p为湖泊网格单元积水面积部分占整个网格单元的面积比例。

(8)若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

低于单元内湖底高程的最低值

时，如图11所示，此时湖泊网格单元积水面积部分湖泊水体渗漏到地下水，且与地下水位无关，该湖泊网格单元的水量平衡项为：

(9)若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

位于湖泊平均水位

和单元内湖底高程的最低值

之间时，如图12所示，此时湖泊水体在该单元处渗漏到地下水，但地下水位以下的积水面积部分按与地下水位有关的公式计算湖泊渗漏流量，地下水位以上的积水面积部分按与地下水位无关的公式计算湖泊渗漏流量，单元处湖泊的总渗漏流量为两者之和，该湖泊网格单元的水量平衡项为：

其中R_a,p1为湖泊网格单元中地下水位以下的积水面积部分占整个网格单元的面积比例，R_a,p2为湖泊网格单元中地下水位以上的积水面积部分占整个网格单元的面积比例。

(10)若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

高于单元内湖底高程的最高值

时，如图13所示，此时未积水面积部分地下水排泄到湖泊，且同时作用有最大潜水蒸发强度，该湖泊网格单元的水量平衡项为：

其中R_a,n为湖泊网格单元未积水面积部分占整个网格单元的面积比例。

(11)若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

低于湖泊平均水位

时，如图14所示，此时未积水面积部分作用降水入渗补给和潜水蒸发，该湖泊网格单元的水量平衡项为：

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

(12)若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

位于湖泊平均水位

和单元内湖底高程的最高值

之间时，如图15所示，此时低于地下水位的未积水面积部分作用地下水排泄和最大潜水蒸发，高于地下水位的未积水面积作用降水入渗补给和潜水蒸发，该湖泊网格单元的水量平衡项为：

其中R_a,n1为湖泊网格单元中低于地下水位的未积水面积部分占整个网格单元的面积比例，R_a,n2为湖泊网格单元中高于地下水位的未积水面积部分占整个网格单元的面积比例；

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

S8、根据各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量，计算公式为：

其中M为湖泊网格单元总数，

为时段内湖泊积水区含水层向湖泊的渗出流量，

为时段内湖泊积水区的渗漏流量，

为时段内湖泊未积水区含水层的渗出流量，

为时段内湖泊未积水区的降水入渗量，

为时段内湖泊未积水区的潜水蒸发量。

S9、根据湖泊与地下水的水分交换统计量和湖泊参数数据计算得到当前迭代下的湖泊蓄水量。

如图16所示，在非稳定流情况下，详细分析湖泊积水区的相关水量平衡项，可以总结出5项补给项和4项排泄项。在任何时段内，根据湖泊水量平衡原理，以下湖泊积水区的水量平衡控制方程成立：

其中Vⁿ为湖泊积水区在时段末的积水量，V^n-1为湖泊积水区在时段初的积水量，Δt为当前计算时段，P为时段内湖泊积水区的水面降水通量，Q_si为时段内湖泊上游河道汇入量，Rnf为时段内湖泊未积水区的产流汇入流量，

为时段内湖泊积水区含水层向湖泊的渗出流量，

为时段内湖泊未积水区含水层的渗出流量，E为时段内湖泊积水区的水面蒸发通量，

为时段内湖泊积水区的渗漏流量，W为时段内湖泊人工取水量，包括生产、生活、生态等用途，Q_so为时段内湖泊下泄量，指从湖泊出水口(闸门、泵站)等下泄的水量。

针对公式(21)，其中Q_si、W以及Q_so三项为用户输入项，需要直接给定，其他各项采用如下公式计算：

其中

为时段内湖泊水位为

时，湖泊积水区的水面降水通量，p为时段内的降水强度，

为时段内湖泊的平均水面面积，计算公式为：

为湖泊水位-水面面积曲线上对应

的水面面积值。

本发明实施例中，假设湖泊水位-水面面积关系曲线是线性连续的，如图17所示，不同湖泊水位对应的湖泊水面面积可通过关系曲线上相邻的两个离散点的线性插值来确定。以图3的湖底平均高程分布为示例，曲线上各离散点的值确定如下：第1个离散点的湖泊水位值为L_b,1-(L_b,2-L_b,1)/2，定义为湖泊水位的最低值，对应的水面面积为0；第2个离散点的湖泊水位值为(L_b,2+L_b,1)/2，对应的水面面积为具有平均湖底高程L_b,1的湖泊网格单元面积之和；第3个离散点的湖泊水位值为(L_b,2+L_b,3)/2，对应的水面面积为具有平均湖底高程L_b,1、L_b,2的湖泊网格单元面积之和。第4～7个离散点的湖泊水位值与水面面积的关系按以上类推。最后一个离散点(第8个)的湖泊水位值为L_b,7+(L_b,7-L_b,6)/2，对应的水面面积为所有湖泊网格单元的总面积。当湖泊水位超过离散点中的最高水位时，认为湖泊水面面积保持最大面积不变。这里所说的最大面积，为所有湖泊网格单元的总面积。

其中

为时段内湖泊水位为

时，湖泊未积水区上的降水产流量，γ为时段内的降水产流系数，

为时段内未积水区的平均面积，A_T为湖泊网格单元总面积。

其中

为时段内湖泊水位为

时，湖泊积水区的水面蒸发通量，e₀为时段内的水面蒸发强度。

结合公式(21)～(25)得到当前迭代下的湖泊蓄水量的计算公式为：

其中

为当前迭代下的湖泊蓄水量，

为前一次迭代下的湖泊蓄水量。

S10、判断当前迭代下湖泊水位是否收敛，收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值，若是则模拟结束，否则返回步骤S3进入下一次迭代。

当前迭代下湖泊水位

可由湖泊蓄水量

通过湖泊水位-蓄水量关系曲线插值得到。如图18所示，类似于湖泊水位-水面面积关系曲线，湖泊水位-蓄量关系曲线也是线性连续的，不同湖泊水位对应的湖泊蓄量可通过关系曲线上相邻的两个离散点的线性插值来确定。各离散点处的湖泊蓄量的计算过程为，先将离散点处的湖泊水位减去湖泊网格单元的湖底平均高程，得到湖泊网格单元上的水深，再将水深乘以湖泊网格单元的面积从而计算湖泊网格单元上的蓄量，再对各湖泊网格单元上的蓄量进行累加，得到湖泊总的蓄量。与湖泊水位-水面面积关系曲线不同的是，当湖泊水位超过离散点中的最高水位时，湖泊蓄水量的增幅将按照湖泊水面的最大面积乘以湖泊水位的增幅来确定。

通过对湖泊刻画精度的增加，传统湖泊地下水相互作用方法需要通过增加含水层的垂向剖分层数来实现，而本发明则不需要；并且随着垂向剖分层数的增加，传统湖泊地下水相互作用方法的计算出现了显著的不稳定性。由此可以看出，与传统方法相比本发明具有显著的优势。

本发明与传统湖泊地下水相互作用方法根据各自方法属性，对湖泊网格的剖分如图19～21所示，图19中A和B为模拟过程中的观测点，阴影部分表示湖泊单元，图21中阴影部分表示湖泊单元。其他所有参数和模型驱动数据均保持一致，对比两种方法的模拟计算结果，如图22所示，可以看出两种方法模拟结果差别不大，从而证明了本发明方法的合理性。

以上计算结果可以看到，两种方法虽然差异不大，但在模拟结果上仍存在一定范围内的差异，我们也分析得出这些差异是不同的湖泊离散格式造成的。继续进行对比模拟计算，增加湖泊的刻画精度，两种方法对湖泊单元的刻画，如图23-25所示，模拟结果如图26所示。图23中A和B为模拟过程中的观测点，阴影部分表示湖泊单元，图25中阴影部分表示湖泊单元，图26中LAK3表示传统湖泊地下水相互作用方法，SLM表示本发明提供的方法。可以明显看出，对湖泊网格单元进行细化之后，传统湖泊地下水相互作用方法的模拟结果表现出了不稳定性(即取不同的时间步长，得到的模拟结果不同)。

通过对比计算可以看出，与传统湖泊地下水相互作用方法相比，本发明在相同的湖泊刻画精度和参数设置下，能够得出和传统湖泊地下水相互作用方法相类似的模拟结果，由此可以看出本发明的合理性。通过增加对湖泊的刻画精度，可以看出，随着湖泊刻画精度的增加，传统湖泊地下水相互作用方法含水层垂向必然要增加，从而导致了计算的不稳定性出现。而本发明湖泊刻画精度不受含水层垂向剖分的影响，依然能够正常的计算，从而得到了，与传统湖泊地下水相互作用方法相比，本发明在适用性和稳定性方面具有显著优势。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (4)

1.一种湖泊与地下水非稳定流作用模拟方法，其特征在于，包括以下步骤：

S1、对所有湖泊进行离散化处理及边界连续性处理，得到所有湖泊的各湖泊网格单元；

S2、采集并获取所有湖泊的参数数据；

S3、计算当前湖泊平均水位；

S4、根据当前湖泊平均水位判断当前湖泊的各湖泊网格单元所处的积水状态；

S5、不考虑源汇项建立地下水数值计算矩阵方程；

S6、求解地下水数值计算矩阵方程，得到当前迭代下的地下水位；

S7、根据当前湖泊的各湖泊网格单元所处的积水状态，结合当前迭代下的地下水位计算得到当前湖泊的各湖泊网格单元的水量平衡项；

S8、根据当前湖泊的各湖泊网格单元的水量平衡项计算得到湖泊与地下水的水分交换统计量；

S9、根据湖泊与地下水的水分交换统计量和湖泊参数数据计算得到当前迭代下的湖泊蓄水量；

S10、判断当前迭代下湖泊水位是否收敛，收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值，若是则模拟结束，否则返回步骤S3进入下一次迭代；

所述步骤S2中的所有湖泊的参数数据包括：时段内的降水强度p、时段内的降水产流系数γ、降水入渗补给系数k、时段内的水面蒸发强度e₀、湖泊上游河道汇入量Q_si、湖泊人工取水量W以及湖泊下泄量Q_so；

所述步骤S3中当前湖泊平均水位的计算公式为：

其中

为当前湖泊平均水位，η∈[0,1]为隐式加权因子，

为时段初，即上一时段末的湖泊水位，

为时段末的湖泊水位；

所述步骤S4具体为：

若

则该湖泊网格单元处于完全积水状态；若

则该湖泊网格单元处于完全未积水状态；若

则该湖泊网格单元处于部分积水状态；其中

为单元内湖底高程的最高值，

为单元内湖底高程的最低值；

所述步骤S5中地下水数值计算矩阵方程为：

[A]{h}＝{q}

其中[A]为系数矩阵，{h}为地下水数值计算矩阵，{q}为所有常数项和已知项集合；

所述步骤S6具体为：

将系数-C_m加入到系数矩阵[A]主对角线系数中，将

加入到矩阵方程右端项{q}中，得到当前迭代下的地下水位

其中C_m为该湖泊网格单元处湖底与含水层之间的综合水力传导系数；

所述步骤S7具体为：

若湖泊网格单元处于完全积水状态，当地下水位

高于单元内湖底高程的最高值

时，该湖泊网格单元的水量平衡项为：

其中

和

分别为对应湖泊平均水位

第m个湖泊网格单元积水面积部分地下水渗出排泄到湖泊的流量和湖泊水体渗漏到地下水的流量；

若湖泊网格单元处于完全积水状态，当地下水位

低于单元内湖底高程的最低值

时，该湖泊网格单元的水量平衡项为：

若湖泊网格单元处于完全积水状态，当地下水位

位于单元内湖底高程的最高值

和单元内湖底高程的最低值

之间时，该湖泊网格单元的水量平衡项为：

其中R_a为湖泊网格单元处地下水位之下面积部分所占的面积比例，

为地下水位之下面积部分的湖泊渗漏流量，

为地下水位之上面积部分的湖泊渗漏流量；

若湖泊网格单元处于完全未积水状态，当地下水位

高于单元内湖底高程的最高值

时，该湖泊网格单元的水量平衡项为：

其中

和

分别为对应湖泊平均水位

第m个湖泊网格单元未积水面积部分地下水渗出排泄到湖泊的流量和降水入渗量，E₀为单元上作用的最大潜水蒸发强度，

为该单元的面积；

若湖泊网格单元处于完全未积水状态，当地下水位

低于单元内湖底高程的最低值

时，该湖泊网格单元的水量平衡项为：

其中

为对应湖泊平均水位

第m个湖泊网格单元未积水面积部分的潜水蒸发量，p为时段内的降水强度，k为降水入渗补给系数，E_p为单元上作用的潜水蒸发强度，计算公式为：

其中D_M为潜水蒸发极限埋深，D为实际水位埋深，此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

若湖泊网格单元处于完全未积水状态，当地下水位

位于单元内湖底高程的最高值

和单元内湖底高程的最低值

之间时，该湖泊网格单元的水量平衡项为：

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

高于湖泊平均水位

时，该湖泊网格单元的水量平衡项为：

其中R_a,p为湖泊网格单元积水面积部分占整个网格单元的面积比例；

若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

低于单元内湖底高程的最低值

时，该湖泊网格单元的水量平衡项为：

若湖泊网格单元处于部分积水状态，在积水面积部分，当地下水位

位于湖泊平均水位

和单元内湖底高程的最低值

之间时，该湖泊网格单元的水量平衡项为：

其中R_a,p1为湖泊网格单元中地下水位以下的积水面积部分占整个网格单元的面积比例，R_a,p2为湖泊网格单元中地下水位以上的积水面积部分占整个网格单元的面积比例；

若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

高于单元内湖底高程的最高值

时，该湖泊网格单元的水量平衡项为：

其中R_a,n为湖泊网格单元未积水面积部分占整个网格单元的面积比例；

若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

低于湖泊平均水位

时，该湖泊网格单元的水量平衡项为：

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

若湖泊网格单元处于部分积水状态，在未积水面积部分，当地下水位

位于湖泊平均水位

和单元内湖底高程的最高值

之间时，该湖泊网格单元的水量平衡项为：

其中R_a,n1为湖泊网格单元中低于地下水位的未积水面积部分占整个网格单元的面积比例，R_a,n2为湖泊网格单元中高于地下水位的未积水面积部分占整个网格单元的面积比例；

此时用于计算潜水蒸发强度E_p时的实际水位埋深D取值为：

2.根据权利要求1所述的湖泊与地下水非稳定流作用模拟方法，其特征在于，所述步骤S1具体包括：

S1-1、对所有湖泊进行离散化处理：将湖底所在的含水层网格单元标识为湖泊网格单元，并将湖底高程数据以离散的方式赋予每个湖泊网格单元，使得每个湖泊网格单元都具有其单元面积范围内的平均湖底高程值；如果某个湖泊网格单元不位于第一层，则该湖泊网格单元上方的所有网格单元都将被定义为无效单元；

S1-2、对所有湖泊进行边界连续性处理：假设湖底高程在湖泊网格单元内是倾斜的，具有单元内的最低值和最高值；对于非最上级和非最下级的湖泊网格单元，其湖底高程最低值为本单元湖底平均高程与下一级单元湖底平均高程的中间位置，最高值为本单元湖底平均高程与上一级单元湖底平均高程的中间位置；对于最下级的湖泊网格单元，其湖底高程最低值为其平均湖底高程减去其与上一级单元平均湖底高程差的一半；对于最上级湖泊网格单元，其湖底高程最高值为其平均湖底高程加上其与下一级单元平均湖底高程差的一半。

3.根据权利要求1所述的湖泊与地下水非稳定流作用模拟方法，其特征在于，所述步骤S8中湖泊与地下水的水分交换统计量的计算公式为：

其中M为湖泊网格单元总数，

为时段内湖泊积水区含水层向湖泊的渗出流量，

为时段内湖泊积水区的渗漏流量，

为时段内湖泊未积水区含水层的渗出流量，

为时段内湖泊未积水区的降水入渗量，

为时段内湖泊未积水区的潜水蒸发量。

4.根据权利要求3所述的湖泊与地下水非稳定流作用模拟方法，其特征在于，所述步骤S9中当前迭代下的湖泊蓄水量的计算公式为：

其中

为当前迭代下的湖泊蓄水量，

为前一次迭代下的湖泊蓄水量，p为时段内的降水强度，γ为时段内的降水产流系数，e₀为时段内的水面蒸发强度，

为时段内湖泊的平均水面面积，Δt为当前计算时段，A_T为湖泊网格单元总面积，Q_si为湖泊上游河道汇入量，W为湖泊人工取水量，Q_so为湖泊下泄量。

Images