[发明专利]一种基于分层土壤的水沙过程计算方法

权利要求书

1.一种基于分层土壤的水沙过程计算方法，其特征在于，土壤进行改良，对改良后的土壤的水流进行汇流计算，对改良后的泥沙进行汇沙计算，来分析水循环过程和泥沙侵蚀过程；

所述水流的汇流计算包括单元内水分通量的计算：具体为水文气象数据展布的计算、植被截流的计算、枯枝落叶储流的计算、膨胀性土壤土石混合介质层水分运动的计算、沙沟及其以下土层水分运动的计算、地下水水层水分运动的计算、坡面汇流过程的计算和土壤蒸散发过程的计算；

所述泥沙的汇沙计算包括：土壤侵蚀的计算和河沟道泥沙输移过程的计算；所述土壤侵蚀的计算包括雨滴溅蚀的计算、薄层水流侵蚀的计算、股流侵蚀的计算和重力侵蚀的计算；

河沟道泥沙输移过程的计算包括沟道水流挟沙能力的计算和河道水流挟沙能力的计算；

所述水文气象数据展布的计算包括：

1)水文气象过程空间尺度展布

降水是驱动因素，而其他气象要素风速、相对湿度、大气温度、日照是计算蒸散发所必须的参数，以日为时间尺度进行计算，需要逐日气象数据作为输入；

气象数据作为输入值，采用外部气象展布算法将时间序列内的气象数据展布到计算单元内，实际运行期间直接读取对应数据；中采用泰森多边形法和反距离加权平均法进行流域内气象数据的展布，计算公式如下：

其中，D表示待插值点估计值，mm；D_pi表示第pi个参证站点数据，mm；m表示参证站点个数；λ_pi表示第pi个参证站点数据权重；d_pi表示第pi个参证站点同待插值点的距离，km；n表示权重指数，n＝0时，式(1)退化为算术平均法；n＝1时，式(1)为简单距离反比法；n＝2时式(1)为广泛应用的距离平方反比法；

2)降雨数据的时间降尺度展布

降雨期内，采用小时尺度进行计算，具体公式如下：

S＝a·P+b (5)

其中，i表示历时t内最大降水平均雨强；S表示暴雨参数或称之为雨力，等于单位时间内最大平均雨强；t表示时段；n表示暴雨衰减系数，与气候区有关，可用实测资料率定获取；P表示日降雨量；T表示日降雨总历时；a,b表示参数；

所述植被截流的计算：

忽略降雨期间蒸散发量条件下，植被截留采用下式计算：

Wr_max＝0.2·Veg·LAI (8)

其中，Veg表示植被的面积、盖度；Wr表示植被截留水量，mm；Wr_max表示最大植被截留水量，mm；P表示降雨量，mm；Rr表示植被冠层流出水量，即超出最大植被截留水量的部分，mm；LAI表示叶面积指数；

枯枝落叶储流的计算：

枯枝落叶储留量消耗于后期蒸散发过程，枯枝落叶储留采用下式计算：

H_umax＝α_maxG (9)

式中：G为枯枝落叶干重；α_max为枯枝落叶最大持水系数；

膨胀性粘土土石混合介质层的计算：

膨胀性粘土土石混合介质层水分运动改进本人提出的膨胀性土壤水分运动模型，该模型在分析土壤干湿受力变形的基础上，基于质量守恒原理得到了膨胀性土壤土石混合介质水分运动方程；

1)水分运动控制方程

基质区水分运动过程方程为：

上边界：

下边界：

裂隙优先流区水分运动过程方程为：

上边界：

下边界：

其中，Rv为碎石质量比系数；θ为土壤体积含水量，cm³/cm³；θ₁为时段初的土壤含水量，cm³/cm³；e为孔隙度，cm³/cm³；e₁为时段初的孔隙度，cm³/cm³；K_e(ψ)为土壤非饱和导水系数，cm/min；ψ为土壤水吸力，cm；S为源汇项，min^-1；W_e为两流区水量交换量，min^-1；Φ为土壤总水势，cm；q为水分通量，cm/min；上述字母的下标f、j分别表示裂隙优先流区、基质区；t为时间，min；z表示坐标轴z轴方向；w_f为裂隙优先流区面积比例，w_j为基质区面积比例；

2)土壤水分运动参数

