Finite-element modeling of spontaneous potential in an axisymmetric reservoir model with account of its shale content

: (1) The article is devoted to the development of a theoretical and algorithmic basis for numerical modeling of the spontaneous potential method (SP) as applied to the study of sandy-argillaceous reservoirs. (2) In terms of coupled flows, we consider a physical-mathematical model of SP signals from an electrochemical source, with regard to the case of fluid-saturated shaly sandstone. (3) An algorithm for 2D finite-element modeling of SP signals was developed and implemented in software, along with its internal and external testing with analytical solutions. The numerical SP modeling was carried out, with determining the dependences on the reservoir thickness and porosity, the amount of argillaceous material and the type of minerals. We performed a comparative analysis of the simulated and field SP data, using the results of laboratory core examina-tions taken from wells in a number of fields in the Latitudinal Ob Region of Western Siberia. (4) The results of the study may be used either for the development of the existing SP techniques, by provid-ing them with a consistent computational model, or for the design of new experimental approaches.


Introduction
Over the last decade, there has been a significant technological development of logging equipment for geophysical surveys in oil and gas wells. High-precision devices for synchronous wireline logging in one trip in vertical wells were created, as well as autonomous complexes for logging on drill pipes in highly deviated wells. Among a large number of geophysical methods implemented in modern equipment, one of the most demanded in the study of geological sections and widely used in all drilled wells is the method of spontaneous polarization potentials (SP).
SP logging consists in the recording along the wellbore of electrical signals arising naturally due to physical and chemical processes occurring in the near-wellbore space [1][2][3]. SP sources include electrokinetic [4][5][6], electrochemical [7,8] and redox processes [9,10]. Despite the long history of the SP method, it still remains the key method for studying geological sections exposed by oil and gas wells in order to identify reservoirs, quantify their filtration-capacitative properties, namely porosity, sandiness, salinity of formation water, etc. [11][12][13][14].
With all the development of computational technology, correction charts remain the main tool for analyzing SP logs. These are calculated using equivalent electrical circuits and analytical solutions in simplified formulations [1,15,16]. Correction charts allows one to determine the value of static spontaneous potential (SSP) as an electromotive force, under the action of which electric currents circulating in and near the borehole arise. It is assumed that the region of the borehole adjacent to permeable formation is electrically isolated from the rest of the borehole. As a result, SSP depends only on the ratio of salinity of the drilling mud filtrate and of the adjacent formations (and, analogously, the ratio of the resistivities of the drilling mud filtrate and formation water). However, this closedcircuit approach does not fully reflect the physical and chemical processes leading to SP formation [17,18]. To fully describe the physical and chemical processes occurring in the near-wellbore space, it is preferable to use the theory of coupled fluxes and nonequilibrium thermodynamics [17,19] and apply numerical methods to solve differential equations describing the formation of SP signals [20][21][22]. In this connection, one of the most important tasks today is the application of modern methods of computational mathematics and the development of efficient computational algorithms and computer programs allowing the simulation of spontaneous polarization potentials in a class of realistic models that adequately reflect the experimental conditions. To date, a number of approaches for the simulation of SP logging data based on the finite element method (FEM) have been proposed [23][24][25].
We propose an algorithm for numerical simulation of SP logging data in an axisymmetric model of a fluid-saturated sandstone reservoir in-between shale sediments. A physical-mathematical model of a clayey sandstone reservoir is used, which takes into account both the clay material content and the type of clay minerals. Numerical solution of the forward SP logging problem is performed in a two-dimensional formulation and is based on the finite element method with approximation of electrophysical parameters by sigmoidal functions. The developed algorithm is tested with successful reproduction of analytical results obtained in approximation of infinitely extended dipole layers [2]. The computational algorithm is used to numerically simulate SP signals with the establishment of their dependencies on model parameters. The capabilities of the algorithm are demonstrated on field SP logging data using the results of laboratory studies of core from intervals of the Lower Cretaceous deposits, penetrated by wells in a number of Western Siberian fields.

