黄土高原是中国发生滑坡等地质灾害频繁的地区之一^[1]。据统计，中国每年约30%的地质灾害发生在占国土面积6%的黄土高原上^[2]。受季风影响，降雨是黄土高原滑坡的主要触发因素^[3]。据统计，降水诱发黄土高原滑坡的比例占滑坡灾害的70%以上^[4-5]。2013年，甘肃天水地区“7·25”暴雨引发地质灾害达700多处，其中浅层滑坡灾害点占58%^[6]。2018年7月，连续17 d的降雨造成甘肃庆阳市东庄沟小流域发生35场浅层滑坡，导致土体下移至坡脚，大量泥沙冲入河道，而这些滑坡均发生在有植被覆盖的坡面上^[7]。

对于黄土坡面，短历时强降雨或连续降雨造成土体基质吸力减小，以及孔隙水压力增加，使得土体力学强度降低，进而影响坡体稳定性^[8]。在这一过程中，土体的纵向渗透与侧向排水能力决定了坡体是否因水分累积而形成局部超载失稳现象^[9]。土壤颗粒的物理性质，如孔隙度、颗粒组成、饱和导水率等影响了土体的水分渗透和排水能力。浅层土体的物理性质会受到植被根系发展和生物活动等影响，产生异质性^[10]。随着1999年退耕还林工程的实施，2000—2020年黄土高原的植被覆盖度由0.39提高到0.61，植被恢复区域占黄土高原总面积的21.74%^[11-13]。大规模退耕还林使浅层土壤物理性质产生了极大的变化^[14-15]，除了减水减蚀^[16]，还改变了山坡-河谷的水文过程^[17]。植被恢复导致近地表土壤孔隙度和饱和导水率的增加^[18]，提高了降雨入渗能力，降低了地表产流能力。最新研究显示^[19-21]，植被恢复有可能改变黄土高原小流域的产流机制，使其由超渗产流机制为主导转变为由上层滞水引起的浅层蓄满产流机制为主导。这种改变的主要原因是植被恢复引起的土壤饱和导水率纵向异质性。研究发现，在植被恢复的土地上，土壤饱和导水率存在显著的纵向异质性，近地表层土壤的饱和导水率通常最大^[14]。同时，土壤饱和导水率纵向异质性的特征有利于层间水分积累^[17]。从坡体稳定性的角度看，入渗能力的增强和近地表层间水分的累积均暂时增大了坡体的自重，在一定程度上会增加坡体失稳的风险。Wang等^[22]通过对延安市“2013·7”暴雨后的8 135个浅层滑坡的采样分析，发现黄土滑动面深度与入渗深度、坡度和降雨强度等因素紧密相关，其中入渗深度(一般不超过2 m)与植被覆盖类型有关。许阳光等^[23]通过野外原位模拟降雨试验发现，在极端降雨条件下，草被沟坡易发生浅层滑坡，滑坡面上、下层土壤的物理性质存在显著差异，特别是上层饱和导水率可达下层饱和导水率的8.2~311.0倍。Guo等^[7]通过对黄土高原南小河沟流域的35个浅层滑坡野外调查，发现滑动面上下层的土壤饱和导水率差异巨大，进而导致下渗受阻。在这个过程中，植被可能对浅层滑坡的发生起了不可忽视的促进作用，应当引起重视。一些学者认为在植被恢复的过程中，根系会增加土壤含水率，在黄土高原降雨过程中易诱发滑坡的发生^[7]。在强降雨作用下，植被及根系会加剧滑坡侵蚀^[24]。相比于无植被覆盖的边坡，黄土高原有植被覆盖的坡面浅层滑坡总量增加153%^[25]，而传统国内外研究认为，植被及其根系对土体稳定性具有明显的促进作用^{[5, 26]}。以上研究通过实验和野外调查的方法初步揭示了黄土高原植被恢复与浅层滑坡的潜在相关性，然而尚未揭示植被恢复类型及程度对黄土高原浅层滑坡的影响机理，特别是植被恢复造成的浅层土壤纵向异质性对滑坡的发生、滑动面深度的影响。

