基于GMS的水力控制措施下污染物遷移預測分析

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

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

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

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

因此，本次研究基于GMS地下水數值模擬軟件，結合工程實踐案例，預測了在有局部地下水力控制措施的情況下，滲漏的污染物在天津南部平原區淺部第四系孔隙含水層中的遷移特征，并分析了假設污染滲漏情景對地下水環境可能產生的影響。

1 研究區概況

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

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

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

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

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

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

2 污染遷移模型

2.1 數學計算模型

1) 地下水水流模型

(1) 控制方程

(1)

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

(2) 初始條件

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

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

(2)

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

(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類邊界上的已知水位函數；τ1為I類邊界條件。

② 第II類邊界

(4)

③ 第III類邊界

(5)

2) 地下水溶質運移模型

(1)控制方程

(6)

(2) 初始條件

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

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

(7)

式中：Ω為模擬預測的區域；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)為該邊界上已知的彌散通量函數。

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

(10)

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

2.2 初始與邊界條件

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

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

通過現場踏勘及訪談調查了解，擬建區域現狀的坑塘內水位隨季節有所變化，但變化幅度較小，坑塘水位變幅痕跡垂向距離小于0.2 m，在識別驗證期內，模擬計算范圍內的坑塘水位高程按實測平均值-0.56 m設定為一類水頭邊界；在預測期內按照地下水位導排位置高程-6.70 m設定為一類水頭邊界。

2.3 污染預測情景分析

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

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

本次研究主要考慮該廠區生產區域的滲濾液調節池等地下池體，其污染物濃度高、滲濾液存量大、且污染泄漏不易發現(位于地下、半地下無檢漏設備設施)，存在防滲措施失效等導致的滲漏風險，故將其概化為給定污染物濃度持續輸入的面狀污染源。

2.4 校準后參數

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

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

由于本次模擬研究為假設污染情景的預測，實際上并未發生污染形成污染暈，因此本次驗證僅針對水流模型開展，暫且認為校準之后的水流模型用于污染遷移預測是可行的。

圖1 研究區地下水流場擬合結果

預測期采用研究區降雨和蒸發多年平均值進行計算，多年平均降雨量為526.6 mm，多年平均水面蒸發量1 830.3 mm。

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

表1 校準后的入滲及蒸發參數

滲透系數校正：模型初始值采用幾何平均值0.21 m/d，調參中依據現場抽水試驗成果，潛水含水層滲透系數介于0.17 m/d～0.26 m/d，參考滲透系數等值線插值結果，將滲透系數按照相應范圍進行參數分區，進一步優化調整了參數結構。

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

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

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

彌散度校正：模型初始值統一以經驗公式估算αm=5.3 m。彌散度主要與分選程度有關、均質程度、孔徑迂曲度等有關，分選好的介質彌散度小，分選差的介質彌散度大；均質介質彌散度小，非均質介質彌散度大。因此按各概化層先估算質點遷移距離Ls，再結合成因、顆粒分析，及模型運行結果進行調整。經模型校準后的遷移參數調整如表2所示。

表2 校準后的水流及遷移參數

3 預測結果

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

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

圖2 預測期地下水流場

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

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

第一承壓含水層滲透性明顯高于潛水含水層，污染暈雖然時間上滯后進入，但是進入后的擴展速度相對較快，對流作用對污染暈的分布起明顯主導作用，導致了該層"超標污染區域"范圍小，而“影響區域”范圍大的現象。

4 討論與結果

4.1 討 論

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

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

4.2 結 果

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

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

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

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