河床透-阻型岩溶塌陷形成机理

余政兴,金福喜，段选亮

(1.中南大学地球科学与信息物理学院，湖南 长沙 410083； 2.有色金属成矿预测与地质环境监测教育部重点实验室(中南大学)，湖南 长沙 410083)

0 引言

河床岩溶塌陷的情况时有发生。1984年，铜陵市新桥矿区的疏干排水试验，造成圣冲河新、老河道以及附近水塘内出现30多处塌陷，河水水位大幅下降、水塘漏失，随后附近房屋大量开裂[1]。1995年，临沂市兰山药材批发市场小涑河河床内出现多个圆形塌坑，造成小涑河污水随塌陷坑倒灌污染地下水，调查发现塌陷东边和南边各有一眼抽水井，出水量120 m3/h[2-3]。河床塌陷造成河水倒灌，地下水污染，河水断流，影响附近人们的正常生产生活。同时也威胁水利设施，存在严重的安全隐患。鉴此，对河床岩溶塌陷情况应予以充分重视，应查明河床岩溶塌陷形成原因与塌陷机理，为地面塌陷地质灾害防治及规划工作提供依据。

长期以来，大量学者在岩溶地面塌陷的成因机理与预测评价等方面进行了大量研究。程星等[4]依据覆盖层性质将岩溶塌陷分为单一阻水型、单一透水型、无盖层型、阻—透型、透—阻型、透—阻—透型、阻—透—阻型7种概化模型，并对每一种概化模型的致塌机制进行了分析和讨论；康彦仁[5]按塌陷产生的主导因素和受力状态，将岩溶地面塌陷划分为：重力致塌、潜蚀致塌、气爆致塌、真空吸蚀致塌、振动致塌、荷载致塌、溶蚀致塌和根蚀致塌，并提出了相应8种致塌模式。大量调查发现，岩溶塌陷主要诱发因素是降雨、地下水位波动。 XIAO H等[6]评价了降雨、地下水渗漏和地下水位差对覆盖型岩溶塌陷影响，发现佛罗里达州中部塌陷主要为矿山抽排水造成的地下水位差导致；苏添金等[7]建立了岩溶塌陷力学模型，并采用FLAC3D模拟了地下水位下降条件下岩溶土洞的发育过程，得出了岩溶塌陷系数K与地下水位差呈线性关系。吴庆华等[8]基于物理模型试验定量研究了砂土覆盖岩溶塌陷过程与机理。

本文在前人研究的基础上，对研究区岩溶塌陷机理进行定性分析，再从渗透力学、土力学等理论出发，进行相关公式推导与数值模拟，解释岩溶水位下降对河床透—阻型盖层岩溶土洞稳定性的影响，以期为区域覆盖型岩溶塌陷的预防工作提供参考。

1 研究区概况

1.1 岩溶塌陷基本情况

2012年11—12月，某地区河床分别发生2处塌陷(图1)，造成河水倒灌，塌陷区影响范围约5 000 m2，塌陷坑附近4户民房严重开裂，直接经济损失约300万元。2013年河床围堰建成后塌陷坑已回填。2013年9—10月，距河床塌陷仅200 m的河流阶地出现2处塌坑，2014年11月塌坑再次活化和发展，同时附近居民区房屋开始出现大面积开裂。截止2019年6月，该区域共开裂房屋105户，威胁人口达504人。由于区域内河床塌陷具有优先性，故研究河床塌陷的成因机理对区域内地质灾害防治工作布置具有重要意义。

图1 塌陷坑及部分河床Fig.1 Karst collapse and partial riverbed

1.2 地质环境概况

