一種用于突發(fā)性場地污染模擬的解析模型

陳家軍，王 浩，全向春，文 軍

（1.北京師范大學(xué)水沙科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100875；2.北京師范大學(xué)環(huán)境學(xué)院，北京 100875；3.新疆水利水電勘測設(shè)計(jì)院，新疆烏魯木齊 830091）

0 引 言

隨著化學(xué)品使用范圍日益擴(kuò)大，在其生產(chǎn)、使用、銷售和儲運(yùn)過程中發(fā)生泄漏、爆炸和火災(zāi)等突發(fā)性環(huán)境污染事故也相應(yīng)增多。這類事故發(fā)生突然，破壞性大，直接關(guān)系到人民群眾生命財產(chǎn)的安全，同時也會使人類賴以生存的生態(tài)環(huán)境遭到破壞。做好突發(fā)性環(huán)境污染事故的預(yù)防、監(jiān)測、應(yīng)急處置及善后工作，建立運(yùn)行有效、行動快速的事故應(yīng)急管理系統(tǒng)是最大限度減少損失的關(guān)鍵［1］。突發(fā)性污染造成的污染場地是潛在的污染源，一旦發(fā)生污染事故，不僅后果嚴(yán)重，往往會造成惡劣的社會影響［2］。如美國Love運(yùn)河污染事件是一起典型的污染場地引起的公害事件，致使大量人群的生命健康受到損害，并引起巨大經(jīng)濟(jì)損失［3］。隨著中國經(jīng)濟(jì)的快速發(fā)展與經(jīng)濟(jì)結(jié)構(gòu)的調(diào)整，污染場地問題備受關(guān)注。2004年4月，北京地鐵5號線工程宋家莊附近挖掘出含六六六粉、滴滴涕的污泥超過10 000t，后查明此地是北京農(nóng)藥廠原址［4］。

為了更有效地處理各種突發(fā)性環(huán)境污染事故，減少事故對環(huán)境的損害，確保環(huán)境盡可能地安全，有必要開發(fā)適用于污染場地污染物遷移轉(zhuǎn)化的數(shù)學(xué)模型?；谀Ｐ皖A(yù)測分析一定時間內(nèi)污染物在包氣帶及地下水中的遷移動向，開展相應(yīng)技術(shù)評估，采取合理的管理、應(yīng)對措施并選取適宜的場地污染修復(fù)技術(shù)，為應(yīng)急監(jiān)測與處置工作提供決策支持，確保突發(fā)性污染事故的應(yīng)急監(jiān)測與處置工作快速有效地進(jìn)行。目前，國內(nèi)外已經(jīng)針對場地污染建立了一些污染物運(yùn)移的數(shù)學(xué)模型，其中數(shù)值模型［5－8］適用條件廣，更為有效，但在應(yīng)用中對水動力彌散參數(shù)和幾何數(shù)據(jù)等要求極高（污染場地范圍需有很高的研究程度才能達(dá)到要求），操作較為復(fù)雜，從而限制了其實(shí)際應(yīng)用。解析模型雖然適用的場地條件不如數(shù)值模型廣，但所要求的參數(shù)較少，模擬計(jì)算工作量小，操作使用簡單易行，更加符合突發(fā)性場地污染模擬預(yù)測的要求［9－11］。Beltman等研究了土壤吸附、微生物降解的空間變異性對農(nóng)藥滲漏量的干擾，建立相應(yīng)的數(shù)學(xué)模型并預(yù)測出這兩種因素對農(nóng)藥遷移、降解的影響規(guī)律；此外，該模型還考慮了植被對污染物的被動吸附作用［9］。Van Der Zee等根據(jù)前人研究建立了飽和帶和非飽和帶的農(nóng)藥遷移解析模型［12］。但這些研究大多針對難揮發(fā)性污染物，模型建立過程中一般都忽略了污染物的揮發(fā)性，因此在應(yīng)用時具有一定的局限性。另外，在突發(fā)性污染事故發(fā)生初期，若沒有遇到降雨，污染物在重力作用下會向地下自由入滲，這種情況也應(yīng)納入場地突發(fā)性污染模擬的范疇。中國目前針對此類突發(fā)性場地污染的解析模型研究還比較匱乏，相關(guān)的文獻(xiàn)資料非常少見。筆者將建立能夠應(yīng)用于一般性污染場地適用于各種污染物，并且方便使用的包氣帶、飽和帶污染物遷移數(shù)學(xué)模型并得到相應(yīng)的解析解。

