濰河下游典型鹵水開采區(qū)海咸水入侵趨勢研究

蒙永輝 王集寧 張麗霞 羅梅

摘要：為了掌握濰河下游鹵水開采對海咸水入侵的影響，在地質(zhì)、水文地質(zhì)、地下水位和地下水水化學(xué)等實測資料的基礎(chǔ)上，基于濰河下游典型鹵水貯存、運動特征和開采現(xiàn)狀，建立了濰河下游典型鹵水開采區(qū)的地下水流系統(tǒng)數(shù)值模型，并對鹵水開采區(qū)不同開采條件下海水入侵的趨勢進行了預(yù)測。結(jié)果表明：在保持目前地下水和鹵水資源開采狀態(tài)的情況下，海咸水入侵鋒線向鹵水區(qū)推進，且鹵水區(qū)地下水礦化度大幅度降低；若保持地下水開采、停止鹵水資源開采，則鹵水區(qū)地下水位逐漸恢復(fù)，鹵水區(qū)地下水礦化度有降低現(xiàn)象；海咸水入侵鋒線處于相對穩(wěn)定狀態(tài)，原因是其處于地下水降落漏斗中心處，南部低于鋒線值和北部高于鋒線值地下水均向漏斗區(qū)匯聚，使得海咸水入侵鋒線處于相對穩(wěn)定狀態(tài)。

關(guān)鍵詞：海咸水入侵；典型鹵水開采區(qū)；數(shù)值模擬；濰河下游

中圖分類號：P314.1 文獻標(biāo)志碼：A doi：10.3969/j.issn.1000-1379.2018.01.015

地下鹵水資源具有很高的經(jīng)濟價值，可用于軍工、鹽化工、電子、制藥等領(lǐng)域[1]。

濰河下游位于萊州灣中部，具有得天獨厚的地下鹵水資源，是全國最大的海鹽及鹽化工生產(chǎn)基地，截至2006年年底，濰坊北部沿海地區(qū)共有鹵水開發(fā)企業(yè)4000多家[2]。近年，鹵水資源開采導(dǎo)致了一系列的生態(tài)環(huán)境問題，包括水位下降、咸淡水分界線向內(nèi)陸遷移、海水入侵、海洋生態(tài)環(huán)境破壞、鹵水質(zhì)量濃度和有益組分含量降低等[3]，其中海水入侵的危害尤為突出，目前該地區(qū)海水入侵面積已超過4000km2，入侵陸地的速度達500～1000m/a[4]。日益加劇的海水人侵嚴重污染了地下水環(huán)境，使土壤次生鹽漬化，對海岸區(qū)域的工農(nóng)業(yè)生產(chǎn)和生態(tài)環(huán)境造成嚴重的影響。因此，研究鹵水資源開采與海水入侵的關(guān)系、評估海水人侵的趨勢，具有重要的理論和實際意義。國內(nèi)外學(xué)者對此做了很多研究工作，較為常用的研究方法包括數(shù)值模擬、水文地球化學(xué)和同位素方法等。Oahman等[5]利用二維有限差分模型SEAWAT模擬了Gaza海水入侵狀況，認為地下水開采使海水入侵狀況惡化；Karahanoglu等[6]利用有限元方法模擬了土耳其Kocaeli-Darica采石場海岸含水層海水入侵狀況，研究了海水入侵行為、地下水位變化規(guī)律和鹽分濃度的關(guān)系；Mahlknecht等[7]利用地下水中主要組分、Sr和B同位素，評價了人類活動影響下干旱地區(qū)海岸含水層的海水入侵狀況；Cimino等[8]利用地球物理和地球化學(xué)方法，評價了西西里島北部Acquedolci海岸含水層海水入侵狀況，并提出了咸淡水混合模型；我國學(xué)者林文盤等[9]探討了萊州灣海水入侵災(zāi)害防治的進展，提出海水入侵的機理本質(zhì)上是海岸帶咸淡水界面在地下水位變化情況下水動力平衡被破壞；章斌等[10]根據(jù)秦皇島洋戴河平原地下潛水和海水的水化學(xué)特征，運用數(shù)理統(tǒng)計和模糊數(shù)學(xué)方法研究了該地區(qū)潛水含水層的海水入侵程度；李國敏等[11]根據(jù)海水入侵研究現(xiàn)狀，認為數(shù)值方法是現(xiàn)今模擬和求解海水入侵問題的最有力工具；段梅[12]綜合分析硇洲島的地質(zhì)和水文地質(zhì)條件，建立了海水入侵數(shù)值模型，對硇洲島在不同開采條件下的海水入侵趨勢進行了模擬預(yù)測。但是，上述研究僅局限于單一影響因素，對海平面水位、風(fēng)暴潮等影響因素鮮有考慮。

