劉 娟，楊國(guó)勇，蔣兆華，王萬(wàn)林，徐 文

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

1 FEFLOW 簡(jiǎn)介［1］

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

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

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

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

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

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

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

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

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

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

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

2.1 污染源

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

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

近年來(lái)，由于污染源(某農(nóng)藥廠)生產(chǎn)周期變化和有時(shí)故障停產(chǎn)，其外排廢水四氯化碳濃度變化劇烈，最高達(dá)1 390.7，最低為3.1。近年來(lái)由于地下水污染引起重視，外排廢水經(jīng)過(guò)一定的處理，污染源四氯化碳濃度整體有了一個(gè)下降。

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

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

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

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

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

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

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

先建立水流數(shù)學(xué)模型，在水流數(shù)學(xué)模型的基礎(chǔ)上建立溶質(zhì)運(yùn)移模型。

3.1 水流模型

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

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

3.2 溶質(zhì)運(yùn)移模型

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

邊界和初始條件:

式中:C為地下水中的污染物濃度，量綱為ML－3;為水動(dòng)力彌散系數(shù)，量綱為L(zhǎng)2T－1;為地下水的實(shí)際平均流速，量綱為L(zhǎng)T－1;為源匯項(xiàng)所在單位體積含水層的水流通量，量綱為T(mén)－1;為源匯項(xiàng)的污染物濃度，量綱為MT－3;為有效孔隙度，無(wú)量綱;為模擬區(qū);為第一類(lèi)濃度邊界;為平均水力傳導(dǎo)系數(shù)，量綱為L(zhǎng)T－1;

4 模型的識(shí)別與檢驗(yàn)

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

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

圖2 2007年8月實(shí)測(cè)地下水等水位線

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

圖3 2007年8月地下水位模擬結(jié)果

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

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

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

溶質(zhì)運(yùn)移模擬選取時(shí)間步長(zhǎng)為2190天，為6個(gè)完整的日歷年，即模擬始于2001年8月，終止于2007年8月，圖6是兩者的染羽擬合圖。

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

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

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

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

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

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

5 結(jié)論及建議

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

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

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

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

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

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

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

［4］李旭東，曹玉清，胡寬瑢.水文地質(zhì)單元內(nèi)水化學(xué)類(lèi)型形成某些機(jī)制問(wèn)題的探討－以辛安泉域潞安礦區(qū)為例［J］.地球科學(xué)－中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2000，25(2)

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

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

［7］薛禹群，謝春紅。地下水?dāng)?shù)值模擬。科學(xué)出版社。2007

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