考慮不同海岸類型的深圳市東海岸海水入侵多尺度耦合模擬

宋紹宇,成建梅*,程天舜,盧薇艷,職敬儒

(1.中國地質(zhì)大學(xué)(武漢) 環(huán)境學(xué)院,湖北 武漢 430078;2.深圳市水務(wù)規(guī)劃設(shè)計院股份有限公司,廣東 深圳 518023;3.深圳市地質(zhì)局,廣東 深圳 518023)

全世界許多國家和地區(qū)都出現(xiàn)了海水入侵問題,對海水入侵的研究也引起了眾多學(xué)者的關(guān)注[1-4]。數(shù)值模擬技術(shù)是定量研究海水入侵發(fā)展演化的重要手段,常見的海水入侵?jǐn)?shù)值模擬模型包括突變界面模型和變密度過渡帶模型[5-7]。目前常采用變密度過渡帶模型來模擬并解決實際海水入侵問題,主要包括基于單個組分或總礦化度的考慮變密度海水入侵三維溶質(zhì)運移模型[8-9]、考慮海咸水運移過程中多個水化學(xué)組分的三維水流-水質(zhì)模型[10-11]和考慮耦合水化學(xué)反應(yīng)過程的高濃度海水入侵三維水質(zhì)模型[12]。濱海含水系統(tǒng)非均質(zhì)結(jié)構(gòu)與海水入侵優(yōu)勢通道的刻畫、不同類型臨海邊界的處理、多尺度的耦合是實現(xiàn)高精度海水入侵過程模擬的關(guān)鍵[4,13-15]。面對不同的研究目標(biāo),可選擇的海水入侵?jǐn)?shù)值模擬模型類型不同。吳吉春等[16]提出要以不同類型海岸帶的海水入侵發(fā)生、演化、防治為主線,探究不同類型海岸帶海水入侵發(fā)生的機理與作用機制,并探求高效的海水入侵問題管理對策。此外,針對濱海含水層的復(fù)雜性和海岸類型的差別,一些學(xué)者對單一類型海岸帶(如砂質(zhì)海岸、裂隙基巖海岸或灰?guī)r海岸等)海水入侵特征開展了數(shù)值模擬研究[1,3,17-20],但極少的模擬案例能涉及多種海岸類型[21],無法精確表征實際復(fù)雜的海岸帶結(jié)構(gòu)特點和海水入侵演化規(guī)律。

20世紀(jì)80年代以來,隨著深圳市的改革開放和城市迅速發(fā)展,地下水過量開采導(dǎo)致深圳海岸帶發(fā)生了嚴(yán)重的海水入侵,加之受城鎮(zhèn)化建設(shè)、填海造陸、海水養(yǎng)殖等一系列海岸帶人類工程活動的影響,使得該地區(qū)海水入侵的發(fā)展演化機制更為復(fù)雜。目前已有的研究多集中在利用水化學(xué)、同位素及電阻率等指標(biāo)對深圳海岸帶海水入侵過程進行識別和評估,涉及數(shù)值模擬模型的研究極少[22-23],且研究區(qū)僅限于深圳西海岸。而深圳市東部出露有河口海岸、基巖海岸、砂質(zhì)海岸等多種海岸帶類型,海水入侵方式更為復(fù)雜多樣。為此,本文以深圳市東海岸為研究區(qū),充分發(fā)揮區(qū)域滲流對局域滲流的總體控制優(yōu)勢,采用多尺度耦合模擬方法,建立了研究區(qū)區(qū)域尺度及典型海岸帶局域尺度的三維變密度溶質(zhì)運移模型,并通過數(shù)值模擬方法探究了研究區(qū)不同類型海岸帶的海水入侵機理與地下水水質(zhì)演化規(guī)律。

1 研究區(qū)背景條件

