基于小尺度數(shù)值模型的海底滑坡運動敏感性分析

修宗祥，劉樂軍，解秋紅，李西雙，胡光海，李家鋼，趙強

（1.國家海洋局第一海洋研究所，山東　青島　266061；2.海岸和近海工程國家重點實驗室，遼寧　大連　116024；3.中海油研究總院，北京　100027）

基于小尺度數(shù)值模型的海底滑坡運動敏感性分析

修宗祥1，2，劉樂軍1，解秋紅1，李西雙1，胡光海1，李家鋼3，趙強1

（1.國家海洋局第一海洋研究所，山東青島266061；2.海岸和近海工程國家重點實驗室，遼寧大連116024；3.中海油研究總院，北京100027）

開展海底滑坡運動特性研究是深水陸坡區(qū)滑坡地質災害認識與防治的基礎，建立了基于非牛頓流體歐拉-歐拉兩相流理論的小尺度海底滑坡數(shù)值模型。在與實驗數(shù)據(jù)和BING程序結果對比驗證的基礎上，模擬分析海底滑坡的一般運動規(guī)律及特性，并同無水條件下的滑坡模擬結果進行了對比。結果表明：環(huán)境水的存在可引發(fā)“滑水”現(xiàn)象，延長滑坡運動時間，增加運動距離，但端部最大峰值速度相對無水條件時較小；滑坡體物質組成、地形坡度、初始速度、初始厚度等因素，對最終的運動距離有較大的影響；滑坡體在運動過程中因擾動、混水而導致的屈服強度和粘滯系數(shù)的不斷降低是海底滑坡長距離運動的主要原因。

海底滑坡；歐拉-歐拉兩相流模型；Herschel-Bulkley模型；滑水；軟化

深水陸坡區(qū)海底滑坡產生的高速滑動，可能給海底管道等設施造成破壞性影響 （Locat et al，2002；Nadim et al，2005；Mosher et al，2010）。由于陸坡區(qū)水深較大，海底地貌崎嶇復雜，海底管道一旦破壞，不但維修難度和費用較大，還會帶來嚴重的環(huán)境問題。因此，海底管道設計與路由選擇必須考慮海底滑坡地質災害的影響（Parker et al，2008；Randolph et al，2010；Yuan et al，2015）。而開展海底滑坡運動規(guī)律特性研究，分析土的組成、地形坡度、初始速度、初始厚度等因素與運動速度、距離的關系，有利于提高對陸坡區(qū)海底滑坡地質災害認識。

限于觀測條件限制，室內物理實驗與數(shù)值模擬成為研究海底滑坡運動特性的主要手段（Locat et al，2002；Wright et al，2007）。而數(shù)值模擬技術作為一種快速有效的方法，近年來已被廣泛采用。Imran等（2001a）采用Herschel-Bulkley和雙線性流變模型模擬海底滑坡運動，該模型未能考慮“滑水”現(xiàn)象與Mohrig等（1998）的部分實驗結果不符。Harbitz等（2003）基于潤滑理論提出了穩(wěn)定“滑水”狀態(tài)下的一維模型，模型忽略了前部、尾部及頂部的動水壓作用，但研究表明動水壓作用是十分重要的 （Wright et al，2007；Ilstad et al，2004a）。De Blasio等（2005）在Imran等人的基礎上將滑坡滑動過程分為初始流動、產生水楔、“滑水”、“滑水”停止的四階段模型，該模型與實驗符合較好，但臨界Froude數(shù)、初始水楔形狀、滑坡土體與海床間的粗糙度高度參數(shù)等均為人為假定。多相流理論近年來被引入到海底滑坡模擬中（Gauer et al，2006；Xiu et al，2015），室內實驗對比證明該方法可以有效模擬“滑水”現(xiàn)象，為海底滑坡機理研究提供了新的途徑。

目前國內海底滑坡相關研究主要集中在前期的參數(shù)識別、觸發(fā)機制與穩(wěn)定性研究方面（劉保華等，2005；吳時國等，2008；劉杜娟等，2010；胡光海，2011；彭俊等，2014；劉樂軍等，2014），針對滑后的運動特性研究相對較少。考慮實際海底滑坡的復雜性和多樣性，本文通過小尺度歐拉-歐拉兩相流數(shù)值模型探討海底滑坡運動的一般特性，并分析了多種工況參數(shù)對滑坡滑動影響的敏感性。研究結果可為我國深水陸坡海底滑坡地質災害認識與防治提供參考。

1　歐拉-歐拉兩相流模型

兩相流模型中環(huán)境水采用牛頓流體模擬，滑坡體可以看作粘塑性流體（Imran et al.，2001a），其流變特性可以由Herschel-Bulkley模型表達

