充气法治理承压含水层海水入侵的数值模拟方法

摘要

本发明涉及防止海水入侵技术领域。为治理承压含水层海水入侵，本发明采用的技术方案是：充气法治理承压含水层海水入侵的数值模拟方法，包括如下步骤：步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考虑温度对系统的影响；步骤二：模型求解；步骤三：模型边界条件确定：步骤四：模型验证：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证；步骤五：分析充气法治理承压含水层海水入侵的效果。本发明主要应用于防止海水入侵。

说明书

技术领域

本发明涉及防止海水入侵技术领域，具体讲，涉及充气法治理承压含水层海水入侵的数 值模拟方法。

背景技术

近年来，随着沿海地区经济的快速发展，地下水超采越来越严重，造成地下水位持续下 降，海水侵入地下含水层，形成海水入侵现象。我国海岸线长达1.8×104Km，出现海水入 侵的城市有十几座，入侵面积己超过10000Km²，由此引发的耕地土壤盐碱化、地下设备腐蚀、 水质恶化等问题严重阻碍了沿海地区的经济和社会发展。因此，开展海水入侵深入研究，提 出防治和减轻海水入侵的措施，具有重要的科学价值和社会意义。

目前，防治海水入侵的措施主要有：限制地下水开采量、迁移抽水井、人工回灌补给、 人工抽取咸水、修建地下帷幕。限制地下水开采量并不能从根本上解决海水入侵，有可能会 加剧供需水矛盾。迁移抽水井的成本较高，也可能面临其它问题，如建筑层或含水层的尺寸 条件。人工回灌需要大量的淡水资源，对于缺水地区，需兴建大量引水工程，势必增加成本 费用。人工抽取咸水能够减少海水入侵的程度，但主要问题在于如何处理抽取的咸水。地下 帷幕包括实体帷幕、水力屏障和地下充气墙等，实体帷幕和水力屏障的造价昂贵，可能会引 起水质恶化和污染，而修建地下充气墙的方法造价低，无需注水或泥浆即可形成挡水帷幕， 也不会引起二次污染，较为理想。然而，尽管充气法防治海水入侵在理论上是可行的，但针 对充气法防治效果展开试验的研究费用较为昂贵，且耗时耗力。

发明内容

为克服现有技术的不足，治理承压含水层海水入侵，为达到上述目的，本发明采用的技 术方案是：充气法治理承压含水层海水入侵的数值模拟方法，包括如下步骤：

步骤一：建立地下水气一液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考 虑温度对系统的影响，具体如下：

模型的基本控制方程为：

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

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

步骤二：模型求解：以TOUGH2/EOS7为工具，空间上采用积分形式的有限差分方法IFDM 进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用 Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：

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

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

$\frac{d}{\mathrm{dt}}{\int }_{V_n}{M}^{\kappa }{\mathrm{dV}}_n={\int }_{{\Gamma }_n}{F}^{\kappa }·{\mathrm{nd\Gamma }}_n+{\int }_{V_n}{q}^{\kappa }{\mathrm{dV}}_n---\left(\text{2}\right)$

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

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

${\int }_{V_n}{M}^{\kappa }{\mathrm{dV}}_n=V_n{M}_n^{\kappa }---\left(3\right)$

${\int }_{V_n}{q}^{\kappa }{\mathrm{dV}}_n=q_n^{\kappa }{V}_n---\left(4\right)$

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

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

${\int }_{{\Gamma }_n}{F}^{\kappa }·{\mathrm{nd\Gamma }}_n=\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }---\left(5\right)$

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

将式(3)、(4)和(5)代入到式(2)中，得到一组关于时间的一阶微分方程组

$\frac{{\mathrm{dM}}_n^{\kappa }}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }+q_n^{\kappa }---\left(6\right)$

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

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

$R_n^{\kappa ,\kappa +1}=M_n^{\kappa ,k+1}-M_n^{\kappa ,k}-\frac{\mathrm{\Delta t}}{V_n}\{\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa ,k+1}+V_n{q}_n^{\kappa ,k+1}\}---\left(7\right)$

