中国石化涪陵页岩气田是中国首个实现商业化开发的页岩气田，也是最大的页岩气田^[1-3]，截至2020 年探明储量为 7 926×10⁸ m^3[4]。涪陵页岩气田气井生产主要采用定产量生产方式，气井生产可以分为稳产阶段和产量递减阶段。现有资料显示，涪陵页岩气田老井第 1 年递减率约为 60%，气井平均稳产期为 2.9 a ^[3]。目前，涪陵老井基本已进入产量递减阶段。涪陵页岩气田初始开发方案采用500～600 m井距，水平井穿行气藏下部层系，矿场生产数据显示，上部气层和下部气层现有井的井间储量动用不充分^[3]。页岩气水平井立体开发模式就是针对储层多个层系甜点分别部署井网，有效挖掘未动用储量^[5]。因此，为了实现涪陵页岩气田稳产接替， 2020年以来现场开展了在上部、中部及下部气层现有井的井间部署调整井的立体开发试验^[6-7]，试验结果显示，立体开发模式可以提高区块整体的动用率和采收率。

为提升页岩气水平井立体开发效果，可采用数值模拟技术，定量确定储量动用和剩余储量分布状况，并以此制定合理的加密井方案。利用数值模拟评估立体开发效果的关键是精确刻画人工裂缝及天然裂缝分布，模拟两者并存条件下的渗流过程，可用于页岩气数值模拟的模型主要包括等效连续介质模型、离散裂缝网络模型和混合模型^[8]。常用的商业模拟软件（如Eclipse，tNavigator等）主要采用等效连续介质模型。该模型简单易用，但当裂缝非均质性较强时，模拟会产生较大误差。因此，并不适用于刻画多甜点层水平井立体开发的储量动用状况。离散裂缝网络模型可以细分为 2 类：①基于非结构化网格的离散裂缝模型（DFM）。它以裂缝面为约束，在裂缝面附近剖分非结构化网格，计算相应传导率。该方法可以较为清晰地刻画每一条裂缝的展布规律、开度等信息，因此近几年在页岩气领域逐渐得到应用^[9-12]。但该方法基于非结构化网格，剖分困难、计算量大，难以在实际矿场模型中得到大规模应用。②嵌入式离散裂缝模型（ED⁃ FM）。它以结构化网格为背景网格，将裂缝嵌入背景网格中，以非相邻连接将二者相连，可以在保证裂缝刻画精度的条件下大幅提高模拟效率^[13-16]。近年来，嵌入式离散裂缝模型也逐渐在页岩气开发中得到应用^[17-19]。鉴于页岩气水平井立体开发优化要求，模拟结果必须能够准确地反映不同工艺压裂效果、井间干扰等因素对井组产能的影响，综合考虑，采用嵌入式离散裂缝模型是较优的选择。为此，笔者提出了一套基于嵌入式离散裂缝模拟的页岩气水平井立体开发优化设计方法，并以涪陵江东区块为例，验证了该方法在实际矿场中的应用效果。模拟结果揭示了现有储层的储量动用情况和剩余气分布，模拟了立体开发的效果，提出调整方案部署建议。

1 研究方法

1.1 数值模拟渗流控制方程

以嵌入式离散裂缝模型模拟页岩气渗流过程，储层中的渗流过程细分为裂缝中的渗流和基质中的渗流两部分，分别由不同的控制方程进行描述^[20-21]。

在基质中，气相和水相的控制方程分别为：

在裂缝中，气相和水相控制方程分别为：

与裂缝型气藏不同，页岩气藏必须考虑基质吸附气解吸的流量。在控制方程中，气相解吸扩散流量可以表示为：

通过实验数据分析，涪陵页岩气平衡吸附浓度满足Langmuir等温吸附公式^[22]，即：

基质与裂缝之间的水相流量的计算公式为：

基质与裂缝之间的气相流量的计算公式为：

与常规气藏基质与裂缝均能向井筒供气情况不同的是，由于页岩基质渗透率极低，因此在模拟过程中，气体只能通过裂缝进入井筒，井筒流量的计算公式为：

通过 C++语言，已将上述理论在中国石化石油勘探开发研究院自主研发的数值模拟软件 COM⁃ PASS 中加以实现^[21]，数值模拟工作均是基于 COM⁃ PASS来完成的。

1.2 优化设计工作流程

