引言

近年来，储层工程实践经验^[1]和科学研究^[2-5]表明，近远场岩土力学效应（如地表隆起、沉陷，井内位移等），能够有效反映储层的体积膨胀与压缩、岩体开裂、流体迁移、流体泄漏等储层动态（reservoir behavior）。地表变形监测结合岩土力学数值模拟，作为洞察储层动态的重要方法^[6-9]之一，受到广泛关注和持续研究^{[4, 10-12]}。尤其是针对阿尔及利亚Insalah咸水层二氧化碳埋存工程，所监测到的KB502井附近迥异地表变形的研究^[13-19]，颇具代表性。众多研究团队使用不同方法，尝试拟合和解释这一现象，取得很多有价值的成果，但同时也暴露出一些问题。比如，至今鲜有通过地表变形监测信息，反演储层应力、应变及孔隙压力分布等变化的尝试，并且对地表变形监测结果的历史拟合效果，也远不能令人满意^{[10, 20-21]}。

储层注气工程的目标储层埋深一般在500~5 000 m，自储层到地表即包含成十上百不同岩性的地层。地表变形数值模拟，模型中须包含储层下方相当范围的台基部分（达到回弹效应可以忽略的深度）直至地表的所有岩层。对于尺寸如此庞大的模型，若如传统储层流动模拟一样，对整个模型细划网格，模拟运行将笨重不堪。现有研究^{[2, 13-14]}一般通过细划储层近场网格，粗划远场范围，甚至将储层上方至地表的成十上百不同岩性参数的凹凸地层，简化为数层水平均质地层，未能取得理想的计算效果。Morris等^[14]建议当盖层参数可以得到时，模型应包含尽量详细的盖层参数，强调了建立复杂地质模型的必要性。Eiken等^[10]也认为，地质模型的建立和完善还有很大的提升空间。

本文通过自主研发大规模网格离散方法，建立了计入详细地层拓扑形状和岩性参数，并能够随注气扩散过程调整网格粗细分布的储层全场地质模型(full-field reservoir model)，利用Eclipse与FLAC$^{\rm{3D}}$实现储层流动与岩石力学单向耦合，模拟了某大型咸水层二氧化碳埋存工程注气一定时间后的地表变形，研究了常规模拟流程在远场微小变形分析问题中出现计算精度不足等局限性的原因，并给出了解决方案，最后探讨了地质模型中所包含的储层下伏岩层的厚度对模拟精度的影响。

1 建模方法介绍

本文建模的主要思路如图 1所示，利用自主研发的大规模网格化离散方法，将储层流动模拟软件Eclipse所建立的储层近场地质模型转换并加入远场岩层地质信息，使之延伸为FLAC$^{\rm{3D}}$全场地质模型；并转换Eclipse储层流动模拟结果，使之为FLAC$^{\rm{3D}}$岩石力学模拟所用。进而实现储层流动-岩石力学效应单向耦合分析。若进一步改进该方法，使之不仅能实现上述过程，而且能实现虚线所示的逆过程，则能够利用Eclipse与FLAC$^{\rm{3D}}$两软件实现方便高效的储层流动与岩石力学效应的全耦合分析。

1.1 建模流程

本文的建模流程如图 2所示。

（1）设定模型所含网格层数（工程场地深度方向），以及模型每层网格核心区域和放射区域的节点数目。按照图 3a所示的规律分别为每层网格核心区域和放射区域的节点赋予节点编号，为每个单元赋予单元编号，并指定组成该单元所需要节点的编号。

(2) 设定模型实际长宽尺寸，以及模型核心区域占整个模型长宽的比例。通过数值离散，按照一定的顺序为各层网格核心区域和放射区域的每个节点分别赋予横坐标(x)和纵坐标(y)。

