海陸過渡帶氯離子遷移規(guī)律分析研究

賈超，朱恒華，虞未江，衛(wèi)政潤(rùn)，劉治政，李雙

(1.山東大學(xué)海洋研究院，山東 青島 266237;2.山東省地質(zhì)調(diào)查院，山東 濟(jì)南 250014)

0 引言

隨著經(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í)有效的建議。

1 理論方法

1.1 海陸過渡帶模型

沿海地區(qū)地下含水層中的水力梯度基本上是傾向海洋，內(nèi)陸地下含水層中的大部分淡水最終流向海洋。因?yàn)楹Ｆ矫嬉韵碌牡叵潞畬又写嬖谟泻Ｏ趟?，所以流向海洋的淡水與海咸水之間存在一條接觸帶。此外，海咸水與淡水是兩種可以互溶的流體，在流體動(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.2 地下水運(yùn)動(dò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)

2 研究區(qū)海陸過渡帶遷移規(guī)律研究

2.1 研究區(qū)概況

研究區(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 萊州灣及附近海域夏季洋流分布示意圖

2.2 地下水流場(chǎng)計(jì)算

對(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 海陸過渡帶遷移規(guī)律研究

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)意義。

3 結(jié)論

該文選擇海陸過渡帶模型，結(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)格控制地下水開采。