CN115983158B - 一种地下水模型与二维水动力模型松散耦合的方法 - Google Patents

Abstract

本发明涉及一种地下水模型与二维水动力模型松散耦合的方法，具体步骤包括：1、构建二维水动力模型；2、在云平台使用超级计算机进行水动力模型的并行模拟计算；3、构建地下水模型；4、根据水动力模型结果计算洪水地下水补给量；5、根据地下水模型与水动力模型的时间、空间耦合关系，将洪水地下水补给量输入地下水模型；6、补充地下水模型的其他补给项、排泄项，完成地下水耦合模拟。本发明克服了地下水模型与水动力模型耦合中存在的模拟效率不同、时间步长不一致的问题，建立了水动力‑地下水耦合模型。本发明建立的耦合模型可充分体现在洪水演进过程中地表水对地下水的动态影响，实现对地下水位的准确模拟。

Description

一种地下水模型与二维水动力模型松散耦合的方法

技术领域

本发明属于水利数值模拟领域，具体涉及一种地下水模型与二维水动力模型松散耦合的方法。

背景技术

地下水主要接受地表水补给，地表水的范围、水深常处于充满不确定性的变化过程中，如果忽略这种变化，直接进行地下水模拟，会导致模拟结果与实际情况存在比较大的偏差，因此，需进行地表水与地下水耦合模拟，以提高地下水模拟的准确性。

目前，模拟地表水变化的常用数值模型有水文模型、水动力模型；为提高地表水模型的准确性，水文水动力耦合模型、一维二维水动力耦合模型得到了广泛应用，已有较多技术成果；为提高地下水模拟的准确性，也有许多水文模型与地下水模型耦合的技术成果，但水动力模型与地下水模型进行耦合的相关技术较少。

地下水模型与水动力模型耦合的实现存在两大难点，一是两个模型的模拟效率不同，水动力模型模拟耗时通常以天为单位，而地下水模型模拟耗时通常只需要几分钟，水动力模型完成地下水模型同样时长的模拟耗时太长，阻碍模型耦合的实现；二是两个模型的时间步长不一致，水动力模型的时间步长通常为几秒，难以与地下水模型以天为单位的时间步长进行耦合。

发明内容

针对地下水模型与水动力模型耦合中存在的模拟效率不同、时间步长不一致的问题，本发明的目的在于提供一种地下水模型与二维水动力模型松散耦合的方法，以提高地下水模拟的准确性。

为解决现有技术问题，本发明采取的技术方案为：

一种地下水模型与二维水动力模型松散耦合的方法，包括以下步骤：

S1、构建二维水动力模型，用于模拟地表二维水流运动，反映地表水(洪水)的存在区域及在各个区域的停留时间，即洪水演进过程；

S2、设置水动力模型的并行数，在云平台使用超级计算机进行水动力模型的并行模拟计算；

S3、构建地下水数值模型，用于模拟地下水的变化；

S4、根据地下水模型与水动力模型的时空耦合关系构建耦合模型，并完成对地下水位的准确模拟。

作为改进的是，步骤S1中所述的二维水动力模型是指采用非结构化三角形或四边形网格的水动力模型，具体步骤如下：

S1.1根据实测或遥感数据，分析地表高程及地表形态，综合确定模型的模拟区域；

S1.2将模型的模拟区域剖分成非结构化网格，并将河道处的网格加密；

S1.3将地表高程数据及模型的边界条件输入模型，即可。

作为改进的是，步骤S2中在云平台使用超级计算机进行水动力模型的并行模拟计算，具体步骤如下：

S2.1通过远程控制软件，连接云平台的超级计算机，并在超级计算机上安装水动力模拟软件构建二维水动力模型；

S2.2根据云平台的计算机性能以及模型的网格数量，综合确定模型的并行数；

S2.3将并行数输入水动力模型后，在超级计算机上完成水动力模型的并行模拟计算。

作为改进的是，步骤S3中构建地下水数值模型的步骤如下：

S3.1根据地形、地质及地下水流场数据综合确定模型的侧向边界及各边界的类型；

S3.2根据模拟区域的水文地质图、钻孔资料对含水层进行概化；

S3.3根据气象资料、水文地质资料、灌溉资料对模型的补给项、排泄项进行概化，补给项、排泄项的最小时间步长为天。

作为改进的是，步骤S4构建耦合模型的具体操作如下：

S4.1分析洪水下渗过程，将下渗过程简化，构建由洪水淹没时长、淹没面积及下渗率计算洪水地下水补给量的计算公式：

W_补给G＝f(t,f_c(t),S(t))