(3) 设定每层岩层所要划分的网格层数，根据Eclipse储层近场地质模型以及其余远场岩层的地质信息，利用二元线性或非线性插值(根据节点x、y坐标对节点z坐标进行线性或非线性插值)，为每层网格核心区域和放射区域各节点分别赋予竖坐标(z)。至此可以输出适应FLAC$^{\rm{3D}}$格式的建模文件，导入FLAC$^{\rm{3D}}$后可建立计入详细地层拓扑形状和岩性参数的储层全场地质模型。

通过以上步骤所建立的简单模型的例子如图 3b所示，均质凹凸起伏岩体，共划分两层网格，每层网格的核心区域划分为4×4个单元，放射区域共划分为4×3×4个单元，核心区域节点均布，放射区域在环向上节点均布，在径向上节点按1.5的比例发散分布。

(4) 计算FLAC$^{\rm{3D}}$模型各单元中心点坐标。

(5) 根据Eclipse储层流动分析结果，计算得到对应于Eclipse模型各单元由注气渗流过程引起的孔隙压力净增量数据，并计算Eclipse模型各单元中心点坐标，输出对应于Eclipse模型各单元中心点坐标的孔隙压力净增量数据。

(6) 根据FLAC$^{\rm{3D}}$模型各单元对应的中心点坐标和Eclipse模型各单元中心点坐标及其对应的孔隙压力净增量数据，通过插值方法，获得FLAC$^{\rm{3D}}$模型中心点坐标对应的孔隙压力净增量，由此得到FLAC$^{\rm{3D}}$模型各单元对应的孔隙压力净增量数据。

1.2 本文建模方法的特点

(1) 能够对结构复杂储层场地建立全场地质模型并划分网格，较精确地计入实际地层拓扑形状，而不必简化成直线和平面。

(2) 所构建的三维网格主要由正交单元垂直堆叠而成，网格划分情况对数值分析结果的影响较小。

(3) 在获取地质信息后，可快速简便地构建储层全场地质模型，且可以快速重构现有地质模型，是一种对现有地质模型进行补充和改进的良好工具。

(4) 在保持确定的地层拓扑形状的条件下，可快速地自由调整网格划分情况，具有良好的用户可控性。

(5) 针对储层注气工程的特点，将模型设计为由核心区域和放射区域两部分组成。方便使模型的核心区域随着注气过程精细地模拟注入气体在储层内的扩散动态，在保证计算精度的前提下，减少了模型总的节点数和单元数，实现兼顾计算精度与计算效率的数值分析。提供了在耦合分析中动态聚焦注气扩散区域的初步探索。

2 储层注气工程地表变形数值模拟 2.1 计算原理

根据多孔介质力学有效应力原理，有效应力等于总应力减去孔隙压力。具体到本文研究的问题，储层注气前后(如图 4所示)，有

$ \sigma _{0ij}^{\rm{e}} = \sigma _{0ij}^{\rm{t}} - p_0 $

(1)

$ \sigma _{1ij}^{\rm{e}} = \sigma _{1ij}^{\rm{t}} - p_1 $

(2)

储层注气过程中，上覆岩层自重不变，因此储层岩体总应力不变，即

$ \sigma _{1ij}^{\rm{t}}=\sigma _{0ij}^{\rm{t}} $

(3)

由于注入储层气体的扩散，导致储层孔隙压力增大，即

$ \Delta p = p_1 - p_0 > 0 $

(4)

综合以上各式，有

$ \Delta \sigma _{ij}^{\rm{e}} = \sigma _{1ij}^{\rm{e}} - \sigma _{0ij}^{\rm{e}} = - \Delta p < 0 $

(5)

式中：

$\sigma _{0ij}^{\rm{e}}$─注气前储层岩体有效应力，MPa；

$\sigma _{0ij}^{\rm{t}}$─注气前储层岩体总应力，MPa；

$p_0$─注气前储层岩体孔隙压力，MPa；

$\sigma _{1ij}^{\rm{e}}$─注气后储层岩体有效应力，MPa；

