在自然环境下，滨海含水层咸淡水处于一种动态平衡状态。然而，海平面的上升以及许多沿海地区地下水的过渡开采，导致旧的平衡被打破，海水向内陆含水层迁移，造成海水入侵，含水层水质恶化^[1-2]。深入研究海水入侵过程的盐度分布变化规律，对于科学管理、利用滨海地区地下水资源、防范滨海地区海水入侵危害有着重要意义。

海水入侵最关键的问题是研究咸淡水界面形态、位置及发展规律。早期的研究通常忽略潮汐等引起的海水面波动的影响, 并假设海相边界水位为平均海平面的固定边界，如Ghyben和Herzberg早在100多年前就推导了用以确定咸淡水突变界面位置的Ghyben-Herzberg公式^[3]。Henry^[4]根据咸淡水混溶原理，将变密度流与盐度输移方程耦合求解，提出了垂直海岸线剖面上盐度分布的解析解。

而对于滨海含水层，地下水流模式及盐度分布规律通常受到潮汐作用的影响^[5]。潮汐作用会引起含水层中地下水位的波动，这种现象已经得到了较为广泛的关注^[6-8]。地下水位的波动必然会对盐度输移过程造成影响，目前越来越多的国内外研究者已经开始研究并考虑海水入侵过程中的潮汐作用，Robinson等^[9]采用有限元法模拟了地下河口地下水流动和盐分运移过程，并将模拟结果与实验数据进行对比分析。Boufadel等^[10]通过实验和数值模拟两种方法研究了稳态与瞬态下溶质的运移机理。Li等^[11]考虑了浸润面及毛细压力等多种因素，研究了地下水与海水的循环模式。Vandenbohede等^[12]通过野外观测和地下水流模拟，研究了斜坡海岸潮汐作用下地下水流模式以及盐度运移规律。

潮汐波动动态地改变了海相边界条件，尤其当海滩为小坡度倾斜海岸时，在潮汐落潮阶段海滩上会形成渗出面，非线性变化的边界使海相边界条件变得更为复杂。近年来，国内外学者对于滨海潜水层海水入侵的数值模拟过程中常用的一些数值模拟软件，已经不足以准确模拟上述情形。在这方面，前人已作了一些探索，如荷兰学者Post^[13]在模拟滨海含水层海水入侵问题时, 编写了周期性边界条件程序包, 来模拟潮汐周期性波动边界的影响, 并将其嵌入MODFLOW与SEAWAT中，模拟结果与前人数据吻合较好。

本文基于OpenGeoSys科学模拟软件，研究潮汐对潜水含水层海水入侵机制的影响，海相边界分别设定为固定边界和移动边界两种不同情形，建立斜坡海滩潜水层海水入侵的二维剖面数值模型。模拟结果表明，移动边界能够很好的体现出海相边界随潮汐变化的动态变化特征。在此基础上探讨了潜水含水层受潮汐影响下的水位波动及盐度运移规律。

1 移动边界处理方法 1.1 移动边界特征

以往的海水入侵研究常将海相边界设置为固定边界条件，即设为定水头和定浓度边界。然而海相边界常受到潮汐波动的影响，潮汐引起的海水位变化会对海水入侵过程造成影响，将海相边界设置为固定边界显然是不合理的。为考虑潮汐波动作用，后来研究者们常将海相边界概化为“L”型的垂直边界，即假设海岸线水平(水陆交界面)位置保持不变，潮汐波动仅引起海水位在边界位置上垂向变化，这种简化对于倾角较大的岩质海岸是合理的。

然而对于倾角较小的砂质海岸(见图 1)，由于海水水位的波动不仅会造成水位垂向的变化，而且会造成海岸线位置沿水平方向移动，因此将海相边界概化为“L”型边界也不太合理。这种情况下，海相边界会在垂直和水平方向都存在与潮汐相关的周期性变化，我们将这种边界称为移动边界。移动边界条件下的海水入侵问题求解更为复杂，解析法很难求出其解析解，通常采用数值模拟的方法。

1.2 移动边界判定

潮汐引起的海水水位变化可以分解为多个规则的简谐振动，每个简谐振动称为一个分潮，实际运用中通常选取11个较大的分潮，就可以得到较为准确的海面水位变化规律，即

$ {H_{\rm{s}}}(t) = {H_0} + {A_i}\sin \left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{{T_i}}}t - {\varphi _i}} \right)。$

(1)

式中：H_s(t)表示海水水位；H₀表示平均海平面位置；A_i、T_i与φ_i分别表示各个分潮的振幅、周期和迟角。

