基于FEFLOW的三維土壤-地下水耦合鉻污染數(shù)值模擬研究

劉 玲 ，陳 堅(jiān) ，牛浩博 ，李 璐 ，殷樂(lè)宜 ，魏亞強(qiáng),3

（1.生態(tài)環(huán)境部環(huán)境規(guī)劃院/長(zhǎng)江經(jīng)濟(jì)帶生態(tài)環(huán)境聯(lián)合研究中心，北京 100012；2.中國(guó)地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；3.上海交通大學(xué)環(huán)境科學(xué)與工程學(xué)院，上海 200240）

污染物下滲及其遷移轉(zhuǎn)化過(guò)程是地下水污染防治工作和地下水環(huán)境研究的焦點(diǎn)，已引起國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[1]。已有學(xué)者將各土壤水運(yùn)動(dòng)模擬軟件與地下水模擬軟件MODFLOW相結(jié)合建立了土壤-地下水水流運(yùn)動(dòng)耦合模型[2-8]。Yakirevich等[9]采用有限差分法建立了擬三維飽和-非飽和水流運(yùn)動(dòng)和溶質(zhì)運(yùn)移耦合模型，有效模擬場(chǎng)地尺度的水流運(yùn)動(dòng)及溶質(zhì)運(yùn)移。林琳等[10]將三維Richards方程分成一維非飽和方程和三維飽和方程，并將飽和模型與非飽和模型在潛水面處耦合建立了基于迦遼金有限單元法的擬三維數(shù)學(xué)模型，為研究大范圍飽和-非飽和水分運(yùn)動(dòng)提供了有效工具。查元源等[11-12]將一維非飽和水流模型同三維地下水模型結(jié)合，模擬了包含大氣邊界的區(qū)域地下水運(yùn)動(dòng)問(wèn)題。并在2014年全面分析了在不同水流條件下各迭代數(shù)值模型的穩(wěn)定性、質(zhì)量守恒、計(jì)算成本以及適用范圍。并釆用線性化的方法，將各種類(lèi)型方程迭代模型轉(zhuǎn)化為非迭代模型，提高了模型的穩(wěn)定性和計(jì)算效率。

由于數(shù)據(jù)采集難度大并且模擬計(jì)算量受限，現(xiàn)有研究多是用非飽和帶模型計(jì)算飽和帶的補(bǔ)給流量，將其代入飽和帶方程計(jì)算新的地下水位，更新非飽和帶的下邊界。這種方式不利于模型之間的計(jì)算反饋，易出現(xiàn)計(jì)算誤差[13]。因此，針對(duì)大范圍場(chǎng)地尺度的三維土壤-地下水水流及溶質(zhì)運(yùn)移模型的構(gòu)建和應(yīng)用顯得尤為迫切。

FEFLOW軟件基于變飽和的控制方程，可以實(shí)現(xiàn)完全三維的飽和-非飽和數(shù)值模擬，但是現(xiàn)有研究多是將FEFLOW應(yīng)用于飽和帶[14-16]。為了精確刻畫(huà)實(shí)際場(chǎng)地土壤-地下水系統(tǒng)中污染物遷移規(guī)律，揭示變飽和反應(yīng)溶質(zhì)遷移模型的參數(shù)敏感性，本文以某鉻污染場(chǎng)地為例，建立了完全三維土壤-地下水水流及溶質(zhì)運(yùn)移模型，采用統(tǒng)一的變飽和水流及溶質(zhì)運(yùn)移控制方程，以水頭作為變量，實(shí)現(xiàn)了大范圍場(chǎng)地的土壤-地下水耦合模擬。

1 研究區(qū)概況

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

研究區(qū)位于湘鄉(xiāng)市丘陵地區(qū)（圖1），多年平均降雨量為1 484 mm，平均蒸發(fā)量為1 142 mm。隸屬華南加里東褶皺帶，地殼運(yùn)動(dòng)包括雪峰、加里東、印支——燕山等。

圖1 研究區(qū)位置圖Fig.1 Location of the study area

