污染物在填埋場底部非飽和土中的一維擴散解析解

關馳，謝海建，，陳云敏，唐曉武，樓章華

（1.浙江大學 建筑工程學院，浙江 杭州 310058；2.浙江大學 軟弱土與環境土工教育部重點實驗室，浙江 杭州 310058）

為防止和減緩降雨入滲和垃圾自身降解等產生的高污染物濃度填埋場滲濾液對地下水造成的污染，必須對滲濾液擴散引起的地下水污染問題進行研究與評價[1－7]。雖然填埋場滲濾液對地下水的污染問題可以通過數值模型分析和求解[8－10]，但解析解在許多方面是代替數值方法的一種經濟、有效的方法[11－12]。一方面解析解是數值解的理論基礎[13]；另一方面，數值方法的復雜程度往往較高（尤其在填埋場的初步設計階段），采用詳盡的數值模擬及其必要的發散性檢查所需要的代價（按工時）往往不合算[14]。

Liu and Ball[15]和Foose等[16]獲得了半無限空間污染物在土中的一維擴散解析解。Shackelford and Lee[17]，陳云敏等[18]，謝海建等[19]和Li and Cleall[20]分別針對不同的邊界條件和土層的分層性等因素，得到了有限厚度邊界時污染物在土中的一維擴散解析解。但這些解析解均基于飽和土層，許多情形下填埋場底部的防滲層是非飽和的[21－22]。Fityus等[21]對已有填埋場的調查發現，土層體積含水量達到穩態需要的時間遠小于填埋場污染物的運移時間；因此假設填埋場底部土層中體積含水量穩態分布是合理的，且由此導致的誤差可以忽略。對每一單層土而言，線性體積含水量分布是分析非飽和土中體積含水量分布的實用和簡便的方法[21]。以往學者主要針對污染物在非飽和含水層中的運移問題進行研究，水力梯度引起的對流作用和機械彌散作用往往是這類問題的主導作用，而分子擴散作用則往往被忽略[23]。因此，基于此得到的污染物在非飽和含水層中運移的解析解不能用于分析污染物在填埋場底部粘性土中的運移問題。

本文在假設填埋場底部非飽和土層內體積含水量線性分布的條件下，建立了污染物在填埋場底部非飽和土中的擴散模型，利用分離變量法結合廣義正交理論得到了該模型的解析解。通過與多物理場耦合分析軟件COMSOL Multiphysics 4.2[24]的計算結果比較，驗證了本文解析解的正確性和有效性。基于該解析解，研究了土層體積含水量變化對污染物擊穿非飽和土的影響。該解析解可用于數值模型的驗證和試驗數據的擬合等。

1 數學模型

1.1 基本假設

假設滲濾液中污染物濃度保持不變，為C0。底部的污染物濃度為0或污染物通量為0。并假設豎直向下為z的正方向，且z軸的原點在非飽和土頂端，土層厚度為L（見圖1）。

圖1 污染物在非飽和土中運移模型

該模型的主要假設有：（1）污染物在非飽和土層中運移是一維的；（2）土層是均質的；（3）僅考慮污染物在非飽和土中的擴散作用，不考慮對流作用的影響；及（4）土層體積含水量是深度的線性函數。

1.2 控制方程及邊界條件

在上述假設條件下，污染物在非飽和土中的一維擴散控制方程為[21，25]：

其中C（z，t）為土層孔隙水中污染物濃度（mg／L），隨深度z和時間t變化；D＊為污染物在非飽和土中的有效擴散系數（m2／s）；θ（z）為體積含水量，即單位體積土體中水的含量；這里假設其為深度的線性函數[21]

其中A、B為常數。

將式（2）代入式（1），得到：

式（3）即為假設體積含水量線性變化時污染物在非飽和土中的擴散方程。

當上部污染物濃度C0不變時，上邊界條件為：

若下邊界保持污染物濃度為0，則下邊界條件為：

當土層下方含水層的水流流速較大以致能及時帶走滲漏出的污染物時，適宜采用上述零濃度邊界條件。

若土層底部通量為零，則下邊界條件為：

當粘土襯墊的下臥含水層流速很小，或其下臥層是不透水的基巖時，可采用這種邊界條件。

在初始時刻，土層中污染物濃度假設為深度的某一函數：

其中f（z）為初始時刻非飽和土中污染物的濃度分布函數。

2 模型的解

在滿足邊界條件式（4）～（6）和初始條件式（7）下，控制方程（3）的解為（詳細的求解過程見附錄）：

J0為0階第一類貝塞爾函數，Y0為0階第二類貝塞爾函數，λm為方程的第m個特征根。

