劉 娟，楊國勇，蔣兆華，王萬林，徐 文

(1.華東冶金地質勘查研究院，安徽 合肥230088;2.中國礦業(yè)大學，江蘇 徐州221008;3.江蘇省工程勘測研究院有限責任公司，江蘇 揚州225009)

1 FEFLOW 簡介［1］

FEFLOW是迄今為止功能最為齊全的地下水水量及水質計算機模擬軟件系統(tǒng)。該軟件具有圖形人機對話、地理信息系統(tǒng)數(shù)據接口，自動產生空間多種有限單元網、空間參數(shù)區(qū)域化以及快速精確的數(shù)值算法和先進的圖形視覺化技術等特點。

在FEFLOW系統(tǒng)中，用戶可以很方便迅速地產生空間有限單元網格，設置模型參數(shù)和定義邊界條件，運行數(shù)值模擬以及實時圖形顯示結果與成圖。

它可以模擬污染物在地下水中遷移過程及其時間空間分布規(guī)律(分析和評價工業(yè)污染物及城市廢物堆放對地下水資源和生態(tài)環(huán)境的影響，研究最優(yōu)治理方案和對策)。

通過標準數(shù)據輸入接口用戶既能直接利用已有的 空間多邊形數(shù)據生成有限單元網格。也可以用鼠標設計和調整網格幾何形狀，增加和放疏網格密度。在建立水流場和遷移模型時，用戶不僅能夠視具體情況定義第一、第二和第三類邊界，而且可以對邊界條件增加特定的限制條件，以避免非現(xiàn)實的數(shù)值解。用戶也能夠直接定義多含水層中的抽水和注水井邊界條件。所有邊界及附加條件既可設置為常數(shù)，也能定義為隨時間變的函數(shù)。已知的邊界及模型參數(shù)可按點、線或面的形式直接輸入。對散的空間抽樣數(shù)據進行內插或外數(shù)據區(qū)域化 ，提供克里格(Kriging)、阿基瑪(Akima)和距離反比加權法(IDW)。輸入數(shù)據格式既可以是ASCⅡ碼文件，也可以是GIS地理信系統(tǒng)文件。FEFLOW支持ARC/INFO點、線、面的廣義數(shù)據格式。

FEFLOW具有齊全的地下水模擬功能，表現(xiàn)在:

(1)三維空間模型，二維平面，二維剖面或者軸對稱二維模型。

(2)非穩(wěn)定流或穩(wěn)定流模擬。

(3)多層自由表面含水系模擬，包括滯水模擬。

(4)化學物質遷移及熱傳遞模擬，包括溫度鹽分 遷移模擬。

FEFLOW的計算結果既有水位，污染物濃度及溫度等標量數(shù)據，也包括流速，流線和流徑線等向量數(shù)據。模型參數(shù)和計算結果既能按 ASCII碼文件，GIS地理信息系統(tǒng)文件，DXF或 HPGL文件輸出，又能在 FEFLOW系統(tǒng)中直接顯示。FEFLOW提供了其他任何地下水模擬軟件都無法比擬的，豐富實用的圖形顯示和數(shù)據結果分析工具。其先進的圖形可視化及數(shù)據分析技術表現(xiàn)在:①有限單元網格、邊界條件和模型參數(shù)的三維可視化;三維彩色等勢面顯示以及二維平面彩色或等值線顯示;②三維地下水流線追蹤，流動時間及流速動畫顯示;③三維交叉斷面圖、剖面圖與切片圖顯示;④三維圖形的交互旋轉、放大或縮小;⑤模型整體和局部水量均衡分析(包括任意幾何多邊形內的水流通量分析);⑥各種邊界條件上的水流通量、物質通量以及在特定時間、空間內的水流積分量都可以由模型算出并以圖形顯示出來。

2 四氯化碳在巖溶含水層中運移的數(shù)值模擬

2.1 污染源

我國北方某缺水城市，其南郊水源地位于一北東向山間盆地中，主要含水層為奧陶系灰?guī)r，基巖埋深0－30 m，主要含水層段埋深70－210 m，近幾年來的檢測發(fā)現(xiàn)該水源地收到四氯化碳的污染，通過污染源調查判定污染源位于水源地補給區(qū)的一農藥廠，該農藥廠生產過程中使用了大量四氯化碳作為溶劑［2］，

