宽级配砾石土，也称为泥石流物源土，通常发育于泥石流沟谷两侧的斜坡，是由地震或滑坡等作用形成的第四纪松散的土体物质^[1–2]。宽级配砾石土的级配较宽，粒径跨越几个数量级，在降雨入渗的条件下易于产生优先流现象，即由于水在较大的孔隙内的入渗速率远大于在小孔隙中的入渗速率而使大孔隙部分成为了水量运输的高速通道，宏观上提高了整体的渗透性能。

降雨入渗的计算模型根据其适用条件可以分为两大类—均匀流和优先流。均匀流模型是建立在均匀介质假设条件下的，如Horton入渗模型、Philip入渗模型、Smith模型、Green-Ampt模型等。大尺度的土体中普遍存在优先流现象，与均匀流模型相比反映了土体具有更复杂的内部结构。大量研究证实了优先流现象的存在^[3–8]。描述优先流的传统模型主要是两域模型，两域模型中人为划分了基质区和优先流区，不能很好地体现土体自身的结构性。宽级配砾石土是一种多孔隙材料，大量研究表明天然岩土体等多孔材料的颗粒和孔隙都具有统计分形的特征，分形理论能够从整体上反映宽级配砾石土体多尺度的结构特性^[9–11]。Katz等^[12]采用扫描电镜测量一组砂岩样品，研究了孔隙率与分维值之间的关系。Oades等^[13]研究表明土具有明显的分级特征。Tyler等^[14]利用Sierpinski地毯模型研究了土体孔隙的分布，Bra-kensiek等^[15]研究土体的分形特性，得到了大孔隙的分布以幂函数表示及大孔隙的面孔隙度。Pachepsky等^[16]推导了分形多孔介质中的水分运动方程。Yu等^[17–18]研究孔隙个数满足的分形关系，推导了分维值的表达式。Yu等^[19]基于分形毛细管束模型提出非饱和多孔介质的相对渗透率的表达式。Guarracino^[20]采用Sierpinski地毯模型研究了裂隙化岩体的非饱和水力传导性。Xu等^[21]采用分形方法研究了各向同性多孔介质的Kozeny-Carman常数。Xu等^[22]以分形毛细管束模型为基础，考虑毛细管压力对多孔介质中水相与气体相的划分影响，推导出非饱和多孔介质相对渗透率。

作者采用分形理论描述宽级配砾石土的多尺度结构，建立一个能够反映优先流特征的降雨入渗模型，研究降雨作用下，入渗水量的分布规律及雨强对入渗特征的影响。

1 孔隙的结构和渗透性 1.1 孔隙的分形结构

在一个分形单元内，孔隙数量满足式（1）^[17]：

${N(\lambda) } = {\left( {\frac{{{\lambda _{\max }}}}{\lambda }} \right)^{{D_\text{f}}}}$

(1)

式中：λ为度量直径，在分形理论中，度量直径作为量测孔隙数量的尺子，即孔隙数量不是定值而是度量直径的函数；N(λ)为度量直径为λ时的孔隙数量；λ_max为最大的孔隙直径；D_f为孔隙面积的分维值。

孔隙的总数量为^[17]：

$N = {\left(\frac{{{\lambda _{\max }}}}{{{\lambda _{\min }}}}\right)^{{D_\text{f}}}}$

(2)

式中，λ_min为孔隙的最小直径。

由于孔隙数量巨大，式（1）可看做是连续可微的函数，在式（1）两边对λ取微分，得 ${\rm{ - d}}N \! = \! {D_\text{f}}\lambda _{\max }^{{D_\text{f}}}{\lambda ^{ - ({D_\text{f}} + 1)}}\text{d}\lambda $ 。孔隙分布的概率密度函数可以写成 $f(\lambda ) = {\rm{ }}\displaystyle\frac{{ - \text{d}N}}{{N\text{d}\lambda }} = $ ${D_\text{f}}\lambda _{\min }^{{D_\text{f}}}{\lambda ^{ - ({D_\text{f}} + 1)}}$ ，根据概率密度函数性质 $\displaystyle\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {f(\lambda )\text{d}\lambda } = $ $1 - {\displaystyle(\frac{{{\lambda _{\min }}}}{{{\lambda _{\max }}}})^{{D_\text{f}}}} \equiv 1$ ，得^[17]：

${\left(\frac{{{\lambda _{\min }}}}{{{\lambda _{\max }}}}\right)^{{D_\text{f}}}} = 0$

(（3）)

式（3）为分形理论是否可用的判别式。当最小直径与最大直径之间相差几个数量级时，可以认为式（3）近似满足，可以采用分形理论处理。