其中，W_补给G为整场洪水持续期间的洪水地下水补给量，t为时间，f_c(t)为随时间变化的下渗率，S(t)为随时间变化的洪水淹没面积；

S4.2分析二维水动力模型的结果文件，直接从中读取洪水持续总天数K_max；

S4.3分析二维水动力模型的结果文件，按海伦公式计算洪水淹没面积，当二维水动力模型的网格为非结构化四边形时，网格面积计算公式如下：

若有2个顶点坐标的水深大于0，则四边形网格处于淹没状态；

当二维水动力模型的网格为非结构化三角形网格时，面积计算公式如下：

对于三角形网格，若有1个顶点坐标的水深大于0，则认为该网格处于淹没状态；其中，i表示一个被淹没的网格，s_i(t)为t时刻第i个网格的面积，p_i(t)、a_i(t)、b_i(t)、c_i(t)、d_i(t)为第i个网格t时刻的中间量，

为t时刻第i个网格的顶点坐标；

t时刻的总淹没面积由下式计算：

其中，N为t时刻处于淹没状态的网格总数；

S4.4以某一时刻T的淹没情况代表该天的淹没情况，则可由步骤S4.1、步骤S4.2、步骤S4.3得到下式，以计算模拟期内全域逐日的洪水地下水补给量，建立地下水模型与二维水动力模型的时间耦合关系：

其中，K表示洪水持续天数的顺序，

为第K天的洪水地下水补给量；f_c ^K(T)为第K天T时刻的下渗率，单位为mm/d；S^K(T)为第K天的洪水淹没面积；

S4.5对于某一模拟区域，f_c(T)可由经验值或现场试验求得，则由步骤S4.1、步骤S4.2、步骤S4.4得到整场洪水持续期间的洪水地下水补给量计算公式：

其中，W_补给G为整场洪水持续期间的洪水地下水补给量，K_max为洪水持续总天数，

为洪水地下水补给量的最大值，等于包气带的体积乘以给水度；

S4.6根据步骤S1中的模型，选来水量最大的几年按步骤S2中的方法进行洪水模拟，根据模拟结果，选出淹没范围最广的一年，分析其结果文件；提取结果文件中所有节点每天T时刻(与S4.4步骤中所取时刻相同)的水深数据，将水深大于零、小于零的节点淹没属性分别设置为1、0，用ArcGIS的空间插值工具，按节点的淹没属性进行插值，得到K_max个逐日淹没范围栅格文件；用空间分析工具中的Raster calculator工具，将K_max个淹没范围栅格文件相加，得到一个淹没天数栅格文件；

S4.7使用空间分析工具中的Reclassify工具对淹没天数栅格文件进行栅格重分类，将二维水动力模型分为n个区域，并使用ArcGIS统计不同区域的面积，分区方法如下：

按淹没天数占总淹没天数的比例，共设置n个区间，各区间的比例分别为0～5g₁％、5g₁％～5g₂％、...、5g_n-2％～5g_n-1％、5g_n-1％～100％，与之对应的水动力模型的各区域面积分别为S₁、S₂、…、S_n-1、S_n；

其中，n＝2,…,20；g₀＝0；g₁＝1,2,…,19；g₂＝g₁+1,g₁+2,…,19；g_n＝g_n-1+1,g_n-1+2,…,19；若g_n取值加1，而S_n增幅小于5％，则g_n取增加前的值；

S4.8根据步骤S4.7中各分区的面积及淹没天数比例，计算各分区的平均淹没天数，公式如下：

其中，K_n为第n个区域的平均淹没天数，p为第n个区域淹没天数占总天数的比例，L、H分别为第n个区域p的下限、上限，S_n(p)为第n个区域中淹没天数占总天数比例为p的面积，q_n为第n个区域平均淹没天数占全域总淹没天数的比例，

S4.9根据步骤S4.5中的整场洪水持续期间的洪水地下水补给量、步骤S4.7中各区域的面积、步骤S4.8中各区域平均淹没天数，得到各区域洪水地下水补给量的计算公式：

其中，W_n为第n个区域的洪水地下水补给量；

S4.10根据地下水模型与二维水动力模型的坐标信息建立地下水模型与水动力模型的空间耦合关系；

S4.11将S4.9步骤的结果转化成与地下水模型时间步长一致的洪水地下水补给强度，以其作为纽带，按S4.10步骤的空间耦合关系将其输入地下水模型，实现地下水模型与水动力模型的耦合；

S4.12除步骤S4.11中的补给量外，将其它补给项、排泄项输入地下水模型，完成地下水位的模拟。

