徐鐵兵, 夏凡, 閆佰忠, 王宏亮, 張艷

(1.河北省生態(tài)環(huán)境科學研究院, 石家莊 050031; 2.河北省高校生態(tài)環(huán)境地質應用技術研發(fā)中心, 石家莊 050031;3.河北省水資源可持續(xù)利用與產業(yè)結構優(yōu)化協(xié)同創(chuàng)新中心, 石家莊 050031; 4.河北省水資源可持續(xù)利用與開發(fā)重點實驗室, 石家莊 050031)

近年來,隨著中國經濟的高速發(fā)展,對水資源的需求不斷增加,開采量逐漸增大,地下水危機日漸凸顯。河北省屬于中國北方典型資源性缺水地區(qū),受地下水持續(xù)超采影響,地下水位呈不斷下降趨勢,并因此誘發(fā)了含水層疏干、水位降落漏斗、水質惡化、地面沉降等一系列地質環(huán)境問題[1-2],已成為影響和制約當?shù)亟洕鐣沙掷m(xù)發(fā)展的重要瓶頸。政策性人工干預作為地下水生態(tài)恢復一項重要途徑,成為中國華北地區(qū)地下水超采治理的重要手段。

目前,已有諸多學者對地下水動態(tài)及人工干預治理理論進行了研究。研究方法有神經網(wǎng)絡法[3-5]、多元回歸分析法[6-7]、數(shù)值模擬法[8],其中數(shù)值模擬法運用最為廣泛。賈金生等[9]利用數(shù)值模擬軟件構建了河北省欒城縣三維地下水數(shù)值模型,模擬了不同開采強度情境下的地下水動態(tài)響應。韓宇平等[10]在建立華北平原人民勝利渠灌區(qū)地下水運動數(shù)值模擬模型的基礎上,進一步預測了不同氣候情境下地下水位變化。汪麗芳等[11]采用Visual MODFLOW數(shù)值模擬模型預測了邯鄲東部平原多個開采情境下的地下水位動態(tài),分析了開采量與地下水位的關系,并建立了二者之間的多元回歸方程。這些研究所選典型區(qū)地下水類型多以孔隙水、裂隙水為主,對于巖溶區(qū)地下水人工調控研究較少。

涉縣盆地作為河北省典型巖溶區(qū),巖溶水發(fā)育較好,但由于地下水的開采,導致泉流量減少,甚至產生斷流[12]。現(xiàn)利用涉縣2011—2021年地下水位動態(tài)監(jiān)測資料,系統(tǒng)分析由于人類活動影響產生的地下水環(huán)境變化特征,并在此基礎上以GMS(groundwater modeling system)數(shù)值模擬為平臺,考慮地下水壓采的實際情況,對地下水壓采人工干預調控下地下水環(huán)境的動態(tài)響應進行研究,分析未來地下水環(huán)境的動態(tài)變化,以期為河北省典型巖溶區(qū)規(guī)劃和使用地下水資源提供科學合理的依據(jù)。

1 研究區(qū)概況

東風湖泉域位于河北省涉縣西南部,西部與陜西省省黎城縣交界,泉域面積924 km2。所屬涉縣地處113°26′E～114°00′E,36°17′N～36°55′N。年平均氣溫10.7～14.2 ℃,年平均降水量為571.7 mm,年內降水分配不均,多集中在7—9月,約占全年70%,年平均蒸發(fā)量為1 732.1 mm,屬暖溫帶半濕潤大陸性季風氣候。境內水系主要有漳衛(wèi)南運河水系的清漳河、濁漳河和漳河,屬山溪河流。地勢東北低西南高,差距大,海拔相差1 300 m,地貌復雜。

東風湖泉域主要發(fā)育地層包括中元古界長城系的大紅峪組的石英砂巖,分布于清漳河兩岸河谷地帶;古生界寒武系和奧陶系白云巖、灰?guī)r,在全區(qū)廣泛分布;以及新生界第四系松散堆積物,主要分布于東南部、河谷兩岸。含水層類型主要包括松散巖類孔隙含水巖組、碳酸巖類巖溶裂隙含水巖組和巖漿巖、砂巖風化裂隙含水巖組。碳酸巖類巖溶裂隙含水巖組富水性較好,是區(qū)內主要取水層。

