余政興,金福喜，段選亮

(1.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083； 2.有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室(中南大學(xué))，湖南 長(zhǎng)沙 410083)

0 引言

河床巖溶塌陷的情況時(shí)有發(fā)生。1984年，銅陵市新橋礦區(qū)的疏干排水試驗(yàn)，造成圣沖河新、老河道以及附近水塘內(nèi)出現(xiàn)30多處塌陷，河水水位大幅下降、水塘漏失，隨后附近房屋大量開裂[1]。1995年，臨沂市蘭山藥材批發(fā)市場(chǎng)小涑河河床內(nèi)出現(xiàn)多個(gè)圓形塌坑，造成小涑河污水隨塌陷坑倒灌污染地下水，調(diào)查發(fā)現(xiàn)塌陷東邊和南邊各有一眼抽水井，出水量120 m3/h[2-3]。河床塌陷造成河水倒灌，地下水污染，河水?dāng)嗔?，影響附近人們的正常生產(chǎn)生活。同時(shí)也威脅水利設(shè)施，存在嚴(yán)重的安全隱患。鑒此，對(duì)河床巖溶塌陷情況應(yīng)予以充分重視，應(yīng)查明河床巖溶塌陷形成原因與塌陷機(jī)理，為地面塌陷地質(zhì)災(zāi)害防治及規(guī)劃工作提供依據(jù)。

長(zhǎng)期以來(lái)，大量學(xué)者在巖溶地面塌陷的成因機(jī)理與預(yù)測(cè)評(píng)價(jià)等方面進(jìn)行了大量研究。程星等[4]依據(jù)覆蓋層性質(zhì)將巖溶塌陷分為單一阻水型、單一透水型、無(wú)蓋層型、阻—透型、透—阻型、透—阻—透型、阻—透—阻型7種概化模型，并對(duì)每一種概化模型的致塌機(jī)制進(jìn)行了分析和討論；康彥仁[5]按塌陷產(chǎn)生的主導(dǎo)因素和受力狀態(tài)，將巖溶地面塌陷劃分為：重力致塌、潛蝕致塌、氣爆致塌、真空吸蝕致塌、振動(dòng)致塌、荷載致塌、溶蝕致塌和根蝕致塌，并提出了相應(yīng)8種致塌模式。大量調(diào)查發(fā)現(xiàn)，巖溶塌陷主要誘發(fā)因素是降雨、地下水位波動(dòng)。 XIAO H等[6]評(píng)價(jià)了降雨、地下水滲漏和地下水位差對(duì)覆蓋型巖溶塌陷影響，發(fā)現(xiàn)佛羅里達(dá)州中部塌陷主要為礦山抽排水造成的地下水位差導(dǎo)致；蘇添金等[7]建立了巖溶塌陷力學(xué)模型，并采用FLAC3D模擬了地下水位下降條件下巖溶土洞的發(fā)育過(guò)程，得出了巖溶塌陷系數(shù)K與地下水位差呈線性關(guān)系。吳慶華等[8]基于物理模型試驗(yàn)定量研究了砂土覆蓋巖溶塌陷過(guò)程與機(jī)理。

本文在前人研究的基礎(chǔ)上，對(duì)研究區(qū)巖溶塌陷機(jī)理進(jìn)行定性分析，再?gòu)臐B透力學(xué)、土力學(xué)等理論出發(fā)，進(jìn)行相關(guān)公式推導(dǎo)與數(shù)值模擬，解釋巖溶水位下降對(duì)河床透—阻型蓋層巖溶土洞穩(wěn)定性的影響，以期為區(qū)域覆蓋型巖溶塌陷的預(yù)防工作提供參考。

1 研究區(qū)概況

1.1 巖溶塌陷基本情況

2012年11—12月，某地區(qū)河床分別發(fā)生2處塌陷(圖1)，造成河水倒灌，塌陷區(qū)影響范圍約5 000 m2，塌陷坑附近4戶民房嚴(yán)重開裂，直接經(jīng)濟(jì)損失約300萬(wàn)元。2013年河床圍堰建成后塌陷坑已回填。2013年9—10月，距河床塌陷僅200 m的河流階地出現(xiàn)2處塌坑，2014年11月塌坑再次活化和發(fā)展，同時(shí)附近居民區(qū)房屋開始出現(xiàn)大面積開裂。截止2019年6月，該區(qū)域共開裂房屋105戶，威脅人口達(dá)504人。由于區(qū)域內(nèi)河床塌陷具有優(yōu)先性，故研究河床塌陷的成因機(jī)理對(duì)區(qū)域內(nèi)地質(zhì)災(zāi)害防治工作布置具有重要意義。

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

