基于Sobol法对辽宁南部海水入侵敏感性分析探究

赵小龙

(辽宁省营口水文局，辽宁 营口 115003)

海水入侵是一个全球性的环境问题，沿海省市均被其困扰。海水入侵使地下淡水体咸化，进而导致淡水资源减少，土壤盐渍化，工业机械腐蚀，甚至危害人体健康，严重制约沿海城市的经济发展[1]。海水入侵的诱发因素有两类，一类是自然因素，包含海平面上升、干旱、潮汐等；另一类是人为因素，包含地下水开采、海水养殖等。WHO指出，若淡水混入约1%体积的海水(氯离子浓度大于250mg/L)则不能直接饮用。由于沿海城市经济发展对地下水的需求以及全球变暖引起的海平面上升等问题的存在，关于海水入侵的研究目前仍是一个热点问题。

1 研究区概述

1.1 自然地理

研究区位于辽东湾东侧的辽宁省营口市,大清河流域望宝山水文监测站的下游，内有长大铁路通过，交通便利，总面积为170km2(陆地面积150km2，潮汐波动带20km2)，东西方向约24km，南北方向2～13km，海岸线约6km。属大陆性季风气候，四季分明，雨热同季。1991—2016年多年平均降水量为600～800mm，年平均日照时数为2600～2880h，多年平均蒸发量1000～1200mm，蒸发极限深度为4m。多年平均气温9～10℃，冬季1月气温最低，平均为-9～-10℃，夏季7月气温最高，平均为24.5～25.0℃，无霜期为180～210天。

1.2 地质条件

营口市属华北地台辽东台背斜营口至宽甸隆起的南翼，受燕山运动的影响形成的千山余脉呈东—西向纵贯本地区。地形自东向西逐渐由高变低。受构造、岩性与新构造运动的影响，地貌形态由东向西呈规律性变化，即低山—高丘陵—低丘陵—滨海平原。

研究区是发源于营口东部山间的大清河入海冲洪积形成的山间河谷平原，主要由一级阶地和二级阶地组成。辽宁省水文地质大队曾对大清河中下游河谷平原进行地质勘查，顺河流而下，共勘探了5个垂直于河流的剖面和一个平行于河流的剖面。根据勘测所得地层剖面，概化出研究区含水介质概念模型,并对研究区地下水含水层特征进行分析。

区内第四纪覆盖层受下伏基岩面起伏影响厚度变化大，变化范围20～65m，整体上呈现出由河谷向两侧丘陵逐渐变薄，北侧丘陵覆盖层由于侧向山谷冲洪积作用靠山一侧覆盖层较厚。由于河流流量历史上随季节气候变化，因此，由河流冲洪积形成的平原含水层性质也具有很大差异，自上而下可以分为5层：第一层是以亚砂土、黏土为主的第四纪覆盖层，含水层类型为潜水含水层，厚度5～10m，渗透能力相对较弱，特别是潮汐波动带的淤泥质海岸，在海潮波动下饱和，涨潮时被淹没，退潮时水分疏干速度慢，加之此处为半日潮，退潮到下次涨潮之间淤泥质海岸中水位变幅基本可以忽略；第二层以砂卵砾石为主，局部沉积亚黏土、亚砂土等，为主要含水层，以潜水为主，局部微承压，该层介质粒径较大，孔隙度较高，渗透系数20～100m/d，自河谷向两侧逐渐减小；第三层以亚砂土为主，局部含有黏土，厚度为0～8m，透水性稍弱；第四层为砾卵石，厚度约5～15m；第五层为黏土层，厚度约3～10m。

1.3 水文地质条件

区内主要的地表水体为大清河地表水系，河床宽度约20～300m。大清河发源于营口市东部山区，流域面积1468km2。为监测大清河的水位、流量、流速等，1959年在大清河的中下游修建了望宝山水文站。根据望宝山水文站多年观测，大清河最大径流量为57m3/s，最小径流量为0.316m3/s，径流随季节变化明显。为调节大清河径流量随季节的变化，在其上游修建了石门水库。此外，在大清河下游修建了1座集蓄水、灌溉、挡潮于一体的拦河闸。

区内潜水含水层的补给来源主要有大气降水入渗、河流侧向补给、山前地下水径流补给、下伏碳酸盐岩岩溶水的顶托补给、山前冲洪积扇的侧向径流补给和开采条件下的越流补给。含水层接受河流侧向补给、大气降水入渗补给等多项补给后，沿地势自东向西流，径流速度随着含水层厚度、透水性和地形的变化而变化[1]，到滨海地带含水介质的颗粒逐渐变小，透水性相应减小，地下水径流也随之减缓。天然状态下，大清河流域在上游河段地下水向河流排泄，在下游河段枯水期受地下水补给，丰水期河水补给地下水。随着工农业发展，用水需求增加，大清河流域先后修建了4个水源地(井深较深，分布在第四层中)，自上游至下游依次是团甸水源地(19口井，其中研究区涵盖4口)、化纤水源地(8口井)、盖州二三水源地(13口井)、永安水源地(18口井)，以及大量农业机电井(井深较浅，分布在第二层)来满足供水需求。自水源地开采井和农业用井修建以来，人工开采已经成为地下水最主要的排泄方式。由于农业灌溉用井和水源地开采井对地下水大量开采，地下水位低于河水水位，研究区地下水不再向河流排泄，而是长时间受河流补给，由于河流补给速度和补给量有限，在水源地开采井和农业灌溉井广泛分布的永安水源地附近形成降落漏斗。

2 典型区选取

