一种考虑空气阻力作用的土壤降雨入渗测定系统和方法与流程

本发明涉及一种考虑空气阻力作用的土壤降雨入渗测定系统和方法，是一种水文领域中，针对降雨在土壤中入渗过程的模拟系统和方法。

背景技术：

分析流域降雨-入渗-产流间的内在关系，实现流域产汇流过程的模拟，在水资源评价和管理、防洪减灾等领域应用广泛。降水进入土壤的过程中，水分不断进入土壤，土壤中的空气很难及时排出土体，空气被禁锢在土壤内，被禁锢在土壤中的空气阻碍了降雨在土壤中的入渗。禁锢空气的存在，导致湿润区的土壤孔隙一部分还残留着空气，土壤实际含水量和土壤实际导水系数均小于土壤饱和含水量和饱和导水系数；另外，禁锢空气产生的压力也影响了湿润锋处实际水头，湿润锋处实际水头为禁锢空气压力和土壤水吸力之差。当前，大部分降雨入渗产流模型均认为降雨过程中雨强恒定、土壤均质和无空气阻力，导致模型出现模拟失真、精度不高等问题。

技术实现要素：

为了克服现有技术的问题，本发明提出了一种考虑空气阻力作用的土壤降雨入渗测定系统和方法。所述系统可以更加真实地模拟流域降雨入渗产流过程，测定结果可为空气阻力对降雨入渗产流影响机理研究提供支撑。所述方法修正了green-ampt(格林-安姆)模型，方法中引入了饱和度系数(sθ、sk)量化了滞留空气对土壤含水量和土壤导水系数的影响，引入了进水值(hab)和进水值(hwb)量化了禁锢空气产生的空气压力对土壤水势的影响。同时，所述系统和方法充分考虑了土壤分层特征和降雨非稳定性，改善了模拟失真问题，提高了模拟精度。

本发明的目的是这样实现的：一种考虑空气阻力作用的土壤降雨入渗产流测定系统，包括马氏瓶供水装置，所述的马氏瓶供水装置与降雨器管路连接，所述的降雨器下部为有机玻璃土柱，所述有机玻璃土柱为竖直的柱形体，下端密闭，内填充土样，降雨器与土样之间为降雨空间，所述的有机玻璃土柱侧面设置土壤水分监测装置，所述的降雨器和有机玻璃土柱之间密封连接，使有机玻璃土柱和降雨器形成完整的密封体；用以收集密封体内所有空气和水分的集水集气管设置在降雨空间中的土样表面处，所述的集水集气管与水气收集装置管路连接，所述的有机玻璃土柱侧面按照不同的高度设置测定空气压力的u型管。

进一步的，所述马氏瓶供水装置包括：封闭的圆筒形有机玻璃的瓶体，所述的瓶体下侧设置供水口，瓶体顶部设置插入瓶内水中并带有调节阀的导气管，和用于注水的导水管，所述的导水管上设有开关阀，所述瓶体上设有校准瓶体水平的水准泡。

进一步的，所述的水气收集装置包括：与所述罐体降雨空间部分连通的水气集管，所述的水气集管插入密闭的集水瓶底部，所述的集水瓶顶部设置u型导气管，所述u型导气管与密闭的集气瓶下半部分连通，所述的集气瓶底部设有用于排水的u型导水管，所述u型导水管出口与所述u型导气管在集气瓶上的进口保持在同一水平面内。

一种使用上述系统的考虑空气阻力作用的土壤降雨入渗测定方法，所述方法的步骤如下：

土壤分层的步骤：用于使用不同质地的土壤分层填充所述的有机玻璃土柱，形成剖面分为多层的土壤土样，同一层内土壤水分特征参数相同，所述的土壤水分特征参数包括：土壤初始含水量、饱和导水系数、土壤饱和含水量、田间持水量、进气值和进水值；

降雨过程分段的步骤：用于在使用所述的降雨器进行模拟降雨，并将降雨过程分为连续的多个相等时段，每个时段内降雨强度恒定；

确定湿润锋处实际水头的步骤：用于考虑空气阻力作用并基于对每层土壤理化性质的分析，引入各层土壤的土壤进气值hab和土壤进水值hwb：

确定实际湿润锋处土壤水的水头swwf：

swwf＝swcf－swaf

式中：swcf为湿润锋处土壤水吸力，swcf＝hab；swaf为湿润锋处禁锢空气压力，swaf＝0.5(hab+hwb)；

可得：swwf＝hab-0.5(hab+hwb)；

进一步确定为湿润锋处总水头hl：

hl＝－l－(swcf－swaf)