在井台地质及动态资料、分析化验、测试及试采资料分析的基础上，首先建立压前三维地质模型，并根据动态分析初步确定缝网分布及属性，基于嵌入式离散裂缝模型建立压后三维数值模型；然后，通过历史拟合确定各井的裂缝分布校正模型参数；最后，在获得校正的数值模拟模型上开展储量动用状况及剩余气分布、井组立体开发影响因素分析，对比评估多甜点层水平井立体开发效果，评价开发技术政策，从而实现立体开发方案优化（图1）。

2 压前三维地质模型的构建

2.1 构造建模

涪陵页岩气田江东区块构造上处于焦石坝背斜带、乌江背斜带和江东向斜带的结合部位，发育2 组不同走向的断层：西南部乌江断裂带附近断层呈近南北向，其他断裂多为北东走向，东北部浅，西南部深，埋深为 2 500～3 500 m，工区总面积为 93 km²。A 平台位于江东区块的东北部，相对构造平缓，小断层不发育，工区总面积为 16.6 km²，总井数为 8口，被近南北向断层分割为 2个区域，断层以西井号分别为 A-1HF，A-2HF，A-3HF，A-4HF，断层以东井号分别为 A-5HF，A-6HF，A-7HF，A-8HF。将区块划分为 25 m×25 m的网格，纵向上分为 9层，共计网格为128×365×9=42×10⁴。

2.2 属性建模

江东区块上部气层页岩品质与焦石坝区块大体相当。基于测井解释结果建立属性模型。①— ⑨小层平均厚度为5.3～12 m，其中①—⑤小层为下部气层，⑥—⑨小层为上部气层。①—⑤小层 TOC 值平均大于3%，孔隙度为5.3%～5.4%，含气饱和度为63%。⑥—⑨小层TOC值平均大于1.6%，孔隙度为 4.6%～4.8%，含气饱和度为 57.7%。下部气层储层压力系数为 1.49～1.74，平均为 1.54，属于高压系统；上部气层储层压力系数为 1.28～1.34，平均为 1.30。根据区块内 8 口水平井的静压测试数据确定，断层分为 2 个压力系统，断层东北地层压力为 51 MPa，断层西南地层压力为48 MPa。根据吸附解吸实验结果，下部气层Langmuir体积平均为1.57 m³ / t，上部气层Langmuir体积平均为0.65 m³ /t。

根据地质模型统计显示，该区块储层物性较好，下部气层储量为100.26×10⁸ m³，丰度为 6.04×10⁸ m³ /km²；上部气层储量为 73.87×10⁸ m³，丰度为 4.45× 10⁸ m³ / km²；总储量为 174.13 × 10⁸ m³，总丰度为 10.49×10⁸ m³ /km²。

3 压后三维数值模型构建

压后三维数值模型构建是页岩气数值模拟的重点和难点。综合运用微地震监测、压裂施工数据及气藏工程解释建立了8口井压后非均匀缝网的初始分布，为下一步数值模拟和历史拟合提供基础。根据微地震监测信息显示，A-1HF井微地震解释平均缝高为39 m，裂缝主要在①—⑤小层内延伸。A-1HF 井微地震监测 46 段，波及体裂缝半长平均为 190 m。

在微地震监测基础上，利用动态储量分析、特征线分析及RTA分析等气藏工程分析方法，确定各井的解释裂缝半长（表1）。其中，改造地层系数是指裂缝总流动面积和地层有效渗透率的负二分之一次幂的乘积，该系数可以在一定程度上反映压裂改造效果和产能大小^[15]。

以解释裂缝半长为压裂改造区的半长，改造区内部随机生成次级裂缝，以微地震响应裂缝半长为各段主裂缝半长，生成该井的压裂缝网。从A-1HF 井压后缝网模型（图2）中可以看出，压裂改造区 （SRV）区域宽度为 123 m，主裂缝半长约为 190 m。该压后缝网可以初步反映压裂和气藏工程的信息，具体参数将在数值模拟中进一步加以调整确定。

4 立体开发方案数值模拟与优化

4.1 数值模拟模型构建及生产数据历史拟合

在地质模型基础上建立 A 平台的数值模拟模型。主裂缝采用嵌入式离散裂缝模拟，分支缝采用粗化算法作为SRV的等效渗透率。

A 平台生产井采用放压生产，中途无明显调配产以及关井开井等制度变化，产量曲线明显为递减特征。通过对 A 平台 8 口井进行生产历史拟合，得到每口井的裂缝参数，主要包括每段裂缝半长、裂缝导流能力和 SRV 平均渗透率。导流能力为 0.61～1.74 mD·m，SRV 平均渗透率为 0.58×10^-4～1.2×10^-4 mD。

