0 引言

对于水驱气藏而言，在其开发过程中随着气藏压力的下降，边（底）水逐渐推进至井底，导致气井见水，产水量迅速增加，产气量大幅下降，甚至部分气井出现水淹的状况，从而影响了气藏开发效果^[1-2]。水驱气藏在我国气藏类型中占有相当大的比例（约占65%），正确认识水驱气藏的水侵动态变化规律及控制因素，是成功开发该类气藏的关键^[3-4]。目前，采用水驱特征曲线^[5-6]、物质平衡法^[7-8]、经验公式^[9]等方法研究水驱气藏各参数的动态大小及变化规律已有较多研究成果，但这些方法仅能笼统分析水侵量、累计产气量和累计产水量随时间的变化^[10-13]，难以体现复杂的地下渗流及驱替特征，也无法了解不同时刻地层水的推进情况及水侵动态的控制机制。一些研究人员^[14-15]通过气藏数值模拟手段，虽然可以分析不同开发阶段的水线推进情况来了解气藏的水侵动态变化，但是存在复杂渗流机理认识不深入的缺点，分析结果存在较多模糊性。

基于分流理论，提出一种新的改进分流方程。首先推导考虑气体高速非达西渗流的分流方程，以四川盆地某边水气藏为例，分析典型因素对边水气藏动态变化的影响规律；其次，结合正交实验设计与熵权法，明确边水气藏动态变化的主控因素；然后，利用新的改进分流方程和传统的基于达西定律的分流方程，分别计算多口产水气井的理论见边水时间，并将计算结果分别与矿场统计值进行对比与分析；最后，基于研究区水侵主控因素的分析结果，探讨延缓边水气藏水侵的技术思路与对策，以期为水侵动态分析、治水稳产方案的编制提供一定参考。

1 边水气藏分流方程推导 1.1 经典分流理论的适用性

Buckley等^[16]基于一维水驱油线性系统，提出了经典分流理论（图 1a），描述了油水两相非混相驱替过程，其中油水两相流体渗流过程满足达西定律。经典分流理论包含2个重要方程，即分流率方程和等饱和度推进方程。

驱替相分流率方程^[12]：

$f_{\mathrm{D}}=\frac{q_{\mathrm{D}}}{q_{\mathrm{T}}}=\frac{q_{\mathrm{D}}}{q_{\mathrm{D}}+q_{\mathrm{o}}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{D}}+v_{\mathrm{o}}}=\frac{1}{1+\frac{K_{\mathrm{o}} \mu_{\mathrm{D}}}{K_{\mathrm{D}} \mu_{\mathrm{o}}}}$

(1)

式中：f_D为驱替相分流率；q_D，q_o，q_T分别为驱替相、被驱替相及两相总体积流量，m³/s；v_D，v_o分别为驱替相与被驱替相渗流速度，m/s；K_D，K_o分别为驱替相与被驱替相有效（相）渗透率，mD；μ_D，μ_o分别为驱替相与被驱替相黏度，mPa·s。

等饱和度推进方程（或推进速度方程）^[16]：

$\Delta u=\frac{q_{\mathrm{T}}}{\varphi A}\left(\frac{\mathrm{d} f_{\mathrm{D}}}{\mathrm{d} S_{\mathrm{D}}}\right) \Delta \theta$

(2)

$\left(\frac{\partial u}{\partial \theta}\right)_{S_{\mathrm{D}}}=\frac{q_{\mathrm{T}}}{\varphi A}\left(\frac{\partial f_{\mathrm{D}}}{\partial S_{\mathrm{D}}}\right)_\theta$

(3)

式中：u为驱替相沿着渗流方向的推进距离，m；S_D为驱替相饱和度，%；θ为时间，ks；A为渗流截面积，m²；φ为孔隙度，%。

