基于GMS模擬預測云南保山盆地孔隙水中氨氮的遷移規(guī)律

楊加明,徐世光,2,李勝,張佳慧

(1.昆明理工大學國土資源工程學院,昆明 650093;2.云南地礦工程勘察集團公司,昆明 650011;3.云南南方地勘工程有限公司,大理 671000)

1 引言

水資源是人類生存生活及生產(chǎn)所必備的條件,尤其是地下水資源對人類的生存和發(fā)展具有助力作用,地下水是我們?nèi)祟惣安糠治锓N賴以生存和整個社會發(fā)展的前提基礎(chǔ),是經(jīng)濟社會可持續(xù)發(fā)展和生態(tài)系統(tǒng)循環(huán)的最基礎(chǔ)的物質(zhì)[1-2]。近些年地下水污染問題顯得日漸突出,為了促進經(jīng)濟的發(fā)展,工業(yè)的飛速發(fā)展,造成地下水污染現(xiàn)象明顯[3-4]。地下水開發(fā)利用過程中存在正常水位出現(xiàn)明顯下降、地面沉降、水質(zhì)狀況惡化等重要問題[5]。王攀[6]等人對永城市淺層地下水中污染物的來源作了解析,導致該區(qū)域中硝酸鹽和硫酸鹽超標的主要來源是工業(yè)污水排放,農(nóng)業(yè)農(nóng)藥和化肥的過度使用;李政紅[7]等人探究人類活動對內(nèi)蒙古托克縣淺層地下水中的硝酸鹽的驅(qū)動作用,結(jié)果顯示人類污水灌溉、糞便堆放、化肥施用對該區(qū)域地下水中的硝酸鹽具有很高的貢獻率;郇環(huán)[8]等人利用普通克里金插值方法對松花江吉林段沿岸淺層地下水中硝酸鹽的氮濃度分布特征作了分析,結(jié)果顯示常數(shù)趨勢球型模型適用于描述研究區(qū)地下水NO3-N濃度分布空間變異性,NO3-N濃度在一定范圍內(nèi)存在空間相關(guān)性,空間相關(guān)距離為相對較大,前人對淺層地下水中的污染因子研究較多的是其污染來源、空間分布變化及水質(zhì)評價。本篇論文是采用地下水數(shù)值模擬軟件對研究區(qū)中的氨氮在現(xiàn)狀污染的基礎(chǔ)上,預測其未來的污染變化趨勢,對保護地下水資源具有一定的指導意義,數(shù)值模擬方法靈活,適用于各種復雜的水文地質(zhì)條件,在處理地下水流、地下水溶質(zhì)運移和水資源環(huán)境分析等方面有良好的效果[9-11]。

本文以云南保山盆地地下水污染為例,借助水文地質(zhì)相關(guān)資料,結(jié)合研究區(qū)的實際核心特征,利用GMS專業(yè)軟件對污染物進行仿真模擬,建立研究區(qū)內(nèi)流場的數(shù)學研究模型,選用氨氮作為特征因子,預測特征因子的運移規(guī)律對地下水造成的環(huán)境影響。數(shù)值模擬對地下水開采、治理、環(huán)境問題和資源的評價開發(fā)等提供參考意義[12-13]。

經(jīng)研究,該區(qū)域氨氮污染的主要來源是保山地區(qū)動物的禽畜糞便的排放,經(jīng)相關(guān)資料統(tǒng)計出該區(qū)域的動物的數(shù)量,折合成標準豬,每年每頭豬的糞尿排放量為1.2 t,糞尿中的2%排入到土壤中,經(jīng)降雨入滲排的氨氮含量每年是25.01 t;農(nóng)業(yè)過度施肥也是造成地下水中氨氮超標的主要來源,保山的第一產(chǎn)業(yè)是農(nóng)業(yè),大面積農(nóng)作物需要施肥,導致土壤中含有大量的含氮化肥,經(jīng)計算,保山每年經(jīng)降雨淋濾進入地下水中的總氮為212.4 t;工業(yè)污水的排放與泄露,經(jīng)調(diào)查研究,研究區(qū)內(nèi)共有92座工業(yè)企業(yè)的廢水直接排入到東河內(nèi),每年排放污水的總氮含量是29.28 t。若不及時治理,地下水中的氨氮會在特定環(huán)境下轉(zhuǎn)化為亞硝酸鹽,亞硝酸鹽含量的升高會增加人體患癌癥的幾率;當然地下水中的氨氮超標對生態(tài)環(huán)境也有極大的危害,會造成魚類的死亡,破壞生態(tài)平衡。