现假设海平面水位为0，存在一个振幅为A、周期为T、迟角为0的潮汐波。

在$0 < t < \frac{1}{4}T$时，处于涨潮阶段，海水水位上升，海岸线向内陆方向移动。在此过程的各个时刻中处于海岸线下方(被海水淹没)的海相边界为定水头和定浓度边界，水头等于该时刻海水面高程，盐度等于海水盐度；处于海岸线上方(未被海水淹没)的海相边界为零流量边界。

在$\frac{1}{4}T < t < \frac{3}{4}T$时，处于落潮阶段，海水水位下降，海岸线向海侧移动。在此过程中各个时刻处于海岸线下方的海相边界为定水头和定浓度边界，水头等于该时刻海水面高程，盐度等于海水盐度；由于在落潮时地下水位回落速度比海水位的变化速度慢一些，因此会在斜坡海岸上形成地下水流渗出面，渗出面上水头为该位置的高程值，盐度等于流出的地下水中盐度值；在渗出面上方的边界为零流量边界。

在$\frac{3}{4}T < t < T$时，再次处于涨潮阶段。在此过程中，若海岸上存在渗出面，则边界条件与$0 < t < \frac{1}{4}T$时的边界条件设置相同，若由于海水位的上升导致渗出面完全被海水淹没，则边界条件与$\frac{1}{4}T < t < \frac{3}{4}T$时的边界条件设置相同。

通过进一步分析发现，在一个潮汐波动周期内，海岸线的水平位置变化幅度为$\frac{{2A}}{{\tan \alpha }}$。即海岸坡度越小，海岸线在水平方向变动范围越大，这与江峰的研究结果一致^[14]。这种情况下，采用移动边界可以更好的刻画海相边界，减小数值模拟过程中产生的误差。

1.3 移动边界条件模拟方法

本文的研究工作基于德国亥姆霍兹环境研究中心开发的OpenGeoSys科学模拟软件，它是一个由C++实现、面向对象的建模源程序工具包，数值计算采用有限元方法，能够模拟多孔介质或裂隙介质中单个或者耦合的传热-流动-力学-化学过程。

在模拟潮汐作用下的海水入侵过程时，采用以下方法来设定海相移动边界条件：

(1) 每个时间步开始计算时，通过方程(1)动态地计算每个时间步的海水水位值H_s。

(2) 将海相边界上所有离散点的初始水头都设置为H_s，盐度设为海水盐度C_s。

(3) 运行OpenGeoSys科学模拟软件，计算得到所有海相边界离散点的压强p与浓度C。

(4) 对于海相边界上所有的离散点，若该点高程z < H_s，则设为定水头和定浓度边界。

(5) 若该点高程z > H_s且压强p > 0，则该点处在渗出面上，即令该点压强值p=0。

(6) 将与条件(2)、(3)都不满足的海相边界上其它点设置为零流量边界。

(7) 重复步骤(3)~(6)，直至计算收敛，进入下一时间步。

(8) 重复步骤(2)~(7)，直至计算完成。

图 2给出了在OpenGeoSys中实现动边界处理的每个时间步长的数值模拟流程图。

2 海水入侵数值模型 2.1 模型构建

本文在Kuan^[15]等物理模型试验的基础上，建立如图 3所示的实验室尺度下的滨海潜水含水层海水入侵的数值模型。其中，潜水含水层长为3.4 m，厚度为0.7 m，假设潜水含水层介质各向同性，底板水平隔水，不考虑降雨入渗等影响。左侧为内陆淡水边界，有定流量的淡水补给。右侧为坡度为1:3.8的斜坡海岸，海岸水平长度为1.482 m，海岸最右端端点处高度为0.27 m。在不考虑潮汐作用时，右侧海水位维持在平均海平面位置保持不变，在考虑潮汐作用时，右侧海水位会随潮汐波动而上下波动，则此时海相边界应设置为移动边界。

在OpenGeoSys软件中，建立变密度流情况下的渗流-弥散模型，即水流控制方程选取用来描述非饱和流动的Richards Flow方程，盐度扩散过程选取溶质输移方程，不考虑含水层中的挥发、吸附、生物及化学反应等作用。

水力特征函数采用van Genuchten水力特征函数，即

$ S = {S_{\rm{r}}} + \left( {1 - {S_{\rm{r}}}} \right){\left\{ {\frac{1}{{\left[ {1 + {{\left( {a{p_{\rm{c}}}} \right)}^m}} \right]}}} \right\}^{1 - \frac{1}{m}}}。$

(2)