本研究利用濰河下游典型鹵水開采區(qū)地下水位、水質(zhì)長期動態(tài)系列資料，以及海平面水位、風(fēng)暴潮侵入次數(shù)和距離等數(shù)據(jù)，建立以地下水礦化度為判斷因子，以地下水對流遷移為主導(dǎo)，并考慮分子擴散、吸附作用、密度變化因素的變密度三維地下水流模型和溶質(zhì)運移模型，對未來海咸水入侵的發(fā)展趨勢進行預(yù)測，旨在為研究區(qū)鹵水開發(fā)利用和保護提供技術(shù)參考。

1 研究區(qū)概況

研究區(qū)位于萊州灣南岸濱海平原東部、濰河下游，面積約946km2。氣候?qū)倥瘻貛О霛駶櫞箨懶詺夂騾^(qū)，四季分明，氣候溫和，多年平均氣溫為11.9℃，多年平均年降水量560mm，其中6-9月降水量占年降水量的55%～77%，年均水面蒸發(fā)量為1788mm。地貌類型由南至北可分為南部山前沖洪積平原、中部沖積海積平原和北部海積平原。地層主要為第四系更新統(tǒng)至全新統(tǒng)沖洪積、海積沉積層，巖性以粉細沙、中粗沙和礫沙為主。地下水主要賦存于第四系沙礫石層等含水介質(zhì)中，黏性土作為相對隔水層，形成松散巖類孔隙式的多層含水結(jié)構(gòu)，地下水類型主要為松散巖類孔隙水，垂向分為淺層含水巖組和深層含水巖組。淺層含水巖組主要指深度小于120m的淺層潛水、微承壓水，含水層巖性從山前平原中粗沙、礫石變?yōu)闉I海區(qū)的粉沙、細沙，分為淺層淡水、淺層微咸水和淺層咸水，其中：淺層淡水主要賦存于山前沖洪積扇、濰河及其古河道堆積形成的河谷、階地含水層中，礦化度<1g/L，單井涌水量為500～3000m3/d，水化學(xué)類型為HCO3·Cl-Na；淺層微咸水分布于柳瞳以北，礦化度為1～3g/L，單井涌水量1000～3000m3 /d，水化學(xué)類型為HCO3·Cl-Na；淺層咸水分布于第四系海相地層的松散沉積物中，礦化度>3g/L，水化學(xué)類型為Cl-Na。深層含水巖組主要分布于海積、沖洪積海陸交互堆積平原，巖性以中沙和細沙為主，單井涌水量為1000～3000m3/d，礦化度為1～2g/L，水化學(xué)類型為HCO3·Cl-Na。

2 數(shù)據(jù)來源和研究方法

2.1 數(shù)據(jù)來源

研究所用數(shù)據(jù)來自山東省地質(zhì)環(huán)境監(jiān)測總站和濰坊市礦產(chǎn)資源管理中心，包括氣象、水文、地形地貌、地質(zhì)、水文地質(zhì)、地下水位、地下水水化學(xué)和地下水開采量等數(shù)據(jù)。目前，研究區(qū)鹵水資源年開采量為4800萬m3。

2.2 水文地質(zhì)概念模型

2.2.1 地下水系統(tǒng)邊界條件的概化

研究區(qū)處于濰河下游，南部以濰河沖洪積扇頂部為界，北部到海岸線，東部以膠萊河為界，西部以濰河與白浪河沖積地層的交接處為界。南部、東部和西部邊界均為流量邊界；北部海岸線邊界為給定水位邊界，即按照渤海灣海平面的觀測，并考慮海水潮位變化，結(jié)合實際觀測的海水潮位變化給定水位邊界。垂向上，由于很多區(qū)域淺部鹵水含水層（潛水層）地下水位已低于含水層底板，含水層被疏干，因此將淺部鹵水含水層與第一承壓（微承壓）鹵水含水層概化為一層，下部各承壓鹵水含水層概化為一層，咸水層以下深層淡水概化為一層。按照地質(zhì)剖面地層結(jié)構(gòu)，從北部鹵水區(qū)向南將整個含水層結(jié)構(gòu)概化為3層含水層和2層相對隔水層（弱透水層），把地下水系統(tǒng)概化為非均質(zhì)各向異性三維非穩(wěn)定流。

2.2.2 模型初始值確定