1 模型假定及概化

1.1 假設(shè)條件

若沒有降雨發(fā)生，污染物發(fā)生泄漏前后主要通過自由入滲方式向地下遷移，可采用自由入滲模型模擬污染物的垂向遷移。自由入滲模型中，污染物主要借助重力作用向地下遷移，由于此模型中污染物本身構(gòu)成了流動相，遷移過程中忽略分子擴(kuò)散作用，選用一維垂直入滲模型進(jìn)行模擬。若污染發(fā)生時有降雨發(fā)生，則污染物隨雨水向下遷移。降雨入滲模型中，污染物在包氣帶中主要借助重力作用向下遷移，遷移過程中伴隨植物根系吸收、污染物揮發(fā)、分子擴(kuò)散等多種機(jī)制共同作用，因此污染物在包氣帶中的遷移可使用剖面二維非飽和淋溶模型進(jìn)行模擬，水平方向僅存在分子擴(kuò)散作用。污染物進(jìn)入地下水后隨地下水流運(yùn)移，選用平面二維飽和溶質(zhì)運(yùn)移模型進(jìn)行模擬，除根系吸收外，主要考慮機(jī)制與包氣帶相同。另外，該模型還包括以下假設(shè)條件：①污染物在包氣帶中的運(yùn)動視為重力作用下的遷移，地表坡度為0；②自由入滲中，僅考慮污染物在包氣帶的運(yùn)移，忽略植物吸收作用及生物降解作用，污染物在土壤中的最大飽和度為土壤的有效孔隙度；③若有降雨發(fā)生，降雨量平均分配至每天，分子擴(kuò)散作用在所有區(qū)域都發(fā)生；④不同包氣帶位置的土壤質(zhì)地、水動力參數(shù)以及有機(jī)物成分與種類均相同；⑤污染物在土壤中的吸附為線性吸附；⑥根系對污染物的吸附為被動式吸收，選用根系被動吸附模型進(jìn)行模擬，揮發(fā)作用只通過孔隙氣體進(jìn)行。

1.2 模型概化

現(xiàn)實(shí)中污染源的污染物含量分布不均，污染源也以不規(guī)則形狀存在。本研究旨在建立一種能夠適用各種污染場地模型，因此假定污染源都位于地表位置，形狀為一矩形區(qū)域，厚度較?。ü什豢紤]污染源的厚度）。污染物在污染區(qū)域內(nèi)平均分布，各位置污染物含量（質(zhì)量分?jǐn)?shù)）相同，污染物含量可由污染物泄漏總量計(jì)算得到

式中：C為污染物平均含量；M為污染物總量；l1×l2為地面污染源的初始尺寸，l1為初始長度，l2為初始寬度；h為污染物在土壤中的厚度；ρs為土壤密度。

自由入滲模型對自由入滲的概化較為簡單，僅考慮垂直一維重力條件下的污染物遷移，以便于污染發(fā)生后迅速對污染物的入滲情況進(jìn)行模擬和預(yù)測。

圖1 包氣帶剖面污染物遷移概化Fig.1 Generalization of sectional contaminant migration in unsaturated zone

