一种适用于任意初始含水量分布下的边坡降雨入渗-径流计算方法

技术领域

本发明涉及到一种适用于任意初始含水量分布下的边坡降雨入渗-径流计算方法，属于土壤物理和山坡水文学技术领域。

背景技术

降雨滴落至土壤表层时，会入渗进入土壤，直到降雨强度超过土壤有限的入渗能力，此时土壤表层产生积水，径流和侵蚀也随之发生。准确的估计积水发生时间以及径流产生量对于水文预报、农业灌溉设计、水土保持至关重要。

降雨入渗是一个复杂的动态过程，主要的影响因素包括降雨强度随时间的变化、土壤剖面的水力特性参数、剖面上初始含水量的分布状态以及坡度等。当前主要采用有限元和有限差分法数值求解Richards方程来估计积水发生时间和径流产生量，但数值模拟方法往往需要在时间和空间上进行高度的离散化来适应上边界复杂的降雨强度变化。此外求解一个高度非线性的Richards方程需要消耗大量的算力，因此这类方法的模拟过程往往面临着数值运算不稳定、计算成本较高等问题。作为一种改进方法，经验、解析、半解析的入渗方程因其形式简单，不需要时空离散化，和保证了水量平衡而被发现，并应用于入渗-径流过程的模拟。但由于降雨入渗问题的复杂性，现存的入渗方程要么忽略了降雨强度随时间变化这一自然特征(CN 105547957；CN 112685874；CN 109898489)，亦或是本身具有的初始含水量均匀分布的假定无法满足田间初始含水量呈非均匀分布的条件(Chen and Young,2006；Wang et.al,2018),均难以真正应用于真实的野外环境。故目前仍缺乏一种可在野外使用的，能够模拟非稳态降雨下初始含水量呈非均匀分布的边坡入渗-径流计算方法，这将极大地降低水文预报的时效性与准确性。

发明内容

本发明的目的是克服现有技术的不足，提供一种适用于任意初始含水量分布下的边坡降雨入渗-径流计算方法。

本计算方法所预先需要测定的模型参数包括：1、土壤的初始含水量分布剖面θ

步骤1，对于已知的将土壤的初始含水量分布剖面θ

式中，Z

步骤2，根据步骤1的分割剖面，土壤累积入渗量I被表示为黎曼和的形式：

式(13)中，J

利用一阶泰勒级数展开式和中心差分方法，θ

结合公式(13)和公式(14)计算得到ΔI

对于已知的ΔI

对于已知的Δt

进而，对于给定的土壤的初始含水量分布剖面θ

步骤3，对于步骤2所得到的离散的土壤入渗能力曲线，采用Philip二项边坡入渗公式进行参数α、β的拟合，进而获得一条连续的表层土壤入渗率曲线J

J

式中，γ为坡度，根据质量守恒原理，一条连续的累积入渗量曲线I(t)表示为，

步骤4，根据降雨监测分辨率的大小，将非稳态降雨事件划分为分段式的稳态降雨区间，在第n段稳态降雨区间内的降雨速率r

r

式中，t

步骤5，对于第一段稳态降雨，降雨速率为r

此时的土壤入渗能力通过如下公式描述，

J

