基于Sobol法對遼寧南部海水入侵敏感性分析探究

趙小龍

(遼寧省營口水文局，遼寧 營口 115003)

海水入侵是一個(gè)全球性的環(huán)境問題，沿海省市均被其困擾。海水入侵使地下淡水體咸化，進(jìn)而導(dǎo)致淡水資源減少，土壤鹽漬化，工業(yè)機(jī)械腐蝕，甚至危害人體健康，嚴(yán)重制約沿海城市的經(jīng)濟(jì)發(fā)展[1]。海水入侵的誘發(fā)因素有兩類，一類是自然因素，包含海平面上升、干旱、潮汐等；另一類是人為因素，包含地下水開采、海水養(yǎng)殖等。WHO指出，若淡水混入約1%體積的海水(氯離子濃度大于250mg/L)則不能直接飲用。由于沿海城市經(jīng)濟(jì)發(fā)展對地下水的需求以及全球變暖引起的海平面上升等問題的存在，關(guān)于海水入侵的研究目前仍是一個(gè)熱點(diǎn)問題。

1 研究區(qū)概述

1.1 自然地理

研究區(qū)位于遼東灣東側(cè)的遼寧省營口市,大清河流域望寶山水文監(jiān)測站的下游，內(nèi)有長大鐵路通過，交通便利，總面積為170km2(陸地面積150km2，潮汐波動帶20km2)，東西方向約24km，南北方向2～13km，海岸線約6km。屬大陸性季風(fēng)氣候，四季分明，雨熱同季。1991—2016年多年平均降水量為600～800mm，年平均日照時(shí)數(shù)為2600～2880h，多年平均蒸發(fā)量1000～1200mm，蒸發(fā)極限深度為4m。多年平均氣溫9～10℃，冬季1月氣溫最低，平均為-9～-10℃，夏季7月氣溫最高，平均為24.5～25.0℃，無霜期為180～210天。

1.2 地質(zhì)條件

營口市屬華北地臺遼東臺背斜營口至寬甸隆起的南翼，受燕山運(yùn)動的影響形成的千山余脈呈東—西向縱貫本地區(qū)。地形自東向西逐漸由高變低。受構(gòu)造、巖性與新構(gòu)造運(yùn)動的影響，地貌形態(tài)由東向西呈規(guī)律性變化，即低山—高丘陵—低丘陵—濱海平原。

研究區(qū)是發(fā)源于營口東部山間的大清河入海沖洪積形成的山間河谷平原，主要由一級階地和二級階地組成。遼寧省水文地質(zhì)大隊(duì)曾對大清河中下游河谷平原進(jìn)行地質(zhì)勘查，順河流而下，共勘探了5個(gè)垂直于河流的剖面和一個(gè)平行于河流的剖面。根據(jù)勘測所得地層剖面，概化出研究區(qū)含水介質(zhì)概念模型,并對研究區(qū)地下水含水層特征進(jìn)行分析。

區(qū)內(nèi)第四紀(jì)覆蓋層受下伏基巖面起伏影響厚度變化大，變化范圍20～65m，整體上呈現(xiàn)出由河谷向兩側(cè)丘陵逐漸變薄，北側(cè)丘陵覆蓋層由于側(cè)向山谷沖洪積作用靠山一側(cè)覆蓋層較厚。由于河流流量歷史上隨季節(jié)氣候變化，因此，由河流沖洪積形成的平原含水層性質(zhì)也具有很大差異，自上而下可以分為5層：第一層是以亞砂土、黏土為主的第四紀(jì)覆蓋層，含水層類型為潛水含水層，厚度5～10m，滲透能力相對較弱，特別是潮汐波動帶的淤泥質(zhì)海岸，在海潮波動下飽和，漲潮時(shí)被淹沒，退潮時(shí)水分疏干速度慢，加之此處為半日潮，退潮到下次漲潮之間淤泥質(zhì)海岸中水位變幅基本可以忽略；第二層以砂卵礫石為主，局部沉積亞黏土、亞砂土等，為主要含水層，以潛水為主，局部微承壓，該層介質(zhì)粒徑較大，孔隙度較高，滲透系數(shù)20～100m/d，自河谷向兩側(cè)逐漸減小；第三層以亞砂土為主，局部含有黏土，厚度為0～8m，透水性稍弱；第四層為礫卵石，厚度約5～15m；第五層為黏土層，厚度約3～10m。

1.3 水文地質(zhì)條件

區(qū)內(nèi)主要的地表水體為大清河地表水系，河床寬度約20～300m。大清河發(fā)源于營口市東部山區(qū)，流域面積1468km2。為監(jiān)測大清河的水位、流量、流速等，1959年在大清河的中下游修建了望寶山水文站。根據(jù)望寶山水文站多年觀測，大清河最大徑流量為57m3/s，最小徑流量為0.316m3/s，徑流隨季節(jié)變化明顯。為調(diào)節(jié)大清河徑流量隨季節(jié)的變化，在其上游修建了石門水庫。此外，在大清河下游修建了1座集蓄水、灌溉、擋潮于一體的攔河閘。