降雨入滲包氣帶以及飽和帶的模型概化見圖1，其中，O為原點(diǎn)，z軸為垂直向下方向。污染物隨降雨入滲進(jìn)入地下土壤，污染物遷移過程中存在于水相，部分污染物被根系吸收進(jìn)入植物體內(nèi)，另一部分在重力以及分子擴(kuò)散條件下發(fā)生遷移，遷移過程中不斷伴隨揮發(fā)、吸附作用及微生物降解。殘留的污染物繼續(xù)向地下遷移進(jìn)入飽和帶（圖2），污染物擴(kuò)散至飽和帶時初始的污染源尺寸擴(kuò)大為l′1×l′2矩形區(qū)域（l′1和l′2分別為污染羽到達(dá)飽和帶的長度和寬度），包氣帶與飽和帶界面的污染物濃度作為飽和帶的污染源。污染物在飽和帶（含水層）的遷移概化見圖2，假定含水層厚度較小，飽和帶中重點(diǎn)考慮污染物在平面的遷移和擴(kuò)散，主要考慮的機(jī)制有介質(zhì)吸附、微生物降解、對流彌散作用，擴(kuò)散作用中包氣帶只考慮分子擴(kuò)散，飽和帶水流方向取機(jī)械彌散，橫向則只考慮分子擴(kuò)散。

圖2 飽和帶剖面和平面污染物遷移概化Fig.2 Generalization of sectional and planar contaminant migration in saturated zone

2 自由入滲模型

若污染事故發(fā)生后短時間內(nèi)無降雨發(fā)生，污染物在重力作用下會以自由態(tài)形式緩慢滲透進(jìn)入土壤，此時可將污染物在土壤中的遷移簡化為單相流體的流動，選用流體一維垂向流模型模擬污染物的流動，其基本控制方程為

式中：θc為污染物在土壤中的飽和度；Kc（θc）為污染物相的非飽和滲透系數(shù)；Dc（θc）為污染物擴(kuò)散率；t為擴(kuò)散時間。

定解條件為

式中：θs為污染物在土壤中的最大飽和度。式（3）的解析解為式中：θc（z，t）為z方向t時刻的污染物飽和度；K－c為土壤滲透系數(shù)；D－c為污染物擴(kuò)散系數(shù)。式（4）即為無降雨時突發(fā)事故場地污染的自由入滲解析解。

3 降雨入滲模型

3.1 水分運(yùn)動方程

由降雨強(qiáng)度及根系吸附作用得到土壤包氣帶孔隙水實(shí)際平均流速υu為

式中：ξ為降雨入滲系數(shù)；Jw為降雨入滲速度；CET為植物的蒸騰系數(shù)；θ為包氣帶平均體積含水率；qu為包氣帶滲流速度。

根據(jù)均質(zhì)各向同性條件得到飽和帶孔隙水實(shí)際平均流速υs為

式中：H2、H1為2號和1號地下水井的水頭；L為兩口水井間在地下水流向上的距離；Ks為飽和帶的飽和滲透系數(shù)；qs為飽和帶滲流速度；n為飽和帶孔隙度。

3.2 溶質(zhì)運(yùn)移方程

考慮生物降解以及土壤對污染物的線性吸附機(jī)制，液相污染物的遷移可表示為

式中：C1為土壤污染物含量；ku為污染物在包氣帶的生物降解系數(shù)；Ru為包氣帶阻滯因子；Dzu和Dxu分別為包氣帶中z、x方向水動力彌散系數(shù)。

3.2.1 包氣帶模型

（1）根系吸附。若污染場地存在植被覆蓋，根系對污染物的吸收速率Rr可以表示為［10］

式中：F為蒸騰流濃度因子；S為植物吸收水分的速率常數(shù)。

這樣根系區(qū)的溶質(zhì)運(yùn)移方程可以寫成

（2）揮發(fā)。揮發(fā)途徑主要考慮游離態(tài)或土壤吸附態(tài)的污染物進(jìn)入空氣，揮發(fā)速度Rv可以表示為［10］

式中：Dg為土壤－氣體擴(kuò)散系數(shù)；Kdu為包氣帶中污染物土水分配系數(shù)。

式（9）可以寫為

假設(shè)土壤體積含水率不隨時間和空間變化，則式（10）可寫為

式中：A為總源匯系數(shù)；Cu為包氣帶污染物含量。

式（11）邊界條件為

式中：Cu（x，z，t）為包氣帶點(diǎn)（x，z）在t時刻的污染物含量；C0為土壤初始污染物含量；a為污染源衰減系數(shù)。

3.2.2 飽和帶方程

飽和帶污染物運(yùn)移的二維方程（平面x，y方向）可以寫為

