非饱和带溶质迁移离心模拟的可行性研究

秦欢欢， 喻立珊， 郑春苗

(1.东华理工大学省部共建核资源与环境教育部重点实验室培育基地，江西 南昌 330013；2 北京大学工学院水资源研究中心，北京 100871；3.南方科技大学环境科学与工程学院，广东 深圳 518055)

非饱和带是大气水和地表水同地下水发生联系并进行水分交换的地带，由于污染物大多是随水扩散而发生迁移，因此它是地表污染物向下迁移污染地下水的通道。随着中国经济的发展，一味追求经济增长而不重视环保的弊端逐步显现，各种污染事件时有发生，如2013年山东潍坊地下水污染事件就一度引起公众的关注(Zheng et al., 2013)；山西省太原市工业结构以冶金、化工、煤炭、电力、机械为主，人口密集，人类活动频繁，城市建设发展迅速，由此导致的地下水污染非常严重(蒋方媛等，2009)；近几年来随着大同县社会经济的快速发展，污染物的大量排放，地下水水质呈现出逐渐恶化的趋势(张哲等，2012)；工业固体废物填埋所产生的渗滤液等污水对地下水环境及人类健康存在潜在的威胁(程晓燕等，2016；杨文琳等，2017)。据《全国地下水污染防治规划2011—2020》，全国地下水符合Ⅰ—Ⅲ类水质标准的占63%，Ⅳ—Ⅴ类的占37%，总体情况算好，但对北京、辽宁、吉林、上海、江苏、海南、宁夏和广东等8个省(区、市)641眼井的分析中水质为Ⅳ—Ⅴ类的占73.8%，可见在经济发达区地下水污染状况堪忧。因此，开展非饱和带污染物迁移转化研究，深入掌握其规律，有利于针对性地制定污染预防和治理措施。

非饱和带是岩土颗粒、水、空气三相共存的复杂系统，研究困难极大。近年来，国外对这方面开展研究并做了大量工作，把它当作为研究地下水污染问题的突破口(刘仁冬，2005)，相比之下我国这方面研究开展得较晚。迄今为止，针对溶质在土壤、非饱和带中迁移转化规律进行了大量研究，取得了大量实验资料和理论进展，主要采用三种方法：室内土柱实验、现场实验和数值分析法。但是，这三者都有一定的局限性。土柱实验不能模拟实际土体的重力效应，也不适于做污染物迁移的长期观测，多数土柱实验中装置较小，缺乏与实际情况的类比性；现场实验由于历时较长，受外界因素影响较大，实验本身会造成环境污染，所用示踪剂受到限制，各种边界条件也不易控制(Mattson et al., 2010)；由于缺少足够的数据，数值模型的检验和校正变得非常困难，也很难对污染扩散过程做出准确预测。

1980年代，离心机开始应用于污染物在土壤中的迁移研究。目前，离心模拟技术已被广泛应用于非饱和带溶质迁移研究(Timms et al., 2009; 张建红等，2006a，2006b; Ataie-Ashtiani et al., 2003; Kumar, 2006)。总体来说，离心机可分为土工离心机(大中型)和小型离心机。前者因其操作空间更大，可控制变量更多并能通过技术手段实现动态数据采集，得到了更广泛的应用，模拟的对象也从惰性金属离子(Mckinley et al., 1998)扩展到重金属(张建红等，2004)和非水相流体(胡黎明等，2003；Pantazidou et al., 2000)。近来则主要集中在利用离心机研究非饱和带非均质性对溶质迁移的影响(Menezes et al., 2011)、粘土弱透水层中溶质迁移(Timms et al., 2009)和污染物的去除(Timms et al., 2009; 郝荣福等，2004)。

研究非饱和带溶质迁移对地下水污染防治具有重要意义，而离心模拟技术在低渗透土壤长期污染物迁移及获取数据进行模型验证等方面体现了优越性。然而，目前仍缺乏将离心模拟技术应用于非饱和带溶质迁移规律研究是否可行的系统、全面的论述。因此，本文以非饱和带溶质迁移离心模拟的可行性为目标，通过建立数值模拟，试图在理论上分析、总结该技术在非饱和带各类溶质迁移规律研究中的可行性及适用条件，为离心模拟技术的发展及地下水污染物的防治提供科学的支持。

1 材料与方法

