天然气水合物开采的土力学问题：现状与挑战

韦昌富 颜荣涛 田慧会 周家作 李文涛 马田田 陈 盼

1.岩土力学与工程国家重点实验室·中国科学院武汉岩土力学研究所 2.广西岩土力学与工程重点实验室·桂林理工大学

0 引言

天然气水合物（以下简称水合物）可存在于海洋大陆架边缘地层或寒区永久冻土层中[1]，本文主要讨论海底水合物储层的情况。海底水合物储层不同于冻土储层，除了环境温度在0 ℃以上以及孔隙中不含冰外，所处水深大、埋深浅，骨架固结程度弱或没有固结，工程地质条件复杂。此外，海底水合物储层不仅具有普通深海沉积物所具有的内部结构与工程地质特性，而且还因水合物的可相变性，其内部结构及工程特性随着水合物开采过程均在不断地演变中。可相变性指在合适温度、压力及水化学条件下，沉积物孔隙中固体水合物可以分解成液态水和天然气；反之，孔隙中液态水和天然气也可生成水合物。正是这种可相变性，使得含水合物沉积物在环境荷载作用下表现出比普通海底沉积物更为复杂的工程力学行为特征，并伴有显著的多相、多场、多过程耦合效应。特别地，在水合物开采过程中，水合物分解产生天然气和液态水并吸收周围的热量，从而引起渗流、扩散、热传导等耦合过程；同时，随着水合物的分解，水合物的胶结作用逐渐受到破坏，沉积物骨架抗剪强度减弱，压缩性增大，从而导致沉积物变形或骨架坍塌，最终引发流泥、流沙、井壁破坏、储层沉降或失稳、大范围海底滑坡等重大工程或地质灾害问题[2-3]。因此，有效确定含水合物沉积物的工程力学特性并厘清其随着水合物分解过程的变化规律，是实现水合物安全、高效开采的前提保证。

笔者围绕着水合物开采中所涉及的关键土力学问题，从含水合物沉积物的物理力学特性的关键测试技术、含水合物沉积物的变形与强度特性、含水合物沉积物的物理力学模型、水合物开采过程的数值模拟方法等几个方面，分析和评述国内外相关研究现状和进展，并在此基础上探讨这些方向的未来发展趋势及所面临的挑战。从岩土工程的角度来看，水合物开采还涉及许多其他的重要方面，如原位取样、原位土工测试、原位监测、海底稳定性分析、深海基础工程（如锚固基础）等。毫无疑问，这些方面的工作也是不可或缺的，但本文除了在讨论未来发展趋势时有所涉及外，重点只放在含水合物储层的工程力学特性方面。本文中，“含水合物沉积物”和“水合物储层”这对名词可以互用。在一些文献中[4]，含水合物沉积物也被称为“能源土”，但本文将采用国际比较流行的说法，即“含水合物沉积物”（hydratebearing sediment）[2]。

1 物理力学特性的测试方法和技术

1.1 含水合物沉积物试样的制备方法

由于水合物稳定条件较为苛刻，钻取和运输过程中压力和温度的变化会导致水合物分解，使得原状含水合物沉积物试样的获取与运输极为困难，因此，实现人工制备含水合物沉积物试样显得尤为重要。在此之前，了解水合物于沉积物中的赋存模式则是人工制备含水合物沉积物试样的前提条件。

自然界中，天然气水合物主要存在扩散型、渗漏型及复合型等3种成藏模式[5-6]。在扩散型成藏模式下，水量丰富，天然气渗漏通量低，称为富水环境，气体只能通过扩散进入水合物稳定带与水结合形成水合物；而在渗漏型成藏模式下，构造裂隙发育，天然气渗漏量高，称为富气环境，下部游离气体沿着构造裂隙往上运移至水合物稳定区和当地的孔隙水结合形成水合物[7-8]；复合型兼具扩散型和渗漏型两者的成藏特征。因此，在这3种典型环境中，水合物在粗粒土中主要以4种模式存在：①孔隙填充模式。水合物颗粒主要悬浮于孔隙水中，与附近土颗粒无实质接触[9]。②持力体模式。这种模式由孔隙填充模式发展而来，当水合物含量达到一定量后部分水合物起支撑作用并桥接相邻的颗粒，从而成为承重骨架的一部分。一般地，当水合物饱和度超过25%～40%时，填充型水合物会变成持力型水合物[9-12]。③胶结模式。水合物优先在多孔介质颗粒接触处及附近生成，使土颗粒接触并黏合在一起，从而胶结成一个固体骨架[13-14]。④包裹模式。水合物包裹于土颗粒周围。在细粒沉积物中，水合物一般以分散的结核状、透镜体或脉状块体形式存在[15-18]。

为了模拟原位的水合物赋存情况可在实验室制得含水合物沉积物试样，目前实验室制备试样使用的主要有饱和法、非饱和法、冰种子法、混合法等。饱和法是指在一定压力条件下制备气体的超饱和溶液，将其注入沉积物后，通过降低温度在沉积物孔隙中生成水合物[19]。这种方法可以很好地制得孔隙填充模式和持力体模式试样，但由于甲烷等天然气体的溶解度较低，该方法制备的水合物饱和度较低且耗时。为了克服这一困难，有学者尝试利用添加剂来增大甲烷气体溶解度[20-21]。非饱和法利用沉积物在非饱和情况下具有较高气—水接触面积的特点，使过量气体与孔隙水在合适温度、压力条件下充分反应，从而在沉积物中生成较高饱和度的水合物，该方法已被证实能很好地模拟胶结模式和包裹模式的试样[22-24]。冰种子方法是将冰粉与土混合制样或将试样预先冻结，在供气条件下逐渐升温，孔隙冰融化后与气体结合形成水合物[25-26]。该方法可以加快水合物的生成过程，但在制样过程中存在冰粉压融现象，如果预先冻结试样，则可能因为“冻吸力”引起水分迁移而造成水合物分布不均。上述3种方法均是先制备沉积物，再通过控制温度和压力条件在沉积物孔隙中形成水合物。为了进一步提高水合物制备效率及水合物饱和度，Hyodo等[27]通过先制备水合物粉，然后与沉积物按一定比例混合，压制成含水合物沉积物试样，这类方法称为混合制样法。由于水合物的形成条件较为苛刻，对试验条件要求较高，应用这种方法的关键在于抑制制备过程中水合物粉的分解，避免水合物含量计算不准确。