四氯化碳(CCl4，Carbon tetrachloride)，熔點 22.96 oC，沸點 76.8 oC，比重(20 oC)1.595，無色有芳香氣味的透明液體，極易揮發(fā)，不燃燒，遇火可分解為二氧化碳、氯化氫、光氣及氯氣。四氯化碳化學性質相當穩(wěn)定，很難分解。在地下水中，揮發(fā)擴散極少，也幾乎沒有生物分解，很難自凈消除。四氯化碳在環(huán)境中具有持久性、長期殘留性和生物蓄積性，可經呼吸道、皮膚、消化道進入人體。它是一種具有“三致”作用的有毒有害物，能傷害肝臟，因而受到嚴格的控制［3］。

近年來，由于污染源(某農藥廠)生產周期變化和有時故障停產，其外排廢水四氯化碳濃度變化劇烈，最高達1 390.7，最低為3.1。近年來由于地下水污染引起重視，外排廢水經過一定的處理，污染源四氯化碳濃度整體有了一個下降。

2.2 水文地質概念模型

本區(qū)含水系統(tǒng)主要包括第四系松散沉積物孔隙含水層和裂隙巖溶含水層［4］。裂隙巖溶含水層，一部分是裸露的巖溶丘陵，另一部分為淺埋于沖洪積群之下的巖溶區(qū)。巖溶含水層在裸露半裸露區(qū)直接或間接接受大氣降水或地表水體的入滲補給，在隱伏區(qū)接受上覆孔隙水的越流補給。由于在巖溶含水層中存在非可溶巖相對隔水地層的分隔，盡管存在裂隙、斷層切割，各層間還有一定的水力聯(lián)系，但水力聯(lián)系較弱。所以仍可將巖溶含水層系統(tǒng)概化為單一的巖溶地下水系統(tǒng)處理。該地下水系統(tǒng)在裸露或半裸露區(qū)為潛水，在覆蓋區(qū)為承壓水。由于研究區(qū)巖溶地層巖性及構造控水特征明顯，巖溶發(fā)育程度及富水程度極不均勻，滲透特征在不同方向變化較大，巖溶發(fā)育往往呈帶狀分布，而在廢黃河斷裂帶中巖溶發(fā)育及滲透特征同時受地層巖性及斷裂構造控制，在綜合作用的情況下，表現(xiàn)為有時沿地層走向發(fā)育，有時沿斷裂方向發(fā)育，有時則方向性不明顯，兩者接近。同時，計算區(qū)巖溶裂隙水一般富集在淺部(在斷裂帶或巖體接觸帶可達200 m以下)，開采井基本也在此深度［5］。因此，本區(qū)巖溶水文地質系統(tǒng)的結構及水動力條件可概化為非均質各向異性的承壓—無壓平面二維流。

2.3 研究區(qū)邊界的確定

研究區(qū)屬于賈汪復式向斜水文地質區(qū)中的七里溝裂隙巖溶水亞區(qū)的一部分，位于徐州市的東南郊。根據含水巖組性質，確定模擬范圍。在平面上，總面積103.29 km，北以廢黃河故道為邊界，南以銀山—三堡一線為邊界，東以邵樓—徐村一線為邊界，西以云龍山—泉山—龍腰山為邊界。根據研究區(qū)地形圖，運用Sufer軟件繪制模擬范圍。

北部的廢黃河故道，由于富水性較好，導水系數(shù)高達9 500 m/d，設置為第一類水頭邊界;南部由于向斜樞紐抬高部位而形成地下水分水嶺，設置為第二類邊界，水量交換為零;東側為北東向山鏈頂部分水嶺，大致沿邵樓斷層走向，設置為第二類零流量邊界;西部以云龍山—泉山—龍腰山一線七里溝盆地西部分水嶺為邊界，主要地層為毛莊組—徐莊組多由砂頁巖組成，其透水性極差，單孔涌水量小于50 m/d，起阻水作用，可作為隔水邊界，設置為第二類零流量邊界;各抽水井作為模型的第四類邊界條件，賦抽水量值［6］，邊界條件約束見圖1。

圖1 研究區(qū)邊界條件示意圖

3 數(shù)學模型建立及求解

先建立水流數(shù)學模型，在水流數(shù)學模型的基礎上建立溶質運移模型。

3.1 水流模型

根據地下水的水均衡原理和達西定律，結合對研究區(qū)的水文地質條件以及初始水頭、邊界條件(東、西兩側為隔水邊界;南部為地下水分水嶺;北部為定水頭邊界)，概化為非均質各向異性介質中的二維非穩(wěn)定流問題，可建立下列承壓含水層地下水水流的二維數(shù)學模型［7］:

