在天然条件下, 滨海含水层海岸带咸水和淡水之间始终保持着一种动态的平衡状态。由于全球平均海平面上升和区域性降水量在某些年份内偏小, 特别是不合理开采海岸带地下水资源造成滨海含水层中地下水位下降, 海岸带原有的咸淡水之间的平衡关系被打破, 致使咸淡水界面向内陆方向运移。海水入侵是指在开发利用地下水资源过程中, 造成海水向内陆滨海含水层迁移的现象。国内外许多沿海地区都出现了程度不同的海水入侵现象, 不合理的开采井布局和地下水的过量开采通常是导致海水入侵的直接原因, 此外, 连续多年的干旱少雨、上游兴建地表水拦蓄工程、地表建筑物覆盖面积扩大导致滨海含水层地下水补给量减少等因素也加剧海岸带的海水入侵^[1]。

滨海含水层海岸带咸水与淡水之间界面的位置、形状、运移及其预测是研究海水入侵的核心内容。咸淡水之间关系的研究模型可根据是否考虑其间的混溶问题分为两类: 过渡带模型和突变界面模型。前者考虑到海岸带中咸、淡水之间的混溶, 可将其视为一个水动力弥散问题来研究。后者是在过渡带相对于整个含水层的厚度很窄的前提下, 将咸淡水之间的接触带近似地处理为一个突变界面^[2]。运用咸淡水突变界面模型或者过渡带模型定量研究海水入侵问题, 对不同条件下海岸带咸淡水之间突变界面或过渡带的发展趋势进行预测, 可以为海水入侵的防治与合理开采利用和保护滨海含水层地下水资源提供科学依据^[3]。

1.   咸淡水界面的解析研究

海岸带咸淡水界面的解析研究已有100多年的历史。Ghyben和Herzberg分别在1889和1901年各自独立提出经典的Ghyben-Herzberg公式, 用来确定咸淡水之间突变界面位置。假定海岸带咸淡水之间处于静水平衡状态, 在海岸处潜水面与咸淡水界面相交, 且不存在淡水出口。根据咸淡水界面上同一点受到的淡水压力和咸水压力相同, 可建立如下关系:

式中: M——咸淡水界面上任意点距平均海平面的深度;

h_f——在同一垂直线上的潜水面相对于平均海平面的高度;

ρ_s和ρ_f——咸水和淡水的密度;

g——重力加速度。

由式(1)可得Ghyben-Herzberg公式^[4]:

式(2)仅依赖潜水位h_f来确定界面的深度, 简单易行。当h_f=0时M=0, 即式(2)描述的咸淡水界面在海岸处不存在淡水出口, 这种情形与实际海岸带不相符。即便这样, Ghyben-Herzberg公式在区域性研究中仍不失为一种简单的方法^[1]。Hubbert对存在突变界面, 且密度也不同的两种流体的水头与观测点之间的位置关系进行了研究, 推导出两种流体水头和该点位置之间关系的Hubbert公式^[5], Hubbert公式是正确的严格的描述界面位置的数学公式, 然而由于界面的位置本身并不能事先知道, 从而使得Hubbert公式受到限制而难以在实际中得到应用^[1]。

