截滲工程條件下濱海地下水源地地下水補(bǔ)給源解析

黃俞銘， 閻 妮，2，3， 鄭西來，2，3， 劉貫群，2，3

(1. 中國海洋大學(xué)環(huán)境科學(xué)與工程學(xué)院， 山東 青島 266100；2. 中國海洋大學(xué)山東省海洋環(huán)境地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室， 山東 青島 266100；3. 中國海洋大學(xué)海洋環(huán)境與生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室， 山東 青島 266100)

地下水是寶貴的淡水資源，為全世界大多數(shù)的城市提供生活、農(nóng)業(yè)和工業(yè)用水。在沿海地區(qū)，地下水需求量的增加引起地下水過量開采，從而導(dǎo)致地面沉降、海水入侵等環(huán)境地質(zhì)問題，海水入侵引起的地下水水質(zhì)惡化已經(jīng)成為了全球性的水文地質(zhì)問題[1]。為了防治海水入侵，在濱海地區(qū)修建了各類防治海水入侵的工程措施。濱海含水層本身就具有地層結(jié)構(gòu)復(fù)雜、水利條件脆弱等特點(diǎn)[2]，在防治工程措辭的影響下，濱海含水層變成了一個(gè)更為復(fù)雜的地質(zhì)體，地下水來源解析也因此變得十分困難。

國內(nèi)外研究者已經(jīng)對(duì)濱海含水層地下水來源展開了多方面的研究。地下水來源解析方法主要分為水化學(xué)法和模型計(jì)算法。水化學(xué)方法包括統(tǒng)計(jì)分析法[3]、離子比值法[4-5]、Piper圖解法[6]和吉布斯圖解法[7]等。水化學(xué)法已經(jīng)被廣泛應(yīng)用于不同海水入侵區(qū)地下水來源分析的研究，但是水化學(xué)法作為一種定性分析方法，無法定量解析出地下水的來源，特別是化學(xué)成分相似的補(bǔ)給源，僅計(jì)算1～2個(gè)保守物種來確定混合比例。同時(shí)研究者發(fā)現(xiàn)，對(duì)于同一地下水樣品，當(dāng)選擇不同的保守組分進(jìn)行混合比計(jì)算時(shí)，計(jì)算出的混合比可能會(huì)有很大的不同[8-9]。

基于這種情況，為了對(duì)地下水來源進(jìn)行定量解析，研究者建立了多種結(jié)合多元統(tǒng)計(jì)分析、水化學(xué)分析和數(shù)學(xué)計(jì)算在內(nèi)的模型計(jì)算方法。最小二乘法首先被廣泛應(yīng)用于地下水來源解析的研究[8]，其主要目的是為了克服分析誤差，將樣品點(diǎn)投影到不同端元的混合線上。

但是最小二乘法也僅僅使用了2～3種指標(biāo)進(jìn)行來源解析，缺少對(duì)樣品的整體認(rèn)識(shí)。在最小二乘法的基礎(chǔ)上，Hooper等[10]結(jié)合主成分分析(PCA)提出了端元混合分析(End-Member Mixing Analysis, EMMA)，該方法綜合選取地下水樣品中的保守和非保守組分同時(shí)作為參數(shù)，期望用最少的端元數(shù)量來解釋特定區(qū)域的地下水來源，該方法被廣泛應(yīng)用于地下水、地表水等不同水體以及沉積物、大氣污染物等研究領(lǐng)域[11]。但是，端元混合分析法不能計(jì)算混合比，因此，需要結(jié)合其他方法。Carrera等[12]基于最大似然法提出了混合計(jì)算(MIX calculations)程序，該方法應(yīng)用最大似然法計(jì)算計(jì)算混合比和端元組成的方法。該方法承認(rèn)端元組分的不確定性，通過非線性優(yōu)化考慮所有混合物，選擇合適的選擇示蹤劑來計(jì)算潛在的補(bǔ)給源和混合比例。驗(yàn)證表明，相比水化學(xué)法和其他模型計(jì)算法，混合計(jì)算程序的計(jì)算結(jié)果更接近真實(shí)情況。但是模型計(jì)算法簡(jiǎn)化了地下水流動(dòng)過程和其中可能存在的水文地球化學(xué)反應(yīng)，假定了地下水流動(dòng)的連續(xù)性，并在此基礎(chǔ)上建立溶質(zhì)和水的質(zhì)量平衡計(jì)算，可能無法準(zhǔn)確反映研究區(qū)復(fù)雜的水文地質(zhì)條件。因此，對(duì)特定研究區(qū)進(jìn)行長(zhǎng)期水文地質(zhì)監(jiān)測(cè)，為模型計(jì)算提供依據(jù)是地下水源解析過程中必不可少的步驟。

