一种基于甲烷碳同位素的页岩含气量评估方法与流程

本发明涉及石油工程领域，特别是涉及一种评估页岩含气量时利用甲烷碳同位素的分馏表征解析过程的方法。

背景技术

页岩气是滞留于烃源岩中聚集形成的天然气藏，具有典型的“自生自储”的特点，主要以吸附和游离的方式赋存于页岩中。其中吸附气含量介于20～85％之间(平均为50％)，是页岩气的重要组成部分，在很大程度上决定着页岩中的含气量。虽然富有机质泥页岩地层层系中可能广泛含气，但只有地层中的含气量达到一定水平(如美国的含气量底限为0.5m³/t-1.0m³/t)并形成相对富集时，才具有工业开发价值。如果泥页岩地层中的含气量太低，达不到一定水平，那么在目前的经济技术条件下可能就不具备工业开发的条件，对这部分页岩气开展的资源评价就无实际意义了。页岩气的勘探开发已经引起国内外广泛地关注和重视，北美页岩气勘探开发的巨大成功更是掀起了页岩气的热潮。

与美国等发达国家相比，中国的页岩气勘探开发研究尚处于快速起步阶段，尤其在页岩含气量实验方法和评价技术方面研究还比较薄弱。页岩含气量是页岩气资源评价、有利区优选、勘探开发部署等的关键参数和重要依据，针对页岩含气量分析的实验和理论研究迫在眉睫。

现有对含气量评估方法主要包括直接法和间接法两种。

直接法是对页岩样品进行解析测量，然后通过计算去拟合实验结果得到含气量。直接法主要包括usbm法、acf法和mcf法；

(1)usbm方法借鉴煤层气的评价方法，解析量可以通过实验相对准确获取，研究含气量的关键在于估算损失气量。当采收率小于0.5的时候解吸气量与时间的1/2次方成正比。由于样品在最初的几个小时里面温度与储层温度不同，采用初始的解析数据计算会带来误差，不能反映真实的气体扩散特征；

(2)acf方法可以使用所有的可用解析数据，acf方法考虑从零时刻到解吸结束这个过程是连续的，可以用统一扩散模型描述；

(3)在mcf中，区分宏孔和微孔，建立双孔扩散模型。

但上述方法在估算损失气量时存在不准确的问题。损失气量是指在钻井、提钻直到密封装罐之前逸散的天然气量，通常为游离气及靠近岩心表面的吸附气。损失气量无法通过实验测量，一般通过计算求取。在计算过程中，岩心解吸测定前的暴露时间这一参数的确定尤为重要。由于早期扩散速度快，损失气量是非常大的一部分，需要对这部分进行准确的定量评估，但直接法对损失气量的计算方式非常粗糙，同时其残余气量测量方式也不够准确。

间接法分别测量自由气量和吸附气量，然后求出总含气量。其测自由气量依赖于孔隙度的测试，而吸附气量则利用高压甲烷吸附实验得到。对于游离气部分，通常借助于地质模型来描述气藏的几何形态，通过对气藏厚度、孔隙度、含水饱和度及储层的平面展布进行评估来确定模型所需参数，将这些参数输入到地质模型，从而确定气藏的体积。结合气藏压力、温度条件下的流体性质，评估出气藏中单位岩石游离气的体积。对于吸附气部分，页岩的吸附能力可以通过甲烷等温吸附实验，建立吉布斯模型、势差模型或兰格缪尔模型来计算得到。间接法存在问题：实测吸附解吸时，很难保证与实际地层条件一致，而且测试时将岩芯磨成粉末，与实际情况同样存在差异。

技术实现要素：

本发明的目的是要提供一种涉评估页岩含气量时利用甲烷碳同位素的分馏表征解析过程的方法。

特别地，本发明提供一种基于甲烷碳同位素的页岩含气量评估方法，包括：

步骤100，取一规则的岩屑，给定岩屑的内部压力值、外部压力值、吸附量和孔隙度参数；

步骤200，利用前述参数分别计算出两种同位素甲烷由岩屑内部扩散出来的气体扩散通量、表面覆盖率和吸附量，根据控制体积内的自由气量和吸附气量的变化之和等于从边界流入的气量，建立质量守恒方程；

步骤300，将质量守恒方程转换为对两种同位素甲烷的分压求解，计算出两种同位素甲烷在任一时刻的分压，然后求出任一时刻两种同位素甲烷累积放气量，计算累积同位素值，与岩屑解析实验进行匹配，如果误差大则返回步骤100重新给定参数，如果误差小则进入下一步骤；

