秦歡歡， 喻立珊， 鄭春苗

(1.東華理工大學(xué)省部共建核資源與環(huán)境教育部重點實驗室培育基地，江西 南昌 330013；2 北京大學(xué)工學(xué)院水資源研究中心，北京 100871；3.南方科技大學(xué)環(huán)境科學(xué)與工程學(xué)院，廣東 深圳 518055)

非飽和帶是大氣水和地表水同地下水發(fā)生聯(lián)系并進行水分交換的地帶，由于污染物大多是隨水?dāng)U散而發(fā)生遷移，因此它是地表污染物向下遷移污染地下水的通道。隨著中國經(jīng)濟的發(fā)展，一味追求經(jīng)濟增長而不重視環(huán)保的弊端逐步顯現(xiàn)，各種污染事件時有發(fā)生，如2013年山東濰坊地下水污染事件就一度引起公眾的關(guān)注(Zheng et al., 2013)；山西省太原市工業(yè)結(jié)構(gòu)以冶金、化工、煤炭、電力、機械為主，人口密集，人類活動頻繁，城市建設(shè)發(fā)展迅速，由此導(dǎo)致的地下水污染非常嚴重(蔣方媛等，2009)；近幾年來隨著大同縣社會經(jīng)濟的快速發(fā)展，污染物的大量排放，地下水水質(zhì)呈現(xiàn)出逐漸惡化的趨勢(張哲等，2012)；工業(yè)固體廢物填埋所產(chǎn)生的滲濾液等污水對地下水環(huán)境及人類健康存在潛在的威脅(程曉燕等，2016；楊文琳等，2017)。據(jù)《全國地下水污染防治規(guī)劃2011—2020》，全國地下水符合Ⅰ—Ⅲ類水質(zhì)標(biāo)準的占63%，Ⅳ—Ⅴ類的占37%，總體情況算好，但對北京、遼寧、吉林、上海、江蘇、海南、寧夏和廣東等8個省(區(qū)、市)641眼井的分析中水質(zhì)為Ⅳ—Ⅴ類的占73.8%，可見在經(jīng)濟發(fā)達區(qū)地下水污染狀況堪憂。因此，開展非飽和帶污染物遷移轉(zhuǎn)化研究，深入掌握其規(guī)律，有利于針對性地制定污染預(yù)防和治理措施。

非飽和帶是巖土顆粒、水、空氣三相共存的復(fù)雜系統(tǒng)，研究困難極大。近年來，國外對這方面開展研究并做了大量工作，把它當(dāng)作為研究地下水污染問題的突破口(劉仁冬，2005)，相比之下我國這方面研究開展得較晚。迄今為止，針對溶質(zhì)在土壤、非飽和帶中遷移轉(zhuǎn)化規(guī)律進行了大量研究，取得了大量實驗資料和理論進展，主要采用三種方法：室內(nèi)土柱實驗、現(xiàn)場實驗和數(shù)值分析法。但是，這三者都有一定的局限性。土柱實驗不能模擬實際土體的重力效應(yīng)，也不適于做污染物遷移的長期觀測，多數(shù)土柱實驗中裝置較小，缺乏與實際情況的類比性；現(xiàn)場實驗由于歷時較長，受外界因素影響較大，實驗本身會造成環(huán)境污染，所用示蹤劑受到限制，各種邊界條件也不易控制(Mattson et al., 2010)；由于缺少足夠的數(shù)據(jù)，數(shù)值模型的檢驗和校正變得非常困難，也很難對污染擴散過程做出準確預(yù)測。

