河岸潛流帶中非穩定地下水流的溶質運移模擬研究

郭 昆，李 寧，熊啟華，王芮瓊，廖春來

(1.湖北省地質環境總站，湖北 武漢 430034; 2.資源與生態環境地質湖北省重點實驗室，湖北 武漢 430034；3.中國地質大學(武漢)環境學院，湖北 武漢 430078)

在工農業密集的沖洪積平原，賦存著豐富的孔隙水，且地表水系發育，通常地表水和地下水具有一定的水力聯系和相互作用[1]，這一相互作用的關鍵帶稱為潛流帶，它對局部地區的生態環境和水資源起著重要作用。潛流帶的地下水往往隨著地表水水位的波動存在著與地表水互補的現象，形成非穩定流，因此對于潛流帶中的溶質遷移研究一直是個熱點，也是難點。杜堯等[2]系統梳理了潛流帶的水文和生物地球化學過程的相關原理和研究方法，此外有學者通過實驗尺度的數值模擬和室內實驗等方法，探究了溶質在潛流帶的遷移規律[3-5]。但是這種實驗尺度的研究較為理想化，對于實際復雜的水文地質條件的適用性還有待檢驗，另外以往的研究對象多為河床潛流帶，而對范圍較大的河岸潛流帶的研究較為少見。

地下水溶質運移模型是預測溶質在地下水中運移的有力手段，主要利用數學公式和一些基本條件刻畫溶質在地下水中的運動、降解等物理化學過程。模型的主體部分包括控制方程、初始條件和邊界條件，它們決定了溶質在研究區的初始狀態及其與外界的轉換方式。目前溶質運移模型已成功應用到地下水環境影響評價和水污染防治等工作和研究中[6-9]。本文基于FEFLOW軟件的地下水溶質運移模型，并以武漢漢江邊某污染場地為研究對象，對河岸潛流帶非穩定地下水流的溶質運移進行模擬研究。該實例模型探討地表水—地下水相互補給條件下污染物在承壓含水層中的遷移規律，旨在為場地水土環境質量評價及地下水污染治理提供科學依據。

1 研究區概況

污染場地位于武漢市硚口區，處于漢江一級階地上，地勢平坦開闊(圖1)。階地由第四系全新統粉質粘土、粉細砂及含礫中粗砂組成，主要地下水類型為第四系松散巖類孔隙承壓水，承壓含水層厚約為30 m，其上覆相對隔水層為粉質粘土，平均厚度為11.5 m，下伏隔水層為泥巖地層。區內地下水水位隨漢江豐枯水期變化而變化。在豐水期，地表水水位均高于承壓水測壓水位，在水動力驅動下，一方面漢江水向河岸潛流帶補給，使得承壓水壓力增大，測壓水位升高；另一方面存在部分向上越流補給潛水，抬高了潛水位。在枯水期，漢江水位下降幅度大，此時的承壓含水層水位高于江水水位，潛水越流補給承壓含水層水進而又向漢江排泄[10](圖2)。

圖1 研究區平面示意圖

圖2 研究區剖面示意圖及地表水—地下水水位動態變化圖

研究區原為一制藥廠，場地主要污染物包括揮發性有機化合物、總石油烴類，以及硝酸鹽等無機污染物。通過對廠區鉆孔地下水采樣分析檢測顯示，地下水中的揮發性酚類、氨氮、硫酸鹽、氯化物等指標較高，達到了Ⅴ類水標準；在場區中心取樣進行淋濾實驗，結果顯示，表層土壤濾出液的硫酸鹽、硝酸鹽、氯化物的濃度均最高，分別為2 007 mg/L、1 263 mg/L、448 mg/L，均超過地下水質量Ⅲ類標準，說明該場地的土壤包氣帶和地下含水層均受到了一定程度的污染，影響了城市供水安全和周邊的環境健康。

2 地下水流模型

溶質運移的載體是地下水流模型，因此正確的地下水流模型是使溶質運移模型結果可靠、合理的關鍵性因素。非穩定流溶質運移模擬是在一定的初始流場條件下開始進行的，該初始流場通常取自代表該模擬區地下水系統的穩定模擬結果，它可以來自多個鉆孔地下水位的插值，也可以是人為校正后的天然流場，將穩定流模型的邊界賦上隨時間變化的序列值(本文是漢江水位)，運行后便是一種非穩定流。地下水流模型一般包含水文地質概念模型的建立、邊界條件和模型校正等過程。

2.1 水文地質概念模型