涉縣東風湖泉域地下水監(jiān)測工作起步較早,區(qū)內監(jiān)測點位較多,本次工作選擇區(qū)內32眼地下水監(jiān)測井2011—2021年的動態(tài)資料(圖1),并結合對應年份的氣象和開采量等資料進行綜合分析,地下水位數(shù)據(jù)等資料主要通過實際調查測量,此外部分收集自邯鄲水文水資源勘測局。

圖1 研究區(qū)位置及地下水監(jiān)測點分布圖Fig.1 Location of the study area and distribution map of groundwater monitoring points

2 地下水環(huán)境動態(tài)特征

2.1 年際動態(tài)變化特征

研究區(qū)巖溶地下水年際動態(tài)呈現(xiàn)一定的周期性與滯后性,從時間上可分為:2011—2013年、2014—2017年、2018—2020和2020—2021年4個時段。其中2011—2021年水位整體呈上升趨勢,地下水位平均變幅15 m左右;2011—2013年地下水位由于開采量增加整體呈下降趨勢,包氣帶厚度變厚,水位變幅10 m左右;2014—2018年水位呈波動回升狀態(tài),特別在2013年后由于降雨量增加,開采量減小,地下水位抬升明顯;2018—2020年水位呈下降狀態(tài);2020年由于降雨量較大,水位上升幅度較大,較前一年同期地下水位上漲27 m;可見區(qū)內巖溶地下水動態(tài)受降雨和開采的雙重影響(圖2和圖3)。

圖2 巖溶地下水位與降雨量關系圖Fig.2 Relationship between Karst groundwater level and rainfall

圖3 巖溶地下水位與開采量關系圖Fig.3 Fig.2 Relationship between Karst groundwater level and rainfall

東風湖泉域地下水位空間變幅不同,其中G701監(jiān)測井位于上游偏城村,2011—2021年水位整體呈上升趨勢,地下水位平均變幅15 m左右,最高水位出現(xiàn)在2020年10月為705.18 m,最低水位出現(xiàn)在2014年7月為646.76 m。受開采影響,2011—2014年水位整體呈下降趨勢,水位變幅8 m左右;隨著地下水壓采2014—2020年水位呈波動狀態(tài);2020年由于年降水量多達803.7 mm,地下水水位上升幅度較大,上漲45 m,年末有下降趨勢(圖4)。

圖4 上游偏城村G701監(jiān)測孔地下水位變化曲線Fig.4 Groundwater level variation curve of G701 in Piancheng

SD01監(jiān)測井位于下游河南店鎮(zhèn),2011—2021年水位也整體呈上升趨勢,地下水位平均變幅20 m左右,最高水位出現(xiàn)在2017年12月為432.23 m,最低水位出現(xiàn)在2013年12月為406.89 m。受開采影響,2011—2013年水位整體呈下降趨勢,水位變幅11 m左右;隨著地下水壓采2014—2020年水位呈波動上升狀態(tài),且水位基本可以維持高位,如圖5所示。可見,相對于上游地區(qū),下游區(qū)地下水位能夠得到上游補給,地下水位回升后基本能夠維持。

圖5 下游河南店鎮(zhèn)SD01監(jiān)測孔地下水位變化曲線Fig.5 Groundwater level variation curve of SD01 in Henandian

由圖6可知,研究區(qū)地下水位空間分布特征不同,其中大部分地區(qū)水位為上升區(qū)或者基本穩(wěn)定區(qū)。上升區(qū)主要分布在研究區(qū)的下游地區(qū),包括井店、神頭、涉城、索堡、河南店等一帶,其中井店鎮(zhèn)和南溝一帶地下水水位上升幅度較大,上升幅度基本在10 m以上;穩(wěn)定區(qū)基本分布于研究區(qū)的中部條帶區(qū),地下水位變幅在-2～2 m;在研究區(qū)的東北部和西北區(qū)的上游區(qū),由于地下水得到補給來源較為單一,此外地下水還需補給下游地區(qū),地下水位基本呈現(xiàn)下降的趨勢,特別是在東北部的偏城鎮(zhèn)-王家莊-雞鳴鋪一帶地下水位下降幅度在15 m左右。由于上游處于補給區(qū),地下水位較高,包氣帶厚度較薄,得到降雨入滲的能力較強,在上游補給區(qū)地下水水質狀況較好,水化學類型主要為HCO3-Ca型。水中TDS整體上表現(xiàn)出西北高、東南低的特點。