其中：S为饱和度；S_r表示孔隙水残余饱和度；p_c表示毛细压强，与饱和度S度相关；a表示进气压强水头比例因子，反应岩土空隙的大小；m为指数参数，表示孔隙尺寸分布参数。

变密度流密度ρ与浓度的关系按下式计算，即

$ \rho (c) = {\rho _0}\left[ {1 + {\beta _{\rm{c}}}\left( {c - {c_0}} \right)} \right]。$

(3)

其中：c表示溶质质量浓度；ρ₀表示淡水密度；c₀表示淡水浓度；β_c表示溶质膨胀系数。

本文数值模型采用有限元法进行模拟研究，将整个模拟区域通过三角形网格进行划分，共有5 670个节点，包含三角形网格10 944个。为了更加精确体现海水侵入潜水含水层的整个过程，对海岸线潮间带及潮间带以下区域的网格进行局部加密，网格密度增大为其他区域网格密度的一倍。数值模型中采用的相关参数取值见表 1。

2.2 模型验证

依据物理模型试验分别模拟有、无潮汐两种条件下海水入侵过程。平均海平面为0.398 m，潮汐周期T为62 s，振幅A为3.2 cm，海相边界设为移动边界，内陆淡水边界设为定流量边界，流量为20 mL/min。并将数值模拟结果与试验结果进行比较，以验证数值模拟结果的有效性。为更清晰刻画盐度输移结果，选取0.5盐度等值线作为咸淡水的分界线，模拟结果如图 4所示。

模拟结果表明，不考虑潮汐作用时，海水从海相边界处开始入侵潜水含水层，在含水层中形成单一的海水入侵盐水楔。而考虑潮汐作用时，由于潮汐波动引起的潮间带区域海水-地下水交换，因此会在潮间带区域内的潜水含水层中形成半椭圆形状的高盐度区，称为上盐水楔。上盐水楔的形成改变了含水层中盐度的分布规律，从图 4中可以看出，潮汐作用减小了盐水楔向内陆延伸的距离，这在一定程度上减缓了海水入侵的程度。数值模拟结果与Kuan等^[15]的实验结果吻合较好，表明所建模型以及移动边界的处理方法能够准确模拟潮汐波动动态边界条件下的滨海潜水含水层海水入侵过程。

3 潮汐对海水入侵的影响分析

为进一步探讨潮汐波动对潜水含水层海水入侵规律的影响，引入不同的淡水边界条件和潮汐条件，分析潜水含水层淡水边界和潮汐条件对海水入侵过程的影响。不同模型条件参数取值见表 2。

3.1 潮汐对海水入侵过程的影响

潮汐作用使得咸淡水交换过程更为频繁，在涨潮时，海水由高潮线附近的海滩渗入潜水含水层，落潮时，由低潮线附近的海滩出流。与无潮条件相比，这种海水-地下水循环模式以及地下水位的波动必然会对海水入侵过程产生影响，从而影响潜水含水层中的海水入侵过程。

3.1.1 潜水含水层中的潮汐效应

选取Q1A2为研究对象，在模型中距离淡水边界1.6、1.9和2.2 m位置处选取三点，来获取地下水位的变化情况，不同位置处水头变化幅度随时间的关系曲线如图 5所示。

模拟结果表明，周期性的潮汐波动会引起潜水含水层中水位的周期性波动，且潜水含水层中水位波动周期与潮汐周期一致。位置离海岸线最近，则水位波动振幅最大，水位波动振幅随着与海岸线距离的增大而依次减小，表明潮汐的影响在向内陆的传播过程中不断衰减，越靠近内陆，水位受潮汐的影响越小，水位波动最终会衰减为零。同时，地下水波动的相位与潮汐波动相比存在滞后性，滞后时间随着离岸的水平距离的增加而不断增大。

另外，地下水波动的平均水位要高于初始水位值，即潮汐作用会引起地下水位超高。Li等^[7]总结得到水位超高是由于沙滩倾角导致的边界条件的不对称、浸润面的形成以及控制方程的非线性共同作用下的结果。

3.1.2 地下淡水排放通道变化

地下淡水排放通道为地下淡水通过含水层边界流入海洋的区域。选取Q1A2、Q1A0为研究对象，模拟潮汐波动对地下淡水排放通道的影响，模拟结果如图 6所示。

在不考虑潮汐作用时，海平面和海岸的交点与盐水楔和海岸的交点两点之间的海岸为潜水含水层中地下淡水排放通道。潮汐作用不仅会改变含水层中盐度的分布，同时也会改变含水层中地下水淡水的排放通道。