（1）初始流場和初始濃度場。將2015年12月的地下水流場作為模型的初始流場（水位等值線），見圖1。把地下水礦化度作為海咸水入侵的鑒別因子，淺部地下水初始濃度場：南部以2015年年底地下水礦化度實測值為初始濃度值；北部鹵水區(qū)根據(jù)2007年鹵水勘察的鹵水波美度值，依據(jù)鹵水波美度與濃度的轉(zhuǎn)化關(guān)系（表1）和研究區(qū)鹵水主要陰陽離子比例分布，將鹵水波美度轉(zhuǎn)化為濃度，近似作為鹵水礦化度值，并將其作為第一含水層和第二含水層的礦化度初值，形成第一、第二含水層模擬的初始濃度場（圖2）。依據(jù)南部深層水界線上的礦化度情況，設(shè)定深層含水層（第三含水層）地下水礦化度均為2000mg/L。

（2）源匯項處理。源匯項主要為大氣降水入滲、潛水蒸發(fā)和地下水人工開采。2016年研究區(qū)地下水資源開采量為5 177.5萬m3，鹵水資源開采量為4800萬m3，地下水和鹵水資源開采區(qū)均位于研究區(qū)北部的漏斗區(qū)，在模型中均以面狀開采形式體現(xiàn)。

2.3 數(shù)學(xué)模型

根據(jù)上述概化的水文地質(zhì)概念模型，建立研究區(qū)地下水概念模型和地下水溶質(zhì)運移模型，其數(shù)學(xué)模型（控制方程）和定解條件為式中：H為地下水位（相對于淡水）；Kij為滲透系數(shù)張量（i，j=1，2，3）；η為密度耦合系數(shù)，η=ε/Cs，ε為密度差率，ε=（ρs-ρ0）/ρ0，Cs為與流體最大密度ρs對應(yīng)的濃度，ρ0為參考密度（淡水密度）；ρ為混合溶液的密度；ρ*為現(xiàn)狀條件下混合液密度；C為溶液濃度；ej為重力方向單位矢量第j個分量；AR為貯水率；φ為孔隙率；q為單位體積孔隙介質(zhì)源（或匯）流量；xi、xj為笛卡兒坐標(biāo)（i，j=1，2，3）；H0為地下水初始水位；HB為邊界Г1上的給定地下水位；qB2為邊界Г2上補排強度；t為時間；W'為潛水面邊界Г3上的補排強度；H*為潛水面邊界Г3上各點地下水位；He為潛水面邊界Г3上各點的高程；ni為邊界上外法線單位矢量；D為水動力彌散系數(shù)張量；μi為地下水實際流速在xi方向的分量；C*為抽出或注入液體的濃度；C0為初始濃度；CB為邊界Г1上的濃度；C'為降水入滲的濃度。

2.4 模型識別及結(jié)果

利用2015—2016年的地下水動態(tài)資料對模型進行調(diào)試與識別。把2015年年底實測的地下水流場作為初始流場，對典型觀測孔實測的地下水位動態(tài)曲線和2016年年底的地下水流場、地下水濃度（礦化度）場進行擬合。圖3、圖4、圖5分別為長觀孔地下水位、地下水流場和地下水濃度場擬合結(jié)果。經(jīng)模型識別及檢驗，將全區(qū)分為5個降雨入滲系數(shù)分區(qū)，將含水層和弱透水層分別分為12、8個水文地質(zhì)參數(shù)分區(qū)，見表2～表4和圖6～圖8。由于V區(qū)為曬鹽池，且已作防滲處理，因此該區(qū)大氣降水入滲系數(shù)為0?？紤]到第四系含水層中夾有黏性土以及垂向上的水流沉積作用，模型中第四系水平滲透系數(shù)取垂向滲透系數(shù)的10倍。

從典型觀測孔模擬水位和實測水位的擬合曲線可看出，水位模擬誤差較小，所建模型能夠反映研究區(qū)地下水流場的動態(tài)變化，可以預(yù)測不同開采條件下地下水的動態(tài)變化。

2.5 預(yù)測情景

海咸水入侵預(yù)報設(shè)置了兩種情景：一是保持目前地下水開采和鹵水開采狀態(tài)，預(yù)測未來5、10、15a咸淡水界面變化情況；二是保持目前地下水開采狀態(tài)，鹵水開采在5a后停止，預(yù)測未來10、15a海咸水入侵情況。預(yù)報過程中，保持研究區(qū)水文地質(zhì)參數(shù)不變，大氣降水量按照2001-2015年的實際降水序列加載到模型中。在上述條件下，利用驗證后的模型進行海咸水入侵預(yù)測。

3 結(jié)果與分析

預(yù)測情景一和情景二的前5a是相同的，均保持目前地下水資源和鹵水資源的開采量，預(yù)測結(jié)果見圖9；圖10，圖11分別為按情景一和情景二預(yù)測的10、15a海咸水入侵情況（圖中地下水位單位為m，礦化度單位為mg/L）。

