非饱和土壤紊散动力学系数测定方法及装置与流程

本发明属于土壤质量检测技术领域，具体的说是一种非饱和土壤紊散动力学系数测定方法及装置。

背景技术：

随着农业生产的发展和科学技术水平的提高，要求农田合理施肥，经济灌溉和有效预测、控制土壤盐碱化，防治农药，化肥及废水灌溉引起的土壤和地下水污染，以及土壤溶质随地表径流迁移，引起水环境污染和土壤生产力降低，已经成为当今世界普遍关注的问题。对此研究涉及到的对土壤溶质的研究。对土壤溶质研究最广泛的是利用数学模型，模拟土壤溶质运动，其中，水动力弥散系数是不可或缺的参数之一。

溶解于土壤水的物质(化肥、农药、盐分、污染物等)在土壤空隙中运移，称之为土壤溶质运移。促使可溶性物质在土壤中运移的机制有：质流(也称对流)，水动力弥散。水动力弥散系数是表征一定流速下，多孔介质对某种污染物质弥散能力的参数。在宏观上反映了多孔介质中地下水流动过程和空隙结构特征对溶质运移过程的影响。包括机械弥散系数与分子扩散系数。是土壤水盐运移的一个重要参数。

非饱和土壤紊散动力学系数是化肥、农药在农田的移动规律、盐碱地水盐运动检测、地下水资源保护中不可缺少的参数。目前，虽然有关饱和条件下土壤弥散系数溶质运移各类参数的测定已较为成熟，但非饱和条件下由于土壤中实际孔隙流速随含水率的变化表现出数量级的变化，以及非饱和条件下溶质迁移非线性驱动机制，由于各类问题本身的复杂特性，使得对水动力紊散系数参数的测定相当困难，需要发展测试理论和测试技术。

国内外测定非饱和土壤弥散动力学系数的方法有：公式法、水平土柱吸渗法、瞬时剖面法。这些方法的不同在于测试的边界条件和控制过程的差异，以及与边界条件和控制过程对应的解析方法的差异，然而就测试理念来说，并没有差异，都是测定能够表述溶质迁移所表现出的宏观特性的等效参数，而不是正真意义的参数。而对等效参数，仅适合在与测试条件下相似的情况下，使用参数描述土壤中的溶质迁移，而当实际土壤条件与测试条件有较大差异的情况下，参数的适用性则显著降低。而本项研究提出的非饱和土壤水动力紊散参数测定方法则能够有效的解决这一问题。

技术实现要素：

为了克服上述现有技术的不足，本发明提供了一种非饱和土壤紊散动力学系数测定的方法，具有操作简便、费用低、测定周期短、误差小的特点。

本发明的方法所采用的技术方案是：

一种非饱和土壤紊散动力学系数测定的方法，包括以下步骤：

步骤1：进行土样装填，在装填过程中，通过土壤初始含水率控制，土壤容重控制和土壤非饱和渗透性控制消除土壤分层和各断面不均匀现象；

步骤2：根据土样的质地，以及水流在土壤中的流动速度，按照一定间隔(5-10cm)布设监测位点，各监测位点均布设土壤基质势监测系统和土壤示踪溶质浓度监测系统，用于对各监测位点的土壤含水率变化以及土壤中示踪溶质浓度变化进行连续监测；

步骤3：基于各监测位点进行断面划分，各断面位于相邻两个监测位点中间位置，将各断面相邻的上、下监测位点所在断面分别作为各断面的上、下边界，通过各断面上、下边界的土壤基质势监测系统测定的土壤基质势确定水势差(负压坡降)，然后采用达西定律确定各断面的水分运动通量；

步骤4：在任一时刻t，根据各断面水分运动通量，计算该时刻平均断面水分通量；

步骤5：基于步骤2各监测位点测定的土壤负压，根据土壤基质势与土壤含水率之间土壤水分特征曲线关系，将土壤负压转化为土壤含水率，进一步耦合步骤4测定的断面平均水分通量，采用如下公式求得各测试位点的溶质迁移的紊散动力学系数：

