梅 凱
(河北省煜環(huán)環(huán)保科技有限公司,河北石家莊050000)
根據水文地質調查,按照地下水的賦存條件、水理性質及水力特征,將研究區(qū)地下水劃分為第四系松散巖類孔隙水和奧陶系石灰?guī)r溶裂隙水。兩含水巖組無穩(wěn)定的隔水層,水力聯系密切,具有統(tǒng)一的水位,同時具有潛水特征,可視為同一地下水系統(tǒng)。
1)第四系松散巖類孔隙水。松散巖類孔隙水含水巖組,巖性主要為第四系上更新統(tǒng)的粉土、粗、中砂,底板埋深約50~100m,地下水的主要補給方式為大氣降水入滲、地表水下滲補給及側向徑流補給,主要排泄方式為向下越流、側向徑流。地下水徑流方向基本與地形坡向一致,總體徑流方向自西向東。
2)巖溶裂隙水。巖溶裂隙水含水巖組,分布于場地及其周圍的地面以下的奧陶系石灰?guī)r巖溶裂隙中,單位涌水量100~250m3/h·m,地下水的主要補給方式為上覆第四系孔隙水入滲越流補給及側向徑流,主要排泄方式為人工開采。
評價區(qū)的包氣帶巖性主要為第四系松散巖類,包氣帶厚度變化較大,松散巖類分布區(qū)包氣帶厚度般為小于40m。從地表看包氣帶巖性,第四系松散巖類分布于評價區(qū)大部分區(qū)域,主要巖性為粉土,其次為河道附近的細砂。
1)地下水的補給、徑流、排泄條件。區(qū)內地下水的補給來源主要為大氣降水入滲補給與側向徑流補給。松散巖類孔隙水:主要補給來源有兩種;其一是來自同一水文地質單元的上游水平方向地下水徑流補給及側向徑流補給;其二是地表水入滲補給。巖溶裂隙水:主要補給來源有兩種;其一是來自同一水文地質單元的上游水平方向地下水徑流補給及不同水文地質單元的側向補給;其二是松散巖類孔隙水的垂直滲入補給。兩種類型地下水在不同地段相互補給,具有統(tǒng)一的地下水位,聯系密切。
地下水徑流方向基本與地形坡向一致,評價范圍的內徑流方向是由西向東,水力坡度一般為5‰~25‰。且呈現從上游至下游水力坡度度逐漸減緩的趨勢。
排泄方式主要為工業(yè)開采利用,其次為向下越流、側向徑流及潛水蒸發(fā)。本區(qū)各類型地下水逕流距離短,具有就近補給、當地排泄的特點。
2)地下水水位動態(tài)特征。本區(qū)地下水位動態(tài)變化與補徑排條件密切相關,其年變化過程大致可分為開采下降、補給回升及相對穩(wěn)定三個階段。
水位下降期:地下水位下降主要受開采影響。歷年的4-5月份降水量比較少,又是農業(yè)春灌集中開采期,加之河道徑流量較小,造成水位緩慢下降,至5月底出現最低水位。水位下降幅度在1~3m以內。
水位回升期:6月份雨季來臨,河道徑流量增加,河道滲漏補給和大氣降水入滲補給明顯增多,加之農灌停采,區(qū)域地下水位普遍回升,回升幅度1~3m,至9月底出現最高水位。
水位穩(wěn)定期:進入10月份,受農業(yè)季節(jié)性開采影響,水位又轉為緩慢下降階段,持續(xù)到11月初,12月初至次年2月,地下水位進入相對平穩(wěn)階段,水位升降變化幅度一般小于1m。
研究區(qū)地下水動態(tài)特征為降水入滲補給——開采排泄型,本類型地下水受大氣降水垂直滲透補給和開采消耗影響明顯,水交替作用強烈,水位升降幅度較大,尤其是近年來隨著區(qū)域工業(yè)發(fā)展和新民居小區(qū)的建設,地下水開采量影響地下水動態(tài)變化特征越來越明顯。
地下水動態(tài)變化規(guī)律是:主要隨大氣降水補給量而發(fā)生變化,同時區(qū)域地下水開采量又決定地下水水位變化幅度。通過本次研究發(fā)現,研究區(qū)內除了工業(yè)用水外,大部分為鄉(xiāng)鎮(zhèn)居住小區(qū)和村莊生活用水,農業(yè)用水極少。隨著降雨量的變化,在6月份以后大氣降雨量驟然增大,但由于天氣變暖工業(yè)及生活開采量增加,地下水位并未上升速度變緩,最高水位與最低水位變幅差值一般低于3m。
由前面的論述可知,松散巖類孔隙水與巖溶裂隙水是評價區(qū)的兩個主要的地下水類型,兩類地下水之間沒有穩(wěn)定的隔水層,水力聯系密切,并且地下水位埋藏較深,決定了這兩類地下水成為本區(qū)的主要開采層。
本區(qū)的地表水體是評價區(qū)中部和南部的沙河,調查期間該河流有少量的地表徑流。從地下水位等值線圖可以看出,沙河水位高于地下水位,河流滲漏補給地下水。
采用地下水數值法進行地下水資源評價首先建立地下水系統(tǒng)的概念模型。在建立地下水系統(tǒng)概念模型的基礎上再建立地下水流動模型。
水文地質概念模型是地下水系統(tǒng)的一種近似的形象化表示,其目的是為了簡化野外實際問題,便于對該地下水系統(tǒng)進行分析,建立數學模型,組織有關數據。水文地質概念模型的概化主要包括計算范圍和邊界條件的概化,含水層結構和水文地質特征的概化等[1]。
1)計算區(qū)范圍。根據評價區(qū)水文地質條件特點,結合現有各類資料確定研究區(qū)模擬范圍如下:模型北部邊界是以山區(qū)與平原銜接處,作為弱透水的流量邊界;南部邊界為沙河,定為River邊界;西部以唐家莊——林東一線的地下水分水嶺為界,從歷年地下水流場圖上看,此處等水位線的形狀變化不大,一直為地下水分水嶺位置,因此西部邊界可作為零流量邊界處理;東部邊界參照評價區(qū)多年平均的流場圖定為通用水頭邊界;模擬區(qū)總面積約90 km2,如圖1所示。
圖1 模擬計算區(qū)范圍示意圖
2)含水層的概化。松散巖類孔隙水與巖溶裂隙水是評價區(qū)的兩個主要的地下水類型,兩類地下水之間沒有穩(wěn)定的隔水層,水力聯系密切,可以將水流模型概化為一層,作為非均質各向同性潛水二維平面流處理。實際上對于區(qū)域水流來講,垂向上水力聯系比較密切,多年地下水動態(tài)資料也沒有分層檢測,在垂向上沒有水壓力分布資料驗證模型,只能暫時將各層的水力特征認為是一樣的(表現為等效的水力特征,具有同樣的水文地質參數,統(tǒng)一的地下水流場)。
綜上所述,將評價區(qū)的地下水流系統(tǒng)概化為非均質,各向同性,二維非穩(wěn)定流的地下水流系統(tǒng)。其數學模型為[2]:
式中:Ω-滲流區(qū)域;
H-地下水水位標高(m);
K-含水層在水平方向上的滲透系數(m/d);
ε-含水層的源匯項(m/d);
H0-初始流場(m);
Γ2—滲流區(qū)域的二類邊界;
n-邊界面的法線方向;
q-Γ2邊界上的單寬流量(m2/d),流入為正,流出為負;
Z(x,y)—含水層底板高程。
1)區(qū)域剖分。對整個區(qū)域模型采用矩形網格剖分,剖分為80行100列,邊長為0.1 km×0.1 km,共剖分矩形網格單元8000個,其中有效單元4664個,計算節(jié)點位于單元中心。模擬區(qū)網格剖分見圖2。
圖2 模擬區(qū)網格剖分圖
2)源匯項處理。
(1)大氣降水入滲補給量。潛水含水層通過包氣帶接受大氣降水入滲補給。降水入滲補給條件的不均勻性用入滲分區(qū)概化處理。依據有關降水入滲資料,并參考包氣帶巖性、潛水位埋深、地形、植被等因素,繪出全區(qū)降水入滲系數分區(qū)圖,分別給出各區(qū)降水入滲系數平均值,加在模型對應的剖分網格單元上。根據各區(qū)面積、降水量、降水入滲系數來計算降水入滲補給量。
當降水量較小時,難以補給地下水,所以當月降水量小于10 mm時,不計入有效降水量。
(2)地下水開采量。研究區(qū)內潛水主要是用于居民生活用水和企業(yè)用水。
研究區(qū)內地下水開采有集中開采、農村生活用水分散開采。集中開采為卑家店水廠供水水源地和各個企業(yè)自備井開采,其開采量按實際調查的開采量加在水源地對應的網格節(jié)點上。農村生活用水分散開采,按開采強度進行分區(qū)概化,分別給出各區(qū)開采強度,加在模型對應的剖分網格單元上。
(3)蒸發(fā)。潛水蒸發(fā)是指潛水(埋深小于4米時)在毛細管力的作用下向上運動,最終以參加陸面蒸散發(fā)形式散逸到大氣中的水分損失量。評價期內潛水埋深均超過了4米,潛水蒸發(fā)量按零計[3]。
(4)河流入滲[4]。模型中將沙河作為第三類邊界條件——混合邊界(River)處理。從河流橫剖面示意圖(圖3)上看出,河流通過低滲透性底積物與含水層發(fā)生水力聯系,根據達西定律,橫剖面法線上河流和含水層的交換量示,所以整個河流和含水層的交換量可用公式(1)表示(圖4):
式中:QR為河流側滲量;
hR為河流水位(m);
h為地下水位(m);
W為河流寬度(m);
M為河床底積層厚度(m);
KC為河床底積層厚度滲透系數m/d;
L為河流長度。
因此,河流和含水層之間的交換量由河水水位hR、潛水水位h、底積物的滲透系數Kc、底積物的厚度M以及河流的長度L寬度M共同控制。把各段水位和參數寫入Modflow中的River程序包,便可計算沙河滲漏補給。
圖3 河流橫剖面示意圖
圖4 河床水力傳導的理想化模型示意圖
3)初始條件設置及模擬識別
根據所掌握的資料,本次模擬識別期選為2012年9月1日到2013年3月1日,應力期以月為單位,共劃分為6個應力期,每個應力期又包括若干個時間步長,時間步長為模型自動控制,嚴格控制每次的迭代誤差,在同一應力期內地下水補排項不變。
4)模擬識別
模型的識別過程是整個模擬中極為重要的一步工作,通常要在反復修改參數和調整某些源匯項基礎上才能達到較為理想的擬合結果。此模型的識別過程采用的方法稱為試估—校正法,屬于反求參數的間接方法之一。
為了確保模型求解的唯一性,在模型調試過程中充分利用各種定解條件,也就是用那些靠得住的實測資料來約束模型對原形的擬合。在模型調試過程中,還充分利用水文地質調查中獲得的有關信息及計算者對水文地質條件的認識,來約束模型的調試和識別。
本次模擬首先進行了穩(wěn)定流計算,以便擬合潛水的初始流場(見圖5),這樣做避免了直接建立非穩(wěn)定模型多參數識別的不便,通過建立相對于非穩(wěn)定流模型輸入輸出簡單的穩(wěn)定流模型,運用了模型反求參的方法獲得含水層滲透系數。另外,概化的含水層的結構也在建立穩(wěn)定流模型時確定下來,直接運用于非穩(wěn)定流模型。這樣非穩(wěn)定流模型的參數識別過程就可以只確定給水度的大小,因此增加了此次模型的可信性[5]。
接著用穩(wěn)定流擬合的潛水的初始流場(2012年9月1日流場)作為非穩(wěn)定流模擬的初始值(和實測的初始等水位線比起來,穩(wěn)定流模擬計算得出的流場能更明顯地表現出評價區(qū)的水文地質條件),運行計算程序,通過擬合同時期的流場,識別水文地質參數、邊界值和其它均衡項,使建立的模型更加符合模擬區(qū)的水文地質條件。
模型的識別和驗證主要遵循以下原則:(1)模擬的地下水流場要與實際地下水流場基本一致,即要求地下水模擬等值線與實測地下水位等值線形狀相似;(2)從均衡的角度出發(fā),模擬的地下水均衡變化與實際要基本相符;(3)識別的水文地質參數要符合實際水文地質條件。根據以上三個原則,對模擬區(qū)地下水系統(tǒng)進行了識別和驗證。通過反復模擬、識別后的水文地質參數較好的刻劃了地下水系統(tǒng)的水文地質特征,基本反映了地下水隨時間和空間的變化規(guī)律,使水位擬合誤差較小,達到預期效果。驗證后的平面流場和優(yōu)化調整后確定的參數分區(qū)分別見圖6至圖8,分區(qū)參數值見表1和2。
圖5 2012年9月1日潛水等水位線擬合圖(初始流場擬合)
圖6 2013年3月1日淺層水等水位線擬合圖(驗證)
圖7 包氣帶降水入滲系數分區(qū)圖
圖8 含水層參數分區(qū)圖
表1 潛水含水層水文地質參數分區(qū)表
表2 降水入滲系數參數分區(qū)表
廠區(qū)新鮮供水擬由位于廠區(qū)東北的卑家店水廠水源井(卑家店水廠水源井位于后巍峰山東與沙河之間)來供給,新增開采量規(guī)模為1 235m3/d。
圖9 開采1年后等降深預測圖
為了預測該項目增采地下水對水源地附近區(qū)域流場的影響,利用上述識別和驗證的模型按照既定開采方案,分別對水源地連續(xù)增采20年后的地下水等降深場進行了預測。結果表明:增采1年后水源地漏斗中心降深達1.32米,3年后水源地漏斗中心降深達1.43米,5年后水源地漏斗中心降深達1.46米,且基本達到穩(wěn)定,其降深等值線圖分別見圖9至11。由水源地漏斗中心降深歷時曲線(圖12)與降深等值線圖(圖9至11)分析可知,增采所引起的區(qū)域水位下降并不大,擴展范圍相對比較小。但是鑒于研究區(qū)目前由于多個大型企業(yè)以及居民生活用水對地下水的開采,已經造成了該區(qū)域地下水位的明顯下降,含水系統(tǒng)遭到了一定的破壞,使得該區(qū)域的地下水資源量正在逐步減少,因此建議政府水務行政職能部門嚴格規(guī)范控制該地區(qū)各行各業(yè)的地下水開采量,開展一水多用,開源節(jié)流等節(jié)水措施,并逐步減少對地下水開采,維護地下水采補平衡,同時建議進一步加大外來水的供應,減小水資源供需矛盾。
圖10 開采3年后等降深預測圖
圖11 開采5年后等降深預測圖
圖12 水源地漏斗中心降深歷時曲線
利用通用的地下水模型軟件 Visual Modflow建立了研究區(qū)的地下水流模擬模型。此次數值模擬過程,筆者首先建立了穩(wěn)定流模型,然后以此為基礎建立非穩(wěn)定流模型,這在建立地下水模型時避免了直接建立非穩(wěn)定模型多參數識別(滲透系數、給水度和儲水系數一起識別)的不便。建立的地下水非穩(wěn)定流模型可以有效預測地下水動態(tài),從而對該廠區(qū)開采地下水進行了可靠的預測評價。
[1]薛禹群,李同斌,賈貴庭.地下水動力學[M].北京:地質出版社,1997.
[2]王旭升,萬力.地下水運動方程[M].北京:地質出版社,2009.
[3]遲寶明,王文科,宮輝力,等.專門水文地質學[M].北京:地質出版社,2009.
[4]郭衛(wèi)星,盧國平.模塊化三維有限差分地下水流動模型[M].南京:南京大學出版社,1999:111-113
[5]朱奎.大區(qū)域地下水數值模擬[D].武漢:武漢大學,2004:46-47