研究區(qū)地下水主要為孔隙水和巖溶水，以及少量裂隙水（圖2）。地下水補(bǔ)給來(lái)源有大氣降水、地表水，孔隙水徑流途徑短，大部分地區(qū)就地排泄；裂隙水沿巖石裂隙、孔隙徑流，向地勢(shì)低洼的溝谷排泄；巖溶水沿地下管道、裂隙徑流，向地勢(shì)低洼的溝谷排泄。

圖2 研究區(qū)水文地質(zhì)平面圖（本圖來(lái)源為全國(guó)地質(zhì)資料館1∶20萬(wàn)水文地質(zhì)圖G4902幅）Fig.2 Hydrogeological map of the study area (The source of this figure is G4902 1∶200 000 hydrogeological figure of the national geological data center)

研究場(chǎng)區(qū)內(nèi)地層由上至下共分為5層：①雜填土（Q4ml）：由粉質(zhì)黏土、碎石構(gòu)成，厚3～4 m，孔隙度較大，透水性佳；②粉質(zhì)黏土（Q4al）：以黏土為主，含少量砂，厚2～4 m，透水性差，含水量較少，為相對(duì)阻水層；③中-粉砂（Q4al）：由砂和黏土構(gòu)成，孔隙發(fā)育，厚1 m左右，含水量較少；④圓礫（Q4al）：由砂礫構(gòu)成，厚3～4 m，孔隙水發(fā)育，以側(cè)向補(bǔ)給為主；本層為模型的主要含水層；⑤泥質(zhì)粉砂巖（Kdnd）：由砂泥構(gòu)成，深度約9～10 m，弱透水，為隔水底板（圖3）。

圖3 場(chǎng)地水文地質(zhì)剖面圖Fig.3 Hydrogeological profiles of the site

1.2 場(chǎng)地污染特征

本研究結(jié)合省地方標(biāo)準(zhǔn)與場(chǎng)地的實(shí)際情況，采用重金屬污染場(chǎng)地修復(fù)標(biāo)準(zhǔn)DB43/T1165-2016的工業(yè)用地標(biāo)準(zhǔn)作為土壤污染物濃度的參考標(biāo)準(zhǔn)值。按照調(diào)查技術(shù)規(guī)范要求，調(diào)查共布設(shè)297個(gè)土壤采樣點(diǎn)、45個(gè)地下水采樣點(diǎn)。場(chǎng)地污染情況如表1所示。

表1 六價(jià)鉻檢測(cè)結(jié)果Table 1 Statistics of soil test results

2 研究方法

2.1 數(shù)學(xué)模型構(gòu)建

2.1.1 水流數(shù)學(xué)模型

以水頭h作為主要變量，三維變飽和多孔介質(zhì)中的水流運(yùn)動(dòng)控制方程為[17]：

式中：S0——貯水系數(shù)/m-1；

h——水力水頭/m；

?h——水力梯度，無(wú)量綱；

Sw—— 飽和度，非飽和帶為水頭h的 函數(shù)，0＜Sw≤1，飽和帶Sw=1，無(wú)量綱；

krw—— 相對(duì)滲透率，非飽和帶krw是飽和度的函數(shù)，由Van Genuchten-Mualem模型計(jì)算可得；飽和帶為1，無(wú)量綱；

K——介質(zhì)滲透系數(shù)/（m·s-1）；

t——時(shí)間/s；

V——達(dá)西流速/（m·s-1）；

φ——孔隙度，無(wú)量綱；

e——單位法向量，無(wú)量綱；

Q——流量系統(tǒng)的源匯項(xiàng)/s-1；

χ——浮力系數(shù)，表示流體的密度效應(yīng)，無(wú)量綱；

?V——流速梯度/（s-1）。

2.1.2 溶質(zhì)運(yùn)移數(shù)學(xué)模型

三維變飽和溶質(zhì)運(yùn)移控制方程為[18]：

式中：c——污染物源濃度/（mg·L-1）；

c*——滲入地下水中的污染物濃度/（mg·L-1）；