2 研究區(qū)概況與污染現(xiàn)狀

2.1 概況

以云南保山盆地為研究區(qū),其區(qū)域內(nèi)水文地質(zhì)條件稍顯復雜,地下水類型劃分較多,盆地內(nèi)的第四系松散層孔隙水與地貌類型密不可分,因此該區(qū)內(nèi)的孔隙水在外在形式上表現(xiàn)出多種多樣,沖洪積扇孔隙水具有較強的富水性,主要以鵝卵石為主,鉆孔單位涌水量很大;沖湖積平原孔隙水因地勢較為平坦,地下水位埋深較淺,等水位線走勢常年保持基本不變;沖洪積平原孔隙水,富水性中等,盆地內(nèi)的地下水具有各項異性的特征。本次重點研究的工作是第四系松散層內(nèi)的孔隙水的污染核心特征,該地下水類型主要接受大氣降雨的自然補給,還有部分接受周圍具有水力聯(lián)系的基巖裂隙水的徑流補給,最終均以泉的形式以及下游向侵蝕基準面較低的東河排泄。研究區(qū)內(nèi)的地下水取樣點見圖1。

圖1 研究區(qū)地下水取樣點圖

2.2 污染現(xiàn)狀

根據(jù)2019年采取盆地內(nèi)35個淺層地下水樣品的檢測結(jié)果顯示,氨氮的濃度在0.02～17.9 mg/l,取樣點中漢莊鎮(zhèn)914處、盆地東側(cè)青陽村916處、板橋鎮(zhèn)920處及河圖鎮(zhèn)924處的氨氮濃度高于地下水Ⅲ類標準(<0.5 mg/ l),且最高濃度在916處,研究區(qū)內(nèi)的氨氮主要集中在東河區(qū)域附近,越靠近東河,污染越嚴重,由此可見工業(yè)廢水的排放對東河附近的地下水污染是最嚴重的,若不及時處理,東河附近的居民較多,不知情的情況下,飲用氨氮超標的地表水和地下水,均會影響居民的身體健康。氨氮的污染現(xiàn)狀圖見圖2。由氨氮的污染現(xiàn)狀圖顯示,盆地內(nèi)孔隙水中的氨氮污染區(qū)域較大,主要集中在盆地的上游,下游盆地內(nèi)地下水中的氨氮并未污染。

圖2 氨氮污染現(xiàn)狀圖

3 水文地質(zhì)特征及概念模型

3.1 研究區(qū)區(qū)域地質(zhì)條件

水文地質(zhì)概念模型可以很好地反應(yīng)研究區(qū)內(nèi)地下水的流向,可以將復雜的問題簡單化,在計算機上清晰地反應(yīng)出來。

研究區(qū)的水文地質(zhì)條件較為復雜,盆地內(nèi)主要是第四系覆蓋層,該層主要是富水性較強的孔隙水含水層組,該層地下水因地貌類型不同,主要有第四系沖洪積扇孔隙水,沖湖積平原孔隙水,沖洪積平原孔隙水,孔隙水中的隔水層主要是上部粘土層,而含水層主要是下部的帶有卵礫石的泥砂層。主要分布在盆地內(nèi)及盆地邊緣與巖層的交界處。

裂隙水含水層主要分布在盆地西部和東部的基巖中,因基巖區(qū)中的裂隙發(fā)育多,地下水賦存于此,主要的含水層是粉砂巖,而主要的隔水層則是板巖。

巖溶水含水層主要是研究區(qū)內(nèi)碳酸鹽巖夾碎屑巖裂隙溶洞水,主要分布在東側(cè)及盆地下游側(cè),研究區(qū)內(nèi)的水文地質(zhì)圖見圖3。從研究區(qū)的剖面圖可以看出地下水含水層的補給關(guān)系,盆地內(nèi)的第四系孔隙水,主要接受碎屑巖裂隙水及碳酸鹽巖巖溶水的補給,見圖4。

圖3 研究區(qū)水文地質(zhì)圖

圖4 研究區(qū)水文地質(zhì)剖面圖

3.2 研究區(qū)概念模型