The description of the physical model
The SP problem is considered in a two-dimensional geological model, including a permeable sandstone bed and adjacent impermeable shale deposits (Fig. 1). A constant element of the model is a well, the boundary of which is described by a cylindrical surface. When the permeable formation is penetrated, fresh water-based mud filtrate penetrates the medium, forming a zone with changed fluid saturation outside the well -the zone of drilling mud filtrate invasion [26,27]. The correct description of the SP phenomenon from the basic physical principles is carried out with the help of the coupled fluxes formalism. The coupled fluxes are understood as a combination of conduction current, caused by the presence of an electric field in the medium, and charged particle current caused by physical and chemical processes of different nature occurring in the near-wellbore area [18,19,25]. The density of the total volumetric current flowing in the medium is represented as a sum of two components: where is volumetric conduction current and is volumetric current of SP source. These source currents are caused by the presence of ion activity gradients, filtration of drilling mud and formation water, as well as redox processes. In this case, the divergence of the total volumetric current density is assumed to be equal to zero, which is equivalent to the requirement of conservation of the total charge in the medium: Expressing the conduction current density through the electric field in the medium ( = ) and introducing its scalar potential ( = −∇ ), we obtain the Poisson equation of the SP problem: where is electrical conductivity of the medium, is scalar potential. The resulting equation defines the electrostatic potential distribution for an arbitrary source of SP that induces an external current .
In this paper, the contribution of ion activity gradients, also called the electrochemical current source, is considered as such SP source. Since the other two contributions to the formation of SP, namely, caused by drilling mud/ formation water filtration and redox reactions, are much smaller than the contribution of ion diffusion, in most cases they can be neglected. The expression for the current density of an electrochemical source can be derived from the basic principles of ionic diffusion in charged porous media. Here we utilize the approach presented in [2], according to which the current density of an electrochemical source can be determined as: where and are diffusion and diffusion-adsorption EMF coefficients, respectively, is diffusion-adsorption activity of rock, is electrical conductivity of the fluid, depending on the concentration of ions at a given temperature. Given the definition of the source current density (4), the Poisson equation (3) takes the following form: where = and = are the corresponding values of electrical resistivity.

Accounting for shale content
It is well established that in the case of pure sandstone with a non-conductive matrix, the resistivity of the formation and the resistivity of the fluid saturating it are related by the Dakhnov-Archie relation [28,29]. However, this empirical model does not take into account a number of factors, including shale material, the content and type of which significantly affects the overall resistivity of the medium. Partial consideration of shalyness using the description of the problem (5) is possible only with the involvement of data from laboratory core studies, which will be shown below.
To correctly account for the effect of clay material content and type on the SP signal in the case of a clayey reservoir (Fig. 1), it is necessary to use the parameterization, in which the values of and in equation (5) are redefined as follows. Thus, the diffusion-adsorption EMF coefficient can be represented in the following form: where ( ) is macroscopic Hittorff number (transfer number for cations), is Boltzmann constant, is temperature, is elementary charge. The macroscopic Hittorff number for cations reflects the fraction of electric current in the porous medium carried as a result of cation transport, and can be expressed through the petrophysical properties of the formation as follows [17,25]: where is porosity, is water saturation, is effective mobility of cations on the mineral surface, is solid phase density, is cation exchange capacity, ( ) is microscopic Hittorff number, determined by the relative mobility of cations in the free electrolyte solution. Using these values, we can also obtain an expression for the reservoir resistivity in the following form: where is cementation exponent and is saturation exponent. This expression takes into account the contribution to conductivity from both formation water resistivity (Dakhnov-Archie law) and the effect of surface charges of clay minerals present in the formation (Waxman-Smits model) [30]. The magnitude of the latter is determined mainly by the cation exchange capacity. If the sandy reservoir contains clay minerals, it is necessary to use the average cation exchange capacity [17]: where is clayeness, is cation exchange capacity of sand and is cation exchange capacity of shale.

