Detector response model development

To obtain the model form of the detector response and ensure that the accuracy of the FFCM meets the requirements, a numerical model of the detector response was obtained based on the perturbation theory [23]. Perturbation theory is a mathematical method for determining the approximate solution of a problem under a small perturbation based on the exact solution of an existing relevant problem. Its underlying formula can be written as $N_{\text{P}}=N_{\text{B}}+\text{Δ}N\left(\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}\right)=N_{\text{B}}\left(1+\frac{\text{Δ}N\left(\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}\right)}{N_{\text{B}}}\right),$ (1) where $N_{\text{P}}$, $N_{\text{B}}$ are the detector counts under perturbation $P$ and base case $B$, respectively. $\text{Δ}N$ is the change in counts caused by the environmental perturbation $\text{Δ}\sigma$. According to perturbation theory, as $\text{Δ}N\to 0$, $\text{Δ}\sigma \to 0$ and $N^{\prime }=\underset{\text{Δ}\sigma \to 0}{\lim }\frac{\text{Δ}N}{\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}}$.

According to the Rytov approximation [24], Eq. (1) can be expressed as $N_{\text{P}}=N_{\text{B}}\text{Δ}\cdot e^{\left(\frac{\text{Δ}N\left(\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}\right)}{N_{\text{B}}}\right)}.$ (2)

Eq. (2) is expanded into a Taylor series at $N=N_{\text{B}}$ (corresponding to the environmental parameter ${\sigma }_{\text{B}}$) as shown in Eq. (3), with the expansion term $\sigma$. The first-order Taylor expansion is given by Eq. (4). $\begin{array}{c}{N}_{\text{P}} =N_{\text{B}}\left[1+{\displaystyle \sum _{m=1}^{\infty }\frac{1}{m!}}\frac{d^m{e}^{\left(\frac{\text{Δ}N}{N_{\text{B}}}\right)}}{d{\sigma }^m}\text{Δ}{\left(\sigma -{\sigma }_{\text{B}}\right)}^m\right]\\ =N_{\text{B}}\left[1+{\displaystyle \sum _{m=1}^{\infty }\frac{1}{m!}}\frac{d^m{e}^{\left(\frac{\text{Δ}N}{N_{\text{B}}}\right)}}{d{\sigma }^m}\text{Δ}{\left(\text{Δ}\sigma \right)}^m\right]\end{array}$ (3) $N_{\text{P}}=N_{\text{B}}+\frac{d{e}^{\left(\frac{\text{Δ}N}{N_{\text{B}}}\right)}}{d\sigma }\cdot \text{Δ}\sigma =N_{\text{B}}+N^{\prime }\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}$ (4)

N’ in Eq. (4) is obtained based on nuclear physics.

The detector count $N$ can be expressed as the sum of the particles in space from the source weighted by their contribution to the detector count, as shown in Eq. (5) [25–27]. $N=\int \text{d}r\int \text{d}E\text{Φ}\left(r,E,\text{Ω|}{r}_{\text{S}}\right){\text{Φ}}^{+}\left(r,E,\text{Ω|}{r}_{\text{R}}\right)=\int \text{d}r\int \text{d}EFSF,$ (5) where $\text{Φ}$ and ${\text{Φ}}^{+}$ are the forward flux and adjoint flux, respectively. $r_{\text{S}}$ and $r_{\text{R}}$ are the source and adjoint source locations, respectively. The flux sensitivity function is defined as $FSF=\text{Φ}\left(r,E,\text{Ω|}{r}_{\text{S}}\right){\text{Φ}}^{+}\left(r,E,\text{Ω|}{r}_{\text{R}}\right)$.

Based on perturbation theory, the FSF of the perturbation model is approximately equal to that of the base case, as shown in Eq. (6). $N^{\prime }=FSF\approx FS{F}_{\text{B}}$ (6)

According to Eqs. (4) and (6), the difference between the base and perturbed counts at a certain spatial position due to the single neutron energy is given by $\text{Δ}{N}_{{\text{P}}_{\text{SinglePosition}}} =N^{\prime }\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}=FS{F}_{\text{B}}\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}.$ (7)