此后, 海岸带咸淡水界面的解析模型的研究取得了较快的发展。Henry发表了正交于海岸线的承压稳定流动条件下垂直断面上盐分浓度分布的解析解(后来称为Henry模型), 并注意到海水的环流现象^[6]。Bear等对承压含水层中淡水与咸水之间的界面变化进行了探讨, 给出了界面的解析解, 并基于此假设, 给出了抽水导致界面升锥问题的解析解, 结果表明, 用近似解描述界面运动具有足够的精度^[7]。Dagan等基于小扰动法的线性化近似解, 提出了解决移动界面的问题所需的精确方程, 给出了存在小偏差的初始或平均稳定界面的解析解。同时考虑了二维和三维两种情况, 通过沙箱模型中的实验验证, 其分析结果均在有效的范围内^[8]。Youngs基于Ghyben-Herzberg公式, 给出了以淡水边界为定水头边界的潜水含水层模型以及圆岛模型中咸淡水界面上部揭穿整个淡水带的完整井抽水的最大抽水量计算公式^[9]。Bear在《多孔介质流体动力学》和《地下水水力学》中对界面在稳定与移动情况下的近似解进行了详细的论述, 分析了井在界面上部抽水时所引起的升锥问题, 并指出应用Ghyben-Herzberg公式估算出的咸淡水界面位置要比实际深度偏小, 且越靠近海岸其误差越是明显^[10~11]。Mercer等采用数学模型并以Dupuit假设为基础对界面的移动进行了模拟。结果表明, 该模型能与解析解较好地吻合, 但在复杂的条件下, 利用Dupuit近似求解可能不太合适^[12]。Izuka等基于垂直水头梯度提出一种运用测量界面上方非完整井中的钻井水头来估算界面深度的方法。虽然数据收集和含水层的非均匀性等方面会带来一些困难, 但这些困难都可以被最小化^[13]。Kim等利用单个钻孔中获得的淡水区和咸水区两组压力数据估算咸淡水界面的深度, 结果与实测的电导率剖面有良好的一致性^[14]。Elena应用水化学相演化图对西班牙东海岸的Vinaroz沿海平原海岸带咸淡水界面的动态进行了研究, 并结合其在GIS地图上空间变化的量化结果, 对含水层中海水入侵的状态进行了确定, 为识别沿海冲积含水层的海水入侵的状态和演变提供了一个简单有用的方法^[15]。魏玲娜等利用由Strack推导的地下水流动的饱和深度公式, 联立饱和达西定律, 结合对应势函数之间的关系式, 运用势函数法建立突变界面模型, 并对咸淡水界面的动态变化进行了研究。但因其忽略了咸淡水之间的水动力弥散作用, 太过于理想化, 且又是概化模型, 并不适用于过渡带较宽的地区^[16]。上述研究在考虑到流量不同、渗透率不同、咸淡水界面倾斜等其他影响因素的基础上, 进一步对咸淡水界面问题进行了探讨, 但所考虑的实际条件大大简化, 多属理论上的探讨, 算例过于理想化, 即使应用于实例计算, 所得的解析解也只是近似解, 与实际情况有所差别, 且具有一定局限性。

根据海岸带水头或压力来确定咸淡水界面的位置仍然是界面问题的最简单的方法, 近几年来又取得了一些新进展。周训等考察同一垂直线上咸淡水界面之上和之下任意点的水头对滨海均质各向同性潜水含水层淡水带地下水流动满足Dupuit假设时咸淡水界面的位置进行了确定, 并指出Ghyben-Herzberg公式是这一方法在咸水带水头与平均海平面相同情形下的特例, 描述界面位置的Hubbert公式也是该方法的一个特例^{[1, 4]}。其计算公式如下:

式中: z_i——咸淡水界面上任意一点的标高;

h_s和h_f——与该点在同一垂直线上的咸水带任意点咸水水头和淡水带任意点淡水水头。

周训等还提出了用压力代替水头的描述界面位置的公式^{[1, 4]}:

式中: z_s和z_f——同一垂直线上咸水带任意点和淡水带任意点的标高;

P_s和P_f——z_s和z_f点的压力。

利用该方法测算广西北海市西北海岸带咸淡水界面位置。该方法不要求海岸带处于绝对的静水状态, 使用条件较为宽泛。

2.   咸淡水界面的数值研究

与解析法相比, 采用数值模型研究海岸带的咸淡水界面具有可以刻画诸多复杂的水文地质条件和人为因素影响下的过渡带或咸淡水界面运移规律的优点。数值模型方法已成为模拟和求解海水入侵问题的有力工具^[17]。