研究區(qū)位于廣東省深圳市東部沿海區(qū)(圖1),北與惠州市、坪山區(qū)相接,以山前為邊界,西至沙頭角河與深圳河,東部和南部靠海,與大鵬灣和大亞灣相望,面積約為361 km2。區(qū)內(nèi)地勢總體北高南低,發(fā)育有山地和丘陵,南側(cè)有海積平原。區(qū)內(nèi)較大的河流有鹽田河、葵涌河、王母河和東涌河等,河流兩側(cè)及濱海地帶的海積平原發(fā)育有第四系松散地層,巖性以中細(xì)砂、粗砂及卵石等為主,賦存豐富的孔隙水,單井出水量為100～150 m3/d。區(qū)內(nèi)發(fā)育泥盆系中統(tǒng)鼎壺山群砂巖、白堊系-侏羅系花崗巖、侏羅系上統(tǒng)高基坪群酸性火山巖、石炭系下統(tǒng)石磴子組灰?guī)r等基巖地層,賦存風(fēng)化裂隙水和構(gòu)造裂隙水,單井出水量為50 m3/d,水量相對貧乏。區(qū)內(nèi)第四系潛水含水層分布不規(guī)律且不連續(xù),而其下的基巖裂隙含水層分布廣泛且具有統(tǒng)一性,其隔水頂板多在海平面之下。區(qū)內(nèi)地下水補給以大氣降水、河流側(cè)滲和上游水庫遠源補給為主,地下水徑流方向主要受地形控制,總體由北至南徑流,最終向海排泄。

圖1 研究區(qū)含水巖組結(jié)構(gòu)簡圖Fig.1 Structure diagram of water-bearing rock group of the study area

深圳市東部海岸線長約123 km,海岸帶類型復(fù)雜。其中,砂質(zhì)海岸主要分布在東涌、西涌地區(qū)以及大鵬澳沿岸,地形平緩,地貌上屬于沖洪積平原及海積平原,如圖1中c-d、e-f剖面所示,巖性以第四系砂礫石層為主,透水性較強;基巖海岸主要分布在東部大鵬灣和大亞灣西側(cè)的大鵬半島沿岸,長度約為81.9 km,地形上以山地、丘陵為主,基巖裂隙相對發(fā)育;河口海岸依托研究區(qū)內(nèi)河流發(fā)育,如圖1中a-b剖面所示,葵涌河河口與海洋相接,且入海口花崗巖基巖出露,葵涌河成為內(nèi)陸地下水與海水的重要交換媒介。深圳市東部海岸帶的海水入侵咸淡水界面如圖1所示,海水在不同類型海岸段入侵程度不同,其中基巖海岸的海水入侵距離以300 m到1 000 m為主,海水入侵主要受到地形、巖性和裂隙構(gòu)造發(fā)育等因素影響;而砂質(zhì)海岸和河口海岸的海水入侵范圍較大,海水入侵距離受第四系分布和感潮河流水位變幅的影響較大。

2 研究區(qū)海水入侵多尺度耦合模型構(gòu)建

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

根據(jù)深圳市東海岸涉及的地層巖性特征及透水性強弱,將表層覆蓋的人工填土、海陸交互成因的淤泥和黏土以及全風(fēng)化的基巖概化為弱透水層,將第四系松散沉積物砂、卵石以及強風(fēng)化、中風(fēng)化基巖概化為含水層,將研究區(qū)地層概化為9層42區(qū),體現(xiàn)為“三弱三含”的地層結(jié)構(gòu)特征。在區(qū)域模型的邊界概化上,研究區(qū)向西延伸至梧桐山,以相鄰深圳河、沙頭角河以及小型的地下分水嶺為邊界;研究區(qū)東部和南部靠海,在模型中處理為定水頭邊界和定濃度邊界;研究區(qū)北部地勢較高,巖性出露為侏羅系下統(tǒng)橋源組砂巖、石英砂巖和花崗巖,此區(qū)域巖石的滲透系數(shù)較小,選取地下分水嶺為模型北部的隔水邊界;模型上邊界為潛水面,接受大氣降水的補給,模型底部為侏羅系-白堊系微風(fēng)化花崗巖,由于其風(fēng)化程度較低,滲透性較差,處理為零通量邊界。

