初始土體含水率對泥石流方量影響的數值模擬

謝鑫鵬,楊 龍,張 明,焦偉之,李成文,謝駿錦

(中國地質大學(武漢)工程學院，湖北 武漢 430074)

我國特殊的地質構造和復雜的地理環境，導致各種地質災害頻發，其中以泥石流地質災害最為嚴重。泥石流地質災害不僅會對生態環境造成破壞，而且嚴重威脅著人民的生命和財產安全[1]。泥石流是在降雨條件下，于溝谷中形成的攜帶大量泥砂、石塊等固體物質的復雜混合流體[2-4]。目前，關于降雨觸發泥石流地質災害的研究較多。如：Wilson等[5]和Miller等[6]將降雨誘發的泥石流地質災害歸納為兩種類型，一種是持續時間較短且短時間內降雨強度較大誘發的泥石流，另一種是持續時間較長且長時間內降雨強度較小誘發的泥石流；Rudolph等[7]分析了前期土體含水率和降雨類型對微地貌變化的響應；Coe等[8]分析了泥石流啟動前后土體的含水量變化特征，認為徑流對沖溝中堆積物的沖蝕作用是泥石流啟動的主要原因;Berti等[9]對泥石流溝道堆積物的滲透性進行了研究，發現降雨條件下溝道堆積物的滲透性越差土體越易達到飽和，并得出水動力條件是泥石流啟動的主要原因之一;宋兵等[10]通過對魯家溝泥石流進行詳細的現場調查，指出魯家溝泥石流是在暴雨作用下順層斜坡失穩破壞，堰塞溝道并潰決形成的；毛碩等[11]研究認為野牛溝泥石流的啟動主要是因為強降雨條件下導致的雨水快速匯集效應造成的；蘇遠東等[12]通過對熱河溝泥石流的形成過程進行詳細分析，得出熱河溝泥石流的形成機理為“暴雨-潰決”型；王毅等[13]以毛家灣隧道口溝泥石流為研究對象，對該泥石流的運動形式與強降雨之間的關系進行了深入研究；王征強等[14]以后溝泥石流為研究對象，通過對其形成機制進行分析得出充足的水源是該泥石流形成的重要條件之一；楊敏[15]以小秦嶺金礦區礦渣泥石流為例，通過對其形成的主要影響因素進行分析，得出該泥石流為暴雨型水石流，其啟動方式主要包含“河谷起動型”和“溝道坡面滑塌—溝道堵塞—形成潰決”兩種。

近年來，隨著計算機技術的革新與發展，越來越多的學者采用數值模擬的方法來研究泥石流的啟動過程和危險性。如：O’Brien等[16]通過對不同類型的泥石流進行分析研究，建立了FLO2D模型，為泥石流啟動過程的數值模擬提供了依據；胡明鑒等[17]采用二維顆粒流程序(PFC2D)分析了降雨作用下松散碎屑物質啟動形成泥石流的過程及其與土體含水率的關系；劉鑫磊等[18]為了預測潰決型泥石流沖出量特征值與威脅范圍，基于FLO2D軟件對不同降雨頻率下牛圈溝泥石流的危險性進行了數值模擬研究；胡卸文等[19]利用基于有限體積法的CFX軟件，對江口溝泥石流危險區范圍和速度場分布情況進行了數值模擬研究；左正金等[20]利用PFC2D軟件對欒川縣柿樹溝泥石流的失穩啟動過程進行了數值模擬研究，結果表明強降雨是柿樹溝溝內松散堆積體啟動形成泥石流的觸發因素。

目前的研究主要關注降雨對泥石流啟動及運移過程的影響，而針對泥石流失穩方量與土體含水率之間關系的研究較少。因此，本文以新疆維吾爾自治區蘇約克河流域內典型的泥石流溝——喬庫而泥石流為研究對象，通過分析其工程地質條件，建立了離散元數值模型，并利用PFC2D軟件對該泥石流物源區堆積物的失穩方量與土體含水率之間的關系進行了數值模擬研究，以為泥石流地質災害的預警預報提供科學依據[21-22]。

1 喬庫而泥石流概況