总体上，虽然目前的实验室制样技术能够初步满足含水合物沉积物物理力学试验要求，但仍有以下几个关键问题亟待解决：①制备大尺度的孔隙填充模式和持力体模式的含水合物试样，以供开展物理模型试验；②克服制样过程中的重力和吸力的影响[28]，制备物理力学试验所需的均匀含水合物沉积物试样；③在水合物形成过程中实现含量的精准控制与测量。

1.2 水合物含量的确定方法

含水合物沉积物的物理力学特性与水合物含量密切相关[12,23-24,27]，因此准确地确定水合物含量是研究含水合物沉积物物理力学特性的关键问题之一。目前常用的水合物含量的确定方法主要包括[29-31]：压降法计算、含水量消耗计算、时域反射技术、核磁共振技术等。压降法是指通过水合物生成前后气压变化计算出气体消耗量从而计算得到水合物的生成量。在非饱和法制样中，客体气体分子一直保持过量，通过假定沉积物中水分完全参与水合反应生成水合物，即可利用所消耗的水量来计算出水合物饱和度。而时域反射（TDR）技术则是通过测定水合物生成前后孔隙水介电常数变化来反算水合物饱和度。近年来，田慧会等[31]利用低场核磁共振技术确定了天然气水合物的水合数，进而求得了水合物含量大小。现有的各种水合物含量的确定方法均存在一定局限性，如测量精度难以确定、不同方法间缺乏一致性等，研发高精度的沉积物水合物含量测量方法仍是一项亟需开展的任务[32]。

1.3 微细观组构探测技术

现有的沉积物组构的探测技术主要包括：①声波技术[33-35]。即根据实测的声波参数如波速、衰减及其频散特性反演水合物沉积物的孔隙度、水合物饱和度、弹性模量等参数，并分析水合物在沉积物中的赋存状态。②电阻率法[36-38]。基于沉积物电阻率与复电阻率频散特性随着水合物含量及赋存模式的变化规律，测定沉积物中水合物的含量与分布。③TDR技术[29,39]。根据介电常数与含水量的关系，推算出沉积物试样中水合物含量的变化。④计算机辅助成像（CT）技术[40-42]。利用含水合物沉积物不同组分对X射线衍射的吸收与透过率不同而形成含水合物沉积物的微观图像，并由此观测和表征水合物体积以及剪切过程中微细观结构变化。⑤核磁共振（NMR）技术[31,43-44]。利用孔隙水核磁信号与弛豫时间探测含水合物沉积物中水合物含量及赋存状态。

在以上5种探测技术中，前3种探测技术主要根据含水合物沉积物的波速、电阻及介电常数随沉积物各组分和组分含量的变化关系来反演其内部微细观组构，此类办法能定性反映体系中水合物相变过程以及微细观结构变化，但不适用于在微细观层面进行定量探测和分析。与其他方法相比，CT技术和NMR技术是两种极具潜力的组构探测及参数测定技术，其最大的优势在于它们是直接、定量、快捷的测试方法。但是，X射线衍射存在辐射效应需要特殊材料进行试验防护，NMR技术则必须保证相关试验设备对含水合物沉积物无核磁信号干扰。另外，昂贵的价格也限制了CT技术和NMR技术在水合物研究领域的应用。

1.4 水力特性的测试方法

含水合物沉积物的水力特性包括渗透性和持水曲线。传统上，沉积物渗透率的测定方法主要有非稳态法和稳态法。稳态法是在稳态渗透状态下通过测定试样两端的压力差及孔隙水的渗流量，再由达西定律求得渗透系数。在含水合物沉积物渗透系数的测定试验中，稳态法比较流行[45-47]。然而，稳态法的测试时间比较长，在测试过程中水合物易发生分解，从而影响到测试精度。相比之下，非稳态法所需时间较短，只需要测定水流量（或孔隙压力）随着时间变化关系，然后利用数学模型来反演得到渗透率。目前，非稳态法已开始应用于含水合物沉积物渗透率的测定中，例如：Jaiswal[48]采用JBN非稳态法测定了含水合物沉积物的有效渗透率；张宏源等[49]基于瞬态压力脉冲法测定了含水合物粉细砂的渗透率。除了上述两种方法，一些间接的渗透系数方法，例如，Kleinberg[50]利用T2分布特征采用Kenyon公式计算了含水合物沉积物中的相对渗透率，这种方法在石油领域测井中广泛应用；Kneafsey等[51]基于CT技术监测沉积物中不同组分的密度变化，推算出了含水合物砂土的渗透系数。总之，含水合物沉积物渗透性的测试难点在于水渗流过程会引起水合物的溶解或分解，减小水合物含量。因此，如何有效避免或减小水合物含量的变化是优化含水合物沉积物渗透性测试的关键。