非饱和带溶质迁移和水分迁移密不可分，溶质须溶于水才具备可迁移性，非饱和带水分迁移模型在很大程度上能帮助理解污染物的迁移行为，但不能直接对某种类型污染物的扩散迁移进行分析。地下水中常见污染物有保守性盐类(如NaCl)、挥发性有机物(如BTEX)、非水相流体(石油类)、重金属(如Cd)和放射性核素(如铀)等，理化过程包括微生物过程、吸附、解吸、气化、溶解、化学反应、放射衰变等，涉及水、气、固三相，非常复杂。采用常规实验方法进行研究具有相当大的困难，而且类似NAPLs和重金属等迁移缓慢物质很难在实验室内开展长期观测，离心模拟技术在这类物质的迁移研究中具有较大的优越性。笔者在喻立珊(2015)建立的非饱和带水分迁移离心模型的基础上开发非饱和带溶质迁移离心模拟理论及数值模型，评价将该技术应用于非饱和带溶质迁移研究的可行性。鉴于篇幅，本文只把重点放在溶质型盐类污染物上，考虑的理化过程包括吸附和简单的化学反应等几种。有关离心模拟相关的基础理论，如离心相似比、模拟溶质迁移的相似基础、非饱和带水分迁移离心模拟理论模型等，详见喻立珊(2015)，喻立珊等(2014)等文献。

1.1 非饱和带溶质迁移模型

非饱和带溶质迁移控制方程有多种，主要和目标溶质A的迁移特性有关。

(1)当A为保守性物质，忽略它和土壤之间的吸附与解吸过程，控制方程为：

(1)

式中，cm是溶液浓度，tm是模拟时间，vm是渗流速度，r为半径，Dm是水动力弥散系数。

(2)当A为保守性物质不参与化学反应，假设A与土壤之间的吸附与解吸瞬时达到平衡，控制方程为：

(2)

式中，Rm=1+Kdρb/θm为A的迟滞系数，Kd是吸附系数，ρb为土壤容重，θm为土壤体积含水率，饱和时为土壤孔隙度。

(3)当A为反应性物质时，考虑吸附与解吸作用，假设能瞬时达到平衡，控制方程为：

(3)

式中，cAm，cBm分别为溶质A和产物溶质B的浓度，ka，kb分别为反应速率常数，DAm，DBm分别为水动力弥散系数，RAm和RBm分别为迟滞系数。

1.2 数值计算方法

1.2.1 数值化方案

离心状态下非饱和带溶质迁移的控制方程和常规情况一致，不同之处体现在渗流速度vm。三种情景下非饱和带溶质迁移过程分别由(1)—(3)描述，数值化方案可分为两类：(1)和(2)为一类，(1)可看成是A的迟滞系数Rm为1时(2)的特殊情况；(3)为一类。本文使用算子分裂法对微分方程进行数值化。算子分裂法是众多数值化方法中的一种，是分步骤数值化的方法，Yanenko(1971)将其进行推广应用。它将对流弥散方程分为对流算子和弥散算子，将反应迁移方程分为迁移算子和反应算子。

(1)方程(2)的数值化。 先将方程分裂成对流算子和弥散算子：

(4)

将对流算子显式向前差分：

(5)

(6)

进一步整理：

(7)

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

(2)方程(3)的数值化。将方程分裂成迁移算子和反应算子：

(8)

将迁移算子用隐式差分进行离散：

(9)

(10)

式中

(11)

1.2.2 初始和边界条件

鉴于本文研究地下水污染，故假定离心模型和原型的土在初始时刻不含溶质，即初始条件为：

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

(12)

考虑到现实情况下溶质污染地下水的一般场景为：含有污染物的垃圾散布在地表，在一个雨天，垃圾中的溶质淋溶而出并随雨水下渗污染土壤及地下水。因此，本文考虑的边界条件为：一维土柱顶部提供一个短暂的溶液供给，溶质A的浓度先维持不变后转为0 g/m3；上边界水力条件从固定基质吸力为0 kPa转变为水流通量为0 m/s的边界；模型底部为自由流边界，基质吸力固定为0 kPa，迁移模拟过程中处理成弥散通量固定的边界条件。用数学表达为：

(13)

式中，c0为初始浓度，tpulse为降雨事件时长。方程(13)中第三式可离散为：

(14)

2 结果与分析

为方便和原型进行比较，笔者利用Visual Basic编写了自然条件下非饱和带溶质迁移数值模拟代码，参照前文阐述的数值化方案，对自然条件及离心状态下的溶质迁移过程进行计算。

图2 非饱和带保守性溶质迁移离心模拟可行性比较(图中标注均为原型时间)Fig.2 Feasibility comparison of centrifugal modeling of nonreactive solute transport in unsaturated zone (time in the figure is the prototype time)

2.1 保守性溶质

由于保守性溶质既不发生吸附也不发生化学反应，从离心相似理论上看，离心模型能够正确反映原型下溶质迁移行为而被广泛应用于非饱和带溶质迁移机理的研究(Nakajima et al., 1998; Kumar, 2006)。此处利用数值模拟结果的比较考察到底是否可行。数值模拟的固定参数及值为：θ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模型参数(喻立珊，2015)，Ks是土壤的饱和水力传导系数，rb为离心模型底部半径。离心模型和原型空间离散节点数n= 101，A浓度为100 g/m3。原型的时长共100天，下雨时长tpulse为75天，原型非饱和带厚度为2 m，其它变动参数的取值列于表1，离心模型和原型模拟结果对比如图2所示。

