重庆缙云山森林火灾干扰边坡的滑坡易发期预测

李 通，王云琦 ，祁子寒，何相昌，李克文

（1. 北京林业大学水土保持学院重庆三峡库区森林生态系统教育部野外科学观测研究站,100083,北京；2. 北京林业大学水土保持学院重庆缙云山三峡库区森林生态系统国家定位观测研究站,100083,北京）

0 引言

森林自然火灾常由夏秋季节地表层干旱和积温造成，过火林地受高温灼烧影响，一方面地表失去植被屏障，土壤性质发生明显变化，另一方面地下植物根系消亡，固土能力丧失，残留根孔，火灾迹地成为浅层滑坡和泥石流等次生灾害的敏感区域。据统计火灾迹地上泥石流和滑坡发生的概率高达40%[1-2]。火后森林根系固土护坡研究对火灾迹地次生滑坡灾害的预测对未来森林管护和防灾具有重要意义。

关于森林火灾后次生地质灾害和土壤侵蚀的研究最早见于国外文献统计，如1934 年的加州洛杉矶火后泥石流摧毁483 座房舍，30 人遇难[3]；2002 年科罗拉多州传教士岭的火灾损毁 29 185 hm2森林，火后泥石流冲毁了桥梁公路及车辆[4]。国内对火灾后浅层边坡稳定的关注也在近十年内逐渐增多，最有名为2014 年甘孜仁额拥沟流域的森林火灾及灾后泥石流事件[4-6]和2020 年凉山森林火灾及泥石流事件，两例事件均造成了大量公共设施和生命财产损失。多项研究表明火灾后的1～2 a 内易发生地表径流型泥石流[1]，而在5～10 a 内则容易发生浅层滑坡[7]。前者常因过火林地表层渗透性降低甚至形成斥水层[8]，由边坡表层超渗产流引起，而后者则多因为植被根系逐渐消亡和固土能力丧失而导致[9]。植物根系的固土强度是动态的，该动态变化引起的土体强度在坡体稳定性中发挥重要作用，根系除了因地上部分移除而腐烂造成活根数量变化外，也包括植物重植和生长造成的根系纤维自身强度的衰减-增强变化，因此根系对土体的净增加强度，即根系固土能力在灾害干扰后会发生衰减-回弹的现象，从而导致了一个滑坡易发期[10]。

当前数值模拟手段已成为分析植被边坡防护效益的主流方法，MAO 等[11]通过“林岛-黏聚力增量模型”模拟了三维边坡形态，根层厚度，林岛位置对稳定系数的影响，指出未来研究中对根系固土能力时间效应的迫切需求；李通等[12]通过有限元理论模拟了外力荷载下根土体的变形和整体稳定过程；稽晓雷等[13-15]通过强度折减法模拟了不同根系分布和生长期限对边坡位移的限制作用。然而当前关于火灾后植物根系死亡-再生过程中根系固土强度的时间动态研究成果极少，暂无可靠的土体强度动态量化方法，少有的研究案例[14-17]只提供了短期的测量数据案例，很难考虑灾后长期根系强度全生命周期动态，这成为了火灾迹地边坡滑动易发期预测研究的最大障碍。SIDLE 等[18-19]对人工砍伐森林的根系进行了长期监测研究，发现森林砍伐后根系强度衰减和恢复的时间动态可分别描述为反曲函数和指数函数，并通过实际监测的曲线拐点值、最大最小强度值回归了根系固土强度模型，该模型在后来的俄勒冈州艾里奥特公园森林两次滑坡事件中得到验证[10]。虽然人工砍伐和火灾有着本质的区别，但同为土壤内部根系固土能力的衰减恢复问题，因此在长时间序列监测数据缺乏的情况下，LEE 等[20]对韩国首尔等地2002—2020 年森林火灾后滑坡研究中发现SIDLE 曲线和滑坡频率具有高度一致性，本研究将SIDLE 模型应用于本区未来浅层滑坡灾害预测。