式中:T=KM為稱為導水系數(shù)，量綱為L2T－1;物理含義是水力坡度等于1時，通過整個含水層厚度上的單寬流量。W=Wp+Ws－Wg－Q為源、匯項。其中為 Wp降水入滲補給強度;Ws為灌溉入滲補給強度;Wg為蒸發(fā)強度;Q為開采強度，量綱為T－1;μ為多孔介質的貯水系數(shù)(或釋水率)，表示當水頭下降一個單位時，從單位體積孔隙介質中貯存(或釋放)的水量，量綱為 T－1;H0(x，y)為水頭，量綱為L;φ(x，y)為第一類邊界上水頭，量綱為L;Ω為研究的空間區(qū)域。

3.2 溶質運移模型

大量研究表明四氯化碳在地下水中化學性質相當穩(wěn)定，揮發(fā)擴散作用較弱，在自然條件下幾乎沒有生物分解，很難自凈消除［8］，因此可以不考慮化學反應項。非均質介質中水流和溶質運移，對流為主要作用，彌散為次要作用。常用的是對流—彌散方程，該方程如下表示［5］:

邊界和初始條件:

式中:C為地下水中的污染物濃度，量綱為ML－3;為水動力彌散系數(shù)，量綱為L2T－1;為地下水的實際平均流速，量綱為LT－1;為源匯項所在單位體積含水層的水流通量，量綱為T－1;為源匯項的污染物濃度，量綱為MT－3;為有效孔隙度，無量綱;為模擬區(qū);為第一類濃度邊界;為平均水力傳導系數(shù)，量綱為LT－1;

4 模型的識別與檢驗

模型參數(shù)的識別采用間接法。將巖溶含水層按照水文地質參數(shù)分成5個區(qū)，假定每個分區(qū)具有相同的參數(shù)，如各向滲透系數(shù)、貯水系數(shù)、補給模數(shù)等水文地質參數(shù)，先給每個分別賦參數(shù)初值，代入有限元模型中，通過調整參數(shù)值，將數(shù)值模擬的等水位線與實測等水位線擬合，進行模型校正。

本次模擬用實測數(shù)據繪制的2007年8月等水位線來擬合模擬結果，用于水流模型的校正，即選取模擬的時間步長為730天，為2個完整的日歷年，模擬開始時間為2005年8月，終止時間為2007年8月，模擬結束時巖溶水水位見(圖2)。由2007年實測地下水等水位線(圖3)可知研究區(qū)南部與北部的水位幾乎達到同等高度，東西兩側水位低于南北兩端，這顯然與實際的水位不符，而且東西兩側是山地，地下水等水位線走向應該與山的等高線基本一致，因此水流模型需要進一步修正。模型將研究區(qū)南端地下水分水嶺設為第二類邊界條件無流量交換，經過分析相關資料得出在南端有部分邊界存在流量交換，因此在南端打開部分邊界設為第二類水頭邊界條件，對水頭進行賦值。其它具體參數(shù)的小幅調整在此不再冗述。

圖2 2007年8月實測地下水等水位線

將修改條件后模擬出的2007年8月地下水等水位線與2007年8月實測數(shù)據繪制的地下水等水位線進行擬合，結果見圖4。從擬合圖上可以看出，中心漏斗區(qū)范圍基本一致，等水位線走向也基本一致，中間盆地區(qū)域水位較低，東西兩側水位高且等水位線與山的等高線走向一致，符合研究區(qū)的地下水水文地質條件。

圖3 2007年8月地下水位模擬結果

溶質運移模型是在水流模型的基礎上進行的擬合，通過初始條件和污染物邊界、污染物遷移的物質參數(shù)確定。研究區(qū)內較大范圍的巖溶地下水已監(jiān)測到四氯化碳污染。因此污染背景值不能假定為零。本次模擬由于缺乏初始污染資料，這里采用2001年8月豐水期四氯化碳濃度等值線為背景值，即以2001年8月為模擬的開始時間。這里將研究區(qū)2001年8月四氯化碳濃度等值線作為背景值插入，插入結果見圖 5。

圖4 2007年8月地下水等水位線擬合圖

圖5 研究區(qū)CCl4背景濃度圖

溶質運移模擬選取時間步長為2190天，為6個完整的日歷年，即模擬始于2001年8月，終止于2007年8月，圖6是兩者的染羽擬合圖。