研究区主要为溶蚀侵蚀河谷地貌，由于下游有一拦水坝，塌陷处河面高程常年稳定在86.5 m。区内沟谷宽缓，河谷平坦开阔，地形坡度一般小于5°。河床基底为二叠系下统栖霞组(P1q)灰岩、泥质灰岩。据勘察孔资料，河床内第四系厚度在4.80～12.6 m，主要为卵石—粉质黏土二元结构土体(图2)。区内构造以断层为主，褶曲次之。主要断层逆断层F5位于研究区南部，长20 km，走向近东西；隐伏断层F8为F5断层北盘分布的北西向次级断层，位于区内的中部，走向北西，走向延展长1 260 m，该断层控制了区内冲沟的发育方向和地下岩溶的发育走向，导致灰岩溶洞水溃入附近的露采矿坑，是区内主要含水层向采坑直接充水的通道。隐伏断层F9为F5断层北盘分布的东西向次级断层，位于河流的南侧部，走向近东西，走向延展长约1 500 m，河床内两处塌陷均位于该断层附近。该断层基本控制了河流南岸附近地下水的流向，使地下水向隐伏断层F8汇集流动，最终河水与露采矿坑贯通。

图2 塌陷坑T1地质剖面图Fig.2 The geological profile of karst collapse 1-淤泥；2-素填土；3-卵石；4-粉质黏土；5-灰岩；6-隐伏断层

1.3 人类工程活动概况

距塌陷点直线距离1.2 km处有一露天采矿场，矿场自1980年开采至今已形成一个近东西向的长约1 600 m，宽约 600 m的大采坑。采场底部已延深至±0 m标高，垂直坑深达104 m，2012～2015年，每天抽水量在4×104～6×104m3。矿山开采前，地下水流向总体受地形控制，即由高处向低处汇流，由山坡向溪沟，再经溪沟汇入河流。开采后，地下水尤其是下层岩溶水开始由四周向矿坑汇集。在矿山开采的条件下，采坑排水形成了降落漏斗，据有关勘察单位抽水试验结果计算，降落漏斗半径为2 200 m。塌陷处河水与矿坑坑底水位相差86.5 m，水平距离约1 200 m，水力坡降达72.1‰。同时勘察发现F8附近岩溶水位沿矿坑方向逐渐降低，说明河水通过F8断层大量排泄到矿坑。水动力条件的改变，使地下水流场变得复杂，以矿坑为中心汇集的岩溶地下水被人工大量排出，导致区域内发生岩溶塌陷。由于大量抽排水造成附近多处塌陷，该矿场已于2016年1月停产，目前由于地下水补给，矿坑形成巨大水库，水库水位基本与河水位平齐。2016—2019年，由于地下水位回升，周围未新增塌陷。塌陷坑、断裂带、矿坑相对位置关系见图3。

2 研究区河床岩溶塌陷致塌机理分析

2.1 塌陷分布特征与发育模式

在时间分布上，两处河床岩溶塌陷均发生在2012年11—12月，两处阶地塌陷发生于2013年9—10月，复活于2014年11月。故区内10—12月为岩溶塌陷高发期，且阶地塌陷明显滞后于河床塌陷。根据对阶地塌陷坑T4附近的岩溶水位动态监测数据分析(图4)，在两处阶地塌陷复活的2014年11月，地下水位监测钻孔K1测得岩溶水水位为7.30～11.40 m，而该处土层厚度约9 m。可知岩溶水在基岩附近波动时期更易发生地面塌陷。该地区岩溶水位则受降雨以及矿山抽排水共同控制。

空间分布上，4处岩溶塌陷均分布于二叠系下统栖霞组(P1q)质纯层厚碳酸盐岩分布区，且沿断裂带F8、F9分布；在地形上，区内塌陷与房屋开裂均发育在相对低洼阶地、谷地，标高85～98 m；同时，塌陷均发生于第四系相对较薄处，覆盖层厚度为6～10 m；区内塌陷的分布受地下水径流作用的控制明显，塌陷与房屋开裂均位于磷矿抽排水降落漏斗影响范围内。