式中，引入了组分κ=w，b，a的余量Δt为时间步长，上标k和k+l表示两相邻的 时间步长指标。其中，分别表示k、k+1时刻M^κ在单元体积V_n上的平均值，表示k+1时刻F^κ在面A_nm上的内法线方向的平均值，表示k+1时刻q^κ在单元体积V_n上的平 均值；右端的流量项和源汇项均采用新的时间步长值。

(3)Newton-Raphson迭代方法

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

$-\underset{i}{\Sigma }\frac{\partial R_n^{\kappa ,k+1}}{\partial x_i}{|}_P(x_{i,p+1}-x_{i,p})=R_n^{\kappa ,k+1}\left(x_{i,p}\right)---\left(8\right)$

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

(1)Dirichlet边界条件

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

①对于空气边界，其边界条件单元的相态为仅有气相状态，主要变量为孔隙气压力p_g、 参考盐水占相态的质量百分数X_b、空气占气相的质量分数和温度T，对于与大气接触的 边界上p_g＝p_atm，对于有空气超压作用的边界上，如人工充气墙边界上，p_g＝p_atm+Δp，其 中p_atm为大气压力，Δp为充气压力；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边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设 置与之相邻的边界条件单元，则与不透水面单元的流量交换为零；

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

步骤五：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件， 施加充气作用，分析充气法治理承压含水层海水入侵的效果：

(1)盐度变化：充气结束时，相对盐度等值线与不透水底板的交点向海水方向撤退，海 水入侵范围减小；

(2)水、气相压力及流场变化：水相压力和气相压力的变化规律基本一致，分为三个阶 段：迅速增大至峰值、由峰值迅速减小、缓慢减小趋近于相对稳定值；另外，由于毛细压力 的作用，观察点的水相压力略小于气相压力；在注入区及含水层顶部，水、气相压力增大的 较多，形成指向海水侧的水力梯度，从而驱替入侵的咸水退出含水层；

(3)空气损失变化：充气初始时期的空气损失较小，然后迅速上升，之后逐渐趋于平稳， 达到相对稳定值。

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

$M^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }{\mathrm{\phi S}}_{\beta }{\rho }_{\beta }---\left(9\right)$

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

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

$F_{\mathrm{adv}}^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }{F}_{\beta }---\left(10\right)$

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

$F_{\beta }=-k\frac{{\rho }_{\beta }{k}_{\mathrm{\gamma \beta }}\left(S_w\right)}{{\mu }_{\beta }}({▿p}_{\beta }-{\rho }_{\beta }g)---\left(11\right)$

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

液相的孔隙压力p_l为气相压力p_g和毛细压力p_c(负值)之和：

p_l＝p_g+p_c      (12)

关于毛细压力～饱和度的关系，采用Leverett模型来表征：

p_c＝p₀γ[1.417(1-S_e)-2.120(1-S_e)²+1.263(1-S_e)³]      (13)

式中，p₀为进气值，γ为表面张力，有效水饱和度S_e的表达式为 S_e＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为液相剩余饱和度，S_ls为液相最大饱和 度；

关于相对渗透率～饱和度的关系吗，采用Fatt and Klikoff模型来表征：

$k_{\mathrm{rl}}=S_e^3$(液相)      (14)

k_rg＝(1-S_e)³(气相)      (15)

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

$F_{\mathrm{diff}}^{\kappa }={-\mathrm{\phi \tau }}_0\underset{\beta }{\Sigma }{\tau }_{\beta }\left(S_{\beta }\right){\rho }_{\beta }{d}_{\beta }^{\kappa }{▿X}_{\beta }^{\kappa }---\left(16\right)$

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

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

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

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

relative permeability model：τ₀τ_β＝τ₀k_rβ(S_β)

Millington-Quirk Model：

constant diffusivity：τ₀τ_β＝τ₀S_β      (18)

Relative permeability model为相对渗透率模型，Millington-Quirk Model是由 Millington-Quirk提出的非饱和土壤导水率模型，constant diffusivity为常扩散系数模型。

本发明可带来如下效果：

