黑龍洞泉域巖溶地下水水質(zhì)運(yùn)移模擬分析

馬志敬， 邊 凱， 龐 宇， 劉 博

(1.河北工程大學(xué)水利水電學(xué)院, 邯鄲 056038； 2. 中國煤炭地質(zhì)總局水文地質(zhì)局， 邯鄲 056004; 3.河北工程大學(xué)地球科學(xué)與工程學(xué)院， 邯鄲 056038)

黑龍洞泉域是邯鄲市的主要飲用水源地之一[1-2]，由于長期的不合理開發(fā)利用，造成該泉域的地下水水位下降，形成了地下水位降落漏斗，水質(zhì)變差，出現(xiàn)了許多地質(zhì)問題[3]。地下水污染不同于地表水污染，當(dāng)發(fā)現(xiàn)地下水被污染時其污染程度已經(jīng)很高，很難進(jìn)行恢復(fù),由于人類活動對區(qū)域水資源可持續(xù)利用的影響會隨著人類社會的發(fā)展和進(jìn)步繼續(xù)擴(kuò)大[4]，所以合理開發(fā)利用和有效保護(hù)地下水資源就顯得尤為重要。

建立地下水水源地保護(hù)區(qū)是保護(hù)地下水水質(zhì)安全的有力手段[5]，雖然國內(nèi)對地下水水源地保護(hù)區(qū)的研究主要借鑒國外的3級劃分法[6]。但是隨著研究的不斷深入，中國水源地保護(hù)區(qū)的劃分方法主要形成了經(jīng)驗(yàn)值法、公式計(jì)算法、解析模型法、數(shù)值模擬法等[7]。近年來，巖溶水越來越受到學(xué)者的關(guān)注，常用的巖溶水?dāng)?shù)值模擬的方法主要有等效多孔介質(zhì)法、雙重連續(xù)介質(zhì)法和三重介質(zhì)法[8-9]。地下水?dāng)?shù)值模擬方法正逐漸取代傳統(tǒng)技術(shù),已成為現(xiàn)代研究地下水系統(tǒng)各種問題的重要技術(shù)手段之一[10-11],地下水?dāng)?shù)值模擬方法的各種模型及研究范圍已經(jīng)基本滿足了國民經(jīng)濟(jì)發(fā)展建設(shè)的需要[12]。

利用地下水流數(shù)值模型分析評價地下水資源量的研究有許多，但針對地下水系統(tǒng)動態(tài)演變規(guī)律的區(qū)域研究就相對匱乏。白曉等[13]采用Visual MODFLOW建立峰峰礦區(qū)巖溶地下水流數(shù)值模型，在目前開采條件下，對研究區(qū)地下水資源量及均衡狀態(tài)進(jìn)行分析，并結(jié)合ARIMA降水量預(yù)測模型，預(yù)測了2018—2025年地下水位動態(tài)變化，為峰峰礦區(qū)地下水資源合理開發(fā)利用提供理論依據(jù)和參考；秦歡歡[14]通過建立北京市平原區(qū)三維地下水流穩(wěn)態(tài)及非穩(wěn)態(tài)數(shù)值模型，采用情景分析法模擬了北京平原2015—2030年的地下水流動及水均衡狀態(tài)，對北京平原地下水儲量、水均衡及可持續(xù)利用狀況進(jìn)行一定分析，制訂出最優(yōu)的地下水開采情景，以保證北京地下水資源的合理、可持續(xù)開發(fā)利用，控制地下水環(huán)境的繼續(xù)惡化，為北京地下水的可持續(xù)利用和發(fā)展提供科學(xué)的指導(dǎo)。

史啟朋等[15]和楊釗等[16]都利用Visual MODFLOW軟件建立了污染物運(yùn)移的數(shù)值模型，以模擬和預(yù)測研究區(qū)污染物運(yùn)移過程與趨勢特征，為地下水環(huán)境評價和治理提供了依據(jù)；張文喆等[17]將FEFLOW應(yīng)用于徐州市奎河地下水流數(shù)值模擬，分析了河水與地下水的補(bǔ)排關(guān)系，為研究地下水污染提供了依據(jù)。