Because the energy of the americium-beryllium neutron source is a fixed energy spectrum, it can be considered as a linear superposition of neutron sources with multiple single energies. The detector has a certain volume; therefore, it must be integrated over space, as shown in Eq. (8). $\text{Δ}{N}_{\text{P}}=\int \text{d}r\int \text{d}E\text{Δ}{N}_{{\text{P}}_{\text{SinglePosition}}}=\int \text{d}r\int \text{d}EFS{F}_{\text{B}}\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}$ (8)

According to Eqs. (4) and (8), the detector count in the case of perturbation is given by Eq. (9). $\begin{array}{c}{N}_{\text{P}}=N_{\text{B}}+\int \text{d}r\int \text{d}EFS{F}_{\text{B}}\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}=N_{\text{B}}\\ +\int \text{d}r\int \text{d}E{\text{Φ}}_{\text{B}}\left(r,E,\text{Ω|}{r}_{\text{S}}\right){\text{Φ}}_{\text{B}}^{+}\left(r,E,\text{Ω|}{r}_{\text{R}}\right)\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}\end{array}$ (9) where $\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}$ is derived based on the volumetric detection constraints method, which is described in detail in Section 2.2; the base case’s forward flux ${\text{Φ}}_{\text{B}}$ and adjoint flux ${\text{Φ}}_{\text{B}}^{+}$ are solved based on the time-independent Boltzmann equation, as shown in Eq. (10) [28, 29], which describes the entire process of neutrons attenuating in the formation and finally reaching the detector under neutron porosity logging using a radioactive neutron source. $\begin{array}{l}\left[\text{Ω}\cdot \nabla +{\mu }_{\text{t}}\left(r,E\right)\right]\text{Φ}\left(r,E,\text{Ω}\right)\\ ={\displaystyle \int \text{d}{E}^{\prime }}{\displaystyle \int \text{d}{\text{Ω}}^{\prime }}{\mu }_{\text{s}}\left(r,E^{\prime }\to E,{\text{Ω}}^{\prime }\to \text{Ω}\right)\text{Φ}\left(r,E^{\prime },{\text{Ω}}^{\prime }\right)+S\end{array}$ (10) where $\text{Φ}$ is the neutron flux at position $r$ in the angular direction $\text{Ω}$ and with energy $E$; ${\text{μ}}_{\text{t}}$ is the total reaction cross section, which can be expressed as the inverse of the migration length $L_{\text{m}}$; ${\mu }_{\text{s}}\left(r,E^{\prime }\to E,{\text{Ω}}^{\prime }\to \text{Ω}\right)$ is the scattering cross section from energy $E^{\prime }$ to $E$ and angular direction from ${\text{Ω}}^{\prime }$ to $\text{Ω}$; and $\nabla$ is a Laplacian operator.

Considering only the monoenergetic neutron group, $\Phi \left(r,E,\text{Ω}\right)$ can be written as $\Phi \left(r,\text{Ω}\right)$. Based on the simple diffusion approximation, i.e., no dependence of the angular flux on the direction $\text{Ω}$, Eq. (10) can be simplified to a neutron diffusion equation, as shown in Eq. (11) [30]. $\left[D{\nabla }^2-{\mu }_{\text{t}}\right]\text{Φ}\left(r\right)=-S$ (11) $D=\frac{1}{3{\mu }_{\text{t}}}$ (12) where $D$ is the diffusion coefficient, given by Eq. (12).

Based on the one-group neutron diffusion model [31] and a neutron source located at $r_{\text{S}}$, Eq. (11) can be expressed as $\left[D{\nabla }^2-\sigma \right]\text{Φ}\left(r\text{|}{r}_{\text{S}}\right)=-S$ (13) where σ is the perturbation parameter with $\sigma =\frac{1}{L_{\text{m}}}$ [32]. $\text{Φ}\left(r\text{|}{r}_{\text{S}}\right)$ is the flux field existing at $r$ generated by a positive point source of $S$ intensity located at $r_{\text{S}}$, i.e., the spatially diffuse scalar flux at observation position $r$.