从图 6中可以看出，在考虑潮汐作用时，由于上盐水楔的存在，上盐水楔与下盐水楔之间区域成为淡水排放的通道，与无潮条件下相比，潮汐作用使下盐水楔和海岸的交点向低潮线方向移动，淡水排放通道也同时向低潮线方向移动。

3.2 潮汐振幅对海水入侵的影响

分别选取潮汐振幅不同的三种工况Q1A1、Q1A2、Q1A3为研究对象，探讨潮汐振幅对海水入侵的影响，不同潮汐振幅条件下海水入侵数值模拟结果如图 7所示。可以看出，潮汐振幅大小对潜水含水层海水入侵过程有着显著的影响。潮汐振幅越小，上盐水楔的分布范围越小，潜水含水层中地下水淡水排放通道向低潮线偏移距离越短。这是因为潮汐是引起潮间带海岸区域水体循环的主要动力因素，潮汐振幅的减小减弱了潮间带水体循环的动力因素，从而使得入侵淡水形成的上盐水楔范围随之减小。同时，由于上盐水楔的形成会限制下盐水楔的入侵距离，上盐水楔的减小会令下盐水楔向内陆移动，入侵潜水含水层的距离增大。

因此，陆地边界流量相同时，潮汐振幅的减小使得上盐水楔收缩，下盐水楔向陆地扩张。

3.3 陆源流量对海水入侵的影响

从图 8可以看出，在相同潮汐边界条件下，内陆淡水边界的流量越小，则潜水含水层中海水入侵淡水的距离越远。同时，内陆边界流量为20 mL/min时，潮汐对盐水楔成长的影响明显小于内陆边界流量为15 mL/min时，这表明在潮汐振幅相同时，陆地边界流量越小受潮汐影响越大。这是因为当淡水流量减小时，在潮汐作用下形成的上盐水楔会进一步在潜水含水层中扩散，而分布范围更大的上盐水楔会再次限制下盐水楔的入侵距离。

因此，潮汐条件相同时，内陆淡水边界流量越小，海水入侵距离越远，并且受潮汐的影响越大。

3.4 多个分潮对海水入侵的影响

潮汐引起的海水位波动是多个分潮共同作用下的结果，为进一步研究多个分潮的影响，假设向海一侧同时存在两个分潮，探讨两个分潮共同作用下的海水入侵作用，两个分潮其中一个潮周期为62 s，振幅为4.2 cm，另一个潮周期为124 s，振幅为2.1 cm，参数取值如表 2所示。

图 9给出了分别在Q1A2、Q1A4以及两个分潮共同作用下的潜水含水层海水入侵模拟结果。当海侧潮汐振幅为0.021 m时，由于潮汐动力因素太弱，因此在模拟结果中并未在潜水含水层中形成上盐水楔，此时的盐度分布基本上与稳态时保持一致。但在两个分潮共同作用下，上盐水楔的扩散范围比只考虑单一分潮的二者都要大，这是由于多个分潮的存在增强了海洋中的水动力作用，促进了潮间带海水向潜水含水层的扩散，使得上盐水楔进一步成长，下盐水楔入侵距离减小。

4 结论

本文建立了潮汐影响下复杂海相边界海水入侵的数值模型，实现了OpenGeoSys模拟海水入侵时移动边界的动态处理，得到了与以固定边界为海相边界的传统海水入侵模型不同的盐度分布规律，重点探讨了潮汐波动对潜水含水层海水入侵规律的影响，达到了深入研究海水入侵机理的目的。本文得到的主要结论有以下几个方面：

(1) 潮汐波动会引起潜水含水层地下水水位超高；潮汐引起的水位波幅随着离岸的水平距离的增加而不断衰减，相位与潮汐波相比存在滞后，滞后时间随着离岸的水平距离的增加而不断增大。

(2) 斜坡海滩上的潮汐波动会引起潮间带区域相对较快的水体循环，这种海水-淡水循环促使了上盐水楔的形成，上盐水楔的存在改变了海水入侵过程中盐度的分布及淡水排放通道，并且潮汐的存在会在一定程度上减缓海水入侵程度。

(3) 陆地边界流量相同时，潮汐振幅的减小使得上盐水楔收缩，下盐水楔向陆地扩张。

(4) 潮汐条件相同时，内陆淡水边界流量越小，海水入侵距离越远，并且受潮汐的影响越大。

(5) 多个分潮的存在使得上盐水楔加剧扩散，下盐水楔向海收缩。

本文只进行了实验室尺度下物理模型的数值模拟分析，对于实际尺度下的应用，需要考虑更为复杂的降雨入渗、地下水开采、各向异性介质等问题，还有待进一步研究。