CN109376433B - 基于土壤非饱和水和地下水耦合的区域水流运动模拟方法 - Google Patents

Abstract

本发明提供一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，包括以下步骤：对于研究区域范围内每个有限元节点P，均采用以下方式处理，由此得到与每个有限元节点P在空间上对应的栅格；计算得到有限元节点P的降雨入渗量；建立一维土壤非饱和流模拟模型；计算有限元节点P的影响面积A_P；建立三维地下水流模型，实现对地下水流运动状态的模拟。优点为：(1)将研究区域的遥感栅格数据和有限元节点实现了空间上的关联，保证了后续采用有限元方法进行地下水流模拟的精度；(2)一维土壤非饱和流模拟模型和三维地下水流模型进行了耦合，保证了对地下水流运动状态的模拟精度。

Description

基于土壤非饱和水和地下水耦合的区域水流运动模拟方法

技术领域

本发明属于区域水流运动模拟技术领域，具体涉及一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法。

背景技术

降雨、蒸发、土壤入渗等地表水文过程和地下水流运动在一定地形、地质和气候条件下是相互作用、具有内在联系的有机整体。因此，在对地下水流运动数值模拟时，需要以地表水文数据作为输入驱动，才能够保证地下水流运动数值模拟的精确性。

然而，在实现本发明的过程中，发明人发现，现有技术至少存在以下问题：

由于地表水文数据主要通过卫星遥感获得，而卫星遥感数据是以栅格形式存储的，因此，对于采用基于有限元方法的地下水模型，当其进行地下水流运动数值模拟时，栅格形式的卫星遥感数据无法直接输入到地下水模型中，从而不利于对地下水流运动进行数值模拟。

发明内容

针对现有技术存在的缺陷，本发明提供一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，包括以下步骤：

步骤1，获取研究区域的遥感栅格数据；其中，所述研究区域的遥感栅格数据对应M行N列个栅格；每个栅格的边长为s；基于所述遥感栅格数据，得到每个栅格顶点的降雨入渗量值；

步骤2，对研究区域的土壤和地下水层进行有限元网格剖分，水平方向上得到多个三角形形状的有限元网格，进而也得到多个有限元节点，垂直方向上按不同地层介质的厚度分层；

步骤3，将所述研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中；对于每个有限元节点P，均采用以下方式处理，由此得到与每个有限元节点P在空间上对应的栅格，其中，有限元节点P在空间上对应的栅格含义为：当将研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中时，有限元节点P落入到某个栅格中，即为与有限元节点P在空间上对应的栅格；

步骤3.1，在x-y坐标系中，对于M行N列的栅格，其左下角的坐标为(x₀，y₀)；

假设某个有限元节点P的坐标为P(x_P，y_P),则其在空间上对应的栅格W的行数h₁和列数h₂通过以下公式得出：

步骤3.2，在定位到与有限元节点P在空间上对应的栅格W后，栅格W的左下角、右下角、右上角和左上角这四个顶点分别记为：W₁、W₂、W₃和W₄；通过下式分别计算得到W₁、W₂、W₃和W₄的坐标：

W₁的坐标为：(x₀+s(h₁-1)，y₀+s(h₂-1))；

W₂的坐标为：(x₀+sh₁)，y₀+s(h₂-1))；

W₃的坐标为：(x₀+sh₁，y₀+sh₂))；

W₄的坐标为：(x₀+s(h₁-1)，y₀+sh₂))；

步骤4，采用下述方法计算得到有限元节点P的降雨入渗量q^(p)；

步骤4.1，在得到有限元节点P在空间上对应的栅格W的四个顶点W₁、W₂、W₃和W₄的坐标后，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的降雨入渗量分别为q₁，q₂，q₃，q₄；

步骤4.2，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的基函数分别为：N₁(x_P，y_P)、N₂(x_P，y_P)、N₃(x_P，y_P)和N₄(x_P，y_P)；即：栅格顶点的基函数与有限元节点P的坐标相关；

