疏干条件下隧洞涌水量及孔压分布解析研究

吴建; 周志芳; 庄超

doi:10.11779/CJGE202010015

疏干条件下隧洞涌水量及孔压分布解析研究 English Version

河海大学地球科学与工程学院,江苏　南京　211100

基金项目:

国家重大研究计划重点支持项目 91747204

详细信息

作者简介:
吴建(1994—),男,博士研究生,主要从事地下水数值模拟相关研究。E-mail：wujian1994@hhu.edu.cn

通讯作者:
周志芳, E-mail：zhouzf@hhu.edu.cn

中图分类号: U45
计量
- 文章访问数: 287
- HTML全文浏览量: 45
- PDF下载量: 152
出版历程
- 收稿日期: 2019-01-06
- 网络出版日期: 2022-12-07
- 刊出日期: 2020-09-30

Analytical study on water inflow and pore water pressure of tunnels under drainage

School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China

摘要

摘要: 隧洞涌水和孔隙水压力计算一直是地下工程领域最为关注的焦点问题,然而对于疏干条件下隧洞涌水量及衬孔隙水压力的变化规律仍未得到很好的解决。针对浅埋、深埋及衬砌隧洞分别进行解析分析,然后利用数值计算方法反算出关键参数（影响距离、等效地下水位及衬砌等效地下水位）的计算公式,进一步提出了隧洞涌水量的解析公式,探讨了孔隙水压力的分布规律。结果表明,疏干条件下隧洞涌水量取决于隧洞半径、初始地下水位高度及围岩与衬砌相对渗透性,其中围岩与衬砌相对渗透性起到决定性因素；衬砌圈显著改变了隧洞涌水量大小与外水压力分布规律,随着相对渗透性的增大,孔压解析解与数值解的拟合程度逐渐提高；通过与传统解析公式及数值的相互对比,提出的解析公式真实刻画了浅埋、深埋和衬砌隧洞渗流场特征,计算结果更加合理和准确。
- 隧洞 /
- 涌水量 /
- 孔隙水压力 /
- 解析公式 /
- 数值模拟
Abstract: The calculation of water inflow and pore water pressure of tunnels is always the focus of underground engineering. However, the effect of water inflow and pore water pressure of tunnels under drainage has not been well solved. The analytical analysis of shallow, deep and lined tunnels is carried out respectively. Then, the expressions for the key parameters (influencing distance, equivalent underground water level and equivalent underground water level of linings) are obtained by using the numerical method, and the analytical formulae for water inflow into tunnels are put forward. Finally, the distribution rules of pore water pressure are discussed. The results show that the water inflow into tunnels depends on the tunnel radius, the initial groundwater level and the relative permeability of the surrounding rock and linings, among which the relative permeability of the surrounding rock and linings plays a key role. The lining ring significantly changes the amount of water inflow and distribution of pore water pressure of tunnels. With the increase of the relative permeability, the fitting degree between the analytical solution and the numerical solution of pore water pressure increases gradually. Compared with traditional analytical solution and numerical solution, the proposed analytical formula more truly describes the seepage field characteristics of shallow, deep and lined tunnels, and the calculated results are more reasonable and accurate.
- tunnel /
- water inflow /
- pore water pressure /
- analytical formula /
- numerical simulation

HTML全文

0. 引言

天然气水合物是一种笼形晶状化合物（图 1），具有能量密度高、储量巨大等特点，其中97%以上分布在海底沉积物中^[1-2]。据估算，全球天然气水合物中含碳量相当于全球已探明传统化石燃料总碳量的两倍，被普遍认为是21世纪最具潜力的新型清洁能源之一^[2]。目前世界很多国家都开展了水合物勘探和试采工作，例如美国JIP计划^[3]、日本MHR & DP项目^[4]，挪威GANS项目^[5]。另外，俄罗斯、德国、韩国、印度等国家都相继开展了水合物相关的研究项目和开发计划^[6]。中国在水合物领域的研究起步较晚，但取得了重要成果。在南海北部陆坡圈定了6个水合物成矿远景区，19个成矿区带，25个有利区块和24个钻探目标区^[7-8]。特别是在神狐、东沙等海域钻探获取了水合物实物样品，发现了多种类型的高品位水合物储层，锁定了试采目标^[9]。并于2017年和2020年成功完成了两次水合物试采工作，取得了突破性进展^[10-11]。

图 1 水合物实物与分子结构图^[2]

Figure 1. Physical and molecular structural diagrams of hydrate

下载: 全尺寸图片幻灯片

目前水合物领域的研究主要集中在资源储量勘探评估、水合物开发技术和水合物地层稳定性评价3个方面。其中水合物地层稳定性评价是水合物安全可持续生产的前提条件，特别是考虑开采影响下水合物储层的动态响应及海底斜坡稳定性问题。通过海洋工程地质调查、室内模型试验和数值模拟研究发现，含水合物储层的海底斜坡失稳破坏，往往与水合物分解后产生的超孔隙压力有关^[12-14]。因此，近年来专家学者对与水合物分解相关的超孔压计算模型进行了深入研究。例如，Nixon等^[15]提出一个数学模型，通过假定水合物的分解比例，量化了不排水条件下水合物分解产生的超孔压。在此基础上，Xiong等^[16]考虑甲烷的溶解度对孔隙压力的影响，对模型进行了改进。Zhang等^[17]通过RK方程计算了真实甲烷气体的压缩系数，并将此系数引入到上述模型。Chen等^[18]利用此模型，通过考虑水合物分解过程中压力和温度的平衡，计算了水合物分解的最大超孔隙压力。此外，刘锋等^[19]、Grozic^[20]、Zhu等^[21]进行了相关研究。然而，这些研究中所采用的孔隙压力模型很少考虑水合物储层的热-水-化(THC)耦合作用，导致无法准确获取瞬态超孔隙压力及其演化特征，因此很难准确评价海底斜坡的实时动态稳定性。

