海陆过渡带氯离子迁移规律分析研究

贾超，朱恒华，虞未江，卫政润，刘治政，李双

(1.山东大学海洋研究院，山东 青岛 266237;2.山东省地质调查院，山东 济南 250014)

0 引言

随着经济社会的发展，人类生活和工业生产需水量激增，许多城市都开采地下水来满足用水量，尤其沿海城市因为过度开采地下水，导致内陆地下水位大幅度下降，破坏了海咸水与内陆地下淡水之间的动态平衡，从而引发海岸带咸淡水界面向内陆迁移，即海水入侵[1-2]。美国、荷兰、西班牙等国家的滨海城市均出现海水入侵问题，中国辽东湾、山东半岛、上海市等地区也发生海水入侵，其中山东莱州湾地区是海水入侵较为严重的地区之一，海咸水入侵将污染内陆地下淡水资源，破坏地下水环境，造成沿海地区土地固结化和盐渍化，严重阻碍当地经济社会发展。

海水入侵现象在全球范围内普遍存在[3-9]，国内外众多学者对海水入侵问题进行深入研究，在海水入侵模型、基本理论、咸淡水界面形状和迁移机理等方面取得许多重要成果。目前，主要有两类模型对海水入侵进行刻画，即咸淡水界面突变模型和海陆过渡带模型。海陆过渡带模型因其在刻画复杂水文地质条件和人类活动条件等方面具有优势，模型的仿真性和可靠性较高，而被国内外普遍采用[10-15]。

国内对海水入侵已有较多研究，该文在相关研究的基础上[16-18]，选择海陆过渡带模型，结合地下水运动理论，选取受海水入侵危害较严重的山东省龙口、招远滨海区域作为研究区，研究该地区海陆过渡带迁移的规律。在对该区域水文地质基础资料进行研究分析的基础上，构建研究区水文地质概念模型，建立地下水流场数值模型，并对数值模型进行识别和检验。选择作为模拟因子，在地下水流场的基础上，对研究区海陆过渡带的迁移进行数值模拟并研究内陆地下水位对海陆过渡带迁移的影响，总结分析研究区海陆过渡带迁移的规律，并为该地区合理开发利用沿海地下淡水资源和防控海咸水入侵提供切实有效的建议。

1 理论方法

1.1 海陆过渡带模型

沿海地区地下含水层中的水力梯度基本上是倾向海洋，内陆地下含水层中的大部分淡水最终流向海洋。因为海平面以下的地下含水层中存在有海咸水，所以流向海洋的淡水与海咸水之间存在一条接触带。此外，海咸水与淡水是两种可以互溶的流体，在流体动力弥散作用下，该接触带将以过渡带形式存在。海陆过渡带模型及咸淡水流动如图1所示。

图1 地下含水层中海咸水和淡水流动过程示意图

海陆过渡带模型由两个数学偏微分方程描述，一个是描述密度不断改变的混合液体流动方程，另一个是描述海水中盐分运移的对流—弥散方程，通过两个数学方程，把流体的密度、浓度和水位联系在一起，得到海陆过渡带的分布范围、浓度值和水位值。

(1)描述密度不断改变的混合液体流动方程

(1)

该式是在假设密度随浓度呈线性变化和忽略液体压强对密度影响的基础上推导出来[12]。

(2)描述海水中盐分运移的对流—弥散方程

(2)

式中Dij—弥散系数张量，(i,j)=(1,2,3)；u=(u1,u2,u3)—地下水实际平均流速矢量；c*—注入或者抽出液体中相应组分的浓度。

1.2 地下水运动理论

(1)地下水运动微分方程：在达西定律和质量守恒定律基础上，推导出非均质各向异性含水层中承压水运动的基本微分方程：

(3)

(2)溶质运移模型的控制方程：描述模拟因子k在三维地下水流系统中的迁移和去向：

(4)

2 研究区海陆过渡带迁移规律研究

2.1 研究区概况

研究区龙口、招远滨海区域位于山东省胶东半岛西北部，渤海湾南岸，西与莱州湾相接，东北与天津、大连等城市以及朝鲜半岛隔海相望，地理坐标为东经120°13′38″～120°34′12″，北纬37°45′16″～37°29′40″之间，研究区东西最大横距约为23km，南北最大纵距约为25km，海岸线长约为55km，研究区总面积约为424km2，地理位置如图2所示。

图2 研究区地理位置示意图