4.2 储量动用分析

如图3所示，由于水平井轨迹主要穿越③小层，导致临近的①—⑤小层储量被动用，动用部位压力由初始的 52 MPa 下降到 2~37 MPa。⑥—⑨小层距离水平井较远，除⑥小层由于局部微裂缝沟通导致储量有少量动用外，其余大部分区域压力仍然保持在42 MPa左右的原始状态，储量基本未动用。

选取 3 个剖面分析平面纵向储量动用状况（图4），A-5HF 和 A-6HF 井之间改造较好、未动用区宽度相对较小（约为 200 m），其他井间未动用区宽度约为 270～320 m；A 平台平面上的动用面积为 6.5 km²，动用储量为39.26×10⁸ m³。

如图5所示，按照目前井网预测30 a后，单井压降范围可增加 50～65 m，井间仍然存在较大的未有效动用区域。根据井组数值模拟评价的储量动用状况，截至 2021 年 1 月，①—⑤小层储量未动用区域宽度为240~380 m，⑥—⑨小层储量基本无动用，具有加密调整潜力。

4.3 上、下部气层立体开发调整及效果分析

根据储量动用模拟结果，进行立体开发调整，其中①—⑤小层在老井间（井距为 500 m）加密 1 口井，水平井井轨迹穿行③小层底部；在⑦小层底部，如图6 所示位置处部署 3 口调整井，井距为 250 m。根据此方案，3个井台共设计①—⑤小层中加密4口井，在⑥—⑨小层中调整9口井。

①—⑤小层加密的 4 口井，水平井井段长度为 2 000 m，穿行③小层的下部；在⑥—⑨小层调整的9 口井，水平井井段长度为 2 000 m，穿行⑦小层的下部。上、下部调整/加密井参考试采井动态分析结果，均采用密切割新工艺压裂，考虑风险及工程技术进步，考虑 3 种改造情况（表2）。其中，情况 2 的参数取值为老井气藏工程解释的平均值，即加密井裂缝半长为20～100 m，平均为60 m，①—⑤小层加密井的 SRV 渗透率为 4×10^-5 mD。鉴于⑥—⑨小层可压性劣于①—⑤小层，导致⑥—⑨小层预期压裂效果较①—⑤小层差，⑥—⑨小层调整井的SRV渗透率设为2×10^-5 mD。在相同的压裂规模下，一般压裂施工会出现形成长主缝和网状缝2种情况。与情况 2相比，情况 1考虑的是主裂缝更短，改造区渗透率更高的情况。相反，情况3则考虑主裂缝更长、改造区渗透率略低的情况。经过模拟预测可知，不进行加密，全区块30 a末累积产气量为8.59×10⁸ m³，控制储量的采收率（30 a末的累积产量/（上部井总长× 井距×上部丰度+下部井总长×井距×下部丰度））为 14.2%；若只部署下部井网（包含老井和加密井），3 种情况的全区块30 a末的累积产气量分别为12.52× 10⁸，13.48×10⁸和 14.08×10⁸ m³，平均单井产气量为 1.04×10⁸，1.12×10⁸ 和 1.17×10⁸ m³，控制储量的采收率分别为 20.7%，22.4% 和 23.4%；若只部署上部井网，上部储层动用面积为 4.5 km²，动用储量为 20.03×10⁸ m³，3 种情况下的全区块 30 a 末的累积产气量分别为 6.21×10⁸，8.10×10⁸ 和 8.82×10⁸ m³，平均单井产气量分别为0.69×10⁸，0.9×10⁸ 和0.98×10⁸ m³，相应的控制储量的采收率分别为 10.3%，13.4% 和 14.6%，总体可见上部井比下部井开发效果差；若采用立体开发调整方案，总动用储量为 59.29×10⁸ m³，较单层加密增加了 51%，3 种情况下的全区块 30 a 末的累积产气量分别为 18.73 × 10⁸，21.58 × 10⁸ 和 22.90×10⁸ m³，上部气层调整井与下部加密井间不存在明显干扰。3种情况的立体开发控制储量的采收率分别为31.0%，35.8%和38.0%。

4.4 三层立体开发调整及效果分析

