苏北沿海三市三维地下水流数值模拟

陈 雄，张 岩，王艺伟，叶淑君，吴吉春，于 军，龚绪龙

1.南京大学地球科学与工程学院，南京 210023 2.国土资源部地裂缝地质灾害重点实验室(江苏省地质调查研究院)，南京 210018

0 引言

地下水是水资源的重要组成部分，对国民经济的发展起着举足轻重的作用，然而不合理的开采[1]和过量抽取[2]地下水资源会引起地面沉降[3-7]、地裂缝[8-11]、海(咸)水入侵[12-14]、土地荒漠化[15]及地下水污染等一系列环境地质问题[16]。自20世纪80年代起，随着社会经济的快速发展，苏北沿海地区对地下水资源的需求量与日俱增[17]，三市深层地下水(第Ⅱ、Ⅲ、Ⅳ承压含水层)的大规模开采造成部分地区水位持续下降，引发了地面沉降、水质咸化等环境地质问题，严重影响当地社会、经济的可持续发展[18-22]。

凌家荣等[23]、骆祖江等[24-26]、单卫华[21]、徐玉琳[27]、安晓宇等[28-29]先后尝试用建立数学模型的方法来研究部分区域的地下水资源及地下水环境问题，认识到过度开采是导致地下水位下降的主要原因。随着地下水大量开采，地下水流场发生巨大改变，区域漏斗连成一片，建立统一的地下水流模型非常必要且更为合理。本文将苏北沿海地区三市的区域地下水流系统作为研究对象(图1)，建立三维地下水流数值模型，进行研究区三维地下水流的模拟和预报，以期为合理开采地下水资源、防治地面沉降、解决水质咸化等环境地质问题提供科学支撑。

1 研究区水文地质条件

研究区自新生代以来受到河流的搬运与海洋沉积作用的影响，堆积了厚达200～1 600 m的松散沉积层，形成了一套巨厚的松散孔隙含水系统[20,23]。根据岩性和时代，自上而下划分为：全新统孔隙潜水含水层，上更新统第Ⅰ承压含水层(Ⅰ含)，中更新统第Ⅱ承压含水层(Ⅱ含)，下更新统第Ⅲ承压含水层(Ⅲ含)，中、上更新统第Ⅳ承压含水层[21](Ⅳ含)。各个含水层间以弱透水的黏土层相分隔(图2)。由于潜水及Ⅰ含中的地下水主要为水质较差的咸水和微咸水，长久以来该地区的地下水开采主要集中在Ⅱ含、Ⅲ含、Ⅳ含[18,26,30]。所以本次研究模拟Ⅱ含、Ⅲ含、Ⅳ含的深层地下水。

钻孔及水工实验资料显示，地下水系统为非均质各向异性的多层结构，各个含水层通过其间的弱透水层发生水力联系；故本文将含水层与弱透水层作为统一的复合含水系统来处理，将其概化为三维空间上的非均质各向异性介质，含水层之间的弱透水层也按独立的层位参与计算，整个地下水流概化为三维非稳定流[31]。

Ⅰ--Ⅰ′为水文地质剖面线。图1 研究区地理位置图Fig.1 Location map of the study area

Ⅰ含、Ⅱ含、Ⅲ含、Ⅳ含分别表示第Ⅰ、第Ⅱ、第Ⅲ、第Ⅳ承压含水层;Ⅰ弱、Ⅱ弱、Ⅲ弱、Ⅳ弱分别表示第Ⅰ、第Ⅱ、第Ⅲ、第Ⅳ弱透水层。图2 研究区水文地质剖面图Fig.2 Hydrogeological profile of the study area

