水力屏障和截渗墙在海水入侵防治中的数值模拟研究

吕盼盼，宋 健，吴剑锋，吴吉春

（南京大学地球科学与工程学院/表生地球化学教育部重点实验室，江苏 南京 210023）

海水入侵是由淡水和海水之间的密度差异驱动的咸-淡水界面向陆地方向迁移并达到平衡的自然现象[1−2]。但是，随着社会经济的发展以及气候变化等外部因素，滨海地区的地下水超采、土地利用类型的改变以及海平面上升等现象破坏了原有的咸淡水界面的动态平衡，进一步加剧了海水入侵的程度[3−7]。海水入侵已在全球范围内发展，造成了滨海地区可利用淡水资源的减少以及开采井的报废等不利影响，极大地制约了滨海地区的社会经济发展[8−9]。在以往研究中，依据滨海含水层中过渡带的类型可将海水入侵研究分为突变界面模型和过渡带模型，其中突变界面模型假设海水和淡水为非混溶流体，分界面固定，但是在实际场地条件下海水和淡水具有可混溶性，因此过渡带模型更符合海水入侵的规律[10−11]。在人类活动引发海水入侵的滨海地区，需要耗费大量经济成本构建复杂的工程措施防止滨海含水层中可利用淡水资源的进一步污染，并且治理的周期十分漫长。

当前，国内外学者已提出多种方法应对海水入侵问题，如控制地下水开采量、人工补给地下水、建立地下水截渗墙等[12−16]。如Ebeling等[17]利用SEAWAT和FloPy建立了一个二维变密度海水入侵模型，研究混合水力屏障方案的可行性和最优管理策略，并与单个抽注水力屏障进行对比分析。Kaleris等[18]利用数值模拟方法研究了截渗墙的深度、距海岸线的距离、截渗墙渗透率以及含水介质的各向异性对海水入侵的影响，并评估了其对海岸附近地下水开采的保护作用。Abdoulhalik等[19]提出了结合不透水截渗墙和半透水地下坝的混合物理屏障作为控制海水入侵的新屏障系统。该系统可引起淡水推动咸水朝海岸线方向向上抬升，可显著减小滨海含水层咸水体的入侵范围。

本文从建立水力屏障和物理屏障（截渗墙）2个角度，构建基于SEAWAT-2000程序的室内二维砂箱试验数值模型，并从物理屏障（截渗墙）的角度构建基于山东龙口地区某典型剖面的海水入侵数值模型。结果揭示了不同水力或物理屏障布设方式对咸淡水界面运移规律的影响，可为场地条件下滨海含水层海水入侵防治中的水力屏障和截渗墙布设提出合理化建议与决策。

1 研究方法

1.1 海水入侵模型

基于SEAWAT-2000程序构建数值模型模拟二维砂箱试验以及山东龙口地区典型剖面在不同注水井、截渗墙布设方案下的咸水入侵过程[20]。SEAWAT是由美国地质调查局开发的用于模拟变密度地下水水流和盐分迁移的模型。水流方程与溶质运移方程为：

式中：ρ—流体密度；

Kfi—主轴方向的等效淡水渗透系数；

ρf—淡水密度；

Z—计算点高程；

hf—等效淡水水头；

Sf—等效淡水单位储水系数；

t—时间；

θ—有效孔隙度；

C—溶质浓度；

ρs—源（汇）流体密度；

qs—源（汇）单位体积流量；

Ck—溶质k的溶解浓度；

Dij—水动力弥散张量；

vi—渗流速度；

—源（汇）流中溶质k的浓度；

Rn—化学反应项。

该程序设计的基本原理是：以MODFLOW模拟地下水流过程，MT3DMS模拟溶质运移过程，并将两者耦合求解变密度条件下地下水流和溶质运移的耦合模型。SEAWAT已成为目前研究变密度地下水问题最常用的模型之一，并表现出较好的稳定性与有效性[21−25]。