The neutron flux equation can be solved using Eq. (13) as expressed in Eq. (14). $\text{Φ}\left(r\text{|}{r}_{\text{S}}\right)=\frac{S}{4\pi D\left|r-r_{\text{S}}\right|}{e}^{-\sigma \left|r-r_{\text{S}}\right|}$ (14)

Combining Eqs. (5) and (14), when the position and energy are fixed, the detector count $N$ is related to the perturbation $\sigma =\frac{1}{L_{\text{m}}}$ of the environment.

The numerical model of the detector count can be expressed as $N_{\text{P}} =N_{\text{B}}+\int \text{d}r\int \text{d}E\frac{S^{\text{B}}}{4\pi D^{\text{B}}\left|r-r_{{\text{S}}^{\text{B}}}\right|}{e}^{-\frac{\left|r-r_{{\text{S}}^{\text{B}}}\right|}{L_{\text{m}}^{\text{B}}}}\cdot \frac{S^{\text{B}+}}{4\pi D^{\text{B}+}\left|r-r_{{\text{S}}^{\text{B}+}}\right|}{e}^{--\frac{\left|r-r_{{\text{S}}^{\text{B}+}}\right|}{L_{\text{m}}^{\text{B}}}}\frac{L_{\text{m}}^{\text{B}}}{\text{Δ}{L}_{\text{m}}},$ (15) where $S^{\text{B}}$ and $S^{\text{B}+}$ are the source and adjoint source of base case, respectively. $L_{\text{m}}^{\text{B}}$ is the migration length of the standard model. $\frac{L_{\text{m}}^{\text{B}}}{\text{Δ}{L}_{\text{m}}}$ is derived based on the detection volumetric constraints method, which is described in detail in next section.

Detection volumetric constraints method

The detection volume contains important information such as the tool position, tool direction, and angle of the wells. For example, in the case of a high-angle well, the detection volume of the tool is inclined and the proportion of each layer in the detection volume differs from that of vertical wells. Therefore, the detection volume can be introduced to constrain the effect of perturbation to obtain $\frac{L_{\text{m}}^{\text{B}}}{\text{Δ}{L}_{\text{m}}}$ in Eq. (15) and optimize the detector response model to achieve a high-accuracy calculation.

An open-hole well consists mainly of two parts: a borehole filled with water and a formation. The hydrogen content of the formation rock mainly originates from the water in the pores; therefore, the hydrogen content of the borehole fluid is much larger than that of the formation rock, and the borehole fluid has a greater ability to slow down neutrons. Therefore, the particle transport path of neutrons in the borehole is shorter than that in the formation. According to Eq. (14), it can also be proved that the particle transport trace in the borehole is shorter when the neutron flux decays to the same value. This implies that the detection volume can be approximated as an ellipsoid, which is further demonstrated using the actual tool in Sect. 3.1.

Assuming that the well slope angles in the X, Y, and Z directions are ${\theta }_x$, ${\theta }_y$, and ${\theta }_z$, respectively, the detection volume equation is given by Eq. (16). $V_{\text{Detection}}=\frac{{\left(x^{\prime }-C_x\right)}^2}{a^2}+\frac{{\left(y^{\prime }-C_y\right)}^2}{b^2}+\frac{{\left(z^{\prime }-C_z\right)}^2}{c^2}$ (16) $\begin{array}{c}\left(x^{\prime },y^{\prime },z^{\prime }\right)=\left(x,y,z\right)\left[\begin{array}{lll}1\hfill & 0\hfill & 0\hfill \\ 0\hfill & \cos {\theta }_x\hfill & -\sin {\theta }_x\hfill \\ 0\hfill & \sin {\theta }_x\hfill & \cos {\theta }_x\hfill \end{array}\right]\\ \left[\begin{array}{lll}\cos {\theta }_y\hfill & 0\hfill & \sin {\theta }_y\hfill \\ 0\hfill & 1\hfill & 0\hfill \\ -\sin {\theta }_y\hfill & 0\hfill & \cos {\theta }_y\hfill \end{array}\right]\left[\begin{array}{lll}\cos {\theta }_z\hfill & -\sin {\theta }_z\hfill & 0\hfill \\ \sin {\theta }_z\hfill & \cos {\theta }_z\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill \end{array}\right]\end{array}$ (17) where $V_{\text{Detection}}$ is the detection volume; $\left(C_x,C_y,C_z\right)$ are the center coordinates of the detection volume; $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ are the coordinates of the detection volume equation after rotation; $a$, $b$, and $c$ denote half of the length of the ellipsoid along the three coordinate axes.