其中，d为溶质迁移的紊散动力系数，j为溶质的质量通量，θ为土壤体积含水率，c为土壤示踪溶质浓度，t为时间，z为垂直坐标，q为水分运动通量，相邻监测位置节点深度差值δz＝zi+1-zi，相邻监测时间节点差值δt＝tj+1-tj，i表示监测位置节点，i＝1,2,3…,j表示监测时间节点，j＝1，2，3…，zi和zi+1分别为i和i+1监测位置节点的深度，1/2表示监测位置或时间节点的中点位置，tj和tj+1分别表示j和j+1时刻。

上述步骤1中，土壤初始含水率控制：土壤在105℃条件下完全烘干后，用孔径为2mm的震动筛去除影响土壤结构的大孔径土粒后，根据土壤重量，确定添加水量，采用表面喷洒拌合的方式实现土壤初始含水率，要求黏土达到10～15％的饱和含水率，要求壤土达到20～30％的饱和含水率。

上述步骤1中，土壤容重控制：根据所设定的土壤容重，确定分层装填的土壤重量，采用气动压缩机根据设定容重所确定的断面位置进行土样装填。

上述步骤1中，为控制装填土壤的质量，在10⁶pa的压强下进行压气实验，装填土壤控制在12小时，压缩位移和非饱和渗透性指标需要达到以下控制要求：黏土：位移小于装填土壤长度的1％，非饱和渗透性小于10^-5m/s；壤土：位移小于装填土壤长度的5％，非饱和渗透性小于10^-4m/s。

本发明还提供基于上述方法的一种非饱和土壤紊散动力学系数测定装置，所述装置包括边界控制系统和监测系统，边界控制系统用于控制实验边界条件以及实现土壤的示踪溶质的控制注入，监测系统包括土壤基质势监测系统和土壤示踪溶质浓度监测系统；所述边界控制系统能够实现多种条件下的多参数测试控制，包括：①定水头、定浓度示踪溶质连续注入，②定水头、脉冲示踪溶质注入，③定流量、定浓度示踪溶质连续注入，以及④定流量、脉冲示踪溶质注入；所述监测系统对示踪溶质迁移过程中的土壤含水率变化以及土壤中示踪溶质浓度变化进行全过程连续监测。

作为优选，所述土壤基质势监测系统选用负压计，所述土壤示踪溶质浓度监测系统选用盐桥传感器。

定流量边界，根据达西定律通过动态控制外压水头边界实现，土壤中的流动符合达西定律，通量与压力正比，根据土壤负压计测定的土壤基质，控制外压水头，实现定流量过程入渗。

定水头直接通过供平水装置(马氏瓶)控制外压水头即可。

根据测定的土壤含水率，以及示踪溶质(盐分)的分布，能够基于式(13)采用有限差分法确定土壤紊散动力学系数。通过对不同的工况条件下的溶质迁移动力学过程进行测试，各监测位置的土壤含水率以及示踪溶质不断发生变化，而这种变化，是由于对流和弥散作用所造成的，而紊散动力学系数的测定需要通过对土柱的所有检测点的连续过程监测，确定弥散作用的贡献，进而确定参数。

作为优选，本方法为了更加准确的得到剖面含水率，计算每一时段土柱含水量的增值δw：

式中，θ为计算时段第i点的含水率，xi为第i点与第i+1点的间距，t和t+1分别表示前后两个时刻。

土柱含水量的增值δw应等于边界处在这一时段的供水量。如两者相差悬殊，则必须检查是否计算有误，在相对误差大于5％的情况下则应修正含水率剖面，直至计算值和测定值相对误差在5％以内。

本发明方法包括测试方法，以及数据解析方法两个相互联系的部分，根据测定的土壤含水率，以及示踪溶质(盐分)的分布，基于式(13)采用有限差分法确定土壤紊散动力学系数，与现有技术相比，本发明具有以下有益效果：

(1)边界条件和控制过程对应的解析方法的差异，然而就测试理念来说，并没有差异，都是测定能够表述溶质迁移所表现出的宏观特性的等效参数，而不是正真意义的参数。等效参数所表现出的宏观特性在很大程度上依赖于实验条件，因此，当求解问题与参数的测定条件有显著差异的情况，等效参数的有效性显著降低。而本发明的方法则在程度上克服了这一问题。