區(qū)內(nèi)潛水含水層的補(bǔ)給來源主要有大氣降水入滲、河流側(cè)向補(bǔ)給、山前地下水徑流補(bǔ)給、下伏碳酸鹽巖巖溶水的頂托補(bǔ)給、山前沖洪積扇的側(cè)向徑流補(bǔ)給和開采條件下的越流補(bǔ)給。含水層接受河流側(cè)向補(bǔ)給、大氣降水入滲補(bǔ)給等多項(xiàng)補(bǔ)給后，沿地勢自東向西流，徑流速度隨著含水層厚度、透水性和地形的變化而變化[1]，到濱海地帶含水介質(zhì)的顆粒逐漸變小，透水性相應(yīng)減小，地下水徑流也隨之減緩。天然狀態(tài)下，大清河流域在上游河段地下水向河流排泄，在下游河段枯水期受地下水補(bǔ)給，豐水期河水補(bǔ)給地下水。隨著工農(nóng)業(yè)發(fā)展，用水需求增加，大清河流域先后修建了4個(gè)水源地(井深較深，分布在第四層中)，自上游至下游依次是團(tuán)甸水源地(19口井，其中研究區(qū)涵蓋4口)、化纖水源地(8口井)、蓋州二三水源地(13口井)、永安水源地(18口井)，以及大量農(nóng)業(yè)機(jī)電井(井深較淺，分布在第二層)來滿足供水需求。自水源地開采井和農(nóng)業(yè)用井修建以來，人工開采已經(jīng)成為地下水最主要的排泄方式。由于農(nóng)業(yè)灌溉用井和水源地開采井對地下水大量開采，地下水位低于河水水位，研究區(qū)地下水不再向河流排泄，而是長時(shí)間受河流補(bǔ)給，由于河流補(bǔ)給速度和補(bǔ)給量有限，在水源地開采井和農(nóng)業(yè)灌溉井廣泛分布的永安水源地附近形成降落漏斗。

2 典型區(qū)選取

營口市海岸屬于淤泥質(zhì)(較細(xì)粒底質(zhì))海岸，海岸線較長，進(jìn)行現(xiàn)場調(diào)查時(shí)主要沿大清河流域開展工作，由上游石門水庫開始，沿途搜集相關(guān)資料，并現(xiàn)場觀測部分?jǐn)?shù)據(jù)，直至下游入海處。大清河中游修建有望寶山水文站，此處大清河河床寬約17m，河漫灘兩側(cè)為低山丘陵。團(tuán)甸水源地在望寶山水文監(jiān)測站附近，共19口開采井(下游4口)，原水源地總開采量為3.0萬m3/d，經(jīng)2013年壓采后的開采量減少至0.7萬m3/d。沿大清河向下游方向南側(cè)有蓋州二水源地與蓋州三水源地(簡稱蓋州二三水源地)，共有11口開采井，全部位于大清河南岸，原開采量為3.0萬m3/d，壓采后現(xiàn)開采量為1.2萬～2.0萬m3/d。

繼續(xù)沿岸踏勘可見化纖水源地，共8口開采井，位于大清河流域的南岸，原開采量為1.5萬m3/d，壓采后現(xiàn)開采量減少為0.5萬m3/d。其下游修建有永安水源地，共19口開采井，原開采量為5.0萬m3/d，壓采后現(xiàn)開采量為1.5萬m3/d。此處有微咸水。

3 敏感性分析

3.1 分析方法

在復(fù)雜非線性模型中存在較多參數(shù)，使得模型計(jì)算難以進(jìn)行，因此需要分析出模型中較為敏感和重要的參數(shù)。Sobol法是一種基于方差的全局敏感性分析方法[2]，它可以測試整個(gè)輸入空間的敏感性，即通過全局方法獲取非線性響應(yīng)和測試參數(shù)之間相互作用的影響。

3.1.1 基于水文地質(zhì)資料建立地下水?dāng)?shù)值模型

把該模型看作一個(gè)函數(shù)Y=f(X)，其中，X是一個(gè)n個(gè)不確定模型輸入{X1，X2，…，Xd}的向量[3]，如降雨入滲率、滲透系數(shù)、儲水率等，Y是一個(gè)單值模型輸出水頭、濃度。同時(shí)得到X的先驗(yàn)分布。

3.1.2 產(chǎn)生參數(shù)序列

利用蒙特卡羅方法中的拉丁超立方抽樣算法，根據(jù)模型參數(shù)先驗(yàn)信息，抽樣得到一個(gè)N×2d維的矩陣[AB]，即每行是一個(gè)2d維的采樣點(diǎn)。

3.1.3 參數(shù)點(diǎn)構(gòu)造

3.1.4 模型計(jì)算

3.1.5 計(jì)算敏感性系數(shù)

