楊元旭，柴建華，徐世光,

(1.昆明理工大學(xué)國(guó)土資源工程學(xué)院，昆明 650093；2.云南地礦工程勘察集團(tuán)公司，昆明 650041)

我國(guó)地下水污染正面臨著由點(diǎn)到面、由淺到深、由城市到農(nóng)村不斷擴(kuò)展和污染程度日益嚴(yán)重的趨勢(shì),地下水污染防治迫在眉睫[1]。對(duì)于地下水污染問(wèn)題的研究，目前主流的方法是通過(guò)計(jì)算機(jī)軟件建立地質(zhì)模型，以此來(lái)模擬預(yù)測(cè)地下水污染狀況[2]。GMS軟件在地下水溶質(zhì)運(yùn)移模擬方面應(yīng)用已比較成熟[3]，可用于處理本次研究問(wèn)題。

本文以某污染廠區(qū)受污染地下水為研究對(duì)象，利用GMS軟件建立了廠區(qū)區(qū)塊的地下水流場(chǎng)，并在分析廠區(qū)污染源的基礎(chǔ)上模擬驗(yàn)證了廠區(qū)地下水污染過(guò)程，最終依據(jù)驗(yàn)證的溶質(zhì)運(yùn)移模型分析研究了污染場(chǎng)地特征污染物在松散層地下水中的污染規(guī)律。為其他類(lèi)似性質(zhì)的場(chǎng)地在工業(yè)選址污染防控以及對(duì)已污染場(chǎng)地的地下水修復(fù)提供可借鑒的研究方法和參考依據(jù)。

1 研究區(qū)水文地質(zhì)條件

研究區(qū)位于一個(gè)斷陷盆地內(nèi)，海拔在1 280～1 350 m間，盆地內(nèi)地勢(shì)相對(duì)平坦，年平均降雨量815.8 mm，屬于亞熱帶季風(fēng)氣候。盆地內(nèi)主要沉積第四系地層。區(qū)域水文地質(zhì)圖如圖1所示。

研究區(qū)水文地質(zhì)條件較為復(fù)雜，含水層(組)類(lèi)型展布與區(qū)域主體構(gòu)造密切相關(guān)，而次一級(jí)構(gòu)造主要控制著地下水的運(yùn)移和富集[4]。盆地周邊及盆地基底主要發(fā)育碳酸鹽巖，且?guī)r溶十分發(fā)育，富水性較好，是區(qū)內(nèi)主要的含水層。

第四系松散層主要接受大氣降雨補(bǔ)給，同時(shí)盆地周邊的巖溶水及裂隙水會(huì)通過(guò)巖溶通道及裂隙通道側(cè)向徑流補(bǔ)給盆地第四系孔隙水。地下水主要以徑流方式向盆地內(nèi)最低點(diǎn)長(zhǎng)橋海和大屯海排泄；其次，地下水還通過(guò)地表蒸發(fā)進(jìn)行排泄，但是盆地地下水主要為承壓水，其蒸發(fā)強(qiáng)度較弱。

圖1 區(qū)域水文地質(zhì)圖

2 第四系地下水流場(chǎng)特征

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

通過(guò)場(chǎng)地水文地質(zhì)鉆探、民井調(diào)查以及已有資料分析，本次研究將盆地內(nèi)第四系地層從上至下依次分為粉質(zhì)黏土層(弱透水層)，中層砂礫石層(含水層)，下層粘土層(相對(duì)隔水層)。