1980年代，離心機開始應(yīng)用于污染物在土壤中的遷移研究。目前，離心模擬技術(shù)已被廣泛應(yīng)用于非飽和帶溶質(zhì)遷移研究(Timms et al., 2009; 張建紅等，2006a，2006b; Ataie-Ashtiani et al., 2003; Kumar, 2006)?？傮w來說，離心機可分為土工離心機(大中型)和小型離心機。前者因其操作空間更大，可控制變量更多并能通過技術(shù)手段實現(xiàn)動態(tài)數(shù)據(jù)采集，得到了更廣泛的應(yīng)用，模擬的對象也從惰性金屬離子(Mckinley et al., 1998)擴展到重金屬(張建紅等，2004)和非水相流體(胡黎明等，2003；Pantazidou et al., 2000)。近來則主要集中在利用離心機研究非飽和帶非均質(zhì)性對溶質(zhì)遷移的影響(Menezes et al., 2011)、粘土弱透水層中溶質(zhì)遷移(Timms et al., 2009)和污染物的去除(Timms et al., 2009; 郝榮福等，2004)。

研究非飽和帶溶質(zhì)遷移對地下水污染防治具有重要意義，而離心模擬技術(shù)在低滲透土壤長期污染物遷移及獲取數(shù)據(jù)進行模型驗證等方面體現(xiàn)了優(yōu)越性。然而，目前仍缺乏將離心模擬技術(shù)應(yīng)用于非飽和帶溶質(zhì)遷移規(guī)律研究是否可行的系統(tǒng)、全面的論述。因此，本文以非飽和帶溶質(zhì)遷移離心模擬的可行性為目標(biāo)，通過建立數(shù)值模擬，試圖在理論上分析、總結(jié)該技術(shù)在非飽和帶各類溶質(zhì)遷移規(guī)律研究中的可行性及適用條件，為離心模擬技術(shù)的發(fā)展及地下水污染物的防治提供科學(xué)的支持。

1 材料與方法

非飽和帶溶質(zhì)遷移和水分遷移密不可分，溶質(zhì)須溶于水才具備可遷移性，非飽和帶水分遷移模型在很大程度上能幫助理解污染物的遷移行為，但不能直接對某種類型污染物的擴散遷移進行分析。地下水中常見污染物有保守性鹽類(如NaCl)、揮發(fā)性有機物(如BTEX)、非水相流體(石油類)、重金屬(如Cd)和放射性核素(如鈾)等，理化過程包括微生物過程、吸附、解吸、氣化、溶解、化學(xué)反應(yīng)、放射衰變等，涉及水、氣、固三相，非常復(fù)雜。采用常規(guī)實驗方法進行研究具有相當(dāng)大的困難，而且類似NAPLs和重金屬等遷移緩慢物質(zhì)很難在實驗室內(nèi)開展長期觀測，離心模擬技術(shù)在這類物質(zhì)的遷移研究中具有較大的優(yōu)越性。筆者在喻立珊(2015)建立的非飽和帶水分遷移離心模型的基礎(chǔ)上開發(fā)非飽和帶溶質(zhì)遷移離心模擬理論及數(shù)值模型，評價將該技術(shù)應(yīng)用于非飽和帶溶質(zhì)遷移研究的可行性。鑒于篇幅，本文只把重點放在溶質(zhì)型鹽類污染物上，考慮的理化過程包括吸附和簡單的化學(xué)反應(yīng)等幾種。有關(guān)離心模擬相關(guān)的基礎(chǔ)理論，如離心相似比、模擬溶質(zhì)遷移的相似基礎(chǔ)、非飽和帶水分遷移離心模擬理論模型等，詳見喻立珊(2015)，喻立珊等(2014)等文獻。

1.1 非飽和帶溶質(zhì)遷移模型

非飽和帶溶質(zhì)遷移控制方程有多種，主要和目標(biāo)溶質(zhì)A的遷移特性有關(guān)。

(1)當(dāng)A為保守性物質(zhì)，忽略它和土壤之間的吸附與解吸過程，控制方程為：

(1)

式中，cm是溶液濃度，tm是模擬時間，vm是滲流速度，r為半徑，Dm是水動力彌散系數(shù)。