步骤4.3，采用下述的线性插值函数计算得到有限元节点P的降雨入渗量q^(p)的值：

q^(p)＝q₁·N₁(x_P，y_P)+q₂·N₂(x_P，y_P)+q₃·N₃(x_P，y_P)+q₄·N₄(x_P，y_P)

；

步骤5，基于以下的一维土壤非饱和流运动方程，建立一维土壤非饱和流模拟模型：

其中：

n_r为土壤孔隙度，S为土壤水的饱和度；

t为时间；

z为坐标，代表土壤深度；

q_s为源项，等于步骤4得到的每个有限元节点P的降雨入渗量；

q_z为达西流速，为待求值；

运行所述一维土壤非饱和流模拟模型，计算得到与每个有限元节点P对应的达西流速值；

步骤6，假设有限元节点P为u个有限元网格的共享节点，u个有限元网格分别为V₁、V₂,…,Vu；则采用下式计算有限元节点P的影响面积A_P:

其中：N^K为有限元节点P在有限元网格V^K上的基函数，K＝1、2,…,u；

Ω_K为有限元节点P在有限元网格V^K上的积分空间；

采用下式计算得到有限元节点P对应的地下水补给量Q：

Q＝q_z*A_P

步骤7，基于以下的三维地下水流运动方程，建立三维地下水流模型：

其中：

h为地下水水头，为因变量；

k_x、k_y、k_z分别为x，y，z方向的水力传导系数；

S₀为含水层储水系数；

t为时间；

W_i为源汇项，等于步骤7计算得到的地下水补给量Q；

运行所述三维地下水流模型，实现对地下水流运动状态的模拟。

优选的，步骤1中，采用以下方法得到每个栅格的顶点的降雨入渗量值：

通过所述遥感栅格数据，得到与栅格顶点位置对应的降雨量P_a、太阳辐射能R_a、大气的平均气温T_a，最高气温T_max和最低气温T_min，则采用下式得到栅格的顶点的降雨入渗量值q：

q＝P_a-E；

其中：

E为蒸发潜能；是大气温度和太阳辐射能的函数；

κ，λ，β均为经验公式常数，分别取值2.3*10^-3，2.45和17.8。

本发明提供的基于土壤非饱和水和地下水耦合的区域水流运动模拟方法具有以下优点：

(1)将研究区域的遥感栅格数据和有限元节点实现了空间上的关联，并通过本发明提供的映射投影方式，能够快速搜索定位到与有限元节点在空间上对应的栅格，从而将栅格遥感数据传递到有限元节点，作为有限元节点的输入驱动数据，保证了后续采用有限元方法进行地下水流模拟的精度；

(2)一维土壤非饱和流模拟模型和三维地下水流模型进行了耦合，一维土壤非饱和流模拟模型得到的达西流速与影响面积的乘积，得到有限元节点P对应的地下水补给量Q，直接作为三维地下水流模型的上边界条件，保证了对地下水流的地下水水头运动状态的模拟精度。

附图说明

图1为本发明提供的基于土壤非饱和水和地下水耦合的区域水流运动模拟方法的流程示意图；

图2为本发明提供的定位到与有限元节点P对应的栅格的原理图；

图3为本发明提供的计算有限元节点P的影响面积A_P的原理图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

参考图1，本发明提供一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，包括以下步骤：

步骤1，获取研究区域的遥感栅格数据；其中，所述研究区域的遥感栅格数据是指通过卫星遥感获得的数据，遥感栅格数据对应M行N列个栅格，即：研究区域栅格化为M行N列；每个栅格的边长为s；基于所述遥感栅格数据，得到每个栅格顶点的降雨入渗量值；

其中，步骤1中，采用以下方法得到每个栅格的顶点的降雨入渗量值：

通过所述遥感栅格数据，得到与栅格顶点位置对应的降雨量P_a、太阳辐射能R_a、大气的平均气温T_a，最高气温T_max和最低气温T_min，则采用下式得到栅格的顶点的降雨入渗量值q：

q＝P_a-E；

其中：

