李白業(yè)，李雪妮

(新疆地礦局第一水文工程地質(zhì)大隊(duì)，新疆烏魯木齊830091)

近幾年來庫車縣工業(yè)迅速發(fā)展，當(dāng)?shù)亟?jīng)濟(jì)持續(xù)增長，地表水與地下水的需水量大幅度增加，廢污水的排放及處理面臨巨大的壓力。要控制地下水的污染，保護(hù)地下水資源，就必須了解地下水現(xiàn)狀污染的程度，預(yù)測將來可能影響的范圍。在現(xiàn)狀污染調(diào)查與評價基礎(chǔ)上對研究區(qū)典型地段地下水中的特定污染組分濃度隨空間與時間的變化進(jìn)行計算預(yù)測，以掌握地下水污染的空間變化規(guī)律及發(fā)展趨勢，以備及時有效地采取防控措施，有預(yù)見地減少污染危害。

1 預(yù)測方法的選用

地下水污染物遷移(見圖1)模型研究是在現(xiàn)狀污染調(diào)查與評價基礎(chǔ)上對研究區(qū)典型地段地下水中的特定污染組分濃度隨空間與時間的變化進(jìn)行計算預(yù)測，以掌握地下水污染的空間變化規(guī)律及發(fā)展趨勢，以備及時有效地采取防控措施，有預(yù)見地減少污染危害。本次計算典型地段選擇新老城區(qū)及工業(yè)園區(qū)兩個區(qū)段;預(yù)測單源污染液連續(xù)滲入在飽水帶形成的污染暈(也稱為彌散帶范圍)發(fā)展趨勢及在未來30年對下游地段地下水的影響。

圖1 地下水中污染物運(yùn)移示意圖

通常，地下水污染的預(yù)報分析可細(xì)分為定性分析、半定量分析及定量分析三種類型，以后者的分析精度最高，但分析精度越高對資料序列的要求也相應(yīng)越高，野外投入也越大。水質(zhì)污染定量預(yù)報模型可分為黑箱模型和灰箱模型?；蚁淠Ｐ鸵缘叵滤畡恿W(xué)方法為主，方法主要包括簡單的解析模型、半解析模型、數(shù)值模型。其中，數(shù)值模型法要求有長序列的特定組分監(jiān)測資料和足夠的野外彌散試驗(yàn)參數(shù)值，加上計算工作量繁重，多在環(huán)評一級或復(fù)雜條件的環(huán)評二級中使用。

地下水流態(tài)的非均質(zhì)性，約束了常規(guī)解析法的應(yīng)用范疇，盲目使用會產(chǎn)生很大誤差[1]。計算機(jī)技術(shù)的廣泛應(yīng)用為解析法的改進(jìn)提供了有力的條件，借鑒地下水水量數(shù)值模型中潛水的非線性問題(導(dǎo)水系數(shù)T隨含水層厚度變化，呈現(xiàn)動態(tài)參數(shù)特征)的迭代處理辦法，在解析計算中增加一個外部迭代過程，采用追趕法計算地下水非均質(zhì)參數(shù)(稱動態(tài)參數(shù)，如滲透系數(shù) k=f[kx，x])，使參數(shù)域無限逼近污染暈鋒面，這種方法可稱為改進(jìn)的解析模型方法。

綜上，本次工作采用:定量分析方法→灰箱模型→地下水水動力學(xué)方法→改進(jìn)解析法→適用于非均質(zhì)條件的:改進(jìn)的一維解析模型方法。

2 數(shù)學(xué)方程

有限寬度、無限長度，存在源匯項(xiàng)及綜合作用，無邊界的一維均勻穩(wěn)定流孔隙滲流溶質(zhì)運(yùn)移偏微分方程或泛定方程式[2]通常可以表示為:

式中:C為溶質(zhì)濃度，[ML－3];v為在滲流方向的平均地下水孔隙流速，[LT－1];D為縱向的彌散系數(shù)張量，[L2T－1];C'為源匯項(xiàng)中的溶質(zhì)濃度，[ML－3];W'為源匯項(xiàng)流量，[M3T－1];n為有效孔隙度，[無量綱];K為水力傳導(dǎo)系數(shù)張量，[無量綱];Rk為離子交換反應(yīng)、分解反應(yīng)、化合反應(yīng)等的溶質(zhì)生成量，[ML－3];x為笛卡爾坐標(biāo)中滲流方向，[M];i，m為溶質(zhì)綜合反應(yīng)的序次及總數(shù)。

公式(1)中的第一項(xiàng)代表水動力彌散分量，第二項(xiàng)代表對流分量，第三部分多代表(源)匯項(xiàng)作用分量，第四項(xiàng)代表其它綜合物理化學(xué)(降解)作用量如:離子交換、分解作用、沉淀作用及放射性衰變等分量。