(2)當(dāng)A為保守性物質(zhì)不參與化學(xué)反應(yīng)，假設(shè)A與土壤之間的吸附與解吸瞬時達到平衡，控制方程為：

(2)

式中，Rm=1+Kdρb/θm為A的遲滯系數(shù)，Kd是吸附系數(shù)，ρb為土壤容重，θm為土壤體積含水率，飽和時為土壤孔隙度。

(3)當(dāng)A為反應(yīng)性物質(zhì)時，考慮吸附與解吸作用，假設(shè)能瞬時達到平衡，控制方程為：

(3)

式中，cAm，cBm分別為溶質(zhì)A和產(chǎn)物溶質(zhì)B的濃度，ka，kb分別為反應(yīng)速率常數(shù)，DAm，DBm分別為水動力彌散系數(shù)，RAm和RBm分別為遲滯系數(shù)。

1.2 數(shù)值計算方法

1.2.1 數(shù)值化方案

離心狀態(tài)下非飽和帶溶質(zhì)遷移的控制方程和常規(guī)情況一致，不同之處體現(xiàn)在滲流速度vm。三種情景下非飽和帶溶質(zhì)遷移過程分別由(1)—(3)描述，數(shù)值化方案可分為兩類：(1)和(2)為一類，(1)可看成是A的遲滯系數(shù)Rm為1時(2)的特殊情況；(3)為一類。本文使用算子分裂法對微分方程進行數(shù)值化。算子分裂法是眾多數(shù)值化方法中的一種，是分步驟數(shù)值化的方法，Yanenko(1971)將其進行推廣應(yīng)用。它將對流彌散方程分為對流算子和彌散算子，將反應(yīng)遷移方程分為遷移算子和反應(yīng)算子。

(1)方程(2)的數(shù)值化。 先將方程分裂成對流算子和彌散算子：

(4)

將對流算子顯式向前差分：

(5)

(6)

進一步整理：

(7)

圖1 算子分裂迭代法解對流彌散方程思路Fig.1 Process of the operator-split method solving convection diffusion equation

(2)方程(3)的數(shù)值化。將方程分裂成遷移算子和反應(yīng)算子：

(8)

將遷移算子用隱式差分進行離散：

(9)

(10)

式中

(11)

1.2.2 初始和邊界條件

鑒于本文研究地下水污染，故假定離心模型和原型的土在初始時刻不含溶質(zhì)，即初始條件為：

cm(rt≤r≤rb,tm=0)=0

(12)

考慮到現(xiàn)實情況下溶質(zhì)污染地下水的一般場景為：含有污染物的垃圾散布在地表，在一個雨天，垃圾中的溶質(zhì)淋溶而出并隨雨水下滲污染土壤及地下水。因此，本文考慮的邊界條件為：一維土柱頂部提供一個短暫的溶液供給，溶質(zhì)A的濃度先維持不變后轉(zhuǎn)為0 g/m3；上邊界水力條件從固定基質(zhì)吸力為0 kPa轉(zhuǎn)變?yōu)樗魍繛? m/s的邊界；模型底部為自由流邊界，基質(zhì)吸力固定為0 kPa，遷移模擬過程中處理成彌散通量固定的邊界條件。用數(shù)學(xué)表達為：

(13)

式中，c0為初始濃度，tpulse為降雨事件時長。方程(13)中第三式可離散為：

(14)

2 結(jié)果與分析

為方便和原型進行比較，筆者利用Visual Basic編寫了自然條件下非飽和帶溶質(zhì)遷移數(shù)值模擬代碼，參照前文闡述的數(shù)值化方案，對自然條件及離心狀態(tài)下的溶質(zhì)遷移過程進行計算。

圖2 非飽和帶保守性溶質(zhì)遷移離心模擬可行性比較(圖中標(biāo)注均為原型時間)Fig.2 Feasibility comparison of centrifugal modeling of nonreactive solute transport in unsaturated zone (time in the figure is the prototype time)

