針對城市污水管網滲漏預測的數字模擬研究

宋雷震

(淮南聯合大學 智能制造學院，安徽 淮南 232038)

城市化進程給人們生活帶來便利的同時，也引發城市污水管網外滲污染問題，同時建筑工程開發過程中的基坑降水也導致了地面沉降和地下水降落漏斗[1]等問題。地下水模擬系統(Groundwater Modeling System，GMS)能夠對地下水位進行數值模擬[2]。GMS軟件是在現有地下水模型的基礎上建立的全面綜合性的圖形界面軟件，是目前已經被廣泛應用于世界范圍內的地下水模擬方法[3]，它具備強大的數據分析處理功能、方便快捷的操作界面、三維可視化的效果。本文利用GMS軟件對城市污水管網外滲污染地下水進行地下水流和水質的數值模擬分析，并運用溶質運移模型預測排管泄漏的情況。

1 管網外滲污染地下水的數值模型

1.1 水文地質條件和地下水環境

本研究選取某市區域對城市污水管網外滲污染地下水情況進行數值模擬和預測。該地區氣候為溫帶濕潤季風氣候，海拔高度約為60 m，相對高度差為8，全年平均氣溫為6.3 ℃，最低氣溫和最高氣溫分別為-31 ℃和35 ℃，無霜期和封凍期均約為150 d，全年平均日照時間約為2 700 h[4]。采集氣象站該區域2005—2018年的年降水量進行處理分析，分別得到偏態系數、變異系數、降水量均值3個統計參數，數值分別為0.30，0.75，656.68 mm。2008—2018年系列降水頻率曲線如圖1所示。從圖1可以看出，降水頻率范圍為3%～96%且集中分布在15%～85%。當降水頻率為75%時，相應年平均降水量為524.7 mm。年降水量和頻率可以擬合為一條光滑的曲線。

圖1 2005—2018年系列降水頻率曲線

研究區域的地下水廣泛存在于圓礫土層、粉細砂、中粗砂、粉細砂中。自然狀態下水位年幅度變化為3 m，同時基坑降水產生降落漏斗，該漏斗以城區為中心，面積約為16 km2。整個研究時間段內，地下水埋深為1.49～10.25 m，溫度約為8～12 ℃，水質主要含有重碳酸鈣鈉和重碳酸鈣，礦化度為0.5 g/L，鐵離子含量豐富，質量濃度通常為2～11 mg/L，部分地段濃度更高。研究區域處于建設階段，排泄方式為基坑降水、外圍機井灌溉、植物蒸騰作用。地下水通過徑流直接流入研究區域，補給源為側向徑流、水田區灌溉回滲水、降雨[5]。城區內地下水類型為人工開采型，城區外為氣候型。利用滲水試驗得到該區域包氣帶的滲透系數范圍為0.08～1.44 m/d，滲透深度范圍為40～140 m。含水層水文地質參數值如表1所示，展示9種鉆孔編號在不同降深下的不水文參數值。通過附加人工流場下穩定流抽水的野外彌散試驗可知，縱向和橫向彌散系數分別為0.5 m2/d和0.015 m2/d，縱向彌散度和橫向彌散度分別為0.212 8 m和0.006 4 m。

表1 水文地質參數情況

依據GB/T 14848—2017《地下水質量標準》對研究區域的地下水水質進行評價，通過指標值的限值范圍得到地下水質量情況，檢測得到地表水質量情況如表2所示。因此，研究區域菌群總數、總大腸菌群、鋁、錳、鐵等指標均超出集中式飲用水源限定標準，該區域的地下水質情況被人類生產生活影響。通過地下水水位等值線圖，分析得到該區域形成以2個施工工地為中心的沉降漏斗，中心水位埋深為9.14～10.20 m。但區域歷史資料顯示漏斗中心下降的深度約為6 m，地下水流場從天然狀態變成從4個方向流向新城。

表2 地表水質情況

1.2 地下水流和水質數學模型

