目前已知的海底滑坡多数由火山活动或者地震等自然地质因素所引发。但近几年的研究表明，赋存于海底低温高压环境下的天然气水合物由于其对于外界环境的敏感性，在外界环境变化时容易产生分解，分解过程中产生的大量水和天然气导致地层的孔隙压力上升，同时在此过程中水合物储层的黏聚力和土体颗粒的胶结程度也会降低，两者的同时作用也是导致斜坡稳定性降低的重要因素，严重时可引起海底滑坡^[1-2]。

国内外学者对合物分解诱发的海底滑坡有诸多研究。SULTANA^[3]通过软件在考虑温度、压力、孔隙水化学和平均孔径分布情况下，对地层中天然气水合物的热力学化学平衡进行了理论研究，并提出Storegga滑坡的产生根源在于水合物的分解。杨晓云^[4]讨论了海底边坡的坡角、水合物的饱和度以及分解速率对于滑坡产生的影响，并且建立了模拟水合物开采的边坡模型对滑坡产生的各项因素进行了讨论。刘峰等^[5]通过实验室三轴试验结合数值模拟，对水合物分解产生的孔隙压力以及分解过程中土体强度的变化进行了模拟，并结合相关因素对于海底边坡的稳定性影响进行了讨论。周丹^[6]分析认为水合物分解时导致的储层强度减弱是产生海底滑坡的重要因素，并且上覆土体越厚，海底边坡的稳定性越高。马云等^[7]采用ABAQUS软件建立了包含水合物层的海底边坡模型，利用数值模拟评估了边坡稳定性，给出了潜在滑坡的塑性贯通区，以及水合物分解过程中土体强度变化对于边坡稳定性的影响。GUPTA等^[8]建立了热、化学、力学多场耦合模型，模拟了水合物的分解过程，讨论了水合物的分解过程对于地层形变的影响。LU等^[9]对各因素控制的的地层变形和应力分布进行了数值模拟，并对边坡产生的位移进行了讨论。从国外发展现状可以看到，专家学者主要对于水合物分解诱发海底滑坡的产生机理、边坡失稳的影响因素等进行了理论研究，但是缺少与实际工区地质情况的结合，且一些地区数值模拟过程中的力学等参数并非实际地质参数，因此，模拟得到结果与实际情况存在差别。

针对当前研究的不足之处，本文利用神狐海域已有的地震成果数据，通过Gassmann方程计算出储层的孔隙度和饱和度^[10]，加之反演得到的其他岩石物理参数，基于海底斜坡的地质力学模型和强度折减法，揭示了应力的变化过程，并利用数值模拟，得到水合物分解程度对海底边坡稳定性的影响，实现水合物分解诱发海底滑坡致灾风险的初步评价。

为直观表达出海底滑坡致灾的风险程度，本文采用了安全系数法来评价海底边坡的稳定性，根据安全系数是否<1（即抗剪切力是否小于剪切力）来判断边坡是否失稳产生破坏。采用安全系数法进行安全风险判别，可以直观、精确地评估出工区面临的灾害风险，判断边坡是否失稳产生破坏。安全系数的定义如式（1）

式中：$\phi $为土体的内摩擦角，（°）；

c为土体的黏聚力，Pa；

τ为有效应力法向分量，Pa；

$ \mathrm{\alpha } $为有效应力剪切分量，Pa。

忽略其他环境因素变化对于海底斜坡影响的情况下，海底边坡所受的应力主要为静水压力、海底土体的自身重力、以及水合物分解过程中产生的孔隙压力，据此建立了海底斜坡力学模型（图1），用于计算潜在滑动面上有效应力的法向分量和剪切分量。

图  1  简化海底斜坡力学模型

Figure 1.  Simplified submarine slope mechanical model

式中：${\alpha _{\rm{w}}}$为沉积物自身重力，N；

${\gamma _{\rm{s}}}$为单位体积沉积物重量，N/m³；

${\alpha _{{\rm{hp}}}}$为静水压力，Pa；

${\gamma _{\rm{w}}}$为单位体积水的重量，N/m³；

${\alpha _{\rm{p}}}$为地层本身孔隙水压力，Pa；

${\alpha _{{\rm{ep}}}}$为超孔隙水压力，Pa；

$\alpha $有效应力的剪切分量，Pa；

$\beta $为海底边坡的剪切面倾角，（°）；

$\tau $为有效应力的法向分量，Pa。

最终通过三轴实验获得的抗剪切强度参数和由孔隙压力平衡方程计算得到的孔隙压力^[5]，联合实现安全系数的求取。