E为蒸发潜能；是大气温度和太阳辐射能的函数；

κ，λ，β均为经验公式常数，分别取值2.3*10^-3，2.45和17.8。

步骤2，对研究区域的土壤和地下水层进行有限元网格剖分，水平方向上得到多个三角形形状的有限元网格，进而也得到多个有限元节点，垂直方向上按不同地层介质的厚度分层；有限元节点就是指有限元网格的顶点。

步骤3，将所述研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中；对于每个有限元节点P，均采用以下方式处理，由此得到与每个有限元节点P在空间上对应的栅格，其中，有限元节点P在空间上对应的栅格含义为：当将研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中时，有限元节点P落入到某个栅格中，即为与有限元节点P在空间上对应的栅格；

对于待研究的流域，有限元节点P及其对应的栅格实际即为：沿有限元节点P做一条垂直线，有限元节点P及其对应的栅格即位于同一垂直线上。

步骤3.1，在x-y坐标系中，对于M行N列的栅格，其左下角的坐标为(x₀，y₀)；参考图2，为定位到与有限元节点P对应的栅格的原理图；在图2中，其左下角的坐标即为点D的坐标。

假设某个有限元节点P的坐标为P(x_P，y_P),则其在空间上对应的栅格W的行数h₁和列数h₂通过以下公式得出：

参考图2，有限元节点P对应的栅格W的行数h₁＝1+1＝2；有限元节点P对应的栅格W的列数h₂＝1+1＝2。也就是说，以栅格中的点D为参考时，有限元节点P落入第二行第二列的栅格中。

步骤3.2，在定位到与有限元节点P在空间上对应的栅格W后，栅格W的左下角、右下角、右上角和左上角这四个顶点分别记为：W₁、W₂、W₃和W₄；W₁、W₂、W₃和W₄的相对位置参考图2。通过下式分别计算得到W₁、W₂、W₃和W₄的坐标：

W₁的坐标为：(x₀+s(h₁-1)，y₀+s(h₂-1))；

W₂的坐标为：(x₀+sh₁)，y₀+s(h₂-1))；

W₃的坐标为：(x₀+sh₁，y₀+sh₂))；

W₄的坐标为：(x₀+s(h₁-1)，y₀+sh₂)；

步骤4，采用下述方法计算得到有限元节点P的降雨入渗量q^(p)；

步骤4.1，在得到有限元节点P在空间上对应的栅格W的四个顶点W₁、W₂、W₃和W₄的坐标后，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的降雨入渗量分别为q₁，q₂，q₃，q₄；

步骤4.2，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的基函数分别为：N₁(x_P，y_P)、N₂(x_P，y_P)、N₃(x_P，y_P)和N₄(x_P，y_P)；即：栅格顶点的基函数与有限元节点P的坐标相关；

步骤4.3，采用下述的线性插值函数计算得到有限元节点P的降雨入渗量q^(p)的值：

q^(p)＝q₁·N₁(x_P，y_P)+q₂·N₂(x_P，y_P)+q₃·N₃(x_P，y_P)+q₄·N₄(x_P，y_P)

由于栅格也具有一定的长度范围，如果简单的将栅格的中心点坐标作为落入该栅格的有限元节点P的降雨入渗量q^(p)的值，会带来以下问题：如果有两个不同的有限元节点P落入某个栅格的不同位置，将会导致这两个不同的有限元节点P的降雨入渗量的值相同，这显然与实际情况不符，降低了有限元节点P的降雨入渗量的计算精度。因此，本发明引入栅格顶点的基函数，在计算有限元节点P的降雨入渗量时，通过引入的基函数，考虑了有限元节点P在其所落入的栅格的位置，从而提高了有限元节点P的降雨入渗量的计算精度。

步骤5，基于以下的一维土壤非饱和流运动方程，建立一维土壤非饱和流模拟模型：

其中：

n_r为土壤孔隙度，S为土壤水的饱和度；

t为时间；

z为坐标，代表土壤深度；

q_s为源项，等于步骤4得到的每个有限元节点P的降雨入渗量；

