宋壽鵬,楊 濤,康衛(wèi)東

(1.同濟(jì)大學(xué)地下建筑與工程系,上海 200092;2.石家莊經(jīng)濟(jì)學(xué)院,石家莊050031;3.西北大學(xué) 地質(zhì)學(xué)系,西安710029)

靖邊縣地處陜北黃土高原北部與毛烏素沙漠南緣過渡地帶[1],境內(nèi)石油、天然氣礦產(chǎn)資源十分豐富,是國家能源重化工榆林基地的重要組成部分。四柏樹水源地和煉油廠水源地是縣內(nèi)主要的供水水源,且兩者處于同一水文地質(zhì)單元,地下水聯(lián)系比較緊密。在兩水源地已開展的地水文地質(zhì)研究工作中,往往只研究單一水源地,而忽略兩水源地間的相互影響[2]。

本次工作依托前人的水文地質(zhì)成果[2],建立聯(lián)合兩水源地的區(qū)域地下水流模型,模擬不同開采條件,以期查明兩水源地間的相互聯(lián)系與影響,這對區(qū)域水資源的開發(fā)及生態(tài)環(huán)境的保護(hù)具有重要的意義。

1 研究區(qū)概況

靖邊縣境內(nèi)的四柏樹水源地和煉油廠水源地分別位于蘆河的北南兩側(cè),相距約6 km,屬干旱-半干旱內(nèi)陸性季風(fēng)氣候區(qū),多年平均降水量396.1 mm,降水量年際間呈豐、平、枯交替的周期性變化規(guī)律,6-9月降水約占全年降水量的71%,且多以暴雨形式降落。

研究區(qū)地處祁-呂-賀山字型構(gòu)造體系伊陜盾地與新華夏系第三沉降帶-陜甘寧盆地的復(fù)合部位的北部,為一向西緩傾的單斜構(gòu)造,地層有老到新依次為：侏羅系中統(tǒng)安定組(J2a)、白堊系下統(tǒng)洛河組(K1L)、環(huán)河組(K1h)、第三系上新統(tǒng)泥巖(N2)和第四系地層(Q)。

2 三維地下水流數(shù)值模擬

2.1 水文地質(zhì)概念模型

本次模擬計(jì)算區(qū)面積154 km2;平面上包含兩水源地分布區(qū);垂向上為地面以下侏羅系頂面以上之間的部分,厚約400～500 m。

現(xiàn)狀條件下,南邊界和西邊界為補(bǔ)給邊界,北邊界和東邊界為排泄邊界;出于安全供水考慮,計(jì)算時(shí)對各含水層的南邊界處理為定流量邊界,西邊界徑流補(bǔ)給量較小,計(jì)算時(shí)大部分地段處理為零流量邊界,北邊界和東邊界處理為變流量邊界。

模擬區(qū)含水層結(jié)構(gòu)可分為5層：第1層為第四系潛水含水層,北部的砂層為透水性較好的非均質(zhì)、各向同性含水層,南部的黃土層為透水性較差的均質(zhì)、水平與垂直各向異性的含水層;第2層為弱透水的粉質(zhì)黏土層和黏性土層,視為均質(zhì)、各向同性介質(zhì);第3層為白堊系環(huán)河組砂泥巖承壓水含水層,第4、5層是按開采層與非開采層人為劃分的白堊系洛河組厚層交錯(cuò)砂巖承壓水含水層。第3層與第4,5層之間水力聯(lián)系密切,實(shí)際上就是一個(gè)垂向上呈現(xiàn)非均質(zhì)的雙層結(jié)構(gòu)承壓水含水層;第1層與第3,4層之間通過透水性較弱的第2層發(fā)生水力聯(lián)系,第2層是上部潛水與下部承壓水發(fā)生水量交換的必經(jīng)“通道”。

模擬區(qū)南、西邊界接受地下水側(cè)向徑流補(bǔ)給,潛水含水層還接受降水入滲補(bǔ)給,北部存在潛水蒸發(fā)排泄與灌溉水入滲補(bǔ)給。模擬區(qū)的蘆河在上段補(bǔ)給地下水,在下段排泄地下水;因河床底部存在一層淤泥質(zhì)粉質(zhì)黏土層,故上段河水下滲補(bǔ)給地下水的水量較小(約900 m3/d),從水源地安全考慮可忽略不計(jì)。

2.2 數(shù)學(xué)模型

綜上所述,將模擬區(qū)地下水流系統(tǒng)可用下面偏微分方程的定解問題來描述[3]：

式中：H——地下水位標(biāo)高(m);K——滲透系數(shù)(m/d);μ——給水度;Ss——彈性釋水率 ;x,y,z——坐標(biāo)變量(m);t——時(shí)間變量(d);ε——潛水面垂向交換量(入為正、出為負(fù))(m3/d?m2);h0——地下水初始水位標(biāo)高(m);q——第2類邊界流量(m3/d?m2);q0——第3類邊界處初始時(shí)刻流量(m3/d?m2);α——第3類邊界系數(shù)(1/d);Qi——第 i眼開采井開采量(m3/d);r,θ,z′——輔助柱坐標(biāo)變量 ;r′,θ′,φ′——輔助球坐標(biāo)變量;n——邊界外法線方向;Γ2——二類邊界;Γ3——三類邊界 ;wi——第 i眼開采井的井點(diǎn)位置;Ω——模擬區(qū)范圍。