持水曲线是指沉积物中吸力与液态水含量的关系，其受到水合物含量的影响。Mahabadi等[52]在传统的非饱和土持水曲线测试装置上增加了控温系统，测试了含THF水合物砂土在不同水合物饱和度条件下的持水特征曲线；Ghezzehei等[53]对压力室进行了改造，使试样处于高压环境，并且提高溢出水端口压力保证气—水压力差在高进气值陶土板可控范围内，同时增加了控温装置来控制系统温度，在此基础上测得了含水合物沉积物的脱湿和吸湿曲线。由于持水特征曲线测试非常耗时，其中水渗出过程可能导致水合物分解，而施加的气体压力则有可能导致水合物进一步形成，从而影响到测试结果的精度。为了克服这一问题，有学者采用物理仿真模拟来研究水合物赋存对含水合物沉积物持水特征曲线的影响[54-55]。目前有关含水合物沉积物持水曲线的直接测试工作非常稀缺，主要问题在于如何通过控制温度、压力以及其他手段来保证测试过程中水合物不产生二次形成或分解，以保证水合物含量的恒定不变，从而提高测试结果的精度。

1.5 力学特性的测试技术

力学特性包括应力应变关系、压缩特性、抗剪强度特性等，可通过直剪仪、三轴剪力仪等仪器来测定。尽管直剪仪操作简便，适合于工程应用，但目前利用直剪仪来测试含水合物沉积物抗剪强度的报道非常少。中国科学院武汉岩土力学研究所是最早开展这方面研究的单位之一，研制了国内首台含水合物土直剪仪（图1）。该设备将剪切盒安装于反应釜内，通过气瓶提供生成水合物所需的高压气体，用恒温箱控制反应釜内试样的温度，最终形成水合物沉积物试样；通过活塞式加载杆实现法向力和切向力的加载与卸载，用气瓶输出压力提供法向压力而用恒流泵提供切向推力。Liu等[56-57]采用该设备开展了一系列含CO2水合物沉积物的直剪试验，揭示了含水合物沉积物的抗剪强度随水合物饱和度的变化规律。

三轴剪切试验是测试含水合物沉积物力学特性最常用的方法，可以测得应力—应变关系以及变形和强度特性。表1列出了国内外各主要研究单位的含水合物沉积物三轴剪切系统的基本信息。各种三轴仪的主要差别在于所采用的试样制备、水合物含量确定、试样体变测量等方法上，这些正是含水合物沉积物三轴试验的难点。在大多数三轴试验中，试样均采用非饱和法制备，通过配置不同初始含水量来控制水合物饱和度。假设试样中孔隙水完全转化成水合物，Ghiassian等[24]及Hyodo等[66]通过确定试验结束后水合物分解所释放的气体量反算得到水合物含量大小。目前在三轴试验中体变测定主要有以下3种方法：①在恒定围压时，通过围压泵驱动液体积变化换算获得[24,65]；②类似于非饱和土试样体积量测的双压力室法[66]；③采用径向应变传感器进行量测[67]。为了评价现场所取的含水合物沉积物试样的力学特性，Yoneda等[68]开发了一套新型的三轴试验系统（TACTT），可将原位取得含水合物试样直接嵌入三轴压力室进行试验，并保证整个过程中试样处于稳定状态。这套系统采用了图像处理技术来监测三轴剪切过程中试样的体积变化。

图1 含水合物土直剪仪的结构示意图[56-57]

表1 各种含水合物沉积物三轴试验系统对比表

目前有关水合物沉积物的工程力学特性试验仍然存在许多问题，如测试方法与标准没有统一、各种试验数据之间缺乏可比性、试样制备方法单一、试验结果不能真实反映原位沉积物的行为特性等。这些均为目前亟待解决的难题，是弄清含水合物沉积物工程力学特性，揭示刚度、强度等物理力学参数的变化规律以及识别其关键影响因素所必须解决的关键问题。

2 含水合物沉积物的变形与强度特性

2.1 含水合物沉积物的抗剪强度特性

含水合物沉积物的抗剪强度是指沉积物达到破坏时所能承受的剪应力。Winters等[69]通过试验证明了水合物沉积物的抗剪强度取决于孔隙中水合物、液体水、天然气、冰等填充物的相对数量和空间分布形态。李洋辉等[70]将水合物颗粒与土颗粒混合制成试样并开展了三轴试验，发现水合物含量增加显著提升了试样的抗剪强度，且水合物的主要贡献在于增加黏聚力，而对内摩擦角影响较小。Masui等[25]通过开展三轴剪切试验发现，冰—砂混合物制得试样只有在水合物饱和度超过25%时才能提升其抗剪强度，而水—砂混合物制得试样的抗剪强度随水合物饱和度的增大而显著上升。Miyazaki等[71]的试验结果表明，含水合物砂的抗剪强度和刚度随水合物饱和度和有效围压的增加而增加，而割线泊松比随有效围压增加而降低。颜荣涛等[23]分别用非饱和法、饱和法制备水合物沉积物试样并开展三轴试验，发现非饱和法制作的水合物试样抗剪强度随饱和度增大而显著提高，而饱和法制作试样的力学性质受水合物含量的影响不大；这种差异性可归因于水合物赋存模式的不同，非饱和法制备的试样中水合物以胶结形式存在，而饱和法产生填充型水合物。Liu等[57]采用直剪仪研究了CO2水合物沉积物的抗剪强度特性，发现峰值抗剪强度和残余抗剪强度均随水合物饱和度和法向应力的增加而提高。综上所述，水合物含量和赋存模式对含水合物沉积物强度都会产生明显影响，并且水合物主要贡献于增加黏聚力，而对内摩擦角影响较小[12,23,25,71]。

2.2 含水合物沉积物的剪胀/剪缩特性