1.2 地質(zhì)環(huán)境概況

研究區(qū)主要為溶蝕侵蝕河谷地貌，由于下游有一攔水壩，塌陷處河面高程常年穩(wěn)定在86.5 m。區(qū)內(nèi)溝谷寬緩，河谷平坦開闊，地形坡度一般小于5°。河床基底為二疊系下統(tǒng)棲霞組(P1q)灰?guī)r、泥質(zhì)灰?guī)r。據(jù)勘察孔資料，河床內(nèi)第四系厚度在4.80～12.6 m，主要為卵石—粉質(zhì)黏土二元結(jié)構(gòu)土體(圖2)。區(qū)內(nèi)構(gòu)造以斷層為主，褶曲次之。主要斷層逆斷層F5位于研究區(qū)南部，長(zhǎng)20 km，走向近東西；隱伏斷層F8為F5斷層北盤分布的北西向次級(jí)斷層，位于區(qū)內(nèi)的中部，走向北西，走向延展長(zhǎng)1 260 m，該斷層控制了區(qū)內(nèi)沖溝的發(fā)育方向和地下巖溶的發(fā)育走向，導(dǎo)致灰?guī)r溶洞水潰入附近的露采礦坑，是區(qū)內(nèi)主要含水層向采坑直接充水的通道。隱伏斷層F9為F5斷層北盤分布的東西向次級(jí)斷層，位于河流的南側(cè)部，走向近東西，走向延展長(zhǎng)約1 500 m，河床內(nèi)兩處塌陷均位于該斷層附近。該斷層基本控制了河流南岸附近地下水的流向，使地下水向隱伏斷層F8匯集流動(dòng)，最終河水與露采礦坑貫通。

圖2 塌陷坑T1地質(zhì)剖面圖Fig.2 The geological profile of karst collapse 1-淤泥；2-素填土；3-卵石；4-粉質(zhì)黏土；5-灰?guī)r；6-隱伏斷層

1.3 人類工程活動(dòng)概況

距塌陷點(diǎn)直線距離1.2 km處有一露天采礦場(chǎng)，礦場(chǎng)自1980年開采至今已形成一個(gè)近東西向的長(zhǎng)約1 600 m，寬約 600 m的大采坑。采場(chǎng)底部已延深至±0 m標(biāo)高，垂直坑深達(dá)104 m，2012～2015年，每天抽水量在4×104～6×104m3。礦山開采前，地下水流向總體受地形控制，即由高處向低處匯流，由山坡向溪溝，再經(jīng)溪溝匯入河流。開采后，地下水尤其是下層巖溶水開始由四周向礦坑匯集。在礦山開采的條件下，采坑排水形成了降落漏斗，據(jù)有關(guān)勘察單位抽水試驗(yàn)結(jié)果計(jì)算，降落漏斗半徑為2 200 m。塌陷處河水與礦坑坑底水位相差86.5 m，水平距離約1 200 m，水力坡降達(dá)72.1‰。同時(shí)勘察發(fā)現(xiàn)F8附近巖溶水位沿礦坑方向逐漸降低，說(shuō)明河水通過(guò)F8斷層大量排泄到礦坑。水動(dòng)力條件的改變，使地下水流場(chǎng)變得復(fù)雜，以礦坑為中心匯集的巖溶地下水被人工大量排出，導(dǎo)致區(qū)域內(nèi)發(fā)生巖溶塌陷。由于大量抽排水造成附近多處塌陷，該礦場(chǎng)已于2016年1月停產(chǎn)，目前由于地下水補(bǔ)給，礦坑形成巨大水庫(kù)，水庫(kù)水位基本與河水位平齊。2016—2019年，由于地下水位回升，周圍未新增塌陷。塌陷坑、斷裂帶、礦坑相對(duì)位置關(guān)系見圖3。

2 研究區(qū)河床巖溶塌陷致塌機(jī)理分析

