考虑不同海岸类型的深圳市东海岸海水入侵多尺度耦合模拟

宋绍宇,成建梅*,程天舜,卢薇艳,职敬儒

(1.中国地质大学(武汉) 环境学院,湖北 武汉 430078;2.深圳市水务规划设计院股份有限公司,广东 深圳 518023;3.深圳市地质局,广东 深圳 518023)

全世界许多国家和地区都出现了海水入侵问题,对海水入侵的研究也引起了众多学者的关注[1-4]。数值模拟技术是定量研究海水入侵发展演化的重要手段,常见的海水入侵数值模拟模型包括突变界面模型和变密度过渡带模型[5-7]。目前常采用变密度过渡带模型来模拟并解决实际海水入侵问题,主要包括基于单个组分或总矿化度的考虑变密度海水入侵三维溶质运移模型[8-9]、考虑海咸水运移过程中多个水化学组分的三维水流-水质模型[10-11]和考虑耦合水化学反应过程的高浓度海水入侵三维水质模型[12]。滨海含水系统非均质结构与海水入侵优势通道的刻画、不同类型临海边界的处理、多尺度的耦合是实现高精度海水入侵过程模拟的关键[4,13-15]。面对不同的研究目标,可选择的海水入侵数值模拟模型类型不同。吴吉春等[16]提出要以不同类型海岸带的海水入侵发生、演化、防治为主线,探究不同类型海岸带海水入侵发生的机理与作用机制,并探求高效的海水入侵问题管理对策。此外,针对滨海含水层的复杂性和海岸类型的差别,一些学者对单一类型海岸带(如砂质海岸、裂隙基岩海岸或灰岩海岸等)海水入侵特征开展了数值模拟研究[1,3,17-20],但极少的模拟案例能涉及多种海岸类型[21],无法精确表征实际复杂的海岸带结构特点和海水入侵演化规律。

20世纪80年代以来,随着深圳市的改革开放和城市迅速发展,地下水过量开采导致深圳海岸带发生了严重的海水入侵,加之受城镇化建设、填海造陆、海水养殖等一系列海岸带人类工程活动的影响,使得该地区海水入侵的发展演化机制更为复杂。目前已有的研究多集中在利用水化学、同位素及电阻率等指标对深圳海岸带海水入侵过程进行识别和评估,涉及数值模拟模型的研究极少[22-23],且研究区仅限于深圳西海岸。而深圳市东部出露有河口海岸、基岩海岸、砂质海岸等多种海岸带类型,海水入侵方式更为复杂多样。为此,本文以深圳市东海岸为研究区,充分发挥区域渗流对局域渗流的总体控制优势,采用多尺度耦合模拟方法,建立了研究区区域尺度及典型海岸带局域尺度的三维变密度溶质运移模型,并通过数值模拟方法探究了研究区不同类型海岸带的海水入侵机理与地下水水质演化规律。

1 研究区背景条件

研究区位于广东省深圳市东部沿海区(图1),北与惠州市、坪山区相接,以山前为边界,西至沙头角河与深圳河,东部和南部靠海,与大鹏湾和大亚湾相望,面积约为361 km2。区内地势总体北高南低,发育有山地和丘陵,南侧有海积平原。区内较大的河流有盐田河、葵涌河、王母河和东涌河等,河流两侧及滨海地带的海积平原发育有第四系松散地层,岩性以中细砂、粗砂及卵石等为主,赋存丰富的孔隙水,单井出水量为100～150 m3/d。区内发育泥盆系中统鼎壶山群砂岩、白垩系-侏罗系花岗岩、侏罗系上统高基坪群酸性火山岩、石炭系下统石磴子组灰岩等基岩地层,赋存风化裂隙水和构造裂隙水,单井出水量为50 m3/d,水量相对贫乏。区内第四系潜水含水层分布不规律且不连续,而其下的基岩裂隙含水层分布广泛且具有统一性,其隔水顶板多在海平面之下。区内地下水补给以大气降水、河流侧渗和上游水库远源补给为主,地下水径流方向主要受地形控制,总体由北至南径流,最终向海排泄。

