宋雷震

(淮南聯(lián)合大學(xué) 智能制造學(xué)院，安徽 淮南 232038)

城市化進(jìn)程給人們生活帶來(lái)便利的同時(shí)，也引發(fā)城市污水管網(wǎng)外滲污染問(wèn)題，同時(shí)建筑工程開(kāi)發(fā)過(guò)程中的基坑降水也導(dǎo)致了地面沉降和地下水降落漏斗[1]等問(wèn)題。地下水模擬系統(tǒng)(Groundwater Modeling System，GMS)能夠?qū)Φ叵滤贿M(jìn)行數(shù)值模擬[2]。GMS軟件是在現(xiàn)有地下水模型的基礎(chǔ)上建立的全面綜合性的圖形界面軟件，是目前已經(jīng)被廣泛應(yīng)用于世界范圍內(nèi)的地下水模擬方法[3]，它具備強(qiáng)大的數(shù)據(jù)分析處理功能、方便快捷的操作界面、三維可視化的效果。本文利用GMS軟件對(duì)城市污水管網(wǎng)外滲污染地下水進(jìn)行地下水流和水質(zhì)的數(shù)值模擬分析，并運(yùn)用溶質(zhì)運(yùn)移模型預(yù)測(cè)排管泄漏的情況。

1 管網(wǎng)外滲污染地下水的數(shù)值模型

1.1 水文地質(zhì)條件和地下水環(huán)境

本研究選取某市區(qū)域?qū)Τ鞘形鬯芫W(wǎng)外滲污染地下水情況進(jìn)行數(shù)值模擬和預(yù)測(cè)。該地區(qū)氣候?yàn)闇貛駶?rùn)季風(fēng)氣候，海拔高度約為60 m，相對(duì)高度差為8，全年平均氣溫為6.3 ℃，最低氣溫和最高氣溫分別為-31 ℃和35 ℃，無(wú)霜期和封凍期均約為150 d，全年平均日照時(shí)間約為2 700 h[4]。采集氣象站該區(qū)域2005—2018年的年降水量進(jìn)行處理分析，分別得到偏態(tài)系數(shù)、變異系數(shù)、降水量均值3個(gè)統(tǒng)計(jì)參數(shù)，數(shù)值分別為0.30，0.75，656.68 mm。2008—2018年系列降水頻率曲線(xiàn)如圖1所示。從圖1可以看出，降水頻率范圍為3%～96%且集中分布在15%～85%。當(dāng)降水頻率為75%時(shí)，相應(yīng)年平均降水量為524.7 mm。年降水量和頻率可以擬合為一條光滑的曲線(xiàn)。

圖1 2005—2018年系列降水頻率曲線(xiàn)

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

表1 水文地質(zhì)參數(shù)情況

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

表2 地表水質(zhì)情況

1.2 地下水流和水質(zhì)數(shù)學(xué)模型

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

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

結(jié)合前述分析，研究建立非勻質(zhì)各向同性松散巖類(lèi)孔隙微承壓水和潛水的地下水流模擬數(shù)值模型?？紫稘撍涂紫段⒊袎核臄?shù)學(xué)模型計(jì)算公式分別如式(1)和式(2)所示。

(1)

(2)

式中：K、T分別表示潛水含水層的滲透系數(shù)和承壓含水層導(dǎo)水系數(shù)；μ*和μ表示承壓含水層彈性釋水系數(shù)和潛水含水層給水度；H、H0、h分別表示承壓水水頭、潛水水位、地下水位；B表示含水層底板標(biāo)高，Qd和Qr表示排泄強(qiáng)度和入滲補(bǔ)給強(qiáng)度；Qi表示大井開(kāi)采量；h0和h1表示初始水位和一類(lèi)邊界點(diǎn)的水位；q表示二類(lèi)邊界單寬流量；x，y為坐標(biāo)；D表示計(jì)算區(qū)范圍；Γ1和Γ2分別表示一類(lèi)邊界和二類(lèi)邊界；t表示時(shí)間；n表示開(kāi)采井總數(shù)。

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

Qpr=0.1×α×P×F。

(3)

式中：α和P分別表示降水量和降水入滲補(bǔ)給系數(shù)；F表示一分區(qū)面積減去水面的面積。側(cè)向徑流補(bǔ)給量為：

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

(4)

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

Qe=c·F·E601。

(5)

