畢桂真 孟令杰 韋愛菊

（山東省菏澤市水文局 菏澤 274008）

1 地下水流數(shù)值模型的建立及求解

1.1 水文地質(zhì)概念模型

從含水層結(jié)構(gòu)上看，研究區(qū)為多層結(jié)構(gòu)區(qū)。據(jù)實(shí)測資料，園區(qū)內(nèi)潛水埋深一般為1.51～4.73m，含水層厚度為14.6m，潛水層之間夾有平均厚度在5～10m的粘質(zhì)砂土、粉質(zhì)粘土弱透水層。此次研究把潛水作為計(jì)算目的層，通過分析鉆孔資料和水文地質(zhì)調(diào)查資料，將含水層概化為非均質(zhì)、各向同性含水層。所以模擬區(qū)地下水的水動(dòng)力條件概化為穩(wěn)定的三維流。

1.2 地下水流數(shù)學(xué)模型

根據(jù)前述的水文地質(zhì)概念模型，研究區(qū)地下水運(yùn)動(dòng)數(shù)學(xué)模型如下：

式中：H—地下水位標(biāo)高（m）；H0—第一類邊界水位標(biāo)高（m）；t—時(shí)間（t）；A—模擬區(qū)四周邊界；K—滲透系數(shù)（m/d）；Z—潛水面水頭分布值（m）；x，y，z—坐標(biāo)變量（m）；rw—抽水井的半徑（m）；SS—含水層的釋水率（L-1）；Qi—第 i口井的抽水流量（m3/d）；W—研究區(qū)三維流的源匯項(xiàng)（L/d）；μ—給水度；U—滲流研究區(qū)域；f（x，y，z）—初始水位分布（m）；（r，θ）—輔助柱坐標(biāo)變量。

1.3 模型離散化及基礎(chǔ)資料的給定

1.3.1 空間的離散化

研究范圍是一個(gè)規(guī)則的矩形區(qū)域。將研究區(qū)垂向上剖分為4層，其中潛水含水層1層，弱透水層1層，微承壓含水層1層和隔水層1層，由于研究區(qū)地處平原，地形平坦，面積較小，地形高程和各層厚度以平面的形式輸入到模型中。

1.3.2 時(shí)間的離散化

選取2017年4月1日為模擬的起始時(shí)間，7300日后為模擬的終止時(shí)間。最大時(shí)間步長為200d，嚴(yán)格控制每次迭代的誤差。在每個(gè)應(yīng)力期保持含水層補(bǔ)給和排泄強(qiáng)度不變。

1.3.3 初始條件的輸入

初始水位采用2017年4月統(tǒng)測的分層觀測水位，將實(shí)測水位輸入到模型中從而得到研究區(qū)數(shù)值模擬的初始流場，其中弱透水層的初始水位取上下兩個(gè)含水層水位的算術(shù)平均值。

1.3.4 邊界條件的輸入

四周邊界概化為第一類邊界條件，每月1日水位值作為水位觀測值輸入到Visual MODFLOW中。

1.3.5 水文地質(zhì)參數(shù)的初值

根據(jù)2017年3月16日穩(wěn)定流抽水試驗(yàn)資料和水文地質(zhì)勘察室內(nèi)土工試驗(yàn)資料，結(jié)合研究區(qū)地質(zhì)、水文地質(zhì)條件，將以上資料所得參數(shù)輸入到地下水流模型中。

1.3.6 源匯項(xiàng)的處理

潛水的補(bǔ)給來源包括降雨入滲補(bǔ)給、灌溉入滲補(bǔ)給和渠系滲漏補(bǔ)給，其中降水入滲補(bǔ)給占主導(dǎo)地位。地下水排泄方式包括側(cè)向徑流排泄、潛水蒸發(fā)排泄和人工開采。補(bǔ)給和排泄按補(bǔ)給強(qiáng)度和開采強(qiáng)度處理，人工開采按雙井開采流量計(jì)算。降水采用鄄城縣多年平均降水量609.7mm，蒸發(fā)強(qiáng)度采用魏樓閘水文站多年平均水面蒸發(fā)量，潛水蒸發(fā)極限埋深參考鄧集試驗(yàn)站成果取4.0m。

1.4 數(shù)學(xué)模型的識(shí)別、驗(yàn)證

1.4.1 模型識(shí)別

