CN106503463A - 模拟海平面上升情况下海水入侵内陆边界的处理方法 - Google Patents

Abstract

本发明所属技术领域为“海水入侵内陆边界控制方法”领域，为建立地下水气‑液二相流及溶质运移模型，来模拟分析不同内陆边界条件对沿海地区海平面上升引起的海水入侵的效果。本发明采用的技术方案是，模拟海平面上升情况下海水入侵内陆边界的处理方法，步骤如下：步骤一：建立地下水气‑液二相流及溶质运移模型，包括基本控制方程及辅助方程，假定步骤二：模型求解；步骤三：模型边界条件确定：模型计算中的边界条件包括狄利克雷Dirichlet边界条件和诺伊曼Neumann边界条件两种，其数学处理方法如下：本发明主要应用于防止海水入侵场合。

Description

模拟海平面上升情况下海水入侵内陆边界的处理方法

技术领域

本发明所属技术领域为“海水入侵内陆边界控制方法”领域。具体讲,涉及模拟海平面上升情况下海水入侵内陆边界的处理方法。

背景技术

世界上大部分的人口分布在沿海地区，近年来,受地壳运动、自然压密、气候变暖引起的海平面上升、大量抽用地下水等自然及人为因素的影响，致使沿海各地区已不同程度地遭到海水入侵,严重制约沿海地区的社会经济发展，引起了世界各沿海国家的充分重视。最近，由于海平面上升引起的海水入侵问题的研究也受到越来越多的关注，其中，内陆边界条件对模拟海平面上升引起的海水入侵有至关重要的影响。因此，开展内陆边界条件对模拟海平面上升引起的海水入侵影响的深入研究，寻求一种合适的内陆边界条件的处理方法，具有重要的科学价值和社会意义。

虽然流量控制系统边界和水头控制系统边界在数值模拟中是最常应用的两种边界条件，但事实上，以上两种边界都有一定的局限性，因为随着海平面的上升，真实的模拟区域内陆边界既不是水头恒定，也不是流量恒定，而是一个间接的过程，因此有人提出了第三种边界条件——依靠水头的流量边界条件，这种边界条件在一些特地场地也进行了海平面上升对海水入侵的影响评估，而且也有人提出了解析解，但是却没有基于假定的理想的概念模型的一般性的研究，更缺乏基于理想的数值模拟在预测海平面上升对海水入侵的三种不同内陆边界条件的对比并得出合适的模拟海平面上升情况下海水入侵内陆边界的处理方法。

发明内容

为克服现有技术的不足，本发明旨在在深入分析地下多相流运移机理的基础上，建立地下水气-液二相流及溶质运移模型，来模拟分析不同内陆边界条件对沿海地区海平面上升引起的海水入侵的效果。本发明采用的技术方案是，模拟海平面上升情况下海水入侵内陆边界的处理方法，步骤如下：

步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，假定系统一直处于恒温状态。具体如下：

模型的基本控制方程为：

式中，M^κ表示κ组分的累积质量密度，组分是指纯水w、参考盐水b和空气a，F^κ为κ组分的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项；

(1)累积质量密度等于β相中组分κ的质量之和，其表达式为：

式中，β表示流体相，流体相为气相或液相，液相为纯水和参考盐水的混合物，表示孔隙率，S_β为β相的饱和度，ρ_β为β相的密度，p_atm为κ组分占β相的质量百分数；

(2)平流流量p_g的表达式为：

式中，F_β为相的平流流量，符合达西定律，其多相流形式的表达式为：

其中，k为固有渗透系数，k_rβ为β相的相对渗透系数，μ_β为β相的粘滞性系数，为梯度符号，p_β为β相的孔隙压力，g为重力加速度矢量；

(3)扩散流量的表达式为：

式中，为κ组分在β相中的分子扩散系数；τ₀τ_β(S_β)为迂曲度，τ₀是多孔介质特性参数，τ_β＝τ_β(S_β)是饱和度的函数；

分子扩散系数随压力p和温度T变化的计算式为：

式中，θ为温度系数；

非饱和地下水流及热流传输TOUGH的第二个版本TOUGH2中提供了三种迂曲度模型：

τ₀τ_β＝τ₀k_rβ(S_β) 相对渗透率模型

米林顿夸克模型

τ₀τ_β＝τ₀S_β 常量扩散 (7)

步骤二：模型求解：以TOUGH2的水、盐水和空气的混合模块EOS7为工具，空间上采用积分形式的有限差分方法IFDM进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉弗森Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；

步骤三：模型边界条件确定：模型计算中的边界条件包括狄利克雷Dirichlet边界条件和诺伊曼Neumann边界条件两种，其数学处理方法如下：

(1)Dirichlet边界条件

Dirichlet边界条件上，边界条件单元的主要变量在计算过程中保持不变：

①对于空气边界，其边界条件单元的相态为仅有气相状态,变量为p_g,X_b,T，对于与大气接触的边界上p_g＝p_atm，X_b＝0.0；T为气温；