地形、土层厚度及水力性质、降雨强度及持续时间、植被特征、前期水分条件等会对近地表水文响应产生综合影响^[9]。当入渗和侧向流动累积超过排水能力时，基质吸力和孔隙压力的骤变会引发边坡的破坏。滑坡的水文过程变量难以通过事后野外调查或原位实验获得，因而现象与机理之间存在难以逾越的鸿沟^[27]。为了解决这个问题，本文研究使用基于物理的建模方法，地下水流动物理方程的数值解近似地反映了孔隙压力在空间和时间上对降雨的响应，将地下水流动物理方程的数值解与边坡稳定分析法相结合，以评估研究区域的安全系数或边坡失效的可能性，并揭示水文响应过程与边坡失稳的内在联系。降雨型浅层滑坡作为一种水分累积运动导致坡体超载失稳的滑坡类型，对其水文响应过程的描述越精确，滑坡预测效果越好^[28]。基于此，本文通过将基于物理的分布式综合水文模型(integrated hydrology model, InHM)与无限边坡稳定分析(infinite slope stability analysis, ISA)模型进行耦合，基于不同植被恢复条件下的饱和导水率实测数据，研究黄土高原植被恢复背景下饱和导水率纵向异质性对浅层坡体稳定性的影响，为黄土高原退耕还林还草的水土保持效益提供建议，也为未来黄土高原地区的植被恢复和水土保持工作提供理论参考。

1 研究区域 1.1 流域概况

本文研究选取黄土高原多沙粗沙区的佘家沟流域(109°58′ E, 37°42′ N)作为研究对象。佘家沟流域位于陕西省榆林市子洲县，属于无定河岔巴沟(180 km²)的子流域，如图 1a所示。研究区处于吴起—志丹—延安—安塞滑坡灾害发育区，受青藏高原的影响，极易形成气旋，造成该地区多暴雨，地质灾害多为多日降水诱发型的揭皮式滑坡^{[1, 23]}。佘家沟流域面积为4.26 km²，海拔范围为920~1 143 m，河道平均比降约达1.15%。基于1959—2006年的观测数据，该流域年平均降水和潜在蒸发分别为480和1 570 mm，短历时高强度降雨集中发生在汛期，约占全年降水量的70%^[29]。研究区位置、地形和站点分布，如图 1所示。图中DEM(digital elevation model) 指数字高程模型。佘家沟流域汛期主要产流模式为超渗产流，流域在短历时、高强度暴雨条件下快速响应，产生大量地表径流^[30]。根据USDA(U.S. Department of Agriculture)土壤分类，该流域的土壤主要为砂壤土，易受降雨侵蚀和水力侵蚀的影响。在大规模开展植被恢复以前，该流域年均输沙量为每年12 800 t/km²。其中，大部分输沙来自汛期暴雨引发的降雨侵蚀和黄土滑坡^[29]。与黄土高原众多小流域类似，退耕还林工程在佘家沟流域已经开展了30余年，流域内主要的土地利用类型已由20世纪70年代的裸土和耕地(占比63.20%)转变为如今的草地和林地(占比85.80%)^[19]，是一个典型的退耕还林治理流域。

1.2 数据来源

本文研究使用了包括地形、土壤、降雨-径流事件和土地利用类型等数据，用于建模、率定验证和工况设计工作，如表 1所示。表中，K_s为土壤饱和导水率。土壤水分特征曲线和K_s均为室内测试。

佘家沟流域2010年30 m精度的数字高程模型来源于国家测绘局，对比20世纪60年代黄河水利委员会的实地测绘地图，主要地形变化来自主河道淤地坝修建导致的河道淤积。流域内30 m精度的土地利用数据来源于黄土高原科学数据共享平台，展示了1975—2020年流域的土地利用类型转变^[31]。1959—1969年子洲试验站在佘家沟设置了水文站和雨量站(见图 1a)，开展暴雨-径流-泥沙观测工作，其间共发生31次有纪录的暴雨洪水过程^[31]。考虑观测时间精度和洪水代表性，本文研究从中筛选出12场具有代表性的典型暴雨-洪水事件，用于模型第1阶段的参数率定和模型效能验证，所选场次在洪水涨水-峰值-退水阶段的观测时间精度为5~60 min，在采样频次上保障了对事件的观测精度，同时涵盖了大、中、小型洪水，保障了流域内径流过程的代表性。

2021年8月对土壤含水量进行测定，在佘家沟左侧三号支沟坡面上选取坡顶、坡中和坡脚3个位置埋设TDR (time domain reflectometry)，进行0.1~1.5 m深度的土壤含水量观测(见图 1b)。土壤测定期间，共获得4场降雨-径流事件后的土壤含水量观测数据，用于模型第2阶段的参数率定和模型效能验证工作。

土壤水分特征曲线表征土壤在非饱和条件下渗透能力和土壤含水量的非线性关系，是描述非饱和区土壤水分运动的重要条件^[32]。本文研究于2014年开展了黄土高原小流域土壤采样-测试工作，采用高速离心机法测定了裸土和不同植被覆盖下浅层土壤的水分特征曲线、孔隙度和容重的基础参数。数据源于佘家沟的邻近流域(王茂沟)，二者在地形、气候、土壤性质和退耕还林还草的进程上基本相似，可代表研究区的土壤水分特性^[19]。

