数学建模论文——重金属污染

城市表层土壤重金属污染分析模型

摘要

本文针对城市地质演变受人类活动影响所引发的土壤重金属污染问题，依据附件中给出的数据，先用surfer软件绘制出各主要重金属元素的空间分布三维曲面图，接着通过单因子指数法、内梅罗综合污染指数法与潜在生态危害指数法，对不同区域内的重金属污染程度进行分析；再通过变异系数法，结合实际，找出污染产生的主要原因；然后结合土壤溶质运移的基本理论，运用多元非线性拟合，建立模型，确定出污染源的位置及其坐标；最后将模型由二维推广至三维。

针对问题1，先运用surfer软件绘制出8种主要重金属元素在该城区的空间分布三维曲面图，再从污染指数角度，通过单因子指数法反映出不同区域内的重金属污染程度，并用内梅罗综合污染指数法，全面分析各污染物对土壤的不同作用，突出高浓度污染物

害等级：

区自身的特点，比较分析，得出重金属污染的主要原因为：工厂生产与交通污染。

针对问题3，先由土壤溶质运移的基本理论，分析重金属污染物的传播特性，由此确立传播系数，运用多元非线性拟合，建立模型，结合surfer软件绘图，从而使用matlab 编程确定出污染源的位置及其坐标为：

元素污染源坐标

As （18367,10475）、（11965,2756.3）

Cd （21391,11613）、（2100.6,2714.9）、（6423.1,1913.3）

Cr （4353.1,65.015）、（4474.2,4358.4）

Cu （4248.1,179.43）、（1838.1,3317.3）

Hg （14901,9379.7）、（13501,2361.6）、（2766.5,2469.8）、（1777,2662.2）

Ni （6425.5,1914.3）

Pb （2433.7,1348.5）、（4506.3,5270.5）

Zn （13576,9598.3）、（9754.5,4686.8）

针对问题4，评价了问题3中模型的优缺点后，确定出还应收集的信息包括土壤密度、土壤含水率、土壤的横纵向弥散系数等参数，利用土壤溶质运移方程建模求解，最终得到三维的土壤溶质运移方程。

关键词：重金属surfer软件污染指数变异系数溶质运移

一、问题重述

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。

现要求通过数学建模来完成以下任务：