如圖1所示，在新疆維吾爾自治區蘇約克河流域內19 km長的河谷中，野外調查發現發育小型(V<2萬方)、中型(2萬方≤V<小于20萬方)、大型(20萬方≤V<50萬方)和巨型(V≥50萬方)泥石流溝共計39處，且該區段發育3條斷層，巖性以砂巖、泥巖、玄武巖為主，沿斷層兩側巖體較為破碎，經風化作用坡腳處堆積大量的風化物，為泥石流提供了豐富的物源，每逢雨季該流域內泥石流地質災害頻發。為了保證當地人民的生命和財產安全，亟待對該流域內泥石流失穩方量與土體含水率的關系進行研究。

本文選取蘇約克河流域內典型的泥石流溝喬庫而泥石流為研究對象，該泥石流溝位于蘇約克河的西岸(見圖1中的102號泥石流)，主要威脅對象為河流漫灘上的牧場。喬庫而泥石流的巖性主要為第三系的紅色砂泥巖和白堊系的黑色砂巖、礫巖，該泥石流區被一逆沖斷層橫切(見圖1)，斷層下盤為紅色的砂泥巖，巖體表面較為破碎，受風化作用在其橫穿的斜坡坡腳處往往堆積有小規模的細粒崩坡積層；上盤為黑色的砂巖、礫巖，在斷層面附近巖體極為破碎，形成的崩坡積物是喬庫而泥石流的主要物源。

喬庫而泥石流溝的整個流域面積為4.42 km2，其物源區、流通區和堆積區如圖2所示。其中，物源區主要為臨近斷層面的上盤、溝谷南側破碎山體失穩堆積在坡腳處的滑坡和崩塌堆積物，堆積物的方量約為120萬方；流通區分布在斷層下盤的溝谷底部，長度約為1 500 m，兩側為紅色砂泥巖組成的溝壁；堆積區面積約為0.11 km2，堆積物呈扇形堆積，平均厚度約為4 m，體積約為44萬方。根據這三個區域得到的喬庫而泥石流典型的剖面如圖3所示。

圖1 新疆維吾爾自治區蘇約克河流域內泥石流分布圖

圖2 喬庫而泥石流的平面圖

圖3 喬庫而泥石流典型的剖面圖

2 泥石流物源區堆積物失穩方量與土體含水率之間關系的數值模擬研究

本次針對喬庫而泥石流物源區堆積物失穩方量與土體含水率之間關系的數值模擬研究的具體思路為：首先，根據遙感影像和現場調查，估算喬庫爾泥石流的全部物源量以及威脅牧場和牧民的最小泥石流沖出量；然后，利用PFC離散元數值方法反演喬庫而泥石流在相當于最小泥石流沖出量的土體失穩啟動和物源區堆積物土體全部啟動并正常運行至河谷兩種工況下，土體所對應數值模型中顆粒的細觀參數；隨后，利用雙軸壓縮試驗獲取上述細觀參數對應的土體宏觀剪切強度參數；最后，結合直剪試驗所得到的泥石流物源區堆積物土體剪切強度參數隨土體含水率變化的曲線，得到土體相應的含水率。

喬庫而泥石流物源區堆積物土體為黃褐色的粉質黏土，硬塑，結構較為密實，土體粒徑較小，土體表面分布塊石，塊石粒徑約5～10 cm。如圖2所示，在喬庫而泥石流物源區堆積物土體(B01)中取了1個樣品，利用樣品在含水率分別為13%(天然)、20%、27%和35%(飽和)條件下進行了直接剪切試驗，樣品每種含水率對應100 kPa、200 kPa、300 kPa和400 kPa 4種不同的法向荷載,最終獲取的泥石流物源區土體(B01)剪切強度參數隨含水率的變化曲線，見圖4。

圖4 喬庫而泥石流物源區土體剪切強度參數隨含水率 的變化曲線

2. 1 數值模型建立

根據如圖3所示的喬庫而泥石流典型剖面建立其數值模型。在數值模擬過程中，為了清晰地觀察泥石流的啟動和運移過程，將物源區的粉質黏土劃分為4個部分，分別用顆粒填充。考慮到計算速度與結果可靠性之間的平衡，本次模擬采用的顆粒最小半徑為0.25 m、最大半徑為1 m，顆粒總數為16 785個，采用高斯分布。同時，為了更好地模擬喬庫而泥石流物源區散體物質的特征，顆粒即球體之間的接觸模型選用接觸黏結模型。除了物源區的土體，其余的邊界都用墻(Wall)來代替。在不影響結果可靠性的前提下為了提高數值計算效率，最終建立的喬庫而泥石流數值模型，見圖5。