Yu等^[17]以Sierpinski模型为基础，得到分维值D_f与孔隙率 $\phi $ 的关系：

$\displaystyle{D_\text{f}} = \displaystyle2 - \frac{{\ln \phi }}{\displaystyle{\ln \frac{{{\lambda _{\min }}}}{{{\lambda _{\max }}}}}}$

(4)

作者以上述分形理论为基础，建立宽级配砾石土的结构。其中，最大直径与最小直径为分形结构的上下限直径，且假定在这个范围内分维值为常数。

1.2 渗透特征

以地面一点O为原点建立直角坐标系，见图1。竖直向下方向为z轴正方向；孔隙管道的实际长度为L_t，其在竖直方向上的投影长为L₀。入渗的深度越大，L_t越接近L₀；为简化模型，假定L_t=L₀，即所有孔隙管道都呈直线。

单个孔隙管道的流量采用Poiseuille等式（5）表示^[18]：

$q(\lambda ) = \frac{{\text{π} {\lambda ^4}\Delta p}}{{128\mu {L_0}}} = \frac{\text{π} }{{128}}\frac{\gamma }{\mu }{\lambda ^4}$

(5)

式中：λ为孔隙直径；μ为水的动力黏滞系数，在文中取温度20℃时的值1×10^-3 Pa·s ；γ为水的重度； $\displaystyle\frac{{\Delta p}}{{{L_0}}}$ 为水力坡降， $\displaystyle\frac{{\Delta p}}{{{L_0}}}$ 水在降雨入渗条件下仅考虑重力作用，取 $\displaystyle\frac{{\Delta p}}{{{L_0}}} = \gamma $ 。

一个分形单元的总流量可表达为：^[18]

$Q \! = \! \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} \! {q(\lambda )( - \text{d}N)} \! = \! \frac{\text{π} }{{128}}\frac{\gamma }{\mu }\displaystyle \frac{{{D_\text{f}}\lambda _{\max }^4}}{{4 - {D_\text{f}}}}\left[1 \! - \! {\left(\frac{{{\lambda _{\min }}}}{{{\lambda _{\max }}}}\right)^{4 - {D_\text{f}}}}\right]\! $

(6)

λ_min与λ_max相差几个数量级，所以式（6）可简化为^[18]：

$Q = \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {q(\lambda )( - \text{d}N)} = \frac{\text{π} }{{128}}\frac{\gamma }{\mu }\frac{{{D_\text{f}}\lambda _{\max }^4}}{{4 - {D_\text{f}}}}$

(7)

对于非饱和土，有部分孔隙管道中同时存在水和空气两相介质。将单个孔隙管道中水相部分所占的面积换算为等效的圆形面积，其直径用λ_ω表示， ${\lambda _\text{ω} } \le \lambda $ ，那么不完全充水的单个孔隙管道的流量仍用式（5）表示，其中孔隙直径λ取为 ${\lambda _\text{ω} }$ ^[19]。

2 孔隙结构分区

孔隙直径跨越多个量级，具有多尺度的特征，这使得降雨总入渗量不均匀地分散在各孔隙管道上。小尺度的孔隙管道入渗量小，大尺度的孔隙管道入渗量大。图2为降雨入渗示意图，降雨强度R为常数。在降雨过程中，存在流量特征直径λ_R( ${\lambda _{\min }} \le {\lambda _{\rm R}} \le {\lambda _{\max }}$ )。λ_R为完全充满水与不完全充满水的孔隙管道的分界直径，小于该特征直径的孔隙管道均饱和，流量均为各孔隙管道的最大流量，但各不相等；大于该特征直径的孔隙管道呈非饱和状态，这些孔隙管道的流量都相等，都等于直径为λ_R的孔隙管道的流量。从水量平衡的角度，小直径的孔隙表面产生剩余流量，直径相对较大的孔隙不完全被水充满，小孔隙的剩余流量首先流入大孔隙；当雨强大于土体饱和入渗速率时，所有大孔隙都被充满，此时 ${\lambda _{\rm R}} = {\lambda _{\max }}$ ，剩余流量在土体表面转化为超渗产流。

宏观土体面积可划分为若干个分形单元，1个分形单元的孔隙面积为：

$S = \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\frac{\text{π} }{4}{\lambda ^2}} ( - \text{d}N) = \frac{\text{π} }{4}\frac{{{D_\text{f}}\lambda _{\max }^2}}{{2 - {D_\text{f}}}}(1 - \phi )$

(8)

宏观上单位土体面积包含分形单元的个数n为：