滲透性強的砂礫石層和裂隙發(fā)育的基巖破碎帶是海岸帶海水入侵內(nèi)陸含水層的重要通道,從而加劇海水的入侵[24]。海岸帶依據(jù)海岸帶底質(zhì)和空間形態(tài)的不同劃分為基巖海岸、砂質(zhì)海岸、淤泥質(zhì)海岸、生物海岸和河口海岸[25]。因此,為了更好地探究研究區(qū)不同類型海岸帶的海水入侵發(fā)展趨勢,在區(qū)域尺度模型的基礎(chǔ)上選擇砂質(zhì)海岸、基巖海岸和河口海岸3種不同類型海岸帶建立局域尺度模型,如圖2所示。

在局域尺度模型的邊界概化上,基巖海岸帶模型范圍選擇研究區(qū)大鵬半島的東北部基巖海岸帶,區(qū)內(nèi)地勢較高,主要出露有泥盆系雙頭群和鼎湖山群的石英砂巖和粉砂巖,構(gòu)造裂隙發(fā)育,地下水類型為層狀基巖裂隙水;選擇研究區(qū)葵涌河地區(qū)建立河口海岸帶模型,區(qū)內(nèi)主要出露有石磴子組結(jié)晶灰?guī)r、白云質(zhì)灰?guī)r及少量壺天群大理巖,發(fā)育數(shù)量較多、規(guī)模較大的巖溶裂隙和溶洞,成為地下水理想的儲運空間,且沿河流兩側(cè)發(fā)育有第四系沉積物,在入海口處有花崗巖出露,隨著潮汐波動,葵涌河成為感潮河流,海水在大潮時期入侵內(nèi)河,并侵染灰?guī)r含水層;選擇研究區(qū)東涌和西涌地區(qū)的海積平原區(qū)構(gòu)建砂質(zhì)海岸帶模型,區(qū)內(nèi)主要分布有第四系海相沉積物、鼎湖山群石英砂巖和白堊系花崗巖,地下水類型為第四系孔隙水和基巖裂隙水。

2.2 數(shù)學(xué)模型建立

過渡帶海水入侵的數(shù)學(xué)模型由兩個偏微分方程及其定解條件來描述,其中一個偏微分方程用于描述受密度變化影響的混合液體(咸淡水混合物)的流動,另一個偏微分方程用于描述混合液體中鹽分的運移[26]。考慮到研究區(qū)不同海岸帶含水層分布的特點,這里采用考慮含水層各向異性的變密度三維水流模型(模型I)和變密度三維對流-彌散模型(模型Ⅱ)來描述海水入侵過程,模型I和模型Ⅱ分別表示如下。

模型Ⅰ:

模型Ⅱ:

式中:Kij為含水層滲透系數(shù)(m/d);xi、xj為坐標(biāo)系(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個分量;φ為孔隙度;η為密度耦合系數(shù);μs為比儲水系數(shù);μd為重力給水度;ρ0為淡水密度(kg/m3);ρ為混合液體密度(kg/m3);vi為地下水滲透速度在xi上的分量(m/d);q為單位體積多孔介質(zhì)中抽取或注入的水量(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為彌散系數(shù)張量;CB為邊界Γ1上的濃度(kg/m3);C2為降水或灌溉用水等入滲水的濃度(kg/m3);C3為定通量邊界流入或流出水的濃度(kg/m3)。

2.3 數(shù)值建模

為了探究研究區(qū)不同類型海岸帶的海水入侵演化規(guī)律,先對研究區(qū)整體區(qū)域進行數(shù)值建模,在完成區(qū)域水流-水質(zhì)模型的識別與驗證后,采取參數(shù)套用的方式初步建立研究區(qū)不同類型海岸帶的局域尺度模型;然后根據(jù)區(qū)域尺度模型模擬過程中的流場和濃度場計算出局域尺度模型的邊界條件,并根據(jù)不同類型海岸帶的特點對局域尺度模型海岸帶的參數(shù)進行精細(xì)化處理,并調(diào)整穩(wěn)定流模型使其水均衡穩(wěn)定,最終完成研究區(qū)不同類型海岸帶局域尺度模型的建立。

2.3.1 區(qū)域尺度模型建立

在濱海含水系統(tǒng)的模型概化中,臨海邊界條件的正確處理直接影響到海水入侵模型的仿真程度。目前絕大多數(shù)海水入侵模型將平面圖的海岸線作為含水系統(tǒng)的邊界,但實際上海岸近海地形一般為緩坡向外海傾斜,延伸的距離至少有幾十米甚至超過百米。在本模型中,通過查閱前人資料,并綜合潮汐數(shù)據(jù)選定距海岸線150 m處作為模型向海一側(cè)外邊界,且處理為定水頭邊界。

依據(jù)數(shù)學(xué)模型,利用Feflow軟件建立數(shù)值模型并求解。使用Triangle算法對研究區(qū)進行三角網(wǎng)格剖分,對研究區(qū)重點河流、斷層、長期觀測孔進行剖分單元加密,并且保證觀測孔位置在結(jié)點上。根據(jù)概念模型中對整個研究區(qū)的垂向分層及橫向分區(qū)(圖3),考慮含水層的滲透系數(shù)K、重力給水度、比儲水系數(shù)、縱/橫向彌散度、擴散系數(shù)等水文地質(zhì)參數(shù),用來刻畫各向異性含水層的變密度三維溶質(zhì)運移過程。

圖3 研究區(qū)部分含水層參數(shù)分區(qū)圖Fig.3 Hydrogeological parameter zones in some aquifers in the study area

在模型的源匯項處理方面,根據(jù)土地利用類型將研究區(qū)劃分為填海區(qū)、建筑區(qū)、植被區(qū)、水域、灌木和耕地6種類型的降雨入滲分區(qū)[圖3(a)],并收集研究區(qū)19個雨量站的日降雨量數(shù)據(jù),采用泰森多邊形與降雨入滲補給系數(shù)疊加的方式賦于模型;將模型中河流和水庫概化為柯西邊界(第三類邊界Fluid-transfer邊界),根據(jù)深圳市海水入侵地質(zhì)災(zāi)害調(diào)查資料[27],設(shè)定地表水與地下水水力交換注入傳輸率Kin為0.05/d,抽取傳輸率Kout為0.10/d,并收集研究區(qū)內(nèi)羅屋田、鐵扇關(guān)門、楓木浪等水庫的日水位變化數(shù)據(jù),通過時間序列的方式賦于模型,同時收集研究區(qū)內(nèi)鹽田河、葵涌河、王母河等河流河口及中上游的日水位數(shù)據(jù),通過反距離插值法且參考河流地表高程生成河流各結(jié)點水位。

潮汐作用會使得咸淡水交換過程更為頻繁,在漲潮時海水由高潮線附近的海灘滲入潛水含水層,落潮時由低潮線附近的海灘出流,這種潮汐水位的波動引起的海水-地下水循環(huán)模式會對海水入侵過程產(chǎn)生重要的影響。本文選取大鵬灣區(qū)域潮汐觀測站每日6:00、12:00和18:00的潮汐水位和大亞灣區(qū)域潮汐觀測站每日12:00和18:00的潮汐水位以線性時間序列的方式作為模型的海洋邊界。模型的初始流場和初始濃度場通過Surfer軟件對收集到的實際數(shù)據(jù)進行克里金插值,并與穩(wěn)定流模型結(jié)果進行對比后確定出更接近實際情況的初始流場和初始濃度場。

模型的識別與驗證是地下水?dāng)?shù)值模擬中極為重要的工作,為了驗證區(qū)域尺度模型的參數(shù)設(shè)置是否合理,本次選取的模型識別期為2019年10月31日至2020年10月31日,收集研究區(qū)內(nèi)6個長期觀測孔的實測水位數(shù)據(jù),觀測間隔為5 d,并收集2020年10月下旬的114個觀測點的實測Cl-濃度值進行內(nèi)插,得到識別期末區(qū)域?qū)崪yCl-濃度場。在模型的識別過程中,對研究區(qū)的水文地質(zhì)參數(shù)、邊界條件以及源匯項等進行反復(fù)調(diào)整,并將研究區(qū)長期觀測孔的模擬水位與實測水位變化進行了比較,同時將識別期末研究區(qū)模擬Cl-濃度場與實測Cl-濃度場進行了對比(圖4和圖5),當(dāng)兩者的變化趨勢一致時,則認(rèn)為區(qū)域尺度模型完成識別驗證, 說明獲得的研究區(qū)水文地質(zhì)參數(shù)能夠基本反映研究區(qū)的水文地質(zhì)條件。