(1)建立了地下水气-液二相流及溶质运移模型，可以用来模拟充气作用下，地下水中气 -液二相流动过程及其对咸-淡水交界面运移的影响。

(2)本发明避免了由于开展实验性研究的高费用和长工期问题，模拟结果合理可信，为 进一步展开充气法治理海水入侵的试验性研究提供了科学指导和分析依据。

附图说明

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

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

图3稳态下海水入侵的盐度分布。

图4稳态下地下水流速分布，水平流速为零的等值线。

图5稳态下的水压力分布。

图6稳态下的孔隙水压力分布。

图7充气作用的网格剖分及充气位置。

图8相对盐度为0.5的等值线与不透水底板的交点随时间的变化情况。

图9相对盐度为0.25、0.5和0.75的等值线与不透水底板的交点随时间的变化情况。

图10含水层中总盐度随时间的变化情况。

图11充气结束时水相和气相压力、水流和气流及气体饱和度的分布情况。

图12含水层中各观察点在充气过程中的水相和气相压力变化情况。

图13充气过程中的空气损失率变化情况。

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

具体实施方式

本发明采用的技术方案是：

步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考 虑温度对系统的影响。具体如下：

模型的基本控制方程为：

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

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

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

$M^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }\phi S_{\beta }{\rho }_{\beta }---\left(2\right)$

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

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

$F_{\mathrm{adv}}^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }{F}_{\beta }---\left(3\right)$

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

$F_{\beta }=-k\frac{{\rho }_{\beta }{k}_{\mathrm{\gamma \beta }}\left(S_w\right)}{{\mu }_{\beta }}({▿\rho }_{\beta }-{\rho }_{\beta }g)---\left(4\right)$

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

液相的孔隙压力p_l为气相压力p_g和毛细压力p_c(负值)之和：

p_l=p_g+p_c      (5)

关于毛细压力～饱和度的关系，采用Leverett模型来表征：

$p_c=p_0\gamma [1.417(1-S_e)-2.120{(1-S_e)}^2+1.263{(1-S_e)}^3]---\left(6\right)$

式中，p₀为进气值，γ为表面张力，有效水饱和度S_e的表达式为S_e=(S_l-S_lr)/(S_ls-S_lr)， S_l为液相饱和度，S_lr为液相剩余饱和度，S_ls为液相最大饱和度。

关于相对渗透率～饱和度的关系吗，采用Fatt and Klikoff模型来表征：

$k_{\mathrm{rl}}=S_e^3$(液相)           (7)

k_rg=(1-S_e)³(气相)      (8)

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

$F_{\mathrm{diff}}^{\kappa }=-{\mathrm{\phi \tau }}_0\underset{\beta }{\Sigma }{\tau }_{\beta }\left(S_{\beta }\right){\rho }_{\beta }{d}_{\beta }^{\kappa }▿X_{\beta }^{\kappa }---\left(9\right)$

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

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

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

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

τ₀τ_β=τ₀k_rβ(S_β)(relative permeability model))

(Millington-Quirk Model)

τ₀τ_β=τ₀S_β(constant diffusivity)         (11)

Relative permeability model为相对渗透率模型，Millington-Quirk Model是由 Millington-Quirk提出的非饱和土壤导水率模型，constant diffusivity为常扩散系数模型。

步骤二：模型求解：以TOUGH2/EOS7为工具，空间上采用积分形式的有限差分方法 (IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用 Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：

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

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

$\frac{d}{\mathrm{dt}}{{\int }_V}_n{M}^{\kappa }{\mathrm{dV}}_n={{\int }_{\Gamma }}_n{F}^{\kappa }·\mathrm{nd}{\Gamma }_n+{{\int }_V}_n{q}^{\kappa }{\mathrm{dV}}_n---\left(12\right)$

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

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

${{\int }_V}_n{M}^{\kappa }{\mathrm{dV}}_n=V_n{M}_n^{\kappa }---\left(13\right)$

${{\int }_V}_n{q}^{\kappa }{\mathrm{dV}}_n=q_n^{\kappa }{V}_n---\left(14\right)$

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

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

${{\int }_{\Gamma }}_n{F}^{\kappa }·{\mathrm{nd\Gamma }}_n=\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }---\left(15\right)$

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