肖杰等[18]通過應(yīng)用數(shù)值模擬的方法對崇州市城區(qū)飲用水源地進(jìn)行了各級水源地保護(hù)區(qū)的劃分；吳樂等[19]通過對北京西山地區(qū)進(jìn)行區(qū)域含水層系統(tǒng)數(shù)值模擬，分析了在不同的開采條件下含水層系統(tǒng)響應(yīng)特征，為地下水資源持續(xù)利用提供依據(jù)；孟茜等[20]通過采用地下水?dāng)?shù)值模擬的方法，對準(zhǔn)格爾地區(qū)地下巖溶水資源量變化進(jìn)行了預(yù)測評價，為本區(qū)的地下水合理開采提供了依據(jù)，對類似條件的地下水?dāng)?shù)值模型建立提供了一定的借鑒。研究表明，地下水?dāng)?shù)值模擬法在實(shí)際應(yīng)用中結(jié)果準(zhǔn)確可靠，具有廣泛的實(shí)用性。但在地下水?dāng)?shù)值模擬中結(jié)合遙感影像來確定水文地質(zhì)參數(shù)分區(qū)的研究相對匱乏。

現(xiàn)以黑龍洞泉域巖溶地下水為例，在充分了解地下水?dāng)?shù)值模擬法的理論和解決地下水問題的步驟后[21]，對研究區(qū)建立概念模型和數(shù)學(xué)模型，在已有資料的基礎(chǔ)上結(jié)合遙感影像來確定含水層參數(shù)，對其數(shù)學(xué)模型進(jìn)行識別和驗(yàn)證,通過模擬獲取地下水流場和以開采井為中心溶質(zhì)質(zhì)點(diǎn)向外遷移的運(yùn)動軌跡，按照《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》(HJ 338—2018)[22]要求劃分各級水源地保護(hù)區(qū)，進(jìn)而為有關(guān)部門提供合理的劃分結(jié)果。

1 區(qū)域概況

黑龍洞泉域在河北省邯鄲市峰峰礦區(qū)，西邊緊鄰太行山，華北平原在泉域的東面，總體地勢趨于西北、西部高，東南、東部低。

泉域內(nèi)從老到新發(fā)育的地層依次為：震旦系、寒武系、奧陶系、石炭系、二疊系、三疊系、古近系、新近系和第四系。其中震旦系和三疊系地層僅在局部出漏分布。泉域內(nèi)以新華夏系為主要構(gòu)造體系，構(gòu)造方向主要是北東、北西及東西向，構(gòu)造形態(tài)主要有褶皺和斷裂。

2 地下水流數(shù)值模型構(gòu)建

2.1 水文地質(zhì)概念模型

泉域內(nèi)主要含水巖系有寒武系、奧陶系碳酸鹽巖裂隙巖溶含水巖系，石炭系、二疊系、三疊系碎屑巖夾碳酸鹽巖裂隙含水巖系，第三系、第四系松散巖類孔隙含水巖系三種，其中寒武-奧陶系灰?guī)r作為研究區(qū)的主要含水巖組。含水層中的巖溶裂隙非常發(fā)育，且都為非均質(zhì)各向異性，由西北向東南逐漸變薄，厚度在109～230 m。上覆由石炭系本溪組的鋁土質(zhì)泥巖構(gòu)成了相對穩(wěn)定的隔水頂板，使地下水處于承壓狀態(tài)，但由于局部地區(qū)的人工開采，造成水位下降，地下水已由承壓狀態(tài)轉(zhuǎn)化為無壓狀態(tài)。下伏的寒武系灰?guī)r由于埋藏位置較深，節(jié)理裂隙極不發(fā)育且多被充填，因此形成了相對穩(wěn)定的隔水底板。雖然它們之間都有相對隔水層,但是還會通過各種構(gòu)造來發(fā)生水力聯(lián)系,因此在區(qū)域上可以當(dāng)成是統(tǒng)一的含水巖體。在研究區(qū)內(nèi)除了有集中式開采的羊角鋪水源地外，還有工業(yè)用水井、礦井疏干排水井、分散式開采的鄉(xiāng)鎮(zhèn)工農(nóng)業(yè)及生活用水水源井。

