閆正和 石軍太 秦 峰 洪舒娜 白美麗
(1. 中海石油深海開發有限公司 廣東珠海 518000; 2. 中國石油大學(北京)油氣資源與探測國家重點實驗室 北京 102249;3. 石油工程教育部重點實驗室(中國石油大學) 北京 102249 )
國內大多數常規氣藏都屬于不同程度的水驅氣藏,邊底水活躍的氣藏占40%~50%[1]。在水驅氣藏的開發過程中,隨著氣體的產出和地層壓力降低,水體逐漸侵入氣藏。水體的侵入,一方面可以有效補充地層能量,另一方面也會提高氣藏廢棄壓力、降低氣藏采收率。因此,水驅氣藏水侵量及儲量計算是氣藏合理高效開發的基礎,對水驅氣藏后期生產調整具有重要意義。
常規氣驅氣藏的動態儲量方法[2-5]已相對成熟,但在水驅氣藏儲量方面適應性較差,如應用常規物質平衡方法計算水驅氣藏動態儲量往往會導致計算結果偏大[6]。國內外學者針對水驅氣藏水侵量及動態儲量計算開展了大量研究,推導了若干水侵量計算方法,主要可分為穩態水侵模型[7]和非穩態水侵模型[8-11]兩類。但是上述模型假設條件較多,計算過程相對繁瑣,且計算所需的水體參數如水侵角、水體內外邊界半徑等往往難以準確確定,不利于實際應用。因此,一些學者基于物質平衡原理和最優化思想,提出了水侵量計算方法與數學模型[12-19]。如唐圣來 等[12]應用水驅氣藏物質平衡方程和變流壓下產量響應函數,通過擬合產量、流壓、靜壓等生產動態數據,經過多次循環迭代得出強水驅氣藏的動態儲量和水侵量;王怒濤 等[13]結合數值反演法,提出了一種自適應遺傳算法用于求解地質儲量、水侵量、水侵系數以及水體大小;杜凌云 等[14]針對常規物質平衡法中的差值法、圖版法、視地質儲量法進行分析,建立了綜合考慮這3種計算方法的水侵量最優化模型。也有一些學者提出了水驅氣藏儲量和水侵量簡易計算方法[20-22],基于氣井二項式產能方程、試井滲透率、氣水相滲曲線的水驅氣藏水侵量和儲量計算方法[23],基于水驅特征曲線的儲量計算方法[24]。但是,以上最優化計算方法通常需要利用計算機編程來實現,現場應用較為困難,而簡易方法大多需要數據擬合得出經驗相關式,理論基礎薄弱。
因此,筆者基于氣藏壓力下降引起水體壓力下降、水體孔隙體積收縮和水體中水的彈性膨脹導致水侵的機理,應用水體物質平衡原理,推導出了水侵量計算模型,并結合氣藏物質平衡方程,建立了水驅氣藏動態儲量和水侵量計算新方法,并進行了驗證和實際應用。該方法可同時計算出氣藏的動態儲量和水侵量,也可判斷產氣驅動能量占比的變化,為氣藏防水治水提供指導。
水侵量的計算是水驅氣藏動態儲量計算的先決條件,即在計算水驅氣藏動態儲量之前,須提前用其他方法計算水驅氣藏的水侵量。關于水侵量的計算,一直是水驅氣藏動態儲量計算的難點。本文將應用物質平衡原理,以水驅氣藏的水體為研究對象,將有限水體的彈性膨脹量與水侵量關聯起來(圖1),結合水驅氣藏的物質平衡方程,形成水體-氣藏耦合物質平衡方程。