在2018—2019年，对土壤K_s进行采样测定，在实验室内采用恒定水头法测定^[20]。采样测定期间，共测得黄土高原10个典型流域185个采样点的1 021个土壤样品K_s。利用SPSS软件进行变异系数、Kolmogorov-Smirnov检验和单因素方差分析，确定土壤在不同深度和植被覆盖下K_s的显著性差异，获得黄河中游主要支流(皇甫川、孤山川、窟野河、佳芦河、无定河和延河等)流域的裸土、草地、灌木和林地在10、20和30 a恢复期后不同深度的K_s，选取佘家沟所在支流的K_s数据用于模拟工况设计。

2 研究方法 2.1 模型结构

InHM模型是一种基于物理的综合水文模型，已成功用于流域-场次尺度的降雨径流模拟^[32]，水文响应-地貌演变模拟^[30]、浅层滑坡预报^{[9, 27, 33]}等研究。降雨型浅层滑坡深度小，破坏面常常平行于坡面，基于此，ISA模型被广泛应用于山坡和流域尺度的稳定性评估和边坡失稳的机理研究中^{[27, 34-35]}。因此，在模拟浅层滑坡方面，InHM基于物理过程的水文响应模拟结果耦合无限边坡模型是可靠的。

2.1.1 InHM模型

本文研究运用基于物理过程的InHM，模拟土壤中的三维变饱和水流运动^{[27, 32-33]}，三维土壤水流运动控制方程表示如下：

$ \boldsymbol{\nabla} \cdot f_{\mathrm{a}} \boldsymbol{q} \pm q_{\mathrm{b}} \pm q_{\mathrm{e}}=f_{\mathrm{v}} \frac{\partial \lambda S_{\mathrm{w}}}{\partial t} . $

(1)

其中：q为土壤水达西流量, m/s；q_b为土壤水特定源汇项，s^-1；q_e为土壤孔隙与地表的水量交换，s^-1；λ为孔隙度；S_w为土壤饱和度；f_a为连续介质中单元体的相关面积分数；f_v为连续介质中单元体的相关体积分数；t为时间，s。土壤水达西流量，表示如下：

$ \boldsymbol{q}=-k_{\mathrm{rw}} \frac{\rho_{\mathrm{w}} g}{\mu_{\mathrm{w}}} \boldsymbol{k} \boldsymbol{\nabla}\left(\psi+z_{1}\right). $

(2)

其中：k_rw为土壤相对渗透系数，ρ_w为水的密度，kg/m³；g为重力加速度，m/s²，μ_w为水的动力黏度，kg/(m·s)；k为土壤固有渗透向量；z_l为高程，m；ψ为压力水头, m。土壤固有渗透向量与土壤饱和导水率K_s的关系表示如下：

$ K_{\mathrm{s}}=|\boldsymbol{k}| \frac{\rho_{\mathrm{w}} g}{\mu_{\mathrm{w}}} . $

(3)

土壤中压力水头ψ和孔隙水压力p, kg/(m·s²)，表示如下：

$ p=\rho_{\mathrm{w}} g \psi. $

(4)

InHM运用控制体积有限元法来求解三维土壤水流运动控制方程式(1)，求解过程中时间步长可随计算量大小进行自动调整，即自适应。结合式(1)—(4)可知，通过InHM求解土壤中水流运动过程，可获得任意时刻、任意土体单元的土壤饱和度和孔隙水压力。

2.1.2 ISA模型

在ISA模型中用安全系数(factor of safety, FS)表征任一土壤单元在任意时刻的坡面稳定性。FS大小与土体抗滑力和滑动力的大小有关。FS＞1代表坡体处于稳定状态，FS＜1代表坡体处于不稳定状态，存在滑坡风险。安全系数可表示如下：

$ \mathrm{FS}=\tau_{2} / \tau_{1}, $

(5)

$ \tau_{2}=c+(\sigma-p) \tan \varphi. $

(6)

其中：τ₂为变饱和土体的土体抗滑力，τ₁为变饱和土体的滑动力；c为土壤的有效黏聚力，kPa；σ为土体正应力，kPa；φ为土体有效内摩擦角, (°)；p为土体任意位置、任意时刻的孔隙水压力，可根据其与压力水头的关系式(4)，由InHM求得。土体正应力可表示如下：

$ \sigma=\cos ^{2} \theta \int_{0}^{z_{2}} \rho g \mathrm{~d} z_{2}. $

(7)

其中：θ为土体单元的坡度；z₂为土体单元距离地表的深度，m；ρ为土壤容重，kg/m³。变饱和土体容重随土壤饱和度的变化而变化，表示如下：

$ \rho=\lambda \cdot S_{\mathrm{w}} \cdot \rho_{\mathrm{w}}+(1-\lambda) \cdot d_{\mathrm{s}} \cdot \rho_{\mathrm{w}} . $

(8)

其中d_s为土壤颗粒比密度，为无量纲。