表1 保守性溶质情况下数值模拟参数

注：D为水动力弥散系数；θ0是初始含水率；N是重力水平；表示离心模型是在N个g的加速度环境下开展实验；L是模型高度；Δt是时间步长；LX表示离心模型；YX表示原型。

从图2a可看到，离心状态下非饱和带保守性溶质的迁移和原型基本一致，但在靠近底部有一定偏差，此处模型中的迁移要稍滞后于原型。由于离心加速度存在径向分布，模型2/3高度处至底部的这一部分相似比被低估，空间尺寸被缩小，导致水分及溶质迁移均有一定的滞后。另外，N为40比20时底部的偏差要小，说明同一设备使用更大的N值有利于减轻这种偏差。图2b展现了分子弥散系数对溶质迁移的影响：分子弥散系数越大，溶质的污染范围更大，但最高浓度峰值会变小。这是因为分子弥散系数越大，它对溶质浓度剖面曲线的“削峰填谷”的能力越强。图2c比较了不同初始含水率对溶质迁移的影响，初始含水率越低，土壤对溶质的迁移阻碍越大。由于孔隙度和渗透系数之间关系不太明确且太过复杂，故没有考察土壤压实度对溶质迁移的影响。总体来看，离心模拟能较好地反映保守性溶质在原型下的迁移行为，而要减小离心加速度分布不均产生的误差，可通过严格计算不同部位的重力水平来正确处理模型尺寸的缩放从而避免这种误差，或者使用更大的N值来尽力减小这种误差。

图3 非饱和带吸附性溶质迁移离心实验模拟可行性比较(图中标注均为原型时间)Fig.3 Feasibility comparison of centrifugal modeling of adsorption solute transport in unsaturated zone (time in the figure is the prototype time)

2.2 吸附性溶质

在非饱和带保守性溶质迁移基础之上，考虑吸附作用的影响。数值模拟的固定参数及值为：θ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，离心模型和原型空间离散节点数n=101，初始含水率为0.2(均匀分布)，A浓度为100 g/m3。原型的时长共100天，下雨时长tpulse为75天，原型非饱和带厚度为2 m，其它变动参数的取值列于表2，离心模型和原型模拟结果对比如图3所示。

表2 吸附性溶质情况下数值模拟参数

注：N是重力水平，表示离心模型是在N个g的加速度环境下开展实验；L是模型高度；Δt是时间步长；Kd是吸附系数；LX表示离心模型；YX表示原型。

从图3a可以看出，离心模拟的结果和原型吻合度较好，说明该技术用于非饱和带瞬时平衡吸附性溶质的迁移研究是可行的。和图2a相比较，可知吸附作用的存在会明显阻滞溶质的迁移扩散，对于控制污染的范围具有积极意义。但溶质会长期残留在土壤中，在农业上对作物有影响。图3b比较了分配系数对溶质迁移的影响，分配系数越大，土壤对溶质的迟滞作用越强，导致溶质迁移得更加缓慢，污染范围也更加集中。

比较遗憾，本文只讨论了平衡吸附，而实际情况下非平衡吸附更加普遍。对于非平衡吸附，目前常用双区模型来描述其迁移过程，如Li等(1993)，离心模拟技术在非平衡吸附下非饱和带溶质的迁移研究中也是可行的。

图4 非饱和带反应性溶质迁移离心实验模拟可行性比较(图中标注均为原型时间)Fig.4 Feasibility comparison of centrifugal modeling of reactive solute transport in unsaturated zone (time in the figure is the prototype time)

2.3 反应性溶质

研究非饱和带伴随产物溶质B出现的反应性溶质A的迁移过程，着重考察B的性质对A迁移过程的影响。数值模拟的固定参数及值为：θ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，离心模型和原型空间离散节点数n=101，初始含水率为0.2(均匀分布)，A浓度为100 g/m3。原型的时长共100天，下雨时长tpulse为75天，原型非饱和带厚度为2 m，其它变动参数的取值列于表3，离心模型和原型模拟结果对比如图4。

表3 反应性溶质情况下数值模拟参数

注：N是重力水平，表示离心模型是在N个g的加速度环境下开展实验；L是模型高度；Δt是时间步长；ka和kb为A和B的反应速率常数，YXF代表考虑反应性溶质的原型，LXF代表考虑反应性溶质的离心模型。

