基于GMS的長春齊家地下水水源地二級保護區劃分

鮑新華, 于瀚博, 計 量

1.吉林大學新能源與環境學院，長春 130021 2.中國電力工程顧問集團東北電力設計院有限公司，長春 130021

0 引言

水源地是指提供城鎮居民生活及公共服務用水取水工程的水源地域，分為地表水水源地和地下水水源地，包括河流、湖泊、水庫、地下水等。為保護我國居民正常生產生活用水安全，建立水源地保護區是一種成本低廉、效果顯著的措施[1-4 ]。常用的水源地保護區劃分方法有經驗值法、經驗公式法及數值法等[5-8]。經驗值法和經驗公式法總體而言成本較低、操作便捷，但由于這兩種方法一般采用固定的參數值模糊化處理水文地質特征，導致保護區劃分得不夠精確；而數值法擁有結果精確、能夠對含水層進行細致刻畫的優勢，因此本文主要采用數值法，通過建立地下水流及水質模型進行保護區劃分。現有的采用數值法對地下水保護區進行劃分的研究主要包括數值法建立地下水水流模型、利用粒子示蹤法來劃分保護區范圍[9-11]等，如韓京龍等[12]采用數值法對磐石市下水水源地保護區進行了劃分，其模型能夠較為準確地描述研究區含水層，提出的保護區劃分方案實用性較高[13-14]。

齊家水源地是長春市重要的后備水源地和當地農業生活用水的主要來源。由于其臨近石油開采區，隨著油田前期勘探試采工作的進行，原有的水源地一級保護區已經無法解決油田發展與用水安全之間的矛盾。為了保護周邊居民用水環境不受石油污染[1-4]的同時能兼顧油田開采工作的持續有效進行，在一級保護區的基礎上劃分范圍更大的二級保護區就顯得尤為重要。由于研究區臨近石油開采區，傳統的利用粒子示蹤法建立的保護區并不適用；且由于污染質受分子擴散作用影響，在其隨地下水滲流遷移過程中遷移速度要快于地下水流速，因此僅采用水流數值模擬并不足以準確描述污染物的運移范圍。有別于前人以水源地為中心進行溶質運移模擬確定保護區邊界會帶來的不夠細化等問題，本文在地下水水流模型的基礎上，在本研究區采取假定污染源位置反向模擬溶質到達水源井的時間，以期更精準地確定保護區邊界[15-21]。

1 研究區概況

齊家飲用水水源地研究區面積約260 km2，位于長春市雙陽區東北部，北臨石頭口門水庫，南距雙陽城區3 km，西北距長春市區45 km，雙陽河自南向北流經研究區，如圖1所示。

研究區北部平原區地下水資源較為豐富，南部低山丘陵區水資源較貧乏。根據地下水賦存條件和水力特點，區內地下水分為5種類型：第四系松散巖類孔隙潛水，主要沿飲馬河、雙陽河分布；新近系碎屑巖類裂隙孔隙承壓水，主要分布在雙陽區北部平原區；白堊系碎屑巖類裂隙孔隙承壓水，主要分布在飲馬河、雙陽河河谷地帶；碳酸鹽溶隙水，主要分布在低山丘陵地區；基巖裂隙水，主要分布在雙陽盆地東西兩側及南側。

研究區共有集中式飲用水水源地供水井30眼，井深200～220 m，一般單井出水量約40 m3/h，總計服務人口約10萬人。作為長春市第二大地下水水源地，齊家水源地主要為周圍居民提供生活用水及作為長春市后備應急供水水源地，可保證長春及周邊地區供水安全[22]。為應對緩解我國能源危機，1981年對研究區周邊進行全區勘察工作，探明齊家鎮周邊地區石油地質儲量1 017×104t[23]，具有一定的開采價值。然而石油開采井的建設會對周圍地下水水質安全形成一定威脅，保障居民用水安全與石油開采之間的矛盾十分突出，因此科學建立地下水水源地保護區的需求日益迫切[24-26]。2017年5月，長春市環保局雙陽分局根據《飲用水水源保護區劃分技術規范》(HJ/T 338—2007)[22]中的相關技術指標對研究區進行了水源地一級保護區的劃分工作，保護區范圍如圖2所示。新規范(HJ/T 338—2018)[6]施行后，關于一級保護區的規定與(HJ/T 338—2007)相同。