$\sigma _{1ij}^{\rm{t}}$─注气后储层岩体总应力，MPa；

$p_1$─注气后储层岩体孔隙压力，MPa；

$\Delta p$─注气过程储层孔隙压力增量，MPa；

$\Delta \sigma _{ij}^{\rm{e}}$─注气过程储层有效应力增量，MPa。

即由注气引起的储层岩体有效应力的减少值等于储层岩体孔隙压力的增加值。注气引起了储层岩体有效应力的减少，导致储层岩体膨胀，进而导致地表隆升等远场岩石力学效应。

根据以上分析，利用自主研发的建模方法建立FLAC$^{\rm{3D}}$格式的储层全场地质模型，且获得模型各单元对应的孔隙压力净增量数据后，将模型中各单元初始有效应力减去其对应的孔隙压力增量，计算平衡后，即可得到由注气引起的全场地质模型的近远场变形及位移情况。

2.2 计算模型

根据某大型咸水层二氧化碳埋存工程的实际地质信息资料，及其Eclipse储层流动模型和Eclipse模拟分析结果，利用自主研发的建模方法，建立计入详细岩层拓扑形状和物理力学参数的全场地质模型，并得到该模型各个单元所对应的孔隙压力净增量数据。

所用Eclipse储层流动分析结果为，模拟该工程自注气起始之日，每天向储层注入301 t CO$_2$，持续至第1 096 d时的储层状态。

2.2.1 模型尺寸及网格划分情况

如图 5，该咸水层二氧化碳埋存工程的目标储层埋深平均约为3 000 m，储层厚度约为200 m。图 6为该工程场地地表等高线图，地势偏低的蓝色部分为一条河流流经区域。

本文所建立的储层全场地质模型，其水平方向尺寸为14 418 m×14 418 m。由于测井信息只到二氧化碳埋存目标储层的下一层岩层，因此模型深度方向取储层下一层岩层直至地表的范围(本文第2.2.2节和第3.2节介绍了模型中所包含的下伏岩层的厚度对地表变形模拟结果的影响研究)，平均约为3 400 m。模型中共包含25种不同物理力学参数的岩层，划分为68层网格，自地表至储层方向各层网格分别编号为第1层，第2层…第68层。目标储层的编号是第52~59层。模型网格划分如图 7所示，核心区域(6 488 m×6 488 m)单元划分较密，放射区域距离核心越远处单元划分越稀疏。

对比图 6和图 7，可以看出，自主研发的大规模网格化离散方法通过自适应插值较好地处理了岩层拓扑形状。

2.2.2 物理力学参数

地质模型物理力学参数均按实测地质信息赋值。地质信息来自3D地震测量、岩石物理岩芯数据及综合的测井数据等。按照由场地表面到储层方向的顺序，各岩层弹性模量、泊松比和密度的取值如表 1所示。

2.2.3 本构模型

借鉴InSalah项目的研究经验^{[9, 21]}，使用线弹性本构模型分析该问题。

2.2.4 边界条件

模型底部边界施加x、y、z 3方向位移约束，对四周边界分别施加法向位移约束。模型上部边界自由。

计算过程中当最大不平衡力与典型内力比值小于1.0×10$^{-7}$时，即认为达到平衡。

2.3 研究问题及方法 2.3.1 地表变形模拟与局部误差的消除

岩土工程数值模拟，往往采取先使模型在自重作用下计算平衡，获得初始地应力状态，再模拟具体工程问题的计算流程。本文研究过程中发现，使用这一常规计算流程(下文方案1) 模拟储层注气引起的地表变形时，在贴近地表的数层岩层处产生了较大的位移误差，难以取得理想的地表变形模拟结果。为解释这一现象和提供更适合的分析方法，采用3种不同的数值模拟方案分析由储层注气引起的场地地表变形。

由于采用完全弹性本构模型，根据弹性力学叠加原理，3个方案理论上应取得相同的变形计算结果。

