林廣宇，符亞兵，焦志亮， 梁玉凱(天津市勘察院， 天津 300191)

天津南部平原區(qū)淺部第四系孔隙介質(zhì)一般以水平向成層規(guī)律明顯的粉質(zhì)黏土、黏土、淤泥質(zhì)土、粉土、粉砂互層組成，以全新統(tǒng)下組沼澤相沉積層(Q41h)粉質(zhì)黏土以上劃為潛水含水層的劃分方法為例，一般在地下水環(huán)境影響評價(jià)的預(yù)測過程中多以“解析法”進(jìn)行溶質(zhì)運(yùn)移的預(yù)測，在這種預(yù)測過程中不考慮含水介質(zhì)的差異性，只能在水平方向上簡單的按平均滲透系數(shù)、縱向彌散系數(shù)進(jìn)行計(jì)算水平一維或二維污染遷移情況，這種方法固然便捷，但其應(yīng)用具有一定限制性，無法獲得溶質(zhì)在垂向的運(yùn)移趨勢、不能反映復(fù)雜參數(shù)條件下的溶質(zhì)遷移特征和規(guī)律[1]。例如本次研究案例，由于所在區(qū)域地下水位埋深較淺，擬建項(xiàng)目運(yùn)營期間為保證地下水面與填埋坑區(qū)域底部第二防滲層保持一定安全距離，在底部第二防滲層的下部設(shè)置地下水導(dǎo)排層，形成了局部人工水力控制措施。水力梯度的改變將會對溶質(zhì)運(yùn)移產(chǎn)生影響[2-3]。

隨著計(jì)算機(jī)軟硬件技術(shù)的普及和進(jìn)步，地下水?dāng)?shù)值模擬技術(shù)已被廣泛地應(yīng)用于地下水量計(jì)算、地下水污染溯源、地下水環(huán)境調(diào)查、評估、修復(fù)等領(lǐng)域，并且獲得了長足的發(fā)展[4]。尤其是近年來，地下水?dāng)?shù)值模擬技術(shù)逐漸成為了解決復(fù)雜地下水環(huán)境問題時(shí)不可或缺的技術(shù)工具[5-6]。國內(nèi)外現(xiàn)今地下水?dāng)?shù)值模擬工作中應(yīng)用的主流商業(yè)軟件包括Visual MODFLOW、Processing MODFLOW、FEFLOW以及GMS等[7]。

地下水?dāng)?shù)值模擬系統(tǒng)(Groundwater Modeling System,GMS)是由美國猶他州楊百翰大學(xué)環(huán)境模型研究實(shí)驗(yàn)室和美國軍隊(duì)排水工程試驗(yàn)工作站在綜合MODFLOW、FEMWATER、MT3DMS、RT3D、MODPATH、UTCHEM 等多個(gè)已有地下水模擬計(jì)算模塊的基礎(chǔ)上，聯(lián)合開發(fā)的一個(gè)綜合性的、用于地下水?dāng)?shù)值模擬的圖形界面軟件[8]。該軟件交互界面優(yōu)良，數(shù)據(jù)處理功能強(qiáng)大，與其他專業(yè)軟件接口較好，可直接利用鉆孔數(shù)據(jù)建立地質(zhì)實(shí)體，提升工作效率與便捷性，滿足研究的要求，故本次選擇應(yīng)用GMS軟件對地下水流場和溶質(zhì)運(yùn)移進(jìn)行數(shù)值模擬研究。

因此，本次研究基于GMS地下水?dāng)?shù)值模擬軟件，結(jié)合工程實(shí)踐案例，預(yù)測了在有局部地下水力控制措施的情況下，滲漏的污染物在天津南部平原區(qū)淺部第四系孔隙含水層中的遷移特征，并分析了假設(shè)污染滲漏情景對地下水環(huán)境可能產(chǎn)生的影響。

1 研究區(qū)概況

天津平原區(qū)潛水在自然條件下側(cè)向徑流較為緩慢，局部受地形、地勢、地表水影響較為明顯；垂向上主要由大氣降水補(bǔ)給、以蒸發(fā)形式排泄，體現(xiàn)為入滲-蒸發(fā)動態(tài)類型。研究區(qū)埋深30.00 m以淺按埋藏條件主要劃分為包氣帶、潛水、第一承壓水及其隔水層。研究區(qū)水文地質(zhì)剖面描述如下：