根據實際的地下水流場并結合場地地質條件，圈劃出模擬區面積約21 km2(圖1)。模擬區地層由上到下可概化為：雜填土、粉質粘土、粉細砂和泥巖等。雜填土、粉質粘土滲透性較小，模型上將它們歸為一層；粉細砂層作為主要的含水層，為第二層；底部的泥巖可不予模擬，系統自動將其作為不透水地質體。最后將模擬區概化成非均質、各向異性、三維非穩定流的水文地質模型。

2.2 邊界條件與水文地質參數

模擬區南邊以漢江作為給定水頭邊界，北邊作為一個地下水水位基本穩定的邊界，概化為15 m水頭的定水頭邊界，東、西兩側作為流量邊界處理[11-12](圖1)，即Q=KIA，由于軟件可自動計算A，因此只需將KI輸入模型即可，本次模型取I為0.01。上邊界主要接受大氣降水的入滲補給，根據武漢市多年平均降雨量及入滲系數，模擬區地下水補給凈輸入設置為12 mm/a，下邊界泥巖作為底部相對隔水邊界。模擬區初始流場的南部水頭邊界為漢江2015年7月5日的水位數據，運行后形成穩定流，并利用模擬區內鉆孔水位不斷擬合校正，得到了該區校正后的相關水文地質參數：粉質粘土層的Kxx、Kyy、Kzz、給水度、降雨入滲系數分別為0.04 m/d、0.04 m/d、0.01 m/d、0.05、0.01；砂層的Kxx、Kyy、Kzz、給水度分別為5 m/d、5 m/d、0.8 m/d、0.21。

2.3 模型驗證

為了驗證地下水流模型的合理性和有效性，利用正演的方法進行檢驗。首先在穩定流的基礎上，將南邊水頭邊界賦為變化的河流時間序列水位值，并利用2015年8月—2016年1月的場地內CG-1鉆孔的30組地下水位實測值與模擬值進行對比(圖3)?？梢娪^測孔的計算水位與實測水位的誤差在1 m之內，整體變化趨勢基本一致，表明最終選取的參數基本能代表整個含水層的水文地質特征。

3 溶質運移模型

3.1 控制方程和內部參數

污染物在地下水中的運移過程是對流、分子擴散和機械彌散綜合作用的結果，分子擴散和機械彌散對溶質的遷移作用通常難以厘定，習慣上將它們統稱為彌散作用[13]，因此有以下溶質運移模擬的控制方程：

圖3 CG-1孔模擬水位與實際水位對比圖

(1)

式中：等號右端前三項為彌散項，后三項為對流項，最后一項為由于化學反應或吸附解析所產生的溶質增量；Dxx、Dyy、Dzz分別為x、y、z三個主方向的彌散系數；μx，μy、μz分別為x、y、z方向的實際水流速度；c為溶質濃度。

在溶質運移方程中，彌散系數和有效孔隙度是兩個重要的參數，彌散系數除了和含水介質本身性質有關，還依賴于地下水流速等因素，是個多指標變量，于是引入彌散度來反映溶質在地下水中的彌散能力，它和彌散系數之間的關系如下：

(2)

式中：Dd為有效分子擴散系數；δij為Kroneker數；|v|為流速向量v的模；αL、αT分別為縱向彌散度和橫向彌散度；vi、vj為沿i、j方向的流速。在彌散度等參數的選取上，該模型的相關參數采用前人經驗值：粉質粘土層與砂層的彌散度分別為0.006 m、0.5 m，有效孔隙度分別為0.08、0.2。

3.2 初始條件和邊界條件

溶質運移模型的初始條件是描述地下水溶質濃度非穩定變化的初始狀態，是后續模擬進行迭代計算的初始值，因此在賦值時必須加以明確，即使初始溶質濃度為零濃度。邊界條件是模型內部和外界進行溶質交換的一種通道，它可以接受外界的溶質輸入，也可以向外界輸出溶質，因此邊界條件是一個過程量，它伴隨著模擬從開始到結束。在本模型中由于污染場地工廠已搬遷，除原有污染物外無其他污染物輸入，因此假設污染場地均被污染，且濃度分布均勻。根據前期淋濾實驗結果，以硝酸鹽作為典型污染物，以淋濾實驗獲取的最大硝酸鹽淋濾濃度400 mg/L作為初始條件輸入，假設模擬區四周均無污染，并且假設地表水補給時不帶來污染，因此邊界條件不設置濃度邊界。

4 模擬結果及分析

