趙小龍
(遼寧省營(yíng)口水文局,遼寧 營(yíng)口 115003)
海水入侵是一個(gè)全球性的環(huán)境問(wèn)題,沿海省市均被其困擾。海水入侵使地下淡水體咸化,進(jìn)而導(dǎo)致淡水資源減少,土壤鹽漬化,工業(yè)機(jī)械腐蝕,甚至危害人體健康,嚴(yán)重制約沿海城市的經(jīng)濟(jì)發(fā)展[1]。海水入侵的誘發(fā)因素有兩類(lèi),一類(lèi)是自然因素,包含海平面上升、干旱、潮汐等;另一類(lèi)是人為因素,包含地下水開(kāi)采、海水養(yǎng)殖等。WHO指出,若淡水混入約1%體積的海水(氯離子濃度大于250mg/L)則不能直接飲用。由于沿海城市經(jīng)濟(jì)發(fā)展對(duì)地下水的需求以及全球變暖引起的海平面上升等問(wèn)題的存在,關(guān)于海水入侵的研究目前仍是一個(gè)熱點(diǎn)問(wèn)題。
研究區(qū)位于遼東灣東側(cè)的遼寧省營(yíng)口市,大清河流域望寶山水文監(jiān)測(cè)站的下游,內(nèi)有長(zhǎng)大鐵路通過(guò),交通便利,總面積為170km2(陸地面積150km2,潮汐波動(dòng)帶20km2),東西方向約24km,南北方向2~13km,海岸線(xiàn)約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℃,無(wú)霜期為180~210天。
營(yíng)口市屬華北地臺(tái)遼東臺(tái)背斜營(yíng)口至寬甸隆起的南翼,受燕山運(yùn)動(dòng)的影響形成的千山余脈呈東—西向縱貫本地區(qū)。地形自東向西逐漸由高變低。受構(gòu)造、巖性與新構(gòu)造運(yùn)動(dòng)的影響,地貌形態(tài)由東向西呈規(guī)律性變化,即低山—高丘陵—低丘陵—濱海平原。
研究區(qū)是發(fā)源于營(yíng)口東部山間的大清河入海沖洪積形成的山間河谷平原,主要由一級(jí)階地和二級(jí)階地組成。遼寧省水文地質(zhì)大隊(duì)曾對(duì)大清河中下游河谷平原進(jìn)行地質(zhì)勘查,順河流而下,共勘探了5個(gè)垂直于河流的剖面和一個(gè)平行于河流的剖面。根據(jù)勘測(cè)所得地層剖面,概化出研究區(qū)含水介質(zhì)概念模型,并對(duì)研究區(qū)地下水含水層特征進(jìn)行分析。
區(qū)內(nèi)第四紀(jì)覆蓋層受下伏基巖面起伏影響厚度變化大,變化范圍20~65m,整體上呈現(xiàn)出由河谷向兩側(cè)丘陵逐漸變薄,北側(cè)丘陵覆蓋層由于側(cè)向山谷沖洪積作用靠山一側(cè)覆蓋層較厚。由于河流流量歷史上隨季節(jié)氣候變化,因此,由河流沖洪積形成的平原含水層性質(zhì)也具有很大差異,自上而下可以分為5層:第一層是以亞砂土、黏土為主的第四紀(jì)覆蓋層,含水層類(lèi)型為潛水含水層,厚度5~10m,滲透能力相對(duì)較弱,特別是潮汐波動(dòng)帶的淤泥質(zhì)海岸,在海潮波動(dòng)下飽和,漲潮時(shí)被淹沒(méi),退潮時(shí)水分疏干速度慢,加之此處為半日潮,退潮到下次漲潮之間淤泥質(zhì)海岸中水位變幅基本可以忽略;第二層以砂卵礫石為主,局部沉積亞黏土、亞砂土等,為主要含水層,以潛水為主,局部微承壓,該層介質(zhì)粒徑較大,孔隙度較高,滲透系數(shù)20~100m/d,自河谷向兩側(cè)逐漸減?。坏谌龑右詠喩巴翞橹?,局部含有黏土,厚度為0~8m,透水性稍弱;第四層為礫卵石,厚度約5~15m;第五層為黏土層,厚度約3~10m。
區(qū)內(nèi)主要的地表水體為大清河地表水系,河床寬度約20~300m。大清河發(fā)源于營(yíng)口市東部山區(qū),流域面積1468km2。為監(jiān)測(cè)大清河的水位、流量、流速等,1959年在大清河的中下游修建了望寶山水文站。根據(jù)望寶山水文站多年觀測(cè),大清河最大徑流量為57m3/s,最小徑流量為0.316m3/s,徑流隨季節(jié)變化明顯。為調(diào)節(jié)大清河徑流量隨季節(jié)的變化,在其上游修建了石門(mén)水庫(kù)。此外,在大清河下游修建了1座集蓄水、灌溉、擋潮于一體的攔河閘。
區(qū)內(nèi)潛水含水層的補(bǔ)給來(lái)源主要有大氣降水入滲、河流側(cè)向補(bǔ)給、山前地下水徑流補(bǔ)給、下伏碳酸鹽巖巖溶水的頂托補(bǔ)給、山前沖洪積扇的側(cè)向徑流補(bǔ)給和開(kāi)采條件下的越流補(bǔ)給。含水層接受河流側(cè)向補(bǔ)給、大氣降水入滲補(bǔ)給等多項(xiàng)補(bǔ)給后,沿地勢(shì)自東向西流,徑流速度隨著含水層厚度、透水性和地形的變化而變化[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ī)電井(井深較淺,分布在第二層)來(lái)滿(mǎn)足供水需求。自水源地開(kāi)采井和農(nóng)業(yè)用井修建以來(lái),人工開(kāi)采已經(jīng)成為地下水最主要的排泄方式。由于農(nóng)業(yè)灌溉用井和水源地開(kāi)采井對(duì)地下水大量開(kāi)采,地下水位低于河水水位,研究區(qū)地下水不再向河流排泄,而是長(zhǎng)時(shí)間受河流補(bǔ)給,由于河流補(bǔ)給速度和補(bǔ)給量有限,在水源地開(kāi)采井和農(nóng)業(yè)灌溉井廣泛分布的永安水源地附近形成降落漏斗。
營(yíng)口市海岸屬于淤泥質(zhì)(較細(xì)粒底質(zhì))海岸,海岸線(xiàn)較長(zhǎng),進(jìn)行現(xiàn)場(chǎng)調(diào)查時(shí)主要沿大清河流域開(kāi)展工作,由上游石門(mén)水庫(kù)開(kāi)始,沿途搜集相關(guān)資料,并現(xiàn)場(chǎng)觀測(cè)部分?jǐn)?shù)據(jù),直至下游入海處。大清河中游修建有望寶山水文站,此處大清河河床寬約17m,河漫灘兩側(cè)為低山丘陵。團(tuán)甸水源地在望寶山水文監(jiān)測(cè)站附近,共19口開(kāi)采井(下游4口),原水源地總開(kāi)采量為3.0萬(wàn)m3/d,經(jīng)2013年壓采后的開(kāi)采量減少至0.7萬(wàn)m3/d。沿大清河向下游方向南側(cè)有蓋州二水源地與蓋州三水源地(簡(jiǎn)稱(chēng)蓋州二三水源地),共有11口開(kāi)采井,全部位于大清河南岸,原開(kāi)采量為3.0萬(wàn)m3/d,壓采后現(xiàn)開(kāi)采量為1.2萬(wàn)~2.0萬(wàn)m3/d。
繼續(xù)沿岸踏勘可見(jiàn)化纖水源地,共8口開(kāi)采井,位于大清河流域的南岸,原開(kāi)采量為1.5萬(wàn)m3/d,壓采后現(xiàn)開(kāi)采量減少為0.5萬(wàn)m3/d。其下游修建有永安水源地,共19口開(kāi)采井,原開(kāi)采量為5.0萬(wàn)m3/d,壓采后現(xiàn)開(kāi)采量為1.5萬(wàn)m3/d。此處有微咸水。
在復(fù)雜非線(xiàn)性模型中存在較多參數(shù),使得模型計(jì)算難以進(jìn)行,因此需要分析出模型中較為敏感和重要的參數(shù)。Sobol法是一種基于方差的全局敏感性分析方法[2],它可以測(cè)試整個(gè)輸入空間的敏感性,即通過(guò)全局方法獲取非線(xiàn)性響應(yīng)和測(cè)試參數(shù)之間相互作用的影響。
3.1.1 基于水文地質(zhì)資料建立地下水?dāng)?shù)值模型
把該模型看作一個(gè)函數(shù)Y=f(X),其中,X是一個(gè)n個(gè)不確定模型輸入{X1,X2,…,Xd}的向量[3],如降雨入滲率、滲透系數(shù)、儲(chǔ)水率等,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ù)
通過(guò)對(duì)Y方差分解,敏感性指數(shù)可以劃分為一階敏感性系數(shù)和全階敏感性系數(shù)[6]。對(duì)于一階敏感性系數(shù),它主要是由Xi所引起的對(duì)輸出方差的貢獻(xiàn),所以它主要量測(cè)了Xi單獨(dú)變化的影響,可寫(xiě)為
(1)
式中:Xi為第i個(gè)參數(shù);X~i為除去Xi以外的其他參數(shù);Var為輸出的方差;E為期望;Si為參數(shù)Xi的一階敏感性系數(shù)。
全階敏感性系數(shù)用以量測(cè)Xi對(duì)輸出方差的貢獻(xiàn),包括由Xi與其他輸入變量間相互作用所造成的所有任意階方差[7-9],可寫(xiě)為
(2)
(3)
(4)
Var(Y)=Var(YA,YB)
(5)
估計(jì)值的精度取決于N。N值的選取可先順序加點(diǎn),然后計(jì)算指數(shù)直到估計(jì)值達(dá)到一定程度可接受的收斂[10-11]。
在進(jìn)行敏感性分析之前,將當(dāng)?shù)啬杲邓亢腿斯尉_(kāi)采量(農(nóng)用井Q1、永安水源地Q2、化纖水源地Q3、二三水源地Q4)作為中間值,同時(shí)將中間值減少和增加50%之后的范圍作為參數(shù)的敏感性分析分布范圍,另外給定平均海水位一個(gè)波動(dòng)范圍。第二、第四含水層各參數(shù)的中間值及分布范圍設(shè)置見(jiàn)表1、表2。
表1 含水層第二層Sobol敏感性分析參數(shù)設(shè)置
表2 含水層第四層Sobol敏感性分析參數(shù)設(shè)置
本次研究中,第二層主要選取平均海水位、農(nóng)田灌溉井開(kāi)采量、降水量這3個(gè)相關(guān)性比較大的參數(shù)進(jìn)行敏感性分析;第四層主要選取下游3個(gè)水源地各自開(kāi)采量作為參數(shù)進(jìn)行敏感性分析。
對(duì)第二、第四含水層各自3個(gè)參數(shù)進(jìn)行敏感性分析,分別以不同的觀測(cè)點(diǎn)的Cl-濃度值作為模型輸出的Y。在第二層和第四層已入侵部分均勻地選擇6個(gè)觀測(cè)點(diǎn)的Cl-濃度。通過(guò)Sobol法進(jìn)行敏感性分析時(shí),將研究區(qū)海水入侵模型運(yùn)行一年,根據(jù)6個(gè)觀測(cè)點(diǎn)的Cl-濃度可以分析得出不同位置處各參數(shù)全階敏感性系數(shù)STi,見(jiàn)圖1、圖2。
圖1 第二層O1~O6觀測(cè)點(diǎn)的全階敏感性系數(shù)
圖2 第四層O1~O6觀測(cè)點(diǎn)的全階敏感性系數(shù)
由圖1可知:平均海水位對(duì)第二層各點(diǎn)位的濃度變化均有影響,O4點(diǎn)距離大清河最近,相應(yīng)受平均海水位的影響最大。由圖2可知:農(nóng)用井開(kāi)采量對(duì)O4~O6點(diǎn)影響較小,對(duì)O1~O3點(diǎn)影響較大。由于研究區(qū)大清河北岸平原面積小,農(nóng)用井?dāng)?shù)量少,大清河南岸河谷平原較平坦寬闊,農(nóng)用井分布多,總開(kāi)采量較大,兩岸之間又有大清河相隔,一定程度上干擾了水力聯(lián)系。因此,南岸觀測(cè)點(diǎn)相對(duì)敏感。
本文應(yīng)用Sobol法建立地下水?dāng)?shù)值模型,對(duì)遼寧南部大清河海水入侵敏感性進(jìn)行分析。通過(guò)模型結(jié)果分析可知:該區(qū)發(fā)生海水入侵的主要原因來(lái)自?xún)煞矫?,一方面是潮汐作用下的海水沿河道上溯,并向四周補(bǔ)給低水位地下水所產(chǎn)生的入侵;另一方面是由于地下水的大量開(kāi)采形成降落漏斗,從而導(dǎo)致降落漏斗中心水位遠(yuǎn)低于海平面的平均水力坡降,海水沿水力坡降方向入侵。水源地進(jìn)行壓釆后,農(nóng)業(yè)用水成為研究區(qū)地下水的主要排泄方式,為加速海水入侵回退速率,應(yīng)對(duì)研究區(qū)農(nóng)業(yè)灌溉用井的開(kāi)采量加以控制??傮w上,研究區(qū)進(jìn)行壓釆后地下水位抬升,原來(lái)指向內(nèi)陸的水力梯度重新指向海洋,海水入侵回退,回退速率由水力梯度大小決定,回退比例與回退時(shí)間呈線(xiàn)性相關(guān)。大量人工開(kāi)采是該區(qū)產(chǎn)生海水入侵的主要原因,控制開(kāi)采量是回退海水入侵最直接的辦法。研究區(qū)淤泥質(zhì)海岸由于其弱透水性的特性,不利于產(chǎn)生有利海侵回退的水力梯度。若想海水入侵快速回退,應(yīng)在近海岸處進(jìn)行處理(如大量抽水使水位下降等),使得由內(nèi)陸指向海洋的水力梯度變大。