基于GIS的地下水模型邊界地下水流估算方法

于國寶

(遼寧省朝陽水文局,遼寧 朝陽 122000)

1 概述

地下水和河流系統(tǒng)的生態(tài)完整性經(jīng)常受到人類活動的威脅[1]。社會也面臨日益嚴(yán)重的水資源管理問題[2]。通常地下水建模是確定地下水管理策略的重要工具，然而，為了開發(fā)一個可靠的地下水模型，研究人員需要克服許多方法上的挑戰(zhàn)[3]。本研究論文集中于一種基于地理信息系統(tǒng)(GIS)的新方法，用于在地下水流模擬中估算研究區(qū)域邊界上的地下水流。

作為開發(fā)地下水流動模型的第一步，需要確定研究區(qū)域和模型邊界。模型邊界是研究區(qū)域和周圍環(huán)境之間的界面[4]。阿慧娟等[5]證明了分配正確邊界條件的重要性，并表明邊界位置的不準(zhǔn)確性會導(dǎo)致校準(zhǔn)過程中引入不必要的不均勻性。同時河流、湖泊和濕地與地下水的相互作用受水體相對于地下水流動系統(tǒng)的位置、河床的地質(zhì)特征及其氣候環(huán)境的控制[6]。陳張建[7]全面總結(jié)了河流和地下水之間的相互作用。并認(rèn)為在降雨量少的情況下，許多溪流中的基流大部分時間是來自于地下水的輸入。

GIS是為任何地下水流和污染物遷移建模工作開發(fā)地下水模型的重要工具。GIS提供數(shù)據(jù)管理和空間分析能力，可用于地下水建模。它提供數(shù)據(jù)插值、系統(tǒng)的模型參數(shù)分配、空間統(tǒng)計生成和模型結(jié)果的可視化顯示[8]。唐敏等[9]設(shè)計了一個GIS數(shù)據(jù)庫，研究了地下水脆弱性分析和水文地質(zhì)建模。GIS在地下水流系統(tǒng)的概念化、特征化和數(shù)值模擬中也具有廣泛應(yīng)用。龐園等[10]應(yīng)用GIS技術(shù)對堅硬巖石地形的地下水補(bǔ)給進(jìn)行定量建模。而且GIS還提供了許多工具，如創(chuàng)建剖面圖的工具、提取某一點(diǎn)的值、創(chuàng)建緩沖區(qū)的工具、插值工具等。因此，它可以和達(dá)西定律一起用來估算通過給定邊界的平均流量。

2 研究區(qū)

本研究選取遼寧省渾河的某支流作為研究對象。研究區(qū)的主要含水地層是沖積層，約占總面積的60%，其次是花崗巖和千枚巖。研究區(qū)的地下水位深度從平均海平面(MSL)的136～380m不等。MSL的地面高度從147～446m不等。該地區(qū)的平均降雨量為452mm。大部分降雨發(fā)生在6月至9月的季風(fēng)季節(jié)。其余時間被認(rèn)為是非季風(fēng)季節(jié)。本文從水文、水文地質(zhì)、地下水流方向和形態(tài)特征等方面對該區(qū)域進(jìn)行了深入研究，并在此基礎(chǔ)上確定了地下水流的邊界條件。在研究區(qū)域的東側(cè)，是傾向NNE-SSW的山脈?？梢詫⑵湟暈橐粋€重要的無流量邊界。在西側(cè)剩余的小部分以封閉研究區(qū)域的邊界。研究區(qū)域如圖1所示。選擇了研究區(qū)域內(nèi)及其周圍的108個探井，如圖2所示，根據(jù)這些探井的歷史水位數(shù)據(jù)，用逆距離加權(quán)(IDW)插值法計算了1999—2003年季風(fēng)和非季風(fēng)季節(jié)的水位面。插值時選擇360m×360m的像元尺寸作為輸出圖像。圖2還顯示了2003年季風(fēng)前地下水位的分布圖像。

圖1 研究區(qū)域

圖2 測點(diǎn)位置

根據(jù)獲得的28個探井記錄，繪制了研究區(qū)的水文地層圖。建造了3個沖積層、花崗巖和千枚巖的水文地層，并計算出所有三層的底部高程。ASTER高程值用于開發(fā)數(shù)字地形模型。研究區(qū)部分探井的詳細(xì)信息以及所有三層的井底高程見表1。利用鉆孔測井信息和IDW插值工具，為沖積層、花崗巖和千枚巖的水文地層底部高程準(zhǔn)備表層。還計算了監(jiān)測井所有季節(jié)的平均水位深度，并通過插值將其近似轉(zhuǎn)換為1個層。