圖1 研究區地理位置圖

圖2 研究區水源地一級保護區范圍

一級保護區的確定，為齊家水源地水質保護提供了基本的技術支撐，但由于其范圍過小，不足以更全面地保護水源地。本文擬在此基礎上，進一步對二級保護區進行研究劃分。

2 水文地質條件概化

通過研究區的地層結構、水動力場、水化學場分析，對研究區進行水文地質條件概化。從采油區到水源地附近的區域為第四系潛水和新近系承壓水的三維流場(圖3)。新近系地層中有多套含水層與相對弱透水層互層分布，由于地層情況較為復雜，已有水文地質資料不足，將其概化為一套含水層。模擬中滲透系數采用了砂礫石參數作為概化后的參數，概化后局部地區會出現溶質沿滲透性較差的介質運移距離較近的情況，在劃分保護區時應將砂礫石參數考慮在內進行適當調整；但宏觀層面上對整體保護區劃分影響不大，因此可將新近系地層中含水層概化為一套含水層。

2.1 確定模型模擬區

模擬區應盡量包含油田開發可能影響到的區域。根據研究區地下水的流向，確定可能影響的區域為下游齊家水源地的準保護區。

本次模擬中，區內河流較多，容易通過地表高程尋找相對獨立的地下水分水嶺。現有水位觀測點繪制的地下水流場分布圖表明了地表與地下水分水嶺具有一致性；因此北部和西北部邊界擴大到伊通河與雙陽河分水嶺，西南部以雙陽河為界，東部以雙陽河與飲馬河的分水嶺為界。如圖4所示。

2.2 含水層概化

研究區地下水類型主要為第四系松散巖類孔隙潛水和新近系碎屑巖類裂隙孔隙承壓水。

第四系松散巖類孔隙潛水含水層包括河谷一級和二級階地含水巖組，地勢低平。含水層巖性具有典型的二元結構特征，上部顆粒較細，主要為亞砂土、亞黏土等，下部顆粒較粗，主要為含砂礫卵石等，厚度總體由南向北漸增，顆粒由粗漸細，水位埋深由深漸淺。其含水層總體埋藏淺，厚度大，透水性較好，有利于鄰區地下水的徑流補給。地表巖性多為亞黏土，局部為亞砂土。區內水稻田多，灌溉回滲水可直接補給地下水。河流直接與第四系潛水含水層相接，地下水與地表水水力聯系密切，水交替頻繁，河水常年或間歇性補給地下水[27]。

第四系之下為新近系碎屑巖類裂隙孔隙承壓水。含水層巖性為粉細砂巖、中粗砂巖、砂礫巖等，頂板埋深20～50 m。承壓水接收上部孔隙潛水垂向和相鄰地下水的側向徑流補給。

通過對研究區地質和水文地質條件分析，將研究區概化為非均質各向同性含水系統，滲流按三維達西流研究。

3 地下水流數值模擬

3.1 地下水流數學模型

根據研究區地質條件與水文地質條件，地下水流數學模型為：

(1)

根據文獻[25]修改。

圖4 模擬區范圍及剖分圖

式中：Kx，Ky，Kz分別為x,y,z方向滲透系數(m/d)；h為地下水水頭(m)；W為源匯項(1/d)；μ為貯水率(1/ m)；t為時間(d)；h0為初始水頭(m)；h1為一類邊界水頭(m)；Σ1為一類邊界；D為研究區。

3.2 研究區剖分

根據研究區水文地質條件和滲流場特征，將模擬區垂向上剖分為3層，分別為第四系潛水含水層、粉質黏土相對弱透水層和新近系承壓含水層。水平方向的網格大小為300 m×300 m，在油井區和水源地等主要模擬區附近加密到50 m×50 m。如圖4所示。