根據(jù)前述并結(jié)合羊角鋪水源地周邊水文地質(zhì)條件進(jìn)行合理概化，東部邊界：以岳城、新坡、中史村一線為邊界，灰?guī)r的埋深大多在標(biāo)高-800 m以下，巖溶已經(jīng)極不發(fā)育，地下水循環(huán)變得滯緩，所以設(shè)為阻水邊界；西部邊界：以由許多相對阻水及導(dǎo)水?dāng)鄬咏M成的賈壁斷裂帶為邊界，該邊界的西面有大面積的碳酸鹽巖裸露區(qū)，在雨季時接受大氣降水補(bǔ)給后，以側(cè)向徑流的形式通過導(dǎo)水?dāng)鄬恿魅朐撓到y(tǒng)，因此設(shè)置為補(bǔ)給邊界;南部邊界：以老爺山背斜及陡貢-楊花莊斷裂隆起帶為邊界，是一個地下水分水嶺，阻水性能良好，將此邊界設(shè)為隔水邊界；北部邊界：以八特至康二城西一線為邊界，其中西段八特至楊二莊一線，由于苑城地塹使奧陶系灰?guī)r埋藏加深，巖溶極不發(fā)育，為阻水邊界，中段由于鼓山斷層將斷層兩側(cè)含水層錯開，同樣為阻水邊界，東段奧陶系埋深較深，且北部有紫山阻水巖體存在，也可視為阻水邊界，總體來說北部邊界為阻水邊界。

將研究區(qū)的地下水流概化為三維非均質(zhì)各向異性的承壓-無壓非穩(wěn)定流，地下水溶質(zhì)通過對流-彌散隨地下水流一起運(yùn)移。水文地質(zhì)概念模型如圖1所示。

圖1 水文地質(zhì)概念模型圖Fig.1 Hydrogeological conceptual model diagram

2.2 數(shù)學(xué)模型

選取北東18°主構(gòu)造線方向-主滲透方向?yàn)閤軸，根據(jù)上述水文地質(zhì)概念模型，建立下列與之相適應(yīng)的數(shù)學(xué)模型。

2.2.1 地下水流數(shù)學(xué)模型

根據(jù)已經(jīng)建立的水文地質(zhì)概念模型，用以下的數(shù)學(xué)模型進(jìn)行描述：

(1)

H(x,y,z)t=0=H0(x,y,z,t0), (x,y,z)∈Ω

(2)

(3)

(4)

式中：kxx、kyy、kzz為含水層各向異性主方向滲透系數(shù)，m/d;H為點(diǎn)(x,y,z)在t時刻的水頭值，m;H0(x,y,z,t0)為點(diǎn)(x,y,z)處的初始水位，m;Ss為貯水率，m-1;W為源匯項(xiàng)，d-1;t為時間，d；Ω為計(jì)算區(qū)；qw為自由面單位面積上大氣降水入滲補(bǔ)給量與地下水蒸發(fā)量之和，m/d;cos(n,x)、cos(n,y)、cos(n,z)分別為流量邊界外法線方向與坐標(biāo)軸方向夾角的余弦；μ為飽和差(當(dāng)自由面上升時)或給水度(當(dāng)自由面下降時)；q(x,y,z,t)為第二類邊界上單位面積的補(bǔ)給量，m/d。

2.2.2 溶質(zhì)運(yùn)移數(shù)學(xué)模型

(5)

C(x,y,z,t)t=0=C0(x,y,z,t0)， (x,y,z)∈Ω,t≥0

(6)

(7)