步骤400，对两种同位素分压取初始状态，则将两种同位素量相加可得到岩屑的总含气量，同时可根据前述获取不同时刻的数据计算出自由气量、吸附气量、损失气量、解析气量和残余气量。

在本发明的一个实施方式中，所述步骤200中，计算气体扩散通量的公式如下：

式中，ni为气体组分i的扩散通量，i值取1和2并对应c13甲烷和c12甲烷，pi为气体组分i的分压，φf为孔隙度，τ为迂曲度，r为纳米孔道半径，μ为气体粘度，r为普适气体常数，t为温度(k)，dk,i为气体i的knudsen扩散系数，对于纳米管而言，为

mi为气体i的摩尔质量。

在本发明的一个实施方式中，所述步骤200中的表面覆盖率计算公式如下：

式中，θi为组分i在流动通道壁面的表面覆盖率bi通过下式进行计算：

bi＝k/p⁰(4)

式中p⁰是标准压力，k为langmuir系数，可以表示为

q＝ea-ed是吸附热，δs⁰＝rln(aa/ad)是吸附的标准熵。

在本发明的一个实施方式中，所述步骤200中的吸附量计算公式如下：

vi,stp＝vlθi(6)

vi,stp为组分i的吸附量，vl为极限吸附量；

单位体积岩屑吸附气中组分i的摩尔数：

式中，ρshale为岩屑密度，ρgi,stp为组分i在标况下的密度,mi为组分i的摩尔质量。

在本发明的一个实施方式中，建立质量守恒方程的过程如下：

根据控制体积内的自由气变化量和吸附气变化量之和等于从边界流入的气量，得到质量守恒方程

式中，v为控制体积，φr为存储气空间孔隙度，ρg1,stp和ρg2,stp分别两种同位素在标况下的密度，和为时间偏导。为控制体积内的自由气变化量，为控制体积内吸附气变化量，等于从控制体积边界流入的气量之和。

在本发明的一个实施方式中，所述步骤300中将质量守恒方程转换为对两种同位素的分压求解过程如下：

对质量守恒方程中的进行转换，化为用气体分压表示的形式：

最终方程化为用分压p1和p2为未知数的方程:

上式中的矩阵a为：

首先，根据高斯公式将式(10)右边的积分化为从控制体积的边界面流入控制体积的扩散通量之和，并进行离散化，即：

式中，ω指的是控制体积的边界；

采用显式迭代方法，以p1、p2为未知数，建立离散化的代数方程，式(10)变为：

进一步转换后得到任一时间时两种同位素的分压：

式中，δt为时间步长，p1ⁿ⁺¹和p2ⁿ⁺¹分别代表第n+1时间步时两种同位素的分压，p1ⁿ和p2ⁿ分别代表第n时间步时两种同位素的分压，在显式迭代中，矩阵a和扩散通量ni都采用第n时间步的参数进行计算。

在本发明的一个实施方式中，所述步骤300中计算累积同位素值的过程如下：

首先根据t时刻的压力，计算两种同位素t时刻的剩余气量：

式中，g1(t)和g1(t)分别代表岩屑内两种同位素的剩余气量；g1k(t)和g2k(t)分别代表网格k内两种同位素的剩余气量，计算时，右边的参数都取t时刻的值；

累积放气量q1(t)和q2(t)＝t0时刻的剩余气量-t时刻的剩余气量：

式中，δcum(t)为t时刻的累积同位素值。

在本发明的一个实施方式中，所述总含气量的计算过程如下：

首先取两种气体w1、w2初始条件下的分压，以确定未放气时的储气量，再通过下式计算出t＝0时两种气体w1、w2的气体量：

式中，右边的参数都取t＝0时刻的；

总含气量w为：w＝w1+w2。

在本发明的一个实施方式中，初始状态下的所述自由气量w自采用下式得到：

在本发明的一个实施方式中，初始状态下的所述吸附气量w吸通过下式得到：

在本发明的一个实施方式中，初始状态下的所述损失气量w损通过下式得到：

设t＝0～t1时间内的气量损失，那么损失气量等于含气量减去t1时刻的剩余气量，计算公式为：

在本发明的一个实施方式中，初始状态下的所述残余气量为：

设放气进行到t2时刻停止，那么能够收集到的解析气量w析为：

则残余气量w残为：

本发明利用页岩气解析放气过程中甲烷碳同位素值的变化，与放气过程中的气体采收率建立对应关系，实现总气量的计算。本发明利用同位素分馏提供的估计含气量方式更重要的是，一方面可以将损失气、解吸气、残余气都计算出来；另一方面又可以将自由气、吸附气都算出来。

附图说明