图1 研究区含水岩组结构简图Fig.1 Structure diagram of water-bearing rock group of the study area

深圳市东部海岸线长约123 km,海岸带类型复杂。其中,砂质海岸主要分布在东涌、西涌地区以及大鹏澳沿岸,地形平缓,地貌上属于冲洪积平原及海积平原,如图1中c-d、e-f剖面所示,岩性以第四系砂砾石层为主,透水性较强;基岩海岸主要分布在东部大鹏湾和大亚湾西侧的大鹏半岛沿岸,长度约为81.9 km,地形上以山地、丘陵为主,基岩裂隙相对发育;河口海岸依托研究区内河流发育,如图1中a-b剖面所示,葵涌河河口与海洋相接,且入海口花岗岩基岩出露,葵涌河成为内陆地下水与海水的重要交换媒介。深圳市东部海岸带的海水入侵咸淡水界面如图1所示,海水在不同类型海岸段入侵程度不同,其中基岩海岸的海水入侵距离以300 m到1 000 m为主,海水入侵主要受到地形、岩性和裂隙构造发育等因素影响;而砂质海岸和河口海岸的海水入侵范围较大,海水入侵距离受第四系分布和感潮河流水位变幅的影响较大。

2 研究区海水入侵多尺度耦合模型构建

2.1 多尺度耦合概念模型建立

根据深圳市东海岸涉及的地层岩性特征及透水性强弱,将表层覆盖的人工填土、海陆交互成因的淤泥和黏土以及全风化的基岩概化为弱透水层,将第四系松散沉积物砂、卵石以及强风化、中风化基岩概化为含水层,将研究区地层概化为9层42区,体现为“三弱三含”的地层结构特征。在区域模型的边界概化上,研究区向西延伸至梧桐山,以相邻深圳河、沙头角河以及小型的地下分水岭为边界;研究区东部和南部靠海,在模型中处理为定水头边界和定浓度边界;研究区北部地势较高,岩性出露为侏罗系下统桥源组砂岩、石英砂岩和花岗岩,此区域岩石的渗透系数较小,选取地下分水岭为模型北部的隔水边界;模型上边界为潜水面,接受大气降水的补给,模型底部为侏罗系-白垩系微风化花岗岩,由于其风化程度较低,渗透性较差,处理为零通量边界。

渗透性强的砂砾石层和裂隙发育的基岩破碎带是海岸带海水入侵内陆含水层的重要通道,从而加剧海水的入侵[24]。海岸带依据海岸带底质和空间形态的不同划分为基岩海岸、砂质海岸、淤泥质海岸、生物海岸和河口海岸[25]。因此,为了更好地探究研究区不同类型海岸带的海水入侵发展趋势,在区域尺度模型的基础上选择砂质海岸、基岩海岸和河口海岸3种不同类型海岸带建立局域尺度模型,如图2所示。

在局域尺度模型的边界概化上,基岩海岸带模型范围选择研究区大鹏半岛的东北部基岩海岸带,区内地势较高,主要出露有泥盆系双头群和鼎湖山群的石英砂岩和粉砂岩,构造裂隙发育,地下水类型为层状基岩裂隙水;选择研究区葵涌河地区建立河口海岸带模型,区内主要出露有石磴子组结晶灰岩、白云质灰岩及少量壶天群大理岩,发育数量较多、规模较大的岩溶裂隙和溶洞,成为地下水理想的储运空间,且沿河流两侧发育有第四系沉积物,在入海口处有花岗岩出露,随着潮汐波动,葵涌河成为感潮河流,海水在大潮时期入侵内河,并侵染灰岩含水层;选择研究区东涌和西涌地区的海积平原区构建砂质海岸带模型,区内主要分布有第四系海相沉积物、鼎湖山群石英砂岩和白垩系花岗岩,地下水类型为第四系孔隙水和基岩裂隙水。

2.2 数学模型建立