根据水体与天然气的分布可将水驱气藏分为边水气藏和底水气藏。边水气藏的流体性质和渗流规律与水驱油藏具有明显差异：①水驱油藏与边水气藏的高压物性参数不同，导致二者的驱替特征不同；②对于水驱油藏，渗流速度通常达不到紊流发生的速度界限，而对于边水气藏，当其整体采气强度或气井采气量提高时，气井容易在近井地带产生高速非达西效应，随之产生由流量变化引起的附加压降，流体的渗流规律偏离达西定律；③边水气藏径向系统下的基本几何流动模型与水驱油藏不同，径向系统下边水气藏流体的流动方向由供气边界流向井筒，基本几何流动模型为“点汇”模型（图 1b），而常规水驱油径向基本几何流动模型为“点源”模型，即流体的流动方向由注入端向外延伸（图 1c）。以上原因导致常规水驱油的分流方程难以应用于边水气藏。

1.2 改进的分流率方程

基于达西定律的水相运动方程^[15-17]（不考虑毛管力、重力）为

$v_{\mathrm{w}}=\frac{k k_{\mathrm{rw}}}{\mu_{\mathrm{w}}} \frac{\mathrm{d} p_{\mathrm{w}}}{\mathrm{d} r}$

(4)

式中：v_w为水相渗流速度，m/s；k为绝对渗透率，mD；k_rw为水相相对渗透率；p_w为水相流体压力，MPa；μ_w为水相黏度，mPa·s；r为径向距离，m。

考虑到径向地层气井高速生产过程中容易产生高速非达西效应，根据Forchheimer^[18]提出的二次方程，对于气相有

$\frac{\mathrm{d} p_{\mathrm{g}}}{\mathrm{d} r}=\frac{\mu_{\mathrm{g}} v_{\mathrm{g}}}{k k_{\mathrm{rg}}}+\beta \rho_{\mathrm{g}} v_{\mathrm{g}}{}^2$

(5)

对Forchheimer二次方程求根，物理上方程仅存在唯一解，可将式（5）进一步表示为

$v_{\mathrm{g}}=\frac{-\frac{\mu_{\mathrm{g}}}{k k_{\mathrm{rg}}}+\left[\frac{\mu_{\mathrm{g}}^2}{k^2 k_{\mathrm{rg}}^2}+4 \rho_{\mathrm{g}} \beta \frac{\mathrm{d} p_{\mathrm{g}}}{\mathrm{d} r}\right]^{1 / 2}}{2 \rho_{\mathrm{g}} \beta}$

(6)

式中：v_g为气相渗流速度，m/s；μ_g为气相黏度，mPa·s；k_rg为气相相对渗透率；ρ_g为气相密度，kg/m³；β为非达西流系数，m^-1；p_g为气相流体压力，MPa。

Evans等^[19]的研究表明，如果利用单相条件下的非达西流系数与绝对渗透率、孔隙度构建的经验方程，来预测气水两相渗流系统下的非达西流系数，会导致预测结果与实验结果之间出现数量级的偏差。由此提出采用气相有效渗透率来替换绝对渗透率，采用气相饱和度来修正孔隙度，所获经验关系更能反映气水两相渗流系统下的渗流规律及由惯性效应主导的物理规律。Evans等^[19]、LIU等^[20]将两相渗流系统下的非达西流系数称为“有效非达西流系数”（本文用“β_eff”表示），关系式为

$\beta_{\mathrm{eff}}=\frac{\mathrm{C}_\beta}{\left.\left(k k_{\mathrm{rg}}\right)^{5 / 4}\left[\varphi\left(S_{\mathrm{g}}-S_{\mathrm{gr}}\right)\right]\right]^{3 / 4}}$

(7)

式中：β_eff为有效非达西流系数，m^-1；C_β为非达西流常数，m^3/2，砂岩储层取值高于碳酸盐岩储层1~2个数量级^[20-21]，SI单位制下取值3.2×10^-7~3.2×10^-5 m^3/2；S_g，S_gr分别为含气饱和度和残余气饱和度，%。

计算气体高速非达西渗流条件下的各相分流率时，首先需要获得各相渗流速度。假设不考虑毛管压力，气相与水相具有相同的压力梯度，式（5）可变形为