盆地西部、南部、東部地區(qū)主要發(fā)育碳酸鹽巖類(lèi)灰?guī)r、白云巖，富水性較強(qiáng)。這些碳酸鹽巖地層巖溶裂隙都十分發(fā)育，且這些地層直接與上覆第四系地層接觸，盆地周邊的巖溶裂隙水會(huì)通過(guò)巖溶裂隙通道補(bǔ)給盆地第四系地層的孔隙水，故本次將盆地周邊碳酸鹽巖發(fā)育地區(qū)概化為向盆地內(nèi)第四系補(bǔ)給的定流量邊界；盆地東南部及東部部分地區(qū)發(fā)育有砂巖地層，其富水性較弱，本次研究將其概化為隔水邊界；大屯海和長(zhǎng)橋海兩個(gè)天然湖泊為盆地內(nèi)的最低水平面，盆地內(nèi)第四系地下水在接受大氣降雨補(bǔ)給和盆地周邊基巖區(qū)側(cè)向補(bǔ)給后，以徑流方式最終排泄于兩個(gè)天然湖泊內(nèi)。因此將盆地內(nèi)兩個(gè)天然湖泊概化為定水頭邊界(詳見(jiàn)圖1)。

2.2 研究區(qū)數(shù)學(xué)模型

分析整個(gè)盆地的區(qū)域水文地質(zhì)條件，本次研究將研究區(qū)第四系地下水流概化為穩(wěn)定的均質(zhì)各項(xiàng)異性的地下水流系統(tǒng)，并選擇三維承壓水含水層地下水穩(wěn)定流運(yùn)動(dòng)方程進(jìn)行概化分析，見(jiàn)式(1)。

(1)

其中，Kxx為x方向滲透系數(shù)主值，m/d；Kyy為y方向滲透系數(shù)主值，m/d；Kzz為z方向滲透系數(shù)主值，m/d；Us為彈性釋水系數(shù)；W為垂直方向補(bǔ)給強(qiáng)度。

求解上述地下水?dāng)?shù)學(xué)模型偏微分方程需要初始條件和邊界條件等一系列條件，初始條件為零時(shí)刻的邊界條件。

(1) 第一類(lèi)邊界(已知水頭邊界)

(2)

其中，H(x，y，z，t)表示三維條件下邊界段Si上點(diǎn)(x，y，z)在t時(shí)刻的水頭，是Si上的已知函數(shù)。

(2) 第二類(lèi)邊界(已知流量邊界)

(3)

其中，n為邊界Si的外法線方向。qi為已知函數(shù)，表示Si上單位面積的側(cè)向補(bǔ)給量。

(3) 第三類(lèi)邊界(混合邊界)

(4)

其中，α、β為已知函數(shù)；H為潛水含水層水頭，m。

2.3 源匯項(xiàng)

(1) 補(bǔ)給項(xiàng)

盆地內(nèi)第四系地下水主要接受大氣降雨補(bǔ)給及盆地周邊山區(qū)基巖裂隙巖溶水側(cè)向補(bǔ)給。據(jù)統(tǒng)計(jì)，當(dāng)?shù)囟嗄昶骄涤炅?15.8 mm，降雨入滲系數(shù)為0.03～0.25?；诒敬我巴夤ぷ髌陂g屬旱季-平季水位，故本次選擇降雨入滲補(bǔ)給量稍小，調(diào)整為0.002 6 m/d。

側(cè)向徑流補(bǔ)給量采用達(dá)西公式(式5)進(jìn)行計(jì)算。

Q=KIHL

(5)

式中,Q為側(cè)向徑流補(bǔ)給量，m3；K為剖面含水層的平均滲透系數(shù)，m/d；I為地下水的平均水力坡度，無(wú)因次；H為含水層厚度，m；L為垂直于地下水流方向的斷面長(zhǎng)度，m。

依據(jù)公式(式5)，計(jì)算可得盆地側(cè)向凈補(bǔ)給總量為931.7 m3/d。

(2) 排泄項(xiàng)

盆地內(nèi)第四系地下水，主要以地下徑流集中向盆地北部大屯海及長(zhǎng)橋海排泄?；谂璧氐叵滤饕獮槌袎核?，其蒸發(fā)強(qiáng)度較弱，故本次模型不考慮地下水的蒸發(fā)排泄方式。

