1 INTRODUCTION

The calculation of oil saturation in hydrocarbon reservoirs and the judgment of flooded level in flooded zones depend on the accurate acquisition of connate water resistivity which has always troubled the petrophysicists and log analysts to a great extent^[1]. Spontaneous Potential (SP) logging is an important method which can provide the information of formation permeability as well as initial water resistivity in place. There are four resources contributing to the SP,including the diffusion potential,the membrane or Nernst potential,the electro-kinetic or streaming potential,and the oxidation/reduction potential,of which the latter two are relatively weak and can be ignored. Thus the SP in permeable formation is primarily the sum of the membrane and diffusion potentials which are mainly determined by the ratio of the mud filtrate to formation connate water salinity^[2]. Since the membrane and diffusion potential are mainly generated in the interface of different fluids,it will be helpful to analyze and identify the radial distribution of the fluid salinity during the process of mud filtrate invasion.

Nowadays,array induction logging tools capable of providing signals with different probing depths and same resolution,are widely used in China’s oilfields,including array Induction Imaging logging tool (AIT) produced by Schlumberger,high resolution induction logging tool (HDIL) produced by Baker Hughes and high resolution array induction tool (HRAI) developed by Halliburton^{[3, 4, 5, 6, 7]}. Based on the abundant information provided by array induction logging data,the proper formation model and reasonable algorithm,the formation resistivity distribution profile can be well reconstructed,meanwhile the generated position of the SP can also be determined. Thus,combined with SP logging data,formation water resistivity could be inverted effectively. In this work,through numerical simulation of mud invasion and analysis of SP generation position,a new “five-parameter” formation model is proposed. Then based on the joint inversion of array induction logging and SP,reconstruction of formation resistivity distribution profile and accurate calculation of the formation water resistivity are accomplished.

2 NUMERICAL SIMULATION OF MUD FILTRATE INVASION AND DETERMINATION

OF GENERATED POSITION OF SP

2.1 Numerical Simulation of Mud Filtrate Invasion

The process of mud filtrate invasion can be divided into two individual processes: the displacement of formation fluids by mud filtrate which results in the change of saturation in formation radial direction,and the diffusion of mud filtrate and formation water icon which contributes to the change of salinity of formation^{[8, 9]}. Assume that mud filtrate and formation fluids in the displacement process are not miscible; reservoir rock and fluids are both slightly compressible; fluid transfusion complies with the law of Darcy. Then according to the B-L’s water oil displacement theory,the water-oil flow equation can be depicted as follows:

where,P_w and P_o are water phase pressure and oil phase pressure,respectively; S_oand S_ware oil saturation and water saturation,respectively; μ_o and μ_w are oil viscosity and water viscosity,respectively; K_ro and K_rw are the relative permeability of oil and water,respectively; φ is formation porosity; h is formation thickness; g is gravity acceleration; and ρ_w and ρ_o are water density and oil density,respectively.

Then the diffusion equation can be written as

where C_wis the formation water salinity or mud filtrate salinity. Then consider the Archie’s equation where a,b,m and n are Archie parameters,and R_wis the resistivity of formation or mud filtrate.

Numerical simulation of reservoir saturation,salinity,resistivity distribution during mud invasion process can be achieved with an appropriate formation model and finite difference method (FDM).

2.2 Determination of Generated Position of SP

In order to study the distribution characteristics of fluid and resistivity around the borehole during the mud invasion,typical oil layer,water layer and oil-water layer models were established. Reservoir saturation,salinity and resistivity profiles of different models were simulated (Fig.1). Model parameters are: formation porosity is 10%,permeability is 100 μm²,water salinity is 80000 mg/L,mud filtrate salinity is 44000 mg/L,residual hydrocarbon saturation is 20%,irreducible water saturation is 20%; Oil saturation: 75% in the oil layer,25% in water layer,and 50% in the oil-water layer.

The following conclusions can be obtained from the distribution patterns of reservoir saturation,salinity,resistivity with different fluid types:

① The saturation and salinity vary monotonically in the radial direction during the mud filtrate invasion. The non-monotonic variation of radial resistivity is due to the speed difference of two effects;

②Saturation profile change is largely affected by original water saturation,as original water saturation increases,the displacement depth reduces gradually;

③Salinity change is less affected by original water saturation. Compared with the slowly changed saturation profiles,steep stepwise change appeared in salinity profiles;

④The SP is generated at the ladder of salinity mutation. In the resistivity profile,for the oil layer the position is at the rising part of the ladder; for the water layer lay is at the declining part of the ladder; and for the oil-water layer the position is at the declining part of “low resistivity annulus”.

3 “FIVE-PARAMETER” FORMATION MODEL

According to the resistivity distribution characteristic of mud filtrate invasion^[10],the approximate resistivity distribution model around borehole is established as shown in Fig.2.

