塔河油田奥陶系油藏是典型的碳酸盐岩缝洞型油藏，裂缝、溶洞非常发育，流体储存于缝洞中，基质基本不具备储渗能力^[1,2,3]。由于缝洞型油藏的储集层具有非均质性强、空间结构复杂、缝洞随机分布、溶洞大小不一、裂缝倾角较大等特点，导致碳酸盐岩缝洞型油藏的开发在世界范围内仍处于探索阶段^[4,5,6]。

与常规碎屑岩油藏相比，空间离散性是缝洞型油藏的突出特点^[7,8]。同时，储集体空间分布又具有极强的随机性，这种离散性和随机性导致缝洞型油藏在地质建模、油藏数值模拟方面面临巨大的困难^[9,10]。现有的缝洞型油藏数值模型包括多重介质等效连续模型^[11,12]、耦合流动模型^[13,14]、离散介质试井模型^[15,16]等。多重介质模型是将溶洞作为一种连续介质处理，但实际上溶洞在空间上离散分布，其表征单元体并不存在，将其等效为连续介质是不合适的。耦合流动模型考虑了缝洞介质的离散特征，但溶洞内流体的流动用N-S方程描述，导致计算成本过高，难以用于油藏尺度的模拟。离散介质试井模型针对几个具体组合模式，应用范围有限。为此，笔者基于缝洞型油藏缝洞分布的离散性和随机性，从缝洞型油藏的溶洞、裂缝物性参数统计特征入手，基于地球物理识别结果，建立了规整化的离散缝洞分布模式，引入了等效连通系数和弹性储集系数，建立了缝洞单元离散数值模拟方法。

缝洞型油藏储集体的分布极为复杂，且具有埋藏较深、尺度结构差异较大、溶洞充填特征多样化等特征，很难精确获得缝洞储集体的分布特征。通常，井点周围缝洞分布的特征依据岩心和钻井信息确定，而非井点缝洞分布则在地球物理预测的基础上，结合地质统计信息确定。根据塔河油田9口井的缝洞网状分布状况，建立10 m×10 m×10 m的规则网格系统，将裂缝和溶洞归位到网格中并进行统计，结果见表 1。

分析5口井的成像测井资料，得到裂缝开度、裂缝倾角和裂缝倾向的频率和累计频率分布,如图 1所示。

图 1 裂缝特征参数统计分布特征 Fig.1 Statistical distribution characteristics of characteristic parameters for fractures

图选项

从图 2可以看出：裂缝开度呈现对数正态分布，裂缝开度最频值分布于0.000 1～0.010 0 cm；裂缝倾角变化呈现正态分布，裂缝倾角的最频值分布于65°～90°，中值为80°，即绝大多数裂缝为高角度缝；裂缝倾向在360°范围内均匀分布。

图 2 规整化的离散缝洞分布模式 Fig.2 Regularized discrete fracture-karst distribution pattern

图选项

为了便于研究，基于地球物理解释和地质统计获得的缝洞分布信息，建立了规整化的离散网络缝洞分布模型。其中，贯穿井点的裂缝和溶洞分布依据钻井过程中的钻柱放空、钻井液漏失和岩心等信息确定。井点裂缝、溶洞的位置及相对大小根据地球物理识别预测结果确定，但各条裂缝的开度、倾角、走向根据地质统计信息赋值，各个溶洞的控制储量根据单元地质储量及相对大小分别赋值，如图 2所示。

流体在裂缝中的流动主要受裂缝特征、流体特征和节点两端压差的影响，因此引入裂缝连通系数表征流体在裂缝中的流动能力。单位时间内，流经节点1和节点2之间裂缝的流量为：

式中：J_L12为节点1和节点2间缝内的液相连通流动系数，m³/(MPa·d)；q_L12为经过节点1和节点2之间裂缝的流量，m³/d；pⁿ⁺¹₁和pⁿ⁺¹₂为当前时刻节点1和节点2的压力，MPa。

高速流动条件下，裂缝中流体流动可能表现出非线性渗流，流动规律符合高速非达西表达式。

式中：μ为黏度，mPa·s；v为流速，cm/s；K₀为绝对渗透率，D；ρ为流体密度，kg/m³；p为压力，MPa。

满足广义达西定律的等效渗透率为：

式中：K_e为等效渗透率，D；β为非达西系数。

考虑高速非达西效应的等效渗透率与绝对渗透率的比值F_J为：

则考虑高速非达西效应的流量为：

由于裂缝渗透率与开度有关，该模型采用内蕴裂缝渗透率计算方法。内蕴裂缝渗透率是指流体通过单一裂缝或裂缝网的传导能力，与周围岩石(基质)无关，内蕴裂缝渗透率可表示为：

式中：b为裂缝开度，μm。

