摘要：本文以时域反射仪(Time Domain Reflectometry，TDR)作为溶质锋的探测手段，对利用边界层方法确定的参数和利用穿透曲线拟合法确定的参数进行了比较，结果表明边界层方法确定的扩散弥散系数(D)与穿透曲线拟合法确定的扩散弥散系数基本相近。但，因为TDR灵敏度的限制，边界层方法确定的延迟因子(R)大于穿透曲线拟合法确定的延迟因子。应用边界层方法确定的参数，比较了用边界层方法和精确方法预测的浓度剖面图，表明边界层方法在一定时间内可以精确地预测污染物的动态浓度分布。实验结果不仅说明了TDR作为溶质锋探测手段的可行性，同时说明了边界层方法在确定溶质迁移参数方面有一定的可靠性。
关键词：多孔介质 溶质迁移 边界层理论 溶质锋
CDE模型是描述化学物质迁移规律的数学模型，该模型的应用需要两个基本的参数：水动力弥散系数D和延迟因子R，因此准确确定化学物质在特定多孔介质中迁移的模型参数D和R是模型得以应用的前提。关于模型参数的估计，已提出了许多方法[1～4]，而这些方法存在着收敛性和参数唯一性的问题。边界层方法是利用溶质锋迁移位置与时间的关系进行CDE方程参数的估计，可以省去做穿透曲线的麻烦。但此方法还仅仅停留在理论层次上，尚未从实践的角度加以验证，况且如何准确确定溶质锋的位置还是一个有待解决的问题。本文通过对边界层方法和传统的穿透曲线拟合方法进行比较，初步确定利用TDR探测确定溶质锋位置具有可行性。
1 理论基础
邵明安等(1998)通过对CDE方程的研究，在均质稳态条件下，结合边界层理论，提出了确定CDE方程中参数的边界层理论模型[5]。该理论通过对农业污染物在多孔介质中迁移过程中浓度锋面的描述和概化，得到了下述表达式：
|
|
(1) |
式中：d(t)表示t时间溶质锋的位置。D是水动力弥散系数，R是延迟因子，V是平均孔隙流速。
土壤中的驻留浓度：
|
cr(x，t)=vd(t)c0/vd(t)+3D(1-x/d(t))3 |
(2) |
对于半无限土柱或田间土壤剖面，相应于通量型入流边界条件，土壤溶质迁移过程可以利用边界层理论进行描述。
由方程(1)可知，只要求得时间t时的溶质锋的深度d(t)，我们就可确定溶质迁移参数D、R。为了便于对数据t和d(t)进行处理，下面对方程(1)进行进一步简化。两边同时平方，得
|
d2(t)=4vt/Rd(t)+12Dt/R |
(2) |
方程两边同除以t2，得
|
1/t=R/12Dd2(t)/t2-v/3Dd(t)/t |
(3) |
2 材料与方法
表1 实验设置及土柱物理参数
|
| ||||||
|
土样 |
高度/cm |
ρs/(g/cm3) |
ρv/(g/cm3) |
Vp/(cm/h) |
探测点个数 |
水头/cm |
|
| ||||||
|
黄绵土 |
46 |
2.74 |
1.30 |
2.14 |
4 |
8.0 |
|
沙壤土 |
68 |
2.70 |
1.34 |
1.05 |
6 |
5.5 |
|
娄 土 |
48 |
2.65 |
1.27 |
1.25 |
4 |
9.0 |
|
| ||||||
供试土样分别为榆林沙壤土耕层、安塞黄绵土耕层及杨凌娄土耕层，过2mm筛，实验所用的有机玻璃管高为100cm，管内径为14cm.玻璃管壁上每隔10cm钻内径为2cm的小孔，以便插TDR(时域反射仪)探头，共有6个。TDR为波兰EASY TEST公司生产的FOM/mts型，能够同时测定土壤中的水分、盐度和土壤温度。探头选择10cm长的双针式探头。控制一定的容重，采用分层填装法填装土柱，采用马氏瓶控制水头。根据边界层方法，实验在稳态饱和条件下进行。由于边界层方法是在均质稳态不考虑源汇项的条件下建立的，所以本实验中选择化学性质不活泼的Cl-为迁移离子。微机自动控制TDR，连续测定探测点电解质的浓度变化。下端用100ml的容量瓶接取出流液，获得该土柱Cl-的穿透曲线。AgNO3滴定法测定出流液中Cl-含量。
3 结果与讨论
3.1 溶质锋位置的判断 边界层方法的应用首先需要确定溶质锋的位置，而溶质锋是一个迁移物质在迁移过程中从无到有的界面，而所采用的设备手段总是有一定的灵敏度，所以准确探测溶质锋有一定的难度。常用的监测溶质迁移的技术设备一般有溶液采样器和电阻感应探头，而TDR则是另一种新的用于溶质迁移研究的设备手段[6～8]。本实验所采用的TDR的具体性能参数见表2。实验中采用长度为10cm的双针式探头，它的感测范围是一个内径为5cm，长度为13cm圆柱体，与实验所采用的土柱相适宜。
表2 FOM/mts型TDR性能参数
|
| |||
|
|
测量范围 |
测量精度 |
分辨率 |
|
| |||
|
土壤含水量 |
0%～100% |
±2 |
%0.1% |
|
总盐/(s/m) |
0.000～1 |
±10% |
0.0001 |
|
土壤温度(℃) |
-20℃～+50℃ |
±0.8% |
0.1℃ |
|
| |||
表3 边界层测量数据
|
| ||||||
|
|
沙壤土 |
黄绵土 |
娄土 | |||
|
探测点 |
| |||||
|
|
d(t)/cm |
t/h |
d(t)/cm |
t/h |
d(t)/cm |
t/h |
|
| ||||||
|
1st |
11 |
4.28 |
11 |
2.63 |
11 |
3.45 |
|
2nd |
21 |
9.62 |
21 |
6.32 |
21 |
8.78 |
|
3rd |
31 |
13.48 |
31 |
9.95 |
31 |
13.58 |
|
4th |
41 |
19.17 |
41 |
13.65 |
41 |
18.67 |
|
5th |
51 |
26.22 |
|
|
|
|
|
6th |
61 |
36.18 |
|
|
|
|
|
| ||||||
由于实验处理为饱和稳态条件，在溶质锋未到达探测点前，TDR所探测到电解质的浓度为一恒定值。当该值发生变化时，也就是电解质浓度开始增大时，就断定在这个时刻，溶质迁移到了该点，这样就得到了溶质锋迁移到第一个探测点时的时间。随后将TDR移向下一个探测点。各探测点位置及溶质锋到达的时间见表3.
图1为沙壤土土柱中TDR在前5个点监测电解质浓度变化的部分过程图，图2为TDR在第6个点的完整过程图。从图1可以看出，每个监测点的电解质浓度变化曲线从开始的一段时间内是一条直线，浓度没有变化，此时我们认为溶质锋还没有到达该点。随时间的推移，浓度开始逐步增加，我们认为浓度开始增加的那个点对应的时间是溶质锋到达该点的时间。从图1还可以看出，随着探测点的下移，电解质浓度变化的曲线逐渐变缓，也就是迁移物质的扩散距离逐渐增大，浓度梯度逐渐变小，这与实际情况相符。从图2可以看出，TDR记录的完整的电解质浓度变化过程图的形状跟穿透曲线的形状极其相似，这也说明了用TDR探测溶质锋的可行性。
|
图1 沙壤土前5个探测点浓度变化过程 |
图2 沙壤土第6探测点浓度变化过程 |
|
3.2 穿透曲线方法求D、R参数 穿透曲线形状不仅可以反映溶质迁移的机理和溶质与土壤之间的作用，而且可以把穿透曲线的数据借助CXTFIT程序用最小二乘法进行拟合，得到溶质迁移的两个重要参数扩散弥散系数D和延迟因子R.因为Cl-的化学性质不活泼且带负电荷，故选用CXTFIT程序中的第二个数学模型即线形平衡吸附模型进行拟合。其控制方程为： |
|
|
|
(4) |
式中：cf为通量型溶质浓度。求得的3种土壤的Cl-迁移参数，见表4.
表4 两种方法求得D、R参数对照
|
| ||||||
|
|
沙壤土 |
黄绵土 |
娄土 | |||
|
参数 |
| |||||
|
|
BTC |
BLTM |
BTC |
BLTM |
BTC |
BLTM |
|
| ||||||
|
D |
2.09 |
2.26 |
4.44 |
5.13 |
1.47 |
3.39 |
|
R |
0.93 |
2.44 |
1.16 |
3.34 |
0.75 |
2.88 |
|
| ||||||
3.3 边界层方法求D、R参数 对表3中的数据按式(3)的形式再进一步处理，进行二次多项式拟合，分别得到Cl-溶质锋在3种土壤中运动的数学表达式(5)～(7)及图3.
沙壤土：
|
y=0.0902x2-0.1548x |
(5) |
|
r2=0.8091 |
黄绵土：
|
y=0.0552x2-0.1391x |
(6) |
|
r2=0.9971 |
娄土：
|
y=0.0633x2-0.1103x |
(7) |
|
r2=0.9896 |
拟合的相关系数较高，拟合的结果具有可信性。根据公式各项系数可求得溶质Cl-迁移参数D、R，见表4.
3.4 两种方法的比较与分析 从表4可以看出，两种方法得到的水动力弥散系数D值除娄土的相差比较大外，其他两种土样的D值比较相近，而用边界层方法得到的3种土样的延迟因子R值都比用穿透曲线方法获得的R值大。其原因可能是由于TDR所断定的溶质锋的位置比实际溶质锋的位置靠后所致。因为TDR有一个灵敏度，只有当电解质的浓度变化量达到仪器的灵敏度时，浓度的变化才被指示出来。而理论上边界层是迁移物质从无到有的一个层面，层面上的浓度变化很小时，TDR探测不到，所以等TDR的探测值发生变化时，溶质锋已经过了探测点的位置。溶质锋位置确定的延迟在表观上就等同于迁移离子的延迟因子R增大，使得溶质锋向前推移的速度变慢。从式(2)中也可以看出，在其他参数不变的前提下，随d(t)增大，R也随之增大，这正与我们的上述分析相符。
为检验边界层方法获得参数的可应用性，我们可以用边界层方法拟合的驻留浓度剖面与用CDE方程的精确解拟合的浓度剖面进行比较。根据方程(2)
|
cr(x，t)=vd(t)c0/vd(t)+3D(1-x/d(t))3 |
(8) |
在d(t)、D已知的情况下可求得从溶质锋到土柱上边界的浓度剖面图。而CDE方程的解析解的表达式为[9]：
|
cr(x，t)/c0=1/2erfc[Rx-vt/2(DRt)0.5]+(v2t/πDR)1/2exp[-(Rx-vt)2/4DRt]-f(x，t) |
(9) |
|
f(x，t)=1/2(1+vx/D+v2t/DR)exp(vx/D)erfc[Rx+vt/2(DRt)0.5] |
(10) |
图4～图93种土样以边界层方法获得D、R参数代入方程(8)和方程(9)计算的浓度剖面对比图。
|
图4 沙壤土4.28h浓度剖面 |
图5 沙壤土26.2h浓度剖面 |
|
图6 黄绵土2.62h浓度剖面 |
图7 黄绵土13.65h浓度剖面 |
|
图8 娄土3.45h浓度剖面 |
图9 娄土18.67h浓度剖面 |
由图4～图9可以看出，在一定时间内，边界层计算的结果同精确方法计算的结果基本一致。但随时间的增加，边界层方法计算的结果与精确方法计算的结果之间的误差会逐渐增大，说明用边界层方法在给定迁移参数情况下预测迁移物质浓度分布时，给定的时间不能过长，即不能预测日期太久后的浓度分布。
4 结 论
研究结果表明，TDR可以用来监测溶质锋的迁移，与拟合穿透曲线法相比，边界层方法确定的参数具有一定可靠性。TDR因受其灵敏度的制约，所判断的溶质锋比实际的位置靠后，导致边界层方法确定的延迟因子R大于拟合穿透曲线方法得到的R，但两种方法得到的D值相差不大。从两种方法预测的在不同时段内的迁移物质动态分布来看，在一定的时间内，用边界层方法预测迁移物质的浓度分布剖面具有较高的精度。
参 考 文 献：
[1] Rifai M N E，Kaufman W J，Todd D K.Dispersion phenomena in laminar flow through porous media[J]。Inst.of Eng.Res.1956，93(2).Sanitary Eng.Lab.Univ.of Calif.，Berkeley.
[2] Elprince A M，Day P R.Fitting solute breakthrough equations to data using two adjustable parameters[J]。Soil Sci.Soc.Am.J.，1977，41:39～41.
[3] Kool J B，Park J C，Van Genuchten M T H.Parameter estimation for unsaturated flow and transport model-a review[J]。Hydrol.1987，91:255～293.
[4] Buthter B，Hiniz C，Flury M，Fluhler H.Heterogeneous Flow and solute transport in an unsaturated stony soil monolith[J]。Soil Sci.Soc.Am.J.，1995，59:14～21.
[5] Shao Mingan et al.An approximate solution to the convection-dispersion equation of solute transport in soil[J]。Soil Science，1998，163:339～345.
[6] Bear J.Some experiments in Dispersion[J]。Geophys Res.1961，66:2455～2467.
[7] Mallants D，M Vanclooster， N Toride，J Vanderborght，J Feyen.Comparison of three methods to calibrate TDR for monitoring solute movement in undisturbed soil[J]。Soil Sci.Soc.Am.J.1996，60:747～754.
[8] Ridiclife D E，S M Gupte，J E Box Jr.Solutetransport at the pedon and polypedon scale[J]。Nutr.Cycl.Agroecosyst，1998，50:77～84.
[9] Van Genuchten M T H，Park J C.Boundary conditions for displacement experiment through short laboratory soil columns[J]。Soil Sci.Soc.Am.J.，1984，48:703～708.