注水开发油藏、储层含水饱和度和井口含水率是评价油藏储量、油层水淹程度，进行油藏监测和管理的重要指标^[1-4]。

含水饱和度获取方法主要有密闭取芯资料法^[5-6]、地球物理测井法^[7]、油藏数值模拟法和单井化学示踪剂法^[8-11]。密闭取芯费用昂贵，且数据点有限；地球物理测井法通常采用经验关系式或者统计回归技术，通过测井信息来预测，计算精度和效果难以满足生产实践要求；油藏数值模拟工作量大，计算周期长；而单井化学示踪测试确定含水饱和度，是利用示踪剂在油层中的色谱分离原理，取决于示踪剂在平衡状态下的色谱效应，测试结果精度高，可成本也较高。

含水率预测方法目前主要有数值模拟法、经验公式法和解析法。含水率预测困难，受储层物性及流体等多因素影响^{[1, 12]}。数值模拟方法是一种较好的预测方法，但需花费很多时间^[13]。经验公式法只是对数据进行数学方法上的处理，缺乏具体物理意义，包括结合井的动、静态资料，利用神经网络法^[14]，应用Usher模型^[15]，Weibull预测模型、广义翁氏预测模型和瑞利预测模型^[16]等建立含水率预测模型；利用水驱油实验数据拟合含水率变化曲线^[17]，建立含水率预测模型；基于实际聚合物驱含水率变化曲线的特征及定量表征模型的建立原则，建立含水率变化定量表征模型^[18]。解析法包括基于分流量理论，结合相渗曲线预测含水率^[19-22]；1999年蒋明^[23]、2007年Kazeem ^[3]和2011年Li等^[2]对油藏做一定假设后，建立含水率与时间的数学关系模型；2010年许寒冰等^[24]对注采井网内复势的研究和流线分布的刻画建立含水率解析解。这些方法均未考虑分层注水开发条件下油藏含水率的预测。

笔者应用井温曲线确定注水剖面，应用Poisson过程分析法，结合概率论方法，计算分层注水开发油藏任意时刻各小层中含水饱和度、见水时间和含水率。为油藏开发采取分层堵水措施等提供理论依据。

1 模型建立及求解 1.1 注水剖面确定

分层注水开发油藏，可通过注水井的温度剖面，依据注入水、储层岩石的传质传热基本方程，确定注水剖面，即确定每一小层的注水量^[25-26]。井眼质量和热量传递示意图见图 1。

1) 水相传质传热方程。

井眼中水温变化：

$ \begin{array}{l} {\rho _{\rm{w}}}{C_{\rm{w}}}\frac{{\partial {T_{\rm{w}}}}}{{\partial t}} = {\rho _{\rm{w}}}{C_{\rm{w}}}v\frac{{\partial {T_{\rm{w}}}}}{{\partial z}} + {k_{\rm{w}}}\frac{{{\partial ^2}{T_{\rm{w}}}}}{{\partial {z^2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{2h}}{{{r_{\rm{w}}}}}\left[ {\left( {{r_{\rm{w}}},z,t} \right) - {T_{\rm{w}}}\left( {z,t} \right)} \right]。\end{array} $

(1)

初始条件：

$ {T_{\rm{w}}}\left( {z,0} \right) = a + bz, $

(2)

边界条件：

${T_{\rm{w}}}\left( {{z_{\max }},t} \right) = {T_{\rm{s}}}\left( t \right),0 < t \le {t_{si}};$

(3)

$ {k_{\rm{w}}}\frac{{\partial {T_{\rm{w}}}\left( {{z_{\max }},t} \right)}}{{\partial z}} = 0,{t_{si}} \le t; $

(4)

$ {k_{\rm{w}}}\frac{{\partial {T_{\rm{w}}}\left( {0,t} \right)}}{{\partial z}} = 0,0 < t。$

(5)

其中：

$ h = \alpha {\left[ {T\left( {{r_{{\rm{in}}}},z,t} \right) - {T_{\rm{w}}}\left( {z,t} \right)} \right]^{\frac{1}{4}}}。$

方程(1)是描述井眼中水相流量与温度变化的方程。当对流项${\rho _{\rm{w}}}{C_{\rm{w}}}v\frac{{\partial {T_{\rm{w}}}}}{{\partial z}} = 0$时，也适用于关井后的情况。

2) 岩层传质传热方程。

储层温度变化：

$ \begin{array}{l} {\rho _{\rm{e}}}{C_{\rm{e}}}\frac{{\partial T}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{k_{fe}}\frac{{\partial T}}{{\partial r}}} \right) + \frac{\partial }{{\partial z}}\left( {{k_{fe}}\frac{{\partial T}}{{\partial z}}} \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{q_{\rm{w}}}{\rho _{\rm{w}}}{C_{\rm{w}}}}}{{2{\rm{ \mathsf{ π} }}r}}\frac{{\partial T}}{{\partial r}} + \frac{{{q_{\rm{w}}}\left( {{p_{\rm{w}}} - {p_{\rm{e}}}} \right)}}{{2{\rm{ \mathsf{ π} }}{r^2}{J_c}\ln \left( {{r_{\rm{e}}}/{r_{\rm{w}}}} \right)}}。\end{array} $