自然状态下，深层地下水自西向东缓慢流动，接受西部地下水流的补给，并向东部海洋流动排泄；在大规模开采地下水后，人为开采成为主要排泄方式，地下水由四周向开采强烈的降落漏斗区流动。研究区北部基岩已侵入至Ⅱ含，故将Ⅱ含概化为隔水边界，其余侧向边界在有观测孔的部分设定为一类给定水头边界，其他设定为第二类流量边界(图3a)。此次研究中：将Ⅰ含作为模拟的三维含水系统的上部边界，由于Ⅰ含水质差，开采利用少，观测数据显示该层水位年变化幅度很小，所以Ⅰ含可概化为一类定水头边界；含水系统底部即Ⅳ含底部为封闭性良好的黏土层或基岩，故Ⅳ含概化为隔水边界。大规模的人工开采为该区地下水的主要排泄方式，将各乡镇的开采强度概化为大井[23]。各个乡镇的年开采量由实际统计获得，Ⅲ含开采井的分布如图3b所示。

2 数学模型

以连续性原理及达西定律为基础，结合研究区的水文地质条件，建立相应的三维非稳定流数学模型[32-34]：

式中：Ω为渗流区域；t为时间(s)；h为地下水位(m)；Kxx,Kyy为水平渗透系数(m/d)，Kzz为垂向渗透系数(m/d)；W为源汇项(d-1)；h0为初始水位(m)；h1为一类边界上的水位(m)；SS为贮水率(m-1)；Γ1为一类边界；Γ2为二类边界；Kn为二类边界法线方向渗透系数(m/d)；q为二类边界上的流量(m/d)。

a．边界条件； b．开采井分布。图3 研究区Ⅲ含边界条件及开采井分布Fig.3 Boundary conditions and wells of the third confined aquifer in the study area

3 数值模拟

研究区使用地下水模拟系统(groundwater modeling system，GMS)中的MODFLOW模块来建立地下水数值模型。GMS是目前国际上通用的地下水模拟软件，可用来模拟饱和、非饱和流条件下的地下水流和溶质运移等问题[35]。本次数值模型采用矩形剖分，垂向上剖分为6层，分别对应含水系统的Ⅱ弱、Ⅱ含、Ⅲ弱、Ⅲ含、Ⅳ弱、Ⅳ含；平面上每层剖分成300行和300列(网格长1 292 m，宽1 146 m)，单元总数540 000，其中活动单元数128 510个(图4)。模型采用有限差分法求解，并采用预处理共轭梯度迭代法(PCG)联立迭代求解代数方程组。时间上，模型模拟时间为2005-01-01—2013-12-31，其中，2005-01-01—2010-12-31为模型识别期，2011-01-01—2013-12-31为模型验证期。一个季度为一个应力期，共36个应力期，每个应力期90 d。各含水层的初始水位采用2004年末的实测水位，弱透水层的初始水位取上下相邻含水层水位插值获得[24]，初始流场见图5，边界条件设置见图3(以Ⅲ含为例)。苏北沿海三市地下水开采经历了开采、超采、压缩3个阶段，2000年开始全面实施地下水限采， 2005—2013年各个深层承压水的开采量逐年减少，从2005年的2.03亿m3减少至2013年的1.56亿m3左右。

模型采用“试错法”(trial-and-error)来间接率定水流模型的相关参数[36]。对含水层和弱透水层进行参数分区，给定初始值，然后根据实测水位与模拟计算水位的对比，不断调整参数，直到误差达到允许范围。各含水层均具有一定数量的水位观测孔来进行水位拟合，其中Ⅱ含有72个水位观测孔、Ⅲ含有122个水位观测孔、Ⅳ含有67个水位观测孔，基本控制全区(图6)。除了观测孔水位的拟合，还对识别期末和验证期末各含水层流场进行了拟合。表1给出各含水层与弱透水层参数值范围，表2、表3给出了第25应力期(2011年)各含水层的流入量与流出量。各含水层部分观测孔的水位拟合情况见图7；所有观测孔的模拟水位与观测水位的对比关系见图8。识别期末与验证期末各含水层的模拟流场与实测流场对比如图9所示，识别及验证结果均表明，模型计算值与观测值拟合良好，水位绝对误差不超过1.00 m，且模型计算所得的含水系统的地下水流入量与流出量基本相等(表2、表3)，说明所建地下水值数学模型能够刻画研究区的地下水流问题。