有益效果：

与现有技术相比，本发明一种地下水模型与二维水动力模型松散耦合的方法，通过在云平台的超级计算机上进行水动力模型的并行模拟计算，解决了水动力模型计算效率低的问题，为完成长时间的洪水演进模拟提供了可能，最终为地下水模型提供充足且符合实际时空变化的洪水地下水补给量打下了基础；

通过本发明提供的耦合方法，实现了地下水模型与水动力模型的耦合计算，提高了地下水模拟的准确性；

地下水的准确模拟可为预测植被生长情况提供理论依据，进而提高水资源综合管理能力、为生态环境保护提供有效支撑。

附图说明

图1为地下水模型与二维水动力模型的耦合原理；

图2为实际模拟区域(额尔齐斯河中游河谷区)的二维水动力模型计算范围、网格剖分情况及模型边界条件；

图3为实际模拟区域(额尔齐斯河中游河谷区)的地下水数值模型模拟范围及边界类型；

图4为二维水动力模型某个时刻的可视化模拟结果；

图5为实际模拟区域(额尔齐斯河中游河谷区)二维水动力模型的分区结果；

图6为实际模拟区域(额尔齐斯河中游河谷区)地下水模型与二维水动力模型空间耦合示意；

图7为实际模拟区域(额尔齐斯河中游河谷区)地下水模型的率定结果；

图8(a)为实际模拟区域(额尔齐斯河中游河谷区)受洪水淹没的部分地下水观测井的模拟值与观测值对比；

图8(b)为实际模拟区域(额尔齐斯河中游河谷区)不受洪水淹没的部分地下水观测井的模拟值与观测值对比；

图9为实际模拟区域(额尔齐斯河中游河谷区)在给定未来10年来水条件下的地下水埋深情况。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。

图1是本发明一种地下水模型与二维水动力模型的松散耦合流程图：

S1，构建二维水动力模型，用于模拟地表二维水流运动，具体用于模拟洪水的淹没范围及持续时间，反应洪水的存在区域以及在各个区域的停留时间，即洪水的演进过程；

S2，设置水动力模型的并行数，在云平台的计算机上进行水动力模型的并行模拟计算；

S3，构建地下水数值模型，用于模拟地下水的变化；

S4、根据地下水模型与水动力模型的时空耦合关系构建耦合模型，并完成对地下水位的准确模拟。

实施例

计算区域以额尔齐斯河中游河漫滩区域为例，二维水动力模型采用非结构化的三角形网格。

图2是本发明的实际二维水动力模型的计算区域示意图及水动力模型的网格剖分情况。

图3为实际模拟区域(额尔齐斯河中游河谷区)的地下水数值模型模拟范围及边界类型。

图4是本发明的二维水动力模型某个时刻的可视化模拟结果。据其可以获取各时刻的洪水淹没范围及洪水的持续时间：

S4.1，额尔齐斯河中游河漫滩地区位于我国西北干旱地区，结合该区域的气象、土壤、洪水特性分析该区域的洪水下渗过程，可将该区域的洪水下渗过程视为供水强度充分大条件下的下渗。即下渗的洪水，一部分被土壤吸收，另一部分补给地下水，并且当土壤吸水饱和后，土壤不再吸水，洪水将按照固定的下渗率补给地下水，此时的下渗率称为稳定下渗率。将本实施例的洪水下渗过程进行简化，即在洪水消失以前，只要地面积水存在，下渗率不发生改变，均为稳定下渗率，且从洪水下渗开始，便以稳定下渗率补给地下水。被土壤所吸收的下渗的洪水与本发明无关，因此只考虑洪水地下水补给量。则整个洪水持续过程的洪水地下水补给量为：

其中，W_补给G为洪水地下水补给量，f_c为稳定下渗率，t为时间，t₁、t₂分别为洪水下渗开始时刻、结束时刻，S(t)为随时间变化的洪水淹没面积。

S4.2，分析二维水动力模型的结果文件，从模拟结果中直接读取洪水持续总天数K_max。

S4.3，分析二维水动力模型的结果文件，按海伦公式计算洪水淹没面积。

如图2所示，实施本发明的二维水动力模型采用非结构化的三角形网格，由海伦公式计算处于淹没状态的某个三角形网格面积，计算公式如下所示：

其中，i表示一个被淹没的网格，s_i(t)为t时刻第i个网格的面积，p_i(t)、a_i(t)、b_i(t)、c_i(t)为第i个网格t时刻的中间量，

为t时刻第i个网格的顶点坐标；

t时刻的总淹没面积由下式计算：