3.3 邊界處理

根據前期收集的水文地質資料，將研究區西北部伊通河與雙陽河分水嶺[27]概化為零流量邊界，將西南部雙陽河概化為水頭邊界，將東部雙陽河與飲馬河分水嶺概化為零流量邊界。

3.4 源匯項處理

模擬區潛水含水層主要補給為降水入滲補給、灌溉回滲水補給和河流側向補給，主要排泄有蒸發排泄、人工開采和側向排泄。承壓水的補給主要為上部孔隙潛水的垂向補給，排泄主要為人工開采。降水入滲系數取0.17～0.23，稻田灌溉回滲強度取0.000 8 mm/d，最大蒸發強度埋深為1.9 m，蒸發極限埋深為4.25 m。

3.5 水文地質參數

模擬用到的參數主要有滲透系數、潛水含水層的給水度、承壓含水層的貯水系數和降雨入滲系數等。各參數初值先通過抽水實驗、土工試驗、水文地質手冊和經驗值綜合給定，再通過模型參數調整，最終獲得模擬采用的水文地質參數。

3.6 模型的識別和驗證

使用數值法軟件GMS(groundwater model system)中的MODFLOW模塊建立地下水流數值模型后，利用給定水文地質參數、邊界和各均衡項數值初值進行模擬，得到概化后的水文地質模型地下水流場空間分布。通過2017—2018年豐水期、枯水期的流場，識別驗證水文地質參數、邊界值和各均衡項，以便建立的模型更好地擬合研究區的水文地質條件。模型識別擬合主要原則為：識別后的水文地質參數與其他資料獲取的水文地質參數基本一致；模擬的地下水流場與流場基本一致；邊界條件、各水均衡項與實際情況基本一致。

模擬采用2017—2018年豐水期的觀測水位進行模型識別，取枯水期的觀測水位進行驗證。通過反復識別和驗證，最終確定上部第四系潛水含水層的滲透系數為6.0～18.0 m/d、給水度為0.1～0.2，中間層粉質黏土相對弱透水層的滲透系數為0.1～0.2 m/d，下部新近系承壓含水層的滲透系數為2.0～14.0 m/d、彈性釋水系數為0.003～0.008，其中滲透系數(K)分區見圖5。圖5中參數區范圍與圖4不完全一致(西北側)，主要是參數區為水源地所在區域，資料較多；而西北側資料有限，且考慮采用伊通河與雙陽河分水嶺邊界較好，故西北側所用參數參照鄰區給出。圖6為水源地所在區域上模型計算豐水期流場與實測流場的對比，擬合情況較好，水位相對誤差為5%～10%，說明模型建立基本合理。

采用2017—2018年枯水期的觀測水位進行模型驗證，圖7為模型計算的水源地所在區域上枯水期流場與實測流場的對比。計算結果表明：相對誤差在1%～5%之間，擬合情況較好；地下水流動方向基本與河流方向一致，并且在主要的飲用水開采區形成明顯的降落漏斗，地下水流速明顯加快，可能會提高污染物遷移速率，使其迅速進入飲水井，從而對周邊居民飲用水安全造成威脅。選取2018年1月—2019年12月，對代表性觀測井W18井和齊家觀測井潛水水位變化計算值與實測值進行擬合觀察，上述流場擬合檢驗結果如圖8所示，結果表明水位擬合相對誤差均不超過5%。經過模型驗證，表明本次模擬具有較高的精度和可信性。

圖5 模擬區潛水含水層(a)和承壓含水層(b)滲透系數分區

4 溶質運移模擬

以上述地下水流數值模型為基礎，通過使用GMS中MT3DMS模塊建立溶質運移模型來進行飲用水水源地保護區圈定。

4.1 地下水溶質運移模型

4.1.1 地下水溶質運移數學模型

地下水溶質運移的對流-彌散方程為

