0.   引言

中国西南岩溶地区土壤大量流失及其所引发的石漠化问题严重制约着区域经济社会发展^[1-2]，岩溶地区特殊的地表-地下二元水文结构导致了该地区不仅有土壤地表流失，同时存在土壤地下漏失^[3]。岩溶地区的土壤地下漏失是指土壤以滑动、沉积、伴随水力侵蚀等方式流失^[4]。目前，岩溶地区土壤漏失的研究方法主要有野外径流小区监测^[5]、人工模拟降雨^[6]、USLE/RUSLE模型^[7]、核素示踪^[8]以及探地雷达技术^[9]等。这些传统的方法中在各自情况内具有一定优势，但也有很多的局限性，比如较难应用于岩溶地区土壤地下漏失过程的观测与分析。岩溶坡面土壤的地下漏失过程是一种典型的非连续介质问题，离散单元法（DEM）可以很好地捕捉细观颗粒的运移特征，在模拟非连续介质中具有较大的优势。目前离散单元法已广泛应用于边坡失稳、碎屑流运移、土颗粒潜蚀等问题的研究^[10-11]。其中，CFD-DEM耦合计算模拟被逐渐应用于涉及水土相互作用的颗粒运移问题的研究。离散元法作为一种研究非连续介质的方法，通过研究颗粒的细观力学特性可以较好地揭示土壤地下漏失机理，但是目前将该方法用于研究岩溶地区流固耦合作用下坡面土壤颗粒地下漏失的成果比较少，有必要对该问题进行深入研究。

本文通过室内模拟试验与CFD-DEM数值模拟相结合的方式，探究石漠化地区土壤地下漏失过程规律及细观机理。系统揭示岩溶高原峡谷地区二元水文地质结构下土壤的侵蚀特征、驱动及漏失过程机制，为岩溶石漠化地区水土流失特别是土壤漏失的精准防治提供技术支撑。

1.   研究方法

1.1   岩溶坡面土壤地下漏失物理模拟试验

利用自制的表层岩溶裂隙带土壤地表流失和地下漏失模拟装置进行室内模拟试验，试验具体装置、坡面布设与试验方法见已发表成果^[12]，试验条件设置如下：降雨强度为60 mm/h，坡面坡度为20°，岩溶裂隙管道宽度为20 mm、裂隙走向为0°。

1.2   基于CFD-DEM的土壤地下漏失离散元数值模拟

（1）CFD-DEM流固耦合模型建立

为分析典型岩溶含裂隙坡面土壤地表流失、地下漏失特征，需要使建立的数值模型中颗粒的细观参数表现的宏观力学特性与实际颗粒的力学特性一致，以此提高数值模拟的可靠性，本文选取操作简便的直剪快剪试验获取土的强度指标。固相模型的生成参照室内模型试验实际情况设置，模型箱设置为0.5 m×0.25 m×0.5 m（长×宽×高），固体颗粒及墙体参数的选择按照PFC^3D直剪试验数值模拟标定的结果选取，粒径范围为1~3 mm，服从随机分布特征。土壤颗粒与墙体间的接触为线性模型，固相颗粒间接触为线型接触黏结模型。模型采用砂雨法成样，以更好地模拟边坡的地应力特征。在离散元软件中直接制作含裂隙墙体较为复杂，本文采用Rhino 5软件制作含裂隙坡面模型，再导入PFC软件中模拟墙体。

在PFC^3D中运用Python模块实现流固耦合，渗流条件设置方法为在边界设定压力，产生压力梯度，求得流场速度，最终实现CFD-DEM耦合作用，通过改变入水口和出水口的压力值实现坡面流体流速的变化，进一步实现降雨强度大小的调节，本研究主要考虑降雨条件下坡面土壤流失过程，对不同降雨强度的影响不做具体工况分析。入口压力值恒为400 Pa，地下流出面压力值设置为200 Pa，地表出水口压力值恒为0 Pa，本研究构建的岩溶含裂隙边坡土壤侵蚀流固耦合模型如图 1所示。

图  1  建立的流固耦合模型

Figure  1.  Established fluid-solid coupling model

下载: 全尺寸图片 幻灯片

（2）模型验证

对不同时刻数值试验和物理试验地下产流量及地下产沙量结果进行对比，由于本模型是通过设置恒定出口压力的方式所确定的模型边界条件，受到模型渗透系数、渗透率变化的影响地表出水口流速一直在变化，故地下产流量的统计方法为：监测每一时刻地下出水口流速，将该值乘以模型地表出水口面积得到出水口流量模拟值。将物理和数值模拟试验的地下产流量、地下产沙量及试验时间T以最大值为基准进行归一化处理，绘制对比曲线如图 2（a），（b）所示。通过综合对比分析物理试验与数值模型的地下产流产沙曲线可以看出，采用本文所创建的CFD-DEM耦合数值模型可以对岩溶坡面水土漏失问题进行有效模拟。

图  2  地下产流产沙量物理试验与数值试验的对比

Figure  2.  Comparison between physical and numerical experiments on underground runoff and sediment yield

下载: 全尺寸图片 幻灯片

2.   土壤漏失过程分析

2.1   坡面不同层位土壤颗粒运移过程