2.1 保守性溶質(zhì)

由于保守性溶質(zhì)既不發(fā)生吸附也不發(fā)生化學(xué)反應(yīng)，從離心相似理論上看，離心模型能夠正確反映原型下溶質(zhì)遷移行為而被廣泛應(yīng)用于非飽和帶溶質(zhì)遷移機理的研究(Nakajima et al., 1998; Kumar, 2006)。此處利用數(shù)值模擬結(jié)果的比較考察到底是否可行。數(shù)值模擬的固定參數(shù)及值為：θr=0.034，θs=0.46，nv=1.37，α=0.016 cm-1，Ks=1×10-7m/s，rb=2 m，其中θr是殘余含水率，θs是飽和含水率，nv和α是描述土水特征曲線的Van Genuchten模型參數(shù)(喻立珊，2015)，Ks是土壤的飽和水力傳導(dǎo)系數(shù)，rb為離心模型底部半徑。離心模型和原型空間離散節(jié)點數(shù)n= 101，A濃度為100 g/m3。原型的時長共100天，下雨時長tpulse為75天，原型非飽和帶厚度為2 m，其它變動參數(shù)的取值列于表1，離心模型和原型模擬結(jié)果對比如圖2所示。

表1 保守性溶質(zhì)情況下數(shù)值模擬參數(shù)

注：D為水動力彌散系數(shù)；θ0是初始含水率；N是重力水平；表示離心模型是在N個g的加速度環(huán)境下開展實驗；L是模型高度；Δt是時間步長；LX表示離心模型；YX表示原型。

從圖2a可看到，離心狀態(tài)下非飽和帶保守性溶質(zhì)的遷移和原型基本一致，但在靠近底部有一定偏差，此處模型中的遷移要稍滯后于原型。由于離心加速度存在徑向分布，模型2/3高度處至底部的這一部分相似比被低估，空間尺寸被縮小，導(dǎo)致水分及溶質(zhì)遷移均有一定的滯后。另外，N為40比20時底部的偏差要小，說明同一設(shè)備使用更大的N值有利于減輕這種偏差。圖2b展現(xiàn)了分子彌散系數(shù)對溶質(zhì)遷移的影響：分子彌散系數(shù)越大，溶質(zhì)的污染范圍更大，但最高濃度峰值會變小。這是因為分子彌散系數(shù)越大，它對溶質(zhì)濃度剖面曲線的“削峰填谷”的能力越強。圖2c比較了不同初始含水率對溶質(zhì)遷移的影響，初始含水率越低，土壤對溶質(zhì)的遷移阻礙越大。由于孔隙度和滲透系數(shù)之間關(guān)系不太明確且太過復(fù)雜，故沒有考察土壤壓實度對溶質(zhì)遷移的影響。總體來看，離心模擬能較好地反映保守性溶質(zhì)在原型下的遷移行為，而要減小離心加速度分布不均產(chǎn)生的誤差，可通過嚴格計算不同部位的重力水平來正確處理模型尺寸的縮放從而避免這種誤差，或者使用更大的N值來盡力減小這種誤差。

圖3 非飽和帶吸附性溶質(zhì)遷移離心實驗?zāi)M可行性比較(圖中標(biāo)注均為原型時間)Fig.3 Feasibility comparison of centrifugal modeling of adsorption solute transport in unsaturated zone (time in the figure is the prototype time)

2.2 吸附性溶質(zhì)

在非飽和帶保守性溶質(zhì)遷移基礎(chǔ)之上，考慮吸附作用的影響。數(shù)值模擬的固定參數(shù)及值為：θr=0.034，θs=0.46，nv=1.37，α=0.016 cm-1，Ks=1×10-7m/s，rb=2 m，Dm=20 cm/d，ρb=1.7 g/cm3，離心模型和原型空間離散節(jié)點數(shù)n=101，初始含水率為0.2(均勻分布)，A濃度為100 g/m3。原型的時長共100天，下雨時長tpulse為75天，原型非飽和帶厚度為2 m，其它變動參數(shù)的取值列于表2，離心模型和原型模擬結(jié)果對比如圖3所示。