由于區內地下水主要為承壓含水層，且其與地表江水聯系十分密切，因此整個模擬針對承壓含水層而言，即假設硝酸鹽污染在承壓含水層頂部，模擬其未來的變化情況。硝酸根離子易遷移，故整個模擬過程不考慮其被固相吸附[14]。本次模擬考慮兩種情景：非穩定流和穩定流(流場見圖4)。非穩定流符合地下水實際流場，地下水位隨漢江水位的變化而變化；穩定流是假設的情景，其完全由江水補給地下水，南北水頭值分別設置為22.9 m(漢江豐水期平均水位)和15 m。兩種情景均設置了5個監測井用以了解研究區不同部位的污染物遷移分布情況(圖4)，兩種情景的監測井位置完全相同，其中非穩定流井編號為1、2、3、4、5，穩定流編號為1’、2’、3’、4’、5’。1、2、3(1’、2’、3’)號井在場地內部，4(4’)號和5(5’)號在場地外部東北方向。模擬采用自動步長，設置總模擬時長為50年，并單獨對第1、5、10、30和50年后的污染物濃度變化及分布進行討論分析。

圖4 穩定流流場(a)和非穩定流初始流場(b)

4.1 污染物平面分布

非穩定流模擬區污染暈的擴散方向整體呈北東向，和地下水流向基本一致(圖5，圖中為了更好地展現模擬效果，對模擬區進行了適當剪裁，裁剪范圍見圖1)，同時污染物還存在向四周彌散的現象。在模擬時長為1年時，污染物有了一定距離的擴散，約為9 m，但擴散出的濃度并不大，為80 mg/L左右。在模擬5年之后，污染物向外擴散加快，污染場地的濃度在減小，場地西南處的濃度減小速度比場地東北處更快。在模擬期為10年、30年和50年時，污染物向場地東北方向擴散的距離在逐漸增大，分別為74 m、167 m和274 m，與此同時場地內部污染暈中心的污染物濃度和污染暈高值區面積在逐漸減小(圖5-b)。模擬結束時承壓含水層表面污染物的最大濃度為208.71 mg/L，與模擬初始時相比減小了約50%。

圖5 穩定流與非穩定流溶質運移模擬結果對比

與非穩定流的溶質遷移模擬相比較，污染物在穩定流中的遷移速度和濃度減小速度更快，模擬1年的時間里便已經開始出現了較大擴散，擴散距離達到35 m；至50年模擬期結束時，污染物最大濃度為198.48 mg/L，最遠擴散距離約為470 m(圖5-a)。兩種情景在不同模擬時間的最大擴散距離和濃度值見表1。

表1 不同模擬期的最大污染物濃度及擴散距離

非穩定流下幾個監測井的模擬結果顯示(圖6-a)，污染場地內部的井點濃度呈鋸齒狀逐年減小，其中1號井污染物濃度下降最快；2號和3號井下降速度大致相當，并且在50年時趨于平穩；4號井在模擬至2 700 d 時開始監測到污染物，且在10 000 d之前污染物濃度呈現上升的特點，之后增速極為緩慢，說明10 000 d 時污染暈的高值區經過了4號井，但是至模擬期結束，其污染物濃度到達了187 mg/L，且沒有出現下降趨勢，說明這段時間內污染物濃度高值區并沒有完全遷移出4號井；5號井點由于距離污染場地較遠，在模擬期內未檢測出濃度。

穩定流的監測井污染物濃度曲線沒有出現鋸齒狀變化，而是平滑地升高或降低。場地內的3個監測井同樣出現下降，但是它們的下降速度比非穩定流更快，其中1’號井最為明顯，且在11 000 d左右濃度已趨近于零，說明該廠區的污染范圍整體上都進行了一定距離的遷移(圖6-b)；4’號井在模擬1 020 d時監測到了污染物，這一時間比非穩定流模擬提前了近1 700 d，穩定流模擬中4’號井的污染物濃度同樣存在一個先快速上升后穩定的過程，但是在15 000 d時則開始緩慢下降，說明這一模擬污染暈濃度高值區已經遷移出4’號監測井；5’號井在9 000 d左右時也已監測到污染物，隨后污染物濃度逐漸增大，至模擬期結束其濃度達到100 mg/L。

圖6 監測井污染物濃度變化曲線