2.1 塌陷分布特征與發(fā)育模式

在時(shí)間分布上，兩處河床巖溶塌陷均發(fā)生在2012年11—12月，兩處階地塌陷發(fā)生于2013年9—10月，復(fù)活于2014年11月。故區(qū)內(nèi)10—12月為巖溶塌陷高發(fā)期，且階地塌陷明顯滯后于河床塌陷。根據(jù)對(duì)階地塌陷坑T4附近的巖溶水位動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)分析(圖4)，在兩處階地塌陷復(fù)活的2014年11月，地下水位監(jiān)測(cè)鉆孔K1測(cè)得巖溶水水位為7.30～11.40 m，而該處土層厚度約9 m?？芍獛r溶水在基巖附近波動(dòng)時(shí)期更易發(fā)生地面塌陷。該地區(qū)巖溶水位則受降雨以及礦山抽排水共同控制。

空間分布上，4處巖溶塌陷均分布于二疊系下統(tǒng)棲霞組(P1q)質(zhì)純層厚碳酸鹽巖分布區(qū)，且沿?cái)嗔褞8、F9分布；在地形上，區(qū)內(nèi)塌陷與房屋開裂均發(fā)育在相對(duì)低洼階地、谷地，標(biāo)高85～98 m；同時(shí)，塌陷均發(fā)生于第四系相對(duì)較薄處，覆蓋層厚度為6～10 m；區(qū)內(nèi)塌陷的分布受地下水徑流作用的控制明顯，塌陷與房屋開裂均位于磷礦抽排水降落漏斗影響范圍內(nèi)。

綜合上述分析，該地區(qū)河床巖溶塌陷發(fā)生的主要因素有：(1)區(qū)內(nèi)碳酸鹽質(zhì)純層厚，巖溶發(fā)育強(qiáng)烈。F8、F9斷裂為導(dǎo)水性斷層，斷層附近巖溶更為發(fā)育，致使河水與巖溶地下水有了直接的水力聯(lián)系塌陷，并且塌陷處有開口型溶洞，存在巖溶塌陷發(fā)育的條件；(2)礦山長(zhǎng)期疏排巖溶水導(dǎo)致區(qū)域水文地質(zhì)條件發(fā)生變化，塌陷區(qū)由天然狀態(tài)下的巖溶水排泄區(qū)變成了河水補(bǔ)給巖溶水的補(bǔ)給區(qū)(圖5)。在采坑大量抽排水的情況下，河水與采坑之間形成較大水力梯度，造成河水大量滲流補(bǔ)給巖溶水，滲透力增大，加速塌陷產(chǎn)生；(3)由于河水對(duì)土體的侵蝕、潛蝕作用，河床內(nèi)覆蓋層較薄，且土體在河水的浸潤(rùn)下力學(xué)性質(zhì)變差，易于塌陷發(fā)生。在綜合因素的作用下，可溶巖上覆土體的力學(xué)平衡遭到破壞，最終導(dǎo)致溶洞上覆土體的抗塌力小于下塌力，土體失穩(wěn)陷落而形成巖溶塌陷。

圖3 研究區(qū)巖溶塌陷分布圖Fig.3 The distribution of karst collapse in research area

圖4 2014年巖溶水位動(dòng)態(tài)監(jiān)測(cè)曲線Fig.4 Dynamic monitoring curve of karst water level in 2014

圖5 研究區(qū)水文地質(zhì)剖面圖Fig.5 Hydrogeologic profile of research area 1-第四系覆蓋層；2-二疊系下統(tǒng)棲霞組第三段；3-二疊系下統(tǒng)棲霞組 第二段；4-二疊系下統(tǒng)棲霞組第一段；5-震旦系上統(tǒng)；6-斷裂帶

該地區(qū)河床巖溶塌陷的破壞模式為“潛蝕—失托—重力致塌”型塌陷。即第四系覆蓋層在巖溶水位下降后浮托力減小或消失，在地下水流的滲透作用下，使土層產(chǎn)生潛蝕或流失，形成土洞，土層厚度減小。同時(shí)上層河水不斷補(bǔ)給巖溶水并沖刷土體，致使溶洞上覆的土體塌落，并在自重作用下形成巖溶地面塌陷。

