1、第 29 卷第 9 期 岩 土 力 学 Vol.29 No. 9 2008 年 9 月 Rock and Soil Mechanics Sep. 2008 收稿日期： 2006-12-26 基金项目：国家自然科学基金（ No. 50374041； 50309015； 50574048） ；农业部农业环境与农产品安全重点开放实验室开放基金（ No. 2004A03） 。 作者简介：梁冰，女， 1962 年生，教授，博士生导师，主要从事环境流体力方面的研究。 E-mail: 文章编号： 1000 7598 (2008) 09 2373 05 土壤中污染物迁移转化规律的数值模拟研究 梁 冰1，赵
2、颖1,2，薛 强2，刘 磊1,2 （ 1. 辽宁工程技术大学 力学与工程科学系，辽宁 阜新 123000； 2. 中国科学院武汉岩土力学研究所，武汉 430071） 摘 要：考虑到田间土壤的结构性以及运移其中的污染物的化学反应特征，建立了两区双点平衡 /动力学吸附溶质运移模型。用拟交替显隐式有限差分格式对模型进行了离散，通过编制的相应程序，对一施药期的水田进行了模拟。研究了土壤中污染物的迁移转化规律，分析了不可动水和非平衡吸附相对浓度分布的影响。结果表明：不可动区水相和非平衡吸附相的浓度具有的滞后现象，在污染前期一定程度上起到了延缓污染作用，而在污染后期又在一定程度上加重了污染。二者的存在均对
3、溶质的迁移转化具有重要的影响，进行数值模拟时应充分考虑二者的存在。模拟及分析结果为土壤环境质量评价提供了定量依据，同时对环境污染的治理与控制具有理论和现实意义。 关 键 词：污染物；土壤 -水环境；平衡 /非平衡双点吸附 /解吸；两区模型；数值模拟 中图分类号：TU 441 文献标识码：A Research on numerical simulation of contaminant transportation and transformation in soil LIANG Bing1, ZHAO Ying1,2, XUE Qiang2, LIU Lei1,2(1. Department
4、of Mechanics and Engineering Sciences， Liaoning Technical University, Fuxin 123000, China; 2. Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China) Abstract: Considering the structural characteristics of soil and the chemical reactivity of contaminants, the kinetic
5、model of contaminant transportation in the soil-water system is established, in which the immobile water and equilibrium/non-equilibrium adsorption-desorption are considered. The model is solved using the quasi-ADI method; and the corresponding program is used to simulate a paddy field, which is in
6、the pesticide application period. The law of the contaminant transportation and transformation in soil is studied; and the effect of immobile water and non-equilibrium on concentration distribution is analyzed. The result from numerical simulation shows that the concentration hysteresis exists both
7、in immobile water phase and in non-equilibrium phase; so it retards the pollution in prior period, and increases the pollution in latter period. The existences of the two phases both have important effect on the contaminant transportation. So they must be sufficiently considered in numerical simulat
8、ion. The results of numerical simulation and analysis can provide quantitative basis for environmental evaluation, and have important theoretical and practical significance for controlling of contamination. Key words: contaminant; soil-water system; equilibrium/nonequilibrium adsorption/desorption;
9、two-region model; numerical simulation 1 引 言 随着经济的发展和人口的增多， 土壤 -水环境的污染也在不断加重。 为满足经济建设和生活的需要，人们通过施药施肥、填埋垃圾、排放污水、堆放废弃物质等方式将大量污染物质投入到土壤中，造成了土壤 -水环境的严重污染。 目前全球每年农药总投入量已超过 1.80 106t1。大量事实表明，农药的使用已给环境尤其是土壤和地下水造成了极其严重的污染。如：经欧美各国广泛检验，已有 60 多种农药在地下水中被检出。据有关部门统计， 2000 年我国城市垃圾的产量已达 1.5 108t， 处置不当的垃圾长期受雨水淋溶，其淋滤
10、液造成了土壤和地下水的严重污染。此外，煤矿开采造成的土壤和水环境污染也是目前遇到的一个重要的环境问题2。我国仅国有重点煤矿现有矸石山就达 1 000 多座，堆积量约 3 109t，占地约 5.5 107m2，今后我国每年排放矸石约 1.7 108t， 约占土地面积 （ 3 4） 106m23。岩 土 力 学 2008 年 煤矸石被日晒雨淋，经过空气和水的综合作用，发生一系列物理、化学和生物化学变化，其淋滤液严重污染了土壤及其水环境。 人类各种不恰当的活动已经造成了土壤 -水环境的严重污染。避免和治理上述污染的一个重要前提是了解污染物质在土壤 -水环境中运移规律，因此，分析污染物在土壤 -水环境
11、中运移的机制，并用数学模型直观、定量地模拟出来显得尤为重要2。 实践表明，土壤中不可动水和非平衡吸附相是影响土壤 -水环境中污染物浓度分布的重要因素。 因此， 有必要深入研究二者对污染物分布影响的规律。本文在前人研究的基础上，在综合考虑土壤结构、平衡和非平衡吸附解吸以及微生物降解等条件下，建立二维情况下污染物在土壤 -水环境中迁移转化的动力学模型，用拟交替显隐式有限差分法对模型进行离散后，通过编制的程序，以农药为例，定量化模拟污染物浓度分布规律，研究了不可动水和非平衡吸附相对污染物浓度分布的影响。这对于预测预报污染物质迁移归宿、评价环境质量以及环境污染的治理与控制具有理论和现实意义。 2 数学
12、模型的建立 2.1 污染物在土壤中迁移转化的控制方程 污染物在土壤中迁移转化是一个复杂的过程，不同土壤特性和污染物特性都会导致不同运移结果4。实际上土壤具有一定的结构性，且其中运移的农药等溶质一般具有化学反应的特征，势必造成结构土壤中化学反应性溶质运移行为十分复杂。研究表明，不可动水和土壤吸附相的存在是影响田间土壤中污染物迁移转化的重要因素。因此，采用同时考虑可动水和不可动水和吸附 /解吸行为的两区模型进行模拟所得结果将更为精确。本文在该模型中对可动水区考虑非平衡吸附模式，对不可动水仅考虑平衡吸附模式，模型建立如下5 8。 二维情况下考虑非平衡吸附和一阶微生物降解的可动水相的浓度控制方程为 m
13、m m mmmmmmmmmmimmdmmmcmm()()()xx xyyy yxyxCCCDDtx x yCCDDyyqCqCCCxykf K C S C =+ + （ 1）不可动水区的水分虽然不直接参与淋溶，但与可动水区发生溶质交换，并伴随生物降解现象。因此，不可动区总控制方程为 im d immimimcimim dsimim(1)()(1 )KfCCCtCfKC + =（ 2）由质量守恒原理和考虑一级动力学吸附过程可得可动水区非平衡吸附相浓度控制方程为 mmdmm smm()SkKC S St=（ 3） 不可动水区平衡吸附相浓度控制方程为 im d imSKC= （ 4） 式中：mC 和
14、imC 分别为污染物在可动水相和不可动水相中的浓度；mS 和imS 分别为污染物在可动和不可动吸附相上的吸附浓度；m 和im 分别为可动水相和不可动水相的体积含水率；xxD ,xyD ,yyD 和yxD均为弥散系数分量；xq 和yq 为达西流速分量； f 为平衡吸附交换点位占总吸附点位的比例； 为土壤体积密度；dK 为土壤 -水分配系数； k 为一阶吸附解吸速率常数； 为可动水区与不可动水区之间的质量转移系数；cm ，cim ，sm ，sim 分别为可动水相、不可动水相、可动吸附相和不可动吸附相一阶微生物降解速率系数。 2.2 定解条件 溶质运移模型定解条件包括初始条件和边界条件，其常用形式如
15、下： （ 1）初始条件 已知计算区域内的浓度初始分布，有： 0(, ,0) (, ,)Cxyz C xyz= （ 5） （ 2）边界条件 Dirichlet 边界条件： 111(,) (,), ( , )iiiCxt C xt x xyz =（ 6） Neuman 边界条件： 222(,), ( , )ij i i ijCDn qxtxxyzx= = （ 7） Cauchy 边界条件： 33(,)ij i i ijCD vC n q x tx=,3(, )ixxyz=（ 8）3 数学模型的有限差分解 为便于方程的离散，采用如下形式对方程进行简写：记mCu= ，imCv= ，mSS= ，ma =
16、 ，m xxD b = ，mmxy yxD Dc= = ，m yyD d = ，xqe= ，2374第 9 期 梁 冰等：土壤中污染物迁移转化规律的数值模拟研究 yqg= ，mdmcmkf K h +=，mkf p= ， im(1 )kfq+=，im cim sim(1 )fk + =，mdkK = ，msmk +=。 方程（ 1）（ 3）分别简记为 au u u ubdctxxyyxyueuquchvpSyx x y=+ + （ 9） qvuvt =（ 10） Suvt =（ 11） 本文交替显隐式有限差分格式是求解二维和三维方程的一种常用方法，该方法可以将二维和三维方程离散为 2 组或 3
17、 组 3 对角方程组，从而采用追赶法进行最终求解， 避免求解 5 对角或 9 对角矩阵。由于溶质运移方程中混合偏导项的存在，本文在求解模型时对交替显隐式有限差分格式稍加修改，即对混合偏导项始终采用显式格式，给出求解二维溶质运移方程的“拟交替显隐式有限差分格式” ，则方程离散过程如下： 在 到 12 + 时刻对方程（ 9）采用下列差分形式： () ()() ()1/2,14,12 12 12 121, , , 1,12 1212, 12,22,1 , , ,1,12 ,12221, 1 , 1 1, 112,(/2)ij ijijij ij ij ijij ijij ij ij ijij iji
18、j ij ijijuuatuu uubbxxuu uuddyyuuuc + + + +=+,1,1 1,1 ,1 1,112,1, 1 1, 1, 1 1,121, 1, 1 1, 1, 1,1212 121, 1,121,444422ijij i j ij i jijij ij ij ijijij ij ij ijijij ijij iuxyuu uucxyuuuucxyuu uucxyeeuuxx + + + + + +121,1 ,1,1 ,1, , ,22jij ijij ijij ij ij ij ij ijqquuyyhu v pS +（ 12）在 12 + 1 + 时刻对方程（
19、9）采用下列差分形式： () ()() ()11/2,34,12 12 12 121, , , 1,12 1212, 12,2211 11,1 , , ,1,12 ,1222121, 11212,(/2)ij ijijij ij ij ijij ijij ij ij ijij ijij iijuuatuu uubbxxuu uuddyyuuc + + += + +12 12 12,1 1,1 ,112 12 12 12,1 1,1 ,1 1,11212,12 12 12 121, 1 1, 1, 1 1,12,1212 121, 1, 112,12444jijijij i j ij i jij
20、ij ij ij ijijij ijijuuxyuu uucxyuuuucxyuuc + + + + + + + + +12 121, 1, 112 121, 1,12 121, 1,11,1 ,111,1 ,112 12 12 12 12 12, , ,42222ij ijij ijij ijij ijij ijij ij ij ij ij ijuuxyeeuuxxqquuyyhu v pS + + +（ 13）将式（ 10）和（ 11）中 1 + 时段分成同样的 2 个时段进行求解，即 ()()12,14 12 12,12 12, , ,1(/2) 212ij ijij ij ij ij
21、ijij ij ij ijvvquutvv += + （ 14） ()112,34 12 12 1 1,12 12 1 1, ,1(/2) 212ij ijij ij ij ij ijij ij ij ijvvquutvv + += +（ 15） ()()12, 12 12, , ,12 12, , ,1/2 212ij ijij ij ij ijij ij ij ijSSuutSS += +（ 16） ()()112, 12 12 1 1, ,12 12 1 1, ,1/2 212ij ijij ij ij ijij ij ij ijSSuutSS + + += +（ 17） 由式（ 14）
22、和式（ 16）分别整理出1/2,nijv+和1/2,nijS+的表达式，代入式 (12)，求得1/2,niju+，再由1/2,nijv+和1/2,nijS+表达式求得1/2,nijv+和1/2,nijS+。由式（ 15） 、式（ 16）和式（ 13）采用同样2375岩 土 力 学 2008 年 的顺序求出1,niju+，1,nijv+，1,nijS+。循环上述过程即可求得各时刻各相浓度场分布。 4 算 例 本文通过模拟一水田来研究土壤中农药阿特拉津的分布规律以及土壤各相间相互影响关系。研究区域为水平宽度 8 mX = 、垂直深度 4 mY = 的矩形垂直剖面，且处于稳态的含水率分布。该区域初始
23、浓度为 0， 0t = d，在上边界 1m 7mx 处施有浓度301 000 g/mC = 的农药，116t = d 后换以清水。参数取值见表 1，部分模拟结果绘于图 1 9。 表 1 模拟参数取值 Table 1 The values of parameters in the model 模型参数 取值 可动区含水率 m / % 0.250不可动区含水率 im % .050达西流速 yq /-1md 0.063土壤密度 /-3gcm 1.560 可动水降解系数 cm /-1d 0.008不可动水降解系数 cm /-1d .004可动区吸附相降解系数sm/-1d 0.004 不可动区吸附相降解
24、系数 sim /-1d .004纵向弥散度 L /2-1md 0.080横向弥散度 T /2-1md .040一阶吸附解吸速率常数 k /-1d 0.060 土壤 -水分配系数dK/3-1cm g 1.200质量转移系数 0.015平衡吸附交换点位占总吸附点位的比例 f / % 0.200 图 1 8 d 时土壤可动水相浓度等值线（单位 : g/m3） Fig.1 The concentration contours of mobile water phase after being polluted 8 d (unit: g/m3) 0 -1-2-3-4X/m 21 3 4 5 678495
25、 480 460440 420 400 380 365 35030010020020 1 2 5 350300200100201 2 5 Y/ m图 2 8 d 时土壤不可动水相浓度等值线（单位 : g/m3） Fig.2 The concentration contours of immobile water phase after 8 d (unit: g/m3) 图 3 8 d 时土壤可动水区吸附相浓度等值线 （单位 : mg/kg） Fig.3 The concentration contours of adsorption phase in mobile water after be
26、ing polluted 8 d (unit: mg/kg) 0-1-2-3-412 3456 7 8980965950 935920900 880850800400104220060050850400104260020080050980 Y/ mX/ m 图 4 16 d 时土壤可动水相浓度等值线（单位 : g/m3） Fig.4 The concentration contours of mobile water phase after being polluted 16 d (unit: g/m3) 0-1-2-3-457356 78720 700 675 625 65060055040
27、0105120050600550510504002001Y/ m X/ m 1 2 3 4 图 5 16 d 时土壤不可动水相浓度等值线（单位 : g/m3） Fig.5 The concentration contours of immobile water phase after being polluted 16 d (unit: g/m3) 0-1-2-3-4234 5 6 7 8235230225220 215 210 205200 195 19018518015020510.51005050100150180510.5Y/ m1X/ m 图 6 16 d 时土壤可动水区吸附相浓度等
28、值线 （单位 : mg/kg） Fig.6 The concentration contours of adsorption phase in mobile water after being polluted 16 d (unit: mg/kg) 分析各图可知：前 16 d，由于污染物不断从上边界输入，土壤中各点处各相浓度均随时间增加而增加，且污染浓度沿深度方向呈递减变化，即上游各相浓度大于下游各相浓度； 16 d 后，撤去污染源，0-1-2-3-412 3 4 5 6 7 8Y/ m 125120115 110 105 100 95 90 85 60 20 0.5 1 10 5 85 60
29、 20 80 40 10 5 1 0.5 0-1 -2 -3 -4 X/ m Y/ m 1 6 43 5 82 7960 940 920900 875 850 825 800 780 7506001047400960 940 400300750100210 X/ m 2376 第 9 期 梁 冰等：土壤中污染物迁移转化规律的数值模拟研究 污染物随水分向下游淋溶。 由图 7 和图 8 可以看出，淋溶作用使 24 d时下游水相浓度要比上游水相浓度高些，且上、下游浓度均较 16 d 时有所减少，因为可动水相直接参与淋溶，因此，要比不可动水相减少得更快些。由于上游浓度的降低，使得可动区吸附相有部分污染
30、物从土壤中解吸出来，吸附相的浓度较 16 d 时略有减小。而对于下游，虽然水相浓度也比 16 d 时减少了许多， 但因未降低到使污染物从吸附相解吸出来的浓度，因此，下游吸附相的浓度比 16 d 时反而略高些。 Y/ m30 40 5060708090 110130 0-1 -2 -3 -4 1 2 34 5678102050.5122520101120 100 X/ m 图 7 24 d 时土壤可动水相浓度等值线（单位 : g/m3） Fig.7 The concentration contours of mobile water phase after being polluted 24 d
31、 (unit: g/m3) 0-1 -2-3 -435030010050200101210020010250350300360 370380 385 Y/ m1 X/ m 2 3 4 5 6 78图 8 24 d 时土壤不可动水相浓度等值线（单位 : g/m3） Fig.8 The concentration contours of immobile water phase after being polluted 24 d (unit: g/m3) 0-1-2 -3-4 2 3 4 5678209 205 200195 190180520501501002050.5 1 Y/ m1X/ m
32、图 9 24 d时土壤可动水区吸附相浓度等值线 （单位 :mg/kg） Fig.9 The concentration contours of adsorption phase in mobile water after being polluted 24 d (unit: mg/kg) 由上述分析表明：不可动区水相和非平衡吸附相浓度具有滞后现象，在污染前期一定程度上起到了延缓污染作用，而污染后期又在一定程度上加重污染。二者的存在对溶质迁移转化具有重要影响。 5 结 论 本文采用同时考虑可动水和不可动水以及平衡吸附和非平衡吸附的溶质运移模型进行模拟，用拟交替显隐式有限差分格式对模型进行了求解，
33、通过编制的相应程序用该模型对一施药期的水田进行模拟。模拟结果表明：不可动区水相和非平衡吸附相浓度具有滞后现象，在污染前期一定程度上起到了延缓污染作用，而污染后期又在一定程度上加重污染。 二者的存在均对溶质的迁移转化具有重要影响，在模拟中应充分考虑二者的存在。 本文所建立的模型和所编制的程序不仅适用于农药、化肥等污染物迁移规律的数值模拟研究，对于煤矸石淋滤液、垃圾淋滤液等污染物中某些物质在土壤中运移规律的研究也同样适用。 参 考 文 献 1 杨永岗 , 彭钰欣 . 农药与未来农业的发展 J. 环境导论 , 1997, (6): 35 38. YANG Yong-gang, PENG Yu-xin
34、g. Modern agriculture pesticide environment developmentJ. Environment Herald, 1997, (6): 35 38. 2 XUE QIANG, LIANG BING, WANG HUIYUN, et al. Effects of trace elements leaching from coal mine spoil on soil water environment pollution J. Journal of Trace Elements in Biology and Medicine, 2006, 20 (2):
35、 97104. 3 李树志 , 周锦华 , 张怀新 . 矿区生态破坏防治技术 M. 北京 :煤炭工业出版社 , 1998: 46 86. 4 XUE QIANG, LIANG BING, WANG QINXIN. Environmental prediction models for transport of chemicals and pesticides in soils J. Journal of Exper- imental Botany, 2003,54 (Supp. 1): 59 59. 5 刘晓丽 , 梁冰 , 薛强 . 地下水环境中有机污染物迁移转化动力学模型的研究 J. 工程
36、勘测 , 2003, (1): 24 28. LIU Xiao-li, LIANG Bing, XUE Qiang. Organic chemical mathematical model sensitive analysis for parameters groundwater environmentJ. Geotechnical Investigation & Surveying, 2003, (1): 24 28. 6 薛强 , 梁冰 , 刘晓丽 . 有机污染物在土壤中迁移转化的研究进展 J. 土壤与环境 , 2002, 11(1): 90 93. XUE Qiang, LIANG Bing, LIU Xiao-li. Progress on organic contaminant transport and transform in soilJ. Soil and Environmental Sciences, 2002, 11(1): 90 93. 7 李韵珠 , 李保国 . 土壤溶质运移 M. 北京 : 科学出版社 , 1998. 2377