a.平面剖分；b.横向剖面(第165行)剖分；c.纵向剖面(第142列)剖分。图4 研究区单元剖分Fig.4 Units division in the study area

图5 研究区Ⅲ含、Ⅳ含初始水位Fig.5 Initial water levels of the confined aquifers Ⅲ and Ⅳ in the study area

a．所有观测孔；b．本文中列举的观测孔。图6 研究区含水层水位观测孔分布Fig.6 Distribution of all water level observation hole in the study area

含水层/弱透水层Kxx /(m/d)Kyy /(m/d)Kzz /(m/d)Ss /m-1Ⅱ弱Ⅱ含3.00×10-7～8.00×10-50.08～40.003.00×10-7～8.00×10-50.08～40.001.00×10-8～8.50×10-60.01～4.001.00×10-5～2.50×10-32.00×10-6～8.50×10-4Ⅲ弱5.00×10-8～7.00×10-55.00×10-8～7.00×10-55.00×10-8～7.00×10-51.00×10-5～9.00×10-4Ⅲ含1.00～40.001.00～40.000.10～4.008.50×10-6～9.85×10-4Ⅳ弱1.00×10-8～1.00×10-51.00×10-8～1.00×10-51.00×10-8～1.00×10-51.00×10-5～8.00×10-4Ⅳ含0.45～45.000.45～45.000.05～4.505.50×10-7～8.60×10-4

表2 研究区2011年各层平均流入量

表3 研究区2011年各层平均流出量

识别期末及验证期末各个含水层整体拟合效果较好，2005—2013年研究区深层地下水开采量总体处于下降趋势，Ⅱ含水位为-20.00～-5.00 m。由表2、表3经计算可知：Ⅱ含接受的流入量小于流出量，储存量在2011年日平均净减少速率为48 419 m3/d；2005—2013年区域水位平均下降2.55 m，下降速率为0.28 m/a(表4)。

因Ⅱ含不是南通的目标开采层，故南通大部分区域的水位在-10.00 m左右，水位变化幅度较小；但是南通地区Ⅱ含越流补给大规模开采的Ⅲ含，启东、海门处的水位有一定程度的下降。在模型识别期，盐城市北部水位下降5.00 m左右，部分地区的水位降幅接近10.00 m，如滨海县0229观测孔处水位从2005年初的-21.94 m下降至2010年底的-27.23 m，降幅为5.29 m(图7b)；盐城南部由于开采强度减少，水位回升，如盐都区0221观测孔处水位从2005年初的-32.73 m回升至2010年底的-28.70 m，水位回升4.03 m(图7c)。模型验证期，北部水位持续下降，但降幅明显减小；0229观测孔处水位从2011年初的-27.23 m下降至2013年底的-28.30 m，降幅为1.07 m(图7b)，滨海县-25.00 m漏斗中心范围扩大较为明显(图9a、9b)；南部水位持续回升且升幅趋于平缓，0221观测孔处的水位从2011年初的-28.90 m回升至2013年底的-28.70 m，升幅为0.20 m(图7c)，-20.00 m水位漏斗消失，水位维持在-15.00～-20.00 m之间；连云港市东部燕尾港在识别期末形成-25.00 m降落漏斗，验证期漏斗中心水位持续下降，漏斗范围扩大(图9a、9b)。