Finite element method computation
The model of the geological medium shown in Fig. 1 is axisymmetric, considered in the cylindrical coordinate system ( , ). In general, the model may include a bundle of formations with plane-parallel horizontal boundaries, with several radial zones in the near-wellbore area of each permeable formation separated from the well and from each other by coaxial-cylindrical boundaries (Fig. 2).
In equation (5), for each value of , the diffusion-adsorption EMF coefficient is a constant value, while and are functions of and take different values in the well and formation. The values of the functions and may be different in each of the radial zones of the near-wellbore area, while in the borehole they are taken the same and equal to the drilling fluid resistivity . The solution of the forward SP problem is based on the calculation of the electric potential distribution in the investigated area using equation (5) with a given distribution of functions , and parameter . A numerical algorithm based on FEM is proposed for modeling SP logging data. The use of FEM makes it possible to effectively calculate the electric potential distribution in logging problems both in the case of axial symmetry and three-dimensional geometry [31]. Here the FEM is used to solve the partial differential equation (5) within an axially symmetric region. The boundary conditions and the right-hand side of the equation are set and, importantly for the calculations, the parameters of the grid of the computational domain are determined and selected.
The electrode used to record signals during SP logging can be considered a sphere of small radius. The size of the electrode is much smaller than the borehole radius and its influence on the values of recorded signals is insignificant. The first boundary condition is that the normal component of the electric field on the electrode surface is zero, which is equivalent to the vanishing of the potential derivative along the normal: where B 1 is electrode surface. Also, a necessary boundary condition is that the potential at the distance from the current sources is equal to zero: where B 0 is the outer boundary of the modelling area. The choice of the size of the outer boundary is dictated by the requirement of zero potential and at the same time the need to provide an optimal calculation time.
Consider a two-dimensional region Ω, in the general case with inhomogeneous distribution of physical properties, on which we define the following functional spaces: where L (Ω) is a Lebesgue space for the elements of which scalar product can be introduced: The equation (5) in the general form is, as follows: In order to use FEM to solve equation (14) with boundary conditions (10) and (11), we turn to a variational formulation of the boundary problem. This requires finding a function φ ∈ H (Ω) for which any ν ∈ H (Ω) holds: This formulation of the problem automatically leads to the fulfillment of boundary conditions (10) and (11). To find an approximate solution φ of the variational problem, we will use the basis function expansion, the set of which forms a finite-dimensional subspace H (Ω) of the space H (Ω): where ψ is i-th basis function and is its weight in decomposition. For an approximate solution, the variational problem becomes discrete and takes the following form: At the same time, since any function ν ∈ H (Ω) can be represented as a basis function expansion of a given subspace, it is sufficient to require (17) only for all given basis functions, which leads to a set of N equalities of the following form: This expression can be rewritten in matrix form: where = ( , , … , ) is a vector of weights of basis functions uniquely determining the required approximate solution. Thus, the problem of finding the electric potential is reduced to the search for the solution of the system of linear equations (19), the coefficients of which are defined by the following expressions: Since the SP problem is two-dimensional in the axial symmetry approximation, the bilinear functions defined on a rectangular grid were chosen as basis functions. To find the solution of (19) the Kholetsky decomposition was used. This is due to the fact that the resulting matrix is positively definite, and can be easily reduced to the symmetric form.