研究结果表明，水合物的形成会影响含水合物沉积物的剪胀/剪缩特性。Miyazaki等[71]用人工合成水合物砂试样进行三轴试验，在轴向荷载逐渐增加，试样逐渐从剪缩变为剪胀，偏应力也在到达最高点后迅速下降呈现应变软化现象。Hyodo等[66]的试验也证实了水合物沉积物的剪胀性随着水合物饱和度的升高而明显增大。Masui等[72]利用水合物沉积物的钻井岩心进行三轴试验，发现与室内合成的水合物砂相比，剪胀性弱于人工合成水合物试样。从微观层面，Pinkert等[73]认为在剪切作用下水合物颗粒破碎损伤引起了试样应变软化的出现，破碎后的水合物颗粒与砂颗粒一起运动，相互翻转跨越导致试验出现剪胀变形。显然，水合物生成使含水合物沉积物试样从开始阶段的应变硬化向应变软化转变，同时伴随着剪胀现象。随着水合物饱和度的增加，应变软化和体积剪胀性会更加明显[66,74-77]。另外，增大的围压会限制颗粒之间的相互翻转跨越运动，使沉积物的剪胀性变弱，甚至出现剪缩现象[78-79]。

水合物沉积物的剪胀机理可用图2来进行说明。剪切面两侧土颗粒受到法向力FN和剪切力FT，作用在每个颗粒上。当颗粒松散排列时切向变形使其排列更密集，从而使土体收缩并显示出应变硬化行为（图2a→b）；相反，当颗粒处于密实状态时，在剪切过程中则可能发生体积膨胀即剪胀（图2b→a）。一旦水合物形成，孔隙空间就会被水合物部分填充，从而使土体的固体骨架变得更加致密，因此剪胀性随水合物饱和度增大得以提高（图2-c）。如图2-d所示，取剪切面上局部土颗粒和水合物颗粒对其进行受力分析，FN和FT的合力与土颗粒和水合物颗粒接触面构成夹角α。在较低围压时FN较小，夹角α较小，土颗粒和水合物颗粒发生滑动跨越导致土体膨胀增大而发生剪胀，同时由于胶结结构丧失导致应变软化变形特征（图2-e）；当围压增大时，FN增大导致夹角α增大，当α大于土颗粒与水合物颗粒之间的摩擦角时，土颗粒与水合物颗粒组成的系统将产生自锁效应而无法沿接触面滑动，水合物颗粒被不断增大的FT压碎，从而形成更紧密的结构表现为剪缩（图2-f），在此过程中水合物沉积物应变软化特征减弱甚至转变为应变硬化。

2.3 含水合物沉积物的小应变刚度

小应变刚度指的是应变为10－5～10－8情况下试样的变形模量或剪切模量，声波探测是研究试样小应变刚度的有效方法。Lee等[81]在砂土、粉土和黏土中形成THF水合物，通过测试试样的P波和S波研究了小应变刚度特性。其试验结果表明，含水合物沉积物的P波和S波受水合物饱和度、有效应力和沉积物比表面积的影响显著。由于水合物的胶结和填充作用，含水合物沉积物的小应变体积模量和剪切模量会随着水合物增大呈现明显的上升趋势[82-84]。Priest等[84]采用“富气法”和“富水法”分别制备含水合物砂土试样，发现两种试样显示明显不同的波速响应。对于“富气法”形成试样，水合物主要起到胶结作用，少量的水合物就可以明显提高含水合物沉积物的刚度；而“富水法”形成试样，水合物主要填充于孔隙中，只有当水合物达到一定含量时，水合物赋存模式从填充型转变为持力体型，含水合物沉积物的刚度才能得以提升。这一结论在Choi等[34,82]所进行的实验室中也得到了验证。可见，水合物含量及赋存模式对沉积物的小应变（10－5～10－8）刚度具有显著的影响[9-11,85]。

图2 水合物沉积物剪切变形机理示意图

基于水合物含量及赋存模式对含水合物沉积物试样的小应变刚度的影响机制，研究人员针对于不同水合物赋存模式建立了预测含水合物沉积物波速的理论模型，通过测定剪切波速、压缩波速与水合物饱和度关系来综合判别沉积物中水合物的赋存模式[86-87]。这种水合物赋存模式的识别方法目前已得到一定程度的应用[34,88-90]。

综合上述分析可知，水合物含量及其赋存模式对沉积物强度和变形特性的影响已经成为当前重点研究的问题。尽管现有的研究成果已初步揭示了水合物含量和赋存模式对沉积物物理力学特性的影响规律，但囿于制样方法与实验技术，目前对含水合物沉积物物理力学特性的研究仍然存在以下几个问题：①对在水合物含量变化下填充型向持力体型转变过程中沉积物的物理力学变化规律缺乏深入系统的认识；②对水合物分解过程中因水合物局部分解而导致的含水合物沉积物抗剪强度的衰减及压缩性的增加规律几乎没有涉及；③流泥流沙情况下含水合物沉积物骨架的物理力学特性衰化规律仍然没有厘清。另外，目前有关沉积物物理力学特性的试验数据均在特定试验条件下获得，缺乏一致性、可比性和系统性，需要针对典型水合物储层介质的物理力学特性开展系统深入的试验研究工作。

3 含水合物沉积物的物理力学模型

3.1 沉积物体系中水合物相平衡模型

沉积物中水合物相平衡模型描述了水合物在平衡时温度、压力与其他状态变量（如水合物含量、液态水含量、含盐量等）之间的关系，是含水合物沉积物最基础的本构关系之一。最早的水合物相平衡模型是由 van der Waals 等[91]所建立的，主要预测非沉积物中纯水合物的相平衡压力和温度条件。在沉积物中，由于沉积物颗粒表面物理吸附和孔隙毛细管效应的影响，孔隙水活度降低影响了相平衡条件，因此相平衡模型应该考虑沉积物孔隙大小及分布特征的影响。Clennell等[92-93]将Gibbs-Thomson方程引入vdWP模型，从而考虑了孔径引起的毛细管效应，建立了沉积物中水合物相平衡模型。该模型只适合于描述具有单一或窄孔径分布的沉积物中水合物相平衡条件。为突破这一局限性，颜荣涛等[94]引入了孔径分布函数，并根据水合物分解的孔径大小顺序，确定了孔径大小分布与水合物饱和度的关系，从而得到了p—T—Sh三维相平衡关系。最近，利用沉积物孔径分布与持水特性的内在关联性，颜荣涛等[95]通过引入持水曲线，建立了沉积物中水合物的宏观相平衡模型。