通過對Y方差分解，敏感性指數(shù)可以劃分為一階敏感性系數(shù)和全階敏感性系數(shù)[6]。對于一階敏感性系數(shù)，它主要是由Xi所引起的對輸出方差的貢獻(xiàn)，所以它主要量測了Xi單獨(dú)變化的影響，可寫為

(1)

式中：Xi為第i個(gè)參數(shù)；X～i為除去Xi以外的其他參數(shù)；Var為輸出的方差；E為期望；Si為參數(shù)Xi的一階敏感性系數(shù)。

全階敏感性系數(shù)用以量測Xi對輸出方差的貢獻(xiàn)，包括由Xi與其他輸入變量間相互作用所造成的所有任意階方差[7-9]，可寫為

(2)

(3)

(4)

Var(Y)=Var(YA,YB)

(5)

估計(jì)值的精度取決于N。N值的選取可先順序加點(diǎn)，然后計(jì)算指數(shù)直到估計(jì)值達(dá)到一定程度可接受的收斂[10-11]。

3.2 敏感性分析過程和結(jié)果

在進(jìn)行敏感性分析之前，將當(dāng)?shù)啬杲邓亢腿斯尉_采量(農(nóng)用井Q1、永安水源地Q2、化纖水源地Q3、二三水源地Q4)作為中間值，同時(shí)將中間值減少和增加50%之后的范圍作為參數(shù)的敏感性分析分布范圍，另外給定平均海水位一個(gè)波動范圍。第二、第四含水層各參數(shù)的中間值及分布范圍設(shè)置見表1、表2。

表1 含水層第二層Sobol敏感性分析參數(shù)設(shè)置

表2 含水層第四層Sobol敏感性分析參數(shù)設(shè)置

本次研究中，第二層主要選取平均海水位、農(nóng)田灌溉井開采量、降水量這3個(gè)相關(guān)性比較大的參數(shù)進(jìn)行敏感性分析；第四層主要選取下游3個(gè)水源地各自開采量作為參數(shù)進(jìn)行敏感性分析。

對第二、第四含水層各自3個(gè)參數(shù)進(jìn)行敏感性分析，分別以不同的觀測點(diǎn)的Cl-濃度值作為模型輸出的Y。在第二層和第四層已入侵部分均勻地選擇6個(gè)觀測點(diǎn)的Cl-濃度。通過Sobol法進(jìn)行敏感性分析時(shí)，將研究區(qū)海水入侵模型運(yùn)行一年，根據(jù)6個(gè)觀測點(diǎn)的Cl-濃度可以分析得出不同位置處各參數(shù)全階敏感性系數(shù)STi，見圖1、圖2。

圖1 第二層O1～O6觀測點(diǎn)的全階敏感性系數(shù)

圖2 第四層O1～O6觀測點(diǎn)的全階敏感性系數(shù)

由圖1可知：平均海水位對第二層各點(diǎn)位的濃度變化均有影響，O4點(diǎn)距離大清河最近，相應(yīng)受平均海水位的影響最大。由圖2可知：農(nóng)用井開采量對O4～O6點(diǎn)影響較小，對O1～O3點(diǎn)影響較大。由于研究區(qū)大清河北岸平原面積小，農(nóng)用井?dāng)?shù)量少，大清河南岸河谷平原較平坦寬闊，農(nóng)用井分布多，總開采量較大，兩岸之間又有大清河相隔，一定程度上干擾了水力聯(lián)系。因此，南岸觀測點(diǎn)相對敏感。

4 結(jié) 語

本文應(yīng)用Sobol法建立地下水?dāng)?shù)值模型，對遼寧南部大清河海水入侵敏感性進(jìn)行分析。通過模型結(jié)果分析可知：該區(qū)發(fā)生海水入侵的主要原因來自兩方面，一方面是潮汐作用下的海水沿河道上溯，并向四周補(bǔ)給低水位地下水所產(chǎn)生的入侵；另一方面是由于地下水的大量開采形成降落漏斗，從而導(dǎo)致降落漏斗中心水位遠(yuǎn)低于海平面的平均水力坡降，海水沿水力坡降方向入侵。水源地進(jìn)行壓釆后，農(nóng)業(yè)用水成為研究區(qū)地下水的主要排泄方式，為加速海水入侵回退速率，應(yīng)對研究區(qū)農(nóng)業(yè)灌溉用井的開采量加以控制?？傮w上，研究區(qū)進(jìn)行壓釆后地下水位抬升，原來指向內(nèi)陸的水力梯度重新指向海洋，海水入侵回退，回退速率由水力梯度大小決定，回退比例與回退時(shí)間呈線性相關(guān)。大量人工開采是該區(qū)產(chǎn)生海水入侵的主要原因，控制開采量是回退海水入侵最直接的辦法。研究區(qū)淤泥質(zhì)海岸由于其弱透水性的特性，不利于產(chǎn)生有利海侵回退的水力梯度。若想海水入侵快速回退，應(yīng)在近海岸處進(jìn)行處理(如大量抽水使水位下降等)，使得由內(nèi)陸指向海洋的水力梯度變大。