其巖溶地面塌陷可以概括為如下地質(zhì)過(guò)程：(1)初始狀態(tài)：可溶巖中巖溶通道與上覆土體連通，但由于巖溶地下水位較高，與河水未形成相對(duì)滲流，巖溶通道附近土體保持穩(wěn)定。(2)塌陷孕育：由于礦山抽排水導(dǎo)致巖溶水位下降，巖溶地下水與河水間形成相對(duì)滲流，巖溶通道開口處土體逐漸被水流沖刷帶走，土洞初步形成。同時(shí)在土巖界面可溶巖開口處的集中滲流效應(yīng)使得土洞邊界不斷擴(kuò)大，土洞開始向上方延伸。(3)塌陷形成：隨著土洞不斷擴(kuò)大，土洞頂板不能承受自重與河水重量，開始向下迅速塌落，塌落土體堆積在土洞中并最終形成塌陷坑。(4)塌陷擴(kuò)大：塌陷坑形成后，河水與巖溶通道完全貫通，河水不斷向巖溶通道內(nèi)補(bǔ)給，同時(shí)沖刷并帶走通道附近土體，塌陷范圍不斷擴(kuò)大。

2.2 塌陷力學(xué)分析

該地區(qū)塌陷點(diǎn)附近覆蓋層厚度較薄，為卵石—粉質(zhì)黏土二元結(jié)構(gòu)土體。根據(jù)勘察資料，河床塌陷T1、T2上部卵石層(淤泥層為后期河水沖刷產(chǎn)生)厚度分別為2 m、2.5 m，下部粉質(zhì)黏土層厚度分別為7 m、7.5 m，符合程星等[4]提出的7種模型中的透—阻型蓋層地質(zhì)概化模型。對(duì)于該類地質(zhì)模型，邢宇健[9]將土洞形態(tài)概化為半球體，考慮了水土自重，側(cè)壁摩擦力等，根據(jù)極限平衡理論對(duì)某一時(shí)刻土洞的受力狀態(tài)進(jìn)行分析，但尚未考慮地表水滲透對(duì)土洞產(chǎn)生的作用力。而陶小虎等[10]比較了潛水位上升與承壓水位下降對(duì)阻水層滲透穩(wěn)定性的影響后發(fā)現(xiàn)，水對(duì)飽和土壤中土體的作用力主要受滲透坡降的影響。由于礦山開采時(shí)間較長(zhǎng)(20世紀(jì)80年代開始開采)，巖溶地下水下降速度較慢，土洞內(nèi)負(fù)壓有足夠時(shí)間消散，故不考慮真空吸蝕作用。為方便力學(xué)分析，作出以下假設(shè)：(1)河水僅對(duì)土體有向下方向的滲透力，不產(chǎn)生水平方向的沖刷力；(2)土洞形態(tài)為理想半球體；(3)土洞為完全空腔，內(nèi)部氣壓與大氣壓力相等；(4)不考慮“土拱”作用。綜合考慮土體自重、河水荷載、側(cè)壁摩擦力、滲透力等因素，可得該區(qū)域巖溶地面塌陷的力學(xué)模型(圖6)。

圖6 河床巖溶地面塌陷力學(xué)模型Fig.6 The mechanical model of karst collapse in riverbed

假設(shè)土洞為半球形，水、兩層土重度分別為γ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。

河水對(duì)土洞產(chǎn)生的滲透力包括兩部分，首先由于土洞內(nèi)水體排出空氣進(jìn)入，浮力減小，這一過(guò)程相對(duì)于給土洞頂板施加一個(gè)與浮托力大小相同，方向相反的力，其次水流下滲對(duì)土體僅產(chǎn)生同方向的拖拽力[9-10]，根據(jù)周建等[11]推導(dǎo)結(jié)果，該力可由式(2)表示：

(2)

式中：uw——滲透流速，m/s，通常小于10-4m/s，故第二部分計(jì)算中可忽略不計(jì)。

土洞頂板所受向上荷載主要是土洞頂板側(cè)面產(chǎn)生的靜摩擦力，該力主要由土的側(cè)壓力產(chǎn)生，借鑒庫(kù)倫土壓力理論，假設(shè)外部土體對(duì)影響區(qū)域內(nèi)的土體作用處于主動(dòng)極限平衡狀態(tài)，則兩層土水平側(cè)壓力分別為：