q_z为达西流速(Darcy’s flow flux)，为待求值；

运行所述一维土壤非饱和流模拟模型，计算得到与每个有限元节点P对应的达西流速值；

步骤6，假设有限元节点P为u个有限元网格的共享节点，u个有限元网格分别为V₁、V₂,…,Vu；则采用下式计算有限元节点P的影响面积A_P:

其中：N^K为有限元节点P在有限元网格V^K上的基函数，K＝1、2,…,u；

Ω_K为有限元节点P在有限元网格V^K上的积分空间；

参考图3，为计算有限元节点P的影响面积A_P的原理图。在图3中，有限元节点P为5个有限元网格的共享节点，5个有限元网格分别为V₁、V₂,…,V₅；

采用下式计算得到有限元节点P对应的地下水补给量Q：

Q＝q_z*A_P

(2)

该方法基于基函数定义了有限元节点在研究范围内的影响面积，使地下水补给量在空间上的表达式(2)连续、可导，有利于地下水流微分方程(3)的求解。

本步骤的主要原理为：

在步骤5确定一维土壤非饱和流模拟模型时，由于垂向上的土壤入渗作用远远大于水平方向的流动，因此，步骤5中的一维土壤非饱和流运动方程为一维方程，即：只是土壤深度方向的一维方程，用于表征垂向上的土壤入渗作用的运动过程。

而地下水在水平方向上的运动速度远远大于垂向上的流速，因此地下水流模型需要建立三维模型。任意一维土壤非饱和流在有限元节点P的达西流速与其影响面积A_P相乘，即可得到对地下水的补给量，用于三维地下水流模型的源汇项W_i，由此实现了土壤非饱和水和地下水的耦合。

步骤7，基于以下的三维地下水流运动方程，建立三维地下水流模型：

其中：

h为地下水水头，为因变量；

k_x、k_y、k_z分别为x，y，z方向的水力传导系数；

S₀为含水层储水系数；

t为时间；

W_i为源汇项，等于步骤7计算得到的地下水补给量Q；

运行所述三维地下水流模型，实现对地下水流运动状态的模拟。

传统地下水流模拟软件，例如MODFLOW不计算非饱和带水分运动，而是通过蒸发量及入渗量的线性插值表达非饱和带对饱和带的影响。本发明方法通过将非饱和带底部通量，即地下水补给量Q作为饱和带上边界条件分别求解方程(1)和(3)并迭代循环，实现了非饱和带垂向水分分布特征的详细刻画，并能准确描述非饱和带与饱和带的水力联系。

本发明提供一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，具有以下优点：

(1)将研究区域的遥感栅格数据和有限元节点实现了空间上的关联，并通过本发明提供的映射投影方式，能够快速搜索定位到与有限元节点在空间上对应的栅格，从而将栅格遥感数据传递到有限元节点，作为有限元节点的输入驱动数据，保证了后续采用有限元方法进行地下水流模拟的精度；

(2)一维土壤非饱和流模拟模型和三维地下水流模型进行了耦合，一维土壤非饱和流模拟模型得到的达西流速与影响面积的乘积，得到有限元节点P对应的地下水补给量Q，直接作为三维地下水流模型的上边界条件，保证了对地下水流的地下水水头运动状态的模拟精度。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (2)

1.一种基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，其特征在于，包括以下步骤：

步骤1，获取研究区域的遥感栅格数据；其中，所述研究区域的遥感栅格数据对应M行N列个栅格；每个栅格的边长为s；基于所述遥感栅格数据，得到每个栅格顶点的降雨入渗量值；

步骤2，对研究区域的土壤和地下水层进行有限元网格剖分，水平方向上得到多个三角形形状的有限元网格，进而也得到多个有限元节点，垂直方向上按不同地层介质的厚度分层；

步骤3，将所述研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中；对于每个有限元节点P，均采用以下方式处理，由此得到与每个有限元节点P在空间上对应的栅格，其中，有限元节点P在空间上对应的栅格含义为：当将研究区域的遥感栅格数据和所述有限元网格投影到同一个平面x-y坐标系中时，有限元节点P落入到某个栅格中，即为与有限元节点P在空间上对应的栅格；