Shamirh等假设咸淡水界面符合线性方程, 忽略淡水带地下水垂向流速, 建立了滨海含水层一维模型, 并给出了数值解^[18]。Segol等建立了剖面二维有限元模型, 模型考虑了流体压力、速度和浓度等因素, 并对美国佛罗里达州南部的某个海岸带垂直剖面上非稳定流条件下的海水入侵进行了计算。然而, 为了确保通过边界的流速连续分布, 模型必须同时解决流体的速度和压力, 但经验表明, 只有当含水层的渗透性逐级变化时, 才能使用渗透系数的函数表示^[19~20]。Wilson等提出采用固定网格的间接数值程序和特殊的尖端追踪算法追踪移动边界, 运用有限元方法对突变界面的运移进行了模拟^[21]。该模型在利用区域基本水平流动模型预测地下水含水层的海水入侵上, 比等效网格和移动网格体系更经济和准确。Diersch等通过FEFLOW建立了有关界面升锥问题的有限元模型, 结合水动力弥散作用的分析, 给出了界面的最大升锥高度, 并计算了地下淡水的流量及含盐量。模型结果与Dagan和Bear的预测结果基本符合^[22]。Essaid建立了通用的准三维有限差分数值模型, 用来探讨多层滨海含水层中咸淡水突变界面向内陆的移动问题, 并运用该模型研究了美国加利福尼亚州Soquel-Aptos盆地的海水入侵, 结果表明, 即使深层含水层的地下水位低于海平面, 浅层含水层的淡水仍流向海洋^[23]。Lin等运用SEAWAT模块对美国阿拉巴马州海湾未来可能发生的海水入侵范围进行了预测, 结果显示海水入侵会进一步加剧, 且海水入侵在深层含水层的侵入状况比浅层含水层更为严重^[24]。Choquet等引入一个考虑突变和扩散界面的相场模型, 将三维问题简化为二维模型, 描述了咸淡水突变界面和潜水面两个自由表面之间深度的演化, 并指出该模型在初始条件和边界条件下存在弱解, 同时还证明了两个界面的深度满足耦合极大值原理^[25]。

国内学者对咸淡水界面数值模拟也进行了有关的探索。吕贤弼应用二维边界数值模拟方法和狭缝槽模型, 对影响咸淡水界面的因素进行了研究, 总结了界面变化的主要特征, 认为突变界面向内陆方向的推进比界面向海洋方向的消退更为缓慢, 为利用和保护海岸带淡水资源提供了依据^[26]。薛禹群等建立了三维特征有限元模型, 对山东龙口-莱州滨海地区含水层中的海水入侵进行了研究, 并与实际观测资料拟合较好, 是国内比较早的研究咸淡水界面运移规律的三维数值模型^[27~29]。艾康洪等给出了剖面二维非均质流体的对流-弥散海水入侵数值模型, 利用上游加权有限元方法, 考虑到水头、浓度与密度相互作用, 进而分析广西漫尾岛咸淡水界面的运移规律, 模型中边界的处理方法较以往有所改进, 经与观测孔实测水位的拟合来识别模型参数并探讨了漫尾岛咸淡水界面的成因^[30]。姜效典等提出用B样条函数方法求解山东省莱州市滨海地区咸淡水突变界面的位置来确定海水的入侵范围, 结果表明该方法原理简明、计算可靠、易于实现, 且比采用有限单元方法更为有效^[31]。唐心强等考虑到地下水的渗流特征, 利用数理分析的方法, 建立了模拟滨海地区含水层中咸淡水界面运移规律的二维仿真数值模型, 与以往不同的是, 该模型对含水层中底部隔水层为不平坦的不透水层时的咸淡水界面随着潮汐波动而变化的规律进行了探讨。与实验结果比较, 数学模型具有能够客观反映海岸带含水层中因海水入侵而引起的咸淡水界面的变化规律, 且精度也比较高^[32~33]。