圖5 喬庫而泥石流數值模型

2. 2 土體物理力學參數確定

本次數值模擬采用墻(Wall)和球(Ball)的計算參數，如表1所示。墻的參數包括法向剛度、切向剛度、摩擦系數，均取常用值，在數值模擬中采用固定參數。球體用接觸黏結模型，模型中需要8個參數，但是多參數反演不僅會大大增加計算量而且會導致多解情況的加劇。為了盡量避免這一情況，首先利用工程地質類比法[23]確定密度、法向剛度、切向剛度、半徑系數等對泥石流動力學過程影響較小的參數；然后，通過數值反演獲取法向抗拉強度、切向抗拉強度、黏聚力和內摩擦角等對泥石流動力學過程影響較大的參數。

表1 喬庫而泥石流物源區土體的細觀參數

2. 3 模擬工況設置

根據遙感影像和現場調查，假設當泥石流堆積物約為1 m厚時，就會對牧場和牧民產生威脅；而喬庫而泥石流堆積區的面積為11萬方，因此估算當泥石流一次性沖出11萬方的沖積物時，就會對泥石流扇前方和側緣的牧場和牧民產生威脅。同時根據估算，喬庫而泥石流物源區的全部物源為53萬方。

本次數值模擬反演了兩種工況：①工況1下，喬庫而泥石流物源區至少有11萬方物質失穩，這種工況下反演出的土體剪切強度參數和含水率是泥石流啟動所需要的土體最低剪切強度參數和含水率；②工況2下，喬庫而泥石流物源區的物質全部啟動，形成的泥石流到達河谷，這種工況下反演出的土體剪切強度參數和含水率是泥石流全部失穩啟動所需要的土體剪切強度參數和含水率。

2. 4 模擬結果分析

2.4.1 喬庫而泥石流啟動和運移過程反演結果分析

經過多次反演計算，當球體的法向抗拉強度、切向抗拉強度、黏聚力和內摩擦角分別等于20.4 kPa、22.7 kPa、20 kPa和10°時，模擬所得到的喬庫而泥石流啟動和運移過程的結果與工況1要求最為符合。圖6為工況1下喬庫而泥石流啟動和運移過程的數值模擬結果，數值反演計算得到的該泥石流堆積體的分布特征如圖6中150×104steps所示。

由圖6可見，雖然該泥石流物源區整體失穩，但是最終大部分土體(深藍色、綠色和紅色部分)穩定下來，幾何形態沒有發生變化，泥石流前緣藍色部分土體發生解體并有較大的位移，這部分土體的方量約為10.1萬方，與估算的最小威脅泥石流的土體方量接近。

經過多次反演計算，當球體的法向抗拉強度、切向抗拉強度、黏聚力和內摩擦角分別為12.8 kPa、15.4 kPa、15 kPa 和8°時，模擬得到的喬庫而泥石流啟動過程的結果與工況2要求最為接近，即此時喬庫而泥石流的流動距離能達到河谷。圖7為工況2下喬庫而泥石流啟動和運移過程的數值模擬結果，直觀地顯示了該泥石流物源區土體從失穩、運移到動能耗盡并最終停積下來的全過程，數值反演計算得到的該泥石流堆積體的分布特征如圖7中150×104steps所示。

圖6 工況1下喬庫而泥石流啟動和運移過程的數值 模擬結果

圖7 工況2下喬庫而泥石流啟動和運移過程的數值 模擬結果

由圖7可見，由物源區運移到流通區和堆積區的土體方量約為48萬方，該泥石流的最大運動距離約為2 450 m，且部分物源運移至蘇約克河道內。

2.4.2 喬庫爾泥石流物源區土體雙軸壓縮試驗結果分析

繼續利用上述工況1和工況2反演所得到的球體細觀參數包括法向抗拉強度、切向抗拉強度、黏聚力、內摩擦角和表1中其他參數，通過建立雙軸壓縮試驗數值模型(見圖8)，對喬庫而泥石流物源區土體進行雙軸壓縮試驗，可得到兩種工況下喬庫而泥石流啟動時物源區土體的剪切強度參數[24]，其試驗結果見圖9和圖10。

圖8 雙軸壓縮試驗數值模型

圖9 工況1下喬庫而泥石流物源區土體雙軸壓縮試驗結果