(6)

初始条件：

$ T\left( {r,z,0} \right) = a + bz, $

(7)

边界条件：

$ T\left( {{r_{\rm{e}}},z,t} \right) = a + bz, $

(8)

$ T\left( {r,0,t} \right) = a, $

(9)

$ T\left( {r,{z_{\max }},t} \right) = a + b{z_{\max }}。$

(10)

其中：${k_{{\rm{fe}}}}\frac{{\partial T}}{{\partial r}}\left( {{r_{\rm{w}}},z,t} \right) = h\left[ {{T_{\rm{w}}}\left( {z,t} \right) - T\left( {{r_{\rm{w}}},z,t} \right)} \right]$，(0≤z≤z₁，z₂≤z≤z_max，0＜t；z₁＜z＜z₂，t_si＜t)；T(r_w，z，t)=T_w(z，t)，(z₁＜z＜z₂，0＜t≤t_si)；J_c为功热换算系数，427 m·kg/kJ。

方程(6)是描述井眼周围储层岩石的热量平衡，其适用于注水时的注水层；当径向对流$\frac{{{q_{\rm{w}}}{\rho _{\rm{w}}}{C_{\rm{w}}}}}{{2{\rm{ \mathit{ π} }}r}}\frac{{\partial T}}{{\partial r}} = 0$和摩擦热$\frac{{{q_{\rm{w}}}\left( {{p_{\rm{w}}} - {p_{\rm{e}}}} \right)}}{{2{\rm{ \mathit{ π} }}{r^2}{J_c}\ln \left( {{r_{\rm{e}}}/{r_{\rm{w}}}} \right)}} = 0$时，也适用于注水层外的其他岩层及关井以后的注水层。

结合水相和岩层初始条件和边界条件，联立式(1)和式(6)，就可通过有限差分计算出任意时刻各小层注水量q_w (1，j，t)。具体求解过程：先给定注水剖面q_w(1，j，t)的数值量，通过模型求出T_1，j^t数值解，然后分析T_1，j^t数值解和井温T_1，j^t测量值之间的差别，逐步调整注水剖面q_w (1，j，t)的数值量，使T_1，j^t数值解逼近测量值，这时注水剖面q_w (1，j，t)的数值量就可认为t时刻各小层的实际注水量。其中q_w (1，j，t)为第t时刻径向第1网格纵向第j网格注水量，m³/d；T_1，j^t为第t时刻径向第1网格纵向第j网格井温，℃。

1.2 水流线速度推导

