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

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

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

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

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

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

研究區(qū)位于廣東省深圳市東部沿海區(qū)(圖1),北與惠州市、坪山區(qū)相接,以山前為邊界,西至沙頭角河與深圳河,東部和南部靠海,與大鵬灣和大亞灣相望,面積約為361 km2。區(qū)內(nèi)地勢(shì)總體北高南低,發(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,水量相對(duì)貧乏。區(qū)內(nèi)第四系潛水含水層分布不規(guī)律且不連續(xù),而其下的基巖裂隙含水層分布廣泛且具有統(tǒng)一性,其隔水頂板多在海平面之下。區(qū)內(nèi)地下水補(bǔ)給以大氣降水、河流側(cè)滲和上游水庫(kù)遠(yuǎn)源補(bǔ)給為主,地下水徑流方向主要受地形控制,總體由北至南徑流,最終向海排泄。

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

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

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

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

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

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

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

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

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

模型Ⅰ:

模型Ⅱ:

式中:Kij為含水層滲透系數(shù)(m/d);xi、xj為坐標(biāo)系(i、j=1,2,3);C為液體濃度(kg/m3);C0為液體初始濃度(kg/m3);C1為源(匯)項(xiàng)濃度(kg/m3);H為求解空間的位置水頭(m);t為時(shí)間(d);H0為初始水頭(m);HB為邊界Γ1給定水頭(m);ej為重力方向單位矢量第j個(gè)分量;φ為孔隙度;η為密度耦合系數(shù);μs為比儲(chǔ)水系數(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為潛水面變動(dòng)帶地下水密度(kg/m3);H1為潛水面Γ2-1上的參考水頭(m);x3為潛水面Γ2-1上的水頭高程;W為潛水面單位面積上的降水入滲補(bǔ)給量(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ī)律,先對(duì)研究區(qū)整體區(qū)域進(jìn)行數(shù)值建模,在完成區(qū)域水流-水質(zhì)模型的識(shí)別與驗(yàn)證后,采取參數(shù)套用的方式初步建立研究區(qū)不同類型海岸帶的局域尺度模型;然后根據(jù)區(qū)域尺度模型模擬過(guò)程中的流場(chǎng)和濃度場(chǎng)計(jì)算出局域尺度模型的邊界條件,并根據(jù)不同類型海岸帶的特點(diǎn)對(duì)局域尺度模型海岸帶的參數(shù)進(jìn)行精細(xì)化處理,并調(diào)整穩(wěn)定流模型使其水均衡穩(wěn)定,最終完成研究區(qū)不同類型海岸帶局域尺度模型的建立。

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

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

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

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

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

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

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

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

圖5 識(shí)別期末研究區(qū)模擬Cl-濃度場(chǎng)與實(shí)測(cè)Cl-濃度場(chǎng)的對(duì)比Fig.5 Comparison of simulated and measured Cl- concentration fields in the study area at the end of the identification period

由圖4和圖5可以看出:研究區(qū)各長(zhǎng)期觀測(cè)孔的模擬水位曲線與實(shí)測(cè)水位曲線的變化趨勢(shì)一致,且識(shí)別期末的模擬濃度場(chǎng)與實(shí)測(cè)濃度場(chǎng)的分布基本一致,說(shuō)明所建立的區(qū)域尺度模型基本可靠,能夠用于下一步局域尺度模型的建立。

2.3.2 局域尺度模型

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

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

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

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

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

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

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

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

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

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

兩種預(yù)測(cè)方案下研究區(qū)區(qū)域尺度模型以及不同類型海岸帶局域尺度模型Cl-濃度的預(yù)測(cè)結(jié)果,如圖8所示。其中,圖8(a)表示3種類型海岸帶的Cl-初始濃度場(chǎng);圖8(b)表示采用方案一在深圳市東海岸禁止開(kāi)采地下水政策條件下模型運(yùn)行5年Cl-濃度的預(yù)測(cè)結(jié)果;圖8(c)表示采用方案二在深圳市東海岸在678萬(wàn)m3/a地下水開(kāi)采條件下模型運(yùn)行5年Cl-濃度的預(yù)測(cè)結(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為無(wú)入侵,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ù)測(cè)方案下研究區(qū)不同類型海岸帶Cl-濃度分級(jí)圖Fig.8 Classification of Cl- concentrations in different types of coastal zones in the study area under two prediction scenarios

由圖8可以看出:

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

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

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

4 結(jié) 論

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

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

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

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