基于改進積分型Richards方程的地下水開采方案環(huán)境影響評估研究

孔令堯

(遼寧省阜新水文局，遼寧 阜新 123000)

對于干旱半干旱區(qū)域而言，地下水占據(jù)水資源的比重均在50%以上，是區(qū)域水資源總量的重要組成。隨著當前水資源矛盾的日益緊張，對于區(qū)域地下水的開采量也逐步較大，而由于不當?shù)牡叵滤_采，對區(qū)域地下水環(huán)境產(chǎn)生的影響也越來越顯著。為此對不同開采方案下的地下水環(huán)境進行有效評估，確定適宜的地下水開采方案對于區(qū)域地下水環(huán)境保護十分重要。近些年來，對于地下水開采下環(huán)境影響評估的研究也逐步成為國內(nèi)學者的研究熱點[1- 5]。采用的方法主要為兩種，一種是觀測試驗方式，這種方式較為直接，但這種試驗方式需要大量的人力和物力；另一種方式為采用數(shù)值模型，對區(qū)域地下水進行數(shù)值模擬，分析不同地下水開采方案下的地下水模擬結(jié)果，來對環(huán)境進行評估[6- 10]。這種方式較為快捷，這種方式在地下水模擬試驗中應(yīng)用較多。遼西為遼寧省地下水開采量最大區(qū)域，年地下水開采量在900萬m3左右，而地下水開采量的加劇，勢必造成區(qū)域地下水環(huán)境的影響，為對區(qū)域地下水環(huán)境進行保護，本文結(jié)合改進的積分型Richards方程對遼寧西部某區(qū)域不同地下水開采方案下的環(huán)境進行評估，提出適宜的地下水開采方案。

1 改進的積分型Richards方程計算原理

改進的積分型Richards方程計算每個計算節(jié)點的控制水量，各節(jié)點控制水量的計算方程為：

(1)

式中，ΔWi—第i個節(jié)點控制水量，m3；Bβ—β單元的控制面積，m；Aβ—β單元的控制深度，m；ΔHi—第i個節(jié)點水頭，m；Ss—節(jié)點給水率。在進行數(shù)值求解時，需要對能量方程的源匯項進行確定，源匯項方程為：

(2)

式中，Qsi—源匯等效流量，m3/s；Si—控制單元i的源匯項。在進行源匯項計算后，對控制單元的垂向流量進行計算，方程為：

(3)

(4)

(5)

(6)

式中，Δβ—控制單元β的平均水力傳導(dǎo)度；Hs—節(jié)點s的水頭，m；bi、bs、ci、cs—水力傳導(dǎo)系數(shù)。此外，在控制單元不同方向流量計算的基礎(chǔ)上，對各節(jié)點進行均衡分析，分析方程為：

(7)

式中，各變量同上述方程中的變量含義，Qsi和Qvi采用隱式方程進行求解。

2 實例應(yīng)用

2.1 區(qū)域概況

本文以遼寧西部某區(qū)域為研究區(qū)域，區(qū)域地下基巖含水層的富集程度分布十分不均勻，地下水主要分布在風化巖裂隙以及線性充水帶中，上、中、下更新層的砂礫石分布不均，水文地質(zhì)較為復(fù)雜，區(qū)域單井涌量變化較為穩(wěn)定，穩(wěn)定在50～150m3/d，整個地下系統(tǒng)主要為沖擊巖組成。受地下水開采量逐年增加的影響，區(qū)域地下水環(huán)境造成一定程度的影響。為對區(qū)域不同地下水開采方案進行環(huán)境影響評估，對區(qū)域內(nèi)布設(shè)6組地下水觀測井。各觀測點的布設(shè)位置如圖1所示。

圖1 觀測井分布位置

2.2 模型參數(shù)動態(tài)識別及驗證

為對模型進行參數(shù)動態(tài)識別并對模型進行驗證，結(jié)合6個區(qū)域觀測點的組水文地質(zhì)參數(shù)進行設(shè)定，并對比分析改進前后模型的計算質(zhì)量及時間，結(jié)果見表1、2及圖1。

表1 模型參數(shù)動態(tài)識別結(jié)果

表2 模型計算質(zhì)量及計算時間對比

*Em表示計算質(zhì)量；tc表示網(wǎng)格控制單元計算時間。

表1中給出了改進的積分型Richards方程下區(qū)域不同水文地質(zhì)下的主要敏感參數(shù)。結(jié)合參數(shù)設(shè)置結(jié)果對比分析改進前后模型的計算質(zhì)量和計算時間，從6組樣本序列計算質(zhì)量和計算時間的對比結(jié)果可以看出，相比于傳統(tǒng)模型，改進模型下各組樣本序列計算時間縮短約11.5倍，格網(wǎng)點計算質(zhì)量整體提高22倍。

2.3 不同開采方案下的潛水位動態(tài)變幅分析