近年来，Wei[96]基于多孔材料的连续介质理论建立了能考虑骨架与孔隙溶液之间复杂物理化学作用的多相岩土介质化学—力学耦合理论，并提出了能考虑渗透、毛细管、吸附等效应的孔隙溶液各组分化学势的一般数学表达式。在此基础上，Zhou等[97]建立了能够考虑复杂物理化学效应的水合物相平衡方程，该方程给出了平衡温度偏移量与毛细管吸力和盐溶液浓度的关系，揭示了沉积物持水特性与相平衡条件的内在联系。结合所提出的相平衡方程和沉积物的持水特征曲线，进而提出了水合特征曲线（SHCC）的概念，即含水合物沉积物中液态水含量与平衡温度偏移量的关系曲线。在p—T—w三维空间中（w为液态水含量），SHCC描述了沉积物中水合物的相边界。如图3所示，对于任意指定的水合物沉积物，在平衡时含水合物沉积物中的温度、压力和液态水之间存在唯一关系，这与实验观测结果是完全一致的。该模型克服了传统水合物相平衡模型不能预测任意平衡的温度和压力下沉积物中液态水含量问题，不仅形式简单（没有涉及微观参数），而且还能综合考虑盐分、毛细管以及物理化学吸附的影响。

图3 在p—T—w空间中含CO2水合物沉积物的相平衡面及其在各平面上的投影[97]

3.2 含水合物沉积物水力模型

含水合物沉积物的水力模型包括渗透率函数模型和持水曲线模型，其中，渗透率函数不是严格意义上的本构关系，而只是一条辅助的经验关系式。在处理土体渗透性时，通常的做法是根据孔隙分布特征把孔隙简化成不同孔径大小的平行的毛细管，再利用流体力学的原理来确定其渗透性[98]。基于这种思想，考虑到孔隙包裹型和孔隙填充型对平行毛细管特征的影响，Kleinberg[50]分别建立了能够考虑孔隙包裹型与填充型水合物赋存模式的水合物沉积物渗透率模型。然而，这类模型没有考虑土体中渗流路径的曲折性，也没有考虑水合物形成后对孔隙联通性的改变，因此该模型的预测结果与实际情况有一些差距。在自然界中，根据所处环境条件与水合物含量，包裹型和填充型水合物可同时出现于海底沉积物中，并且同时对渗透率存在影响，只是程度不同而已。针对于这种情况，Delli等[99]提出了一个混合式模型，该模型可以考虑包裹型和填充型两种水合物模式共存的情况。此外，对于水合物分布不均匀的情况，Dai等[100]修正了Kozeny-Carman模型，较好地模拟了渗透率随水合物饱和度的变化。

含水合物沉积物的持水特性控制了气—水渗流过程，并影响着沉积物的其他工程力学特性。在有关水合物开采的边界值问题分析中[101-102]，一般均采用非饱和土的van Genuchten模型来描述含水合物沉积物的持水特性[103]。然而，由于水合物的生成/分解，含水合物沉积物的孔隙大小及分布都在不断变化中，因此沉积物持水特性也随之发生变化。Gupta等[104]通过引入修正系数来考虑水合物饱和度和孔隙率对毛细管吸力的影响。该模型主要是从模仿孔隙中微生物生长对骨架产生的影响修正而来，并没有真正考虑水合物生成对土体结构的影响。Yan等[105]分析了水合物形成对含水合物土毛细管吸力的影响机理，考虑了水合物填充的效应，建立了含水合物沉积物的土水特征曲线模型。虽然该模型很好地描述了含水合物沉积物的持水特征关系，特别是水合物对进气值和残余饱和度的影响，但是其有效性仍需要进一步验证。

总体上，现有的多数含水合物沉积物的水力模型（渗透率函数和持水曲线模型）都是基于高度简化的含水合物沉积物结构，对现有的水力模型修正来建立的，所以其模拟精度及适用性受到限制。因此，后期发展的水力模型应该能够考虑水合物含量及赋存模式对孔隙水/气迁移过程的影响，在大量试验结果基础上，建立能统一考虑水合物含量及赋存模式影响的含水合物沉积物的渗透率函数和持水曲线模型。

3.3 应力应变本构关系模型

含水合物沉积物的应力应变本构关系大体上可分为非线性弹性模型和弹塑性模型两大类。前者主要包括Duncan-Chang类[106-109]和弹性损伤模型[110-113]等。Miyazaki等[106]建立了含水合物沉积物强度和刚度与饱和度的经验关系式，修正Duncan-Chang模型来模拟了含水合物沉积物的应力应变关系。Yu等[107]把含水合物沉积物在受荷过程中的破坏分成结构性快速破坏和结构性完全破坏（即屈服）两个阶段，并且在初始切线模量和极限偏应力中考虑了温度的影响，建立了包含温度、围压、应变速率等不同参数的修正Duncan-Chang本构模型。这类模型的优势在于形式简单，参数易于测定，易于在数值模型中进行求解。然而，该类模型虽然考虑了水合物饱和度对试样强度和刚度的增强效应，但不能模拟应变软化及剪胀行为。下一步工作中，可以重点考虑对Duncan-Chang本构模型进行修正，使其能反映含水合物沉积物的应变软化或剪胀行为。