包氣帶：主要指地下水位以上的雜填土(地層編號①1)、素填土(地層編號①2)、新近沖積層(Q43Nal)黏土為主(地層編號③1)局部夾粉土(地層編號③1-1)透鏡體組成，厚度一般為1.284 m～4.058 m。

潛水含水層：主要由地下水位以下的新近沖積層(Q43Nal)黏土為主(地層編號③1)局部夾粉土(地層編號③1-1)透鏡體、全新統(tǒng)上組陸相沖積層(Q43al)粉質(zhì)黏土(地層編號④1)、全新統(tǒng)中組海相沖積層(Q42m)粉質(zhì)黏土(地層編號⑥1)、粉質(zhì)黏土(砂性大)夾粉土(地層編號⑥3)組成，場地范圍內(nèi)厚度一般為8.05 m～8.95 m。潛水含水層滲透系數(shù)介于0.17 m/d～0.26 m/d，水位埋深一般介于1.519 m～3.954 m。

潛水隔水層：由全新統(tǒng)下組沼澤相沉積層(Q41h)粉質(zhì)黏土(地層編號⑦)、全新統(tǒng)下組陸相沖積層(Q41al)粉質(zhì)黏土(地層編號⑧1)組成，場地范圍厚度一般約為2.50 m～4.90 m。室內(nèi)土工滲透試驗(yàn)結(jié)果介于3.83×10-6cm/s ～3.98×10-6cm/s。

第一承壓水含水層：由全新統(tǒng)下組陸相沖積層(Q41al)粉土(地層編號⑧2)組成，場地范圍厚度一般約為2.80 m～4.70 m。該層滲透系數(shù)一般介于1.12 m/d～2.94 m/d，勘察資料顯示水位埋深一般約3.445 m。

第一承壓水隔水層：由上更新統(tǒng)第五組陸相沖積層(Q3eal)由粉質(zhì)黏土(地層編號⑨1)及上更新統(tǒng)第三組陸相沖積層(Q3cal)粉質(zhì)黏土(地層編號1)組成，場地范圍厚度一般約為10.70 m～15.70 m。室內(nèi)土工滲透試驗(yàn)結(jié)果介于4.78×10-6cm/s ～5.03×10-6cm/s。

2 污染遷移模型

2.1 數(shù)學(xué)計(jì)算模型

1) 地下水水流模型

(1) 控制方程

(1)

式中：W為水流源匯項(xiàng),m3/d；t為模擬預(yù)測的時(shí)間,d；h為地下水水位值,m；μs為含水層儲水率,1/m；Kx,Ky,Kz為不同軸向上的滲透系數(shù)值,m/d。

(2) 初始條件

h(x,y,z,t)=h0(x,y,z)，

(x,y,z)∈Ω,t=0

(2)

式中：h0(x,y,z)為已知水位分布；Ω為模型模擬區(qū)。

(3) 邊界條件

① 第I類邊界

h(x,y,z,t)|τ1=h(x,y,z,t)，

(x,y,z)∈τ1,t≥0

(3)

式中：h(x,y,z,t)為I類邊界上的已知水位函數(shù)；τ1為I類邊界條件。

② 第II類邊界

(4)

③ 第III類邊界

(5)

2) 地下水溶質(zhì)運(yùn)移模型

(1)控制方程

(6)

(2) 初始條件

C(x,y,z,t)=C0(x,y,z),

(x,y,z)∈Ω1,t=0

(7)

式中：Ω為模擬預(yù)測的區(qū)域；C0(x,y,z)為已知的污染物濃度分布情況。

(3) 定解條件

① 第I類邊界條件(給出濃度)

C(x,y,z,t)τ1=c(x,y,z,t),

(x,y,z)∈τ1,t≥0

(8)

式中：τ1為含義為給出濃度的邊界；c(x,y,z,t)為該邊界上的濃度分布情況。

② 第II類邊界條件(給出彌散通量)

(9)

式中：τ2含義為彌散通量邊界；fi(x,y,z,t)為該邊界上已知的彌散通量函數(shù)。