2.2 年內動態(tài)變化特征

由圖7和圖8可知,地下水動態(tài)類型為降水入滲-開采型,主要影響因素為降水補給和人工開采,水位升降變幅明顯。地下水位年內動態(tài)變化規(guī)律為:年初的1月份春灌開始后,地下水水位由上升變?yōu)橄陆?3—6月隨著對地下水開采量的增加,地下水水位下降速率變快,在汛期到來之前的6月末降至年內最低水位(開采量大者會延遲到9月)。進入雨季后,降雨補給增加、農業(yè)開采相應減小,地下水水位開始逐漸回升,直至年末或來年1月份達到最高水位。

圖7 上游偏城村G701監(jiān)測孔2021年年內地下水位變化曲線Fig.7 Groundwater level variation curve of G701 in Piancheng during 2021

圖8 下游河南店鎮(zhèn)SD01監(jiān)測孔2021年年內地下水位變化曲線Fig.8 Groundwater level variation curve of SD01 in Henandian during 2021

3 地下水數(shù)值模型構建

3.1 水文地質概念模型

研究模型范圍為:北部東段以大巖村—土林背村—王家莊一帶存在的地下分水嶺為界;西段以震旦系、寒武系下統(tǒng)為界,概化為二類流量邊界;西部為震旦系區(qū)域隔水底板構成的確定性隔水邊界,概化為零流量邊界;東部北段以西戌鎮(zhèn)以北的長亭阻水斷層為界,南段以西戌鎮(zhèn)以南的涉縣阻水斷裂為界,概化為零流量邊界;南部東段以涉縣阻水斷裂為界,西段以黎城—源泉地下分水嶺為界,概化為隔水邊界;總面積為92 400 hm2,寒武系、奧陶系灰?guī)r裸露面積65 900 hm2。

對于垂向含水層結構的概化,中上寒武統(tǒng)、奧陶統(tǒng)作為研究的目的層。目的層上覆地下水類型為第四系孔隙水,該含水層薄且不穩(wěn)定,富水性較差,其下部為連續(xù)巨厚的黏土層,隔水性較好,故可概化為目的層的隔水頂板,其下伏的寒武系下統(tǒng)與其下的震旦系,因富水性較差構成模型區(qū)域性隔水底板。

平面上根據(jù)地層巖性結構及水動力條件,劃分為補給區(qū)、徑流區(qū)、排泄區(qū)。其中,研究區(qū)西、北部存在較大面積灰?guī)r,透水性好,故該區(qū)域為模型目的層的主要補給區(qū)、涉縣向斜盆地覆蓋區(qū)為模型目的層的徑流區(qū),東風湖泉群為排泄區(qū)。區(qū)內含水介質為巖溶裂隙發(fā)育,具有較好的連通性,同時因其具有非均質各向異性特征,故模型可概化為各向異性非均質二維地下水系統(tǒng)模擬。

3.2 地下水流數(shù)學模型

基于研究區(qū)水文地質概念模型,數(shù)值模型地下水數(shù)學模型如式(1)～式(4)所示。

(x,y)∈Ω,t≥0

(1)

(2)

(3)

(4)

式中:x、y為水流方向;t為時間;Ω為計算區(qū)域;H為地下水流系統(tǒng)水位標高;K為滲透系數(shù);Kn為二類邊界法線方向上n的滲透系數(shù);Ke為零通量邊界法線方向e上的滲透系數(shù);M為含水層厚度(L);ε為單元體內的源匯項;H0為初始水位分布;μ為含水層給水度;q(x,y,t)為二類邊界的單位面積流量;Γ2為滲流區(qū)域二類邊界。

3.3 模型識別與驗證