综合上述分析，该地区河床岩溶塌陷发生的主要因素有：(1)区内碳酸盐质纯层厚，岩溶发育强烈。F8、F9断裂为导水性断层，断层附近岩溶更为发育，致使河水与岩溶地下水有了直接的水力联系塌陷，并且塌陷处有开口型溶洞，存在岩溶塌陷发育的条件；(2)矿山长期疏排岩溶水导致区域水文地质条件发生变化，塌陷区由天然状态下的岩溶水排泄区变成了河水补给岩溶水的补给区(图5)。在采坑大量抽排水的情况下，河水与采坑之间形成较大水力梯度，造成河水大量渗流补给岩溶水，渗透力增大，加速塌陷产生；(3)由于河水对土体的侵蚀、潜蚀作用，河床内覆盖层较薄，且土体在河水的浸润下力学性质变差，易于塌陷发生。在综合因素的作用下，可溶岩上覆土体的力学平衡遭到破坏，最终导致溶洞上覆土体的抗塌力小于下塌力，土体失稳陷落而形成岩溶塌陷。

图3 研究区岩溶塌陷分布图Fig.3 The distribution of karst collapse in research area

图4 2014年岩溶水位动态监测曲线Fig.4 Dynamic monitoring curve of karst water level in 2014

图5 研究区水文地质剖面图Fig.5 Hydrogeologic profile of research area 1-第四系覆盖层；2-二叠系下统栖霞组第三段；3-二叠系下统栖霞组 第二段；4-二叠系下统栖霞组第一段；5-震旦系上统；6-断裂带

该地区河床岩溶塌陷的破坏模式为“潜蚀—失托—重力致塌”型塌陷。即第四系覆盖层在岩溶水位下降后浮托力减小或消失，在地下水流的渗透作用下，使土层产生潜蚀或流失，形成土洞，土层厚度减小。同时上层河水不断补给岩溶水并冲刷土体，致使溶洞上覆的土体塌落，并在自重作用下形成岩溶地面塌陷。

其岩溶地面塌陷可以概括为如下地质过程：(1)初始状态：可溶岩中岩溶通道与上覆土体连通，但由于岩溶地下水位较高，与河水未形成相对渗流，岩溶通道附近土体保持稳定。(2)塌陷孕育：由于矿山抽排水导致岩溶水位下降，岩溶地下水与河水间形成相对渗流，岩溶通道开口处土体逐渐被水流冲刷带走，土洞初步形成。同时在土岩界面可溶岩开口处的集中渗流效应使得土洞边界不断扩大，土洞开始向上方延伸。(3)塌陷形成：随着土洞不断扩大，土洞顶板不能承受自重与河水重量，开始向下迅速塌落，塌落土体堆积在土洞中并最终形成塌陷坑。(4)塌陷扩大：塌陷坑形成后，河水与岩溶通道完全贯通，河水不断向岩溶通道内补给，同时冲刷并带走通道附近土体，塌陷范围不断扩大。

2.2 塌陷力学分析

该地区塌陷点附近覆盖层厚度较薄，为卵石—粉质黏土二元结构土体。根据勘察资料，河床塌陷T1、T2上部卵石层(淤泥层为后期河水冲刷产生)厚度分别为2 m、2.5 m，下部粉质黏土层厚度分别为7 m、7.5 m，符合程星等[4]提出的7种模型中的透—阻型盖层地质概化模型。对于该类地质模型，邢宇健[9]将土洞形态概化为半球体，考虑了水土自重，侧壁摩擦力等，根据极限平衡理论对某一时刻土洞的受力状态进行分析，但尚未考虑地表水渗透对土洞产生的作用力。而陶小虎等[10]比较了潜水位上升与承压水位下降对阻水层渗透稳定性的影响后发现，水对饱和土壤中土体的作用力主要受渗透坡降的影响。由于矿山开采时间较长(20世纪80年代开始开采)，岩溶地下水下降速度较慢，土洞内负压有足够时间消散，故不考虑真空吸蚀作用。为方便力学分析，作出以下假设：(1)河水仅对土体有向下方向的渗透力，不产生水平方向的冲刷力；(2)土洞形态为理想半球体；(3)土洞为完全空腔，内部气压与大气压力相等；(4)不考虑“土拱”作用。综合考虑土体自重、河水荷载、侧壁摩擦力、渗透力等因素，可得该区域岩溶地面塌陷的力学模型(图6)。