式中：R為遲滯系數(shù)；θ為介質(zhì)孔隙度；C為組分的濃度，ML-3;T為時間，d;x,y,z為空間位置坐標(biāo)，L；Dij為水動力彌散系數(shù)張量，L2/T;Vi為地下水滲流速度張量，L/T;qs為源和匯，T-1；Cs為源或匯水流中組分的濃度，ML-3;C0(x,y,z,t0)為已知濃度分布；Ω為模型模擬區(qū)域；Γ2為通量邊界；fi(x,y,z,t)為邊界上Γ2已知的彌散通量函數(shù)。

2.3 數(shù)值模型

2.3.1 研究區(qū)空間網(wǎng)格剖分

研究區(qū)空間離散采用等距或不等距的正交長方體對研究區(qū)進(jìn)行網(wǎng)格剖分，從而得到了研究區(qū)的空間離散單元，對于需要提高模擬精度的區(qū)域進(jìn)行局部加密剖分，最終在研究區(qū)平面上剖分成110×70的矩形網(wǎng)格單元，共計(jì)7 700個網(wǎng)格單元，其中有效單元為4 962個，無效單元格2 738個，平面剖分網(wǎng)格和垂向剖分網(wǎng)格如圖2～圖4所示。

圖2 研究區(qū)平面剖分圖Fig.2 Plane section map of the study area

圖3 第32列剖面圖Fig.3 Column 32 profile

圖4 第51行剖面圖Fig.4 Line 51 profile

2.3.2 源匯項(xiàng)處理

根據(jù)研究區(qū)地下水系統(tǒng)的補(bǔ)給、徑流及排泄條件，最終確定大氣降水入滲補(bǔ)給和側(cè)向徑流補(bǔ)給為研究區(qū)的地下水補(bǔ)給來源。大氣降水補(bǔ)給主要分布在中部鼓山和西部廣大的碳酸鹽巖裸露區(qū)，入滲系數(shù)取0.5(根據(jù)河北省邯鄲市水資源研究)。側(cè)向徑流補(bǔ)給主要在研究區(qū)西部邊界通過導(dǎo)水?dāng)鄬恿魅朐撓到y(tǒng)。

地下水排泄以人工開采井和黑龍洞泉群自然排泄為主要途徑，排泄項(xiàng)采用抽水井的方式表示。根據(jù)自來水公司及峰峰礦區(qū)水利局、邯鄲市水利局收集2018年地下水排泄量數(shù)據(jù)，如表1所示。

表1 地下水排泄量統(tǒng)計(jì)表

2.3.3 水文地質(zhì)參數(shù)分區(qū)

已建立的數(shù)值模型把滲透系數(shù)、儲水系數(shù)、給水度等作為主要的水文地質(zhì)參數(shù)。早在20世紀(jì)70年代，專家學(xué)者就已對黑龍洞泉城研究做了大量的工作，如1976—2000年，煤炭部水文勘探指揮部、河北省電力勘測設(shè)計(jì)研究院等對峰峰礦區(qū)奧灰地下水疏放降壓進(jìn)行了水文地質(zhì)勘探，并在羊角鋪和王鳳等地進(jìn)行了多次抽水試驗(yàn)。通過以前多次的水文地質(zhì)勘探和抽水試驗(yàn)，對黑龍洞泉域內(nèi)的水文地質(zhì)條件有了充分了解。故奧陶系灰?guī)r含水層參數(shù)主要依據(jù)20世紀(jì)70年代以來多次水文地質(zhì)勘探成果資料進(jìn)行分區(qū)并賦值。根據(jù)黑龍洞泉域內(nèi)水文地質(zhì)條件差異，對該區(qū)域進(jìn)行參數(shù)分區(qū)劃分，如圖5所示，共分為了14個參數(shù)分區(qū)，并賦予每個參數(shù)分區(qū)不同的給水度初值、儲水系數(shù)Ss和滲透系數(shù)k，來作為模型模擬計(jì)算的基礎(chǔ)。

