Regular Article
Thermodynamically consistent Darcy–Brinkman–Forchheimer framework in matrix acidization^{★}
^{1}
College of Mathematics and Statistics, Shenzhen University, Shenzhen, Guangdong 518048, PR China
^{2}
School of Civil Engineering, Shaoxing University, Shaoxing, Zhejiang 312000, PR China
^{3}
School of Mathematics and Statistics, Hubei Engineering University, Xiaogan, Hubei 432000, PR China
^{4}
Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, King Abdullah University of Science and Technology, Thuwal 239556900, Saudi Arabia
^{5}
Department of Petroleum Engineering, Colorado School of Mines, 1600 Arapahoe Street, Golden, CO 80401, USA
^{*} Corresponding author: shuyu.sun@kaust.edu.sa
Received:
29
May
2020
Accepted:
17
November
2020
Matrix acidization is an important technique used to enhance oil production at the tertiary recovery stage, but its numerical simulation has never been verified. From one of the earliest models, i.e., the twoscale model (Darcy framework), the Darcy–Brinkman–Forchheimer (DBF) framework is developed by adding the Brinkman term and Forchheimer term to the momentum conservation equation. However, in the momentum conservation equation of the DBF framework, porosity is placed outside of the time derivation term, which prevents a good description of the change in porosity. Thus, this work changes the expression so that the modified momentum conservation equation can satisfy Newton’s second law. This modified framework is called the improved DBF framework. Furthermore, based on the improved DBF framework, a thermal DBF framework is given by introducing an energy balance equation to the improved DBF framework. Both of these frameworks are verified by former works through numerical experiments and chemical experiments in labs. Parallelization to the complicated framework codes is also realized, and good scalability can be achieved.
© Y. Wu et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Acidization is a useful technique for promoting or restoring oil production in reservoirs and can be classified as either fracture acidization or matrix acidization. In fracture acidization, a highly pressurized acid flow is injected into a well to physically enlarge the fractures and chemically dissolve the deposits that inhibit permeability. However, in matrix acidization, the pressure of the acid flow is not high enough to destroy the fractures, and thus, the acid flow can only enlarge the natural pores of the matrix. Both kinds of acidization attempt to enlarge the voids in reservoirs and ease the outflow of hydrocarbons from the subsurface matrix. Many previously published papers [1–3] have studied fracture acidization, but this work pays more attention to matrix acidization.
Theoretically, matrix acidization is a chemical dissolution–front instability problem. Many studies have focused on the factors that influence matrix acidization, such as the mineral reactive surface area [4], mineral dissolution ratio [5], and solute dispersion [6]. Numerically, four major models have been proposed to investigate matrix acidization, including the capillary tube model [7], the network model [8–10], the dimensionless model [11, 12] and the twoscale model [13–15]. Since the twoscale model can better predict dissolution patterns and more accurately capture the formation of wormholes, this work focuses on the twoscale model. The two scales indicate the Darcy scale and pore scale. In each scale, there are a series of equations used to describe the progress of matrix acidization. In earlier literature, the momentum conservation equation in the Darcy scale was described by Darcy’s law, so the twoscaled model was also named the Darcy framework in these works [15–18]. Later, due to the nature of matrix acidization, Wu et al. [19–22] concluded that the Darcy framework was insufficiently accurate to describe matrix acidization and provided the Darcy–Brinkman–Forchheimer (DBF) framework by adding the Brinkman term and Forchheimer term to the momentum conservation equation. Li et al. [23–25] further analyzed the numerical stability and accuracy of the DBF framework. However, the DBF framework still has a defect when processing the momentum conservation equation, which degrades its reliability. This work reviews the DBF framework in more detail than [19] and focuses on this defect, which is covered in Section 2. Moreover, in the DBF framework, a pseudo parameter ε is introduced into the mass conservation equation [26] to solve the linear system by an iterative solver HYPRE [27], which also degrades the reliability of the DBF framework. However, by replacing HYPRE with a direct solver MUMPS [28, 29], the introduction of ε is not necessary, which is realized in this work. Furthermore, the flowchart of the simulation should also be changed accordingly. Thus, the new framework provided in this work is called the improved twoscale model based on the DBF framework, or the improved DBF framework for short.
However, the frameworks above only consider the mass conservation law and momentum conservation law, and the energy conservation law is not included, which is a major drawback. Temperature is a key variable in the energy conservation law and has a significant influence on thermodynamic parameters such as the surface reaction rate and the molecular diffusion coefficient. However, these thermodynamic parameters are deemed constants in those frameworks. Moreover, in real applications, matrix acidization is performed in a subsurface environment where the matrix is warmed by terrestrial heat. Thus, the temperature should be considered as a major factor in matrix acidization, which promotes the necessity and reasonability of introducing the energy conservation equation into these frameworks.
In addition to this work, many researchers have acknowledged the temperature issue and upgraded the twoscale model in their own ways. For example, Li et al. [30] introduced a heat transmission equation in the form of radial flows to the twoscale model and produced simulation results near the wellbore. Ma et al. [2] also developed a temperatureinfluenced model based on the twoscale model and used it to simulate matrix acidization in fractured carbonate rocks. Kalia and Glasbergen [31] studied cases when the temperatures of the acid fluid and the matrix were different and simulated acidization in the matrix under both adiabatic and nonadiabatic conditions. They concluded that the fluid temperature could be designed as a parameter to control matrix acidization. Although these endeavors attempted to consider the thermal effect on matrix acidization, all of them were based on the Darcy framework, which is not accurate enough to simulate matrix acidization, as mentioned above. As a result, their reliability was degraded. Therefore, this work provides a heat transfer model as an expansion of a more reasonable improved DBF framework and aims to output more reliable results. The new model is a thermodynamically consistent DBF framework, which can be called the thermal DBF framework for short. Different from [30], this work studies matrix acidization in the form of linear flows. Meanwhile, fractures in the matrix are not considered, and a general matrix is acidized, which is different from [2]. Inspired by [31], two cases are studied – when the temperatures of the acid flow and the matrix are the same and when they are different – and the results are verified against [31] and [32], respectively.
The work in [19] only realized a 2D parallel code of the DBF framework; the present study further realizes a 3D parallel code of the DBF framework. In addition, the 2D and 3D parallel codes of the improved DBF framework and thermal DBF framework are also realized in the present work.
In the following discussions, the improved DBF framework is developed, and then the thermal DBF framework is provided. In the model verification section, the correctness of the improved DBF framework is checked by comparing its 2D and 3D results with existing works. Only when the correctness of the improved DBF framework is guaranteed can the reasonability of the thermal DBF framework be assured. After this, a series of thermal experiments are performed to investigate the temperature effect on matrix acidization. The performance of the 3D parallel code is evaluated at the end of this work.
2 Improved DBF framework and its solution scheme
The meanings of all the notations in the statements below are listed in Table 1. The Darcy framework has been widely used to simulate the matrix acidization procedure. At the pore scale, a group of semiempirical equations is provided to describe the relationship of parameters at the pore scale, such as porosity, permeability [33], and local masstransfer coefficient. These equations have seen few changes during the study of matrix acidization simulation. However, at the Darcy scale, which is the other scale of the Darcy framework, these equations have seen many changes with continued study. Generally, there are three kinds of equations on the Darcy scale: momentum conservation equations, mass conservation equations, and concentration balance equations. The changes have mainly occurred in the momentum conservation equations, while the other two kinds of equations have remained more or less the same, although some small changes have been made according to the needs of different cases. Initially, the momentum conservation equation is represented by Darcy’s law:
assuming that the permeability is homogeneous and isotropic. However, the applicability of this equation is limited to conditions where the Reynolds number Re < 1 [34] and the Darcy number Da ≪ 1. At the beginning of matrix acidization, the porosity in the porous medium is not high, and Darcy’s law can be leveraged to properly describe the flow of fluid in the porous medium. However, with the propagation of channels due to matrix acidization, the areas eroded by the channels become enlarged. In these channels, the porosity can grow and even approach the value of one, which creates high permeability in these areas. From the definition of the Darcy number, it can be determined that the Darcy number is much higher than one as a result of this change. Moreover, in highpermeability areas, the fluid velocity increases, which may lead to a high Reynolds number that can exceed the value of one, assuming that the viscosity and mass density of the fluid remain more or less the same and that the particle size remains constant. These reasons indicate that Darcy’s law is not suitable for use in matrix acidization simulations.
Nomenclature.
To address the issue, some corrections have to be made to Darcy’s law to cope with the conditions of high permeability and high Reynolds numbers. The first correction is called the Brinkman correction. According to Darcy’s law, the uniform velocity in a crosssectional direction can be seen when the permeability is low. However, in a porous medium with high permeability, a noslip condition should be considered instead of uniform velocity. The Brinkman correction introduces a viscous shear stress term to Darcy’s law that can describe the noslip well. The Brinkmancorrected Darcy’s law can be expressed as,
in which ∇^{2} u is called the Brinkman term. In addition, under the condition of high Reynolds numbers, form drag can be much larger than viscous drag, which can be suitably described by a Forchheimer correction. In this correction, a term called the Forchheimer term, which is expressed as , is added to Darcy’s law. Combining both corrections together, Darcy’s law can be modified as,
with as the Forchheimer coefficient [35]. In the equation above, the righthand side equals zero, which means that the sum of all the external forces imposed on the fluid is zero. However, this is not a general case. For cases in which the sum of all the external forces is not zero, the righthand side of the equation should equal the product of mass density and acceleration. With a Eulerian expression of acceleration, the righthand side of the equation can be written as,
If the gravity effect is considered, the momentum conservation equation in its final form can be written as,
Equation (1) is the momentum conservation equation used in our previous work [19]. Because this kind of momentum conservation equation introduces the Brinkman term and Forchheimer term, this model is a twoscale model based on the Darcy–Brinkman–Forchheimer (DBF) framework, or the DBF framework for short. More details can be found in [19]. However, the framework cannot meet Newton’s second law, since porosity ϕ is changed during the simulation procedure. It is noted that in the first term on the lefthand side of equation (1), ϕ is outside of the time derivative, which indicates that ϕ does not change with time. This contradicts the true physical observation, and thus Newton’s second law is also violated. Thus, ϕ should be placed inside the time derivative. After this operation, can be deemed a new variable – the effective velocity. Accordingly, the second term on the lefthand side of equation (1) and the third term on the righthand side of equation (1) are changed, which can be expressed as,
It is noted that the lefthand side of equation (2) is in reality the material derivative,
The discussion above, taken from the theories of fluid dynamics, demonstrates that the new momentum conservation equation should describe matrix acidization more reasonably, and thus the new model is an “improved” twoscale model based on the DBF framework, or the improved DBF framework for short.
Since the flow in matrix acidization is assumed to be incompressible, the mass conservation equation can be expressed as,
However, considering the local volume change in the matrix acidization procedure, the mass conservation equation should be modified as,
The concentration balance equation can be derived from the principle of species balance during matrix acidization. The balance of species can be achieved by accumulation, advection, diffusion, and reaction effects, which bring about a straightforward expression of the concentration balance equation,
In the equation, D_{e} is a function of u,
In the 3D condition,
u_{x}, u_{y}, and u_{z} represent the xdirection, ydirection, and zdirection velocities, respectively. The 2D conditions are similar. In fact, the concentration balance equation is the other expression of the mass conservation law, in addition to the mass conservation equation.
In equations (2)–(4), the velocity vector u, pressure p, and the cupmixing concentration of the acid C_{f} are deemed unknowns to be solved for. However, the three equations cannot accomplish this since the values of the other variables are unknown. Thus, more auxiliary equations and other necessary assumptions are needed. Three additional equations are given as follows:
These equations are derived from mathematical deduction and chemical experiments, and more details can be found in [15]. It should be emphasized that equation (7) can be put in equation (4) to substitute C_{s} when solving C_{f}. Moreover, since the mass density ρ_{f} and viscosity μ of the fluid will not change much during the matrix acidization procedure, they should be given constants for simplicity. It is also easy to understand that the dissolving power of the acid α and the mass density of the solid phase ρ_{s} can also be deemed as given constants. For the reason mentioned later, the surface reaction rate k_{s} is no longer deemed as a constant, as shown in [19]. Instead, it is a function of the temperature T in this work. As a result, the reaction rate also becomes a function of T and can be rewritten as R(C_{s}, T), which is different from [19].
In addition, it is noted that some variables at the pore scale, such as porosity, permeability, the interfacial surface area per unit volume and the local masstransfer coefficient, appear in the equations at the Darcy scale. Thus, their values should be known before we attempt to solve for the unknowns on the Darcy scale with the help of a series of equations on the pore scale. It is stipulated that the subscript 0 represents the initial value or reference value of the corresponding variable in the following equations, and all the initial values are known. First, three equations called the Carman–Kozeny correlation are provided as follows:
From these equations, it can be seen how the permeability, pore radius, and interfacial surface area per unit volume change with porosity. Thus, as long as the porosity is known, the values of these three variables can be computed from the Carman–Kozeny correlation. At this point, we can compute porosity. From equations (5) to (7), the equation below can be derived:
The lefthand side of equation (11) describes the change in porosity with time, and its righthand side includes many variables; except for a_{v} and k_{c}, all the other variables have no direct relationship with porosity. By using equations (8)–(10), a_{v} can be expressed as a function of ϕ,
Equation (11) can then be changed as,
Moreover, the local masstransfer coefficient k_{c} can be calculated from the expression of the Sherwood number Sh, which is a dimensionless masstransfer coefficient. The expression is given as follows:
On the righthand side of equation (14), Sh_{∞} is a given constant, and the Reynolds number Re can be expressed as,
The Schmidt number Sc is expressed as,
Thus, by leveraging equations (14)–(16) together, k_{c} can be calculated as,
Equation (17) includes the variable r_{p}, which is a function of porosity from equations (8) and (9),
(18)which means that k_{c} is in fact a function of porosity and velocity. Thus, equation (17) can substitute for k_{c} in equation (13), and then a new equation with porosity being the unknown can be derived. However, the new equation is too complex to derive an analytic formula of porosity, and therefore, the value of porosity has to be computed by a numerical scheme. Under that condition, a semiimplicit scheme is applied. In the following statements, the superscripts of the notations represent the time step. The porosity at time step τ is used in equation (18) to calculate the pore radius at time step τ, which is then put into equation (17) to compute k_{c} at time step τ. With the semiimplicit scheme, equation (13) can be rewritten as,
(19)by which ϕ^{τ+1} can be computed easily.
From the discussion above, it can be seen that the main unknown to be solved at the pore scale is porosity ϕ. As long as the porosity value is obtained, the values of the other variables at the pore scale can be derived from a series of porescale equations. Furthermore, with the values of all the porescale variables, the main unknowns in the Darcy scale can be calculated by the Darcyscale equations. In addition, from the computation of porosity, it can be seen that the variables u and C_{f} in the Darcy scale will affect the porosity value. Therefore, the computations of the Darcy scale and pore scale are coupled with each other, with the Darcyscale variables u and C_{f} and the porescale variable ϕ being their interaction media, which can be shown in Figure 1.
Fig. 1 The interaction of the pore scale and Darcy scale. The upper arrow indicates that the pore scale will affect Darcy scale by the porescale variable ϕ. As long as ϕ is changed, the porescale variables used in the Darcy scale are changed, which brings about the changes of the main variables u, p, and C_{f} in the Darcy scale. The lower arrow indicates that the Darcy scale will affect the pore scale by the Darcyscale variables C_{f} and u, since the two variables will affect the value of ϕ directly. Once ϕ is changed, all the other variables in the pore scale will be changed accordingly. 
The derivations of all the equations at both the Darcy scale and pore scale have been achieved in the former discussions, and then numerical schemes are applied to these derivations to solve for the variables. Since porosity ϕ plays a central role in the improved DBF framework, it should be computed from equation (19) with the semiimplicit scheme first. Then, with the computed porosity, permeability K and the interfacial surface area per unit volume a_{v} can be computed with equations (8) and (12), respectively. Next, equations (2) and (3) are combined together as a linear system and solved for the velocity u and pressure p with the semiimplicit scheme. It is emphasized that the term in equation (3) is substituted by equation (6) when the linear system is solved. Moreover, with the update of porosity ϕ, the local masstransfer coefficient is also updated by equations (17) to before solving the linear system. With the semiimplicit scheme, equations (2) and (3) can be rewritten as follows:
After that, since the velocity u is updated, the local masstransfer coefficient can be updated again by equations (17) to . Last, the semiimplicit scheme is used to solve equation (4) for concentration C_{f}, and then another linear system is formed. It is emphasized that equation (7) is put into equation (4) when the linear system is solved, which is shown as,
In brief, the solution procedure can be described as a flowchart, which is shown in Figure 2.
Fig. 2 The flowchart of the improved DBF framework. 
From the flowchart, another difference of this work from our former work [19] can be seen. In the former work, for each iteration, the simulation begins with the computation of the variables in the Darcy scale, such as pressure, velocity, and concentration, and ends with the computation of the variables in the pore scale, such as porosity, permeability, and the interfacial surface area per unit volume. However, in the present work, the computation of the variables at the pore scale is performed ahead of the computation of the variables at the Darcy scale. The flowchart of this work is more reasonable, which can be demonstrated by the following statements. From equation (20), it is observed that to obtain , we have to know ϕ^{τ+1}, and therefore the computation of ϕ^{τ+1} should occur before the computation of . However, in our former work, the computation of ϕ^{τ+1} occurs after the computation of , which is not reasonable. With this philosophy, in the simulation of the improved DBF framework, the variables are first computed at the pore scale and then at the Darcy scale.
3 Thermal DBF framework and its solution scheme
Based on the improved DBF framework, a heat transfer model that considers the heat transmission process in matrix acidization is developed. The model is composed of the improved DBF framework and the energy conservation equation, which can be expressed as the governing equation of the temperature T,
H_{r}(T) is the reaction heat [30]. λ_{f} and λ_{s} are the heat conduction coefficients of the fluid and solid phase, respectively, and thus, λ is the average heat conduction coefficient between the two phases. θ_{f} and θ_{s} are the heat capacities of the fluid and solid phase, respectively. λ_{f}, λ_{s}, θ_{f}, and θ_{s} are deemed constants in this work. ϑ_{f} and ϑ_{s} are the amounts of heat per unit volume of the fluid and solid phase, respectively, and thus, ϑ is the total amount of heat per unit volume. The total energy is a sum of the fluid energy and solid medium energy, which may vary with time. The energy transportation caused by the fluid flow occurs only between fluids in different spatial positions due to fluid flow. The heat conduction in the interiors of both fluids and solids and with each other are considered in the first term of the righthand side of (21). The effects of the pressure and viscosity force of the fluid are described in the second and third terms of the righthand side of (21). The effects of friction forces between fluids and solids are described in the fourth and fifth terms of the righthand side of (21) based on the Darcy–Forchheimer framework. The chemical reaction can produce heat, which is considered in the last term of (21). All the terms constitute the source or sink of energy. It is emphasized that the temperature of the acid and matrix is assumed to become the same immediately when acid is injected into the matrix, since the speed of heat transfer is much faster than the fluid speed. Thus, the temperature of the acid and matrix can be represented by a single notation T. In fact, differentiating between the acid temperature and matrix temperature brings challenges to theories and applications, and the details of the heat transfer between the acid and matrix must be studied thoroughly. Relevant work can be left to the future. Furthermore, the surface reaction rate k_{s} is deemed a function of the temperature T, which can be expressed as,
(22)in which k_{s0} is the surface reaction rate at temperature T_{0}, E_{g} is the activation energy, and R_{g} is the molar gas constant [31]. As a result, the reaction rate also becomes a function of T, and its expression is rewritten as R(C_{s}, T) in (21). Moreover, the molecular diffusion coefficient d_{m} is also affected by the temperature T, which can be expressed as [30],
(23)d_{m0} is the molecular diffusion coefficient at temperature T_{0}. The heat transfer model induced from the improved DBF framework is called the thermal DBF framework for short.
After solving a series of equations in the improved DBF framework, the semiimplicit scheme can be used to solve equation (21) for the temperature T, which can be expressed as,
and the third linear system is formed. It is noted that equation (7) is put into (21) in the above expression. When using the thermal DBF framework to simulate matrix acidization, the flowchart is shown in Figure 3. It can be seen in the figure that the molecular diffusion coefficient d_{m} and the surface reaction rate k_{s} are calculated first, followed by a series of computations in the improved DBF framework, and the computation of the temperature is performed last. It can also be seen that by using the variables ϕ, u, p, and C_{f}, the improved DBF framework changes the main variable T in the energy conservation equation, while the energy conservation equation also changes the variables of the improved DBF framework by T, or k_{s} and d_{m}, which can be seen in Figure 4. From these results, it can be determined that the thermal DBF framework considers all three kinds of conservation laws: mass, momentum and energy, and should simulate matrix acidization more reasonably.
Fig. 3 The flowchart of the thermal DBF framework. 
Fig. 4 The interaction of the improved DBF framework and the energy conservation equation. 
At the end of this section, the relationships among the Darcy framework, the DBF framework, the improved DBF framework and the thermal DBF framework can be summarized. Initially, the Darcy framework is provided to simulate matrix acidization and achieves great success. In the framework, the communication between the pore scale and Darcy scale advanced the progress of simulation. Then, considering the clear fluid area that cannot be described accurately by Darcy’s law, the Brinkman term and Forchheimer term are introduced to the momentum conservation equation, and the DBF framework is developed. However, this framework cannot obey Newton’s second law, and a modification is made, thus suggesting the improved DBF framework. Until now, all the frameworks have only considered two kinds of conservations, i.e., mass conservation and the momentum conservation, which is insufficient. Thus, based on the improved DBF framework, a third kind of conservation, energy conservation, is introduced in the thermal DBF framework. Their relationships can also be seen in Figure 5.
Fig. 5 The relationship among the Darcy framework, the DBF framework, the improved DBF framework, and the thermal DBF framework. 
4 Discretization and parallelization
Possible discretization methods include the multipoint flux approximation method and the hybrid finite volume method. The multipoint flux approximation method is “designed to give a correct discretization of the flow equations for general nonorthogonal grids as well as for general orientation of the principal directions of the permeability tensor” [36], while the hybrid finite volume method is “the ideal method for computing discontinuous solutions arising in compressible flows” [37]. Since this work considers orthogonal grids and incompressible flows, the finite difference method should be a better fit from all possible discretization methods. In the following discussion, the finite difference method is used to discretize the model, and the experimenting field approach [38–42] is used to compute the coefficients in the two linear systems. Although these approaches have been used in the 2D simulation of matrix acidization in our previous work [19], it is necessary to expand them to the 3D simulation, which is a focus of this work.
The equations used in the thermal DBF framework are discretized one by one according to the flowchart shown in Figure 3. For simplicity, suppose there is a 3D Cartesian grid. For equations (8), (12), (19), (22), and (23), every variable is imposed at the center of the cube. For equation (17), except the variable u, the other variables are imposed at the center of the cube. However, should also be imposed at the center of the cube. Generally, u = (u_{x}, u_{y}, u_{z}) and its three components are imposed on the faces of the cube. In other words, the xdirection velocity u_{x} is imposed on the xdirection face, which is vertical to the xaxis. The processes of the ydirection velocity u_{y} and zdirection velocity u_{z} are similar. Thus, for a cube with its xcoordinate ranging from i to i + 1, ycoordinate ranging from j to j + 1 and zcoordinate ranging from k to k + 1, there is,
with the subscript representing the coordinate. Equation (2) is discretized on the faces of the cube, with the xdirection momentum equation discretized on the xdirection face, the ydirection momentum equation discretized on the ydirection face and the zdirection momentum equation discretized on the zdirection face. Thus, the porosity and permeability on the faces should be known. However, from the discussion above, it can be seen that the porosity and permeability are imposed at the center of the cube. Thus, the harmonic method must be applied to obtain their values on the faces. It is emphasized that the advection term in equation (2) is discretized with the upwind scheme. After performing all the operations, the xdirection momentum equation imposed on the xdirection face with its xcoordinate being i, ycoordinate ranging from j to j + 1 and zcoordinate ranging from k to k + 1 can be discretized as below:
in which and represent the ydirection average velocity and zdirection average velocity on the face, respectively. is computed as the average of the ydirection velocities on the four ydirection faces adjacent to the face. The computation of is similar. g_{x} is the xdirection component of g. It is noted that in the above equation, it is assumed that,
The discretization of the ydirection and zdirection momentum equations is similar. Under the Neumann boundary condition for pressure, the momentum conservation equation on the boundary degenerates to,
in which u_{B} is the boundary normal velocity. The discretization of equation (3) is straightforward. Equation (4) is discretized at the center of the cube. For a cube with its xcoordinate being from i to i + 1, ycoordinate being from j to j + 1 and zcoordinate being from k to k + 1, the lefthand side of equation (4) is discretized as,
with the upwind scheme being used to obtain the concentration value on the face of the cube, since the computed concentration is imposed at the center of the cube. If D_{ e } is written as,
the first term of the righthand side of equation (4) can be discretized as,
The notation stands for the average value of on the xdirection face with its xcoordinate being i + 1, which is calculated by the four on the four ydirection faces adjacent to the xdirection face. The meanings of the other similar notations are analogous. The discretization of the second term of the righthand side of equation (4) is trivial and thus not given here. It is easy to see that the stencil pattern of T is the same as C_{f}; thus, the discretization of equation (21) holds the same philosophy as that of equation (4), and the details are no longer given.
To use the experimenting field approach to compute the coefficients of the three linear systems, the unknowns to be computed can be divided into four fields: the velocity field, the pressure field, the concentration field, and the temperature field. If there is a 3D domain and it can be divided into eight cubes as shown in Figure 6, then the velocity field can be represented as arrows on the faces of the cubes, and the pressure field, concentration field and temperature field can be represented as points at the centers of the cubes. Each xmomentum, ymomentum, and zmomentum conservation equation can be discretized on each xdirection, ydirection, and zdirection face, respectively. Each mass conservation equation, concentration balance equation, and energy conservation equation can be discretized at each center of the cube. This kind of grid is called a staggered grid in CFD. The experimenting field approach used in the 2D simulation [19] is expanded to the 3D case directly, and the details are no longer given in this work.
Fig. 6 Staggered grid. 
To capture the details of the configuration of the matrix after acidization, a fine 3D grid is needed in the simulation, which brings about a large number of cells in the 3D grid, and as a result, it is necessary to introduce parallelization in the simulation. In the first step, domain decomposition must be performed on the 3D domain. The main purpose of domain decomposition is to allocate the discretized equations to the processors. Suppose there is a 3D Cartesian grid with nx, ny, and nz cubes in the x, y, and zdirections, respectively, and there are npx, npy, and npz processors in the x, y, and zdirections, respectively. nx, ny, and nz are supposed to be divisible by npx, npy, and npz, respectively. Furthermore, it is stipulated that , , and . For the processor with the coordinate (I,J,K), the following equations discretized at the centers of the cubes with the coordinate (i, j, k) are allocated to it,
and the equations discretized on the xdirection faces with the coordinate (i, j, k), which are the xmomentum conservation equations, are allocated to it as,
The allocation strategy of the ymomentum conservation equations and the zmomentum conservation equations is similar. After the allocation of the discretized equations, the variables that are needed by the equations should also be allocated to the processors, with some variables being communicated among the processors. By the domain decomposition strategy, each discretized equation can be allocated to only one processor, with the benefits that it can keep load balance of the processors and reduce the communication cost among the processors.
After domain decomposition, a suitable parallel solver can be leveraged to solve the three linear systems. In the 2D parallel simulation of the work [19], the parallel solver HYPRE is used. There are many numerical algorithms that can be used in HYPRE, such as the Generalized Minimal RESidual method (GMRES) and the Algebraic Multigrid Method (AMG). However, few simple cases can be solved by these algorithms. As a result, more complicated cases have to be solved by the direct solver UMFPACK [43] in a serial code, which limits the application of the parallel code. The main issue of these complicated cases is that their condition number is too large. In this work, another parallel solver called MUMPS is used, which can solve complicated cases that HYPRE cannot solve. With the help of MUMPS, most of the cases for the improved DBF framework and the thermal DBF framework can be run in parallel, which makes the fine 3D simulation feasible. In MUMPS, a variant of the Gaussian elimination method – the multifrontal method – is used. This method can solve a large sparse system of equations in parallel: the equation system is first divided into independent subsets, which are called fronts, and then the fronts are processed on different processors simultaneously. It is emphasized that unlike HYPRE, which is an iterative solver, MUMPS is a direct solver, which makes the time step larger in the simulation. Moreover, the direct solver can solve the linear system directly, with no need to add a pseudo parameter ε in the mass conservation equation, which is the case in the work [19]. There, ε is introduced to ensure a linear system with an invertible coefficient matrix; otherwise, the iterative solver HYPRE cannot solve it. However, the introduction of ε changes the attribute of the flow in matrix acidization from incompressible to slightly compressible, which contradicts the real case and makes the DBF framework less reliable. Last, with the help of FORTRAN 90 and MPI, a series of 2D and 3D parallel codes are developed. In the following sections, these codes are used to run a series of numerical experiments on the Shaheen supercomputer [44].
5 Verification of the improved DBF framework
5.1 3D sheardriven cavity flows
It can be seen that the model of sheardriven cavity flows [45] is in fact a reduction of the improved DBF framework, since equations (2) and (3) reduce to the following two equations:
which describes sheardriven cavity flows. Therefore, the 3D code of the improved DBF framework can be used to simulate 3D sheardriven cavity flows identically as long as some parameters are simplified, such as ignoring the concentration balance equation. In terms of 3D sheardriven cavity flows, a laminar incompressible flow is inside a unit cube cavity whose ydirection top surface is moved by an xdirection uniform velocity of 1 m/s, as shown in Figure 7. The Reynolds number (Re) is 100. The gravity effect is ignored. The grid has 20^{3} cubes. The simulation results of stable flows are shown in Figures 8 and 9, respectively. The two figures display the velocity profiles of the xdirection component on the vertical centerline and the ydirection component on the horizontal centerline of the plane z = 0.5. The simulation results can be compared with Figure 6 in [45]. To the naked eye, we cannot clarify their differences, which proves the correctness of the 3D code of the improved DBF framework to some extent.
Fig. 8 Velocity profile of the xdirection component (u) on the vertical centerline of the plane z = 0.5. 
Fig. 9 Velocity profile of the ydirection component (v) on the horizontal centerline of the plane z = 0.5. 
5.2 2D linear flows
The 2D linear flows in a previous work [17] are simulated again by the improved DBF framework, with more or less the same experimental parameters, which are shown in Table 2. In [17], the flows are simulated with the Darcy framework [15], but NavierStokes fluid dynamics are considered. Thus, among all the stateoftheart models developed from the Darcy framework, the model is close to the improved DBF framework, and its results can be compared with the results from the improved DBF framework. It is noted that the subscript 0 represents the initial value. represents the initial average porosity in the medium, with a heterogeneity magnitude of 0.03. The gravity effect is ignored. In the 2D simulation, there is a rectangular matrix of 0.1m length (xdirection) and 0.04m width (ydirection). Acid flow is injected into the matrix from the left boundary and goes out of the matrix from the right boundary, which means that the injected velocity is imposed on the left boundary, and the Dirichlet boundary condition for pressure is imposed on the right boundary. It is stipulated that the pressure imposed on the right boundary is the same as the initial pressure in the matrix. For concentration, the Dirichlet boundary condition is imposed on the left boundary, and the noflux boundary condition is imposed on the right boundary. The upper and lower boundaries are closed for both pressure and concentration, which means that noflow, noflux boundary conditions are imposed. The acid concentration is initially zero in the matrix. The injected velocity of the acid flow of 0.5 M hydrogen chloride (HCl) on the left boundary is changed in the simulation, which leads to different configurations of the matrix after acidization. The grid has 180 cells in the xdirection and 72 cells in the ydirection, which is the same size as that in [17]. Since the author of [17] declared that their grid is fine enough to describe matrix acidization, our grid is also capable of doing so.
Experimental parameters.
The Pore Volumes to BreakThrough (PVBT) of different injected velocities are given in Figure 10, and the curve in the figure is called the acidefficiency curve. Breakthrough is defined as the moment when the pressure drop across the medium drops to 1% of its initial value [46]. In Figure 10, it can be seen that the acidefficiency curve in this work matches the corresponding curve in Figure 8 of [17] well when injected velocity ranges from 4.17 × 10^{−7} m/s to 1.67 × 10^{−4} m/s. When the injected velocity is 1.67 × 10^{−7} m/s, its corresponding PVBT is 6.62, which is not reasonable. The injected velocity of 1.67 × 10^{−7} m/s is very slow, which may indicate a face dissolution pattern. According to [16], under such conditions, approximately 100 million grid cells are needed to capture the face dissolution pattern accurately, which indicates that our grid is not fine enough to output accurate results. Unfortunately, due to the limits of supercomputing power, a simulation based on 100 million grid cells cannot be finished in a reasonable time, so it is not performed in this work. The minimum PVBT is 4.54, which is achieved at an injected velocity of 4.17 × 10^{−6} m/s.
Fig. 10 Acidefficiency curve of the 2D linear flows. The numbers beside the points represent the values of PVBT. 
A fixed time step is assumed for the simulations, and the time steps corresponding to the injected velocities are shown in Table 3. It is emphasized that all the time steps ensure that the Courant number is less than one. A sensitivity test is performed when each of the time steps is increased by two times to assure that the Courant number is less than one and all the values of PVBT are the same, which demonstrates that the simulation results in Figure 10 can be deemed the true results.
Time steps for the injected velocities of the 2D linear flows.
The porosity profiles at breakthrough corresponding to five different injected velocities are given in Figure 11. In the figure, it can be seen that five dissolution patterns appear in their turns when the injected velocity increases.
Fig. 11 Porosity profiles at breakthrough for five different injected velocities of the 2D linear flows. (a) u_{x} = 4.17 × 10^{−7} m/s, face dissolution. (b) u_{x} = 1.67 × 10^{−6} m/s, conical wormhole. (c) u_{x} = 4.17 × 10^{−6} m/s, dominant wormhole. (d) u_{x} = 7.17 × 10^{−6} m/s, ramified wormhole. (e) u_{x} = 1.67 × 10^{−5} m/s, uniform dissolution. 
5.3 3D linear flows
The simulation of 2D linear flows has achieved reasonable results, and it can be expanded to the simulation of 3D linear flows by adding another dimension to the matrix, with a length of 0.04 m. After this addition, a 3D matrix is created with a 0.1m length in the zdirection and 0.04m length in the x and ydirections, respectively. According to [17], dissolution patterns from a conical wormhole to uniform dissolution can be captured accurately only when the grid has at least 180 cells in the zdirection and 72 cells in the xdirection and ydirection. However, the face dissolution pattern requires a finer grid. To ensure that the Courant number is less than one, the number of iteration steps can be large, which brings about a long simulation period. Even though the code runs on Shaheen, at least one month is needed for the fastest case to achieve breakthrough. Thus, that kind of grid currently exceeds computing capacity, and a coarser grid is thus given in this work. The number of cells is divided by two in each dimension, and a coarser grid, with 90 cells in the zdirection and 36 cells in the xdirection and ydirection, is used to simulate the 3D linear flows. Although the coarser grid is not fine enough to capture all dissolution patterns accurately, the simulation results can still be used to verify the correctness of the 3D code of the improved DBF framework. The acid flow is injected into the matrix along the zdirection. The other experimental parameters and boundary conditions are the same as the 2D simulation.
The time steps for different injected velocities are shown in Table 4. There are two groups of time steps, with the values of the second column being two times the corresponding values of the third column, the purpose of which is to test the convergence of the results. All the time steps guarantee that the Courant number is less than one. It is emphasized that due to the limits of Shaheen, a code can run on the supercomputer for three days at most. Therefore, simulations for velocities of 3.04 × 10^{−7} m/s and 1.04 × 10^{−4} m/s at fine time steps are not performed since their simulation time to achieve breakthrough is beyond three days. The acidefficiency curves for both groups of time steps are shown in Figure 12. The numbers beside the points represent the values of PVBT, with the blue points coming from coarse time steps and the red points coming from fine time steps. In the figure, it can be seen that the values of PVBT from coarse time steps are very close to those from fine time steps, which demonstrates the convergence of the results. Moreover, the values of PVBT tend to decrease with finer time steps. It can be expected that with finer grids, smaller PVBT values can be achieved. The minimum PVBT is 3.567, which is achieved at the injected velocity of 1.04 × 10^{−6} m/s. This coincides with the work [17], where the minimum PVBT is achieved at the injected velocity of 0.1 cm^{3}/min, which is the same as 1.04 × 10^{−6} m/s. However, in [17], the minimum PVBT is approximately two, which is smaller than ours due to the finer grid. Except for the points at the injected velocity of 1.04 × 10^{−7} m/s, the shape of the acidefficiency curves matches the corresponding shape in the [17] well. The drop in PVBT values at that velocity is due to the inaccuracy of the simulation when the grid is not fine enough, which can also be seen in the simulation of 2D linear flows above. Furthermore, the minimum PVBT of 2D simulations is larger than that of 3D simulations, and the injected velocity of 2D simulations at which the minimum PVBT is achieved is also larger than that of 3D simulations, which conforms to the qualitative trends in [15].
Fig. 12 Acidefficiency curves of the 3D linear flows. The numbers beside the points represent the values of PVBT. 
Time steps for the injected velocities of the 3D linear flows.
These effects of the injected velocities on dissolution patterns are shown in Figure 13. It can be seen in the figure that five different dissolution patterns can be simulated.
Fig. 13 Porosity isosurfaces at breakthrough for five different injected velocities of the 3D linear flows. (a) u_{z} = 3.04 × 10^{−7} m/s, face dissolution. (b) u_{z} = 7.04 × 10^{−7} m/s, conical wormhole. (c) u_{z} = 1.04 × 10^{−6} m/s, dominant wormhole. (d) u_{z} = 3.04 × 10^{−6} m/s, ramified wormhole. (e) u_{z} = 1.04 × 10^{−5} m/s, uniform dissolution. 
6 Verification of the thermal DBF framework
6.1 Isothermal conditions
The correctness of the improved DBF framework is a major premise of the thermal DBF framework, which has been verified in the last section. The correctness of the thermal DBF framework is discussed in this section. First, an experiment is carried out with isothermal conditions in which the injected acid temperature and the initial matrix temperature are the same. Since 2D experiments are eligible to verify the correctness of the model, the grid of the 2D linear flows is used again, with 180 × 72 cells in total. 3D experiments are left to future work. To compare the numerical results with the chemical results of the “effects of temperature” experiment in [32], three temperatures are chosen: 295 K, 323 K, and 353 K, which correspond to 22 °C, 50 °C, and 80 °C in [32]. The boundary conditions and initial conditions for pressure and concentration are the same as those in the 2D linear flow experiment above. In addition, for temperature, adiabatic conditions are applied, which means that except for the acid injection boundary (left boundary), all the other boundaries are adiabatic boundaries. The experimental parameters can be seen in Table 2. It is noted that the values of d_{m} and k_{s} in Table 2 are not used in the experiments of this section, since they are variables in the thermal DBF framework. All the parameters are more or less the same as those in the “effects of temperature” experiment of [32].
The values of PVBT for different injected velocities and temperatures are shown in Table 5. It can be seen in the table that when the temperature is 295 K, the minimal PVBT is 4.350, which is achieved at the optimal injected velocity of 2.67 × 10^{−6} m/s; when the temperature is 323 K, the minimal PVBT is 4.362, which is achieved at the optimal injected velocity of 9.17 × 10^{−6} m/s; when the temperature is 353 K, the minimal PVBT is 4.416, which is achieved at the optimal injected velocity of 4.17 × 10^{−5} m/s. To verify the convergence of the results, both coarsetimestep and finetimestep results are computed. The coarse time step is two times the fine time step for every injected velocity. All the time steps can guarantee that the Courant number is less than one. It can be seen in the table that the differences of the coarsetimestep and finetimestep results are very small, which means convergence is achieved. The values of PVBT from fine time steps constitute the acidefficiency curves of different temperatures in Figure 14. Since these experimental parameters hint that the injected acid is 0.5 M HCl, the matrix is limestone; Figure 14 is compared with Figure 6 in [32]. In the two figures, it is evident that both the minimal PVBT and the optimal injected velocity increase when the temperature increases, which means that the numerical simulation results can be observed in labs. It is noted that when the injected velocity is below approximately 4.17 × 10^{−6} m/s, the values of PVBT increase with increasing temperature; when the injected velocity is above approximately 4.17 × 10^{−5} m/s, the values of PVBT decrease with increasing temperature, which means that the former is a masstransfer controlled regime and the latter is a kinetically controlled regime. The porosity profiles at breakthrough in the optimal injected velocity are given in Figure 15 for the three different temperatures. It can be seen in the figure that with the increase of the temperature, the diameter of the wormhole also increases, which matches the observation in Figure 7 of [32]. This also explains why the minimal PVBT value increases with increasing temperature. In fact, the transferring efficiency of HCl decreases due to increased acid consumption on the walls of the wormhole, which brings about the phenomenon above.
Fig. 14 Acid efficiency curves of different temperatures in isothermal conditions. 
Fig. 15 Porosity profiles at breakthrough in the optimal injected velocity for three different temperatures. (a) 295 K, u_{x} = 2.67 × 10^{−6} m/s. (b) 323 K, u_{x} = 9.17 × 10^{−6} m/s. (c) 353 K, u_{x} = 4.17 × 10^{−5} m/s. 
Values of PVBT. The first row represents injected velocities, and their unit is m/s.
6.2 Nonisothermal conditions
Isothermal conditions are common in labs. However, in field cases, the injected acid temperature is often different from the initial matrix temperature. Due to the geothermal factor, the initial matrix temperature may be higher than the injected acid temperature. Thus, to expand the simulations from labs to fields, nonisothermal conditions are considered. In the experiment, four cases with different combinations of the injected acid temperature and initial matrix temperature are simulated, with an injected velocity of 9.17 × 10^{−6} m/s, and the values of PVBT are shown in Table 6. The other experimental parameters are the same as the isothermal conditions. All the results are from the finetimestep simulations. It can be seen in Table 6 that when the injected acid temperature is fixed at 323 K, the changed initial matrix temperatures will not change the values of PVBT. For three different initial matrix temperatures: 295 K, 323 K, and 353 K, the values of PVBT are the same as 4.362. However, when the injected acid temperature is changed, even though the initial matrix temperature is unchanged, the values of PVBT will be changed, which is evident from the first two rows in Table 6. From the two rows, it can be further seen that the two PVBT values are nearly the same as the corresponding values in the isothermal conditions, which also demonstrates that the injected acid temperature, instead of the initial matrix temperature, has an effect on the PVBT value. Thus, the injected acid temperature governs the PVBT value and can be a design parameter in matrix acidization, which can also be concluded from [31]. The porosity profiles at breakthrough for three different initial matrix temperatures are given in Figure 16, where the same kind of porosity profiles can be seen clearly. Thus, Figure 16 demonstrates that the temperature of the initial matrix is not a key factor affecting matrix acidization once again.
Fig. 16 Porosity profiles at breakthrough for three different initial matrix temperatures. The injected acid temperature is 323 K, and its injected velocity is 9.17 × 10^{−6} m/s. (a) 295 K. (b) 323 K (c) 353 K. 
Values of PVBT for nonisothermal conditions.
To learn the reason why the injected acid temperature has such a significant effect on matrix acidization, the change in the average matrix temperature with time is investigated. The initial matrix temperature is set as 323 K, and the injected acid temperatures are 295 K and 353 K, respectively, which are also the cases represented by the first two rows in Figure 6. The history of the average matrix temperature from the beginning to breakthrough is shown in Figure 17. It can be seen in the figure that the matrix temperature becomes more or less the same as the injected acid temperature immediately after acidization begins and continues to be similar until breakthrough. This explains why the initial matrix temperature has almost no effect on matrix acidization.
Fig. 17 Average matrix temperature against time. The solid curve represents the case where the temperature of the injected acid is 295 K, and the dashed curve represents the case where the temperature of the injected acid is 353 K. The initial matrix temperature is 323 K. 
After the discussion above, it is concluded that the thermal DBF framework can simulate reasonable numerical results, which are verified by the numerical and chemical experiments of other works. This indicates that the thermal DBF framework can be an effective tool in the field of matrix acidization.
7 Performance evaluation
The performance of the 2D parallel code was evaluated in [19], and this work tries to evaluate the performance of the 3D parallel code. The test is performed on the 3D grid used above. Moreover, the number of iterations is set to 100 to save supercomputing resources. Meanwhile, since the sparsity pattern of the coefficient matrix of the linear system of the energy conservation equation is the same as that of the concentration balance equation, evaluating both of the linear systems is a redundancy. Therefore, only the linear system of the concentration balance equation is evaluated, which means the performance of the 3D parallel code of the improved DBF framework is evaluated in this section. This evaluation result can foresee the performance of the 3D parallel code of the thermal DBF framework.
The experiment of the improved DBF framework in the verification section is performed again, and the performance results are shown in Figure 18. It can be seen in the figure that the solver time takes up most of the run time (>98%), which is reasonable. Except for the solver code, the other parts of the code only perform simple operations such as filling arrays and communicating data among processors, which will not cost much time. Moreover, the run time decreases with the increase in the number of processors, which means a certain level speedup can be achieved. However, when the number of processors increases to 144, further speedup seems impossible. The similar scalability of the MUMPS solver can be seen in [47].
Fig. 18 Performance of the 3D parallel code. The numbers above the bars represent the run time, which is the sum of the solver time and the other time. It can be seen that for each bar, the percentage of solver time over run time >98%. 
8 Conclusion and future work
Since Wu et al. [19] contributed the DBF framework to the field of matrix acidization, improvements to the framework have been ongoing. This work is such as endeavor and tries to correct a defect in the momentum conservation equation of the DBF framework and maintain the momentum conservation equation when the porosity is changed. Furthermore, by introducing a direct solver called MUMPS, the pseudo parameter ε in the mass conservation equation can be deleted, which keeps the incompressible attribute of the acid flow in matrix acidization and thus makes the framework more reasonable. In addition, the simulation flowchart is also changed in this work, which is another correction to the DBF framework. After these revisions, the new framework can be called the improved DBF framework for short. The improved DBF framework is realized by 2D and 3D parallel codes with the help of MPI and FORTRAN 90 and verified by comparison with a series of previous works. It is emphasized that the 3D simulation results of the improved DBF framework are given for the first time in this work. The improved DBF framework can simulate similar numerical results with [17] and [45], which demonstrates its reliability.
The correctness of the improved DBF framework makes it feasible to develop a thermal DBF framework based on this verification. In addition to the mass conservation law and momentum conservation law, which are included in the improved DBF framework, the thermal DBF framework also considers the energy conservation law and thus introduces the energy balance equation to the improved DBF framework. Verification to the thermal DBF framework is done under isothermal conditions and nonisothermal conditions, and the numerical simulation results match the conclusions from other chemical and numerical experiments such as [31] and [32]. Therefore, the thermal DBF framework is reasonable and trustable.
Since the accuracy of matrix acidization simulation is highly dependent on the size of the grid, very fine grids are required for trustable results, which brings about the need to develop parallel codes to finish simulations in reasonable time. However, parallelizing the improved DBF framework and thermal DBF framework is not an easy task due to their complex equation systems. With the help of MPI and FORTRAN 90 and the experimenting field approach, this work overcomes these difficulties and develops scalable parallel codes, which is another large contribution to the field of matrix acidization.
With the reliable improved DBF framework and thermal DBF framework, a series of numerical investigations on matrix acidization can be carried out in the future, and more reasonable results are expected.
Acknowledgments
For computer time, this research used the resources of the Supercomputing Laboratory at King Abdullah University of Science & Technology (KAUST) in Thuwal, Saudi Arabia.
References
 Liu P., Yao J., Couples G.D., Ma J., Huang Z., Sun H. (2017) Modeling and simulation of wormhole formation during acidization of fractured carbonate rocks, J. Pet. Sci. Eng. 154, 284–301. [Google Scholar]
 Ma G., Yun C., Yan J., Huidong W. (2018) Modelling temperatureinfluenced acidizing process in fractured carbonate rocks, Int. J. Rock Mech. Min. Sci. 105, 73–84. [CrossRef] [Google Scholar]
 Szymczak P., Kwiatkowski K., Jarosinskí M., Kwiatkowski T., Osselin F. (2019) Wormhole formation during acidizing of calcitecemented fractures in gasbearing shales, SPE J. 24, 795–810. [CrossRef] [Google Scholar]
 Zhao C., Hobbs B.E., Ord A., Hornby P., Peng S. (2008) Effect of reactive surface areas associated with different particle shapes on chemicaldissolution front instability in fluidsaturated porous rocks, Transp. Porous Media 73, 75–94. [Google Scholar]
 Zhao C., Hobbs B.E., Ord A., Peng S. (2010) Effects of mineral dissolution ratios on chemicaldissolution front instability in fluidsaturated porous media, Transp. Porous Media 82, 317–335. [Google Scholar]
 Zhao C., Hobbs B.E., Ord A. (2010) Theoretical analyses of the effects of solute dispersion on chemicaldissolution front instability in fluidsaturated porous media, Transp. Porous Media 84, 629–653. [Google Scholar]
 Buijse M.A. (2000) Understanding wormholing mechanisms can improve acid treatments in carbonate formations, SPE Prod. Facil. 15, 168–175. [CrossRef] [Google Scholar]
 Hoefnger M.L., Fogler H.S. (1988) Pore evolution and channel formation during flow and reaction in porous media, AIChE J. 34, 45–54. [Google Scholar]
 Fredd C.N., Fogler H.S. (1998) Influence of transport and reaction on wormhole formation in carbonate porous media, AIChE J. 44, 1933–1949. [Google Scholar]
 Daccord G., Touboul E., Lenormand R. (1989) Carbonate acidizing: Toward a quantitative model of the wormholing phenomenon, SPE Prod. Eng. 4, 63–68. [CrossRef] [Google Scholar]
 Daccord G., Lenormand R., Lietard O. (1993) Chemical dissolution of a porous medium by a reactive fluid – I. Model for the “wormholing” phenomenon, Chem. Eng. Sci. 48, 1, 169–178. [Google Scholar]
 Daccord G., Lietard O., Lenormand R. (1993) Chemical dissolution of a porous medium by a reactive fluid – II. Convection vs. reaction, behavior diagram. Chem. Eng. Sci. 48, 1, 179–186. [Google Scholar]
 Liu X., Ormond A., Bartko K., Ortoleva P. (1997) A geochemical reactiontransport simulator for matrix acidizing analysis and design, J. Pet. Sci. Eng. 17, 181–197. [Google Scholar]
 Golfier F., Zarcone C., Bazin B., Lenormand R., Lasseux D., Quintard M. (2002) On the ability of a Darcyscale model to capture wormhole formation during the dissolution of a porous medium, J. Fluid. Mech. 457, 213–254. [Google Scholar]
 Panga M.K.R., Ziauddin M., Balakotaiah V. (2005) Twoscale continuum model for simulation of wormholes in carbonate acidization, AIChE J. 51, 3231–3248. [Google Scholar]
 Maheshwari P., Balakotaiah V. (2013) Comparison of carbonate HCl acidizing experiments with 3D simulations, SPE Prod. Oper. 28, 4, 402–413. [Google Scholar]
 Akanni O.O., NasrElDin H.A., Gusain D. (2017) A computational NavierStokes fluiddynamicssimulation study of wormhole propagation in carbonatematrix acidizing and analysis of factors influencing the dissolution process, SPE J. 22, 6, 2049. [CrossRef] [Google Scholar]
 De Oliveira T.J.L., De Melo A.R., Oliveira J.A.A., Pereira A.Z.I. (2012) Numerical simulation of the acidizing process and PVBT extraction methodology including porosity/permeability and mineralogy heterogeneity, in: SPE International Symposium and Exhibition on Formation Damage Control, Society of Petroleum Engineers. [Google Scholar]
 Wu Y., Salama A., Sun S. (2015) Parallel simulation of wormhole propagation with the Darcy–Brinkman–Forchheimer framework, Comput. Geotech. 69, 564–577. [Google Scholar]
 Wu Y. (2015) Parallel reservoir simulations with sparse grid techniques and applications to wormhole propagation. Diss. [Google Scholar]
 Kou J., Sun S., Wu Y. (2016) Mixed finite elementbased fully conservative methods for simulating wormhole propagation, Comput. Meth. Appl. Mech. Eng. 298, 279–302. [CrossRef] [Google Scholar]
 Kou J., Sun S., Wu Y. (2019) A semianalytic porosity evolution scheme for simulating wormhole propagation with the Darcy–Brinkman–Forchheimer model, J. Comput. Appl. Math. 348, 401–420. [Google Scholar]
 Li X., Rui H. (2019) Superconvergence of a fully conservative finite difference method on nonuniform staggered grids for simulating wormhole propagation with the Darcy–Brinkman–Forchheimer framework, J. Fluid Mech. 872, 438–471. [Google Scholar]
 Li X., Rui H., Chen S. (2019) A fully conservative blockcentered finite difference method for simulating DarcyForchheimer compressible wormhole propagation, Numer. Algorithms 82, 2, 451–478. [Google Scholar]
 Li X., Rui H. (2020) A fully conservative blockcentered finite difference method for DarcyForchheimer incompressible miscible displacement problem, Numer. Meth. Partial Differ. Equ. 36, 1, 66–85. [CrossRef] [Google Scholar]
 Wang Y., Sun S., Gong L., Yu B. (2018) A globally massconservative method for dualcontinuum gas reservoir simulation, J. Nat. Gas Sci. Eng. 53, 301–316. [Google Scholar]
 Falgout R.D., Yang U.M. (2002) hypre: A library of high performance preconditioners, Springer, Berlin, Heidelberg. [Google Scholar]
 Amestoy P.R., Duff I.S., L’Excellent J.Y., KosterJ. (2001) A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM J. Mat. Anal. Appl. 23, 1, 15–41. [CrossRef] [MathSciNet] [Google Scholar]
 Amestoy P.R., Guermouche A., Lexcellent J.Y., Pralet S. (2006) Hybrid scheduling for the parallel solution of linear systems, Parallel Comput. 32, 2, 136–156. [Google Scholar]
 Li Y., Liao Y., Zhao J., Peng Y., Pu X. (2017) Simulation and analysis of wormhole formation in carbonate rocks considering heat transmission process, J. Nat. Gas Sci. Eng. 42, 120–132. [Google Scholar]
 Kalia N., Glasbergen G. (2010) Fluid temperature as a design parameter in carbonate matrix acidizing, in: SPE Production and Operations Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Fredd C.N., Fogler H.S. (1999) Optimum conditions for wormhole formation in carbonate porous media: Influence of transport and reaction, SPE J. 4, 3, 196–205. [CrossRef] [Google Scholar]
 Wang Y., Sun S. (2016) Direct calculation of permeability by highaccurate finite difference and numerical integration methods, Commun. Comput. Phys. 20, 2, 405–440. [Google Scholar]
 Wang Y. (2019) Reynolds stress model for viscoelastic dragreducing flow induced by polymer solution, Polymers 11, 10, 1659. [Google Scholar]
 Whitaker S. (1996) The Forchheimer equation: A theoretical development, Transp. Porous Media 25, 1, 27–61. [Google Scholar]
 Aavatsmark I. (2002) An introduction to multipoint flux approximations for quadrilateral grids, Comput. Geosci. 6, 3, 405–432. [Google Scholar]
 Šekutkovski B., Kostić I., Stefanović Z., Simonović A., Kostić O. (2015) A hybrid RANSLES method with compressible komegaSSTSAS turbulence model for high Reynolds number flow applications/Hibridna RANSLES metoda s kompresibilnim komegaSSTSAS turbulentnim modelom namjenjena analizi strujanja pri velikim Reynoldsovim brojevima, Tehnicki VjesnikTechnical Gazette 22, 5, 1237–1246. [Google Scholar]
 Sun S., Salama A., ElAmin M.F. (2012) An equationtype approach for the numerical solution of the partial differential equations governing transport phenomena in porous media, in: The International Conference on Computational Science, ICCS, June 4–6, Omaha, NE. [Google Scholar]
 Salama A., Sun S., El Amin M.F. (2013) A multipoint flux approximation of the steady state heat conduction equation in anisotropic media, ASME J. Heat Trans. 135, 1–6. [CrossRef] [Google Scholar]
 Salama A., Li W., Sun S. (2013) Finite volume approximation of the three dimensional flow equation in axisymmetric, heterogeneous porous media based on local analytical solution, J. Hydrol. 501, 45–55. [CrossRef] [Google Scholar]
 Salama A., Sun S., Wheeler M. (2014) Solving global problem by considering multitude of local problems: Application to flow in anisotropic porous media using the multipoint flux approximation, J. Comput. Appl. Math. 267, 117–130. [Google Scholar]
 Salama A., Sun S., El Amin M.F. (2015) Investigation of thermal energy transport from an anisotropic central heating element to the adjacent channels: A multipoint flux approximation, Ann. Nucl. Energy 76, 100–112. [Google Scholar]
 Davis T.A. (2004) Algorithm 832: UMFPACK V4. 3 – an unsymmetricpattern multifrontal method, ACM Trans. Math. Softw. 30, 196–199. [Google Scholar]
 Hadri B., Kortas S., Feki S., Khurram R., Newby G. (2015) Overview of the KAUST’s Cray X40 system–Shaheen II, in: Proceedings of the 2015 Cray User Group. [Google Scholar]
 Ku H.C., Hirsh R.S., Taylor T.D. (1987) A pseudospectral method for solution of the threedimensional incompressible NavierStokes equations, J. Comput. Phys. 70, 2, 439–462. [Google Scholar]
 Kalia N., Balakotaiah V. (2009) Effect of medium heterogeneities on reactive dissolution of carbonates, Chem. Eng. Sci. 64, 2, 376–390. [Google Scholar]
 Raju M.P. (2009) Parallel computation of finite element NavierStokes codes using MUMPS solver, Int. J. Comput. Sci. Issu. 4, 2, 20–24. [Google Scholar]