從圖9可以看出，若保持目前地下水開采和鹵水資源開發(fā)狀態(tài)，則鹵水區(qū)水位降速大于地下水資源開采區(qū)的水位降速，海咸水入侵鋒線向鹵水區(qū)推進。

從圖10可以看出，按照情景一保持目前地下水開采和鹵水資源開發(fā)狀態(tài)，未來10、15a鹵水區(qū)水位降速仍大于地下水資源開采區(qū)的水位降速，海咸水入侵鋒線向鹵水區(qū)推進，且鹵水區(qū)地下水礦化度大幅度降低。

從圖11可以看出，按照情景二保持目前地下水開采和鹵水資源開發(fā)5a后，停止鹵水資源的開發(fā)，仍保持地下水開采，則地下水資源開采區(qū)水位不斷下降、鹵水區(qū)地下水位逐漸恢復(fù)，由于補給水源礦化度相對較低，因此鹵水區(qū)地下水礦化度有所降低，但是之前鹵水區(qū)礦化度太高，所以其變化不太明顯。海咸水入侵鋒線與前5a基本一致，主要原因是海咸水入侵鋒線處于地下水降落漏斗中心處，南部低于鋒線值和北部高于鋒線值均向漏斗區(qū)匯聚，使得海咸水入侵鋒線處于相對穩(wěn)定狀態(tài)。

4 結(jié)論

（1）濰河下游典型鹵水開采區(qū)保持目前地下水和鹵水資源開采狀態(tài)時，鹵水區(qū)水位降速大于地下水資源開采區(qū)的水位降速，海咸水入侵鋒線向鹵水區(qū)推進，且鹵水區(qū)地下水礦化度大幅度降低；若保持地下水開采、停止鹵水資源的開采，則鹵水區(qū)地下水位逐漸恢復(fù)、地下水礦化度有所降低。

（2）研究區(qū)海咸水入侵鋒線處于相對穩(wěn)定狀態(tài)，主要原因是其處于地下水降落漏斗中心處，南部低于鋒線值和北部高于鋒線值均向漏斗區(qū)匯聚，使得其處于相對穩(wěn)定狀態(tài)。

參考文獻：

[1]鄒祖光，張東生，譚志容.山東省地下鹵水資源及開發(fā)利用現(xiàn)狀分析[J].地質(zhì)調(diào)查與研究，2008，31（3）：214-221.

[2]林存菊，姚英強，付娟.黃河三角洲高效生態(tài)經(jīng)濟區(qū)鹵水資源開采潛力評價[J].山東國土資源，2014，30（9）；48-52.

[3]陳廣泉，徐興永，彭昌盛，等.萊州灣地區(qū)海水入侵災(zāi)害風(fēng)險評價研究[J].自然災(zāi)害學(xué)報，2010，19（2）：103-112.

[4]苗青，陳廣泉，劉文全，等.萊州灣地區(qū)海水入侵災(zāi)害演化過程及成因[J].海岸工程，2013，32（2）：69-78.

[5]OAHMAN K，LARARI A.Evaluation and Numerical Modelingof Seawater Intrusion in the Gaza Aquifer（Palestine）[J].Hydrogeology Journal，2006，14（5）：713-728.

[6]KARAI4ANOGLU N，DOYURAN V.Finite Element Simula-tion of Seawater Intrusion into a Quarry-Site Coastal Aquifer，Kocaeli-Darica，Turkey[J].Environmental Geology，2003，44（4）：456-466.

[7]MAHLKNECHTJ，MERCHAN D，ROSNER M，et al.As-sessing Seawater Intrusion in an Arid Coastal Aquifer UnderHigh Anthropogenic Influence Using Major Constituents，Srand B Isotopes in Groundwater[J].Science of the Total En-vironment，2017，587：282-295.

[8]CIMINO A，COSENTINO C，OIENI A，et al.A Geophysicaland Geochemical Approach for Seawater Intrusion Assessmentin the Acquedolci Coastal Aquifer（Northern Sicily）[J].Envi-ronmental Geology，2008，55（7）：1473-1482.

[9]林文盤，尹澤生.萊州灣海水入侵災(zāi)害防治研究進展[J].中國減災(zāi)，1992，2（1）：31-34.

[10]章斌，宋獻方，韓冬梅，等.運用數(shù)理統(tǒng)計和模糊數(shù)學(xué)評價秦皇島洋戴河平原的海水入侵程度[J].地理科學(xué)，2013，33（3）：342-348.

[11]李國敏，陳崇希.海水入侵研究現(xiàn)狀與展望[J].地學(xué)前緣，1996，3（1）：161-168.

[12]段梅.硇洲島海水入侵數(shù)值模擬[J].地質(zhì)災(zāi)害與環(huán)境保護，2016，27（1）：108-112.