随着大气中CO₂浓度不断攀升，全球平均温度较1850—1900年间已升高近1.5℃，由此导致了极端气候频发、粮食产量减少、海平面上升等问题^[1]。为了应对此类现象，全球各国陆续公布了碳中和目标。CO₂捕集、利用和封存技术（Carbon Capture, Utilization and Storage简称CCUS），被认为是实现减排目标的关键技术^[2-3]。其中，地质封存是实现CO₂永久封存最有效的方式^[4]。地质封存的主要场所是海底深部的咸水层或枯竭的油气藏，其中深部咸水层因其巨大的封存潜力以及不具备其他经济价值而被视为最理想的CO₂地质封存场所^[5-6]，据统计，咸水层CO₂封存潜力占总封存潜力的98%以上^[4]。国外目前正在运行的比较成功的咸水层封存项目包括挪威北海的Sleipner项目、Snøhvit项目，加拿大的Gorgon项目等，其中Sleipner作为世界首个海上咸水层CO₂封存示范项目，以运行时间长、封存量大的特点著称，并部署了各种封存监测手段。自1996年起开始注入，至今已封存超过2400万吨CO₂，是研究CO₂咸水层封存的示范案例^[7-8]。在国内，鄂尔多斯盆地的神华咸水层封存项目正在持续运行，中国海油建设的中国首个海上咸水层CO₂封存项目恩平15-1油田也已正式开始CO₂注入，但目前这些项目仍未部署CO₂监测。

针对CO₂注入后的运移和泄露状况的动态监测问题^[9]，地球物理被认为是一种有效的监测手段，用于确认地质封存的安全性、可靠性^[10]。其中，与重力、测井、微地震等其他方法相比，四维地震在监测CO₂体积和波及范围方面具有独特优势，经过先进方法^[11-12]处理后能够提供随时间变化的地下动态^[13]、高分辨率成像^[14]。Huang等^[15]对德国的咸水层CO₂封存项目通过四维地震的振幅差异监测注入的CO₂波及范围的变化，Chadwick和Noy^[16]对Sleipner四维地震数据定性解释来匹配CO₂羽流迁移的历史，Roach和White^[17]使用了类似的方法来描述加拿大的Aquistore封存点的CO₂羽流迁移，Fawad^[18]为提高监测准确性利用四维地震和电磁测量相结合的方法来描绘CO₂羽流。

为了建立CO₂在海底咸水层封存过程中四维地震属性对CO₂饱和度变化的定量描述关系，提高监测CO₂封存波及范围的准确性。本研究结合Sleipner咸水层CO₂封存项目所采集的四维地震监测资料和测井资料，通过Gassmann流体替换方程和Xu-White模型建立不同CO₂饱和度下的速度模型，应用井控地震属性分析技术研究CO₂咸水层封存过程中CO₂-盐水两相介质变化引起的各向异性地震响应特征，提取并分析对CO₂饱和度变化敏感的地震属性。同时，对Sleipner采集的7次时延地震数据进行解释，提取沿层地震属性，监测CO₂封存波及范围的时空变化。

1.   研究区概况

1.1   Sleipner CO₂封存项目概况

Sleipner是Equinor（原挪威国家石油公司）运营的一个重要天然气田，位于挪威北海。Sleipner CCS项目是其商业化的CCS站点，项目将Sleipner East气田的天然气生产过程中所产生的CO₂压缩并通过约2.3 km的斜井注入高孔隙度Utsira咸水层，以防止CO₂释放到大气中^[19-21]（图1）。自1996年9月15日开始，Sleipner平台通过注入15/9-A-16井开始注入CO₂，注入点位于海平面以下1012 m的咸水层Utsira底部，距离储层顶部约200 m，注入射孔长度为38 m^[8]。为了达到所需的井口压力，CO₂需要先经过4个压缩阶段，在注入井井口CO₂处于气液两相流动的相变阶段。最初几年的年注入率约为0.9 Mt，计划注入200年^[22]。由于Sleipner East气田的天然气流量减少，后期略有下降，截至2010年，累计注入量和注入速度统计见图2。

图  1  Sleipner项目CO₂封存点及Utsira地层注入示意图^[21]

Figure  1.  Schematic representation of the Sleipner project CO₂ storage site and injection in the Utsira Formation^[21]