表1 部分探井?dāng)?shù)據(jù)

3 研究方法與結(jié)果討論

本文采用了以下方法對邊界進(jìn)行計算。

(1)在研究區(qū)的邊界上創(chuàng)建了一個5公里的緩沖區(qū)。邊界被分成22個橫截面，如圖3所示。西側(cè)的三個橫截面編號為F1、F2和F3。上側(cè)的部分從U1到U10順時針編號，底端的部分從L1到L9順時針編號。圖3顯示了其中一些部分的編號。在遇到邊界急劇彎曲的地方，應(yīng)采用更多的橫截面。

圖3 橫截面示意圖

(2)利用三維分析工具，從1999—2003年的不同季風(fēng)前和季風(fēng)后水面剖面中提取了沿每個橫截面區(qū)域大約30—35個均勻分布點(diǎn)的地下水位。

(3)繪制這些點(diǎn)，并將四階多項(xiàng)式擬合到每條曲線，以確定精確邊界點(diǎn)的梯度，如圖4所示。觀察到四階多項(xiàng)式的適合度在所有情況下都非常好。2002年后季風(fēng)的U4截面擬合相關(guān)性為0.998。從四階方程出發(fā)，通過簡單的微分，可以更容易地得到任意點(diǎn)的梯度(dy/dx)。因此，在每個這樣的輪廓的邊界點(diǎn)計算梯度。

圖4 2002年季風(fēng)后沿U4斷面的水位變化

(4)盡管不同剖面之間的梯度值不同，但觀察到季風(fēng)季節(jié)前后和不同年份之間給定剖面的梯度計算值沒有顯著變化。典型線段幾個橫截面不同年份的季風(fēng)和非季風(fēng)季節(jié)的梯度計算值見表2。流入研究區(qū)域邊界的梯度值保持為正，流出為負(fù)?？梢悦黠@的發(fā)現(xiàn)，在季風(fēng)前后各年份和各位置地下水位并沒有很明顯的變化，特別是對于L9而言，水位梯度變化值小于0.00001。

(5)因此，可以通過取所有季節(jié)所有梯度計算

表2 典型線段的梯度計算值

值的平均值來計算每個橫截面的梯度平均值。然后對所有橫截面獲得的邊界段的平均梯度值進(jìn)行分析，并識別和合并具有幾乎相同梯度值的相鄰邊界段。通過這種方法，整個邊界被重新分組為8個邊界段，并計算梯度的平均值。對于這8個邊界段，基于單個線段的平均值。這些分段和梯度平均值如圖5所示。東部邊界由于靠近山脈，因此顯示值為“0”(無流量邊界)。

(6)在每個邊界段選擇5個點(diǎn)，并從相應(yīng)的層中提取這些點(diǎn)的地面高度、水文地層底部高度和平均水深值。因此，在8個邊界段上選擇的40個邊界點(diǎn)的細(xì)節(jié)制成表格。然后用這些值計算三層含水層系統(tǒng)在每個點(diǎn)的等效橫向滲透系數(shù)，然后用達(dá)西定律計算各點(diǎn)邊界每米長度的凈流量，再計算每個區(qū)段的平均流量值，并乘以區(qū)段長度，計算每個區(qū)段的總凈流量，見表3，大部分的梯度值很低。梯度最大值(0.0024)發(fā)生在在4號段。用達(dá)西定律計算得到的流量結(jié)果表明單位長度下的凈流量比從-0.5040到0.9778之間變化，在第2段流入量最大，為35725m3/s。在第6段流出量最大，為29000m3/s。

表3 每個邊界段的凈流量

4 結(jié)論

本文在無法確定研究區(qū)域附近無流動邊界的情況下，構(gòu)建了一種基于地理信息系統(tǒng)計算研究區(qū)域邊界的凈流量的方法。并認(rèn)為地下水位梯度在不同季節(jié)和不同年份之間不會發(fā)生顯著變化，此假設(shè)也在地下水模型中得到驗(yàn)證。在本研究中，研究區(qū)域邊界分為9個區(qū)段，并計算每個邊界區(qū)段的流量梯度值。通常，梯度值很低。梯度最大值(0.0024)發(fā)生在在4號段。然后使用達(dá)西定律來計算穿過這些邊界的流量。結(jié)果表明：在第2段流入量最大，為35725m3/s；在第6段流出量最大，為29000m3/s。本文提供了一種計算邊界平均凈流量的方法，以期對地下水污染研究和地下水模型建立提供一定的參考。

圖5 不同邊界段的最終水位梯度