Ⅲ含是深层地下水开采的主要含水层，始终处于超采状态，储存量在2011年的日平均净减少速率为240 500.6 m3/d(表2、表3)，因此研究区大部分地区水位处于持续下降状态，2005—2013年区域水位平均下降4.14 m，下降速率为0.46 m/a(表4)。模型识别期：在盐城市滨海县、射阳县、建湖县及盐都区形成4处-30.00 m水位漏斗，其中以建湖县水位漏斗的降幅最大，水位从2005年初的-10.25 m下降至2010年末的-32.82 m，降幅达22.57 m(图1，图9c)；-25.00 m水位等值线覆盖盐城市大部分区域并向外扩展与连云港市灌南县相连，灌南县大部分地区水位位于-20.00 m至-25.00 m之间(图9c)；南通市如东漏斗处的平均水位从2005年初的-43.85 m回升至2010年底的-40.25 m，其余地区水位均在下降，水位平均下降5.10 m(图9c)；海门市三厂镇原有的-35.00 m水位漏斗发展形成-40.00 m的水位漏斗，漏斗中心水位从2005年初的-37.37 m下降至2010年末的-44.80 m，降幅为7.43 m(图7i)，漏斗的范围明显扩大(图9c)。验证期开采量虽有下降但开采量仍然巨大：灌南县大部分地区水位已低于-25.00 m(图9d)；盐城市各漏斗中心水位继续下降，漏斗范围继续扩大，盐都区和建湖县-30.00 m水位漏斗已扩大连成一片(图9d)；滨海县漏斗范围扩大；南通如东漏斗平均水位继续回升至2013年底的-38.85 m，海门市漏斗中心2013年底的水位为-44.84 m(图7i)，变化不大，但两处漏斗的范围均有不同程度的扩展。

a．东台市Ⅱ含0113观测孔；b．滨海县Ⅱ含0229观测孔；c．盐都区Ⅱ含0221观测孔；d．射阳县Ⅱ含0232观测孔；e．大丰市Ⅲ含0226观测孔；f．滨海县Ⅲ含0217观测孔；g．阜宁县Ⅲ含0090观测孔；h．南通市Ⅲ含7020观测孔；i．三厂镇Ⅲ含3002观测孔；j．东台市Ⅳ含0187观测孔；k．盐城市Ⅳ含0197观测孔；l．射阳县Ⅳ含0224观测孔。图7 研究区部分观测孔水位拟合图Fig.7 Fitting graphs of partial water level observation holes in the study area

图8 研究区观测孔水位模拟值与观测值对比图Fig.8 Comparison graph of simulated and observed values of water level in observation holes in the study area

Ⅳ含水位的模拟值基本上与观测值一致。该层在模拟期开采总量变化不大，2011年储存量日平均净减少速率为92 224.7 m3/d(表2、表3)，部分地区发生向上越流补给第Ⅲ承压含水层，地下水位在-40.00～-20.00 m(图9e、f)，除了原有的-40.00 m水位漏斗中心附近的水位从2005年初-49.81 m上升至2013年底-40.53 m(图7k)，其余地区的水位均处于下降状态，水位平均下降2.60 m，下降速率为0.29 m/a(表4)。

表4研究区2005--2013年各含水层区域水位平均变化量

Table4Averagewaterlevelvariationofeveryconfinedaquiferfrom2005to2013inthestudyarea

含水层平均下降量/m下降速率/(m/a)Ⅱ含2.550.28Ⅲ含4.140.46Ⅳ含2.600.29

4 模型预测

建立地下水模型的目的是为了预报在拟定的地下水开采方案下地下水位的发展，通过预报，可以了解未来地下水开采方案是否合理以及不同开采方案下地下水位的响应情况。苏北沿海三市从2000年开始全面实施地下水限采，并逐步实施地下水开采计划制度，如2005年开采约2.03亿m3，2013年已减少至1.56亿m3(图10)。运用上文经过识别验证后建立的研究区三维地下水流数值模型建立预测模型，初始条件为2013年末的地下流场，边界条件则均为二类边界，预报时间为2013-12-30—2020-12-30，共分28个应力期。

采用两种开采方案预测。

第1种方案，即现状开采方案：在现状开采条件下，保持现有2013年开采井的数量、开采量及位置不变，每年维持总开采量1.56亿m3。

第2种方案，即限制开采方案：在苏北三市2013年地下水开采方案的基础上，规划至2020年地下水压缩开采方案(表5)(此方案由江苏省水利厅及江苏省地调院规划制定)。其中：盐城市2013年深层地下水开采量为8 718.39万m3， 2020年开采目标为4 365.78万m3，减少开采4 352.61万m3，减少49.92%；南通市2013年深层地下水开采量为