注水开发油藏，注入水在储层中的推进速度与孔隙结构密切相关。以岩心做注水驱油实验，如图 2，岩心孔隙度为φ，横截面为A，单位时间注水量为Q，压力梯度为dp/dl，单位时间内水淹岩心体积为Q/φ，则水在岩心中推进速度为Q/(Aφ)。

设岩心孔隙直径分为N级：γ₁，γ₂，…，γ_N，所占比例分别为p(γ_j)=p_j，其中j=1，2，…，N；p₁+p₂+…+p_N=1。假设平行于岩心轴线的任一直线穿过直径为γ_j的孔隙个数满足强度λ_j的Poisson过程^[27]，即假定岩心长为L，任一长度l：0≤l≤L，穿过直径为γ_j的孔隙个数为随机变量ζ_j(l)，其概率分布为

$ P\left( {{\zeta _j}\left( l \right) = k} \right) = \frac{{{{\left( {{\lambda _j}l} \right)}^k}}}{{k!}}{{\rm{e}}^{{\lambda _j}l}},k = 0,1, \cdots 。$

(11)

穿过长度l各级直径的孔隙总个数为随机变量$\eta \left( l \right) = \sum\limits_{j = 1}^N {{\xi _j}\left( l \right)} $，由Poisson再生性知，η(l)仍满足强度$\lambda = \sum\limits_{j = 1}^N {{\lambda _j}} $的Poisson过程，其概率分布为

$ \begin{array}{l} P\left( {\eta \left( l \right) = k} \right) = \frac{{{{\left( {\lambda l} \right)}^k}}}{{k!}}{{\rm{e}}^{\lambda l}} = \\ \frac{{{{\left[ {\left( {{\lambda _1} + {\lambda _2} + \cdots + {\lambda _N}} \right)l} \right]}^k}}}{{k!}}{{\rm{e}}^{\left( {{\lambda _1} + {\lambda _2} + \cdots + {\lambda _N}} \right)l}}。\end{array} $

(12)

强度参数λ为平行于岩心轴线的任一直线上单位距离所穿过的孔隙平均总个数，可通过岩心切片或电镜照片观测统计得到，也可通过数学推导方法确定。

$ \tau = \frac{{{l_{孔}}}}{L} = \sqrt[3]{\varphi }, $

(13)

其中τ为线孔隙度，l_孔为岩心长度L孔隙所占的总长度。

孔隙平均直径为

$ \gamma = \sum\limits_{j = 1}^N {{\gamma _j}{P_j}} 。$

(14)

单位长度穿过孔隙的平均个数λ为

$ \lambda = \frac{\tau }{\gamma } = \frac{{\sqrt[3]{\varphi }}}{\gamma }。$

(15)

水驱岩心实验中单位时间内水淹孔隙的平均个数为

$ n = \lambda \cdot \frac{Q}{{A\varphi }}。$

(16)

这些孔隙所连成的平均距离为

$ l = n \cdot \gamma = \frac{{\lambda Q\gamma }}{{A\varphi }}。$

(17)

则水在岩心推进的平均速度可表示成

$ \bar \upsilon = \frac{{\frac{{\lambda Q\gamma }}{{A\varphi }}}}{{{\varphi ^{\frac{1}{3}}}}} = \frac{{\lambda Q\gamma }}{{A{\varphi ^{\frac{4}{3}}}}}。$

(18)

在实际油藏开发中，注入水沿各流线的推进速度并不都是υ，而是随机变量U_w。同样，水流线上单位时间内水淹孔隙的总个数也为随机变量，记为ζ_w，都服从Poisson分布，即

$ P\left( {{\zeta _{\rm{w}}} = k} \right) = \frac{{{{\left[ {\frac{{\lambda Q}}{{A\varphi }}} \right]}^k}}}{{k!}}{{\rm{e}}^{\frac{{\lambda Q}}{{A\varphi }}}}。$

(19)

并且有