在式（8）濃度解的基礎上，可通過下式得到任意時刻t和任意位置z處污染物的通量J：

當下邊界為零濃度時：

或當下邊界為零通量時：

由式（8）可知，當t→∞時，即在穩態時，污染物的濃度分布即為φ（z）。

3 與數值解及飽和時解的比較

為了驗證本文解析解的有效性和合理性，比較了本文解和多物理場耦合分析軟件 COMSOL Multiphysics 4.2[24]的計算結果。計算比較時模型相關參數參考文獻[21]，取土層頂部體積含水量為0.36，底部為0.3795，土層厚度0.75m，擴散系數D＊為5×10－10m2／s。如圖2所示，經過10年時間的運移，由本文解得到的污染物濃度結果和對應數值計算結果符合得很好，這說明了本文解的正確性。

圖2 本文解與數值解的比較

圖3為在底部邊界為零濃度情況下改變上下體積含水量分別為0.3和0.6時，土層污染物相對濃度變化曲線本文解及飽和情況解比較，發現本文非飽和時解在同一深度較飽和時的解結果偏小，說明相對飽和土，污染物在非飽和土中隨深度增加其相對濃度變化較慢。

圖3 非飽和與飽和土中污染物相對濃度曲線

圖4為底部零濃度邊界條件下本文非飽和土的解與飽和時的解在土底端通量上的比較，非飽和土時頂端體積含水量為0.3，底端體積含水量為0.6，飽和時的孔隙度分別取為0.3和0.6。非飽和土底端污染物通量擊穿曲線介于飽和度0.3和0.6得到的通量擊穿曲線之間。100年時非飽和土底端污染物通量為91g·ha－1·a－1，是孔隙度為0.3時飽和土底端通量（63g·ha－1·a－1）的1.4倍，是孔隙度為0.6時通量（126g·ha－1·a－1）的0.7倍。這也說明非飽和土的體積含水量對污染物的通量擊穿曲線影響較大。

圖4 飽和與非飽和土層底端污染物通量

4 參數分析

為了從理論上考察土層不同體積含水量分布對污染物擴散的影響，本文選取了四種不同體積含水量分布。粘土孔隙率變化范圍為0.1－0.8[27]。本文選取各工況體積含水量取值見表1。工況1為土層頂部和底部體積含水量分別為0.3和0.6。工況2與工況1相反，即體積含水量在頂部和底部分別為0.6和0.3。工況3和4與工況1和2類似，僅土層上下部的體積含水量分別為0.5和0.55。

表1 四種工況體積含水量取值

圖5為四種不同工況下土層底部污染物相對濃度隨時間的變化曲線，即濃度擊穿曲線的變化。分析表明工況3和4結果介于工況1和2之間，曲線形狀非常類似。當非飽和土頂端體積含水量小于底端體積含水量時（工況1和3），經過5年時間的運移，達到非飽和土底端的污染物相對濃度較小；相反地，當非飽和底端體積含水量大于底端體積含水量時（工況2和4），經過相同時間，達到非飽和土底端污染物相對濃度較大。經過10年時間的運移，工況1和3達到非飽和土底端的相對濃度較小，分別為0.29和0.33。工況2和4對應的相對濃度較大，分別為0.44和0.40。這說明當非飽和土頂端體積含水量小于底端體積含水量時，污染物在該種條件下擴散運移地較慢；相反，在非飽和土頂端體積含水量大于底端體積含水量時，污染物運移較快。

圖5 不同工況污染物濃度擊穿曲線比較

圖6為不同工況下非飽和土底部污染物通量隨時間的變化曲線。工況1和2的底部通量差別不大。工況3和4的通量擊穿曲線在10年之前有所差別，但10年后逐漸趨于一致。不同于濃度擊穿曲線，土層上下端體積含水量變化較小的工況，即工況3和4的底部污染物通量較大，100年時其通量是工況1和工況2的1.09倍。

圖6 不同工況非飽和土底部污染物通量擊穿曲線比較

5 結 論

本文在假設土層體積含水量隨深度呈線性分布的基礎上，根據污染物在非飽和土中的擴散模型，利用分離變量法結合及廣義正交理論得到了該模型的解析解。通過與數值解法的比較驗證本文模型解的合理性。基于該解析解本文分析了土層不同體積含水量分布時污染物通過非飽和土的濃度和通量擊穿曲線，得到的主要結論有：

（1）本文解的計算結果和飽和情形的解的結果比較發現體積含水量變化對污染物運移影響較大；且隨著時間的增大，本文解得到的污染物濃度曲線逐漸接近穩態的情形，說明了本文解的正確性和合理性。