表2 吸附性溶質(zhì)情況下數(shù)值模擬參數(shù)

注：N是重力水平，表示離心模型是在N個g的加速度環(huán)境下開展實驗；L是模型高度；Δt是時間步長；Kd是吸附系數(shù)；LX表示離心模型；YX表示原型。

從圖3a可以看出，離心模擬的結(jié)果和原型吻合度較好，說明該技術(shù)用于非飽和帶瞬時平衡吸附性溶質(zhì)的遷移研究是可行的。和圖2a相比較，可知吸附作用的存在會明顯阻滯溶質(zhì)的遷移擴散，對于控制污染的范圍具有積極意義。但溶質(zhì)會長期殘留在土壤中，在農(nóng)業(yè)上對作物有影響。圖3b比較了分配系數(shù)對溶質(zhì)遷移的影響，分配系數(shù)越大，土壤對溶質(zhì)的遲滯作用越強，導(dǎo)致溶質(zhì)遷移得更加緩慢，污染范圍也更加集中。

比較遺憾，本文只討論了平衡吸附，而實際情況下非平衡吸附更加普遍。對于非平衡吸附，目前常用雙區(qū)模型來描述其遷移過程，如Li等(1993)，離心模擬技術(shù)在非平衡吸附下非飽和帶溶質(zhì)的遷移研究中也是可行的。

圖4 非飽和帶反應(yīng)性溶質(zhì)遷移離心實驗?zāi)M可行性比較(圖中標(biāo)注均為原型時間)Fig.4 Feasibility comparison of centrifugal modeling of reactive solute transport in unsaturated zone (time in the figure is the prototype time)

2.3 反應(yīng)性溶質(zhì)

研究非飽和帶伴隨產(chǎn)物溶質(zhì)B出現(xiàn)的反應(yīng)性溶質(zhì)A的遷移過程，著重考察B的性質(zhì)對A遷移過程的影響。數(shù)值模擬的固定參數(shù)及值為：θr=0.034，θs=0.46，nv=1.37，α=0.016 cm-1，Ks=1×10-7m/s，rb=2 m，DAm= 20 cm/d，DBm=20 cm/d，ρb=1.7 g/cm3，KAd=0.000 1 m3/kg，KBd=0.000 1 m3/kg，離心模型和原型空間離散節(jié)點數(shù)n=101，初始含水率為0.2(均勻分布)，A濃度為100 g/m3。原型的時長共100天，下雨時長tpulse為75天，原型非飽和帶厚度為2 m，其它變動參數(shù)的取值列于表3，離心模型和原型模擬結(jié)果對比如圖4。

表3 反應(yīng)性溶質(zhì)情況下數(shù)值模擬參數(shù)

注：N是重力水平，表示離心模型是在N個g的加速度環(huán)境下開展實驗；L是模型高度；Δt是時間步長；ka和kb為A和B的反應(yīng)速率常數(shù)，YXF代表考慮反應(yīng)性溶質(zhì)的原型，LXF代表考慮反應(yīng)性溶質(zhì)的離心模型。

從圖4a，e可看出，化學(xué)反應(yīng)的存在會使得離心模擬的結(jié)果偏離原型。和放射性核素一樣，當(dāng)反應(yīng)進行得比較慢時，離心模擬A的遷移過程能較好地和原型匹配；但若反應(yīng)較為迅速，則誤差偏大，往往會高估A的遷移擴散能力和污染區(qū)A的濃度。再分別對比圖4b，c，f，g，可知對于產(chǎn)物溶質(zhì)B，離心模擬的結(jié)果和原型相差很大，濃度的預(yù)測值誤差能達到2個數(shù)量級。進一步仔細觀察可發(fā)現(xiàn)原型中土壤頂部附近區(qū)域A的濃度分布有明顯的陡降特征，而B在頂部區(qū)域分布集中，這些特征在離心模型中無法體現(xiàn)，表明離心模擬技術(shù)用于反應(yīng)較快的非飽和帶溶質(zhì)遷移研究時很難抓住過程的細節(jié)。