研究区属暖温带半湿润季风型大陆性气候，年平均气温为11.6℃，多年平均降水量为592.7mm，主要集中在夏季的7，8月份，约占全年降水量的50%，该区域多年平均地表蒸发量为1479mm。研究区地势由东南向西北逐渐降低，沿海地势平缓，海拔高程多在0～15m之间，属于滨海平原区。该区属于龙口-威海侵入岩变质岩冲洪积海积工程地质亚区，含水层呈多元结构，一般发育有3～4层，且其表层主要以砂性土为主，砂性土以下依次以粉土和粉质黏土为主，区域内地下水从东南流向西北的海域，主要是第四系的松散岩类孔隙水，水化学类型为氯化物型，水交替强烈。

该区地下水的赋存条件差，富水性弱。仅在局部大理岩分布区及断裂构造发育带，为地下水的赋存创造小面积的或点状、线状的富集环境，形成小区域大理岩岩溶水富水区及变质岩、侵入岩裂隙水点或线状较富水带。区内河流发育短小，第四系松散层多呈条带状分布于现代河床及两侧附近，厚度小，富水性弱，部分河流谷地较开阔，形成小型河谷平原及沿岸山前倾斜平原孔隙水较富水地段。区内基岩裂隙水广泛分布，可分为层状岩类裂隙水、块状岩类裂隙水及喷出岩孔洞裂隙水，大部分属风化裂隙水。一般富水性较弱，个别地段孔洞裂隙水水量较大，如蓬莱市南王镇一带，形成小的富水地段。山区由于侵蚀切割作用，形成水系较多，基岩裂隙水主要靠大气降水补给。碳酸盐岩类岩溶裂隙水主要分布在胶东半岛的大理岩分布区，含水层主要为张格庄组大理岩，另外在蓬莱南王镇一带隐伏的张格庄组大理岩和牟平养马岛野头组大理岩也形成了富水地段。区内主要分布寿光水源地、潍寒水源地、潍黄旗堡水源地、朱里水源地、穆村水源地、昌邑水源地等龙口、招远地区城市供水近70%来自地下水，第四系松散层中的孔隙潜水是主要的开采对象，地下水水资源的超量开采使该地区遭受海水入侵，地下水被污染。

该区海域洋流水动力特征季节性显著，洋流分布特性如图3、图4所示。

图3 莱州湾及附近海域冬季洋流分布示意图

图4 莱州湾及附近海域夏季洋流分布示意图

2.2 地下水流场计算

对研究区域地形地貌、地质构造以及地层结构进行调查和分析，认为该区域地质结构完整，各含水层结构基本连续，可利用GMS软件MODFLOW模块对研究区地下水流场进行数值模拟。

2.2.1 水文地质概念模型

建立该区域的水文地质概念模型，主要包括含水层系统、边界条件、源汇项等3个方面。

(1)含水层系统：龙口、招远滨海研究区内地下各含水层的分布较连续，垂向上呈明显的层状分布，根据研究区地质资料，可将研究区地下水含水层概化为3层，由地表向下，含水层岩性依次以砂性土、粉土和粉质黏土为主。

(2)边界条件：垂向上，研究区含水层上部主要接受大气降水补给，为水量交换边界，下部与相对隔水的黏土层、基岩接触，为隔水边界；水平向上，模拟区北侧、西侧和西北侧与渤海相接，东北角和西南角分别以黄水河、界河为边界，都概化为一类定水头边界，南侧边界为龙口内陆低山丘陵区，概化为二类定流量边界，研究区边界条件及地下水位观测点如图5所示。

图5 研究区边界条件及地下水位观测点图

(3)源汇项：研究区内地下水的补给以降水入渗和山丘区基岩裂隙水的侧向径流为主；研究区濒临渤海湾，属滨海平原地带，区域内的地下水位埋深较浅，普遍小于5m，需要考虑地下水的蒸发量。同时，孔隙潜水的开采也是该区域地下水排泄的一个重要途径，根据研究区的统计资料，确定开采井的位置、开采深度和具体的开采量。

2.2.2 地下水流场数学模型建立及求解

在对区域水文地质概念模型分析研究的基础上，对地下水水流进行概化，建立非均质、各向同性的二维非稳定地下水流数学模型。

H(x,y,t)|t=0=H0(x,y,t) (x,y)∈D

H(x,y,t)|Γ1=H1(x,y,t) (x,y)∈Γ1,t>0

(5)