When the borehole diameter and lithology are fixed, the variation in counts is caused by formation porosity. Therefore, the corrected migration length value $L_{{\text{m}}_{\text{final}}}$ and the correction value $\text{Δ}{L}_{\text{m}}$ can be obtained using Eqs. (18) and (19), respectively. Owing to the difference in the mechanical structure, the correction values are different for various detectors of different tools. $L_{{\text{m}}_{\text{final}}} =\frac{N_{\text{P}}-N_{\text{B}}}{{{\displaystyle \int }}_{\text{r}\in {\text{V}}_{\text{formation}}}^{ }\text{d}r\int \text{d}EFS{F}_{\text{B}}}+{\sigma }_{\text{B}}$ (18) $\text{Δ}{L}_{\text{m}} =L_{{\text{m}}_{\text{final}}}-L_{\text{m}}^0$ (19) where $L_{\text{m}}^0$ is the initial migration length calculated in an infinite medium with a point source.

Based on Eq. (19), the migration length correction formula $f$ for various detectors in each base case can be fitted using Eq. (20), and the fitted data are the eigenvalues and adjacent data of the standard model. $L_{{\text{m}}_{\text{final}}}=L_{\text{m}}^0+\text{Δ}{L}_{\text{m}}=f\left(L_{\text{m}}\right)$ (20)

When the vertical resolution of the logging tool is not sufficiently high, or when thin layers are present, the tool detects signals from multiple geological layers simultaneously, as shown in Fig. 2, which is commonly referred to as a pile-up. This may result in logging data that do not accurately represent the characteristics of individual layers but rather a combination of their characteristics. Therefore, to accurately calculate the detector responses, it is necessary to consider the influence weight of the different geological layers detected by the detector in the response, which is achieved by weighting the migration length. Within the effective detection volume，let the detected volume of the $i^{\text{th}}$ rock layer be $V_{{\text{for}}_i} \left(i=1,2,\cdots ,n\right)$, the corresponding migration length of the $i^{\text{th}}$ rock layer be $L_{{\text{m}}_{{\text{for}}_i}} \left(i=1,2,\cdots ,n\right)$, and the detected volume of the borehole liquid be $V_{\text{fluid}}$. A two-dimensional schematic of the detection volume for multiple rock layers is shown in Fig. 2.

Fig. 2Full-screen View

(Color online) Two-dimensional detection volume schematic

The detection volumes for different rock layers and borehole fluids can be calculated using Eqs. (21) and (22), respectively. $V_{{\text{for}}_i}={\displaystyle \underset{V}{\iiint }\left\{\left[{x^{\prime }}^2+{y^{\prime }}^2-{\left(\frac{\text{BD}}{2}\right)}^2\le 0\right]\cap V_{\text{Investigation}}\right\} \text{d}V}$ (21) $V_{\text{fluid}}={\displaystyle \underset{V}{\iiint }\left\{\left[{x^{\prime }}^2+{y^{\prime }}^2-{\left(\frac{\text{BD}}{2}\right)}^2\ge 0\right]\cap V_{\text{Investigation}}\right\} \text{d}V}$ (22) where BD denotes the borehole diameter and $\cap$ denotes the intersection of the detection volume with the formation or borehole.

