天然气水合物热分解引起多场响应的特征时间分析1)

张旭辉 鲁晓兵

(中国科学院力学研究所流固耦合系统力学重点实验室，北京 100190)

(中国科学院大学工程科学学院，北京 100049)

天然气水合物（以下简称水合物）是资源量丰富的高效清洁能源。水合物赋存区广泛分布在海洋大陆边缘、陆地永久冻土带和部分内陆深水湖泊。据估计，水合物中甲烷资源量约为2.0×1016m3，其资源量约是当前已探明化石燃料（煤、石油和天然气）总量的两倍[1-2]。海域水合物资源量主要分布在水深800 m以上、海底以下数十米到百米的沉积层中[3-5]。我国海域水合物勘探的远景资源量高达800亿吨油当量，中国地质调查局组织实施的两轮试验性开采，深化了对水合物开采工程实践的认识[4,6]。

要实现水合物的商业化开采，必须解决两个方面的问题：灾害控制及环境保护、经济高效的开采方法[7]。尤其在全球气候变暖环境下，水合物开采相关的地质灾害和环境安全问题引起了世界各国政府和学者的极大关注，对水合物开采均采取了谨慎态度。我国南海的水合物赋存于浅层未固结岩土层中，致其开采难度大，亟需回答水合物开采会不会使得原来稳定的海水-上覆盖层-储层-下伏层四层系统内发生失稳、什么条件失稳、如何防控等问题[8]。

水合物开采与常规油气的根本区别在于前者涉及相变，生成天然气和水，一方面使含水合物地层中原有土颗粒之间的胶结和微结构受到大幅扰动和破坏；另一方面随天然气超过当地饱和度而产生自由气，使原地层成为一种强度很低的非饱和欠固结土层[9-10]。相变区域土体强度的急剧衰减，尤其在热分解相变条件下，若孔隙流体压力的缓慢耗散会导致超静孔隙压力上升与有效应力降低，以及地层承载性能很弱，则可能演化为海底滑坡、海床隆起或沉陷、盖层的裂隙发展[8,11-12]。这些灾害在四层系统环境一般不是单独发生，常常以灾害链的形式发生，如滑塌转化为泥流、浊流，甚至剧烈的海流，冲击海洋平台、油气管道及海底电缆等设施；地层开裂、滑塌将水合物储层与海水连通，海水入侵引起大面积的储层中水合物分解，导致甲烷气体泄漏，带来海洋生态环境以及全球气候平衡等的破坏灾难[8]。

以水合物热分解后的岩土体中的层裂和喷发为案例[11-12]。岩土体中的层裂是地层因高压气体或水而导致断裂，断裂的地层被水和气填充形成。压力足够大且能量足够大时则产生喷发。喷发是水、气和岩土颗粒混合物从地层中类似火山的强烈地喷出。层裂和喷发这两种破坏皆可能导致水合物大面积暴露于海水，进而产生大面积的水合物分解。人们在西伯利亚发现了多处直径20～100 m、深度约50 m的巨大天坑，认为是平均温度逐渐升高，水合物分解生成高压气区，且冻土区退化，上覆强度减弱而形成[13-15]。在百慕大地区发现的数量巨大的海底坑，呈椭圆形、水平尺度300～1 000 m、深度达30 m、侧边坡度达50°。人们认为这是由于从下部的油气藏通过断层和裂隙迁移到浅层来的天然气在上一次冰期形成水合物；在冰退期，温度升高导致水合物分解产生大量气体，进而引起气体伴随沉积物的突然喷发，形成巨坑[16-17]。一旦溢出到海床面上，甲烷水合物或气泡可能溶解于水中或上升到海面，如果大量的水合物释放到海水中，则可在约500 m水深处形成大而集中的泡状羽流，进而导致海面上的喷发[18-19]。

水合物热分解引起地层软化与破坏的关键科学问题归根结底是深海浅层含相变的多过程、多时间空间尺度的岩土力学问题，即多孔介质内的传热、水合物相变的演化、孔隙内渗流和地层应力传播，甚至岩土破坏以及固体破坏向多相流动转化的一系列过程，是一个包含热、物理化学和力学的多个物理效应问题[7,20-22]。

这些过程均遵循自身特有的特征时间，可能具有多个时间尺度，同时相变导致的土体物性的不均匀、不连续，使得地层在空间上具有多个尺度。对这几个物理过程组合在一起构成的极其复杂的问题，不论是理论分析，还是现场实验、模型实验的分析难度均极大。首先要判断这些效应是否耦合，即紧密配合并产生相互作用？还是相互之间具有时间尺度的量级差别，从而可以进行解耦分析？解耦分析的思想是通过量纲分析得到各个物理效应的特征时间，问题中若多个物理场的自身特征时间之间的差异具有2～3个数量级以上时，进行先慢后快的简化分析，更明确地体现出问题中的主次和各个物理过程之间的关联，找到整个物理过程的主要控制因素，弄清整个物理过程的因果关系，便于简化实验设计和理论分析，比如，先分析慢的过程，以慢过程的结果作为稳态或者准静态条件，分析快的过程，依次类推，这与渐近方法中的多重尺度分析方法是一致的。同理，在建立数值方法时，各物理场计算的时间步长可以根据特征时间差异给出，提高计算效率且减少系统误差[2,23]。