圖4 研究區(qū)各長期觀測孔水位模擬值與水位實測值的對比曲線Fig.4 Simulated and measured water levels in long-term observation holes in the study area

圖5 識別期末研究區(qū)模擬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可以看出:研究區(qū)各長期觀測孔的模擬水位曲線與實測水位曲線的變化趨勢一致,且識別期末的模擬濃度場與實測濃度場的分布基本一致,說明所建立的區(qū)域尺度模型基本可靠,能夠用于下一步局域尺度模型的建立。

2.3.2 局域尺度模型

為了提高局域尺度模型的計算精度,將研究區(qū)內(nèi)3個局域海岸帶模型沿海岸線1 000 m范圍內(nèi)的網(wǎng)格進行加密處理,以保證局域尺度模型能夠精細(xì)刻畫模型內(nèi)部含水層的分布和結(jié)構(gòu)特點。局域尺度模型的邊界以第一類定水頭和定濃度邊界為主,根據(jù)區(qū)域尺度模型的流場和濃度場運用達西定律和水均衡模型計算研究區(qū)第四層、第七層和第八層3個含水層與外界交換的流量和通量值構(gòu)成第二類邊界(圖6)。此外,針對河口海岸,由于感潮河流是海水和地下水的重要紐帶,其水位受海洋潮汐水位和地下水水位共同的影響,因此選取潮汐觀測站的小時數(shù)據(jù)作為模型的海洋邊界取值。

圖6 研究區(qū)局域尺度模型含水層參數(shù)分區(qū)圖Fig.6 Hydrogeological parameters of shallow aquifer in local model in the study area

為了充分刻畫不同海岸帶含水層結(jié)構(gòu)和滲透系數(shù)的各向異性特征,在基巖海岸帶局域尺度模型中,采用離散單元來刻畫斷層的水力傳導(dǎo)屬性,同時采用多個有不同方向賦值的滲透系數(shù)分區(qū)來刻畫基巖裂隙含水層的各向異性特征;同樣,對于砂質(zhì)海岸,考慮到第四系沉積物在空間上的巖性變化,可以通過調(diào)整其滲透系數(shù)參數(shù)值的大小分布來刻畫孔隙含水介質(zhì)的滲透性和運移性能變化。局域尺度模型的初始參數(shù)和初始條件均套用區(qū)域尺度模型,并調(diào)整局域尺度模型穩(wěn)定流的地下水水均衡穩(wěn)定,最終構(gòu)建局域尺度模型。

收集研究區(qū)內(nèi)6個長期觀測孔的實測Cl-濃度數(shù)據(jù),實測Cl-濃度監(jiān)測間隔為1～3個月不等,將研究區(qū)各長期觀測孔內(nèi)模擬Cl-濃度與實測Cl-濃度變化進行了對比,其對比曲線如圖7所示。

圖7 研究區(qū)各長期觀測孔內(nèi)Cl-濃度模擬值與實測值的對比曲線Fig.7 Comparison curves of simulated and measured Cl- concentrations in long-term observation holes in the study area

由圖7可以看出:研究區(qū)內(nèi)JC002號長期觀測孔位于咸淡水分界面附近,Cl-濃度值變化幅度較大;其余長期觀測孔在淡水含水層,Cl-濃度值波動平緩;研究區(qū)內(nèi)所有長期觀測孔的模擬Cl-濃度值與實測Cl-濃度值基本吻合且兩者的變化趨勢一致,認(rèn)為局域尺度模型完成識別,可用于下一步模型的預(yù)測。

3 海水入侵預(yù)測方案與結(jié)果分析