式中：Cs為飽和帶中污染物含量；Dxs為飽和帶中x方向水動力彌散系數(shù)；Dys為飽和帶中y方向水動力彌散系數(shù)；Rs為飽和帶阻滯因子；ks為飽和帶污染物生物降解系數(shù)。

飽和區(qū)溶質(zhì)運(yùn)移方程的邊界條件為

式中：Cs（x，y，τ）為飽和帶點(diǎn)（x，y）在τ時刻的污染物含量；Cu（x，d，τ）為包氣帶點(diǎn)（x，d）在τ時刻的污染物含量，d為飽和帶埋深；C－u為污染羽由包氣帶到達(dá)飽和帶時的平均含量。飽和帶初始條件取τ時刻當(dāng)污染羽到達(dá)飽和帶中的平均含量（圖1）。

3.3 溶質(zhì)運(yùn)移解析解

3.3.1 包氣帶模型的解析解

根據(jù)包氣帶的溶質(zhì)運(yùn)移模型及邊界條件可以求得解析解，其表達(dá)式為［13］式中：αzu為包氣帶z方向水動力彌散度，αzu＝Dzu／υu。

3.3.2 飽和帶方程的解析解

飽和帶方程的解析解C（x，y，t）為

式中：D′xs＝Dxs／Rs，D′ys＝Dys／Rs，α′xs＝D′xs／υ′s，υ′s＝ υ／R，k′＝k／R，均為中間變量；

sssss為污染物在飽和帶介質(zhì)的土水分配系數(shù)，ρb為多孔介質(zhì)密度。

4 案例分析

本研究使用Matlab7.0語言編輯計(jì)算機(jī)程序，實(shí)現(xiàn)整套模型的求解。

4.1 自由入滲

自由入滲模型的建立主要應(yīng)用于以下兩種情形：①突發(fā)性污染事故發(fā)生后的短期模擬；②中國北方大部分干旱及半干旱地區(qū)的場地污染模擬。自由入滲模擬過程中，污染物的最大飽和度為土壤的有效孔隙度（0.3）。選取甲醇作為污染物，模擬計(jì)算主要參數(shù)見表1。模擬結(jié)果以污染物飽和度在深度上的變化規(guī)律表示（圖3），模擬結(jié)果顯示甲醇的自由入滲速度非?？?，在污染發(fā)生21d后污染物就遷移至地表下5.5m處，但自由入滲受污染物自身性質(zhì)影響較大，如污染物的黏度會直接影響滲透系數(shù)，從而改變污染物向下遷移的速度。

表1 自由入滲模擬參數(shù)Tab.1 Parameters of free infiltration simulation

4.2 降雨入滲

降雨入滲一般應(yīng)用于長期模擬，或用于污染事故發(fā)生時有降雨條件下的短期模擬。污染發(fā)生后首先對包氣帶污染物的遷移進(jìn)行模擬，污染源假定為1 000cm的長條區(qū)域，重點(diǎn)模擬剖面（x為－1 000～1 000cm，z為－800～0cm）污染物的遷移情況，地下水位埋深假定為800cm，詳細(xì)模擬參數(shù)見表2。表2中a反映了污染發(fā)生后，一些自然或人為因素引起的污染物總量減少或含量變化情況，如污染發(fā)生后的降雨會引起部分污染物隨地表徑流遷移至其他地點(diǎn)，同時污染物含量會因降雨而下降，本案例未考慮這些因素。模擬時間為14d，模擬結(jié)果以阿特拉津含量隨時間的變化等值線表示（圖4）。

圖3 甲醇自由入滲模擬計(jì)算結(jié)果Fig.3 Results of free infiltration simulation for methanol

包氣帶阿特拉津含量模擬結(jié)果顯示，污染物主要在重力作用下向下遷移，同時由于分子擴(kuò)散作用，污染羽在水平方向也逐漸擴(kuò)大。由圖4可以看出，污染發(fā)生14d后污染羽達(dá)到飽和帶，污染羽的前端含量為29.94×10－6。

圖4 包氣帶降雨入滲阿特拉津含量等值線模擬結(jié)果Fig.4 Results of isolines of atrazine contents with rainfall infiltration simulation for unsaturated zone