其中，τC為屈服應力，γ˙為剪切速率，K為粘度系數(shù)，n為流動指數(shù)。滑坡體與水的相互作用可采用歐拉-歐拉兩相流模型模擬。忽略滑坡運動運動過程中的熱傳導作用，流體基本控制方程如下

質量守恒方程：

式中，▽為梯度算子，ρd，ρw，fd，fw，Vd，Vw分別為滑坡體與水的密度、體積分數(shù)以及速度向量，其中fw=1-fd。滑坡體與水的動量平衡方程如下：

其中，P為壓力，g為重力加速度，μd、μw分別為滑坡體和水的有效粘度，Dwd與Ddw為滑坡體與水在交界面處的相互作用力。滑坡體采用層流模型，環(huán)境水采用k-ε湍流模型模擬。

2　兩相流模型驗證

為了驗證兩相流模型的有效性，分別與文獻（Imran et al，2001 b）中的實驗數(shù)據(jù)及BING程序結果（其中筆者對BING模型進行簡單修改，以模擬初始形狀為矩形的情況）進行了對比。實驗中將膨潤土合成泥漿從容器倒入3.44°傾角的矩形槽中。泥漿初始長度L=1.8 m，初始厚度為H=0.3 m，泥漿密度ρd=1 073 kg/m3，屈服強度τy=42.5 Pa，參考應變速率γr=193.2 s-1。限于尺寸，僅給出兩相流模型局部示意圖，見圖1。其中，泥漿采用Bingham流體，環(huán)境空氣采用牛頓流體，流體單元總數(shù)164 360個，計算步長0.002 s。模型頂部為壓力開放邊界，前后兩側采用對稱邊界，左右兩側采用自由滑動邊界。BING程序單元數(shù)取21，時間步長取1×10-6。計算結果可以看出（圖2），兩相流模型模擬結果，與BING模型結果接近，與實驗值誤差約為6%。考慮到數(shù)值模型因在左側設置墻邊界以阻止泥漿向后流動，可能會導致計算結果稍大。認為兩相流模型的結果是可靠的。

由于目前公開文獻中尚未有提供水環(huán)境下數(shù)值模擬所需的全部實驗初始參數(shù)。這里僅對兩相流模型與BING模型結果進行比較。BING忽略了環(huán)境水的動力作用，僅考慮靜水壓的影響。通過對比可體現(xiàn)出兩相流模型的優(yōu)越性。根據(jù)文獻（Ilstad et al，2004 b）數(shù)據(jù)，設泥漿為Bingham流體，環(huán)境水采用牛頓流體，泥漿密度ρd=1 600 kg/m3，屈服強度τy=60 Pa，粘度系數(shù)K=0.035 Pas。斜坡角度6°，泥漿初始形狀采用矩形條帶，長度L=2 m，初始厚度分別為0.3 m和0.5 m兩種情況。

圖1　局部兩相流模型

圖2　最終形態(tài)對比

圖3為兩種方式下泥漿端部到達不同位置時對應的速度值。可以看出，兩相流模型所得到的峰值速度小于BING模型模擬值，這主要是由于BING模型沒有考慮水的阻力，而兩相流模型考慮了泥漿與水的動態(tài)相互作用。同時由圖4和圖5可以看出，雖然存在水的阻力作用，但兩相流模型的最終滑動距離卻并未減小，這主要是由于泥漿滑動過程中發(fā)生“滑水”效應所致。圖3中兩相流模型的泥漿端部速度變化趨勢也與之相符，即在前半部分保持在一定的速度，之后慢慢衰減。圖4中紅圈位置拉伸作用明顯，也主要是由于泥漿端部因滑水導致其速度相對后部較快產生的。Morrig等根據(jù)實驗曾給出的滑水啟動的臨界速度條件

圖3　端部速度變化對比

圖4　兩相流模型與BING最終形態(tài)對比

圖5　水槽實驗設計圖

其中，Vc為滑水啟動時泥漿前端的臨界速度；h為端部的厚度；φ為斜坡角度；k為密度傅汝德數(shù)，約為0.3～0.4。若取保守值k=1，則根據(jù)公式（6），泥漿端部厚度范圍0.15～0.3之間的滑水臨界速度區(qū)間為0.94～1.32，這也與兩相流模型的速度曲線相符合。由上述對比可以看出兩相流模型適合計算水與土體的相互作用。

3　滑坡運動特性及參數(shù)敏感性分析

3.1模型設計