本次雙軸壓縮試驗選用的顆粒尺寸和黏結模型(接觸黏結模型)與上述泥石流數值模擬相同。在進行雙軸壓縮試驗時，為了保證試樣具有足夠的顆粒并消除試樣的尺寸效應，將雙軸壓縮試驗的數值模型尺寸設定為240 m×480 m，顆粒的數量為8 920個(見圖8)。根據泥石流數值模擬得到的細觀參數，設置雙軸壓縮試驗顆粒的細觀參數，并分別在圍壓為50 kPa、100 kPa、150 kPa條件下對喬庫爾泥石流物源區土體進行雙軸壓縮試驗。

由圖9(a)、(b)、(c)可知，通過采用工況1中球體的細觀參數在圍壓分別為50 kPa、100 kPa、150 kPa條件下對喬庫而泥石流物源區土體進行雙軸壓縮試驗獲取的土體應力-應變曲線，可得到該土體的峰值剪切強度分別為164 kPa、246 kPa、304 kPa；由圖9(d)可知，根據采用3種圍壓及對應的峰值剪切強度繪制的土體摩爾破壞包絡線，可得到喬庫而泥石流物源區土體的黏聚力為19.2 kPa、內摩擦角為21.5°。

由室內直剪試驗得到的喬庫而泥石流物源區土體黏聚力、內摩擦角與土體含水率的關系曲線(見圖4)可知：當土體黏聚力為19.2 kPa時，對應的土體含水率為33%時，當土體內摩擦角為21.5°時，對應的土體含水率約為31%，說明當物源區土體含水率為31%～33%時，即土體快達到飽和狀態時(土體飽和含水率為35%)，喬庫而泥石流物源區才約有10.1萬方的土體失穩啟動并有一定的運動位移。

圖10 工況2下喬庫而泥石流物源區土體雙軸壓縮試驗結果

由圖10(a)、(b)、(c)可知，根據采用工況2中球體的細觀參數在圍壓分別為50 kPa、100 kPa、150 kPa條件下對喬庫而泥石流物源區土體進行雙軸壓縮試驗獲取的土體應力-應變曲線，可得到該土體的峰值剪切強度分別為115 kPa、170 kPa、230 kPa；由圖9(d)可知，根據采用3種圍壓及對應的峰值剪切強度繪制的土體摩爾破壞包絡線，可得到喬庫而泥石流物源區土體的黏聚力為12 kPa、內摩擦角為9°。

由室內直剪試驗得到的喬庫而泥石流物源區土體黏聚力、內摩擦角與土體含水率的關系曲線(見圖4)可知：當土體黏聚力為12 kPa、內摩擦角為9°時，該泥石流物源區土體含水率已經超過其飽和含水率，說明喬庫而泥石流物源區的堆積物要全部轉化為泥石流并達到蘇約克河河谷，土體含水率要超過其飽和含水率。

由于喬庫而泥石流物源區的地形較高較陡，除了物源區本身接受的雨水外，其余匯水區匯集的雨水很難滲入到物源區土體中，加上研究區年降雨量少，可以推測喬庫而泥石流物源區土體的含水率很難高于其飽和含水率。因此，本研究認為喬庫而泥石流物源區土體在降雨作用下僅能失穩，如果沒有其他因素導致這些土體在啟動和運移過程中剪切強度降低，它們將很難形成遠程運動的泥石流。故喬庫而泥石流的運移機理是值得進一步深入研究的問題。

3 結 論

本文通過對新疆烏恰縣托云鄉喬庫而泥石流分別處于最小泥石流沖出量和物源區土體全部失穩啟動兩種工況下土體的含水率特征進行了數值模擬分析，得到如下結論：

(1) 喬庫而泥石流達到最小沖出量時(10.1萬方)，物源區土體含水率為31%～33%，即土體快達到飽和狀態(土體飽和含水率為35%)，此時物源區土體的黏聚力為19.2 kPa、內摩擦角為21.5°。

(2) 喬庫而泥石流物源區的土體全部失穩啟動時，物源區土體含水率需要超過其飽和含水率，此時物源區土體的黏聚力和內摩擦角分別為12 kPa和9°。

(3) 喬庫而泥石流物源區土體在降雨作用下僅能失穩，如果沒有其他因素導致這些土體在啟動和運移過程中剪切強度降低，它們將很難形成遠程運動的泥石流。因此，有必要進一步深入研究喬庫而泥石流的運移機理。