以2015年12月的地下水位作為模型初始水位,通過反復調整模型參數(shù)和源匯項變量進行地下水流場和觀測井孔的水位擬合,并利用2020年12月地下水模型流場圖與實際進行擬合對比,整體對模型進行識別驗證,使得模擬的地下水均衡變化與計算均衡量基本相符,且識別的水文地質參數(shù)要符合實際水文地質條件。經過參數(shù)調試,監(jiān)測孔模擬水位和觀測水位相差較小,如圖9和圖10所示,水位變化符合年際水位變化特點,說明數(shù)值模擬對含水層的結構、邊界條件的概化,水文地質參數(shù)的選取是合理、可靠的,可以比較真實地反映研究區(qū)地下水流系統(tǒng)特征,識別后的模型參數(shù)分區(qū)值及水均衡分析如圖11、表1和表2所示。

表1 模擬區(qū)識別驗證后水文地質參數(shù)Table 1 Hydrogeological parameters after identification and verification of simulated areas

表2 模擬期地下水量均衡表Table 2 Groundwater balance table during simulation period

圖9 監(jiān)測孔模擬地下水水位和觀測地下水水位Fig.9 Monitoring hole simulation of groundwater level and observation of groundwater level

圖10 驗證期地下水實測流場和計算流場擬合圖(2020年12月)Fig.10 Validation period groundwater measured flow field and calculated flow field fitting diagram (December 2020)

圖11 模擬區(qū)水文地質參數(shù)分區(qū)圖Fig.11 Simulation area hydrogeological parameter zoning map

4 人工干預下地下水環(huán)境響應特征

根據(jù)研究區(qū)開采現(xiàn)狀及華北平原地下水超采區(qū)綜合治理方案,利用識別后的數(shù)值模型,對未來地下水壓采條件下地下水流場變化規(guī)律進行預測。其中源匯項中降雨量按照豐、平、枯規(guī)律選擇歷史數(shù)據(jù)給入,地下水開采量在2021年基礎上壓采20%后為2 590×104m3/a,其他源匯項利用2011—2021年平均值給入模型。預測期至2029年。

將2029年地下水位預測結果(圖12)與2020年12月相比,得出變差圖(圖13)。由圖12和圖13可知,經過進一步地下水壓采,研究區(qū)地下水位都有不同程度的回升,起到了積極良好的水位恢復效果,區(qū)內地下水位平均響應特征值達到了5～10 m,其中,位于研究區(qū)南部的鐵西水源地、泉域區(qū)、下游斷裂邊界位置處淺層水位響應特征值超過15 m,壓采后水位上升幅度較大;研究區(qū)中游的天平莊、后池耳、鳳崗一帶地下水位響應特征值達到了5～10 m,水位回升明顯;研究區(qū)的上游地區(qū),即偏城鎮(zhèn)、西廟灣等,地下水位響應特征值為3～5 m,總體變化不大;整體來看,研究區(qū)大部分地區(qū)地下水位響應特征值為5～10 m,下游地下水位響應要整體大于上游,分析其原因是地下水位響應在空間分布上呈現(xiàn)不同的規(guī)律。包氣帶厚度在地下水位回升較大地區(qū)逐漸變薄,而且主要集中于研究區(qū)的中下游?？梢?地下水壓采對于中下游地下水水環(huán)境改善有較大促進作用。

圖12 人工壓采下2029年預測地下水流場圖Fig.12 Prediction of groundwater flow field under artificial extraction in 2029

圖13 2029年12月相對2020年12月地下水水位變幅Fig.13 Change in groundwater level in December 2029 compared to December 2020

5 結論

結合研究區(qū)地質、水文地質資料,在系統(tǒng)分析區(qū)內地下水動態(tài)特征的基礎上,構建了東風湖泉域地下水數(shù)值模擬模型,并結合實際壓采措施,模擬了未來在壓采情景下的地下水動態(tài)變化過程,得出以下結論。

(1)研究區(qū)地下水動態(tài)變化明顯,具有周期性、滯后性,具體表現(xiàn)為4個變化階段。地下水類型為降水入滲-開采型,低水位期基本在5—6月,高水位期在1—2月。

(2)研究區(qū)內巖溶地下水整體呈現(xiàn)回升的狀態(tài),但是空間分布特征不同,上升區(qū)主要分布于下游地區(qū)。

(3)在干預條件下,區(qū)內地下水位有著不同程度的回升,南部(下游)淺層水位響應特征值超過15 m,中游地下水位響應特征值達到了8～10 m,北部(上游)地下水位響應特征值為3～5 m;上游和中游地下水位影響整體小于下游,下游水位響應尤為突出。