因本篇論文主要研究的第四系孔隙水,研究區(qū)內(nèi)地下水流場模型的范圍是第四系含水層,地層中的黏土層是主要的隔水層,而主要的含水層是卵礫石層。盆地中的卵礫石層與碳酸鹽巖相交的部分,接受基巖裂隙水的補給,因而可以看作是向研究區(qū)內(nèi)松散孔隙水的定流量邊界的界限,而基巖覆蓋層中的黏土與紅黏土,因透水性差,故而可以看作是隔水邊界。而主要的含水層的補給途徑主要是大氣降水以及基巖的橫向補給,最終向低洼地區(qū)的東河流域排泄,該流域的整體流向貫穿整個盆地的南北向,由地下水與地表水的水力聯(lián)系可看出,該河流可作為定水頭邊界。依據(jù)上述假設(shè)的邊界條件的界限以及該區(qū)域內(nèi)的水文地質(zhì)特征,建立研究區(qū)的水文地質(zhì)概念模型見圖5所示。

圖5 研究區(qū)水文地質(zhì)概念模型圖

3.3 數(shù)學模型

對于新建的地下水流場模型使用偏微分方程數(shù)值化。因研究區(qū)內(nèi)的地下水是承壓含水層,且具有非均質(zhì)各項異性的特性,故使用以下的微分方程來求解:

(1)

其中,Kxx,Kyy,Kzz為各主方向上的滲透系數(shù),m/d;H為承壓含水層水頭,m;Ss為儲水率,單位m-1;W為源匯項,m3;t為時間,單位d。

求解該方程需要的初始條件為零時刻的邊界條件。

H(x,y,z,t)Si=φi(x,y,z,t)(x,y,z)∈Si

(2)

因地下水流場中的偏微分方程具有較多的未知數(shù),也缺少很多條件,因此很難解得具體的解析解,為了解決這一難題,引進數(shù)值模擬的辦法來求得上述微分方程的近似解。

數(shù)值模擬過程主要是利用有限差分的方式將模擬的區(qū)域劃分為若干的小網(wǎng)格,研究每個小網(wǎng)格中的流場特征,最終以替代的思想將研究區(qū)的整體流場特征用每個小網(wǎng)格的模擬水位來代替,進而實現(xiàn)離散的目的,將不能求解的復雜微分方程轉(zhuǎn)化為可求解的簡單數(shù)學問題,忽略微小的誤差。

3.4 數(shù)值模型建立

3.4.1 建立模型

根據(jù)該區(qū)的巖性特征分析,以卵礫石含水層為主要的研究對象,含水層上下部的黏土可概括為隔水層,也就是對應(yīng)的儲水界面,按模型面積劃分成200 200×400×3的網(wǎng)狀結(jié)構(gòu)。利用以上假定的邊界條件,將研究區(qū)內(nèi)的水文地質(zhì)概念模型轉(zhuǎn)化為數(shù)值研究模型,參見附圖6.

圖6 研究區(qū)數(shù)學模型三維地質(zhì)網(wǎng)格圖

3.4.2 源匯項

(1) 補給項

該區(qū)的卵礫石含水層是由大氣降水和基巖裂隙水橫向供給的。根據(jù)相關(guān)的水文氣象數(shù)據(jù),得出了以卵礫巖為基礎(chǔ)的降雨入滲系數(shù)為0.24的降水量,并據(jù)此推算出了降水對地下水的補充量為0.000 692 3 m/d。

而四周基巖裂隙水的補給量可根據(jù)達西公式計算:

Q=KAJ

(3)

式中,Q為周圍基巖區(qū)側(cè)向徑流孔隙水的補給量,m3;K為滲透系數(shù),m/d;J為水力梯度;A為過水斷面的面積,m2。

由研究區(qū)內(nèi)巖性出露的條件及地下水位的關(guān)系看出,與該區(qū)內(nèi)的孔隙水具有水力聯(lián)系的基巖裂隙水主要是西部胡家壩、莊家坡、龍王廟,北部沙壩及東部青陽村、大關(guān)廟等地。據(jù)此計算出了側(cè)向補給總量為2 096.83 m3/d。

(2) 排泄項

因主要的含水層是承壓水且上覆土層具有隔水性質(zhì),故模型中不考慮蒸發(fā)排泄。研究區(qū)的孔隙水主要通過周圍地勢較低的河流、泉及少量的民井排泄。