1.2 海水入侵管理措施

为探究注水井的井位、流量及截渗墙的位置、贯穿深度等参数对咸淡水界面运移规律的影响，本文采用SEAWAT模型，模拟咸淡水界面达到自然条件下的稳定状态，设置注水井或截渗墙的不同工况。注水井布设于自然稳定状态下盐水楔的外部，为仅底部渗水非完整井的点源补给。截渗墙沿地表向地下贯穿且假设为不透水介质，模拟均为单井（或单个截渗墙）。注水井系列模拟方案的变量仅为注水井位置或补给强度，截渗墙系列模拟中变量仅为墙体水平位置或贯穿深度。在加入注水井或截渗墙后，咸淡水界面达到新的稳定状态，此时记录盐水楔前锋到达的位置，并对比初始状态下盐水楔前锋位置，计算表征盐水楔驱退程度的海水入侵回退系数（R）[26]，即为注水井和截渗墙管理措施的有效性。回退系数（R）为：

式中：L0—初始稳定条件下盐水楔前锋到海岸线的 距离；

L—设置注水井或截渗墙措施后再次达到稳定时盐水楔前锋到海岸线的距离。

2 二维砂箱数值模型

2.1 试验模型

模拟区为二维砂箱，可概化为二维、均质、各向同性的潜水含水层（图1）。模拟范围长90.0 cm，宽1 cm，高41.6 cm，为典型的二维剖面。将无外界干预下咸淡水界面达到自然平衡的状态设为初始条件。模拟区域顶部为零通量边界；底部为隔水边界；左侧设置定水头、定浓度边界，边界浓度以TDS表征设为33.6 g/L；右侧设置定水头边界，初始浓度为0 g/L。理想砂箱实验的概念模型见图1。采用SEAWAT-2000程序对模型求解，利用有限差分方法将其在空间上剖分成41列和91行的正方形网格，差分网格间隔为Δx=Δz=1 cm。在天然咸淡水界面达到稳定状态后，设定不同工程措施的管理期为0.5 h。该理想算例主要参数取自文献[26]，见表1。

图1 理想砂箱概念模型示意图Fig.1 Conceptual model of the ideal sandbox

表1 理想砂箱试验地下水数值模型主要参数[26]Table 1 Input parameters for the numerical model of the ideal sandbox

2.2 模拟结果与分析

2.2.1 注水井工程措施评价

（1）注水井流量

为探究注水井流量对咸淡水界面的影响，在初始咸水楔前锋附近设置注水井（图1，假设仅在注水井底部进水），并改变流量进行不同补给情景下的数值模拟。稳定后记录盐水楔前锋位置L，根据式（3）计算不同补给情景下的海水入侵回退系数。海水入侵回退系数（R）与补给流量（Q）的关系见图2。

图2 海水入侵回退系数R与补给流量Q的关系Fig.2 Repulsion ratio with respect to recharge quantity

模拟结果表明，在相同补给点位下，海水入侵回退程度随补给强度的增大而增大，回退系数R与补给流量Q呈非线性正相关关系。当Q大于0.005 m3/d时，随着Q的增大R产生的增量明显减小，且缓慢趋近于1。研究结果表明随着注水井流量的增大，淡水体水头不断抬升，造成滨海含水层中地下水向海洋的排泄量增加，引起咸淡水界面向海岸线方向移动，实现了驱退海水入侵前锋的目的。

（2）注水井位置

为探究不同注水井位置对咸淡水界面的影响，识别最大化海水入侵回退距离的最优注水井位置，在初始咸水楔外部不同位置以20 cm×10 cm的矩形网格节点形式均匀设置一系列相同补给强度（4.2×10−5m3/d）注水井，并在初始盐水楔前锋附近局部加密了注水井分布。图3表示海水入侵回退系数R与注水井（仅底部进水）位置的关系图。