图6 河床岩溶地面塌陷力学模型Fig.6 The mechanical model of karst collapse in riverbed

假设土洞为半球形，水、两层土重度分别为γw、γ1、γ2(其值等于饱和重度)，土洞所受水土自重为：

(1)

式中：G——土洞上部水土体总重/kg；

G1——土洞上部水体重量/kg；

G2——土洞上部第一层土重量/kg；

G3——土洞上部第二层土重量/kg；

D——土洞直径/m；

γw——河水重度/(kg·m-1)；

γ1——第一层土体重度/(kg·m-1)；

γ2——第二层土体重度/(kg·m-1)；

h——河水厚度/m；

H1——第一层土厚度/m；

H2——第二次土厚度/m。

河水对土洞产生的渗透力包括两部分，首先由于土洞内水体排出空气进入，浮力减小，这一过程相对于给土洞顶板施加一个与浮托力大小相同，方向相反的力，其次水流下渗对土体仅产生同方向的拖拽力[9-10]，根据周建等[11]推导结果，该力可由式(2)表示：

(2)

式中：uw——渗透流速，m/s，通常小于10-4m/s，故第二部分计算中可忽略不计。

土洞顶板所受向上荷载主要是土洞顶板侧面产生的静摩擦力，该力主要由土的侧压力产生，借鉴库伦土压力理论，假设外部土体对影响区域内的土体作用处于主动极限平衡状态，则两层土水平侧压力分别为：

(3)

式中：σx1——第一层土水平侧压力；

Ka1——第一层层土主动土压力系数，Ka1=tan2(45°-φ1/2)；

C1——第一层土黏聚力/kPa；

φ1——第一层土内摩擦角/(°)。

σx2=Ka2(γwH2h+γ1H1H2)+(1/2)Ka2γ2H22-

(4)

式中：σx2——第二层土水平侧压力；

Ka2——第二层层土主动土压力系数，Ka2=tan2(45°-φ2/2)；

C2——第二层土黏聚力/kPa；

φ2——第二层土内摩擦角/(°)。

最大静摩擦力为：

f=f1+f2=πDH1(σx1tanφ1+C1)+

πDH2(σx2tanφ2+C2)

(5)

式中：f——最大静摩擦力/kN；

设Ks为稳定系数，则：

(6)

式中：Ks——稳定系数。

若第一层土为砂卵石，则C1=0。将(1)(2)(5)式带入(6)式可求得Ks。当稳定系数Ks<1时，则可能产生地面塌陷，稳定系数Ks=1时处于塌陷临界状态，Ks>1时土洞稳定。

将勘察报告提供的塌陷坑T1岩土参数表1带入(6)式。在塌陷临界状态下，即Ks=1时，得D=4.2 m，即该塌陷坑临界土洞高度仅为2.1 m，说明在土洞形成初期就发生塌陷，这也与实际情况相符合。但由于未考虑土拱效应等影响，因此结果比较保守，需要用数值模拟来进一步研究。

表1 水土参数取值

3 基于数值模拟的岩溶土洞塌陷过程分析

对于岩溶塌陷的数值模拟研究，多采用FLAC3D有限差分数值软件[12-14]。针对研究区河床岩溶塌陷，利用FLAC3D进行数值模拟和预测评价，主要模拟地下水对覆盖层潜蚀作用引起的岩溶土洞扩大的过程，分析土洞大小变化以及岩溶水位下降对地面稳定性的影响，进而进行岩溶塌陷发展预测与评价。

3.1 数值模型的建立