圖1 水驅氣藏示意圖Fig .1 Schematic diagram of gas reservoir with water influx
水驅氣藏物質平衡方程的通式為
GBgi=(G-Gp)Bg+We-WpBw+ΔVp+ΔVw
(1)
式(1)中:G為水驅氣藏天然氣原始地質儲量,108m3;Bgi為原始地層壓力下天然氣體積系數,無因次;Gp為累計產氣量,108m3;Bg為當前平均地層壓力下天然氣體積系數,無因次;We為水侵量,108m3;Wp為累計產水量,108m3;Bw為水的體積系數,無因次;ΔVp和ΔVw表示孔隙體積的縮小量和氣藏中孔隙水的彈性膨脹量,108m3。此處體積單位均采用億方,便于數值計算。
ΔVp和ΔVw通常采用平均地層壓力降幅的簡化表達式,這種簡化在孔隙壓縮系數和水的等溫壓縮系數數值較小的情況下,可以滿足計算要求,但是當其數值較大,簡化表達式將造成較大誤差,本文將應用孔隙壓縮系數和水的等溫壓縮系數的定義式,嚴格推導ΔVp和ΔVw的精確表達式。
應用孔隙壓縮系數的定義式
(2)
式(2)中:Cp為氣藏孔隙壓縮系數,MPa-1;p為平均地層壓力,MPa;Vp為孔隙體積,108m3。
對式(2)從原始地層壓力到目前平均地層壓力進行積分得
Vp=VpieCp(p-pi)
(3)
壓力從原始地層壓力pi降為p,孔隙體積縮小量ΔVp可表示為
ΔVp=Vpi[1-eCp(p-pi)]
(4)
同理,應用水的等溫壓縮系數定義式

(5)
式(5)中:Cw為水的等溫壓縮系數,MPa-1;p為平均地層壓力,MPa;Vw為孔隙水體積,108m3。
對式(5)從原始地層壓力到目前平均地層壓力進行積分得
Vw=VwieCw(pi-p)
(6)
壓力從pi降為p,孔隙中原始含水的體積膨脹量ΔVw可表示為
ΔVw=Vwi[eCw(pi-p)-1]
(7)
假設水體孔隙體積為Ve,水體孔隙中充滿水,因此水體中水的體積也為Ve,水侵量We為水體壓力從pi降為p時擠出的水量。在氣藏流體產出過程中,氣藏和水體間為一非平衡系統,水體壓力高于氣藏壓力,即水體壓降小于氣藏壓降。假設水體壓降與氣藏壓降的比值為β,其數值與氣井的平均產氣量大小有關。平均產氣量越小,平衡條件越容易滿足,β越接近1;平均產氣量越大,非平衡程度越嚴重,則β越小(圖2)。β的取值大于0而小于等于1。β=1,表明水體和氣藏達到平衡,即水體平均地層壓力和氣藏平均地層壓力同步下降。