以上方程中僅考慮污染物直接作用于飽水帶的情況，實(shí)際上污染物在通過包氣帶時，由于吸附等溫作用影響使溶質(zhì)的遷移與飽水帶水的滲流速度產(chǎn)生差異，出現(xiàn)遲滯現(xiàn)象。根據(jù)Grove[1976]的意見，考慮平衡控制的離子交換反應(yīng)，將上述通用表達(dá)式變形后得到考慮吸附等溫影響滯后作用，存在其它作用(包括分解，化合及放射性衰變等)但無源匯項(xiàng)時的一維流水質(zhì)運(yùn)移偏微方程式為:

其中:R為阻滯系數(shù)，[無量綱];μ為綜合作用一階系數(shù)，當(dāng)存在某種作用時反演求得;γ(x)為綜合作用零階系數(shù)，與空間坐標(biāo)有關(guān)，亦在反演時求得。

3 建模步驟與解算方法

3.1 建模步驟

建模步驟見右側(cè)的流程框圖(見圖2)。

圖2 解析模型法建模步驟流程框圖

3.2 解算方法

本次工作經(jīng)篩選采用中美環(huán)保科技交流推薦的軟件包STANMOD 2.1(STudio of ANalytical MODels，1999 年)，這是一個WINDOWS－XP圖形界面軟件包，整合了包括CXTFIT[Toride et al.，1995]在內(nèi)不同解法的幾個一維流及二維流解算模型，本次計算選用的CXTFIT模型算法[3]是由美國鹽漬土實(shí)驗(yàn)室研制的，逆計算基于列文伯格－馬夸爾特法(Levenberg－Marquardt算法[4]，簡稱LM算法)對均質(zhì)一維穩(wěn)定流動中溶質(zhì)易混合置換實(shí)驗(yàn)所得到的穿透曲線(Breakthrough curve，BTC)進(jìn)行非線性最小二乘擬合對偏微分方程進(jìn)行解算的。該算法是高斯－牛頓法和梯度下降法的結(jié)合，既有高斯牛頓法的快速、也有梯度下降法的全局搜索特性[5]。

4 模型參數(shù)

模型參數(shù)包括實(shí)測參數(shù)和經(jīng)驗(yàn)參數(shù)兩部分。主要實(shí)測參數(shù)包括:污染源下游帶水流方向的含水層參數(shù)(K－平均滲透系數(shù)、n－平均有效孔隙度、I－平均水力坡度)，污染源示蹤組分的遷移速率Vc及平均濃度C0;經(jīng)驗(yàn)參數(shù)主要包括:DL－縱向彌散度，R－阻滯系數(shù)。

其中，地下水參數(shù)宜采用動態(tài)平均加權(quán)值，權(quán)重采用污染暈過境參數(shù)區(qū)的值徑向長度(注:二維流一般采用的是參數(shù)分區(qū)面積)。參數(shù)分區(qū)則采用本次水量數(shù)值模型分區(qū)成果。

地下水動態(tài)平均參數(shù)計算公式為:

4.1 平均滲透系數(shù)(K)

平均滲透系數(shù)為動態(tài)地下水實(shí)測參數(shù)，參數(shù)分區(qū)及參數(shù)值采用本次水量數(shù)值模型的參數(shù)分區(qū)成果;權(quán)值采用平水位期的地下水等水位線圖進(jìn)行圖解量算。根據(jù)公式(3)，參數(shù)計算結(jié)果見表1。

表1 平均滲透系數(shù)K(m/d)計算一覽表

4.2 平均水力坡度(I)

在地下水等水位線圖上，對徑向的水力坡度進(jìn)行分段圖量計算，徑向水力坡度計算的起點(diǎn)為污染源處，終點(diǎn)為污染暈鋒面位置。計算公式為I=△H/L，計算結(jié)果見表2。

4.3 有效孔隙度(n)

在沒有進(jìn)行彌散試驗(yàn)時，有效孔隙度主要用來將參與地下水實(shí)際流速的計算，由于地下水孔隙度變化多依賴于含水層巖性，變化幅度也較K值小的多，因此，采用在水量數(shù)值模型計算中其分區(qū)范圍與數(shù)值。延用動態(tài)滲透系數(shù)計算辦法，根據(jù)公式(3)計算的各平均有效孔隙度值見表3。

表2 平均水力坡度I計算一覽表

表3 平均有效孔隙度n計算一覽表

4.4 縱向彌散度(DL)