hl＝-l-hab+0.5(hab+hwb)

式中：l为湿润锋到土壤表层距离；

引入饱和度系数(sθ,n、sk,n)，确定湿润锋以上各层土壤实际含水量θs,n和实际导水系数ka,n：

θa,n＝sθ,n·θs,n

ka,n＝sk,n·ks,n

式中：ks,n为饱和导水系数；θs,n为土壤饱和含水量；sθ,n、sk,n分别为土壤饱和含水量、土壤饱和导水系数所对应的饱和度系数；n为湿润锋以上土壤层数；

降雨入渗产流过程分析的步骤：根据时段初积水深度、本时段内降雨强度、本时段积水入渗率，分入渗过程为积水入渗过程和非积水入渗过程进行分析，分时段、分情景考虑降雨入渗产流过程，所述的情景包括情景a、情景b、情景c、情景d四种：

情景a：h0＝0，

其中，h0为时段初积水深度，i为雨强，为湿润区实际平均导水系数，fpt为积水时段的土壤入渗率，fnpt为非积水时段的土壤入渗率；

此时，地表开始积水入渗，tp为地表积水开始时刻，入渗过程分为积水时段和非积水时段两个阶段，两个时段以tp为分界时间点；

情景b：h0>0，

其中：p’为对应深度的雨强，

由于积水深度h0与雨强间单位不统一，计算时，先将h0除以时段长度△t，转换为雨强p’,即：p’＝h0/△t，△t＝tx-tx-1；tx、tx-1分别为第x时段和第x-1时段；

此时，地表积水逐渐消失；

情景c：h0>0，

此时，下渗过程一直处于积水入渗过程；

情景d：h0＝0，此时，入渗过程一直处于不积水状态；

降雨入渗产流过程计算的步骤：包括积水入渗过程计算、产流过程计算、积水与非积水转换时刻计算：

所述的积水入渗过程计算包括：

湿润锋在第1层时土壤入渗过程计算：

a＝[swcf,1-0.5(hab,1+hwb,1)]δθ1

式中：为土壤累积入渗量；为积水发生时刻土壤累计入渗量；ka，1为第1层的实际导水系数；swcf，1为第1层中湿润锋处的土壤水吸力；hab，1为第1层的土壤进气值；hwb，1为第1层土壤进水值；a为参数；dq为土壤含水变化量，加脚标表示某一层的土壤含水变化量。

其中：dq1＝qa，1－q0，1，qa，1为湿润区第1层土壤实际含水量，q0，1为第1层初始土壤含水量；

湿润锋在第m(2≤m≤n)层时，积水入渗过程计算：

式中，ka，m为第m层土壤的实际导水系数；bm－1、cm－1为参数；fm－1为m-1层以上的土壤累积入渗量；tm－1为湿润锋进入m-1和m层界面时的时间；swcf，i为第i层中湿润锋处的土壤水吸力；hab，i为第i层土壤的进气值；hwb，i为第i层土壤的进水值；li为第i层土壤的厚度；

其中：dqi＝qa，i－q0，i，qa，i为湿润区第i层土壤实际含水量，q0，i为第i层土壤初始含水量；i和m为土层数；

时段产流过程计算：

式中：为时段tx的累积地表产流；为时段tx的累积降雨；h0为时段tx的地表积水深；为时段tx的累积入渗量；

非积水时段入渗过程计算：

积水与非积水转换时刻tp计算：

当湿润锋在第1层时：

当湿润锋在第m层时：

式中，tx-－1为时段x-1的区间时间；为时段tx－1时的累积入渗量；tp为土壤表层积水发生或消失时间；为土壤表层积水发生或消失时刻所在时段的降雨强度，fp为非积水时段与积水时段转换时刻的土壤入渗率。

进一步的，所述的swaf湿润锋处空气压力的最大值为：hb＝h0+l0+hab；

最小值为：hc＝h0+l0+hwb；

式中：hb为突破压力值；hc为闭合压力值；l0为饱和区厚度。

由于饱和区的厚度很小，远小于hab和hwb，可以忽略不计，因此，在初始阶段和过渡阶段hb≈h0+hab，hc≈h0+hwb；稳定阶段时hb≈hc≈h0+0.5(hab+hwb)，假设初始阶段和过渡阶段，空气压力从突破压力到闭合压力的变化呈线性，则平均空气压力可以表示为：

swaf＝h0+0.5(hab+hwb)

进一步的，所述的土壤进气值hab和土壤进水值hwb，土壤进气值hab＝1/αd，土壤进水值hwb＝hab/2－δ

式中：αd为brooks-corey模型中的参数，δ为经验值。