1 基于量纲分析的特征时间推导

水合物沉积层因水合物热分解引起的软化和破坏的物理过程可以描述为：随着热量的不断输入，一方面提供水合物相变热量，促使水合物相变，一方面提高沉积物骨架的温度以及分解产生的水气的温度，分解区域逐步向前扩展，形成水合物沉积层中水合物热分解演化过程；水合物热分解产生的液气封闭在分解区域内，可以形成超静孔隙压力；水合物含量的减少也增加了孔隙流动的有效空间，使得分解区域的渗透率提高，水合物热分解后沉积物孔隙内的压差推动气体和水在沉积层分解区域中渗流；水合物沉积层因水合物热分解后处于欠固结状态，超静孔隙压力引起沉积层和覆盖层的应力状态重新分布，引起沉积层和覆盖层的变形。随着沉积层中水合物分解区域的扩展，水合物分解相变产生的气体压力可以克服上覆阻力时，地层破坏发生，如层裂、剧烈的喷发等[21]。 基于上述物理过程，以一维问题为例，进行特征时间的分析。

由于含水合物沉积物是一种多相材料，包括土、水、气、水合物。各组分的体积分数之和恒定，即

其中，ε表示体积分数，下标h，w，g和s分别表示水合物、水、气和骨架。当然，初始时刻各相体积分数也满足εs0+εw0+εg0+εh0=1 。

简便起见，基于混合物理论对多相体的热物性（比热、密度、热传导系数）进行平均化处理，比热C和热传导系数λ按照各组分含量平均加和模式进行计算

假定水、骨架、水合物的密度为常数；不考虑水合物颗粒与细土颗粒的运移；甲烷气体是理想气体；气体与水在骨架渗流暂认为满足达西定律，且不考虑多孔介质中流体流动产生的对流传热的影响。

1.1 热传导特征时间

热传导特征时间考虑水合物地层边界一端恒定加热引起的一维热传导问题[23-24]

地层的特征长度为L，则热传导特征时间

一般地，热传导特征时间可以近似[23-25]

1.2 相变特征时间

甲烷水合物稳定存在的相平衡条件：pe=pe(T)。pe的数值由拟合经验公式计算[21]

其中p0，T0，a，b，c分别为 1 atm(=101.325 kPa)，273.15 K，9.34，37.3，1.48。

当温度超过相平衡温度时，水合物发生分解相变，多孔介质内单位体积水合物分解速率满足[21,26]

其中kd为多孔介质内甲烷水合物分解速率常数，kd=4.4×10-16kg/(m2·Pa·s)；As取 决 于 有 效√流体孔隙度和渗透率K两个物理量，即为As=为有效流体孔隙度，ε=(1-εs)(1-εh)；f为实验或者实际下甲烷气体逸度；fe为三相平衡条件下甲烷气体逸度。

一般地，用地层孔隙压力pg和甲烷水合物相平衡压力pe分别代替f和fe，保持恒定条件下分解，且在相变过程对有效孔隙度和渗透率平均化处理，对式（9）两边积分可得

从式（10）可以看出，相变特征时间不同于热传导特征时间，它与地层的特征长度L无关，而取决于热分解条件下地层孔隙内分解动力学演化过程。

1.3 渗流的特征时间

渗流的特征时间是压力扰动从多孔介质的一端传递到另一端的时间，是由地层本身的物理性质所决定的，而不是直接通过达西定律计算流体质点从一端流动到另一端的时间[27-28]。考虑简单的一维单相渗流问题，基于孔祥言[29]提出的渗流方程，推导出渗流特征时间。 根据质量守恒和达西定律，并且不考虑流体可压缩性和重力作用，渗流方程为

其中，p为流体压力，µ为流体黏性，ε为有效流体孔隙度，渗透率为K，ct为地层的压缩系数。

地层的特征长度为L，则渗流特征时间

与热传导类似，同为抛物线型偏微分方程，将特征时间可以近似简化为

1.4 应力重分布特征时间

应力重分布的特征时间为水合物热分解相变引起的局部应力应变扰动以弹性波的形式向周围传播特征尺度所需的时间。弹性波的传播速度可以用表示。水合物地层的应力状态受到弹性波的扰动以后恢复到准静态的应力分布的时间，大致与弹性波传播到边界的时间属于同一数量级[23-24]。地层的特征尺度为L，则特征时间可以近似表达为

