当溶质在土／水界面达到平衡之后，会存在以下关系：

Ft(j)?F?(j)?e?(j)/?0 （18）

式中：Ft(j) — 土壤中溶质的总量在单位质量土粒上的质量分数，包括土粒

上吸附溶质和孔隙溶液中的溶质；

F?(j) — 土颗粒表面所吸附的溶质（j）的质量分数；

?(j) — 孔隙中的溶质的质量浓度（ｇ／Ｌ）； ｅ — 土壤孔隙比；

?0 — 土样干密度（g/cm3）

当溶质在土粒上的吸附呈线性状况时，沉淀满足Henry定律：

Ft(j)?vd?(j)． （19）

式中：vd — 溶质在土水界面的分配系数．此时，满足下式：

Ft(j)?(A?B)?(j) （20）

式中：Ａ、Ｂ 均为模型参数，可表示为A?(e?b??s)/s?b,B???b

? — 土壤空隙率；

S — 土壤容重． ．

对于溶质在土壤介质上的吸附，其分配系数通常情况下并非常数，而是随溶质浓度变化而变化的参数．为了便于计算，本文对溶质在土柱中土水间的分配关系做出如下简化假定：

Ft(j)??(?(j))n， （21）

传统的运移控制方程如下：

??(j)???(j)??(j) （22） Rd?(D(?(j)))?v?t?x?x?x

式中：ｔ为时间，Rd为阻滞因子，D(?(j)）表示随土壤孔隙中溶质的质量浓度 ?(j)而变化的非线性的扩散系数，v为孔隙流体的流速．

式（22）是以土壤孔隙中溶质的平衡浓度为基础建立起来的．假定扩散系数Ｄ 为常数，考虑吸附的非线性特性，将式（21）代入式（22），可得改进的非线性重金属运移控制方程为

1?n?F(j)1?n?F(j)?F(j)?(?bn(?)1n/e)t?(DFt(j)nt)?vFt(j)nt （23）

?t?x?x?x

13

?Ft0(j)??(?0(j))n其中：???Ft(j)?0??xFt0(j)?(上边界)(下边界)

基于土粒上溶质质量分数的运移参数提取方法，根据有限差分法，可将运移方程（22）

写成如下差分格式：

RddFi(j)F(j)?Fi(j)F?Fi?1vFi?1(j)?Fi(j)2D?(Fi?1/2(j)i?1?Fi?1/2(j)i)?(dtxi?1?xi?1xi?1?xixi?xi?12xi?1?xiF(j)?Fi?1(j)?i)xi?xi?1 （23）

式中：xi — Ｘ 方向被差分之后的差分网格节点的空间位置；

Fi(j)≈F(j)，ｉ＝１，…，Ｎ－１；F0(j)、FN(j)为上下边界处溶质的质量分

数；Fi0(j)为0时刻土柱内溶质的质量分数．

上述差分格式可以利用Matlab软件进行求解，结合最小二乘法，利用内嵌的最优化工具箱求解运移参数．对比拟合结果与试验得到的浓度剖面数据，当方差最小时即认为此时的参数即为有代表性的典型参数．

由模型可得重金属在土壤中的传播模型如下图所示：

在图中可以清晰地看出重金属在土壤中的传播随时间的增减而减缓，而且速度较水体传播慢。

3. 大气扩散模型的建立及求解：

大气中重金属主要来源于工业生产、汽车尾气排放、汽车轮胎磨损产生的大量含有重金属的有害气体和粉尘。它们主要分布在工矿的周围及公路、铁路两侧。

重金属在大气中的传播一部分通过自然沉降和雨淋沉降进入土壤圈和水体中，一部

14

分，直接污染空间。

在这里因为数据不足更为了简化模型，我们忽略雨淋沉降，将通过大气沉降地表的重金属看做由自然沉降产生的。而自然沉降的主要作用是重力作用。

在垂直方向，考虑粒子受到重力、浮力、流体阻力和随机力的作用，建立粒子的运动方程如下：

dwG?F浮?F阻???(t) （24） dtm其中w为粒子垂直方向速度，F浮为粒子受到的浮力，F阻为粒子受到的流体阻力，

?(t)为随机力及大气流动引起的随机加速度，满足如下条件:

?(t)?0,?(t)?(t')?K?(t?t?) (25)

粒子受到的阻力为F阻?6????，其中?为动力粘性系数，r为粒子半径。 由式（24）得：

?gdw9??g(1?)?w??(t) 其中?g和?p分别为气体和粒2dt?p2?pr2?pr29?，由上式得到：

?t子密度，通常?g《?p。令??dw1 ?g?w??(t) 解此式可得：

dt?ts?tw(t)?g?(1?e?)?e?t?w(0)??e0??(s)ds （26）

将 w(t)分解成为平均量和脉动量之和： w(t)?w(t)?w?(t)

由这两式可得：

w(t)?g?(1?e?t?)?e??t?w(0) （27）

ts?tw?(t)?e?tw(0)??e0??(s)ds （28）

由上式求出w(t)的方差w?2(t)，由于w(t)和?(t)不相关，则，

w?2(t)?w?2(0)e?2t?)?2K?(1?e?t?2t?) （29）

由（28）式可得： w?(0)w?(t)?w?2(0)e?

2假设w(t)是平稳随机过程，则：w?2(0)?w?2(t)??w 2?w由上面三式可得：2K??? 即：K?

2?2w由（26）式可得：w(t??t)?g?(1?e?(t??t)?)?e?(t??t)?w(0)??(t??t)0es?t??t??(s)ds

15

?tts?t w(0)?et??w(t)?g?(1?e?)??e??(s)ds?

??0?? 讲下式代入上式得：w(t??t)?g?(1?e?(t??t)?)?e?(t??t)?w(t)??(t??t)0es?t??t??(s)ds

将最后一项离散化得：w(t??t)?g?(1?e令t=0得： w(?t)?g?(1?e 方差为：w?2(?t)?e?2?t??t??t?)?e??t?w(t)??? （30）

?)?e??t?w(0)???

?w?2(0)??2 与（29）式比较可得：

??2K?(1?e2?2?t?)??(1?e2w?2?t?) )?e??t代入式（30）得到粒子运动的离散方程：w(t??t)?g?(1?e??t??w(t)??w(1?e??t?0.5)?

在水平x、y方向不考虑重力作用，则在平均风速为 (u，v，w)下粒子运动的离散型方程为：

??? ????u(t??t)?u?ev(t??t)?v?e??t??t??u(t)??u(1?ev(t)??v(1?e??t??t??t?0.5)?

??t??t?0.5)??0.5w(t??t)?w?g?(1?e?)?e?w(t)??w(1?e)?有题可知，结合计算精度和时间的考虑，时间步长?t?(3—5)?r?0.1?m 由公式：

VD?C(t1)HIn,其中，C(t1)为t1时刻的平均浓度，C(t2)为t2时刻的平均浓度。 t2?t1C(t2)由于，在空气中传播的影响因素很多，而且其速度比水体快。

结合上面三个模型可知，在同样的条件下，金属污染源在空气中传播最快，其次是在水中传播，最后在土壤中传播。而且这三种传播状态是相互转化的，对于同一污染源来说，对其进行时间和空间的拟合得出如下传播范围

（图五）所示：

由图可知，污染源一定在三者叠加的正中心处，再结合八种重金属元素在该区的分布图，找出各金属的污染源。 如下图

16