2022 年8 月21—26 日缙云山林火事件引起了社会广泛关注，其中北碚区虎头村上方的针阔混交林影响最为严重，林缘大面积灌草被毁。本研究对该区土壤根系做了实地调查和室内试验，并采用数值模拟手段分析区域可能的滑坡成因，同时对已有文献中各地区植物损毁后固土强度的恢复和衰减数据进行分析，采用SIDLE 根系动态固土强度预测模型对灾后滑坡易发时间做保守估计，以期为灾后森林养护和类似的次生灾害防治提供参考。

1 理论和方法

1.1 根土试样采集与力学试验

选取重庆缙云山北碚区虎头村上方的针阔混交林为研究采样地（图1 a），于严重过火范围的下坡至坡脚处选取3 个采样点，采样点垂直坡向一字排列，间隔4 m。每采样点移除倒木和表层灰烬，采用（图1b）所示的汽油土钻原地钻孔（钻头装载环刀直径Φ=61.8 mm，高h=100 mm），钻头每入土20 cm 取出环刀，现场破碎并收集原状土样中的根系，用于根系量统计，收集钻孔深度为40、140、230 cm 土层深度的原状土柱，每土层至少3 个100 mm 土柱。同时，在取样点附近对林缘优势灌木大头茶（Gordonia acuminata）和四川山矾（Symplocos setchuensis Brand）根系进行挖掘和采集带回实验室。室内先按照TST-55 变水头试验仪环刀的尺寸裁切为Φ=61.8 mm，h=50 mm 的试样（每土柱可裁2 个试样）进行渗透试验，测定饱和渗透系数，然后将渗透完的饱和土样按照Z-J四联直剪仪环刀尺寸二次裁切为Φ=61.8 mm，h=20 mm的土饼（每试样可裁2 个土饼），筛选无根系残留的素土土饼（较细根系拔除）在0.8 mm/min 速率下进行快剪（图1c）。筛选无损伤根样，统一裁剪为100 mm，两端采用树脂处理[20]后进行单根抗拉测试（图1d）并记录断口直径。根系抗拉强度采用式（1）计算，并按照根系直径和抗拉强度负幂函数曲线拟合：

图1 火灾林地调查及试验Fig.1 Field survey and tests on burnt forestland

式中Tr为根系抗拉强度，MPa；F为根系抗拉力，N；d为根系直径，mm；α和β为拟合参数。

1.2 火灾后根系固土能力预测(SIDLE 曲线模型)

植物根系抗拉强度Tr在土体承受剪切时充分激活，产生抵抗剪切作用的效果，即根系固土能力，当前所采用的根系固土能力测试和量化方法多引自文献[21-22]的系列成果，文献[21]中WU 等计算模型为

式中ΔS为根系固土强度，kPa,As为根-土复合体剪切面面积，m2；k为模型中根系系数，取1.2[22]；n为根系径级数；N为各个径级根系数量；RAR为根系截面积比，1；为第i径级根系的平均抗拉强度，MPa；N为 第i径级根系的总数。根系直径分级为：1～2、＞2～3、＞3～5、＞5～10 mm。

将根系固土强度量化为独立于土体正应力的形式纳入Mohr-Coulomb 框架[23]：

式中c′ 为有效黏聚力，kPa；φ′ 为有效内摩擦角，（°）；(ua-uw)为基质吸力，kPa；χ为基质吸力系数，以饱和度取值；σ为正压力，kPa；ua为孔隙气压力，kPa；uw为孔隙水压力，kPa。

火灾发生后，一方面过火范围内地被物被大量烧毁，残留根系因腐烂而减少，强度逐渐丧失，另一方面，部分残存树木和次生灌木逐渐萌蘖，近表层根系强度又逐渐恢复。SIDLE 等[18-19]基于加利福尼亚州某针阔混交林砍伐25 a 后根系强度变化，并结合O’LOUGHLIN 等[24-26]的试验数据，提出了根系恢复系数模型R和衰减系数模型D：

式中R和D分别为根系固土能力恢复系数和衰减系数，分别为根系生长和死亡条件下当前固土能力与最大固土能力的比值，无量纲；t为森林砍伐后的时间，a；d1、d2、d3、d4、d5和d6为拟合参数，通过监测的根系强度-时间数据拟合得到。