$ {U_{\rm{w}}} = {\zeta _{\rm{w}}}\frac{\gamma }{{\sqrt[3]{\varphi }}}, $

(20)

$ P\left( {{U_{\rm{w}}} = k\frac{\gamma }{{\sqrt[3]{\varphi }}}} \right) = \frac{{{{\left[ {\frac{{\lambda Q}}{{A\varphi }}} \right]}^k}}}{{k!}}{{\rm{e}}^{\frac{{\lambda Q}}{{A\varphi }}}}。$

(21)

式(21)描述岩心各处水推进速度的随机性速度，其平均速度为

$ \begin{array}{l} {\rm{E}}\left( {{U_{\rm{w}}}} \right) = \sum\limits_{k = 0}^ \propto {k\frac{\gamma }{{\sqrt[3]{\varphi }}}P\left( {{U_{\rm{w}}} = k\frac{\gamma }{{\sqrt[3]{\varphi }}}} \right)} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\gamma }{{\sqrt[3]{\varphi }}}\frac{{\lambda Q}}{{A\varphi }} = \frac{{\lambda Q\gamma }}{{A{\varphi ^{\frac{4}{3}}}}} = \frac{Q}{{A\varphi }}。\end{array} $

(22)

这与实验中水在岩心中宏观推进速度Q/(Aφ)的推导结果是一致的。

对应方差为

$ {\rm{D}}\left( {{U_{\rm{w}}}} \right) = \frac{{{\gamma ^2}}}{{{\varphi ^{\frac{2}{3}}}}}\frac{{\lambda Q}}{{A\varphi }} = \frac{{\lambda Q{\gamma ^2}}}{{A{\varphi ^{\frac{5}{3}}}}}。$

(23)

假设任一平行于岩心轴线的直线l所穿过的孔隙总个数及其单侧的岩石骨架的长度之和为随机变量Y，其数学期望和方差分别为E(Y)、D(Y)，以Y_i表流线上第i个孔隙与其一侧的岩石骨架长度之和，Y_i与Y相互独立且相同分布，则水淹速度U_w可表示为

$ {U_{\rm{w}}} = \sum\limits_{i = 1}^{\left[ {\frac{{\lambda Q}}{{A\varphi }}} \right]} {{Y_i}} 。$

(24)

U_w均值为

$ {\rm{E}}\left( {{U_{\rm{w}}}} \right) = \frac{{\lambda Q}}{{A\varphi }}{\rm{E}}\left( Y \right)。$

(25)

对应方差为

$ {\rm{D}}\left( {{U_{\rm{w}}}} \right) = \frac{{\lambda Q}}{{A\varphi }}{\rm{D}}\left( Y \right)。$

(26)

对任一时间段t内的水淹距离

$ U\left( t \right) = \sum\limits_{i = 1}^{\left[ {\frac{{\lambda Q}}{{A\varphi }}t} \right]} {{Y_i}} , $

(27)

这里$\left[ {\frac{{\lambda Q}}{{A\varphi }}} \right]$、$\left[ {\frac{{\lambda Q}}{{A\varphi }}t} \right]$是指与${\frac{{\lambda Q}}{{A\varphi }}}$、 最接近的整数。

则U(t)均值为

$ {\rm{E}}\left[ {U\left( t \right)} \right] = \frac{{\lambda Q}}{{A\varphi }} \cdot t \cdot {\rm{E}}\left( Y \right)。$

(28)

对应方差为

$ {\rm{D}}\left[ {U\left( t \right)} \right] = \frac{{\lambda Q}}{{A\varphi }} \cdot t \cdot {\rm{D}}\left( Y \right)。$

(29)

1.3 含水饱和度计算

假设某油藏以五点法注采井网开采，储层中水从注入井到生产井推进过程如图 3，设A为注入井，B为生产井，建立坐标系0~x~y。由对称性，仅分析正方形对角线OB下方的三角形区域。