When the borehole diameter is not equal to the standard borehole diameter, the borehole fluid volume around the tool changes, affecting the neutron transport and causing a variation in the detector count. The standard model is based on a fixed borehole diameter and does not consider the influence of the borehole fluid. When the borehole diameter changes, the weight of the borehole fluid influence is no longer equal to zero. A migration length weighting method is given by Eq. (23). $L_{\text{m}}^1={\omega }_{\text{fluid}}\cdot \left(L_{{\text{m}}_{\text{for}}}^0-L_{{\text{m}}_{\text{fluid}}}^0\right)+L_{{\text{m}}_{\text{for}}}^0$ (23) ${\omega }_{\text{fluid}}=g\left(\frac{V_{{\text{fluid}}_{\text{B}}}-V_{\text{fluid}}}{V_{{\text{fluid}}_{\text{B}}}},L_{{\text{m}}_{\text{for}}}^0\right)$ (24) where $L_{\text{m}}^1$ is the weighted migration length of the impact from the borehole diameter constrained by the detection volume. ${\omega }_{\text{fluid}}$ is the weight reflecting the variation of borehole diameter and can be calculated using Eq. (24). The function $g$ can be obtained by fitting the relative variation of the borehole fluid volume and the formation migration length. $\left(L_{{\text{m}}_{\text{for}}}^0-L_{{\text{m}}_{\text{fluid}}}^0\right)$ represents the effect of borehole fluid and formation replacing each other.

When the detector detects multiple layers of rock, the migration length must be weighted based on the detection volume of the different rock layers and $L_{\text{m}}^1$. The weighting equation is given by Eq. (25). $L_{\text{m}}^2={\displaystyle \sum _{i=1}^n{\omega }_{{\text{for}}_i}} L_{{\text{m}}_i}^1$ (25) ${\omega }_{{\text{for}}_i}=\frac{V_{{\text{for}}_i}}{{\displaystyle {\sum }_{i=1}^n{V}_{{\text{for}}_i}}}$ (26) where $L_{\text{m}}^2$ is the weighted migration length of multi-formation impact. $L_{{\text{m}}_i}^1$ is the weighted migration length of borehole diameter impact of the $i^{\text{th}}$ formation. ${\omega }_{{\text{for}}_i}$ is the weight of the $i^{\text{th}}$ formation.

Combining Eqs. (20), (23), and (25), the final migration length can be obtained as $L_{{\text{m}}_{\text{final}}}=f\left({\displaystyle \sum _{i=1}^n{\omega }_{{\text{for}}_i}} \left[g\left(\frac{V_{{\text{fluid}}_{\text{B}}}-V_{\text{fluid}}}{V_{{\text{fluid}}_{\text{B}}}},L_{{\text{m}}_{\text{for}}}^0\right)\cdot \left(L_{{\text{m}}_{\text{for}}}^0-L_{{\text{m}}_{\text{fluid}}}^0\right)+L_{{\text{m}}_{\text{for}}}^0\right]\right)$ (27)

To select the base case that best matches the environmental perturbation, the distance method shown in Eq. (28) is used. $\text{Basecase}=\arg \min \left(L_{\text{m}}^{\text{B}}-L_{{\text{m}}_{\text{final}}}\right)$ (28)

According to Eqs. (9) and (27), the detector response using volumetric detection constraints is given by $\begin{array}{c}{N}_{\text{P}}=N_{\text{B}}+\int \text{d}r\int \text{d}EFS{F}_{\text{B}}\frac{\text{Δ}\sigma }{{\sigma }_{\text{B}}}=N_{\text{B}}+\int \text{d}r\int \text{d}EFS{F}_{\text{B}}\frac{L_{\text{m}}^{\text{B}}}{L_{\text{m}}^{\text{B}}-L_{{\text{m}}_{\text{final}}}}\\ =N_{\text{B}}+\int \text{d}r\int \text{d}EFS{F}_{\text{B}}\frac{L_{\text{m}}^B}{L_{\text{m}}^{\text{B}}-f\left({\displaystyle {\sum }_{i=1}^n{\omega }_{{\text{for}}_i}} \left[g\left(\frac{V_{{\text{fluid}}_{\text{B}}}-V_{\text{fluid}}}{V_{{\text{fluid}}_{\text{B}}}},L_{{\text{m}}_{{\text{for}}_i}}^0\right)\cdot \left(L_{m_{\text{for}}}^0-L_{{\text{m}}_{\text{fluid}}}^0\right)+L_{{\text{m}}_{{\text{for}}_i}}^0\right]\right)}\end{array}$ (29)

Implementation in a neutron porosity tool