$v_{\mathrm{g}}=\frac{k k_{\mathrm{rg}}}{\mu_{\mathrm{g}}} \frac{\mathrm{d} p}{\mathrm{~d} r}\left(1+\frac{\beta_{\mathrm{eff}} \rho_{\mathrm{g}} k k_{\mathrm{rg}} v_{\mathrm{g}}}{\mu_{\mathrm{g}}}\right)^{-1}$

(8)

令

$\delta=\frac{\beta_{\mathrm{eff}} \rho_{\mathrm{g}} k k_{\mathrm{rg}}}{\mu_{\mathrm{g}}}, \chi=\left(1+\delta v_{\mathrm{g}}\right)^{-1}$

(9)

联立式（4）、式（6）、式（7）—（9），考虑气相高速非达西渗流的分流率方程可写为

$f_{\mathrm{w}}^{\mathrm{h}}=\frac{v_{\mathrm{w}}}{v_{\mathrm{w}}+v_{\mathrm{g}}}=\frac{1}{1+\frac{k_{\mathrm{rg}}}{k_{\mathrm{rw}}} \frac{\mu_{\mathrm{w}}}{\mu_{\mathrm{g}}} \chi}$

(10)

$f_{\mathrm{g}}^{\mathrm{h}}=\frac{v_{\mathrm{g}}}{v_{\mathrm{w}}+v_{\mathrm{g}}}=\frac{1}{1+\frac{k_{\mathrm{rw}}}{k_{\mathrm{rg}}} \frac{\mu_{\mathrm{g}}}{\mu_{\mathrm{w}}} \chi^{-1}}$

(11)

式中：f_w^h，f_g^h分别表示考虑气体高速非达西渗流的水相与气相分流率。

1.3 饱和度推进方程、水侵前缘及见边水时间

由图 2可看出，在边水侵入的初始时刻（t=0），地层总体积流量（Q_t）等于水相体积流量（Q_w）；在边水逐渐侵入过程中，相邻微元体中的流体在某一时刻（t或某一时间段Δt内）满足质量守恒定律。

微元体满足如下质量守恒关系：

$\begin{gathered}{\left[\left(Q_{\mathrm{w}} \rho_{\mathrm{w}}\right)_{r+\Delta r}-\left(Q_{\mathrm{w}} \rho_{\mathrm{w}}\right)_r\right] \Delta t=\pi h \varphi\left\{\left[r_{\mathrm{e}}-(r+\Delta r)\right]^2-\right.} \\ \left.\left(r_{\mathrm{e}}-r\right)^2\right\}\left[\left(S_{\mathrm{w}} \rho_{\mathrm{w}}\right)_{t+\Delta t}-\left(S_{\mathrm{w}} \rho_{\mathrm{w}}\right)_t\right]\end{gathered}$

(12)

式中：Q_w为水相体积流量，m³/s；S_w为含水饱和度，%；r_e为供给边界，m；r为某一特定径向位置或推进距离，m；Δr为微元体径向距离变化量，m；ρ_w为水相密度，kg/m³；t为边水推进时间，s；Δt为边水推进Δr所需时间，s；h为储层有效厚度，m。

对式（12）化简可得

$\begin{gathered}{\left[\left(Q_{\mathrm{w}} \rho_{\mathrm{w}}\right)_{r+\Delta r}-\left(Q_{\mathrm{w}} \rho_{\mathrm{w}}\right)_r\right] \Delta t=\pi h \varphi\left(\Delta r^2-2 r_{\mathrm{e}} \Delta r+2 r \Delta r\right)} \\ {\left[\left(S_{\mathrm{w}} \rho_{\mathrm{w}}\right)_{t+\Delta t}-\left(S_{\mathrm{w}} \rho_{\mathrm{w}}\right)_t\right]}\end{gathered}$

(13)

对于微元体，由于Δr → 0，Δt → 0，有下式满足近似相等关系：

$\Delta r^2-2 r_{\mathrm{e}} \Delta r+2 r \Delta r \approx 2 r \Delta r-2 r_{\mathrm{e}} \Delta r$

(14)

将式（14）右项代入式（13）。假设水相密度不变，根据偏微分的定义，式（13）可写为