从图4a，e可看出，化学反应的存在会使得离心模拟的结果偏离原型。和放射性核素一样，当反应进行得比较慢时，离心模拟A的迁移过程能较好地和原型匹配；但若反应较为迅速，则误差偏大，往往会高估A的迁移扩散能力和污染区A的浓度。再分别对比图4b，c，f，g，可知对于产物溶质B，离心模拟的结果和原型相差很大，浓度的预测值误差能达到2个数量级。进一步仔细观察可发现原型中土壤顶部附近区域A的浓度分布有明显的陡降特征，而B在顶部区域分布集中，这些特征在离心模型中无法体现，表明离心模拟技术用于反应较快的非饱和带溶质迁移研究时很难抓住过程的细节。

从模拟结果可以发现，离心模拟技术对于反应过程相对较慢的溶质A的迁移过程预测能够较好地和原型匹配，但当反应过程较快时，离心模拟将产生较大的误差，对于迁移过程的细节捕捉能力也不强。此外，对B的预测，离心模型和原型相差甚远，浓度误差能达到2个数量级，表明离心模拟技术在对B生成和迁移过程的预测上是失败的。若B也是研究的重点，那么离心模拟技术是不值得推荐的。反之，若B在研究中不那么重要，那么在适当的条件(反应进程不快)下，离心模拟技术可行。

3 结论

在水分迁移模型的基础上，建立了离心场内一维非饱和带溶质迁移模型及数值化方案，对离心模拟技术应用于三种类型溶质在非饱和带迁移研究的可行性进行了研究，结果表明离心模拟技术在适当条件下用于非饱和带溶质迁移规律研究是可行的。具体结论如下：

(1)离心模拟可用于保守性溶质的迁移研究。离心模型能准确反映原型下这类溶质的迁移过程，但离心加速度的径向分布会使得靠近底部区域的溶质迁移稍滞后于原型，使用更大的重力水平N值，能使得加速度分布更加均匀从而减小误差。

(2)分子弥散系数和土壤的初始含水率对溶质迁移有影响。分子弥散系数越大，对溶质浓度剖面曲线的“削峰填谷”能力越强，使得溶质分布更加均匀，导致污染范围更大，而浓度峰值会变小。土壤的初始含水率越低，溶质的迁移扩散越慢。

(3)离心模拟适用于吸附瞬时平衡溶质迁移规律研究。吸附与解吸作用的存在对溶质迁移有明显影响，分配系数越大，土壤对溶质的迟滞作用越明显，对于控制污染范围具有积极意义。

(4)对于有产物溶质B生成的反应性溶质A在非饱和带的迁移，当反应较慢时，离心模型预测的A的迁移过程能较好地匹配原型，但对B的预测误差较大，浓度误差可达2个数量级。若B是研究重点，则不推荐离心模拟技术；若B不那么重要，则在适当条件(反应进程不快)下，离心模拟技术是可行的。

参考文献

程晓燕，江宝锋，刘传杰. 2016.危险废物填埋场对地下水环境影响分析[J].东华理工大学学报：自然科学版，39(增刊)：186-189.

郝荣福，胡黎明，邢巍巍.2004.土壤中可挥发性污染物清除的离心试验研究[J].岩土力学， 25(7)：1037-1040.

胡黎明，郝荣福，殷昆亭，等.2003.BTEX在非饱和土和地下水系统中迁移的试验研究[J]．清华大学学报：自然科学版， 43(11)：1546-1549.

蒋方媛，郭清海. 2009.城市地下水污染成因分析——以山西省太原市为例[J].东华理工大学学报：自然科学版，32 (1)：82-88.

刘仁冬. 2005.非饱和带中水分迁移和污染物运移的数值模拟[D].沈阳：东北大学.

杨文琳，张静，袁野，等. 2017.基于有限差分法的地下水溶质运移模拟：以某垃圾填埋场为例[J]. 环境工程，35(12)：30-35.

喻立珊，曹国亮，许模，等. 2014.离心实验在污染物迁移研究中的应用[J].地球科学进展， 29(2)：227-237.

喻立珊. 2015.非饱和带溶质迁移离心实验的数值模拟研究[D].北京：北京大学.

张建红，胡黎明. 2006a.重金属离子和LNAPLs在非饱和土中的运移规律研究[J].岩土工程学报， 28(2)：276-280.

张建红，吕禾，王文成. 2006b.铜离子在非饱和土中迁移的离心模型试验研究[J].岩土力学， 27(11)：1885-18 90.

张建红，严冬. 2004.非饱和粉质砂土中铜离子迁移的离心模型试验研究[J].岩土工程学报， 26(6)：792-797.

张哲，郑秀清，陈军峰，等. 2012.利用灰色聚类法评价大同县地下水水质[J].东华理工大学学报：自然科学版，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.