Neutron porosity logging is an important technique for evaluating reservoirs. This method is applicable to any wireline porosity tool, and is discussed in this study based on a wireline neutron porosity tool. The neutron porosity tool was offset to obtain more information about the formation during logging. As shown in Fig. 3, the tool with a ²⁴¹Am-Be source is placed in a fluid-filled borehole and attached to the borehole wall [33]. Two thermal neutron detectors filled with ³He are placed 380.2 mm and 618 mm from the source, and a shield made of a tungsten-nickel alloy is placed between the source and detectors to shield the axial gamma rays. The detailed parameters of the neutron porosity tool are listed in Table 1. During logging, neutrons are emitted into the formation and reach the detector for collection through processes such as elastic scattering, inelastic scattering, and capture. To simulate deep penetration problems, the Geant4 physics list called “shielding” is used, which is specifically designed for this purpose.

Fig. 3Full-screen View

(Color online) Neutron porosity logging model

To ensure the accuracy of the model in a short computation time, the neutron logging model is described by a two-dimensional Cartesian coordinate system [34], that is, integration over X, where the formation is vertical to the Z axis and the neutron porosity tool is symmetrical about the Z-Y plane. The space around the tool is divided into grids with axial and radial boundaries to collect neutron flux. To build an FFCM covering multiple environments, a data library containing various parameters of borehole diameter, borehole fluid salinity, and porosity was constructed, as listed in Table 2. The corresponding migration length values were calculated for different formation porosities. In this case, the standard borehole diameter for the base case was 311.2 mm and other environmental parameters of the base case are listed in Table 2. The total computational time for the database was 210 h. It should be noted that as the number of environmental factors increases, larger databases are required to ensure computational accuracy.

In this study, Geant4 was used to construct a neutron porosity logging model and obtain a data library. The conversion equation for the detector ratio and porosity is given by Eq. (30) [35, 36] and can be obtained from the detector ratio data in the data library. ${\varphi }_{\text{f}}=1.16\cdot R^2+2.05\cdot R-8.11$ (30) where ${\varphi }_{\text{f}}$ is the measured porosity and $R$ is the near-to-far detector ratio.

The established data library contains the response of the model in different environments, including the forward flux, adjoint flux, and FSFs. The two-dimensional distributions of the forward and adjoint fluxes were obtained using Eq. (14). To better show the distribution, the normalized flux distribution and FSF are shown in Figs. 4, 5, 6, 7, 8. As shown in Fig. 4, the forward neutron flux gradually decreases in all directions centered at the location of the source, and is more widely distributed in the formation direction because the slowing ability of the formation to neutrons is weaker than that of the borehole fluid. The adjoint flux is acquired by setting the adjoint source on the detector; thus, the distribution of the flux centered on the detector and decreasing in all directions is presented in Figs. 5 and 6. Because the increased porosity of the formation increases the hydrogen content of the formation and thus enhances the scattering effect on neutrons, regardless of the forward flux or adjoint flux, the normalized neutron flux distribution is more concentrated at the source center or adjoint source center with increasing porosity of the formation. As shown in Fig. 7 and 8, greater flux sensitivity values are available in the source and detector regions because these regions have the highest neutron flux recorded in the detector, which reflects the contribution of neutrons to the detector in each region. The flux sensitivity functions reflect the common characteristics of the forward and adjoint fluxes; that is, the sensitivity values are more concentrated near the source and detector as the porosity of the formation increases, and the sensitivity values are not concentrated in the part of the tool between the source and detector because of the shielding effect inside the tool [37].

Fig. 4Full-screen View

(Color online) Normalized forward flux distribution at 311.2 mm borehole diameter with different formation porosity

Fig. 5Full-screen View

(Color online) Normalized adjoint flux distribution of Near Detector at 311.2 mm borehole diameter with different formation porosity

Fig. 6Full-screen View

(Color online) Normalized adjoint flux distribution of Far Detector at 311.2 mm borehole diameter with different formation porosity

Fig. 7Full-screen View

(Color online) Normalized flux sensitivity function of Near Detector at 311.2 mm borehole diameter with different formation porosity

Fig. 8Full-screen View

(Color online) Normalized flux sensitivity function of Far Detector at 311.2 mm borehole diameter with different formation porosity