南海神狐海域主要为泥质沉积岩。由于该类岩石的孔隙度和渗透率较小，因此对于外力影响的响应比较明显，内部压力变化也不容易释放到外部。在这种地质条件下，地层中孔隙压力的变化可以等效为地层有效应力的变化，运用天然气水合物分解引起孔隙压力变化的模型来计算有效应力变化^[11-12]，最终求取地层的有效应力。

假设1 m³的天然气水合物完全分解后产生164 m³的气体和0.87 m³的水，并且储层中的固体相和水相处于不可压缩状态，气体体积和压力变化符合波义耳定律。则孔隙压力变化可以表示为公式（4），

式中：${\sigma '}$为水合物分解后产生的孔隙压力，在计算中等效为地层的有效应力u；

$e$为模型的孔隙比；

$\varphi $为孔隙度；

${S_{\rm{w}}}$为含水饱和度；

${C_{\rm{s}}}$为模型的膨胀指数；

$latm$为标准大气压（1.013$ \times $10⁵ Pa）；

${P_2}$为水合物处于平衡状态时的压力，Pa；

${T_2}$为水合物处于平衡状态时的温度，℃；

${P_1}$为水合物未分解时的压力，Pa；

${V_{{\rm{GH1}}}}$为水合物未分解时体积，m³。

由于泥岩的渗透率较低，可以假定水合物的分解过程中模型地层的总体积没有改变，游离气体积不变，公式（4）进一步表示为，

在计算孔隙度和饱和度时，由于速度或密度的变化同时受到孔隙度和饱和度的影响，采用Gassmann方程来求取孔隙度和饱和度的计算精度更高。因此，本文Gassmann方程计算了储层的孔隙度$\varphi $和含水饱和度${S_{\rm{w}}}$^[13]，计算过程如公式（6）所示，

式中：${V_{\rm{P}}}$，${V_{\rm{s}}}$为纵横波速度，m/s，

$\mu $为切变模量，Pa；

$k$为体积模量，Pa；

$\rho $为密度，kg/m³；

$\sigma $为泊松比。

角标f为孔隙中的流体，b为岩石，s为岩石基质，w为水，h为油气类流体。

M为应力与应变的比值，用以表示土体的压缩性高低，由式（7）表示，

所以，

其中膨胀指数${C_{\rm{s}}}$无法直接获取，为了得到膨胀指数，引入了压缩因子的概念${C_{\rm{c}}}$，可以表示为，

上式利用了塑性系数$PI$和比重${G_{\rm{s}}}$表示压缩因子。WROTH等^[14]之前的研究得到：

所以公式（8）可以化简为，

将公式（11）进一步化简，

求解可得，

式中：${C_1}$与${C_2}$为常数，${C_1}$可以通过代入空率、含水饱和度，平衡温度、平衡压力以及膨胀因子求取，${C_2}$可通过模型的边界条件以及模型的初始条件求取，求解可得地层的孔隙压力${\sigma '}$，根据上述假设该孔隙压力可以等效为地层的有效应力u。

根据建立的海底斜坡地质力学模型，安全系数计算公式如式（14）所示，

式中：$C$为土体黏聚力，Pa；

${\gamma _{\rm{w}}}$单位体积水的重量，N/m³；

${\gamma _{\rm{s}}}$为单位体积沉积物的重量，N/m³；

$u$为有效应力，Pa；

$H$为海面到边坡水深，m；

$h$为潜在滑动面的深度，m；

$\phi $为内摩擦角，（°）；

$\beta $为海底边坡的坡角，（°）。

强度折减法的定义为：保持外载荷不变的条件下，边坡土体所能提供的最大抗剪强度与实际受到剪应力之比。边坡在破坏的临界状态时${K^{{\rm{cr}}}}$=1，可得到原始边坡的的安全系数F_s为：

式中：${C^0}$为边坡初始黏聚力；

${\phi ^{\rm{0}}}$初始内摩擦角，（°）

$K$为折减系数；

${C^{{\rm{cr}}}}$为边坡实际黏聚力，Pa；

${\phi ^{{\rm{cr}}}}$为边坡实际内摩擦角，（°）。

结合南海神狐海域实际的地球物理资料，对横向1 500 m、纵向1 200 m范围内、水深1 000 m的海底陆坡深层地层信息进行敏感参数叠前反演，获得不同地层的密度、纵横波速度等，计算了不同地层的孔隙度和含水饱和度，获得了沉积物层的有效应力。

地层弹性参数^[15-16]联合三轴试验获得的抗剪参数（表1）、与经过时深转换后的地震剖面（图2），采取Flac3D软件建立了海底陆坡地质模型（图3），以海面为零点，研究区水合物埋深为190～220 m左右，用于后续诱发海底滑坡影响因素的力学分析与致灾风险评价。

图  2  海底陆坡地震剖面

Figure 2.  Submarine slope seismic profile

图  3  建立的海底陆坡地质模型

Figure 3.  Submarine continental slope model

参照上述地质条件，整个模型设置大小为 1500 m×1 500 m×2 200 m（长×宽×高），在数值模拟中进行网格剖分，模型共划分10 000 个单元。

本文数值模拟采用的边界条件为:

在x = 0 的位置在对 x 进行约束，其他方向不加约束；在 x = 1 500 的位置对x进行约束，其他方向不加约束；在 y 方向上采用固定约束；与x相同在z=0和2 200的位置对z进行约束，其他方向不进行约束；模型的底部边界进行固定约束。

利用Flac3D软件进行了数值模拟，获得了综合考虑初始水合物分解量、以及试采过程中的水合物分解总量等因素变化引起的水合物储层变化的力学响应剖面和对应的安全系数。

研究天然气水合物初始分解量的影响时，假设含水合物储层渗透率很低，则底层内部的压力向外部传递有限，分解产生的压力被约束在内部空间中。分别模拟了初始天然气水合物分解量为0%、10%、20%、30%以及50%情况下的有效压力和安全系数，获得了相应的安全系数云图（图4）以及有效应力分布图（图5）。

图  4  初始天然气水合物分解量对应的安全系数

Figure 4.  The safety factor corresponding to the initial gas hydrate decomposition amount

图  5  初始天然气水合物分解量对应有效应力分布

Figure 5.  Effective stress distribution corresponding to initial gas hydrate decomposition amount

当初始天然气分解量达到30%时（图4d），图中开始出现了安全系数<1的区域，所以当初始天然气分解量达到30%时为该模型出现破坏失稳的临界值。随着初始天然气水合物分解量的增加，最终完全分解后（图4e），在斜坡东、西两侧均出现了安全因子<1的区域。两处海底滑坡表明，随着天然气水合物的不断分解海底边坡的稳定性不断降低，即初始分解量越大，海底越不稳定。

初始水合物分解量为10%时，水合物顶界面处的滑动面上最大有效应力为1.96 MPa（图5a）；随着初始水合物分解量的增大到30%时，滑动面上的最大有效应力降低为1.39 MPa（图5b）；最后，初始水合物全部分解时，滑动面上的最大有效应力变为1.22 MPa（图5c）。有效应力分布表明分解产生的孔隙压力在垂直方向上与模型本身的外部荷载进行了抵消。通过分析可以发现，随着天然气水合物初始分解量的百分比增大，分解产生的孔隙压力也在不断上升，导致天然气水合物完全分解后的有效应力减小，因此，海底边坡的稳定性降低，这与图4的安全系数图正好吻合。

此外，对于天然气水合物分解总量诱发海底滑坡进行了研究。假设初始天然气水合物分解量为10%，计算了开采状态下天然气水合物总量分解达到20%、40%、60%、80%时，安全系数的变化情况，如图6所示。

图  6  水合物分解总量对应的安全系数分布

Figure 6.  The distribution of safety factors corresponding to the total amount of hydrate decomposition

随着开采的进行，水合物持续分解，总量分解到40%时出现了安全系数<1的失稳区域（图6c），此时达到该模型海底滑坡的临界。由图6整体对比分析出，不论是由人工开采还是天然构造运动造成的天然气水合物分解，其总量越大，导致海底斜坡稳定性越差，出现滑坡灾害的程度就越严重。天然气水合物填充着固体胶结物的孔隙增加了沉积物的抗剪切强度，随着水合物分解而不断产生游离气体和水，一方面增大了地层的孔隙压力，另一面降低了储层的胶结程度，从而使得水合物储层的抗剪切能力不断下降，导致边坡稳定性降低，进一步诱发了海底产生滑坡。

基于南海北部实际地球物理资料，获得了反映真实地层结构的力学参数，以强度折减法探讨了初始天然气水合物分解量、水合物分解总量等不同因素影响下诱发海底滑坡的致灾风险和致灾机理，结果表明：

（1）初始天然气水合物分解量越大，其初始孔隙压力越大，对应有效应力就越小，此时海底斜坡越不稳定。在初始水合物分解量达到30%后，稳定边坡达到失稳临界，随着分解量进一步增大，海底斜坡会出现大范围滑坡。

（2）在初始水合物分解量为10%的相同条件下，随着分解总量的不断增加，产生游离气的水合物储层其有限空间里积聚的孔压不断增大，同时降低了水合物储层土体颗粒的胶结程度，两者共同作用使得水合物储层的抗剪切能力降低，增加了诱发海底滑坡产生的可能性，试算总分解量达到40%时边坡已经达到了失稳临界状态。

水合物分解诱发海底滑坡的影响因素还有很多，综合分析不同因素诱发海底滑坡的致灾风险和致灾程度，为水合物实际试采潜在地质灾害的预防提供更加有效的指导和帮助。