因此，森林砍伐后一段时间t的净根系固土能力Δτ(t)，即SIDLE 模型计算式为

式中Δ τ(t)为砍伐时间t后的根增抗剪强度，kPa；Δcmax为最大抗剪强度，kPa；取值式（2）量化的ΔS。

当参数取值d1=-0.05，d2=0.95，d3=20.00，d4=0.25，d5=0.50，d6=0.75[10]时，复现文献[10]中的SIDLE曲线如图2。

图2 SIDLE 曲线模型-森林砍伐后根系固土强度Fig.2 The root reinforcement estimation model by SIDLE after forest harvesting

SIDLE 限定了R的变化区间为（0，1），建议R拟合的特征数据为曲线拐点（ti,Δτi）和2 倍时间拐点（2ti,Δτ），即R(t)的二阶导在ti处为0，因此有以下条件：

以下将根据这3 个假定条件探讨参数之间的关系，将式（5）代入式（8）有：

将式（11）进一步化简：

同样，将式（5）代入式（9）有：

式（13）进一步化简：

至此，d2、d3、d4均可用d1来表示，可得仅包含ti和d1的R表达式：

按照以上参数关系，在区域长期连续性监测数据匮乏的情况下，结合文献中已知的乔灌木数据可推算植物固土强度恢复的拐点时间分布，从而据拐点反推其他参数。SIDLE 曲线中“滑坡易发期”的确定依据为固土强度曲线谷值和实际灾害在时间上的高度一致性，但并无具体强度或稳定系数的数值界定。本研究中将参考多个已知植物的拐点时间来确定快速恢复和慢速恢复2 种模式下的曲线谷值，并以2 个谷值间隔作为滑坡易发期来计算稳定性。

1.3 边坡模型及数值模拟试验

参考调查点边坡和土壤条件构建几何模型，坡度取30°，坡高取15 m，坡肩和坡脚延伸10 m 以减小边界影响。边界条件：力学方面，底部为全位移约束，四围为法向位移约束，表面为自由面；渗流方面，上表面为降雨边界，不考虑陡坡面积水，以降雨流量80 mm/d 输入，持续4 d，模拟降雨小于土体入渗能力的完全入渗情景。前后两面以水头控制地下水位于0 m，模拟初始孔隙水压力（pore water pressure,PWP）随高程的线性分布；两侧面为对称边界，法向流量为0。为做位移数据分析，在边坡中段1 m 深处设置监测点，坐标（14,5,7），详见图3。

图3 边坡模型边界条件及网格剖分示意Fig.3 Boundary conditions and mesh physics generalization of slope model

参考钻孔的土层深度数据，将边坡土体以1.5 m 深度分为0～1.5 m 和＞1.5 m 两层，设定土体力学和渗透性参数，由于土体的弹性模量和泊松比仅会影响弹性位移，对稳定系数的影响极小，采用经验赋值，含根土的泊松比和弹性模量分别略小于和大于无根土[27]，其他实测参数赋值见表1。考虑灌木林情境，个别粗壮根系的锚固效应较弱，根系固土能力统一概化至附加黏聚力，并考虑其在垂直方向的变异，某时间节点Δτ的赋值在软件中以插值函数输入，即Δτ=int(Z)，每隔1 a 赋值一次，重复计算至40 a。

表1 边坡模型力学水力参数Table 1 Mechanical and hydraulic parameters of slope model

计算思路和物理场设置分3 步进行（图4）：1）在渗流和力学边界控制的初始条件下进行稳态分析，计算初始孔隙水压力场和应力场，在弹性赋值和塑性赋值的条件下逐步消除应变和位移[28]；2）施加降雨边界，进行降雨过程中孔隙水压力瞬态变化条件下边坡有效应力场和位移场的求解；3）对第2 步中各个时间节点进行强度折减[29]，计算边坡整体稳定性。

图4 COMSOL 中的计算思路及模拟流程Fig.4 Calculation procedure and simulation step diagram in COMSOL