土壤非饱和导水系数：采用改进的vanGenuchten模型计算膨胀性土壤裂隙优先流区、基质区非饱和导水系数，非饱和导水系数计算公式为：

K_e(ψ)＝K_sz(e)S_e^0.5[1-(1-S_e^1/m)^m]² (16)

其中，K_e(ψ)为膨胀性土壤非饱和导水系数，cm/min；K_sz(e)为深度为z时膨胀性土壤饱和导水系数，cm/min；S_e为饱和度；m和n为参数，m＝1-1/n；

受自重应力和膨胀力影响，土壤孔隙度随深度的变化而变化，导致土壤饱和导水系数随深度的变化而变化；土壤饱和导水系数随深度变化关系采用改进的Lambe模型计算：

e₀＝e_0ns·exp(-C_e·Rv) (19)

K₀＝K_ns·exp(-C_k·Rv) (20)

其中，e₀、K₀分别为零压强下的膨胀性土壤土石混合介质孔隙度、饱和导水系数，cm³/cm³、cm/min；e_0ns、K_ns分别为零压强下的无碎石膨胀性土壤孔隙度、饱和导水系数，cm³/cm³、cm/min；C_e为与孔隙度有关的土壤形状系数；C_k为与导水系数有关的土壤形状系数；e_z为深度为z时膨胀性土壤土石混合介质饱和孔隙度，cm³/cm³；m′为与土壤质地有关的参数；e₁为时段初的土壤孔隙度，cm³/cm³；ρ_d为密度，g/cm³；a₃为参数；α₃为土壤膨胀特征曲线斜率；U为质量含水量，g/g；A和B均为参数；γ为土壤湿比重，N/cm³；z为土壤深度，cm；

土壤膨胀特征曲线：膨胀性土壤土石混合介质变形主要受土壤膨胀力和自重应力影响，土壤膨胀力是土壤含水量的函数，采用三直线模型计算：

式中：ν为比容积，是膨胀性土壤土石混合容重的倒数；U为质量含水量；α₁、α₂、α₃为土壤膨胀特征曲线斜率；U_A、U_B、U_S分别为拐点处质量含水量；a、b和c为参数；

土壤应力应变关系曲线：土壤自重应力与土壤变形关系采用对数函数描述：

式中：ρ_s为膨胀性土壤土石混合介质容重；G为自重应力；γ为湿土比重；z为土壤深度；A和B为参数；

土壤饱和含水量：当土壤达到饱和时，土壤饱和含水量等于孔隙度，则任一深度饱和含水量可以表示为：

θ_sz＝e_z (23)

式中：θ_sz为深度为z时土壤饱和含水量；

模型计算中，基质区和大孔隙优先流区土壤水分特征曲线采用相同曲线；考虑到膨胀性土壤土石混合介质中含有大量碎石，改变了土壤水分特征参数；模型计算中，土壤水分特征曲线采用考虑碎石Van Genuchten模型计算：

式中：S_e为饱和度系数；R_v为碎石碎石质量比系数；θ为土石二元混合介质土壤含水量；θ_s为土石二元混合介质土壤饱和含水量；θ_r为土石二元混合介质土壤残余含水量；θ_n为无碎石土壤含水量；θ_ns为无碎石土壤饱和含水量；θ_nr为无碎石土壤残余含水量；α、n和m为参数，m＝1-1/n；h为土壤水吸力；由于土壤膨胀变形，θ_ns随深度的变化而变化，不同深度处θ_ns可以采用θ_sz代替；

3)两区面积比例

优先流区、基质区的面积比例，其计算公式为：

w_f＝de_w (26)

w_j＝1-de_w (27)

其中，e_w为由土壤吸水膨胀变形导致的孔隙度变化量，cm³/cm³，当土壤饱和时，w_f等于0；

4)土壤变形滞后性计算

土壤膨胀变形不仅受膨胀力和自重应力影响，也受变形时间影响，土壤膨胀变形随时间变化关系采用改进的Merchant模型描述：

式中：P为自重应力和膨胀力间的合力；α₁为弹簧a的体积压缩系数；α₂为弹簧b的体积压缩系数；η为Merchant模型中牛顿黏壶的黏滞系数；e₁为时段初孔隙度；e₂为时段末孔隙度；η、α₁和α₂通过土工试验确定；

沙沟及其以下土层水分运动的计算：

1)沙沟及其下包气带均采用二维Richards方程计算：