2.4 參數(shù)分區(qū)與取值

盆地內(nèi)各區(qū)塊的滲透性能表現(xiàn)出一定的差異性，為了提高數(shù)值模擬的準(zhǔn)確性，就需要對(duì)各區(qū)塊賦以不同的滲透系數(shù)。本次研究主要通過(guò)抽水試驗(yàn)來(lái)獲取研究區(qū)的滲透系數(shù)，對(duì)未進(jìn)行抽水試驗(yàn)的民井等根據(jù)含水層的性質(zhì)給以經(jīng)驗(yàn)值。最終將模擬區(qū)砂礫石層含水層分為13個(gè)不同滲透系數(shù)的小區(qū)域(圖2)，滲透系數(shù)分區(qū)賦值見(jiàn)表1。

表1 滲透系數(shù)分區(qū)表

圖2 滲透系數(shù)分區(qū)圖

2.5 流場(chǎng)模型的建立與驗(yàn)證

(1) 模型建立

本次模擬主要以第四系中的砂礫石層含水層為主要研究對(duì)象。利用GMS軟件將整個(gè)盆地第四系概化成一個(gè)18 177 m×17 851 m范圍的模擬區(qū)。首先將其剖分成100×100×3的單元格，其次將地層分層數(shù)據(jù)及野外地質(zhì)調(diào)查點(diǎn)(民井)、水文地質(zhì)鉆探點(diǎn)數(shù)據(jù)導(dǎo)入到軟件中，運(yùn)用反距離差值的方法建立研究區(qū)三維有限差分網(wǎng)格模型，最后將模型邊界條件、水文地質(zhì)參數(shù)導(dǎo)入至三維有限差分網(wǎng)格模型中，從而建立研究區(qū)的地下水三維穩(wěn)定流流場(chǎng)數(shù)值模型。

(2) 模型驗(yàn)證

為了驗(yàn)證建立的模型與研究區(qū)實(shí)際地下水流場(chǎng)的吻合性，需要對(duì)研究區(qū)的的實(shí)際觀測(cè)水位與模擬的計(jì)算水位做擬合驗(yàn)證[5]。其驗(yàn)證原理是在模型內(nèi)的已有觀測(cè)水位點(diǎn)處設(shè)置誤差檢測(cè)棒，同時(shí)設(shè)置誤差控制值，然后將觀測(cè)點(diǎn)的實(shí)際觀測(cè)水位導(dǎo)入模型中，系統(tǒng)就會(huì)根據(jù)模擬水位與實(shí)際水位進(jìn)行擬合比較[6]。若觀測(cè)水位與計(jì)算水位之差在誤差控制值之內(nèi)，誤差棒顯示為綠色；若超出誤差控制值，但小于控制值的200%，誤差棒顯示為黃色；若超出誤差控制值的200%以上，誤差棒顯示為紅色[7]，如圖3。

圖3 誤差棒示意圖

研究區(qū)內(nèi)共布設(shè)有20個(gè)觀測(cè)孔作為本次模擬的驗(yàn)證依據(jù)。本次模擬研究擬合驗(yàn)證誤差控制值選擇2 m，經(jīng)擬合后發(fā)現(xiàn)，20口觀測(cè)民井中有17口民井水位與模型能夠擬合，有3口民井水位超過(guò)2 m的誤差控制值，但控制在200%的控制值(4 m)內(nèi)，擬合結(jié)果見(jiàn)圖4。

由地下水流場(chǎng)模擬圖可知，本次建立的研究區(qū)第四系地下水流場(chǎng)模型與實(shí)際地下水位擬合較好，基本能夠反映盆地第四系地下水的流場(chǎng)情況。從地下水流場(chǎng)模型圖可以看出，盆地內(nèi)第四系地下水呈南高北低趨勢(shì)，地下水整體由東南向西北方向徑流，最終排泄于長(zhǎng)橋海和大屯海。