若应力重分布后，局部地层达到破坏临界条件，则会发生地层层裂、气体喷发等强度破坏，在考虑更为复杂的重力效应时，可能发生塌陷、滑坡等[12, 27-28]。

2 特征时间数量级比较及解耦分析方法

选取特征长度为1 m，根据水合物沉积层的热学参数和力学参数进行估算[30-31]，且考虑初始土骨架和水合物的体积分数为0.6和0.4。

取骨架的热学性质为代表，可以估计热传导特征时间

考虑两相问题，水合物沉积物分解的比表面积

温度为 3 03.15K ，水合物的相平衡压力数值pe=9.04×107Pa ，估计水合物分解后的气体压力pg=8.76×107Pa。则水合物分解相变特征时间

取水合物沉积层的弹性模量E=1×108Pa ，泊松比v取为0.3，沉积层的侧限压缩模量为Es=压缩系数近似为ct=，则气体相渗流特征时间

水合物地层的应力重分布特征时间，按照等同于弹性波传播的特征时间

通过对热传导特征时间tc、相变特征时间td、渗流特征时间ts和应力重分布特征时间te进行估算，可以得到在水合物沉积物分解区域内四个特征时间的数量级分别为： 0.7×106s ， 4×104s ，1.9×102s ， 5×10-3s 。

从这四个物理过程的特征时间之比来看，应力波传播最快，热传导过程最慢。特别指出的是，相变的发生除了到达相平衡温度，水合物分解需要吸收大量的热量，还需要满足供入能量大于等于水合物相变潜热，即CT≥ΔH，因此将热传导和相变两个物理效应作为含相变的热传导效应。渗流过程比热传导过程快约3个数量级，应力传播又比渗流过程快5个数量级。因此，水合物热分解物理过程可以解耦分析[23]。

该问题可以按以下解耦的顺序进行依次分析。

（1）根据沉积层中含水合物相变的热传导过程分析各个时刻的温度场分布，物理上水合物分解相变近似以分解阵面演化推进，得到多相多阵面演化的自相似解，确定分解区域和未分解区域[30]。

（2）根据稳态或者准稳态的温度场分布，结合气液渗流分析得到孔隙压力分布。沉积层中水合物热分解后的渗流包括两个部分，一个是水合物沉积物分解区域的渗流，由于水合物热分解后，沉积物孔隙渗透率提高，渗流属于相对快的过程；一个水合物沉积物未分解区域，这部分沉积层的渗透率很低，分解产生的气体和水渗流十分缓慢，渗流效应的影响又是一个很慢的过程，几乎可以认为不发生渗流，孔隙流体压力可根据封闭体系内的气体状态方程和已求解的温度场得到[12, 32-33]。

（3）基于以上两个过程（水合物分解范围的大小和压力分布）和地层的力学参数，根据有效应力原理进行沉积层和覆盖层应力应变分析，判断是否达到地层临界破坏条件（如层裂、喷发）等[12, 32]。

综上所述，解耦分析方法已得到热分解相变引起多界面演化、地层层裂、喷发和滑塌等的实验结果的验证。特别指出，通过对含水合物相变源项的质量守恒、动量守恒、能量守恒以及物性方程和初边值条件的数学表述进行无量纲化，得到相似的特征时间，进而通过多重时间尺度渐近法进行解耦分析[34-35]。

这里需要特别指出，不同岩土介质的传热和渗流的物性不同，对于非常低渗透的岩土介质，可能会出现渗流比传热慢的情况。在考虑不同岩土储层中这类多场响应问题时，需要首先获取相应参数，具体问题具体分析，判定孰快孰慢以及是否耦合。

3 结论

对水合物热分解引起多场响应进行了物理上的阐述，并通过物理效应分析和量纲分析的方法得到了传热、相变、渗流和应力重分布的特征时间，它们分别是

将这些特征时间赋予水合物沉积层的物理和力学参数，得到特征时间分别约为 0.7×106s ，4×104s， 1.9×102s ， 5×10-3s 。相变与热传导可以看作是含相变热传导的阵面传播效应，且含相变热传导、渗流与地层应力重分布特征时间相差2个数量级以上，可以进行解耦分析，并得到了系列物理模拟实验的验证，同时物理分析得到的特征时间与从数学表述得到的特征时间是一致的。

对于多效应问题，首先判断各个效应之间是否耦合，可以从特征时间的比较上出发，判定是否紧密配合并相互作用。若可以解耦分析，对于问题的解决，无论从模型实验还是从理论求解、数值模拟上均进行了有效的简化，均可以方便深入认识物理机制，并且工程实用。