进一步的，所述的经验值δ的适用范围：砂土层的δ＝0～2cm，壤土中的δ＝2～5cm，黏土中的δ＝8～10cm。

进一步的，针对土壤实际导水系数计算，第一层土壤i＝1中，所述的饱和度系数sk约等于土壤持水量qh与饱和土壤含水量qs之比：sk≈qh/qs；第i层土壤i>1中，所述饱和度系数sk约等于0.5。

进一步的，针对实际土壤含水量计算，将湿润锋以上区分为饱和区和非饱和区，饱和区sθ等于1，即饱和区实际含水量等于饱和含水量qa＝qs；在非饱和区，饱和度等于土壤持水量qh与饱和土壤含水量qs之比：sθ≈qh/qs，即：非饱和区实际含水量等于田间持水量qa＝qh。

本发明产生的有益效果是：基于green-ampt(格林-安姆)模型，通过考虑降雨不稳定、土壤分层和土壤空气禁锢等影响因素，进一步修正了green-ampt模型，使其更加符合空气阻力影响下的分层土壤非稳定降雨入渗产流过程，模型改进提高了入渗产流过程模拟精度。所述方法提供的降雨入渗产流模拟方法，充分考虑了土壤分层、非稳定降雨和空气阻力对降雨入渗产流模拟的影响，改善了模拟失真问题，提高了模拟精度。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的实施例一、二、三、所述系统的结构原理示意图；

图2是本发明的实施例四所述方法的计算流程图；

图3是本发明的实施例四所述方法的分层土壤非稳定降雨入渗产流示意图；

图4是本发明的实施例四所述的入渗产流过程中的情景a情景示意图；

图5是本发明的实施例四所述的入渗产流过程中的情景b情景示意图；

图6是本发明的实施例四所述的入渗产流过程中的情景c情景示意图；

图7是本发明的实施例四所述的入渗产流过程中的情景d情景示意图。

具体实施方式

实施例一：

本实施例是一种考虑空气阻力作用的土壤降雨入渗测定系统，如图1所示。所述系统包括：马氏瓶供水装置1，所述的马氏瓶供水装置与降雨器3管路连接，所述的降雨器下部为有机玻璃土柱2，所述有机玻璃土柱为竖直的柱形体，下端密闭，内填充土样202，降雨器与土样之间为降雨空间201，所述的有机玻璃土柱侧面设置土壤水分监测装置4，所述的降雨器和有机玻璃土柱之间密封连接，使有机玻璃土柱和降雨器形成完整的密封体；用以收集密封体内所有空气和水分的集水集气管设置在降雨空间中的土样表面处，所述的集水集气管与水气收集装置6管路连接，所述的有机玻璃土柱侧面按照不同的高度设置测定空气压力的u型管501，按照土层数量设置的多个u型管组成空气压力测定装置5。

本实施例所述的马氏瓶供水装置用于持续、稳定地向降雨器中供水。马氏瓶是一种可提供稳定供水装置。这种装置简单，使用方便，具有一定的控制精度。

有机玻璃土柱是放置土壤样品的容器，为模拟降雨时土壤中水分由浅入深的变化，有机玻璃土柱一般要有一定的高度，使土壤样品有一定的深度，因此，有机玻璃土柱可以是柱形的罐状，如使用竖直的圆筒(一段管子)，上下两端封闭。为观察方便，圆筒使用有机玻璃管子，作为有机玻璃土柱的材料。有机玻璃土柱中装填土壤样品，模拟入渗产流过程，装填方式可以根据实验需要确定。土壤样品可以从被研究区域的土壤中取样，并在有机玻璃土柱中装填、压实，形成与被研究区域相同的土层样本。取样时最好使用与有机玻璃土柱相同直径的土壤取样器，将土层完整的从原始土壤中取出，保持原土样的构成，然后整体放入有机玻璃土柱中进行实验，以最接近的方式还原被实验土壤的构成。

为实现分层土壤实验，在装填土样的时候，可以根据需要取不同的土样：如砂石土、壤土等，分层装填，形成不同的土层，以便进行实验。有机玻璃土柱的直径和高度试验时根据要求制定，例如可选择直径40cm、高度100cm。

所述的降雨器采用喷头阵列，设置在土柱罐的顶部。喷头阵列是有多个喷头在一个平面内纵横排列而形成。喷头可以采用注射器的针头，并根据实验雨量的需要变换使用不同型号的针头。整个喷头阵列尽可能的保持水平，以便产生十分均匀的模拟降雨。为此，可以在喷头阵列上设置水准泡。

