黄河上游沙漠小流域粗沙岸坡重力侵蚀滑塌过程及型态的数值模拟

王之君, 张 科, 李余杰, 郑 健

(兰州理工大学 能源与动力工程学院, 甘肃 兰州 730050)

沙丘坡沟系统是黄河上游沙漠流域风水复合侵蚀产沙的主要固体物质来源.沙丘岸坡在降雨入渗及坡脚水流侵蚀等因素共同作用下发生重力滑塌,使得沟道水流含沙量剧增,进而演变为高含沙洪水乃至泥流灾害.因此,其侵蚀产沙机制研究可为沙漠流域水土流失综合治理及高含沙洪水灾害防控提供科学依据.相关学者发现:黄河流域坡沟发生重力侵蚀的概率极大,并且重力侵蚀对坡沟产沙过程有显著影响[1],但重力侵蚀由于具有发生突然、滑塌危险性大、诱因多等特点,使得野外定量观测等受到极大限制,成为当前国内外研究的难点.张霞等[2]对流域重力侵蚀破坏概率可靠度进行了分析,同时为建立流域稳定性与重力侵蚀可靠度指标评价体系提供了Rosenblueth法.杨吉山等[3]研究了重力侵蚀发生的特点和影响因素,发现坡面坡度是小型重力侵蚀的主要触发因素.高哲等[4]对运行期淤地坝坝控区域的重力侵蚀来源、诱发因素及预测等进行了分析,发现在黄土丘陵沟壑区淤地坝的库区拦沙来源,既有上游水力侵蚀产沙又有库区自身重力侵蚀.余璐等[5]对影响黄土高原重力侵蚀的地貌等相关因素进行敏感性分析,发现坡度对重力侵蚀总量、滑坡侵蚀量和崩塌侵蚀量的影响较大,且都呈正相关.王光谦等[6-7]、刘清泉等[8]从土壤水动力学角度研究了黄土高原丘陵沟壑区坡面侵蚀机制.然而,上述研究主要集中在黄河中游的黄土高原流域,研究对象多为黄土岸坡,关于黄河上游沙漠流域坡沟重力侵蚀规律与机制的研究,却鲜有报道.笔者研究团队通过前期有关黄河上游“十大孔兑”沙漠流域高含沙洪水侵蚀产沙过程研究[9-11]发现:由沙漠粗沙(粒径0.08～0.25 mm)组成的沙丘岸坡,其重力侵蚀滑塌过程较之黄土岸坡而言,滑塌过程极其迅速,滑塌后的岸坡多为直立、微弯型态,而非黄土岸坡的圆弧型态.笔者认为:上述现象或许与沙漠粗沙颗粒间黏聚力小、降雨入渗系数大等因素有关.因此,本文选取黄河上游“十大孔兑”典型沙漠小流域—苏达尔沟,重点针对入渗条件下其沙漠沟道沙丘岸坡的重力侵蚀滑塌过程和型态开展初步的数值模拟研究,与黄土岸坡进行对照,以期揭示沙漠粗沙岸坡的滑塌机制.研究将为下一步建立符合真实物理图景的沟坡重力侵蚀模型,进而定量计算其侵蚀产沙量奠定理论基础,就维护新时期黄河上游沙漠流域水土生态系统健康稳定而言,亦具有较强的实践指导意义.

1 研究区域与数据来源

1.1 研究区域概况