3.   咸淡水界面的实验研究

通过实验室模拟实验研究滨海含水层海岸带咸淡水界面可以直观地观察界面的形态和移动特点, 但这方面的文献报导相对较少, 并且实验内容和研究目的各有不同。Rumer等通过采用长150 cm、宽15.2 cm、高60 cm的实验装置, 探讨了三种不同含水介质中零淡水通量的情况下咸淡水界面向内陆的运移速率问题^[34]。Kuan等在假设海岸含水层均质的前提下, 采用实验室模拟实验的方法研究了潮汐波动下咸淡水界面的变化, 并认为应充分考虑固定流量和固定水头控制的混合陆地边界与超高的潮汐水位引起的较少淡水排泄的联合影响^[35]。Mehdizadeh等通过采用长117 cm、宽5.2 cm、高60 cm的实验装置, 模拟了存在垂直通量滨海多层含水层系统的海水入侵, 比较分析了弥散模型和砂槽实验, 指出当向上的淡水越流绕过任何上层海水且仅流入上覆含水层中的淡水时, SEAWAT模拟和砂箱实验的结果良好^[36]。

唐杰等考虑到咸淡水边界水位有突变、咸水边界浓度不一、介质水平分层等情况, 通过利用大型S形回旋加长型渗流弥散砂槽模型试验, 应用电导率法对咸淡水之间混合带的动态特征进行了分析观测, 发现了混合带在发展演化的过程中存在回咸现象, 对混合带及其动态特性进行了室内砂槽试验, 给海水入侵研究带来新的活力^[37]。张奇考虑了海水位的潮汐变化对入侵界面的影响, 在渗流槽内填充玻璃珠作为水流介质, 利用着色咸水结合图像处理方法对咸水向含水层中的入侵进行室内实验研究, 该方法优越于传统的水样化验方法, 不干扰流场、且数据在空间连续分布、精度较高等。利用实验对海洋的潮汐条件进行了模拟, 揭示了入侵界面在海水水位的周期性变化下的动态特征, 为数学模型的建立和验证提供了基础。并指出在潮汐作用下, 入侵界面随着海水水位的周期性波动作往返运动中, 其形态和宽度并没有发生明显的改变^[38]。刘晓红等采用长、宽、高分别为120 cm、2 cm、5 cm的砂箱装置, 以白色石英砂为含水层介质, 对海水入侵界面的形成和发展规律进行了探讨, 并对实验数据进行拟合, 求解出海水入侵速度的经验公式^[39]。唐心强等采用实验模拟的方法, 探讨了不同水面坡度条件下咸水的侵入和后退特点, 对不同的抽水位置和水面坡度下咸水楔的运移特征进行了考察, 并结合实验的结果给出了能预测因地下水位的升降而引起的咸水楔前后来回移动的数学模型。由于海岸带含水层地下水时刻受到周期性的潮汐作用影响, 即咸水楔始终处于不停的动态变化状态, 而该实验却是以假定咸水楔可以达到平衡状态为前提进行的。因此, 与实际情况相比该实验结果存在一定的差异^[40]。何丽等通过砂箱模型进行室内试验模拟研究了咸淡水界面的运移速率, 利用着色技术, 获得了咸、淡水之间界面的形态及其运移变化, 为数值模拟打下了基础^[41]。宋超采用实验室模拟实验的方法, 结合数值模拟的方法, 研究了滨海含水层中海水入侵的过程, 固定其他条件不变, 分析依次改变抽水量、咸水密度、水头差条件下的咸淡水界面的运移规律, 对合理开发海岸带地下水具有重要意义。但是该研究忽略了咸水边界处淡水排泄会导致咸水密度降低的问题, 观测结果可信度有待提高。其数值模拟只考虑到二维剖面情况, 在后续研究中可以考虑建立三维海水入侵数值模型^[42]。