從兩者的擬合圖可以明顯地看出:南部污染源區(qū)和北部匯集區(qū)四氯化碳濃度等值線擬合較好，研究區(qū)中部濃度較南區(qū)和北區(qū)偏低，僅有一小部分沒有在實測四氯化碳濃度等值線范圍內，但大體趨勢擬合較好，可以在這個模型的基礎上進行預測研究區(qū)未來四氯化碳濃度分布變化情況。本論文預測了2008年8月、2010年8月、2012年 8月、2020年 8月、2030年8月、2050年8月四氯化碳污染羽分布情況，見圖7－12。有以上各圖可以直觀地總結出以下幾點:

(1)研究區(qū)北部和南部兩個研究區(qū)污染中心的形狀和范圍基本一致，研究初期的污染形狀可以反映出污染濃度中間低，南北兩端高，較好地體現(xiàn)了條帶狀“啞鈴型”分布特征。

(2)巖溶含水層中的四氯化碳污染呈逐年衰減的趨勢，尤其是北區(qū)四氯化碳污染羽在逐漸消失。南區(qū)污染羽也有一定的衰減，但相對于北區(qū)衰減層度稍小。中區(qū)污染羽不明顯。

圖6 2007年8月四氯化碳污染羽擬合圖

圖7 2010－2050年8月巖溶含水層中CCL4污染羽分布圖

(3)污染羽衰減幅度隨著時間的推移在減小，也就是隨著污染羽濃度的衰減，其衰減幅度也在降低，尤其南部污染源區(qū)和北部匯集區(qū)到了最后，其濃度的變化幾乎為零，中區(qū)本身污染較低，其污染羽濃度變化亦不明顯。這種演變規(guī)律正好吻合了四氯化碳的物理化學性質，即當濃度降低到一定程度后，其再降解就會十分困難，這也是地下水四氯化碳污染處理的難點之處。

5 結論及建議

本文的研究成果主要表現(xiàn)為以下幾個方面:

(1)通過對水文地質條件的概化，建立了水文地質概念模型。確定了模擬區(qū)的范圍和計算目的層，概化了含水層水力特征、垂向邊界和側向邊界，闡明了地下水的補給、徑流、排泄規(guī)律和地下水動態(tài)特征。依據水文地質概念模型建立了地下水數(shù)學模型，采用FEFLOW軟件求解確定性數(shù)學模型。

(2)建立七里溝水源地地下水水流數(shù)值模擬模型，利用多年實測資料對模型進行識別和驗證，識別和驗證效果良好，表明含水層結構、水文地質參數(shù)的確定、邊界條件的概化、源匯項的處理是合理的。因此可以在此基礎上進行地下水溶質運移模擬。

(3)在水流模型的基礎上，建立七里溝水源地地下水水質模型，以四氯化碳作為模擬因子，通過研究區(qū)四氯化碳初始濃度、邊界條件、所需參數(shù)的調整進行模擬。結果表明，研究區(qū)巖溶地下水四氯化碳污染羽呈南北濃度高，中間濃度低的“啞鈴型”;研究區(qū)地下水四氯化碳污染濃度隨著時間的推移呈下降趨勢，南區(qū)污染源處濃度降低最明顯，其次是北區(qū)廢黃河處，中區(qū)也有相應的下降，這和研究區(qū)污染特征，地下水文地質結構，污染源的治理和人為因素有關。

［1］谷源澤，張勝紅，郭書英等.FEFLOW有限元地下水流系統(tǒng).中國礦業(yè)大學出版社.2001.8

［2］韓保平，朱雪強等。某巖溶水源地四氯化碳污染途徑研究.中國礦業(yè)大學，2006.1

［3］朱承駐，于勇等。水相中CCL4和CHCL3的紫外光解機理。環(huán)境科學，2001

［4］李旭東，曹玉清，胡寬瑢.水文地質單元內水化學類型形成某些機制問題的探討－以辛安泉域潞安礦區(qū)為例［J］.地球科學－中國地質大學學報，2000，25(2)

［5］徐州市區(qū)地下水水質演化及污染控制研究研究報告.2008.12

［6］徐州市規(guī)劃區(qū)巖溶地下水資源管理研究報告.1999.9

［7］薛禹群，謝春紅。地下水數(shù)值模擬?？茖W出版社。2007

［8］R D Morrison.Application of forensic techniques for age dating and source identification in environmental litigation.Environ.Forensics，2000(1):131－153