海水入侵預(yù)測以識別后的模型為基礎(chǔ),預(yù)測時間選取為2020年10月31日至2025年10月31日共計5年。將模型的參數(shù)保持不變,選擇2020年10月31日流場和濃度場為初始條件;潮汐數(shù)據(jù)和地表水的入滲補給以2019—2020年的水庫和河流數(shù)據(jù)為基準(zhǔn),以循環(huán)時間序列的方式賦值;對常用的降水量預(yù)測模型即馬爾科夫鏈模型進行加權(quán)并采用模糊數(shù)學(xué)的方法進行改進,將深圳市1965年至2020年的降雨數(shù)據(jù)利用均值-標(biāo)準(zhǔn)差法進行枯豐頻率劃分,并考慮近些年的氣候變化,最終采用“豐-偏豐-偏枯-平-枯”降水組合結(jié)合月降水分配系數(shù)進行賦值,見表1。

表1 研究區(qū)平均降雨量年內(nèi)降水分配系數(shù)賦值表

20世紀(jì)末,由于深圳市自來水還未普及,因此該階段地下水開采程度高,開采量較大,2006年深圳市水務(wù)部門發(fā)布了《深圳市取締私采地下水及填埋自備井的通知》,并在2007年末開展了大規(guī)模取締非法開采地下水的行為,全市共查封、填埋非法開采地下水井3 000余口,開采及私采、偷采地下水的現(xiàn)象得到了有效遏制。為了更好地研究深圳市東海岸海水入侵未來的發(fā)展趨勢和不同類型海岸帶的地下水水質(zhì)演化規(guī)律,根據(jù)研究區(qū)地下水開采背景設(shè)定兩種海水入侵預(yù)測方案,利用不同類型海岸帶的局域尺度模型對研究區(qū)不同類型海岸帶的海水入侵過程進行模擬預(yù)測。研究區(qū)海水入侵預(yù)測方案如下:方案一,以當(dāng)前禁止開采地下水的城市政策為條件,預(yù)測5年后深圳市東海岸不同類型海岸帶的海水入侵發(fā)展趨勢;方案二,當(dāng)深圳市東海岸供水不足時,可采取應(yīng)急性地下水開采以保證城市用水。根據(jù)《深圳市水資源公報》得知,深圳市2007年地下水源供水量為3 769萬m3,這其中包括1 548萬m3淺層地下水和2 221萬m3的深層地下水。因此,假定深圳市恢復(fù)2007年地下水開采強度,按照模型面積占比,設(shè)定深圳市東海岸區(qū)域尺度模型的地下水總開采量為678萬m3/a,河口海岸、基巖海岸和砂質(zhì)海岸局域尺度模型的地下水總開采量分別為28.3萬m3/a、33.9萬m3/a和56.5萬m3/a,依據(jù)2007年淺層地下水和深層地下水的供水比例在第四層潛水含水層和第七層承壓含水層中設(shè)置井群抽水,井群的位置選擇在河流下游發(fā)育第四系沉積物的居民區(qū),潛水含水層和承壓含水層交替抽水。

兩種預(yù)測方案下研究區(qū)區(qū)域尺度模型以及不同類型海岸帶局域尺度模型Cl-濃度的預(yù)測結(jié)果,如圖8所示。其中,圖8(a)表示3種類型海岸帶的Cl-初始濃度場;圖8(b)表示采用方案一在深圳市東海岸禁止開采地下水政策條件下模型運行5年Cl-濃度的預(yù)測結(jié)果;圖8(c)表示采用方案二在深圳市東海岸在678萬m3/a地下水開采條件下模型運行5年Cl-濃度的預(yù)測結(jié)果。海水入侵程度由Cl-濃度劃分,大多數(shù)研究中都以250 mg/L Cl-濃度值作為海水入侵的侵染指標(biāo),由于研究區(qū)內(nèi)部分地區(qū)海水入侵程度不大,本文以100 mg/L Cl-濃度作為海水入侵預(yù)警指示指標(biāo),認(rèn)為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為嚴(yán)重入侵。

圖8 兩種預(yù)測方案下研究區(qū)不同類型海岸帶Cl-濃度分級圖Fig.8 Classification of Cl- concentrations in different types of coastal zones in the study area under two prediction scenarios

由圖8可以看出:

1) 根據(jù)初始階段的海水入侵污染狀況來看,河口海岸的海水入侵方式為沿河口呈舌狀入侵并沿河流兩側(cè)階地不斷回退,海水最大入侵距離為1.72 km;基巖海岸的海水入侵主要沿著海岸線通過地質(zhì)構(gòu)造和裂隙入侵,平均入侵距離約為0.68 km;砂質(zhì)海岸海水主要通過第四系沉積物沿咸淡水分界面入侵,海水平均入侵距離為0.84 km[圖8(a)]。

2) 在方案一深圳東海岸禁止地下水開采條件下模型運轉(zhuǎn)5年,3種類型海岸帶海水入侵程度均發(fā)生了不同程度回退。其中,砂質(zhì)海岸過渡帶明顯變窄,回退后海水入侵距離為0.59 km,分析原因為砂質(zhì)海岸地勢平坦,東西側(cè)和北部地區(qū)均為山地,在地勢影響下地下水水力梯度過大,地下水的滲流作用大于彌散作用,地下水流向海水導(dǎo)致海水入侵發(fā)生大幅度回退;基巖海岸回退后海水入侵距離為0.48 km;河口海岸海水入侵回退不大,穩(wěn)定后最大海水入侵距離為1.68 km,分析原因為葵涌河上游受到東南側(cè)羅屋田水庫的補給,水庫水位受人為控制較為穩(wěn)定,進而使得葵涌河水位較為穩(wěn)定,因此河口海岸的海水入侵不會出現(xiàn)較大幅度的回退[圖8(b)]。

3) 在方案二深圳市東海岸以678萬m3/a開采地下水條件下模型運轉(zhuǎn)5年,地下水水位出現(xiàn)明顯下降,這為海水入侵提供了水動力條件。其中,砂質(zhì)海岸海水入侵加劇,海水通過第四系沉積物入侵地下水含水層,海水入侵距離達到1.22 km,海水入侵速率為0.076 km/a,海水入侵面積從3.4 km2增加到5.3 km2,海水入侵速率為0.38 km2/a;河口海岸由于河流下游與海水直接接觸,從而受到海水潮汐作用的影響顯著,上游河水位下降導(dǎo)致海水沿河流不斷上溯,最大海水入侵距離達到2.05 km,海水入侵速率為0.066 km/a;海水入侵面積從3.6 km2增加到4.4 km2;基巖海岸由于遠離海岸線基巖裂隙風(fēng)化程度不斷降低,海水入侵速度較為緩慢,但基巖海岸海水嚴(yán)重入侵帶范圍明顯擴大,海水入侵距離最終穩(wěn)定為0.81 km,海水入侵速率為0.026 km/a,海水入侵面積從6.3 km2增加到6.9 km2,砂質(zhì)海岸的海水入侵速率約為基巖海岸的3～4倍。

4 結(jié) 論

本文通過建立深圳市東海岸區(qū)域尺度及典型海岸帶局域尺度的三維變密度溶質(zhì)運移模型,探究了不同類型海岸帶的海水入侵機理,并討論了多種情景下研究區(qū)地下水水質(zhì)演化規(guī)律,得出了以下結(jié)論:

1) 研究區(qū)不同類型海岸帶的海水入侵機理不同,河口海岸的海水沿河口呈舌狀入侵,最大海水入侵距離可達1.72 km,受潮汐和感潮河流水位變化控制明顯;基巖海岸平均海水入侵距離約為0.68 km,其受控于地質(zhì)構(gòu)造和巖石裂隙發(fā)育程度;砂質(zhì)海岸的海水入侵主要呈面狀入侵,平均海水入侵距離為0.84 km,海水入侵程度主要與第四系沉積物分布和局部流場相關(guān)。

2) 在保持現(xiàn)有的禁止開采地下水政策條件下,未來5年內(nèi)研究區(qū)不同類型海岸帶的海水入侵程度均會發(fā)生不同程度的回退,其中砂質(zhì)海岸回退幅度最大。

3) 在深圳市東海岸以678萬m3/a地下水開采條件下,未來5年內(nèi)河口海岸最大海水入侵距離的增長速率為0.066 km/a,砂質(zhì)海岸海水入侵速率為0.076 km/a,基巖海岸海水入侵速率約為0.026 km/a,且砂質(zhì)海岸的海水入侵速率約為基巖海岸的3～4倍。