过渡带海水入侵的数学模型由两个偏微分方程及其定解条件来描述,其中一个偏微分方程用于描述受密度变化影响的混合液体(咸淡水混合物)的流动,另一个偏微分方程用于描述混合液体中盐分的运移[26]。考虑到研究区不同海岸带含水层分布的特点,这里采用考虑含水层各向异性的变密度三维水流模型(模型I)和变密度三维对流-弥散模型(模型Ⅱ)来描述海水入侵过程,模型I和模型Ⅱ分别表示如下。

模型Ⅰ:

模型Ⅱ:

式中:Kij为含水层渗透系数(m/d);xi、xj为坐标系(i、j=1,2,3);C为液体浓度(kg/m3);C0为液体初始浓度(kg/m3);C1为源(汇)项浓度(kg/m3);H为求解空间的位置水头(m);t为时间(d);H0为初始水头(m);HB为边界Γ1给定水头(m);ej为重力方向单位矢量第j个分量;φ为孔隙度;η为密度耦合系数;μs为比储水系数;μd为重力给水度;ρ0为淡水密度(kg/m3);ρ为混合液体密度(kg/m3);vi为地下水渗透速度在xi上的分量(m/d);q为单位体积多孔介质中抽取或注入的水量(1/d);ni为边界Γ2、Γ2-1、Γ2-2在xi轴方向上的法向单位矢量;ρB2为定通量边界上流入或流出的水的密度(kg/m3);vB2为单位截面积的定通量边界上流入或流出的水量(m3/d);ρ1为潜水面变动带地下水密度(kg/m3);H1为潜水面Γ2-1上的参考水头(m);x3为潜水面Γ2-1上的水头高程;W为潜水面单位面积上的降水入渗补给量(m3/d);Γ1、Γ2、Γ2-1、Γ2-2分别为浓度给定边界、隔水边界、潜水面边界和弥散通量边界;ui为孔隙水平均流速(m/d);Dij为弥散系数张量;CB为边界Γ1上的浓度(kg/m3);C2为降水或灌溉用水等入渗水的浓度(kg/m3);C3为定通量边界流入或流出水的浓度(kg/m3)。

2.3 数值建模

为了探究研究区不同类型海岸带的海水入侵演化规律,先对研究区整体区域进行数值建模,在完成区域水流-水质模型的识别与验证后,采取参数套用的方式初步建立研究区不同类型海岸带的局域尺度模型;然后根据区域尺度模型模拟过程中的流场和浓度场计算出局域尺度模型的边界条件,并根据不同类型海岸带的特点对局域尺度模型海岸带的参数进行精细化处理,并调整稳定流模型使其水均衡稳定,最终完成研究区不同类型海岸带局域尺度模型的建立。

2.3.1 区域尺度模型建立

在滨海含水系统的模型概化中,临海边界条件的正确处理直接影响到海水入侵模型的仿真程度。目前绝大多数海水入侵模型将平面图的海岸线作为含水系统的边界,但实际上海岸近海地形一般为缓坡向外海倾斜,延伸的距离至少有几十米甚至超过百米。在本模型中,通过查阅前人资料,并综合潮汐数据选定距海岸线150 m处作为模型向海一侧外边界,且处理为定水头边界。

依据数学模型,利用Feflow软件建立数值模型并求解。使用Triangle算法对研究区进行三角网格剖分,对研究区重点河流、断层、长期观测孔进行剖分单元加密,并且保证观测孔位置在结点上。根据概念模型中对整个研究区的垂向分层及横向分区(图3),考虑含水层的渗透系数K、重力给水度、比储水系数、纵/横向弥散度、扩散系数等水文地质参数,用来刻画各向异性含水层的变密度三维溶质运移过程。

图3 研究区部分含水层参数分区图Fig.3 Hydrogeological parameter zones in some aquifers in the study area

在模型的源汇项处理方面,根据土地利用类型将研究区划分为填海区、建筑区、植被区、水域、灌木和耕地6种类型的降雨入渗分区[图3(a)],并收集研究区19个雨量站的日降雨量数据,采用泰森多边形与降雨入渗补给系数叠加的方式赋于模型;将模型中河流和水库概化为柯西边界(第三类边界Fluid-transfer边界),根据深圳市海水入侵地质灾害调查资料[27],设定地表水与地下水水力交换注入传输率Kin为0.05/d,抽取传输率Kout为0.10/d,并收集研究区内罗屋田、铁扇关门、枫木浪等水库的日水位变化数据,通过时间序列的方式赋于模型,同时收集研究区内盐田河、葵涌河、王母河等河流河口及中上游的日水位数据,通过反距离插值法且参考河流地表高程生成河流各结点水位。