其中，N为t时刻处于淹没状态的网格总数；

S4.4，将零点的洪水状况作为全天的洪水状况，则该天洪水地下水补给量由下式计算：

其中，K表示洪水持续天数的顺序，

为第K天的洪水地下水补给量；f_c为稳定下渗率，单位为mm/d；S^K(0)为第K天的洪水淹没面积；

整场洪水持续期内的洪水地下水补给总量由下式计算：

其中，W_补给G为整场洪水持续期间的洪水地下水补给量，K_max为洪水持续总天数，

为洪水下渗补给量的最大值，等于包气带的体积乘以给水度；

S5.1,本发明根据淹没时长占总时长的0～5％、5％～80％、80％～100％将二维水动力模型模拟范围分别分为区域1、区域2、区域3，结果如图5所示，三个区域的淹没时长分别为总淹没时长的2％、50％、95％，三个区域的面积分别为S₁、S₂、S₃。

S5.2，由水动力模型各区域的面积及淹没时长比例得到各区域的洪水下渗补给量分别为：

其中，W₁、W₂、W₃分别是区域1、区域2、区域3的洪水地下水补给量。

根据坐标信息，将地下水模型模拟范围划分为与水动力模型三个分区相对应的区域，并将计算的洪水下渗补给量按照空间耦合关系输入地下水模型，实现了地下水模型与二维水动力模型的空间耦合，空间耦合示意图如图6所示。

图7是实施本发明的地下水模型率定结果，图中7zk-1、7zk-2、7zk-3、7zk-4、10zk-1、10zk-2、10zk-3、10zk-4是观测井编号，由图7可知，各观测井地下水位的观测值、模拟值基本吻合，地下水模型模拟精度较高。

图8(a)显示地下水位模拟结果能够较准确模拟7zk-2观测井处受洪水淹没造成的地下水位的上升、下降。

图8(b)显示地下水位模拟结果能够较准确模拟7zk-2观测井处不受洪水淹没时的地下位变化情况。

图9是使用本发明的地下水-水动力耦合模型模拟对给定未来10年来水条件下的地下水埋深情况，从图中可以看出，在该给定来水条件下，研究区地下水埋深为2～4m，该结果可为研究区植被保护、水资源综合管理提供依据。

Claims (4)

1.一种地下水模型与二维水动力模型松散耦合的方法，其特征在于，包括以下步骤：

S1、构建二维水动力模型，用于模拟地表二维水流运动；

S2、设置水动力模型的并行数，在云平台使用超级计算机进行水动力模型的并行模拟计算；

S3、构建地下水数值模型，用于模拟地下水的变化；

S4、根据地下水模型与水动力模型的时空耦合关系构建耦合模型，并完成对地下水位的准确模拟，其中，步骤S4构建耦合模型的具体操作如下：

S4.1分析洪水下渗过程，将下渗过程简化，构建由洪水淹没时长、淹没面积及下渗率计算洪水地下水补给量的计算公式：

W_补给G＝f(t,f_c(t),S(t))

其中，W_补给G为整场洪水持续期间的洪水地下水补给量，t为时间，f_c(t)为随时间变化的下渗率，S(t)为随时间变化的洪水淹没面积；

S4.2分析二维水动力模型的结果文件，直接从中读取洪水持续总天数K_max；

S4.3分析二维水动力模型的结果文件，按海伦公式计算洪水淹没面积，当二维水动力模型的网格为非结构化四边形时，网格面积计算公式如下：

若有2个顶点坐标的水深大于0，则四边形网格处于淹没状态；

当二维水动力模型的网格为非结构化三角形网格时，面积计算公式如下：

对于三角形网格，若有1个顶点坐标的水深大于0，则认为该网格处于淹没状态；其中，i表示一个被淹没的网格，s_i(t)为t时刻第i个网格的面积，p_i(t)、a_i(t)、b_i(t)、c_i(t)、d_i(t)为第i个网格t时刻的中间量，

为t时刻第i个网格的顶点坐标；

t时刻的总淹没面积由下式计算：

其中，N为t时刻处于淹没状态的网格总数；

S4.4以某一时刻T的淹没情况代表该天的淹没情况，则由步骤S4.1、步骤S4.2、步骤S4.3得到下式，以计算模拟期内全域逐日的洪水地下水补给量，建立地下水模型与二维水动力模型的时间耦合关系：

其中，K表示洪水持续天数的顺序，

为第K天的洪水地下水补给量；f_c ^K(T)为第K天T时刻的下渗率，单位为mm/d；S^K(T)为第K天的洪水淹没面积；

