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

孔令堯

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

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

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

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

(1)

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

(2)

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

(3)

(4)

(5)

(6)

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

(7)

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

2 實例應用

2.1 區域概況

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

圖1 觀測井分布位置

2.2 模型參數動態識別及驗證

為對模型進行參數動態識別并對模型進行驗證，結合6個區域觀測點的組水文地質參數進行設定，并對比分析改進前后模型的計算質量及時間，結果見表1、2及圖1。

表1 模型參數動態識別結果

表2 模型計算質量及計算時間對比

*Em表示計算質量；tc表示網格控制單元計算時間。

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

2.3 不同開采方案下的潛水位動態變幅分析

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

表3 不同地下水開采方案

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

圖2 方案1下地下水變幅結果

圖3 方案2下地下水變幅結果

圖4 不同開采方案下的區域地下水鹽堿化風險評估結果

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

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

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

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

3 結論

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

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

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

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

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

[3] Hussein Ibrahim Ahmed. 中國北方大型灌區地下水數值模擬和預測[D]. 西北農林科技大學, 2014．

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

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

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

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

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

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

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