变饱和土体的滑动力主要来自土体在水分入渗条件下不断增加的重力，表示如下：

$ \tau_{1}=\sin \theta \cdot \cos \theta \cdot \int_{0}^{z_{2}} \rho g \mathrm{~d} z_{2}. $

(9)

整合式(5)—(9)可得ISA中边坡稳定性安全系数的计算公式^[9]，表示如下：

$ \mathrm{FS}=\frac{c+\left(\cos ^{2} \theta \cdot \int_{0}^{z_{2}} \rho g \mathrm{~d} z_{2}-\rho_{\mathrm{w}} g \psi\right) \tan \varphi}{\sin \theta \cdot \cos \theta \cdot \int_{0}^{z_{2}} \rho g \mathrm{~d} z_{2}} . $

(10)

ISA模型假定滑移面平行于地表。野外调查显示近几年黄土高原植被恢复区浅层滑坡的滑移面大多平行于地表^[22-23]，因此ISA模型在本文所关注的植被恢复区是适用的。

2.2 模型设定 2.2.1 初始条件及边界条件

模型使用的计算网格基于前期的研究工作^[29]，采用增加层数的方法将地表DEM转化为三维网格，如图 2所示。此网格水平和垂直方向的离散步长分别由水文响应过程非活跃区的180 m(分水岭)和39.3 m(下边界)加密至水文响应活跃区的5 m(河道)和0.10 m(近地表)。此种离散方法既保证了计算网格在水文响应过程活跃的地方(河道、沟道、近地表土壤)有足够的精度，又能通过在不活跃的地方(分水岭、深层土壤)降低网格精度来节省计算资源的损耗^[32]。

雨季实际土壤水分布的初始条件由饱和排水法确定。通过饱和排水法获得佘家沟流域汛前平均流量(0.003 m³/s)，此流量下的流域土壤含水量分布即为模拟的初始条件，代表了佘家沟流域基于物理且有效的连续压力水头分布^[19]。

覆盖全流域的降雨和流域出口的自由出流构成了模型的地表边界条件。不透水边界为侧面的A-E-J-F，E-D-I-J，C-D-I-H，B-C-H-J和底面F-G-H-I-J。此外，A-B-G-F为流域的地下出流边界。在此边界条件下，壤中流以一定的水力梯度向下游源汇点出流^[29]。由于黄土高原的降雨多为短历时且高强度的降雨，因此在研究中对蒸发量忽略不计^[36]。

2.2.2 模型参数化

模型参数组来源于实测、率定和文献资料，如表 2所示。水文模块中土壤孔隙度、残余含水率、容重来源于实测数据，其纵向变幅不大，因此取平均值，取值均在参考范围之内。模型采用V-G(Van-Genuchten)公式^[32]描述非饱和土壤水分特征曲线，采用高速离心机法测得不同含水率下的土壤吸力变化曲线，进而拟合得到V-G公式中的参数α和β^[37]。Manning系数影响水文模块的地表产流过程，其值由率定得到。坡体稳定性模块中，研究区黄土的土体有效内摩擦角和有效黏聚力由文[9, 38]获得。植被恢复前，模型中各层土壤的饱和导水率由第1阶段的参数率定获得；植被恢复后，各层土壤的饱和导水率由第2阶段的参数率定获得，并与2019年的实测数据作对比。植被恢复前后K_s的率定值的取值范围为2.9×10^-6~50.0×10^-6 m/s，表明其在植被恢复的50多年内发生了明显变化。

2.2.3 模型的率定和验证

基于已有的出口流量过程和土壤含水量观测数据，模型的参数率定工作分为2个阶段进行, 如图 3所示，K_{s, i}指待校准的饱和导水率，n_i指待校准的Manning系数。首先，利用1960—1969年12场降雨-径流事件的出口流量观测数据初步校准控制降雨入渗和地表产流的参数组，保证模型能够准确描述研究区在进行植被恢复工作前的降雨再分配(降雨-产流-入渗)过程。然后在已获得参数组的基础上，利用2021年土壤含水量观测数据进一步校准参数组，保障模型能够准确描述降雨入渗后的土壤水分运动过程。在第1阶段的率定过程中，用6场降雨径流数据进行参数率定，另外6场数据用于模型效能验证。在第2阶段的率定过程中，用2场土壤含水量观测数据进行参数率定，另外2场数据用于模型效能验证。在率定验证过程中，采用NSE(Nash-Sutcliffe efficiency)系数和决定系数R²来评价模拟结果的准确度，相关指标可表示如下：

$ {\rm{NSE}} = \frac{{\left[ {\sum\limits_{i = 1}^m {{{\left( {{O_i} - \bar O} \right)}^2}} - \sum\limits_{i = 1}^m {{{\left( {{S_i} - {O_i}} \right)}^2}} } \right]}}{{\sum\limits_{i = 1}^m {{{\left( {{O_i} - \bar O} \right)}^2}} }}, $