4.   咸淡水界面的其他方法研究

除了理论研究和实验研究以外, 运用地球物理探测也可以对咸淡水界面的变化进行分析, 在观测井中测量地下水电导率是常用的方法。张彦增等利用地面电测法对河北省阜城县的深层地下水中的咸淡水界面进行观测, 并结合1975-1994年间深井机井位置、深度、含水层组、咸淡水界面等资料, 绘制出20世纪80年代和90年代两个不同时期的深层咸淡水界面等值线图, 从水文地质条件、成井结构、深井开采强度等影响因素, 分析了咸淡水界面的下移规律, 提出抑制界面下移的有关措施, 但对咸水入侵的机理、危害及更有效防治措施的研究还有待加强^[43]。杨化勇等对山东潍坊市的滨海平原咸水入侵进行了调查, 根据监测井点水样中氯离子含量化验结果, 绘制出咸淡水分界线, 并分析了咸水入侵的成因, 提出了一些防治措施, 为滨海地区海岸带咸水入侵调查研究与防治提供了参考^[44]。唐心强等采用Wenner法对太平洋沿岸某地区滨海含水层的咸淡水界面进行了探查, 实地调查验证了电法探查咸淡水界面的可行性、适用性和可能性, 根据海岸纵深方向的调查井垂直方向的电导率确定了不同类型土质构成的岩土层咸淡水界面的深度, 是一种成本低廉、可靠、简单、快速的咸淡水界面深度的确定方法, 具有较高的灵敏度, 且探查到海岸线附近含水层咸淡水界面因潮汐作用而引起的微小变动, 是监测海水入侵的一种有效手段^[45]。

5.   结论与展望

滨海含水层海岸带咸淡水界面的研究历经百余年的努力, 已取得了许多成就。早期研究工作主要侧重于理论研究, 特别是解析研究。近数十年来, 随着数值方面模拟技术的发展, 数值研究已逐渐取代解析研究, 成为海水入侵问题研究的最有力工具。与此同时, 海岸带过渡带模型的研究也取得了重要的进展。今后可以关注和加强如下几个方面的研究:

(1) 边界条件的确定与概化是否合理, 明显影响过渡带的形状或淡水透镜体的形状和厚度, 已有数值模型的研究多侧重于计算格式和数值方法的讨论, 对边界条件的论证与处理不够细微^[2]。如何准确地确定边界处的水头、浓度或通量是不可忽视的问题。

(2) 通常水动力参数和水质参数是含水层中咸淡水界面模型计算模拟的基础, 其大小会影响过渡带或咸淡水界面的宽度。在滨海地区, 周期性的潮汐波动变化会引起观测孔中的水位周期性变化, 如何能够更加合理地直接利用地下水位潮汐动态求取含水层参数是值得深入研究的方向。

(3) 海水入侵的研究均是在海水入侵危害发生比较严重的地区开展, 然而一旦发生海水入侵, 治理起来就会变得十分困难, 因而应当进行超前的研究和预报, 以防止海水入侵。这就需要分析在不同的天然条件和人类活动影响下咸淡水界面的运移和变化规律, 从而达到减轻和防止海水入侵的目的。

(4) 一般在海岸带布设井孔, 建立地下水三维动态监测网对海水入侵的程度和范围进行观测判别, 耗时多、投入大, 积极探求有效且简便的能显示和追踪过渡带或界面运移变化的方法并进行相关仪器的研制, 是有待努力的方向。

(5) 实验室模拟实验方面, 已有咸淡水界面向内陆延伸距离的计算公式未考虑咸水边界处淡水排泄, 对咸水边界处淡水排泄导致咸水密度降低等问题研究不够, 还应多做研究和改进。

(6) 咸淡水界面的解析研究多以满足Dupuit假设为前提而进行, 如何确定不满足这一假设条件下的咸淡水界面位置值得深入研究。