$2 \pi h \varphi\left(r_\mathrm{w}-r\right) \frac{\partial\left(S_\mathrm{w}\right)}{\partial t}=-\frac{\partial\left(Q_\mathrm{w}\right)}{\partial r}$

(15)

根据分流率基本定义，式（15）可表示为

$2 \pi h \varphi\left(r_{\mathrm{e}}-r\right) \frac{\partial\left(S_{\mathrm{w}}\right)}{\partial t}=-\frac{\partial\left(Q_{\mathrm{t}} f_{\mathrm{w}}^{\mathrm{h}}\right)}{\partial r}$

(16)

式中：Q_t为总体积流量，m³/s。

对于气井某一特定生产时刻，式（16）可进一步化简为

$\frac{2 \pi h \varphi\left(r_{\mathrm{e}}-r\right)}{Q_{\mathrm{t}}} \frac{\partial\left(S_{\mathrm{w}}\right)}{\partial t}=-\frac{\partial f_\mathrm{w}^{\mathrm{h}}}{\partial r}$

(17)

由于水相分流率仅是含水饱和度的函数^[15-16]，对上述偏微分方程利用链式求导法展开，可得到边水气藏径向地层中两相流动的渗流微分方程（或驱替方程）：

$\begin{aligned} \frac{2 \pi h \varphi\left(r_{\mathrm{e}}-r\right)}{Q_{\mathrm{t}}} \frac{\partial S_{\mathrm{w}}}{\partial t}+\frac{\partial f_{\mathrm{w}}^{\mathrm{h}}}{\partial S_{\mathrm{w}}} \frac{\partial S_{\mathrm{w}}}{\partial r}= \\ \frac{2 \pi h \varphi\left(r_{\mathrm{e}}-r\right)}{Q_{\mathrm{t}}} \frac{\partial S_{\mathrm{w}}}{\partial t}+{f_{\mathrm{w}}^{\mathrm{h}}}^{\prime} \frac{\partial S_{\mathrm{w}}}{\partial r}=0\end{aligned}$

(18)

式中：f_w^h′为水相分流率导数。

含水饱和度S_w是时间t和位置r的二元函数，即S_w=f（r，t）。选择地层流体等饱和度作为特征线，在特征线上，下式成立^[15-16]：

$\frac{\mathrm{d} r}{\mathrm{d} t}=-\frac{\frac{\partial S_{\mathrm{w}}}{\partial t}}{\frac{\partial S_{\mathrm{w}}}{\partial r}}$

(19)

将式（19）带入式（18），整理得到径向地层任意水相等饱和度推进速度方程为

$\frac{\mathrm{d}\left[r_{\mathrm{e}}{ }^2-\left(r_{\mathrm{e}}-r\right)^2\right]}{\mathrm{d} t}=\frac{{f_{\mathrm{w}}^{\mathrm{h}}}^{\prime} Q_{\mathrm{t}}}{\pi h \varphi}$

(20)

对式（20）关于时间t和径向距离r积分，得

$r_{\mathrm{e}}{ }^2-\left(r_{\mathrm{e}}-r\right)^2=\frac{Q_{\mathrm{t}} f_{\mathrm{w}}^{\mathrm{h}}{ ^{\prime}}}{\pi h \varphi}\left(t-t_0\right)$

(21)

式中：t₀为水侵初始时刻，s。

显然，式（21）仅存在唯一解满足物理规律，即边水气藏任意水相等饱和度面推进距离为

$r=r_{\mathrm{e}}-\left[r_{\mathrm{e}}{ }^2-\frac{Q_1 {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}}{\pi h \varphi}\left(t-t_0\right)\right]^{1 / 2}$

(22)

由于水侵流量在气藏开发过程中并非常数，平均水侵速度v_t可采用数值积分方法计算，即采用物质平衡法获得不同阶段的水侵流量Q_ti(t)，计算Q_ti(t)与dt围成的面积，并除以累计生产时间t，表达式为

$\overline{v_{\mathrm{t}}}=\frac{\int_0^t Q_{\mathrm{ti}}(t) \mathrm{d} t}{A t}=\frac{\int_0^t Q_{\mathrm{ti}}(t) \mathrm{d} t}{2 \pi r_{\mathrm{e}} h t}$

