張 野

(遼寧省盤錦水文局，遼寧 盤錦 124010)

1 研究區(qū)概況

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

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

2 水文地質(zhì)

大清河流域面積較廣，在其上游地區(qū)河岸階地狹窄，面積較小，望寶山下游河谷平原逐漸開闊且望寶山水文監(jiān)測站水位監(jiān)測資料齊全，作為已知水頭邊界能夠反映出上游河流對下游地區(qū)的影響。因此，研究區(qū)模擬范圍為：南北兩側(cè)以低山丘陵為界，東側(cè)以望寶山水文監(jiān)測站為界，西側(cè)延伸至遼東灣距海岸2km處。模型的各邊界設(shè)定如下：南北兩側(cè)含水層受山前地下水徑流補給、山前沖洪積扇的側(cè)向徑流補給等多項側(cè)向補給，是已知流量邊界；東側(cè)根據(jù)望寶山水文監(jiān)測站多年水位監(jiān)測數(shù)據(jù)可定義為已知水頭邊界；西側(cè)邊界取海岸向海洋延伸2km處，可定義為已知潮汐水頭H(t)(淤泥質(zhì)海岸，水頭值恒大于等于淤泥質(zhì)海岸頂板高程)邊界和已知Cl-濃度20g/L邊界(該區(qū)監(jiān)測資料以Cl-濃度為主)；攔河閘可作為已知水頭邊界；區(qū)域內(nèi)河段在模型中作為河流處理(已知水頭和河床底板高程)；底部隔水巖層可視為零流量邊界；含水層與大氣之間的水分交換發(fā)生在潛水含水層，包括大氣降水(大氣降水量取營口氣象站測定的1991—2015年年平均變化降水量進行賦值)和蒸發(fā)(蒸發(fā)極限深度為4m，研究區(qū)地下水位埋深多大于4m，因此蒸發(fā)量忽略不計)；源匯項主要有農(nóng)業(yè)開采井(在農(nóng)耕區(qū)按照1.5口井/km2進行布設(shè)，共有90口井)和水源地開采井(下游3個水源地的所有開采井和團甸水源地的4口開采井)。

3 數(shù)值模型

數(shù)學模型，變密度水流方程為

(1)

式中：hf為等效淡水水頭,m;Kf為等效淡水滲透系數(shù),m/h;ρf為淡水的密度,kg/m3;ρ為實際液體密度,kg/m3；Z為位置點高程,m;Sf為等效淡水彈性釋(儲)水率,1/m;t為時間,h;θ為有效孔隙度;C為溶質(zhì)濃度,mg/h;ρs為單位體積源(匯)項的密度,kg/m3;qs為單位體積源(匯)項的流量,m3/s。

式(1)加上相應的初始條件和邊界條件，就構(gòu)成了描述地下水流運動的數(shù)學模型[1]。定解條件可表示為

初始條件：

H(x,y,z,0)=H0(x,y,z)

(2)

第一類邊界條件：

H(x,y,z,t)|Γ=H1(x,y,z,t)

(3)

第二類邊界條件：

(4)

式中：H0(x,y,z)為研究區(qū)各層初始水頭值；H1(x,y,z,t)為研究區(qū)各層第一類邊界上的實測水頭值；q(x,y,z,t)為研究區(qū)各層第二類邊界上的單位面積流量[2]。另外，由于不考慮溶質(zhì)運移過程中的化學反應，因此溶質(zhì)運移方程表達式為

(5)

式中：C為溶質(zhì)濃度,mg/L;D為溶質(zhì)擴散系數(shù),m2/h;Cs為單位體積源(匯)項的溶質(zhì)濃度,mg/L。

式(5)加上相應的初始條件，就構(gòu)成了地下水溶質(zhì)運移的數(shù)學模型。定解條件可表示為

初始條件：

C(x,y,z,0)=C0(x,y,z)

(6)

第一類邊界條件：

C(x,y,z,t)|Γ1=C1(x,y,z,t)

(7)

式中：C0(x,y,z)為已知的濃度分布；C1(x,y,z,t)為研究區(qū)各層已知濃度邊界上的實測濃度值。

3.1 模型離散