下载: 全尺寸图片 幻灯片

图  2  CO₂日均注入体积及累计注入量

Figure  2.  CO₂ Average daily injection volume and cumulative injection

下载: 全尺寸图片 幻灯片

迄今为止，Sleipner共进行了10次三维地震调查和4次重力调查，通过这些调查，深入研究了储层中CO₂的运移行为。当CO₂进入储层时，通常处于超临界状态，与初始的咸水层形成了强烈的声速对比，为地震监测提供了有利条件。Sleipner的地震监测分别于1994年（Base line）、1999年、2001年、2002年、2004年、2006年、2008年、2010年、2013年和2016年进行了重复三维拖曳地震勘测，其中1994年的地震监测获得了注入CO₂之前的基线三维地震数据^[23-24]。

1.2   地质概况

Sleipner CO₂封存地点位于Utsira砂层^[25]，这是一个深度约800～1100 m、厚度约200～300 m、面积约26000 km²的咸水层，南北长400多千米、东西长50～100 km^[26]，主要由中新世晚期至上新世早期的砂质单元和部分泥质储层组成^[27]（图3）。

图  3  Utsira储层区域剖面^[27]

Figure  3.  WSW-ENE orientation of the Utsira reservoir regional profile ^[27]

下载: 全尺寸图片 幻灯片

通过伽玛测井曲线以及强地震反射（图4）可以有效地识别Utsira储层中的砂泥岩。声速和密度分别在地层单元顶部和底部急剧下降和增加，构成了明显的阻抗差异，并引起强烈的清晰可识别的地震反射，根据15/9-13的测井曲线，对砂泥岩进行速度分析（图5），结果显示随着泥质含量的提高，声波速度整体呈上升的趋势，泥岩的速度整体要大于砂岩的速度。

图  4  15/9-13井GR测井曲线图（a）及2010年Inline 187地震剖面（b）

数据来自https://CO2datashare.org。

Figure  4.  Well 15/9-13 GR logging profile (a), and the Inline 187 seismic profile in 2010 (b)

Data from https://CO2datashare.org.

下载: 全尺寸图片 幻灯片

图  5  15/9-13井砂泥岩速度分析图

Figure  5.  Velocity analysis of Well 15/9-13 sand mudstone

下载: 全尺寸图片 幻灯片

根据测井曲线和三维地震资料，将目标地层分为9个主要的砂层（图4），Utsira储层之上的Nordland页岩层厚度约200～300 m，构成主要的盖层，从测井曲线可以看出，在8层和9层之间还有一层比较厚、横向分布比较广泛的泥岩，厚度约为5 m（图4a）。储层中的薄泥岩层平均厚度为1～1.5 m，构成了储层砂层内重要的渗透性屏障，并对CO₂在储层内的运移和圈闭产生影响。储层单元的含砂量一般为0.7～1.0，15/9-A-16井获得的不同深度Utsira储层岩芯样品的主要岩石和碎屑成分分析结果见表1。

表  1  15/9-A-16井不同深度岩芯样本主要岩石和碎屑成分百分比

Table  1.  Percentage of major rock and debris compositions in core samples from Wells 15/9-16 at different depths

	碎屑含量/%
	850～860 m	890～900 m	1000～1010 m
石英	50.7	66.7	76.7
长石	7.3	3.7	2.7
方解石	18	17	7.7
页岩	4.3	1	4.7

下载: 导出CSV 

| 显示表格

根据挪威国家石油公司数据，Sleipner地区在1038 m附近的地下温度测量值为37℃。地下温度梯度为33℃/km，因此，Sleipner地区Utsira储层的温度范围为储层顶部的28℃到储层底部的41℃^[27]，具体温度、压力以及CO₂密度关系如图6所示。

图  6  温度、压力、CO₂密度随深度的变化曲线^[27]

Figure  6.  Variation curves of temperature, pressure, and CO₂ density with depth^[27]

下载: 全尺寸图片 幻灯片