图1是本发明一个实施方式的评估方法流程示意图；

图2是岩屑的放气过程模型；

图3是本发明一个实施方式的岩屑放气实验示意图；

图4是岩屑放气实验测得的同位素值随时间变化示意图；

图5是本发明一个实施方式的计算值中气量组成示意图；

图6是本发明一个实施方式的计算值与实验值的匹配示意图：

图7是本发明一个实施方式的计算值中气量组成示意图。

具体实施方式

在下面的说明中，当某个参数后面带i时，表明此时计算的是两种同位素甲烷分别的量，如果未带i则指整体的量。

如图1所示，本发明一个实施例公布一种基于甲烷碳同位素的页岩含气量评估方法，包括如下步骤：

步骤100，取一规则的岩屑，给定岩屑的内部压力值、外部压力值、吸附量和孔隙度参数；

其中的岩屑如图2所示，初始条件：岩屑内含高压气体，压力为储层压力。其中，岩屑是由井底钻井后提取上来的，初始岩层里面含有高压的甲烷，根据页岩储层气体压力可得到岩屑中高压气体的压力值。

边界条件：出口定压边界，在岩屑放气过程中，外界压力为一个大气压，这里以一个大气压作为定压边界。

区分两个孔隙度：岩屑内存在储存空间孔隙度φr和流动通道孔隙度φf。在页岩中，气体主要储存在有机质的孔隙中，主要有吸附气和游离气两种赋存形式。而在气体的扩散过程中，气体流动受狭窄的吼道控制。页岩的储气空间大，而流动通道空间小。具体的岩屑可以直接从获取的岩样图像中从有机质总孔隙度大小和基质中连通的孔隙度大小分别给出储存空间孔隙度φr和流动通道孔隙度φf。在计算的方程中，φr影响气量，φf影响有效扩散系数。两种气体的扩散通量中采用φf，而自由气量采用φr。扩散通量指的是从样品内部扩散出来的气体量，也就是下面公式中ni。自由气量指的是储存在样品孔隙中的自由气体，区别于吸附在岩屑孔隙度表面的气体，φr在质量守恒方程第一项中引入，以控制自由气量的多少。

步骤200，利用前述参数分别计算出两种同位素甲烷由岩屑内部扩散出来的气体扩散通量、表面覆盖率和吸附量，根据控制体积内的自由气量和吸附气量的变化之和等于从边界流入的气量，建立质量守恒方程计算：

在尘气模型中，着重考虑努森扩散和粘性流动，气体扩散通量采用下式得到：

式中，ni为气体组分i的扩散通量，分馏过程中存在两种碳同位素甲烷气体，i值取1和2，对应组分1和组分2，其中气体组分1即c13甲烷，气体组分2即c12甲烷。pi为气体组分i分压，φf为流动通道孔隙度，τ为迂曲度，r为纳米孔道半径，μ为气体粘度，r为普适气体常数，t为温度(k)。dk,i为气体i的knudsen扩散系数，对于纳米管而言，为

mi为气体i摩尔质量。

上述公式得到了两种气体组分的扩散通量ni，为后面的质量守恒方程提供计算基础。

表面覆盖率采用如下方式得到：

在计算时需要考虑竞争吸附，即当管道内有两种气体时，两种气体会在管道壁面发生吸附，由于两种气体的吸附能力不同，会竞争壁面位置，这就是竞争吸附，竞争吸附的结果体现在相同压力下两种气体的吸附量不同。对于页岩，吸附气占了很大的比例，故而引入吸附和解吸是必要的；对于一般的砂岩，可以不考虑吸附解吸的影响。

引入表面覆盖率来表示竞争吸附，表面覆盖率指单位面积表面已吸附分子数与单位面积表面按二维密堆积所覆盖的最大吸附分子数之比。表面覆盖率θ取值在0～1之间，0表示管道壁面没有吸附，1表示管道壁面吸附了最大分子数。存在两种组分气体时，气体i表面覆盖率θi可以用下式计算：

式中，θi为组分i在管道壁面的表面覆盖率用下式进行计算：

bi＝k/p⁰(4)

式中p⁰是标准压力(1bar)，k为langmuir系数，可以表示为

q＝ea-ed是吸附热，δs⁰＝rln(aa/ad)是吸附的标准熵。这里获取的表面覆盖率θi同样是后面质量守恒方程式(8)的计算基础，代表了吸附气在分馏过程中的影响。对于c12甲烷，q＝15.1kj/mol，δs⁰＝-76.0j/mol/k。对于c13甲烷，k^*＝αk，