式(30)函数族的曲线与按等产量—源—汇的流线方程(31)

$ y = {x^\alpha },\left( {0 \le x \le 1,参变量\;\alpha \ge 1} \right), $

(30)

$ {\rm{arctg}}\frac{y}{{x - a}} - {\rm{arctg}}\frac{y}{{x + a}} = c\left( {c\;常数} \right) $

(31)

的流线非常接近，因此将三角形区域内的流线近似地以函数y=x^α表示。当α=1，是正方形对角线，长为$\sqrt 2 $；当α→∞，对应两条直角边连成的折线，长为2，任一值α，所对应流线长度

$ {l_\alpha } = \int_0^1 {\sqrt {1 + {{\left[ {\frac{{{\rm{d}}{x^\alpha }}}{{{\rm{d}}x}}} \right]}^2}} {\rm{d}}x} = \int_0^1 {\sqrt {1 + {\alpha ^2}{x^{2\left( {\alpha - 1} \right)}}} {\rm{d}}x} 。$

(32)

适当选取一组α值α₁＞α₂＞…＞α_n，则可绘出n条流线

$ y = {x^{{\alpha _i}}},i = 1,2, \cdots ,n。$

(33)

当注入井压力P_e和生产井井底流压P_w都稳定时，注入水沿长度为l_α流线的压力梯度为

$ {\rm{grad}}{P_\alpha } = \frac{{{P_{\rm{e}}} - {P_{\rm{w}}}}}{{{l_\alpha }}}。$

(34)

结合达西定律

$ Q = A\frac{K}{\mu }{\rm{grad}}P, $

(35)

将式(35)代入式(28)，则经时间t水淹推进平均距离为

$ {\rm{E}}\left[ {U\left( t \right)} \right] = \frac{{\lambda Kt}}{{\mu \varphi }}{\rm{E}}\left( Y \right) \cdot {\rm{grad}}P。$

(36)

同理，水淹推进距离的方差为

$ {\rm{D}}\left[ {U\left( t \right)} \right] = \frac{{\lambda Kt}}{{\mu \varphi }}{\rm{D}}\left( Y \right) \cdot {\rm{grad}}P。$

(37)

根据式(36)和(37)可计算出任意时刻不同流线的水淹推进距离和对应方差，据此得到t时的各流线水淹前缘位置，依次连接各流线水淹前缘位置就可画出t时油水界面曲线，此油水界面曲线包围的面积与总面积之比就是t时储层含水饱和度。

1.4 含水率计算

对于同一小层，根据达西定律，流量与压力梯度成正比，将各流线压力梯度累积求总和ΔP_k，记为

$ \Delta {P_k} = \sum\limits_{i = 1}^{{n_k}} {{\rm{grad}}{P_k}{\alpha _i}} , $

(38)

其中k为层号，n_k为第k层流线的总条数，gradP_kα_i为第α_i条流线的总压降。

假设在t时刻共有n(t)条流线水驱前缘已到达生产井，此时这些流线已全部为水相，将这些流线上压力梯度累积求和为

$ \Delta {P_{k水}} = \sum\limits_{i = 1}^{{n_k}\left( t \right)} {{\rm{grad}}{P_k}{\alpha _i}} 。$

(39)

则此时第k小层含水率f_kw(t)为

$ {f_{k{\rm{w}}}}\left( t \right) = 100\% \cdot \frac{{\Delta {P_{k水}}}}{{\Delta {P_k}}}。$

(40)

将各小层含水率与对应产量加权平均，得到t时刻的生产井含水率f_w(t)为

$ {f_{\rm{w}}}\left( t \right) = \frac{{\sum\limits_{i = 1}^k {{q_k} \cdot {f_{k{\rm{w}}}}\left( t \right)} }}{{{Q_{\rm{L}}}}}, $

(41)

其中q_k为第k层产液量，Q_L为生产井总产液量。

2 模型应用