(23)

式中：v_t为平均水侵速度，m/s；Q_ti(t)为气井不同阶段生产阶段的水侵流量，m³/s。

联立式（22）、式（23），当t₀=0时，有

$\begin{align} r= & r_{\mathrm{e}}-\left[r_{\mathrm{e}}^2-\frac{\int_0^t Q_{\mathrm{ti}}(t) \mathrm{d} t}{\pi h \varphi} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}\right]^{1 / 2}= \\ & r_{\mathrm{e}}-\left[r_{\mathrm{e}}^2-\frac{2 \overline{v_{\mathrm{t}}} r_{\mathrm{e}} t}{\varphi} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}\right]^{1 / 2}\end{align}$

(24)

则井筒距水侵前缘的距离r_w-f为

$\begin{aligned} r_{\mathrm{w}-\mathrm{f}}= & {\left[r_{\mathrm{e}}^2-\frac{\int_0^t Q_{\mathrm{ti}}(t) \mathrm{d} t}{\pi h \varphi} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}\right]^{1 / 2}=} \\ & {\left[r_{\mathrm{e}}^2-\frac{2 \overline{v_{\mathrm{t}}} r_{\mathrm{e}} t}{\varphi} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}\right]^{1 / 2} }\end{aligned}$

(25)

式中：r_w-f为井筒距水侵前缘的距离，m。

通过以上推导，可得到边水气藏（或气井）径向地层任意水相等饱和度的线移动速度及推进方程[式（20）、式（22）或式（24）]，也可得到井筒距水侵前缘的距离[式（25）]。

当气相区增加的水量等于边水侵入量时有：

$\begin{gathered}\int_0^{r_1}\left(S_{\mathrm{w}}-S_{\mathrm{wi}}\right) \pi h \varphi \mathrm{d}\left[r_{\mathrm{e}}^2-\left(r_{\mathrm{e}}-r\right)^2\right]= \\ Q_{\mathrm{t}}\left(t-t_0\right)\end{gathered}$

(26)

式中：r_f为水侵前缘推进距离，m。

整理式（26）为

$\begin{gathered}\int_0^r\left(S_{\mathrm{w}}-S_{\mathrm{wi}}\right)\left(2 r_{\mathrm{e}}-2 r\right) \pi h \varphi \mathrm{d} r= \\ Q_{\mathrm{t}}\left(t-t_0\right)\end{gathered}$

(27)

对式（21）求微分，可得

$\left(2 r_{\mathrm{e}}-2 r\right) \mathrm{d} r=\frac{Q_{\mathrm{t}}\left(t-t_0\right)}{\pi h \varphi} \mathrm{d} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}$

(28)

将式（28）带入式（27），可得

$\left.\left(S_{\mathrm{w}}-S_{\mathrm{wi}}\right) {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}\right|{\;}_0 ^{r_{\mathrm{f}}}-\int_0^{r_\mathrm{f}} \mathrm{~d} f_{\mathrm{w}}^{\mathrm{h}}=1$

(29)

式中：S_wi为初始地层含水饱和度，%；$f_{\mathrm{wf}}^{\mathrm{h}}, {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}$分别为水侵前缘含水饱和度（S_wf）对应的分流率及分流率导数。

进一步积分，得到边水气藏（或气井）水侵前缘饱和度方程为

${f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}=\frac{f_{\mathrm{wf}}^{\mathrm{h}}}{S_{\mathrm{wf}}-S_{\mathrm{wi}}}$

(30)

式（30）的几何意义是由初始含水饱和度位置作水相分流率曲线的切线，切线的斜率值即${f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}$，切点位置对应的含水饱和度为水侵前缘饱和度。

水侵前缘到达井筒时，r_f=r_e-r_w且r_w²≪r_e²，由式（22）可得径向地层气井见水时间t_bt为