Determining source current density
To determine the function of the right-hand side of equation (14), it is necessary to specify the distribution of current density of external sources, or, similarly, the resistivity, as well as the parameter as a function of spatial coordinates ( , ). Since the functions and take different values inside the well and inside each of the radial zones of the reservoir (Fig. 2), they can be represented as piecewise constant on and stepwise on functions of the following form: where ( ) is a radial function of the resistivity for i-th formation, depicted in fig.  2, is a resistivity of the fluid, , is a resistivity of j-th radial zone of i-th formation, ( ) is Heaviside function, , is a radius of j-th radial zone of i-th formation, -number of radioal zones in i-th formation, is a number of formations. A similar expression can be written for . Since the analytical calculation of the right-hand side of equation (14) requires finding the gradient of as well as finding the gradient and divergence of sequentially, the corresponding functions must be sufficiently smooth. Since the Heaviside function is not smooth, this approximation is not applicable. In this connection, the approximation of the radial resistivity functions with sigmoidal functions of the following form is used: where parameter determines the slope of the sigmoid, or the size of the transition region between adjacent zones. This transition region is a consequence of approximation of the initial piecewise constant function and does not reflect real transients between the well and the formation. In the limit -> ∞, equation (22) turns into (21). An example of the resulting form of the function of the right part for the i-th formation when setting the radial functions and by the method described above is shown in Fig. 3. To calculate the elements of the system of linear equations, it is necessary to calculate integrals from the obtained function on the finite element grid of the form (20). Since the function changes rapidly near the wellbore boundary, the Simpson method with an adaptive integration step is chosen for numerical integration. In this approach, the number of partitions of the integration interval is doubled until the difference of integral values calculated at the current and the previous step becomes less than the specified value. In this case, the function values calculated at the current step are used to calculate the following iterations, which reduces the computational cost. Despite the fact that Simpson's method is usually less efficient than some other methods of numerical integration (e.g., Gauss method), in this case its usage allows to get rid of the procedure of integration grid setting. This simplifies the software implementation of the algorithm, in particular, when varying the mesh of finite elements or values of simulation parameters.

Testing the computational algorithm
To test the developed algorithm, simulation of SP logging data was performed on several typical examples. A particular case of the geoelectric model shown in Fig. 2 includes a well, a formation of clean fluid-saturated sandstone without an invasion zone and the adjacent layers of shale sediments. In the well, the function takes the mud resistivity , while in the permeable reservoir it takes the sandstone reservoir resistivity . The value of the function in the well, as well as for the function , is equal to , and in the reservoir -to the resistivity of formation water . The diffusion-adsorption EMF coefficient in the limiting case of pure sandstone is -11.6 mV, while in shale sediments it takes the value of 58 mV. The simulation parameters used in testing the computational algorithm on all examples are shown in Table 1. By selecting a sandstone formation which thickness is many times greater than the borehole radius, it is possible to obtain a SP logging curve being close to the static SP curve. To match the model case of static SP, the resistivities of the adjacent formations and reservoir are assumed to be equal to the resistivity of the drilling mud. Calculated at different values of SP anomaly amplitude Δ in this approximation agree with the equation describing static SP (Fig. 4): where is the static SP coefficient, in this case equal to 69.6 mV. The influence of geometrical parameters on the magnitude and shape of the SP logging curve can be evaluated within the framework of the considered model. Simulation results were obtained for formations of different thickness and compared with the analytical formula given in [2] for the approximation of infinitely extended dipole layers (Fig. 5). The effect of finite reservoir thickness manifests itself in a decrease in the observed amplitude of the SP anomaly Δ relative to the static spontaneous potential . Figure 5. SP logging curves normalized to the static SP value from the thickness of the permeable formation. The cyphers of the corresponding curves denote the ratio of reservoir thickness to borehole radius. Circles are the results of calculation with the proposed numerical algorithm, the solid line is the calculation according to the analytical formula presented in [2].
In general case, the ratio of the resistivities of both the permeable formation and the adjacent formations to the resistivity of the drilling mud can significantly affect the shape and amplitude of SP signals. Fig. 6 shows the dependence of SP anomaly amplitude, normalized to the value of static SP, on the thickness of the sandstone reservoir for different ratios of reservoir and drilling mud resistivity. In this case, the resistivity of the adjacent formations is assumed to be equal to the resistivity of the reservoir. As can be seen from Fig. 6, the greater the sandstone reservoir resistivity exceeds the mud resistivity, the greater the difference between the SP anomaly amplitude and the static SP at a given value of ℎ / . The result is consistent with the correction charts for SP logging data [32].