把各種水文地質(zhì)資料代入模型，以迭代殘差和最大水頭改變量最小為目標(biāo)。通過調(diào)整分區(qū)參數(shù)值使二者盡量小，并據(jù)此來判斷所用水文地質(zhì)參數(shù)及分區(qū)是否合理。經(jīng)反復(fù)調(diào)整參數(shù)，獲得了較為滿意的水文地質(zhì)參數(shù)。

1.4.2 模型驗(yàn)證

選用2016年（年降水量611.0mm）作為驗(yàn)證年份，以周邊11#、15B#、32A#多年地下水動(dòng)態(tài)監(jiān)測井作為參照井進(jìn)行驗(yàn)證。經(jīng)驗(yàn)證，各監(jiān)測井實(shí)測水位與計(jì)算水位差值的絕對(duì)值小于0.05m，表明各觀測孔的水位計(jì)算結(jié)果與實(shí)測結(jié)果吻合很好，驗(yàn)證了所取參數(shù)的合理性。同時(shí)對(duì)比2017年4月潛水的實(shí)測水位與模型計(jì)算水位，模擬流場與實(shí)際流場的變化趨勢(shì)基本一致，可以運(yùn)用到地下水水質(zhì)模型中。

1.5 地下水流數(shù)學(xué)模型預(yù)測

在考慮各項(xiàng)補(bǔ)給項(xiàng)的條件下，使用建立的地下水流數(shù)學(xué)模型，假定園區(qū)工業(yè)取用水戶雙井穩(wěn)定流開采地下水，單井出水量1200m3/d。連續(xù)抽水48h后實(shí)際觀測記錄與預(yù)測結(jié)果見表1。

2 地下水水質(zhì)模型及預(yù)測

通過Visual MODFLOW中的MT3DMS模塊計(jì)算污染物質(zhì)的運(yùn)移情況，可以求出污染物在地下水系統(tǒng)中的變化規(guī)律，預(yù)測研究區(qū)污染物質(zhì)在不同時(shí)刻、不同的情況下所導(dǎo)致的地下水污染程度。

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

根據(jù)研究區(qū)的具體條件，采用下述的溶質(zhì)運(yùn)移模型：

式中：C—溶解于水中的污染物濃度；n—孔隙度；xi—空間坐標(biāo)；qw—源（正值）或匯（負(fù)值）的單位流量；Dij—水動(dòng)力彌散系數(shù)張量；t—時(shí)間；F—固相表面的溶質(zhì)濃度；C0—源匯項(xiàng)的濃度；Ω—空間區(qū)域；Vi—地下水滲透流速；（x，y，z）—空間位置；r—研究區(qū)邊界。

2.2 水質(zhì)模型的建立

水質(zhì)模型以水流模型為基礎(chǔ)建立，水質(zhì)模型的概化與所建立的水流概念模型相符。水質(zhì)模擬區(qū)范圍、含水層結(jié)構(gòu)、邊界類型劃分、源匯項(xiàng)的概化均與水流概念模型相同，流體概化為不可壓縮的均質(zhì)流體，粘度和密度均為常數(shù)。

2.2.1 模擬因子的選擇

此次研究為穩(wěn)定流情況下預(yù)測研究區(qū)污染物的運(yùn)移情況，考慮污染物質(zhì)在含水層中的線性等溫吸附（平衡），無反應(yīng)項(xiàng)，選擇化學(xué)需氧量（COD）為示蹤劑。

2.2.2 邊界條件、初始條件

研究區(qū)水質(zhì)模型選定2017年4月1日作為初始時(shí)刻，初始時(shí)刻含水層中污染物的濃度為0。

表1 園區(qū)抽水試驗(yàn)觀測記錄與預(yù)測結(jié)果對(duì)照表

2.2.3 模型參數(shù)

溶質(zhì)運(yùn)移模型涉及的參數(shù)中含水介質(zhì)的有效孔隙度（n）由試驗(yàn)所得n=0.19，其他參數(shù)的取值如下：

（1）地下水滲流速度：按照2017年3月16日實(shí)際抽水試驗(yàn)結(jié)果，含水層地下水垂向滲透系數(shù)速度為1.0521×10-4cm/s，由2017年3月26日室內(nèi)土工試驗(yàn)，弱透水層地下水垂向滲透系數(shù)確定為2.73×10-6cm/s，隔水層為9.94×10-8cm/s。