方案1：(1) 将建模文件导入FLAC$^{\rm{3D}}$中，按照表 1所列信息设置岩层物理力学参数，设置边界条件；(2) 使模型在自重作用下计算平衡，并将位移场和速度场清零，获得场地初始地应力场；(3) 读入模型各单元的孔隙压力净增量，根据有效应力原理，编写FLAC$^{\rm{3D}}$ fish命令流，将各单元3个方向正应力减去该单元对应的孔隙压力净增量，并赋值给该单元对应方向的正应力，重新计算平衡，得到由注气引起的变形。

方案2：(1) 将建模文件导入FLAC$^{\rm{3D}}$中，将模型地表部分3层网格的密度设置为零，其余按照表 1所列信息设置岩层物理力学参数，设置边界条件；(2)、(3) 与方案1相同。

方案3：(1) 与方案1相同；(2) 设置重力加速度为零，即不模拟初始地应力场；(3) 与方案1相同

2.3.2 模型包含下伏岩层的厚度对模拟精度的影响

对岩土工程问题的数值模拟，建模时依照经验为研究对象设置一定厚度的基座，而不研究模拟结果对基座厚度取值的敏感性，几乎成为惯例。这对边坡稳定性分析等研究的模拟精度可能影响不大。然而对于储层注气工程地表变形这类远场岩石力学效应的模拟精度则有不容忽视的影响。这是由于，储层注气过程中，下伏岩层受储层膨胀扰动引发的弹性压缩变形，相当于对地表竖向位移起到减小的作用，对地表变形的模拟结果有较大影响。因此地质模型中必须包含储层至储层下部相当厚度的岩层，直至岩层弹性变形可以忽略的深度。

为探讨合理的建模范围，研究了模型中包含下伏岩层的厚度对地表变形分析结果的影响。建立3组对照模型：对于储层部分及储层上部已有详细地质信息的岩层，物理力学参数均按实测值，将已有详细地质信息的下伏岩层(厚约200 m)物理力学参数改为表 2中所列信息，并将模型中储层下伏岩层的厚度分别依次扩展200，400，600，…，5 200 m，形成3组模型。不同组模型的储层下伏岩层物理力学参数不同，同一组模型又由下伏岩层厚度依次增加的若干个模型组成。对这3组模型按照方案3计算由注气引起的地表变形。

3 结果与分析 3.1 地表变形计算结果分析

方案1为常规模拟流程，先使模型获得初始地应力状态，再模拟注气引起的变形；方案2，计算过程中由于将凹凸起伏较大的地表岩层的密度设为零，在重力作用下计算平衡后，模型中地表岩层初始应力为零，其余岩层初始应力状态与方案1相近；方案3，模型中所有岩层初始应力全部为零。前已述及，根据弹性力学叠加原理，3个数值模拟方案理论上应取得相同的变形计算结果。

模拟该工程持续注气1 096 d，共注入约33×10$^4$ t CO$_2$后，所产生的竖向位移如图 8所示。W代表注气井位置，A─A'剖面代表图 6中的x=654 500 m剖面。图 8a、图 8b、图 8c为使用本文第2.2节中介绍的仅包含3 400 m有详细地质信息的岩层的地质模型模拟得到的竖向位移结果；图 8d为使用本文第2.3.2节所述第一组模型中下伏岩层扩展2 400 m的模型计算得到的竖向位移结果。

从图 8a、图 8b中可以看出，经过3年的储层持续注气过程，场地地表产生了毫米尺度的竖向位移，竖向位移在地表以注气井所在位置为中心，向四周逐渐减小，位移等值线呈不规则的同心圆状。

图 8dA─A'剖面位移等值线图清楚显示了注气引起的竖向位移在储层和上下盖层中的分布情况。注气渗流过程引起储层孔隙压力增加，有效应力相应减小，储层发生膨胀，使得储层上半部分岩体和下半部分岩体的位移情况展现出相反的趋势，表现为图中上半部分位移为正，下半部分位移为负。

