CN108763797B - 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 - Google Patents

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所示为湖泊为干涸状态时两种方法湖泊水位模拟结果对比示意图。

具体实施方式

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

本发明实施例提供了一种基于地下水模型的湖泊与地下水稳定流作用模拟方法，包括以下步骤S1-S9：

步骤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、根据上一次迭代下的湖泊水位判断各湖泊网格单元所处的积水状态。

若

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

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

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

为上一次迭代时段末的湖泊水位，

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

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

S4、建立地下水数值计算矩阵方程，本发明实施例中，建立地下水数值计算矩阵方程时不考虑源汇项，具体方程为：

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

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

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

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

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

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

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

其中η∈[0,1]为隐式加权因子，

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

为本次迭代时段末的湖泊水位。

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

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

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

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

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

其中

和

分别为对应湖泊平均水位

第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取值为：

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

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

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

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

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

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

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

S8、根据湖泊与地下水的水分交换统计量和湖泊参数数据更新当前迭代下的湖泊水位。

如图16所示，将湖泊各项源汇项表达为湖泊水位

的函数，可得待求解的湖泊稳定流水量平衡控制方程：

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

其中

为时段内湖泊水位为

时，湖泊积水区的水面降水通量，p为时段内的降水强度，A_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)得到湖泊水量平衡控制方程的具体公式为：

以上湖泊稳定流水量平衡控制方程可采用牛顿迭代法对湖泊水位进行求解：

其中

为上一次迭代的湖泊水位，

为本次迭代的湖泊水位，F′(·)为F(·)的导数，则公式(27)只需确定F′(·)的函数形式即可进行迭代求解。

对公式(26)求导得：

公式(28)中A_P′的确定方法如下：

其中

表示湖泊水位-积水面积线性相关曲线上对应

时的折线斜率。

的计算公式为：

其中，A为完全积水状态的湖泊网格单元中，单元处地下水位高于湖泊水位的湖泊网格单元的总数量；C_a为上述A个湖泊网格单元中第a个处湖底与含水层之间的综合水力传导系数；B为部分积水状态的湖泊网格单元中，单元处地下水位高于湖泊水位的湖泊网格单元的总数量；C_b为上述B个湖泊网格单元中第b个处湖底与含水层之间的综合水力传导系数；

和

分别为第b个湖泊网格单元的单元内湖底最高和最低高程。

的计算公式为：

其中，D为完全积水状态的湖泊网格单元中，单元处地下水位低于湖泊水位的湖泊网格单元的总数量；C_d为上述D个湖泊网格单元中第d个处湖底与含水层之间的综合水力传导系数；E为部分积水状态的湖泊网格单元中，单元处地下水位低于湖泊水位的湖泊网格单元的总数量；C_e为上述E个湖泊网格单元中第e个处湖底与含水层之间的综合水力传导系数；

和

分别为第e个湖泊网格单元的单元内湖底最高和最低高程。

的计算公式为：

综合公式(28)～(32)得：

由于A与D之和为完全积水状态的湖泊网格单元的总数量，B与E之和为部分积水状态的湖泊网格单元的总数量，则公式(33)可进一步化简为：

其中C_s为第s个完全积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数，S为完全积水状态的湖泊网格单元的总数量，C_f为第f个部分积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数，F为部分积水状态的湖泊网格单元的总数量，

和

分别为第f个部分积水状态的湖泊网格单元的单元内湖底最高和最低高程。

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

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

本发明与传统湖泊地下水相互作用方法根据各自方法属性，对湖泊网格的剖分如图18～20所示，图18中A和B为模拟过程中的观测点，阴影部分表示湖泊单元，图20中阴影部分表示湖泊单元。其他所有参数和模型驱动数据均保持一致，考虑到发明中湖泊水面蒸发是影响湖泊水位主要因素，继续尝试分别增加和减少湖泊的水面蒸发速率，以期对湖泊的不同蓄水状态进行全面的模拟和验证，对比两种方法的模拟计算结果(如图21～23所示)，可以看出两种方法模拟结果差别不大，从而证明了本发明方法的合理性。但在模拟过程中发现，LAK3并不是全局收敛，只有给其合适的初始值很靠近真实值时，迭代才会收敛，这一缺陷将很大程度上限制LAK3在实际工作的应用，而本发明则可以全局收敛。

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

Claims (6)

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

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

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

S3、根据上一次迭代下的湖泊水位判断各湖泊网格单元所处的积水状态；

S4、建立地下水数值计算矩阵方程；

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

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

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

S8、根据湖泊与地下水的水分交换统计量和湖泊参数数据更新当前迭代下的湖泊水位；

S9、判断当前迭代下湖泊水位是否收敛，若是则模拟结束，否则返回步骤S3进入下一次迭代；

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

所述步骤S3具体为：

若

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

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

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

为上一次迭代时段末的湖泊水位，

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

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

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

[A]{h}＝{q}

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

所述步骤S5具体为：

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

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

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

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

其中η∈[0,1]为隐式加权因子，

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

为本次迭代时段末的湖泊水位。

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

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

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

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

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

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

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

其中

和

分别为对应湖泊平均水位

第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取值为：

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

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

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

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

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

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

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

5.根据权利要求4所述的湖泊与地下水稳定流作用模拟方法，其特征在于，所述步骤S8中更新当前迭代下的湖泊水位的公式为：

其中

为上一次迭代的湖泊水位，

为本次迭代的湖泊水位，F(·)表示湖泊水量平衡控制方程，具体公式为：

其中，p为时段内的降水强度，γ为时段内的降水产流系数，e₀为时段内的水面蒸发强度，A_P为时段内湖泊的平均水面面积，A_T为湖泊网格单元总面积，Q_si为湖泊上游河道汇入量，W为湖泊人工取水量，Q_so为湖泊下泄量；

F′(·)为F(·)的导数，具体公式为：

其中

表示湖泊水位-积水面积线性相关曲线上对应

时的折线斜率，C_s为第s个完全积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数，S为完全积水状态的湖泊网格单元的总数量，C_f为第f个部分积水状态的湖泊网格单元处湖底与含水层之间的综合水力传导系数，F为部分积水状态的湖泊网格单元的总数量，

和

分别为第f个部分积水状态的湖泊网格单元的单元内湖底最高和最低高程。

6.根据权利要求5所述的湖泊与地下水稳定流作用模拟方法，其特征在于，所述步骤S9中收敛条件为当前迭代下湖泊水位与上一次迭代下湖泊水位结果之差小于设定的收敛阈值。

Images