GMS模擬系統包括UTCHEM、RT3D、SEEP2D、MT3DMS、MODFLOW等計算模塊和Solid、MAP、PEST等輔助模塊。研究依據區域含水層情況，將其設置為松散巖類孔隙微承壓和潛水層含水層。設置地下水水力坡度在2/1 000以內，水流運動方式為平面二維流，水力特征遵從達西定律，同時當前和預期內水位均為非穩定流運動。研究依據建設區階段和建設完成后等水位情況，均將研究區域邊界設置為二類邊界[6-7]。整體區域面積23.6 km2，人為干擾下和天然狀態下邊界條件劃分情況如圖2所示。根據圖2中線條和箭頭所表示的河流邊界及水流方向，將區域劃分為I1和I22個部分。

(a)人為干擾下邊界條件

結合前述分析，研究建立非勻質各向同性松散巖類孔隙微承壓水和潛水的地下水流模擬數值模型。孔隙潛水和孔隙微承壓水的數學模型計算公式分別如式(1)和式(2)所示。

(1)

(2)

式中：K、T分別表示潛水含水層的滲透系數和承壓含水層導水系數；μ*和μ表示承壓含水層彈性釋水系數和潛水含水層給水度；H、H0、h分別表示承壓水水頭、潛水水位、地下水位；B表示含水層底板標高，Qd和Qr表示排泄強度和入滲補給強度；Qi表示大井開采量；h0和h1表示初始水位和一類邊界點的水位；q表示二類邊界單寬流量；x，y為坐標；D表示計算區范圍；Γ1和Γ2分別表示一類邊界和二類邊界；t表示時間；n表示開采井總數。

通過一類、二類邊界條件、初始條件、偏微分方程即可求解問題。具體求解過程為：首先在區域范圍進行線性插值和三角形剖分，利用伽遼金有限單元法離散數學模型為單元方程組，即可通過計算機程序求解問題。全部等值線圖利用Surfer軟件繪制，經分析檢測沒有問題的情況下，通過GMS軟件實現數據自動采集和單元自動剖分。區域剖分使用Modflow模塊，含水層分為2層，每層含有30 000單元，單元平均面積為1 405 m2，活動單元格16 807個，該種求解方法在保證計算準確率的前提下，也能夠提升運行速度[8-9]。鑒于研究區域水位信息統計資料只有2019年1月的數據，因此不能實現模型識別和檢測。研究利用區域平水年的源匯項進行模型識別，假如地下水位通過平水年的源匯項運行呈現穩定的變化狀態，則模型參數設置合理，具有很好的可行性。在此基礎上，初始流場設置為地下水流場，接著將2008年作為模型的源匯項，連續運行10 a。1個應力期為1個月。降水入滲補給量的計算公式為：

Qpr=0.1×α×P×F。

(3)

式中：α和P分別表示降水量和降水入滲補給系數；F表示一分區面積減去水面的面積。側向徑流補給量為：

Qlr=36.5×K×I×M×L。

(4)

式中：K和I表示含水層滲透系數和垂直于剖面方向的水力坡度；M和L分別表示含水層厚度和剖面長度。研究分析水位低于4 m的蒸發排泄情況，潛水蒸發排泄量Qe為：

Qe=c·F·E601。

(5)

式中：c表示潛水蒸發系數；F表示計算面積；E601表示水面蒸發量。在此基礎上，利用MT3DMS模塊構建地下水溶質運移模型，完成對污染的預測，選用的模擬因子為化學耗氧(O2計)[10]。二維非穩定流溶質運移方式下彌散數學模型為：

(6)

求解問題采用3rd-orderTVD方法和有限差分析方法。I1區域預測5 a后完成建設，因此5 a內地下水流場仍然會受到人工因素影響，研究選取2種特征污染物完成管網外滲過程，設置模型運行10 a。

2 外滲污染地下水數值模擬分析

2.1 水質模型識別結果

研究首先沿驗證所提出地下水流數值模型的可行性，2層區域的初始設置值如下：第1層分區I1區域和I2區域面積依次為12.25 m2和11.35 m2，滲透系數依次為14.43 m/d和24.48 m/d，給水度分別為0.10和0.12。第2層分區I1區域和I2區域面積與第1層分區一致，滲透系數依次為80 m/d和100 m/d，給水度分別為0.18和0.20，彈性釋水系數均為0.000 1。經側向徑流補給量計算可知，滲透系數為14.43 m/d，水力坡度為2%，含水層厚度為30 m，側向補給值為688.12萬m3/a。農業灌溉回滲量計算可知，農業灌溉量和回滲量為134.46萬m3/a和20.17萬m3/a。人工開采包括基坑降水和農業灌溉，總值為836.2萬m3/a。