S4.5对于某一模拟区域，f_c(T)由经验值或现场试验求得，则由步骤S4.1、步骤S4.2、步骤S4.4得到整场洪水持续期间的洪水地下水补给量计算公式：

其中，W_补给G为整场洪水持续期间的洪水地下水补给量，K_max为洪水持续总天数，

为洪水地下水补给量的最大值，等于包气带的体积乘以给水度；

S4.6根据步骤S1中的模型，选来水量最大的几年按步骤S2中的方法进行洪水模拟，根据模拟结果，选出淹没范围最广的一年，分析其结果文件；提取结果文件中所有节点每天T时刻的水深数据，其中，每天T时刻与S4.4步骤中所取时刻相同，将水深大于零、小于零的节点淹没属性分别设置为1、0，用ArcGIS的空间插值工具，按节点的淹没属性进行插值，得到K_max个逐日淹没范围栅格文件；用空间分析工具中的Raster calculator工具，将K_max个淹没范围栅格文件相加，得到一个淹没天数栅格文件；

S4.7使用空间分析工具中的Reclassify工具对淹没天数栅格文件进行栅格重分类，将二维水动力模型分为n个区域，并使用ArcGIS统计不同区域的面积，分区方法如下：

按淹没天数占总淹没天数的比例，共设置n个区间，各区间的比例分别为0～5g₁％、5g₁％～5g₂％、...、5g_n-2％～5g_n-1％、5g_n-1％～100％，与之对应的水动力模型的各区域面积分别为S₁、S₂、…、S_n-1、S_n；

其中，n＝2,…,20；g₀＝0；g₁＝1,2,…,19；g₂＝g₁+1,g₁+2,…,19；g_n＝g_n-1+1,g_n-1+2,…,19；若g_n取值加1，而S_n增幅小于5％，则g_n取增加前的值；

S4.8根据步骤S4.7中各分区的面积及淹没天数比例，计算各分区的平均淹没天数，公式如下：

其中，K_n为第n个区域的平均淹没天数，p为第n个区域淹没天数占总天数的比例，L、H分别为第n个区域p的下限、上限，S_n(p)为第n个区域中淹没天数占总天数比例为p的面积，q_n为第n个区域平均淹没天数占全域总淹没天数的比例，

S4.9根据步骤S4.5中的整场洪水持续期间的洪水地下水补给量、步骤S4.7中各区域的面积、步骤S4.8中各区域平均淹没天数，得到各区域洪水地下水补给量的计算公式：

其中，W_n为第n个区域的洪水地下水补给量；

S4.10根据地下水模型与二维水动力模型的坐标信息建立地下水模型与水动力模型的空间耦合关系；

S4.11将S4.9步骤的结果转化成与地下水模型时间步长一致的洪水地下水补给强度，以其作为纽带，按S4.10步骤的空间耦合关系将其输入地下水模型，实现地下水模型与水动力模型的耦合；

S4.12除步骤S4.11中的补给量外，将其它补给项、排泄项输入地下水模型，完成地下水位的模拟。

2.根据权利要求1所述的一种地下水模型与二维水动力模型松散耦合的方法，其特征在于，步骤S1中所述的二维水动力模型是指采用非结构化三角形或四边形网格的水动力模型，具体步骤如下：

S1.1根据实测或遥感数据，分析地表高程及地表形态，综合确定模型的模拟区域；

S1.2将模型的模拟区域剖分成非结构化网格，并将河道处的网格加密；

S1.3将地表高程数据及模型的边界条件输入模型。

3.根据权利要求1所述的一种地下水模型与二维水动力模型松散耦合的方法，其特征在于，步骤S2中在云平台使用超级计算机进行水动力模型的并行模拟计算，具体步骤如下：

S2.1通过远程控制软件，连接云平台的超级计算机，并在超级计算机上安装水动力模拟软件，构建二维水动力模型；

S2.2根据云平台的计算机性能以及模型的网格数量，综合确定模型的并行数；

S2.3将并行数输入水动力模型后，在超级计算机上完成水动力模型的并行模拟计算。

4.根据权利要求1所述的一种地下水模型与二维水动力模型松散耦合的方法，其特征在于，步骤S3中构建地下水数值模型的步骤如下：

S3.1根据地形、地质及地下水流场数据综合确定模型的侧向边界及各边界的类型；

S3.2根据模拟区域的水文地质图、钻孔资料对含水层进行概化；

S3.3根据气象资料、水文地质资料、灌溉资料对模型的补给项、排泄项进行概化，补给项、排泄项的最小时间步长为天。

Images