2.3 數(shù)值模型的建立

2.3.1 網(wǎng)格剖分及參數(shù)分區(qū) 本次研究選用Processing Modflow Pro集成軟件系統(tǒng)中的Modflow模塊進(jìn)行數(shù)值模擬[4-5]。模型區(qū)潛水含水層滲透系數(shù)和給水度有3個(gè)分區(qū),承壓含水層滲透系數(shù)和貯水率有1個(gè)分區(qū)(圖1)。各分區(qū)參數(shù)的初值按抽水試驗(yàn)成果并結(jié)合前人資料給出[2]。潛水和承壓含水層的初始流場由抽水試驗(yàn)前第四系含水層及白堊系含水層各自的井孔水位統(tǒng)測資料獲得,采用數(shù)值模擬方法獲得了天然狀態(tài)下的地下水流場(圖2)。

2.3.2 邊界流量 模擬區(qū)南、西邊界為二類定流量邊界,邊界流量取現(xiàn)狀地下水徑流量。北、東邊界為變流量邊界(即3類邊界),邊界流量為邊界初始流量加上邊界的增減量;邊界初始流量取現(xiàn)狀地下水徑流量,未來邊界流量的增減量按與邊界處水位升降值成正比計(jì)算,比例系數(shù)為現(xiàn)狀地下水徑流量除以潛水含水層厚度或承壓水頂板以上的水頭高度。

3.3.3 模型的識別與檢驗(yàn) 模型校驗(yàn)和識別利用大型抽水試驗(yàn)資料[2],利用抽水階段及恢復(fù)水位階段的試驗(yàn)數(shù)據(jù)用作模型檢驗(yàn)[6-7]。本次群孔抽水試驗(yàn)有7個(gè)觀測孔的實(shí)測資料可用作調(diào)參依據(jù),經(jīng)反復(fù)調(diào)參[8],典型觀測孔的實(shí)測水位曲線與計(jì)算水位曲線擬合誤差絕對值絕大多數(shù)小于0.5 m,由此可見,各觀測孔的水位擬合效果是較好的,從而得到水文地質(zhì)參數(shù)分區(qū)結(jié)果(表1)。模型識別和檢驗(yàn)結(jié)果證明所建立的數(shù)學(xué)模型、邊界條件、水文地質(zhì)參數(shù)和源匯項(xiàng)的確定都是合理的,基本反映了地下水流系統(tǒng)和含水層系統(tǒng)特征,故可利用模型進(jìn)行水源地開采的地下水位預(yù)報(bào)。

圖中坐標(biāo)為6度帶投影公里網(wǎng)圖1 水文地質(zhì)參數(shù)分區(qū)圖

3 數(shù)值評價(jià)結(jié)果分析

3.1 數(shù)值評價(jià)思路介紹

依據(jù)所建模型,預(yù)測兩水源地現(xiàn)狀條件及不同增采條件下開采30 a后的結(jié)果,評價(jià)兩水源地增采對水源地自身及區(qū)域水環(huán)境的影響。

表1 參數(shù)識別結(jié)果

現(xiàn)狀條件下,模擬區(qū)總開采量為農(nóng)灌開采量(0.98萬m3/d)和水源地開采量(四柏樹水源地為0.66萬m3/d,煉油廠水源地為0.54萬m3/d)之和。與現(xiàn)狀條件相比,水源地增采條件下,農(nóng)灌開采量均保持不變;開采井?dāng)?shù)目略有增加,其中,四柏樹水源地開采井由11眼增至14眼,煉油廠水源地由9眼增至12眼,兩水源地開采井布局見圖3。

圖中坐標(biāo)為6度帶投影公里網(wǎng)圖2 模型地下水初始流場圖

圖中坐標(biāo)為6度帶投影公里網(wǎng)圖3 水源地已建與設(shè)計(jì)開采井布局圖

3.2 四柏樹水源地單獨(dú)增采

(1)地下水位時(shí)空變化趨勢。由圖4可知：地下水開采2 a后觀測孔地下水位降速開始變緩,開采20 a后變的更加平緩,但開采30 a也未達(dá)到穩(wěn)定狀態(tài),說明水源地屬開采“消耗”型水源地;隨著四柏樹水源地開采量的增加,開采20 a后煉油廠水源地水位(頭)下降速度由9 mm/a至33 mm/a逐漸增加。

圖 4 地下水位(頭)歷時(shí)曲線

四柏樹水源地增采后并未影響地下水徑流的總方向,地下水仍由南西向北東方向徑流。由圖 5可知：水源地開采后的水位降已擴(kuò)展到四周邊界,隨著開采量的增加水位降在增大,區(qū)域上地下水位(頭)的下降幅度差異不是很大,近于等幅下降,僅在水源地開采區(qū)內(nèi)降幅增大,兩水源地大致以蘆河為界分別形成兩個(gè)地下水降落漏斗,且均以水源地開采區(qū)為中心呈近似圓形的寬淺的封閉降落漏斗,四柏樹水源地漏斗較為明顯。