(2)非饱和流动条件下，溶质迁移通量中对流通量和弥散通量的解析是非常困难的，本发明方法能够同步实现对流通量的解析，所提出的紊散水动力参数则在机理上更为准确。

附图说明

图1为非饱和弥散参数测试试验装置示意图；

图2为定流量入渗条件下土壤水分和盐分分布；

图3为定流量入渗条件下各点弥散系数和对应点的平均孔隙流速关系曲线；

图4为定水头入渗条件下土壤水分和盐分分布；

图5为定水头入渗条件下各点弥散系数和对应点的平均孔隙流速关系曲线。

具体实施方式

下面结合附图及实施例对本发明进行详细说明。以下实施例将有助于本领域技术人员进一步理解和实施本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

步骤1：进行土样装填，实验装置如附图1所示，实验装填时，通过土壤初始含水率控制、容重控制和非饱和渗透性控制，消除由于土壤结构具有较强的空间变异性和各项异性而造成的分层和各断面不均匀等现象：

1.1，土壤壤初始含水率控制：土壤在105℃条件下完全烘干后，用孔径为2mm的震动筛去除影响土壤结构的大孔径土粒后，根据土壤重量，以及控制含水率(黏土：10～15％的饱和含水率，壤土20～30％的饱和含水率)，确定添加水量，采用表面喷洒拌合的方式实现土壤初始含水率。

1.2，容重控制根据设定土壤容重，确定分层装填的土壤重量，采用气动压缩机根据设定容重所确定的断面位置进行土样装填。

1.3，为控制装填土壤的质量，在10⁶pa的压强下进行压气实验，装填土壤在12hr时间，压缩位移和非饱和渗透性指标需要达到控制要求(黏土：位移小于装填土壤长度的1％，非饱和渗透性小于10^-5m/s，壤土位移小于装填土壤长度的5％，非饱和渗透性小于10^-4m/s)

实验装置(附图1)包括(1)边界控制系统，控制实验边界条件，以及实现土壤的示踪溶质的控制注入，(2)监测系统，包括土壤基质势监测系统(如负压计)和土壤中示踪溶质的监测系统(如盐桥传感器)。非饱和条件下，土壤中流速随土壤含水率变化而显著变化，为了使得测定的水动力紊散系数能够覆盖整个非饱和条件下流速变化范围，至少需要对4种边界控制条件下的2种条件进行测试。边界控制系统能够实现：①定水头、定浓度示踪溶质连续注入，②定水头、脉冲示踪溶质注入，③定流量、定浓度示踪溶质连续注入，以及④定流量、脉冲示踪溶质注入等多种条件下的多参数测试控制。监测系统对示踪溶质迁移过程中，土壤中含水率以及示踪溶质浓度变化进行全过程连续监测。

定流量边界，通过动态控制外压水头边界实现，土壤中的流动符合达西定律，通量与压力正比，根据土壤负压计测定的土壤基质，控制外压水头，实现定流量过程入渗。

定水头直接通过供平水装置(马氏瓶)控制外压水头即可。

根据测定的土壤含水率，以及示踪溶质(盐分)的分布，能够基于式(13)采用有限差分法确定土壤紊散动力学系数。

通过对不同的工况条件下的溶质迁移动力学过程进行测试，各监测点的土壤含水率以及示踪溶质不断发生变化，而这种变化，是由于对流和弥散作用所造成的，而紊散动力学系数的测定需要通过对土柱的所有检测点的连续过程监测，确定弥散作用的贡献，进而确定参数。

填装土样的容器为直径为12cm，长100cm的有机玻璃，土壤按照预定的初始含水率和浓度拌匀，按照容重进行装填，在试验土柱中按一定距离埋设负压计和盐分传感器，观测水盐运动过程。

步骤2：根据测试土样的质地以及水流在土壤中的流动速度，确定负压计和盐分传感器的安装间隔；对于黏土，按照5cm间隔要装负压计和盐分传感器，对于壤土，按照10cm间隔要装负压计和盐分传感器。在测试阶段，盐分传感器对各监测点的土壤示踪溶质浓度进行连续监测，负压计对各监测点的土壤水势进行连续监测，土壤基质势根据相应溶质条件下的土壤水势、示踪溶质浓度得到。