(11)

$ {R^2} = \frac{{{{\left[ {\sum\limits_{i = 1}^m {\left( {{O_i} - \bar O} \right)} \left( {{S_i} - \bar S} \right)} \right]}^2}}}{{\sum\limits_{i = 1}^m {{{\left( {{O_i} - \bar O} \right)}^2}} \sum\limits_{i = 1}^m {{{\left( {{S_i} - \bar S} \right)}^2}} }}. $

(12)

其中：O_i和S_i分别为实测值和模拟值，共有m个数据进行率定验证；$ \bar{O}$和$ \bar{S}$分别为实测数据和模拟数据的平均值。NSE的值为对水文模型模拟结果的评价，取值范围为(-∞，1)。NSE值越接近1，代表模拟效果越好，模型可信度越高；NSE值越接近0，代表模拟结果接近观测值的平均水平，总体结果可信，但过程模拟误差较大；当NSE值远小于0时，代表模型不可信。R²值越接近1，代表模拟值与观测值的拟合程度越好。

第1阶段模型参数率定验证结果显示模型对流域内降雨-径流事件的模拟效果较好，率定期NSE系数为0.763~0.960，R²值为0.81~0.99；验证期NSE系数为0.730~0.890，R²值为0.71~0.91，能够较好地捕捉流域出口洪水的涨水过程、峰值流量和落水过程，证明此阶段获得的参数组能够较好地描述佘家沟流域下垫面对暴雨的再分配过程^{[19, 30, 36]}，如表 3所示。第2阶段模型率定验证结果显示模型对坡面降雨入渗和土壤水分运动过程模拟效果较好，如图 4所示。率定期NSE系数为0.730~0.870，R²值为0.78~0.97；验证期NSE系数为0.480~0.830，R²值为0.59~0.87(见表 3)。图 5a显示模型能够较好地展现雨前土壤水分的纵向分布，可以为降雨-滑坡事件的模拟提供一个较为准确的前期条件；图 5c和5d显示模型能够较准确地模拟由雨后水分入渗导致的浅层土壤水分累积过程，对水分累积的最大深度模拟较好。对于0.5 m深度以下的水分消退过程，模型能较好地模拟其趋势，但对深层土壤含水量的波动模拟效果欠佳(见图 5d)，这可能与土壤深层局部区域的异质性、根系和虫洞造成的优先流路径有关。总体而言，经过2个阶段的参数校准后，模型能够较好地模拟不同的土壤深度、坡面位置在不同降雨时刻的土壤含水量变化，为模型对暴雨-滑坡事件过程中土体饱和度和孔隙水压力变化的模拟提供基础。

通过阶段1和阶段2的参数率定，获得2组参数组。对比2个阶段的参数组发现，与20世纪60年代相比，佘家沟流域当前状态的土壤饱和导水率K_s发生了较为明显的变化，与佘家沟流域实测的现状土壤饱和导水率一致，这是由植被恢复改变土壤浅层物理性质导致的。Manning系数n在2个阶段的模拟中未发生明显变化，这是由于在第2阶段模拟中，校准对象并非地表径流的运动过程。

2.3 工况设计

为了研究不同植被恢复条件下K_s纵向异质性对坡体稳定性的影响，本文选取黄土高原3种典型植被恢复类型(草地、灌木和林地)，与无植被恢复条件(裸土)作对比，并考虑植被恢复期(10、20和30 a)的差异。基于1 021个土壤样品测试结果，选取靠近模拟区域的测试数据，获得不同恢复期下各类植被覆盖区不同深度土壤的K_s平均值，具体土层设置为0~10 cm(第1层)，10~20 cm(第2层)，20~50 cm(第3层)，和50 cm以下(第4层)^[39]。从表 4可以看出，植被恢复增大了各层土壤的饱和导水率，第1层的K_s1由裸土的6.34×10^-6 m/s增大至林地恢复期为30 a的78.44×10^-6 m/s；第4层K_s4也由裸土的1.51×10^-6 m/s增大至林地恢复期为30 a的27.88×10^-6 m/s，且恢复期的差异也影响了K_s的变化。然而，在实际情况中植被恢复、土地利用存在空间分布差异，并非整个流域为同一植被恢复类型。为了简化研究问题，本文在工况设置中不考虑植被空间分布的差异性对流域水文过程和坡体稳定性的影响，而是对任意植被恢复工况，将佘家沟流域饱和导水率按空间均一分布设置。基于山顶植树，山坡种草的传统治理模式，对于流域中某一坡面(长200 m，宽50 m)，其产汇流区域内的植被类型可认为是均一的。此外，研究区典型草被、灌木的根系较短，且以散生、水平根型为主，在本文模拟中不考虑根土复合体的加筋锚固作用，在模型中用额外黏聚力Δc表征不同植被对坡体稳定性的力学影响。