為了更精確的反映廠區(qū)及其周邊第四系松散層的地下水流場(chǎng)特征，本次研究將依據(jù)整個(gè)盆地的地下水流場(chǎng)特征再進(jìn)一步模擬廠區(qū)及周邊范圍第四系地下水流場(chǎng)特征。

圖4 地下水流場(chǎng)模擬圖

3 廠區(qū)地下水流場(chǎng)特征

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

第四系地下水流場(chǎng)的建立，為廠區(qū)區(qū)域水文地質(zhì)單元邊界的劃定提供了重要依據(jù)。本次研究將廠區(qū)區(qū)域單獨(dú)劃分成一個(gè)由東南部邊界補(bǔ)給、西北部邊界排泄的定水頭水文地質(zhì)單元。其中，將垂直于等水位線的東西兩側(cè)邊界視為等效隔水邊界。邊界條件概化后的概念模型如圖5。

圖5 廠區(qū)第四系水文地質(zhì)單元概念模型圖

3.2 流場(chǎng)模型的建立與驗(yàn)證

廠區(qū)區(qū)域地下水流場(chǎng)模型建立的過(guò)程與第四系地下水流場(chǎng)模型相似?；诒敬文M區(qū)范圍較小，考慮到精度問(wèn)題，本次模擬研究擬合驗(yàn)證誤差控制值選擇1.5 m，經(jīng)擬合后發(fā)現(xiàn)，本次共設(shè)置的12個(gè)監(jiān)測(cè)點(diǎn)水位中有11個(gè)監(jiān)測(cè)點(diǎn)水位能夠與模型計(jì)算水位相擬合(顯示為綠色)，僅有廠區(qū)內(nèi)1個(gè)孔的監(jiān)測(cè)水位1.5 m的誤差控制值，但只低于計(jì)算水位的1.558 m，在200%的控制值(3 m)內(nèi)(顯示為黃色)，擬合結(jié)果見(jiàn)圖6。

總體上，本次模擬計(jì)算水位與實(shí)際觀測(cè)水位整體擬合較好，基本能夠反映廠區(qū)區(qū)域的地下水流場(chǎng)特征。

圖6 廠區(qū)區(qū)域地下水流場(chǎng)模擬圖

4 地下水污染規(guī)律研究

該廠區(qū)污染方式主要是廢渣污染物伴隨降雨入滲擴(kuò)散。污染物主要通過(guò)上層粉質(zhì)粘土遷移至中層粉細(xì)砂層承壓含水層，最終在承壓含水層中順?biāo)飨虮睌U(kuò)散。本次研究利用GMS軟件的MODFLOW和MT3D模塊，結(jié)合廠區(qū)流場(chǎng)模型建立廠區(qū)區(qū)域溶質(zhì)運(yùn)移模型并對(duì)模型進(jìn)行驗(yàn)證。以驗(yàn)證過(guò)的模型為基礎(chǔ)，深入分析污染物在降雨入滲擴(kuò)散方式下在不同介質(zhì)和不同源強(qiáng)條件下的運(yùn)移規(guī)律。

4.1 高濃度下污染物的遷移規(guī)律

本次選擇在氨氮泄露工況下，模擬研究污染物在降雨入滲方式下在粉質(zhì)粘土層中的遷移規(guī)律。將氨氮污染源強(qiáng)設(shè)置為250 000 mg/l，將其注入ZK3處，任其自由擴(kuò)散，模擬結(jié)果顯示，在這種超高濃度的污染物泄露工況下，污染物僅5 d就能滲透穿過(guò)粉質(zhì)粘土層進(jìn)入第二層含水層，見(jiàn)圖7。據(jù)鉆孔資料顯示，ZK3處粉質(zhì)粘土層層厚約7.8 m，故在這種超高濃度工況下，污染物在粉質(zhì)粘土含水層垂向上的平均遷移速率為1.56 m/d。