考虑到裂缝的各向异性特征，基于三维空间转换原理,得到单元体内任意一条倾角为α_i、方位角为β_i、内蕴渗透率为K₀的裂缝在大地坐标系中的渗透率表达式为：

将裂缝归位后，网格块的渗透率为：

式中：K_x0，K_y0，K_z0为裂缝在x，y和z方向的绝对渗透率，D；α_i为裂缝倾角，(°)；β_i为裂缝方位角，(°)。

离散缝洞网络中，每一个溶洞类似一个源汇相，溶洞中的流体依靠弹性作用流入或者流出溶洞，因此引入弹性储集系数表征单位压力变化时，流入或者流出溶洞流体的量为：

式中：q_L3为3号溶洞依靠弹性流进或流出流体的量，m³；E_L3为3号溶洞弹性储集系数，m³/MPa；pⁿ₃和pⁿ⁺¹₃为3号溶洞相邻时刻内的压力，MPa。

溶洞的弹性储集系数与其中流体和岩体的弹性有关，可以参考多孔介质综合压缩系数的概念进行计算。根据缝洞型油藏储集空间差异分类，无论溶洞大小的差异性如何，流体在其中的分布仍然遵循重力分异原则，因此当有流体流进或流出时，其中的流体会按照密度差异重新分布，如图 3所示。

图 3 溶洞内流体分布状态 Fig.3 Distribution state of fluids in a karst

图选项

在缝洞型油藏中同时存在裂缝-裂缝或裂缝-溶洞连通形式，规整后的连通形式基本不变，如图 4所示。

图 4 缝洞连通形式示意 Fig.4 Schematic diagram for fracture-karst connection form

图选项

如果缝洞连通形式如图 4(a)所示，以溶洞为节点，则根据流体流动的连续性，可以得到下述质量守恒方程：

如果缝洞连通形式如图 4(b)所示，以裂缝为节点，则根据流体流动的连续性，可以得到下述质量守恒方程：

式中：J_Lk→1为节点1与周围节点的液相连通流动系数，m³/(MPa·d)；pⁿ⁺¹_k为当前时刻周围节点的压力，MPa；pⁿ⁺¹₁为当前时刻溶洞和周围节点的压力，MPa；Δtⁿ⁺¹为当前时间步长，d。

对于具有统一水动力系统的缝洞型油藏，缝洞系统具有相同的油水界面。但是，在缝洞型油藏中，由于溶洞是以离散孤立形式存在的，且由于裂缝与溶洞连通位置的结构不同，部分溶洞内存在各自的油水界面，这些溶洞内油水同时存在或全部为水，如图 5所示。

图 5 初始流体分布模式 Fig.5 Distribution mode of initial fluids

图选项

流体流动过程中，缝洞节点的含水率和含油率是不断变化的，这种变化又决定了其中的流体从该节点流入下一节点的动态。周围节点中的流体通过裂缝流入缝洞节点时，缝洞节点的含水率为：

式中：q_ok1为节点k流入节点1的油量，m³/d；J_ok1为节点k与节点1的油相流动系数，m³/(MPa·d)；p_k为节点k处的压力，MPa；q_wk1为节点k流入节点1的水量，m³/d；J_wk1为节点k与节点1的水相流动系数，m³/(MPa·d)；f_w1为节点1的含水率；Q_iw1为节点1作为源汇项时流入的水量，m³/d；Q_io1为节点1作为源汇项时流入的油量，m³/d。

流体在裂缝中的流动能力与流体性质有较大关系，笔者基于理想裂缝油水相渗曲线，并利用含水率、含油率加权得到油水共同流动时，裂缝连通流动系数的表达式：

式中：J_L12为节点1与节点2的液相流动系数，m³/(MPa·d)；f_o12为流出节点1与节点2的含油率；J_io12为节点1与节点2的油相流动系数，m³/(MPa·d)；f_w12为流出节点1与节点2的含水率；J_iw12为节点1与节点2的水相流动系数，m³/(MPa·d)。

同理，可得溶洞节点弹性储集能力的的表达式：

式中：E_io3为3号溶洞全部储集油时的弹性储集系数，m³/MPa；f_o3为3号溶洞的含油率；E_iw3为3号溶洞全部储集水时的弹性储集系数，m³/MPa；f_w3为3号溶洞的含水率。

由于井为流入或流出节点，无论井点贯穿裂缝，还是同时贯穿裂缝和溶洞，都可以将井看作溶洞体，如图 6所示。

图 6 井点虚拟溶洞示意 Fig.6 Schematic diagram for well point virtual karst

图选项

对于井所在节点，根据流体流动的连续性可以得到定液量和定流压生产时的质量守恒方程。

定液量生产时：

定流压生产时：

式中：qⁿ⁺¹_L1为井定压p₁时的井点产量，m³/d；pⁿ⁺¹₁为井定产q_L1时的井点压力，MPa。