式中cm为临界状态黏聚力，kPa；φm为临界状态内摩擦角，（°）；FOS为折减系数，即边坡的整体稳定系数。

降雨入渗时表层孔隙水压力梯度较大，对上层土体加密处理，自上而下渐疏，最小尺寸0.005 m，最大尺寸1 m，共7 2 887 个单元（图3），以保证收敛性和结果精度。其中经步骤1）导入步骤2）边坡初始位移值7.37×10-5mm，远在毫米级以下，说明用于正式分析的初始应力和初始孔隙水压力在给定边界下平衡效果良好，确保了模拟计算的可靠性。

2 结果与分析

2.1 SIDLE 模型曲线拟合参数分析及取值

林火毕竟为低概率事件，区域自1961—2021 年都无火灾事件，期间也无根系的连续性监测数据，不可能实现模型6 个参数的严格拟合。对极少文献资料中长时间序列数据进行了自定义非线性拟合，拟合效果见图5～图6。其中文献[6]为四川凉山苏格兰松（Pinus sylvestris）固土强度衰减数据；文献[18-19]及文献[6]分别为SIDLE 对各地6 种植物和四川凉山苏格兰松的固土强度衰减数据；文献[30]为原作者模拟的9 种乔灌木固土强度增长数据，本研究将其换算为0～1 之间的R值。

图5 9 种植物根系固土能力恢复系数拟合曲线Fig.5 Fitted root reinforcement recovery coefficient curves of 9 plant species

图6 7 种植物根系固土能力衰减系数拟合曲线Fig.6 Fitted root reinforcement decay coefficient curves of 7 plant species

由图5 可知，9 种乔、灌植物的拟合参数差别较大，R2均在0.97 以上，SIDLE 提出的R曲线模型可较好地反映树木移除后的根系固土强度动态，但拐点时间ti不易从图中直接分辨。图6 中SIDLE 统计的根系强度衰减曲线主要为针叶树种，衰减速率差别较大，完全消亡的时间在5～25 a 不等，其中四川凉山研究的苏格兰松（图6f），研究地距本区仅700 km，同属亚尔热带季风气候，本区域森林生境最为接近，因此本研究直接采用了其衰减曲线参数d5=-0.19，d6=1.09。基于图5 中的拟合参数，通过式（15）得到这些植物固土能力恢复的理论拐点时间ti，其分布见图7，可知灌木树种的恢复时间拐点较乔木更靠前，超过该时间后进入快速恢复阶段。考虑本研究植物对象为灌木根系，仅以5 种灌木（图7a）的拐点时间ti范围（4.9～9.7 a）和恢复系数R范围（0.51～0.62）为参考，来确定一般灌木的恢复模式，并应用到本地区。分别取“最长拐点时间+最小恢复系数组合（ti=9.7 a，R=0.51）”和“最短拐点时间+最大恢复系数组合（ti=4.9 a，R=0.62）”作为慢速恢复和快速恢复模式，从而计算其他参数（式（12）～（15））。由此确定本研究区域的SIDLE 曲线及参数见图7b。

图7 由时间拐点估计的根系固土强度曲线Fig.7 Estimated root reinforcement curve according to inflection time points

2.2 林地土体抗剪和渗透特性

图8 中为3 个钻孔至少4 个成功渗透试样的饱和渗透系数（图8a）和抗剪强度参数（图8b）的均值和极值，可知土体渗透系数Ks随深度增大呈降低趋势，表层相对松散的土体测试偏差明显高于深层紧密结构的土体，并且在1.4 m 以下有明显的数据降低，上下层的均值分别5.2×10-6和1.5×10-6m/s，相差3.47 倍。与饱和渗透系数分布相似，表层松散土体具有较大的黏聚力，深层紧密土体具有较大的内摩擦角和较小黏聚力。黏聚力和内摩擦角也在1.4 m 附近存在数据突变，以约1.4 m 为分界，上下层土体黏聚力总均值分别为5.5 和1.8 kPa，上下层内摩擦角总均值分别为 18°和26 °，上下层黏聚力和摩擦角的总均值相差分别达到3.06 倍和0.69 倍。