圖2 氣藏壓降、水體壓降與平均產氣量之間關系的示意圖Fig .2 Schematic diagram for the relationship among reservoir pressure drop,aquifer pressure drop,and average gas production rate
水體中孔隙體積收縮和水體膨脹量可參照式(4)和式(7)的推導過程。假設水體孔隙壓縮系數與氣藏孔隙壓縮系數相同,因為氣藏和水體一般處于同一儲層,因此該條件不難滿足。根據水體孔隙壓縮系數和水體中水的等溫壓縮系數公式,可得水體中孔隙體積收縮量為
ΔVpa=Ve[1-eCpβ(p-pi)]
(8)
水體中水的膨脹量為
ΔVwa=Ve[eCwβ(pi-p)-1]
(9)
那么水侵量為水體中孔隙體積收縮量加上水體中水的膨脹量,表達式為
We=Ve[eCwβ(pi-p)-eCpβ(p-pi)]
(10)
假設水體倍數為n,則Ve=nVpi,則
We=nVpi[eCwβ(pi-p)-eCpβ(p-pi)]
(11)
鑒于水侵初期水體壓降只波及到了部分水體,在水侵的過程中,水體壓降波及是個動態過程,這種動態變化是氣藏中天然氣逐漸被采出所導致的,氣藏原始狀態累計產氣量為0,水體壓降波及體積為0,水侵量為0;隨著累計產氣量的增大,水體波及體積增大,直到全部水體被波及,因此選用Gpb進行修正。選用冪函數的原因在于既要滿足水體波及體積隨著累計產氣量的增大而增大的特征,又要滿足累計產氣量為零時水體波及體積為0。因此在式(11)的基礎上加入對累計產氣量的修正,表達式為
We=GpbnVpi[eCwβ(pi-p)-eCpβ(p-pi)]
(12)
式(12)中:We為水侵量,108m3;Gp為累計產氣量,108m3;b為水體波及體積修正因子,無因次。若b=0,表明水體波及速度極快,即不考慮水體波及的動態過程。若b=1,表明水體波及體積和累計產氣量是正比關系。
將ΔVp的表達式(式(4))、ΔVw的表達式(式(7))和We的表達式(式(12))代入式(1)可得
GBgi=(G-Gp)Bg+GpbnVpi[eCwβ(pi-p)-eCpβ(p-pi)]-WpBw+Vpi[1-eCp(p-pi)]+Vwi[eCw(pi-p)-1]
由于
那么
(13)
式(13)中:Swi為氣藏原始含水飽和度。
整理得
GpBg+WpBw=G{(Bg-Bgi)+
(14)
若令
F=GpBg+WpBw
(15)
Eg=Bg-Bgi
(16)
(17)
(18)
E=Eg+Efw+Eaq
(19)
那么
F=GE
(20)
式(15)~(20)中:F為產出流體的地下體積,108m3;Eg為氣體彈性膨脹能,無因次;Efw為氣藏中孔隙體積收縮和氣藏中水的彈性膨脹能,無因次;Eaq為水體孔隙體積收縮和水體中水的彈性膨脹能,無因次;E為天然氣產出綜合彈性膨脹能,無因次。
應用多個實測壓力點,以F為縱軸,E為橫軸做圖,先設置b=0,試算水體倍數n和β的值,β值通常取1,主要改變n值,使所有數據都近乎落在同一條直線上,且截距盡可能接近0。若截距大于0,則增大n的取值,若截距小于0,則減小n的取值。最終得出最接近零截距的n值即為最終的n值,對應直線斜率值即為該水驅氣藏的動態儲量。
將式(20)變形為
(21)
應用多個實測壓力點,以F/E為縱軸,Gp為橫軸做圖,先設置b=0,試算水體倍數n和β的值,使所有數據都近乎落在同一水平線上,對應的縱軸數值即為該水驅氣藏的動態儲量。若直線向下偏離水平線,則增大n值,若直線向上偏離水平線,則減小n值,直到直線非常接近水平線,得出最終的n值。若在b=0的情況下,無論如何調整n和β,都無法滿足水平線,那么再調整b值。
以上2種確定方法可結合使用,在F與E的直角坐標圖中,使所有數據都近乎落在同一條直線上,且截距盡可能接近零;在F/E與Gp的直角坐標圖中,使所有數據都近乎落在同一水平線上;且F與E的直角坐標圖中直線點的斜率近乎等于F/E與Gp的直角坐標圖中水平線縱軸數值。那么該數值即為水驅氣藏動態儲量G。
根據最終確定出的b、n和β的值,可以計算出水侵量,公式如下:
(22)
計算出動態儲量和水侵量之后,就可以評價天然氣開采過程中各種驅動能量在整個驅動能量中占比的變化,將式(19)變形可得
(23)
對于非異常高壓的水驅氣藏,若Cp和Cw數值很小,約小于5×10-4MPa-1,則式(14)可以簡化為
GpBg+WpBw=G[(Bg-Bgi)+
(24)
令nβ=n′,則最終表達式可以簡化為
GpBg+WpBw=G[(Bg-Bgi)+
(25)
根據式(15)~(19),Efw+Eaq將簡化為
Efw+Eaq=
(26)
式(26)中:n′為等效水體倍數。通過這樣的簡化,將簡化前的3個調整參數b、n和β減少成了簡化后的2個調整的參數b和n′。
根據最終確定出的b和等效水體倍數n′的值,可以計算出水侵量,公式如下:
(27)
首先用理論數值模型對本文方法進行驗證,然后利用文獻中已有的實例進行驗證。
采用CMG軟件構建理論數值模型,水驅氣藏埋深3 000 m,原始地層壓力為30 MPa,儲層溫度為100 ℃,氣層厚度為32 m,孔隙度為0.2,滲透率為8 mD,原始含水飽和度為0.5,孔隙壓縮系數為0.000 32 MPa-1,水的等溫壓縮系數為0.000 445 9 MPa-1,水的體積系數為1.030 15,天然氣相對密度為0.58,氣藏外邊界半徑為400 m,水體模型選用Fetkovitch模型,水體外邊界與氣藏外邊界的比值分別設置為1.7和1.9,那么計算得出對應的水體倍數為1.89和2.61。由于CMG軟件中設定的地表溫度為15.6 ℃,因此計算得出該氣藏天然氣原始體積系數為0.0042 98,采用容積法計算得出該氣藏的天然氣原始地質儲量為3.742 2×108m3。應用建立的CMG理論數值模型,定產氣量5×104m3/d以一個月的時間步長生產了10年,輸出累計產氣量、累計產水量、平均地層壓力、平均視壓力(p/Z)和累計水侵量的動態數據。水體外邊界與氣藏外邊界的比值為1.7和1.9的2個案例,以年為時間步長的動態數據和水侵量計算表如表1和表2所示。

表1 水體外邊界與氣藏外邊界的比值為1.7時動態數據與水侵量計算表Table 1 Calculation table of dynamic data and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.7

表2 水體外邊界與氣藏外邊界的比值為1.9時動態數據與水侵量計算表Table 2 Calculation table of dynamic data and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.9
然后假設天然氣原始地質儲量、水體倍數和水侵量未知,已知原始含水飽和度、孔隙壓縮系數、水的等溫壓縮系數和水的體積系數,并采用CMG數值模擬模型輸出的累計產氣量、累計產水量和平均視壓力(p/Z)數據,應用本文建立的方法,反求該氣藏的天然氣原始地質儲量和水侵量動態。若數據點滿足本文提出的直線關系和水平線關系,且求出的動態儲量與模型真實儲量(容積法計算得出)相差不大,水侵量計算結果與理論數值模型輸出結果接近,則證明本文提出的水驅氣藏動態儲量和水侵量半解析方法是合理和可靠的。
理論數值模型驗證過程中,b=0且β=1就能很好地擬合出直線和水平線,因此對于理論數值模型,可以認為水侵初期就波及到了整個水體,水體和氣藏實時達到平衡狀態。水體邊界與氣藏外邊界的比值為1.7和1.9的擬合圖如圖3、4所示。可以看出這2個例子F與E的直線關系非常好,且幾乎穿過原點(圖3a、4a);F/E與Gp的水平線完全水平(圖3b、4b)。調整確定的水體倍數n分別為1.858 6和2.605 2,計算的天然氣動態儲量分別為3.725×108m3和3.747×108m3,與理論數值模型真實天然氣地質儲量3.742 2×108m3相差不大。如表1、表2、圖3c和圖4c所示,新方法計算的水侵量與理論數值模型輸出的真實水侵量幾乎完全吻合。理論數值模型驗證結果表明,本文建立的水驅氣藏動態儲量和水侵量計算新方法在理論上是正確的。

圖3 水體外邊界與氣藏外邊界比值為1.7時動態儲量與水侵量計算結果Fig .3 Calculation results of reserves and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.7

圖4 水體外邊界與氣藏外邊界比值為1.9時動態儲量與水侵量計算結果Fig .4 Calculation results of reserves and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.9
實例的數據選自文獻[20],氣藏原始地層壓力為78.789 MPa,儲層溫度為111.85 ℃,綜合壓縮系數為0.000 542 69 MPa-1,原始含水飽和度為0.562 2[25],假設水的體積系數為1.071 36,水的等溫壓縮系數為0.000 235 MPa-1,那么計算出孔隙壓縮系數為0.000 105 5 MPa-1,累計產氣量、累計產水量、壓力、視壓力和實際水侵量等動態數據[20]見表3所示。應用本文方法,計算出Eg、Efw和Eaq、E和F的值(表3)。

表3 某實例水驅氣藏動態數據[20]與水侵量計算表Table 3 Dynamic data and water influx calculation table of a water drive gas reservoir
應用本文方法,F與E擬合圖、F/E與Gp擬合圖和水侵量計算結果如圖5所示。可以看出,F與E的直線關系非常好,且幾乎穿過原點(圖5a);F/E與Gp的水平線近乎水平(圖5b)。調整確定的水體倍數n為5.03,水體壓降與氣藏壓降的比值β為1,b值為1.04,求得該氣藏天然氣動態儲量為19.232×108m3,與文獻中該氣藏實際天然氣地質儲量19.321×108m3相差不大。如表3和圖5c所示,新方法計算的水侵量與文獻中該氣藏的實際水侵量吻合程度也很高。該實例驗證結果表明,本文建立的水驅氣藏動態儲量和水侵量計算新方法在實際應用中也是可靠的,可以推廣應用。

圖5 文獻中氣藏實例動態儲量與水侵量計算結果Fig .5 Calculation results of reserves and water influx of the gas reservoir in the literature
番禺某氣田分布著多個中高孔、中特高滲邊底水氣藏,氣藏A被A1井控制,氣藏B被B1井控制,氣藏C被C1井控制。這些井所控制的氣藏的物性參數見表4。單井的累計產氣量、累計產水量、平均地層壓力和平均視壓力數據見表5。
應用本文方法計算出各井的Eg、Efw和Eaq、E和F的值(表5)。F與E擬合圖、F/E與Gp擬合圖和水侵量計算結果如圖6所示。可以看出,3口井F與E的直線關系都非常好,且都幾乎穿過原點(圖6a);3口井F/E與Gp的曲線接近水平線(圖6b)。3口井調整確定的水體倍數n、水體壓降與氣藏壓降的比值β和b值,以及求得的各井控制的動態儲量見表4。本文方法計算的各井水侵量見表5和圖6c。A1、B1和C1井產氣驅動能量占比變化曲線如圖6d、e和f所示。可以看出A1井和B1井天然氣膨脹能和水體能量相當,而C1井水體能量是天然氣生產的主要驅動能量,一直維持在90%左右。此外,A1井和B1井水體能量在開發過程中是變化的,A1井從最初的45%上升為53%,水體能量越來越強;B1井水體能量從最初的48%降低為42%,水體驅動能量逐漸變弱。

表4 番禺某氣田各井物性參數表Table 4 Physical property parameters of three wells in Panyu gas field

表5 番禺某氣田3口井動態Table 5 Production performance of three wells in Panyu gas field

圖6 番禺某氣田3口井動態儲量與水侵量計算結果Fig .6 Calculation results of reserves and water influx of three wells in Panyu gas field
1) 本文建立的水驅氣藏動態儲量和水侵量計算新方法,可同時計算出氣藏的動態儲量和水侵量,適用于同一壓力系統的氣藏、區塊或單井,除原始地層壓力外,還需要至少2個靜壓值才能計算,數據點越多,且采出程度越高,儲量和水侵量計算結果越可靠。
2) 綜合產出流體地下體積F與綜合彈性能E指示曲線和F/E與Gp指示曲線,增加了擬合的精度和解釋結果的唯一性,理論數值模型和實例驗證均表明:F與E指示曲線確實符合穿過原點的直線,斜率即為氣藏的動態儲量,F/E與Gp指示曲線的確為一水平線,縱軸數值即為氣藏的動態儲量,且水侵量計算結果與理論數值模型和實際水侵量吻合程度很高,從理論和實踐均證明所建立的水驅氣藏動態儲量與水侵量計算方法是正確的。
3) 建立的水驅氣藏動態儲量和水侵量計算方法也可判斷天然氣彈性能、氣藏巖石和孔隙水彈性膨脹能、水體能量在整個驅動能量中的占比的變化,評價水驅氣藏水侵能量強弱,為氣藏防水治水措施提供依據。
4) 對于異常高壓水驅氣藏,水體中的溶解氣的彈性膨脹能量不可忽略,考慮水體中溶解氣的影響將是異常高壓水驅氣藏動態儲量和水侵量計算的下步研究方向。