Numerical modelling
The proposed algorithm is used to perform numerical simulations in typical model cases. First of all, we simulated a set of alternating sandstone beds of different thicknesses interlaid with shale sediments, for various ratios of the resistivity of the formation and the drilling fluid (Fig. 7). The model is represented by sandstone layers of different thickness (ℎ = 40 ; ℎ = 10 ; ℎ = 4 ), separated by clay interlayers of equal thickness (ℎ = 2 ), and includes a sandstone layer of 10 , containing in the center a thin clay interlayer (ℎ = 2 ). Modeling was performed for three cases with the ratio of the sandstone reservoir resistivity to the mud resistivity equal to 1, 5 and 20. In the intervals of thick reservoirs, the amplitudes of SP anomalies are close in value at different resistivity ratios. In this case, with decreasing reservoir thickness, a pronounced decrease in the amplitude, accompanied by smoothing of the SP signal diagram in the areas of intersection of reservoir boundaries with increasing permeability ratio of the reservoir and drilling fluid is observed. The amplitude of the SP anomaly of a sandstone reservoir containing a thin clay interlayer is less than the amplitude of the anomaly of a homogeneous reservoir of the same thickness. At the same time, the amplitude of the SP anomaly in the interbedded layer is greater than the amplitude of the anomalies in the constituent sandstone beds, if they are separated by a thick clay interlayer. The above numerical experiments were performed without taking into account the clay content of the reservoir formation. In accordance with section 2.2, simulation of SP logging data in the interval of reservoir formation characterized by both changes in clay content and its type (Fig. 8) was carried out. The clay content is accounted for by calculating the average cation exchange capacity of the reservoir bed, using (9). The calculated cation exchange capacity determines the value of both diffusion-adsorption EMF coefficient of the reservoir, according to (6) and (7), and its resestivity (8). Kaolinite, illite and smectite are considered as clay minerals. Cation exchange capacities of these minerals and other modeling parameters are given in Table 2. Note that an increase in clay material content in the reservoir formation is accompanied by a change in porosity. It is seen that increasing the clay material content in the reservoir reservoir by 30% and decreasing the porosity by 5% relative to the case of pure sandstone leads to a more than twofold decrease in the amplitude of SP anomaly. At fixed values of clay material content and porosity, the amplitude of SP anomaly significantly decreases for sandstone reservoir with clay mineral having higher values of cation exchange capacity. As described in section 2.1, during the process of drilling the filtrate of fresh water clay mud penetrates into the permeable formation ( fig. 1), which leads to a change of the radial profile of the resistivity, characterizing the depth of invasion and the value of change of resistivity of the formed zone of drilling mud filtrate invasion . The latter, according to (Schlumberger, 1989), can be considered equal to 75% of the drilling mud resistivity in the borehole (1.5 Ohm m). Two-dimensional potential distribution in the studied area was analyzed in the presence of the borehole, invasion zone and unchanged part of the permeable formation (Fig. 9). Since the electrochemical source of SP is non-zero only in the area characterized by the salinity gradient, increasing the radius of the invasion zone of the mud filtrate is accompanied by a decrease in the amplitude of the SP anomaly.