③ 第III類邊界(混合邊界)

(10)

式中：τ3為濃度通量混合邊界；gi(x,y,z,t)為已知的總對流-彌散通量函數(shù)[9]。

2.2 初始與邊界條件

研究區(qū)2006年至2018年間地形地貌等均無明顯變化，淺層地下水受蒸發(fā)、降雨等因素影響，處于隨時(shí)間變化的非穩(wěn)定流狀態(tài)，淺層地下水呈現(xiàn)入滲-蒸發(fā)動態(tài)類型，區(qū)域淺層地下水無人工開采，天然狀態(tài)下地下水的側(cè)向徑流條件差、水量交換弱，因此水流模型中可以將降水入滲作為源，蒸散發(fā)作用作為匯。此外，還考慮污染物發(fā)生泄漏的池體以及填埋坑局部地下水人工導(dǎo)排措施作為水質(zhì)及水流的源匯項(xiàng)。

根據(jù)場地所在區(qū)域2018年全年潛水觀測數(shù)據(jù)，項(xiàng)目所在區(qū)域流場方向主要為由西南向東北，區(qū)域水力坡度約0.16‰～0.25‰，近年來總體趨勢隨時(shí)間無明顯變化。因此，在模擬范圍內(nèi)沿地下水流場方向的上、下游邊界(具有水位觀測井的區(qū)域)設(shè)定為一類水頭邊界，邊界水位賦值以該區(qū)域觀測井實(shí)測水位插值結(jié)果確定；模擬范圍內(nèi)垂直于流場方向(即流場法向量方向)考慮其側(cè)向水量交換較小設(shè)定為二類邊界中的零流量邊界[10]。且由于模型邊界距離模擬的泄漏及遷移核心區(qū)域較遠(yuǎn)，最大程度的降低了邊界條件不確定性對模擬的影響，保證了研究區(qū)的模擬精度。