降雨是土壤漏失的主要诱导因素，物理试验中难以监测流体和颗粒间的相互作用，本文所建立的CFD -DEM数值模型为分析土壤漏失过程中流体与颗粒的相互作用提供了便利。模型中颗粒的运动特征受流体作用的控制，因此将不同时刻下流体速度云图与不同位置颗粒速度进行对比（图 3），其中T为模拟总时间。由图 3（a），（b）可知，在降雨初期（t < 0.4T），不同层位的颗粒速度均表现出从坡顶到坡底逐渐增大的趋势，其中当t由0.2T增加到0.4T，地下漏失口的颗粒速度明显增大；此外当t=0.2T时地表流失出口处不同层位的土壤颗粒速度差异明显（表层 > 中层 > 底层），而当t = 0.4T时，流失口处的颗粒速度趋于一致。由图 3（b）~（d）可知，在降雨的中后期（t > 0.4T），距坡顶不同坡面位置的土壤颗粒速度差异较小。t = 0.6T~T时，坡面颗粒速度不断增大。以上结果表明降雨历程中不同层位不同坡面位置的土壤运动规律为：降雨初期漏失口附近的颗粒运动加快，同时坡脚的土壤也发生流失；降雨的中后期坡面土壤流失速度不断加快且不同层位的土壤流失速度接近。

图  3  不同坡面位置土壤颗粒速度变化

Figure  3.  Changes in soil particle velocity at different slope positions

下载: 全尺寸图片 幻灯片

2.2   坡面侵蚀形态变化过程

将t = 0.5T时的物理试验结果（坡面土壤侵蚀形态）和CFD-DEM模拟结果（坡面颗粒位移云图）进行对比得到图 4的坡面侵蚀形态对比，可见物理试验的坡面侵蚀形态与模拟结果较为吻合，均表现为在降雨中期流失口的土壤由角部先发生侵蚀，此时漏失口的土壤颗粒也有了较为明显的位移；在降雨后期，坡面前端的土壤颗粒发生明显的整体滑移流失，此时漏失口附近有明显位移的土壤颗粒范围扩大。即岩溶坡面土壤流失规律表现出由流失口和漏失口率先侵蚀发展到整体滑移的过程。

图  4  试验坡面形态与模拟坡面位移云图对比

Figure  4.  Comparison of experimental slope morphology and simulated slope displacement cloud map

下载: 全尺寸图片 幻灯片

2.3   颗粒-坡面相互作用变化过程

颗粒流失过程中颗粒与坡面的相互作用是土壤侵蚀过程中土壤流失难易程度及是否会发生整体滑移的重要依据。其中，坡面底层颗粒与坡面接触的摩擦力是主要的影响颗粒滑移进程的重要因素。图 5为与坡面接触的颗粒在Y方向（即顺坡方向）与坡面间的接触及摩擦力云图，由图可见，当降雨开始后，坡面颗粒在流体力的作用下接触数有明显的降低，约为初始状态的37%；当t > 0.3T之后，颗粒的接触数小幅增大，对比图 5（b），（c），（d）可以看出，裂隙口下缘的颗粒与坡面接触数逐渐增多，主要原因是随着降雨历程，颗粒逐渐在裂隙口堆积密实，表现出裂隙口下缘的颗粒与坡面接触数增多的趋势。由t=0.9T时颗粒与坡面Y方向接触力分布情况（图 6）可以看出来，在降雨后期大部分颗粒与坡面的接触力依然较小，超过70%颗粒与坡面接触力小于0.1 N，并呈现逐级递减的趋势，表明在降雨后期坡面与颗粒之间的摩擦力较弱，流体拖曳作用控制着坡面土壤流失的发生和发展。

图  5  颗粒与坡面Y方向接触力云图

Figure  5.  Cloud map of contact force between particles and slope in Y-direction

下载: 全尺寸图片 幻灯片

图  6  t=0.9T时颗粒与坡面Y方向接触力分布

Figure  6.  Distribution of contact force between particles and slope in the Y-direction at t=0.9T

下载: 全尺寸图片 幻灯片

3.   结论

通过进行降雨驱动下的岩溶坡面土壤地下漏失物理模拟试验和CFD-DEM数值模拟，通过综合对比验证了所构建的CFD-DEM流固耦合模型的适用性和准确性，模拟得到以下结论：

（1）在降雨初期，坡面不同层位的土壤颗粒速度均表现出从坡顶到坡底逐渐增大的趋势，土壤颗粒速度从剖面上表现为表层 > 中层 > 底层；在降雨的中期和后期，坡面土壤颗粒速度不断加快且不同层位的土壤颗粒速度较为接近。

（2）坡面侵蚀形态变化过程表现为，在降雨初期，地表流失出口和地下漏失口的土壤颗粒先发生侵蚀，随后坡面部分土壤颗粒开始发生侵蚀并形成细沟，降雨中期裂隙被土壤颗粒充满并伴随着坡面侵蚀加剧，降雨后期坡面前端的土壤颗粒发生明显的滑移。

（3）当降雨启动后，坡面颗粒在流体力的作用下与坡面的接触数有显著的降低，随着降雨的继续进行，颗粒与坡面的接触数小幅增大，颗粒在裂隙口逐渐堆积，裂隙口下缘的颗粒与坡面接触数逐渐增大；此外，在降雨过程中，坡面上的颗粒与坡面在Y方向的接触力整体上始终较小，表明降雨过程中颗粒与坡面的相互作用力较弱，流体作用力对土壤流失起主要作用。