潮汐作用会使得咸淡水交换过程更为频繁,在涨潮时海水由高潮线附近的海滩渗入潜水含水层,落潮时由低潮线附近的海滩出流,这种潮汐水位的波动引起的海水-地下水循环模式会对海水入侵过程产生重要的影响。本文选取大鹏湾区域潮汐观测站每日6:00、12:00和18:00的潮汐水位和大亚湾区域潮汐观测站每日12:00和18:00的潮汐水位以线性时间序列的方式作为模型的海洋边界。模型的初始流场和初始浓度场通过Surfer软件对收集到的实际数据进行克里金插值,并与稳定流模型结果进行对比后确定出更接近实际情况的初始流场和初始浓度场。

模型的识别与验证是地下水数值模拟中极为重要的工作,为了验证区域尺度模型的参数设置是否合理,本次选取的模型识别期为2019年10月31日至2020年10月31日,收集研究区内6个长期观测孔的实测水位数据,观测间隔为5 d,并收集2020年10月下旬的114个观测点的实测Cl-浓度值进行内插,得到识别期末区域实测Cl-浓度场。在模型的识别过程中,对研究区的水文地质参数、边界条件以及源汇项等进行反复调整,并将研究区长期观测孔的模拟水位与实测水位变化进行了比较,同时将识别期末研究区模拟Cl-浓度场与实测Cl-浓度场进行了对比(图4和图5),当两者的变化趋势一致时,则认为区域尺度模型完成识别验证, 说明获得的研究区水文地质参数能够基本反映研究区的水文地质条件。

图4 研究区各长期观测孔水位模拟值与水位实测值的对比曲线Fig.4 Simulated and measured water levels in long-term observation holes in the study area

图5 识别期末研究区模拟Cl-浓度场与实测Cl-浓度场的对比Fig.5 Comparison of simulated and measured Cl- concentration fields in the study area at the end of the identification period

由图4和图5可以看出:研究区各长期观测孔的模拟水位曲线与实测水位曲线的变化趋势一致,且识别期末的模拟浓度场与实测浓度场的分布基本一致,说明所建立的区域尺度模型基本可靠,能够用于下一步局域尺度模型的建立。

2.3.2 局域尺度模型

为了提高局域尺度模型的计算精度,将研究区内3个局域海岸带模型沿海岸线1 000 m范围内的网格进行加密处理,以保证局域尺度模型能够精细刻画模型内部含水层的分布和结构特点。局域尺度模型的边界以第一类定水头和定浓度边界为主,根据区域尺度模型的流场和浓度场运用达西定律和水均衡模型计算研究区第四层、第七层和第八层3个含水层与外界交换的流量和通量值构成第二类边界(图6)。此外,针对河口海岸,由于感潮河流是海水和地下水的重要纽带,其水位受海洋潮汐水位和地下水水位共同的影响,因此选取潮汐观测站的小时数据作为模型的海洋边界取值。

图6 研究区局域尺度模型含水层参数分区图Fig.6 Hydrogeological parameters of shallow aquifer in local model in the study area

为了充分刻画不同海岸带含水层结构和渗透系数的各向异性特征,在基岩海岸带局域尺度模型中,采用离散单元来刻画断层的水力传导属性,同时采用多个有不同方向赋值的渗透系数分区来刻画基岩裂隙含水层的各向异性特征;同样,对于砂质海岸,考虑到第四系沉积物在空间上的岩性变化,可以通过调整其渗透系数参数值的大小分布来刻画孔隙含水介质的渗透性和运移性能变化。局域尺度模型的初始参数和初始条件均套用区域尺度模型,并调整局域尺度模型稳定流的地下水水均衡稳定,最终构建局域尺度模型。