在水平方向上, 將該模型離散為100m×100m的矩形網(wǎng)格, 垂向上依照相應的含水層特征進行劃分。模型共分成5層, 激活邊界內(nèi)單元格，有效單元格共有9個。模型模擬的校正期是1991年11月至2015年11月，將兩年作為一個模擬區(qū)間，兩個月作為輸出時間步長。

3.2 模型識別

初始含水層水平滲透系數(shù)根據(jù)遼寧地質(zhì)勘察報告抽水試驗所得結(jié)果進行賦值，假設(shè)水平滲透系數(shù)是垂直滲透系數(shù)的10倍。該區(qū)內(nèi)共有3個水位監(jiān)測井和3個多年濃度(含鹽量)取樣點，將計算出的水位和濃度與監(jiān)測的水位和濃度進行匹配，通過調(diào)整滲透系數(shù)、彈性釋(儲)水率進行標定(即模型識別的主要參數(shù)是滲透系數(shù)和彈性釋水率)，以滿足水流場和濃度場以及水頭和Cl-溶度隨時間的分布規(guī)律[3]。反復標定至3個水位觀測井(W1、W2、W3)的水位觀測值和計算值之間的平均絕對誤差小于觀測值的10%，3個濃度觀測井(W4、W5、W6)的濃度觀測值和計算值之間的平均絕對誤差小于觀測值的10%。使得最后的水流場和濃度場分別與真實水流場濃度場相似，識別參數(shù)結(jié)果見表1，另外，水平滲透系數(shù)是垂直滲透系數(shù)的5倍。模型中剩余的參數(shù)沒有進行識別，根據(jù)已有模擬案例取經(jīng)驗值，其中孔隙度為0.3，縱向彌散度和橫向彌散度分別為500m和50m，擴散系數(shù)設(shè)定為1×10-4m2/d。

表1 研究區(qū)分區(qū)水文地質(zhì)參數(shù)識別結(jié)果

該區(qū)內(nèi)觀測井觀測數(shù)據(jù)與計算值之間的擬合曲線見圖1，2015年計算水位見圖2。從圖2中可以看出，無論是水位觀測井還是濃度觀測井的水位觀測值還是濃度觀測值，與模擬計算所得變化趨勢都是基本相同的。從圖1中可以看出，水位與濃度的計算值均在觀測值之間均勻地波動，造成這種現(xiàn)象的原因可能是模型上邊界大氣降水賦值是按照年平均降水量進行的，在一個模擬區(qū)間內(nèi)，計算水位是該時間區(qū)間內(nèi)水位變化的平均值，監(jiān)測水位是具體時間點的水位。

3.3 模型驗證

為了驗證降水量的時間不均一性是否是產(chǎn)生上述誤差的原因，選取2015年11月至2016年11月作為模型的驗證期。在模型驗證期內(nèi)降水量按照實際月降水量進行賦值，其他參數(shù)不變，將1個月作為一個模擬區(qū)間，3天作為一個輸出時間步長進行運算。將水位監(jiān)測井(W1、W2、W3)的計算水位值分別與實際水位進行對比，見圖3，通過驗證可知，該模型可以對海水入侵未來趨勢作出分析。

4 模型不確定性分析與模擬結(jié)果

4.1 模型參數(shù)不確定性分析

建模中含水層的劃分是按照實際含水層結(jié)構(gòu)和分布確定的。由于河谷平原實際含水層面積自上而下是逐層減少的，為使計算結(jié)果穩(wěn)定、收斂，假設(shè)模型的每一層含水層面積相等，而每一層含水層多出的面積部分賦予較小的滲透參數(shù)值[4]。通過多組計算對比，對計算結(jié)果產(chǎn)生的誤差很小，可以忽略。

該區(qū)降水入滲系數(shù)是按照已有研究的經(jīng)驗數(shù)據(jù)進行賦值的。由于模型校正周期時間較長，在校正過程中直接取年均降水量值進行計算[5]，而在驗證模型過程中按照月均降水量進行賦值。對比同一時間按年均降水量與按月均降水量賦值的計算結(jié)果，誤差較小，可忽略。因此，在后期對未來20年的預測分析中可按照多年降水量的平均值進行賦值。