Dij——水動(dòng)力彌散系數(shù)/（m2·s-1）；

D0——源溶液的分子擴(kuò)散系數(shù)/（m2·s-1）；

|V|——達(dá)西流速的絕對(duì)值/（m·s-1）；

R——阻滯系數(shù)，表征土壤吸附能力，無(wú)量綱；

ρs——土壤顆粒密度/（kg·m-3）；

kd——分配系數(shù)/（L·kg-1）；

λ——反應(yīng)系數(shù)/（10-4·s-1），本文中主要考慮土壤及地下水中鉻的氧化還原反應(yīng)；

αL、αT——縱向和橫向彌散度/m；

τ——彎曲度，由Millington-Quirk方程可得，無(wú)量綱；

qc——邊界濃度通量/（mg·m-2·s-1）；

ni——邊界上的單位法向量，無(wú)量綱；

xi、xj—— 空間坐標(biāo)（xi、xj=x,y,z）/m；

Vi、Vj——xi、xj方向的達(dá)西流速/（m·s-1）。

2.2 耦合方法及數(shù)值求解過(guò)程

2.2.1 耦合方法

由于水流和溶質(zhì)運(yùn)移控制方程均與流體密度相關(guān)密切，因此，可將水流和溶質(zhì)運(yùn)移通過(guò)表示流體密度效應(yīng)的浮力項(xiàng)非線性耦合，利用FEFLOW軟件實(shí)現(xiàn)土壤-地下水耦合數(shù)值模擬。

模型在非飽和帶基于以下假設(shè)[19]：

（1）氣液分離，唯一動(dòng)態(tài)相是液相；

（2）流體運(yùn)動(dòng)符合達(dá)西定律，在模擬中不斷根據(jù)瞬時(shí)條件改變潛水面的位置。

2.2.2 數(shù)值求解過(guò)程

采用Galerkin有限單元（FEM）網(wǎng)格及上游加權(quán)法對(duì)水流及溶質(zhì)運(yùn)移的控制方程進(jìn)行求解，采用預(yù)處理共軛梯度法（PCG）求解水流控制方程，用預(yù)處理正交最小化（ORTHOMIN）法求解溶質(zhì)運(yùn)移控制方程[20]。

2.3 數(shù)值模型構(gòu)建

2.3.1 概念模型

研究區(qū)面積約為11.16 km2，由于研究區(qū)不是完整的水文地質(zhì)單元，因此上游沒(méi)有完整的自然邊界。根據(jù)水文地質(zhì)調(diào)查結(jié)果，將漣水河西北側(cè)約2.5 km處等水位線作為模擬區(qū)上游邊界，下游以漣水河為邊界。垂向上，以地表作為模擬的上邊界，底部泥質(zhì)粉砂巖為下邊界，總深度約為23 m。根據(jù)研究場(chǎng)地滲透系數(shù)不同（表2），將模型概化為5層，由于第一層中碎石土主要存在于上部1 m范圍內(nèi)，導(dǎo)致第一層上部和下部孔隙度不同，因此將模型第一層細(xì)化為2個(gè)亞層，每層厚度由軟件插值可得。區(qū)內(nèi)地下水主要接受降雨及地表水的補(bǔ)給，流向?yàn)樽晕鞅毕驏|南流向漣水河。

表2 研究區(qū)水文地質(zhì)參數(shù)取值表Table 2 Values of hydrogeological parameters in the study area

污染物主要是由于上層渣土長(zhǎng)期淋溶導(dǎo)致六價(jià)鉻不斷進(jìn)入含水層所導(dǎo)致。六價(jià)鉻和三價(jià)鉻在運(yùn)移過(guò)程中可能會(huì)發(fā)生相互轉(zhuǎn)化。本次模擬將場(chǎng)地調(diào)查獲取的濃度數(shù)據(jù)作為模擬的初始條件，在污染源全部清理后，無(wú)人工干預(yù)的情況下，模擬土壤和地下水中六價(jià)鉻的時(shí)空變化趨勢(shì)。