圖7 高濃度氨氮污染物擴(kuò)散圖

在第一層粉質(zhì)粘土含水層中，污染物在水平面上擴(kuò)散緩慢，其擴(kuò)散過(guò)程如圖8。由模擬結(jié)果可知，在這種超高濃度污染泄露工況下，污染物從100 d～1 000 d擴(kuò)散了4.5 m，1 000～2 000 d擴(kuò)散了3.2 m，2 000～3 000 d擴(kuò)散了2.5 m，3 000～4 000 d擴(kuò)散了1.4 m，4 000～5 000 d擴(kuò)散了1 m，平均擴(kuò)散速率0.002 5 m/d。整體看來(lái)，隨著時(shí)間推移和擴(kuò)散范圍的擴(kuò)大，污染物在水平面上的擴(kuò)散速率越來(lái)越慢。

圖8 粉質(zhì)粘土層高濃度氨氮污染物水平擴(kuò)散圖

污染物穿過(guò)第一層粉質(zhì)粘土層后，污染物運(yùn)移受第二層承壓含水層水力作用，運(yùn)移速度加快，主要順?biāo)餮厮μ荻容^低的西北方向呈長(zhǎng)條狀運(yùn)移，如圖9。

圖9 粉細(xì)砂層高濃度氨氮污染物水平擴(kuò)散圖

由模擬結(jié)果可知，污染物從第2 d～500 d運(yùn)移了46 m，500～1 000 d運(yùn)移了37 m，1 000 d～1 500 d運(yùn)移了32 m，從1 500～2 000 d運(yùn)移了35 m，從2 000～2 500 d運(yùn)移了38 m。其擴(kuò)散統(tǒng)計(jì)結(jié)果見(jiàn)圖10。由擴(kuò)散統(tǒng)計(jì)圖可知，該污染物在第二層含水層中大致呈勻速運(yùn)移，平均遷移速率0.075 m/d。

圖10 粉細(xì)砂層高濃度氨氮污染物擴(kuò)散統(tǒng)計(jì)圖

4.2 低濃度下污染物的遷移規(guī)律

本次選擇廠區(qū)南側(cè)廢渣中的錳離子污染物在降雨入滲方式下，分析研究錳離子在粉質(zhì)粘土中的遷移規(guī)律。在廢渣污染物中，本次擬定的錳離子的污染源強(qiáng)為3 mg/l。模型顯示，在這種低濃度工況下，錳離子通過(guò)降雨入滲作用穿過(guò)粘土層大約需要600 d。其遷移情況見(jiàn)圖11。鉆孔資料顯示，廠區(qū)南部的粉質(zhì)粘土層厚約8.5 m左右，則在低濃度下，錳離子在粉質(zhì)粘土垂向上的平均遷移速率為0.014 m/d。

圖11 低濃度錳污染物擴(kuò)散圖

在粉質(zhì)粘土介質(zhì)中，錳污染物在水平面上的擴(kuò)散速度非常慢，經(jīng)過(guò)5 000 d的時(shí)間，污染物僅向北邊擴(kuò)散了接近一個(gè)格子，約8 m左右，平均水平運(yùn)移速率約為0.001 6 m/d，污染物擴(kuò)散情況如圖12。

綜上，在降雨入滲方式下，低濃度的錳污染物在粉質(zhì)粘土介質(zhì)中主要以垂向運(yùn)移為主，垂向運(yùn)移速率約0.014 m/d，水平運(yùn)移速率約0.001 6 m/d，垂向運(yùn)移速率約為水平運(yùn)移速率的10倍。

當(dāng)污染物穿過(guò)第一層粉質(zhì)粘土層進(jìn)入第二層粉細(xì)砂含水層后，污染物運(yùn)移受承壓含水層水力作用，運(yùn)移速度加快，主要順?biāo)餮厮μ荻容^低的西北方向呈長(zhǎng)條狀運(yùn)移，如圖13。