②对于已知水头边界，包括地下淡水水头边界和海水水头边界，其边界条件单元的相态均为仅有液相状态,变量为p_g，X_b，T，由于液相饱和状态下毛细压力为零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l为淡水或海水的密度，h为边界上的水深；地下淡水水头边界上X_b等于零，海水边界上的X_b是根据参考盐水的特性以及海水的盐度计算出来；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，边界条件单元的单位流量对应式中的源汇项，流入为正：

①对于不透水边界，是一类特殊的Neumann边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设置与之相邻的边界条件单元，则与不透水面单元的流量交换为零；

②对于降雨入渗条件的处理方法，则设置一个适当大小的源汇项来实现；

(3)依靠水头的流量边界条件

设区域外的已知水位又称参考水头为H_ref，水流流入或者流出研究区域的流量大小与研究区内的水头H和参考水头H_ref之间的水头差成比例，并且用如下线性方程式表示：

q＝-C(H-H_ref) (15)

其中q是从内陆侧流入研究区的单宽淡水补给量，C是导水系数。

而研究区域淡水侧的等效淡水水头H_f跟淡水补给量q_f在承压含水层和非承压含水层中分别有如下关系：

承压含水层：

非承压含水层：

其中K是含水层的渗透系数，α＝ρ_f/(ρ_s-ρf)是相对密度，ρ_f和ρ_s分别是淡水密度和海水密度，L和B分别为研究区域含水层的长度和厚度，H_s是距离含水层底部的平均海水深度。

把公式(15)分别带入公式(16)和(17)中，并且让H＝H_f和q＝q_f，可得出如下关系式：

承压含水层：

非承压含水层：

模型求解具体步骤如下：

(1)空间上采用积分有限差分法IFDM进行离散

首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散，对于任意单元n，单元体积为V_n，边界面为Γ_n，单元的质量平衡方程的积分形式如下：

式中，n为表面单元dΓ_n的单位法向量，指向控制单元体内为正；

引入体积平均值，有

式中，为M^κ，q^κ在V_n上的平均值。

Γ_n上的面积分可近似为其所包含的各个表面A_nm的面积分的平均值之和，有

式中，m为与单元n相邻的所有单元，A_nm是单元n和m相邻的交界面，是F^κ在面A_nm上的内法线方向的平均值；

将式(9)、(10)和(11)代入到式(8)中，得到一组关于时间的一阶微分方程组

(2)时间上采用一阶向后差分方法进行离散

对式(12)的时间微分采用一阶向后差分方法，得到任意单元的全隐式非线性方程组，见式(13)：

式中，引入了组分κ＝w,g,a的余量△t为时间步长，上标k和k+1表示两相邻的时间步长指标；右端的流量项和源汇项均采用新的时间步长值；

(3)Newton-Raphson迭代方法

对式运用Newton-Raphson迭代方法进行线性化，引入迭代指标p，对式(13)中的余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，得到包含3×NEL个方程的线性方程组，NEL为计算域内单元数，并且以两迭代步的增量为未知量；最后得到大型稀疏系数矩阵的线性方程组，如式(14)：

还包括模型验证步骤：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证；

还包括如下分析步骤：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件，采用依靠水头的流量边界，分析模拟海平面上升情况下海水入侵的效果：

(1)入侵程度变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升过程中，用相对盐度等值线表征的海水入侵范围始终介于流量恒定和水头恒定中间，具体如下：当内陆边界为流量恒定时，随海平面的上升，海水入侵程度基本与初始情况保持不变；而当内陆边界为水头恒定时，随着海平面的上升，海水入侵范围明显增大；当内陆边界为依靠水头的流量边界时，随着海平面的上升，海水入侵范围有所增大，但入侵程度明显比水头恒定情况小；

(2)研究区域左右两侧等效水头差和左侧补给量的变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升，研究区域左右两侧等效水头差和左侧补给量都在逐渐减小，但是内陆边界为流量恒定时，以上两项都基本保持不变，而水头恒定时以上两项减少的幅度要大于依靠水头的流量边界。

本发明的特点及有益效果是：

(1)建立了地下水气-液二相流及溶质运移模型，可以用来模拟不同内陆边界条件对，由海平面上升引起的海水入侵的影响。

(2)本发明采用了基于假定的理想的概念模型来模拟海平面上升下海水入侵的内陆边界处理方法，这种依靠水头的流量边界，更加符合实际，因为它的补给会随着海平面的上升而逐渐减小，它是由参考水头(远处参考的某一水面)与海平面之间的高度差以及含水层的性质共同决定的。内陆边界条件对海水入侵有至关重要的影响，因此提出一种更加合理的模拟海平面上升下海水入侵的内陆边界处理方法，对于合理预测沿海地区海水入侵的位置和移动情况，有效防治地下水污染，保护有限的地下淡水资源具有重大的现实意义，而且对于水资源管理者更好的评估和管理海平面上升对海水入侵潜在的长期影响具有一定的指导意义。