The model is described by

where R_xo = R₁is the resistivity of flushed zone,r₁and r₂(r₂ > r₁) are the radius of flushed zone and transition zone,R_t= R₂is the original formation resistivity,R_iis the resistivity of transition zone which is assumed to present a parabola distribution,that is,R_i= ar² +br+c,where r is the radial depth,a,b and c are coefficients and can be determined by (r₁,R₁),(r₂,R₂) and ((r₁ + r₂)/2,R_1/2),where R_1/2 represents the resistivity of the middle point of the transition zone. Formation parameters (r₁,r₂,R₁,R₂,R_1/2) reflect the property of mud filtrate invasion,and determine the response of the array induction logs in mud-invaded formation.

When the parabola coefficient a = 0,the new model degenerates into a “4-parameter model” in which the transition zone resistivity is assumed as a slope,and when r₁≈ r₂,the model degenerates into the “3-parameter model” in which the formation resis- tivity is described as a piston. Thus various patterns of resistivity distribution including the “low resistivity annulus” can be described by adjusting the parameters.

4 JOINT INVERSION OF ARRAY INDUCTION AND SPONTANEOUS POTENTIAL LOG

GING During the calculation process of formation water resistivity,the joint inversion method in combination with physical parameters is adopted based on the data of array induction and SP logging. First,the resistivity profile is reconstructed and the generation position of SP is determined based on the “five parameters” formation model. Then,take the formation resistivity and the generation position of SP as known parameters,the formation water resistivity is inverted by using SP logging data in conjunction with mud filtrate resistivity. Generally,the inversion starts from non-permeable formation. By using the method of depth domain windowing,multi-layers above and under the target layer are inverted at the same time in order to eliminate the influence of surrounding rocks. The specific process is as follows.

Formation parameters (r₁,r₂,R₁,R₂,R_1/2) determine the response of array induction logs in mud-invaded formation,which can be expressed by the following mathematical function

where x is formation parameter vector and y_i is the response of array induction log or SP of the ith point. By using the optimized least square method,the problem can be written as

To minimize the objective function,the following equation is required

Eq.(8) is transformed into the following equation based on the Taylor formula

Solve the problem using regularization method with a regularization factor

where λ is a regularization factor that can be determined by numerical experiments,and in our case is equal to 2.0.

4.1 Array Induction Inversion Algorithm

Take the High Resolution Induction Logging tool (HDIL) from Baker Hughes for example,which can provide resistivity curves of three resolutions (M1-0.3 m,M2-0.6 m,M4-1.2 m) and six probing depths (R1-0.25 m,R2-0.5 m,R3-0.75 m,R6-1.5 m,R9-2.3 m,RX-3.0 m). Considering the signal processing of array induction signals and the overfocusing phenomenon of HDIL^[11],the curve M2R1 (the resolution is 0.6 m,probing depth is 1.2 m) is abandoned,meanwhile the curves M1R2 (the resolution is 0.3 m,the probing depth is 0.5 m),and M4RX (the resolution is 1.2 m,the probing depth is 3.0 m),which are less processed,are adopted for inversion. Assuming that there are M layers in the depth window,inversion is based on the array induction logging response of N points N(N ≥ M),for which Jacobian matrix is expressed as

where y_n,i denotes the ith induction log of nth depth point,x_m,1,x_m,2,· · · ,x_m,5 are the five parameters (r₁^m ,r₂^m ,R^m ₁ ,R₂^m ,R_1/2^m) of mth layer. During the inversion process,the calculation of Jacobian matrix occupies most of the time. In order to speed up the inversion speed,in this paper the followed method is proposed to calculate inversion matrix A. Assuming that the transverse conductivity of the medium is piecewise continuous,the induction logging apparent conductivity can be expressed as where σ₀(r),σ_M(r),σ_m(r) are the conductivities of each layer and can be determined by the following formula where H(z) is the Heaviside function,r_oh is the borehole radius,r^* _l (l = 2,3,· · · ,L) is the discrete point between r₁and r₂,r^* ₁ = r₁,r^* ₂ = r₂.

Consider

Thus there is

where g(r,z) is 2D geometric factor of HDIL which can be used directly in the process of inversion because the skin effect has been eliminated. This approach is highly computation efficient because the geometric factor is previously known,avoiding repeated calculation of the geometric factor.

4.2 Spontaneous Potential Logging Inversion Algorithm

During the drilling process,due to the differences of drilling fluid and formation water salinity as well as the formation lithology,dipole layer generates on the surface of different salinity solution. From the principle of the SP logging,a dipole source of SP is distributed on the borehole and the boundaries between different formations (Fig.4).