（2）污染物在非飽和土中的擴散較其在飽和土中慢。當土層上下端體積含水量分別為0.3到0.6時，5年后0.3m深處非飽和土中污染物相對濃度為0.38，而飽和土時為0.45。控制填埋場底部土層的飽和度有利于填埋場的防污。

（3）體積含水量的線性變化對非飽和土底部污染物通量變化影響較大，分別為飽和時最大體積含水量和最小體積含水量對應的土層底部污染物通量的0.7和1.4倍。

（4）當非飽和土層頂端體積含水量小于底端體積含水量時，污染物在非飽和土中擴散運移較慢；相反，當頂端體積含水量大于底端體積含水量時，污染物運移較快。

[1]郭永龍，王焰新，蔡鶴生，等.垃圾填埋場滲濾液對地下水環境影響的評價[J].地質科技情報，2002，21（1）：87-92.Guo Y L，Wang Y X，Cai H S，et al.Assessment of environment impact of leakage of landfill site on ground water[J].Geological Science and Technology Information，2002，21（1）：87-92.

[2]張文杰，陳云敏，詹良通.垃圾填埋場滲濾液穿過垂直防滲帷幕的滲漏分析[J].環境科學學報，2008，28（5）：925-929 Zhang W J，Chen Y M，Zhan L T.Transport of leachate through vertical curtain grouting in landfills[J].Acta Scientiae Circumstantiae，2008，28（5）：925-929.

[3]馬志飛，安達，姜永梅，等.某危險廢物填埋場地下水污染預測及控制模擬[J].環境科學，2012，33（1）：64-70.Ma Z F，AN D，Jiang Y M，et al.Simulation on Contamination Forecast and Control of Groundwater in a Certain Hazardous Waste Landfill[J].Chinese Journal of Environmental Science，2012，33（1）：64-70.

[4]馮世進，高麗亞，王印.垃圾填埋場邊坡上土工膜的受力分析[J].巖土工程學報，2008，30（10）：1484-1489.Feng S J，Gao L Y，Wang Y，et al.Analysis of tension of geomembranes placed on landfill slopes[J].Chinese Journal of Geotechnical Engineering，2008，30（10）：1484-1489.

[5]Liu C L，Zhang F E，Zhang Y，et al.Experimental and numerical study of pollution process in an aquifer in relation to a garbage dump field[J].Environmental Geology，2005，48（8）：1107-1115.

[6]Dong J，Zhao Y S，Zhang W H，et al.Laboratory study on sequenced permeable reactive barrier remediation for landfill leachate-contaminated groundwater[J].Journal of Hazardous Materials，2009，161（1）：224-230.

[7]薛強，趙穎，劉磊，等.垃圾填埋場災變過程的溫度-滲流-應力-化學耦合效應研究[J].巖石力學與工程學報，2011，30（10）：1970-1988.Xue Q，Zhao Y，Liu L，et al.Study of thermo-hydro-mechanical-chemical coupling effect of catastrophe process of landfill[J].Chinese Journal of Rock Mechanics and Engineering，2011，30（10）：1970-1988.

[8]Rowe R K，Quigley R M，Brachman R W L，et al.Barrier Systems for Waste Disposal[M].London and New York：Spon Press，2004

[9]黃冠華，黃權中，詹紅兵，等.均質介質中污染物遷移的分數微分對流 彌散模擬[J].中國科學D輯：地球科學，2005，35（增刊I）：249-254.Huang G H，Huang Q Z，Zhan H B，et al.Simulation of contaminants advection-diffusion in homogeneous medium [J].Science in China Series D：Earth Sciences，2005，35（Suppl）：249-254.

[10]何俊，高子坤，詹國才.有機污染物在復合襯墊中的運移分析[J].水文地質工程地質，2009，1：124-127.He J，Gao Z K，Zhan G C.Transport of organic contaminants through a composite liner[J].Hydrogeology and Engineering Geology，2009，1：124-127.

[11]謝海建，陳云敏，樓章華.污染物通過有缺陷膜復合襯墊的一維運移解析解[J].中國科學：技術科學，2010，40（5）：486-495.Xie H J，Chen Y M，Lou Z H.An analytical solution to contaminant transport through composite liners with geomembrane defects[J].Science China Technological Sciences，2010，40（5）：486-495.

[12]謝海建，樓章華，陳云敏，等.污染物通過GCL／AL防滲層的對流－彌散解析解[J].科學通報，2010，55（21）：2148－2155.Xie H J，Lou Z H，Chen Y M，et al.An analytical solution to contaminant advection and dispersion through a GCL／AL liner system [J].Chinese Science Bulletin，2010，55 （21）：2148-2155.