上述参数的值可由实验中测试得到，也可以从现有文献中获取，该上述几个参数先代入式(5)中计算langmuir系数k，然后通过式(4)计算bi,最后通过式(3)计算表面覆盖率，最终放入质量守恒方程式(8)中。

吸附量的计算过程如下：

在标况下，组分i的吸附量(m³/kg)为：(i取1和2，组分1即c13甲烷，组分2即c12甲烷)

vi,stp＝vlθi(6)

vi,stp为组分i的吸附量，vl为极限吸附量(m³/kg)。

单位体积页岩中吸附气中组分i的摩尔数为：

式中，ρshale为页岩密度，ρgi,stp为组分i在标况下的密度,mi为组分i的摩尔质量。

在控制体积v上，根据质量守恒，可知控制体积内的自由气量和吸附气量的变化之和等于从边界流入的气量，建立方程：

式中，v为控制体积，φr为储气空间的孔隙度，ρg1,stp和ρg2,stp分别为c13甲烷和c12甲烷在标况下的密度。和为时间偏导。为控制体积内的自由气变化量，为控制体积内吸附气变化量，等于从控制体积边界流入的气量之和。

步骤300，将质量守恒方程转换为对两种同位素甲烷的分压求解，计算出两种同位素甲烷在任一时刻的分压，然后求出任一时刻两种同位素甲烷累积放气量，计算累积同位素值，与岩屑解析实验进行匹配，如果误差大则返回步骤100重新给定参数，如果误差小则进入下一步骤；

对式(8)中第二项中的表面覆盖率的时间偏导进行转换，化为用气体分压表示的形式：

最终方程化为用分压p1和p2为未知数的方程，即c13和c12甲烷的分压：

上式中的矩阵a为：

首先，根据高斯公式，将式(10)右边的积分化为从控制体积的边界面流入控制体积的扩散通量之和，并进行离散化，即：

式中，ω指的是控制体积的边界。

采用显式迭代方法，以p1、p2为未知数，建立离散化的代数方程，式(10)变为：

式中，δt为时间步长。p1ⁿ⁺¹和p2ⁿ⁺¹分别代表第n+1时间步c13甲烷和c12甲烷的分压，p1ⁿ和p2ⁿ分别代表第n时间步c13甲烷和c12甲烷的分压，在显式迭代中，矩阵a和扩散通量ni都采用第n时间步的参数进行计算。

上述计算过程步骤如下：

step1:已知第n时间步的每个网格p1、p2值，求组分摩尔分数

x1＝p1/pt＝p1/(p1+p2)，x2＝1-x1(14)

step2:计算扩散通量ni，(p,t,x，μ)等参数均采用界面上第n步的值计算，界面上的值由相邻的两个网格中心加权后的压力计算得到。

step3:计算矩阵a，显式迭代中，参数采用网格内第n步的值。

step4:更新压力

step5：得到第n+1步的压力，判断是否达到给定的时间，没有达到则返回step1。

同位素值的定义为：

式中，(¹³c/¹²c)sample是实测的c13甲烷和c12甲烷的浓度比，(¹³c/¹²c)standard是标准中的浓度比，pdb中常用的标准(¹³c/¹²c)standard＝0.011237。

接下来，通过两种气体的放气量计算累积同位素值。根据质量守恒方程式(8)求出t时刻两种气体的累积放气量，c13甲烷的累积放气量为q1(t),c12甲烷的放气量为q2(t)。

首先根据t时刻的压力，计算两种气体t时刻的剩余气量：

式中，g1(t)和g1(t)分别代表样品内c13甲烷和c12甲烷总的剩余气量；g1k(t)和g2k(t)分别代表网格k内c13甲烷和c12甲烷的剩余气量。计算时，右边的参数都取t时刻的值。

累积放气量q1(t)和q2(t)(摩尔数mol)＝t0时刻的剩余气量-t时刻的剩余气量：

计算放气累积同位素值：(积分值，从采集的0时刻到t时刻的累积值，代表收集腔内气体的同位素值)同样，根据方程求出t时刻的累积放气量，q1(t)和q2(t)，然后计算t时刻的累积同位素值，实验中测试的也是累积同位素值，计算值与实验值进行匹配。

式中，q1(t)和q2(t)分别为c13甲烷和c12甲烷的累积放气量，mol；δcum(t)为t时刻的累积同位素值。这个公式的结果是计算出了t时刻的累积同位素值，进而得到同位素值随时间变化的曲线，并与实验结果进行匹配。如果匹配结果误差大，返回步骤100，调整计算的输入参数，以使得误差减小；当结果匹配很好时，那么这些参数也就确定下来，进入步骤400。通过这些确定的参数，计算出整个放气的过程。