2.3.2 初始和邊界條件

初始條件：模型初始水位為48 m，土壤含水量由模型計(jì)算可得。

垂直邊界：模型上邊界概化為降雨入滲、蒸發(fā)邊界；本文中補(bǔ)給由年平均降雨量與蒸發(fā)量差值計(jì)算所得，設(shè)為343 mm/a；下邊界為隔水邊界。

側(cè)向邊界：上游邊界概化為定水頭邊界，水頭最大值為56 m，最小值為54 m；下游（漣水河）概化為定水頭邊界，水頭最大值為45 m，最小值為37 m。

溶質(zhì)邊界：本次模擬中為自由出流邊界。

2.3.3 模型參數(shù)選取

根據(jù)水文地質(zhì)條件將研究區(qū)概化為非均質(zhì)、各向異性三維非穩(wěn)定滲流系統(tǒng)。通過(guò)土壤試驗(yàn)及相關(guān)文獻(xiàn)獲得土壤相關(guān)參數(shù)[21]，含水層的水文地質(zhì)參數(shù)主要通過(guò)前期的水文地質(zhì)調(diào)查獲得（表2）。

溶質(zhì)運(yùn)移考慮對(duì)流、彌散、吸附、反應(yīng)衰減等過(guò)程。彌散系數(shù)與土壤性質(zhì)有關(guān)，通過(guò)試驗(yàn)由“三點(diǎn)公式”求得[22-23]?？v向彌散度為彌散系數(shù)與平均孔隙水流速度之比，橫向彌散度取縱向彌散度的1/5。根據(jù)線性等溫平衡吸附規(guī)律，利用試驗(yàn)確定六價(jià)鉻的分配系數(shù)，阻滯系數(shù)為分配系數(shù)與介質(zhì)密度的乘積（式4）。反應(yīng)衰減主要考慮六價(jià)鉻的氧化還原反應(yīng)，反應(yīng)系數(shù)在野外分析的基礎(chǔ)上結(jié)合了前人的研究資料進(jìn)行整理分析綜合確定[24-29]（表3）。

表3 溶質(zhì)運(yùn)移模型參數(shù)取值表Table 3 Parameter values of the solute transport model

2.4 模型識(shí)別驗(yàn)證

模型參數(shù)校驗(yàn)主要采用“試錯(cuò)法”進(jìn)行調(diào)整。初始流場(chǎng)是將研究區(qū)參數(shù)初始值輸入模型，經(jīng)過(guò)穩(wěn)定流計(jì)算得到天然流場(chǎng)，然后根據(jù)實(shí)際觀測(cè)水位對(duì)天然流場(chǎng)進(jìn)行參數(shù)校正，得到校正后的地下水初始流場(chǎng)（圖4），將模擬值與實(shí)際值做擬合（圖5），得到其相關(guān)系數(shù)為0.999 6，模擬值與實(shí)際值相比，均方差為0.092 1，均方根誤差為0.303 5，表明模型初始流場(chǎng)基本符合研究區(qū)實(shí)際水文地質(zhì)條件，反映了實(shí)際流場(chǎng)特征，故可利用該模型得到的流場(chǎng)作為非穩(wěn)定流的初始流場(chǎng)并以此為基礎(chǔ)進(jìn)行溶質(zhì)運(yùn)移模擬。

圖4 研究區(qū)模擬流場(chǎng)圖與實(shí)際流場(chǎng)圖Fig.4 The simulated flow field and actual flow field in the study area

圖5 模擬值與實(shí)際值擬合圖Fig.5 Fitting diagram of the simulated value and actual value

3 結(jié)果和討論

3.1 模擬預(yù)測(cè)結(jié)果