黑龍洞泉域巖溶水系統(tǒng)補(bǔ)給區(qū)位于西部山區(qū)和中部鼓山灰?guī)r裸露區(qū)，面積約1 262 km2。為減少因模型參數(shù)的不準(zhǔn)確性帶來的誤差，在上述參數(shù)分區(qū)的基礎(chǔ)上結(jié)合遙感影像(圖6)，依據(jù)灰?guī)r裸露區(qū)覆蓋程度對其再進(jìn)行分區(qū)，賦予不同的入滲系數(shù)，如表2所示，最終模擬計(jì)算時采用最準(zhǔn)確的參數(shù)分區(qū)。

2.3.4 模型的識別與驗(yàn)證

模型識別與驗(yàn)證也稱為“參數(shù)反演”，通過野外現(xiàn)場水文地質(zhì)試驗(yàn)獲得的水文地質(zhì)參數(shù)是模型的解析解，與實(shí)際值相差較大，但水位觀測數(shù)據(jù)是比較可靠、準(zhǔn)確的，因此水位觀測數(shù)據(jù)可用于反求水文地質(zhì)參數(shù)。當(dāng)計(jì)算水位與觀測水位趨勢相同時，說明修正后的水文地質(zhì)參數(shù)可反映實(shí)際的水文地質(zhì)條件。

圖5 水文地質(zhì)參數(shù)的分區(qū)圖Fig.5 Zonal map of hydrogeological parameters

Ⅰ區(qū)為灰?guī)r表面覆蓋程度較大區(qū)域； Ⅱ區(qū)為灰?guī)r表面覆蓋程度較小區(qū)域圖6 系統(tǒng)邊界的遙感圖Fig.6 Remote sensing diagram of system boundary

模型選擇2017年1月1日—12月31日為模型的校正、識別期，分為10個制度期，每個制度期又分為10個計(jì)算時間步長。在全區(qū)共建立了20口水位觀測井，其分布如圖7所示。該模型選取了9眼地下水位觀測井，因?yàn)?眼井大致分布均勻，具有長期的水位觀測資料。

通過該模型的運(yùn)行，根據(jù)預(yù)測水頭值和觀測井實(shí)測值的擬合情況，對參數(shù)不斷進(jìn)行試算檢驗(yàn)，直到兩者擬合較好為止，得到各分區(qū)的參數(shù)值。經(jīng)校正、識別，各分區(qū)的參數(shù)識別結(jié)果如表3和表4所示。各地下水位觀測井的水位擬合精度如圖8所示，平均水位誤差在1 m左右。

從上述擬合圖中可以看出，在模型識別的驗(yàn)證期內(nèi)，每層觀測孔的觀測水位基本與計(jì)算出的水位一致，滿足了模型識別和驗(yàn)證的要求。結(jié)果表明，水文地質(zhì)概念模型是合理的，所建立的數(shù)值模型基本反映了模擬區(qū)域地下水運(yùn)動規(guī)律，可用于水源保護(hù)區(qū)的劃分。

表2 泉域補(bǔ)給區(qū)入滲系數(shù)分區(qū)統(tǒng)計(jì)表

圖7 地下水位觀測井分布示意圖Fig.7 Schematic diagram of groundwater level observation well distribution

表3 各分區(qū)參數(shù)表

表4 各分區(qū)參數(shù)數(shù)據(jù)來源表

圖8 不同時間的地下水位擬合圖Fig.8 Groundwater level fitting map of different time

2.4 溶質(zhì)運(yùn)移模型(流線追蹤)

利用上述識別、驗(yàn)證后的數(shù)值模型，以2018年1月1日作為預(yù)測的初始時刻，通過運(yùn)行Visual Modflow模擬，并結(jié)合《中華人民共和國國家環(huán)境保護(hù)標(biāo)準(zhǔn)》關(guān)于地下水源地一、二級保護(hù)區(qū)劃分的原則，模擬得到運(yùn)移100、1 000 d的運(yùn)動軌跡。

受地下水徑流和滲透系數(shù)的影響，質(zhì)點(diǎn)在羊角鋪水源地西部和北部運(yùn)移跡線較長，南部和東部相對較短。在羊角鋪水源地北部，粒子運(yùn)移100 d時北部最大運(yùn)移半徑約為1 120 m，西部約770 m，南部約960 m，東部約618 m；粒子運(yùn)移1 000 d時北部最大運(yùn)移半徑約為3 230 m，西部約3 980 m，南部約2 215 m，東部約950 m。