吴二林等[110-111]较早引入弹性损伤概念来模拟含水合物沉积物的本构行为，后期经过了李彦龙、祝效华、颜荣涛等[112-114]研究者的发展。该模型的主要思想是将沉积物视为由水合物和土体骨架组成的复合材料，利用混合律计算出沉积物的等效弹性模量，将统计损伤、概率等概念考虑进去，使其能反映沉积物受力时水合物对微结构和应力应变的影响，从而建立含水合物沉积物的弹性损伤模型。此类模型较好地反映了剪切过程中含水合物沉积物中的损伤过程，从而有效模拟了含水合物沉积物的应变软化行为，但该模型不能很难描述剪切过程中试样体积的变化。

Pinkert及其合作者主要发展了含水合物沉积物的Mohr-Coulomb类本构模型[73,115-116]，该模型很好地模拟了含水合物沉积物的应变软化效应，采用的Rowe应力剪胀理论对含水合物砂土剪胀行为有较好的描述。尽管如此，Mohr-Coulomb类模型还是没有充分考虑体积屈服效应以及剪切过程中水合物胶结作用衰减或丧失所带来的影响。

为了更全面地描述含水合物沉积物的本构行为，一些学者尝试在临界状态弹塑性理论框架内建立水合物沉积物的本构模型。目前比较常用的方式主要是通过改变屈服面大小来反映水合物的影响，结合硬化规律、屈服函数、流动法则和加卸载准则及弹性行为模拟，给出了一套内容完整、物理意义明确、结构清晰的理论体系。例如，Klar等[102]所建立的理想弹塑性模型、Sultan等[117]和Uchida等[118]提出的修正剑桥模型、Zhang等[119]提出的修正亚塑性模型等、Lin等[120]基于SMP准则的弹塑性本构模型等。在目前发表的临界状态弹塑性模型中，主要关注含水合物沉积物的以下几个典型特征：①强度和刚度随水合物含量的增大；②应变软化效应会随着水合物含量的增加而增强；③随着水合物含量的增大，试样的剪胀效应会明显增大[117-118,120-123]。这些模型主要通过修正原剑桥模型中的屈服函数和剪胀函数来考虑水合物含量对含水合物沉积物力学特性的影响规律。这些屈服函数和剪胀函数的建立主要依赖于理论假设，缺少足够的试验数据支持。此外，一些其他影响因素也在临界状态模型中进行了考虑。例如，Shen等[124]分析了水合物对临界状态以及剪胀函数的作用，建立了依赖密实度的含水合物砂土的临界状态模型；Yan等[121]重新定义了含水合物土的有效应力，并且在屈服函数、流动法则和硬化规律中嵌入了水合物胶结和毛细管吸力的影响，提出了有效考虑水合物饱和度和赋存模式对力学特性、模拟分解诱发的土体变形的含水合物土的本构模型；邹远晶等[125]考虑重点分析了受荷情况下含水合物沉积物胶结结构的破坏情况，基于扰动状态概念建立了描述胶结结构破坏的全过程的含水合物土的应力—应变模型。

近年来，含水合物沉积物的应力应变本构关系模型发展迅速，新的模型变得越来越完备，除了能够考虑一般土体的主要力学行为特征外，还能模拟水合物含量及赋存模式的影响，但代价是模型变得越来越复杂，包含了许多物理意义不明确的且难以获取的模型参数，从而导致了这些新模型很难应用于解决实际工程问题。

4 水合物开采土力学多场耦合过程的数值模拟方法

数值模拟方法为水合物储层的开采潜力评价及开采场地安全性预测提供了一条重要途径[126-127]。然而，针对水合物开采的数值模拟需要求解复杂的非线性多过程耦合问题。由于存在强耦合作用，常常会出现物理上的奇异和数值计算上的不稳定。另外，水合物储层是一种多相孔隙介质，其物理和化学性质较强地依赖于水合物含量、温度、压力等状态变量，强烈的非线性特征使得求解的困难程度大为增加。因此，发展针对水合物开采模拟的稳定高效数值计算方法目前依然是重要的科学难题。

针对于开采潜力评价的数值分析（针对于产能预测），目前国际上已研发出一系列较为成熟的水合物开采数值模拟器，如美国劳伦斯伯克利国家实验室的TOUGH+HYDRATE[128]、日本石油工程公司和东京大学合作开发的MH-21HYDRES[129]和加拿大CMG公司的CMG-STARS[130]等。这些模拟器主要是基于孔隙介质的质量和能量守恒方程，考虑了相变、传热、液气运移和组分扩散过程，忽略了土体变形/破坏对其他物理场的影响。笔者主要针对开采过程中含水合物储层的稳定性和井壁稳定性的研究进行论述。宁伏龙[3]研究了水合物地层钻探过程中出现的井壁稳定性问题，指出实际钻井过程中井壁围岩地层的孔压和应力分布情况是流体流动和固体变形共同作用的结果。程远方等[131]基于有限元法，采用流固耦合模型对降压开采过程中储层的应力状态以及物性参数的动态变化规律进行了研究。然而，该模型忽略了温度场的变化，并采用了解耦方法对渗流场和变形场进行了单向顺次求解。Rutqvist等[132-133]采用单向耦合的方法将TOUGH+HYDRATE和FLAC结合起来，采用前者计算温度、孔压、水合物含量等，模拟水合物储层的传热、渗流以及水合物分解过程，采用后者计算应力应变，模拟分析储层的变形和稳定性问题。两者之间通过专门设计的通信模块来实现数据的单向传递。与此思路相同，Qiu等[134]通过单向耦合MH21-HYDRES和3D FE地质力学模拟器分析评价了日本Nankai海槽水合物藏降压开采过程中的井筒完整性。