(3)

式中：σx1——第一層土水平側(cè)壓力；

Ka1——第一層層土主動(dòng)土壓力系數(shù)，Ka1=tan2(45°-φ1/2)；

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

φ1——第一層土內(nèi)摩擦角/(°)。

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

(4)

式中：σx2——第二層土水平側(cè)壓力；

Ka2——第二層層土主動(dòng)土壓力系數(shù)，Ka2=tan2(45°-φ2/2)；

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

φ2——第二層土內(nèi)摩擦角/(°)。

最大靜摩擦力為：

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

πDH2(σx2tanφ2+C2)

(5)

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

設(shè)Ks為穩(wěn)定系數(shù)，則：

(6)

式中：Ks——穩(wěn)定系數(shù)。

若第一層土為砂卵石，則C1=0。將(1)(2)(5)式帶入(6)式可求得Ks。當(dāng)穩(wěn)定系數(shù)Ks<1時(shí)，則可能產(chǎn)生地面塌陷，穩(wěn)定系數(shù)Ks=1時(shí)處于塌陷臨界狀態(tài)，Ks>1時(shí)土洞穩(wěn)定。

將勘察報(bào)告提供的塌陷坑T1巖土參數(shù)表1帶入(6)式。在塌陷臨界狀態(tài)下，即Ks=1時(shí)，得D=4.2 m，即該塌陷坑臨界土洞高度僅為2.1 m，說(shuō)明在土洞形成初期就發(fā)生塌陷，這也與實(shí)際情況相符合。但由于未考慮土拱效應(yīng)等影響，因此結(jié)果比較保守，需要用數(shù)值模擬來(lái)進(jìn)一步研究。

表1 水土參數(shù)取值

3 基于數(shù)值模擬的巖溶土洞塌陷過(guò)程分析

對(duì)于巖溶塌陷的數(shù)值模擬研究，多采用FLAC3D有限差分?jǐn)?shù)值軟件[12-14]。針對(duì)研究區(qū)河床巖溶塌陷，利用FLAC3D進(jìn)行數(shù)值模擬和預(yù)測(cè)評(píng)價(jià)，主要模擬地下水對(duì)覆蓋層潛蝕作用引起的巖溶土洞擴(kuò)大的過(guò)程，分析土洞大小變化以及巖溶水位下降對(duì)地面穩(wěn)定性的影響，進(jìn)而進(jìn)行巖溶塌陷發(fā)展預(yù)測(cè)與評(píng)價(jià)。

3.1 數(shù)值模型的建立

依照地區(qū)塌陷附近實(shí)際情況設(shè)計(jì)兩組數(shù)值模擬方案，其中模型二為對(duì)照組，土體c、φ值均較大，以模擬階地巖溶塌陷，兩種模型土體參數(shù)均由勘察報(bào)告提供(表2)。

模型一：覆蓋層土體主要分為2層: 卵石2 m厚，粉質(zhì)黏土7 m厚，由此進(jìn)行簡(jiǎn)化建立數(shù)值模型，模型尺寸為18 m×9 m×9 m，上部施加20 kPa的均布荷載，等效于2 m河水的自重荷載。由于河水作用，設(shè)置土體頂面孔隙水壓力為20 kPa，底部孔隙水壓力為110 kPa。

模型二：覆蓋層土體主要分為3層: 素填土1 m厚(等效于2 m河水荷載)，卵石2 m厚，粉質(zhì)黏土7 m厚，模型尺寸為20 m×10 m×10 m。初始狀態(tài)下，假定地下水位在土體頂面，則土體頂面孔隙水壓力為0 kPa，底部孔隙水壓力為100 kPa。

在河水位為2 m情況下，設(shè)置兩者初始土洞高度均為1.5 m，土洞按ΔH=0.5 m速率往各方向等速擴(kuò)大，最終直到土洞破壞，得到臨界土洞大小。由于土體中的滲流為持續(xù)過(guò)程，滲流產(chǎn)生的托拽力僅對(duì)土洞的發(fā)育過(guò)程影響較大，而對(duì)某一時(shí)刻土洞穩(wěn)定性的影響較小，故對(duì)臨界土洞模擬時(shí)關(guān)閉滲流功能，僅考慮自重效應(yīng)以及浮托力消失作用。而在臨界土洞條件下，再開啟FLAC3D滲流計(jì)算功能，對(duì)抽水前后的孔隙水壓力、應(yīng)力、位移等變化情況進(jìn)行模擬分析。