步骤400，对两种同位素分压取初始状态，则将两种同位素量相加可得到岩屑的总含气量，同时可根据前述获取不同时刻的数据计算出自由气量、吸附气量、损失气量、解析气量和残余气量。

含气量计算：气体1和2的分压取为初始条件，此时得到未放气时的储气量。先分别算出两种气体的含气量，即t＝0时的样品中的气体量，如下：

式中，右边的参数都取t＝0时刻的值。

总含气量为：w＝w1+w2。

初始时刻的自由气量为：

吸附气量为：

损失气量也称散失气量或者逸散气量，是在钻井、提钻过程中及装入解析灌之前逸散的天然气数量，通常为游离气及靠近岩心表面的吸附气，需根据散失时间的长短及实测解析气量的变化速率进行理论计算。损失气量计算：设t＝0～t1时间内的气量损失，那么损失气量等于含气量减去t1时刻的剩余气量，计算公式为：

解析气量指将岩心装入解析仪之后一定时间范围内所解析出的天然气数量，主体反映为岩心中的吸附气数量及部分游离气量。设放气进行到t2时刻停止，那么能够收集到的解析气量为：

残余气量是解析过程结束后残余在岩心内部或经过更长时间才能获得的天然气数量，一般不代表岩心含气量的主体。残余气量也是含气量的重要组成部分，残余气量的多少视样品及实验的情况不同而不同，有时候残余气量可以占到相对较大的比例，因此不能忽视残余气量的测量。残余气量为：

页岩含气量由三部分组成，即损失气量、解析气量和残余气量，三者相加之和为页岩含气量。

实施例一：

步骤一：进行岩屑现场解析实验，测试时长72小时，在1小时、6小时、12小时、24小时、48小时和72小时这几个时间点测量解析气的同位素值。获取同位素值随时间变化曲线。

实验过程主要利用页岩岩屑放气过程中同位素的动态变化来反映页岩气的勘探地质特征。如图3所示，将从井筒中反排出来的岩屑用图示的罐子收集起来，以一个连续的时间序列测量罐子顶部空气中ch4的同位素值。因为页岩孔吼细小，在放气的过程中会发生同位素的分馏(c12甲烷和c13甲烷跑的速度不一样)。比如，对于丁页五井其中一块样品的同位素测试结果如图4所示，反映了同位素值随时间的变化。实验过程中，岩屑从储层钻取到反排的过程中存在损失气，最终测量的值是损失之后岩屑内的累积同位素值。

步骤二：给定一系列参数，计算模型的累积同位素值随时间变化结果。给定孔隙半径为2nm，储气孔隙度为5％，岩屑半径为5mm，流动孔隙度为0.125％，vl为0.01m³/kg，损失时间为20分钟。按照前述公式进行计算，得到累积同位素值随时间的变化。

步骤三：计算结果与实验结果进行对比，然后根据与实验结果的参数估计相关参数，在估计参数时考虑参数的取值范围及合理性。

将实验结果和计算结果进行对比，将累积同位素值放入同张图6内，发现计算结果与实验结果有较好的一致性，故而验证了参数的正确性。

步骤四：利用得到的相关参数计算放气过程中采收率的变化，计算含气量。通过前述公式计算损失气量，解析气量，残余气量以及自由气量和吸附气量。

图5、6、7为本发明方法得到的效果示意，其中，图5表明总含气量为14.4m³/t，损失气占比56％，解吸气41％，残余气3％。此外，吸附气为7.9m³/t，自由气为6.5m³/t，吸附相占比为55％。

图6表明累积同位素值随时间变化。其中深度：3808.6m，同位素真值-28.25；计算参数：孔隙半径为2nm，储气孔隙度为5％，岩屑半径为5mm，流动孔隙度为0.125％，vl为0.01m³/kg，损失时间为20分钟，含气量为14.4m³/t。

图7表明总含气量为14.4m³/t，损失气占比56％，解吸气41％，残余气3％。此外，吸附气为7.9m³/t，自由气为6.5m³/t，吸附相占比为55％。

本发明的同位素分馏提供了一种估计含气量的方式。更重要的是，一方面可以将损失气、解吸气、残余气都计算出来；另一方面又可以将自由气、吸附气都算出来。

至此，本领域技术人员应认识到，虽然本文已详尽示出和描述了本发明的多个示例性实施例，但是，在不脱离本发明精神和范围的情况下，仍可根据本发明公开的内容直接确定或推导出符合本发明原理的许多其他变型或修改。因此，本发明的范围应被理解和认定为覆盖了所有这些其他变型或修改。