從模擬結(jié)果可以發(fā)現(xiàn)，離心模擬技術(shù)對于反應(yīng)過程相對較慢的溶質(zhì)A的遷移過程預(yù)測能夠較好地和原型匹配，但當(dāng)反應(yīng)過程較快時，離心模擬將產(chǎn)生較大的誤差，對于遷移過程的細節(jié)捕捉能力也不強。此外，對B的預(yù)測，離心模型和原型相差甚遠，濃度誤差能達到2個數(shù)量級，表明離心模擬技術(shù)在對B生成和遷移過程的預(yù)測上是失敗的。若B也是研究的重點，那么離心模擬技術(shù)是不值得推薦的。反之，若B在研究中不那么重要，那么在適當(dāng)?shù)臈l件(反應(yīng)進程不快)下，離心模擬技術(shù)可行。

3 結(jié)論

在水分遷移模型的基礎(chǔ)上，建立了離心場內(nèi)一維非飽和帶溶質(zhì)遷移模型及數(shù)值化方案，對離心模擬技術(shù)應(yīng)用于三種類型溶質(zhì)在非飽和帶遷移研究的可行性進行了研究，結(jié)果表明離心模擬技術(shù)在適當(dāng)條件下用于非飽和帶溶質(zhì)遷移規(guī)律研究是可行的。具體結(jié)論如下：

(1)離心模擬可用于保守性溶質(zhì)的遷移研究。離心模型能準確反映原型下這類溶質(zhì)的遷移過程，但離心加速度的徑向分布會使得靠近底部區(qū)域的溶質(zhì)遷移稍滯后于原型，使用更大的重力水平N值，能使得加速度分布更加均勻從而減小誤差。

(2)分子彌散系數(shù)和土壤的初始含水率對溶質(zhì)遷移有影響。分子彌散系數(shù)越大，對溶質(zhì)濃度剖面曲線的“削峰填谷”能力越強，使得溶質(zhì)分布更加均勻，導(dǎo)致污染范圍更大，而濃度峰值會變小。土壤的初始含水率越低，溶質(zhì)的遷移擴散越慢。

(3)離心模擬適用于吸附瞬時平衡溶質(zhì)遷移規(guī)律研究。吸附與解吸作用的存在對溶質(zhì)遷移有明顯影響，分配系數(shù)越大，土壤對溶質(zhì)的遲滯作用越明顯，對于控制污染范圍具有積極意義。

(4)對于有產(chǎn)物溶質(zhì)B生成的反應(yīng)性溶質(zhì)A在非飽和帶的遷移，當(dāng)反應(yīng)較慢時，離心模型預(yù)測的A的遷移過程能較好地匹配原型，但對B的預(yù)測誤差較大，濃度誤差可達2個數(shù)量級。若B是研究重點，則不推薦離心模擬技術(shù)；若B不那么重要，則在適當(dāng)條件(反應(yīng)進程不快)下，離心模擬技術(shù)是可行的。

參考文獻

程曉燕，江寶鋒，劉傳杰. 2016.危險廢物填埋場對地下水環(huán)境影響分析[J].東華理工大學(xué)學(xué)報：自然科學(xué)版，39(增刊)：186-189.

郝榮福，胡黎明，邢巍巍.2004.土壤中可揮發(fā)性污染物清除的離心試驗研究[J].巖土力學(xué)， 25(7)：1037-1040.

胡黎明，郝榮福，殷昆亭，等.2003.BTEX在非飽和土和地下水系統(tǒng)中遷移的試驗研究[J]．清華大學(xué)學(xué)報：自然科學(xué)版， 43(11)：1546-1549.