（2）縱向彌散系數(shù)的確定：溶質(zhì)運(yùn)移模型需要的彌散系數(shù)是在結(jié)合抽水試驗(yàn)進(jìn)行的野外彌散實(shí)驗(yàn)，同時(shí)取研究區(qū)的土樣，采用一維土柱彌散實(shí)驗(yàn)法，進(jìn)行室內(nèi)土柱彌散實(shí)驗(yàn)而確定的。經(jīng)過計(jì)算，評(píng)價(jià)區(qū)縱向彌散系數(shù)為1.46m/d。

（3）其他參數(shù)選用經(jīng)驗(yàn)值。

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

假定園區(qū)發(fā)生污染事故，污染物下滲污染當(dāng)?shù)氐叵滤?。以地面表層受到一定濃度COD污染作為初始條件。在假定的園區(qū)工業(yè)取用水戶雙井穩(wěn)定流開采地下水的條件下，通過Visual MODFLOW中的MT3DMS模塊和已建立的溶質(zhì)運(yùn)移模型，預(yù)測不同時(shí)間段內(nèi)污染物垂直和水平運(yùn)移的情況。設(shè)定污染羽外圍污染物的濃度達(dá)到0.5mg/L時(shí)為該污染物的影響范圍。

3.1 不同時(shí)段污染物的垂直運(yùn)移情況

以地面表層受到濃度為5000mg/L的COD污染為初始條件，經(jīng)過365d（1年）的垂直運(yùn)移，地下水已經(jīng)受到了輕度的污染。此時(shí)，在潛水的表層有仍有少許積累外，污染物已經(jīng)到達(dá)潛水含水層的中下部，中心位置的濃度達(dá)到100mg/L。經(jīng)過5475d（15年）的垂直運(yùn)移，污染物在垂向上已經(jīng)能夠到達(dá)隔水層的下部，實(shí)現(xiàn)了對(duì)隔水層的滲透性運(yùn)移。運(yùn)移7300d（20年）后，垂向上的影響范圍沒有顯著的增加，但潛水表層的污染物濃度明顯增加，污染物在橫向和縱向上得到了積累，地下水中的污染物運(yùn)移范圍進(jìn)一步擴(kuò)大，濃度遞增顯著，中心位置處不但濃度高達(dá)1000mg/L，而且范圍擴(kuò)大了近10倍。

3.2 不同時(shí)段污染物的水平運(yùn)移情況

為對(duì)比不同濃度下污染物的水平運(yùn)移情況，分別對(duì)初始濃度5000mg/L、500mg/L的COD做污染物的運(yùn)移模擬，對(duì)比在7300d后的運(yùn)移等值線，5000mg/L濃度COD的10mg/L等值線運(yùn)移了582m，污染羽（濃度為0.5mg/L）運(yùn)移了680m；500mg/L濃度COD的10mg/L等值線運(yùn)移了426m，污染羽運(yùn)移了606m。說明污染物濃度越高，污染羽在相同時(shí)段運(yùn)移距離越遠(yuǎn)速度越快。

4 結(jié)論與建議

根據(jù)鄄城縣工業(yè)園區(qū)水文地質(zhì)條件的調(diào)查，運(yùn)用Visual MODFLOW模型對(duì)工業(yè)園區(qū)受到污染物污染后污染物垂直和水平運(yùn)移的模擬計(jì)算，得出以下結(jié)論：

（1）當(dāng)園區(qū)受到污染物污染后，隨著時(shí)間的推移，園區(qū)地下水從受到輕度污染逐步演變到重度污染。從初期污染物在垂向上到達(dá)隔水層的下部并滲透性運(yùn)移過隔水層，到后期地下水中的污染物運(yùn)移范圍進(jìn)一步擴(kuò)大，濃度遞增顯著，污染物濃度逐年線性增加，若不及時(shí)治理，污染范圍會(huì)進(jìn)一步擴(kuò)大。隨著抽水時(shí)間的延長，地下水埋深逐年增加，地下水位逐年降低，污染羽的運(yùn)移速度在第二含水層顯著增加。

（2）園區(qū)場地基礎(chǔ)之下第一巖土層不能滿足天然防滲層的要求，園區(qū)內(nèi)事故污染下滲會(huì)造成當(dāng)?shù)氐叵滤奈廴?。鑒于園區(qū)地層不能滿足最低滲透性標(biāo)準(zhǔn)，建議采用在園區(qū)四周和底部鋪設(shè)土工膜的方法，進(jìn)行人工水平防滲和垂直防滲，防止由于園區(qū)受到污染物污染對(duì)地下水水質(zhì)造成影響