朱洪生,王繼華,豆敬峰,朱瑞霞
(1.河南省地質(zhì)環(huán)境監(jiān)測院,河南鄭州450016;2.河南省地質(zhì)環(huán)境保護重點實驗室,河南鄭州450016;3.平輿縣楊埠鎮(zhèn)農(nóng)業(yè)服務(wù)中心,河南平輿463400)
近年來,隨著全球氣候變暖,水循環(huán)速度加快,汛期極端強降水事件增多,引起洪澇災(zāi)害頻發(fā)。我國北方地區(qū)受地下水資源開采影響,地下水位多呈下降趨勢。降水作為地下水的主要補給來源,影響著地下水資源的動態(tài)變化特征,極端強降水的輸入導(dǎo)致地下水產(chǎn)生變化。目前,研究地下水動態(tài)與其影響因素的理論方法包括多元回歸分析法[1-2]、時間序列分析預(yù)測模型[3-4]、灰色系統(tǒng)模型[5-6]和地質(zhì)統(tǒng)計學(xué)方法[7-8]等。平建華等[9]總結(jié)回顧了地下水動態(tài)預(yù)測模型,將各種模型分為確定性模型和隨機模型,并指出了各種模型的適用條件;陶虹等[10]研究了關(guān)中城市群50 a地下水動態(tài)變化及其影響因素,指出該地區(qū)地下水位整體呈下降趨勢,地下水開采是其主要影響因素;危潤初等[11]采用標(biāo)準(zhǔn)化降水指數(shù)研究了黑龍江建三江地區(qū)地下水位變化的突變點,指出地下水位下降受水稻種植面積以及降水影響;遲寶明等[12]利用基于遺傳算法的BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測了元寶山露天礦區(qū)地下水動態(tài)特征;張文化等[13]研究了氣候變化與人類活動對石羊河流域地下水動態(tài)變化的影響,指出人類活動在地下水動態(tài)變化中占主導(dǎo)因素。上述主要是對降水的周期性、頻率,以及降水的變化趨勢、空間分配等的研究,而對于典型區(qū)域強降水條件下地下水動態(tài)響應(yīng)的研究較少,單次強降水引起的地下水動態(tài)響應(yīng)和預(yù)測研究更少。
豫北平原位于河南省北部,地下水過量開采嚴(yán)重,造成地下水埋深很深,包氣帶厚度增大。近些年,該地區(qū)汛期強降水頻發(fā),特別是2016年7月,華北地區(qū)出現(xiàn)5次明顯的降水過程,河北北部、北京、河南等地出現(xiàn)了50 mm以上的暴雨天氣。筆者利用豫北地區(qū)水文地質(zhì)、地下水位統(tǒng)測資料和水質(zhì)監(jiān)測資料等,建立研究區(qū)地下水?dāng)?shù)值模型,對強降水條件下地下水動態(tài)響應(yīng)進行研究。
研究區(qū)位于黃河以北的河南省北部。地形西高東低,高程42~1 725 m。研究區(qū)受地質(zhì)構(gòu)造控制明顯,地貌類型分為山高谷深、山勢陡峻雄偉的斷塊山地、堆積地貌,局部發(fā)育侵蝕、剝蝕型地貌。
研究地層主要有太古界、中太古界、古生界、中生界與新生界的古近系、新近系和第四系地層。研究區(qū)地下水類型可分為松散巖類孔隙水、碎屑巖類裂隙孔隙水和碳酸鹽巖類裂隙巖溶水,其中:松散巖類孔隙水在垂向上劃分為淺層地下水系統(tǒng)和深層地下水系統(tǒng),在平面上分為漳衛(wèi)河地下水系統(tǒng)(Ⅰ)、黃河地下水系統(tǒng)(Ⅱ);碎屑巖類裂隙孔隙水可劃分為安陽丘陵區(qū)裂隙孔隙水系統(tǒng)(Ⅰ)和湯陰火龍崗裂隙孔隙水系統(tǒng)(Ⅱ);碳酸鹽巖類裂隙巖溶水系統(tǒng)可分為焦作九里山泉巖溶水(Ⅰ)、輝縣百泉巖溶水系統(tǒng)(Ⅱ)、鶴壁許家溝泉巖溶洞水系統(tǒng)(Ⅲ)和安陽小南?!渲槿獛r溶水系統(tǒng)(Ⅳ)。
地下水位監(jiān)測包括2016年1月—2017年12月豫北地區(qū)淺層地下水動態(tài)監(jiān)測點34個,深層地下水監(jiān)測點4個;水質(zhì)監(jiān)測點41個,2016年4月共采集41組水樣,2016年7月、9月10個加密點各采集1次共20組水樣。
為了研究2016年7—9月強降水對地下水位的影響,建立地下水?dāng)?shù)值模擬模型,對地下水位變化趨勢進行研究。
2.2.1 水文地質(zhì)概念模型
模擬區(qū)為整個豫北平原,總面積約21 460 km2。模擬計算的目的層為第四系全新統(tǒng)-下更新統(tǒng)孔隙含水層,含水層概化為2層,包括潛水含水層和承壓含水層。模型上部邊界為潛水面,為水量交換邊界,包括降水入滲、潛水蒸發(fā)、河渠入滲和灌溉水回滲等,下部多為第三系碎屑巖類,處理為隔水邊界。模型西部以山區(qū)與平原自然分界線為邊界,為地下水流入邊界。模型北側(cè)以行政區(qū)邊界為流量邊界,側(cè)向流入流出量根據(jù)達西定律計算得出。模型其他邊界為河流邊界,主要由黃河、金堤河和徒駭河組成。地下水接受大氣降水和地下水側(cè)向補給,整體流向為由西南流向東北,以蒸發(fā)、人工開采和地下水側(cè)向流出為主要排泄方式。
2.2.2 地下水流數(shù)學(xué)模型
模擬區(qū)第四系松散沉積物巖性在水平、垂直方向上都有較大變化,在同一點,水平方向滲透系數(shù)在各個方向大小一致,而垂直方向滲透系數(shù)小于水平方向滲透系數(shù),因此將計算區(qū)域含水層概化為非均質(zhì)各向異性介質(zhì)。模擬區(qū)地下水分為潛水和深層承壓水,因此將研究區(qū)地下水流系統(tǒng)概化為非均質(zhì)各向異性準(zhǔn)三維非穩(wěn)定地下水流系統(tǒng),可用下面偏微分方程的定解問題來描述:
式中:K為滲透系數(shù);S為自由面以下含水層的貯水率;W為源匯項;μ為潛水含水層在潛水面上的重力給水度;P 為潛水面上的降水入滲和蒸發(fā)量等;h(x,y,z,t)為水頭;h0(x,y,z)為初始壓力水頭;h1(x,y,z,t)為第一類邊界上的壓力水頭;f1(x,y,z,t)為第二類邊界上的水分通量;Kn為邊界法線方向的滲透系數(shù);t為時間;Ω 為研究范圍;→n 為有效孔隙度;Γ0為潛水面;Γ1為第一類邊界;Γ2為第二類邊界。
2.2.3 模型的識別和驗證
模擬軟件采用Modflow軟件,初始(2016年6月)流場見圖1,第一含水層滲透系數(shù)、給水度分區(qū)和第二含水層滲透系數(shù)、貯水率分區(qū)見圖2。根據(jù)模擬區(qū)地下水位資料,將2016年6—12月作為模型識別驗證期。驗證期末刻(2016年12月)地下水流場擬合見圖3。驗證期模擬水位與觀測水位擬合較好,說明含水層結(jié)構(gòu)、邊界條件的概化、水文地質(zhì)參數(shù)的選取是合理的,所建立的數(shù)學(xué)模型能較真實地刻畫研究區(qū)地下水系統(tǒng)特征,識別后含水層參數(shù)見表1。
圖1 地下水初始流場(單位:m)
圖2 模擬區(qū)分區(qū)
圖3 2016年12月地下水流場擬合(單位:m)
表1 含水層參數(shù)
2.2.4 預(yù)測情景
2016年保持降水量、蒸發(fā)量、開采量等源匯項不變;其他年份保持降水量總量不變,7—9月增加降水量30%,其他月份降水量減少30%,其他源匯項保持不變。預(yù)測未來10、20 a地下水位變化情況,以及典型漏斗區(qū)地下水位變化情況。
利用6—8月、6—12月地下水位等值線得到地下水位變幅,將研究區(qū)淺層地下水劃分為基本平衡區(qū)(-0.5~0.5 m)、緩慢下降區(qū)(-0.5~-2.0 m)、急劇下降區(qū)(<-2.0 m)、緩慢上升區(qū)(+0.5~2.0 m)和急劇上升區(qū)(>2.0 m)5個水位變幅分區(qū)。
2016年6月淺層地下水平均水位為61.11 m,平均變幅0.75 m,最大降幅(-12.50 m)位于濮陽市胡村鄉(xiāng)北豆門村,最大升幅(13.81 m)位于鶴壁市??h淇門水文站。2016年6—8月豫北地區(qū)淺層地下水位變幅分區(qū)分布見圖4,其中:急劇上升區(qū)面積2 121.89 km2,占10.65%;緩慢上升區(qū)面積6 891.28 km2,占34.58%;基本平衡區(qū)面積9 505.54 km2,占47.69%,緩慢下降區(qū)面積682.53 km2,占2.90%;急劇下降區(qū)面積577.63 km2,占4.18%。與2016年6—12月豫北地區(qū)地下水位變幅分區(qū)(見圖5)相比較,急劇上升區(qū)面積減少200.6 km2,緩慢上升區(qū)面積增大1 373.97 km2,基本平衡區(qū)面積減少905.65 km2。2016年6—8月淺層地下水位急劇上升區(qū)面積與緩慢上升區(qū)面積之和比2016年6—12月淺層地下水位急劇上升區(qū)面積與緩慢上升區(qū)面積之和多1 173.37 km2,說明7—8月強降水過后水位埋深明顯減小。
圖4 豫北地區(qū)淺層地下水位變幅分區(qū)(2016年6—8月)
圖5 豫北地區(qū)淺層地下水位變幅分區(qū)(2016年6—12月)
淺層地下水各單項組分變化情況見表2。根據(jù)2016年4—7月淺層地下水水質(zhì)數(shù)據(jù)分析水化學(xué)類型和單個化學(xué)組分的變化情況。2016年4月HCO3-Ca·Mg·Na型水占80%;2016年7月HCO3-Ca·Mg·Na型水占70%,且與2016年4月相比井位有所改變。2016年4月、7月HCO3·SO4-Ca·Mg·Na型水都占20%,但2016年7月井位有所改變。SO4·Cl-Ca·Mg·Na型水僅在2016年7月測試中有1眼井為此類型。
單個化學(xué)組分變化情況:pH值為 7.3~8.3,與2016年4月相比,7月pH值增高井?dāng)?shù)為5眼,減少井?dāng)?shù)為2眼,不變井?dāng)?shù)為3眼。溶解性總固體含量7月低于4月的有9眼。與地下水質(zhì)量標(biāo)準(zhǔn)Ⅲ類水相比,2016年4月和7月,氯化物含量均不超標(biāo),最高含量為231.49 mg/L;硫酸鹽超標(biāo)率均為10%;硝酸鹽超標(biāo)率均為50%,并且是相同井;氨氮均不超標(biāo);鐵含量超標(biāo)率均為60%,最高含量為5.09 mg/L;耗氧量4月超標(biāo)率為10%,7月不超標(biāo)。與2016年4月相比,2016年7月指標(biāo)含量減少的占比相對較多,其原因是,4月為相對枯水期,到7月受降水量增加影響,淺層地下水接受降水補給,地下水受到稀釋,指標(biāo)含量減小。指標(biāo)含量增加或不變的淺層地下水,可能是淺層地下水埋深大,降水入滲量小所致。
表2 淺層地下水各單項組分變化情況
通過預(yù)測得出10 a和20 a后不同含水層地下水流場,見圖6。研究區(qū)整個地下水流場發(fā)生了一定程度的變化,整體地下水位有所下降,不同地區(qū)水位下降程度不同,但整體地下水流向未發(fā)生明顯變化,尤其是深層地下水。對于淺層地下水,在安陽市西北側(cè),10 a后出現(xiàn)一定面積的疏干,但深層水未疏干;在安陽市東側(cè)和濮陽市北部,水位下降比較明顯,形成較大的地下水漏斗區(qū);在模擬區(qū)南部,水位下降不明顯。研究區(qū)整體地下水位的下降說明地下水的補采不平衡,開采量大于補給量。對于當(dāng)前存在的降落漏斗,如果按照當(dāng)前開采量繼續(xù)開采,則20 a后漏斗范圍將持續(xù)擴大,地下水埋深將越來越大。
圖6 含水層地下水水位等值線(單位:m)
為了詳細(xì)分析漏斗區(qū)地下水位變化情況,將安陽市東側(cè)和濮陽市大部分區(qū)域的漏斗區(qū)作為典型漏斗區(qū)進行分析。將典型漏斗區(qū)作為一個均衡區(qū),同時在漏斗區(qū)的上、中、下3個地方安置水位觀察點,分別位于內(nèi)黃縣、濮陽市和留固鎮(zhèn),通過分析地下水位隨時間的變化了解漏斗區(qū)水位的變化。由圖7可知,不同觀察點淺層地下水位和深層地下水位均出現(xiàn)不同程度的下降,且下降速率不同。濮陽市淺層地下水位20 a下降了1.6 m,下降速率比較平穩(wěn);深層地下水位下降了3.3 m,前2 a下降速率較快,后期下降速率比較穩(wěn)定。內(nèi)黃縣淺層地下水位20 a變化不大,總體下降幅度不超過1.5 m;深層地下水位下降3.6 m,前期水位下降迅速,后期較為穩(wěn)定。留固鎮(zhèn)淺層地下水位20 a下降了1.8 m,深層地下水位下降了3.7 m。按照當(dāng)前開采速率,典型漏斗區(qū)淺層和深層地下水位將持續(xù)下降。
圖7 不同觀察點地下水位隨時間變化曲線
受強降水影響,研究區(qū)地下水位變化可劃分為基本平衡區(qū)、緩慢下降區(qū)、急劇下降區(qū)、緩慢上升區(qū)和急劇上升區(qū)5個水位變幅分區(qū)。受強降水影響,地下水化學(xué)類型和水化學(xué)組分均發(fā)生了變化,水質(zhì)超標(biāo)率減小。保持當(dāng)前地下水開采量不變,降水量總量保持不變,7—9月增加降水量30%,其他月份減少降水量30%情景下,研究區(qū)地下水位將下降,不同地區(qū)水位下降程度不同,但地下水流向未發(fā)生明顯變化,尤其是深層地下水;在滲透性較好的地區(qū),10 a后將出現(xiàn)一定面積的疏干,但深層地下水未疏干,在安陽市東側(cè)和濮陽市北部形成地下水漏斗區(qū)。典型漏斗區(qū)保持現(xiàn)狀開采,20 a后淺層地下水位和深層地下水位均將出現(xiàn)不同程度的下降,濮陽市淺層地下水位下降1.6 m,深層地下水位下降3.3 m;內(nèi)黃縣淺層地下水位下降幅度不超過1.5 m,深層地下水位下降3.6 m;留固鎮(zhèn)淺層地下水位下降1.8 m,深層地下水位下降3.7 m。