由模型可知，在第二層粉細(xì)砂含水層中，污染物從1 000 d～2 000 d運(yùn)移了43 m，從2 000 d～30 00 d運(yùn)移了49 m，從3 000 d～4 000 d運(yùn)移了60 m，從4 000 d～5 000 d運(yùn)移了63 m，其擴(kuò)散統(tǒng)計(jì)結(jié)果見(jiàn)圖14。由統(tǒng)計(jì)結(jié)果折線圖可知，該污染物在第二層含水層中大概呈勻速擴(kuò)散趨勢(shì)，平均擴(kuò)散速率0.072 m/d。

綜上，本次模擬研究了在降雨入滲擴(kuò)散方式下，污染物以不同污染源強(qiáng)在粉質(zhì)粘土和粉細(xì)砂兩種介質(zhì)地下水中的遷移規(guī)律，統(tǒng)計(jì)結(jié)果見(jiàn)表2。

表2 污染物擴(kuò)散速率統(tǒng)計(jì)表

由統(tǒng)計(jì)結(jié)果可知，降雨入滲擴(kuò)散方式下污染物在地下水中的遷移速率與污染物的濃度及地下水介質(zhì)密切相關(guān)，總結(jié)規(guī)律如下：

(1) 在粉質(zhì)粘土潛水含水層中污染物主要以垂向運(yùn)移為主，并且受污染源強(qiáng)影響較大，污染物濃度越高，其垂向運(yùn)移速率越快。

(2) 在垂向上，污染物擴(kuò)散速率還與污染物載體有關(guān)。

(3) 在粉質(zhì)粘土潛水含水層中，污染物在水平面上的擴(kuò)散速率受污染源強(qiáng)的影響較小。

(4) 在污染物穿過(guò)上層粉質(zhì)粘土層進(jìn)入粉細(xì)砂承壓含水層后，污染物主要順?biāo)飨蛩μ荻容^低的方向運(yùn)移，其運(yùn)移速率受濃度影響較小，主要以機(jī)械彌散為主。

圖12 粉質(zhì)粘土層低濃度錳污染物水平擴(kuò)散圖

圖13 粉細(xì)砂層低濃度錳污染物水平擴(kuò)散圖

圖14 粉細(xì)砂層低濃度錳污染物水平擴(kuò)散圖

5 結(jié)語(yǔ)

(1) 將該盆地第四系地層作為一個(gè)獨(dú)立的水文地質(zhì)單元，第四系含水層地下水主要接受大氣降雨補(bǔ)給以及東部、南部及西部的山區(qū)基巖碳酸鹽巖裂隙巖溶水的側(cè)向補(bǔ)給，經(jīng)地下徑流最終向北部長(zhǎng)橋海和大屯海排泄。

(2) 在盆地區(qū)域地下水流場(chǎng)的基礎(chǔ)上圈定了盆地內(nèi)廠區(qū)地塊地下水流場(chǎng)邊界，建立了污染廠區(qū)區(qū)塊地下水流場(chǎng)，整體擬合結(jié)果較好。

(3) 通過(guò)驗(yàn)證的溶質(zhì)運(yùn)移模型分析研究了污染物在歷史污染過(guò)程中不同工況下在地下水中的遷移規(guī)律。降雨入滲擴(kuò)散方式下，污染物在地下水中的遷移速率與污染物的濃度及地下水介質(zhì)密切相關(guān)。在潛水含水層中，污染物主要以垂向運(yùn)移為主，并且受污染源強(qiáng)影響較大，污染物濃度越高，其垂向運(yùn)移速率越快，污染物在水平面上的擴(kuò)散速率受污染源強(qiáng)的影響較??；在承壓含水層中，污染物沿水流方向遷移，并且呈現(xiàn)勻速擴(kuò)散趨勢(shì)。