上下盖层岩体通过自身压缩变形不断减弱储层膨胀引起的扰动位移。对上覆盖层而言，这一过程造成岩层隆升，且距离储层越远隆升幅度越小，表现为图中上覆岩层竖向位移为正，且距离储层越远竖向位移越小；对下伏岩层而言，这一过程造成岩层沉陷，距离储层越远沉陷幅度越小，表现为图中下伏岩层竖向位移为负值，且距离储层越远处，竖向位移的绝对值越小。由于场地地表自由，而储层下伏岩层受到下部岩体的约束，储层上覆岩层的隆升幅度要普遍大于储层下伏岩层的沉陷幅度，表现为图中上覆盖层的位移绝对值普遍大于下伏岩层。

图 8c位移分布规律与图 8d大致相同，但是由于模型中包含的下伏岩层厚度不够，储层下伏岩层的压缩变形在模拟中计入不足，导致上覆盖层竖向位移结果偏大，这也证明了开展模型中所包含下伏岩层的厚度对地表变形模拟精度的影响这一研究的必要性。

图 9为采用方案2计算得到的竖向位移云图。方案2与方案3模拟所得结果基本相同，不再累述。

3.2 常规模拟流程的局限性及解决方案

图 10为采用岩土工程领域的常规模拟流程(先模拟初始地应力场，再计算具体工程问题)，即方案1，计算得到的模型竖向位移云图。按照弹性力学理论，仅出现在模型局部的较大位移，必然只能由该局部岩体内的较大应力改变引起，而实际工程中，距离储层几千米远的地表岩层不存在因为储层注气而产生局部较大应力改变的可能，因此图中仅出现在贴近地表数层岩层的局部较大位移(图 10中圈出部分)，属于计算误差。

此计算误差仅出现在贴近地表的岩层中，且对比图 10与图 6，发现计算误差出现区域与地表起伏较大区域相吻合。由图 8、图 9可知，采用方案3、方案2模拟该工程问题，计算结果并未出现此计算误差。可见此计算误差的出现需要两个条件：凹凸起伏较大的地表岩层、计算流程中模拟了贴近地表附近岩层的初始地应力。据此，有以下分析。常规计算流程，模型在自重下逐渐沉降平衡模拟初始地应力场的计算过程中，模型地表附近岩层节点力较小，模型下部节点力较大。所有网格点力的平均值为典型内力，FLAC$^{\rm{3D}}$认为当最大不平衡力与典型内力的比值小于某定值时，计算达到平衡，此时，受力大的模型下部岩层计算精度较高，而受力小的地表处单元应力计算精度难以保证。若再加上地表局部高低起伏较大，网格刻画不能十分平滑，该处的单元应力则可能存在较大的计算误差。另一方面，储层注气对数千米外地表岩层的变形扰动十分微小，这一微小变形扰动与相对较大且精度不高的初始地应力场相互作用，产生了较大的误差，以致出现了图 10中的地表岩层凹凸起伏较大处的局部异常位移。

详细的误差产生原因还有待进一步的研究。但可以肯定的是，对储层注气引起地表变形等远场微小变形的数值模拟，计算初始地应力场近乎相当于画蛇添足，因而不适宜采用常规模拟流程。本文方案2和方案3的对地表变形的良好模拟效果证明，当采用完全弹性本构模型时，计算远场微小变形问题可以不计算初始地应力，而直接模拟具体工程问题；当储层部分采用弹塑性本构模型，不得不计算场地初始地应力时，可以采用忽略地表附近岩层自重的方法改善地表变形模拟效果。

3.3 建模范围对地表变形模拟精度的影响

采用本文第2.3.2节所述的3组模型，利用方案3的模拟流程，得到图 11所示的3组结果，其中横坐标代表每个模型下伏岩层所扩展的厚度，纵坐标代表该模型W点竖向位移计算结果。