$\begin{gathered}t_{\mathrm{bt}}=t_0+\frac{\pi\left(r_{\mathrm{e}}{ }^2-r_{\mathrm{w}}{ }^2\right) h \varphi}{Q_{\mathrm{t}} f_{\mathrm{wf}}^{\mathrm{h}}} \approx \\ t_0+\frac{\pi r_{\mathrm{e}}{ }^2 h \varphi}{Q_{\mathrm{t}} {f_{\mathrm{w}}^{\mathrm{h}}}^{\prime}}\end{gathered}$

(31)

式中：t_bt为径向地层气井见水时间，s；r_w为井筒半径，m。

2 水侵动态影响因素分析与水侵对策研究 2.1 水侵动态影响因素分析

边水气藏水侵动态变化的影响因素较多，选取部分典型因素（如非达西流常数、润湿性、水侵流量）分析其对水侵动态参数的影响规律。

当非达西流常数（C_β）取值不同时，有效非达西流系数（β_eff）差距较大。在相同的水侵流量下（Q_t=0.40 m³/s），对C_β取不同数值，计算川中L02-X1井在稳定水侵流量下的水相分流率与水侵前缘推进情况，并与达西流条件下的曲线进行对比（图 3、图 4）。由图 3、图 4可知，在水侵流量相同的情况下，非达西流常数C_β越大，惯性效应越严重，气藏的“非活塞”驱替特征就越显著；水侵前缘饱和度越小，水侵速度越快，水侵前缘推进距离就越远。

对于某一特定气藏，润湿性/饱和历史对气水相渗曲线也具有显著影响，进而影响气藏水侵动态^[18-21]。计算不同润湿性（驱替与渗吸）条件下的分流率曲线与前缘推进曲线如5、图 6所示。

在不同润湿性条件下，多相流体的运动机制会发生变化。当储层润湿相饱和度较高时，非润湿相逐渐取代润湿相，驱替作用占主导；当储层润湿相饱和度较低时，润湿相逐渐取代非润湿相，渗吸作用占主导。由图 5、图 6可看出，不同润湿性条件下水相分流率曲线、水侵前缘推进曲线具有明显差异，对气藏的水侵动态具有一定影响。

以L01-1井为例，计算不同水侵流量下水相分流率及水侵前缘推进曲线如图 7、图 8所示。

由图 7、图 8可看出，水侵流量越大，水侵前缘饱和度越小，边水推进速度就越快。当水侵流量达到2.5 m³/s时，边水早已推进至井筒，且该水侵流量下的水侵前缘饱和度（18.82%）仅略大于初始地层含水饱和度，反映了更为显著的“非活塞”驱替特征（图 8）。

2.2 水侵动态影响因素敏感性分析

边水气藏水侵动态的影响因素较多，通过边水气藏水侵动态因素敏感性分析，可进一步明确边水气藏动态变化的主控因素，探究延缓水侵的技术思路与对策。

2.2.1 水侵动态影响因素与评价指标的选择及量化

选择储层孔隙度、渗透率、水侵流量、有效厚度、供气半径、非达西流常数、相渗曲线特征参数等作为水侵动态的影响因素，将水侵前缘推进速度与见边水时间作为水侵特征的评价指标。由于相渗曲线特征参数较多，若将表征同一类因素的所有特征参数纳入，会因多因素分析而难以形成规律性认识，由此提出采用相渗曲线“共渗三角区面积占比”来表征各类储层渗流能力的大小，并作为反映相渗特征的唯一因素纳入量表（表 1）。“共渗三角区面积占比”为相对渗透率曲线与饱和度轴围成的双曲边三角形闭合面积占总相渗曲线区域面积的比例，可结合改进Corey幂律模型^[19]，利用定积分计算，表达式为

$S_{\Delta}=\int_{S_{\mathrm{w}}}^{S_{\mathrm{w}}^{\mathrm{e}}} k_{\mathrm{rw}}\left(S_{\mathrm{w}}\right) \mathrm{d}\left(S_{\mathrm{w}}\right)+\int_{S_{\mathrm{w}}^{\mathrm{c}}}^{1-S_{\mathrm{gr}}} k_{\mathrm{rg}}\left(S_{\mathrm{w}}\right) \mathrm{d}\left(S_{\mathrm{w}}\right)$

