土壤中污染物迁移转化规律的数值模拟研究

1、第29卷第9期 岩 土 力 学 V ol.29 No. 9 2008年9月 Rock and Soil Mechanics Sep. 2008收稿日期:2006-12-26文章编号:1000-7598-(2008 09-2373-05土壤中污染物迁移转化规律的数值模拟研究梁 冰1,赵 颖1,2,薛 强2,刘 磊1,2(1. 辽宁工程技术大学 力学与工程科学系,辽宁 阜新123000;2. 中国科学院武汉岩土力学研究所,武汉 430071摘 要:考虑到田间土壤的结构性以及运移其中的污染物的化学反应特征,建立了两区双点平衡/动力学吸附溶质运移模型。用拟交替显隐式有限差分格式对模型进行了离散,通过编

2、制的相应程序,对一施药期的水田进行了模拟。研究了土壤中污染物的迁移转化规律,分析了不可动水和非平衡吸附相对浓度分布的影响。结果表明:不可动区水相和非平衡吸附相的浓度具有的滞后现象,在污染前期一定程度上起到了延缓污染作用,而在污染后期又在一定程度上加重了污染。二者的存在均对溶质的迁移转化具有重要的影响,进行数值模拟时应充分考虑二者的存在。模拟及分析结果为土壤环境质量评价提供了定量依据,同时对环境污染的治理与控制具有理论和现实意义。关 键 词:污染物;土壤-水环境;平衡/非平衡双点吸附/解吸;两区模型;数值模拟 中图分类号:TU 441 文献标识码:AResearch on numerical s

3、imulation of contaminanttransportation and transformation in soilLIANG Bing 1, ZHAO Ying 1,2, XUE Qiang 2, LIU Lei 1,2(1. Department of Mechanics and Engineering Sciences ,Liaoning Technical University, Fuxin 123000, China;2. Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 4

4、30071, ChinaAbstract: Considering the structural characteristics of soil and the chemical reactivity of contaminants, the kinetic model of contaminant transportation in the soil-water system is established, in which the immobile water and equilibrium/non-equilibrium adsorption-desorption are conside

5、red. The model is solved using the quasi-ADI method; and the corresponding program is used to simulate a paddy field, which is in 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

6、 concentration distribution is analyzed. The result from numerical simulation shows that the concentration hysteresis exists both 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

7、 phases both have important effect on the contaminant transportation. So they must be sufficiently considered in numerical simulation. The results of numerical simulation and analysis can provide quantitative basis for environmental evaluation, and have important theoretical and practical significan

8、ce for controlling of contamination.Key words: contaminant; soil-water system; equilibrium/nonequilibrium adsorption/desorption; two-region model; numerical simulation1 引 言随着经济的发展和人口的增多,土壤-水环境的污染也在不断加重。为满足经济建设和生活的需要,人们通过施药施肥、填埋垃圾、排放污水、堆放废弃物质等方式将大量污染物质投入到土壤中,造成了土壤-水环境的严重污染。目前全球每年农药总投入量已超过1.80 ×

9、106 t 1。大量事实表明,农药的使用已给环境尤其是土壤和地下水造成了极其严重的污染。如:经欧美各国广泛检验,已有60多种农药在地下水中被检出。据有关部门统计,2000年我国城市垃圾的产量已达1.5 × 108 t ,处置不当的垃圾长期受雨水淋溶,其淋滤液造成了土壤和地下水的严重污染。此外,煤矿开采造成的土壤和水环境污染也是目前遇到的一个重要的环境问题2。我国仅国有重点煤矿现有矸石山就达1 000多座,堆积量约3 × 109 t ,占地约 5.5 × 107 m 2,今后我国每年排放矸石约1.7× 108 t ,约占土地面积(34×106 m

10、 2 3。岩 土 力 学 2008年煤矸石被日晒雨淋,经过空气和水的综合作用,发生一系列物理、化学和生物化学变化,其淋滤液严重污染了土壤及其水环境。人类各种不恰当的活动已经造成了土壤-水环境的严重污染。避免和治理上述污染的一个重要前提是了解污染物质在土壤-水环境中运移规律,因此,分析污染物在土壤-水环境中运移的机制,并用数学模型直观、定量地模拟出来显得尤为重要2。实践表明,土壤中不可动水和非平衡吸附相是影响土壤-水环境中污染物浓度分布的重要因素。因此,有必要深入研究二者对污染物分布影响的规律。本文在前人研究的基础上,在综合考虑土壤结构、平衡和非平衡吸附解吸以及微生物降解等条件下,建立二维情况下

11、污染物在土壤-水环境中迁移转化的动力学模型,用拟交替显隐式有限差分法对模型进行离散后,通过编制的程序,以农药为例,定量化模拟污染物浓度分布规律,研究了不可动水和非平衡吸附相对污染物浓度分布的影响。这对于预测预报污染物质迁移归宿、评价环境质量以及环境污染的治理与控制具有理论和现实意义。2 数学模型的建立2.1 污染物在土壤中迁移转化的控制方程污染物在土壤中迁移转化是一个复杂的过程,不同土壤特性和污染物特性都会导致不同运移结果4。实际上土壤具有一定的结构性,且其中运移的农药等溶质一般具有化学反应的特征,势必造成结构土壤中化学反应性溶质运移行为十分复杂。研究表明,不可动水和土壤吸附相的存在是影响田间

12、土壤中污染物迁移转化的重要因素。因此,采用同时考虑可动水和不可动水和吸附/解吸行为的两区模型进行模拟所得结果将更为精确。本文在该模型中对可动水区考虑非平衡吸附模式,对不可动水仅考虑平衡吸附模式,模型建立如下5-8。二维情况下考虑非平衡吸附和一阶微生物降解的可动水相的浓度控制方程为m m m m m m m m m m mm m im m d m m m cm m (xx xy yy yxy x C C C D D t x x y C C D D y y x q C q C C C x ykf K C S C µ=+ (1 不可动水区的水分虽然不直接参与淋溶,但与可动水区发生溶质交换,

13、并伴随生物降解现象。因此,不可动区总控制方程为im d imm im im cim im d sim im (1(1K f C C C tC f K C µµ+= (2由质量守恒原理和考虑一级动力学吸附过程可得可动水区非平衡吸附相浓度控制方程为mm d m m sm m (S k K C S S tµ= (3 不可动水区平衡吸附相浓度控制方程为im d im S K C = (4式中:m C 和im C 分别为污染物在可动水相和不可动水相中的浓度;m S 和im S 分别为污染物在可动和不可动吸附相上的吸附浓度;m 和im 分别为可动水相和不可动水相的体积含水率;

14、xx D ,xy D ,yy D 和yx D 均为弥散系数分量;x q 和y q 为达西流速分量;f 为平衡吸附交换点位占总吸附点位的比例;为土壤体积密度;d K 为土壤-水分配系数;k 为一阶吸附解吸速率常数;为可动水区与不可动水区之间的质量转移系数;cm µ,cim µ,sm µ,sim µ分别为可动水相、不可动水相、可动吸附相和不可动吸附相一阶微生物降解速率系数。 2.2 定解条件溶质运移模型定解条件包括初始条件和边界条件,其常用形式如下:(1初始条件已知计算区域内的浓度初始分布,有:0(,0(,C x y z C x y z = (5(2边界条件

15、Dirichlet 边界条件:111(,(, (,i i i C x t C x t x x y z =(6Neuman 边界条件:222(, (,ij ii i j C D n q x t x x y z x =(7Cauchy 边界条件:33(,iji i i j CD v C n q x t x =,3(,i x x y z =(83 数学模型的有限差分解为便于方程的离散,采用如下形式对方程进行简写:记m C u =,im C v =,m S S =,m a =,m xx D b =,m m xy yx D D c =,m yy D d =,x q e =,2374第9期 梁 冰等:土壤

16、中污染物迁移转化规律的数值模拟研究y q g =,m d m cm kf K h µ+=,m kf p =,im (1k f q +=,im cim sim (1f k µµ+=,m d kK =,m sm k µ+=。方程(1(3分别简记为au u u u b d c t x x y y x y u eu quc hu v pS y x x y=+(9qvu v t= (10 Su v t= (11 本文交替显隐式有限差分格式是求解二维和三维方程的一种常用方法,该方法可以将二维和三维方程离散为2组或3组3对角方程组,从而采用追赶法进行最终求解,避免求解

17、5对角或9对角矩阵。由于溶质运移方程中混合偏导项的存在,本文在求解模型时对交替显隐式有限差分格式稍加修改,即对混合偏导项始终采用显式格式,给出求解二维溶质运移方程的“拟交替显隐式有限差分格式”,则方程离散过程如下:在到12+时刻对方程(9采用下列差分形式:(1/2,14,121212121,1,121212,12,22,1,1,12,12221,1,11,112,(/2i j i ji j i j i j i ji j i ji ji j i j i j i j i j i j i j i j i j i j u u a t u u u u b b x x u u u u d d y y u

18、u u c +=+,1,11,1,11,112,1,11,1,11,121,1,11,1,1,1212121,1,121,444422i j i j i j i j i j i j i j i j i j i ji j i j i j i j i j i j i j i j i j i u x yu u u u c x yu u u u c x yu u u u c x ye e u u x x +121,1,1,1,1,22j i j i j i j i j i j i j i j i j i j i jq q u u y y h u v p S +(12在12+1+时刻对方程(9采用下列差

19、分形式: (11/2,34,121212121,1,121212,12,221111,1,1,12,1222121,11212,(/2i j i ji ji j i j i j i j i j i j i j i j i j i j i j i j i j i i j u u a t u u u u b b x x u u u u d d y y u u c +=+121212,11,1,112121212,11,1,11,11212,121212121,11,1,11,12,1212121,1,112,12444j i j i j i j i j i j i j i j i j i j i

20、j i ji j i j i j i j u u x yu u u u c x yu u u u c x yu u c +12121,1,112121,1,12121,1,11,1,111i j u u x ye e u u x x q q u u y y h u v p S +(13将式(10和(11中1+时段分成同样的2个时段进行求解,即(12,141212,1212,1(/221 2i ji j i ji j i j i j i ji j i j i j i jv v q u u t v v +=+ (14 (112,34121211,121211,1(/221 2i j i ji ji

21、 j i j i j i j i j i j i j i j v v q u u t v v +=+ (15 (12,1212,1212,1/221 2i ji j i j i j i j i j i j i j i j i j S S u u t S S +=+ (16(112,121211,121211,1/221 2i j i ji j i j i j i j i j i j i j i j S S u u t S S +=+ (17由式(14和式(16分别整理出1/2,n i jv +和1/2,n i j S +的表达式,代入式(12,求得1/2,n i j u +,再由1/2,n i

22、 j v +和1/2,n i j S +表达式求得1/2,n i j v +和1/2,n i jS +。由式(15、式(16和式(13采用同样2375岩 土 力 学 2008年的顺序求出1,n i j u +,1,n i j v +,1,n i j S +。循环上述过程即可求得各时刻各相浓度场分布。4 算 例本文通过模拟一水田来研究土壤中农药阿特拉津的分布规律以及土壤各相间相互影响关系。研究区域为水平宽度8 m X =、垂直深度 4 m Y =的矩形垂直剖面,且处于稳态的含水率分布。该区域初始浓度为0,0t = d ,在上边界1m 7m x <<处施有浓度30 1 000 g/m

23、C =的农药,116t =d 后换以清水。参数取值见表1,部分模拟结果绘于图19。表1 模拟参数取值Table 1 The values of parameters in the model模型参数取值可动区含水率m /% 0.250不可动区含水率im % 0.050 达西流速y q /-1m d 0.063 土壤密度/-3g cm 1.560 可动水降解系数cm µ/-1d 0.008不可动水降解系数cm µ/-1d 0.004可动区吸附相降解系数sm µ/-1d 0.004不可动区吸附相降解系数sim µ/-1d 0.004 纵向弥散度L /2-1m

24、 d 0.080横向弥散度T /2-1m d 0.040一阶吸附解吸速率常数k /-1d 0.060土壤-水分配系数d K /3-1cm g 1.200质量转移系数 0.015 平衡吸附交换点位占总吸附点位的比例f / % 0.200图1 8 d 时土壤可动水相浓度等值线(单位: g/m 3Fig.1 The concentration contours of mobile water phase after being polluted 8 d (unit: g/m 30 -1-2-3-4X /m21 3 4 5 678495 480460440420400 380365350300 100

25、2002012 5 350300200 100 20125 Y / m图2 8 d 时土壤不可动水相浓度等值线(单位: g/m 3Fig.2 The concentration contours of immobilewater phase after 8 d (unit: g/m 3图3 8 d 时土壤可动水区吸附相浓度等值线(单位: mg/kg Fig.3 The concentration contours of adsorption phase in mobile water after being polluted 8 d (unit: mg/kg0-1-2-3-4123456789

26、80965950935920900880850800400104220060050850400104260020080050980Y / mX / m 图4 16 d 时土壤可动水相浓度等值线(单位: g/m 3Fig.4 The concentration contours of mobile water phase after being polluted 16 d (unit: g/m 3 0-1-2-3-45735678720700 675625 650600550400105120050600550510504002001Y / mX / m 123 4图5 16 d 时土壤不可动水

27、相浓度等值线(单位: g/m 3Fig.5 The concentration contours of immobile water phase after being polluted 16 d (unit: g/m 30-1-2-3-42345678235230225220215 210 205200 19519018518015020510.51005050100150180510.5Y / m1X / m图6 16 d 时土壤可动水区吸附相浓度等值线(单位: mg/kg Fig.6 The concentration contours of adsorption phase in mo

28、bile water after being polluted 16 d (unit: mg/kg分析各图可知:前16 d ,由于污染物不断从上边界输入,土壤中各点处各相浓度均随时间增加而增加,且污染浓度沿深度方向呈递减变化,即上游各相浓度大于下游各相浓度;16 d 后,撤去污染源,0-1-2-3-4123 4 5 6 78Y / m125120115110105 100 95908560 200.5 110 585 60 20 80 40 10 5 1 0.50-1 -2-3-4 X /mY / m1 6 43 5 82 7960 940 920900875850 825800 780 75

29、0600 10 47400960940 400300 750 100 210 X / m2376第9期 梁 冰等:土壤中污染物迁移转化规律的数值模拟研究污染物随水分向下游淋溶。由图7和图8可以看出,淋溶作用使24 d 时下游水相浓度要比上游水相浓度高些,且上、下游浓度均较16 d 时有所减少,因为可动水相直接参与淋溶,因此,要比不可动水相减少得更快些。由于上游浓度的降低,使得可动区吸附相有部分污染物从土壤中解吸出来,吸附相的浓度较16 d 时略有减小。而对于下游,虽然水相浓度也比16 d 时减少了许多,但因未降低到使污染物从吸附相解吸出来的浓度,因此,下游吸附相的浓度比16 d 时反而略高些。

30、Y / m30405060708090 1101300-1-2 -3 -41 2 345678102050.5122520101120100 X / m图7 24 d 时土壤可动水相浓度等值线(单位: g/m 3Fig.7 The concentration contours of mobile water phaseafter being polluted 24 d (unit: g/m 30-1 -2-3-435030010050200101210020010250350300360370380385Y / m1 X / m 2345 6 78图8 24 d 时土壤不可动水相浓度等值线(单

31、位: g/m 3Fig.8 The concentration contours of immobile water phaseafter being polluted 24 d (unit: g/m 30-1-2 -3-42 3 4 5678209205200195190180520501501002050.51 Y / m1X / m图9 24 d 时土壤可动水区吸附相浓度等值线(单位:mg/kg Fig.9 The concentration contours of adsorption phase in mobile water after being polluted 24 d (u

32、nit: mg/kg由上述分析表明:不可动区水相和非平衡吸附相浓度具有滞后现象,在污染前期一定程度上起到了延缓污染作用,而污染后期又在一定程度上加重污染。二者的存在对溶质迁移转化具有重要影响。 5 结 论本文采用同时考虑可动水和不可动水以及平衡吸附和非平衡吸附的溶质运移模型进行模拟,用拟交替显隐式有限差分格式对模型进行了求解,通过编制的相应程序用该模型对一施药期的水田进行模拟。模拟结果表明:不可动区水相和非平衡吸附相浓度具有滞后现象,在污染前期一定程度上起到了延缓污染作用,而污染后期又在一定程度上加重污染。二者的存在均对溶质的迁移转化具有重要影响,在模拟中应充分考虑二者的存在。本文所建立的模型和所编制的程序不仅适用于农药、化肥等污染物迁移规律的数值模拟研究,对于煤矸石淋滤液、垃圾淋滤液等污染物中某些物质在土壤中运移规律的研究也同样适用。参 考 文 献1 杨永岗, 彭钰欣. 农药与未来农业的发展J. 环境导论,1997, (6: 35-38.YANG Y