由于該區(qū)海岸帶地勢平坦，灘涂可向海岸帶延伸約2km，因此，在建立模型邊界范圍時向海岸帶延伸了2km，并將灘涂范圍設(shè)定為潮汐水頭。對海岸帶延伸距離進行敏感性分析，分別假定海岸帶延伸距離為1km和3km。模擬結(jié)果表明，潮汐作用條件下，海岸帶延伸距離對結(jié)果的影響較小，邊界延伸距離的變化對模擬結(jié)果影響并不十分敏感。

4.2 模擬結(jié)果

該區(qū)1991—2015年的模擬結(jié)果顯示，地下水的過度開采(研究區(qū)降雨補給平均約1700萬m3/a，側(cè)向補給400萬m3/a，河流流量約2500萬m3/a，水源地和農(nóng)業(yè)井的最大開采量約5500萬m3/a，補給量小于開采量)在大清河的南側(cè)永安水源地附近形成了一個降落漏斗，見圖2，降落漏斗中心水位最低時約-16m，濃度等值線(見圖4)則由于密度差異、分層水文地質(zhì)參數(shù)不同、分層開采量不同等原因呈現(xiàn)出圖4(a)～圖4(e)中的分層差異，分別代表研究區(qū)第1～5層的濃度分布。從圖4(a)中可以看出，第1層海水入侵主要受潮汐波動的影響，攔河閘下游的大清河河道在潮汐作用下，高潮位海水沿河道上溯，出現(xiàn)高濃度入侵鋒；從圖4(b)中可以看出，第2層由于大量農(nóng)業(yè)灌溉開采地下水，第1層地下水下滲補給時，攜帶Cl-向降落漏斗中心入侵，由于北岸河谷平原面積較小，農(nóng)業(yè)井較少，大清河北岸海水入侵程度明顯低于南側(cè)。圖4(c)和圖4(d)則表明第3、第4層含水層海水入侵基本不受潮汐作用的影響，主要是由于第4層內(nèi)水源地開采井大量開采地下水，海水沿指向內(nèi)陸的水力坡降向內(nèi)陸運移。從圖4(d)和圖4(e)中可以看出，第5層與第4層濃度分布幾乎無差異，其原因可能是第5層滲透系數(shù)相對較小，該層入侵主要由第4層已入侵部分在重力作用下滲透。因此，該區(qū)發(fā)生海水入侵的主要原因有兩方面，一方面是潮汐作用下的海水沿河道上溯，并向四周補給低水位地下水所產(chǎn)生的入侵；另一方面是由于地下水的大量開采形成降落漏斗，降落漏斗中心水位遠低于海平面，產(chǎn)生指向內(nèi)陸的水力坡降，海水沿水力坡降方向入侵[5]。

設(shè)定濃度小于250mg/L的地區(qū)是未入侵區(qū)，天然無開采條件下研究區(qū)入侵面積定義為A0，某一時間節(jié)點對應的入侵面積定義為At，則某一時間節(jié)點對應的海水入侵面積為ΔAt=At-A0[6-9]。研究區(qū)2012年開始實施壓采方案對各水源地開采量進行控制，2012—2013年初步壓采(永安水源地和團甸水源地開采量減半)，2014—2015年再次壓采(各個水源地的單井開采量均不超過1000m3/d), 但由于對農(nóng)業(yè)灌溉井并未進行約束，加之降落漏斗中心水位未恢復至海平面以上，研究區(qū)整體水力梯度仍由海洋指向內(nèi)陸[10]。計算所得研究區(qū)實施水源地壓采之后逐年分層入侵面積變化情況見表2，可以看出，壓采方案實施后，250mg/L的入侵等值線仍舊向內(nèi)陸侵移，入侵面積加大。除地表受潮汐作用影響較大外，下部地層入侵速率雖有波動，但整體上呈減緩趨勢。

表2 壓采后逐年分層入侵面積變化情況

5 結(jié) 論

從以上分析可以看出監(jiān)測數(shù)據(jù)與計算數(shù)據(jù)變化趨勢一致，且計算值與觀測值之間的絕對誤差小于5%。從監(jiān)測井數(shù)據(jù)與計算數(shù)據(jù)的水位對比曲線來看，模型識別所得參數(shù)基本符合研究區(qū)實際情況，可以用來預測研究區(qū)未來海水入侵的趨勢。由此可知，采用變密度地下水流動數(shù)學模型研究區(qū)域海水入侵過程，將會是未來發(fā)展的方向。