式中h为土壤水势，cm；S为源汇项，根系吸水项；下标i为土壤剖面分层数；x、z为坐标轴；Sr为源汇项，min^-1；K(h)为导水系数，cm/min；t为时间，min；

2)边界条件

上下边界均为通量边界：

左右为零通量边界：

式中：q为通量，cm/min；

地下水含水层水分运动的计算：

地下水运动相关计算公式如下：

其中，h表示地下水位，潜水层，m或水头，承压层，m；E表示蒸发蒸腾量，m；RG表示地下水出流量，m；C表示储留系数；k表示导水系数，m/day；z表示潜水层底部高程，m；GWP表示地下水开采量，m；Q₄表示来自不饱和土壤层的渗透量，m；

坡面汇流过程的计算，包括坡面汇流和河道汇流：

汇流过程采用运动波方程进行模拟，公式如下：

S_f＝S₀(运动方程) (34)

式中：Q表示过流断面流量，cm³/s；A表示过流断面面积，cm²；q_L表示单宽入流量，计算单元或河道所有流入的水量，cm³/s/m；S_f表示摩擦坡降，比例系数，tan(α)，α表示坡面和水平面的夹角；S₀表示计算单元平均地面坡降或河道坡降，比例系数；R表示过流断面水力半径，cm；n表示曼宁糙率系数；

土壤蒸散发的计算：

计算单元内的土壤蒸散发量是植被截留蒸发、土壤蒸发和植被蒸腾的之和，计算公式如下：

E＝E₁+E₂+E₃ (36)

式中：E表示计算单元总蒸散发(mm)；下标1表示植被截留蒸发；下标2表示植被蒸腾；下标3表示裸土蒸发；

土壤潜在蒸发能力由Penman公式计算，最大蒸发强度：

式中：RN为净辐射量；G为传入水中的热通量；Δ为饱和水汽压对温度的导数；ρ_a为空气密度；C_p为空气的定压比热；δ_e为实际水蒸气压与饱和水蒸气压的差值；r_a为蒸发表面空气动力学阻抗；λ为水的气化潜热；γ为空气湿度常数；PR为大气压；

空气动力学阻抗计算：其计算公式如下：

式中：r_a为空气动力学阻抗，s/m；hz为气象站观测点离地面的高度，m；hd为置换高度，m；z_om表为水蒸气紊流扩散对应的地表粗度，m；z_ox为地表粗度；κ为von Karman常数；U为风速；h_c为植被高度；

植被截留蒸发量E₁使用Noilhan-Planton模型计算：

E₁＝Veg·δ·E_p (41)

δ＝(Wr/Wr_max)^2/3 (44)

Wr_max＝0.2·Veg·LAI (45)

式中：Veg为裸地-植被域中植被的面积占计算单元的面积比例；δ为湿润叶面占植被叶面的面积比例；E_p为潜在蒸发量，即最大蒸发量；Wr为植被截留量；Wr_max为植被最大截留水量；I为雨强；Rr为植被冠层流出水量，即超出最大植被截留水量的部分，mm；LAI表示叶面积指数；

植被蒸腾量E₂采用Penman-Monteith公式计算，具体公式如下：

E₂＝Veg·(1-δ)·E_PM (46)

式中：E₂为实际植被蒸腾量，mm；E_PM为采用Penman-Monteith公式计算的潜在蒸腾量；G为热通量；r_c为植被阻抗，s/m；

水流控制方程中的源汇项S_r表示植物根系在单位时间内从土壤中吸收的水分，采用Feddes模型计算：

S_r(h)＝α(h)·E₂

E₂＝b(x,z)·E_PM (48)

其中：

式中：α(h)为土壤水势的指定响应函数，0≤α≤1；E₂为潜在吸水速率，min^-1；b(x,z)为标准化二维根系吸水分布，cm^-2；E_PM为潜在蒸腾速率，cm/min；X_m、Z_m分别为x方向和z方向最大根系长度，cm；p_z、p_r、x^*和z^*为经验参数；

植被群落阻抗计算：采用Dickinson提出公式计算植被群落阻抗：

σ₁^-1＝1-0.0016(25-T_a)² (51)

σ₂^-1＝1-VPD/VPD_c (52)