步骤3：基于传感器布设进行断面划分，断面位于相邻两个传感器中间位置，根据断面上、下边界位置所布设的负压计的所测定的土壤基质势确定水势差(负压坡降)后，采用达西定律确定断面的水分运动通量。

步骤4：在任一时刻t，根据各断面水分运动通量，计算该时刻断面平均水分运动通量。

步骤5：基于步骤2测定的各监测点的土壤水势，根据土壤基质势与土壤体积含水率之间土壤水分特征曲线关系，将土壤水势转化为土壤体积含水率，如附图2所示，确定所有传感器布设节点位置的土壤含水率，进一步考核土壤测定电导率关系确定土壤溶液浓度。在不同的控制边界条件下，土壤含水率以及浓度的变化范围应覆盖土壤非标和区域(附图2和附图4)。

根据所确定的土壤含水率和孔隙率，进一步耦合断面平均水分运动通量，确定断面平均流速v：

其中，a为断面面积，n为土壤孔隙率。其他符号意义同上。

采用式(13)，采用有限差分法，求得各监测点紊散动力学系数；

其中，d为溶质迁移的紊散动力系数，j为溶质的质量通量，θ为土壤体积含水率，c为土壤示踪溶质浓度，t为时间，z为垂直坐标，q为水分运动通量，相邻监测位置节点深度差值δz＝zi+1-zi，相邻监测时间节点差值δt＝tj+1-tj，i表示监测位置节点，i＝1,2,3…,j表示监测时间节点，j＝1，2，3…，zi和zi+1分别为i和i+1监测位置节点的深度，1/2表示监测位置或时间节点的中点位置，tj和tj+1分别表示j和j+1时刻。在此基础上，根据各个节点的土壤物理状态以及d的测定(附图3和附图5)关系，采用解析法或瞬时剖面法确立紊散动力学函数关系。

本实验采用了解析法和瞬时剖面法，分别对应实验条件为连续溶质注入和脉冲溶质注入的条件下确定紊散动力学系数。

(一)解析法

解析法是根据影响非饱和紊散动力系数的主要因素，提出的确定紊散动力学系数的直接方法。水动力弥散系数是综合反映溶质、土壤特性的参数,既和多孔介质的状况及溶质的性质有关,又受含水量和孔隙水流速度的影响。水动力弥散系数的一般表达式为:

dsh(θ,v)＝ds(θ)+dh(v)式(1)

式中：ds(θ)为分子扩散系数；dh(v)为机械弥散系数。

情况1：土壤介质中，当流速对应的雷诺数超过3.0后，当对流速度相当大时，机械弥散作用会大大超过分子扩散作用，以致水动力弥散中只考虑机械弥散作用，此时紊散动力学系数的表达式为：

式中αl、αt为非饱和土壤的纵向、横向弥散度，；vi、vj，i、j两个方向的孔隙水流速度；δij为delta函数，当两个下标相同时，取1，相异时取值为零。说明：两个方向为纵向和横向(纵向指流动方向，横向为垂直于流动方向，两个下标表示方向，相同即表示方向相同，相异即表示方向不同)

在进行数值计算，为减少数值弥散的影响，常用d′sh作为弥散系数：

d′sh＝dsh-vδx/2式(3)

式中，v＝q/θ，(见式2)，δx为渗流长度。

当土壤溶液静止时，则机械弥散完全不起作用，只剩下分子扩散作用，此时，dsh可以用以下表达较好：

dsh＝d0ae^bθ

式中，d0为溶质在自由溶液中的扩散系数，取值与溶质分子类型有关；θ为土壤体积的含水量；a，b为经验系数，a＝0.001～0.005，b＝10。

(二)瞬时剖面法

用瞬时剖面法求解紊散动力学系数，根据水分及溶质质量守恒原理,便可建立起剖面上各点的水分、溶质均衡方程,从而求出剖面上各点的水动力弥散系数。