All Tables
Values of PVBT. The first row represents injected velocities, and their unit is m/s.
All Figures
Fig. 1 The interaction of the pore scale and Darcy scale. The upper arrow indicates that the pore scale will affect Darcy scale by the porescale variable ϕ. As long as ϕ is changed, the porescale variables used in the Darcy scale are changed, which brings about the changes of the main variables u, p, and C_{f} in the Darcy scale. The lower arrow indicates that the Darcy scale will affect the pore scale by the Darcyscale variables C_{f} and u, since the two variables will affect the value of ϕ directly. Once ϕ is changed, all the other variables in the pore scale will be changed accordingly. 

In the text 
Fig. 2 The flowchart of the improved DBF framework. 

In the text 
Fig. 3 The flowchart of the thermal DBF framework. 

In the text 
Fig. 4 The interaction of the improved DBF framework and the energy conservation equation. 

In the text 
Fig. 5 The relationship among the Darcy framework, the DBF framework, the improved DBF framework, and the thermal DBF framework. 

In the text 
Fig. 6 Staggered grid. 

In the text 
Fig. 7 3D sheardriven cavity flow configuration and coordinate system [45]. 

In the text 
Fig. 8 Velocity profile of the xdirection component (u) on the vertical centerline of the plane z = 0.5. 