有关山坡水文学的研究表明，前期条件对降雨入渗、坡面产流和土壤水侧向运动过程有不可忽视的影响^[40]。滑坡发生与前期降雨过程也有显著相关性^[41]，连续的降雨过程和有时间间隔跨度的降雨过程对坡体稳定性的影响可能不同^[42]。通过设置不同间隔期的连续降雨，本文探讨不同前期条件下异质性土层的水分运动和累积对降雨的响应过程，以及对坡体稳定性的影响。以历时2 h、单场强度为87.98 mm/h的暴雨为基础，设置间隔期分别为1、3、7 d的3场连续暴雨工况，如表 5所示。降雨工况的设置基于流域相邻雨量站的历史实测数据，降雨强度取20年一遇的特征雨强87.98 mm/h，降雨时长定为2 h，一方面是因为在统计过程中发现70%的最大年降雨强度对应的降雨时长为2 h^[20]，另一方面是确保总降雨量在黄土高原实际降雨范围内^[38]。将土壤参数与降雨类型组合后得到40种工况，各模拟工况独立施加到佘家沟流域上，得到40种模拟结果。

3 模型结果 3.1 场次暴雨条件

在场次暴雨条件下，土壤层K_s异质性显著影响坡体稳定性。在不同恢复期，降雨历时2 h时，裸地、草地、灌木和林地的距地表深度与安全系数FS关系，如图 6所示，在场次暴雨条件下，当所有工况坡体稳定性安全系数FS均出现低于1.00时，均存在滑坡风险。对于无植被恢复的裸土坡面，坡体失稳深度范围为0~0.24 m，且失稳过程随深度向下逐渐发展。在再植草地、灌木和林地恢复条件下，具有失稳风险的深度范围分别为0~0.57 m、0~0.67 m和0~0.70 m，与裸土坡面相比，土地利用程度提高会导致滑动面下移，滑坡厚度增加。这是由于植被根系深度和微生物活动范围不同，导致的K_s影响深度的差异。K_s的纵向异质性导致层间水分积累，从而增加非饱和土孔隙压力，促使滑坡的发生。在模拟工况中，孔隙水压力在降雨结束时刻达到最大值。相比于未恢复的土地情况，土地利用率的提高，使得孔隙水压力在0~0.50 m土层内呈垂直变化。在不同恢复期，降雨历时2 h时，裸地、草地、灌木和林地的地表距离与孔隙水压力关系，如图 7所示。土地利用率越高，表层孔隙水压力的垂直变化的深度愈大，孔隙水压力在深层土壤中变化越平缓。

在不同的植被恢复期条件下，K_s异质性会显著影响坡体滑动面深度、降雨后失稳时刻和降雨持续时间。在强度为87.98 mm/h的场次暴雨条件下，植被恢复期越长，草地和灌木发生滑坡时滑动面深度越大，但林地相反(见图 6)。草地恢复期为10、20和30 a的失稳深度范围分别为0~0.34 m、0~0.45 m和0~0.57 m；灌木恢复期为10、20和30 a的失稳深度范围分别为0~0.42 m、0~0.56 m和0~0.67 m；林地恢复期为10、20和30 a的失稳深度范围分别为0~0.68 m、0~0.67 m和0~0.57 m。植被恢复期的差异影响降雨后失稳时刻，相比于裸地条件，植被恢复期越长，则失稳时刻越滞后。植被恢复期的差异显著影响滑坡的持续失稳时间，相比于裸地条件，再植灌木与草地恢复的时间越久，坡体持续失稳时间越长，然而在再植林地中结果相反。z₂表示距地表深度，如图 8所示，在场次暴雨条件下，草地恢复期为10、20和30 a的失稳时长约为0.55、1.19和1.37 h；灌木恢复期为10、20和30 a的失稳时长约为1.39、1.54和1.67 h；林地恢复期为10、20和30 a的失稳时长约为1.94、0.85和0.36 h。

坡体失稳持续时间与该土层的土壤水饱和度变化趋势一致，如图 8所示。图中，彩色实线为草地、灌木和林地对应的安全系数；彩色虚线为草地、灌木和林地对应的饱和度。饱和度骤增的时段对应安全系数骤降的时段，当土体饱和时，坡体失稳；饱和度下降的时段对应安全系数回升的时段，此时坡体稳定性逐步回升。在安全系数回升的4 h内，3种恢复期的草地工况土壤饱和度由0.96、0.98、1.00降至0.86左右，灌木工况由0.96、0.99、1.00降至0.85左右，而林地工况由0.92、0.96、1.00降至0.86、0.80、0.76。这表明，林地对0.5 m以上的土层具有更强的排水能力，可以避免水分过度累积导致的坡体失稳。此外，对比林地20、30 a与其他7种恢复工况可以发现，林地的恢复期越长，纵向和侧向排水能力越强，进而失稳时间越短，坡体恢复稳定的速度越快。