參考Ilstad等（2004b）實驗（圖5），設計數(shù)值模型如圖6。假定初始滑坡體為矩形并具有一定初始速度。由于實際海底滑坡初始形狀接近矩形（L'Heureux et al，2013），且地震等因素觸發(fā)的滑坡會具有一定的初始動能，所以這種假設具有一定的實際意義。兩相流數(shù)值模型長為16 m，同時設計了1°、3°和6°3種坡度以考慮地形坡度的影響。為了對比水下滑坡與陸地滑坡的不同，同時對無水條件下的運動進行了模擬。為估算網(wǎng)格數(shù)量對模擬結果的影響，計算了85 000和170 000兩種單元數(shù)目下的最終運動距離，兩者基本相同，因此本文計算均選擇85 000網(wǎng)格數(shù)進行計算。

圖6　模型圖

數(shù)值模型中模擬滑坡體的力學參數(shù)，根據(jù)實驗中的4種工況選取。不同粘土含量下的工況定性代表了不同物質組成的滑坡體。其中，模擬滑坡的泥漿由高嶺土（ρ=2 650 kg/m3）、石英砂（ρ=2 750 kg/m3）以及水按質量分數(shù)比例混合而成，混合后密度ρ= 1 690 kg/m3。其流變特性由流變儀測出，表1為4種不同粘土含量下泥漿的Herschel-Bulkley模型擬和參數(shù)，計算過程中將公式（1）轉化為表觀粘度表達式，即μ=τC/Kγ˙+Kγ˙n-1。同時，針對每種材料工況分別進行了不同初始厚度與初始速度下的模擬對比，以討論不同因素對海底滑坡滑動過程的影響。初始厚度根據(jù)實驗觀測值分別選取0.1 m、0.2 m、0.3 m、0.4 m、0.5 m。初始速度分別選取了0.5 m/s、1.0 m/s、1.38 m/s、1.5 m/s、2 m/s，中間工況1.38 m/s為根據(jù)實驗漿體總體積、出口面積以及灌入時間求出的平均值。

表1　不同粘土含量泥漿參數(shù)

3.2參數(shù)變化的影響

圖7為初始模型高為0.5m，初始速度為1.0 m/s時，針對4種不同粘土含量漿體最大滑動距離。可以看出，粘土含量對滑動距離影響明顯，隨著粘土含量的增加，最大運移距離顯著降低，且接近指數(shù)衰減。說明不穩(wěn)定海床的物質成分組成，決定著海床土的力學性質，是影響海底滑坡滑動距離的關鍵因素。不考慮水動力條件下的模擬值，明顯小于考慮水動力時的模擬值，粘度含量越低、地形坡度越大，兩者差異越明顯。

圖7　最大運移距離隨粘土含量變化曲線

圖8為粘土含量為28.7%、初始高度為0.5 m時，不同初始速度下滑坡的最大運移距離。可以看出考慮水動力影響的最大運移距離仍然明顯高于不考慮水動力時的運移距離。隨著初始速度的增大，最大運移距離也相應較大。表明能夠給滑坡體帶來一定初始動能的觸發(fā)機制（如地震等）可能會造成較大的滑動距離。

圖8　最大距離隨初始速度變化曲線（粘土含量28.7%）

圖9為初始速度為1.38 m/s時，不同初始厚度下對應的最大運動距離。可以看出，不管是有水還是無水工況下，最大滑動距離都隨著初始厚度的增加而增大，接近于冪增長函數(shù)。說明相同條件下，滑坡厚度越大，滑動距離往往越大。相比之下，考慮水動力作用時，厚度變化對最大運移距離影響較不考慮水動力作用時更敏感。

圖9　最大距離隨初始厚度變化曲線（粘土含量28.7%，初始速度1.38 m/s）

3.3水動力作用

3.3.1水動力作用

由前文計算結果可以看出，水動力作用能對滑坡運動距離產生較大的影響。水對滑坡體主要影響方式有兩種：形成水楔而啟動“滑水”現(xiàn)象與混入降低強度。“滑水現(xiàn)象”主要是由于滑坡土體快速運動時，在其底部形成類似水楔的水層。圖10為28.7%粘土含量，在6°斜坡、初始速度1 m/s工況下，1.5 s時的體積分數(shù)云圖，此時最大端部速度為1.02 m/s，大于公式（6）計算的滑水臨界速度。其他工況計算結果也表明，滑動過程中出現(xiàn)過滑水現(xiàn)象。這也是海底滑坡土體可以在較小地形角度下滑動較大距離的原因之一。

圖10　T=1.5 s時體積分數(shù)云圖（粘土含量28.7%，初始速度1 m/s）

3.3.2混水軟化效應