通過現(xiàn)場踏勘及訪談?wù){(diào)查了解，擬建區(qū)域現(xiàn)狀的坑塘內(nèi)水位隨季節(jié)有所變化，但變化幅度較小，坑塘水位變幅痕跡垂向距離小于0.2 m，在識別驗(yàn)證期內(nèi)，模擬計(jì)算范圍內(nèi)的坑塘水位高程按實(shí)測平均值-0.56 m設(shè)定為一類水頭邊界；在預(yù)測期內(nèi)按照地下水位導(dǎo)排位置高程-6.70 m設(shè)定為一類水頭邊界。

2.3 污染預(yù)測情景分析

項(xiàng)目擬建區(qū)域現(xiàn)狀為未利用地，局部因挖土歷史形成地勢低洼的坑塘，項(xiàng)目建設(shè)方基于節(jié)約建設(shè)資金角度考慮，計(jì)劃利用現(xiàn)狀坑塘經(jīng)改造后作為生活垃圾填埋坑，以減少土方施工量。

由于所在區(qū)域地下水位埋深較淺，擬建項(xiàng)目運(yùn)營期間為保證地下水面與填埋坑區(qū)域底部第二防滲層保持一定安全距離，在底部第二防滲層的下部設(shè)置地下水導(dǎo)排層，填埋坑底設(shè)計(jì)高程約-5.5 m，地下水導(dǎo)排低點(diǎn)高程約-6.7 m，以保證填埋區(qū)域地下水位不高于防滲層，形成了局部人工水力控制措施。

本次研究主要考慮該廠區(qū)生產(chǎn)區(qū)域的滲濾液調(diào)節(jié)池等地下池體，其污染物濃度高、滲濾液存量大、且污染泄漏不易發(fā)現(xiàn)(位于地下、半地下無檢漏設(shè)備設(shè)施)，存在防滲措施失效等導(dǎo)致的滲漏風(fēng)險(xiǎn)，故將其概化為給定污染物濃度持續(xù)輸入的面狀污染源。

2.4 校準(zhǔn)后參數(shù)

建立模型后，通過將各地下水觀測井點(diǎn)的實(shí)際觀測值與模型計(jì)算值進(jìn)行擬合對模型進(jìn)行校正，校正方法為試錯(cuò)法，校正的主要內(nèi)容包括初始條件、邊界條件以及水文地質(zhì)參數(shù)設(shè)置等，以使得校正后的模型更能合理的反映研究區(qū)水文地質(zhì)特征[11]。

擬合校正模型驗(yàn)證期模擬流場與實(shí)測流場擬合如圖1所示。經(jīng)擬合校正，本次模擬的計(jì)算水位與觀測水位隨時(shí)間的變化趨勢基本一致，二者的擬合誤差小于0.5 m(水位變化小于5.0 m時(shí)的允許誤差)，模擬流場形態(tài)與實(shí)際觀測繪制的流場形態(tài)、地下水流向基本一致，滿足模型校準(zhǔn)的誤差要求[12]。

由于本次模擬研究為假設(shè)污染情景的預(yù)測，實(shí)際上并未發(fā)生污染形成污染暈，因此本次驗(yàn)證僅針對水流模型開展，暫且認(rèn)為校準(zhǔn)之后的水流模型用于污染遷移預(yù)測是可行的。

圖1 研究區(qū)地下水流場擬合結(jié)果

預(yù)測期采用研究區(qū)降雨和蒸發(fā)多年平均值進(jìn)行計(jì)算，多年平均降雨量為526.6 mm，多年平均水面蒸發(fā)量1 830.3 mm。

降雨入滲及蒸發(fā)參數(shù)的設(shè)定及調(diào)整參考相關(guān)文獻(xiàn)、技術(shù)手冊，并結(jié)合研究區(qū)的地表植被條件、包氣帶巖性、水位埋深等因素綜合確定[13-14]。由于計(jì)算模擬范圍內(nèi)南部為已建成工業(yè)園區(qū)，園區(qū)內(nèi)大部分進(jìn)行了路面硬化等措施，因此該區(qū)域的降水入滲及蒸發(fā)系數(shù)采用試錯(cuò)法進(jìn)行調(diào)整，最終調(diào)整為工業(yè)園區(qū)區(qū)域按照天然條件下的參數(shù)做一定折減處理后的參數(shù)如表1所示。

表1 校準(zhǔn)后的入滲及蒸發(fā)參數(shù)

滲透系數(shù)校正：模型初始值采用幾何平均值0.21 m/d，調(diào)參中依據(jù)現(xiàn)場抽水試驗(yàn)成果，潛水含水層滲透系數(shù)介于0.17 m/d～0.26 m/d，參考滲透系數(shù)等值線插值結(jié)果，將滲透系數(shù)按照相應(yīng)范圍進(jìn)行參數(shù)分區(qū)，進(jìn)一步優(yōu)化調(diào)整了參數(shù)結(jié)構(gòu)。

給水度校正：給水度僅針對潛水含水層，模型初始值0.05，因潛水含水層為一定程度概化，以粉質(zhì)黏土為主，但夾有粉土互層，調(diào)整后取值0.10[15]。

儲水系數(shù)校正：儲水系數(shù)僅針對具有承壓性質(zhì)的含水介質(zhì)，模型中粉質(zhì)黏土初始值0.002 m-1、粉土初始值0.000 2 m-1，參數(shù)敏感性不高，未作調(diào)整。

有效孔隙度校正：模型初始值以給水度近似替代計(jì)算，潛水含水層調(diào)整同給水度，其余按照土性對應(yīng)經(jīng)驗(yàn)值作為參考。

彌散度校正：模型初始值統(tǒng)一以經(jīng)驗(yàn)公式估算αm=5.3 m。彌散度主要與分選程度有關(guān)、均質(zhì)程度、孔徑迂曲度等有關(guān)，分選好的介質(zhì)彌散度小，分選差的介質(zhì)彌散度大；均質(zhì)介質(zhì)彌散度小，非均質(zhì)介質(zhì)彌散度大。因此按各概化層先估算質(zhì)點(diǎn)遷移距離Ls，再結(jié)合成因、顆粒分析，及模型運(yùn)行結(jié)果進(jìn)行調(diào)整。經(jīng)模型校準(zhǔn)后的遷移參數(shù)調(diào)整如表2所示。

表2 校準(zhǔn)后的水流及遷移參數(shù)

3 預(yù)測結(jié)果

預(yù)測期模擬流場如圖2所示，可以看出，由于建成后填埋坑區(qū)域設(shè)置地下水導(dǎo)排措施，形成了潛水局部人工排泄，導(dǎo)致地下水向該區(qū)域匯集，在此條件下，水流運(yùn)動對污染物的遷移擴(kuò)散將起主導(dǎo)作用。

預(yù)測期7 300 d的各概化層中污染暈分布如圖3所示。預(yù)測期內(nèi)假設(shè)的定濃度持續(xù)補(bǔ)給污染情景對潛水水質(zhì)造成了較為明顯的影響，下游超標(biāo)范圍約150 m(以III類標(biāo)準(zhǔn)值計(jì))，下游影響范圍約180 m(以檢出限計(jì))。污染暈沿主流場方向擴(kuò)展最快，遷移范圍最大，也驗(yàn)證了污染遷移受水力控制為主，受局部導(dǎo)排地下水的措施影響較明顯，而向上游及側(cè)向擴(kuò)展相較于主軸方向并不顯著。隨著各層深度的增加影響依次減弱，濃度和范圍均有所降低。

圖2 預(yù)測期地下水流場

圖3 7 300 d時(shí)地下水污染暈分布情況

為了更好的分析長時(shí)間尺度污染遷移情況，預(yù)測了10 950 d的污染暈垂向等值線剖面如圖4所示。根據(jù)污染暈垂向等值線分布特征，污染物在相對滲透性較差、孔隙度較低的潛水相對隔水層發(fā)生較明顯的阻隔作用，潛水中污染物最大濃度超過250 mg/L，而第一承壓水中污染物濃度降至50 mg/L以下，淺部含水層中污染物垂向濃度梯度大，隨著深度增加而濃度梯度顯著降低。

第一承壓含水層滲透性明顯高于潛水含水層，污染暈雖然時(shí)間上滯后進(jìn)入，但是進(jìn)入后的擴(kuò)展速度相對較快，對流作用對污染暈的分布起明顯主導(dǎo)作用，導(dǎo)致了該層"超標(biāo)污染區(qū)域"范圍小，而“影響區(qū)域”范圍大的現(xiàn)象。

4 討論與結(jié)果

4.1 討 論

本次研究對地下水中對流-彌散作用主導(dǎo)下污染物的遷移規(guī)律進(jìn)行了分析，獲得了較為直觀的結(jié)果表達(dá)，由于本次模擬中并未設(shè)置吸附平衡常數(shù)以及反應(yīng)速率常數(shù)，未考慮含水層介質(zhì)對污染物的吸附作用，未考慮微生物降解以及污染物可能發(fā)生的化學(xué)反應(yīng)，而這些作用在自然界中是客觀存在的，這些都將對污染物遷移轉(zhuǎn)化起到一定影響，使得污染穿透難度增加。但可以預(yù)見，在加入這些作用后，污染暈的范圍(水平和垂向)以及濃度將出現(xiàn)一定程度的縮減。

圖4 10 950 d時(shí)污染暈垂向剖面等值線圖

4.2 結(jié) 果

(1) 由于建成后填埋坑區(qū)域設(shè)置地下水導(dǎo)排措施，形成了潛水局部人工排泄，導(dǎo)致四周地下水向該區(qū)域匯集，可見局部人工水力控制措施對地下水流場起到了明顯控制作用。

(2) 污染遷移水平方向上受水力控制明顯，沿主流場方向擴(kuò)展最快，遷移范圍最大，污染暈20年下游擴(kuò)展最大范圍約150 m。

(3) 污染遷移垂向上受透性相對較差、孔隙度相對較低的隔水層阻隔作用明顯，第一承壓水中污染物濃度降至50 mg/L以下。潛水中污染物濃度高、影響范圍相對較小，而第一承壓水中污染物濃度低，影響范圍相對較大。

(4) 地下水的水力控制措施對污染物遷移起到了較明顯的控制作用，在開展地下水污染管控及修復(fù)治理過程中，可結(jié)合水文地質(zhì)條件采用一定的水力控制措施、或利用相對隔水層的天然屏障，優(yōu)化實(shí)施方案，以更科學(xué)、經(jīng)濟(jì)、合理的應(yīng)對地下水環(huán)境問題。