用来模拟自然降雨的降雨器，可以使用喷头阵列。喷头阵列有多种形式可以选择，如采用排列的水管，在水管侧面成排的打孔形成成排的喷口，并使用大小不同直径的喷口形成各种不同的雨强。

降雨器包括：喷头阵列，所述的喷头阵列中的各个喷头301通过管路互相连接，并与马氏瓶供水装置的供水口管路连接，所述的喷头阵列上设有校准喷到阵列水平的水准泡302，如图1所示。

也可以使用纵横排列的注射器针头，形成注射器针头阵列。多个注射器针头均匀的排布在一个水平面内，各个针头通过管路互相连通，并与马氏供水装置连通，使水流可以顺畅的流入各个喷头。针头的排列可以采用纵横、交错和放射性排列等方式。为控制雨量的强度，可以采用各种大小不同型号的注射器针头，如：采用5号针头控制的降雨器模拟0.10～0.32mm/min强度的降雨，7号注射器的针头模拟0.32～1.90mm/min强度的降雨，12号注射器的针头模拟大于1.90mm/min强度的降雨。

在同一直径的喷口或同一型号的针头下，通过调整马氏瓶发泡点位置改变降雨器内的水位也可以达到调节降雨强度的目的。

与马氏瓶供水装置类似，降雨器上部边缘以同样的方式呈三角形安放个水准泡，用于调平降雨器，使得降雨器保持水平。

降雨器和土柱罐间采用密封接口，不透气。传统的入渗产流实验，由于没有考虑空气阻力作用，一般忽略入渗时土壤中的空气释放出来，因此，采用敞开的方式，入渗发生时的空气完全自然的释放到空气中了。本实施例将降雨器与土样之间完全密封，将入渗发生时，土壤中释放的空气完全收集起来，进行计量，充分考虑空气的作用。

有机玻璃土柱侧面分别布设土壤水分监测装置和空气压力测定装置。土壤水分监测装置由竖排的土壤水分传感器和数据采集器构成，通过这些设施确定土壤中的水分特征参数。

所述土壤水分监测装置可以这样设置：按照土柱的高度由上到下均匀分布的土壤水分传感器阵列401，所述的土壤水分传感器阵列中的各个土壤水分传感器分别与数据采集器402电连接，如图1所示。

传感器阵列采用的是在土柱罐侧壁，从上到下设置一列传感器，这些传感器分别检测从上到下的不同土层。传感器的个数根据实验时对土层分层的要求，每个土层可以设置1到2个传感器。各个传感器通过导线与数据采集器连接，而数据采集器可以通过直接与电脑连接，或采用存储装置存储所采集的数据，最终输入电脑进行处理。

空气压力测定装置可以使用多个u型管，这些u型管竖直排列在土柱罐的一侧。u型气压检测管的个数也与土层的分层有关，如：每个土层设置一个u型气压检测管。

u型气压检测管设有水柱和刻度，通过水柱的位置的变化，通过刻度直接读取土壤中气压的变化。也可以在u型气压检测管上设置传感器，将气压变化直接输入数据采集器中。

有机玻璃土柱在土层表面布设集气集水口，连接塑料软管作为集水集气管，用于连接积水瓶和有机玻璃土柱，将土壤表面积留的水，以及入渗过程中土壤释放出来的空气收集起来，进行测定。为此，集水集气口可以设置在土柱与降雨空间衔接的位置，即在土柱的顶端，暴露在降雨空间中，以便土样顶部积留的水可以顺畅的流入集水集气口。

水气收集装置可以采用集水瓶和集气瓶串联的方式，在集水瓶和集气瓶上设置刻度，通过读取刻度而检测水、空气的排出量。

实施例二：

本实施例是实施例一的改进，是实施例一关于马氏瓶供水装置的细化。本实施例所述马氏瓶供水装置包括：封闭的圆筒形有机玻璃的瓶体，所述的瓶体下侧设置供水口105，瓶体顶部设置插入瓶内水中并带有调节阀102的导气管101，和用于注水的导水管103，所述的导水管上设有开关阀104，所述瓶体上还设有校准瓶体水平位置的水准泡106，如图1所示。

为使供水装置尽可能水平，可以在马氏瓶供水装置上部边缘呈三角形安放个水准泡，使水准泡均匀分布于马氏瓶供水箱的三个方向，以更好地实现马氏瓶的调平，提高在读取马氏瓶水位变化时的准确性。为使供水量数字化，可以在瓶中设置水位传感器，或者简单在瓶壁上设置刻度，直接读取水位数据。

实施例三：