营口市海岸属于淤泥质(较细粒底质)海岸，海岸线较长，进行现场调查时主要沿大清河流域开展工作，由上游石门水库开始，沿途搜集相关资料，并现场观测部分数据，直至下游入海处。大清河中游修建有望宝山水文站，此处大清河河床宽约17m，河漫滩两侧为低山丘陵。团甸水源地在望宝山水文监测站附近，共19口开采井(下游4口)，原水源地总开采量为3.0万m3/d，经2013年压采后的开采量减少至0.7万m3/d。沿大清河向下游方向南侧有盖州二水源地与盖州三水源地(简称盖州二三水源地)，共有11口开采井，全部位于大清河南岸，原开采量为3.0万m3/d，压采后现开采量为1.2万～2.0万m3/d。

继续沿岸踏勘可见化纤水源地，共8口开采井，位于大清河流域的南岸，原开采量为1.5万m3/d，压采后现开采量减少为0.5万m3/d。其下游修建有永安水源地，共19口开采井，原开采量为5.0万m3/d，压采后现开采量为1.5万m3/d。此处有微咸水。

3 敏感性分析

3.1 分析方法

在复杂非线性模型中存在较多参数，使得模型计算难以进行，因此需要分析出模型中较为敏感和重要的参数。Sobol法是一种基于方差的全局敏感性分析方法[2]，它可以测试整个输入空间的敏感性，即通过全局方法获取非线性响应和测试参数之间相互作用的影响。

3.1.1 基于水文地质资料建立地下水数值模型

把该模型看作一个函数Y=f(X)，其中，X是一个n个不确定模型输入{X1，X2，…，Xd}的向量[3]，如降雨入渗率、渗透系数、储水率等，Y是一个单值模型输出水头、浓度。同时得到X的先验分布。

3.1.2 产生参数序列

利用蒙特卡罗方法中的拉丁超立方抽样算法，根据模型参数先验信息，抽样得到一个N×2d维的矩阵[AB]，即每行是一个2d维的采样点。

3.1.3 参数点构造

3.1.4 模型计算

3.1.5 计算敏感性系数

通过对Y方差分解，敏感性指数可以划分为一阶敏感性系数和全阶敏感性系数[6]。对于一阶敏感性系数，它主要是由Xi所引起的对输出方差的贡献，所以它主要量测了Xi单独变化的影响，可写为

(1)

式中：Xi为第i个参数；X～i为除去Xi以外的其他参数；Var为输出的方差；E为期望；Si为参数Xi的一阶敏感性系数。

全阶敏感性系数用以量测Xi对输出方差的贡献，包括由Xi与其他输入变量间相互作用所造成的所有任意阶方差[7-9]，可写为

(2)

(3)

(4)

Var(Y)=Var(YA,YB)

(5)

估计值的精度取决于N。N值的选取可先顺序加点，然后计算指数直到估计值达到一定程度可接受的收敛[10-11]。

3.2 敏感性分析过程和结果

在进行敏感性分析之前，将当地年降水量和人工单井开采量(农用井Q1、永安水源地Q2、化纤水源地Q3、二三水源地Q4)作为中间值，同时将中间值减少和增加50%之后的范围作为参数的敏感性分析分布范围，另外给定平均海水位一个波动范围。第二、第四含水层各参数的中间值及分布范围设置见表1、表2。

表1 含水层第二层Sobol敏感性分析参数设置

表2 含水层第四层Sobol敏感性分析参数设置

本次研究中，第二层主要选取平均海水位、农田灌溉井开采量、降水量这3个相关性比较大的参数进行敏感性分析；第四层主要选取下游3个水源地各自开采量作为参数进行敏感性分析。

对第二、第四含水层各自3个参数进行敏感性分析，分别以不同的观测点的Cl-浓度值作为模型输出的Y。在第二层和第四层已入侵部分均匀地选择6个观测点的Cl-浓度。通过Sobol法进行敏感性分析时，将研究区海水入侵模型运行一年，根据6个观测点的Cl-浓度可以分析得出不同位置处各参数全阶敏感性系数STi，见图1、图2。

图1 第二层O1～O6观测点的全阶敏感性系数

图2 第四层O1～O6观测点的全阶敏感性系数

由图1可知：平均海水位对第二层各点位的浓度变化均有影响，O4点距离大清河最近，相应受平均海水位的影响最大。由图2可知：农用井开采量对O4～O6点影响较小，对O1～O3点影响较大。由于研究区大清河北岸平原面积小，农用井数量少，大清河南岸河谷平原较平坦宽阔，农用井分布多，总开采量较大，两岸之间又有大清河相隔，一定程度上干扰了水力联系。因此，南岸观测点相对敏感。

4 结 语

本文应用Sobol法建立地下水数值模型，对辽宁南部大清河海水入侵敏感性进行分析。通过模型结果分析可知：该区发生海水入侵的主要原因来自两方面，一方面是潮汐作用下的海水沿河道上溯，并向四周补给低水位地下水所产生的入侵；另一方面是由于地下水的大量开采形成降落漏斗，从而导致降落漏斗中心水位远低于海平面的平均水力坡降，海水沿水力坡降方向入侵。水源地进行压釆后，农业用水成为研究区地下水的主要排泄方式，为加速海水入侵回退速率，应对研究区农业灌溉用井的开采量加以控制。总体上，研究区进行压釆后地下水位抬升，原来指向内陆的水力梯度重新指向海洋，海水入侵回退，回退速率由水力梯度大小决定，回退比例与回退时间呈线性相关。大量人工开采是该区产生海水入侵的主要原因，控制开采量是回退海水入侵最直接的办法。研究区淤泥质海岸由于其弱透水性的特性，不利于产生有利海侵回退的水力梯度。若想海水入侵快速回退，应在近海岸处进行处理(如大量抽水使水位下降等)，使得由内陆指向海洋的水力梯度变大。