圖5 30年末I-I'剖面地下水位

(2)兩水源地相互影響分析。根據(jù)圖6可得：隨著開采量的增加,兩水源地水位(頭)降均呈逐漸增大的趨勢;四柏樹水源地的水位(頭)降最多至8.09 m,煉油廠水源地的水位(頭)降最多至3.58 m;四柏樹水源地增采使得煉油廠水源地30 a末水位(頭)降均以1.07的比例同趨勢逐漸增大,但兩水源地水位(頭)降均在可控范圍以內(nèi)。

3.3 煉油廠水源地單獨(dú)增采

(1)地下水位時(shí)空變化趨勢。對比圖7與圖4、圖8與圖5可得：地下水位下降規(guī)律、地下水流場與降深場變化規(guī)律與四柏樹水源地單獨(dú)開采情況基本一致;區(qū)別僅在于煉油廠水源地的漏斗隨著開采量的增加越來越明顯。隨著四柏樹水源地開采量的增加,開采20 a后煉油廠水源地水位(頭)下降速度由6 mm/a至28 mm/a逐漸增加。

(2)兩水源地相互影響分析。由圖9可知：隨著開采量的增加,兩水源地水位(頭)降均呈逐漸增大的趨勢;四柏樹水源地的水位(頭)降最多至2.60 m,煉油廠水源地的水位(頭)降最多至10.88 m;煉油廠水源地增采使得四柏樹水源地30 a末水位(頭)降均以0.97的比例同趨勢逐漸增大,但兩水源地水位(頭)降均在可控范圍。

圖6 水源地井點(diǎn)單元最大水位降深曲線

3.4 區(qū)域水資源影響對比分析

水源地開采對區(qū)域水資源的影響主要體現(xiàn)在地下水位(頭)。無論四柏樹水源地還是煉油廠增采都會引起水源地水位(頭)下降速度的成倍增加,而同樣兩水源地的增采也會引起水源地30 a末地下水位(頭)降深以某一比例增加,這說明兩水源地間的相互影響較明顯,在水源地實(shí)際增采時(shí)需統(tǒng)籌。

圖 7 地下水位(頭)歷時(shí)曲線

圖8 30年末I-I′剖面地下水位(頭)

圖9 水源地井點(diǎn)單元最大水位降深曲線

4 結(jié)論

通過數(shù)值評價(jià),綜合分析了兩水源地的運(yùn)行對區(qū)域水環(huán)境的影響。研究表明：

(1)兩水源地均屬開采“消耗”型水源地。隨著開采量的增加,水位(頭)降在增大,區(qū)域上地下水位(頭)的下降幅度差異不是很大,近于等幅下降,僅在水源地開采區(qū)內(nèi)降幅增大,兩水源地大致以蘆河為界分別形成兩個(gè)地下水降落漏斗,且均以水源地開采區(qū)為中心呈近似圓形寬淺的封閉降落漏斗。

(2)兩水源地在現(xiàn)狀條件下單獨(dú)增采時(shí),隨開采量的增加,兩水源地開采后相對穩(wěn)定的水位(頭)下降速度及30 a末水位(頭)降深均呈逐漸增大的趨勢,相互的影響也越來越顯著,在水源地實(shí)際增采時(shí)需統(tǒng)籌安排。

(3)水資源數(shù)值評價(jià)是預(yù)測開采條件下遠(yuǎn)期年區(qū)域水資源狀況的有效手段,對區(qū)域生態(tài)環(huán)境的保護(hù)與水資源合理開發(fā)利用有著積極的作用。

[1] 王瑋,馬思錦.靖邊平原地下水開發(fā)利用模式初探[C]//中國西部環(huán)境問題與可持續(xù)發(fā)展國際學(xué)術(shù)研討會論文集.2004：588-592.

[2] 王瑋,馬思錦,郭洪鈞.四柏樹水源地三維水文地質(zhì)數(shù)值模擬中的參數(shù)最優(yōu)估計(jì)[J].地下水,2003,25(3)：141-146.

[3] 孫訥正.地下水流的數(shù)學(xué)模型和數(shù)值方法[M].北京：地質(zhì)出版社,1981：25-51.

[4] 李俊亭.地下水流數(shù)值模擬[M].北京：地質(zhì)出版社,1992：44-52.

[5] Chiang W H,Kinzelbach W.3D-groundwater modeling with PMWIN a simulation system for modeling groundwater flow and pollution[M].New York：Springer-Verlag Berlin Heidelberg,2005.

[6] 李佩成.地下水動力學(xué)[M].北京：農(nóng)業(yè)出版社,1993：183-201.

[7] 陳雨孫,顏明志.抽水試驗(yàn)原理與參數(shù)測定[M].北京：水利電力出版社,1985.

[8] 鄒正盛,鄭清潔.求水文地質(zhì)參數(shù)的計(jì)算機(jī)配線法[J].工程勘察,2001(6)：30-33.