步骤3.1，在x-y坐标系中，对于M行N列的栅格，其左下角的坐标为(x₀，y₀)；

假设某个有限元节点P的坐标为P(x_P，y_P),则其在空间上对应的栅格W的行数h₁和列数h₂通过以下公式得出：

步骤3.2，在定位到与有限元节点P在空间上对应的栅格W后，栅格W的左下角、右下角、右上角和左上角这四个顶点分别记为：W₁、W₂、W₃和W₄；通过下式分别计算得到W₁、W₂、W₃和W₄的坐标：

W₁的坐标为：(x₀+s(h₁-1)，y₀+s(h₂-1))；

W₂的坐标为：(x₀+sh₁)，y₀+s(h₂-1))；

W₃的坐标为：(x₀+sh₁，y₀+sh₂))；

W₄的坐标为：(x₀+s(h₁-1)，y₀+sh₂))；

步骤4，采用下述方法计算得到有限元节点P的降雨入渗量q^(p)；

步骤4.1，在得到有限元节点P在空间上对应的栅格W的四个顶点W₁、W₂、W₃和W₄的坐标后，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的降雨入渗量分别为q₁,q₂,q₃,q₄；

步骤4.2，查找所述遥感栅格数据，获得顶点W₁、W₂、W₃和W₄的基函数分别为：N₁(x_P，y_P)、N₂(x_P，y_P)、N₃(x_P，y_P)和N₄(x_P，y_P)；即：栅格顶点的基函数与有限元节点P的坐标相关；

步骤4.3，采用下述的线性插值函数计算得到有限元节点P的降雨入渗量q^(p)的值：

q^(p)＝q₁·N₁(x_P，y_P)+q₂·N₂(x_P，y_P)+q₃·N₃(x_P，y_P)+q₄·N₄(x_P，y_P)；

步骤5，基于以下的一维土壤非饱和流运动方程，建立一维土壤非饱和流模拟模型：

其中：

n_r为土壤孔隙度，S为土壤水的饱和度；

t为时间；

z为坐标，代表土壤深度；

q_s为源项，等于步骤4得到的每个有限元节点P的降雨入渗量；

q_z为达西流速，为待求值；

运行所述一维土壤非饱和流模拟模型，计算得到与每个有限元节点P对应的达西流速值；

步骤6，假设有限元节点P为u个有限元网格的共享节点，u个有限元网格分别为V₁、V₂,…,Vu；则采用下式计算有限元节点P的影响面积A_P:

其中：N^K为有限元节点P在有限元网格V^K上的基函数，K＝1、2,…,u；

Ω_K为有限元节点P在有限元网格V^K上的积分空间；

采用下式计算得到有限元节点P对应的地下水补给量Q：

Q＝q_z*A_P

步骤7，基于以下的三维地下水流运动方程，建立三维地下水流模型：

其中：

h为地下水水头，为因变量；

k_x、k_y、k_z分别为x，y，z方向的水力传导系数；

S₀为含水层储水系数；

t为时间；

W_i为源汇项，等于步骤7计算得到的地下水补给量Q；

运行所述三维地下水流模型，实现对地下水流运动状态的模拟。

2.根据权利要求1所述的基于土壤非饱和水和地下水耦合的区域水流运动模拟方法，其特征在于，步骤1中，采用以下方法得到每个栅格的顶点的降雨入渗量值：

通过所述遥感栅格数据，得到与栅格顶点位置对应的降雨量P_a、太阳辐射能R_a、大气的平均气温T_a，最高气温T_max和最低气温T_min，则采用下式得到栅格的顶点的降雨入渗量值q：

q＝P_a-E；

其中：

E为蒸发潜能；是大气温度和太阳辐射能的函数；

κ，λ，β均为经验公式常数，分别取值2.3*10^-3，2.45和17.8。

Images