I1區域的降水入滲系數弱于I2區域，這是由于I1區域比I2區域早建設完成，大面積硬化路面覆蓋導致降水入滲情況較差，而I2區域由于還存在大規模的耕地，因此降水入滲情況更好。2個區域的降水量補給計算結果如圖3(a)所示。I1區域和I2區域降水入滲補給量分別為38.7萬m3/a和107.58萬m3/a。潛水蒸發排泄量情況如圖3(b)所示，I1和I2的面積分別為1.8 km2和0.11 km2，地下水蒸發量分別為21.78 萬m3/a和1.33 萬m3/a，水面蒸發量均為1 210 mm，蒸發系數均為0.1。最終地下水補給資源總量和地下水總排泄量達到動態平衡，補給誤差僅為4.74 萬m3/a。

(a)降水入滲補給資源量

漏斗中心和QT05標記處的地下水位變化情況分別如圖4(a)和4(b)所示。從2008—2018年整體來看，2個地段的地下水位表現出動態平衡的現象，同時實際動態情況也和變化幅度接近，因此模型的參數設計較為合理，所建立的數學模型和水文地質模型可靠性較高。

(a)漏斗中心

2.2 溶質運移模型結果

通過GMS軟件中的MT3DMS模塊完成三維地下水流數值模擬。結合尺度效應影響，潛水層縱向和橫向彌散系數初始值設置為3.98 m2/h和0.313 m2/h，橫縱彌散度比值為0.079。微承壓層的層縱向和橫向彌散系數初始值設置為27.34 m2/h和2.875 m2/h，橫縱彌散度比值為0.105。區域內潛水含水層滲透系數為9.8 m/d，有效孔隙度為0.3。在人工加速流場下，縱向和橫向彌散系數初始值設置為0.004 3 m2/h和0.001 5 m2/h，縱向和橫向彌散度分別為0.004 6 m和0.001 6 m。區域內承壓，含水層滲透系數為47 m/d，有效孔隙度為0.192。在人工加速流場下，縱向和橫向彌散系數初始值設置為0.5 m2/h和0.015 m2/h，縱向和橫向彌散度分別為0.212 8 m和0.006 4 m。

前述可知11 a降水系列數據可以計算得知降水均值、變態系數、偏態系數。實驗選取豐平枯指標評判降水情況，分別將降水頻率為25%、50%、75%的年份設置為枯水年、平水年、豐水年，計算結果如表3所示?？菟辍⑵剿?、豐水年對應的降水量分別為506.8，621.3，738.5 mm。因此研究選取2008年作為污染物質的預測年份，通過GMS軟件中的MT3DMS模塊實現對城市污水管網外滲情況下點源連續滲漏對地下水污染情況預測，檢測出研究區域有32處污水管道滲漏，模型預測污染物在5 a和10 a內的遷移過程。人為干擾下，10 a內化學耗氧量最大污染范圍為0.29 km2，最大運移距離是1.88 km；天然狀態下，10 a內化學耗氧量最大污染范圍為0.14 km2，最大運移距離是0.99 km，如圖5所示。

表3 豐水年、平水年、枯水年結果

管道泄露5 a后 管道泄露10 a后

3 結論

此次研究以城市污水管網外滲污染地下水源為例，利用滲水試驗、水質分析等步驟完成區域水文水質分析，建立水流和水質數學模型并設置參數進行水源污染模型識別和預測。源匯項處理得到補給資源總量中大氣降水入滲、側向流入、灌溉回滲量分別為146.28萬，688.12萬，20.17萬m3/a，排泄項中基坑降水、蒸發量、農業灌溉分別為701.74萬，23.11萬，134.46萬m3/a，兩者基本達到平衡狀態，水文地質模型具有較高的識別效果?？菟?、平水年、豐水年對應的年份分別為2015年、2008年、2012年，相應的降水量分別為506.8，621.3，738.5 mm。由于所獲取的數據有限，未對水流數學模型進行參數調優，這在后續研究需要完善。