为了提高计算精度，Kim等[135]进一步发展了TOUGH+HYDRATE和FLAC的交错耦合模型，模拟水合物分解过程中的流体渗流和岩土体变形之间的耦合作用。与单向耦合法相比，交错耦合法（即双向耦合法）将信息由单向传递变为双向传递，使计算过程能够符合不同物理场之间是相互影响的这一基本事实。为了克服FLAC作为商业软件在源代码和内存共享等方面的局限，Queiruga等[136]开发了有限元程序Millstone以替换FLAC来计算储层的力学行为，并提出了新的不同网格间的数据传递方法，有效提高了计算效率。需要指出的是，单向耦合法和交错耦合法同属于解耦算法，虽然在计算方法上相对容易实现，可以利用已有的较为成熟的热传导、流体渗流和固体变形等单场求解程序，但未能将各个物理场放在同一地位进行考虑，理论上不能称为严格的耦合算法。更重要的是，在处理像水合物开采这种瞬态非线性问题时，要统一协调不同软件之间的时域积分和非线性迭代，解耦合法或多或少存在一些“相延迟”问题，通常不得不采用较小的时间步长，很难把握计算过程的收敛性，甚至会出现无法模拟某些重要的物理现象的情况，例如Biot流固耦合问题中经典的Mandel-Cryer效应。为此，不少学者开始尝试将全耦合方法引入到水合物开采的多场耦合模拟中。

Salehabadi等[137]利用力学商业软件ABAQUS结合Heriot-Watt大学开发的水合物分解模拟软件HWHYD，研究了水合物分解条件下套管周围应力场的变化，分析了水合物分解对井眼及套管的稳定性的影响。该研究中对流固耦合问题采用了全耦合方法，但对热流耦合和热固耦合问题只进行了单向耦合求解，即只考虑传热对渗流和变形的作用，而忽略了渗流和变形对温度场的影响。Klar等[138]推导了气液两相渗流方程和动量平衡方程的显式有限差分格式，采用了理想弹塑性模型来描述固体骨架的变形，最终通过对FLAC软件的二次开发实现了解耦求解。不久前，Klar等[138]继续补充了热传导方程，进一步考虑了温度场的变化对开采的影响。Uchida等[139]首次将颗粒运移问题耦合到水合物开采过程中，使得该程序可模拟储层出砂过程。模型中考虑了沉积物变形、颗粒分离、颗粒运移、颗粒沉降等过程，通过数值分析评估了出砂相关参数对固体骨架体变和水合物含量的敏感性。最近，Zhou等[140]在该程序中嵌入了修正的临界状态模型来模拟水合物储层的力学行为，并对日本Nankai海槽东部海域的水合物开采过程进行了分析。由于采用了显示算法，该数值模型虽然考虑了众多参数的非线性特征，但并未引入迭代过程，而是直接采用上一时间步的结果对当前时间步的参数进行更新，其计算结果的准确性及其对时间步长的依赖性仍值得进一步研究。

除此之外，Kimoto等[101,141]、Raz[142]、吴二林等[143]也自主发展了水合物储层开采模拟的多场耦合数学模型，并基于实测数据开展了地层力学参数演化及沉降变形预测的研究。与之前类似，上述研究并不完善，仍有部分计算过程采用了解耦策略，忽略或淡化了某些耦合作用，这与实际发生的强耦合物理过程是不相符合的。此外，当采用解耦合方法时，往往是直接调用现有的单场软件。不同商业软件对输入和输出数据的开放程度有很大的不确定性。受到可选用的单场软件的限制，不同物理场的耦合求解经常不得不采用不同的数值方法，如渗流场采用有限体积法、位移场采用有限元法、温度场采用有限差分法。如何在不同类型的计算软件之间协调网格剖分、统一边界条件、设计数据传递渠道，对于解耦方法来说目前仍然困难重重。

5 发展趋势与未来的挑战

5.1 物理力学特性测试方法和技术

目前发展出的多功能探测技术主要以定性为主，围绕沉积物组构的精细定量探测手段相对缺乏。CT、声波和核磁共振技术可以对土组构进行无损、定量的探测，但由于在信号解析上缺乏成熟的理论模型，这些技术在应用中仍存在困难。随着CT、声波技术与核磁共振技术的不断发展，沉积物组构探测会逐渐走向定量化，并将与其他探测技术以及沉积物工程力学特性测试相结合，这是一种发展趋势。在沉积物工程力学特性测试方面，试样制备—组构探测—力学测试一体化、多种测试功能集成化、面向原位取样的测试方法、模拟原位开采工况的大型模型试验等将成为未来的重要研究方向。另外，由于含水合物沉积物的结构易损性，沉积物组构探测与工程特性的原位和船载测试技术会得到进一步发展，除了传统的海洋原位土工试验、地球物理、声呐等技术会进一步完善外，新的测试技术（如水下光纤、机器人等）将越来越多地被应用。同时，室内与原位测试会更加紧密结合，相互取长补短。

5.2 含水合物沉积物的力学特性研究

在含水合物沉积物的力学特性方面，目前的试验数据比较零散且各类数据相互独立，缺乏针对典型水合物储层介质的系统性试验数据。这些问题对全面、准确地认识含水合物沉积物的力学特性及其变化规律是非常不利的。解决这一问题需要针对一些典型水合物储层，系统开展相关沉积物的水力特性、力学特性以及相变特性的研究，建立完备的沉积物储层物理力学特性数据库。另外，随着水合物分解/生成，水合物含量及其赋存模式在不断变化中，使得水合物储层介质具有明显的动态结构性，因此，厘清水合物储层介质的组构特性及其关键影响因素、揭示组构演化对沉积物力学特性的影响机理和规律也将成为未来有关水合物研究的重点方向。

