The flow equation of injected acid in a porous medium is described by Darcy’s law, the formula is written as: v = − K μ ∇ P (1.1) Where, v is the vector of Darcy’s velocity, m/s; K is the tensor of permeability in study region, m²; μ is fluid viscosity, mPa·s; P is fluid pressure, Pa. The continuity equation is derived from the law of mass conservation, the formula is as follows: ∂ ϕ ∂ t + ∇ ⋅ v + Q m f = 0 (1.2) Where, ϕ is the porosity of carbonate rock; t is reaction time, s. Q _mf is the normal flow rate through the interface between the matrix system and the fracture system.

Generally, when the acid is injected into the carbonate rock, the equation of acid transport in the rock is described by the following mass balance equation. ∂ ( ϕ C f ) ∂ t + ∇ ⋅ ( v C f ) = ∇ ⋅ ( ϕ D e ⋅ ∇ C f ) + R + Φ m f (1.3) Where, C _f is acid concentration in a liquid phase, mol/m³; D _e is the tensor of the diffusion coefficient, m/s²; Φ _mf is solute mass transmitted into the matrix system through the interface between matrix medium and fracture medium. R is reaction term, its physical significance is that the amount of acid transferred to the liquid-solid surface is equal to the consumption of surface reaction, the formula of reaction term is expressed as: R = − a v ⋅ k c ( C f − C s ) (1.4) Where, a _v is specific surface area, m⁻¹; k _c is local transport coefficient of acid, m/s; C _s is the acid concentration of liquid-solid surface in core pore, mol/m³, it is calculated by: C s = C f ( 1 + k s k c ) (1.5)

The porosity will change as acid corrodes carbonate rock. The calculation formula of dynamic change of porosity is as follows: ∂ ϕ ∂ t = a v ⋅ α ⋅ k c ( C f − C s ) ρ s (1.6) Where, α is the dissolving power of acid, kg/mol; ρ _s is the density of carbonate rock, kg/m³.