根據(jù)場(chǎng)地污染現(xiàn)狀（圖6），上層土壤污染暈中心濃度最高值為3 410 mg/L，下層土壤為3 430 mg/L（第一層頂?shù)装澹?，地下水（第四層底板）污染暈中心濃度最高值?09 mg/L，并分別在主要含水層底板位于場(chǎng)區(qū)周?chē)鷧^(qū)域、場(chǎng)地下游方向以及漣水河西側(cè)河岸設(shè)置4處觀測(cè)點(diǎn)監(jiān)測(cè)濃度隨時(shí)間的變化，分別為CG-1、CG-2、CG-3、CG-4。

圖6 初始污染羽分布圖Fig.6 Distribution figures of the initial contaminate plume

3.1.1 土壤模擬預(yù)測(cè)結(jié)果

土壤初始污染暈自二分廠向四周呈遞減趨勢(shì)。在降水淋濾作用下，土壤中六價(jià)鉻不斷向下運(yùn)移，并向四周擴(kuò)散。第28天時(shí)模型第三層下部六價(jià)鉻濃度達(dá)到最大值（圖7a），污染暈中心濃度約為3 382.94 mg/L。此后，土壤中六價(jià)鉻濃度逐漸降低。第378天時(shí)土壤中污染暈的分布范圍達(dá)到最大（圖7b），場(chǎng)地東南側(cè)遷移距離最大，約為300 m，北側(cè)最短，約為90 m，此時(shí)，污染暈中心濃度為357.94 mg/L。而后污染暈范圍逐漸減小，在1 900 d時(shí)土壤中污染暈已完全消失。

圖7 土壤中（第三層）六價(jià)鉻的運(yùn)移模擬預(yù)測(cè)結(jié)果Fig.7 Simulation and prediction results of hexavalent chromium transport in soil (the third layer)

3.1.2 飽和地下水模擬預(yù)測(cè)結(jié)果

飽和帶初始濃度中心為東、西渣場(chǎng)。地下水流動(dòng)是六價(jià)鉻在含水層運(yùn)移的主要驅(qū)動(dòng)力，結(jié)合研究區(qū)地下水流場(chǎng)可知污染物主要向污染場(chǎng)地東南側(cè)的漣水河遷移（圖8），第10天時(shí)運(yùn)移至場(chǎng)地南側(cè)邊界GC-2附近，此時(shí)污染暈中心濃度值為53.01 mg/L，第25天時(shí)運(yùn)移至東北側(cè)邊界GC-1附近，污染暈中心濃度值為422.45 mg/L，第49天時(shí)主要含水層六價(jià)鉻濃度達(dá)到最大值，為1 533.65 mg/L。在此之前，由于土壤中六價(jià)鉻不斷下滲，東、西兩渣場(chǎng)以及兩場(chǎng)之間區(qū)域濃度均較高，而后隨著入滲含水層中污染物濃度逐漸降低。在第585天時(shí)污染羽遷移至漣水河，污染暈中心濃度為25.79 mg/L，第917天時(shí)主要含水層污染暈分布范圍最大，此時(shí)污染暈中心位于西渣場(chǎng)，濃度為7.91 mg/L，而后污染暈逐漸減?。ū?）。第2 000 天時(shí)僅在場(chǎng)地邊界處存在小范圍濃度超標(biāo)區(qū)域。

表4 污染暈面積表Table 4 Contaminant halo area

圖8 地下水六價(jià)鉻運(yùn)移模擬預(yù)測(cè)結(jié)果Fig.8 Simulation and prediction results of hexavalent chromium migration in groundwater

土壤中六價(jià)鉻在場(chǎng)地附近以垂向運(yùn)移為主，水平運(yùn)移速率較??；隨著埋深的增加，水平運(yùn)移逐漸加強(qiáng)（圖9）。由于土壤中六價(jià)鉻的下滲，地下水中的污染物濃度先是逐漸增大，當(dāng)污染物全部進(jìn)入含水層后，地下水中污染物濃度開(kāi)始逐漸減小。土壤中的六價(jià)鉻約需10 h下滲至潛水面處（圖10a），全部污染物下滲至隔水底板約需2 d（圖10b）。

圖9 污染物遷移至漣水河時(shí)污染羽三維分布圖Fig.9 3D distribution of the contaminant plume when contaminants migrate to the Lianshui River