(2)

4.1.2 評價目的層

油田開采過程中發生泄漏，發生泄漏的位置可能是油井的各個層位。研究區的主要地下水類型是第四系孔隙潛水和新近系碎屑巖類裂隙孔隙承壓水。如果泄露事故發生在第四系亞黏土層或者新近系泥巖層，因為這兩個層位的滲透性較差，污水的泄漏量將會很小，不會對齊家水源地構成威脅；如果泄露發生在第四系潛水含水層，采油井在一定的壓力下可向潛水含水層釋放一定的污水，污水有可能通過由南向北方向的水平運移，將污染物運移到齊家水源地所開采的位置，再通過垂向滲透對水源地的承壓含水層水源造成影響；如果泄露發生在新近系承壓含水層，采油井在一定的壓力下會向新近系含水層釋放較多的污水，污水有可能通過由南向北方向的水平運移，將污染物運移到齊家水源地所開采的位置，對水源地的新近系承壓水水源造成污染。

綜上，評價目的層擬選擇發生泄漏事故造成污染較為嚴重的第四系潛水含水層和新近系承壓含水層，以齊家水源地開采的新近系承壓含水層為主。

4.1.3 確定污染運移模型參數

根據筆者前期彌散試驗結果得知，研究區縱向彌散度為0.560 m，橫向彌散度為0.056 m ，垂向彌散度為0.056 m。模擬層主要為中砂，取有效孔隙度為0.3。

圖6 模擬區豐水期潛水水位(a)和承壓水位(b)實測值與計算值擬合圖

4.2 模擬過程

石油類中不易揮發的大分子組分與砂的吸附關系多為線性吸附關系；但考慮到研究區主要含水層為砂礫層，為保守起見，模擬預測中并沒有考慮到這種吸附的阻擋作用。石油類液體本身的黏滯系數較大，在含水層介質中不易流動，但考慮到評價的保守性，在模型中可不加以計算。故模型可采用保守型污染質計算方法進行泄漏事故發生后的污染預測。模擬1 000 d后污染物在承壓含水層中分布見圖9。

由圖9可見，如果位于水源地內的采油井發生泄漏，對水源地的影響將是明顯的。被污染的地下水將很快進入水源井污染飲用水水質，開采井附近受開采影響，具有較大的水力坡度，使得地下水流速加快，污染質運移速度也將加快。根據數值模擬結果，一旦發生泄漏事故，石油類污染質在地下水中經過 1 000 d的運移距離在400～800 m之間。

根據圖9，在水源地周邊距水源井400～800 m范圍內布置多個假想的污染點源，污染點源為假想石油開采井，石油開采井穿過承壓含水層，污染物通過套管部分滲漏入含水層。根據研究區石油開采井特點，采用5%作為采出液滲漏比，利用已有的開采井年開采量及采出液含水率資料可計算出模擬年份污染源向地下水滲漏石油污染物速率為0.25 m3/d、質量濃度為6.2×104mg/L，以此為源強進行預測。

圖7 模擬區枯水期潛水水位(a)和承壓水位(b)實測值與計算值擬合圖

圖8 研究區W18井(a)和齊家觀測井(b)水位擬合圖

圖9 石油開采井泄漏污染物1 000 d后在承壓含水層分布

利用先前建立的溶質運移模型進行模擬，預測污染發生1 000 d后的最大污染邊界。根據污染暈最大邊界與水源井之間的位置關系，在先前劃定的范圍內移動污染源的位置，再次進行模擬，直至污染暈在1 000 d的最大邊界剛好到達水源井。由于含水層污染難以治理，石油類污染物將長期被含水層介質吸附并緩慢釋放，因此上述模擬得出的位于水源井上游的污染源位置應再向水源地以外移動100 m左右留作緩沖區，此時污染源的位置即為二級保護區的邊界。污染源最終位置1 000 d污染暈范圍如圖10所示。