式中：K—潜水含水层渗透系数(m/d)；μ—潜水含水层给水度；Z—含水层底板标高(m)；H—潜水水位(m)；H1—一类边界点的水位(m)；q—二类边界单宽流量(m3/d/m)；D—计算区范围；Γ1,Γ2—一、二类边界。

利用GMS软件中的MODFLOW模块计算研究区地下水流场，模拟区域面积约为424km2，在GMS软件3DGRID模块下对研究进行网格剖分，共剖分3384个有效矩形网格单元(图6)。

图6 研究区网格剖分立体示意图

根据区域水文地质条件，研究区内各含水层水文地质参数(给水度、渗透系数、有效孔隙度)初始值如表1所示。

2.2.3 识别和检验地下水流场数值模型

地下水流场数值模型的识别和检验是相互独立的2个阶段，需要根据不同时间段的资料分别完成。

表1 研究区各含水层水文地质参数

该文以2009年7月—2010年6月作为模型识别期，反求水文地质参数，再次对研究区的地下水流场进行数值模拟；以2010年7月—2011年6月作为地下水流场数值模型的检验期，主要检验模拟研究区的地下水水位动态变化是否与实测的地下水位动态变化一致，两个地下水位观测点在检验期内实测地下水位值与计算值的对比如图7、图8所示。

图7 2010年7月--2011年6月海岱镇河南孙家水位观测点地下水位观测值与计算值对比图

图8 2010年7月--2011年6月徐福镇洼西村水位观测点地下水位观测值与计算值对比图

从图7和图8可知，海岱镇河南孙家和徐福镇洼西村水位监测点地下水位观测值与计算值之间的最大差值为0.48m，拟合程度较好，满足《地下水流技术数值模拟技术要求》GW1-D1中“水位变化值较小(小于5m)的情况下，水位拟合误差一般应小于0.5m”的要求。

此外，对地下水流场、源汇项均衡量和水文地质参数分别进行检验，各项与实际情况的拟合程度都较好，认为建立的地下水流场数值模型计算结果达到精度要求，反映了该区地下水系统的动力特征，可以用于后续地下水溶质运移模拟计算。

2.3 海陆过渡带迁移规律研究

2.3.1 地下水超采状态下海陆过渡带迁移

在研究区地下水流场计算的基础上，对该地区海水中Cl-在内陆地下水中的运移进行数值模拟，因为Cl-是海咸水中最稳定和最主要常量组成元素，可以直接反映出海咸水运移的状况。在龙口、招远滨海区域，可以将地下水中Cl-浓度超过200mg/L的地区视为海水入侵区域[11]。由于研究区地下水处于超采状态，故研究海陆过渡带的迁移是以地下水超采作为条件。

对研究区西北部约12km长海岸线上海咸水中的Cl-向内陆地下水中的运移过程进行数值模拟，并给作用于海岸线上的浓度赋初值为3g/L。模拟运移的计算期为7200d，每300d作为一个计算应力期，共分为24个计算应力期。

在内陆地下水中Cl-运移1000d，3000d，6000d后在第一层、第二层含水层中水平向的浓度分布如图9和图10所示。

图9 Cl-在内陆地下水中运移1000d，3000d，6000d后在第一层含水层中水平向浓度分布图

图10 Cl-在内陆地下水中运移1000d，3000d，6000d后在第二层含水层中水平向浓度分布图

研究海咸水中Cl-浓度在含水层中垂向上的变化，在模型中选取4个剖面，提出各剖面Cl-上运移6000d后的垂向浓度分布图。选取的剖面位置如图11所示。

模型I-I剖面、J-J剖面、K-K剖面、L-L剖面上Cl-运移6000d后，在含水层中垂向上浓度变化如图12所示。

由图9、图10、图11，Cl-在含水层中运移若干天后的水平向和垂向浓度分布如下：

(1)当内陆地下水超采，海咸水入侵，海咸水中的Cl-在内陆地下水中的运移主要发生在第一层砂性土层中，第二层粉土层受Cl-的影响较小，第三层粉质黏土层受Cl-的影响，可以忽略不计，即海咸水入侵污染内陆地下水的最大深度为两层含水层的厚度。

图11 Cl-浓度在垂向变化剖面位置图

a—I-I剖面垂向Cl-浓度变化图；b—J-J剖面垂向Cl-浓度变化图；c—K-K剖面垂向Cl-浓度变化图；d—L-L剖面Cl-垂向浓度变化图图12 模型I-I剖面、J-J剖面、K-K剖面、L-L剖面上Cl-运移6000d后，在含水层中垂向上浓度变化图