(32)

式中：S_Δ为共渗三角区面积占比；S_wc为共渗点水相饱和度，%；k_rw（S_w），k_rg（S_w）为水相与气相相对渗透率与含水饱和度函数。

由表 1可以看出，储层渗透率在整个多元复杂系统中提供的信息最为关键，其对水侵动态的影响最为显著，其次为非达西流常数（或非达西流系数）、共渗三角区面积（或相对渗透率）、孔隙度，而有效厚度、供气边界、水侵流量的影响程度相同且较低。

水侵特征评价指标（见边水时间与水侵前缘推进速度）相较于各影响因素，其熵值低，在反映水侵动态的系统中具有较大信息量，在系统中权重系数总和高达49.18%，水侵特征评价指标与影响因素在整个系统中的地位同等重要，说明选用的评价指标可作为水侵变化的重要衡量指标。

依据岩石类型、孔隙特征及储层物性，可将研究区储层分为Ⅰ类（k ≥ 20 mD，φ≥12%）、Ⅱ类（0.1 mD ≤k<20 mD，6%≤φ<12%），Ⅲ类（0.01 mD≤k<0.1 mD，2%≤φ<6%））、Ⅳ类（k<0.01 mD，φ<2%）。由图 9可看出，从Ⅰ类储层到Ⅳ类储层，随着渗透率与孔隙度的减小，束缚水饱和度增大，气水两相区范围变窄，等渗点变低，等渗三角区面积变小，两相渗流能力变弱。

2.2.2 正交设计

基于研究区静态、动态资料，确定影响因素水平，设计L32.4.7/9正交表（32组算例、4个水平、7个影响因素、2个指标），制定标准正交表。

各影响因素的水平取值为：渗透率（40.00 mD，30.00 mD，10.00 mD，0.05 mD）；孔隙度（24%，18%，9%，4%）；有效厚度（100 m，75 m，50 m，25 m）；供气半径（2 000 m，1 500 m，1 000 m，500 m）；Q_t（1 m³/s，0.75 m³/s，0.5 m³/s，0.25 m³/s）；非达西流常数（0.2×10^-5 m^3/2，3.2×10^-5 m^3/2，6.2×10^-5 m^3/2，9.2×10^-5 m^3/2）。对于砂岩储层，水相与气相的Corey拟合指数n_w和n_g均取值为1.5~4.0^[19]。基于数值优化结果，制定相渗曲线综合表征参数（共渗三角区面积占比SΔ%）的水平为6.26%（n_g=n_w=1.5），3.75%（n_g=n_w=2.0），1.66%（n_g=2.0，n_w=4.0）和0.58%（n_g=n_w=4.0）。

2.3 延缓气藏水侵的技术思路与对策

储层渗透率、非达西流常数（或非达西流系数）、相对渗透率（共渗三角区面积）均是反映储层与流体渗流特性的重要影响因素，且各影响因素之间存在一定的相关性^{[15, 18-19, 22-24]}，而有效厚度、供气边界、水侵流量是反映储层与水体性质的宏观度量参数，并非地层流体渗流的决定性因素。

储层渗透率对边水气藏水侵动态的影响最为显著，渗透率发生变化时，非达西流常数与相对渗透率均会发生变化。通常，渗透率越大，边水推进速度越快，非均质性越强，边水非活塞驱替特征就越明显，采收率也就越低。因此，开展低渗致密有效储层识别与筛选，并结合水体分布及其性质认识，充分发掘物性分布相对均匀的低渗致密段储层的开发潜力是提升边水气藏开发效果的关键。

水侵流量尽管对水侵动态无直接相关关系，但是保持合理的储层压力可降低水侵流量，促进水线均匀推进，减少残余气局部富集，从而提高边水气藏最终采收率。因此，从生产角度，制定合理的气藏采气强度（或气井工作制度）是控制水侵速度及控水稳水的重要手段。

3 实例计算