基于以上分析，本文以黃水河濱海地下水源地為研究對(duì)象，通過現(xiàn)場(chǎng)的地下水水質(zhì)的系統(tǒng)監(jiān)測(cè)，在水化學(xué)解析方法的基礎(chǔ)上，采用端元混合分析確定地下水潛在的補(bǔ)給來源，利用混合計(jì)算程序計(jì)算不同地下水來源的混合比，確定研究區(qū)地下水補(bǔ)給來源和混合比，為該區(qū)域水資源開發(fā)利用和地下水污染治理提供科學(xué)依據(jù)和技術(shù)支持。

1 研究區(qū)的環(huán)境背景

1.1 自然地理概況

研究區(qū)位于中國山東省煙臺(tái)市龍口市東北部，坐標(biāo)北緯37°37′—37°43′，東經(jīng)120°30′—120°37′(見圖1)。北部地區(qū)的海拔不到1 m，南部地區(qū)達(dá)到25 m,屬溫帶季風(fēng)氣候區(qū)。根據(jù)1978—2018年的氣象記錄資料，研究區(qū)年平均降雨量大約680.5 mm，多年平均蒸發(fā)強(qiáng)度為1 932 mm。研究區(qū)內(nèi)降雨在年內(nèi)呈現(xiàn)出分配不均的特點(diǎn)，每年6—9月為汛期，降水量占全年降水量的70%～77%。年平均氣溫約11.2～12.4 ℃，1月份平均氣溫最低，約-4.3—-7.0 ℃，歷史最高低溫為-17.9 ℃；7月份平均氣溫最高，約25.5～27.8 ℃，歷史最高氣溫為38.9 ℃[13]。

圖1 研究區(qū)位置與采樣點(diǎn)分布圖Fig.1 Location of study area and distribution of sampling points

黃水河為研究區(qū)的主要河流，發(fā)源于棲霞市，經(jīng)龍口境內(nèi)流入渤海，年平均流量約4.77 m3/s，平均徑流量約1.51 m3/s，流域面積1 066 km2。為充分利用黃水河徑流，在最上游修建了王屋水庫，王屋水庫為大(二)型水庫，總庫容1.21億m3。同時(shí)在河流中下游修建了攔河閘和地下壩等工程，組成了黃水河下游地下水庫，總庫容5 329萬m3，同時(shí)具有儲(chǔ)存地下水與防止海水入侵的作用[14]。

1.2 水文地質(zhì)條件概況

1.2.1 地下水賦存條件 研究區(qū)含水層為第四紀(jì)孔隙潛水含水層，厚度在0～70 m之間，含水層為沖洪積細(xì)砂、礫質(zhì)粗砂、中粗砂和卵礫石[12]。研究區(qū)含水層大致分為四個(gè)部分：① alQ4礫質(zhì)粗砂，為高滲透層，滲透系數(shù)為237.39 m/d；② mQ4中粗砂；③ pl-alQ3粗砂，分布于黃水河Ⅱ級(jí)階地，滲透系數(shù)約為126.5 m/d；④ pl-alQ2中粗砂，局部夾雜卵礫石，該層的滲透系數(shù)為28.09～130.26 m/d(見圖2)[14]。