黄河上游“十大孔兑”沙漠流域位于黄河上游宁蒙河段下段,干流三湖河口～头道拐水文站区间南岸,由穿越库布齐沙漠的十条由南向北贯入黄河的季节性山洪沟组成(图1a,T-T10,孔兑:蒙语山洪沟),自西向东,地表物质组成依次由沙漠向黄土过渡.毛不拉孔兑(T1)中游沙漠沟道区多为流动沙丘,风沙活动尤为剧烈,暴雨山洪诱发的沙丘岸坡重力侵蚀滑塌及高含沙洪水过程最为典型,是一条危害大、治理难的山洪沟.风季,沙丘随风自西向东移动,除部分粗沙堆积在沟道中以外,沙丘整体停靠在沟道岸边,形成左岸沙丘岸坡,雨季,洪水顺沟道倾泻而下,裹挟沟道中堆存的沙漠粗沙的同时,沟道水流不断掏蚀坡脚处,导致坡沟在重力侵蚀与水力侵蚀的复合作用下发生滑塌,滑塌的坡体进一步为洪水补给了沙源,挟沙水流含沙量剧增,进而引发高含沙洪水乃至泥流灾害.苏达尔沟是毛不拉孔兑的一条代表性支沟,因其沟道相对顺直,流域面积相对较小(约60 km2),且笔者研究团队前期已布设了风沙水沙过程观测仪器设备,建立了高含沙洪水综合观测场(图1b),故选择该流域作为研究区域,其岸坡重力侵蚀的物理过程与“十大孔兑”流域一致.

图1 研究区域及高含沙洪水综合观测场仪器设备布置

1.2 观测数据来源

为方便研究,依托前期研究建立的流域高含沙洪水综合观测场,在8 km长的沙漠沟道布设简易断面180个,间距80 m,其中上游丘陵沟壑区100个,中游沙漠沟道区80个(图1b),用以量测洪水前后的断面型态,据此估算风洪产沙比,定量分析泥沙输移和沟道冲淤型态等.2011年～2015年间,共捕捉到6次暴雨洪水事件,其中4次演变成了泥流灾害,同步进行全流域野外调研与岸坡土壤采样,用环刀法及变水头试验渗透试验联合测定土壤容重及渗透系数,用以间接表征降雨条件下的入渗特性.上述工作为本研究提供了丰富的观测数据基础(前期研究详见文献[9-11],此处不再赘述).本文将在野外观测数据基础上,提取研究区域数字高程(DEM)数据,重点开展沙丘岸坡重力侵蚀滑塌过程及型态的数值模拟.

2 数值模拟准备

基于研究区DEM数据、地表物质组成粒径分析数据及土体渗透系数实验数据,运用岩土工程软件GEO5建模、设定材料属性及边界条件、进行有限元模拟,将其结果与ABAQUS软件和野外实测资料进行对比验证.

2.1 数据来源及几何建模

从GoogleEarth中截取“苏达尔沟” 8 km的沙漠沟道区域,精度为30 m×30 m,如图2a所示.用AutoCAD小程序DXF2xyz提取图2a的数字高程点,将其导入Surfer生成“苏达尔沟”流域坡沟高仿真地形图,如图2b所示.在图2b中截取断面No.33处(图1b)一定范围内的沙丘岸坡区域(x坐标为108.855～108.865 m,y坐标为40.22～40.225 m),将所截取的数字高程点导入岩土工程有限元分析软件GEO5中,建立该断面处二维岸坡几何模型,划分网格,如图2c所示.网格边长为0.05 m,生成的网格为6节点三角形,网格节点数为14 988,单元总数为7 851.

图2 原始地形、三维高仿真地形、网格划分

2.2 本构模型及材料参数设置

数值模拟采用GEO5 岩土工程软件中有限元分析模块,其中有限元分析设置为边坡稳定分析.稳定性分析过程中选用理想弹塑性模型作为本构模型[12].模拟过程中理想弹塑性模型采用M-C(Mohr-coulomb)屈服准则,其土体破坏的强度公式为

τ=c+σtanφ

(1)

式中:c为土体颗粒之间黏聚力,kPa;σ为滑移面上土体颗粒所受主应力,kPa;φ为土体的内摩擦角,(°);τ为土体颗粒所受的切应力,kPa.

模拟过程中以有限元强度折减法为主要方法,其黏聚力、内摩擦角折减为

式中:c′为折减后黏聚力,kPa;φ′为折减后内摩擦角,(°);F为岸坡失稳的最小安全系数,即折减系数.

