賈超,朱恒華,虞未江,衛(wèi)政潤(rùn),劉治政,李雙
(1.山東大學(xué)海洋研究院,山東 青島 266237;2.山東省地質(zhì)調(diào)查院,山東 濟(jì)南 250014)
隨著經(jīng)濟(jì)社會(huì)的發(fā)展,人類生活和工業(yè)生產(chǎn)需水量激增,許多城市都開采地下水來滿足用水量,尤其沿海城市因?yàn)檫^度開采地下水,導(dǎo)致內(nèi)陸地下水位大幅度下降,破壞了海咸水與內(nèi)陸地下淡水之間的動(dòng)態(tài)平衡,從而引發(fā)海岸帶咸淡水界面向內(nèi)陸遷移,即海水入侵[1-2]。美國(guó)、荷蘭、西班牙等國(guó)家的濱海城市均出現(xiàn)海水入侵問題,中國(guó)遼東灣、山東半島、上海市等地區(qū)也發(fā)生海水入侵,其中山東萊州灣地區(qū)是海水入侵較為嚴(yán)重的地區(qū)之一,海咸水入侵將污染內(nèi)陸地下淡水資源,破壞地下水環(huán)境,造成沿海地區(qū)土地固結(jié)化和鹽漬化,嚴(yán)重阻礙當(dāng)?shù)亟?jīng)濟(jì)社會(huì)發(fā)展。
海水入侵現(xiàn)象在全球范圍內(nèi)普遍存在[3-9],國(guó)內(nèi)外眾多學(xué)者對(duì)海水入侵問題進(jìn)行深入研究,在海水入侵模型、基本理論、咸淡水界面形狀和遷移機(jī)理等方面取得許多重要成果。目前,主要有兩類模型對(duì)海水入侵進(jìn)行刻畫,即咸淡水界面突變模型和海陸過渡帶模型。海陸過渡帶模型因其在刻畫復(fù)雜水文地質(zhì)條件和人類活動(dòng)條件等方面具有優(yōu)勢(shì),模型的仿真性和可靠性較高,而被國(guó)內(nèi)外普遍采用[10-15]。
國(guó)內(nèi)對(duì)海水入侵已有較多研究,該文在相關(guān)研究的基礎(chǔ)上[16-18],選擇海陸過渡帶模型,結(jié)合地下水運(yùn)動(dòng)理論,選取受海水入侵危害較嚴(yán)重的山東省龍口、招遠(yuǎn)濱海區(qū)域作為研究區(qū),研究該地區(qū)海陸過渡帶遷移的規(guī)律。在對(duì)該區(qū)域水文地質(zhì)基礎(chǔ)資料進(jìn)行研究分析的基礎(chǔ)上,構(gòu)建研究區(qū)水文地質(zhì)概念模型,建立地下水流場(chǎng)數(shù)值模型,并對(duì)數(shù)值模型進(jìn)行識(shí)別和檢驗(yàn)。選擇作為模擬因子,在地下水流場(chǎng)的基礎(chǔ)上,對(duì)研究區(qū)海陸過渡帶的遷移進(jìn)行數(shù)值模擬并研究?jī)?nèi)陸地下水位對(duì)海陸過渡帶遷移的影響,總結(jié)分析研究區(qū)海陸過渡帶遷移的規(guī)律,并為該地區(qū)合理開發(fā)利用沿海地下淡水資源和防控海咸水入侵提供切實(shí)有效的建議。
沿海地區(qū)地下含水層中的水力梯度基本上是傾向海洋,內(nèi)陸地下含水層中的大部分淡水最終流向海洋。因?yàn)楹F矫嬉韵碌牡叵潞畬又写嬖谟泻O趟?,所以流向海洋的淡水與海咸水之間存在一條接觸帶。此外,海咸水與淡水是兩種可以互溶的流體,在流體動(dòng)力彌散作用下,該接觸帶將以過渡帶形式存在。海陸過渡帶模型及咸淡水流動(dòng)如圖1所示。
圖1 地下含水層中海咸水和淡水流動(dòng)過程示意圖
海陸過渡帶模型由兩個(gè)數(shù)學(xué)偏微分方程描述,一個(gè)是描述密度不斷改變的混合液體流動(dòng)方程,另一個(gè)是描述海水中鹽分運(yùn)移的對(duì)流—彌散方程,通過兩個(gè)數(shù)學(xué)方程,把流體的密度、濃度和水位聯(lián)系在一起,得到海陸過渡帶的分布范圍、濃度值和水位值。
(1)描述密度不斷改變的混合液體流動(dòng)方程
(1)
該式是在假設(shè)密度隨濃度呈線性變化和忽略液體壓強(qiáng)對(duì)密度影響的基礎(chǔ)上推導(dǎo)出來[12]。
(2)描述海水中鹽分運(yùn)移的對(duì)流—彌散方程
(2)
式中Dij—彌散系數(shù)張量,(i,j)=(1,2,3);u=(u1,u2,u3)—地下水實(shí)際平均流速矢量;c*—注入或者抽出液體中相應(yīng)組分的濃度。
(1)地下水運(yùn)動(dòng)微分方程:在達(dá)西定律和質(zhì)量守恒定律基礎(chǔ)上,推導(dǎo)出非均質(zhì)各向異性含水層中承壓水運(yùn)動(dòng)的基本微分方程:
(3)
(2)溶質(zhì)運(yùn)移模型的控制方程:描述模擬因子k在三維地下水流系統(tǒng)中的遷移和去向:
(4)
研究區(qū)龍口、招遠(yuǎn)濱海區(qū)域位于山東省膠東半島西北部,渤海灣南岸,西與萊州灣相接,東北與天津、大連等城市以及朝鮮半島隔海相望,地理坐標(biāo)為東經(jīng)120°13′38″~120°34′12″,北緯37°45′16″~37°29′40″之間,研究區(qū)東西最大橫距約為23km,南北最大縱距約為25km,海岸線長(zhǎng)約為55km,研究區(qū)總面積約為424km2,地理位置如圖2所示。
圖2 研究區(qū)地理位置示意圖
研究區(qū)屬暖溫帶半濕潤(rùn)季風(fēng)型大陸性氣候,年平均氣溫為11.6℃,多年平均降水量為592.7mm,主要集中在夏季的7,8月份,約占全年降水量的50%,該區(qū)域多年平均地表蒸發(fā)量為1479mm。研究區(qū)地勢(shì)由東南向西北逐漸降低,沿海地勢(shì)平緩,海拔高程多在0~15m之間,屬于濱海平原區(qū)。該區(qū)屬于龍口-威海侵入巖變質(zhì)巖沖洪積海積工程地質(zhì)亞區(qū),含水層呈多元結(jié)構(gòu),一般發(fā)育有3~4層,且其表層主要以砂性土為主,砂性土以下依次以粉土和粉質(zhì)黏土為主,區(qū)域內(nèi)地下水從東南流向西北的海域,主要是第四系的松散巖類孔隙水,水化學(xué)類型為氯化物型,水交替強(qiáng)烈。
該區(qū)地下水的賦存條件差,富水性弱。僅在局部大理巖分布區(qū)及斷裂構(gòu)造發(fā)育帶,為地下水的賦存創(chuàng)造小面積的或點(diǎn)狀、線狀的富集環(huán)境,形成小區(qū)域大理巖巖溶水富水區(qū)及變質(zhì)巖、侵入巖裂隙水點(diǎn)或線狀較富水帶。區(qū)內(nèi)河流發(fā)育短小,第四系松散層多呈條帶狀分布于現(xiàn)代河床及兩側(cè)附近,厚度小,富水性弱,部分河流谷地較開闊,形成小型河谷平原及沿岸山前傾斜平原孔隙水較富水地段。區(qū)內(nèi)基巖裂隙水廣泛分布,可分為層狀巖類裂隙水、塊狀巖類裂隙水及噴出巖孔洞裂隙水,大部分屬風(fēng)化裂隙水。一般富水性較弱,個(gè)別地段孔洞裂隙水水量較大,如蓬萊市南王鎮(zhèn)一帶,形成小的富水地段。山區(qū)由于侵蝕切割作用,形成水系較多,基巖裂隙水主要靠大氣降水補(bǔ)給。碳酸鹽巖類巖溶裂隙水主要分布在膠東半島的大理巖分布區(qū),含水層主要為張格莊組大理巖,另外在蓬萊南王鎮(zhèn)一帶隱伏的張格莊組大理巖和牟平養(yǎng)馬島野頭組大理巖也形成了富水地段。區(qū)內(nèi)主要分布?jí)酃馑吹?、濰寒水源地、濰黃旗堡水源地、朱里水源地、穆村水源地、昌邑水源地等龍口、招遠(yuǎn)地區(qū)城市供水近70%來自地下水,第四系松散層中的孔隙潛水是主要的開采對(duì)象,地下水水資源的超量開采使該地區(qū)遭受海水入侵,地下水被污染。
該區(qū)海域洋流水動(dòng)力特征季節(jié)性顯著,洋流分布特性如圖3、圖4所示。
圖3 萊州灣及附近海域冬季洋流分布示意圖
圖4 萊州灣及附近海域夏季洋流分布示意圖
對(duì)研究區(qū)域地形地貌、地質(zhì)構(gòu)造以及地層結(jié)構(gòu)進(jìn)行調(diào)查和分析,認(rèn)為該區(qū)域地質(zhì)結(jié)構(gòu)完整,各含水層結(jié)構(gòu)基本連續(xù),可利用GMS軟件MODFLOW模塊對(duì)研究區(qū)地下水流場(chǎng)進(jìn)行數(shù)值模擬。
2.2.1 水文地質(zhì)概念模型
建立該區(qū)域的水文地質(zhì)概念模型,主要包括含水層系統(tǒng)、邊界條件、源匯項(xiàng)等3個(gè)方面。
(1)含水層系統(tǒng):龍口、招遠(yuǎn)濱海研究區(qū)內(nèi)地下各含水層的分布較連續(xù),垂向上呈明顯的層狀分布,根據(jù)研究區(qū)地質(zhì)資料,可將研究區(qū)地下水含水層概化為3層,由地表向下,含水層巖性依次以砂性土、粉土和粉質(zhì)黏土為主。
(2)邊界條件:垂向上,研究區(qū)含水層上部主要接受大氣降水補(bǔ)給,為水量交換邊界,下部與相對(duì)隔水的黏土層、基巖接觸,為隔水邊界;水平向上,模擬區(qū)北側(cè)、西側(cè)和西北側(cè)與渤海相接,東北角和西南角分別以黃水河、界河為邊界,都概化為一類定水頭邊界,南側(cè)邊界為龍口內(nèi)陸低山丘陵區(qū),概化為二類定流量邊界,研究區(qū)邊界條件及地下水位觀測(cè)點(diǎn)如圖5所示。
圖5 研究區(qū)邊界條件及地下水位觀測(cè)點(diǎn)圖
(3)源匯項(xiàng):研究區(qū)內(nèi)地下水的補(bǔ)給以降水入滲和山丘區(qū)基巖裂隙水的側(cè)向徑流為主;研究區(qū)瀕臨渤海灣,屬濱海平原地帶,區(qū)域內(nèi)的地下水位埋深較淺,普遍小于5m,需要考慮地下水的蒸發(fā)量。同時(shí),孔隙潛水的開采也是該區(qū)域地下水排泄的一個(gè)重要途徑,根據(jù)研究區(qū)的統(tǒng)計(jì)資料,確定開采井的位置、開采深度和具體的開采量。
2.2.2 地下水流場(chǎng)數(shù)學(xué)模型建立及求解
在對(duì)區(qū)域水文地質(zhì)概念模型分析研究的基礎(chǔ)上,對(duì)地下水水流進(jìn)行概化,建立非均質(zhì)、各向同性的二維非穩(wěn)定地下水流數(shù)學(xué)模型。
H(x,y,t)|t=0=H0(x,y,t) (x,y)∈D
H(x,y,t)|Γ1=H1(x,y,t) (x,y)∈Γ1,t>0
(5)
式中:K—潛水含水層滲透系數(shù)(m/d);μ—潛水含水層給水度;Z—含水層底板標(biāo)高(m);H—潛水水位(m);H1—一類邊界點(diǎn)的水位(m);q—二類邊界單寬流量(m3/d/m);D—計(jì)算區(qū)范圍;Γ1,Γ2—一、二類邊界。
利用GMS軟件中的MODFLOW模塊計(jì)算研究區(qū)地下水流場(chǎng),模擬區(qū)域面積約為424km2,在GMS軟件3DGRID模塊下對(duì)研究進(jìn)行網(wǎng)格剖分,共剖分3384個(gè)有效矩形網(wǎng)格單元(圖6)。
圖6 研究區(qū)網(wǎng)格剖分立體示意圖
根據(jù)區(qū)域水文地質(zhì)條件,研究區(qū)內(nèi)各含水層水文地質(zhì)參數(shù)(給水度、滲透系數(shù)、有效孔隙度)初始值如表1所示。
2.2.3 識(shí)別和檢驗(yàn)地下水流場(chǎng)數(shù)值模型
地下水流場(chǎng)數(shù)值模型的識(shí)別和檢驗(yàn)是相互獨(dú)立的2個(gè)階段,需要根據(jù)不同時(shí)間段的資料分別完成。
表1 研究區(qū)各含水層水文地質(zhì)參數(shù)
該文以2009年7月—2010年6月作為模型識(shí)別期,反求水文地質(zhì)參數(shù),再次對(duì)研究區(qū)的地下水流場(chǎng)進(jìn)行數(shù)值模擬;以2010年7月—2011年6月作為地下水流場(chǎng)數(shù)值模型的檢驗(yàn)期,主要檢驗(yàn)?zāi)M研究區(qū)的地下水水位動(dòng)態(tài)變化是否與實(shí)測(cè)的地下水位動(dòng)態(tài)變化一致,兩個(gè)地下水位觀測(cè)點(diǎn)在檢驗(yàn)期內(nèi)實(shí)測(cè)地下水位值與計(jì)算值的對(duì)比如圖7、圖8所示。
圖7 2010年7月--2011年6月海岱鎮(zhèn)河南孫家水位觀測(cè)點(diǎn)地下水位觀測(cè)值與計(jì)算值對(duì)比圖
圖8 2010年7月--2011年6月徐福鎮(zhèn)洼西村水位觀測(cè)點(diǎn)地下水位觀測(cè)值與計(jì)算值對(duì)比圖
從圖7和圖8可知,海岱鎮(zhèn)河南孫家和徐福鎮(zhèn)洼西村水位監(jiān)測(cè)點(diǎn)地下水位觀測(cè)值與計(jì)算值之間的最大差值為0.48m,擬合程度較好,滿足《地下水流技術(shù)數(shù)值模擬技術(shù)要求》GW1-D1中“水位變化值較小(小于5m)的情況下,水位擬合誤差一般應(yīng)小于0.5m”的要求。
此外,對(duì)地下水流場(chǎng)、源匯項(xiàng)均衡量和水文地質(zhì)參數(shù)分別進(jìn)行檢驗(yàn),各項(xiàng)與實(shí)際情況的擬合程度都較好,認(rèn)為建立的地下水流場(chǎng)數(shù)值模型計(jì)算結(jié)果達(dá)到精度要求,反映了該區(qū)地下水系統(tǒng)的動(dòng)力特征,可以用于后續(xù)地下水溶質(zhì)運(yùn)移模擬計(jì)算。
2.3.1 地下水超采狀態(tài)下海陸過渡帶遷移
在研究區(qū)地下水流場(chǎng)計(jì)算的基礎(chǔ)上,對(duì)該地區(qū)海水中Cl-在內(nèi)陸地下水中的運(yùn)移進(jìn)行數(shù)值模擬,因?yàn)镃l-是海咸水中最穩(wěn)定和最主要常量組成元素,可以直接反映出海咸水運(yùn)移的狀況。在龍口、招遠(yuǎn)濱海區(qū)域,可以將地下水中Cl-濃度超過200mg/L的地區(qū)視為海水入侵區(qū)域[11]。由于研究區(qū)地下水處于超采狀態(tài),故研究海陸過渡帶的遷移是以地下水超采作為條件。
對(duì)研究區(qū)西北部約12km長(zhǎng)海岸線上海咸水中的Cl-向內(nèi)陸地下水中的運(yùn)移過程進(jìn)行數(shù)值模擬,并給作用于海岸線上的濃度賦初值為3g/L。模擬運(yùn)移的計(jì)算期為7200d,每300d作為一個(gè)計(jì)算應(yīng)力期,共分為24個(gè)計(jì)算應(yīng)力期。
在內(nèi)陸地下水中Cl-運(yùn)移1000d,3000d,6000d后在第一層、第二層含水層中水平向的濃度分布如圖9和圖10所示。
圖9 Cl-在內(nèi)陸地下水中運(yùn)移1000d,3000d,6000d后在第一層含水層中水平向濃度分布圖
圖10 Cl-在內(nèi)陸地下水中運(yùn)移1000d,3000d,6000d后在第二層含水層中水平向濃度分布圖
研究海咸水中Cl-濃度在含水層中垂向上的變化,在模型中選取4個(gè)剖面,提出各剖面Cl-上運(yùn)移6000d后的垂向濃度分布圖。選取的剖面位置如圖11所示。
模型I-I剖面、J-J剖面、K-K剖面、L-L剖面上Cl-運(yùn)移6000d后,在含水層中垂向上濃度變化如圖12所示。
由圖9、圖10、圖11,Cl-在含水層中運(yùn)移若干天后的水平向和垂向濃度分布如下:
(1)當(dāng)內(nèi)陸地下水超采,海咸水入侵,海咸水中的Cl-在內(nèi)陸地下水中的運(yùn)移主要發(fā)生在第一層砂性土層中,第二層粉土層受Cl-的影響較小,第三層粉質(zhì)黏土層受Cl-的影響,可以忽略不計(jì),即海咸水入侵污染內(nèi)陸地下水的最大深度為兩層含水層的厚度。
圖11 Cl-濃度在垂向變化剖面位置圖
a—I-I剖面垂向Cl-濃度變化圖;b—J-J剖面垂向Cl-濃度變化圖;c—K-K剖面垂向Cl-濃度變化圖;d—L-L剖面Cl-垂向濃度變化圖圖12 模型I-I剖面、J-J剖面、K-K剖面、L-L剖面上Cl-運(yùn)移6000d后,在含水層中垂向上濃度變化圖
(2)內(nèi)陸地下水水位下降,海咸水中的Cl-向沿海地區(qū)內(nèi)陸地下水運(yùn)移入侵,主要發(fā)生在前1000d之內(nèi),這期間內(nèi)地下水中Cl-的濃度變化較大,3000d以后,地下水中Cl-運(yùn)移趨于平緩,該Cl-運(yùn)移特性可以作為防控海咸水入侵的依據(jù)。
2.3.2 Cl-濃度隨時(shí)間變化
選取距海岸線不同距離的A,B,C,D,E,F(xiàn)等6個(gè)點(diǎn),并從模型中提取出6個(gè)點(diǎn)處Cl-濃度變化曲線,研究該點(diǎn)處Cl-濃度隨時(shí)間變化的特性。A,B,C,D,E,F(xiàn)距離海岸線的距離分別約為3km,4.5km,8.4km,12.5km,14.8km,16km,各點(diǎn)位置如13所示。
圖13 距海岸線不同距離的A,B,C,D,E,F(xiàn)等6個(gè)點(diǎn)位置示意圖
A,B,C,D,E,F(xiàn)各點(diǎn)處Cl-濃度隨時(shí)間變化如圖14所示。
距海岸線4.5km的B點(diǎn)處濃度運(yùn)移模擬在500d時(shí)達(dá)到200mg/L,故綜合考慮研究區(qū)在第一層含水層中運(yùn)移3000d后的濃度分布圖和以上六點(diǎn)處濃度隨時(shí)間變化曲線圖可知,沿海地區(qū)存在濃度恒大于200mg/L的地區(qū),該類區(qū)域已被海咸水入侵污染,寬度約為1.5~4.5km,即研究區(qū)龍口、招遠(yuǎn)濱海區(qū)域海陸過渡帶的寬度約為1.5~4.5km。南京大學(xué)薛禹群院士團(tuán)隊(duì)曾在20世紀(jì)通過大量的野外觀察,對(duì)該地區(qū)的海水入侵狀況進(jìn)行了研究,認(rèn)為該地區(qū)海陸過渡帶的寬度約為1.5~2.5km。該文通過對(duì)氯離子運(yùn)移進(jìn)行數(shù)值模擬,對(duì)研究區(qū)劃分的海陸過渡帶為1.5~4.5km,比薛禹群院士團(tuán)隊(duì)研究的海陸過渡帶寬度稍寬,說明經(jīng)過20多年,該地區(qū)的海咸水繼續(xù)向內(nèi)陸運(yùn)移,地下水繼續(xù)被超采,海水入侵未得到有效遏制。
2.3.3 內(nèi)陸地下水位對(duì)海陸過渡帶遷移的影響
在目前地下水流場(chǎng)的基礎(chǔ)上,先對(duì)研究區(qū)的地下水位進(jìn)行改變,將地下水位整體降低1.5m和抬高1.5m,再對(duì)研究區(qū)海咸水中Cl-運(yùn)移進(jìn)行數(shù)值模擬,分析內(nèi)陸地下水位對(duì)研究區(qū)海陸過渡帶遷移的影響。地下水位抬高1.5m和降低1.5m后,在第一層含水層中運(yùn)移1000d,3000d,6000d后的濃度分布如圖15所示。
圖14 A,B,C,D,E,F(xiàn)各點(diǎn)Cl-濃度隨時(shí)間變化曲線圖
a—地下水位抬高1.5m后,Cl-在第一層含水層中運(yùn)移1000d,3000d,6000d后的濃度分布圖;b—地下水位降低1.5m后,Cl-在第一層含水層中運(yùn)移1000d,3000d,6000d后的濃度分布圖圖15 地下水位抬高1.5m和降低1.5m后,Cl-在第一層含水層中運(yùn)移1000d,3000d,6000d后的濃度分布圖
由以上不同地下水位條件下,Cl-運(yùn)移相同時(shí)間后的濃度分布圖可知,內(nèi)陸地下水位對(duì)海咸水中Cl-的運(yùn)移有顯著影響,即內(nèi)陸地下水位對(duì)海陸過渡帶的遷移有重要的影響,內(nèi)陸地下水位越高,海陸過渡帶向內(nèi)陸遷移的范圍就越小,速度也越慢,這對(duì)控制沿海地區(qū)海咸水入侵污染內(nèi)陸地下水也具有重要的指導(dǎo)意義。
該文選擇海陸過渡帶模型,結(jié)合地下水運(yùn)動(dòng)理論,選取山東半島藍(lán)色經(jīng)濟(jì)區(qū)龍口、招遠(yuǎn)濱海地區(qū)作為研究區(qū)域,在對(duì)基礎(chǔ)地質(zhì)資料和水文地質(zhì)資料研究分析的基礎(chǔ)上,利用地下水?dāng)?shù)值模型系統(tǒng)GMS計(jì)算該地區(qū)的流場(chǎng),并在內(nèi)陸地下水超采狀態(tài)下,海咸水入侵對(duì)該地區(qū)地下水中Cl-的運(yùn)移進(jìn)行數(shù)值模擬,研究該地區(qū)海陸過渡帶遷移的規(guī)律。
通過研究,該地區(qū)海陸過渡帶的寬度約為1.5~4.5km,與1992年1.5~2.5km寬的海陸過渡帶相比,該地區(qū)的海陸過渡帶明顯變寬,海咸水繼續(xù)向內(nèi)陸運(yùn)移。此外,內(nèi)陸地下水位對(duì)海陸過渡帶的遷移具有重要的影響,內(nèi)陸地下水位越低,海陸過渡帶向內(nèi)陸遷移的范圍越大、速度越快。
因此,為有效遏制海咸水入侵內(nèi)陸地下水,建議采取限采、回灌、地下水庫(kù)等相關(guān)措施抬高該地區(qū)的地下水位,如:在該區(qū)域內(nèi)選擇具有足夠的天然地下存儲(chǔ)空間和充足水源的地區(qū)修建地下水庫(kù);修建攔蓄補(bǔ)源工程,充分利用雨洪水資源,回灌補(bǔ)給地下水;嚴(yán)格控制地下水開采。