图 9展示了在不同恢复条件下失稳时刻与失稳深度的关系。从整体分析坡体失稳的开始时刻，草地在10、20和30 a的恢复期，失稳时刻约为降雨开始后0.55、0.73和0.86 h；灌木在10、20和30 a的恢复期，失稳时刻约为降雨开始后0.72、0.91和1.48 h；林地在10、20和30 a的恢复期，失稳时刻约为降雨开始后1.01、1.59和1.89 h。总体来说，植被恢复能够增加土体失稳深度，且失稳往往首先发生在水分累积层。此累积层的深度与植被类型有关，草地、短恢复期的灌木和林地的深度较浅，而长恢复期的灌木和林地的深度较深。然而，对于30 a恢复期的灌木和20、30 a恢复期的林地，表层土壤未出现失稳状态。在基岩较深的黄土高原坡面，中间土层的短暂失稳不易诱发深层滑坡。因此，这3种恢复条件下，坡体最不容易发生浅层滑坡。

相邻2层K_s的比值Δk与安全系数的关系，如图 10所示。当Δk＜1.5时，异质性Δk与FS并没有规律性变化；当Δk＞1.5时，第j+1层土壤的FS总是大于第j层。因此，认为Δk=1.5为土壤K_s纵向异质性的关键节点。Δk接近1，即土壤层间K_s无明显变化时，第j层土壤安全系数大于1的概率较高。基本认为，植被恢复改变的K_s层间异质性促进了水分的滞留，导致坡体失稳的发生。在其他条件相同的情况下，层间差异性越大，越有利于浅层土壤水分积累，越容易导致坡体失稳。

3.2 连续暴雨条件

本文研究了场次暴雨和连续暴雨饱和导水率纵向异质性对坡体稳定性的影响，并探讨了植被恢复对此的影响。研究发现，相比场次暴雨条件下，连续暴雨会增加失稳深度，且植被恢复期越长，灌木及草地的失稳深度越大。在87.98 mm/h强度暴雨条件下，连续型暴雨使草地、灌木和林地的失稳深度、坡体失稳持续时间显著大于单一场次暴雨，而在再植林地中，土壤层的失稳深度会随着植被恢复期的延长而减小。连续暴雨结束时刻与场次暴雨相比，草地、灌木和林地在不同恢复期距地表深度与安全系数的关系，如图 11所示。间隔1 d的连续暴雨条件下，再植草地恢复期为10、20和30 a的失稳深度范围分别为0~0.75 m、0~1.66 m、0~1.78 m；再植灌木恢复期为10、20和30 a的失稳深度范围分别为0~1.63 m、0~1.77 m、0~2.73 m；再植林地恢复期为10、20和30 a的失稳深度范围分别为0~2.35 m、0~1.20 m、0~0.82 m。场次暴雨与连续暴雨条件下草地、灌木和林地在不同恢复期的失稳持续时间，如图 12所示。与场次暴雨相比，连续降雨使坡体失稳持续时间延长，降雨后失稳时刻提前。在再植灌木与草地中，植被恢复期越长，连续降雨后FS越小，持续失稳时间越长，失稳概率愈大；但在再植林地中，植被恢复期越长，持续失稳时间越小，连续降雨对坡体稳定性影响越小。

恢复期20 a的草地、灌木和林地在间隔1、3和7 d的连续降雨情况下，距地表深度与安全系数，如图 13所示。与间隔1 d的连续降雨相比，连续暴雨间隔期越长，表层土壤失稳风险越低，但是暴雨间隔期的水分入渗会降低深层土壤的安全系数；与间隔1 d的连续型暴雨相比，间隔3、7 d的连续暴雨后，1.5 m深度以上的土壤失稳深度和持续失稳时间均减小，FS显著提高。然而，由于水分的渗透较充分，1.5 m深度以下土壤的安全系数明显减小，降雨间隔时间越长，渗透愈充分，失稳风险越大(见图 13b)。综上可知，前期降雨的入渗使得深层的土壤处于一个相对不稳定的状态，下一场降雨发生容易诱发深度更深、方量更大的滑坡，且林地的高渗透能力也会放大这种影响。例如在草地条件下，水分下渗1 d不足以使0~5.0 m厚的坡体处于FS＜1的状态(见图 13a)，而在林地条件下，水分下渗1 d足以使得0.1~5.0 m厚的坡面处于FS＜1的状态(见图 13c)；倘若有新的暴雨发生，林地覆盖的坡面将有5 m厚的坡体发生整体滑坡的风险。