将式(13)、(14)和(15)代入到式(12)中，得到一组关于时间的一阶微分方程组

$\frac{{\mathrm{dM}}_n^{\kappa }}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }+q_n^{\kappa }---\left(16\right)$

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

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

$R_n^{\kappa ,k+1}=M_n^{\kappa ,k+1}-M_m^{\kappa ,k}-\frac{\mathrm{\Delta t}}{V_n}\{\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa ,k+1}+V_n{q}_n^{\kappa ,k+1}\}---\left(17\right)$

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

(3)Newton-Raphson迭代方法

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

$-\underset{i}{\Sigma }\frac{\partial R_n^{\kappa ,k+1}}{{\partial x}_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\kappa ,k+1}\left(x_{i,p}\right)---\left(18\right)$

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

(1)Dirichlet边界条件

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

①对于空气边界，其边界条件单元的相态为仅有气相状态，主要变量为p_g，X_b，T，对 于与大气接触的边界上p_g＝p_atm，对于有空气超压作用的边界上，如人工充气墙边界上， p_g＝p_atm+Δp；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边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设 置与之相邻的边界条件单元，则与不透水面单元的流量交换为零。

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

步骤五：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件， 施加充气作用，分析充气法治理承压含水层海水入侵的效果：

(1)盐度变化：充气结束时，相对盐度等值线与不透水底板的交点向海水方向撤退，海 水入侵范围减小；

(2)水、气相压力及流场变化：水相压力和气相压力的变化规律基本一致，可分为三个 阶段：迅速增大至峰值、由峰值迅速减小、缓慢减小趋近于相对稳定值。另外，由于毛细压 力的作用，观察点的水相压力略小于气相压力。在注入区及含水层顶部，水、气相压力增大 的较多，形成指向海水侧的水力梯度，从而驱替入侵的咸水退出含水层；

(3)空气损失变化：充气初始时期的空气损失较小，然后迅速上升，之后逐渐趋于平稳， 达到相对稳定值。

下面结合附图，针对充气法治理承压含水层海水入侵的数值模拟方法进行具体说明，其 步骤如下：

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

假定含水层系统处于恒温状态(T＝25℃)。

表1模型相关参数取值

模型的基本控制方程为：

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

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

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

$M^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }{\mathrm{\phi S}}_{\beta }{\rho }_{\beta }---\left(2\right)$

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

②平流流量的表达式为：

$F_{\mathrm{adv}}^{\kappa }=\underset{\beta }{\Sigma }{X}_{\beta }^{\kappa }{F}_{\beta }---\left(3\right)$

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

$F_{\beta }=-k\frac{{\rho }_{\beta }{k}_{\mathrm{\gamma \beta }}\left(S_w\right)}{{\mu }_{\beta }}({▿p}_{\beta }-{\rho }_{\beta }g)---\left(4\right)$

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

液相的孔隙压力p_l为气相压力p_g和毛细压力p_c(负值)之和：

p_l＝p_g+p_c            (5)

关于毛细压力～饱和度的关系，采用Leverett模型来表征：

$p_c=p_0\gamma [1.417(1-S_e)-2.120{(1-S_e)}^2+1.263{(1-S_e)}^3]---\left(6\right)$

式中，p₀为进气值，γ为表面张力，有效水饱和度S_e的表达式为S_e＝(S_l-S_lr)/(S_ls-S_lr)， S_l为液相饱和度，S_lr为液相剩余饱和度，S_ls为液相最大饱和度。

关于相对渗透率～饱和度的关系吗，采用Fatt and Klikoff模型来表征：

$k_{\mathrm{rl}}=S_e^3$(液相)            (7)

k_rg＝(1-S_e)³ (气相)             (8)

③扩散流量的表达式为：

$F_{\mathrm{diff}}^{\kappa }=-{\mathrm{\phi \tau }}_0\underset{\beta }{\Sigma }{\tau }_{\beta }\left(S_{\beta }\right){\rho }_{\beta }{d}_{\beta }^{\kappa }{▿X}_{\beta }^{\kappa }---\left(9\right)$

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

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

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

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