$n = \frac{{\varOmega \phi }}{{\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\displaystyle\frac{\text{π} }{4}{\lambda ^2}( - \text{d}N)} }} = \frac{{4(2 - {D_\text{f}})}}{{\text{π} {D_\text{f}}\lambda _{\max }^2}}\frac{\phi }{{1 - \phi }}\varOmega $

(9)

式中，Ω为单位土体面积。

1个分形单元内孔隙面积（S）被划分为2个区域，基质区（S₁）和优先流区（S₂），见图3。

1）基质区（完全充水区）是指直径λ_min～λ_R的孔隙管道的区域，其面积为：

$\begin{aligned}[b]{S_1} = & \int_{{\lambda _{\min }}}^{{\lambda _{\rm{R}}}} {\frac{\text{π} }{4}{\lambda ^2}} ( - {\rm{d}}N) = \\[3pt] & \frac{\text{π}}{4}\frac{{{D_\text{f}}\lambda _{\max }^2}}{{2 - {D_\text{f}}}}\left[{\left(\frac{{{\lambda _\text{R}}}}{{{\lambda _{\max }}}}\right)^{2 - {D_\text{f}}}} - {\left(\frac{{{\lambda _{\min }}}}{{{\lambda _{\max }}}}\right)^{2 - {D_\text{f}}}}\right]\end{aligned}$

(10)

占孔隙面积的比例为：

$\frac{{{S_1}}}{S} = \frac{{{{\left(\displaystyle\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)}^{2 - {D_{\rm f}}}} - \phi }}{{1 - \phi }}$

(11)

2）优先流区（不完全充水区），指直径大于λ_R的孔隙管道所占的部分。其面积为：

${S_2} = \int_{{\lambda _{\rm R}}}^{{\lambda _{\max }}} {\frac{\text{π} }{4}{\lambda ^2}} ( - {\rm d}N) = \frac{\text{π} }{4} \frac{{{D_{\rm f}}\lambda _{\max }^2}}{{2 - {D_f}}}\left[1 - {\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)^{2 - {D_{\rm f}}}}\right]$

(12)

占孔隙面积的比例：

$\frac{{{S_2}}}{S} = \frac{{1 - {{\left(\displaystyle\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)}^{2 - {D_{\rm f}}}}}}{{1 - \phi }}$

(13)

优先流区内，每个孔隙管道内同时存在水和空气两相物质。其中水占总面积：

${S_{2\text{ω} }} \! = \! \! \int_{{\lambda _\text{R}}}^{{\lambda _{\max }}}\! {\frac{\text{π} }{4}\lambda _\text{R}^2( - \text{d}N)} \! = \! \frac{\text{π} }{4}\lambda _{\max }^2 \left[{\left(\frac{{{\lambda _\text{R}}}}{{{\lambda _{\max }}}}\right)^{2 - {D_\text{f}}}} \! - \! {\left(\frac{{{\lambda _\text{R}}}}{{{\lambda _{\max }}}}\right)^2}\right]\! \! $

(14)

空气部分占总面积：

${S_{2{\rm a}}} = {S_2} - {S_{2\text{ω} }}$

(15)

3 入渗量与入渗率的推导

入渗量是指在单位时间内一定面积上入渗水的方量，入渗率是指单位时间内单位面积上入渗水的深度，即竖直方向上的入渗速率。

3.1 入渗量

1）基质区内各不同直径孔隙管道的入渗量为式（5）。单位土体面积内基质区的总入渗量为：

$\begin{aligned}[b]{I_1} = & n\int_{{\lambda _{\rm{min}}}}^{{\lambda _{\rm R}}} {q(\lambda )( - {\rm d}N)} = \frac{1}{{32}}\frac{{2 - {D_{\rm f}}}}{{4 - {D_{\rm f}}}}\frac{{\phi \varOmega }}{{1 - \phi }}\frac{\gamma }{\mu }\lambda _{\max }^2 \times \\& \left[{\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)^{4 - {D_{\rm f}}}} - {\left(\frac{{{\lambda _{\rm{min}}}}}{{{\lambda _{\max }}}}\right)^{4 - {D_{\rm f}}}}\right] \end{aligned}$

(16)

式中，n为单位土体面积内包含的分形单元个数，Ω为单位土体面积， $q(\lambda )$ 为直径λ的孔隙管道流量。

2）优先流区内不同直径孔隙管道具有相同的入渗量，因此单位土体面积内优先流区的总入渗量为：