After obtaining the FSF, it is necessary to obtain the three parameters, $a$, $b$, and $c$ for the detection volume in Eq. (16). The transport tracks of the neutron particles that arrived at the detector were collected using Geant4. The maximum distance reached by 90% of the detected neutrons was defined as the maximum detection volume. The $a$, $b$, and $c$ of various detectors in different environments are listed in Table 3.

From Table 3, without considering the angle of the well and direction of the tool, the detection volume is mainly related to the location of the detector and formation migration length, which are slightly influenced by the borehole diameter because the tool is attached to the borehole wall. Therefore, the detection volume can be determined based on the migration length of the base case to constrain environmental disturbance. Thus, a preset data library is built, and the FFCM can be computed according to the pseudocode shown in APPENDIX A.

Simulated data verification

To ensure the validity of the FFCM, a test database containing variations in environmental parameters, such as borehole diameter, lithology, and porosity, was established. Setting the number of source particles to $2\times {10}^9$ so that the statistical error in detector response counts is less than 2% and converges [8], and each model requires more than 5 h of computing using 110 threads on Intel® Xeon(R) Platinum 8280 CPU @ 2.60GHz × 224. The formation in the logging model is composed of alternating sandstone or limestone with different porosities, and the borehole is filled with fresh water. The wireline neutron tool moved continuously from bottom to top, sampling the formation at 100 mm intervals to obtain the detector response. The method was validated in two cases: a vertical well (Case 1) and a high-angle well (Case 2).

Case 1 was designed using a model with different lithologies, formation thicknesses, and porosities for a nonstandard borehole diameter of 350 mm. A schematic of the model is shown in Fig. 9, and its specific parameters are listed in Table 4. The calculation results of the FFCM in Case 1 are shown in Fig. 10, where Track 1 and Track 2 are the normalized counts of the near and far detectors, respectively, and the normalization formula is shown in Eq. (31), which represents the probability of a neutron emitted from the source reaching the detector. MC Near Det Counts and FFCM Near Det Counts represent the counts of near detectors from MC and FFCM, respectively, and MC Far Det Counts and FFCM Far Det Counts represent the counts of far detectors from MC and FFCM, respectively. Track 3 shows the formation porosity logs of the MC and FFCM obtained using Eq. (30), which are called the MC Porosity and FFCM Porosity, respectively. Track 4 shows the absolute and relative errors of the FFCM and MC simulation results, where the red lines are the error lines of 1 p.u. and 8%, respectively. $N_{\text{normalized}}=\frac{N}{N_{\text{S}}}$ (31) where $N_{\text{normalized}}$ is the normalized count. $N_{\text{S}}$ and $N$ are the number of source particles and the detector count, respectively.

Fig. 9Full-screen View

(Color online) Case 1 model schematic

Fig. 10Full-screen View

Calculation results of the FFCM in Case 1

As shown in Fig. 10, 50 measurement points were obtained by sampling at a sampling interval of 100 mm in the Case 1 formation model at a depth of 4.9 m. As shown in Tracks 1–3, the results of the FFCM calculations are in good agreement with the MC simulation results in terms of trends and values, both for detector counts and formation porosity. As shown in Track 4, the absolute error of FFCM is mostly below 1 p.u., and the relative error is mostly below 6%, with a maximum absolute error of 1.19 p.u., the maximum relative error of 6.80%, and the root mean square error (RMSE) of 0.62 p.u.

Case 2 is a high-angle well model with the Y-axis angle ${\theta }_y$ being 10°, while the other axis tilt angles ${\theta }_x$ and ${\theta }_z$ are zero and borehole diameter is 311.2 mm. The working state of the wireline neutron log was simulated to study the results of FFCM with different formation thicknesses and porosities, as shown in Fig. 11. A schematic of the Case 2 model is presented in Fig. 11 and Table 5.

Fig. 11Full-screen View

(Color online) Case 2 model schematic