圖2 研究區(qū)I-I’水文地質(zhì)剖面圖Fig.2 I-I’ hydrogeological profile of the study area

1.2.2 地下水的補(bǔ)給、徑流和排泄 研究區(qū)地下水補(bǔ)給來源主要為大氣降水、地下水側(cè)向補(bǔ)給和地表水，其中最主要的補(bǔ)給來源為大氣降水和地表水。地下水排泄途徑主要是人為開采，由于地下水埋深均大于3 m，因此蒸發(fā)較弱。在天然條件下，地下水面形態(tài)與地形基本一致，即南高北低。地下水水力坡度較小，為1～3‰，水流緩慢。

1.2.3 地下水動(dòng)態(tài)變化規(guī)律 研究區(qū)地下水受人類活動(dòng)和降雨量影響較大，因此在年內(nèi)呈現(xiàn)出明顯的季節(jié)變化規(guī)律。在春季，地下水水位顯著下降，5—6月達(dá)到地下水水位最低值；6—9月進(jìn)入汛期，地下水水位開始上升，8—9月達(dá)到地下水水位最高值。地下水水位的年際變化呈現(xiàn)出明顯的分配不均的現(xiàn)象，在枯水年或連續(xù)枯水年，地下水水位大幅度持續(xù)下降；豐水年或連續(xù)豐水年地下水水位又持續(xù)上升。

1.2.4 地下水水化學(xué)特征

1.2.4.1地下水樣品采集 分別在2019年4和9月以及2020年9月對(duì)黃水河地下水源地進(jìn)行3次采樣，每次采集樣品21個(gè)，包括19個(gè)地下水樣品、1個(gè)上游地表水樣品和1個(gè)海水樣品，并在2020年8月采集大氣降水樣品1個(gè)。采樣點(diǎn)選擇的主要依據(jù)包括研究區(qū)水文地質(zhì)條件、土地利用類型、海水入侵程度等。分別選取濱海區(qū)和農(nóng)業(yè)活動(dòng)區(qū)內(nèi)具有代表性的采樣點(diǎn)。采樣井主要為工業(yè)用井、農(nóng)業(yè)用井和已發(fā)表研究成果[15]的監(jiān)測(cè)井。樣品收集方法以及室內(nèi)水質(zhì)監(jiān)測(cè)方法參照《中華人民共和國地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn) 地下水質(zhì)檢驗(yàn)方法》(DZ/T 0064.1～0064.90—93)。對(duì)得到的離子組分進(jìn)行電荷平衡誤差(Charge balance errors, C.B.E.)計(jì)算，保證所有樣品的電荷平衡誤差小于5%。在現(xiàn)場(chǎng)用便攜式Y(jié)SI型多參數(shù)水質(zhì)測(cè)定儀原位測(cè)定地下水的溫度、TDS、pH、EC、DO和ORP指標(biāo)。

1.2.4.2 地下水水化學(xué)特征 表1列出了研究區(qū)所有地下水樣品點(diǎn)的物理指標(biāo)和化學(xué)組成情況，包括每個(gè)指標(biāo)的最大值、最小值、標(biāo)準(zhǔn)偏差和平均值。研究區(qū)多數(shù)地下水樣品沒有表現(xiàn)出明顯受到海水入侵的地下咸水的離子組成特征，水樣總體的TDS值較低，且γ(Na+)的平均值低于γ(Ca2+)，與1994年[15]相比，大部分地下水樣品中Cl-濃度顯著降低。接近海邊的地下水樣品(1、2、3、4和19號(hào))中Cl-濃度較高，體現(xiàn)出了受海水入侵影響的特征。因此根據(jù)采樣點(diǎn)地理位置、水化學(xué)特征和土地利用類型，將研究區(qū)分為濱海區(qū)和農(nóng)業(yè)活動(dòng)區(qū)，其中濱海區(qū)包括的地下水采樣點(diǎn)為1、2、3、4和19號(hào)，其余采樣點(diǎn)位于農(nóng)業(yè)活動(dòng)區(qū)內(nèi)(見圖1)。