依照地区塌陷附近实际情况设计两组数值模拟方案，其中模型二为对照组，土体c、φ值均较大，以模拟阶地岩溶塌陷，两种模型土体参数均由勘察报告提供(表2)。

模型一：覆盖层土体主要分为2层: 卵石2 m厚，粉质黏土7 m厚，由此进行简化建立数值模型，模型尺寸为18 m×9 m×9 m，上部施加20 kPa的均布荷载，等效于2 m河水的自重荷载。由于河水作用，设置土体顶面孔隙水压力为20 kPa，底部孔隙水压力为110 kPa。

模型二：覆盖层土体主要分为3层: 素填土1 m厚(等效于2 m河水荷载)，卵石2 m厚，粉质黏土7 m厚，模型尺寸为20 m×10 m×10 m。初始状态下，假定地下水位在土体顶面，则土体顶面孔隙水压力为0 kPa，底部孔隙水压力为100 kPa。

在河水位为2 m情况下，设置两者初始土洞高度均为1.5 m，土洞按ΔH=0.5 m速率往各方向等速扩大，最终直到土洞破坏，得到临界土洞大小。由于土体中的渗流为持续过程，渗流产生的托拽力仅对土洞的发育过程影响较大，而对某一时刻土洞稳定性的影响较小，故对临界土洞模拟时关闭渗流功能，仅考虑自重效应以及浮托力消失作用。而在临界土洞条件下，再开启FLAC3D渗流计算功能，对抽水前后的孔隙水压力、应力、位移等变化情况进行模拟分析。

本次土体模型选用莫尔-库仑模型，计算满足莫尔-库仑破坏准则，渗流模拟满足达西定律。

表2 数值模拟土体参数取值Table 2 The parameter of soil in numerical simulation

3.2 计算成果分析

首先研究不同土洞高度在岩溶水位下降后的稳定性，其中模型1(图7)，共两层土体，划分网格数11 664个，底部挖空以模拟天然土洞。由于本次计算未限定计算步数，故计算后需使不平衡力达到收敛。最大不平衡力随着迭代步数的增加初始不平衡力逐渐减小(图8)，减小至1×10-6Pa，趋近于0，说明在天然重力场条件下，模型通过自身的调整可以达到自平衡，天然应力场顺利形成。

图7 三维模型网格划分Fig.7 Mesh of the 3D numerical model

图8 最大不平衡力曲线Fig.8 Maximum unbalanced force curve

地面塌陷发生的标志为土洞周围塑性区与地表贯通。在模型1中，截取模型中心竖直方向剖面(图9)，塑性区首先出现在土洞顶板两侧，随后在地表卵石层也有出现，这是由于卵石层塑性较低，下部土体变形后易产生扰动。随着土洞不断扩大，塑性区的范围也越来越大，溶洞顶板处的塑性区也不断地增大。在土洞高度达到2.5 m时，土洞顶板塑性区与地面贯通，塌陷发生。故模型1中临界土洞高度在2～2.5 m，与(6)式计算得到的临界土洞高度2.1 m的数值较吻合。而在模型2中当土洞高度达到3.5 m时塑性区才出现贯通(图10)，这也佐证相似地层条件下，河床岩溶塌陷较陆面塌陷更易发生。

采用模型1处于临界土洞(R=2.5 m)时状态进行渗流模拟，定义初始水位为土体顶部以上2 m，渗流未开始时，不存在超孔隙水压力。未抽水条件下，土体内未产生渗流，最大孔隙水压力出现在模型底部，为110 kPa，受河水作用，顶部孔隙水压力为20 kPa(图11)。

把土洞底部设置为透水边界，以模拟岩溶水位下降到基岩面以下。选取距离模型底部5 m 处剖面研究孔隙水压力的变化，截取运算时步分别为 500、1 000 和2 000 时步的孔隙水压力图(图12)。分析发现：随着渗流进行，中心位置孔隙水压力逐渐降低，且越靠近中心部分，孔隙水压力越小，水位下降越快；同时降落漏斗形成明显，且影响范围逐渐增大。