As shown in Fig. 12, 46 measurement points were obtained by sampling at a sampling interval of 100 mm in the Case 2 formation model, which had a depth of 4.5 m. In Tracks 1–3, the calculated logs of FFCM and MC coincide well and show a good matching performance. According to Track 4, the absolute error of FFCM is mostly below 1 p.u. and the relative error is mostly below 6%, with a maximum absolute error of 1.09 p.u., a maximum relative error of 6.54%, and an RMSE of 0.56 p.u. Based on the validation of two simulation cases, the good agreement between FFCM and simulations verifies the validity and reliability of the novel algorithm.

Fig. 12Full-screen View

(Color online) Calculation results of the FFCM in Case 2

Once the database has been constructed, FFCM takes less than 1 s to execute, which provides a significant advantage in complicated scenarios that typically require heavy amounts of Monte Carlo modeling. For instance, in Case 2, the conventional Monte Carlo method requires the computation of a minimum of 46 models with the correct parameters, which would take 230 h with 110 threads on Intel® Xeon(R) Platinum 8280 CPU @ 2.60 GHz × 224. Considering the additional calculations of other model parameters during the curve-matching process or longer, well-segmented processing requirements, the Monte Carlo method requires more time. However, using FFCM, each model in Cases 1 and 2 can be computed in less than 1 s regardless of the computational hardware parameters, enabling real-time data processing.

Experimental part: field application

To evaluate the performance of FFCM, it was applied to a field well. The experimental well is located in the southern part of Liaodong Bay in the Bohai oil field and is drilled with a 241.3 mm (9.5 inch) bit. During logging, the tool slides against the borehole wall from the bottom to the top and is sampled at 76 mm intervals. The sieve residue log shows that the logging section is dominated by sandy mudstone formations, and borehole fluid sampling shows that borehole fluid has a density of 1.18 g/cm³ with a salinity of 94 kppm. The compositions of the formations with different porosities and borehole fluids were obtained to calculate the initial migration length $L_m^0$. In addition, the experimental well requires real-time measurements of the borehole diameter and angle of the well, which provides more accurate inputs to the method. The pseudocode of the method is shown in APPENDIX A, and its input is a formation model that includes the lithology, borehole parameters, and formation properties. By iteratively updating the input formation model to ensure that the calculated porosity log of the FFCM matches the measured porosity log of the tool as much as possible, the final porosity model with multiple layers can be used to assist in the rapid completion of formation analysis to support the formation of a layered interpretation. By utilizing the FFCM, as described above, the limitations of the baseline calculation method can be overcome, and the pile-up problem can be resolved, resulting in accurate properties of different geological layers. However, although the proposed method may introduce some errors compared to a full Monte Carlo simulation, it remains the optimal choice that balances time and accuracy considerations.

The formation model obtained using the FFCM method in the X475-X515 well section was compared with the interpretation results, as shown in Fig. 13. Track 1 shows the borehole diameter log (CAL), natural gamma log (GR), and DEVOD log, reflecting the angle of the well. RD, RS, and RMSF in Track 2 denote the investigated deep double lateral resistivity log, shallow investigated double lateral resistivity log, and micro-spherically focused resistivity log, respectively. DT24 in Track 3 represents the acoustic slowness. Measured Porosity and Measured Density are the measured values from the neutron porosity and gamma density tools, respectively. Tracks 2 and 3 are important bases for the interpretation. The FFCM-calculated porosity (FFCM Porosity) log in Track 4 was obtained by iteratively and continuously updating the formation model, which matched well with the measured porosity. The FFCM Formation Model of Track 5 is the final formation model corresponding to Track 4, and the Interpretation Model is the formation model obtained by synthesizing and analyzing a variety of measured logs; its detailed parameters are listed in Table 6.

Fig. 13Full-screen View

(Color online) Calculation results of the FFCM in field

The GR log value of Track 1 in Fig. 13 was approximately 75 API, which is consistent with the sieve residue log, indicating that the lithology was mainly composed of sandy mudstone. The borehole size of the measurement section was maintained at 9–13 in because of the enlargement of the borehole, and the angle of the well in this section increased from 25°to 30°, as indicated by the DEVOD log of Track 1. After multiple updates of the formation model, the FFCM-calculated porosity of Track 4 matched the measured porosity well. In addition, the formation porosity model of Track 5 was in good agreement with the model results provided by the interpreters, proving the validity and reliability of the FFCM constrained by the detection volume.