根据表1矿物成分以及Sleipner地区的实际地质条件、储层温度和压力状况得到岩石基质矿物成分和储层内流体的相关模量和密度。其中在储层条件下，盐水的体积模量为2.3 GPa^[28]，在储层温度和压力下，盐水的密度假设为1.03 g/cm^3[29]。大部分CO₂以超临界状态封存，考虑初始孔隙压力为8 MPa^[30]，储层温度范围为28℃（Utsira上部地层）至41℃（注入点附近），选择CO₂参考密度为700 kg/m^3[31]。CO₂的体积模量变化范围为0.02～0.075 GPa^[32]，根据现场压力和温度条件，选择使用0.075 GPa的CO₂体积模量^[31]，具体参量统计见表2。

表  2  储层中矿物成分和流体的弹性模量及密度

Table  2.  Elastic modulus and density of mineral components and fluids in reservoirs

	体积模量/GPa	剪切模量/GPa	密度/(g/cm3)
石英	37.00	44.00	2.65
长石	37.50	15.00	2.70
方解石	76.80	32.00	2.71
盐水	2.30	0	1.03
CO2	0.075	0	0.70

下载: 导出CSV 

| 显示表格

Utsira储层主要是高孔隙度砂岩，通过岩芯样品的液体侵入测量和密度测井得出的孔隙度值范围为36.0%～40.1%，平均38.0%^[27]。从15/9-A-16井的中子孔隙度测井曲线得到孔隙度为35%～39%，平均37%。Williams^[33]对渗透率和电缆测井数据的重新评估认为东西方向渗透率约为1.975 μm²，南北方向渗透率高于东西方向，约为7.896 μm²。

从地震剖面可以看出，CO₂通过几个相对高渗透性“烟囱”垂直向上运移，穿过储层及其地层内的泥岩^[21]（图7）。“烟囱1”和“烟囱3”对应于地层中微妙的不连续面，位于或接近储层顶部（图7a）。“烟囱1”在重复测量中特别突出（图7b），并且可以作为一个独特的垂直圆形反射，半径约为30 m，它被解释为CO₂通过羽流向上运移的主要垂直管道^[34]。

图  7  四维地震剖面对比显示3个可能的CO₂羽流烟囱

a：1994年基准地震剖面，b：2010年时移地震剖面。

Figure  7.  Comparison of 4D seismic profiles showing three possible CO₂ plume chimneys

a: The 1994 reference seismic profile, b: 2010 time-shifted seismic profile.

下载: 全尺寸图片 幻灯片

2.   CO₂封存地震响应特征正演模拟

为了研究CO₂封存过程中随CO₂饱和度变化的四维地震响应特征，需通过测井数据建立不同CO₂饱和度的正演模型，具体包括以下3个方面的内容：① 采用流体替换方法预测不同CO₂饱和度条件下的纵波和横波速度；② 基于流体替换前后的速度和密度曲线构建水平层状模型，正演模拟得到自激自收地震记录；③ 通过对不同CO₂饱和度下的地震记录提取相应的地震属性，分析对于CO₂饱和度变化敏感的属性以用于地震资料解释和CO₂波及范围监测。

2.1   不同CO₂饱和度地层的速度模型

在进行CO₂封存的过程中，CO₂注入导致储层内原有流体被置换，同时岩石框架和饱和岩石体积模量等岩石物理参数发生改变^[35]。通过分析CO₂注入期间地震资料的变化，能够监测到不同阶段CO₂分布状况。因此需要建立CO₂饱和度与地震速度之间的关系，即对注入CO₂后的地震速度进行估算。

Sleipner咸水层CO₂封存项目的储层为砂泥岩互层，因此本文采用Xu-White模型^[36]求取纵横波速度。该模型是在Gassmann流体替换方程的基础之上，综合考虑了黏土影响和孔隙纵横比，结合Kuster-Toksoz模型^[37]和差分有效介质（DEM）理论，提出的一种依据孔隙度和泥质含量来估算泥质砂岩纵波与横波速度的方法。其中纵横波速度计算公式如下：

式（1）、（2）中，V_S为横波速度，V_P为纵波速度，K_sat为饱和岩石的体积模量，μ_sat为饱和岩石的剪切模量，ρ为饱和岩石密度。

Xu-White速度模型的公式如下：