本文考虑水合物注热开采的影响，发展了一种热-水-化（THC）耦合作用下海底斜坡动态稳定性简化评价方法，并自主开发了相应的数值计算程序，建立了水合物饱和度与超孔隙压力间的内在关系；进而获取了多物理场耦合作用下水合物分解引发的瞬态超孔压演化特征及储层的动态响应；并利用改进的极限平衡法对海底斜坡稳定性进行动态评价；最后以南海北部陆坡水合物试采区典型斜坡为目标，研究了水合物注热开采影响下水合物储层动态响应及实时稳定性演化。

1. 考虑THC耦合作用的超孔压确定方法

1.1 THC耦合方程

含水合物沉积物可视为由多组分（土骨架、水、水合物）和多相（固、液、气）组成的动态系统（图 2）。根据地质调查^[22]，海底水合物大多储存于相对封闭的沉积地层中，其形成由沉积物沉积、孔隙流体流动及底部甲烷供给等因素决定，且这类储层中一般缺乏明显的横向变化，因此只需考虑竖向的流体运移和能量迁移。基于此，采用Xu等^[23]提出的一维简化THC耦合模型描述海底水合物地层的赋存特性。该模型通过求解质量、能量和甲烷守恒方程，给出了水合物储层边界的计算公式，并推导了初始地层温度–压力分布和水合物聚集的表达式。

图 2 水合物沉积物多组分多相孔隙介质示意图

Figure 2. Schematic diagram of multi-component multi-phase porous sediments containing gas hydrate

下载: 全尺寸图片幻灯片

首先，根据质量守恒给出一维简化情况下流体运移的控制方程：

，

$- \frac{{k{\rho _{\text{l}}}}}{{{\mu _{\text{l}}}}}\left( {\frac{{{\text{d}}P}}{{{\text{d}}z}} + {\rho _{\text{l}}}g} \right) = {q_{\text{f}}}\text{，}$

(1)

式中，k为渗透率，ρ_l，μ_l分别为液体相的密度和动态黏度，P为地层压力，z为空间坐标，g为重力加速度，q_f为总质量通量。

其次，水合物储存系统的能量迁移可通过能量守恒方程描述，其简化模式为

，

$- \lambda \frac{{{\text{d}}T}}{{{\text{d}}z}}{\text{ + }}{q_{\text{f}}}{C_{{\text{pl}}}}T{\text{ = }}{q_{\text{e}}}\text{，}$

(2)

式中，λ为有效导热系数，T为地层温度，C_pl为水合物的比热容，q_e为总能量通量。

最后，储存系统中的甲烷运移可由甲烷介质的质量守恒方程描述：

，

$- \phi {K_{\text{m}}}{S_{\text{l}}}\frac{{{\text{d}}{M_{\text{l}}}}}{{{\text{d}}z}}{\text{ + }}{q_{\text{f}}}{M_{\text{l}}}{\text{ = }}{q_{\text{m}}}\text{，}$

(3)

式中，ϕ为沉积物孔隙度，K_m为扩散系数，S_l为液体相饱和度，M_l为液体相中甲烷浓度，q_m为与深度相关的甲烷通量。

以上3组基本控制方程间的耦合作用，其具体的建立过程为对式（1），（2）进行积分，可得到关于压力剖面P和地层温度剖面T的函数表达式：

，

$P - {P_0} = - \left( {\frac{{{q_{\text{f}}}{\mu _{\text{l}}}}}{{k{\rho _{\text{l}}}}}{q_{\text{f}}} + {\rho _{\text{l}}}g} \right)z\text{，}$

(4)

$T= \begin{cases}-\left(\xi-T_0\right) \mathrm{e}^{\frac{q_{\mathrm{e}}}{\lambda \cdot \xi} z}+\xi & \left(q_{\mathrm{f}} \neq 0\right) \\ -\frac{q_{\mathrm{e}}}{\lambda} z+T_0 & \left(q_{\mathrm{f}}=0\right)\end{cases}\text{，}$

(5)

式中，ξ为参数，ξ=q_e/(q_fC_pl)，P₀，T₀为海底面处的压力和温度。联立式（4），（5），可以得到任意深度z处的压力与温度之间的关系式。

水合物的形成需要适当的温压条件，地层温度剖面和水合物相平衡曲线共同确定了水合物的理论稳定区间，如图 3所示。Sloan等^[24]通过试验测定了甲烷水合物的相平衡曲线，表达为