如图3所示，位于盐水楔前锋左上角的W3井（x=0.40 m，z= 0.05 m）的回退系数达到最大值，约为21%，表明W3井补给的淡水形成了一个从含水层下部产生的水力屏障，达到了最优的驱退海水入侵前锋的管理效果。盐水楔前锋处W1井的位置达到19%的回退系数，此时产生的补给强度仍然有效，但没有W3井的效果明显，距离W3井右边10 cm的W4井也能达到约14%的回退系数。总体上，注水点距W3井越远，海水入侵回退系数越小。

图3 二维砂箱试验模型海水入侵回退系数R与注水井位置关系Fig.3 Repulsion ratio with respect to recharge well position

模拟结果表明：注水井注入淡水最有效的位置是靠近W3井的位置。如果将淡水补给直接施加在W3井附近，则产生的效果最佳，但当施加的淡水补给远离W3井的位置时效果逐渐减弱；表面补给（如W15、W16、W17、W18）在控制海水入侵方面不如盐水楔前锋附近的注入措施更有效。对于完整井的补给井，同样可以确定为W3、W8、W12、W16构成最优组合，最有利于驱替海水入侵的前锋。由于注水井工程成本随着深度的增加而增大，为节约成本，可适当缩减井的贯穿深度，将注水井布设在靠近W3处。对于表面补给，最有效的位置应在W3井的正上方。

2.2.2 截渗墙工程措施评价

为研究截渗墙对防治海水入侵的影响，模型达到自然条件下的稳定状态后，加入不透水截渗墙并改变其水平位置和贯穿深度完成系列模拟，并计算不同情景模式下的海水入侵回退系数R。贯穿深度以无量纲的贯穿率h/H表示，截渗墙的水平位置以无量纲的1-Xb/L0表示。其中：h表示截渗墙贯穿深度，H表示含水层厚度，Xb表示截渗墙的水平位置，L0为初始状态下盐水楔前锋位置。

（1）截渗墙水平位置的影响

基于不同水平位置截渗墙条件下的模拟结果，计算海水入侵回退系数R。从图4（a）可以看出，当1-Xb/L0<0时，即截渗墙水平位置位于初始盐水楔前锋靠岸一侧时，海水入侵回退系数R值为负数，盐水楔进一步向淡水区推进。当1-Xb/L0=0时，即截渗墙水平位置恰位于初始盐水楔前锋位置时，海水入侵回退系数R值接近0。当1-Xb/L0>0时，即截渗墙水平位置位于初始盐水楔前锋靠海一侧时，海水入侵回退系数R值为正数，盐水楔形体向海岸线方向回退。因此，在不同贯穿率下，截渗墙水平位置越靠近海岸线，海水入侵回退系数越大，截渗墙工程效果越显著。在不同贯穿率下回退系数R与1-Xb/L0呈线性关系。

图4 二维砂箱试验模型回退系数R与截渗墙水平位置1-Xb/L0和截渗墙贯穿率h/H的关系Fig.4 Repulsion ratio（R）with respect to horizontal barrier location expressed in（1-Xb/L0）and respect to penetration ratio（h/H）

（2）截渗墙贯穿深度的影响

基于不同贯穿深度截渗墙条件下的模拟结果，计算海水入侵回退系数R。从图4（b）可以看出，没有加入截渗墙时，盐水楔前锋位置没有变化，初始值R=0。当1-Xb/L0>0时，随着h/H增加，即截渗墙向含水层底部贯穿时，回退系数R增加；当贯穿率h/H接近1时，盐水楔前锋被阻挡在截渗墙位置之外，实现盐水楔最大回退值。当1-Xb/L0<0时，即截渗墙位于盐水楔靠岸一侧时，回退系数R为负值，表明随着截渗墙贯穿深度增加，截渗墙阻止了淡水水流对海水入侵的排斥作用，导致盐水楔前锋进一步向淡水体推进，加剧了海水入侵程度。因此，截渗墙贯穿深度对控制海水入侵是否有效取决于截渗墙水平位置。在不同的截渗墙水平位置下，R与h/H呈非线性关系。