从图 11可以看出，随着模型中下伏岩层扩展厚度的增加，3组结果中，地表竖向位移均逐渐减小，最终趋于稳定(图中相邻两点的纵坐标值差别在0.1%以内，起始点在图中以圆点标示)。这说明，只有当模型中所包含的下伏岩层厚度足够大时才能得到精确的地表位移模拟结果；若深度方向建模范围太小，即模型中所包含的储层下伏岩层太薄，则计算得到的由注气引起的地表位移将偏大。以第一组为例，下伏岩层扩展厚度为0和扩展厚度为3 000 m的两个模型(两个模型实际包含的下伏岩层厚度分别为200 m和3 200 m)W点的竖向位移模拟结果相差40%左右，也即，如果使用仅包含200 m厚的下伏岩层的地质模型计算注气引起的地表竖向位移，其误差将达到40%左右。

3组结果中，第一组和第二组达到稳定时的模型下伏岩层扩展值相近，而第三组则明显偏大。这是由于，当储层下伏岩层刚度较小时，其受储层膨胀扰动产生的弹性变形较大，因此对地表位移等远场岩土力学效应的分析结果影响较大。同时也说明，储层下伏岩层的弹性模量对深度方向建模范围取值的影响较大，而泊松比对其影响较小，在选取模型深度方向建模范围时，应着重参考储层下伏岩层的弹性模量大小。当储层下伏岩层的弹性模量偏小时，为获得精确的地表位移模拟结果，应适当加大模型深度方向的建模范围。

4 结论

（1）建立详细精确的复杂地质模型对改善储层远场岩土力学效应分析效果非常重要。自主研发了用户可控的大规模网格化离散方法，将Eclipse储层模型自动转换并延伸成计入岩层详细拓扑形状和物理力学参数并能兼顾计算效率的FLAC$^{\rm{3D}}$全场地质模型，并将Eclipse储层流动分析结果转换使其为FLAC$^{\rm{3D}}$岩石力学分析所用，从而实现复杂地质结构储层流动与近远场岩石力学效应的单向耦合。

（2）利用自主研发建模方法及已有的储层流动分析结果，模拟了某大型咸水层二氧化碳埋存工程持续注气一定时间后的场地变形情况，并讨论了由注气引起的场地变形规律。

（3）岩土工程常规模拟流程(先模拟初始地应力，再计算具体工程问题)由于在远场的计算结果精度不足，不适宜用于分析储层注气引起的地表变形等远场微小变形问题。当采用完全弹性本构模型时，计算远场微小变形问题可以不计算初始地应力，而直接模拟具体工程问题；当储层部分采用弹塑性本构模型，不得不计算场地初始地应力时，可以采用忽略地表附近数层岩层自重的方法改善地表变形模拟效果。

（4）储层下伏岩层受储层膨胀扰动引起的压缩变形对储层远场岩土力学效应分析结果有一定的影响，盲目选取深度方向建模范围可能对模拟结果造成较大误差。深度方向的建模范围应主要根据储层下伏岩层的弹性模量的大小来取值。

（5）本文对地表变形的模拟，考虑了储层流动与岩石力学的耦合作用，并简化为线弹性问题，模拟结果对储层管理、历史拟合、灾害预测等有一定的参考价值，但仍有很大的进步空间。储层近场岩层在储层工程中经历高压强($\Delta p$)、温度($\Delta T$)、化学($\Delta C$)改变，这些与岩石有效应力改变($\Delta \sigma^{'}$)相互耦合。因而，对于储层工程岩石力学效应的数值模拟，全场模型如能在储层近场计入$\Delta p-\Delta T-\Delta C-\Delta \sigma^{'}$，耦合作用，且考虑近场岩层弹塑性变形，而远场保持弹性分析，不仅能在远场取得更高精度的模拟结果，也能使储层近场的模拟结果符合实际。