式中，ϕ为总孔隙度；ϕ_s为砂岩孔隙度；ϕ_c为泥岩孔隙度；f_s为砂的体积分数；f_c为黏土的体积分数；K_dry、K_ma、K_f分别为干岩石体积模量、组成岩石的固体矿物基质的体积模量以及孔隙流体体积模量。μ_dry、μ_ma、μ_f分别为对应的剪切模量。F(α)、T_ijij(α)为关于孔隙纵横比α的函数。公式（6）和（7）为Kuster-Toksoz方程。

Xu-White模型将岩石视为砂质和泥质的复合体。其假定岩石骨架的泊松比大致为恒定值，强调了孔隙形态对地震波速的显著影响。通过引入孔隙的纵横比以刻画孔隙的几何特性，并将其与孔隙率、岩石基质的剪切及体积模量相结合，从而推导出干岩石体积模量K_dry与剪切模量μ_dry的简化表达式：

式中，参数$ p $与$ q $的表达式分别为：

同时，引入了Voigt-Reuss-Hill（VRH）^[38-39]平均理论计算岩石基质模量，公式如下:

式中，M_VRH为岩石基质的弹性模量；M_ν为Voigt有效弹性模量；M_R为Reuss等效弹性模量。

引入Brie^[40]经验模型计算不同CO₂饱和度下的混合流体模量，公式为：

式中，K_f为混合流体的体积模量，K_w为液相体积模量，K_g为气相体积模量，S_w为液相饱和度，e为经验指数因子。

基于上述模型编程计算得到储层在不同CO₂饱和度状态下的速度曲线，图8分别对比了从0开始注入CO₂，在储层CO₂饱和度为0、20%、40%、60%、80%时的流体替换反射纵横波速度曲线。可以看出，随着储层CO₂饱和度的提高，纵横波速度均持续下降，这是由于随CO₂饱和度的增加，Xu-White模型计算得到的岩石体积模量K_sat和剪切模量μ_sat均减小，因此流体替换得到的速度减小。

图  8  不同CO₂饱和度的纵波（a）及横波（b）速度曲线对比

Figure  8.  Comparison of primary wave velocity profiles (a) and shear wave velocity profiles (b) for different CO₂ saturations

下载: 全尺寸图片 幻灯片

图9对比了不同饱和度下的均方根速度，可以看出，随着CO₂饱和度的升高，纵横波均方根均下降，注入CO₂后到CO₂饱和度为60%之前，横波速度下降幅度略小于纵波速度下降幅度，并且幅度随饱和度增加而减小，当CO₂饱和度大于60%之后，横波速度减小的幅度将超过纵波速度。

图  9  随CO₂饱和度升高纵横波均方根速度变化趋势

Figure  9.  Trend of root-mean-square velocity of longitudinal and transverse waves with increasing CO₂ saturation

下载: 全尺寸图片 幻灯片

2.2   地震正演模拟与敏感属性分析

基于不同CO₂饱和度替换得到的速度和密度曲线建立水平层状模型。对其进行自激自收的地震正演模拟得到储层CO₂饱和度分别为0、20%、40%、60%、80%的地震记录，如图10所示。对比发现，振幅变化最为明显，随CO₂饱和度的增加，盖层振幅保持不变，储层振幅提高。在储盖层交界处，相位发生变化。

图  10  不同CO₂饱和度的地震记录

Figure  10.  Seismic profiles with different CO₂ saturations

下载: 全尺寸图片 幻灯片

在不同饱和度的地震记录中沿剖面提取振幅、频率、相位、吸收衰减地震属性。

（1） 均方根振幅属性

其中均方根（RMS）振幅属性是计算振幅平方和的平方根，并将其除以特定时间窗内的样本数量。使用RMS振幅属性是为了突出显示异常的高振幅，包括负值和正值。RMS振幅属性是一个很好的声学对比指标，可以指示岩性和油气的变化^[41-42]，计算公式如下：

式中，N为采样点总数，a_i为各采样点振幅值。

图11a对比了不同CO₂饱和度的均方根振幅属性，随着饱和度的增加，盖层均方根振幅值不变，储层均方根振幅值变大。通过将不同饱和度的均方根振幅属性做差（图11b）可以看出，开始注入时，振幅变化最明显，随着CO₂饱和度增加，振幅增大的幅度减弱，饱和度越大振幅变化越不明显，但仍能看出变化趋势。