如果已知土壤水、盐运动过程中任意两个时刻的剖面水盐分布，则根据水分以及盐分质量守恒原理，即可建立起剖面上各点的水分及盐分均衡方程，从而求出剖面各点的弥散系数。

根据质量守恒原理，一维非饱和溶质迁移方程为：

式中，j为溶质的质量通量，表示为：

其中，θ为土壤体积含水率，c为溶液浓度，t为时间，d为溶质迁移的紊散动力系数，z为垂直坐标，q为水分运动通量。由水流连续方程得：

式(4)两边积分，并在测点i位置采用有限差分法进行展开，得：

其中，i和j分别表示测点的位置和时间步长节点，zi和zi+1分别为i和i+1位置节点的深度，1/2表示位置或时间节点的重点位置，tj和tj+1分别表示j和j+1时刻。

式(7)左边可表示为：

表示zi处tj时刻的溶质的质量通量。式(7)右边采用一阶积分近似，得：

式(9)为式(7)的进一步推导，为中间步骤。

对(9)时间导数采用中心差分近似，得：

δt＝tj+1-tj，δz＝zi+1-zi

式(10)为式(9)的进一步推导，为中间步骤。

综合式(7)～式(10)得：

由式(5)得：

由式(11)和式(12)得：

式(13)即为求解紊散动力系数的计算公式，其中为：

式(14)为式(13)中的的表述，对应实验中所确定的各个节点的水分通量。

根据步骤1-5得到的土壤示踪溶质浓度土壤体积含水率i＝1,2,3,…m，j＝1,2,3,…n，计算：

土壤体积含水率

土壤示踪溶质浓度

土壤体积含水率

土壤示踪溶质浓度

土壤含水率和示踪剂浓度额中点值可直接采用几何平均方法计算。

因此，在已知上边界条件(zi＝0，tj，i＝1，j＝1,2,3,…n)的溶质的质量通量i＝1，j＝1,2,3,…n的条件下，

将上述数据代入式(11)中，

求得溶质的质量通量i＝1，j＝1,2,3,…n，再将相关数据代入式(11)，可求得溶质的质量通量i＝1，j＝1,2,3,…n，根据该方法可得到i＝1,2,3,…m，j＝1,2,3,…n。

根据步骤1-5得到的土壤体积含水率i＝1,2,3,…m，j＝1,2,3,…n，

因此，在已知上边界条件(zi＝0，tj，i＝1，j＝1,2,3,…n)的水分运动通量i＝1，j＝1,2,3,…n的条件下，

将上述数据代入式(14)中，

求得水分运动通量i＝1，j＝1,2,3,…n，再将数据代入式(14)，可求得水分运动通量i＝1，j＝1,2,3,…n，根据该方法可得到i＝1,2,3,…m，j＝1,2,3,…n。

当保持上边界(zi＝0，i＝1)为定流量、定浓度入渗时，

j(0,t)＝q(0,t)，c(0,t)＝q0c0式(15)

式(15)中的j(0,t)、q(0,t)、c(0,t)、q0、c0，

j(0,t)为上边界位置(zi＝0，i＝1)、测试时间(tj，j＝1,2,3,…n)的溶质的质量通量，

q(0,t)为上边界位置(zi＝0，i＝1)、测试时间(tj，j＝1,2,3,…n)的水分运动通量，

c(0,t)为上边界位置(zi＝0，i＝1)、测试时间(tj，j＝1,2,3,…n)的土壤示踪溶质浓度，

q0和c0分别为边界的初始通量和溶度。

最简单的边界条件为q(0,t)＝0，j(0,t)＝0，即土壤水盐的再分布期间的任意两个时刻的盐分剖面分布均可对紊散动力系数进行求解。土壤水盐的再分布，也就是实验开始后，土壤中各点的水分和溶质(盐分)会不断发生变化。

通量边界和定水头边界条件下，连续注入与脉冲注入的求解方法相同。

①定水头、定浓度示踪溶质连续注入，

②定水头、脉冲示踪溶质注入，

③定流量、定浓度示踪溶质连续注入，以及

④定流量、脉冲示踪溶质注入

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。