附图说明：

图1含水层的一般概念模型：(a)承压含水层，(b)非承压含水层。

图2 Henry’s问题的研究范围与边界条件。

图3稳态下相对盐度为0.5的等值线位置及其与其他文献结果的对比图。

图4流量恒定(FC)和水头恒定(HC)下，达到稳定时相对盐度为0.5的等值线分布。

图5在承压含水层中三种不同内陆边界在不同海平面上升情况下，达到稳态时相对盐度为0.5的等值线分布图：(a)流量恒定系统边界(FC systems),(b)依靠水头的流量边界(GHB,和(c)水头恒定系统边界(HC systems)。

图6承压和非承压含水层中海平面上升前稳态下相对盐度为0.5的等值线分布。

图7在非承压含水层中三种不同内陆边界在不同海平面上升情况下，达到稳态时相对盐度为0.5的等值线分布图：(a)流量恒定系统边界(FC systems),(b)依靠水头的流量边界(GHB,和(c)水头恒定系统边界(HC systems)。

图8本发明的模拟方法流程图。

具体实施方式

本发明步骤如下：

步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，假定系统一直处于恒温状态。具体如下：

模型的基本控制方程为：

式中，M^κ表示κ组分(纯水w、参考盐水b和空气a)的累积质量密度，F^κ为κ组分的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项。

(1)累积质量密度等于β相中组分κ的质量之和，其表达式为：

式中，β表示流体相(气相或液相)，液相为纯水和参考盐水的混合物，表示孔隙率，S_β为β相的饱和度，ρ_β为β相的密度，p_atm为κ组分占β相的质量百分数。

(2)平流流量p_g的表达式为：

式中，F_β为相的平流流量，符合达西定律，其多相流形式的表达式为：

其中，k为固有渗透系数，k_rβ为β相的相对渗透系数，μ_β为β相的粘滞性系数，p_β为β相的孔隙压力，g为重力加速度矢量。

(3)扩散流量的表达式为：

式中，为κ组分在β相中的分子扩散系数；τ₀τ_β(S_β)为迂曲度，τ₀是多孔介质特性参数，τ_β＝τ_β(S_β)是饱和度的函数。

分子扩散系数随压力p和温度T变化的计算式为：

式中，θ为温度系数，通常可取值为1.8。

TOUGH2中提供了三种迂曲度模型：分子扩散系数τ₀一般取1.0

τ₀τ_β＝τ₀k_rβ(S_β) (相对渗透率模型)

(米林顿夸克模型)

τ₀τ_β＝τ₀S_β (常量扩散) (7)

步骤二：模型求解：以非饱和地下水流及热流传输TOUGH的第二个版本TOUGH2的水、盐水和空气的混合模块EOS7为工具，空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉弗森Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：

(1)空间上采用积分有限差分法(IFDM)进行离散

首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散。对于任意单元n，单元体积为V_n，边界面为Γ_n，单元的质量平衡方程的积分形式如下：

式中，n为表面单元dΓ_n的单位法向量，指向控制单元体内为正。

引入适当的体积平均值，有

式中，为M^κ，q^κ在V_n上的平均值。

Γ_n上的面积分可近似为其所包含的各个表面A_nm的面积分的平均值之和，有

式中，m为与单元n相邻的所有单元，A_nm是单元n和m相邻的交界面，是F^κ在面A_nm上的内法线方向的平均值。

将式(9)、(10)和(11)代入到式(8)中，得到一组关于时间的一阶微分方程组

(2)时间上采用一阶向后差分方法进行离散

对式(12)的时间微分采用一阶向后差分方法，得到任意单元的全隐式非线性方程组，见式(13)：

式中，引入了组分κ＝w,g,a的余量△t为时间步长，上标k和k+1表示两相邻的时间步长指标；右端的流量项和源汇项均采用新的时间步长值，这种全隐式方法提高了模型求解的数值稳定性。

(3)Newton-Raphson迭代方法

对式运用Newton-Raphson迭代方法进行线性化。引入迭代指标p，对式(13)中的余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，得到包含3×NEL(计算域内单元数)个方程的线性方程组，并且以两迭代步的增量为未知量；最后得到大型稀疏系数矩阵的线性方程组，如式(14)：

步骤三：模型边界条件确定：模型计算中的边界条件包括狄利克雷Dirichlet边界条件和诺伊曼Neumann边界条件两种，其数学处理方法如下：

(1)Dirichlet边界条件

Dirichlet边界条件上，边界条件单元的主要变量在计算过程中保持不变，为此，设定边界条件单元的体积非常大，如1×10⁴⁰m³。当边界单元的体积相对于土体单元很大时，与土体单元的流量交换将不会改变边界单元的主要变量值。

①对于空气边界，其边界条件单元的相态为仅有气相状态,主要变量为p_g,X_b,T，对于与大气接触的边界上p_g＝p_atm，X_b＝0.0；T为气温。

②对于已知水头边界，包括地下淡水水头边界和海水水头边界，其边界条件单元的相态均为仅有液相状态,主要变量为p_g、X_b、T，由于液相饱和状态下毛细压力为零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l为淡水或海水的密度，h为边界上的水深；地下淡水水头边界上X_b等于零，海水边界上的X_b可根据参考盐水的特性以及海水的盐度计算出来。

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，边界条件单元的单位流量对应式中的源汇项，流入为正，可以是常量，也可以随时间变化。

①对于不透水边界，是一类特殊的Neumann边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设置与之相邻的边界条件单元，则与不透水面单元的流量交换为零。

②对于降雨入渗条件的处理方法，则可设置一个适当大小的源汇项来实现。

(3)依靠水头的流量边界条件

依靠水头的流量边界条件受研究区域外的已知水位的较大地表水体(如湿地、河流、水库等)影响，区域外的已知水位又称参考水头H_ref。水流流入或者流出研究区域的流量大小与研究区内的水头H和参考水头H_ref之间的水头差成比例，并且可以用如下线性方程式表示：

q＝-C(H-H_ref) (15)

其中q是从内陆侧流入研究区的单宽淡水补给量，C是导水系数。

而研究区域淡水侧的等效淡水水头H_f跟淡水补给量q_f在承压含水层和非承压含水层中分别有如下关系：

承压含水层：

非承压含水层：

其中K是含水层的渗透系数，α＝ρ_f/(ρ_s-ρf)是相对密度，ρ_f和ρ_s分别是淡水密度和海水密度，L和B分别为研究区域含水层的长度和厚度，H_s是距离含水层底部的平均海水深度。

把公式(15)分别带入公式(16)和(17)中，并且让H＝H_f和q＝q_f，可得出如下关系式：

承压含水层：

非承压含水层：

步骤四：模型验证：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证；

步骤五：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件，采用依靠水头的流量边界，分析模拟海平面上升情况下海水入侵的效果：

(1)入侵程度变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升过程中，用相对盐度等值线表征的海水入侵范围始终介于流量恒定和水头恒定中间。具体如下：当内陆边界为流量恒定时，随海平面的上升，海水入侵程度基本与初始情况保持不变；而当内陆边界为水头恒定时，随着海平面的上升，海水入侵范围明显增大；当内陆边界为依靠水头的流量边界时，随着海平面的上升，海水入侵范围有所增大，但入侵程度明显比水头恒定情况小。

(2)研究区域左右两侧等效水头差和左侧补给量的变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升，研究区域左右两侧等效水头差和左侧补给量都在逐渐减小，但是内陆边界为流量恒定时，以上两项都基本保持不变，而水头恒定时以上两项减少的幅度要大于依靠水头的流量边界。

下面结合附图，针对模拟海平面上升下海水入侵的内陆边界处理方法进行具体说明，其步骤如下：

(1)建立模拟海水入侵的数学模型：该数学模型描述的是承压含水层未受污染的淡水受海水入侵的情况，即经典的Henry’s问题。模型选取200米长、100米深、1米厚的含水层为研究对象，其上、下两层不透水顶板、底板之间的含水层均质、各向同性，内陆侧边界上有恒定的淡水流入量(或者等效的淡水作用下的静水压力)，海水侧边界上分布着密度较大的海水作用下的静水压力，是一个理想化的准二维模型。其研究范围及边界条件见图2，模型相关参数取值见表1。假定含水层系统处于恒温状态(T＝25℃)。

模型的基本控制方程为：

式中，M^κ表示κ组分(纯水w、参考盐水b和空气a)的累积质量密度，F^κ为组分κ的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项。

①累积质量密度等于β相中组分κ的质量之和，其表达式为：

式中，β表示流体相(气相或液相)，液相为纯水和参考盐水的混合物，表示孔隙率，S_β为β相的饱和度，ρ_β为β相的密度，p_atm为κ组分占β相的质量百分数。

②平流流量p_g的表达式为：

式中，F_β为相的平流流量，符合达西定律，其多相流形式的表达式为：

其中，k为固有渗透系数，k_rβ为β相的相对渗透系数，μ_β为β相的粘滞性系数，p_β为β相的孔隙压力，g为重力加速度矢量。

③扩散流量的表达式为：

式中，为κ组分在β相中的分子扩散系数；τ₀τ_β(S_β)为迂曲度，τ₀是多孔介质特性参数，τ_β＝τ_β(S_β)是饱和度的函数。

分子扩散系数随压力p和温度T变化的计算式为：

式中，θ为温度系数，通常可取值为1.8。

TOUGH2中提供了三种迂曲度模型：分子扩散系数τ₀一般取1.0

τ₀τ_β＝τ₀k_rβ(S_β) (相对渗透率模型)

(米林顿夸克模型)

τ₀τ_β＝τ₀S_β (常量扩散) (7)

(2)模型求解：对计算域进行网格剖分，含水层被划分为5m×5m×1m的8节点的均匀六面体单元，共800个单元，1722个节点。以TOUGH2/EOS7为工具，空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，见式(13)；模型的线性化采用Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组，见式(14)：

式中，引入了组分κ＝w,g,a的余量 表示M^κ在控制体积V_n上的平均值；表示q^κ在控制体积V_n上的平均值；m为与单元n相邻的所有单元；A_nm是单元n和m相邻的交界面的面积；是F^κ在面A_nm上的内法线方向的平均值；△t为时间步长；上标k和k+1表示两相邻的时间步指标。

引入迭代指标p，对式(13)中余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，得到包含3×NEL(计算域内单元数)个方程的线性方程组，并且以两迭代步的增量为未知量：

(3)模型边界条件确定：假定系统处于恒温状态T＝25℃在整个模拟过程中保持不变。在初始条件下，计算域均为液相饱和状态，毛细压力可以忽略，相对渗透系数取值1.0，主要变量为p_g、X_b、T。因此有p_g＝p_atm+ρ_lgh，其中ρ_l为淡水或海水的密度，h为边界上的水深。对于淡水侧(左侧)，作用淡水源汇项，将Q_in＝V_indρ_w(d为含水层深度100m，ρ_w为淡水密度)等于7.639×10^-3kg/s，均匀分布在左侧边界上，主要变量p_g＝p_atm+ρ_wg(100-z)，X_b＝0.0，和T＝25℃；对于海水侧(右侧)，受海水作用，主要变量p_g＝p_atm+ρ_sg(100-z)，X_b＝1.0，和T＝25℃；其中，p_atm为大气压力，等于1.013×10⁵Pa。初次计算的初始条件为：所有单元液相饱和，压力均为大气压力与静水(淡水)压力之和即p_g＝p_atm+ρ_wg(100-z)，X_b＝0.0，

以上实际上属于流量恒定的情况，为了使初始稳态条件相同，当水头恒定时，由公式(16)计算可知，在承压含水层中恒定的等效水头H_f＝102.57米(q_f＝0.66m²/d且H_s＝100m)。因此其内陆边界条件设置为P_l＝P_atm+ρ_wg(102.57-z),X_b＝1.0 and然而依靠水头的流量的边界条件的设置中，假定参考水头H_ref＝105m，因此由公式(15)知，C＝0.2716m/s(q_f＝0.66m²/d且H_f＝102.57m)，且其初始条件的设置与流量恒定完全相同。

(4)模型验证：

在静力平衡情况下，相对(海水)盐度X_b＝0.5的等值线位置及其与已有文献计算结果的对比如图3。由图可见，本模型的计算结果与其他文献的结果基本一致，从而验证了计算模型模拟海水入侵问题的有效性。

(5)根据地下水气-液二相流及溶质运移模型，分析不同内陆边界条件对模拟海平面上升的海水入侵的效果：

以咸-淡水静力平衡情况作为初始条件，由图4可知三种不同内陆边界条件的海水入侵情况是等效的(图4没有标注依靠水头的流量(GHB)边界情况，因为初始情况下，GHB边界与流量边界完全相同)。

为了比较不同内陆边界条件对海平面上升的海水入侵的影响，右侧海水边界的海水位分别设置为100.25,100.50,100.75 and 101.0m(对应于海平面分别上升0.25,0.50,0.75 and 1.0m)。对于流量恒定和水头恒定的内陆边界条件，左侧(内陆边界)补给流量q_f＝0.66m²/d和等效淡水水头H_f＝102.57m分别保持不变。而对于水头依靠的流量内陆边界条件，由公式(19)可知，补给流量q_f随着海水位H_s的增大而减小，具体的q_f分别为0.61,0.57,0.52,和0.48m²/d对应于H_s分别为100.25,100.50,100.75和101.0m。对于水头依靠的流量内陆边界条件的设置跟流量恒定条件的设置类似，除了随海水位变化的q_f。

图5比较了在承压含水层中三种不同内陆边界在不同海平面上升情况下，达到稳态时相对盐度为0.5的等值线分布。结果表明，内陆边界条件对模拟海平面上升引起的海水入侵有至关重要的影响。水头恒定(HC)的内陆边界条件对海平面上升引起的海水入侵最为敏感、脆弱，因为随着海平面的上升，左右两侧的水力梯度逐渐减小；而流量恒定(FC)的内陆边界条件对模拟海平面上升引起的海水入侵却几乎不变，因为海平面上升时左侧淡水水头同等上升，水力梯度几乎保持不变；而依靠水头的流量边界(GHB)条件介于以上两种边界之间，因为它的淡水补给会随着海平面的上升而逐渐减小，它是由参考水头(远处参考的某一水面)与海平面之间的水头差以及含水层的性质共同决定的。

为了研究非承压含水层的模拟内陆边界对海平面上升的海水入侵的情况，首先把Henry问题扩展到非承压区域，保持相同的边界条件和参数，但是上层边界为潜水自由水面且增加了5米，网格离散化范围为0.25到1米不等。因此非承压含水层研究区域的长为200m，厚度为105米。为了更好的研究海平面上升0到1米的对海水入侵的情况，网格在100米上下被细化。

非承压含水层边界条件的设置类似于承压含水层，对于流量恒定的内陆边界条件，左侧单宽补给q_f＝0.66m²/d保持不变。对于水头恒定的内陆边界条件，等效淡水水头由公式(17)可知，H_f＝102.54m(q_f＝0.66m²/d且H_s＝100m)，而对于依靠水头的流量的边界条件，由公式(15)知，C＝0.2681m/s(q_f＝0.66m²/d，H_ref＝105且H_f＝102.57m)。当右侧海水边界的海水位分别设置为100.25,100.50,100.75 and 101.0m(对应于海平面分别上升0.25,0.50,0.75 and 1.0m)时，对于流量恒定和水头恒定的内陆边界条件，左侧(内陆边界)补给流量q_f＝0.66m²/d和等效淡水水头H_f＝102.54m分别保持不变。而对于水头依靠的流量内陆边界条件，由公式(21)可知，补给流量q_f随着海水位H_s的增大而减小，具体的q_f分别为0.62,0.57,0.53,和0.48m²/d对应于H_s分别为100.25,100.50,100.75和101.0m。

由图6承压和非承压含水层中海平面上升前稳态下相对盐度为0.5的等值线分布可知，这四条曲线基本重合，表明海水入侵过程不论在承压还是非承压含水层中，其初始稳态条件下相一致。对于海平面上升情况下，不同内陆边界对海水入侵的影与承压含水层情况基本一致(图7)。

综上所述，对本实施例总结如下：

(1)基于标准的理想的概念模型，采用数值模拟的方法研究了不同内陆边界条件对海平面上升引起的海水入侵效果。首先利用地下水气-液二相流及溶质运移模型模拟分析了经典的Henry问题的咸-淡水静力平衡情况，通过与以前文献结果的对比，验证了所采用模型的有效性。然后基于咸-淡水静力平衡情况，模拟分析了不同内陆边界条件对海平面上升引起的海水入侵效果。

(2)对模拟结果的分析表明，内陆边界条件对模拟海平面上升引起的海水入侵有至关重要的影响。水头恒定(HC)的内陆边界条件对海平面上升引起的海水入侵最为敏感、脆弱，因为随着海平面的上升，左右两侧的水力梯度逐渐减小；而流量恒定(FC)的内陆边界条件对模拟海平面上升引起的海水入侵却几乎不变，因为海平面上升时左侧淡水水头同等上升，水力梯度几乎保持不变；而依靠水头的流量边界(GHB)条件介于以上两种边界之间且更符合实际。因为它的淡水补给会随着海平面的上升而逐渐减小，它是由参考水头(远处参考的某一水面)与海平面之间的水头差以及含水层的性质共同决定的。

(3)依靠水头的流量内陆边界条件是一种更加合理的模拟海平面上升下海水入侵的内陆边界处理方法，这对于合理预测沿海地区海水入侵的位置和移动情况，有效防治地下水污染，保护有限的地下淡水资源具有重大的现实意义，而且对于水资源管理者更好的评估和管理海平面上升对海水入侵潜在的长期影响具有一定的指导意义。

Claims (4)

1.一种模拟海平面上升情况下海水入侵内陆边界的处理方法，其特征是，步骤如下：

步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，假定系统一直处于恒温状态。具体如下：

模型的基本控制方程为：

$\frac{\∂M^{\mathrm{\κ}}}{\∂t}=-div\left(F^{\mathrm{\κ}}\right)+q^{\mathrm{\κ}}---\left(1\right)$

式中，M^κ表示κ组分的累积质量密度，组分是指纯水w、参考盐水b和空气a，F^κ为κ组分的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项；

(1)累积质量密度M^κ等于β相中组分κ的质量之和，其表达式为：

$M^{\mathrm{\κ}}=\underset{\mathrm{\β}}{\Σ}{X}_{\mathrm{\β}}^{\mathrm{\κ}}{\mathrm{\φS}}_{\mathrm{\β}}{\mathrm{\ρ}}_{\mathrm{\β}}---\left(2\right)$

式中，β表示流体相，流体相为气相或液相，液相为纯水和参考盐水的混合物，φ表示孔隙率，S_β为β相的饱和度，ρ_β为β相的密度，为κ组分占β相的质量百分数；

(2)平流流量的表达式为：

$F_{adv}^{\mathrm{\κ}}=\underset{\mathrm{\β}}{\Σ}{X}_{\mathrm{\β}}^{\mathrm{\κ}}{F}_{\mathrm{\β}}---\left(3\right)$

式中，F_β为β相的平流流量，符合达西定律，其多相流形式的表达式为：

$F_{\mathrm{\β}}=-k\frac{{\mathrm{\ρ}}_{\mathrm{\β}}{k}_{\mathrm{\γ}\mathrm{\β}}\left(S_w\right)}{{\mathrm{\μ}}_{\mathrm{\β}}}(\▿p_{\mathrm{\β}}-{\mathrm{\ρ}}_{\mathrm{\β}}g)---\left(4\right)$

其中，k为固有渗透系数，k_rβ为β相的相对渗透系数，μ_β为β相的粘滞性系数，为梯度符号，p_β为β相的孔隙压力，g为重力加速度矢量；

(3)扩散流量的表达式为：

式中，为κ组分在β相中的分子扩散系数；τ₀τ_β(S_β)为迂曲度，τ₀是多孔介质特性参数，τ_β＝τ_β(S_β)是饱和度的函数；

分子扩散系数随压力p和温度T变化的计算式为：

式中，θ为温度系数；

非饱和地下水流及热流传输TOUGH的第二个版本TOUGH2中提供了三种迂曲度模型：

τ₀τ_β＝τ₀k_rβ(S_β)相对渗透率模型

米林顿夸克模型

τ₀τ_β＝τ₀S_β常量扩散 (7)

步骤二：模型求解：以TOUGH2的水、盐水和空气的混合模块EOS7为工具，空间上采用积分形式的有限差分方法IFDM进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉弗森Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；

步骤三：模型边界条件确定：模型计算中的边界条件包括狄利克雷Dirichlet边界条件和诺伊曼Neumann边界条件两种，其数学处理方法如下：

(1)Dirichlet边界条件

Dirichlet边界条件上，边界条件单元的主要变量在计算过程中保持不变：

①对于空气边界，其边界条件单元的相态为仅有气相状态,变量为p_g,X_b,T，对于与大气接触的边界上p_g＝p_atm，X_b＝0.0；T为气温；

②对于已知水头边界，包括地下淡水水头边界和海水水头边界，其边界条件单元的相态均为仅有液相状态,变量为p_g、X_b、T，由于液相饱和状态下毛细压力为零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l为淡水或海水的密度，h为边界上的水深；地下淡水水头边界上X_b等于零，海水边界上的X_b是根据参考盐水的特性以及海水的盐度计算出来；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，边界条件单元的单位流量对应式中的源汇项，流入为正：

①对于不透水边界，是一类特殊的Neumann边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设置与之相邻的边界条件单元，则与不透水面单元的流量交换为零；

②对于降雨入渗条件的处理方法，则设置一个适当大小的源汇项来实现；

(3)依靠水头的流量边界条件

设区域外的已知水位又称参考水头为H_ref，水流流入或者流出研究区域的流量大小与研究区内的水头H和参考水头H_ref之间的水头差成比例，并且用如下线性方程式表示：

q＝-C(H-H_ref) (15)

其中q是从内陆侧流入研究区的单宽淡水补给量，C是导水系数。

而研究区域淡水侧的等效淡水水头H_f跟淡水补给量q_f在承压含水层和非承压含水层中分别有如下关系：

承压含水层：

$q_f=\frac{KB}{L}(H_f-\frac{(1+\mathrm{\α})}{\mathrm{\α}}{H}_s+\frac{B}{2\mathrm{\α}})---\left(16\right)$

非承压含水层：

$q_f=\frac{K}{2L}({H_f}^2-\frac{(1+\mathrm{\α})}{\mathrm{\α}}{H_s}^2)---\left(17\right)$

其中K是含水层的渗透系数，α＝ρ_f/(ρ_s-ρ_f)是相对密度，ρ_f和ρ_s分别是淡水密度和海水密度，L和B分别为研究区域含水层的长度和厚度，H_s是距离含水层底部的平均海水深度。

把公式(15)分别带入公式(16)和(17)中，并且让H＝Hf和q＝q_f，可得出如下关系式：

承压含水层：

$H_f=\frac{\frac{CL}{KB}{H}_{ref}+\frac{(1+\mathrm{\α})}{\mathrm{\α}}{H}_s-\frac{B}{2\mathrm{\α}}}{\frac{CL}{KB}+1}---\left(18\right)$

$q_f=\frac{{\mathrm{CH}}_{ref}-\frac{(1+\mathrm{\α})}{\mathrm{\α}}{\mathrm{CH}}_s+\frac{BC}{2\mathrm{\α}}}{\frac{CL}{KB}+1}---\left(19\right)$

非承压含水层：

$H_f=\frac{-CL+\sqrt{C^2{L}^2+\frac{(1+\mathrm{\α})}{\mathrm{\α}}{K}^2{H_s}^2+2{\mathrm{CKLH}}_{ref}}}{K}---\left(20\right)$

$q_f=C\left(H_{ref}-\frac{-CL+\sqrt{C^2{L}^2+\frac{(1+\mathrm{\α})}{\mathrm{\α}}{K}^2{H_s}^2+2{\mathrm{CKLH}}_{ref}}}{K}\right)---\left(21\right).$

2.如权利要求1所述的模拟海平面上升情况下海水入侵内陆边界的处理方法，其特征是，模型求解具体步骤如下：

(1)空间上采用积分有限差分法IFDM进行离散

首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散，对于任意单元n，单元体积为V_n，边界面为Γ_n，单元的质量平衡方程的积分形式如下：

$\frac{d}{dt}{\∫}_{V_n}{M}^{\mathrm{\κ}}\·{\mathrm{dV}}_n={\∫}_{{\mathrm{\Γ}}_n}{F}^{\mathrm{\κ}}\·{\mathrm{nd\Γ}}_n+{\∫}_{V_n}{q}^{\mathrm{\κ}}{\mathrm{dV}}_n---\left(8\right)$

式中，n为表面单元dΓ_n的单位法向量，指向控制单元体内为正；

引入体积平均值，有

${\∫}_{V_n}{M}^{\mathrm{\κ}}{\mathrm{dV}}_n=V_n{M}_n^{\mathrm{\κ}}---\left(9\right)$

${\∫}_{V_n}{q}^{\mathrm{\κ}}=q_n^{\mathrm{\κ}}{V}_n---\left(10\right)$

式中，为M^κ，q^κ在V_n上的平均值。

Γ_n上的面积分可近似为其所包含的各个表面A_nm的面积分的平均值之和，有

${\∫}_{{\mathrm{\Γ}}_n}{F}^{\mathrm{\κ}}\·{\mathrm{nd\Γ}}_n=\underset{m}{\Σ}{A}_{nm}{F}_{nm}^{\mathrm{\κ}}---\left(11\right)$

式中，m为与单元n相邻的所有单元，A_nm是单元n和m相邻的交界面，是F^κ在面A_nm上的内法线方向的平均值；

将式(9)、(10)和(11)代入到式(8)中，得到一组关于时间的一阶微分方程组

$\frac{{\mathrm{dM}}_n^w}{dt}=\frac{1}{V_n}\underset{m}{\Σ}{A}_{nm}{F}_{nm}^{\mathrm{\κ}}+q_n^{\mathrm{\κ}}---\left(12\right)$

(2)时间上采用一阶向后差分方法进行离散

对式(12)的时间微分采用一阶向后差分方法，得到任意单元的全隐式非线性方程组，见式(13)：

$R_n^{\mathrm{\κ},k+1}=M_n^{\mathrm{\κ},k+1}-M_n^{\mathrm{\κ},k}-\frac{\mathrm{\Δ}t}{V_n}\left\{\underset{m}{\Σ}{A}_{nm}{F}_{nm}^{\mathrm{\κ},k+1}+V_n{q}_n^{\mathrm{\κ},k+1}\right\}---\left(13\right)$

式中，引入了组分κ＝w,g,a的余量△t为时间步长，上标k和k+1表示两相邻的时间步长指标；右端的流量项和源汇项均采用新的时间步长值；

(3)Newton-Raphson迭代方法

对式运用Newton-Raphson迭代方法进行线性化，引入迭代指标p，对式(13)中的余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，得到包含3×NEL个方程的线性方程组，NEL为计算域内单元数，并且以两迭代步的增量为未知量；最后得到大型稀疏系数矩阵的线性方程组，如式(14)：

$-\underset{i}{\Σ}\frac{\∂R_n^{\mathrm{\κ},k+1}}{\∂x_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\mathrm{\κ},k+1}\left(x_{i,p}\right)---\left(14\right).$

3.如权利要求1所述的模拟海平面上升情况下海水入侵内陆边界的处理方法，其特征是，还包括模型验证步骤：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证。

4.如权利要求1所述的模拟海平面上升情况下海水入侵内陆边界的处理方法，其特征是，还包括如下分析步骤：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件，采用依靠水头的流量边界，分析模拟海平面上升情况下海水入侵的效果：

(1)入侵程度变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升过程中，用相对盐度等值线表征的海水入侵范围始终介于流量恒定和水头恒定中间，具体如下：当内陆边界为流量恒定时，随海平面的上升，海水入侵程度基本与初始情况保持不变；而当内陆边界为水头恒定时，随着海平面的上升，海水入侵范围明显增大；当内陆边界为依靠水头的流量边界时，随着海平面的上升，海水入侵范围有所增大，但入侵程度明显比水头恒定情况小；

(2)研究区域左右两侧等效水头差和左侧补给量的变化：内陆边界采用依靠水头的流量边界时，随着海平面的上升，研究区域左右两侧等效水头差和左侧补给量都在逐渐减小，但是内陆边界为流量恒定时，以上两项都基本保持不变，而水头恒定时以上两项减少的幅度要大于依靠水头的流量边界。