表1 研究區(qū)地下水水化學(xué)分析結(jié)果統(tǒng)計(jì)Table 1 The statistical table of groundwater hydrochemical monitoring in the study area

2 研究區(qū)地下水來源解析

2.1 水化學(xué)分析方法

2.1.1 主要陽離子與氯離子的比值 在受到海水入侵的地下水中，通常TDS主要受到海水的影響，因此主要陽離子與氯離子的比值通常接近海水[16]。主要陽離子與氯離子的比值與對(duì)應(yīng)樣品的氯離子濃度相結(jié)合，可以作為判斷地下水中離子的海洋來源的重要依據(jù)。圖3為研究區(qū)不同水體主要陽離子與氯離子的關(guān)系圖。在Na+與Cl-的比值圖以及K+與Cl-的比值圖中，低Cl-濃度的樣品中更接近海水中的比值，這與通常的海水入侵的研究結(jié)果不同。低Cl-濃度的地下水的Na+和K+與Cl-的比值接近海水的主要原因可能是巖石的溶濾風(fēng)化作用導(dǎo)致巖鹽溶解；而濱海區(qū)地下水中Na+與Cl-的比值低與海水的原因是截滲工程使海水入侵程度降低，濱海區(qū)地下水由咸水轉(zhuǎn)變?yōu)槲⑾趟?，Na+和Cl-的濃度降低，同時(shí)發(fā)生了Na+與Ca2+、Mg2+之間的陽離子交換作用。

圖3 研究區(qū)不同地下水中陽離子與氯離子關(guān)系圖Fig.3 Relationships between main cations vs. chloride of groundwater samples in the study area

Ca2+與Cl-的比值以及Mg2+與Cl-的比值隨著Cl-濃度的增加而降低，但是都遠(yuǎn)高于海水中的比值。說明在Cl-濃度較高的地下水中，Cl-與Ca2+、Mg2+的比例更接近海水，但是均高于海水。說明海水入侵一定程度上影響了地下水離子組成，但是研究區(qū)地下水總體表現(xiàn)出低TDS的特征。同時(shí)濱海區(qū)地下水中Mg2+與Cl-的比值比Ca2+與Cl-的比值更接近海水中的比值，主要是由于Ca2+為非保守離子，在地下水形成過程中，參與了多種礦物的沉淀——溶解過程和陽離子交換反應(yīng)。

Ca2+與Cl-的比值以及Mg2+與Cl-的比值隨著Cl-濃度的增加而降低，但是都遠(yuǎn)高于海水中的比值。說明在Cl-濃度較高的地下水中，Cl-與Ca2+、Mg2+的比例更接近海水，但是均高于海水。說明海水入侵一定程度上影響了地下水離子組成，但是研究區(qū)地下水總體仍表現(xiàn)出低TDS的特征。同時(shí)濱海區(qū)地下水中Mg2+與Cl-的比值比Ca2+與Cl-的比值更接近海水中的比值，主要是由于Ca2+為非保守離子，在地下水形成過程中，參與了多種礦物的沉淀-溶解過程和陽離子交換反應(yīng)。

圖4 研究區(qū)不同地下水吉布斯圖Fig.4 Gibbs maps of groundwater samples in study area

2.1.3 優(yōu)化Piper圖解法 為了識(shí)別海水入侵的影響，Gopinath等[3]使用7種不同的分類標(biāo)準(zhǔn)，將傳統(tǒng)的Piper圖法進(jìn)行了優(yōu)化。圖5為使用改良Piper三線圖對(duì)研究區(qū)地下水化學(xué)類型進(jìn)行分析。結(jié)果表明，濱海區(qū)的地下水樣品均位于海水入侵水區(qū)域，農(nóng)業(yè)活動(dòng)區(qū)的地下水樣品主要分布在輕微海水入侵水區(qū)域，但是也有Cl-濃度較低的地下水樣品落在海水入侵水區(qū)，這些地下水樣品中Na+和Cl-離子的毫克當(dāng)量百分比較高，可能是受到巖鹽溶解的影響，而非海水入侵造成的。