收集研究区内6个长期观测孔的实测Cl-浓度数据,实测Cl-浓度监测间隔为1～3个月不等,将研究区各长期观测孔内模拟Cl-浓度与实测Cl-浓度变化进行了对比,其对比曲线如图7所示。

图7 研究区各长期观测孔内Cl-浓度模拟值与实测值的对比曲线Fig.7 Comparison curves of simulated and measured Cl- concentrations in long-term observation holes in the study area

由图7可以看出:研究区内JC002号长期观测孔位于咸淡水分界面附近,Cl-浓度值变化幅度较大;其余长期观测孔在淡水含水层,Cl-浓度值波动平缓;研究区内所有长期观测孔的模拟Cl-浓度值与实测Cl-浓度值基本吻合且两者的变化趋势一致,认为局域尺度模型完成识别,可用于下一步模型的预测。

3 海水入侵预测方案与结果分析

海水入侵预测以识别后的模型为基础,预测时间选取为2020年10月31日至2025年10月31日共计5年。将模型的参数保持不变,选择2020年10月31日流场和浓度场为初始条件;潮汐数据和地表水的入渗补给以2019—2020年的水库和河流数据为基准,以循环时间序列的方式赋值;对常用的降水量预测模型即马尔科夫链模型进行加权并采用模糊数学的方法进行改进,将深圳市1965年至2020年的降雨数据利用均值-标准差法进行枯丰频率划分,并考虑近些年的气候变化,最终采用“丰-偏丰-偏枯-平-枯”降水组合结合月降水分配系数进行赋值,见表1。

表1 研究区平均降雨量年内降水分配系数赋值表

20世纪末,由于深圳市自来水还未普及,因此该阶段地下水开采程度高,开采量较大,2006年深圳市水务部门发布了《深圳市取缔私采地下水及填埋自备井的通知》,并在2007年末开展了大规模取缔非法开采地下水的行为,全市共查封、填埋非法开采地下水井3 000余口,开采及私采、偷采地下水的现象得到了有效遏制。为了更好地研究深圳市东海岸海水入侵未来的发展趋势和不同类型海岸带的地下水水质演化规律,根据研究区地下水开采背景设定两种海水入侵预测方案,利用不同类型海岸带的局域尺度模型对研究区不同类型海岸带的海水入侵过程进行模拟预测。研究区海水入侵预测方案如下:方案一,以当前禁止开采地下水的城市政策为条件,预测5年后深圳市东海岸不同类型海岸带的海水入侵发展趋势;方案二,当深圳市东海岸供水不足时,可采取应急性地下水开采以保证城市用水。根据《深圳市水资源公报》得知,深圳市2007年地下水源供水量为3 769万m3,这其中包括1 548万m3浅层地下水和2 221万m3的深层地下水。因此,假定深圳市恢复2007年地下水开采强度,按照模型面积占比,设定深圳市东海岸区域尺度模型的地下水总开采量为678万m3/a,河口海岸、基岩海岸和砂质海岸局域尺度模型的地下水总开采量分别为28.3万m3/a、33.9万m3/a和56.5万m3/a,依据2007年浅层地下水和深层地下水的供水比例在第四层潜水含水层和第七层承压含水层中设置井群抽水,井群的位置选择在河流下游发育第四系沉积物的居民区,潜水含水层和承压含水层交替抽水。

两种预测方案下研究区区域尺度模型以及不同类型海岸带局域尺度模型Cl-浓度的预测结果,如图8所示。其中,图8(a)表示3种类型海岸带的Cl-初始浓度场;图8(b)表示采用方案一在深圳市东海岸禁止开采地下水政策条件下模型运行5年Cl-浓度的预测结果;图8(c)表示采用方案二在深圳市东海岸在678万m3/a地下水开采条件下模型运行5年Cl-浓度的预测结果。海水入侵程度由Cl-浓度划分,大多数研究中都以250 mg/L Cl-浓度值作为海水入侵的侵染指标,由于研究区内部分地区海水入侵程度不大,本文以100 mg/L Cl-浓度作为海水入侵预警指示指标,认为Cl-浓度≤100 mg/L为无入侵,Cl-浓度为100<～250 mg/L为微弱入侵,Cl-浓度为250<～500 mg/L为轻度入侵,Cl-浓度为500<～1 000 mg/L为中度入侵,Cl-浓度>1 000 mg/L为严重入侵。