(3) 水文地質(zhì)參數(shù)

滲透系數(shù)是反應(yīng)巖土滲透性的重要指標,也是建立流場模型及溶質(zhì)運移模型必備的參數(shù),因研究區(qū)內(nèi)出露的巖性較多較復雜,因此按照鉆孔抽水試驗、地層巖性及地形地貌等因素劃分滲透系數(shù)分區(qū),盆地內(nèi)的第四系覆蓋層厚度變化較大,展現(xiàn)出中間厚、邊緣薄的特點,且沿著沖積扇,卵礫石層與粘土層交替沉積,且厚度不一,因此盆地內(nèi)不同區(qū)域的滲透系數(shù)存在差異,呈現(xiàn)為非均質(zhì)性。本次將研究區(qū)孔隙水含水層劃分為27個子區(qū)域,在建立模型時,對不同的區(qū)域進行賦值,以達到模型的準確性,其參數(shù)分區(qū)平面圖見圖7,分區(qū)表見表1。

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

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

根據(jù)保山的地質(zhì)勘察資料及實際工程經(jīng)驗,研究區(qū)的平均孔隙度設(shè)定為0.28。

(4) 彌散度

根據(jù)研究區(qū)做的彌散試驗,結(jié)合研究區(qū)的研究面積,本次模擬的孔隙水含水層縱向彌散度為0.82 m,橫向彌散度為0.19 m。

4 地下水溶質(zhì)運移模型

4.1 模型建立

對于利用GMS建立污染物的溶質(zhì)運移模型,諸多學者也有許多研究,劉娟等[14]對某市巖溶地下水污染進行研究,運用FEFLOW 軟件對擬建好的數(shù)學研究模型求解,探究研究區(qū)內(nèi)流場的特征,在水流模型的基礎(chǔ)上進行溶質(zhì)運移模擬,模擬了地下水中四氯化碳污染羽的濃度場,并預測了未來濃度場的變化;于翠翠[15]運用地下水模擬軟件 GMS 對山東濟南明水泉域地下水流場進行了數(shù)值模擬,對泉域巖溶地下水資源總量和地下水可采量進行了評價,并對泉水動態(tài)水位進行預測;顏萍[16]利用GMS對保山某尾礦庫的地下水污染物進行污染預測,針對得到的預測結(jié)果,提出治理污染的措施。

該地區(qū)內(nèi)的氨氮的污染來源主要是禽畜糞便、農(nóng)業(yè)過度使用含氮化肥以及工業(yè)污水的隨意排放造成的,上述污染源通過上部地層的孔隙經(jīng)降雨淋濾進入到研究區(qū)內(nèi)的含水層中,最終沿地下水流向方向遷移。本次工作研究的重點是僅考慮地下水的對流及彌散作用運用GMS專業(yè)模擬軟件中的MT3D模塊結(jié)合上述建立的流場數(shù)值模型建立研究區(qū)內(nèi)的溶質(zhì)運移模型,探究預測特征因子在1 000 d、3 000 d、5 000 d后氨氮的地理位置以及中心濃度的運移方向。根據(jù)水動力彌散理論,研究區(qū)內(nèi)地下水溶質(zhì)運移控制方程如下:

(4)

式中,C為地下水中污染物濃度,mg/l;n為孔隙度;Dij為水動力彌散系數(shù),m;vi為地下水滲流速度,m/d;qs為源匯到含水層內(nèi)單位體積的流量;Cs為源匯項中污染物濃度;Rn為化學反應(yīng)項。

初始條件和邊界條件為:

C(x,y,z,t)|t=0=C0(x,y,z)∈Ω

(5)

C(x,y,z,t)|Γ1=CI(x,y,z) (x,y,z)∈Γ1∈Ω

(6)

式中,C0為地下水污染物初始濃度,mg/l;Ω為模擬區(qū)域;CI為邊界污染物濃度,mg/l;Γ1為濃度邊界。

因條件限制,該模型中僅考慮了地下水的對流作用及彌散作用,故模型不考慮污染物的質(zhì)量衰減的因素。利用GMS軟件模擬預測研究區(qū)內(nèi)氨氮的運移規(guī)律。

4.2 模擬預測及結(jié)論