滑坡體在運動過程除了因擾動強度降低之外，環(huán)境水的混入也導致了強度降低、流動性加強。混水主要通過兩種方式，一種是滑水時滑坡體端部下部水楔的不斷混入，另外滑坡體端部受到拉裂等作用也使水不斷混入。目前對于上述軟化過程（真實狀態(tài)下海底滑坡從啟動到最終的停止，更是一個非常復雜的變化過程），還很難進行準確的力學模型描述。為了分析強度降低對海底滑坡運動的影響，這里進行理想化的假設，只考慮屈服抗剪強度的變化對運動距離的影響。以粘土含量為28.7%的工況為例，其初始屈服抗剪強度約為108.7 Pa，對其逐步降低進行求解。圖11為對應不同的最大運移距離。可以看出，強度的降低能夠大幅提升運動距離。而真實狀態(tài)下海底滑坡土體的強度降低，是一個不斷累計發(fā)生的過程，而強度降低導致的流動性增強，又會進一步滿足滑坡體滑水所需的速度條件，這也是海底滑坡大距離運移的主要原因。另外，當弱化后的屈服抗剪強度分別等于粘土含量為15%和20%的屈服抗剪強度時，弱化后的最大運移距離仍然小于他們對應的最大運移距離，粘土含量越低，差值越大。上述表明，采用塑性流體模擬海底滑坡運動時，要同時考慮其屈服強度和粘滯系數(shù)降低的影響，特別是滑動后期滑坡體表現(xiàn)為較明顯的流體特性階段。

圖11　不同屈服強度下泥漿的最大運移距離

4　結論

本文主要基于小尺度的數(shù)值模擬討論海底滑坡的一般性規(guī)律。歐拉-歐拉兩相流模型能夠通過計算滑坡體與環(huán)境水的相互作用實現(xiàn)對“滑水”現(xiàn)象的捕捉，可較好的模擬海底滑坡的運動情況。地形角度、材料組成、初始動能以及滑坡規(guī)模都對最終運動距離有較大影響。通過與無水條件下的模擬結果對比，小角度地形下海底滑坡的運動距離相對較大，同時海底滑坡在運動過程中因擾動和環(huán)境水混入而導致的強度降低能夠顯著提高其運動距離。

實際的海底滑坡過程是一個非常復雜的過程，目前歐拉-歐拉兩相流模型還很難全面考慮海底滑坡在運動過程中因擾動和混水而導致的強度降低情況，同時也不能反映滑坡體內砂的單獨沉積情況。原因除了模型本身的理論限制外，對于上述現(xiàn)象的理論、實驗研究也有待于進一步開展。

De Blasio F V,Elverhoi A,Issler D,et al,2005.On the dynamics of subaqueous clay rich gravity mass flows-the giant Storegga slide, Norway.Marine and Petroleum Geology,22(1):179-186.

De Blasio F V,Engvik L,Harbitz C B,et al,2004.Hydroplaning and submarine debris flows.Journal of Geophysical Research,109, C01002.

Gauer P,Elverh?i A,Issler D,et al,2006.On numerical simulations of subaqueous slides:Back-calculations of laboratory experiments. Norwegian Journal of Geology,86:295-300.

Harbitz C B,Parker G,Elverh?i A,et al,2003.Hydroplaning of subaqueous debris flows and glide blocks analytical solutions and discussion.Journal of Geophysical Research,108(B7):107-137.

Ilstad T,D E Blasio F V,Elverh?i A,et al,2004a.On the frontal dynamics and morphology of submarine debris flows.Marine Geology,213(1-4):481-497.

Ilstad T,Elverhbi A,Issler D,et al,2004b.Subaqueous debris flow behaviour and its dependence on the sand/clay ratio:a laboratory study using particle tracking.Marine Geology,213:415-438.

Imran J,Harff P,Parker G,2001a.A numerical model of submarine debris flow with graphical user interface.Computers Geoscience, 227:721-733.

Imran J,Parker G,Locat J,et al,2001b.1D numerical model of muddy subaqueous and subaerial debris flows.Journal of Hydraulic Engineering,127(11):959-968.

L'Heureux J S,Vanneste M,Rise L,et al,2013.Stability mobility and failure mechanism for landslides at the upper continental slope off Vester?len,Norway,Marine Geology,346:192-207.

Locat J,Lee H J,2002.Submarine landslides:Advances and challenges. Canadian Geotechnical Journal,39(1):193-212.

Mohrig D,Whipple K X,Hondzo M,et al,1998.Hydroplaning of subaqueous debris flows.Geological Society of America Bulletin, 110(3):387-394.