表2 包氣帶降雨入滲模擬參數(shù)Tab.2 Parameters of rainfall infiltration simulation for unsaturated zone

表3 飽和帶參數(shù)Tab.3 Parameters for saturated zone

由于包氣帶的污染物在14d后到達(dá)含水層，飽和帶的模擬從第14天開始，主要模擬參數(shù)見表3。污染發(fā)生后第20、30天飽和帶的模擬結(jié)果見圖5。圖5描述了污染發(fā)生第20、30天時地下水下游100、 300、500、1 000cm處污染物含量的分布情況。圖5中y方向距離即為地下水下游位置。模擬結(jié)果顯示，盡管到達(dá)飽和帶的污染物濃度比包氣帶中小很多，但隨地下水遷移速度非常之快以至于污染物到達(dá)飽和帶后6d（污染發(fā)生20d）即有污染物遷移至下游1 000cm處。按第14天包氣帶污染物遷移至飽和帶的含量為30×10－6進(jìn)行計(jì)算，污染物在第6天（污染發(fā)生后第20天）到達(dá)下游1 000cm處的含量為初始含量的0.057%。而飽和帶遷移16d（污染發(fā)生30d）后，污染物到達(dá)1 000cm處的含量占初始含量的1.73%。污染物進(jìn)入地下水后，遷移的速度主要取決于地下水的流動速度，y方向由于只存在分子擴(kuò)散，污染物遷移速度相對較慢。因此，在突發(fā)性污染事故發(fā)生后應(yīng)及時采取強(qiáng)力有效的措施避免污染物進(jìn)入飽和帶，否則一旦污染物進(jìn)入含水層后將無限加大后續(xù)污染治理的難度。

圖5 飽和帶降雨入滲阿特拉津含量模擬結(jié)果Fig.5 Results of atrazine contents with rainfall infiltration simulation for saturated zone

從以上計(jì)算過程可看出，本模型采用的解析求解比數(shù)值求解已簡單了很多。盡管如此，由于考慮的機(jī)制較多（尤其是對于降雨入滲模型），求解過程中所需要的參數(shù)仍然較繁雜。實(shí)際上對于降雨入滲模型，最重要的參數(shù)是污染物的土水分配系數(shù)，分配系數(shù)可以通過污染物的辛醇－水分配系數(shù)或溶解度估算得到，而對于其他參數(shù)（如揮發(fā)速率因子、生物降解系數(shù)以及蒸騰因子）可較易從相關(guān)文獻(xiàn)或書籍中查得。污染物的水動力彌散系數(shù)可通過經(jīng)驗(yàn)公式計(jì)算得到，而其他水文氣象參數(shù)可直接從當(dāng)?shù)貧庀蟛块T獲取。針對不同的污染物以及地形條件，可以選擇性地忽略模型中的一些機(jī)制（如分子擴(kuò)散、生物降解），這樣可以有效簡化模型的計(jì)算，在突發(fā)性污染發(fā)生后更及時地進(jìn)行模擬預(yù)測。

5 結(jié) 語

（1）針對突發(fā)性場地污染，建立了能夠應(yīng)用于一般性污染場地并且適用于各種污染物的包氣帶、飽和帶污染物遷移數(shù)學(xué)模型，并得到相應(yīng)的解析解。相比于國外建立的數(shù)值模型，該模型具有參數(shù)少、操作簡單方便的優(yōu)點(diǎn)。

（2）該模型分為自由一維自由入滲模型和降雨入滲條件下二維污染物遷移模型，入滲模型主要應(yīng)用于短期模擬或中國干旱、半干旱地區(qū)的污染物遷移模擬；降雨入滲條件下污染物遷移二維模型主要考慮的機(jī)制包括對流和分子擴(kuò)散作用、污染物的揮發(fā)、植物根系吸收、生物降解、土壤吸附，主要用于長期模擬或突發(fā)降雨模擬。

（3）通過案例計(jì)算驗(yàn)證了該模型的有效性，模擬結(jié)果顯示該模型能夠快捷、高效地模擬污染物在包氣帶和飽和帶的遷移轉(zhuǎn)化，能較好地適用于突發(fā)性場地污染的污染物模擬工作。

