Introduction

During the various stages of oilfield exploration and development, due to the influence of multiple internal and external factors, the original physical, chemical, thermodynamic, and hydrodynamic equilibrium states of the reservoir may change, inevitably causing a decrease in reservoir permeability, blocking fluid flow, resulting in reservoir damage and a decline in oil well production, and even “killing” the reservoir1,2,3,4,5. The causes of reservoir damage are diverse and dynamically changing6. Once reservoir damage occurs, it is necessary to adopt appropriate measures to unblock the fluid flow channels according to the reservoir damage situation to restore or improve oil well production and water injection capacity2,7. Therefore, clarifying the types of reservoir damage and the proportion of each type, as well as the spatial distribution and temporal variation of reservoir damage, directly affects the effectiveness of unblocking and increasing production. Accurate and quantitative simulation of reservoir damage types, degrees, and patterns can also improve the accuracy of reservoir numerical simulation, improve the prediction accuracy of reservoir “sweet spots”, achieve scientific well placement, and ultimately improve the drilling success rate and economic benefits.

At present, there are many methods for diagnosing reservoir damage, which can be divided into field methods mainly based on well testing evaluation, logging evaluation, and oil and gas well production decline curve analysis, as well as laboratory evaluation methods mainly based on core flow experiments2,3,4,7. They either use the total skin factor (which equals the sum of the real damage skin factor and the apparent skin factor) that cannot reflect the true damage situation, or characterize the overall damage degree with reservoir permeability, but they cannot accurately determine the types of reservoir damage and the degree of damage of each type on the reservoir, nor can they achieve spatiotemporal evolution simulation. To address this, in 1994, Jiang Guancheng, Yan Jienian, and others2,8,9 first attempted to combine laboratory experimental data with field evaluation to create a mathematical modeling method, realizing quantitative diagnosis of various damage factors, but it has drawbacks such as difficulty in obtaining representative core samples, long experimental cycles, large workload, and poor universality. In 2015, D. Mikhailov and others2,10 also used a combination of mathematical simulation and special experiments to establish a dynamic characteristic model describing reservoir performance changes in the near-wellbore region caused by drilling fluid component invasion, which can more clearly reflect the damage mechanism and dynamic characteristics, but it cannot separate the proportion of drilling fluid damage to the reservoir from the numerous damage factors. In 2022, addressing the challenge of reservoir damage caused by drilling fluid invasion in sandstone formations, Wang Xiang11 established the impact of drilling fluid invasion damage based on two-phase flow equations for oil and water and convection-diffusion equations. Meanwhile, Lin YF et al.12 considered the combined effects of laminar and non-Darcy flow using the Forchheimer equation to establish a testing and analysis model for drilling fluid that takes reservoir damage into account. The limitations of both approaches are their failure to consider the transport of solid phases in drilling fluids and within the reservoir. El-Monier13 employed data mining and machine learning techniques for analyzing reservoir damage caused by hydraulic fracturing, applying image recognition technology at the fracture scale to model reservoir damage. However, it did not account for reservoir damage occurring on the pore-scale. Li JS et al.13 developed a flow model that considers reservoir damage and, based on material balance, established a production model for horizontal wells, modeling reservoir damage and horizontal well productivity. Yet, it focused solely on the reservoir damage generated during the production phase. In recent research, Li J et al.14 established a predictive model for reservoir damage caused by solid phase invasion, with the advantage of considering the effects of pollution time and media tortuosity. Bui D, et al.15 used Computer Modeling Group software to simulate reservoir damage caused by fracturing fluid invasion and well closure. However, their work lacks corresponding mechanistic research. Dargi et al.16 employed artificial neural networks, random forests, and XGBoost combined with genetic algorithms to predict reservoir damage caused by acidizing procedures, achieving high accuracy, although some field parameters are challenging to obtain. Guancheng Jiang et al.17 established numerical models for near-wellbore reservoir damage in horizontal wells using a polar coordinate transformation and a series of convection-diffusion equations, yet the types of reservoir damage considered were not comprehensive. Therefore, a common approach in recent research on reservoir damage diagnosis is to use numerical modeling to describe reservoir damage, although the perspective (physical mechanism perspective, statistical perspective) and methods (artificial intelligence algorithms, numerical simulation software) for building numerical models may vary. However, a limitation of recent research is that it only considers modeling of a single stage, such as drilling or production, without taking into account the entire process from drilling, completion, workover, production, to enhanced recovery.

The numerical simulation method, with its advantages of comprehensively describing reservoir damage mechanisms, various damage characteristics, and dynamic damage simulation, has gained widespread recognition internationally and has become a research hotspot and cutting-edge research direction. The simulation methods have gradually evolved from continuum mechanics (hydraulic equations) to more complex probabilistic process numerical systems18,19. Simulation technology mainly focuses on the fluid flow behavior in porous media and the interaction between fluid and rock surfaces, conducting research through the establishment of mathematical models, with strong theoretical foundations but limited practicality.

Production demands not only require quantitative diagnosis of damage degrees and types but also dynamic simulation and emulation of the damage generation and development process, which were unattainable with previous technologies. Numerical simulation technology involves iterative calculations in time and space, possessing a natural advantage in simulating spatiotemporal change patterns. Therefore, based on the previous research and starting from the microscopic mechanism of reservoir damage, this paper establishes a numerical model containing various damage factors according to the mass conservation principle and diffusion laws in the reservoir damage process. We systematically study the temporal-spatial evolution laws of reservoir damage on a large scale (the spatial area extending from the wellbore center to the reservoir). Specifically, we simulate the interaction between reservoir permeability, rock surface properties, fluid properties, and external factors (such as pressure difference, temperature, and foreign fluids) in space and time dimensions, resulting in changes in reservoir permeability and real damage skin factor. This approach addresses the shortcomings of previous diagnostic technologies that could only provide vague qualitative conclusions and could not provide quantitative information on the damage degree and damage zone radius of a specific damage type. This new technology, which cannot be replaced by field methods and experimental methods, describes reservoir damage situations in a more refined, dynamic, accurate, and practical manner20,21.

Methodology

Establishment of quantitative simulation control equation for spatial-temporal evolution of common types of damage

The damage issues of different reservoirs have both specificity and commonality. Referring to the investigation of reservoir damage situations in more than 4,000 wells by the American Core Laboratories, and based on the proposal that there are usually 8 types of damage during well construction and oil field development stages2,22, this paper combines the specificity of reservoir damage at various operational stages and conducts spatiotemporal evolution quantitative simulation technology research on 12 types of damage.

The governing equations of the 12 reservoir damage models are nonlinear convection-diffusion equations containing source-sink terms, as shown in Eq. 1:

$${d_a}\frac{{\partial u}}{{\partial t}}+\nabla \cdot \left( { - c\nabla u} \right)+\beta \cdot \nabla u=f$$
(1)

\(\nabla =[\frac{\partial }{{\partial x}},\frac{\partial }{{\partial y}},\frac{\partial }{{\partial z}}]\)

Where C is the diffusion coefficient. In the process of model establishment, C takes 1 to represent isotropy. \(\nabla\) is a Laplace operator and a second-order partial differential operator in n-dimensional Euclidean space. Therefore, in the one-dimensional case, the equation can be simplified to the following general form:

$$\frac{{\partial f}}{{\partial t}}+{a_a}\frac{{\partial f}}{{\partial x}}+{b_b}f={c_c}\frac{{{\partial ^2}f}}{{\partial {x^2}}}$$
(2)

Where \(\:{a}_{a},\:\:{b}_{b},\:\:{c}_{c}\) can be constants (such as diffusion coefficient) or functions (such as average reservoir velocity); f can be pressure, material concentration, stress, etc. For the convection-diffusion equation and the finite difference method, please refer to the paper written by Mark and Charles (2004)23

External solid particle damage

The essence of solid phase intrusion blockage is the migration24 and deposition of exogenous solid phase particles after entering the medium. Therefore, the core of the study is to establish a dynamic model of solid phase particle migration and deposition. Due to the advantage of phenomenological models that include many physical and chemical factors25,26,27, this paper establishes a spatiotemporal evolution control phenomenological model of the particle concentration distribution of foreign solid phase particles in the reservoir based on mass conservation and diffusion relationships. This model contains the flowing particle concentration C and the deposited particle concentration Cd, and by combining the relationship between the deposition concentration and the permeability damage, the spatiotemporal field distribution of permeability damage can be diagnosed.

According to mass conservation, the change in the flowing mass is equal to the negative change in the deposited mass, and the mass balance equation of the solid-liquid mixture and particles is obtained:

$$\frac{\partial }{{\partial t}}\left( {\rho \phi } \right)+\frac{\partial }{{\partial x}}\left( {\rho u} \right)= - \dot {m}\\ \frac{\partial }{{\partial t}}\left( {\rho \phi w} \right)+\frac{\partial }{{\partial x}}\left( {\rho uw+j} \right)= - \dot {m}$$
(3)