根據(jù)中國地質(zhì)大學(xué)教授馬騰等的論著[6]，世界上越來越多的室內(nèi)外彌散試驗(yàn)不斷地證實(shí)了空隙介質(zhì)中水動力彌散具有尺度效應(yīng)。一些研究者綜合世界范圍內(nèi)百余個水質(zhì)模型中所使用的縱向彌散度與空間尺度的關(guān)系，認(rèn)為水動力彌散的尺度效應(yīng)具有分形特征，同時，還給出了不同模型縱向彌散度尺度效應(yīng)的分維數(shù)。根據(jù)分維數(shù)和空間尺度可以按經(jīng)驗(yàn)公式計算出適合的彌散度參數(shù)初值(見表4、表5)。

表5 不同巖性條件下的彌散度初值aL計算結(jié)果表

4.5 阻滯系數(shù)(R)

該參數(shù)反映了污染遷移速度與地下水流動速度之間的差異[7][8]，不同污染因子在不同巖性環(huán)境條件下的遷移速度有很大的差別，但當(dāng)污染因子與污染物整體運(yùn)移差異不大時，這種滯后作用就不明顯了。一般試驗(yàn)測定[8]的粘土類的阻滯系數(shù)為2，～5砂土類多在1.0～1.2之間，砂礫石類就更小了，因此本次工作采用1.0～1.5的經(jīng)驗(yàn)值范圍作為阻滯系數(shù)初值，如有可能就在模型參數(shù)識別時根據(jù)現(xiàn)狀污染趨勢進(jìn)行適當(dāng)?shù)男拚?/p>

5 計算結(jié)果

預(yù)報的起始位置在依巴格－烏恰鎮(zhèn)的鐵路沿線附近，這里也是城市污水處理場所在的位置。預(yù)報模型的時間以5年為間隔，共預(yù)報30年。根據(jù)曲線判別出每年污染暈移動的位置。新老城區(qū)30年末污染暈峰面到達(dá)下游8.0 km外的牛場－喀魯烏古一帶南3km農(nóng)業(yè)區(qū)，前峰帶寬度近400 m，累計污染物容量為8 026.17 C0×B×M(B為污染帶寬度，M為污染含水層厚度)。工業(yè)園區(qū)污染峰面30年末到達(dá)5.7 km的西水東調(diào)干渠一線，前峰帶寬度近300 m，5 154.84 C0×B×M。

6 結(jié)語

本次工作利用國際標(biāo)準(zhǔn)化一維流水質(zhì)運(yùn)移軟件包，采用改進(jìn)的解析法對下游地下水污染遷移問題進(jìn)行研究，在新疆地下水研究中尚屬首次。改進(jìn)的一維流解析法提高了計算精度，完全適用于非均質(zhì)流條件，選用方法符合環(huán)評二級要求。因監(jiān)測資料缺乏，模型識別略顯粗糙，影響污染暈彌散范圍的預(yù)測精度，但在以對流為主的污染遷移中，對污染暈向下游的運(yùn)移距離影響誤差可以接受。計算結(jié)果僅可用于下游污染影響帶地下水開采規(guī)劃決策參考，不能作為污染治理專項(xiàng)設(shè)計及水事糾紛處理的法律依據(jù)。

[1]薛禹群，謝春紅.地下水?dāng)?shù)值模擬.北京:科學(xué)出版社.2007.

[2]I.Javandel，C.Doughty，等著.林學(xué)鈺，李生彩，等譯.地下水運(yùn)移數(shù)學(xué)模型手冊.長春:吉林科學(xué)技術(shù)出版社.1985.

[3]N.Toride，F(xiàn).J.Leij，and.The CXTFIT code for Estimating Transport Parameters from Laboratory or Field Tracer(Version 2.1).U.S.:SALINITY LABORATORY AGRICULTURAL RESEARCH SERVICE AND DEP ARTMENT OF AGRICULTUGE RIVERSIDE CALIFORNIA.1999.

[4]Marquardt，D.W.，An algorithm for least－ squares estimation of nonlinear parameters ，J.soc.ind.appl.Math.，11，431 －441，1963.

[5]袁亞湘，孫文瑜.最優(yōu)化理論與方法.北京:科學(xué)出版社.1997.

[6]馬騰，王焰新.“地下水污染與防治”課程建設(shè)的思考與實(shí)踐.大學(xué)環(huán)境類課程報告論壇論組委會編.大學(xué)環(huán)境類課程報告論壇論文集(2006).北京:高等教育出版社，163－166.

[7]石振華，李傳堯，等.城市地下水工程與管理手冊.北京:中國建筑工業(yè)出版社.1993.

[8]Bruce E.Rittmann，Perry L.Mccarty，等著.文湘華，王建龍譯.環(huán)境生物技術(shù)原理與應(yīng)用.北京:清華大學(xué)出版社.2004.