The numerical simulation of SP responses is adopted to calculate the Jacobian matrix in this paper. Assuming that the dipoles are distributed on the interface uniformly,the current potential can be used to implement the simulation of SP^[12],and the electric field equation is

where V is the current potential; "ssp is the static spontaneous potential. For M layers,the spontaneous potential response of the nth depth point is the response superposition of M static spontaneous potential sources,which could be expressed as where S_n^P m is the spontaneous potential response of nth point generated by the mth static spontaneous potential source.

Define the functional

Eq.(21) could be solved quickly by using the finite element method to obtain SPⁿ.

4.3 Formation Water Resistivity (Rw) Acquisition

For pure sandstone formation,the formation water is mainly NaCl solution,the maximum formation spontaneous potential is determined by the ratio of the mud filtrate to formation water resistivity,which can be expressed as^[13]

where T_F is Fahrenheit temperature,R_mfe is the equivalent mud filtrate resistivity,and R_we is the equivalent formation water resistivity which can be transformed into formation water resistivity through the equation^{[14, 15]}

5 FIELD EXAMPLE

The inversion results of two sections from an oilfield in East China were selected and analyzed. The depth is 3880~3900 m (Fig.5),borehole diameter is 0.2 m,mud resistivity is 0.154 m. First,the activity layering method was adopted to implement the automatic zonation and characteristic value acquisition of the formation (in the third panel). The radial resistivity profile can be reconstructed based on the “five-parameters” model,which can reflect the radial variation of formation resistivity intuitively. The resistivity profiles of No.3 layer and 3897.38~3898.56 m of No.4 layer are shown in Fig.6,in which the SP generated positions are represented by the red lines. The inverted static spontaneous potentials (ε_ssp) of No.2 layer and No.4 layer are both -44 mV indicating both the two layers have good property with the same connate water resistivity. The potential difference between No.2 layer (-34.1 mV) and No.4 layer (-43.3 mV) is due to the resistivity difference of the two layers. The inverted connate water salinity (64282 mg/L) agrees well with chemical analysis of formation water (64443.13 mg/L). The reconstructed SP curve and HDIL curves are displayed in the first and fifth panels,which show good agreement with measured curves which proves the inversion algorithm is stable and reliable.

Meanwhile,the reconstruction of formation radial resistivity through inversion aided to the identification of formation fluid types. A_sshown in Fig.5,though No.3 layer shows similar apparent resistivity as No.2 layer,its inverted resistivity profile exhibites a “low resistivity annulus” (R_1/2 < R_t < R_xo) and its original formation resistivity is significantly lower than No.2 layer which indicates a typical oil-water characteristic,therefore is interpreted as an oil-water layer. The inversion results of No.2 layer present a typical oil layer with conductive invasion,while No.4 layer shows a typical resistive invasion characteristic which indicates a water layer. From the oil testing conclusions,No.2 layer is demonstrated to be an oil-bearing layer while No.4 layer is demonstrated to be a water layer,confirming the inversion results.

To further verify the correctness and applicability of the algorithm,18 sandstone layers in an oilfield with pure lithology were inverted with the proposed method. Part of the inversion results of water salinity and the chemical analysis results is compared in Table1. Connate water salinities in this area were obtained based on the joint inversion of HDIL and SP curves,and the results all show good agreement with the chemical analysis results. Inverted results show that the connate water salinity is 85000~120130212 mg/L in this area.

The joint inversion of HDIL and SP is a reliable method to estimate connate water resistivity,which ensures the accurate calculation of reservoir saturation. It provides a feasible method of formation evaluation,especially for the (natural) water flooded layers and zones lack of chemical analysis of formation water.

6 CONCLUSIONS

The non-monotonic variation of radial resistivity is due to the speed difference of displacement and diffusion effects during mud filtrate invasion. The SP is generated at the rising part of the ladder for the oil layer,the declining part of the ladder for the water layer lay,and the declining part of “low resistivity annulus” for the oil-water layer.

A variety of patterns of resistivity distribution including the “low resistivity annulus” can be described by the “five-parameter” model. The joint inversion of HDIL and SP can reconstruct the radial resistivity,identify the SP generated position,and calculate the connate water resistivity precisely,which can provide evidence for reservoir evaluation and fluid property identification.

The proposed joint inversion algorithm,which is stable and reliable,can provide a feasible method of formation evaluation for the (natural) water flooded layers and zones lack of chemical analysis of formation water. However,additional analysis is needed for the cases of the formation with high shale content or water with ion components other than those of NaCl. In addition,when the invasion depth is too large and exceeds the investigation depth of the array induction logs,the inversion results are open to doubt.

7 ACKNOWLEDGMENTS

The research was supported by the Project of National Natural Science Foundation of China (41174099),CNPC National Major Projects (2011ZX05020-008),and CNPC Science and Technology Development Projects (2011D 4101).