Where \(j= - \phi {\rho _L}D\frac{{\partial w}}{{\partial x}}\), \(D=\alpha v\), \(\dot {m} \equiv \frac{{\partial m}}{{\partial t}}\).

The relationship between the deposition volume concentration and the flow volume concentration, that is, the deposition equation ( connection condition )28,29 can be expressed as:

$$\frac{{\partial ({\rho _p}{C_d})}}{{\partial t}}=k(\rho uw+j)$$
(4)

Where \(k={k_0}{G_1}({C_d}){F_1}(T)\).

The relationship between mass concentration and volume concentration of flowing particles is:

$$w=\frac{{{\rho _p}}}{{{\rho _L}}}C$$
(5)

It is usually assumed that the fluid-particle mixing density ρ ≈ fluid density ρL.

Combining (3 ~ 5) and simplifying, the control equations of flow particle concentration C and deposition particle concentration Cd are obtained respectively:

$$\frac{{\partial C}}{{\partial t}}+\frac{v}{\tau }\left[ {1 - \left( {1 - \frac{{{\rho _p}}}{{{\rho _L}}}C} \right)k\alpha \tau } \right]\frac{{\partial C}}{{\partial x}}+\frac{{kv}}{\tau }\left( {1 - \frac{{{\rho _p}}}{{{\rho _L}}}C} \right)C=\alpha v\frac{{{\partial ^2}C}}{{\partial {x^2}}}$$
(6)
$$\frac{{\partial {C_d}}}{{\partial t}}=\frac{{vk\phi }}{\tau }[C - \alpha \tau \frac{{\partial C}}{{\partial x}}]$$
(7)

The relationship between porosity and permeability, and the relationship between penetration coefficient and porosity can be expressed as follows:

$$\frac{K}{{{K_0}}}={\left( {1 - \frac{{{C_d}}}{{{\phi _0}}}} \right)^{{m_k}}}\\ \frac{k}{{{k_0}}}={(1 - \frac{{{C_d}}}{{{C_{d\hbox{max} }}}})^{{m_k}}}$$
(8)

Pressure is the driving force that drives the solid-liquid mixture to continuously invade the reservoir from the wellbore. The pressure conduction equation and Darcy’s law describing the fluid entering the reservoir are:

$${\nabla ^2}P=\frac{{\phi \mu {c_t}}}{K}\frac{{\partial P}}{{\partial t}}\\ v= - \frac{{\tau K}}{{\mu \phi }}\nabla P$$
(9)

Equation (6 ~ 9) constitute the governing equation for the diagnosis of external solid blockage damage. The main physical quantities C、\({C_d}\)、V、P、K in the equation are time-space field functions.

In Eq. (2), when f is the flow particle concentration C, \(\:{a}_{a}=\frac{v}{\tau\:}\left[1-\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right)k\alpha\:\tau\:\right],\:{b}_{b}f=\frac{kv}{\tau\:}\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right),\:{c}_{c}=\alpha\:v\), this equation is the governing equation of the flow particle concentration C, that is, Eq. (6); when f is the concentration of deposited particles, \(\:{a}_{a}=vk\alpha\:\varphi\:\), \(\:{b}_{b}=-\frac{vk\varphi\:}{\tau\:}\), \(\:{c}_{c}=0\), this equation is the governing equation of the concentration of deposited particles, that is, Eq. (7). Limited to space, the control equations established below are described in this form.

Particle migration damage

When establishing a model for particle migration damage, the critical velocity for particle mobilization is first considered. Subsequently, the manner in which these internal particle sources alter the solid-liquid flow and deposition equation described in Eq. (3) is examined. The force and torque acting on elastic solid particles with a radius of \(\:{r}_{s}\) at the rough inner surface of rock pores are illustrated in Fig. 1. These particles, under the influence of fluid scouring and interactions with the rock surface, are subject to forces and moments. The critical velocity corresponds to the fluid velocity at which these forces and moments are precisely balanced.

Fig. 1
figure 1

Moment balance of particles in formation pore fluid.

Full size image

As shown in Fig. 1, the particles are subjected to a drag force \(\:{F}_{d}\) in the same direction as the flow velocity \(\:{v}_{t}\), as well as gravity \(\:{F}_{g}\), electrostatic force \(\:{F}_{e}\), and lifting force \(\:{F}_{l}\). The corresponding flow velocity when the moment is balanced is the critical flow velocity.

$$\:{F}_{d}\cdot\:{l}_{d}=\left({F}_{e}-{F}_{l}+{F}_{g}\right)\cdot\:{l}_{n}$$
(10)

Only when the actual velocity of the fluid in the reservoir exceeds the critical velocity \(\:{v}_{cr}\), the reservoir particles will migrate and become the source of migration material. In general, the closer to the center of the wellbore, the greater the fluid velocity is, and the particle migration area should be a ring belt near the wellbore.

According to the mass equation, let the amount of particles released by the particle source be \(\:{q}_{s}\), then \(\:{q}_{s}\) has the following properties:

$${q}_{s}:\left\{\begin{array}{c}>0,\quad v\ge\:{v}_{cr}\\\:=0,\quad v<{v}_{cr}\end{array}\right.$$
(11)

The Eq. (11) shows that only when the fluid flow rate exceeds the critical flow rate, the particles will start and join the fluid to participate in the migration, thereby increasing the mass of the fluid-solid mixture. Therefore, when there is particle migration, the Eq. (3) is rewritten as follows:

$$\:\frac{\partial\:}{\partial\:t}\left(\rho\:\varnothing\:\right)+\frac{\partial\:}{\partial\:x}\left(\rho\:u\right)=-\dot{m}+{q}_{s}$$
(12)
$$\:\frac{\partial\:}{\partial\:t}\left(\rho\:\varnothing\:w\right)+\frac{\partial\:}{\partial\:x}\left(\rho\:uw+j\right)=-\dot{m}+{q}_{s}$$
(13)

The simultaneous Eqs. (12) and (13) can be obtained:

$$\:{q}_{s}:\left\{\begin{array}{c}>0,\quad v\ge\:{v}_{cr}\\\:=0,\quad v<{v}_{cr}\end{array}\right.$$
(14)

Then the above geometric model is integrated into the form of Eq. (2). In Eq. (2), when f is the concentration of flowing particles C, \(\:{a}_{a}=\frac{v}{\tau\:}\left[1-\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right)k\alpha\:\tau\:\right],\:{b}_{b}f=(\frac{kv}{\tau\:}C-\frac{{q}_{s}}{{\rho\:}_{p}\varphi\:})\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right),\:{c}_{c}=\alpha\:v\), it means that when the fluid velocity exceeds the critical velocity, the particles will start and join the fluid to participate in the migration, thus increasing the mass of the fluid-solid mixture. The determination of the critical velocity can be referred to the method in References8,9,30,31,32,33,34. Combined with the formula (7 ~ 9), the numerical model of particle migration damage is obtained.

Inorganic precipitation damage

The common inorganic precipitations in the drilling process are mainly calcium carbonate, calcium sulfate, barium sulfate, strontium sulfate, etc. The mathematical model in this study mainly considers these four types of inorganic precipitations. When the foreign fluid is incompatible with the reservoir fluid, inorganic precipitation may occur, blocking the fluid flow channels (interconnected pores in rocks)35 and causing reservoir damage. Under certain temperature, pressure, and ion concentration, the occurrence of a precipitation reaction in the reservoir solution is usually judged by the scaling index, the expression for which is36:

$${I_S}\left( {P,T,{S_i}} \right)=\log \left( {\frac{{\left[ {{M_e}} \right]\left[ {{A_n}} \right]}}{{{K_c}\left( {P,T,{S_i}} \right)}}} \right)$$
(15)

Common inorganic precipitates include: calcium carbonate (CaCO3), calcium sulfate (CaSO4), strontium sulfate (SrSO4), and barium sulfate (BaSO4). The scaling indices of these four inorganic substances can be calculated according to the Tomson-Oddo calculation method36.

In Eq. (2), when f is the concentration of the i-th ion \(\:{C}_{i},\:\:{c}_{c}=\frac{D}{\varphi\:}\), \(\:{a}_{a}=\frac{u}{\varphi\:}\), \(\:{b}_{b}=\frac{\partial\:u}{\varphi\:\partial\:x}\), the equation right end plus ion concentration loss(\(\:-{k}_{i}{C}_{i}^{{\beta\:}_{i}}{C}_{j}^{{\beta\:}_{j}}\)), the equation of the i-th ion transport can be obtained37,38. For each time of the migration equation, it is necessary to judge whether the scaling index of the grid node satisfies the condition of Eq. (16):

