基于變密度地下水流模型分析預測沿海地區海水入侵

張 野

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

1 研究區概況

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

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

2 水文地質

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

3 數值模型

數學模型，變密度水流方程為

(1)

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

式(1)加上相應的初始條件和邊界條件，就構成了描述地下水流運動的數學模型[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)為研究區各層初始水頭值；H1(x,y,z,t)為研究區各層第一類邊界上的實測水頭值；q(x,y,z,t)為研究區各層第二類邊界上的單位面積流量[2]。另外，由于不考慮溶質運移過程中的化學反應，因此溶質運移方程表達式為

(5)

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

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

初始條件：

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)為研究區各層已知濃度邊界上的實測濃度值。

3.1 模型離散

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

3.2 模型識別

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

表1 研究區分區水文地質參數識別結果

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

3.3 模型驗證

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

4 模型不確定性分析與模擬結果

4.1 模型參數不確定性分析

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

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

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

4.2 模擬結果

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

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

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

5 結 論

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