${\tau }_0{\tau }_{\beta }={\tau }_0{k}_{\mathrm{r\beta }}\left(S_{\beta }\right)\left(\mathrm{relativepermeability}\mathrm{mode}l\right)$

${\tau }_0{\tau }_{\beta }={\phi }^{1/3}{S}_{\beta }^{10/3}(\mathrm{Millington}-\mathrm{QuirkModel})$

${\tau }_0{\tau }_{\beta }={\tau }_0{S}_{\beta }\left(\mathrm{cons}\tan \mathrm{tdiffusivity}\right)---\left(11\right)$

Relative permeability model为相对渗透率模型，Millington-Quirk Model是由 Millington-Quirk提出的非饱和土壤导水率模型，constant diffusivity为常扩散系数模型。

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

$R_n^{\kappa ,k+1}=M_n^{\kappa ,k+1}-M_n^{\kappa ,k}-\frac{\mathrm{\Delta t}}{V_n}\{\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa ,k+1}+V_n{q}_n^{\kappa ,k+1}\}---\left(12\right)$

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

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

$-\underset{i}{\Sigma }\frac{{\partial R}_n^{\kappa ,k+1}}{{\partial x}_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\kappa ,k+1}\left(x_{i,p}\right)---\left(13\right)$

(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℃；其中，ρ_s为海水密度，p_atm为 大气压力，等于1.013×10⁵Pa。初次计算的初始条件为：所有单元液相饱和，压力均为大气 压力与静水(淡水)压力之和即p_g＝p_atm+ρ_wg(100-z)，X_b＝0.0，

以咸、淡水静力平衡为初始条件，施加充气作用后，系统的右侧、顶部和底部的边界条 件同稳态状态下，对于左侧边界，p_g＝ρ_wg(102.5-z)，X_b＝0.0，和T＝25 ℃；对于充气井边界，主要变量p_g＝p_atm+9.0×10⁵，X_b＝0.0，和T＝25℃。

(4)模型验证：

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

模拟过程中，可得到静力平衡情况下滨海含水层的相对盐度、地下水流速、水压力水头 和孔隙水压力水头的分布情况，详见图3～图6。由图可见，含水层在淡水和咸水的共同作用 下达到静力平衡，海水沿临海(右)边界上的下半部分大约0～48m的范围内流入含水层，与 地下淡水混合形成楔形的咸水入侵体，并在淡水压力作用下沿临海边界上的上半部分大约 49～100m的范围流出，形成一个海水入侵环流。在含水层的右上角处地下水流速最大，约为 1.59×10^-6m/s。地下水流速为零的等值线反映了海水的入侵范围，其与不透水底板的交点约 为x＝92.0m，即海水入侵含水层的最远距离约为108.0m。

(5)根据地下水气-液二相流及溶质运移模型，分析充气法治理承压含水层海水入侵的 效果：

以咸-淡水静力平衡情况作为初始条件，利用钻孔施加充气，当土体水力传导率的范围在 10^-8～10^-4m/s时，可以采用充气法来阻止地下水侵入含水层，有表1数据显示，计算域满足条 件。

根据初始的盐度分布情况(图3)，海水入侵淡水层的最远距离约为108.0m，将充气井布 置在距离海水边界的水平距离为120m处，直径取1.0m。由于空气的密度小于水的密度，空 气进入含水层后将上升，因此将注气位置布置在含水层的下半部，选取x＝80.5m，15≤z≤30m， 如图7。根据压气法的相关经验，只有当钻孔的注气压力大于注气范围内的最大孔隙水压力 时，采用充气法驱替咸水才是完全可能的。有图6可知，注气区域的最大孔隙水压力为 84.6mH₂O(≈8.3×10²kpa)，因此模型采用的注气压力为9.0×10²kpa。注气区的钻孔位置、 网格单元如图7所示，其顶部是不透水边界，不透水区域以上的部分属性与含水层相同，充 气时间设置为一年(365天)。从三维模型看，钻孔占地下水气流的空间很小，产生的阻碍 可以忽略。计算过程中，液相剩余饱和度S_lr取值0.15，表面张力γ在25℃条件下为0.072N/m， 进气值$p_0=\sqrt{\phi /k}=540\mathrm{Kpa}.$

为了研究充气过程中水相和气相压力随时间的变化情况，取5个观察点进行研究，如图 7。其中①、②、③三个点的水平位置相同，用于研究各相压力随垂向位置的变化规律，③、 ④和⑤的垂向坐标相等，用于研究各相压力随水平位置的变化规律。

①盐度变化

图8为相对盐度0.5的等值线与不透水底板的交点在充气过程中随时间的变化过程。图9 为相对盐度为0.25、0.5和0.75的等值线与不透水底板的交点在充气过程中的瞬态变化，由 图可见，相对盐度为0.25、0.5和0.75的等值线分别向海水方向迁移了53.1m、51.0m、43.1m， 其平均速度为0.145m/d、0.140m/d、0.118m/d。含水层的整体盐度变化如图10，在整个时间 段内从176.5mg/l减少到了47.2mg/l。结果表明充气作用下，海水入侵的范围逐渐减少。

②水、气相压力及流场变化：

图11a-b为充气结束时，水相和气相压力及水流和气流的分布情况，由图可见，由于空 气的注入引起了各相压力的增大，在含水层顶部和空气注入区附近增大较多；由于空气的密 度较小，空气进入含水层后先垂直向上流动，在含水层不透水顶板处聚集，并向两侧扩散； 水相压力的变化形成了由空气注入区指向海水边界的水力梯度，使得含水层中的地下水向海 水边界流动，从而驱替入侵的咸水退出含水层。图11c所示为充气结束时，气体饱和度的分 布情况。另外，由于空气的注入，土体由饱和过渡为非饱和产生了毛细压力，造成水相压力 增大值均略小于气相压力增大值。

图12为各个观察点在充气过程中水相和气相压力的变化情况。由图可见，水相压力和气 相压力的变化规律基本一致，可分为三个阶段：迅速增大至峰值、由峰值迅速减小、缓慢减 小趋近于相对稳定值。图12a所示的点①、②、③的水平坐标相同，气相和水相压力增大值 的峰值随着垂向位置的升高依次增大，且含水层顶部的压力增大值的相对稳定值较大。图12b 所示的点③、④和⑤的纵坐标相同，均紧靠近含水层不透水顶板，随着距离充气井的水平距 离的增大，各相压力增大值的峰值和相对稳定值逐渐减小。

③空气损失变化

如图13所示，空气损失速率由初始时的q＝0.23m³/min逐渐增大，最终在t＝90d时达到稳 定为q＝0.824m³/min；在注气时间(一年)内，总空气损失量大约为Q＝2.61×10⁷m³。然而充 气结束后，作为水源地，需从含水层底部抽水，如此势必打破新建立的气体平衡，导致空气 损失和注气成本增加。因此，开采井应布置在空气扩散的外围。

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

(1)针对滨海地区淡水资源比较匮乏的特点，采用数值模拟的方法研究了充气法治理承 压含水层海水入侵的效果。首先利用地下水气-液二相流及溶质运移模型模拟分析了经典的 Henry问题的咸-淡水静力平衡情况，通过与以前文献结果的对比，验证了所采用模型的有效 性。然后基于咸-淡水静力平衡情况，模拟分析了充气作用对含水层盐度、各相压力及空气损 失的影响。

(2)对模拟结果的分析表明，由于空气的密度小于地下水的密度，空气注入含水层后向 上流动，并聚集在含水层顶板进而向两侧流动，因此水相和气相压力在注入区及含水层顶部 增大的较多，并形成了指向海水侧的水力梯度，驱替入侵的咸水向含水层外流动，从而减小 了海水入侵的范围。

(3)含水层各相压力的增大值及空气损失均在充气作用一段时间后达到某一相对稳定 值，说明充气法防治海水入侵的效果及单位成本是比较稳定的；但是充气法应用于实际工程 还需要很多研究工作，包括充气法的影响因素、成本及其可能存在的风险研究等。