本次土體模型選用莫爾-庫(kù)侖模型，計(jì)算滿足莫爾-庫(kù)侖破壞準(zhǔn)則，滲流模擬滿足達(dá)西定律。

表2 數(shù)值模擬土體參數(shù)取值Table 2 The parameter of soil in numerical simulation

3.2 計(jì)算成果分析

首先研究不同土洞高度在巖溶水位下降后的穩(wěn)定性，其中模型1(圖7)，共兩層土體，劃分網(wǎng)格數(shù)11 664個(gè)，底部挖空以模擬天然土洞。由于本次計(jì)算未限定計(jì)算步數(shù)，故計(jì)算后需使不平衡力達(dá)到收斂。最大不平衡力隨著迭代步數(shù)的增加初始不平衡力逐漸減小(圖8)，減小至1×10-6Pa，趨近于0，說(shuō)明在天然重力場(chǎng)條件下，模型通過(guò)自身的調(diào)整可以達(dá)到自平衡，天然應(yīng)力場(chǎng)順利形成。

圖7 三維模型網(wǎng)格劃分Fig.7 Mesh of the 3D numerical model

圖8 最大不平衡力曲線Fig.8 Maximum unbalanced force curve

地面塌陷發(fā)生的標(biāo)志為土洞周圍塑性區(qū)與地表貫通。在模型1中，截取模型中心豎直方向剖面(圖9)，塑性區(qū)首先出現(xiàn)在土洞頂板兩側(cè)，隨后在地表卵石層也有出現(xiàn)，這是由于卵石層塑性較低，下部土體變形后易產(chǎn)生擾動(dòng)。隨著土洞不斷擴(kuò)大，塑性區(qū)的范圍也越來(lái)越大，溶洞頂板處的塑性區(qū)也不斷地增大。在土洞高度達(dá)到2.5 m時(shí)，土洞頂板塑性區(qū)與地面貫通，塌陷發(fā)生。故模型1中臨界土洞高度在2～2.5 m，與(6)式計(jì)算得到的臨界土洞高度2.1 m的數(shù)值較吻合。而在模型2中當(dāng)土洞高度達(dá)到3.5 m時(shí)塑性區(qū)才出現(xiàn)貫通(圖10)，這也佐證相似地層條件下，河床巖溶塌陷較陸面塌陷更易發(fā)生。

采用模型1處于臨界土洞(R=2.5 m)時(shí)狀態(tài)進(jìn)行滲流模擬，定義初始水位為土體頂部以上2 m，滲流未開始時(shí)，不存在超孔隙水壓力。未抽水條件下，土體內(nèi)未產(chǎn)生滲流，最大孔隙水壓力出現(xiàn)在模型底部，為110 kPa，受河水作用，頂部孔隙水壓力為20 kPa(圖11)。

把土洞底部設(shè)置為透水邊界，以模擬巖溶水位下降到基巖面以下。選取距離模型底部5 m 處剖面研究孔隙水壓力的變化，截取運(yùn)算時(shí)步分別為 500、1 000 和2 000 時(shí)步的孔隙水壓力圖(圖12)。分析發(fā)現(xiàn)：隨著滲流進(jìn)行，中心位置孔隙水壓力逐漸降低，且越靠近中心部分，孔隙水壓力越小，水位下降越快；同時(shí)降落漏斗形成明顯，且影響范圍逐漸增大。

在剖面中心位置設(shè)置孔隙水壓力監(jiān)測(cè)點(diǎn)，得到降落漏斗中心水位降深與運(yùn)算時(shí)步對(duì)應(yīng)曲線(圖13)，可知在計(jì)算時(shí)步達(dá)到2 000后孔隙水壓力變化較小，說(shuō)明土體此時(shí)已發(fā)生塑性破壞，土體內(nèi)部已形成貫通的管流通道。

圖9 模型1土體塑性區(qū)分布Fig.9 Distribution of plastic zone in model 1