本实施例是上述实施例的改进，是上述实施例关于水气收集装置的细化。本实施例所述的水气收集装置包括：与所述罐体降雨空间部分连通的水气集管601，所述的水气集管插入密闭的集水瓶602底部，所述的集水瓶顶部设置u型导气管603，所述u型导气管与密闭的集气瓶604下半部分连通，所述的集气瓶底部设有用于排水的u型导水管605，所述u型导水管出口与所述u型导气管在集气瓶上的进口保持在同一水平面内，如图1所示。

集水瓶是土体排出空气的中转站和水分收集器。试验过程中，集水瓶不仅收集土柱内土壤表面生成的水量，也收集土体内排出的空气，并通过u型管将空气排入集气瓶；集气瓶用来收集从土柱中排出的气体，这些气体来自于土壤空隙。降雨入渗过程中，土体内的空气一部分气体滞留在土壤中，一部分排出土壤，进入集气瓶。试验开始前，集气瓶中装满水分，当空气进入集气瓶中，水分排出，通过测量水分的排出体积就可得知进入空气的体积。

实施例四：

本实施例是一种考虑空气阻力作用的土壤非稳定降雨入渗测定方法，所述方法的流程如图2所示。本实施例所述方法的基本原理是：选定一块研究区域土壤，在整个降雨过程中对土壤的降雨入渗进行测定。本实施例可以使用上述实施例所述的测定系统对被研究区域的土壤进行测定，也可以直接在被研究区域对降雨过程进行测定。图3是本实施例中的分层土壤非稳定降雨入渗产流示意图。与传统的降雨入渗测定方法相比，本实施例提出了：

1.空气阻力计算方法：

由于空气压力和土壤水吸力的存在，湿润区实际的水头swwf，i可以由公式swwf，i＝swcf，i－swaf，i计算，其中swaf，i为湿润锋下方的空气压力，swcf，i为湿润锋处土壤水吸力，其中，脚标i表示的是在第i层的湿润锋处土壤水吸力。

入渗过程中，土壤中空气压力(swaf)是一个动态变化过程。peck推测，土壤中的空气排出土体需承受一定的压力，即突破压力，其压力可能等于静水压力与进气值之和。当土体内空气压力达到突破压力时，土体内产生排气通道，空气通过通道自由排出，土体内空气压力减小。当土体内空气压力减小到一定值后，即闭合压力，由于土体内空气压力的减小，排气通道重新闭合。wang等证实了peck的推测并给出了突破压力和闭合压力表达式，其计算公式分别为：

hb＝h0+l+hab

式中：hb为突破压力值；hab为土壤的进气值；l为湿润锋深度；h0为地表积水深度。

hc＝h0+l+hwb

式中：hc为闭合压力；hwb为进水值。可以看出，当湿润锋足够大时，wang等提出的模型认为土体内空气压力swaf可以趋向无穷，即：这显然与实际情况不符。

李援农通过实验研究了空气阻力影响下的入渗过程，结果表明空气压力的变化过程可以分成初始阶段、过渡阶段和稳定阶段三个阶段。在初始阶段，随着湿润锋的增加，禁锢空气压力值增加，空气从排出到禁锢的周期很短，空气压力在突破压力和闭合压力值之间动态重复变化；在过渡阶段，随着湿润锋的增加，空气从排出到禁锢的周期时间延长；在稳定阶段，空气排出的通道相对固定，土壤中的空气释放和空气禁锢过程一直在进行，此时的禁锢空气压力值相对固定，这个固定的空气压力值约为空气突破压力值和闭合压力值的平均值。湿润区内禁锢空气压力swaf不等于0，且湿润锋上下禁锢空气压力值相等。这说明并不是所有湿润区水分对空气均产生静水压力。

入渗过程中，湿润区可以被分成饱和区、过渡区、传导区和湿润区4部分。除禁锢空气压力值达到突破压力，产生排气通道，空气自由排出时段外，饱和区内土壤孔隙完全被水充满，空气不能自由移动，则饱和区内空气压力swaf＝0；对于过渡区、传导区和湿润区3个区域，孔隙不能完全被水填充，大孔隙间相互连通，土壤中的空气可以随意移动，此时空气压力swaf>0。

土壤中禁锢空气排出土体必须要达到足够高的空气压力来克服静水压力。在入渗过程中，由于过渡区、湿润区和传导区内部分孔隙尚没被水分充满，空气可以自由移动，土体内水分无法对禁锢空气产生静水压力，只有表面积水和饱和区的土壤水会对空气产生静水压力。