In the text 
Fig. 9 Velocity profile of the ydirection component (v) on the horizontal centerline of the plane z = 0.5. 

In the text 
Fig. 10 Acidefficiency curve of the 2D linear flows. The numbers beside the points represent the values of PVBT. 

In the text 
Fig. 11 Porosity profiles at breakthrough for five different injected velocities of the 2D linear flows. (a) u_{x} = 4.17 × 10^{−7} m/s, face dissolution. (b) u_{x} = 1.67 × 10^{−6} m/s, conical wormhole. (c) u_{x} = 4.17 × 10^{−6} m/s, dominant wormhole. (d) u_{x} = 7.17 × 10^{−6} m/s, ramified wormhole. (e) u_{x} = 1.67 × 10^{−5} m/s, uniform dissolution. 

In the text 
Fig. 12 Acidefficiency curves of the 3D linear flows. The numbers beside the points represent the values of PVBT. 

In the text 
Fig. 13 Porosity isosurfaces at breakthrough for five different injected velocities of the 3D linear flows. (a) u_{z} = 3.04 × 10^{−7} m/s, face dissolution. (b) u_{z} = 7.04 × 10^{−7} m/s, conical wormhole. (c) u_{z} = 1.04 × 10^{−6} m/s, dominant wormhole. (d) u_{z} = 3.04 × 10^{−6} m/s, ramified wormhole. (e) u_{z} = 1.04 × 10^{−5} m/s, uniform dissolution. 

In the text 
Fig. 14 Acid efficiency curves of different temperatures in isothermal conditions. 

In the text 
Fig. 15 Porosity profiles at breakthrough in the optimal injected velocity for three different temperatures. (a) 295 K, u_{x} = 2.67 × 10^{−6} m/s. (b) 323 K, u_{x} = 9.17 × 10^{−6} m/s. (c) 353 K, u_{x} = 4.17 × 10^{−5} m/s. 

In the text 
Fig. 16 Porosity profiles at breakthrough for three different initial matrix temperatures. The injected acid temperature is 323 K, and its injected velocity is 9.17 × 10^{−6} m/s. (a) 295 K. (b) 323 K (c) 353 K. 

In the text 
Fig. 17 Average matrix temperature against time. The solid curve represents the case where the temperature of the injected acid is 295 K, and the dashed curve represents the case where the temperature of the injected acid is 353 K. The initial matrix temperature is 323 K. 

In the text 
Fig. 18 Performance of the 3D parallel code. The numbers above the bars represent the run time, which is the sum of the solver time and the other time. It can be seen that for each bar, the percentage of solver time over run time >98%. 

In the text 