a．2010年12月Ⅱ含；b．2013年12月Ⅱ含；c．2010年12月Ⅲ含；d．2013年12月Ⅲ含；e．2010年12月Ⅳ含；f．2013年12月Ⅳ含。图9 研究区含水层水位拟合图Fig.9 Fitting graphs of water level in aquifers in the study area

柱状图表示各个城市3个含水层的地下水开采总量。图10 研究区深层地下水历年开采量图Fig.10 Graph of total exploitation per year of deep groundwater in the study area

104m3

5 508.55万m3，2020年开采目标为4 168.89万m3，减少开采1 339.66万m3，减少24.32%；连云港市2013年深层地下水现状开采量为1 426.60万m3，2020年Ⅱ、Ⅲ含开采目标为613.86万m3，减少开采812.74万m3，减少56.97%。开采总量从2013年的1.56亿m3减少到2020年的0.91亿m3。

现状开采方案下，预报期间各承压含水层最低地下水位变化见表6(均为各年份12月30日水位)，各承压含水层在预测期末的流场如图11所示。

在现状开采条件下，各含水层水位总体继续下降。2014—2020年：Ⅱ含区域水位平均下降1.59 m，下降速率为0.23 m/a(表7)；盐城市滨海县、响水县与连云港灌南县的-25.00 m水位漏斗连成一片，漏斗中心的水位降至-34.70 m，其余多个漏斗都有一定程度的发展，南部大丰市地区的水位有一定程度的回升，南通市在通州区形成-15.00 m的水位漏斗(图11a)。Ⅲ含区域水位平均下降2.31 m，下降速率为0.33 m/a(表7)；原有的响水县-35.00 m水位漏斗向外扩展至灌南县，盐城市区附近的多个水位漏斗连成一片，形成-40.00 m的水位漏斗，漏斗中心的水位降幅超过5.00 m，南通市原有的如东县-35.00 m水位漏斗和海门市-40.00 m水位漏斗因开采强度较历史时期降低，漏斗范围缩小，水位回升(图11b、11d)。Ⅳ含区域水位平均下降1.57 m，下降速率为0.22 m/a(表7)；漏斗中心的最低水位下降至-46.50 m(表6)，各主要漏斗范围均在缓慢增大(图9e,9f,图11c)。

表6 现状及限采方案下研究区各含水层最低地下水水位

a．Ⅱ含水位；b．Ⅲ含水位； c．Ⅳ含水位；d．Ⅲ含水位变幅。图11 现状开采方案下研究区2020年底各含水层等水位线图及水位变幅Fig.11 Water level contour map of every aquifer and the variation range map of every confined aquifer under current exploitation scheme in the end of 2020 in the study area

限制开采方案下，预报期间各承压含水层最低地下水位变化见表5，各承压含水层在预测期末的流场如图12所示，预测期末(2020年)的水均衡分析见表8、表9，此时各含水层的开采量已降至最低。和现状开采方案相比，限采方案效果明显，水位持续下降得到有效控制，多个漏斗区水位有一定的回升，但是继续开采及水流汇聚补给漏斗区导致漏斗相邻近区域水位下降；总体来说，3个含水层的净储存量变化速率仍为负值(表8、表9)，说明地下水资源量处于进一步消耗中。2014—2020年：Ⅱ含区域水位平均下降1.08 m，下降速率为0.15 m/a(表7)，水位降速减少34.78 %；连云港市和盐城市部分地区水位回升比较明显，燕尾港-25.00 m水位漏斗消失，阜宁县、滨海县-25.00 m水位漏斗水位回升明显，盐城市区-20.00 m漏斗消失，南通市因Ⅱ含越流补给Ⅲ含，水位有小幅度的下降(图12a)。Ⅲ含各漏斗区水位回升，原有的水位持续下降、降落漏斗持续扩大等现象得到控制，该含水层区域水位平均下降1.15 m，下降速率为0.16 m/a(表7)，水位降速减少51.52 %；盐城市建湖县、盐都区漏斗中心水位回升超过5.00m，南通市如东漏斗区水位持续回升，海门市-40.00 m水位漏斗范围缩小，漏斗中心水位从-44.84 m回升至-41.60 m，其余地区的水位也有一定程度的回升(图12b、12d)。Ⅳ含总体水位处于下降状态，区域水位平均下降1.07 m，下降速率为0.15 m/a(表7)，水位降速减少31.82 %；主要漏斗区水位有一定的回升，原来的漏斗中心处水位回升至-37.80 m(表6)。总而言之，实施地下水限制开采对减缓了苏北沿海地区的地下水位的持续下降有一定的可行性，但为了彻底扭转该区域的地下水位的下降，还需要进一步加大地下水限制开采的力度。