某注水开发油藏，井深3 980 m，控制半径88 m，原始地层压力P_i=40.8 MPa，生产井中程深度压力为32 MPa，套管压力43.2 MPa，套管直径14.6 cm，井距200 m采用五点法注采井网，注入压力P_e=45 MPa，日注水量150 m³，注入时间5 610 h，注入水温度18.5 ℃，注水层面：L₁为3 562.3~3 609.6 m，L₂为3 609.6~3 668.1 m，L₃为3 668.1~3 726.6 m，L₄为3 726.6~3 785.2 m，L₅为3 785.2~3 852.2 m，注入水的黏度μ_w=1 m·Pa·s，油的黏度μ_o=1.5 m·Pa·s，原始含油饱和度S_oi=0.7，储层平均孔隙度φ=11.25%，渗透率K=32.51×10^-3 μm³。储层孔隙体系分布如表 1。

2.1 注水剖面

依据1.1小节，基于水相和岩石的传质传热方程，给定一初始注水剖面，通过有限差分求解方法，求取井温剖面数值解，通过调整注水井的注水剖面，使得井温剖面数值解接近实测值，最终确定了注水剖面，结果见图 4。

注水井总共分5小层注水，各层注水量分别为Q₁=24 m³，Q₂=28 m³，Q₃=35 m³，Q₄=32 m³，Q₅=31 m³。

2.2 含水饱和度

基于1.1小节已确定各层注水量，根据1.2小节各层水流线速度推导，计算出任意时刻不同流线的水淹推进距离，据此就可画出任意时刻油水界面曲线，从而得到任意时刻各层含水饱和度，见图 5~7。

从图 5可知，含水饱和度上升最快的是第3层，其注水强度最大为0.648 m³/(d·m)，其次为第4层，上升最慢的是第1层，注水强度也最小，为0.508 m³/(d·m)。

从图 6可知，随注水开发进行，第1小层油水界面不断扩展，到了第1 100 d水已在生产井发生突破。

从图 7可知，在第600 d，第3层油水界面推进最快，第1层油水界面推进最慢，与图 6结果一致。

2.3 含水率

基于以上计算结果，根据1.4小节，可计算各层的见水时间及任意时刻的含水率，见图 8~10。

见水最早的是第3层，为779 d，其含水上升率最快；其次是第4层为884 d，最晚见水为第1层，1 078 d，其含水上升率最慢。1 600 d之后，各层含水率都高达80%以上，趋于一致。

见水之后，各层含水率在初期上升较快，之后慢慢减缓。第3层含水率上升起始时刻早，上升快，第1 250 d含水率为98.5%，水淹严重；第1层含水率上升起始时刻晚，上升慢，第1 850 d含水率为99.1%，几乎全部水淹。

按每小层产液量加权平均，计算得到生产井的含水率，将其与井口计量的含水率进行对比，见图 10。初期计算误差较大，后期较小，在5%之内。这是由于见水初期，水量大部分来自裂缝等高渗透通道，建立模型未考虑此影响，当常规层位逐渐见水后，此影响削弱，计算的含水率更接近实际，误差小。

3 结论

1) 建立了分层注水开发油藏水淹动态跟踪模型。基于传质传热原理，根据井温曲线确定注水剖面，应用Poisson过程分析法，结合概率论方法，计算注水开发油藏任一时刻各小层中含水饱和度、见水时间和含水率。

2) 分层注水开发油藏，各小层水淹动态过程不同。注水强度大，存在裂缝等高渗透通道的小层，含水饱和度增大快，见水时间早，且含水率上升快，反之。

3) 在不同开发时期，模型预测精度不同。由于模型考虑了储层的微观孔隙分布，但未考虑微裂缝等高渗透通道，故在见水初期，水主要来自高渗透通道时，计算含水率误差较大，随着常规储层逐步水淹，计算含水率越趋于实测含水率，计算误差小于5%。