Comparative analysis of synthetic and practical data
Using the proposed algorithm, simulation and comparative analysis of synthetic and practical SP logging data were performed. To initialize the model parameters, data on the diffusion-adsorption activity of rocks obtained from laboratory studies of core samples from intervals of Cretaceous deposits uncovered by wells of a number of fields in the Latitudinal Ob Region of Western Siberia are used. Presented results include the simulation of SP logging data for the Lower Cretaceous AS10, BS9-BS13 and BS16-BS20 oilbearing strata of the Achimov formation represented by irregularly interbedded sandstones and siltstones in shale sediments. In the presented interval of well A, the diffusion-adsorption activity of the sandstone is characterized by values from 0.75 to 26.7 mV and that of the siltstone by values from 16.5 to 31.1 mV. In the well B interval, this parameter has values: for sandstone from 6.3 to 21.3 mV, in the well C interval: from 3.65 to 41.7 mV and from 13.45 to 19.4 mV for sandstone and siltstone, respectively. Since these values of the diffusion-adsorption activity of rocks were measured under laboratory conditions at room temperature (25℃), while at the considered depth intervals the temperature reaches about 75℃, a temperature correction is taken into account. In this case, since is linearly dependent on temperature, its correct value can be determined by the formula: The results of numerical simulation of SP logging data using data on diffusion-adsorption activity of rocks obtained from laboratory core studies agree well with the practical data. Thus, the average value of relative discrepancy between synthetic and practical data is: well A -0.24, well B -0.11, well C -0.06.

Conclusions
SP logging is one of the main and demanded in practice methods of geophysical survey of geological sections exposed by oil and gas wells. For quantitative interpretation of SP logging data, correction charts based on description of physical and chemical processes in terms of equivalent electrical schemes are still widely used. However, this approach oversimplifies the origin of SP signals and is poorly applicable in real geological conditions. The most consistent approach to the description of the SP phenomenon at the moment is the formalism of coupled fluxes and nonequilibrium thermodynamics. Further development of the SP method is associated with the development of new ways of data analysis, based on the exact physical principles and implemented using modern computational methods.
We propose a new computational algorithm for modeling SP logging data in the axial symmetry approximation based on the finite element method. The physical-mathematical statement of the SP problem is carried out using the coupled fluxes formalism, which allows considering the effect of the electrochemical source of SP from the fundamental principles. The model of geological medium is presented, with the distribution of its electrophysical parameters approximated by sigmoidal functions that allow simulation of SP signals for a wide range of parameters characterizing real experimental conditions. The developed algorithm is tested by the simulation of SP signals in the model of fluid-saturated sandstone to establish the characteristic effects. They include the dependence of SP anomaly amplitude on formation water salinity, influence of permeable formation thickness on the SP log shape and dependence of SP anomaly amplitude on the relationship between the reservoir and drilling fluid resistivity. The used physical-mathematical model allows the simulation of SP signals in a clayey sandstone formation. The effectiveness of this approach is demonstrated by simulating the SP logging data in the reservoir, taking into account the content of the clay material, as well as its type. The results of simulation of SP signals depending on the depth of invasion of drilling mud filtrate are analyzed. The capabilities of the algorithm are demonstrated by a comparative analysis of synthetic and practical SP logging data, based on the results of laboratory studies of core from intervals of the Lower Cretaceous deposits exposed by wells in a number of Western Siberian fields.
The results obtained aim at the development of the SP method and allow us to formulate relevant directions of practical application. Since SP logging is performed in all wells and the data are widely available, including at the late stage of field development, it is extremely important to search for missed oil-bearing intervals. The solution of this problem becomes possible by re-interpretation of archival materials at a qualitatively new level with the use of modern software. One of the ways to increase informativity of the SP method while studying geological well sections is the use of oxidizers in drilling mud in order to contrast permeable intervals. For this purpose it is necessary to use softwarealgorithmic tools, taking into account full model description of physical and chemical processes, including redox, taking place both in a well and near-wellbore space. Another direction is the creation of realistic models of the geological environment to substantiate new geophysical methods of mineral exploration, including the method of pulsed electromagnetic sounding of the Bazhenov Formation with a system of highly deviated wells. All this is successfully facilitated by the active involvement of modern computational algorithms for solving forward and inverse SP problems in the class of multidimensional media models, including the one developed. Funding: Research work related to the construction of realistic models of the Lower Cretaceous deposits based on SP data and using the developed computational algorithm was carried out with the financial support of the Russian Science Foundation, grant number 19-77-20130.