采用上、下 2 个层系的立体开发井间不存在明显的干扰现象。在上、下 2 套调整/加密井网基础上，在⑤小层部署中部调整井（图7）。同样以 A-1HF和A-2HF井所在区域进行分析，通过3层调整/ 加密方案模拟发现全井台（共 7 口井，下部加密井、中、上部调整井分别为 1，2和 2口，压裂改造参数采用情况 2取值）30 a累积产气量为 5.98×10⁸ m³（平均单井30 a累积产气量为0.85×10⁸ m³），30 a末控制储量的采收率达 48.1%。从表3 可知，A-2HF 井由于井间干扰，累积产气量变化较大；老井 A-1HF 由于空间位置关系，干扰较小，累积产气量变化较小。从图8 结果可知，中部调整井 1 与上部调整井和下部加密井均存在压力连通区域，井间裂缝存在沟通。分析中部调整井层间压力场（图9），从下部② 小层到上部⑧小层均有动用。目前的布井方式下，井间层间动用不均匀，造成局部裂缝干扰严重，局部未动用，影响了整体的开发效果。

注：20～80/40表示裂缝半长范围/平均值。

综合所有模拟结果发现：按照目前的改造工艺技术，加密井及调整井 30 a末累积产气量整体上不如现有老井。而在立体开发中，中部调整井开发效果略低于下部加密井。上部调整井效果不佳，30 a 末累积产气量仅为 0.52×10⁸～0.64×10⁸ m³。这主要是上部储层较下部储层物性差所造成的。因此，为了提高上部调整井产量，有必要进一步提升压裂改造效果。

原始井网下不进行加密控制储量采收率为 14.2%，仅进行下部加密控制储量采收率为 22.4%，上、下部气层立体开发控制储量采收率为 35.8%，上、中、下部气层立体开发控制储量采收率为 48.1%。可见，继续推进上、中、下部气层立体开发部署是提高涪陵老区采收率的重要途径。

5 结论

综合运用微地震数据、气藏工程分析数据和生产历史数据等多类型数据确定主裂缝和压裂改造区分布，采用嵌入式离散裂缝可以准确刻画裂缝的非均匀分布，该数值模拟方法可以提高页岩气数值模拟的准确性，可为页岩气水平井立体开发优化设计提供技术依据。运用该方法，以涪陵江东区块为例，综合分析不同工艺压裂结果、井间干扰等因素对井组产能的影响，评估了立体开发的效果，分析结果显示：立体开发是提高采收率的重要途径，需要进一步加大缝网改造效果，降低单井投资，实现效益开发。下一步可在江东区块开展全面立体开发调整，有效提高气田整体采收率。

符号解释

B_g——气相地层压力系数，无量纲；

B_i ——i相地层压力系数，无量纲；

B_w——水相地层压力系数，无量纲；

C_g——基质当前吸附浓度，kg/m³；

C_g ^∞ ——平衡吸附浓度，kg/m³；

D——海拔高度，m；

g——重力加速度，m/s²；

i——表示气或水；

K——绝对渗透率，mD；

K_rg——气相的相对渗透率，无量纲；

K_ri ——i相的相对渗透率，无量纲；

K_rw——水相的相对渗透率，无量纲；

p_g——气相压力，MPa；

p_g ^m——基质中气相压力，MPa；

p_g ^f ——裂缝中气相压力，MPa；

p_L——Langmuir压力，MPa；

p_i——i相地层压力，MPa；

p_w——水相压力，MPa；

p_w ^f ——裂缝中水相压力，MPa；

p_w ^m——基质中水相压力，MPa；

p^w ——井底压力，MPa；

q_g ^mf ——基质与裂缝之间的气相流量，kg/s；

q_g ^S ——气相解吸扩散流量，kg/s；

q_g ^w ——裂缝与井之间的气相流量，kg/s；

q_i ^w ——裂缝与井之间的i相流量，kg/s；

q_w ^mf ——基质与裂缝之间的水相流量，kg/s；

q_w ^w ——裂缝与井之间的水相流量，kg/s；

S_g——气相饱和度，无量纲；

S_w——水相饱和度，无量纲；

t——时间，s；

T_mf ——基质-裂缝间的传导率，m·mD；

v——扩散系数，s^-1；

V_L——Langmuir体积，m³/t；

V_m——基质体积，m³；

WI_i ——井指数，m·mD；

μ_g——气相黏度，mPa·s；

μ_i——i相黏度，mPa·s；

μ_w——水相黏度，mPa·s；

ρ_g——气相密度，kg/m³；

ρ_i ——i相密度，kg/m³；

ρ_w——水相密度，kg/m³；

ϕ——孔隙度，无量纲。

参考文献