$\begin{aligned}[b]{I_2} = & n\int_{{\lambda _{\rm R}}}^{{\lambda _{\max }}} {q({\lambda _{\rm R}})} ( - {\rm d}N) \! = \! \frac{1}{{32}}\frac{{2 - {D_{\rm f}}}}{{{D_{\rm f}}}}\frac{{\phi \varOmega }}{{1 - \phi }}\frac{\gamma }{\mu }\lambda _{\max }^2 \times \! \\& \left[{\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)^{4 - {D_{\rm f}}}} - {\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)^4}\right]\! \end{aligned}$

(17)

式中， $q({\lambda _{\rm R}})$ 为等效直径为λ_R的孔隙管道流量。

单位土体面积上的总入渗量为：

$I = {I_1} + {I_2}$

(18)

当雨强超过土体饱和入渗速率时，单位土体面积上的超渗产流量为：

$\begin{aligned}[b]{Q_{\text{产流}}} = & R\varOmega - n\int_{{\lambda _{\rm{min}}}}^{{\lambda _{\max }}} {q(\lambda )} ( - {\rm d}N) = R\varOmega - \\& \frac{1}{{32}}\frac{{2 - {D_{\rm f}}}}{{4 - {D_{\rm f}}}}\frac{{\phi \varOmega }}{{1 - \phi }}\frac{\gamma }{\mu }\lambda _{\max }^2 \times \left[1 - {\left(\frac{{{\lambda _{\rm{min}}}}}{{{\lambda _{\max }}}}\right)^{4 - {D_{\rm f}}}} \right]\end{aligned}$

(19)

式中，R为降雨强度。

3.2 入渗率

1）基质区内各不同直径孔隙管道的入渗率为：

${i_1}(\lambda ) = \frac{{4q(\lambda )}}{{\text{π}{\lambda ^2}}} = \frac{1}{{32}}\frac{\gamma }{\mu }{\lambda ^2}$

(20)

2）优先流区内不同直径的孔隙管道具有相同的入渗率，因此该区域的入渗率与流量特征直径λ_R的孔隙管道的入渗率相等，为：

${i_2} = \frac{{\lambda _{\rm R}^2}}{{32}}\frac{\gamma }{\mu }$

(21)

显然，基质区各孔隙管道的入渗率均小于优先流区。因此，整个土体的入渗影响深度取决于优先流区的入渗率。

4 土体孔隙特征参数的确定

在上述分析中，孔隙的最大、最小直径和流量特征直径控制着各分区的比例。

4.1 最大直径和最小直径

最小孔隙通过试验方法很难测量，与最大孔隙相比小几个数量级，但一定不为0。最小直径是分形计算的下限直径，通常取经验值，本文模型中采用Katz等^[12]研究中的值，即2 nm。

孔隙的最大直径与土体渗透性的关系极为密切，但很难通过试验方法直接测量。从流量平衡角度，以土体宏观渗透系数表示最大直径的模型更为合理。仅考虑重力作用，水力坡降取为1，根据达西定律和流量平衡，最大直径满足式（22）^[23]：

$nQ = {K_{\rm{sat}}}\varOmega $

(22)

整理得：

${\lambda _{\max }} = \sqrt {32{K_{\rm{sat}}}\frac{{4 - {D_{\rm f}}}}{{2 - {D_{\rm f}}}}\frac{{1 - \phi }}{\phi }\frac{\mu }{\gamma }} $

(23)

式中， ${K_{\rm{sat}}}$ 为土体饱和渗透系数，Ω为单位土体面积，n为单位土体面积内包含分形单元的个数，Q为一个分形单元的总流量，μ为水的动力黏滞系数，γ为水的重度，D_f孔隙面积的分维值， $\phi $ 为孔隙率。

4.2 流量特征直径

根据流量平衡，流量特征直径λ_R满足：

${I_1} + {I_2} = R\varOmega $

(24)

整理得：

$\begin{aligned}[b]\frac{1}{{32}}\frac{{2 - {D_{\rm f}}}}{{{D_{\rm f}}}}\frac{\phi }{{1 - \phi }}\frac{\gamma }{\mu }\lambda _{\max }^2 \times\left[ {\frac{4}{{4 - {D_{\rm f}}}}{{\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)}^{4 - {D_{\rm f}}}} - } \right.\\\left. {\frac{{{D_{\rm f}}}}{{4 - {D_{\rm f}}}}{{\left(\frac{{{\lambda _{\rm{min}}}}}{{{\lambda _{\max }}}}\right)}^{4 - {D_{\rm f}}}} - {{\left(\frac{{{\lambda _{\rm R}}}}{{{\lambda _{\max }}}}\right)}^4}} \right] = R\end{aligned}$

(25)