(2)内陆地下水水位下降，海咸水中的Cl-向沿海地区内陆地下水运移入侵，主要发生在前1000d之内，这期间内地下水中Cl-的浓度变化较大，3000d以后，地下水中Cl-运移趋于平缓，该Cl-运移特性可以作为防控海咸水入侵的依据。

2.3.2 Cl-浓度随时间变化

选取距海岸线不同距离的A，B，C，D，E，F等6个点，并从模型中提取出6个点处Cl-浓度变化曲线，研究该点处Cl-浓度随时间变化的特性。A，B，C，D，E，F距离海岸线的距离分别约为3km，4.5km，8.4km，12.5km，14.8km，16km，各点位置如13所示。

图13 距海岸线不同距离的A，B，C，D，E，F等6个点位置示意图

A，B，C，D，E，F各点处Cl-浓度随时间变化如图14所示。

距海岸线4.5km的B点处浓度运移模拟在500d时达到200mg/L，故综合考虑研究区在第一层含水层中运移3000d后的浓度分布图和以上六点处浓度随时间变化曲线图可知，沿海地区存在浓度恒大于200mg/L的地区，该类区域已被海咸水入侵污染，宽度约为1.5～4.5km，即研究区龙口、招远滨海区域海陆过渡带的宽度约为1.5～4.5km。南京大学薛禹群院士团队曾在20世纪通过大量的野外观察，对该地区的海水入侵状况进行了研究，认为该地区海陆过渡带的宽度约为1.5～2.5km。该文通过对氯离子运移进行数值模拟，对研究区划分的海陆过渡带为1.5～4.5km，比薛禹群院士团队研究的海陆过渡带宽度稍宽，说明经过20多年，该地区的海咸水继续向内陆运移，地下水继续被超采，海水入侵未得到有效遏制。

2.3.3 内陆地下水位对海陆过渡带迁移的影响

在目前地下水流场的基础上，先对研究区的地下水位进行改变，将地下水位整体降低1.5m和抬高1.5m，再对研究区海咸水中Cl-运移进行数值模拟，分析内陆地下水位对研究区海陆过渡带迁移的影响。地下水位抬高1.5m和降低1.5m后，在第一层含水层中运移1000d，3000d，6000d后的浓度分布如图15所示。

图14 A，B，C，D，E，F各点Cl-浓度随时间变化曲线图

a—地下水位抬高1.5m后，Cl-在第一层含水层中运移1000d，3000d，6000d后的浓度分布图；b—地下水位降低1.5m后，Cl-在第一层含水层中运移1000d，3000d，6000d后的浓度分布图图15 地下水位抬高1.5m和降低1.5m后，Cl-在第一层含水层中运移1000d，3000d，6000d后的浓度分布图

由以上不同地下水位条件下，Cl-运移相同时间后的浓度分布图可知，内陆地下水位对海咸水中Cl-的运移有显著影响，即内陆地下水位对海陆过渡带的迁移有重要的影响，内陆地下水位越高，海陆过渡带向内陆迁移的范围就越小，速度也越慢，这对控制沿海地区海咸水入侵污染内陆地下水也具有重要的指导意义。

3 结论

该文选择海陆过渡带模型，结合地下水运动理论，选取山东半岛蓝色经济区龙口、招远滨海地区作为研究区域，在对基础地质资料和水文地质资料研究分析的基础上，利用地下水数值模型系统GMS计算该地区的流场，并在内陆地下水超采状态下，海咸水入侵对该地区地下水中Cl-的运移进行数值模拟，研究该地区海陆过渡带迁移的规律。

通过研究，该地区海陆过渡带的宽度约为1.5～4.5km，与1992年1.5～2.5km宽的海陆过渡带相比，该地区的海陆过渡带明显变宽，海咸水继续向内陆运移。此外，内陆地下水位对海陆过渡带的迁移具有重要的影响，内陆地下水位越低，海陆过渡带向内陆迁移的范围越大、速度越快。

因此，为有效遏制海咸水入侵内陆地下水，建议采取限采、回灌、地下水库等相关措施抬高该地区的地下水位，如：在该区域内选择具有足够的天然地下存储空间和充足水源的地区修建地下水库；修建拦蓄补源工程，充分利用雨洪水资源，回灌补给地下水；严格控制地下水开采。