蔣方媛，郭清海. 2009.城市地下水污染成因分析——以山西省太原市為例[J].東華理工大學(xué)學(xué)報：自然科學(xué)版，32 (1)：82-88.

劉仁冬. 2005.非飽和帶中水分遷移和污染物運移的數(shù)值模擬[D].沈陽：東北大學(xué).

楊文琳，張靜，袁野，等. 2017.基于有限差分法的地下水溶質(zhì)運移模擬：以某垃圾填埋場為例[J]. 環(huán)境工程，35(12)：30-35.

喻立珊，曹國亮，許模，等. 2014.離心實驗在污染物遷移研究中的應(yīng)用[J].地球科學(xué)進展， 29(2)：227-237.

喻立珊. 2015.非飽和帶溶質(zhì)遷移離心實驗的數(shù)值模擬研究[D].北京：北京大學(xué).

張建紅，胡黎明. 2006a.重金屬離子和LNAPLs在非飽和土中的運移規(guī)律研究[J].巖土工程學(xué)報， 28(2)：276-280.

張建紅，呂禾，王文成. 2006b.銅離子在非飽和土中遷移的離心模型試驗研究[J].巖土力學(xué)， 27(11)：1885-18 90.

張建紅，嚴冬. 2004.非飽和粉質(zhì)砂土中銅離子遷移的離心模型試驗研究[J].巖土工程學(xué)報， 26(6)：792-797.

張哲，鄭秀清，陳軍峰，等. 2012.利用灰色聚類法評價大同縣地下水水質(zhì)[J].東華理工大學(xué)學(xué)報：自然科學(xué)版，35(3)：270-275.

Ataie-Ashtiani B, Hassanizadehb S M, Oungc O, et al. 2003. Numerical modelling of two-phase flow in a geocentrifuge [J]. Environmental Modelling and Software, 18(3):231-241.

Kumar P R. 2006.An experimental methodology for monitoring contaminant transport through geotechnical centrifuge models [J]. Environmental Monitoring & Assessment, 117(1-3):215-233.

Li L, Barry D A, Hensley P J, et al. 1993.Nonreactive chemical transport in structured soil: The potential for centrifuge modelling [A]//Geotechnical Management of Waste and Contamination[C].Rotterdam.

Mattson E D, Palmer C D, Smith R W, et al. 2010. Centrifuge techniques and apparatus for transport experiments in porous media [A]// Springman, Laue, Seward. Physical Modelling in Geotechnics[C]. London: Taylor & Francis Group.

Mckinley J D, Price B A, Lynch R J, et al. 1998. Centrifuge modelling of the transport of a pulse of two contaminants through a clay layer [J]. Geotechnique, 48(3):421-425.

Menezes G B, Ward A, Moo-Young H K. 2011. Unsaturated flow in anisotropic heterogeneous media: a centrifuge study [J]. International Journal of Hydrology Science and Technology, 1(3):147-163.

Nakajima H, Hirooka A, Takemura J, et al. 1998. Centrifuge modeling of one-dimensional subsurface contamination [J]. Journal of American Water Resources Association, 34(6):1415-1425.

Pantazidou M, Abu-Hassanein Z S, Michael F R. 2000. Centrifuge study of DNAPL transport in granular media [J]. Journal of Geotechnical and Geoenvironmental Engineering, 126(2):105-115.

Timms W, Hendry M J, Muise J, et al. 2009. Coupling centrifuge modeling and laser ablation Inductively coupled plasma mass spectrometry to determine contaminant retardation in clays [J]. Environmental Science & Technology, 43 (4):1153-1159.

Yanenko N N. 1971.The method of fractional steps [M].New York: Springer-Verlag.

Zheng C M, Liu J. 2013. China's “Love Canal” moment?[J]. Science, 340(6134):810-811.