累积入渗量I(t

I(t

累积径流量F(t

F(t

步骤6，基于时间压缩方法，计算第n-1段末(n>2)的土壤入渗能力J

I[J

步骤7，根据比较土壤入渗能力J

a).当J

I(t

累积径流量F(t

F(t

b).当J

r

如果t

I(t

F(t

相反的，如果t

J

累积入渗量I(t

I(t

F(t

作为本发明的一种优选技术方案：所述步骤2中，公式(17)中的J

执行不动点迭代，获得湿润区中饱和区所占的长度z

获得J

式中未知参数的定义如下：

作为本发明的一种优选技术方案：所述步骤2中，由于在入渗的初始时刻，公式(17)中表层土壤入渗率J

式中，ΔH＝H

作为本发明的一种优选技术方案：所述步骤6中，根据土壤的饱和导水率K

/>

执行不动点迭代算法获得第n-1段的降雨区间末的虚拟时间t

作为本发明的一种优选技术方案：所述步骤7中，将公式(32)代入公式(21)中，联立计算得到公式：

将上式中的I替换为I(t

作为本发明的一种优选技术方案：所述步骤7中，公式(22)中的参数t

作为本发明的一种优选技术方案：所述步骤1中，N值至少在101以上。

本发明的有益效果如下：

1、与已有技术相比，本发明考虑了降雨强度的非稳态性、土壤初始含水量呈非均匀分布等野外真实条件，基于积水条件下的边坡入渗方程，一种可以模拟任意初始含水量分布条件下的降雨入渗解析计算方法被提出，准确估计了均质土的积水发生时间与径流产生量，这为水文预报、农业灌溉设计、水土保持等领域提供了一个有效的计算工具；

2、计算方法简单。本发明仅需已知待计算土壤剖面的水力学性质、土壤初始含水量的分布以及上边界的降雨强度随时间的变化数据，便可计算得出积水发生时间与径流产生量，其计算过程不需要进行复杂的时间与空间的离散化，计算稳定、快速，且保证水量平衡，是一种简单，准确的计算方法。

附图说明

图1为本发明的计算流程图；

图2为积水条件下的边坡入渗水分剖面概念图；

图3(a)为对于一个任意初始含水量剖面的分割概念图，(b)为虚拟的活塞型剖面在任意初始含水量剖面中的入渗概念图；

图4是在两种典型土壤中，非稳态降雨条件下的边坡入渗-径流模拟结果，图中(a)为干旱、半干旱区典型的V字型初始含水量分布剖面，(b)为在V字型初始含水量分布剖面下的沙壤土降雨入渗对比图，(c)为湿润区典型的存在浅层地下水时的初始含水量分布剖面，(d)为在存在浅层地下水时的初始含水量分布剖面下的粉沙土降雨入渗对比图。

具体实施方式

以下结合附图，进一步说明本发明的内容，但不应理解为对本发明的限制。在不背离本发明精神和实质的情况下，对本发明方法、步骤或条件所作的修改和替换，均属于本发明的范围。若未特别指明，以下实施案例中所用的技术手段为本领域技术人员所熟知的常规手段。图1说明了本发明的计算过程。

实施例1

技术原理：首先分析描述积水条件下的Richards边坡水流控制方程及其边界条件如下：

式中：θ是土壤含水量(cm

并分析Brooks-Corey模型描述土壤非饱和导水率K(h)、土壤基质势和土壤含水量之间关系的表达式分别如下：

式中：S是有效土壤水饱和度，θ

通过最小作用量原理，通量假定以及积分中值定理，获得了描述土壤水分上渗的剖面方程，其中，采用的通量假定与积分中值定理是用以简化剖面方程的函数形式，其中通量假定的表达式为：

式中：J

再基于质量守恒原理与剖面方程，推导得出完整的描述恒定水头边界条件下带有虚拟活塞型剖面的边坡水分入渗的解析模型(5)–(9)(图2)。

/>

式中，

式中，z

基于上述的入渗解析模型，一种计算初始含水量非均匀分布条件下的边坡积水入渗方法被建立。方法成立的关键在于：当虚拟的活塞型的湿润区长度z

实施例2

基于实施例1，一种适用于任意初始含水量分布下的边坡降雨入渗-径流计算方法被提出，具体技术方案主要包括两个部分：1、对于给定的初始含水量非均匀分布的土壤剖面，采用一种计算初始含水量非均匀分布条件下的边坡积水入渗方法得到连续的土壤入渗率曲线J

1、计算初始含水量非均匀分布条件下的边坡积水入渗方法

计算方法所需的模型参数包括：1、土壤的初始含水量分布剖面θ

首先，给予任意分布的土壤初始含水量剖面θ

式中，Z

其次，根据上述的分割剖面，土壤累积入渗量被表示为黎曼和的形式，

式中，J

结合公式(13)和公式(14)可以计算得到ΔI

对于已知的ΔI

式中J

对上式执行不动点迭代，获得获得z

综上，对于一个特定的初始含水量分布剖面θ

最后，对于上述步骤所得到的离散的土壤入渗能力曲线，采用Philip二项边坡入渗公式进行参数拟合(拟合参数α、β)，进而获得一条连续的表层土壤入渗率曲线，

J

根据质量守恒原理，积分得到累积入渗量与时间之间的函数表达式，

2、模拟非稳态降雨条件下的边坡降雨入渗-径流方法

基于上述方法计算得到的累积入渗量曲线I(t)和土壤入渗率曲线J

首先，对于一场稳态降雨事件，降雨速率为rcosγ(cm min

J

式中，t

其次，根据降雨监测分辨率的大小，可以将非稳态降雨事件划分为分段式的稳态降雨区间(Chu,1978)。在第在n段稳态降雨区间内的降雨速率r

r

式中，t

第一段稳态降雨区间内，降雨速率为r

I(t

F(t

基于时间压缩方法，在第n-1段末(n>2)的土壤入渗能力J

I[J

式中，I(t

进而J

最后，根据比较土壤入渗能力J

a).当J

I(t

累积径流量F(t

F(t

b).当J

r

具体计算步骤为，将公式(33)代入公式(20)与(21)中，联立计算得到，

如果t

I(t

F(t

相反的，如果t

I(t

F(t

实施例3

采用以下不同土壤质地(土壤水力学性质见表1)、不同非稳态降雨事件以及两种典型的野外土壤初始含水量分布剖面(见图4a、c)进一步说明本发明公布的一种适用于任意初始含水量分布下的边坡降雨入渗-径流计算方法。

表1两种典型质地土壤的水力学性质

注：表中参数意义同上。

1、图(a)展示了V字型的初始含水量剖面，图(b)展示了非稳态降雨强度随时间的变化数据，土壤类型为砂壤土，坡度设置为30°。

步骤一，按照公式(a)、(b)将V字型的初始含水量剖面进行等距分割，其ΔZ

步骤二：按照公式(c)，(f)，(g)，(h)计算该分割后的水量剖面获得一条离散化的土壤入渗能力曲线J

步骤三：利用公式(i)对J

步骤四：根据降雨监测分辨率的大小，非稳态降雨事件划分为时间分辨率一小时一段的稳态降雨区间，见图4(b)。

步骤5：对于非稳态降雨事件的第一段0-60min，其降雨强度r

对于剩余三段稳态降雨区间，步骤6与步骤7将会被循环执行3次，来计算土壤的累积入渗量I(t)与累积径流量F(t)。

①对于非稳态降雨事件的第二段60-120min，其降雨强度为r

②对于非稳态降雨事件的第三段120-180min，其降雨强度为r

③对于非稳态降雨事件的第四段180-240min，其降雨强度为r

2、图(c)展示了存在浅层地下水时的初始含水量剖面，图(d)展示了非稳态降雨强度随时间的变化数据，土壤类型为粉沙土，坡度设置为30°。

步骤一，按照公式(a)、(b)将V字型的初始含水量剖面进行等距分割，其ΔZ

步骤二：按照公式(c)，(f)，(g)，(h)计算该分割后的水量剖面获得一条离散化的土壤入渗能力曲线J

步骤三：利用公式(i)对J

步骤四：根据降雨监测分辨率的大小，非稳态降雨事件划分为时间分辨率一小时一段的稳态降雨区间，见图4(d)。

步骤5：对于非稳态降雨事件的第一段0-60min，其降雨强度r

对于剩余三段稳态降雨区间，步骤6与步骤7将会被循环执行3次，来计算土壤的累积入渗量I(t)与累积径流量F(t)。

①对于非稳态降雨事件的第二段60-120min，其降雨强度为r

②对于非稳态降雨事件的第三段120-180min，其降雨强度为r

③对于非稳态降雨事件的第四段180-240min，其降雨强度为r

图4分别展示了在两种典型的边坡降雨入渗-径流中解析计算方法模拟的土壤入渗能力曲线与直接对Richards方程进行数值求解的差异比较。从图中可以看出解析解与观测结果几乎是相一致的，二者间的相对误差分别为1.81％和1.65％，因此该边坡降雨入渗-径流计算模型可以准确模拟不同质地、不同降雨强度以及任意初始含水量下的积水发生时间和径流产生量。

最后所应说明的是，以上具体实施方式仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对发明进行了详细说明，本领域的普通技术人员应当理解，可以对发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。