式中：r_c为植被群落阻抗，s/m；r_smin为最小气孔阻抗，s/m；LAI为叶面积指数；σ₁为温度影响函数；σ₂为大气水蒸气压饱和差影响函数；σ₃为光合作用有效放射的影响函数；σ₄为土壤含水量的影响函数；T_a为气温，℃；VPD为饱和水蒸气压同实测值之间的差，kPa；VPD_c为气孔闭合时的VPD值，等于4kPa；r_smax为最大气孔阻抗，5000s/m；PAR为光和作用有效放射，W/m²；PAR_c为PAR的临界值，高植被：30W/m²；低植被：100W/m²；θ为根系层土壤含水量，土壤含水量均指体积含水量；θ_w为植被凋萎时的土壤含水量；θ_s为饱和土壤含水量；

裸地土壤蒸发量E₃由下式计算：

式中：β为土壤湿润函数；θ为表层土壤含水量；θ_m为土壤分子吸力，1000-10000个大气压对应的土壤含水量；θ_h为表层土壤田间持水量；

土壤侵蚀计算主要包括雨滴溅蚀过程、薄层水流侵蚀过程、股流侵蚀过程和重力侵蚀过程的计算；其中面蚀过程和细沟侵蚀适用于薄层水流侵蚀模拟，浅沟侵蚀和切沟侵蚀适用于股流侵蚀；

雨滴溅蚀的计算

雨滴溅蚀是降雨雨滴击打地表，造成土壤颗粒分散并以击溅跃移形式搬运的侵蚀方式，是坡面侵蚀的重要来源；雨滴溅蚀主要受降雨强度、雨滴终点速度、雨滴直径土壤特性、地表覆盖物、地面坡度、风速因素影响；其中，降雨动能是雨滴溅蚀驱动因子；

假定雨滴溅蚀在整个计算单元内均有发生，对于计算单元内不同土地利用方式雨滴溅蚀，分别给出相对于裸地的溅蚀衰减系数，其中裸地土壤雨滴溅蚀采用吴普特模型：

式中：D₁为雨滴击溅侵蚀量，g/m²；I为雨强，mm/min；E为雨滴动能，J/m²；J₁为地表坡度，°；k₁、α₁、β₁为经验系数；式中降雨动能E采用江善忠提出的计算公式：

式中：I为雨强，mm/min；k′₁、α′₁为经验系数；

Guy等研究表明，雨滴击打作用可以增强薄层水流侵蚀能力；由雨滴溅蚀增加的土壤侵蚀能力计算公式为：

式中：q_sl为单宽流量输沙能力，kg/m²；k₂、α₂和β₂为经验系数；其它参数意义同上；

雨滴溅蚀与水深呈负相关关系，当水层增加到一定深度，雨滴动能将全部消耗于击穿水层而不能产生溅蚀；水深超过3倍雨滴直径时，雨蚀作用消失；当水深大于雨滴直径3倍以上时，即水深为0.6cm时，雨滴溅蚀作用消失；

薄层水流侵蚀的计算

本文研究采用Foster模型计算薄层水流侵蚀过程，具体如下：

D_c＝k₃(τ_f-τ_c) (60)

式中：D_c为水流剥离土壤速率，kg/(m²·s)；τ_f为水流对土壤颗粒的剪切应力，kg/m²；τ_c为土壤临界抗剪切应力，kg/m²；k₃为土壤可蚀性系数，1/s；D_r为细沟水流剥蚀率，kg/(m²·s)；q为单宽流量，m³/s；c为水流泥沙含量，kg/m³；T_c为水流挟沙能力；k₄、α₄为经验系数；

股流侵蚀过程的计算

股流侵蚀分为浅沟侵蚀和切沟侵蚀，其中浅沟是细沟至切沟的一种过渡形式；股流侵蚀采用实验概化模型：

式中：D_E为股流侵蚀量，kg；k₅为浅沟水流挟沙能力系数；m为经验系数；Q_E为单宽流量，m³/s；T_SE为股流挟沙能力；D_r为上游来沙量，kg；ω_u水流功率，m/s；

重力侵蚀的计算

重力侵蚀是指坡面岩体、土体在重力作用下，失去平衡而发生位移的过程；重力侵蚀包括泻溜、滑塌、崩塌、泥石流等多种类型；重力侵蚀通常发生在黄土坡面沟谷区，且其发生具有很大的随机性，时空变化很大；降雨产生的径流是重力侵蚀产沙的主要动力；影响重力侵蚀比较大的因素是沟谷长度、沟谷深度、沟坡坡度和土壤含水量；由于重力侵蚀影响因素复杂，发生具有随机性和偶然性，当前对重力侵蚀机理认识尚处于完善阶段，大多数重力侵蚀模型均概化侵蚀过程；采用了龚家国概化的重力侵蚀模型模拟重力侵蚀过程，其具体计算公式如下：