表7现状及限采方案下研究区2014--2020年各含水层区域水位平均变化量

Table7Averagewaterlevelvariationofeveryconfinedaquiferfrom2014to2020undercurrentandlimitexploitationschemeinthestudyarea

含水层平均下降量/m下降速率/(m/a)现状限采现状限采Ⅱ含1.591.080.230.15Ⅲ含2.311.150.330.16Ⅳ含1.571.070.220.15

表8 限制开采方案下研究区2020年各层平均流入量

表9 限制开采方案下2020年各层平均流出量

a．Ⅱ含；b．Ⅲ含；c．Ⅳ含；d．Ⅲ含。图12 限制开采方案下研究区2020年底各含水层等水位线图及水位变幅图Fig.12 Water level contour map of every aquifer and the variation range map of the 3rd confined aquifer under limited exploitation scheme in the end of 2020 in the study area

5 结论

1)在充分研究所收集和整理的地质、水文地质、地下水开采、地下水位监测、抽水试验等资料和数据基础上，建立了包括3个承压含水层和3个弱透水层在内的苏北沿海三市深部地下水系统三维非稳定流数值模型，模型通过识别和验证后，能很好地刻画研究区的Ⅱ含以下深部地下水流系统，能正确模拟研究区地下水流，可用于预报拟定开采方案下地下水流场的变化。

2)模型模拟结果显示：2005—2013年随着总的地下水开采量逐年减少，地下水位下降幅度得到有效控制，地下水位下降速率较历史峰值有明显降低，甚至在局部区域出现水位回升。Ⅱ含区域水位平均下降2.55 m，下降速率为0.28 m/a，盐城市南部地区水位有一定回升；Ⅲ含是深层地下水开采的主要含水层，始终处于超采状态，2005—2013年区域水位平均下降4.14 m，下降速率为0.46 m/a，南通市如东漏斗处的平均水位从2005年初的-43.85 m回升至2013年底的-38.85 m；Ⅳ含区域水位平均下降2.60 m，下降速率为0.29 m/a；主要漏斗区水位则有一定的回升。

3)在现状开采及限制开采方案下，2014—2020年进行了预测。在现状开采条件下：含水层水位继续下降，但下降速率降低；3个含水层的区域水位平均分别下降1.59、2.31及1.57 m，水位下降速率分别为0.23、0.33及0.22 m/a，分别比模拟期减少了17.85 %、28.26 %和24.14 %。在限制开采条件下：随着总的开采量从2013年的1.56亿m3压缩到2020年的0.91亿m3，地下水水位下降趋势明显变缓，3个含水层的区域水位平均下降分别为1.08、1.15及1.07 m，水位下降速率分别为0.15、0.16及0.15 m/a，分别比模拟期减少了46.43%、65.21%和48.28%，水位下降得到有效控制，局部地区水位明显回升，原有的主要水位漏斗均有很好的恢复，部分漏斗消失。但是为了彻底扭转该区域的地下水位的下降，还需要进一步加大地下水限制开采的力度。