結(jié)合不同地下水開采方案(見表3)，應(yīng)用改進的積分型Richards方程模擬分析了不同開采方案下的地下水潛水位的變幅，如圖2、3所示。

表3 不同地下水開采方案

圖2中方案下淺層地下水前期變幅波動性較大，后期區(qū)域穩(wěn)定，各觀測點整體下降率平均為1.2～1.5m，后期區(qū)域地下水變化較為平穩(wěn)，在在深層區(qū)域，由于深層地下水壓采量的減少，深層地

圖2 方案1下地下水變幅結(jié)果

圖3 方案2下地下水變幅結(jié)果

圖4 不同開采方案下的區(qū)域地下水鹽堿化風險評估結(jié)果

下水水位整體抬升。而在方案2下，由于同時開出深層和淺層地下水，淺層地下水水位變幅明顯高于方案1下的水位變幅，各觀測點整體下降水位為2.5～4.8m，而受到深層地下水開采影響，方案2下深層地下水變幅要低于方案1，不同方案下的地下水潛水位降低和抬升幅度主要受到各方案的壓采率密切相關(guān)。

2.4 不同開采方案下的地下水鹽堿化風險評估

結(jié)合改進積分型Richards方程實現(xiàn)了區(qū)域地下水的三維數(shù)值模擬，并結(jié)合地下水鹽堿化風險評估指標對區(qū)域不同開采方案下的地下水鹽堿化風險進行評估，各開采方案下的評估結(jié)果如圖4所示。

從區(qū)域地下水鹽堿化風險評估結(jié)果可看出，發(fā)生地下水鹽堿化風險值較高區(qū)域地下水水位埋深為3.5m，而在潛水埋深2.5～3.5m之間的區(qū)域較易發(fā)生次生的地下水鹽堿化風險。從方案對比結(jié)果可以看出，方案1下區(qū)域發(fā)生地下水鹽堿化風險程度高于方案2，建議區(qū)域主要采用方案2進行地下水的開采。從區(qū)域地下水鹽堿化風險分布可看出，區(qū)域的中南部較易發(fā)生地下水鹽堿化風險，應(yīng)盡量減少該區(qū)域的地下水開采，維持區(qū)域地下水環(huán)境健康。

3 結(jié)論

(1)改進的積分型Richards方程較傳統(tǒng)方法的計算時效性改善程度較高，可用于非飽和水流的數(shù)值模擬。

(2)只開采深層地下水，其水位較同時開采淺層和深層地下水的方案下的水位變幅波動性更大，建議應(yīng)采取淺層和深層地下水同時開采的方案，適當調(diào)整方案的壓采率，保持區(qū)域地下水的平穩(wěn)變化。

(3)方案1下區(qū)域發(fā)生地下水鹽堿化風險程度高于方案2，區(qū)域中南部較易發(fā)生地下水鹽堿化風險，應(yīng)盡量減少該區(qū)域的地下水開采。

[1] 段艷華. 淺層地下水系統(tǒng)中砷富集的季節(jié)性變化與機理研究[D]. 中國地質(zhì)大學, 2016．

[2] A. H. 阿辛頓, 錢卓洲, 付湘寧. 地下水生態(tài)系統(tǒng)的分析與評估[J]. 水利水電快報, 2015, 36(12): 13- 16+26．

[3] Hussein Ibrahim Ahmed. 中國北方大型灌區(qū)地下水數(shù)值模擬和預(yù)測[D]. 西北農(nóng)林科技大學, 2014．

[4] 魏鈺潔. 海水入侵區(qū)的地下水開采控制方法及應(yīng)用研究[D]. 鄭州大學, 2013．

[5] 項國圣. 黑河中游張掖盆地地下水開發(fā)風險評價及調(diào)控[D]. 蘭州大學, 2011．

[6] 苗添升, 盧文喜, 歐陽琦, 等. 地下水數(shù)值模擬的不確定性分析在水質(zhì)預(yù)測中的應(yīng)用[J]. 水電能源科學, 2016, 34(08): 20- 23+44．

[7] 陳南祥, 楊杰, 屈吉鴻. 中牟縣地下水生態(tài)水位研究[J]. 華北水利水電大學學報(自然科學版), 2016, 37(01): 84- 88.

[8] 張麗, 溫斯鈞, 田曉龍. 鞏義市城區(qū)地下水水質(zhì)評價分析[J]. 華北水利水電大學學報(自然科學版), 2016, 37(05): 68- 71.

[9] 張麗, 田曉龍. 鞏義市城區(qū)地下水動態(tài)特征分析[J]. 華北水利水電大學學報(自然科學版), 2015, 36(03): 24- 27.

[10] 嚴鋒, 郭玉法, 劉波, 等. GIS技術(shù)在地下水資源研究中的應(yīng)用現(xiàn)狀及前景[J]. 水利技術(shù)監(jiān)督, 2005(05): 63- 65．