图  11  不同CO₂饱和度均方根振幅属性（a）及均方根振幅属性差（b）

Figure  11.  RMS amplitude attributes (a) and difference in RMS amplitude attributes (b) for different CO₂ saturations

下载: 全尺寸图片 幻灯片

（2） 瞬时频率属性

瞬时频率作为相位随时间变化率的量度，能够反映地层的岩性变化。受地层影响，其变化幅度较大，如当地震波穿过某些特殊的储层时，高频成分将会显著衰减，其中含油气储层经常引起高频衰减^[41]。因此瞬时频率属性可反映地震波衰减特征，以此估计地层的周期性，其计算公式如下：

图12a为不同CO₂饱和度的瞬时频率属性，当注入CO₂之后，储盖层的分界面处以及上下两侧频率发生变化，将不同饱和度的瞬时频率属性做差（图12b）可以发现，在储盖层交界处瞬时频率升高，在上下两侧瞬时频率降低，但随着CO₂饱和度的增大，瞬时频率属性的变化敏感性减弱。饱和度60%以上，频率变化几乎为0。

图  12  不同CO₂饱和度瞬时频率属性（a）及瞬时频率属性差（b）

Figure  12.  Instantaneous frequency attributes (a) and difference in instantaneous frequency attributes (a) for different CO₂ saturations

下载: 全尺寸图片 幻灯片

（3） 瞬时相位属性

瞬时相位属性一般是用来反映同相轴的连续性，不考虑振幅强度，地震波穿越不同的地层会引起地震波的相位变化，因此通常是用来识别断层、尖灭点、河道等^[43]，计算公式如下：

图13a对比了不同饱和度的相位属性，CO₂注入前后几乎看不出变化。将不同饱和度的瞬时相位属性做差可以发现（图13b），当开始注入CO₂时，在储层和盖层交界处出现相位变化，当CO₂饱和度增大到40%以上时，瞬时相位属性基本未发生改变。

图  13  不同CO₂饱和度瞬时相位属性（a）及 瞬时相位属性差（b）

Figure  13.  Instantaneous phase attributes (a) and the difference in instantaneous phase attributes (b) for different CO₂ saturations

下载: 全尺寸图片 幻灯片

（4） 瞬时Q值属性

瞬时Q值属性也叫做瞬时品质因子，是一种吸收衰减属性，可通过低频振幅与高频振幅的比值得到。其主要用于描述地震波在介质中传播时的衰减特性。瞬时Q值越大，说明低频能量与高频能量的比值越大，高频的衰减越剧烈。流体饱和度的变化会影响介质的衰减特性，因此可以用来评估地层中的流体饱和度。

图14a为不同饱和度CO₂的瞬时Q值属性图。当CO₂注入后，在储盖层交界处周围发生了较明显的衰减变化。将不同饱和度的瞬时Q值属性做差可以发现（图14b），注入CO₂后，在储层周围出现正负相间的衰减变化，随CO₂饱和度增长，衰减变化值变小，但仍具有敏感性。

图  14  不同CO₂饱和度瞬时Q值属性（a）及瞬时相位属性差（b）

Figure  14.  Instantaneous Q attributes (a) and the difference in instantaneous phase attributes (b) for different CO₂ saturations

下载: 全尺寸图片 幻灯片

3.   CO₂咸水层封存波及范围地震监测

为了研究CO₂羽流在平面以及垂向上的发育情况，沿着泥岩层反射分别解释了上（层9）、中（层5）、下(层1)3个层位(图15)，网格间距为16×16，解释结果如图16所示。根据正演结果选择均方根振幅属性来描述CO₂波及范围。沿解释层位，创建时窗并提取均方根振幅属性。图17显示，以中间层第5层1999年、2006年、2010年3次时延地震资料为例，图中绿色和黄色为振幅高值，蓝和紫色为振幅低值，通过同一层位几次不同时间的属性图对比可以明显看出整体CO₂波及范围随着时间增长，此外，随CO₂饱和度增加，均方根振幅值升高。周围储层的均方根振幅值为0.1左右，CO₂饱和度最高值处振幅在2左右，CO₂羽流波及前缘处振幅值大于0.5，对波及范围边界的识别效果较好，同时与地震正演结果相符。通过均方根振幅属性开展不同层位、不同时间的CO₂波及范围对比，分上（第1）、中（第5）、下（第9层）3个层位提取2010—1999年6次四维地震数据的均方根振幅属性（图18），对比不同时间，上、中、下三层平面上CO₂羽流的纵向扩散速度均大于横向扩散速度，整体均呈SSW-NNE方向展布。