根据上述研究可以确定：swaf的最大值为hb＝h0+l0+hab，最小值为hc＝h0+l0+hwb，这样计算方法得到的结果更准确。其中，l0为饱和区厚度，h0为积水深度，hab为土壤进气值，hwb为土壤进水值。

由于饱和区的厚度很小，远小于hab和hwb，可以忽略不计。因此在初始阶段和过渡阶段hb≈h0+hab，hc≈h0+hwb；稳定阶段时hb≈hc≈h0+0.5(hab+hwb)。假设初始阶段和过渡阶段，空气压力从突破压力到闭合压力的变化呈线性，则平均空气压力可以表示为swaf＝h0+0.5(hab+hwb)。

其中，hab＝1/αd，hwb＝hab/2－δ，δ是经验常数，单位为cm，在不同的土壤中有不同值：在砂土层的δ＝0～2cm，壤土中的δ＝2～5cm，黏土中的δ＝8～10cm。αd为brooks-corey模型中的参数。

2.参数计算方法：

mgam模型关键参数包括饱和度系数(sk、sθ)、土壤实际水含量(qa)、湿润锋处土壤水吸力(swcf)、进气值(hab)和进水值(hwb)，合理确定这些参数是mgam模型成功应用的关键。

湿润区的土壤孔隙内包括残留空气和水分。由于滞留空气的存在，对土壤实际含水量来说，饱和度系数sθ反映了湿润区的饱和程度。因此，sθ可以表示为：

对土壤实际导水系数来说，由于滞留空气存在，导致土壤湿润区不能达到饱和，土壤湿润区处于非饱和状态，土壤实际导水系数与土壤含水量呈幂函数关系，则sk可以表示为：

式中：qra为残留的空气含量，cm³/cm³；qrw为残留含水量，cm³/cm³；α为与土壤性质有关的参数。

实际研究中，qra、qrw和α均很难确定。基于综合前人相关研究成果，本文近似给定了饱和度系数(sk、sθ)的计算方法。

在计算土壤实际含水量时，湿润区完全饱和的条件下，湿润区实际含水量等于土壤的总孔隙度，饱和度系数(sθ)等于1。降雨入渗过程中，bodman和colman认为可以将湿润区分为饱和区、传导区、过渡区和湿润区，其中饱和区中土壤达到饱和，传导区中土壤含水量接近饱和，sθ≈1。王文焰等认为湿润区可以分为饱和区和过渡区，其中饱和区占湿润区一半。为了便于计算，mgam模型也将湿润区分为饱和区和过渡区，饱和区sθ＝1，过渡区sθ＝qh/qs。

在计算土壤实际导水系数时，其饱和度系数(sk)与土壤分层有关(具体分层根据土质特性确定)，针对第1层土壤，sk＝qh/qs；针对第i(i≥2)层，饱和度系数sk约为0.5(sk≈0.5)。通过研究发现，这种计算方法对沙土、壤土和粘土均较适合。