5.3 含水合物沉积物的物理力学模型

开采过程中水合物分解会导致储层介质会发生变形沉降，分析水合物储层沉降对开采系统（特别是海底管线）所造成的影响，需要进一步发展能够描述水合物开采过程的储层介质多场、多过程耦合理论，建立能有效模拟在压力、温度、化学等复杂荷载作用下含水合物沉积物力学行为的实用本构模型。由于含水合物沉积物具有显著的相变性和动态结构性（其结构随着水合物分解、流泥流沙过程而不断发生变化），需要对水合物沉积物组构进行定量表征，并在构建强度准则和本构模型时充分考虑组构演化的影响。尽管开采场地的沉降问题需要考虑，但相对而言，储层破坏（如海底滑坡）、井壁失稳等破坏问题更加突出。针对这些问题，需要发展相应的计算模型和分析方法，包括：井壁稳定性分析、海底稳定性分析、锚固基础抗拔承载力计算等，因此，含水合物沉积物的强度特性及其影响因素和变化规律将是未来需要重点关注的问题。

5.4 水合物开采土力学多场耦合数值模拟

水合物开采中产能预测以及开采场地力学安全性分析涉及到复杂的非线性问题求解。随着水合物储层特征刻画的精细化、沉积物多过程耦合理论的进一步发展以及数值方法和计算技术的不断进步，这一方向的发展趋势具有以下两种特征：①从分析尺度上看，逐渐实现从储层介质的微细结构演变到宏观场地尺度的变形破坏的跨尺度模拟分析；②从数值求解方法看，正在从传统的单向耦合、交错耦合等解耦算法转到全耦合算法上，逐渐实现水合物开采边值问题数值解法的根本性转变。

5.5 含水合物沉积物力学研究面临的挑战

经过近20多年来的发展，有关水合物开采的理论、方法和技术均取得了显著的进展，但从岩土力学的角度看，目前仍然面临诸多挑战。

5.5.1 大尺度含水合物沉积物试样的人工制备

开展室内试验和物理模型试验需要在实验室中制备大尺度均匀的含水合物沉积物试样。由于含水合物沉积物物理力学特性强烈依赖于沉积物的结构，因此所制备试样必须能够真实反映原位含水合物沉积物的结构特征，特别是孔隙水合物的赋存模式应与原位一致。迄今为止，这仍然是没有完全解决的难题。

5.5.2 含水合物沉积物微细观组构的精细探测与定量表征

含水合物沉积物具有显著的结构性，而且水合物分解/生成使得沉积物的微细观结构处于动态演化中。因此，厘清含水合物沉积物的物理力学特性的变化规律必须能够对其微细观结构进行精细探测与定量表征，这正是目前有关含水合物沉积物力学研究中面临的挑战之一。

5.5.3 组构变化对含水合物沉积物力学特性的影响机理与规律

水合物分解会导致其胶结作用破坏，从而引起沉积物的强度和变形特性的变化。因此，阐明水合物含量、水合物赋存模式等组构因素变化对沉积物宏观力学特性的影响机理与规律，是有效评价水合物储层力学稳定性的重要前提。由于缺乏相应的试验数据，这一关键问题目前尚未解决。

5.5.4 开采扰动下流沙发生条件及其对储层稳定性的影响机制

实现水合物高效开采需要对开采扰动下流（出）沙过程进行有效控制，因此需要阐明流沙过程的发生机理和触发条件。另外，流沙过程也会导致储层或开采井壁破坏，从而影响到开采过程。目前对水合物开采过程的流沙问题认识有限，研究手段单一（主要基于离散元法），没有建立有效的流沙判别准则，缺乏考虑流沙影响的储层稳定性分析方法。

5.5.5 开采扰动下多相多组分含水合物储层的多过程耦合问题

含水合物储层是一种多相多组分孔隙介质，在开采扰动下，会发生热传导、气—液流动、物质扩散、细粒流失（流泥流沙）、骨架变形、水合物生成/分解等耦合的多物理化学过程。虽然传统理论很好模拟热传导—渗流/扩散—骨架变形等耦合过程，但对存在相变等化学过程情况仍存在明显的局限性。另外，现有的多孔介质连续介质力学模型没有充分考虑流沙问题的影响。

5.5.6 水合物分解/生成条件下含水合物沉积物的本构响应

在水合物分解或生成条件下，含水合物沉积物的内部组构可以在胶结模式—填充模式—持力体模式之间相互转化。目前囿于试验技术的限制，对这类组构不断变化的土体的本构响应特征缺乏了解。已有的本构模型无法准确模拟组构变化的影响，而近年来发展起来的能考虑赋存模式影响的沉积物本构关系模型不仅结构复杂，而且涉及到模式判别问题，这限制了新模型的应用。因此，发展能有效模拟水合物分解/生成条件下含水合物沉积物本构关系模型仍然是目前亟需解决的难题。

5.5.7 有关水合物开采过程数学模型的适定性问题

目前模拟水合物开采过程的数学模型都是通过对原本只适用于描述单一过程的数学模型的简单组装而得到，如热传导采用的是经典的傅里叶法则，渗流过程采用的是达西定律，扩散过程则是菲克斯第一法则等。为了有效模拟水合物开采过程，必须首先回答该类数学模型解的存在性、稳定性和唯一性问题。

5.5.8 高效稳定的水合物开采过程的数值模型的全耦合解法

有关水合物开采过程模拟所采用的解耦或半解耦解法与实际物理过程严重不符，不仅精度低而且难以控制。因此，发展全耦合计算方法将是一种理想选择。但是，采用全耦合计算方法除了需要从头开始进行软件设计开发，还需要发展高效稳定的全耦合数值算法，这是一项亟待解决的重大科学问题。