图8 土壤剖面上饱和渗透系数和抗剪强度指标随深度变化Fig.8 Change of saturated conductivity and shear strength parameters with soil depth

2.3 林地根系分布特征与固土强度

图9 为研究点采集单根抗拉强度随直径变化曲线，可知2 种优势植物根系的抗拉强度均随直径呈负相关关系，细根抗拉强度明显高于粗根，大头茶和山矾的拟合方程分别为Tr=25.65d-0.47（R2=0.712）和Tr=29.18d-0.53，（R2=0.638）。2 种植物总体强度-直径曲线相差不大，考虑2 种植物在土壤中的生长交叉，采样时颜色均为灰褐色不易区分，为方便进一步的固土能力估计，将2 种植物的数据混合后拟合如实线所示，拟合参数为Tr=26.12d-0.49（R2=0.699）。

图9 根系抗拉强度和直径拟合曲线Fig.9 The fitted curves of root tensile strength versus diameters

钻孔环刀内的根系数量随土层深度分布（每个深度5 环刀根系数量总和）如图10a 所示，根系数量整体上随土层深度增大而减小，在第2 层具有最大根系数量，1 mm 细根系主要分布与0～0.4 m 土层，＞5～10 mm 较粗根系分布于0.4 m，0.8 m 深度以下几乎未统计到根系，说明调查点的灌木扎根较浅。结合图9 根系强度参数，采用式(2)对土层的根系固土能力计算结果进整理如图10b，可知根系对土体强度的贡献主要在表层砂质壤土层1 m深度范围，与根系数量分布趋势接近，根系截面积比与附加黏聚力基本表现出正相关的特点，说明根系固土能力由含根量主导，根系固土能力在0.4 m 为最大，可达到18 kPa，1 m 以下根系对土体强度无贡献。

图10 土壤剖面上根系数量、截面积比和根系固土能力Fig.10 The root number,root area ratio and root reinforcement along soil profile

2.4 降雨过程中边坡渗流与变形过程

图11 为降雨1～4 d 过程中边坡渗流场变化，可知渗流场的变化分为3 个阶段：1）表层土体湿润，随深度线性分布的初始PWP 在边坡表层有明显升高，湿润锋快速向下推移，深层土PWP 尚无变化（图11a～11b）；2）湿润锋到达土层结构面，雨水聚集形成局部饱和带(白色0 kPa 等值线)，并伴随着顺坡渗流现象（图11c）；3）饱和带水位线逐渐升高向上坡位蔓延，整个结构面出现正PWP（图11d）)，最大约10 kPa，坡脚积水逐渐透过结构面向下层土体渗流，但地下水位仍无升高现象。

图11 降雨过程4 个时间节点边坡孔隙水压力变化云图Fig.11 Temporal pore water pressure (PWP) field of slope at 4 time-nodes during precipitation

图12 为降雨1.0～4.0 d 过程中边坡位移场云图及位移矢量标注，降雨初期，最大位移出现在初始基质吸力最高的坡肩处(图12a)，表现为轻微膨胀，随着降雨进行，最大位移带逐渐向下坡位置移动(图12b)，当结构面出现积水后，下坡位的水平位移分量迅速增大(图12c)，边坡逐渐呈现出沿结构面的平直型滑动 (图12d)。图13 为对1.0～4.0 d 的4 个时间节点的强度折减结果，计算暂不考虑根系固土作用，降雨初期，边坡具有较高的稳定系数FOS=1.91，塑性贯通带位于边坡深层3.9 m，呈现圆弧状（图13a），随着降雨持续，深层的塑性贯通带逐渐变浅（图13b），FOS也逐渐降低1.72，当结构面出现积水后，塑性贯通带开始转移到结构面处（图13c～13d），FOS迅速降低至1.13，濒临失稳（图13d），说明边坡深层土体强度具有较高的安全储备，而表层土体受降雨影响，有效应力降低幅度远大于深层，最初的深层弱滑动趋势转化为浅表层的强滑动趋势。考虑根系固土作用时，边坡特征点位移及FOS与裸坡对比如图14，由图14a 可知坡体中下部位在降雨前期仅发生较小位移，结合图12 可知该位移为表层土体湿润后基质吸力快速降低发生轻微膨胀（0.25 mm），由于土壤容重增大该变形在数小时内恢复为0，在18 h 开始产生阶跃状蠕动变形，至3.5 d 积水带逐渐形成，出现正PWP，进入浅层土体的加速形阶段。有根系加固的植被边坡和无根裸坡在降雨前期位移接近，进入加速阶段后，根系对位移的限制作用逐渐凸显，与裸坡相比在降雨末期根系对总位移的减小达到30%。图14b为植被边坡和裸坡的FOS对比，二者具有相同的时程变化，均在约3.0 d后开始骤降，整个降雨过程中根系对FOS的贡献率由最初6.8%变化至19.8%。