3 保護(hù)區(qū)范圍的初步確定

保護(hù)區(qū)范圍初步確定其拐點(diǎn)坐標(biāo)系采用2000國家大地坐標(biāo)系(CGCS2000)。

3.1 一級保護(hù)區(qū)范圍的確定

根據(jù)《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》(HJ 338—2018)有關(guān)規(guī)定，利用上述識別、驗(yàn)證后的數(shù)學(xué)模型，獲取以羊角鋪地下水開采井為中心溶質(zhì)質(zhì)點(diǎn)向外遷移的運(yùn)動軌跡,根據(jù)溶質(zhì)質(zhì)點(diǎn)運(yùn)移100 d的距離確定地下水水源地一級保護(hù)區(qū)范圍。是由下列地理坐標(biāo)點(diǎn)所連曲線包圍的范圍，面積2.52 km2，如表5和圖9所示。即羊角鋪水源地地下水開采井為中心距粒子運(yùn)移100 d時北部最大運(yùn)移半徑約為1 120 m，西部約770 m，南部約960 m，東部約618 m，包括羊角鋪村和三河底村轄區(qū)范圍。

表5 一級保護(hù)區(qū)邊界拐點(diǎn)坐標(biāo)

圖9 一級保護(hù)區(qū)示蹤粒子運(yùn)移軌跡及定界圖Fig.9 Tracer particle migration trajectory and boundary diagram in the first class reserve

3.2 二級保護(hù)區(qū)范圍的確定

根據(jù)數(shù)值模型預(yù)測，在劃分完一級保護(hù)區(qū)后，根據(jù)溶質(zhì)質(zhì)點(diǎn)運(yùn)移1 000 d的距離確定地下水水源地二級保護(hù)區(qū)范圍。粒子運(yùn)移1 000 d時北部最大運(yùn)移半徑約為3 230 m，西部約3 980 m，南部約2 215 m，東部約950 m。是由下列地理坐標(biāo)點(diǎn)所連曲線包圍的范圍，面積為26.28 km2，北至山底村北300 m，南至豆腐溝村南220 m，西至宿鳳村西560 m，東至S222省道，如表6和圖10所示。

表6 二級保護(hù)區(qū)邊界拐點(diǎn)坐標(biāo)

圖10 二級保護(hù)區(qū)示蹤粒子運(yùn)移軌跡及定界圖Fig.10 Tracer particle migration trajectory and boundary diagram in the second class reserve

4 結(jié)論

(1)通過運(yùn)用數(shù)值模擬的方法，建立研究區(qū)的水文地質(zhì)概念模型、數(shù)學(xué)模型和數(shù)值模型，為減少因模型參數(shù)的不準(zhǔn)確性帶來的誤差，當(dāng)在模型水文地質(zhì)參數(shù)分區(qū)時結(jié)合遙感影像來確定含水層參數(shù)，在借助Visual MODFLOW軟件模擬地下水流場和溶質(zhì)質(zhì)點(diǎn)的運(yùn)動軌跡，并根據(jù)在不同的時間標(biāo)準(zhǔn)內(nèi)所遷移的距離來界定地下水水源地各級保護(hù)區(qū)的范圍。

(2)以羊角鋪水源地為實(shí)例，進(jìn)行了數(shù)值模擬法的應(yīng)用研究。在地下水水流模型的基礎(chǔ)上，通過流線追蹤技術(shù)，以溶質(zhì)質(zhì)點(diǎn)遷移100 d和1 000 d的距離分別作為一級保護(hù)區(qū)和二級保護(hù)區(qū)的邊界，劃定了水源地保護(hù)區(qū)。

(3)在最終確定保護(hù)區(qū)的邊界時，應(yīng)充分考慮水源地的實(shí)際地形情況，以便確定明顯的標(biāo)志物，使相關(guān)地下水管理部門更容易操作和制訂管理措施。