目前,就非饱和土壤水动力弥散系数的测定来看,还没有公认而成熟的方法和规范可寻,当然,国内外一些学者在这方面也做了不少探索和研究。Yule和Gardner(1978)在假设弥散系数与速度成比例关系以及含水量均匀的前提下,进行室内短柱试验求得非饱和纵向和横向弥散系数,但由于假设偏于理想化,求得的参数难以体现实际情况。Smiles和Philip(1978),Smiles等(1978)求得水平吸水过程中溶质运移问题的半解析解,通过一维水平吸水实验,认为弥散系数仅为含水量的函数,与流速无关。De Smedt和Wierenga(1979,1984)在长30cm的一维垂直土柱中对两种不同粒径的玻璃球进行实验,认为弥散系数与平均孔隙流速呈线性关系。Jones和Watson(1982)用沙进行一维吸水实验,通过计算结果分析,当取弥散系数与平均流速呈线性关系时,计算结果完全落在实验结果的范围之内。杨金忠(1986)利用水平土柱试验,由数值方法反求参数,是个较有效的方法,但是求解非线性水流方程和对流弥散方程的复杂性,使之难以推广应用。黄康乐(1987,1988)基于质量守恒原理,借鉴求解水力传导度的瞬时剖面方法,提出了一种在实验和计算上都较为简单的室内和野外试验方法,并通过室内、野外试验证明该方法是较有效、精确的。石元春、李韵珠和陆锦文等(1986)以及清华大学的谢森传、杨诗秀和雷志栋(1989)进行了水平土柱的入渗试验,并根据试验结果求得了以含水率为变量的水动力弥散系数。张瑜芳、张蔚榛和沈荣开等(1997)提出,若已知土壤水、盐运动过程中某两个时刻的剖面分布,从质量守恒原理建立起剖面上各点的水分及盐分均衡方程,从而求出剖面上各点的弥散系数,此结果与根据实验用数值方法反求参数的结果相一致。
图2.3.3 扩散度拟合曲线
目前,对水动力弥散系数的结构形式的认识尚不统一(王亚东、胡毓骐,1992)。从理论上讲,水动力弥散系数Dsh为分子扩散系数Ds和机械弥散系数Dh之和。一般将溶质在土壤中的分子扩散系数仅表示为含水率的函数,而与溶质的浓度无关,常用经验公式来表示(雷志栋,杨诗秀,谢森传,1988)。用经验公式表示的分子扩散系数Ds为:
Ds=D0αebθ (2.3.55)
式中:Ds分子扩散系数(cm2/min);D0溶质在自由水体中的扩散系数(cm2/min);θ土壤含水率(cm3/cm3);α、b均为经验常数。
据文献介绍(Olsen 和Kemper,1968),当土壤水吸力在0.3~15atm 的范围内变化时,上述经验公式中b=10 比较适合,α的变化范围为0.005~0.001(沙壤土-粘土),土壤粘性愈大,α值愈小。
一般认为,一维流情况下,机械弥散系数 Dh与平均孔隙流速 υ 的一次方成正比(Bear,1972)
Dh=α|υ| (2.3.56)
式中:Dh机械弥散系数(cm2/min);υ平均孔隙流速(cm/min);α弥散度(为经验常数)(cm)。
综上所述,弥散系数Dsh表示为分子扩散系数Ds和机械弥散系数Dh之和,即
Dsh=D0αebθ+α|υ| (2.3.57)
当对流速度相当大时,机械弥散的作用会大大超过分子扩散作用,以致于水动力弥散中只需考虑机械弥散作用;反之,当土壤溶液静止时,则机械弥散完全不起作用,而只剩下分子扩散了。一般情况下,土壤中的溶质运移,都同时存在分子扩散和机械弥散作用,但实际上很难区分开来,因此,将分子扩散和机械弥散综合统称为水动力弥散。实际应用中,有的学者将水动力弥散系数表示为形如分子扩散系数形式的指数函数,如 Smiles 和 Philip(1978),谢森传、杨诗秀和雷志栋(1989),认为纵向弥散系数对孔隙水流速不敏感,因此,Dsh可以单独作为含水率的函数来对待。但从文献资料看,目前不少学者将水动力弥散系数表示为形如机械弥散系数形式的线性函数,认为Dsh与平均孔隙流速υ的一次方成正比。本文所测定的水动力弥散系数取前一种形式。测定方法有水平土柱法和垂直土柱法。
图2.3.4 垂直土柱试验装置示意图
(一)垂直土柱法
试验装置如图2.3.4 所示,土柱上装有负压计和盐分传感器以测定土壤负压和土壤溶液浓度,供试溶液由马氏瓶从底部进入土柱。为了计算水动力弥散系数(Hydrodynamic Dispersion Coefficiet),首先计算水分通量,然后计算盐分通量,最后由水分通量和盐分通量计算水动力弥散系数。
1.水分通量
若已知溶液从底部补给土柱的水量,以及不同时刻剖面含水率的分布,则由水量均衡原理,土柱上任一截面z处的水分通量qz可表示为:
土壤水盐运移数值模拟
即
土壤水盐运移数值模拟
式中:qz为任一截面z处的水分通量(cm/d),q0为土柱底部的进水量(cm/d),θ为体积含水率(cm3/cm3),Δt=t2-t1为时段(d)。
上式(2.3.59)写为离散格式:
土壤水盐运移数值模拟
式中:k为时段数。
2.盐分通量
若已知土柱底部溶质通量,以及不同时刻剖面含水率和溶质浓度的分布,则任一截面z处的溶质通量Jz由质量守恒原理得:
土壤水盐运移数值模拟
土壤水盐运移数值模拟
式中:Jz为任一截面 z 处的溶质通量(g/cm2·d);J0为土柱底部的溶质通量(g/cm2·d);c为土壤溶质浓度(g/cm3);θ为体积含水率(cm3/cm3);Δt=t2-t1为时段(d)。
式(2.3.62)写为离散格式为:
土壤水盐运移数值模拟
3.水动力弥散系数
根据水动力弥散原理,溶质通量等于水动力弥散通量与对流通量之和,即:
土壤水盐运移数值模拟
土壤水盐运移数值模拟
式中:J 为溶质通量(g/cm2·d);Dsh为水动力弥散系数(cm2/d);c 为溶质浓度(g/cm3);θ为体积含水率(cm3/cm3);Δt=t2-t1为时段(d)。
式(2.3.65)写为离散格式为:
土壤水盐运移数值模拟
将前面计算出的
(二)水平土柱吸渗法
试验装置如图2.3.5所示,溶液由马氏瓶从土柱一端水平渗入,土柱为初始含水率和盐分含量均匀一致的半无限土柱,这个问题可以用如下的水盐运移方程进行描述。
图2.3.5 水平土柱试验装置示意图
水分方程:
基本方程,
土壤水盐运移数值模拟
式中:D(θ)为水分扩散度(cm2/min);θ为与输入端(进水边界)的水平距离为x处的体积含水率(cm3/cm3)。
定解条件,
土壤水盐运移数值模拟
式中:θi为初始体积含水率(cm3/cm3);θs饱和体积含水率(开始试验后在边界处瞬时形成)。
盐分运移方程:
基本方程,
土壤水盐运移数值模拟
式中:Dsh为水动力弥散系数(cm2/min);c为与输入端(进水边界)的水平距离为x处的溶质浓度(g/cm3);q为水流通量(cm/min);θ为体积含水率(cm3/cm3)。
定解条件,
土壤水盐运移数值模拟
式中:ci为初始土壤溶液浓度(g/cm3);c0为所供给溶液浓度(g/cm3)。
由水分方程可以解出扩散度:
土壤水盐运移数值模拟
由盐分方程可以解出水动力弥散系数,由于,
土壤水盐运移数值模拟
所以盐分运移的基本方程式(2.3.69)可以展成:
土壤水盐运移数值模拟
采用 Boltzmann 变换,将上述偏微分方程化为常微分方程,令
土壤水盐运移数值模拟
将
土壤水盐运移数值模拟
令
土壤水盐运移数值模拟
将式(2.3.71)代入式(2.3.76)得:
土壤水盐运移数值模拟
将上式写为离散格式为:
土壤水盐运移数值模拟
式(2.3.75)可写为:
土壤水盐运移数值模拟
Boltzmann变换后盐分运移问题的定解条件变为:
土壤水盐运移数值模拟
将上式(2.3.79)两边在区间[c,ci]上积分,求出水动力弥散系数:
土壤水盐运移数值模拟
写为离散格式为:
土壤水盐运移数值模拟
根据试验数据用式(2.3.82)即可计算水动力弥散系数Dsh。
(三)水动力弥散系数测定结果
本书采用水平土柱吸渗法进行水动力弥散试验。由于不同溶质在土壤中的弥散系数基本相同(张瑜芳、张蔚榛和沈荣开等,1997);通过不同浓度的入渗试验证明,入渗溶液浓度和初始含水量对Dsh影响不明显(石元春、李韵珠和陆锦文等,1986);理论分析和实验证明,入渗溶液的浓度对土壤水分的运动影响很小(谢森传、杨诗秀和雷志栋,1989)。因此,本书选用氯化钠溶液作为供水水源进行弥散试验。
试验装置为分节的有机玻璃圆柱(图 2.3.5),柱长 70cm,内径 2.5cm,每节长3.5cm,节与节之间为钟罩式连接,柱的一端装有多孔板,供水装置为马氏瓶。测试土样同前,为寅阳1#粉砂壤土,大兴2#粉砂壤土,兴隆沙1#粉质粘壤土,土壤含盐量及离子组成见表2.3.5。其中寅阳1#砂壤土,兴隆沙1#粉质粘壤土的土壤盐分均以氯化钠为主,Cl-和Na+的含量占绝对优势,而大兴2#砂壤土离子含量则以
表2.3.5 土样含盐量及离子组成
试验结束后,迅速将土柱按节拆开取样。土壤含水率采用烘干法测定,土壤含盐量采用电导率仪测定。通过实验数据拟合的电导率与土壤含盐量的换算关系为
s=2.8882Ec+ 0.1016 (2.3.83)
式中:s为土壤含盐量(单位质量干土所含盐分的质量(g/kg));Ec为电导率(土水比为1:5的浸提液,标准为103档下的读数(mS/cm))。
土壤溶液浓度c与土壤含盐量s的换算关系为:
θc=γs (2.3.84)
式中:c为土壤溶液浓度(g/L);θ土壤含水率(cm3/cm3);γ 为干土容重(g/cm3);s土壤含盐量(g/kg)。
根据试验的实测数据,按照上述算法进行计算。拟合的水动力弥散系数的经验公式如下:
寅阳1#(相关系数R=0.987)
Dsh(θ)=8×10-6e30.187θ (2.3.85)
大兴2#(相关系数R=0.981)
Dsh(θ)=4×10-8e47.965θ (2.3.86)
兴隆沙1#(相关系数R=0.993)
Dsh(θ)=0.0061e12.448θ (2.3.87)
主要计算图件及拟合曲线见图2.3.6至图2.3.8。
图2.3.6 寅阳1#曲线图
图2.3.7 大兴2#曲线图
图2.3.8 兴隆沙1#曲线图