$T_{\mathrm{diss}}=\left\{\begin{array}{ll} 9.6349 \ln P+197.65 & (P \geqslant 2576.74 \mathrm{kPa}) \\ 0.0158 P+232.7 & (P<2576.74 \mathrm{kPa}) \end{array},\right.$

(6)

图 3 水合物理论稳定区域与生成区域的空间关系示意图^[26]

Figure 3. Schematic diagram of spatial relationship between hydrate stability zone and hydrate formation zone

下载: 全尺寸图片幻灯片

式中，T_diss为水合物分解的相平衡温度。

将式（4）~（6）联立，可得到水合物理论稳定区域的上边界z_st和下边界z_sb：

${z_{{\text{st}}}} = 0 \text{，}$

(7)

$z_{\mathrm{sb}}=\left\{\begin{array}{ll} -\frac{\lambda}{q_{\mathrm{f}} C_{\mathrm{pl}}} \ln \left(\frac{q_{\mathrm{e}}-q_{\mathrm{f}} C_{\mathrm{pl}} T_0}{q_{\mathrm{e}}-q_{\mathrm{f}} C_{\mathrm{pl}} T_{\mathrm{diss}}}\right) & \left(q_{\mathrm{f}} \neq 0\right) \\ -\frac{\lambda}{q_{\mathrm{e}}}\left(T_{\mathrm{diss}}-T_0\right) & \left(q_{\mathrm{f}}=0\right) \end{array},\right.$

(8)

式中，上边界z_st和下边界z_sb见图 3右侧标注。另外，水合物的生成不仅需要适宜的温度和压力，还取决于孔隙中甲烷的可用量。只有当孔隙中甲烷含量大于或等于甲烷溶解度时，才可能在孔隙中形成水合物，如图 3曲线（c）和曲线（d）所示。在水合物生成区域内部，多余的甲烷会与水结合形成水合物，因此孔隙甲烷含量等于甲烷溶解度，可由下式描述：

${M_{{\text{sl}}}}(T) = {M_{{\text{sl,diss}}}}({T_{{\text{diss}}}},P,S) \cdot \exp \left( {\frac{{T - {T_{{\text{diss}}}}}}{\beta }} \right)\text{，}$

(9)

式中，M_sl为任意温度、压力及孔隙中盐浓度影响下的甲烷溶解度，M_{sl, diss}(T_diss, P, S)为甲烷溶解度关于相平衡温度T_diss、压力P和盐度S的函数，由试验测定^[25]，β=14.4 ℃为经验常数。

而在水合物生成区域上部，孔隙中甲烷含量小于其溶解度，甲烷介质以溶质赋存于孔隙中，可由式（3）描述。为确定水合物生成区域的上边界z_ot和下边界z_ob(图 3右侧标注)，在0→z_ot区间内对式（3）进行积分，并辅以式(4)，(5)，(9)可得

$z_{\mathrm{ot}}= \begin{cases}-\frac{\phi K_{\mathrm{m}}}{q_{\mathrm{f}}} \ln \left(\frac{q_{\mathrm{mt}}-q_{\mathrm{f}} M_0}{q_{\mathrm{mt}}-q_{\mathrm{f}} M_{\mathrm{sl}}\left(z_{\mathrm{ot}}\right)}\right) & \left(q_{\mathrm{f}} \neq 0\right) \\ -\frac{\phi K_{\mathrm{m}}}{q_{\mathrm{mt}}}\left(M_{\mathrm{sl}}\left(z_{\mathrm{ot}}\right)-M_0\right) & \left(q_{\mathrm{f}}=0\right)\end{cases} \text{，}$

(10)

${z_{{\text{ob}}}}: - \phi {K_{\text{m}}}{\left( {\frac{{{\text{d}}{M_{{\text{sl}}}}}}{{{\text{d}}z}}} \right)_{z = {z_{{\text{ob}}}}}}{\text{ + }}{q_{\text{f}}}{M_{{\text{sl}}}}({z_{{\text{ob}}}}){\text{ = }}{q_{\text{m}}} \text{，}$

(11)

式中，M₀为海底面处的甲烷含量，q_mt为z_ot处的甲烷通量。

综合考虑式（7），（8），（10），（11），同时满足水合物理论稳定区域条件和水合物生成区域条件即可得出水合物真实赋存区域，其上边界THZ和下边界BHZ如图 3所示。

进一步，考虑水合物形成的时间尺度和空间分布与甲烷供应有关，可给出确定整个水合物储层初始水合物饱和度的表达式^[23]。甲烷水合物填充孔隙度ϕ到任意平均饱和度S_h所需的时间t₀近似为

${t_0} \approx {S_{\text{h}}}{\left( {\frac{{\partial {S_{\text{h}}}}}{{\partial t}}} \right)^{ - 1}}\text{，}$

(12)

式中，t为时间，∂S_h/∂t为水合物的积累速率，

$\frac{{\partial {S_{\text{h}}}}}{{\partial t}} = \frac{1}{{\phi ({\rho _{\text{h}}}{M_{\text{h}}} - {\rho _{\text{l}}}{M_{{\text{sl}}}})}}\left( {\phi {K_{\text{m}}}\frac{{{{\text{d}}^2}{M_{{\text{sl}}}}}}{{{\text{d}}{z^2}}} - {q_{\text{f}}}\frac{{{\text{d}}{M_{{\text{sl}}}}}}{{{\text{d}}z}}} \right) \text{，}$

(13)

其中，ρ_h，M_h为水合物密度和水合物中甲烷的质量分数。综上，由式（1）~（13）构建了一套完整的控制方程，可用于计算地层初始状态下的温度、压力剖面和水合物饱和度剖面。

1.2 超孔压计算

水合物储层注热开采会破坏地层中温度的动态平衡，一旦地层中的温度T超过水合物的相平衡温度T_diss，水合物会产生分解。地层温度剖面T的演化可以用一维热传导方程描述：

$\frac{{\partial T(z,t)}}{{\partial t}} = \kappa \frac{{{\partial ^2}T}}{{\partial {z^2}}} + s(z,t)\text{，}$

(14)

式中，κ为地层有效热扩散系数，假设为常量，这是因为水合物的导热系数接近于水（相差 < 10%）^[27]，s(z，t)为源项，用以考虑可能的热量增加和消耗。需要注意的是，在地层升温过程中（见图 4（a），A→B），一旦温度超出水合物的稳定区域，水合物会发生分解。这其中包含水合物的相变过程，因此需要消耗一定的潜热。这意味着当地层温度T_{i, j}（下标i代表时间步，j代表某地层深度）达到T_diss后，随着热量的增加，温度会保持T_diss不变，这会延迟温度在地层中的传播（见图 4（b））。

图 4 注热开采导致的地层温度变化路径示意图

Figure 4. Schematic diagram of change path of stratum temperature caused by hydrate heat injection exploitation

下载: 全尺寸图片幻灯片

为了计算潜热的消耗，采用热积分法^[28]进行简化。具体实施过程：首先利用有限差分法求解式(14)，计算在任意时间步i时地层深度j的温度T_{i, j}。在水合物分解区域将此温度修正为T_diss，进而会产生一个温度增量dQ_{i, j}=C_ph(T_ij-T_diss)，此温度增量会根据时间步和深度更新。将每一时间步所产生的热量累计求和，可以计算出地层深度j增加的总增量∑dQ_{i, j}。一旦累积的总热量超过水合物的分解潜热，则意味着相变完成，地层温度重新按照式（14）升高。在此过程中，位于深度j处，时刻i时的水合物饱和度S_{h i, j}为

${S_{{\text{h}} i,j}} = {S_{{\text{h}} j}} \times \frac{{{L_{\text{h}}} - \Sigma {\text{d}}{Q_{i,j}}}}{{{L_{\text{h}}}}}\text{，}$

(15)

式中，S_{h, j}为深度j处的初始饱和度，L_h为水合物的分解潜热。

根据式（15）可得出在地层温度变化过程中水合物的分解量，

$\Delta {S_{{\text{h }}i{\text{,}}j}} = \Sigma ({S_{{\text{h }}i{\text{,}}j}} - {S_{{\text{h }}i - 1{\text{,}}j}})。$

(16)

最后，利用气体非逃逸方法计算温度升高过程中瞬态超孔隙压力及其演化。此方法由Nixon等^[15]首先提出，后续不少学者对该方法进行了完善^[16-21]。这里在前人工作的基础上，充分考虑水合物赋存系统的THC耦合作用，获取水合物实时分解量(式(16))，并引入到该方法中，实现升温后地层中的真实超孔隙压力演化。其中做了如下假设：①依据中国南海水合物试采区的地质调查结果，其水合物储层沉积物以细粒为主，渗透系数较小^[9]，因此忽略水合物分解后气体的渗流；②水合物赋存区的孔隙水在水合物分解前已处于甲烷饱和状态，因此忽略分解后生成甲烷在孔隙水中的溶解；③相对于气体的压缩量其余介质的压缩量很小，忽略不计。

基于以上假设，单位水合物分解后沉积物的体积变化等于生成气体的体积∆V_unit，即

$\Delta {V_{{\text{unit}}}} = {V_{\text{g}}} + {V_{{\text{wh}}}} - 1\text{，}$

(17)

式中，V_g为水合物分解产生的气体体积，V_g=164.6P_atmT_nat/(T_atmP_nat)，P，T分别为地层压力和温度，下标atm和nat分别为标准条件和地层真实条件，V_wh为水合物分解产生水的体积。故沉积物总体积变化为

$\Delta V = \phi \Delta {S_{{\text{h}}i,j}}V(\frac{{164.6{p_{{\text{atm}}}}{T_{{\text{nat}}}}}}{{{T_{{\text{atm}}}}{p_{{\text{nat}}}}}} + {V_{{\text{wh}}\theta }} - 1)。$

(18)

根据有效应力原理，超孔隙压力与沉积物变形有关，

$\begin{aligned} u_{\mathrm{e} i, j} & =-E_{\mathrm{s}} \Delta \varepsilon_{\mathrm{v}}=E_{\mathrm{s}} \frac{\Delta V}{V} \\ & =E_{\mathrm{s}} \phi \Delta S_{\mathrm{h} i, j}\left(\frac{164.6 p_{\text {atm }} T_{\text {nat }}}{T_{\mathrm{atm}} p_{\text {nat }}}+V_{\mathrm{wh}}-1\right), \end{aligned}$

(19)

式中，u_{e i, j}为时刻i时深度j处的超孔隙压力，ε_v为沉积物体应变，E_s为压缩模量，其余参数的物理意义同前。综上，可以得到地层升温后考虑THC耦合的超孔隙压力及其演化特征。

需要指出的是，在超孔压计算中做了一些合理性假设，忽略了一部分气体渗流和其余介质的压缩量，故计算出的超孔压会出现大于实际值的情况。另外，水合物分解产生的孔压会增大地层压力P_nat，进而反向影响水合物的分解，使水合物分解减缓，上述超孔压计算方法为单向耦合过程，即不考虑这种超孔压的反向耦合效应，使得海底斜坡的稳定性评价偏于安全。

1.3 方法验证

为了检验上述超孔压计算方法的正确性，复刻了Reagan等^[29]基于TOUGH+HYDRATE的模拟案例。TOUGH+HYDRATE为目前关于水合物分解模拟较为普遍的商业软件之一，并得到了广泛应用^[30-31]。此案例描述了浅层水合物赋存的典型场景，其中水深为320 m，海底温度为0.4℃，地温梯度为3℃/100m，模拟了海底水合物层100 a内，以每年0.03℃上升工况下水合物饱和度的动态响应。

对比文献[29]中设定升温工况下的计算结果，可见不同时刻（0，10，100 a）水合物饱和度的演化曲线具有较好的一致性，如图 5所示。根据式（19）可知，超孔压计算结果的合理性取决于水合物饱和度变化量∆S_{hi, j}的正确性，上述饱和度演化曲线的一致性很好的检验了本文提出的超孔压计算方法的正确性。

图 5 水合物饱和度剖面演化曲线对比图

Figure 5. Comparison of evolution curves of hydrate saturation profile

下载: 全尺寸图片幻灯片

2. 海底斜坡动态稳定性评价方法

极限平衡方法(LEM)被广泛用于海底斜坡的稳定性评价，尤其对于长度远大于沉积层厚度的无限坡。海洋工程地质调查表明，大多数水合物储层与海底面平行分布，贮藏于一定厚度的不透水覆盖层下方^[10]，因此将海底斜坡概化为由上覆层、水合物储层和下卧层组成的无限斜坡（图 6（a））。注热开采条件下，水平井附近的水合物开始分解，将水合物的分解前缘视为潜在滑动面，并以此对海底斜坡进行稳定性评价，如图 6所示。

图 6 水合物开采影响下典型海底斜坡稳定性评价示意图^[32]

Figure 6. Schematic diagram of stability evaluation of typical submarine slope under influences of hydrate exploitation

下载: 全尺寸图片幻灯片

根据极限平衡方法，海底斜坡的稳定性安全系数等于滑动面上的抗剪强度和剪应力的比值，

${F_{\text{s}}} = \frac{{\int {{{\vec \tau }_{\text{r}}}{\text{d}}l} }}{{\int {{{\vec \tau }_{\text{d}}}{\text{d}}l} }}\text{，}$

(20)

式中， ${\vec \tau _{\text{r}}}$ ， ${\vec \tau _{\text{d}}}$ 分别为滑动面处沉积物的抗剪强度和剪应力，dl为单元滑动面长度，通过建立垂直于l轴的力系平衡，可以导出：

${\vec \tau _{\text{r}}} = c' + {\sigma '_0}\tan \varphi '{\text{ = }}c' + (\sigma - {u_{\text{e}}})\tan \varphi ' \text{，}$

(21)

${\vec \tau _{\text{d}}} = \gamma 'H\sin \alpha \cos \alpha \text{，}$

(22)

式中，φ′为有效内摩擦角， ${\sigma '_0}$ 为上覆有效应力 ${\sigma '_0}$ = $\sigma - {u_{\text{e}}}$ ， $\sigma$ 为总上覆压力 $\sigma$ = $\gamma '$ Hcos² $\alpha$ ， $\gamma '$ 为沉积物有效重度，H为滑面深度，c′为沉积物有效黏聚力。水合物分解会引起沉积物胶结强度弱化，进而导致黏聚力降低^{[31, 33]}，这里假设黏聚力随着水合物的分解呈线性折减：

${c'_i} = - ({c'_1} - {c'_2})\frac{{\Delta {S_{{{\text{h}}^{}}i,j}}}}{{{S_{{{\text{h}}^{}}i,j}}}}{\text{ + }}{c'_1} \text{，}$

(23)

式中， ${c'_i}$ 为i时刻的黏聚力， ${c'_1}$ ， ${c'_2}$ 分别为初始黏聚力和最终黏聚力。

将式（21），（22）代入式（20）中，可得出海底斜坡不同深度、不同时刻的瞬态稳定性安全系数F_{s i, j}为

${F_{{{\text{s}}_{}}i,j}} = \frac{{{{c'}_i}}}{{\gamma 'H\sin \alpha \cos \alpha }} + \left( {1 - \frac{{{u_{{{\text{e}}_{}}i,j}}}}{\sigma }} \right)\frac{{\tan \varphi '}}{{\tan \alpha '}} 。$

(24)

3. 注热开采下海底斜坡稳定性实例分析

3.1 研究区域

选取南海北部陆坡白云凹陷作为研究区，区域整体位于珠江口盆地珠二坳陷内，走向近东西，面积大于2.0×10⁴ km²，是南海北部最大的一个深水凹陷，其北部与番禺低隆起相接，南端是珠江口盆地的南部隆起带^[34]。区域内发育有大量现代海底峡谷（见图 7），沉积物失稳频发。中国地质调查局于2019年10月—2020年4月在此区域进行了第二次天然气水合物试采，实现了连续产气30 d，总产气量86.14×10⁴ m³，大大提高了日产气量和产气总量^[11]。此次试采首次提出了利用水平井开采水合物的工作思路，可有效克服直井开采大压差导致二次水合物和冰生成的风险。由于复杂地质环境影响，开采引起的水合物分解可能在海底斜坡内部产生潜在滑动面，即使在坡度很小的情况下也会引发斜坡破坏。特别是考虑到海底浅表层沉积物强度低^[35-36]，易在外部作用干扰下发生土体失稳，进而可能导致海底滑坡等地质灾害，威胁海洋开采平台、锚固基础和海底管线等水下设施的安全^[37-39]。这里主要探讨水平井注热开采工况下水合物储层的超孔压演化过程及海底斜坡动态稳定性评价。

图 7 研究区位置及区域海底地形条件^{[11, 40]}

Figure 7. Location of study area and topography of regional seabed

下载: 全尺寸图片幻灯片

为服务于第二次水合物试采井的轨迹优选，在试采水平井靶点位置附近进行了先导井测试工作，并利用测井得到的P波速度和S波速度预测了水合物饱和度S_h和游离气饱和度S_g分布（图 8）^[41]。分析可见，该区域水合物储层位于海底面以下深度207.8~253.4 m，厚度为45.6 m，平均孔隙度为37.3%，平均水合物饱和度为31.0%^[41]，全井段岩性主要是泥质粉砂，渗透系数较小。另外，海洋地质调查显示此区域水深1000~1500 m，平均地形坡度为3°，海底面温度为3.3~3.7℃，地热梯度为45~67 ℃/km^[9]。基于开采井原位测试数据，构建二维海底斜坡模型，由上覆层、水合物层和下卧层组成，如图 6（a）所示。为达到最大传热效率，将水平井设置在水合物储层中间位置，即H=230.6 m处。设定由于注热影响，地层温度以0.1 K/d的速度匀速增加，总升高温度为10 K，当温度达到设定值后保持不变。模型所用具体参数见表 1。

图 8 先导井测井响应柱状图^[41]

Figure 8. Histogram logging response of pilot well

下载: 全尺寸图片幻灯片

表 1 模型参数^{[9, 11, 26, 41-42]}

Table 1. Model parameters

参数		取值	单位
坡体参数	z (深度)	—	m
	α (坡度)	3	°
	H(水平井位置)	230.6	m
地层参数	T₀(海底温度)	3.5	℃
	G/(地热梯度)	56	K/km
	S_h (初始饱和度)	31	%
	k (渗透率)	2.38	mD
	κ (有效热扩散系数)	4.5×10⁻⁷	m²/s
	λ (有效导热系数)	3.66	W/(m·K)
	c₁′ (初始黏聚力)	800	kPa
	c₂′ (最终黏聚力)	30	kPa
	φ′ (内摩擦角)	8.03	°
	ρ (体密度)	1650	kg/m³
	ϕ (孔隙度)	0.37	m³/m³
通用参数	ρ_l/(液相密度)	1024	kg/m³
	ρ_h (水合物密度)	930	kg/m³
	C_pl (液体比热)	4.18×10³	J/(kg·K)
	C_ph/(水合物比热)	2.16×10³	J/(kg·K)
	μ_l (液体动态黏度)	8.87×10⁻⁴	kg/(m·s)
	S (液体盐度)	3.5	%
	L_h/(水合物潜热)	56552-16.8/T	J/kg

下载: 导出CSV

| 显示表格

3.2 水合物储层动态响应分析

图 9为模拟工况下地层温度剖面的演化曲线，其中THZ和BHZ分别表示水合物储层的上顶面和下底面。初始地温曲线（0 a）和水合物相平衡曲线（T_diss）所围成的阴影部分为水合物的理论稳定区间（下底面为z_sb）。可以看到，由于水平井中热量的注入，水合物储层中部的温度开始升高，并逐渐向上下扩散，一旦温度超过水合物的相平衡温度，水合物开始分解（图 10）。可以观察到，在曲线T_diss附近地层的温度出现局部的降低，这是因为水合物分解是吸热过程，消耗了一定量的潜热，潜热的消耗会造成温度传导的延迟。

图 9 地层温度剖面演化曲线

Figure 9. Evolution curves of stratum temperature profile

下载: 全尺寸图片幻灯片

图 10 升温过程中水合物储层动态响应

Figure 10. Dynamic responses of hydrate reservoir during heating process

下载: 全尺寸图片幻灯片

图 10（a）为水合物饱和度的演化曲线。

可以看到，随着热量的注入储层中部水合物最先开始分解，饱和度逐渐降为0，在地层中形成一条分解带，且随着热量的传递，上下两侧的水合物也开始分解，分解带不断延伸。图 10（b）为水合物分解前缘的演化曲线，可以更加直观地看到水合物的分解及传播过程，其中红色曲线代表水合物分解前缘，黑色曲线为水合物储层的上下界面。可以看到，分解前缘最先产生于储层中部，并不断向两侧移动，直到与水合物储层上下界面重合，意味着水合物全部分解；在储层中部，水合物分解前缘的演化非常迅速，而在储层两侧则相对缓慢。这是因为储层中部靠近水平井热源，此处的水合物分解反应更加剧烈，而如上所述，水合物分解会导致地层中热量的延迟传导，故储层两侧的水合物分解需要更长的时间。

3.3 海底斜坡动态稳定性分析

图 11显示海底斜坡安全系数沿深度的变化曲线。可以看到，初始状态下海底斜坡处于稳定状态，且由于水合物的胶结作用，增加了沉积物的强度，使得初始状态下水合物储层处的安全系数明显高于地层其它位置。但随着水合物分解，海底斜坡的安全系数迅速减小，并达到临界状态。一方面是由于胶结作用的丧失，沉积物强度迅速降低，另一方面水合物分解产生的超孔隙压力减小了地层的有效应力，随着水合物分解区域的不断延伸，不稳定区域也不断扩大。

图 11 海底斜坡安全系数沿深度的变化

Figure 11. Variation of safety factor of submarine slope along depth

下载: 全尺寸图片幻灯片

图 12为海底斜坡最小安全系数的演化曲线。可以看到，经过96 d的升温后，地层中的温度达到分解条件，水合物开始分解，此时海底斜坡安全系数开始下降，在开采358 d后安全系数降为1.0，意味着海底斜坡达到临界状态。此时为了后期的持续安全开采，需要调整开采方案。例如：调整施工参数，配合降压法，既可以有效降低地层中的孔隙压力，又可以提高开采效率。

图 12 海底斜坡最小安全系数随时间的变化曲线

Figure 12. Variation of minimum safety factor of submarine slope with time

下载: 全尺寸图片幻灯片

4. 结论

（1）考虑注热开采影响下海底水合物储层的THC耦合作用，推导了水合物分解引起的瞬态超孔隙压力求解格式，并基于极限平衡分析理论建立了海底斜坡动态稳定性简化评价方法。

（2）水平井注热开采影响下，水合物分解前缘最先形成于储层中部，并不断向上、下两侧延伸；由于水合物潜热消耗的影响，导致热量延迟传导，使储层上下两侧的水合物分解需要更长时间。

（3）水合物试采区典型海底斜坡初步分析表明，由于水合物的胶结作用，初始状态下水合物储层处的安全系数偏高；升温96 d后水合物开始分解，海底斜坡安全系数迅速降低，并于358 d后到达临界状态；为保证水合物的持续安全开采，提出调整开采方案。

（4）本文基于无限坡理论分析了水合物注热开采影响下海底斜坡的动态稳定性问题，而实际含水合物层的海底斜坡失稳破坏或滑坡模式，可能表现为弧形、盘式或拉裂刺穿等破坏形式^[13]，未来应加强不同破坏模式的海底斜坡动态稳定性评价方法的深入研究。

图 1 浅埋隧洞渗流场示意图

Figure 1. Schematic representation of shallow tunnel

下载: 全尺寸图片幻灯片

图 2 深埋隧洞渗流场示意图

Figure 2. Schematic representation of deep tunnel

下载: 全尺寸图片幻灯片

图 3 衬砌隧洞渗流场示意图

Figure 3. Schematic representation of lined tunnel

下载: 全尺寸图片幻灯片

图 4 有限元模型网格剖分示意图

Figure 4. Finite element mesh of models

下载: 全尺寸图片幻灯片

图 5 隧洞涌水量数值计算结果

Figure 5. Water inflows of tunnels in different cases

下载: 全尺寸图片幻灯片

图 6 不同初始地下水位高度下疏干隧洞水平向与垂向的孔隙水压力分布

Figure 6. Distribution of pore water pressure along horizontal and vertical directions of drained tunnel under various initial groundwater levels

下载: 全尺寸图片幻灯片

图 7 不同围岩与衬砌相对渗透性下衬砌隧洞水平向与垂向的.孔隙水压力分布

Figure 7. Distribution of pore water pressure along horizontal and vertical directions of lined tunnel with different k_m/k_l

下载: 全尺寸图片幻灯片

图 8 d/D关于r/h的散点分布及拟合曲线

Figure 8. Scatter distribution and fitting curves of d/D with respect to r/h

下载: 全尺寸图片幻灯片

图 9 $h^{'} / h$ 关于 $r / h$ 的散点分布及拟合曲线

Figure 9. Scatter distribution and fitting curves of $h^{'} / h$ with respect to $r / h$

下载: 全尺寸图片幻灯片

图 10 $h^{'} / h$ 关于 $r / h$ 的散点分布及拟合曲线

Figure 10. Scatter distribution and fitting curves of $h^{'} / h$ with respect to $r / h$

下载: 全尺寸图片幻灯片

图 11 不同计算方法下隧洞中心右侧水平向与顶拱垂向的水压力分布曲线

Figure 11. Distribution of pore water pressure along horizontal and vertical directions of tunnel by different methods

下载: 全尺寸图片幻灯片

图 12 折减系数与相对渗透性之间的关系曲线

Figure 12. Relationship between discount coefficient and k_m/k_l

下载: 全尺寸图片幻灯片

图 13 不同计算方法隧洞涌水量结果对比

Figure 13. Comparison of water inflows into tunnel by different methods

下载: 全尺寸图片幻灯片

参考文献(17)

[1]	刘志春, 王梦恕. 隧道工程因素对地下水环境影响研究[J]. 岩土力学, 2015, 36(增刊2): 281-288. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2038.htm LIU Zhi-chun, WANG Meng-shu. Research on impact of tunnel engineering factors on groundwater environment[J]. Rock and Soil Mechanics, 2015, 36(S2): 281-288. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2038.htm
[2]	潘东东, 李术才, 许振浩, 等. 岩溶隧道承压隐伏溶洞突水模型试验与数值分析[J]. 岩土工程学报, 2018, 40(5): 828-836. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805009.htm PAN Dong-dong, LI Shu-cai, XU Zhen-hao, et al. Model tests and numerical analysis for water inrush caused by karst caves filled with confined water in tunnels Chinese[J]. Journal of Geotechnical Engineering, 2018, 40(5): 828-836. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201805009.htm
[3]	HARR M E, Groundwater and Seepage[M]. New York: Dover Publications, 1962: 249-264.
[4]	吴建, 周志芳, 李鸣威, 等. 隧洞涌水量预测计算方法研究进展[J]. 工程地质学报, 2019, 27(4): 890-902. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201904023.htm WU Jian, ZHOU Zhi-fang, LI Ming-wei, et al. Advance on the methods for predicting water inflow into tunnels[J]. Journal of Engineering Geology, 2019, 27(4): 890-902. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201904023.htm
[5]	GOODMAN R E, MOYE D G, VAN SCHALKWYK A, et al. Ground Water Inflows During Tunnel Driving[M]. California: College of Engineering, University of California, 1964.
[6]	LEI S. An analytical solution for steady flow into a Ttonnel[J]. Groundwater, 1999, 37(1): 23-26. doi: 10.1111/j.1745-6584.1999.tb00953.x
[7]	RAYMER J H. Predicting groundwater inflow into hard-rock tunnels: estimating the high-end of the permeability distribution[C]//2001 Rapid Excavation and Tunneling Conference, 2001, Littleton: 1027-1038.
[8]	EL TANI M. Circular tunnel in a semi-infinite aquifer[J]. Tunnelling and Underground Space Technology, 2003, 18(1): 49-55. doi: 10.1016/S0886-7798(02)00102-5
[9]	MOON J, FERNANDEZ G. Effect of excavation-induced groundwater level drawdown on tunnel inflow in a jointed rock mass[J]. Engineering Geology, 2010, 110(3/4): 33-42.
[10]	苏凯, 陈锐, 周亚峰, 等. 圆形隧洞渗流量计算的解析方法及其应用[J]. 岩石力学与工程学报, 2017, 36(增刊1): 3332-3341. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2017S1024.htm SU Kai, CHEN Rui, ZHOU Ya-feng, et al. Analytical method for calculating the groundwater inflow of circular tunnel and its application[J]. Rock and Soil Mechanics, 2017, 36(S1): 3332-3341. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2017S1024.htm
[11]	EL TANI M, KAMALI A, GHOLAMI M A. Analytic assessment of the water table drawdown, seepage, and back pressure at rudbar PSPP[J]. Rock Mechanics and Rock Engineering, 2019, 52(7): 2227-2243.
[12]	HOOGHOUDT S B. Bijdragen Tot De Kennis Van Eenige Natuurkundige Grootheden Van Den Grond: Algemeene Beschouwing Van Het Probleem Van De Detailontwatering En De Infiltratie Door Middel Van Parallel Loopende Drains, Greppels, Slooten En Kanalen[R]. Hague: Algemeene Landsdrukkerij, 1940: 515-707.
[13]	FERNANDEZ G, MOON J. Excavation-induced hydraulic conductivity reduction around a tunnel-part 1: Guideline for estimate of ground water inflow rate[J]. Tunnelling and Underground Space Technology, 2010, 25(5): 560-566.
[14]	FERNANDEZ G, MOON J. Excavation-induced hydraulic conductivity reduction around a tunnel-part 2: Verification of proposed method using numerical modeling[J]. Tunnelling and Underground Space Technology, 2010, 25(5): 567-574.
[15]	周亚峰, 苏凯, 伍鹤皋. 水工隧洞钢筋混凝土衬砌外水压力取值方法研究[J]. 岩土力学, 2014, 35(增刊2): 198-203. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2014S2027.htm ZHOU Ya-feng, SU Kai, WU He-gao. Study of external water pressure estimation method for reinforced concrete lining of hydraulic tunnels[J]. Rock and Soil Mechanics, 2014, 35(S2): 198-203. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2014S2027.htm
[16]	DAI Y, ZHOU Z. Steady seepage simulation of underground oil storage caverns based on Signorini type variational inequality formulation[J]. Geosciences Journal, 2015, 19(2): 341-355.
[17]	郑宏, 刘德富, 李焯, 等. 一个新的有自由面渗流问题的变分不等式提法[J]. 应用数学和力学, 2005, 26(3): 363-371. https://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200503016.htm ZHENG Hong, LIU De-fu, LEE C F, et al. New variational inequality formulation for seepage problems with free surfaces[J]. Applied Mathematics and Mechanics, 2005, 26(3): 363-371. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200503016.htm

施引文献(2)

期刊类型引用(2)

1.	刘姝，李文杰，李莉佳，罗永江，陶瑞，李晓璇，杨亚会. 天然气水合物开采方法研究现状及展望. 钻探工程. 2024(05): 12-23 . 百度学术
2.	于倩男，唐慧敏，李承龙，梁爽，梁家修，陈志静，张琨. 天然气水合物注热分解渗流特征及数值模拟. 东北石油大学学报. 2023(06): 38-54+127-128 . 百度学术