該模擬中若將初始的硝酸鹽污染物看成是穩定示蹤劑，以場地外的井點為污染物接收點，則污染物的遷移可簡單處理為一個示蹤過程，并由此產生穿透曲線。穿透曲線可以反映投放的示蹤劑在含水介質中的遷移特征。本次模擬單獨對4(4’)號井點的硝酸鹽濃度曲線進行分析(圖6-c)，比較穩定流和非穩定流兩種情況的溶質穿透過程。由于模擬年限和硝酸鹽污染物初始量等原因，穿透曲線未能展示溶質濃度衰減的過程，從穿透曲線中得出了穿透時間等數據(表2)?？梢钥闯?，穩定流的溶質穿透時間更短，且其平均速度比非穩定流的最大遷移速度更大，說明穩定流具有更強的溶質運移能力，溶質容易在含水通道里擴散，而河岸潛流帶中非穩定流對溶質的遷移起到了一定的阻滯作用。和穩定流相比，非穩定流的最大遷移速度和平均速度相差較小，這一現象也可以說明溶質在非穩定流中受到了地下水年際補排關系變化的影響。通過FEFLOW軟件獲取的地下水流速來看，地下水流速度比污染物遷移的平均速度小將近一個量級，表明地下水的滲透流速要遠小于溶質的實際遷移速度，同時也反映了分子擴散在溶質遷移過程中發揮了較大作用。

表2 穩定流和非穩定流條件下穿透曲線的相關參數對比

4.2 污染物垂向分布

為了解硝酸鹽在場地垂向上的遷移過程，對5個不同模擬時間的硝酸鹽濃度剖面圖進行了分析(圖7)，剖面長度1.4 km，剖面方向30°(大致為地下水流方向)，剖面線垂直穿過了整個污染場地。兩種流場的污染物在模擬1年時間里基本到達了含水層底部，但其濃度較小，約為75 mg/L。隨著模擬的進行，污染暈的濃度高值區逐漸向下移動，模擬期結束時，高濃度區均已運移至含水層底部。因此在垂向上，硝酸鹽污染物的遷移總體表現為隨著深度的增加，在模擬初期污染物濃度逐漸減小，而在模擬后期，濃度逐漸增大。模擬期結束時，非穩定流下污染物高濃度區的面積相對較大，說明在垂向上遷移也受到了地下水流的影響，即當污染物在水平方向運移受阻時，便會在垂向上有更多的聚積。因此在沿剖面方向上，非穩定流模擬的污染物擴散距離也小于穩定流情況下的擴散距離，比如模擬30年和50年時，非穩定流的污染物在含水層底部擴散最大距離分別為263 m和384 m，而穩定流情況下的擴散距離為462 m和745 m。

圖7 不同模擬時間的污染物濃度剖面圖(左：非穩定流；右：穩定流)

濃度剖面圖(圖7)顯示隨著深度的增加，污染物沿剖面方向的擴散距離增大，這一現象說明垂向上的溶質鋒面往往最先受水動力影響，發生遷移，使得距離更大。此外也可能有地下水流速差異這一方面的原因，因為在河岸潛流帶內，含水層底部的水頭可能因為深部江水的高壓強而增大，造成了含水層底部水流速度大，淺層流速小，這一解釋在第10年、30年和50年的濃度剖面圖可以很好地體現。當模擬至第10年時，污染物已經到達了含水層底部，且此時底部的擴散距離與淺層相比并不大，也就是說此時污染物失去了溶質鋒面最先被遷移的優勢，深部和淺層的污染物大致在“同一起跑線上”，但是在第30年和50年時，深部和淺層的擴散距離差明顯加大。以非穩定流模擬情況為例，在第30年時，硝酸鹽污染物在含水層頂部和底部的最大擴散距離分別為167 m和269 m，兩者相差102 m；而在第50年時，底部和頂部的污染物最大擴散距離為123 m，比模擬30年時有所增大，表明溶質在垂向上的水平擴散很可能受到了地下水流速的影響。

5 結論

(1)河岸潛流帶的承壓含水層與漢江地表水存在著相互補給的關系，形成非穩定流。在平面上，非穩定流對潛流帶的溶質遷移起到了一定的阻礙作用，比穩定流情況下的溶質更難遷移。模擬時間為50年時，污染物向東北方向的最大遷移距離為274 m，而相同時間下的穩定流模擬為470 m，溶質遷移主要受彌散作用的控制。

(2)在垂向上，污染物的遷移總體表現為隨著深度的增加，在模擬初期濃度逐漸減小，而在模擬后期，污染物濃度逐漸增大。從剖面上來看，污染物在承壓含水層底部的擴散距離比淺層遠，一方面是溶質鋒面最先被運移的結果，另一方面是含水層底部的地下水流速較大帶來的流速差異而造成的。

(3)研究成果可為季節性河流或者水位波動較大的河流階地地下水污染防控提供科學參考，研究地下水和地表水交互補給引起的潛流帶內溶質運移的規律，可以更科學地闡釋河流、湖泊潛流帶中溶質的交換特征，為土壤和地下水修復提供科學依據。