式中：c表示潛水蒸發(fā)系數(shù)；F表示計(jì)算面積；E601表示水面蒸發(fā)量。在此基礎(chǔ)上，利用MT3DMS模塊構(gòu)建地下水溶質(zhì)運(yùn)移模型，完成對(duì)污染的預(yù)測(cè)，選用的模擬因子為化學(xué)耗氧(O2計(jì))[10]。二維非穩(wěn)定流溶質(zhì)運(yùn)移方式下彌散數(shù)學(xué)模型為：

(6)

求解問(wèn)題采用3rd-orderTVD方法和有限差分析方法。I1區(qū)域預(yù)測(cè)5 a后完成建設(shè)，因此5 a內(nèi)地下水流場(chǎng)仍然會(huì)受到人工因素影響，研究選取2種特征污染物完成管網(wǎng)外滲過(guò)程，設(shè)置模型運(yùn)行10 a。

2 外滲污染地下水?dāng)?shù)值模擬分析

2.1 水質(zhì)模型識(shí)別結(jié)果

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

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

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

漏斗中心和QT05標(biāo)記處的地下水位變化情況分別如圖4(a)和4(b)所示。從2008—2018年整體來(lái)看，2個(gè)地段的地下水位表現(xiàn)出動(dòng)態(tài)平衡的現(xiàn)象，同時(shí)實(shí)際動(dòng)態(tài)情況也和變化幅度接近，因此模型的參數(shù)設(shè)計(jì)較為合理，所建立的數(shù)學(xué)模型和水文地質(zhì)模型可靠性較高。

(a)漏斗中心

2.2 溶質(zhì)運(yùn)移模型結(jié)果

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

前述可知11 a降水系列數(shù)據(jù)可以計(jì)算得知降水均值、變態(tài)系數(shù)、偏態(tài)系數(shù)。實(shí)驗(yàn)選取豐平枯指標(biāo)評(píng)判降水情況，分別將降水頻率為25%、50%、75%的年份設(shè)置為枯水年、平水年、豐水年，計(jì)算結(jié)果如表3所示。枯水年、平水年、豐水年對(duì)應(yīng)的降水量分別為506.8，621.3，738.5 mm。因此研究選取2008年作為污染物質(zhì)的預(yù)測(cè)年份，通過(guò)GMS軟件中的MT3DMS模塊實(shí)現(xiàn)對(duì)城市污水管網(wǎng)外滲情況下點(diǎn)源連續(xù)滲漏對(duì)地下水污染情況預(yù)測(cè)，檢測(cè)出研究區(qū)域有32處污水管道滲漏，模型預(yù)測(cè)污染物在5 a和10 a內(nèi)的遷移過(guò)程。人為干擾下，10 a內(nèi)化學(xué)耗氧量最大污染范圍為0.29 km2，最大運(yùn)移距離是1.88 km；天然狀態(tài)下，10 a內(nèi)化學(xué)耗氧量最大污染范圍為0.14 km2，最大運(yùn)移距離是0.99 km，如圖5所示。

表3 豐水年、平水年、枯水年結(jié)果

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

3 結(jié)論

此次研究以城市污水管網(wǎng)外滲污染地下水源為例，利用滲水試驗(yàn)、水質(zhì)分析等步驟完成區(qū)域水文水質(zhì)分析，建立水流和水質(zhì)數(shù)學(xué)模型并設(shè)置參數(shù)進(jìn)行水源污染模型識(shí)別和預(yù)測(cè)。源匯項(xiàng)處理得到補(bǔ)給資源總量中大氣降水入滲、側(cè)向流入、灌溉回滲量分別為146.28萬(wàn)，688.12萬(wàn)，20.17萬(wàn)m3/a，排泄項(xiàng)中基坑降水、蒸發(fā)量、農(nóng)業(yè)灌溉分別為701.74萬(wàn)，23.11萬(wàn)，134.46萬(wàn)m3/a，兩者基本達(dá)到平衡狀態(tài)，水文地質(zhì)模型具有較高的識(shí)別效果?？菟辍⑵剿?、豐水年對(duì)應(yīng)的年份分別為2015年、2008年、2012年，相應(yīng)的降水量分別為506.8，621.3，738.5 mm。由于所獲取的數(shù)據(jù)有限，未對(duì)水流數(shù)學(xué)模型進(jìn)行參數(shù)調(diào)優(yōu)，這在后續(xù)研究需要完善。