[13]Park E，Zhan H B.Analytical solutions of contaminant transport from finite one-，two-，and three-dimensional sources in a finite-thickness aquifer[J].Journal of Contaminant Hydrology，2001，53（1-2）：41-61

[14]Rowe R K，Nadarajah P.An analytical method for predicting the velocity field beneath landfills[J].Canadian Geotechnical Journal，1997，34：264-282.

[15]Liu C X，Ball W P.Analytical modeling of diffusion-limited contamination and decontamination in two-layer porous medium[J].Advances in Water Resources，1998，21：297-313.

[16]Foose G J，Benson C H，Edil T B.Analytical equations for predicting concentration and mass flux from composite liners[J].Geosynthetics International，2001，8（6）：551-575.

[17]Shackelford C D，Lee J M.Analyzing diffusion by analogy with consolidation[J].Journal of Geotechnical and Geoenvironmental，2005，131（11）：1345-1359.

[18]陳云敏，謝海建，柯瀚，等.層狀土中污染物的一維擴散解析解[J].巖土工程學報，2006（4）：521-524.Chen Y M，Xie H J，Ke H，et al.Analytical solution of contaminant diffusion through multi-layered soils [J].Chinese Journal of Geotechnical Engineering，2006（4）：521-524.

[19]謝海建，唐曉武，陳云敏，等.原始土層影響下成層介質污染物一維擴散模型[J].浙江大學學報（工學版），2006（12）：2191-2195.Xie H J，Tang X W，Chen Y M，et al.One-dimensional model for contaminant diffusion through layered media[J].Journal of Zhejiang University （Engineering Science），2006 （12）：2191-2195.

[20]Li Y C，Cleall P J.Analytical solutions for contaminant diffusion in double-layered porous media[J].Journal of Geotechnical and Geoenvironmental，2010，136（11）：1542-1554

[21]Fityus S G，Smith D W，Booker J R.Contaminant transport through an unsaturated soil liner beneath a landfill[J].Canadian Geotechnical Journal，1999，36：330-354

[22]Bonaparte R，Gross B A.Field behavior of double liner systems.In：Bonaparte R eds.Waste containment systems：construction，regulation and performance[J].New York：ASCE Geotechnical Special Publication，1990.52-83.

[23]Sander G C，Braddock R D.Analytical solutions to the transient，unsaturated transport of water and contaminants through horizontal porous media [J].Advances in Water Resources，2005，28（10）：1102-1111.

[24]COMSOL Multiphysics：Engineering simulation software.Version 4.2.Burlington：COMSOL，Inc，2011.

[25]席永慧.鋅和鎳在粉煤灰－粘土介質中擴散系數的測定[J].同濟大學學報：自然科學版，2007，35（11）：1520-1524.Xi Y H.Determination of Diffusion Coefficients of Zn2＋and Ni2＋in saturated fly ash-clay medium[J].Journal of Tongji U-niversity：Natural Science，2007，35（11）：1520-1524.

[26]zisik M N.Heat conduction[M].2nd edition.New York：John Wiley and sons，1993.

[27]袁聚云.土質學與土力學[M].人民交通出版社.北京，2009.

附 錄

令x＝A＋Bz，則z＝（x－A）／B，從而將控制方程及邊界條件轉化為：

為了將非齊次邊界條件（A2）轉化為齊次邊界條件，根據疊加原理，可假設C（x，t）的解具有的形式為[26]

對零通量下邊界條件（A15），控制方程（A12）解為

關于u（x，t）的求解可采用分離變量法，設控制方程（A7）的解具有的形式為

其中Γ（t）和R（x）分別為時間t和空間x的函數。

由控制方程（A7）可知，Γ（t）由下式確定

空間變量x的函數R（x）滿足如下特征根方程：

其中αm為待定參數。

將式（A24）代入初始條件（A11）得到

方程（A20）－（A25）是Sturn-Liouville問題的一個特例，因此R（λm，x）具有如下的正交性質[22]：

其中N（λm）由下式確定

Rm（x）是常微分方程（A20）的解，其具有的形式為

設a＝A，b＝A＋BL；將上式代入上邊界條件（A21）和下邊界條件（A22）、（A23）得到：

式中J1（x），Y1（x）分別為x的一階第一類貝塞爾函數和一階第二類貝塞爾函數。

因為β1m和β2m是不為零的參數，則由式（A29）和 （A30）或（A31）可得確定特征根λm的方程如下：

將式（A37）代入式（A24）可得：

兩邊同乘x[qmJ0（λmx）＋Y0（λmx）]再積分并由正交關系式（A26）、（A27）可以確定ηm為

從而控制方程（A1）的解為：

將x＝A＋Bz代入上式即可得到控制方程（1）的解為式（8）。