圖10 場(chǎng)地六價(jià)鉻垂向濃度模擬結(jié)果Fig.10 Simulation results of chromium vertical concentrations at the site

3.2 參數(shù)敏感性分析

土壤是地下水污染的重要媒介，污染物主要是通過(guò)大氣降水、地表水或灌溉水的入滲淋濾從土壤進(jìn)入地下污染地下水。因此本文將溶質(zhì)運(yùn)移模型簡(jiǎn)化，考慮初始污染物僅存在于上層土壤的情況對(duì)參數(shù)敏感性進(jìn)行分析。

已有研究[30]表明降雨量、阻滯系數(shù)以及反應(yīng)常數(shù)會(huì)對(duì)污染物運(yùn)移產(chǎn)生影響。本文分別改變降雨量、阻滯系數(shù)和反應(yīng)常數(shù)研究其影響。

3.2.1 降雨量對(duì)污染物運(yùn)移的影響

為了討論降雨入滲變化對(duì)污染物運(yùn)移的影響，本文將模型降雨入滲條件設(shè)置為時(shí)變的情況。依據(jù)研究區(qū)近年的降雨量和蒸發(fā)量統(tǒng)計(jì)資料，將豐水期（每年的7—8月）入滲補(bǔ)給量設(shè)為200 mm/a，枯水期（每年的9月至次年6月）設(shè)為15 mm/a。

結(jié)果顯示，每一個(gè)豐水期過(guò)后約15 d，非飽和帶土壤底部的污染物高濃度區(qū)域面積均會(huì)有所回升，上升約8 d后，含水層達(dá)到新的收支平衡，此后非飽和帶中污染物濃度再次開(kāi)始下降（圖11）。主要是由于降雨首先要補(bǔ)給包氣帶，地下水對(duì)降雨補(bǔ)給的響應(yīng)存在滯后性。當(dāng)補(bǔ)給到達(dá)含水層之后，地下水位上升，非飽和帶中污染物濃度隨之增大，隨著污染物的不斷下滲以及降雨補(bǔ)給的持續(xù)淋洗，非飽和帶和飽和帶中的污染物濃度均逐漸降低。

圖11 雨季后土壤中污染物濃度變化Fig.11 Changes in concentrations of pollutants in soil after rainy season

3.2.2 阻滯系數(shù)對(duì)溶質(zhì)運(yùn)移的影響

本文中阻滯系數(shù)是實(shí)驗(yàn)室測(cè)定的六價(jià)鉻的分配系數(shù)與固體體積百分?jǐn)?shù)的乘積，無(wú)量綱，與含水率密切相關(guān)。飽和含水率為定值，而非飽和帶含水率變化較大，因此本文僅考慮非飽和帶土壤阻滯系數(shù)對(duì)溶質(zhì)運(yùn)移的影響。

由圖12、圖13可以看出，當(dāng)阻滯系數(shù)相差兩個(gè)數(shù)量級(jí)時(shí)，曲線變化明顯，污染物運(yùn)移速度變緩，同一距離處污染物的濃度降低。阻滯系數(shù)的大小表征土壤顆粒對(duì)污染物吸附能力的強(qiáng)弱，阻滯系數(shù)越大，吸附作用對(duì)污染物的阻滯作用就越大，污染物遷移越慢。阻滯系數(shù)對(duì)污染物運(yùn)移的影響可能由于研究區(qū)土壤中含有黏性土，黏性土顆粒高度分散，電荷不均衡，表面能較大，可以吸附水中的六價(jià)鉻[31]。另外，在較大時(shí)空的模擬中，阻滯系數(shù)一般是由實(shí)驗(yàn)室測(cè)定的分配系數(shù)通過(guò)溶質(zhì)運(yùn)移控制方程中的阻滯系數(shù)計(jì)算公式求得，然而通過(guò)這樣計(jì)算求得的阻滯系數(shù)值往往偏高，可能會(huì)導(dǎo)致模擬結(jié)果不夠精準(zhǔn)，因此在有條件的情況下，應(yīng)盡可能采用試驗(yàn)直接測(cè)定阻滯系數(shù)的方式進(jìn)行模擬[32]。