Mosher D C,Moscardelli L,Shipp R C,et al,2010.Submarine mass movements and their consequences.Mosher D C,Moscardelli L, Shipp R C,et al,eds,Submarine mass movements and their consequences,Springer,Berlin,1-8.

Nadim F,Locat J,2005.Risk assessment for submarine slides.Hunger O, Fell R,Couture R,et al,eds,International conference for landslide risk,Aa Balkema:Vancouver,321-334.

Parker E J,Traverso C,Moore R,et al,2008.Evaluation of landslide impact on deepwater submarine pipelines.OTC.Waves of change. Curran Associates,New York,OTC19459.

Randolph M,Seo D,White D J,2010.Parametric Solutions for Slide Impact on Pipelines,Journal of Geotechnical and Geoenvironmental Engineering,136(7):940-949.

Wright S G,Hu Hongrui.2007.Risk assessment for submarine slope stability-hydroplaning.Final project report prepared for the Minerals Management Service under the MMS/OTRC cooperative research agreement 1435-01-04-CA-35515 task order 39323. Austin:The University of Texas.

Xiu Z X,Liu L J,Xie Q H,et al,2015.Runout prediction and dynamic characteristic analysis of potential submarine landslide in Liwan 3-1 gas field.Acta oceanologica sinica,34(7):116-122.

Yuan F,Li L L,Guo Z,et al,2015.Landslide impact on submarine pipelines:analytical and numerical analysis,Journal of Engineering Mechanics,ASCE,141(2):04014109.

胡光海，2011.東海陸坡海底滑坡識別及致滑因素影響研究.青島：中國海洋大學.

劉保華，李西雙，趙月霞，等，2005.沖繩海槽西部陸坡碎屑沉積物的搬運方式：滑塌和重力流.海洋與湖沼，36（1）：1-9.

劉杜娟，潘國富，葉銀燦，2010.東海陸架典型海洋災害地質因素及其聲反射特征.海洋通報，29（6）：664-668.

劉樂軍，傅命佐，李家鋼，等，2014.荔灣3-1氣田海底管道深水段地質災害特征.海洋科學進展，32（2）：162-174.

彭俊，陳沈良，陳一強，等，2014.黃河三角洲侵蝕性岸段水下岸坡地質災害及其空間分布.海洋通報，33（1）：1-6.

吳時國，陳珊珊，王志君，等，2008.大陸邊緣深水區(qū)海底滑坡及其不穩(wěn)定性風險評估.現(xiàn)代地質，22（3）：430-437.

（本文編輯：袁澤軼）

Sensitivity analysis of submarine landslide mass movement based on the small-scale numerical model

XIU Zong-xiang1,2,LIU Le-jun1,XIE Qiu-hong1,LI Xi-shuang1, HU Guang-hai1,LI Jia-gang3,ZHAO Qiang1
（1.First Institute of Oceanography,SOA,Qingdao 266061,China;2.State Key Laboratory of Coastal and OffshoreEngineering,Dalian 116024,China;3.CNOOC Research Institute,Beijing,100027,China）

Study of submarine landslide movement characteristics is the basis for the understanding and prevention of the landslide geological hazard in the continental slope area.The small-scale numerical model of submarine landslide is established based on the Non-Newtonian Eulerian-Eulerian two-phase model.The validity of the established model is verified through the comparison with experiment data and the BING simulated result.The general movement behaviors of submarine landslide are then studied,and the simulated results are also compared with those without water.The result shows that the ambient water can lead the occurrence of hydroplaning phenomenon,which mainly can extend the running time of submarine landslide,and make a larger run-out distance.But the calculated maximum front peak velocity of submarine landslide is smaller than that without water.Soil composition,topographic slope,initial velocity and initial thickness play important roles in the submarine landslide run-out distance.The continuing decrease of yield strength and viscosity resulted by the soil disturbance and the incorporation of water during the movement process are the main reason for the larger run-out distance of submarine landslide.

submarine landslide;Eulerian-Eulerian two-phase model;Herschel-Bulkley model;hydroplaning;soften

P694

A

1001-6932（2016）04-0380-06

10.11840/j.issn.1001-6392.2016.04.003

2015-04-02；

2015-08-27

國家自然科學基金青年基金項目（41206058）；國家科技重大專項子課題（2011ZX05056-001-02）；中央級公益性科研院所基本科研業(yè)務費專項資金項目（2013a25）；大連理工大學海岸和近海工程國家重點實驗室開放基金項目（LP1514）

修宗祥（1982-），男，博士，副研究員，主要從事海洋工程地質災害風險評價研究。電子郵箱：xiuzongxiang@163.com

解秋紅（1982-），女，博士，助理研究員。電子郵箱：xqh@fio.org.cn。