τ_c＝c+σtanφ+mq_stanφ (67)

式中：T_g为重力侵蚀量，g/m²；Δh为土壤重力侵蚀深度，m，切沟沟坡取10m，浅沟沟坡取0.5m，细沟沟坡取0.1m；θ为土壤天然休止角；γ′_s为土壤湿容重，g/cm³；k为发生重力侵蚀沟长系数；T_Gx为土体下滑剪切力，kg/m²；τ_c土体抗剪强度，kg/m²；c为土壤原始凝聚力，kg/m²；σ为土体剪切面法向应力，kg/m²；φ为土体内摩擦角，°；m为不稳定凝聚力与结构强度的比值；q_s为非饱和黄土结构强度,a、λ为经验系数；W为土壤含水量，cm³/cm³；L_gully为沟道长度，m；

河沟道泥沙输移过程的计算：

河沟道泥沙过程采用一维恒定水流泥沙扩散方程：

式中：S_x为断面平均含沙量,kg/m3；为断面水流挟沙能力；q为单宽流量，m³/s；α为恢复饱和系数；ω_s为河流泥沙沉速，m/s；S为断面出口平均含沙量，kg/m³；S₀为出口断面平均含沙量，kg/m³；S_0*为进口断面水流挟沙能力；S_*为断面出口水流挟沙能力；L为河段长度，m；d₉₀为泥沙上限粒径，μm；ρ_m为含沙水流密度，kg/m³；μ为浑水粘度；S_vm为基线体积比含沙量，m³/m³；S_v为体积比含沙量；μ₀为水的粘度；k₈为泥沙固体浓度修正系数；d_i为某一粒径级级平均直径，μm；P_i为某一粒径级重量百分比；

河道断面水流挟沙能力计算分为沟道和河道两种分别计算：

(1)沟道水流挟沙能力

沟道水流挟沙能力采用费祥俊模型：

式中：S_v为沟道水流含沙量，kg/m³；U为断面平均流速，m³/s；f为达西系数；R为水力半径，m；ω₉₀为上限粒径在一定浓度下的沉速，m/s；k_s为沟道粗糙度，取d₉₀；R_e为雷诺数；γ_m为含沙水流容重，kg/m³；

g为重力加速度

(2)河道水流挟沙能力

河道水流挟沙能力采用张洪武模型：

式中：S_*为水流挟沙能力；S_v为体积含沙量，kg/m³；v为断面流速，m³/s；k为浑水卡门常数，取0.4；γ_s为泥沙容重，kg/m³；γ_m为浑水容重，kg/m³；g为重力加速度，n/kg；h为断面平均水深，m；D₅₀为泥沙级配中值粒径，μm；土壤进行改良方法具体包括在原始土壤上挖沙沟，沙沟内填沙，填完沙后的沙沟上覆盖膨胀性土壤土石混合介质层；

所述膨胀性土壤土石混合介质层包括膨胀性粘土和碎石；

所述的沙沟由沙土组成，所述沙沟横截面呈椭圆形，长直径为0.5-1m，短直径为0.3-0.5m，沙沟的纵截面为弧形；

两个相邻沙沟形成波浪线型，两个相邻沙沟之间的原始土壤上种植有植被。

2.根据权利要求1所述的一种基于分层土壤的水沙过程计算方法，其特征在于，膨胀性土壤土石混合介质层的制备方法包括：直接就地取土，粉碎后过筛，获得粘土，然后与膨润土按比例混合，并搅拌均匀，最后加入碎石，获得含碎石的膨胀性粘土；所述膨胀性土壤土石混合介质层的厚度为10cm；所述的碎石粒径大于5mm，小于50mm；所述胀性粘土过筛后获得的粘土小于2μm；所述粘土与膨润土的质量比为0-1:5；所述碎石的质量占膨胀性土壤土石混合介质层总质量的20-40％；所述膨润土中的矿物成分为蒙脱石或蛭石；所述沙沟和膨胀性土壤土石混合介质层80-90％区域的总厚度超过15cm，沙土为河道采集的细沙或人工制备的细沙。