根据质量守恒方程，可以建立缝洞系统压力节点数值模型，然后采用预处理共轭梯度方法求解高阶系数矩阵。预处理共轭梯度方法把矩阵的不完全分解和共轭梯度法结合起来，计算速度快。

笔者利用TK472C单井单元和TK630多井单元生产动态验证模型的可靠性。首先，根据TK472C单井单元和TK630多井单元钻井放空信息、钻井液漏失信息、地球物理识别预测结果以及单元地质储量，确定裂缝、溶洞的位置以及溶洞的大小；然后，根据井点裂缝统计信息，利用蒙特卡洛方法，确定裂缝的长度、开度、走向等参数，从而生成多个满足地球物理信息以及井点地质统计信息的地质模型；最后，利用上述数学模型计算生产不同时间的含水率，并与实际生产数据对比，选择拟合效果较好的缝洞地质模型。图 7为TK472C单井单元含水率拟合曲线。图 8为TK630多井单元的含水率拟合曲线。从图 7和图 8可以看出，所建模型的计算结果与生产数据拟合较好。

图 7 TK472C井含水率拟合曲线 Fig.7 Fitted curve for water cut of Well TK472C

图选项

图 8 TK630多井单元含水率拟合曲线 Fig.8 Fitted curve for water cut of TK630 Multi-well Unit

图选项

笔者应用上面历史拟合得到的缝洞地质模型分析注水参数对缝洞型油藏开发效果的影响。

首先，以TK472C单井模型为基础分析注采比对缝洞型油藏开发效果的影响。TK472C单井模型有强底水侵入，弹性能量较强，其在不同注采比(R_IP)下的含水率变化如图 9所示。从图 9可以看出，注采比越高含水率上升越快，可见，对于具有强底水的缝洞型油藏，一方面注入水抑制了强底水能量的发挥，另一方面注入水的作用范围通常比底水驱动小，尚未发挥作用便被采出。因此，对于强底水缝洞型油藏依靠底水驱动开发效果更好。

图 9 强底水缝洞油藏不同注采比下的含水率变化 Fig.9 Change of water cut under different injection-production ratios of fracture-karst reservoir with strong bottom water

图选项

仍以TK472C单井缝洞模型为基础，假设油藏有弱底水侵入，能量为强底水模型的四分之一，其在不同注采比下的含水率变化如图 10所示。从图 10可以看出，早期阶段，注采比为0时含水率较低，但是注水生产一段时间后的含水率高于注水比不为0的情况，说明弱底水油藏中，注入水的作用范围增大，适当注入一定量的水可以增加地层能量，能在一定程度上抑制底水锥进。因此，注水对于弱底水缝洞型油藏是有效的。

图 10 弱底水缝洞油藏不同注采比下的含水率变化 Fig.10 Change of water cut under different injection-production ratios of fracture-karst reservoir with weak bottom water

图选项

以TK630多井缝洞单元模型为基础，缝洞单元无底水，弹性能量较强，分析注水时机、采液速度对开发效果的影响。图 11为不同采出程度(R)下转注水时的含水率变化曲线。从图 11可以看出：随着注水时间的推迟，含水上升速度逐渐增大，当采出程度为8%时转注水，含水率上升速度早期较快但后期变慢；当采出程度为12%时转注水，含水率上升速度始终比较高。考虑到经济效益，对于TK630多井单元，当采出程度达到8%~10%时转注水的开发效果较好。

图 11 不同转注时机下的含水率变化 Fig.11 Change curve for water cut under different water injection occasions

图选项

采液速度为年采液量占地质储量的百分数。图 12为不同采液速度下的含水率变化曲线。

图 12 不同采液速度下的含水率变化 Fig.12 Water cut change under different fluid production rates

图选项

从图 12可以看出，随着采液速度的增大，含水率上升速度逐渐加快；采液速度低于3%时，含水率上升速度差异不大，采液速度超过3%以后，含水率上升速度加快。因此，对于TK630多井单元，以3%的采液速度生产，注水开发效果较好。

1) 以缝洞型油藏地球物理识别预测结果以及溶洞、裂缝物性参数统计特征为依据，建立规整化的离散缝洞分布模式，通过引入等效连通系数、弹性储集系数建立了缝洞单元离散数值模拟方法。该方法可用于碳酸盐岩缝洞型油藏的数值模拟。

2 ) 强底水缝洞单元应依靠底水驱动方式进行开采，弱底水缝洞单元可以通过注入一定量的水增加地层能量，在一定程度上抑制底水锥进，从而取得较好的开发效果。

3) 对于TK630多井缝洞单元，采出程度达到8%~10%时转注水的开发效果较好；以3%的采液速度生产，注水开发效果较好。