（3）截渗墙水平位置与贯穿深度综合影响

根据多组模拟结果绘制在不同水平位置及贯穿深度下回退系数（R）变化的等值线图，见图5。

从图5可以看出，截渗墙水平位置越靠近海岸线一侧，即1-Xb/L0越大，海水入侵回退系数越大。且当1-Xb/L0大于0时，截渗墙贯穿深度越大，海水入侵回退系数越大。当1-Xb/L0小于0时，截渗墙贯穿深度越大，海水入侵回退系数越小，且为负值。因此，当截渗墙水平位置靠近海岸线、贯穿深度较大时达到抵挡海水入侵最佳效果。

图5 二维砂箱试验模型不同水平位置与贯穿深度截渗墙下的回退系数R的分布Fig.5 Repulsion ratio distribution based on location and penetration of the barrier

3 龙口滨海含水层截渗墙剖面模型

3.1 海水入侵数值模型

选取位于山东龙口市黄水河附近靠近渤海湾的某典型剖面AB为研究对象（图6）。该研究区受海水入侵影响且在距海岸线1 km处建有地下截渗墙，实际条件与前期研究的理想算例较为符合。选取剖面长度为1 500 m，厚度为20 m。含水层岩性主要为第四系松散沉积物，由黄水河河相沉积物组成。第四系下伏古近系含煤地层，岩性为砂砾岩、泥灰岩、油页岩、黏土岩及褐煤，透水性及富水性差，构成研究区相对隔水层。研究区农业灌溉以开采地下水为主要水资源利用模式之一，同时也是影响该区域地下水动态的主要因素。

图6 山东龙口地区研究区域示意图Fig.6 Location of the study area in Longkou of Shandong Province

将研究区含水层概化为垂向二维、均质、各向同性的潜水含水层（图7）。将未构建截渗墙时滨海含水层中咸淡水界面达到的稳定状态作为初始条件。含水层上部边界为潜水面；下部边界为隔水边界并且相对隔水底板存在一定的地形坡度；北部海岸线边界设为定水头、定浓度边界；南部边界设置定水头边界。在研究区内距海1.1 km处设置一口底部进水部位长5 m的非完整抽水井，模拟农业灌溉抽水对滨海地下水系统的影响，抽水量设为50 m3/d。采用SEAWAT-2000程序对模型进行求解，利用有限差分方法将其在空间上剖分成151列、1行、20层，差分网格的间隔为Δx=10 m、Δy=Δz=1 m，模拟期为40 d。模型主要参数见表2。

图7 研究区典型剖面概念模型示意图Fig.7 Conceptual model of the typical profile in Longkou

表2 山东龙口地区典型剖面地下水数值模型参数Table 2 Input parameters used in the numerical model for the typical profile in Longkou

3.2 模拟结果与分析

3.2.1 截渗墙工程措施评价

（1）截渗墙水平位置的影响

基于不同水平位置截渗墙条件下的模型结果计算海水入侵回退系数R。回退系数R与截渗墙水平位置1-X/Xw的关系图，见图8（a）。截渗墙的水平位置以无量纲的1-X/Xw表示，贯穿深度以无量纲的贯穿率h/H表示。其中：X表示截渗墙的水平位置，Xw表示抽水井的水平位置，h表示截渗墙贯穿深度，H表示含水层厚度。

模拟结果表明：海水入侵的回退系数与截渗墙的水平位置呈正相关关系，截渗墙水平位置越靠近海岸线一侧效果越好。当截渗墙水平位置位于抽水井靠海岸线一侧时（1-X/Xw>0），海水入侵回退系数为正，截渗墙对盐水楔有驱退效果。当截渗墙水平位置位于抽水井靠岸的一侧时（1-X/Xw>0），海水入侵回退系数为负数，截渗墙起反作用，海水入侵进一步加剧。