其中，I₁、I₂分别为单位土体面积内基质区和优先流区的总流量，R为雨强，Ω为单位土体面积。

式（24）中，等式左边代表了从孔隙尺度上计算单位土体面积上的孔隙允许通过的流量，等式的右边代表了宏观尺度上降雨在单位土体面积上提供的流量。整理得到式（25），该式中仅包含流量特征直径一个未知量，可以通过迭代的方法求解。

5 雨强对渗透性的影响

现以四川理县滑坡滑体宽级配砾石土为例，研究雨强对入渗量和入渗率的影响。土体孔隙率0.3，饱和渗透系数15 mm/h。

雨强对特征直径的影响如图4所示。流量特征直径随雨强增大而增大，当雨强达到土体的饱和入渗速率时，流量特征直径达到最大值，即为孔隙最大直径，雨强再增大，流量特征直径保持不变。

雨强对孔隙各区域比例的影响如图5所示。随着雨强的增大，基质区比例逐渐增大，优先流区比例逐渐减小，整个土体饱和度增大。当雨强达到土体饱和入渗速率时，基质区面积达到最大，即为整个孔隙的面积，优先流区消失。雨强继续增强，基质区面积不变。

1 m²土体面积内总入渗量的变化情况如图6所示。随着雨强增大，基质区总入渗量逐渐增大，优先流区总入渗量随雨强的增大先逐渐增大；当达到一定值（11.2 mm/h）时，达到最大，再逐渐减小；直到雨强达到土体饱和入渗速率时，基质区入渗量达到最大，即为土体的饱和入渗量，同时优先流区消失。雨强再增大，剩余的流量转化为土体表面上的超渗产流。

不同雨强条件下，各孔隙管道的入渗量分布如图7所示。

不同雨强条件下，各孔隙管道的入渗率分布如图8所示。

6 降雨入渗不平衡性的机制

考虑孔隙的多尺度性，将土体孔隙划分为基质区和优先流区，从结构上反映了入渗不平衡性分布产生的根源及孔隙内部不同区域的入渗性质。均匀流模型将孔隙看作为均匀的入渗介质，将流量在全断面上平均，掩盖了孔隙结构的影响。

以Green-Ampt模型^[24]与本文优先流模型进行对比。从流量守恒的角度，两个模型的总入渗量相同，两模型的区别在于水量的分布形式。Green-Ampt模型是建立在均匀介质假设和有压入渗边界条件下的，其入渗存在明显的湿润锋面，入渗量在整个入渗断面上均匀分布，分布形式不随雨强变化而变化；优先流模型是建立在分形多尺度的孔隙结构和无压力自由入渗边界条件下，不存在湿润锋面，反映了入渗量按照孔隙直径的量级差异分布，并随雨强变化而变化，如图7所示。Green-Ampt模型的入渗率是全部水量在全断面上的平均入渗率，形式如下：

$i = \left\{ \begin{array}{l}\!\! \!R\text{，} R \le {v_{\rm s}}\text{；}\\\!\!\! {v_{\rm s}}\text{，}R > {v_{\rm s}}\end{array} \right.$

(26)

式中，i为入渗率，R为雨强，v_s为饱和入渗速率。

优先流模型中各孔隙管道的入渗率各不相同，优先流区的入渗率远大于基质区的入渗率，其分布如图8所示，且随雨强变化而变化。优先流模型入渗的影响深度取决于优先流区的入渗率。

2个模型中最大入渗率的差异如图9所示。

优先流模型的最大入渗率即为优先流区的入渗率。显然优先流模型计算的入渗影响深度远大于Green-Ampt均质模型。

优先流入渗模型体现了孔隙内部不同区域的入渗性质的差异，从土体结构上解释了降雨入渗不平衡性的分布机制。优先流区的入渗率决定了渗流能够到达深度的量级。对于解决降雨能否起动深层滑坡或农药能否快速污染地下水等问题，能够提供合理的入渗模型。

7 结 论

在土体分形孔隙结构的基础上，研究了降雨条件下宽级配砾石土的渗透特性，建立降雨入渗模型。

1）宽级配砾石土孔隙结构具有多尺度特性，在降雨作用下，不同尺度的孔隙具有不同渗流特征，以流量特征直径为界将孔隙部分划分为基质区和优先流区两个区域，流量特征直径随雨强的变化而变化。

2）基质区入渗率较小，优先流区高速入渗，体现了不平衡流特点，优先流区的入渗率决定了整个土体入渗的影响深度。

3）优先流现象的产生是由土体多尺度的孔隙结构和无压力的入渗边界条件决定的，其入渗特性主要受土体的饱和渗透系数和雨强的影响。