(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。

(2) 通过数据分析，说明重金属污染的主要原因。

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。

(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

二、问题分析

人类活动对城市环境质量的影响日益突出，使得对人类活动影响下城市地质环境演变模式的调查，显得日益重要。

题中给出了采样点的位置、海拔、空间分布信息，8种主要重金属元素在采样点出的浓度以及各自的背景值。要求给出各重金属元素在该城区的空间分布，结合数据分析不同区域的重金属污染程度以及重金属污染的主要原因，并建立模型分析污染源的传播特征，确定污染源的位置，最终对模型加以改进与推广，以更好地研究城市地质环境的演变模式。

首先，为了更加直观地表示出8种重金属元素在该城区各自的空间分布，可以运用surfer软件绘制各元素对应的三维曲面图。由污染指数角度看污染程度，通过单因子指数法与内梅罗（N C.Nemerow）综合污染指数法，对不同区域内的重金属污染程度进行分析。由生态危害看污染程度，运用潜在生态危害指数法，得出各功能区内各重金属元素及综合的生态危害等级。

其次，为说明重金属污染的主要原因，可以通过变异系数法求出各重金属元素与不同功能区之间的变异系数，以说明人为活动的干扰作用及污染程度的严重性。并通过绘制各重金属元素在不同功能区内浓度的一组柱状图，联系生活实际以及一些常见的污染

源，结合计算结果分析得出主要原因。

再次，由土壤溶质运移的基本理论，分析重金属污染物的传播特性，由此确立传播系数，找出各个量之间的关系，运用多元非线性拟合，建立模型，结合surfer 软件绘图，从而确定出污染源的位置及其坐标。

最后，考虑到模型求解过程中，一部分工作是人为预测和手工计算所得，带有一定的主观因素，与实际情况不完全相符。另外，模型中对实际情况做出了较为理想的简化，使问题复杂度大大降低，但同时也造成与实际情况的偏差。在此，需要收集包括土壤密度、土壤含水率、土壤的横纵向弥散系数等各种参数，运用土壤溶质运移方程建模求解。

三、模型假设及符号说明

3.1模型假设

1.假设重金属污染传播与海拔无关。

2.假设传播系数λ是与土壤本身性质有关的参数。

3.假设重金属污染只在水平面上传播，不考虑纵向传播。

4.假设污染物在传播过程中只扩散，不衰减。 3.2符号说明

变量和缩略表示 意义

i P 土壤中污染物i 的环境质量指数

i C 各采样点污染物i 的实测值（1-?kg mg ） i S 各功能区污染物i 的背景值（1-?kg mg ） 综P 各功能区的综合污染指数

i f C 单个重金属的污染系数

i r E 单个重金属的潜在生态危害系数 RI 多种重金属的潜在生态危害系数 i r T 各重金属元素的毒性系数

CV 各重金属元素浓度的标准差系数 SD 各重金属元素浓度的标准差

ave C 各重金属元素浓度的均值 （1-?kg mg ） λ 重金属污染物在土壤中的传播系数

k C 污染源附近（污染范围内）任一点的土壤容质浓度（1-?kg mg ）

C ? 污染源附近相邻两点的土壤容质浓度差（1-?kg mg ）

r ? 污染源附近相邻两点的距离（m ） ),(y x 各取样点坐标 ()00,x y 重金属污染源坐标

t 时间 b ρ 土壤的密度 θ 土壤的含水率

d k 分配系数（溶质在土壤有机质和土壤液相中的浓度比） x v 横向对流速度 y v 横向对流速度 z v 纵向对流速度

D横向弥散系数

x

D横向弥散系数

y

D纵向弥散系数

z

四、模型建立及求解

4.1问题1

4.1.1各主要重金属元素在该城区的空间分布

为了给出直观的各重金属元素浓度空间分布图，根据附件1和附件2中的数据，运用surfer软件，根据采集点的x、y轴坐标和其重金属浓度绘制8种重金属元素的三维曲面图如下。

图1 As在该城区的浓度空间分布图

图2 Cd在该城区的浓度空间分布图

图3 Cr在该城区的浓度空间分布图

图4 Cu在该城区的浓度空间分布图

图5 Hg在该城区的浓度空间分布图

图6 Ni在该城区的浓度空间分布图

图7 Pb在该城区的浓度空间分布图

图8 Zn在该城区的浓度空间分布图

4.1.2该城区内不同区域重金属的污染程度

1）由污染指数看污染程度——单因子污染指数法与内梅罗综合指数法

土壤质量评价一般以单因子污染指数为主，通过单因子评价，可以确定出主要的污

染物及污染程度，它一般以污染指数来表示，借以消除量纲，便于个污染物之间的比较分析。

数据中给出了土壤中重金属i 的实测数据i C ，也给出了重金属的背景值（评价标准）

i S ，根据单因子污染指数法的计算公式：

i

i i S C

P = （1）

得出各功能区土壤中重金属i 的单项污染指数i P 。

由表1重金属元素，而同一重金属元素污染等级的最高值，集中出现在第二区（工业区）内。

单因子污染指数只能反别反映各个污染物的的污染程度，不能全面、综合地反映土壤的污染程度，故当土壤同时被多种重金属元素污染时，需将单因子污染指数按一定方法综合起来进行评测，即应用综合污染指数法评价。内梅罗（N C.Nemerow ）综合指数法就是综合污染指数法的一种。它同时兼顾了单因子污染指数的平均值和最高值，可以突出污染较重的污染物的作用，给较严重的污染物以较大的权值，也能较全面地反映土壤环境的总体质量，从而更客观地对土壤环境质量进行评价。

根据表1种求得的各功能区土壤中重金属i 的单项污染指数i P ，得到各功能区中所有重金属的单项污染指数最大值max P 和平均值ave P 。运用内梅罗综合污染指数法，代入公式：

2

22m a x

a v e P P P +=综 （2）

得出各功能区的综合污染指数，再采用国家土壤环境二级标准（GB15618-1995）对土壤污染水平分级。

合污染指数而言，第2区的污染程度比第1、4区还要严重很多。

2）由生态危害看污染程度——潜在生态危害指数法

瑞典著名地球化学家Hakanson 提出的潜在生态危害指数法是目前最为常用的重金属污染程度的方法之一。

单个重金属的污染系数i f C ：

i

i

i f S C C =

（3） 单个重金属的潜在生态危害系数i r E ：

i f i r i r C T E ?= （4） 多种重金属的潜在生态危害系数RI ：

i

i

i r i f i r i r S C T C T E RI ?

=?==∑∑∑ （5） 按照Hakanson 制定的标准化重金属毒性系数为评价依据，As 、Cd 、Cr 、Cu 、Hg 、Ni 、Pb 、Zn 的毒性系数i r T 分别取值为10、30、2、5、40、5、5、1。附件数据给出了土壤中重金属i 的实测数据i C ，也给出了重金属的背景值（评价标准）i S ，分别代入潜

在生态危害指数法的计算公式，求解出相应的i r E 、RI 。

表4 各功能区中单个重金属的潜在生态危害系数i

金属元素，而Cr 、Zn 的潜在生态危害系数普遍较小。这也一定程度上体现了这几种重金属元素的毒性大小：Hg 元素的毒性要远大于Cr 、Zn 。

RI

重金属的污染生态危害系数和生态危害指数分级标准列于表6。

表6 重金属的污染生态危害系数和生态危害指数分级标准

指数 类型

所处范围

污染程度

指数 类型

所处范围

污染程度

潜在 生态 危害 系数 i r E

40

潜在 生态

危害 指数

RI

150

300<≤RI 强生态危害

320160<≤i

r E 很强生态危害 320≥i r E

极强生态危害 600≥RI 很强生态危害

综合表4、表5、表6，得到各功能区的生态危害等级。

第1区第2区第3区第4区第5区

As 轻微轻微轻微轻微轻微

Cd 中等较强轻微较强中等

Cr 轻微轻微轻微轻微轻微

Cu 轻微中等轻微轻微轻微

Hg 较强极强中等极强较强

Ni 轻微轻微轻微轻微轻微

Pb 轻微轻微轻微轻微轻微

Zn 轻微轻微轻微轻微轻微

综合中等极强轻微极强中等由表7

区内重金属污染的综合生态危害等级的影响，也是最大的。

4.2问题2

4.2.1运用变异系数法判断重金属污染的主要原因

根据附件2中给出的数据，绘制出各重金属元素在不同功能区内浓度的一组柱状图，如图9，横坐标为功能区划分，纵坐标表示重金属元素的浓度。

As Cd

Cr Cu

Hg Ni

Pb Zn

图9 各重金属元素在不同功能区内的浓度大小柱状图

由图9显而易见，在不同功能区的土壤内，同一重金属元素的含量存在显著差异。结合问题1中求得的各区污染程度数据（表3、表7），总体重金属含量：工业区>主干道路区>生活区>公园绿地区>山区。

8种重金属元素中，除了Cr 元素以外，其他元素的最大值均出现在工业区。工业生产过程中产生 “三废”——废水、废气、废渣，其中的重金属广泛分布于大气、水体、土壤与生物体中，而底泥往往是重金属最后的归宿。冶金、机械、电镀等工业排放中富含Pb 、Ni 、Hg 、Cu 、Zn 等重金属元素，由此判断工业生产对土壤中的各种重金属污染影响都很大。

通过求变异系数来判断重金属污染的来源主要是自然因素还是人为因素，也可以在一定程度上体现污染程度的大小。变异系数越大，说明人类活动的干扰作用影响越大或污染程度越严重。

变异系数又称“标准差率”，是衡量资料中各观测值变异程度的另一个统计量。为了对比分析不同水平的变量数列之间标志值的变异程度，就必须消除水平高低的影响，这时就要计算变异系数。变异系数是以相对数形式表示的变异指标。它是通过变异指标中的全距、平均差或标准差与平均指标对比得到的。变异系数有全距系数、平均差系数和标准差系数等，常用的是标准差系数CV ，为标准差SD 与均值ave C 的比率，即：

a v e

C SD

CV （6）

根据附件2中给出的8种主要重金属元素的浓度，分别计算出各元素的均值与标准

差，进而得出其对应的变异系数。

由表8可以看出，就此8种重金属元素相比较，Cr、Cu、Hg、Zn这四种重金属元素对本地区土壤污染更为严重，受人类活动影响也更大。其中Hg污染最为突出，其次是Cu、Zn、Cr。

Hg污染主要来源于仪表厂、食盐电解、贵金属冶炼、化妆品、照明用灯、燃煤等，可见工业生产及交通污染对土壤中Hg浓度的影响最大。

Cu元素的主要污染来源是铜锌矿的开采和冶炼、金属加工、机械制造、钢铁生产等，可见工业生产是土壤中Cu元素污染的主要原因。

Zn污染的主要污染源有锌矿开采、冶炼加工、机械制造、燃煤、汽车轮胎磨损等，可见工业生产及交通污染对土壤中Zn浓度的影响最大。

总体看来，本城区土壤的重金属污染程度受工业生产及交通污染情况的影响最大，尤其是工业生产中的金属冶炼与矿产开采，交通污染中的汽车尾气排放与轮胎磨损。

综合以上分析及表1数据可得：

（1）Zn、和Cu元素污染对生活区的影响比较大。Zn的主要来源煤燃烧产生的粉尘、烟尘：Cu的主要来源农药等。

（2）Zn、Cu和Hg元素的污染对工业区的影响比较大。Zn主要来源于冶炼加工、机械制造以及镀锌、仪器仪表、有机会合成和造纸等工业的排放；Cu主要污染来源是铜锌矿的开采和冶炼、金属加工、机械制造、钢铁生产等。冶炼排放的烟尘是大气铜污染的主要来源；Hg主要来源于仪表厂、食盐电解、贵金属冶炼等。

（3）山区的污染是轻微的，其中各元素的浓度都接近背景值。

（4）Cu、Cd、Pb、Zn和Hg元素污染对主干道路区的影响比较大。；Cd可能来源于密集交通带来的尾气和粉尘；Zn来源于轮胎的磨损和油的泄漏；Cu主要来源于密集交通带来的尾气和粉尘以及轮胎的磨损和油的泄漏；Hg主要来源于尾气的排放，油的泄漏。

（5）Hg、Cu、Zn和Cd是公园绿地区的主要污染元素。它们主要来源于农药的喷洒。

4.3问题3

4.3.1重金属污染物的传播特征

重金属污染物属于持久性污染物。土壤溶质（包括污染物）在土壤中随着土壤水分的运移而运动，同时土壤溶质运移的过程中还会受到土壤基质、土壤溶质本身的物理化学性质以及土壤中其他化学成分的影响。土壤溶质运移的物理机理主要包括对流、扩散和机械弥散。

对流是指溶质随着土壤水分的运移过程而进行的运动过程。在单位时间内通过单位土壤横截面积的溶质质量称为溶质通量，溶质在土壤中的溶质通量与土壤水分通量和溶质的浓度有关。

扩散是由溶质分子的热运动引起的，具有从高浓度处向低浓度处运移，以求达到浓度相等的趋势，在自由水溶液中，溶质的扩散遵循Fick 第一定律。

机械弥散是由于土壤中的孔隙大小不一，形状各异，使得水溶液在这些孔隙中运移

时每个孔隙中的流速在大小和方向上都存在差异，溶质分子扩散并且加大运移的范围的一种运移现象。溶质的机械弥散作用是由土壤孔隙中水的微观流速的变化引起的，具体的有以下三个方面的因素：

○

1土壤孔隙的中心和边缘的流速不同，如图10中的a 图 所示； ○

2土壤孔隙的直径大小不一，因此通过的孔隙流速不同，如图10的b 图 所示； ○

3土壤孔隙的弯曲程度不同和封闭孔隙或者团粒内部孔隙水基本不流动而导致微观流速不同，如10图中的c 图 所示。

a b c

图10 引起土壤孔隙中水的微观流速的变化的三个因素

溶质的分子扩散作用和机械弥散作用都会引起溶质在土壤中的运移，但是在实际的运移过程中，这两种作用比较难以区分，而且这两种因素引起的溶质运移通量都有相同的数学表达式形式，所以在实际的测量计算过程中就把这两种因素一起考虑，并且把这两种因素的联合形式称为水动力弥散。

吸附是重金属污染物与土壤相互作用的一个重要过程。土壤的吸附能力越强，重金属污染物的运移难度就越大。同时，不同的重金属污染物，其物理化学性质有所差别，土壤对不同重金属的吸附能力也有所差别，因此，重金属污染物本身的物理化学性质对重金属在土壤中的运移也有影响。

另外，生物吸收以及溶质的化学衰变等因素也对重金属污染物的传播有一定的影响。

综合而言，重金属污染物的传播，整体在地域上近似符合以污染源为顶点的圆锥状分布，越靠近污染源，梯度越大。重金属污染物在其进入土壤的早期主要是受对流作用的影响，但是随着时间的推移，对流作用的影响不是那么明显，使得这些污染物进入土壤的后期主要是扩散和弥散在起作用。经过一定的时间后，重金属污染物在一定范围内的浓度趋于均匀。同时，土壤的吸附作用也对重金属污染物在土壤中的运移过程有着重要的影响，吸附越强，重金属污染物的运移难度就越大。 4.3.2模型的建立，污染源位置的确定

1）金属污染物在传播过程中，由污染源传播到下一个点的污染物的多少是与土壤本身性质有关的，如土壤密度、土壤含水率、对流速度、弥散系数、分配系数等。

为简化计算，不必一一考虑各参数值，而是假设一个传播系数λ，λ是与土壤本身性质有关的所有参数的一个线性组合值。污染源附近相邻两点之间的土壤溶质浓度差正比于此两点之间的距离r ?（假设足够小），则可以得到污染源附近（污染范围内）任一点的土壤容质浓度k C 的一个基本方程：

r C C k k ??-=-+λ1 （7） 即

r C ??-=?λ （8）

对（8）式两边积分：

??=-r

c

dr dc 0

λ （9）

可得：

r C C ?=-λ0 （10）

将r 用坐标x 、y 表示，并对（10）式两边取平方，即可得到污染源周围土壤溶质浓度与坐标之间的关系式：

])()[()(2020220y y x x C C -+-?=-λ （11） 其中，(0x ,0y )为污染源点坐标。

由（11）式可以看出，土壤溶质浓度在污染源周围大致呈正圆锥状数据分布，同时可以认为污染源处于空间分布图圆锥状图形的锥顶附近。

2）模型求解

第一步，运用surfer 软件绘得各重金属浓度空间分布三维曲面图和等高线图，可以由此初步确定各重金属污染源(0x ,0y )的所在区域。又因为污染源附近的重金属浓度往往很高，因而可以排除图中的部分圆锥状图形，不认为是金属污染源，从而进一步缩小范围。

第二步，在初步确定的重金属污染源附近，划定一个尽可能多地包含规则等高线的尽量小的范围，并在附件1中给出的所有采集点坐标数据中进行搜索，找出在这个范围内的所有采集点。

第三步，根据所得到的这些采集点，在MATLAB 中运用nlinfit 函数，以其坐标和重金属浓度为样本数据，对函数（11）进行拟合，可以求得（11）式中各项参数：0C 、λ、

0x 和0y ，进而得到污染源的坐标(0x ,0y )。

求出每种重金属的污染源坐标如表9。

表9 各重金属的污染源坐标

元素 污染源坐标 As （18367,10475）、（11965,2756.3） Cd （21391,11613）、（2100.6,2714.9）、（6423.1,1913.3） Cr （4353.1,65.015）、（4474.2,4358.4） Cu （4248.1,179.43）、（1838.1,3317.3） Hg （14901,9379.7）、（13501,2361.6）、（2766.5,2469.8）、（1777,2662.2） Ni （6425.5,1914.3） Pb （2433.7,1348.5）、（4506.3,5270.5） Zn （13576,9598.3）、（9754.5,4686.8）

由表9可见，Hg 的污染源数目是最多的，一定程度上也可体现出该城区Hg 污染的严重性。

4.4 问题4

4.4.1问题3中模型的优缺点 1）优点

○

1使用专业成图软件，给出更为直观的空间分布图； ○

2对实际情况做出了较为理想的简化，使问题复杂度大大降低； ○

3对同一问题运用多种方法分析，得到相近的结果，更贴合实际，由说服力。 2）缺点

○

1部分结果是人为预测和手工计算所得，带有一定的主观因素，使得模型的解与实际相比存在一定的偏差；

○

2模型中对实际情况做出了较为理想的简化，对结果的准确性有一定影响 4.4.2收集新信息，建立新模型

为了更好地研究城市地质环境的演变模式，可以以问题3中的模型为基础，再加以改进与推广。

问题3中的模型，将实际情况中繁杂的参数用一个简单的传播系数λ替代，λ是与土壤本身性质有关的所有参数的一个线性组合值。这样就使得问题的复杂度大大降低了，但同时，也降低了计算与分析结果的准确性。为了提高模型的准确性，使其更符合实际情况，需要充分考虑重金属污染物在土壤环境系统中传输运移的影响因素，并依据土壤水分和土壤溶质运移的基本理论以及数学梯度的概念，将对地质环境演变模式的研究，由二维提升到了三维，大大提高了模型与实际情况的相符度。

此时需要另外收集的信息有：土壤密度b ρ、土壤含水率θ、分配系数d k 、纵向对流

速度z v 、横向对流速度y x v v =、横向弥散系数y x D D =、纵向弥散系数z D 等。 （1）土壤水分运移的连续方程

根据Darcy 定律和Buckingham —Darcy 定律和质量守恒定律可以推导出土壤水分运移的基本方程。

假设在所研究的土壤水分空间中有任意一点),,(z y x ,并且以该点为中心选取一个无限小的六面体，该六面体的长、宽、高分别为z y x ???,,，根据质量守恒原理可知，在任意一小时间段t ?内，流入土壤六面体的水量等于流出土壤六面体的水量加上该时段内土壤六面体的水分含量的变化量加上该时段由于植物根系的吸收而损失的水量的和。用

式子可以表示为：

t ?时段内进入土体的水量=t ?时段内离开土体的水量+ t ?时段内土体内水分的变化量+由于植物根系的吸收而损失的水量。 若仅考虑一维垂直方向的水分流动，则有：

t ?时段内进入土体的水量= t y x t t z y x J w ????+)2

1

,,,( （12）

在上式中，w J 为通过横截面积为y x ??的z 方向上的水通量，因为时段t ?很小，故以t t ?+21

时刻的通量代表平均通量。

同理：

t ?时段内离开土体的水量= t y x t t z z y x J w ????+?+)2

1

,,,( （13）

t ?时段内土体内的水量变化量= z y x t z z y x z y x t t z z y x ????+-????+?+),2

1

,,(),21,,(θθ

（14）

用w r 表示根系的吸水速率，则有：

t ?时段内植物根系吸收的水量= t z y x r w ???? （15）

依据质量守恒原理，整理各式，当0,0→?→?t z 时，得到土壤水分运移的一维连续方程：

0=+??+

??w w r t

x J θ

（16） 在三维水分运移情况下，连续方程可写为：

0=+??+??+??+??w w w w r t

z J y J x J θ

（17） 由以上式子转化得到土壤水分运移的Richards 方程，再逐步转化为它的基质式形式。然后以含水量形式表示比水容重，引入土壤水分扩散率，进行积分变换得到土壤水分运移的连续方程的含水量方程，也称θ方程。

若仅考虑一维水平的水分运移，则θ方程为：

])([x

D x t ????=??θ

θθ （18） 若是一维垂直方向的水分运移，则θ方程为：

z

K z D x t ??-????=??)(])([θθθθ （19） （2）土壤溶质运移的基本理论（推广）

在问题3的分析中已经概述了土壤溶质运移的基本理论，以及对流、扩散、机械弥散的概念，在此不再赘述。

1）对流

土壤溶质的对流通量c J 与土壤水分通量w J 和溶质的浓度c 有关，它们之间的关系可

用数学表达式表示为：

c J =w J c （20）

将土壤的平均孔隙流速θ

w J

v =代入上式得：

c v J c θ= （21） 2）扩散

自由水溶液中溶质的扩散遵循Fick 第一定律，表示为：

z

c

D J d ??-=0 （22）

式中d J 为溶质在自由水中的分子扩散通量，0D 表示溶质分子在自由水中的扩散系数，z

c

??表示的是溶质的浓度梯度，负号表示溶质分子沿着浓度降低的方向扩散。 溶质分子在土壤中的扩散系数s D 与0D 的关系：

τ0D D s = （23）

式中τ为弯曲系数，是个无纲量的物理量，一般变化范围为0.3—0.7。

3）机械弥散

土壤溶质在运移过程中产生的机械弥散通量用数学表达式可写为：

z

c

D J h h ??-=θ （24）

在上式中，h J 表示的是溶质的机械弥散通量，h D 指的是溶质的机械弥散系数。h D 通常可以表示为孔隙的平均流速的函数，其数学表达式为：

n h v D ||λ= （25） 其中λ指的是弥散度，v 指的是孔隙的平均流速，而n 则是经验常数，一般n 近似等于1。

4）水动力弥散

溶质的分子扩散作用和机械弥散作用都会引起溶质在土壤中的运移，在实际的测量计算过程中就把这两种因素一起考虑，并且把这两种因素的联合形式称为水动力弥散，它所引起的溶质通量dh J 可以用数学表达式表达为：

z

c

v D z c D

J n dh ??+-=??-=)||(0λτθ （26）

在上式中，D 表示的是有效水动力弥散系数。

5）土壤溶质运移的基本方程

综合而言，溶质在土壤中的运移总通量是对流通量和水动力弥散通量的总和，用数学表达式表示为：

z

c

D c J J J J w dh c ??-=+= （27）

根据质量守恒原理，类似于推导土壤水分的连续方程，推导出一维土壤溶质运移的连续方程，把对流——弥散方程写成约旦求和的方式，再一步步转化得到土壤溶质流动是静态流时，一维情况下的方程：

z c

v z

c D t c ??-??=??22 （28）

6）有源汇项的土壤溶质运移方程

考虑到溶质受到土壤的吸附、生物吸收以及溶质的化学衰变，有源汇项的土壤溶质运移方程为：

z c v z

c D t c R

d ??-??=??22 （29） 在上式中 θ

ρd b k

Rd +=1 （30）

推广到三维的土壤溶质运移方程可以写成：

z c

v y c v x c v z

z D y y D x c D t c R z

y x z y x d ??-??-??-??+??+??=??222222 （31） 上式是一个多元高阶偏微分方程可以用差分法来解偏微分方程，首先需要对求解区域进行离散，建立差分网格和差分方程，然后再对建立的差分方程进行求解。所求得的结果是污染物的浓度关于时间和三维空间的改变的分布情况。

五、参考文献

[1]郑永林，何腾兵，等，贵阳市乌当区茶园土壤重金属环境质量状况及其评价，贵州农业科学， 39(4)：103-105，2011。

[2]张一修，王济，张浩，贵阳市区地表灰尘重金属污染分析与评价，生态环境学报，2011,20(1)：169-174,2011。

[3] 覃邦余， 重金属污染物在土壤环境系统中运移的建模与仿真[D]，广西师范大学硕士学位论文，2009。

[4]吴昀昭, 南京城郊农业土壤重金属污染的遥感地球化学基础研究[D], 南京大学博士研究生论文, 2005。

[5] 黄昌勇,土壤学[M],中国农业出版社,1999。 [6] 朱祖祥,土壤学[M],农业出版社, 1983。

[7] 陈同斌, 黄铭洪, 黄涣忠, 等,香港土壤中的重金属含量及其污染现状[J]，地理学报, 52(3): 228 – 236,1997。

[8] 李其林, 黄昀，重庆市近郊区蔬菜地土壤重金属含量变化及污染情况[J]，土壤通报, 33(2): 158 – 160，2002。

[9]周建利,北京市蔬菜及其相应土壤重金属污染评价[D],广西大学硕士生论文, 2005。