圖10 模型2土體塑性區(qū)分布(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 滲流過(guò)程孔隙水壓力剖面Fig.12 Pore water pressure profile during seepage

圖13 中心孔隙水壓力—計(jì)算時(shí)步對(duì)應(yīng)曲線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 土洞正上方地表位移監(jiān)測(cè)曲線Fig.15 Surface displacement monitoring curve

圖14為3 200時(shí)步的位移云圖。可看出土洞正上方土體變形大于周圍土體變形，說(shuō)明該土洞先由中心塌陷，周圍土體在水流與自重的影響下隨后塌陷，以致塌陷范圍不斷擴(kuò)大。為更直觀的展現(xiàn)變形變化規(guī)律，在土洞正上方地表土體位置設(shè)置監(jiān)測(cè)點(diǎn)，監(jiān)測(cè)曲線如圖15所示。隨著計(jì)算步數(shù)增加以及滲流進(jìn)行，地表位移量逐步增大，在500時(shí)步后位移增長(zhǎng)明顯放緩，3 200步時(shí)位移已達(dá)55 cm。

圖16為3 200時(shí)步的的應(yīng)力云圖，由圖16可知，土洞形成后，地下水滲流會(huì)對(duì)周圍土體應(yīng)力場(chǎng)造成影響，導(dǎo)致靠近土洞的土體，應(yīng)力變化更大。這主要受兩個(gè)方面影響：首先是由于土洞的形成和發(fā)展使得周圍土體受到的初始?jí)簯?yīng)力得到釋放，以土洞為中心應(yīng)力重分布。再者由于地下水滲流與重力場(chǎng)共同作用，在土洞頂部形成一個(gè)呈漏斗形的應(yīng)力松弛帶，由內(nèi)向外變化逐漸減小。土洞頂部附近土體由高應(yīng)力轉(zhuǎn)變?yōu)榈蛻?yīng)力或零應(yīng)力，導(dǎo)致土體內(nèi)部發(fā)生剪切破壞，最終塌陷發(fā)生。

圖16 模型1(R=2.5 m)應(yīng)力云圖Fig.16 Stress nephogram in model 1(R=2.5 m)

4 結(jié)論

(1)該地區(qū)河床產(chǎn)生巖溶塌陷主要由于露采礦山大量抽排水使河水與采坑之間水力梯度增大，河水下滲產(chǎn)生強(qiáng)烈潛蝕效應(yīng)，使土洞形成并逐步向上擴(kuò)展，當(dāng)土洞頂板不能支撐自重時(shí)，發(fā)生地面塌陷。

(2)理論計(jì)算表明巖溶塌陷穩(wěn)定系數(shù)Ks主要受土層厚度、土體力學(xué)性質(zhì)、土洞高度、巖溶水下降幅度等影響。計(jì)算得出研究區(qū)河床塌陷臨界土洞高度為2.1 m，與實(shí)際情況較符。

(3)數(shù)值模擬結(jié)果表明該地區(qū)河床塌陷的臨界土洞高度為2～2.5 m，小于該地區(qū)階地巖溶地面塌陷的臨界土洞高度3.5 m。由于缺少上層滯水的滲透作用，階地土洞的擴(kuò)展速度也較慢，能夠解釋河床塌陷發(fā)生后1～2年內(nèi)周圍200 m內(nèi)陸續(xù)出現(xiàn)地面塌陷的原因。

(4)滲流模擬結(jié)果證明了巖溶水下降后潛蝕效應(yīng)明顯，區(qū)內(nèi)巖溶塌陷產(chǎn)生的根本原因是水動(dòng)力條件發(fā)生變化。故在巖溶地區(qū)抽排地下水時(shí)，應(yīng)在控制最大降深，還需要控制抽排水速度，防止出現(xiàn)土體滲透破壞。

(5)由于河床巖溶塌陷可作為大規(guī)模階地塌陷發(fā)生的前兆，故對(duì)類似區(qū)域的河床、水塘發(fā)生的塌陷應(yīng)引起重視。在河床塌陷發(fā)生后應(yīng)及時(shí)對(duì)塌陷坑進(jìn)行注漿或回填，防止巖溶管道被疏通；宜立即停止周圍大規(guī)模抽排水活動(dòng)；需對(duì)周圍地下水變化進(jìn)行長(zhǎng)期監(jiān)測(cè)，地下水在巖土交界面波動(dòng)時(shí)期尤其注意防范地面塌陷。