图  15  解释层位Xline 1164

Figure  15.  Interpreting the Layer Xline 1164

下载: 全尺寸图片 幻灯片

图  16  构造解释结果

Figure  16.  Results of tectonic interpretation

下载: 全尺寸图片 幻灯片

图  17  Sleipner咸水层封存第5层四维地震均方根振幅属性

Figure  17.  4-D seismic RMS amplitude attributes for Layer 5 of the Sleipner saline aquifer storage

下载: 全尺寸图片 幻灯片

图  18  挪威Sleipner咸水层封存项目CO₂平面波及范围四维地震属性预测结果

Figure  18.  Predicted 4D seismic attributes of CO₂ planar spread range for the Sleipner saline aquifer storage project in Norway

下载: 全尺寸图片 幻灯片

图18显示了位于注入点上方的第1层平面上的CO₂波及范围，通过其随时间的变化可以看出，开始注入时CO₂呈一个高饱和度的团状羽流，随着时间的推移，CO₂持续注入，CO₂在平面上不断扩散，并在2004年达到最大横向饱和度，这时储层中注入的CO₂量为6.9 Mt，随后沿纵向继续扩散，整体纵向扩散大于横向，随CO₂羽流扩散，中心的饱和度不断减小，表明CO₂正在不断向上层砂岩运移。

图18展示了中部第5层的CO₂波及范围变化情况，在2008年时第5层达到了最大横向扩散，此时储层中注入的CO₂量为10.2 Mt，在2010年注入12 Mt后整体饱和度达到最大，中心振幅小于两端，表明CO₂在波及范围两端的构造高点处形成高饱和度积聚，中心通过高渗透率“烟囱”持续向上运移。对比发现，1999年时第1层和第5层的CO₂波及范围大小相近，随后CO₂波及范围在第5层的增长速度大于第1层的增长速度，最终的CO₂波及范围第1层也大于第5层。

图18同时显示了储层最顶部砂岩层第9层CO₂波及范围变化情况，可以看出1999年已经有CO₂羽流到达最顶层，这进一步证明了存在CO₂高渗透率“烟囱”。2010年CO₂波及范围达到最大，除主体羽流发育以外，CO₂沿着右侧一条通道向北快速运移，到达一个构造高点后积聚，继续向不同的方向扩散，整体上东侧振幅高于西测，通道处振幅值最高，CO₂到达第9层后，先是平面扩散，随后沿通道运移，并在通道中聚集使饱和度提高。第9层整体上CO₂波及范围小于第5层而大于第1层。

4.   结论与建议

本文基于Sleipner咸水层CO₂封存项目的测井资料、四维地震资料，开展了针对CO₂海底咸水层封存的流体替换和地震正演模拟，以及CO₂波及范围地震监测，得到如下认识：

（1）均方根振幅属性对CO₂饱和度变化最为敏感。随CO₂饱和度的增加，饱和岩石的体积模量、体积密度、纵波速度和横波速度均有所下降。相反，地震记录中总体振幅升高，且随着CO₂注入量的增加，振幅变化幅度减小。

（2）利用均方根振幅属性对Sleipner海底咸水层CO₂波及范围进行描述，在注入期间，CO₂在层内主要沿SSW-NNE方向运移，并在构造高部位聚集。垂向上，CO₂从注入点向上层运移，下层达到最大波及范围的时间早于上层。结合构造解释结果和储层物性发现，CO₂在储层内的波及范围主要受各项异性渗透率和构造高低控制。由于井15/9-A-16缺少实测速度值，本文在速度预测中使用的是固定的孔隙纵横比的Xu-White模型，在今后的研究中，可以采用由实测速度进行矫正的岩石物理模型，提高速度预测精度。CO₂在咸水层封存过程中会改变储层的孔隙度等储层参数，在今后岩石物理建模过程中，可以将CO₂对储层孔渗条件的影响考虑在内，建立动态岩石物理模型。