圖12 不同阻滯系數(shù)下場(chǎng)區(qū)至河流濃度隨距離的變化Fig.12 Variation of concentration with distance under different retardation coefficients from site to river

圖13 不同阻滯系數(shù)河岸處觀測(cè)點(diǎn)濃度隨時(shí)間的變化Fig.13 Variation of concentration with time at the observation point on the bank with different retardation coefficients

3.2.3 反應(yīng)常數(shù)對(duì)污染物運(yùn)移的影響

文中反應(yīng)常數(shù)是指六價(jià)鉻被還原成三價(jià)鉻的速率，單位為10-4s-1。由圖14、圖15可見(jiàn)，當(dāng)反應(yīng)常數(shù)增大兩個(gè)數(shù)量級(jí)時(shí)，曲線變化及其明顯，污染物運(yùn)移速率變緩，同一距離處濃度顯著降低。研究發(fā)現(xiàn)[33]，正常pH值條件下，在天然水體中，三價(jià)鉻和六價(jià)鉻之間可以相互轉(zhuǎn)化，六價(jià)鉻可被某些還原性物質(zhì)還原為三價(jià)，三價(jià)鉻也可被氧化成六價(jià)。六價(jià)鉻極易溶于水，而三價(jià)鉻在水溶液中則以難溶的沉淀形式存在。因此，當(dāng)六價(jià)鉻被還原成三價(jià)鉻時(shí)，形成沉淀后析出，地下水中的六價(jià)鉻濃度降低。

圖14 不同反應(yīng)常數(shù)下場(chǎng)區(qū)至河流濃度隨距離的變化Fig.14 Variation of concentration with distance under different reaction constants from site to river

圖15 不同反應(yīng)常數(shù)河岸處觀測(cè)點(diǎn)濃度隨時(shí)間的變化Fig.15 Variation of concentration with time at the observation point on the river bank with different reaction constants

4 結(jié)論

（1）包氣帶的存在使地下水對(duì)降雨的響應(yīng)存在滯后性。土壤和地下水中的濃度均呈現(xiàn)先升高后降低的趨勢(shì)，場(chǎng)地邊界處土壤中的污染物濃度比地下水中的污染物濃度低。在模擬期內(nèi)的夏季降雨量增加期間地下水中污染物濃度呈現(xiàn)明顯的下降趨勢(shì)，之后緩慢回升；而土壤中的污染物濃度則在降雨量增大時(shí)呈現(xiàn)上升趨勢(shì)，隨后緩慢回落。降雨會(huì)引起潛水面波動(dòng)，可能導(dǎo)致地下水中的污染物進(jìn)入土壤，從而影響土壤中的污染分布。

（2）阻滯系數(shù)和反應(yīng)常數(shù)的影響均較大，當(dāng)變化兩個(gè)數(shù)量級(jí)時(shí)，曲線變化明顯；但當(dāng)反應(yīng)常數(shù)僅增大至10-6時(shí)，遷移出場(chǎng)區(qū)的污染物濃度約減少2 000 mg/L，較難遷移至漣水河，影響最為強(qiáng)烈。但是，在較大時(shí)空的模擬中，往往很難將阻滯系數(shù)和反應(yīng)常數(shù)的影響分開(kāi)，在今后的研究中需要探索更為有效的方式將二者分別討論。

（3）采用基于FEFLOW的數(shù)值模型，能夠解決各系統(tǒng)之間交互性差的問(wèn)題，很好的模擬潛水面在土壤和地下水系統(tǒng)之間的波動(dòng)，以及溶質(zhì)在二者之間的運(yùn)移，提供較為精確的模擬結(jié)果，實(shí)現(xiàn)了大范圍場(chǎng)地的完全三維土壤-地下水耦合模擬，為提出科學(xué)有效的地下水污染防控措施提出定量化的理論依據(jù)。