湿润锋处空气压力swaf从气体闭合压力hc转化成为气体突破压力hb。whisler和bouwer(1970)提出，湿润锋处土壤水吸力swcf可以用brooks—corey模型中的进气值代替，swcf＝hab＝1/αd；湿润锋处空气压力swaf≈h0+0.5((hab+hwb)，其中，h0是积水深度；hab＝1/αd，αd是水分特征曲线参数；hwb＝hab/2-δ，沙土δ＝0～2cm；壤土为2～5cm；粘土为8～10cm。

本实施例所述方法的步骤具体如下：

(一)用于使用不同质地的土壤分层填充所述的有机玻璃土柱，形成剖面分为多层的土壤土样，同一层内土壤水分特征参数相同，所述的土壤水分特征参数包括：土壤初始含水量、饱和导水系数、土壤饱和水含量、田间持水量、进气值和进水值。

根据实验的需要在有机玻璃土柱中分层填充几种不同的土样，将土壤剖面分为多个层次，土壤的层次确定后，确定各层土壤水分特征参数。同一土层内应具有相同的土壤初始含水量、饱和导水系数、土壤饱和水含量、田间持水量、进气值和进水值等土壤水分特征参数。

为计算方便，定义土壤剖面分为n层，第i层土壤厚度为li(i＝1，2，3，…，n)，土壤初始含水量为q0，i，土壤饱和导水系数为ks，i，土壤饱和含水量为qs，i(ks，i和qs，i脚标中的s表示为“饱和”的意思)。

(二)降雨过程分段的步骤：使用所述的降雨器进行模拟降雨，并将降雨过程分为连续的多个相等时段，每个时段内降雨强度恒定。

降雨过程中，影响土壤入渗的主要因素包括降雨不稳定、土壤分层和土壤空气禁锢等，考虑到这些因素的影响，本实施例将降雨历时分为连续的x个时段，每个时段内降雨强度恒定，设为it(x)(其中脚标t(x)为时段x对时间t的函数)，以此将降雨强度变化的因素引入计算中。

(三)确定湿润锋处实际水头的步骤：用于考虑空气阻力作用并基于对每层土壤理化性质的分析，引入各层土壤的土壤进气值hab和土壤进水值hwb：

确定实际湿润锋处土壤水的水头swwf：

swwf＝swcf－swaf

其中：swcf为湿润锋处土壤水吸力，swcf＝hab；swaf为湿润锋处禁锢空气压力，swaf＝0.5(hab+hwb)。

swwf＝hab-0.5(hab+hwb)

进一步确定为湿润锋处总水头hl：

hl＝－l－(swcf－swaf)

hl＝-l-hab+0.5(hab+hwb)

其中：l为湿润锋到土壤表层距离。

为了考虑空气阻力作用，引入饱和度系数sk，i、sθ，i(其中脚标k和θ表示分别表示导水系数和含水量各自对应的饱和度系数)，确定空气阻力影响下的各层土壤的实际土壤含水量和实际导水系数分别为sθ，iqk，i、sa，iks，i。其中sθ，iqk，i、sa，iks，i大小通过分析土壤物理性质确定。引入进水值hwb和进气值hab，确定湿润锋处空气压力swaf＝0.5(hab+hwb)，然后根据湿润锋处土壤水吸力swcf，确定了湿润锋处实际水头。

(四)降雨入渗产流过程分析的步骤：用于根据时段初积水深度、本时段内降雨强度、本时段积水入渗率，分入渗过程为积水入渗过程和非积水入渗过程进行分析，分时段、分情景考虑降雨入渗产流过程。

根据时段初积水深度h0、本时段内降雨强度it(x)、本时段积水入渗率(即本时段潜在入渗率)fpt，确定tx-1～tx时段内的降雨入渗产流过程，根据h0、i和fpt不同，可以分为多个情景进行分析。

计算过程中，由于积水深h0与降雨强度单位不统一，无法直接代入模型计算。计算过程中，首先将积水深h0除以时段△t，转为对应深度的净雨强p’，然后代入模型计算，其中p’＝h0/△t，△t＝tx-tx-1。根据时段内降雨强度、时段初积水深和潜在入渗强度，时段内入渗过程可以分为以下情况。本实施例所述的情景可以分为情景a、情景b、情景c、情景d四种，当然这只是情景分析的一种分法，对情景还可以有更多种类的分法：

情景a(如图4所示)：h0＝0，

其中：h0为时段初积水深度，i为雨强，为湿润区实际平均导水系数，fpt为积水时段的土壤入渗率，fnpt为非积水时段的土壤入渗率，tp为积水开始时刻，入渗过程分为非积水时段(tp之前)和积水时段(tp之后)两个阶段，两个时段以tp为分界。这种情况下，随着降雨的持续，地表开始积水，土壤入渗过程可以分为非积水入渗过程和积水入渗过程。

情景b(如图5所示)：h0>0，

其中：p’为对应深度的雨强，

p’＝h0/△t，△t＝tx-tx-1；tx、tx-1分别为第x时段和第x-1时段。

此时，地表积水逐渐消失。这种情况下，随着入渗过程的进行，土壤积水全部渗入土壤，土壤开始进行非积水入渗过程。土壤入渗过程分为积水入渗过程和非积水入渗过程。

情景c(如图6所示)：h0>0，这种情况下，土壤继续进行积水入渗过程。

情景d(如图7所示)：h0＝0，这种情况下，土壤持续进行非积水入渗过程。其中，tx、tx-1分别为第x时段和第x-1时段，为湿润区实际平均导水系数。由于积水深h0与降雨强度单位不统一，无法直接代入模型计算。计算过程中，首先将积水深h0除以时段△t，转为对应深度的净雨强p’，然后代入模型计算，

其中p’＝h0/△t；△t＝tx-tx-1。

(五)降雨入渗产流过程计算的步骤：用于计算降雨入渗产流过程，包括：积水入渗过程计算、时段产流过程计算、非积水时段入渗过程计算、积水与非积水转换时刻计算。

降雨入渗产流过程计算主要有两种计算方法，分别是基于green-ampt模型和基于richards的方法。richards方法应用需要采用数值解计算，导致模型计算速度较慢，适合于小尺度的计算，如农田尺度。richards方法的计算公式如下：

式中：t为时间；q为水分通量；q为含水量；x、y和z表示坐标轴。

green-ampt模型计算速度快，适用于大尺度地区，如流域尺度，故本实施例采取这种方法，其计算过程如下：

所述的降雨入渗产流过程包括：积水入渗过程计算、时段产流过程计算、非积水时段入渗过程计算、积水与非积水转换时刻计算等四个过程：

(1)计算降雨入渗产流过程计算的步骤：用于计算降雨入渗产流过程，包括：积水入渗过程计算、产流过程计算、积水与非积水转换时刻计算。

所述的积水入渗过程计算包括：

湿润锋在第1层时土壤入渗过程计算：

a＝[swcf,1-0.5(hab,1+hwb,1)]δθ1

式中：为土壤累积入渗量；为积水发生时刻土壤累计入渗量；ka，1为第1层的实际导水系数；swcf，1为第1层中湿润锋处的土壤水吸力；hab，1为第1层的土壤进气值；hwb，1为第1层土壤的进水值；a为参数；dq为土壤含水变化量，加脚标表示某一层的土壤含水变化量。

其中，dq1＝qa，1－q0，1，qa，1为湿润区第1层土壤实际含水量，q0，1为第1层初始土壤含水量；

湿润锋在第m(2≤m≤n)层时，积水入渗过程计算：

式中，ka，m为第m层土壤的实际导水系数；bm－1、cm－1为参数；fm－1为m－1层以上的土壤累积入渗量；tm－1为湿润锋进入m－1和m层界面时的时间；swcf，i为第i层中湿润锋处的土壤水吸力；hab，i为第i层土壤的进气值；hwb，i为第i层土壤的进水值；li为第i层土壤的厚度；

其中，dqi＝qa，i－q0，i，qa，i为湿润区第i层土壤实际含水量，q0，i为第i层初始土壤含水量；i和m为土层数；

时段产流过程计算：

式中，为时段tx的累积地表产流；为时段tx的累积降雨；h0为时段tx的地表积水深；为时段tx的累积入渗量；

非积水时段入渗过程计算：

积水与非积水转换时刻tp计算：

当湿润锋在第1层时：

当湿润锋在第m(2≤m≤n)层时：

式中：tx-－1为时段x－1的区间时间；为时段tx－1时的累积入渗量；tp为土壤表层积水发生或消失时刻；为土壤表层积水发生或消失时刻时段降雨强度，fp为非积水时段与积水时段转换时刻的土壤入渗率。

实施五：

本实施例是实施例四改进，是实施例四关于湿润锋处土壤水吸力范围的细化。本实施例所述的swaf湿润锋处空气压力的最大值为：hb＝h0+l0+hab；

最小值为：hc＝h0+l0+hwb；

式中：hb为突破压力值；hc为闭合压力；l0为饱和区厚度；

由于饱和区的厚度很小，远小于hab和hwb，可以忽略不计，因此，在初始阶段和过渡阶段hb≈h0+hab，hc≈h0+hwb；稳定阶段时hb≈hc≈h0+0.5(hab+hwb)，假设初始阶段和过渡阶段，空气压力从突破压力到闭合压力的变化呈线性，则平均空气压力可以表示为：

swaf＝h0+0.5(hab+hwb)。

实施例六：

本实施例是实施例五的改进，是实施例五关于土壤进气值和土壤进水值的细化。本实施例所述的土壤进气值hab和土壤进水值hwb，土壤进气值hab＝1/αd，土壤进水值hwb＝hab/2－δ

其中：αd为brooks-corey模型中的参数，δ为经验值。

实施例七：

本实施例是实施例六的改进，是实施例六关于经验值δ的细化。本实施例所述的经验值δ的适用范围：砂土层的δ＝0～2cm，壤土中的δ＝2～5cm，黏土中的δ＝8～10cm。

实施例八：

本实施例是实施例四的改进，是实施例四关于饱和度系数sk、sθ的细化。针对土壤实际导水系数计算，第一层土壤i＝1中，所述的饱和度系数sk约等于土壤持水量qh与饱和土壤含水量qs之比：sk≈qh/qs；第i层土壤i>1中，所述饱和度系数sk约等于0.5。针对实际土壤含水量计算，将湿润锋以上区分为饱和区和非饱和区，饱和区sθ等于1，即饱和区实际含水量等于饱和含水量qa＝qs；在非饱和区，饱和度等于土壤持水量qh与饱和土壤含水量qs之比：sθ≈qh/qs，即：非饱和区实际含水量等于田间持水量qa＝qh。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案(比如土壤层的划分方式、各种公式的运用、步骤的先后顺序等)进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。