主要原理是用折减后的参数进行岸坡稳定性分析计算.通过反复试算,直到岸坡失稳.数值模拟过程中优化土体容重及渗透系数两个参数.用环刀法测土样的容重,室内变水头渗透实验测渗透系数[13].试验测得粗沙的干容重为1.23 g/cm3,黄土干容重为1.28 g/cm3.

其中变水头渗透系数为

(4)

式中:kt为水温t℃时的试样渗透系数,cm/s;a为变水头管的内径面积,cm2;2.3为ln和lg的变换因数;L试样高度,cm;t1、t2分别为测读水头的起始和终止时间,s;H1、H2分别为起始和终止水头,cm;A为试样的过水面积,cm2.

标准温度下的渗透系数为

(5)

式中:k20为标准水温(20 ℃)时试样的渗透系数,cm/s;ηt为t℃时水的动力粘滞系数,kPa·s;η20为20 ℃时水的动力粘滞系数,kPa·s;ηt/η20为粘滞系数比.变水头渗透实验数据如表1所示.实验测得粗沙渗透系数为3.11×10-4cm/s,黄土为2.64×10-4cm/s.

表1 变水头渗透实验记录

3 结果与验证

3.1 x方向位移变化

模拟粗沙岸坡、黄土岸坡滑塌过程中x方向位移等值线比较分别如图3a、3b所示.

图3 不同岸坡沿x方向位移等值线比较

由图3可知:粗沙岸坡的滑塌沿浅层滑动面移动,且滑塌主要表现为分层、条状滑塌;黄土岸坡的滑塌表现为团状或块状滑塌.分析认为:岸坡土体稳定性不仅与自重有关,还与土体颗粒之间的黏聚力有关,当坡角一定时,岸坡稳定性对黏聚力的变化最为敏感,黏聚力越大,坡沟系统稳定性越强[14].自然界中粗沙土体颗粒之间没有黏聚力,岸坡土体大多以粗颗粒为单元进行滑塌,且部分滑塌对其他区域影响甚微,出现条状、分层滑塌现象.与之相对应的黄土颗粒之间有黏聚力,岸坡土体颗粒间黏聚力作用往往以颗粒团或者块的形式滑塌,且在黏聚力作用下部分滑塌带动其他区域滑塌,从而出现团状或块状滑塌现象,相关学者也得到了相似的结论[15 -16].

3.2 y方向位移变化

图4a、4b分别为粗沙岸坡、黄土岸坡在滑塌过程中y方向位移变化的等值线比较.由图4可知:粗沙、黄土岸坡在重力作用下坡脚处y方向位移较为明显,因此岸坡土体的失稳破坏节点最有可能出现在坡脚处.分析模拟结果认为:在未考虑坡脚水流冲刷作用时,坡脚处已是失稳破坏最危险的区域,岸坡稳定性对其变化敏感[17].实际情况下,水流侧向掏蚀坡脚加速滑塌过程,一旦坡脚处开始滑塌成为滑塌体,岸坡土体便迅速失稳,进而引起岸坡土体整体重力滑塌.因此,后续研究过程中将重点关注坡脚水流侵蚀作用的影响,同时,在侵蚀破坏防护时,亦应着重关注坡脚处.

图4 不同岸坡沿x方向位移等值线比较

3.3 塑性变形

图5a、5b分别为粗沙、黄土岸坡塑性变形的等值线比较.从图5滑塌体分布的范围可知:粗沙岸坡塑性变形主要集中在坡面附近,同图3a一样呈规则的带状分布,沿浅层滑动,且塑性变形区域呈现坡顶与坡脚贯通趋势.这与文献[14]、文献[16]数值模拟得到的结论相同.而与之相对应的黄土岸坡,塑性变形集中在坡体中部一定范围内.粗沙岸坡塑性变形最大值为2.75,黄土岸坡塑性变形最大值为4.40,黄土岸坡的塑性变形值明显高于粗沙,为粗沙的1.6倍.

图5 不同岸坡塑性变形比较