图12 降雨过程中4 个时间节点边坡总位移变化Fig.12 The temporal development of total displacement at 4 time-nodes during precipitation

图13 降雨过程中4 个时间节点边坡折减的等效塑性应变发展Fig.13 The temporal development of equivalent plastic strain at 4 time-nodes during precipitation by strength reduction

图14 降雨过程中植被边坡与无根裸坡的位移和稳定系数对比Fig.14 Comparison of total displacement and FOS between vegetated slope and bar slope during precipitation

2.5 森林火灾后滑坡易发期和恢复预测

火灾后40 a 内自然状态和极端降雨状态下边坡的稳定性结果如图15 所示，快速恢复和慢速恢复模式对应的稳定性谷值分别为9 和6 a，火灾后边坡在自然状态下Fos变化不大，6～9 a 分别从初始2.05 降低至1.9 和2.0，降低幅度分别为7.3%和2.4%，依然处于稳定水平。这是因为森林土体在较干燥状态下潜在破坏层位于远离植被发育层以下的深处（图13a），表层根系的存在与否，仅影响坡脚剪出口和坡顶剪入口局部强度，不会对深层造成较大影响；而在极端降雨后，边坡潜在滑动面位于浅层（图13c～13d），快速和慢速恢复模式下FOS随灾后年限迅速降低，在第6～9 a 从火灾前的1.41 降低至1.32～1.22，降低幅度分别为6.40%和13.48%，此为滑坡易发期，随后随着根系继续恢复，边坡稳定性逐渐提高，至18～24 a 恢复到火灾前的水平。即使在慢速恢复模式下，火灾后第9 a 边坡遭遇极端降雨下FOS为1.22，依然大于裸坡状态（FOS=1.13）（图13d、图14b），根系发挥了关键的强度补偿作用。

图15 火灾后边坡稳定性预测Fig.15 Slope stability prediction after forest fire

3 讨论

本研究对缙云山过火林地的土壤及根系现状进行了初步调查和室内试验，基于流固耦合理论建立了数值模型，采用SIDLE 曲线对森林火灾后40 a 内边坡在自然状态和极端降雨后稳定性动态做了评估。自然状态下，边坡处于较稳定状态，潜在滑动层位于深层的砂质泥岩层，这得益于下层土体较高的摩擦角和较低的渗透系数；极端降雨入渗使表层砂质壤土层基质吸力丧失，在砂质壤土和泥岩层交界面上形成暂态饱和带[31-32]，与前人的双层土坡结构模拟结果类似。所得降雨中位移的发展符合滑坡演变规律[33]。本研究明确了类似本区域地质条件的边坡失稳机制为：土层结构面上基质吸力转化为正孔隙水压力从而降低了土体的有效应力，这不同于深厚黄土基质吸力降低或者滨水岸坡的坡脚水位升高导致的边坡失稳。

降雨中FOS先缓慢降低后快速降低的结果揭示了植被边坡最危险滑动面由深层向浅层演变的失效过程。黄月华等[34]研究也发现降雨中植被覆盖边坡的FOS随时间变化存在临界时间，超过该时间后FOS开始骤降，且该临界时间与上覆层渗透性和顺层渗流现象密切相关，佐证了本研究所揭示的缙云山滑坡形成机制。