The porosity, radius of the pore, and specific surface area of rock are decided by pore structure. It is a dynamic process and the injected acid changes the pore structure of rock when it continuously erodes porous media. This means that the different kinds of physical parameters of carbonate rock will change dynamically. There are two methods to compute the physical parameters of carbonate rock in reactive flow (Panga et al., 2005; Kalia and Balakotaiah, 2009; Maheshwari et al., 2013; Ghommem et al., 2015; Liu et al., 2020). In the first method, these various physical parameters are computed based on the dynamic change of pore structure obtained through the laboratory. Another method is to calculate the various physical parameters of rock according to empirical or semi-empirical formulas. This paper adopts the modified Carman Kozeny equation to describe the dynamic change of porosity and permeability, the formulas are as follows: { k k 0 = ϕ ϕ 0 ( ϕ ( 1 − ϕ 0 ) ϕ 0 ( 1 − ϕ ) ) 2 β r p r ¯ 0 = k ϕ 0 k 0 ϕ a v a ¯ 0 = ϕ r ¯ 0 ϕ 0 r p (1.7) Where, r _p is pore radius, m; a ¯ 0 is initial specific surface area, m-1; k ₀ is initial permeability, m²; ϕ ₀ is initial porosity; β is a constant that depends on the structure of the medium.

When solute of liquid phase flows in the pores, the velocity of solute in the acid through the pore flows from liquid phase to pore surface and contact with the pore surface is expressed by the rate of the mass transfer. It can be seen from Eq. 1.5 that the effect of the rate of mass transfer on reactive flow cannot be ignored when the acid system is certain, because its magnitude makes a difference to the chemical reaction. (Panga et al., 2010) defined a dimensionless number to the analysis of the associated factors, it is referred to the Sherwood number and represents the dimensionless mass transfer coefficient, it is given by: S h = 2 k c r ¯ p D m = S h ∞ + b Re p 1 2 S c 1 3 (1.8) Where, r _p is the average pore radius of porous media, m; D _m is the molecular diffusion coefficient, m/s²; Sh _∞ is asymptotic Sherwood number of the pore, b is a constant number related to the structure of porous media, b = 0.7/m ^1/2, where, the m is the ratio of pore length to pore radius, Re _p is pore Reynolds number, Re _p = 2ur _p/ν, where, ν is hydrodynamic viscosity; Sc is Schmidt number, Sc = ν _k/D _m, ν _k is hydrodynamic viscosity.

When the acid is transported in the homogeneous and anisotropic matrix, in 2D, the diffusion tensor is characterized by axial diffusion coefficient D _eX and transverse diffusion coefficient D _eT. There is only solute molecular diffusion in the liquid phase and the transverse diffusion coefficient equals the axial diffusion coefficient when the fluid does not flow, the formula is as follows: D e X = D e T = α o s D m (1.9) Where, D _m is the molecular diffusion coefficient, m/s²; α _os is a constant related to the structure of porous media (such as tortuosity or connectivity between pores). The two-dimensional diffusion coefficient of acid relates to the geometry of porous media, the flow pattern of pore scale, and the properties of acid. Researchers usually define a dimensionless Peclet number to describe the contrast between convection and diffusion. The expression formula of the Peclet number is found to be P e p = | v | d h ϕ D m (1.10) Where, | v | describes the magnitude of Darcy’s velocity, m/s; d _h is pore diameter of porous media, m.

The formula of diffusion coefficient derived by Panga (Panga et al., 2010) is used to describe the diffusion of acid in porous media in this paper. the diffusion coefficient can be computed as: D e X = α o s D m + λ X D m P e p (1.11) D e T = α o s D m + λ T D m P e p (1.12) Where, the subscript of X and T denotes the injection direction and vertical transverse direction of acid injection; α _os, λ _T, and λ _x are constants that have typical values of 0.5, 0.5,0.1 for a packed-bed of spheres (Kalia, 2008; Kalia and Balakotaiah, 2009; Liu et al., 2017), respectively.

The pressure at the outlet is constant and the outlet boundary is a constant pressure boundary. The upper and down boundaries are closed boundaries. The specific expression of boundary conditions is given by: P = 0 , C f = 0 ,   a t   t = 0 (1.13) u 0 = − k x v ∂ P ∂ x , ∂ P ∂ y = 0 , C f = C 0 ,   a t   x = 0 (1.14) u 0 = − k x v ∂ P ∂ x , ∂ P ∂ y = 0 , C f = C 0 ,   a t   x = 0 (1.15) ∂ C f ∂ x = 0 , P = P e ,   a t   x = L (1.16) ∂ P ∂ y = 0 , ∂ C f ∂ y = 0 ,   a t   y = 0   a n d   y = L (1.17) Where, u ₀ is initial injection velocity, m/s; C ₀ is a concentration of injected acid at the inlet, mol/m³; P _e is a constant, it is boundary pressure at the outlet, Pa.

A random function f ^ is defined to simulate the heterogeneity of rock. Its value varies from -Δϕ ₀ to +Δϕ ₀ and satisfies the random function of uniform distribution, the generated initial porosity field is shown in Figure 1. The heterogeneity of fracture aperture is defined in this way.

The results of the numerical simulation are shown in Figure 4. There are five different dissolution patterns with Damköhler number changing. They are face dissolution, conical dissolution, wormhole, branch dissolution, and uniform dissolution. Figure 4A is the concentration field, Figure 4B is the porosity field. Column (A) is the concentration field and porosity field in the initial phase. (B) is the concentration field and porosity field when injected acid breaks through the core. It can be seen that there is a big difference among dissolution patterns. The reason is that the dissolution patterns are determined by the combined effect of convection velocity and reaction velocity of the solute in acid. When the convection velocity is smaller than the diffusion velocity and reaction velocity, the reaction between acid and rock is mainly controlled by the diffusion and reaction velocity of solute in acid. Therefore, the hydrogen ions transferred from the center of the pore to the surface of the pore will be completely consumed by the chemical reaction, and the rate of injected acid cannot promote the injected acid to penetrate the core deeper. In this case, the dissolution of acid in the core is named surface dissolution. As the rate of injected acid gets higher, so does the convection velocity of solute in acid, and the convection velocity is dominant for diffusion velocity and reaction velocity.

The solute in the acid flows away before it reacts on the pore surface of carbonate rock. This results in the core pores being filled with acid and leads to a uniform variation in pore porosity. There is no dominant channel in the core like wormhole, and the dissolution pattern is called uniform dissolution.

As can be seen from Eq. 1.4, the reaction rate is determined by the properties of acid and rock, and have nothing to do with the existence of fractures in the matrix. In this case, we can draw an important conclusion: in the case of the same acid-rock reaction system, whether there are fractures in the matrix has no effect on the dissolution patterns of the fractured medium.

In addition, the comparison results of numerical simulation show that the fracture is equivalent to the main advantage channel in the matrix. An injection rate that is too high or too low will not substantially change the dissolution patterns. This is because the injected acid has already reacted with the rock completely when the acid is injected into the core at a low injection rate and there is no excess acid to flow into the fracture. In the case of high injection velocity, the dominant role of fracture in the flow is negligible and may even play a negative role. In the above situations, the existence of fractures does not affect the dissolution patterns, however, the dissolution structure will change when the injection rate of acid is between the two. This is because the flow of injected acid is controlled by fracture conductivity, and the ultimate dissolution pattern is closely related to the distribution of fracture in the fractured medium.

In this section, the effect of fracture orientations on the wormhole structure is studied when Damkohler number is 250. Firstly, a fracture is located in the matrix and then the dissolution structure is obtained by numerical simulation. The results of numerical simulation are compared with the results obtained in a matrix without fracture. The numerical simulation results of a porosity field and concentration field are shown in Figure 5.

(A) is concentration field and porosity field of numerical simulation in the initial phase and (B) is the concentration field and porosity field of numerical simulation when injected acid breaks through the core. In Figure 5, it can be obtained from column (A), the existence and orientations of fracture in the matrix do not affect the formation of the initial wormhole, but the direction of wormhole formation is consistent with the injected port of fractures. (B) shows that the fracture orientations influence the orientation of the wormhole and the fracture gradually becomes the part of the wormhole. The generation trend of the wormhole is random and branching in the matrix without fracture. When the orientation of fracture and injection direction of acid is at an angle of 90°, as shown in row (d) of Figure 5, the formation direction of the wormhole is not controlled by the orientation of the fracture. Instead, it is consistent with the direction of injected acid. Different from the wormhole in the matrix without fractures, the wormhole radius has increased and the direction of the wormhole has also changed in the fractured medium. This shows that the influence mechanism of fracture orientation on the formation direction of the wormhole is like to a high permeability zone.

The breakthrough volume of injected acid in a matrix with different fracture orientations is drawn in Figure 6. It shows that when the orientation of fracture is consistent with the injection direction of acid or the orientation of fracture is similar with horizontal direction, the breakthrough volume of injected acid is smaller than the breakthrough volume in the matrix without fracture, the branching randomness of the wormhole is also getting weaker. When the direction of fracture is bigger than a certain angle, the breakthrough volume of injected acid becomes larger, and so does the consumption of acid. At this time, the effect of fracture inhibits the flow of acid in the fractured medium. The waiting time in the fractured medium becomes larger and then consumes more acid. When the direction of fracture is between 67.5° and 112.5°, the breakthrough volume of injected acid increases first and then decreases. It is because the conductivity of fracture has a smaller impact on the development and formation of the wormhole. When the direction of fracture is bigger than 90°, the breakthrough volume of acid is similar to volume when the angle is less than 90°, the slight difference is due to the heterogeneity of the core.

In this example, the breakthrough volume of injected acid is largest in the core with a single fracture and the orientation of fracture is 67.5. The consumption of injected acid is even larger than the matrix without fracture.

In this section, the effect of fracture aperture on wormhole structure is studied when Damkohler number is 250. This model contains a fracture. The numerical simulation results are compared with the matrix without fracture, and the result of the numerical simulation is shown in Figure 7.

(A) is concentration field and porosity field of numerical simulation in the initial phase and (B) is the concentration field and porosity field of numerical simulation when injected acid breaks through the core in Figure 7. In column (A) of the concentration field and porosity field in Figure 7, whether the fracture exists or not and no matter how the fracture aperture changes, it does not affect the initial wormhole structure.

From column (B), it can be concluded that fracture aperture don't affect the formation trend of the wormhole and the fracture gradually becomes a part of the wormhole. The fracture aperture has no significant effect on the dissolution structure. The relationship between the breakthrough volume of injected acid and the fracture aperture is plotted as a columnar graph shown in Figure 8.

Figure 8 shows that the breakthrough volume of injected acid is the smallest when the fracture aperture is largest. In the mathematical model of this paper, the calculation equation of permeability is obtained by the cubic law: the larger the fracture aperture and the higher the fracture permeability, the conductivity of fracture will become stronger. This means that the fracture has a stronger impact on the formation of the wormhole. As a result, the breakthrough volume of injected acid will become smaller and the consumption of acid will become less. In the same way, the smaller the fracture aperture, the lower the fracture permeability. The efforts of fracture on the formation of the wormhole are weaker and the breakthrough volume and consumption of injected acid are less than before. When the fracture aperture equals 5e-5m, the breakthrough volume of injected acid is the smallest in the study.

The numerical simulation results are shown in Figure 9. With the increase of injected acid and the Damköhler number, there are five different dissolution patterns in fractured medium: face dissolution, conical dissolution, wormhole, branch dissolution, uniform dissolution.

The breakthrough volume curve of injected acid is plotted in Figure 10. It depicts the effect of homogeneity of the initial fracture aperture, the heterogeneity of initial fracture aperture, and the dynamic change of fracture aperture to the breakthrough volume of injected acid.

As shown in Figure 10, The dissolution patterns are unchangeable regardless of whether the initial fracture aperture is homogeneous, heterogeneous or fracture aperture is dynamic change with the acid etching fracture wall, the optimal injection rate of acid solution will not change when Da = 250.