4 讨论 4.1 植被恢复与坡体稳定性

本文研究植被恢复对坡体稳定性的影响，这一问题长久以来受到广泛的关注^{[5, 24-25, 27, 43]}。研究发现，在降雨条件下，不同的植被类型和恢复期对坡体稳定性的影响不同。对于黄土高原坡体，裸土条件下土壤下渗能力很低，水分以超渗产流形式快速从地表排出，对坡体稳定性的影响较小^[20]。植被恢复改善了表层土壤的渗透能力，降低了产流率，但是不同类型和恢复期植被的土壤排水能力存在差异^[23]。草地和灌木恢复条件下土壤水分容易富集在某一土层，造成表层土体滞水超重并失稳，主要受上下层土壤渗透能力的差异和排水能力偏低的综合影响。相比之下，在林地条件下土壤各层排水能力较强，层间渗透差异引起的上层滞水效应不明显，通畅的纵向排水能力保障了浅层坡体的稳定性。

文[7]发现在黄土高原草本和灌木覆盖的区域，浅层滑坡占滑坡总数的83%，远高于乔木和裸土，与本文的研究结果一致。同时，本文研究发现在林地条件下，土壤较强的纵向入渗能力会增大降雨的入渗深度(由1.0 m^[7]增大至5.0 m)，增大下一次降水过程中的滑坡风险和深度，与黄土高原南小河沟的调查结果一致。植被的恢复不仅会影响各层土壤的饱和导水率，还会影响浅层土壤的生物量、孔隙度、容重、非饱和条件下的渗透能力(V-G公式中的参数α和β)和土壤的黏聚力^{[7, 27]}等。同时，植被根系导致的优先流路径也会影响水分的滞流和下渗。因此，不同流域中植被类型和分布对滑坡的响应规律有所不同。

4.2 植被恢复与未来黄土高原水土保持工作

黄土高原进行了退耕还林、水土保持的工作，旨在减少沙漠化和土壤侵蚀，已经持续了数十年并取得了令人瞩目的成果。研究表明，黄河年均输沙量从1987—1999年的8.1亿t减少到2000—2018年的2.5亿t^[13]。然而，野外调查和数值模拟研究表明，黄土高原植被恢复区浅层滑坡的频率和方量有所增加。这些滑坡体会堆积在坡脚或河道中，如果遭遇暴雨或洪水，极有可能被冲走。特别是在极端降雨事件增多^[44]、坝库工程拦沙能力逐年下降^[21]的情况下，这些新的沙源将对黄土高原流域的侵蚀产沙过程产生一定影响。此外，实验研究发现已发生过浅层滑坡的黄土高原坡面，其径流率增加范围是原来的4.0~13.1倍，产沙率增加范围是原来的70~841倍^[23]，对新时期的水土保持和生态恢复工作提出了新的挑战。

5 结论与展望

本文通过基于物理过程的数值模拟，研究了植被恢复对黄土高原降雨诱发型滑坡的影响，主要结论如下：

1) 草地和灌木覆盖的坡体比无植被恢复的坡体更容易失稳。随着植被恢复期的增加，失稳风险和滑坡深度增加，表现为表层土体滞水超重并失稳，失稳最先发生在水分累积层，累积层的深度与土壤饱和导水率纵向异质性有关。

2) 因各层土壤的高饱和导水率不同，林地覆盖的坡体比草地、灌木覆盖的坡体更稳定。在同等场次降雨条件下，林地覆盖的表层土壤不容易失稳，且因其通畅的排水能力，失稳持续时间短，稳定性恢复较快。

3) 不同植被恢复条件下，连续暴雨-入渗造成的高前期含水率显著增大了坡体失稳风险。此外，林地的高渗透能力会放大这种前期条件的影响。将降雨的入渗深度由1.0 m增加至5.0 m，增大下一次降水过程中的滑坡风险和深度。

4) 土壤饱和导水率纵向异质性越大，越有利于土壤层间水分积累，坡体失稳发生的可能性越大。

本文研究采用数值模拟的方法初步探究了植被恢复对土壤饱和导水率纵向差异，以及对场次和连续暴雨下坡体稳定性的影响。未来的研究可以从以下2方面进一步探究植被恢复对浅层滑坡过程的影响：将植被恢复引起的浅层土壤容重、黏聚力、孔隙度、非饱和条件土壤渗透特征参数的变异以及蒸发、截留造成的水分损失纳入研究范围，系统性地探究植被恢复导致的下垫面纵向和侧向异质性对坡体稳定性的影响；系统性地量化黄土高原浅层滑坡对河流输沙量的贡献率，明确浅层滑坡堆积体对山坡和河道水沙连通性的影响。这些工作将对未来黄土高原坡面水土保持工程和退耕还林工作的科学开展有很大帮助。