图8 两种预测方案下研究区不同类型海岸带Cl-浓度分级图Fig.8 Classification of Cl- concentrations in different types of coastal zones in the study area under two prediction scenarios

由图8可以看出:

1) 根据初始阶段的海水入侵污染状况来看,河口海岸的海水入侵方式为沿河口呈舌状入侵并沿河流两侧阶地不断回退,海水最大入侵距离为1.72 km;基岩海岸的海水入侵主要沿着海岸线通过地质构造和裂隙入侵,平均入侵距离约为0.68 km;砂质海岸海水主要通过第四系沉积物沿咸淡水分界面入侵,海水平均入侵距离为0.84 km[图8(a)]。

2) 在方案一深圳东海岸禁止地下水开采条件下模型运转5年,3种类型海岸带海水入侵程度均发生了不同程度回退。其中,砂质海岸过渡带明显变窄,回退后海水入侵距离为0.59 km,分析原因为砂质海岸地势平坦,东西侧和北部地区均为山地,在地势影响下地下水水力梯度过大,地下水的渗流作用大于弥散作用,地下水流向海水导致海水入侵发生大幅度回退;基岩海岸回退后海水入侵距离为0.48 km;河口海岸海水入侵回退不大,稳定后最大海水入侵距离为1.68 km,分析原因为葵涌河上游受到东南侧罗屋田水库的补给,水库水位受人为控制较为稳定,进而使得葵涌河水位较为稳定,因此河口海岸的海水入侵不会出现较大幅度的回退[图8(b)]。

3) 在方案二深圳市东海岸以678万m3/a开采地下水条件下模型运转5年,地下水水位出现明显下降,这为海水入侵提供了水动力条件。其中,砂质海岸海水入侵加剧,海水通过第四系沉积物入侵地下水含水层,海水入侵距离达到1.22 km,海水入侵速率为0.076 km/a,海水入侵面积从3.4 km2增加到5.3 km2,海水入侵速率为0.38 km2/a;河口海岸由于河流下游与海水直接接触,从而受到海水潮汐作用的影响显著,上游河水位下降导致海水沿河流不断上溯,最大海水入侵距离达到2.05 km,海水入侵速率为0.066 km/a;海水入侵面积从3.6 km2增加到4.4 km2;基岩海岸由于远离海岸线基岩裂隙风化程度不断降低,海水入侵速度较为缓慢,但基岩海岸海水严重入侵带范围明显扩大,海水入侵距离最终稳定为0.81 km,海水入侵速率为0.026 km/a,海水入侵面积从6.3 km2增加到6.9 km2,砂质海岸的海水入侵速率约为基岩海岸的3～4倍。

4 结 论

本文通过建立深圳市东海岸区域尺度及典型海岸带局域尺度的三维变密度溶质运移模型,探究了不同类型海岸带的海水入侵机理,并讨论了多种情景下研究区地下水水质演化规律,得出了以下结论:

1) 研究区不同类型海岸带的海水入侵机理不同,河口海岸的海水沿河口呈舌状入侵,最大海水入侵距离可达1.72 km,受潮汐和感潮河流水位变化控制明显;基岩海岸平均海水入侵距离约为0.68 km,其受控于地质构造和岩石裂隙发育程度;砂质海岸的海水入侵主要呈面状入侵,平均海水入侵距离为0.84 km,海水入侵程度主要与第四系沉积物分布和局部流场相关。

2) 在保持现有的禁止开采地下水政策条件下,未来5年内研究区不同类型海岸带的海水入侵程度均会发生不同程度的回退,其中砂质海岸回退幅度最大。

3) 在深圳市东海岸以678万m3/a地下水开采条件下,未来5年内河口海岸最大海水入侵距离的增长速率为0.066 km/a,砂质海岸海水入侵速率为0.076 km/a,基岩海岸海水入侵速率约为0.026 km/a,且砂质海岸的海水入侵速率约为基岩海岸的3～4倍。