在剖面中心位置设置孔隙水压力监测点，得到降落漏斗中心水位降深与运算时步对应曲线(图13)，可知在计算时步达到2 000后孔隙水压力变化较小，说明土体此时已发生塑性破坏，土体内部已形成贯通的管流通道。

图9 模型1土体塑性区分布Fig.9 Distribution of plastic zone in model 1

图10 模型2土体塑性区分布(R=3.5 m)Fig.10 Distribution of plastic zone in model 2 (R=3.5 m)

图11 抽水前孔隙水压力分布Fig.11 Pore water pressure distribution before pumping

图12 渗流过程孔隙水压力剖面Fig.12 Pore water pressure profile during seepage

图13 中心孔隙水压力—计算时步对应曲线Fig.13 Center pore water pressure—time step corresponding curve

图14 模型1(R=2.5 m)位移云图Fig.14 Displacement nephogram in model 1(R=2.5 m)

图15 土洞正上方地表位移监测曲线Fig.15 Surface displacement monitoring curve

图14为3 200时步的位移云图。可看出土洞正上方土体变形大于周围土体变形，说明该土洞先由中心塌陷，周围土体在水流与自重的影响下随后塌陷，以致塌陷范围不断扩大。为更直观的展现变形变化规律，在土洞正上方地表土体位置设置监测点，监测曲线如图15所示。随着计算步数增加以及渗流进行，地表位移量逐步增大，在500时步后位移增长明显放缓，3 200步时位移已达55 cm。

图16为3 200时步的的应力云图，由图16可知，土洞形成后，地下水渗流会对周围土体应力场造成影响，导致靠近土洞的土体，应力变化更大。这主要受两个方面影响：首先是由于土洞的形成和发展使得周围土体受到的初始压应力得到释放，以土洞为中心应力重分布。再者由于地下水渗流与重力场共同作用，在土洞顶部形成一个呈漏斗形的应力松弛带，由内向外变化逐渐减小。土洞顶部附近土体由高应力转变为低应力或零应力，导致土体内部发生剪切破坏，最终塌陷发生。

图16 模型1(R=2.5 m)应力云图Fig.16 Stress nephogram in model 1(R=2.5 m)

4 结论

(1)该地区河床产生岩溶塌陷主要由于露采矿山大量抽排水使河水与采坑之间水力梯度增大，河水下渗产生强烈潜蚀效应，使土洞形成并逐步向上扩展，当土洞顶板不能支撑自重时，发生地面塌陷。

(2)理论计算表明岩溶塌陷稳定系数Ks主要受土层厚度、土体力学性质、土洞高度、岩溶水下降幅度等影响。计算得出研究区河床塌陷临界土洞高度为2.1 m，与实际情况较符。

(3)数值模拟结果表明该地区河床塌陷的临界土洞高度为2～2.5 m，小于该地区阶地岩溶地面塌陷的临界土洞高度3.5 m。由于缺少上层滞水的渗透作用，阶地土洞的扩展速度也较慢，能够解释河床塌陷发生后1～2年内周围200 m内陆续出现地面塌陷的原因。

(4)渗流模拟结果证明了岩溶水下降后潜蚀效应明显，区内岩溶塌陷产生的根本原因是水动力条件发生变化。故在岩溶地区抽排地下水时，应在控制最大降深，还需要控制抽排水速度，防止出现土体渗透破坏。

(5)由于河床岩溶塌陷可作为大规模阶地塌陷发生的前兆，故对类似区域的河床、水塘发生的塌陷应引起重视。在河床塌陷发生后应及时对塌陷坑进行注浆或回填，防止岩溶管道被疏通；宜立即停止周围大规模抽排水活动；需对周围地下水变化进行长期监测，地下水在岩土交界面波动时期尤其注意防范地面塌陷。