圖5 研究區(qū)地下水樣品優(yōu)化Piper圖Fig.5 Optimized Piper diagram of groundwater samples in the study area

2.2 耦合模型計(jì)算法

2.2.1 模型的原理和主要計(jì)算方法 本研究針對(duì)黃水河濱海地下水源地，使用Tubau等[11]提出的一種結(jié)合使用端元混合分析和混合計(jì)算來計(jì)算混合比的方法對(duì)研究區(qū)進(jìn)行地下水來源解析。該方法包括4個(gè)主要步驟：①研究區(qū)水文地質(zhì)條件調(diào)查與水質(zhì)分析；②選擇要在分析中使用的指標(biāo)；③ 確定潛在的地下水補(bǔ)給源；④使用混合計(jì)算程序計(jì)算混合比。其中步驟②和步驟③需要在端元混合分析中進(jìn)行，步驟④在混合計(jì)算程序中進(jìn)行。

2.2.1.1 端元混合分析(EMMA) 端元混合分析提供了一種利用水化學(xué)數(shù)據(jù)的特征向量和殘差分析來估計(jì)端元的數(shù)量的方法，它的目標(biāo)是確定解釋該研究區(qū)地下水補(bǔ)給來源的所需的端元的最小數(shù)目。在端元混合分析過程中，主要遵守的原則和步驟為：①對(duì)所有樣品進(jìn)行主成分分析，該步驟主要是為了消除特征值最大的指標(biāo)；②對(duì)于分析后特征向量的選取，遵循“1規(guī)則”[10,19]，即保留所有特征值大于或等于1的特征向量，數(shù)量為n；③最終選取的端元的數(shù)量為n+1；④最后將每個(gè)樣品點(diǎn)的特征值投影到特征向量所定義的空間中。根據(jù)投影圖對(duì)端元進(jìn)行初步選擇，并對(duì)落在投影定義區(qū)域之外的點(diǎn)進(jìn)行分析。

2.2.1.2 混合計(jì)算(MIX calculations)程序 混合計(jì)算是Carrera等[12]提出的一種對(duì)地質(zhì)年代學(xué)線性擬合方法的擴(kuò)展，其主要步驟是：①初始化，即假設(shè)端元濃度已知，定義初始混合比；②通過極大似然函數(shù)重新評(píng)估端元和地下水樣品的組分；③使用地下水樣品組分的濃度的期望值，利用極大似然函數(shù)獲得最終混合比；④重復(fù)步驟②和步驟③直到收斂。針對(duì)端元組分和混合比構(gòu)建的極大似然函數(shù)分別為：

(1)

式中：yp是第p個(gè)樣品中所有指標(biāo)的向量；Ap是它們的協(xié)方差矩陣；F是選擇的端元的指標(biāo)的ns×ne維矩陣；δp是混合比組成的向量。極大似然函數(shù)的具體求解方法可以參考Carrera等[12]。

圖6 EMMA特征向量組分相對(duì)貢獻(xiàn)系數(shù)圖Fig.6 Relative contribution coefficient of eigenvector components of EMMA

根據(jù)端元混合分析的特征向量投影圖(見圖7)，選取了截滲墻外咸水(19號(hào))、大氣降水(R1)以及農(nóng)業(yè)灌溉地表水(11號(hào))作為端元。

圖7 研究區(qū)地下水樣品特征向量投影圖Fig.7 The projection of the eigenvector of groundwater samples in the study area

2.2.3 不同地下水補(bǔ)給來源的混合比

2.2.3.1 混合計(jì)算輸入數(shù)據(jù) 端元混合分析得到的3個(gè)端元 (19號(hào)、R1和11號(hào))分別定義為:EM1、EM2、EM3。將2019年9月采集的樣品中選定組分的濃度作為端元濃度進(jìn)行混合計(jì)算，由于混合計(jì)算承認(rèn)端元濃度的不確定性，計(jì)算對(duì)每個(gè)選定的端元都賦予了標(biāo)準(zhǔn)偏差，標(biāo)準(zhǔn)偏差根據(jù)3次采樣的結(jié)果進(jìn)行計(jì)算得到。端元的初始濃度及輸入標(biāo)準(zhǔn)差見表2。