［1］ 張嶸嶸.城市突發(fā)性環(huán)境污染事故預(yù)警技術(shù)支持系統(tǒng)的研究［D］.昆明：昆明理工大學(xué)，2007.ZHANG Rong－rong.Study on the warning technical support system of environmental pollution accident［D］.Kunming：Kunming University of Science and Technology，2007.（in Chinese）

［2］ 趙娜娜，黃啟飛，易愛華，等.我國污染場地的管理現(xiàn)狀與環(huán)境對策［J］.環(huán)境科學(xué)與技術(shù)，2006，29（12）：39－40.ZHAO Na－na，HUANG Qi－fei，YI Ai－h(huán)ua，et al.Contaminated sites management in China：current status and control［J］.Environmental Sciences and Technology，2006，29（12）：39－40.（in Chinese）

［3］ ZAREMBAMA.Lovecanal：anintroduction［EB／OL］.（2004－08－07）［2011－03－12］.http：／／www.onlineethics.org／resources／cases／lcanal.aspx.（in Chinese）

［4］ 新京報.北京地鐵5號線掘出有毒氣體，工期影響尚難定論［EB／OL］.（2004－05－01）［2011－03－01］.http：／／news.sina.com.cn／c／2004－05－01／03453180322.shtml.The Beijing News.The toxic gases were found in Line 5of Beijing Subway，so the impact of the gases on construction period was not determined［EB／OL］.（2004－05－01）［2011－03－01］.http：／／news.sina.com.cn／c／2004－05－01／03453180322.shtml.（in Chinese）

［5］ 張 升，盈尚軒.淮河干流突發(fā)性污染數(shù)值模擬的應(yīng)用研究［J］.水利科技與經(jīng)濟(jì)，2010，16（8）：845－846.ZHANG Sheng，YING Shang－xuan.Study on the application of numerical model for the sudden pollution accident in the main stream of Huai River［J］.Water Conservancy Science and Technology and Economy，2010，16（8）：845－846.（in Chinese）

［6］ 張防修，王艷平，劉興盛，等.黃河下游突發(fā)性污染事件數(shù)值模擬［J］.水利學(xué)報，2007，38（增）：613－618.ZHANG Fang－xiu，WANG Yan－ping，LIU Xing－sheng，et al.Numerical model for the sudden pollution accident in the lower of the Yellow River［J］.Journal of Hydraulic Engineering，2007，38（S）：613－618.（in Chinese）

［7］ RATHFELDER K M，LANG J R，ABRIOLA L M.A numeri－cal model（MISER）for the simulation of coupled physical，chemical and biological processes in soil vapor extraction and bioventing systems［J］.Journal of Contaminant Hydrology，2000，43（3／4）：239－270.

［8］ CHENG P，WANG C Y.A multiphase mixture model for multiphase，multicomponent transport in capillary porous media II：numerical simulation of the transport of organic compounds in the subsurface［J］.International Journal of Heat and Mass Transfer，1996，39（17）：3619－3632.

［9］ BELTMAN W H J，BOESTEN J J T I，VAN DER ZEE S E A T M.Analytical modelling of pesticide transport from the soil surface to a drinking water well［J］.Journal of Hydrology，1995，169（1／2）：209－228.

［10］ HANTUSH M M，MARINO M A.An analytical model for the assessment of pesticide exposure levels in soils and groundwater［J］.Environmental Modeling and Assessment，1996，1（4）：263－276.

［11］ JIAO J J，TANG Z H.An analytical solution of groundwater response to tidal fluctuation in a leaky confined aquifer［J］.Water Resources Research，1999，35（3）：747－751.

［12］ VAN DER ZEE S E A T M，BOESTEN J J T I.Effects of soil heterogeneity on pesticide leaching to groundwater［J］.Water Resources Research，1991，27（12）：3051－3063.

［13］ CLEARY R W，UNGS M J.Groundwater pollution and hydrology：mathematical models and computer programs［M］.Princeton：Princeton University Press，1978.