$${k_i}=\left\{ \begin{gathered} 0,\quad _{I_S} \leqslant 0 \hfill \\ {k_{i0}},\quad_{I_S}>0 \hfill \\ \end{gathered} \right.$$
(16)

If \({k_i}=0\), It indicates that no precipitation reaction occurs in situ; if the precipitation reaction (Is > 0) occurs, the cumulative precipitation amount in the dt time step is:

$$d(r,t+dt)=d(r,t)+\frac{1}{\phi }\sum\limits_{i} {c{o_i}\Delta {C_i}(r,t)}$$
(17)

Therefore, the numerical model of inorganic scale damage is obtained by combining the migration equation and formula (9,16,17) of the i-th ion.

Clay hydration expansion damage

It is assumed that the diffusion of water molecules from the liquid phase in the pore to the solid phase through the solid-liquid interface obeys Fick diffusion. Assuming that the initial water volume fraction in the rock solid phase is \({c_0}\), the water volume fraction at a certain moment (dynamic) in the rock solid phase is C, the water volume fraction at a certain moment (dynamic) in the pore is \({c_1}\), the original porosity of the reservoir is \({\phi _0}\), and the irreducible water saturation is \({S_{wc}}\), then \({c_1}(t=0)={\phi _0}{S_{wc}}\). Establish a coordinate system z perpendicular to the solid-liquid interface pointing to the inside of the solid phase. z = 0 at the interface, z > 0 in the solid phase. According to Fick ‘s diffusion law, the diffusion equation of water molecules inside the wall is established31,39 :

$$\frac{{\partial c}}{{\partial t}}={D_w}\frac{{{\partial ^2}c}}{{\partial {z^2}}},0 \leqslant z<\infty ,t \geqslant 0$$
(18)

The initial condition of the equation is \(c(t=0)={c_0}\). The boundary condition is \(\:\dot{S}\equiv\:-{D}_{w}\frac{\partial\:c}{\partial\:z}={k}_{f}\left({c}_{1}-c\right){|}_{z=0}\), \(\mathop {\lim }\limits_{{z \to \infty }} c(z,t)={c_0}\). The expression of water absorption rate is derived by Laplace transform of Eq. (18) and its initial and boundary conditions.

$$\:\dot{S}\equiv\:-{{D}_{w}\frac{\partial\:c}{\partial\:z}|}_{z=0}=\left({c}_{1}-{c}_{0}\right){k}_{f}{e}^{\frac{{{k}_{f}}^{2}}{{D}_{w}}t}\text{e}\text{r}\text{f}\text{c}\left(\sqrt{\frac{{{k}_{f}}^{2}}{{D}_{w}}t}\right)$$
(19)

The change rate of porosity caused by clay expansion can be characterized by the following Eqs. .

$$\:-\frac{\partial\:\varnothing\:}{\partial\:t}={\uplambda\:}\dot{S}$$
(20)

Where λ satisfies: \(\lambda =\frac{{k^{\prime}P{I^{2.44}}C{c^{3.44}}}}{{{{\left( {Cc - 10} \right)}^{2.44}}}}\). If \(PI<1\sim 2\), reservoir rock is brittle rock. \(2<PI<6\), plastic brittle rock. \(PI>6\), plastic packing rock. Under reservoir conditions, the water content in the pores is always greater than or equal to the water content in the solid phase. It is necessary to add a limiting condition to \(\:\dot{S}\) to make it equal to 0 when \({c_1}<{c_0}\). Therefore, \(\:\dot{S}\) is:

$$\dot {S}=\left\{ {\begin{array}{*{20}{c}} {\left( {{c_1} - {c_0}} \right){k_f}\exp \left( {\frac{{{k_f}^{2}}}{{{D_w}}}t} \right){\text{erfc}}\left( {\sqrt {\frac{{{k_f}^{2}}}{{{D_w}}}t} } \right),{c_1}>{c_0}} \\ {0,{c_1} \leqslant {c_0}} \end{array}} \right.$$
(21)

In Eq. (2), when f is the volume fraction of water in the pore at a certain moment (dynamic), \(\:{c}_{c}=\frac{{D}_{w}}{{\varphi\:}_{0}}\), \(\:{a}_{a}=\frac{u}{{\varphi\:}_{0}}\), \(\:{b}_{b}=\frac{\partial\:u}{{\varphi\:}_{0}\partial\:x}\), the expression of reservoir water content distribution is obtained. The numerical model of hydration expansion damage is obtained by combining the formula (9,20,21).

Water blocking effect damage

The water phase movement equation in the water-blocking damage reservoir is consistent with the water phase movement equation in the pore of the clay swelling model. The fractal method is used to describe the pore structure of reservoir rocks41. Assuming that the fractal dimension is Df (0 < Df<3), according to the fractal theory, the number of pores N with pore size greater than \(\:{r}_{pore}\) has the following power function relationship with \(\:{r}_{pore}\):

$$\:N(>{r}_{pore})={\int\:}_{{r}_{pore}}^{{r}_{poremax}}f\left({r}_{pore}\right)d{r}_{pore}=a{{r}_{pore}}^{-{D}_{f}}$$
(22)

The pore size distribution function \(\:f\left({r}_{pore}\right)\) satisfies:

$$\:f\left({r}_{pore}\right)=\frac{dN}{d{r}_{pore}}={a}^{{\prime\:}}{{r}_{pore}}^{-{D}_{f}-1}$$
(23)
$$\:{a}^{{\prime\:}}=-a{D}_{f}$$

For the capillary bundle model, the flow rate Q can be obtained by integrating the pore size distribution function, and the flow rate is expressed as:

$$Q=\frac{\pi\:\varDelta\:p}{8\mu\:l}{\int\:}_{{r}_{pore}}^{{r}_{poremax}}<n{{r}_{pore}}^{4}d{r}_{pore}$$
$$\:n=N{\int\:}_{{r}_{pore}}^{{r}_{poremax}}f\left({r}_{pore}\right)d{r}_{pore}$$
$$\:Q=b\left({r}_{poremax}^{3-{D}_{f}}-{{r}_{pore}}^{3-{D}_{f}}\right)$$
(24)
$$\:b=\frac{N\pi\:{a}^{{\prime\:}}\sigma\:{cos}\theta\:}{4\mu\:l(3-{D}_{f})}$$

The fractal expression of total pore volume is:

$$\:{\varPhi\:}_{max}=A{D}_{f}{\int\:}_{{r}_{poremin}}^{{r}_{poremax}}{{r}_{pore}}^{-1-{D}_{f}}d{r}_{pore}=\frac{A{D}_{f}}{3-{D}_{f}}{r}_{poremax}^{3-{D}_{f}}\left(1-{\left(\frac{{r}_{poremin}}{{r}_{poremax}}\right)}^{3-{D}_{f}}\right)$$
(25)

Higher invaded-fluid saturation will cause more-severe formation damage42. Assuming that the pore throat with a pore diameter greater than \(\:{r}_{pore}\) is completely occupied by the non-aqueous phase (i.e., assuming that the rock is hydrophilic), the non-aqueous phase saturation and the water phase saturation are:

$$\:{S}_{nw}=\frac{\varPhi\:}{{\varPhi\:}_{max}}=\frac{1-{\left({r}_{pore}/{r}_{poremax}\right)}^{3-{D}_{f}}}{1-{\left({r}_{poremin}/{r}_{poremax}\right)}^{3-{D}_{f}}}$$
$$\:{S}_{w}=1-{S}_{nw}=\frac{{\left({r}_{pore}/{r}_{poremax}\right)}^{3-{D}_{f}}-{\left({r}_{poremin}/{r}_{poremax}\right)}^{3-{D}_{f}}}{1-{\left({r}_{poremin}/{r}_{poremxa}\right)}^{3-{D}_{f}}}$$
(26)

From the linear Hagen-Poiseuille viscous flow, the permeability of the capillary bundle model can be expressed as:

$$\:K=\frac{\varphi\:}{8{\tau\:}^{2}}{\sum\:}_{i}{\xi\:}_{i}{r}_{i}^{2}$$
(27)

That means, the reservoir permeability channel is regarded as the accumulation of multiple capillary bundles.

Because of the continuity of pore throat size distribution, Eq. (26) can be changed to integral form and the expression of bulk density function can be written:

$$\:K=\frac{\varphi\:}{8{\tau\:}^{2}}\int\:{{r}_{pore}}^{2}d\xi\:$$
$$\:d\xi\:=\frac{(3-{D}_{f}){{r}_{pore}}^{2-{D}_{f}}}{{r}_{poremax}^{3-{D}_{f}}\left(1-{\left({r}_{poremin}/{r}_{poremax}\right)}^{3-{D}_{f}}\right)}d{r}_{pore}$$
(28)

According to the fractal expression of pore volume and density distribution expression, the permeability expression can be calculated.

$${K_d}({S_w})=\left( {\frac{{1 - {{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{5 - {D_f}}}}}{{1 - {{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{3 - {D_f}}}}}} \right) \cdot \frac{{1 - \left( {1 - {{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{3 - {D_f}}}} \right) \cdot {S_w} - {{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{3 - {D_f}}}}}{{1 - {{\left( {\left( {1 - {{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{3 - {D_f}}}} \right) \cdot {S_w}+{{\left( {\frac{{{r_{pore\hbox{min} }}}}{{{r_{pore\hbox{max} }}}}} \right)}^{3 - {D_f}}}} \right)}^{\left( {\frac{{5 - {D_f}}}{{3 - {D_f}}}} \right)}}}}$$
(29)

Combined with the expression and formula (9,29) of reservoir water content distribution in 2.1.4, the numerical model of water blocking damage can be obtained.

Emulsification blockage damage

In porous media, low interfacial tension and high mechanical shear force between oil and water phases are the main factors causing the formation of emulsions43,44. Therefore, the emulsification damage model must consider the effects of oil-water interfacial tension and flow velocity-induced mechanical shear on emulsification formation. This article mainly simulates the conditions for emulsion droplet formation and the formed emulsion droplet radius. According to the emulsification conditions given in the literature, the expression for the formed emulsion droplet radius can be derived from the theory of interfacial hydrodynamic instability:

$${r_o}(r,i\Delta t)=\left\{ \begin{gathered} {\left( {\frac{{9\pi \sigma \overline {r} }}{{2\left( {\frac{{\phi \mu }}{K}u(r,i\Delta t)+\rho g} \right)}}} \right)^{\frac{1}{3}}},\text{i}\text{f}{\left( {\frac{{9\pi \sigma \overline {r} }}{{2\left( {\frac{{\phi \mu }}{K}u(r,i\Delta t)+\rho g} \right)}}} \right)^{\frac{1}{3}}}<\overline {r} \hfill \\ 0,\text{i}\text{f}{\left( {\frac{{9\pi \sigma \overline {r} }}{{2\left( {\frac{{\phi \mu }}{K}u(r,i\Delta t)+\rho g} \right)}}} \right)^{\frac{1}{3}}} \geqslant \overline {r} \hfill \\ \end{gathered} \right.$$
(30)

After the discretization of the partial differential equation, the grid is formed, so the permeability does not exist in the form of matrix. The dimensionless permeability expression after the damage of emulsification plugging is:

$${K_d}(t=n\Delta t)={\left( {1 - \frac{{\sum\nolimits_{{i=0}}^{n} {\beta \left( {{r_o}(i\Delta t)} \right){r_o}(i\Delta t)} }}{{\overline {r} }}} \right)^{{m_K}}}$$
(31)

The pore size distribution function \(N(r)\) is approximately a log-normal function:

$$N(r,{\mu _s},{\sigma _s})=\frac{1}{{{\sigma _s}\sqrt {2\pi } }}\exp \left[ { - \frac{{{{\left( {\ln r - {\mu _s}} \right)}^2}}}{{2\sigma _{s}^{2}}}} \right]$$
(32)

The mean and standard deviation of the aperture in the formula are \(E(r)={e^{{\mu _s}+\frac{{\sigma _{s}^{2}}}{2}}}\)(that is \(\bar {r}\)), \(SD(r)={e^{{\mu _s}+\frac{{\sigma _{s}^{2}}}{2}}}\sqrt {{e^{\sigma _{s}^{2}}} - 1} =E(r)\sqrt {{e^{\sigma _{s}^{2}}} - 1}\). Then the formula of \({\mu _s},{\sigma _s}\)are:

$${\mu _s}=\ln \left[ {\frac{{\bar {r}}}{{\sqrt {1+{{\left( {SD(r)/\bar {r}} \right)}^2}} }}} \right],_{{}}^{{}}{\sigma _s}=\sqrt {\ln \left[ {1+{{\left( {\frac{{SD(r)}}{{\bar {r}}}} \right)}^2}} \right]}$$
(33)

Then the specific form of pore size distribution function is obtained.

For a droplet with a certain radius (\({r_o}\)), only the part of the pore size smaller than \({r_o}\), the distribution function will be blocked by the droplet. Therefore, the cumulative distribution of the blocking probability (\(\beta\)) from 0 to \({r_o}\) is:

$$\beta ({r_o})=\int_{0}^{{{r_o}}} {N(r;{\mu _s},{\sigma _s})\text{d}r} =\frac{1}{2}\text{e}\text{r}\text{f}\text{c}\left( { - \frac{{\ln {r_o} - {\mu _s}}}{{{\sigma _s}\sqrt 2 }}} \right)$$
(34)

Where erfc() is the famous residual error function: \(\text{e}\text{r}\text{f}\text{c}(x)=\frac{2}{{\sqrt \pi }}\int_{x}^{\infty } {{e^{ - {z^2}}}\text{d}\text{z}}\).

The Eqs. (9,30 ~ 34) can establish a numerical model of emulsification damage.

Wetting reversal damage

The following relationship can better describe the relationship between oil relative permeability and oil saturation2:

$${K_{ro}}={\alpha _0}+{\alpha _1}{S_o}+{\alpha _2}S_{o}^{2}+{\alpha _3}S_{o}^{3}+{\alpha _4}S_{o}^{4},_{{}}^{{}}\left( {0 \leqslant {K_{ro}} \leqslant 1} \right)$$
(35)

In Eq. (2), when f is oil saturation \(\:{S}_{o}\), \(\:{c}_{c}=\frac{{D}_{0}}{\varphi\:}\), \(\:{a}_{a}=\frac{u}{\varphi\:}\), \(\:{b}_{b}=\frac{\partial\:u}{\varphi\:\partial\:x}\), this equation is the convection diffusion law satisfied by oil saturation \(\:{S}_{o}\). Combined with Formula (9,35), a numerical model of wetting reversal damage can be established.

Stress sensitive damage

Fig. 2
figure 2

(a) Micro porous media composed by irregular solid particles clusters. (b) Radial micro deformation of solid particle clusters under stress. (c) Transverse micro deformation of solid particle clusters under stress.

Full size image

Microscopic skeletal particle clusters create a specific permeation pathway, with hydraulic bundles interspersed between these clusters. The hydraulic bundles exhibit a particular tortuosity that results in their actual length being longer than L (as depicted in Fig. 2a). Under the influence of effective stresses, the particle clusters contract (Fig. 2b), and the permeation pathway extends (Fig. 2c).

Based on the fractal geometry theory of rock, the number density of particle cluster diameter (\({\lambda _{pc}}\)), that is, the cumulative distribution of a certain diameter particle in unit volume (\({N_{pc}}\)), satisfies the relationship:

$${N_{pc}}({\lambda _{pc}} \leqslant L)={\left( {\frac{{{\lambda _{pc\hbox{max} }}}}{{{\lambda _{pc}}}}} \right)^{{D_{\text{c}\text{f}}}}}$$
(36)

For two-dimensional rock model, \(0<{D_{\text{c}\text{f}}}<2\); For three-dimensional rock model, \(0<{D_{\text{c}\text{f}}}<3\).

For the length of the seepage channel \({L_c}\) (i.e., the length of the particle cluster) in the cell, there are:

$${L_c}=\lambda _{{pc}}^{{1 - {D_{cT}}}}{L^{{D_{cT}}}}$$
(37)

These are the relationships within a spatial micro-element. There are many micro-elements and many particle clusters in a region. Under the action of effective stress \({\sigma _e}\), there are:

$${\lambda _{pc}}=F(1 - {\sigma _e}/({\nu _{po}}E)){\lambda _{pc0}};_{{}}^{{}}{L_c}=(1 - {\sigma _e}/({\nu _{po}}E)){L_{c0}}$$
(38)

Using the flow formula of the hydraulic pipe and the infinitesimal Darcy’s law45, the expression of reservoir permeability under the influence of effective stress can be obtained:

$$K=\frac{{\left( {2 - {D_{\text{c}\text{f}}}} \right)\pi G{F^4}L_{0}^{{1 - {D_{cT}}}}\lambda _{{pc\hbox{max} 0}}^{{1+{D_{cT}}}}{{(1+{\sigma _e}/E)}^2}}}{{32\sqrt 3 \left( {3+{D_{cT}} - {D_{\text{c}\text{f}}}} \right)\left[ {1 - {\sigma _e}/({\nu _{po}}E)} \right]\left[ {1 - {{\left( {{\lambda _{pc\hbox{min} 0}}/{\lambda _{pc\hbox{max} 0}}} \right)}^{2 - {D_{\text{c}\text{f}}}}}} \right]}}$$
(39)

Therefore, the dimensionless permeability of reservoir damage caused by stress sensitivity (\({K_d}\)) is :

$${K_d}=\frac{{K({\sigma _e})}}{{K({\sigma _e}={\sigma _0})}}=\frac{{\left[ {1 - {\sigma _e}/({\nu _{po}}E)} \right]{{(1+{\sigma _0}/E)}^2}}}{{{{(1+{\sigma _e}/E)}^2}\left[ {1 - {\sigma _0}/({\nu _{po}}E)} \right]}}$$
(40)

In formula (9), the relationship between effective stress, overburden pressure and pore pressure is:\({\sigma _e}={P_o} - P\). Detailed calculation of pore pressure can be found in references15,46,47. The numerical model of stress sensitive damage is obtained by combining Eqs. (9,39).

Sand production damage

According to the sand production mechanism, the key parameters causing reservoir damage due to sand production are the threshold flow velocity and critical flow velocity. Therefore, the core of the quantitative simulation of the spatiotemporal evolution of sand production damage in reservoirs is to solve for the threshold flow velocity and critical flow velocity. The threshold flow velocity is the flow velocity at which sand particles start to move. When the fluid flow velocity is greater than the threshold flow velocity, some sand particles begin to move, and the sand produced at this time is adhered sand. If the fluid flow velocity exceeds the critical flow velocity, the rock skeleton is sheared and damaged, and a large amount of sand is produced, which includes both adhered sand and skeleton sand. The critical flow velocity is related to the critical production pressure difference, so first, a numerical simulation model for the critical production pressure difference is established.

After the formation is drilled, the stress distribution around the wellbore changes, and stress concentration occurs at the wellbore wall. For this, the following assumptions are made before establishing the model: (1) A vertical well penetrates a permeable horizontal oil layer; (2) The oil reservoir rock is an isotropic homogeneous elastic body; (3) The rock voids are filled with liquid; (4) Structural stress is ignored; (5) The rock pores have good connectivity.

According to Mohr-Coulomb rock failure criterion48:

$$\left| \tau \right|={S_{in}}+{\sigma _{zs}}\tan \theta$$
(41)

When the stress circle radius \({r_m} \geqslant d\), rock shear failure occurs, so the critical state is \({r_m}=d\). When \({r_m}=d\), the bottom hole flowing pressure is the critical bottom hole flowing pressure \({P_{er}}\):

$${P_{{\text{er}}}}=\frac{{\sin \theta \left[ {\frac{{2{\nu _{po}}}}{{1 - {\nu _{po}}}}{\sigma _{zo}}+\frac{{1 - 2{\nu _{po}}}}{{1 - {\nu _{po}}}}BP} \right]+2{S_{in}}\cos \theta - \frac{{2{\nu _{po}}}}{{1 - {\nu _{po}}}}{\sigma _{zo}} - \frac{{1 - 2{\nu _{po}}}}{{1 - {\nu _{po}}}}BP}}{{\frac{{1 - 2{\nu _{po}}}}{{1 - {\nu _{po}}}}B+\frac{{\sin \theta }}{{1 - {\nu _{po}}}}B - 2}}$$
(42)

The above is the theoretical calculation formula of critical bottom hole flowing pressure. The critical production pressure difference is: \(\Delta {P_T}=P - {P_{{\text{er}}}}\). After obtaining the critical production pressure difference \(\Delta {P_T}\), according to the Jubi formula, the critical production expression is:

$${q_T}=\frac{{2\pi {K_o}H\Delta {P_T}}}{{{\mu _o}{B_o}(\ln \frac{{{r_e}}}{{{r_w}}})}}$$
(43)

In Eq. (2), when f is the flow particle concentration C, \(\:{a}_{a}=\frac{v}{\tau\:}\left[1-\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right)k\alpha\:\tau\:\right],\)\(\:{b}_{b}f=\frac{kv}{\tau\:}\left(1-\frac{{\rho\:}_{p}}{{\rho\:}_{L}}C\right)+\lambda\:C,\:{c}_{c}=\alpha\:v\), combined with Equation (8,9,30,31), the numerical model of sand damage can be obtained. Sand production is similar to particle migration. The space-time numerical simulation of sand production damage also has intrinsic randomness and complexity. The difference is that the criteria for ' \({q_s}\) is 0 ' are different. The difference here is reflected in the value of λ:

$$\lambda =\left\{ \begin{gathered} 0,(q<{q_T}) \hfill \\ \frac{{vk\phi }}{\tau },(q>{q_T}) \hfill \\ \end{gathered} \right.$$
(44)

10 polymer adsorption damage

The adsorption of polymers on the surface of porous media is mainly divided into two adsorption modes : layer adsorption and bridging adsorption. The polymer density functions of layer adsorption and bridging adsorption are as follows:

$$\begin{gathered} \frac{{\partial {\Gamma _l}}}{{\partial t}}=\kappa _{a}^{l}\left( {\Gamma _{l}^{\infty } - {\Gamma _l}} \right){C_{poly}} - \kappa _{d}^{l}{\Gamma _l} \hfill \\ \frac{{\partial {\Gamma _b}}}{{\partial t}}=\kappa _{a}^{b}\left( {\Gamma _{b}^{\infty } - {\Gamma _b}} \right){C_{poly}} - \kappa _{d}^{b}{\Gamma _b} \hfill \\ \end{gathered}$$
(45)

Similarly, according to the theory of material balance, a one-dimensional model of polymer adsorption49,50 can be established as shown in formula (34).

$$\begin{gathered} \frac{{\partial {C_{poly}}}}{{\partial t}}+\frac{u}{\phi }\frac{{\partial {C_{poly}}}}{{\partial x}}= - {\rho _s}{A_s}\frac{{1 - \phi }}{\phi }\frac{{\partial {\Gamma _l}}}{{\partial t}} - {\rho _s}{A_s}\frac{{1 - \phi }}{\phi }\frac{{\partial {\Gamma _b}}}{{\partial t}}+{D_{poly}}\frac{{{\partial ^2}{C_{poly}}}}{{\partial {x^2}}} \hfill \\ \frac{{\partial {P_{poly}}}}{{\partial t}}+\frac{u}{\phi }\frac{{\partial {P_{poly}}}}{{\partial x}}= - {\chi _q}{\rho _s}{A_s}\frac{{1 - \phi }}{{\phi {C_{poly0}}}}\frac{{\partial {\Gamma _l}}}{{\partial t}} - {\rho _s}{A_s}\frac{{1 - \phi }}{{\phi {C_{poly0}}}}\frac{{\partial {\Gamma _b}}}{{\partial t}} \hfill \\ \end{gathered}$$
(46)

Organic scale damage

During the development of oil fields, changes in reservoir temperature and pressure disrupt the equilibrium state between crude oil, reservoir rock, and reservoir fluid in the oil reservoir. This causes some light components of the crude oil to spill over and some heavy components to precipitate, forming organic scales that block oil and gas flow channels and cause serious damage to the reservoir51.

Generally speaking, organic scales such as asphaltene exist in crude oil in two forms: partially dissolved and partially colloidal suspension. The precipitation of asphaltenes and other organic substances is related to their solubility. The following describes the numerical simulation of organic scale damage based on the changes in the solubility of asphaltenes and other organic substances.

(1) Single component solubility parameter δ.

The solubility parameter of single component is defined as ‘cohesive energy per unit volume’49,52,53. The specific equation is as follows:

$$\delta =\sqrt {\left( {\frac{{\Delta H^\prime - R \cdot T}}{{{V_{{\text{liq}}}}}}} \right)}$$
(47)

(2) Dissolution ability ( the proportion of soluble in oil phase in all organic scale such as asphaltene ).

The following equation can describe the solubility:

$${\left( {{\phi _a}} \right)_{\hbox{max} }}=\exp \left\{ {\frac{{{V_a}}}{{{V_L}}}\left[ {1 - \frac{{{V_L}}}{{{V_a}}} - \frac{{{V_L}}}{{RT}}{{\left( {{\delta _a} - {\delta _L}} \right)}^2}} \right]} \right\}$$
(48)

The behavior of the oil phase molecular cohesive energy per unit volume, that is, the oil phase solubility coefficient \(\:{\delta\:}_{L}\), is the most complex. It is necessary to solve each component separately and then use \(\:{\delta\:}_{L}={\sum\:}_{i}{\varphi\:}_{i}{\delta\:}_{i}\) to calculate. For component i, there is:

$$\:{\delta\:}_{i}={\left(\frac{{\varDelta\:{H}^{{\prime\:}}}_{i}{|}_{T}-R \cdot\:T}{{V}_{L}}\right)}^{0.5}$$
(49)
$$\:{\varDelta\:{H}^{{\prime\:}}}_{i}{|}_{{T}_{bi}}=1.014\times\:{T}_{bi}\times\:\left(8.75+4.571ln\left({T}_{bi}\right)\right)$$
(50)
$$\:{\varDelta\:{H}^{{\prime\:}}}_{i}{|}_{T}={\varDelta\:{H}^{{\prime\:}}}_{i}{|}_{{T}_{bi}}\times\:{\left(\frac{{T}_{ci}-T}{{T}_{ci}-{T}_{bi}}\right)}^{0.38}$$
(51)

he data of \(\:{\varphi\:}_{i}\) can be obtained from the analysis data of oil properties. In general, the above expression is only used to establish the relationship between different pressures and temperatures, rather than directly calculating with the above formula. This requires the \(\:{\delta\:}_{i}^{ref}\) under a certain reference pressure (\(\:{P}_{ref}\)) and reference temperature (\(\:{T}_{ref}\)) measured experimentally, and then calculate the \(\:{\delta\:}_{i}\left(P,T\right)\) value at the required temperature and pressure.

$$\:{\delta\:}_{i}\left(P,T\right)={\delta\:}_{i}^{ref}\cdot\:\sqrt{\frac{\frac{1.014\times\:{T}_{bi}\times\:\left(8.75+4.571ln\left({T}_{bi}\right)\right)\times\:{\left(\frac{{T}_{ci}-T}{{T}_{ci}-{T}_{bi}}\right)}^{0.38}-R\cdot\:T}{{V}_{L}\left(P\right)}}{\frac{{\varDelta\:{H}^{{\prime\:}}}_{i}{|}_{{T}_{ref}}-R\cdot\:{T}_{ref}}{{V}_{L}\left({P}_{ref}\right)}}}$$
(52)

(3) Flow model.

Fig. 3
figure 3

Reservoir cylindrical model.

Full size image

The reservoir is divided into cylindrical shell combinations as shown in Fig. 3. For the i-th molar fluid, the molar number of pitch-containing particles is54:

$$\:{\eta\:}_{i}={\int\:}_{0}^{{t}_{0}}\left({\varphi\:}_{a}-{\left({\varphi\:}_{a}\right)}_{max}\right){f}_{trap}\frac{{\rho\:}_{o}{Q}_{o}\left(t\right)}{{MW}_{o}}dt$$
(53)

Among them, \(\:{f}_{trap}\) is a function related to the particle size distribution of asphalt particles and the average pore size \(\:{r}_{p}\):

$$\:{f}_{trap}={\int\:}_{{r}_{p}}^{{\infty\:}}f\left({r}_{a}\right){dr}_{a}$$
(54)

\(\:f\left({r}_{a}\right)\) is the size distribution of asphalt particles (generally 10−1~10 μm).

Pressure in the ith shell:

$$\:{P}_{i}={P}_{e}-1.8421\times\:\frac{{Q}_{o}\mu\:}{{K}_{0}H}\left(\frac{{r}_{e}}{{r}_{i}}\right)$$
(55)

The change of porosity in the ith shell:

$$\:\varDelta\:{\varphi\:}_{i}=\frac{{\eta\:}_{i}{V}_{a}\times\:{10}^{-6}}{2\pi\:h{\varphi\:}_{0}{r}_{i}dr}$$
(56)

The dimensionless permeability calculation of reservoir damage caused by organic scale blockage:

$$\:{K}_{d}\left({r}_{i}\right)={\left(1-\varDelta\:{\varvec{\varphi\:}}_{\varvec{i}}\right)}^{2}$$
(57)

The skin factor of the reservoir caused by organic scale blockage:

$$\:S\left({r}_{i}\right)=\left(\frac{1}{{K}_{d}\left({r}_{i}\right)}-1\right)ln\left(\frac{{r}_{m}}{{r}_{w}}\right)$$
(58)

The initial condition is the state before the damage, that is \({K_d}(r,t=0)=1\). For a certain time point, there are multiple skin factor values in different spaces. The mean value of the skin factor at the same time point on the spatial grid node is the dimensionless permeability and skin factor of the well to be diagnosed, that is, the numerical model of organic fouling damage is obtained.

Bacterial damage

The impact of reservoir temperature on the maximum bacterial growth rate is significant. As a result, scholars at home and abroad have established general models of bacterial content under varying temperature conditions, as well as square root models containing parameters such as activation energy and frequency factors to describe the maximum bacterial growth rate. However, both types of models have poor applicability in the field. Therefore, under the conditions of sufficient water injection sources and bacterial nutrient sources, this paper uses an empirical formula for the maximum bacterial growth rate with temperature as the main variable55,56,57,58,59 to simulate the spatial distribution of the maximum bacterial growth rate in the reservoir.

$${g_{\hbox{max} }}={[{b_1}(T - {T_{\hbox{min} }})]^2} * \left\{ {1 - \exp [{c_1}(T - {T_{\hbox{max} }})]} \right\}$$
(59)

The main factors affecting the amount of bacteria adhering to the rock surface are the bacterial adsorption and desorption rates, and the formation of biofilm when bacteria adhere to the reservoir rock surface. Therefore, the equation for the amount of bacteria adhering to the rock surface is:

$${C_{{\text{deposition}}}}={k_{{\text{clogging}}}}{C_{{\text{bacteriatran}}}} - {k_{{\text{declogging}}}}{V_{{\text{bacteriatran}}}}$$
(60)

The amount of bacteria on the rock surface is determined by two major factors : bacterial growth, net increase in decay and attachment. The bacteria attached to the rock surface will form a biofilm to reduce reservoir porosity. The equation is expressed as\(\:({V}_{\text{bacteriatran}}=\rho\:{V}_{\text{bacteria}})\)60:

$$\frac{{\partial {V_{{\text{bacteriatran}}}}}}{{\partial t}}=({g_{{\text{actual}}}} - {k_{{\text{decay}}}}){V_{{\text{bacteriatran}}}}+{C_{{\text{deposition}}}}$$
(61)

According to the Gauchani equation, the permeability of the reservoir is related to the specific surface area of the reservoir porosity, tortuosity and rock surface volume. The equation is expressed as:

$$K=\frac{{{\theta ^3}}}{{2{\tau ^2}{S_{vp}}^{2}}}$$
(62)

The migration of bacteria and nutrients in reservoirs is similar to that of porous media particles, which can be described by the same method. Therefore, in Eq. (2), when f is the reservoir temperature distribution T, \(\:{a}_{a}=u\), \(\:{b}_{b}=0,\:{{c}_{c}=D}_{\text{c}\text{o}\text{n}}+{D}_{\text{d}\text{i}\text{s}}\), this equation is the control equation of reservoir temperature distribution; when f is the apparent concentration of bacteria in the reservoir fluid \(\:{C}_{\text{b}\text{a}\text{c}\text{t}\text{e}\text{r}\text{i}\text{a}\text{n}}\), \(\:{a}_{a}=u\), \(\:{b}_{b}={k}_{\text{d}\text{e}\text{c}\text{a}\text{y}}-{g}_{\text{a}\text{c}\text{t}\text{u}\text{a}\text{l}}\), \(\:{{c}_{c}=D}_{\text{s}\text{u}\text{m}}\), this equation is the bacterial concentration distribution equation in the reservoir; when f is the apparent concentration of bacteria in reservoir fluid(\({C_{{\text{nutran}}}}\)), \(\:{a}_{a}=u\), \(\:{b}_{b}=0\), \(\:{{c}_{c}=D}_{\text{n}\text{u}\text{s}\text{u}\text{m}}\), this equation is the distribution equation of nutrient concentration in reservoir. Combined with the formula (59 ~ 62), the numerical model of bacterial damage was obtained.

Establishment of quantitative diagnosis technology for space-time evolution of reservoir damage

In this paper, the skin factor S is used to quantitatively characterize the degree of reservoir damage, and serves as one of the main parameters for quantitatively simulating the spatiotemporal evolution of reservoir damage, as shown in Eq. (63). In actual situations, the damage condition of the well to be diagnosed is often complex, involving multiple types of damage. The degree of reservoir damage \(\:{S}_{it}\) caused by each type at a certain moment t can be obtained according to the numerical models described earlier. Ignoring the mutual influence of various types, the total damage degree \(\:{S}_{\text{a}\text{t}}\) at a certain moment t is equal to the sum of the damages \(\:{S}_{it}\) caused by each type, as shown in Eq. (64). The proportion of each type of damage in the total damage degree, i.e., the contribution rate δit, is calculated using Eq. (65).

$$\:S=\left(\frac{1}{\overline{{K}_{d}}}-1\right){ln}\left(\frac{{r}_{sw}}{{r}_{w}}\right)$$
(63)
$$\:{S}_{\text{a}\text{t}}=\sum\:_{i=1}^{n}{S}_{it}$$
(64)
$$\:{\delta\:}_{it}=\frac{{S}_{it}}{{S}_{\text{a}\text{t}}}\times\:100\%$$
(65)

Of course, the calculated \(\:{S}_{it}\), \(\:{S}_{\text{a}\text{t}}\), and \(\:{\delta\:}_{it}\) at different times are different and change with time. At the same time, due to the non-homogeneity of reservoir properties in different well sections and the variability of internal and external factors of reservoir damage, the calculated \(\:{S}_{it}\), \(\:{S}_{\text{a}\text{t}}\), and \(\:{\delta\:}_{it}\) for different well sections or reservoir depths are also different. In other words, the degree of reservoir damage has a distribution in both time t and space (X, Y, and Z-axis directions), which can be obtained using the numerical simulation methods described in Sect. "Establishment of quantitative simulation control equation for spatial-temporal evolution of common types of damage" of this paper. However, the spatiotemporal evolution numerical simulation diagnosis of reservoir damage types and degrees in time and space not only involves large calculations and multiple iterations but also has a complex calculation process. Therefore, the “Reservoir Damage Types and Degrees Spatiotemporal Evolution Numerical Simulation System” software is used for quantitative simulation. Figure 4 shows the quantitative simulation results of the spatiotemporal evolution of reservoir damage during the drilling process of a vertical well using this software. It can be seen that the spatiotemporal evolution numerical simulation has strong intuitiveness and high visualization (the degree of damage is represented by grayscale values in the figure, the darker the color, the greater the damage, the lighter the color, the smaller the damage).

Fig. 4
figure 4

Distribution of permeability damage near wellbore zone at different times.

Full size image

Theory and calculation

Difference solution of control equation

For Eq. (2), backward difference is used for time and central difference is used for space. The above equation can have the following difference scheme:

$$\:\frac{{f}_{i}^{n}-{f}_{i}^{n-1}}{\varDelta\:t}+{a}_{i}\frac{{f}_{i+1}^{n}-{f}_{i-1}^{n}}{2\varDelta\:x}+{b}_{i}{f}_{i}^{n}={c}_{i}\frac{{f}_{i+1}^{n}-2{f}_{i}^{n}+{f}_{i-1}^{n}}{\varDelta\:{x}^{2}},\:i=\text{1,2},3\dots\:{N}_{i},\:{N}_{i}=\frac{{x}_{max}}{\varDelta\:x}n=\text{1,2},3\dots\:,\:t=n\varDelta\:t$$
(66)

\(\:{N}_{i}\) is the number of discrete space points.

Solve the interval \(\:x\in\:(0,{x}_{max})\), \(\:\varDelta\:x\) and \(\:\varDelta\:t\) are the space and time step. Both initial conditions and boundary conditions are under consideration, \(\:{f}_{i}^{n}{|}_{n=0}={f}_{i}^{0},i=\text{1,2},3\dots\:,{N}_{i}\), \(\:{f}_{i}^{n}{|}_{i=1}={f}_{0},n=\text{1,2},3\dots\:\), \(\:\frac{{f}_{i+1}^{n}-{f}_{i-1}^{n}}{2\varDelta\:x}{|}_{i=N}=0,n=\text{1,2},3\dots\:\) ( central difference and construct a virtual grid \(\:i+1\) )

Firstly, for \(\:i=\text{2,3},\dots\:,{N}_{i}-1\), the above difference scheme is rewritten as:

$$\left(\frac{1}{\varDelta\:t}+{b}_{i}+\frac{2{c}_{i}}{\varDelta\:{x}^{2}}\right){f}_{i}^{n}+\left(\frac{{a}_{i}}{2\varDelta\:x}-\frac{{c}_{i}}{\varDelta\:{x}^{2}}\right){f}_{i+1}^{n}-\left(\frac{{a}_{i}}{2\varDelta\:x}+\frac{{c}_{i}}{\varDelta\:{x}^{2}}\right)=\frac{1}{\varDelta\:t}{f}_{i}^{n-1}\Rightarrow$$
$${f}_{i}^{n}+\frac{\frac{{a}_{i}}{2\varDelta\:x}-\frac{{c}_{i}}{\varDelta\:{x}^{2}}}{\frac{1}{\varDelta\:t}+{b}_{i}+\frac{2{c}_{i}}{\varDelta\:{x}^{2}}}{f}_{i+1}^{n}-\frac{\frac{{a}_{i}}{2\varDelta\:x}+\frac{{c}_{i}}{\varDelta\:{x}^{2}}}{\frac{1}{\varDelta\:t}+{b}_{i}+\frac{2{c}_{i}}{\varDelta\:{x}^{2}}}{f}_{i-1}^{n}=\frac{\frac{1}{\varDelta\:t}}{\frac{1}{\varDelta\:t}+{b}_{i}+\frac{2{c}_{i}}{\varDelta\:{x}^{2}}}{f}_{i}^{n-1}\Rightarrow$$
(67)
$${f}_{i}^{n}+\frac{\frac{{a}_{i}}{2}-\frac{{c}_{i}}{\varDelta\:x}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}}{f}_{i+1}^{n}-\frac{\frac{{a}_{i}}{2}+\frac{{c}_{i}}{\varDelta\:x}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}}{f}_{i-1}^{n}=\frac{\frac{\varDelta\:x}{\varDelta\:t}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}}{f}_{i}^{n-1}\Rightarrow$$
$${f}_{i}^{n}=A{1}_{i}{f}_{i-1}^{n}+A{2}_{i}{f}_{i+1}^{n}+A{3}_{i}{f}_{i}^{n-1}(i=\text{2,3}\dots\:,{N}_{i}-1)$$

where \(\:A{1}_{i}\), \(\:A{2}_{i}\) and \(\:A{3}_{i}\) are

$$\:\frac{\frac{{a}_{i}}{2}+\frac{{c}_{i}}{\varDelta\:x}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}},\:\frac{\frac{{c}_{i}}{\varDelta\:x}-\frac{{a}_{i}}{2}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}},\:\frac{\frac{\varDelta\:x}{\varDelta\:t}}{\frac{\varDelta\:x}{\varDelta\:t}+{b}_{i}\varDelta\:x+\frac{2{c}_{i}}{\varDelta\:x}}$$
(68)

Thus the iterative relationship is obtained:

$$\:{f}_{i}^{n}=A{1}_{i}{f}_{i-1}^{n}+A{2}_{i}{f}_{i+1}^{n}+A{3}_{i}{f}_{i}^{n-1}\:(i=\text{2,3}\dots\:,N-1)$$
(69)

The value of field f can be obtained by iterative calculation using initial and boundary conditions.

Boundary condition difference solution

Iterative relation (69) is suitable for non-boundary grids. For \(\:i=1\) (at the shaft wall), because the point-centered mesh is used and the Dirichlet boundary condition is used, it is directly: \(\:{f}_{1}^{n}={f}_{0}\) (constant), \(\:i=1\).

For \(\:i=N\), it is Neumann boundary condition. Adding a virtual mesh \(\:i={N}_{i}+1\), from \(\:\frac{{f}_{i+1}^{n}-{f}_{i-1}^{n}}{2\varDelta\:x}{|}_{i=N}=0,n=\text{1,2},3\dots\:\) it can be known that \(\:{f}_{i+1}^{n}={f}_{i-1}^{n}\). Substituting Eq. (62), it can be seen that:

$$\:{f}_{i}^{n}=\left(A{1}_{i}+A{2}_{i}\right){f}_{i-1}^{n}+A{3}_{i}{f}_{i}^{n-1}\:\left(i=N\right)$$
(70)

Cylindrical coordinate conversion

The space-time variation of field function f can be solved by difference and iteration. Since the numerical model in Sect. 2.1 in this paper is established for the reservoir near the wellbore, the cylindrical coordinate system r is needed to solve the distribution of a physical quantity f around the well. Therefore, the equation \(\:\frac{\partial\:f}{\partial\:t}+{a}_{a}\frac{\partial\:f}{\partial\:x}+{b}_{b}f={c}_{c}\frac{{\partial\:}^{2}f}{\partial\:{x}^{2}}\) needs to be transformed into \(\:\frac{\partial\:f}{\partial\:t}+{a}_{3}\frac{\partial\:f}{\partial\:r}+{b}_{3}f={c}_{3}\frac{1}{r}\frac{\partial\:}{\partial\:r}\left(r\frac{\partial\:}{\partial\:r}f\right)\). This form is not conducive to equidistant difference, and coordinate transformation can be introduced: \(\:r={r}_{w}{e}^{{x}^{{\prime\:}}}\), where \(\:{r}_{w}\) is the wellbore radius and \(\:{x}^{{\prime\:}}\) is a dimensionless spatial coordinate. Substituting this transformation into the general equation, we can get the equation about \(\:{x}^{{\prime\:}}\):

$$\:\frac{\partial\:f}{\partial\:t}+\frac{a}{{r}_{w}}\cdot\:{e}^{-{x}^{{\prime\:}}}\frac{\partial\:f}{\partial\:{x}^{{\prime\:}}}+bf=\frac{c}{{r}_{w}^{2}}{e}^{-{{x}^{{\prime\:}}}^{2}}\frac{{\partial\:}^{2}f}{\partial\:{{x}^{{\prime\:}}}^{2}}$$
(71)

If \(\:\frac{a}{{r}_{w}}\cdot\:{e}^{-{x}^{{\prime\:}}}\) and \(\:\frac{c}{{r}_{w}^{2}}{e}^{-{{x}^{{\prime\:}}}^{2}}\) are taken as new coefficients, then the above equation is essentially the same as \(\:\frac{\partial\:f}{\partial\:t}+{a}_{a}\frac{\partial\:f}{\partial\:x}+{b}_{b}f={c}_{c}\frac{{\partial\:}^{2}f}{\partial\:{x}^{2}}\). Therefore, it is possible to carry out isometric difference in \(\:{x}^{{\prime\:}}\) coordinates and follow the above iterative format. After calculating the value of f, we can get \(\:f(r,t)\) by mapping the space coordinates from \(\:{x}^{{\prime\:}}\) to r.

The simulated spatial range is 0.1–5 m, and 0.1 m is the simulated wellbore radius. The time step is 0.02 d, and the space step is variable. The grid type is structured grid. After exponential transformation of radial coordinates, the closer to the wellbore centerline, the denser the grid is, and the farther away from the wellbore, the sparser the grid is. There are 100 non-uniform grids in the range of 0.1–5 m. The software for numerical simulation is MATLAB.

Results & discussion

Laboratory verification

The mathematical model for diagnosing six reservoir damage mechanisms during the production stage (stress-sensitive damage, particle migration damage, wettability reversal damage, organic scale damage, emulsification damage, sanding damage) established in this paper was experimentally validated using classical core flow experiments. The experiments were conducted using the “High-Temperature High-Pressure Core Dynamic Damage Evaluation System, Model JDMH-II,” which simulates underground high-temperature and high-pressure conditions, displaces the core, and directly displays instantaneous permeability data, as shown in Fig. 5.

Fig. 5
figure 5

JDMH-II dynamic core damage evaluation system.

Full size image

In actual reservoir environments, multiple types of reservoir damage generally occur simultaneously. By sequentially or simultaneously simulating the presence of several damage types within the same core and measuring the permeability changes of the core under different types of damage, the contribution of each damage type to the permeability damage of the core and the total permeability damage ratio for that core can be determined. By inputting relevant core and experimental parameters into the model, the mathematical model simulates the permeability damage ratio, and a comparison between experimental and simulated results is performed to validate the accuracy of the mathematical model.

To prevent incidental errors from affecting the experimental results, two core samples were selected for concurrent experimental validation. The fundamental parameters of the experimental core samples are provided in Table 1.

Table 1 Basic parameters of core #1 and #2.
Full size table

The core samples underwent a sequence of damage, including wettability alteration damage, particle transport damage, sand production damage, stress sensitivity damage, organic scale damage, and emulsification damage. Permeability changes before and after each damage type were recorded for core samples #1 and #2, and the resulting permeability damage ratios were calculated, as shown in Table 2. In Table 2, K0 refers to the initial permeability, and Kd1-Kd6 represents the measured permeability of the core after each damage.

Table 2 Comprehensive formation damage result of cores.
Full size table

Based on the experimental results, both core samples, after undergoing simulation of the six damage types, experienced a significant decrease in permeability. Permeability decreased from the initial values of 80.31mD and 75.07mD to 15.30mD and 2.23mD, resulting in permeability damage ratios of 80.95% and 93.03%, respectively. Among the damage types, wettability alteration damage and emulsification damage were the most severe, followed by sand production damage, particle transport damage, and organic scale damage. Stress sensitivity damage contributed the least to the permeability damage ratios.

In the model, the basic parameters of the two core samples and experimental parameters were input, and numerical simulation results for dynamic permeability damage ratios and damage factor proportions were obtained, as illustrated in Figs. 6 and 7.

Fig. 6
figure 6

Dynamic permeability damage rate and damage factor ratio of core #1.

Full size image
Fig. 7
figure 7

Dynamic permeability damage rate and damage factor ratio of core #1.

Full size image

According to the simulation results from the reservoir damage mathematical model, the comprehensive permeability damage ratios for the two core samples were approximately 78.38% and 87.01%, respectively. When compared to the experimental results, the accuracy of the model was found to be approximately 90% or higher. Therefore, the model has passed the laboratory evaluation.

Applications in the field

Using the relevant parameters provided by the field, calculate the skin factor of the well to be diagnosed at several times later, and compare it with the measured skin factor to verify the correctness and accuracy of the spatiotemporal evolution numerical simulation technology established in this paper. Since the measured skin factor is the sum of the true skin factor and the pseudo skin factor, it is necessary to decompose the measured skin factor and compare the true skin factor obtained after decomposition with the calculated skin factor to judge the correctness and accuracy of the numerical simulation technology.

The numerical models studied in this paper are all for sandstone reservoirs. Now, let’s take the example of well G01s1 in Suizhong 36 − 1 oilfield whose production method is perforation completion. The relevant parameters and measured skin factors of this well are shown in Table 3. Suizhong 36 − 1 oilfield is a sandstone oil and gas reservoir. The surface crude oil has the characteristics of high density and high viscosity, which belongs to heavy heavy oil. The original density of the ground is 0.941 ~ 0.997 g / cm3, and the viscosity of the ground crude oil is 41.0 ~ 7787.4mPa·s, with an average of 1111.5mPa·s. The viscosity of formation crude oil is 37.4 ~ 154.7mPa·s, the saturation pressure is 10.24 ~ 13.33 MPa, the original dissolved gas-oil ratio is 31m3 / m3, and the volume coefficient is 1.090. The total salinity of formation water is 4481 ~ 7154 mg / L, and the water type is NaHCO3.

Table 3 Basic parameters and skin factor of G01s1 well in Suizhong 36 − 1 oilfield.
Full size table

By inputting the parameters in Table 1 into the software, the changes in the damage zone radius \(\:{r}_{\text{i}\text{t}}\) and damage degree \(\:{S}_{\text{i}\text{t}}\) caused by 6 types of damage, namely “sand production, stress sensitivity, emulsification, particle migration, wettability reversal, and organic scaling” over time (the damage degree caused by other types of damage is extremely slight and negligible) are obtained, as well as the distribution of the total damage zone radius \(\:{r}_{\text{a}\text{t}}\) and total damage degree \(\:{S}_{\text{a}\text{t}}\) after 40 days, as shown in Fig. 8. The evolution of the proportions of each damage type in the total damage over time is shown in Fig. 9.

Fig. 8
figure 8

Radius of Damage Zones (40 Days) and Distribution of Damage Degree with Time.

Full size image
Fig. 9
figure 9

The proportion of each damage factor evolved over time.

Full size image

From the simulation calculation results, it can be seen that on the 40th day, the calculated total skin factor is 35.65, the measured skin factor is 36.5, and the pseudo skin factor caused by well deviation and partial perforation is 0.49, with a relative error of only 0.99% and a very high degree of agreement. From the damage factor analysis, emulsification blockage, particle migration, organic scaling, and sand production are the main types of reservoir damage, while wettability reversal and stress sensitivity damage are small. Of course, according to the needs, the spatiotemporal evolution of reservoir damage at different depth points and directions can also be simulated and calculated.

In addition, the spatiotemporal evolution numerical simulation model established in this paper has been extensively verified and applied in different types of oilfields in China, with an average accuracy rate of more than 95%, fully demonstrating the reliability, accuracy, and practicality of the model. The simulation results comparison of some wells is shown in Table 4.

Table 4 Comparison of simulated and measured skin factor of some Wells.
Full size table

Conclusions

(1) Numerical models for 12 common reservoir damage types during drilling, well workover, oil production, water injection, and polymer injection stages have been established. A spatiotemporal evolution numerical model for reservoir damage types and degrees has been created in simple superposition way, which enables quantitative diagnosis and prediction of reservoir damage. For wells with existing reservoir damage, simulating the spatiotemporal evolution of each damage type and the overall damage degree provides the sensitivity of each damage type to the total damage and their proportion in the total damage degree. This is of significant importance for optimizing plugging measures and improving reservoir numerical simulation accuracy. For wells without reservoir damage, using petrophysical parameters and engineering parameters to be implemented allows quantitative prediction of reservoir damage and spatiotemporal extrapolation of damage patterns, providing scientific guidance for preventing or avoiding reservoir damage, formulating reservoir development plans, and subsequent production increase measures.

(2) Field verification shows that the spatiotemporal evolution numerical simulation technology for reservoir damage developed in this paper, based on the material balance theory and diffusion law during the reservoir damage process, is highly accurate, practical, and operable when used for the spatiotemporal evolution numerical simulation of reservoir damage types and degrees throughout the entire process of oil and gas exploration and development.

(3) The causes of reservoir damage are not only diverse and dynamically changing but also interconnected and mutually constrained. This affects the overall damage degree and spatiotemporal evolution pattern of oil and water wells. In the future, research on the mutual influence and constraint of various damage types should be strengthened to further improve and enhance the accuracy of reservoir damage quantitative simulation.

(4) The numerical model of reservoir damage established in this study has only been validated in sandstone porous reservoirs so far. Further improvement and optimization are needed for the diagnosis of damage in fractured reservoirs. Additionally, the coupling relationships between different types of damage considered in this study are relatively simple. In the future, considering fluid-solid coupling or even fluid-solid-thermal-chemical coupling is an urgent direction for the spatiotemporal evolution of reservoir damage.