四川盆地某边水气藏构造上为一近东西向的低幅度背斜，南北两翼稍陡，整体以薄层为主，平面连续性好，气藏温度为85 ℃，压力为34 MPa，储层孔隙度为4.0%~22.7%，渗透率为0.001~28.210 mD，非均质性强，气层有效厚度约30 m，天然气黏度为0.027 8 mPa·s，地层水黏度为0.357 1mPa·s，天然气密度为239.98 kg/m³，地层水密度为984.69 kg/m³，目前该气藏大量气井已见水，甚至部分气井出现水淹的状况。

利用新的改进分流方程对该气藏进行分析，首先统计该气藏10口气井的静态、动态资料，采用物质平衡法^{[2, 7, 10, 25]}，并结合数值积分法，计算各井的水侵流量；然后，利用Evans公式^[15]与改进Corey模型^[19]，采用数值迭代法计算各气井的分流参数及各井的理论见边水时间，并将新的改进分流方程和传统分流方程的理论见边水时间分别与矿场统计的见边水时间进行对比，进而分析水侵动态指标的差异。

3.1 计算水线推进情况与见边水时间

以L02-X1井为例，该井投产后产量相对稳定，2019年7月产气量迅速下降，产水量迅速上升。井筒测试表明，该气井井筒积液逐渐增多，在已有的工作制度下气井难以携液生产。2021年3月气井停产，通过优化气井工作制度，开展排水采气、油管解堵等措施，气井产能得到一定程度的恢复。目前，该井油压由初期的5.8 MPa下降至3.2 MPa，下降了44.82%。该气井生产曲线与产出水矿化度变化如图 10、图 11所示。

边水气藏开发初期产水主要为凝析水与原生地层水，气井见边水后，产水量急剧增大，产气量急剧减小，L02-X1井符合典型的边水气藏气井生产动态特征^{[2, 4]}（图 10）。由于地层水中通常具有多种离子（如氯离子、钠离子、钙离子、锶离子等），其矿化度远大于凝析水中的矿化度。根据产出水矿化度分析资料（图 11）可知，随着气藏不断被开采，产出水矿化度具有明显的上升趋势，见边水时间与生产曲线产水量迅速增加时间基本吻合。利用新的改进分流方程，计算L02-X1井理论见边水时间为1 100 d，矿场统计见边水时间为1 125 d（表 2），二者差距较小。

3.2 见边水时间计算值与矿场统计值对比

理论见边水时间与矿场统计见边水时间的平均相对误差（Mean Relative Error，MRE%）的计算公式为

$MRE(\%)=\frac{100}{N_i} \sum\limits_{i=1}^{N_i} \frac{\left|t^{\text {calc }}-t^{\text {real }}\right|}{t^{\text {real }}}$

(33)

式中：N_i为统计井数量，口；i为第i个统计数据；t^calc和t^real分别为计算的和矿场统计的见边水时间，d。

通过计算10口井的见边水时间（表 2）可知，新的改进分流方程理论见边水时间与矿场统计见边水时间的相对误差为-17.52%~15.82%，MRE%为9.77%。传统分流方程理论见边水时间与矿场统计见边水时间的相对误差为47.32%~163.62%，MRE%高达78.05%。因此，新提出的改进分流方程能较好地反映地层流体的渗流规律，但该分流方程对于常规储层的适应性更好，而对于非常规储层，由于存在特殊渗流机理，适应性仍需进一步研究。

4 结论

（1）通过扩展经典分流理论，推导了考虑气体高速非达西渗流的边水气藏分流方程，定量研究了不同开发时刻产水气井的边水推进状况及水侵动态参数变的化规律。新的改进分流方程相较于传统的基于达西定律的分流方程，计算结果更加可靠。

（2）渗透率对边水气藏水侵动态的影响最为显著，其次为非达西流常数、相对渗透率（共渗三角区面积）、孔隙度，而有效厚度、供气边界、水侵流量的影响程度均较低。

（3）充分发掘物性分布相对均匀的低渗致密段储层的开发潜力是提升边水气藏开发效果的关键，从生产角度，制定合理的气藏采气强度（或气井工作制度）是主动控水、稳水的重要手段。