粗沙岸坡与黄土岸坡塑性变形值的差异,主要是黏聚力差异所致.粗沙坡沟土体颗粒之间黏聚力近似为零,土体的抗剪强度主要由内摩擦角决定.而黄土岸坡土体颗粒之间有黏聚力,当岸坡土体发生滑塌破坏时,不仅要克服土体的内摩擦角,还要克服颗粒之间的黏聚力,因此同样时间步长内,相较于粗沙岸坡滑塌过程,黄土岸坡滑塌产生的塑性变形尤为显著.

3.4 模型验证与实测对比

3.4.1普通边坡模型验证

为验证黏聚力对岸坡土体滑塌型式的影响,本研究进一步采用ABAQUS软件建立普通模型进行验证,普通岸坡模型参数如图6所示.两种工况黏聚力的强度折减设置如表2所列.模拟过程同样采用强度折减法,强度折减系数以场的形式赋予ABAQUS,以塑性贯通区域作为岸坡失稳的依据.模拟结果如图7所示.两种工况下岸坡土体出现塑性变形及完全贯通所需时间如表3所列.

表2 两种黏聚力工况下的折减

表3 不同工况下岸坡土体出现变形及失稳所需时间

图6 岸坡几何尺寸(m)

图7 两种工况下岸坡塑性变形

由图7可知:黏聚力的工况不同,岸坡土体塑性贯通区域大小也不相同,黏聚力越大,塑性贯通区域越小,但岸坡土体失稳都是从坡脚附近开始.

由表3可知:不同工况下岸坡土体发生塑性变形的时间也不相同.两种工况下土体出现变形的时间分别为0.32、 0.40 h,完全失稳所需时间分别为1.37、1.44 h.可以看出:黏聚力越大,岸坡土体出现失稳现象的时间越迟,完全失稳所需的时间也更久.这两种工况的黏聚力相差较小,岸坡失稳就表现出明显的差距.因此,黏聚力是造成粗沙与黄土岸坡在滑塌过程和型态方面不同的主要原因.

3.4.2野外实测资料的对比

图8分别为2011年8月23日暴雨洪水前后典型断面(No.33)型态实测数据及野外实测与模拟粗沙岸坡、野外实测与模拟黄土岸坡的滑塌型态图.由图8a可知:No.33典型断面暴雨洪水冲刷后,左岸沙丘岸坡出现明显的滑塌,且相应的滑塌型态呈明显的折线型.

图8 野外实测与模拟所得滑塌型态图

由图8b～8e可以进一步看出粗沙岸坡滑塌表现为直立型,与野外观测结果定性吻合,亦与拓万全等[18]在同一区域的研究结果一致;黄土岸坡滑塌型态为圆弧型,两者形成鲜明对比.模拟中黄土岸坡呈现的上缓下陡的圆弧型,亦与当前众多学者研究结果相吻合[19-21].

4 结论

依托黄河上游“十大孔兑”沙漠黄土过渡带地貌景观实际,基于强度折减理论,用岩土工程有限元模拟软件GEO5模拟粗沙岸坡与黄土岸坡滑塌过程和型态,并用ABAQUS软件验证模拟结果的适用性和可靠性,最后与野外实测资料进行对比,初步得到以下结论:

1) 受土体颗粒之间黏聚力差异影响,粗沙岸坡主要表现为分层、条状滑塌,而黄土岸坡多表现为块状或团状滑塌,黄土岸坡土体塑性变形值明显大于粗沙岸坡.

2) 较之黄土岸坡的弧型滑动,粗沙岸坡为折线型滑动、滑塌型态直立,模拟结果与野外实测资料定性吻合.

3) 坡脚水流侧蚀作用,是诱发重力侵蚀滑塌的另一重要因素,研究为下一步耦合降雨入渗及水流侧蚀等多因素的试验与模拟相结合,建立复合真实物理图景的岸坡重力侵蚀模型,定量确定侵蚀滑塌量等奠定了基础,指明了方向.