從模擬結果來看，1 000 d時明顯沿地下水流方向污染暈擴散范圍更大，證明石油類污染物在含水層中運移主要受地下水流的影響。此外污染物同時也向其他方向運移，說明污染物受水動力條件影響外，本身分子擴散也起作用。

最終規劃齊家水源地二級保護區范圍如圖11所示，面積為33.5 km2(含一級保護區面積8.5 km2)。

由最終劃定的水源地二級保護區可見，齊古1、昌25、昌31、昌37、昌88、昌90、昌91等部分規劃石油開采井位于二級保護區之內，會對水源地水質安全構成直接威脅，建議取消規劃、停止建設；昌36-4、昌39、昌84、昌87、昌105-11等規劃石油開采井位于二級保護區邊界附近，存在對水源地水質安全造成破壞的風險，建議將開采井規劃位置向二級保護區外遷移500 m以上距離以降低石油勘探、開發過程中對地下水破壞的風險。

圖10 1 000 d污染暈范圍

5 討論

文中采用的研究方法對劃分一級保護區同樣適用，文中涉及到的一級保護區范圍采用長春市環保局相關資料，僅做模擬計算驗證。本文主要對研究區地下水水源地二級保護區劃分進行研究，研究中有關問題進一步說明或討論如下：

1)主要模擬層為新近系，為多套含水層與相對弱透水層互層，模擬中概化為一套含水層。在分層細化資料缺失時一般采用在模擬中概化為一套含水層[28]。這種概化與實際情況有一定出入。污染物運移時，一般在透水性強的地層中運移快[29]。模擬中滲透系數采用了砂礫石參數作為概化后的參數，計算的模擬污染暈范圍總體較實際結果可能偏大，這考慮到了污染物遷移的最大可能。這種概化在資料有限及解決主要污染物最大遷移問題研究中，精度應該是有保證的[29]。

2)水源地西北側非自然邊界，邊界條件不好給定，研究將模擬區邊界擴展到伊通河與雙陽河分水嶺天然邊界(模擬區邊界見圖4)。但模擬區西北部資料較為有限，因此進行模型運算時相關水文地質參數等資料參照水源地鄰近區給出，這與實際相比會有一定誤差。而模擬的重點區在水源地區域，模型調參與校正主要看水源地所在區域。這樣的處理對水源地區域污染問題研究影響較小。

圖11 水源地二級保護區范圍

根據圖6—圖7流場圖，地下水流方向基本沿主徑流方向由西南向東北，位于水源地上游的石油開采井對水源地影響更大，因此在劃分水源地二級保護區時水源地上游劃分范圍應較下游更大，通過溶質運移模型確定的二級保護區范圍也證實了這一點。本文最終劃定水源地二級保護區符合這一基本判斷。

6 結論與建議

1)在地下水水流模型的基礎上使用MT3DMS模塊建立研究區溶質運移模型，在距水源地400～800 m范圍設定若干污染源，1 000 d污染暈最大邊界距水源井30～80 m的污染源位置定為二級保護區邊界，劃定二級保護區面積為33.5 km2(含一級保護區)。

2)建議齊古1井等規劃建于二級保護區范圍內的7口油田開采井取消建設，昌36-4井等規劃建于二級保護區邊界的5口油田開采井向保護區邊界外遷移500 m以上距離；并建立完善的地下水污染監控系統體系和應急響應措施，最大限度地降低石油勘探、開發過程中風險事故對地下水環境的影響。

3)建立在地下水水流模型基礎上的溶質運移數值模型與解析法相比，能夠更加精細地刻畫溶質運移細節，進行更精細的水源地保護區劃分。本文的方法對其他類似問題的研究有借鑒意義。但同時應注意到，模型精細的刻畫程度主要取決于建模依據的地質-水文地質基礎的正確認識。沒有大量的基礎數據支撐，僅靠模型剖分的細化無助于精度的提高。油田運行過程中，還應注重實際污染物運移監測工作，一方面檢驗本文的計算結果，另一方面在資料進一步積累后，可對模型再進一步修改完善。