根據(jù)研究區(qū)采集到的地下水中氨氮濃度的檢測結(jié)果,利用克里金插值方進行插值計算,可知氨氮在研究區(qū)內(nèi)的現(xiàn)狀污染情況,以及現(xiàn)狀條件下氨氮的實際遷移距離見圖8。結(jié)合研究區(qū)內(nèi)氨氮的污染原因,依據(jù)數(shù)據(jù)資料分析可計算出每天地下水中有0.001 35 mg/l的氨氮進入,這里將氨氮的單天滲入量作為模擬的初始源強,模擬預測氨氮在1 000 d、3 000 d、5 000 d后高濃度地區(qū)的位置、濃度變化以及遷移距離,參見附圖9(a、b、c)。

圖8 氨氮污染現(xiàn)狀及實際遷移距離圖

圖9(b) 氨氮遷移 3 000 d濃度變化及理論距離圖

圖9(c) 氨氮遷移5 000 d濃度變化及理論距離圖

由上述模擬預測圖表明,僅考慮對流及彌散作用時,模擬的3個時間段內(nèi)氨氮的理論遷移距離與實際遷移距離相比有明顯的變化。以上述計算出單日氨氮進入地下水的含量作為源強,對比氨氮現(xiàn)狀污染的遷移距離,1 000 d后,東側(cè)部分的的氨氮向盆地中心處的東河運移了14.9 m,而西側(cè)向盆地中心處的東河運移25.7 m;3 000 d后,東側(cè)部分的氨氮向東河運移了49.7 m,西側(cè)向東河運移80.5 m;5 000 d后,東側(cè)部分的氨氮向東河運移了85.7 m,西側(cè)向東河運移了136.3 m。由此可以看出模擬預測的3個時間段內(nèi),氨氮兩側(cè)的濃度均向盆地中心的東河移動,并且理論遷移的距離與時間呈現(xiàn)正相關(guān)性。由前述可知,東河是盆地內(nèi)的低洼地區(qū),研究區(qū)內(nèi)地下水最終排泄到東河內(nèi),根據(jù)模擬結(jié)果也可以看出(表2),氨氮都是沿著地下水流方向遷移的,并且時間越長,最高濃度在緩慢變小。利用氨氮實際運移距離與理論運移距離的比值可以看出,西側(cè)的比值比東側(cè)的大,說明西側(cè)的氨氮更易沿水流方向向東河遷移。

表2 氨氮污染實際遷移距離與理論遷移距離統(tǒng)計表

根據(jù)上述模擬得到氨氮運移的理論距離,繪制研究時間與污染物遷移距離的關(guān)系圖見圖10,根據(jù)關(guān)系圖來分析東西向污染物遷移的變化,為了時間與探究出二者的關(guān)系,將得到的遷移距離與時間進行散點擬合分析得到式(7)、式(8)。

圖10 氨氮不同時間與理論運移距離關(guān)系圖

S東=1×10-7t2+0.016 7t-1.69

(7)

S西=2×10-7t2+0.026 5t-0.78

(8)

由上述得到擬合的關(guān)系式看出,氨氮的理論遷移距離與時間并不是簡單的正相關(guān)性,而是呈現(xiàn)開口向上的二次函數(shù)形式,且不同東西方向的相關(guān)系數(shù)R2分別為0.994和0.992,說明擬合效果特別好,但二次項的系數(shù)很小,表明氨氮雖然沿著地下水流向運移,但運移速度很慢。

5 結(jié)論

(1) 由氨氮在不同方向運移的距離可知,西側(cè)的遷移距離明顯要比東側(cè)遷移的距離大,說明西側(cè)的氨氮更易沿水流方向遷移,并且西側(cè)的水力梯度明顯高于東側(cè)的。

(2) 研究區(qū)內(nèi)的氨氮污染物是沿著地下水流方向向東河遷移的,說明氨氮的運移主要是由對流作用控制的,隨著時間的變化,西側(cè)遷移的理論距離明顯高于東側(cè)的,對流作用控制氨氮運移的影響越大。

(3) 氨氮的遷移距離與時間呈現(xiàn)二次函數(shù)關(guān)系,但二次項的系數(shù)很小,說明氨氮在地下水水中是緩慢移動的。

(4) 根據(jù)模擬得到的結(jié)果,時間越長,東河流域的水質(zhì)越易受到氨氮的污染,因此要提出合理的治理措施,否則會對居住在東河流域的居民生活有很大的影響。