表2 輸入端元初始濃度及標(biāo)準(zhǔn)差Table 2 Initial parameters of MIX program

2.2.3.2 混合比計(jì)算結(jié)果 3個(gè)端元對(duì)每組地下水的貢獻(xiàn)如表3所示。對(duì)于研究區(qū)整體而言，EM3對(duì)研究區(qū)的貢獻(xiàn)最大，混合比平均值為51.17%；EM1和EM2的混合比分別為10.22%和38.60%。

表3 研究區(qū)地下水混合比Table 3 Mixing ratios of groundwater samples in the study area %

截滲墻下游的地下咸水對(duì)濱海區(qū)地下水的影響十分明顯，平均混合比為34.03%。截滲墻下游地下咸水對(duì)3號(hào)樣品點(diǎn)影響最為明顯，混合比為59.40%。大氣降水和灌溉地表水對(duì)濱海區(qū)地下水的平均混合比分別為29.70%和36.27%。表明濱海地區(qū)地下水雖然受到海水入侵的影響，但是海水入侵程度已經(jīng)逐漸減弱，大氣降水和地表水的貢獻(xiàn)增加。

對(duì)于農(nóng)業(yè)活動(dòng)區(qū)地下水，貢獻(xiàn)率最高的是灌溉地表水，平均混合比為55.75%。截滲墻下游的地下咸水對(duì)其的影響十分微弱，平均混合比僅為2.91%，其中混合比最高的是8號(hào)樣品，混合比為11.30%。同時(shí)發(fā)現(xiàn)大氣降水混合大于50%的樣品點(diǎn)，除5號(hào)樣品外，全部分布在黃水河河道兩側(cè)，說明大氣降水通過進(jìn)入河道后入滲補(bǔ)給地下水。

3 結(jié)論

本研究利用水化學(xué)和數(shù)學(xué)計(jì)算方法來探究黃水河流域地下水源地特殊水文地質(zhì)條件下的地下水補(bǔ)給源問題。主要結(jié)論為:

(1)端元混合分析與混合計(jì)算耦合的數(shù)學(xué)方法在判斷特定研究區(qū)的地下水補(bǔ)給來源和計(jì)算每個(gè)補(bǔ)給來源的混合比方面，可以得到傳統(tǒng)的源解析方法無法得到的定量計(jì)算結(jié)果，并具有更高的準(zhǔn)確性，具有良好的應(yīng)用價(jià)值。

(2)水化學(xué)初步判斷可以得出，研究區(qū)地下水鹽分形成過程中海水入侵、礦物溶解—沉淀和陽離子交換均發(fā)揮了不同程度的作用。海水入侵影響了濱海區(qū)地下水的離子組分；礦物溶解——沉淀主要是巖鹽的溶解是農(nóng)業(yè)區(qū)地下水中TDS的主要影響因素；陽離子交換則主要發(fā)生在濱海區(qū)，隨著海水入侵過程的加劇，影響了地下水中Ca2+和Na+的濃度。

(3)根據(jù)采樣點(diǎn)地理位置、水化學(xué)組分和土地利用類型，可以將研究區(qū)分為濱海區(qū)與農(nóng)業(yè)活動(dòng)區(qū)，其中截滲墻下游的地下咸水、大氣降水和灌溉的地表水對(duì)濱海區(qū)地下水的平均混合比分別為34.03%、29.70%和36.27%，對(duì)農(nóng)業(yè)活動(dòng)區(qū)地下水的平均混合比分別為6.69%、49.03%和44.28%。對(duì)于研究區(qū)整體而言，3個(gè)補(bǔ)給源的平均混合比分別為10.22%、38.60%和51.17%。