（2）截渗墙贯穿深度的影响

基于不同贯穿深度截渗墙条件下的模型结果计算了盐水侵入楔的回退系数R，如图8（b）所示。

模拟结果表明：当截渗墙水平位置位于抽水井靠海一侧时（1-X/Xw>0），截渗墙贯穿深度越大，海水入侵回退系数越大，截渗墙效果越好。此时截渗墙有效阻挡了抽水井抽水引起的淡水水头降低所带来的影响。当截渗墙水平位置位于抽水井的靠岸的一侧时（1-X/Xw>0），截渗墙阻挡了内陆淡水区的补给，因此截渗墙贯穿深度越大，海水入侵回退系数越小，截渗墙效果越差，甚至起到反作用。

（3）截渗墙水平位置与贯穿深度综合影响

根据多组模拟数据，得到截渗墙不同水平位置（X）及不同贯穿深度（Z）对海水入侵回退系数（R）的等值线图，见图9。当截渗墙位于注水井靠海一侧，且离海较近、贯穿深度较大时，驱退海水入侵的效果最佳。

图9 龙口滨海含水层截渗墙情况下R值等值线图Fig.9 Repulsion ratio distribution based on location and penetration of the barrier

3.2.2 实际管理措施分析

山东龙口地区所建截渗墙位于距海岸线1.0 km处，贯穿深度为18 m。依据实际截渗墙工程情况进行海水入侵的数值模拟，并与未设置截渗墙的初始状态进行对比，模拟结果如图7。研究表明龙口地区现有截渗墙相对初始状态可以有效降低海水入侵程度，海水入侵前锋向海岸线方向后退了约230 m，回退系数达23.7%。因此，现有截渗墙工程措施对保护黄水河滨海地区的淡水资源具有重要的作用。但是，对该地区截渗墙工程措施的评价分析表明，为进一步提升截渗墙工程措施的治理效果，在后续滨海地下水管理中应重新考虑优化截渗墙的水平位置与贯穿深度，适当向海岸线方向扩展以达到更好的海水入侵回退效果。

4 结论

（1）注水井工程管理措施对海水入侵的影响与补给流量及井位有关。海水入侵的回退系数与注水井流量在一定范围内成正相关关系。注入井位在靠近盐水楔的位置为最佳淡水注入点，并形成以该点为中心的最优注入井布局。表面补给最有效的位置是最佳注水点正上方，滨海含水层的表面补给在控制海水入侵时，不如注入井更有效。

（2）截渗墙工程管理措施对海水入侵的影响与水平位置及贯穿深度有关。截渗墙水平位置越靠近海岸线一侧，海水入侵回退系数越大，截渗墙效果越好。一般情况下截渗墙贯穿深度越大，海水入侵回退系数越大，海水入侵防治效果最优。

（3）当截渗墙水平位置超过某一界限后（理想算例中该界限为盐水楔前锋，山东龙口实例中该界限为抽水井位置），截渗墙贯穿深度越大，海水入侵回退系数越小且为负值，截渗墙对海水入侵防治起反作用。当截渗墙水平位置靠近海岸线、贯穿深度较大时，该工程措施可以达到治理海水入侵最优效果。

（4）山东龙口黄水河地区的截渗墙工程管理措施对防治当地海水入侵具有重要的作用。但是数值模拟结果分析表明，在后续滨海地下水管理中应进一步优化截渗墙空间布局。

研究结果可为场地条件下滨海地下水管理提供实际的管理措施建议。但本文实例模型未考虑含水介质非均质性以及降雨入渗等边界条件的动态变化对管理措施的影响。未来应构建流域尺度下龙口地区的海水入侵数值模型，基于数值模型并采用优化技术对海水入侵严重区进行注水井、截渗墙等海水入侵防治工程的优化管理，提出经济、高效的滨海地下水管理建议。