火灾干扰对土体和根系性质的影响机制复杂，力学方面，除了土体渗透性和抗剪强度的降低外，还有根系自身强度的衰减[16-19]，考虑这些效应在短期的根系更新和土体沉降后得到恢复[17]，因此在在长期预测问题中仅考虑根系数量变化，图10b 的结果也反映了根系固土能力主要取决于根系数量；水文方面，有研究者强调了植物消亡后因蒸散发减少而导致的基质吸力丧失、根系死亡腐烂造成的优先渗流路径等水文效应对土坡稳定产生不利影响[35-37]，本研究中未考虑该效应可能低估渗流饱和的时间，因此所估计的单场降雨失稳时间节点4 d 略迟于实际，而临界状态结构面上的有效应力主要受控土体饱和后的孔隙水压力大小，因此所计算的FOS大小不受该效应影响。

对火灾后数十年内的稳定性变化的模拟依赖于精准的根系监测数据，国内对火灾和边坡稳定的专项研究极少，如雷鸣宇等[6]通过连续调查数据，采用无限斜坡模型分析了四川省凉山州木里县的火后9 a 内的浅层滑坡稳定性动态，并给出了根系数量和稳定性单调递减的指数型预测函数，预测的滑坡期约为火后8 a；LEE 等[20]近期研究中也发现韩国首尔等地的森林火灾后滑坡在6～10 a 频率最高，且恢复期大于20 a，本研究的预测结果与之接近。与文献[10]的研究相比，SIDLE 模型考虑了根系恢复过程，更严谨。有限元流固耦合的方法较无限斜坡模型也具有进步性，所模拟的火后稳定性动态反映了植被死亡-恢复过程中FOS快速降低后缓慢恢复的趋势。缙云山1989 年6 月有特大风灾历史，当时包括本研究点在内的60 hm2林木连根毁坏，11 a 后形成垂直分化结构，16 a 后演替为如今的常绿阔叶林-灌木林结构[38]。这与本研究估计的18～24 a 完全恢复边坡稳定较接近，说明采用多植物恢复拐点时间估计的SIDLE 模型参数基本可信。

该研究仅分析了30°边坡情境，显然区域内更陡边坡或土体结构强度更差的边坡未来则面临更大的滑坡风险。本例为对未然事件的理论预测，无论是上文中SIDLE 曲线参数的取值策略、或是SIDLE 的曲线参数取值[17]，还是LEE 等[20]基于历史数据对SIDLE 曲线的验证研究，都存在多植物数据混合借鉴的情况，这是监测数据匮乏条件下迫切应用和科学严谨性的权衡和妥协，可见，未来在森林生态研究中还应加强对地下植物根系连续监测和调查工作。本研究所给出的结果可为近期火后森林边坡防护提供初步借鉴，有待未来数十年连续调查研究的检验。

4 结论

1）缙云山林缘边坡为约1.5 m 厚的砂质壤土-砂质泥岩双层土坡结构，上下层渗透性差异悬殊，上层饱和渗透系数为5.2×10-6m/s，是下层的3.47 倍，上层黏聚力均值为5 kPa，摩擦角为18°，分别为下层的3.06 和0.69 倍。

2）区域的林缘边坡在正常情况下处于稳定状态，潜在滑动面位于深层（稳定系数FOS＞2）,极端降雨（80 mm/d×4 d）后，在土层结构面上容易出现带状积水带，积水孔隙水压力约10 kPa，局部的正孔隙水压力减小了有效应力导致浅层滑动。极端降雨后裸坡临界稳定系数FOS为1.13，考虑当前根系固土的稳定系数FOS为1.41，当前根系护坡贡献为19.8%。

3）火灾发生后数年内，根系消亡对正常情况下深层边坡稳定性无明显影响，而极端降雨后浅层边坡稳定性随火灾后的年动态呈现快速降低后缓慢恢复的曲线趋势，估计在火灾后6～9 a 间出现滑坡易发期，在火灾后18～24 a 才可完全恢复火灾前的稳定状态，该结果可作为火灾后治理参考。