Analytical Solutions of Carbonate Acidizing in Radial Flow

The flow of reactive fluids into porous media, a phenomenon known as reactive infiltration, is important in natural and engineered systems. While most of the studies in this area cover theoretical and experimental analyses in linear acid flow, the present work concentrates on radial flow conditions from a wellbore in the field and on finding exact analytical solutions to moving boundary problems of the uniform dissolution front. Closed-form solutions are obtained for the transient convection–diffusion which allow the demarcation of the range of applicability of the quasi-static limit. The fluid velocity dependency of the diffusion–dispersion coefficient is also examined by comparing results from analytical solutions from constant and velocity-dependent coefficients. These contributions form the basis for linear stability analyses to describe acid fingering encountered in reservoir stimulation.


Introduction
Reactive infiltration is the phenomenon where fluids enter the free space of a porous solid, while their components react with the solid matrix and result in changes to its porosity and permeability. Reactive infiltration arises both naturally and in engineered systems. Main applications regarding of the latter are acidizing of wells in hydrocarbon reservoirs (Economides et al. 2013;Fredd and Fogler 1998), and chemical reactions during CO 2 subsurface storage (Liu et al. 2011). Most of the research focuses on analyzing either theoretically or experimentally the pattern formation in rocks created by the reactive infiltration instability, a self-organizing system where the interplay between convective and diffusive fluxes sculpts the shape of the dissolution profile in the rock (e.g., Ortoleva et al. 1987;Papamichos et al. 2020).

3
The present work focuses on a specific application of reactive infiltration that of carbonate acidizing which is a common technique to enhance hydrocarbon production in carbonate reservoirs. In this technique, an aqueous solution of hydrochloric acid is injected under pressure from the well into the near wellbore area of the formation. The calcite is dissolved by the injected acid, its porous network is altered with new pathways emerging, and its porosity and permeability are enhanced. Most works in reactive infiltration in general and in carbonate acidizing in particular consider linear flow of the injected acid (Chadam et al. 1986;Hinch and Bhatt 1990;Szymczak and Ladd 2011;Wangen 2013;Zhao et al. 2013). The present study considers radial flow and thus simulates the geometry in applications where reactive fluids are injected through a cylindrical well and the relevant geometry can be assumed to be axisymmetric. Radial flow acidizing has been studied theoretically by Grodzki and Szymczak (2019), and experimentally in Indiana limestone and a dolomite by McDuff et al. (2010), and in Mons chalk by Walle and Papamichos (2015).
This study can be considered as the first but essential step toward a linear stability analysis of reactive infiltration in radial flow conditions. Such an analysis requires a base solution which corresponds to a uniform, axisymmetric dissolution pattern around the well. In linear flow, the base solution corresponds to a planar dissolution front and it is simpler to obtain because the front travels at a constant velocity. This is not the case for radial flow where the base solution corresponds to the axisymmetric propagation of the dissolution front that travels with a decreasing velocity due to the diverging geometry. Grodzki and Szymczak (2019) obtained such a base solution but their derivation assumed first that the small acid capacity number leads to a quasi-static solution for the acid concentration profile, and second that the diffusion-dispersion coefficient is constant. They also neglected the finite radius of the wellbore and assumed point acid injection. In this study, we first address the validity of the assumption of quasi-static approximation by deriving the transient solution and comparing with the quasi-static one to obtain the limits of its applicability. The transient solution requires the assumption of point injection to obtain self-similar formulations that can be solved analytically. We then remove the simplifying assumptions of constant diffusion-dispersion coefficient and point injection and solve the quasi-static problem in a more general framework for more realistic predictions. The finite wellbore radius becomes especially important in linear stability analyses of reactive infiltration as it provides the necessary scale to upscale the problem in the field and lead to wavenumber selection in fingering.
Section 2 introduces the governing equations of reactive infiltration in polar coordinates, describes the associated geometry and introduces the important assumption of thin front approximation in carbonate acidizing which allows the formulation of various Stefan-type moving boundary problems that can be solved analytically in the direction of previous works (Grodzki and Szymczak 2019;Papamichos et al. 2020). Section 3 derives the transient solution of radial reactive infiltration assuming constant in space and time diffusion-dispersion coefficient and compares the results with those obtained at the quasi-static limit of small acid capacity number. The development of this more general closed-form solution allows the analytical validation of the approximate steadystate solution and the demarcation of its range of applicability. Section 4 adopts the small acid capacity number assumption and solves analytically for a finite radius wellbore the quasi-static problem of a velocity-dependent diffusion-dispersion coefficient. It also examines the effect of this dependence on the propagation velocity of the reaction front and the fluid flux and acid concentration profiles. Finally, Section 5 presents the conclusions.

Governing Equations
In a porous medium, the representative volume V consists of the volume of solids V s and the volume of voids V v . In the presence of a fully saturating fluid, the volume of fluid V f = V v . The porosity φ of the porous medium is defined as the ratio of the volume of voids V v to the total volume V. In acidizing, an aqueous acidic solution is injected radially from a wellbore of radius r 0 into the formation at a prescribed fluid flux u 0 [m/s], and prescribed concentration c 0 [mol/m 3 ], as shown in Fig. 1. For carbonates, the acidic solution is usually an aqueous hydrochloric acid (HCl) solution which reacts with the calcium carbonate (CaCO 3 ) according to the chemical reaction The porous solid has an initial porosity φ 0 . After acidizing, it is assumed that its porosity reaches its final, maximum porosity f ≤ 1 . The solids have molar density ρ s defined as the ratio of number of moles of solids per volume of solids V s . In carbonate acidizing, the porous solid consists almost entirely of calcite (CaCO 3 ) with ρ s = 270 mol/L. The acid concentration c [mol/m 3 ] is defined as where n acid is the number of moles of acid per fluid volume V f .
The injected acid solution dissolves the porous solid forming areas where the rock has reached its maximum porosity. The medium is divided into two domains, the upstream where the rock has reached its final porosity value φ f , and the downstream where calcite has not been fully dissolved (Fig. 1). The conservation laws of porous media are applied to both domains to extract the governing equations. Conservation of fluid mass leads to the differential equation for the Darcy fluid flux u [m/s] In the upstream domain, the rock has reached its final porosity and therefore the porosity is everywhere φ f . In the downstream domain, the evolution of porosity is described by the differential equation where ν [−] is the stoichiometric coefficient for the acid in the chemical reaction Eq. (1), and Σ(c) is the reaction term that describes the rate of dissolution per unit volume and time [mol/(m 3 s)]. The reaction term is assumed to be linearly dependent on the acid concentration and on the difference between the maximum and the current porosity (e.g., Zhao et al. 2013) where k is the reaction rate constant [m/s] and s the available surface for reaction per unit volume of material [m −1 ]. In this way, the reaction stops either when the porosity reaches its maximum value or when no acid is present in the fluid. Conservation of acid molar mass yields a convection-diffusion reaction equation (Bear 1972) It contains a convective term, a diffusion-dispersion term, and the reaction term Σ(c). The diffusion-dispersion term involves the diffusion-dispersion coefficient D [m 2 /s].
The governing equations can be expressed in the polar coordinate system (r, θ, z) associated with the wellbore where the z-axis coincides with the axis of the wellbore. It is assumed that there is no vertical dependency of the solution because gravity effects are negligible compared to the convective and diffusive forces developing during acid injection, and because the injection takes place simultaneously and uniformly on the whole height of the formation. Additionally, axisymmetric conditions are assumed because we are after the uniform dissolution pattern of the base solution. Under these considerations, and considering that the porosity upstream is constant and equal to the final porosity φ f , the governing equations upstream can be written as and downstream as where the first equation expresses the conservation of fluid mass, the second the conservation of acid molar mass, and the third in the downstream domain describes the evolution of porosity due to acid reaction and dissolution. The unknowns of the problem are the fluid flux u, the acid molar concentration c, and the porosity φ downstream.
(4) s t = Σ(c) For analytical solutions, the important assumption of the thin front approximation is introduced to transform the problem to a moving boundary problem with a sharp front. Following Ladd and Szymczak (2017), the thin front approximation states that if the reaction rate k is fast enough, the penetration distance of the acid into the unreacted domain is negligible leading to a sharp dissolution front. In other words, the reaction takes place only in a thin front which separates the fully dissolved from the intact rock. This assumption is validated by experimental data that show that in carbonates no transition zone with respect to the change in rock strength is observed around wormholes in acidizing experiments indicating localized thin dissolution fronts (Papamichos et al. 2020). Mathematically, it means that the porosity is not only constant upstream but also downstream where it is equal to the initial porosity φ 0 . The presence of a thin front turns the model into a moving boundary problem analogous to the classical Stefan problem that describes phase-change processes (Stefan 1889). In acidizing, one phase is the fully dissolved rock at final porosity φ f and the other phase is the intact rock at initial porosity φ 0 . Thus, the governing equations upstream remain the same (Eq.(7)), while downstream, only the fluid flux needs to be solved since the porosity and concentration become known and given as The boundary conditions for the solution of the differential equations Eqs. (7) and (9) are the acid concentration c 0 and the fluid flux u 0 at the wellbore radius r 0 , i.e., At the interface r = R(t) between the fully dissolved rock and the intact rock, i.e., between the upstream and downstream domains, fluid volume conservation and continuity of acid concentration holds, i.e., The derivation of the condition for the fluid flux is given in Appendix A. The second term in the right-hand side of this equation is a storage term that expresses the variation of fluid flux in the porous medium due to the change in porosity across the interface. A minus or plus superscript indicates the upstream or downstream side of the interface, respectively.
The initial conditions are everywhere, the acid concentration is zero, and the interface between the two domains is at its initial position Figure 2 shows a schematic of a cross section of the acidizing problem in a wellbore indicating the boundary and continuity conditions at the reaction front under the thin front approximation.
The mathematical closure of the Stefan problem requires an additional equation at the interface. This equation in the classical Stefan problem (Stefan 1889) is called the Stefan condition, and it is extracted through the law of conservation of energy while in this case expresses the conservation of molar mass of calcite. For reactive infiltration in polar coordinates and for axisymmetric front movement, it is written as The Stefan condition provides the solution for the propagation velocity of the dissolution front by solving for dR∕dt . Its derivation is given in Appendix B.
Despite the fluid flux upstream being coupled to the acid concentration in the convection-diffusion equation in Eq. (7), uncoupling of the equations is possible since the fluid mass conservation equation contains only the fluid flux as dependent variable. Downstream the concentration is zero. Thus, the mass conservation equation can be solved explicitly both upstream and downstream as where upstream the boundary condition Eq. (11) for the fluid flux was used to find the integration constant, while downstream A remains to be solved from the continuity condition Eq. (12) for the fluid flux.

Solution for Transient Convection-Diffusion
This section compares the analytical solution for the transient problem to the solution for the quasi-static limit approximation. In this way, the small acid capacity approximation of the quasi-static problem can be studied and the limits to its applicability can be identified. The transient problem can be solved analytically under the assumption of constant diffusion-dispersion coefficient D = D 0 and point injection, i.e., prescribed acid concentration c 0 at r = 0, and not at the wellbore radius r = r 0 . These assumptions were also adopted by Grodzki and Szymczak (2019) for their quasi-static solution. The convection-diffusion equation for the concentration in Eq. (7) simplifies for constant D to Wellbore cross section with the reacted (upstream) and intact (downstream) domains, the boundary conditions at the wellbore radius r = r 0 and the continuity conditions at the reaction front r = R(t) which is solved under the interface conditions Eqs. (12) and (14)  where a c is a dimensionless number called acid capacity number that is defined as The acid capacity number is an important parameter that expresses the potency of the aqueous solution to dissolve the mineral. With increasing a c , less volume of acid solution is necessary to dissolve the same volume of rock. Substituting u from Eq.(15) into Eq. (17), the convection-diffusion equation becomes Using the Boltzmann's transformation (Boltzmann 1894) the partial differential Eq. (22) transforms into an ordinary differential equation with ξ the independent variable Solving the ordinary differential equation and imposing the boundary Eq. (18) and continuity condition Εq. (12) Equation (25) also contains the dimensionless number Pe defined as It is a Peclet-type control parameter which expresses the interplay between the convective and diffusive forces in reactive infiltration. Large Pe implies convection dominated flow as compared to the diffusive-dispersive. For Pe = 0, the convection-diffusion equation turns into the classical diffusion equation for a cylinder (Carlsaw and Jaeger 1959), since the only force present to change the concentration field is molecular diffusion.
The constant λ can be calculated by applying the solution Eq. (25)  The solution can be illustrated by considering a typical carbonate acidizing scenario as the base case. Table 1 lists the boundary values of the problem as well as the geometrical and material parameters obtained from laboratory experiments in chalk (Papamichos et al. 2020). The value of acid concentration c 0 corresponds to 15% v/v HCl acid solution. The table also lists parameters like the acid capacity number a c and the Peclet number Pe that are calculated from the input values according to their definitions Eqs. (21) and (29), respectively. The value for the constant diffusion-dispersion coefficient D 0 was estimated in previous analyses on Mons chalk acidizing (Papamichos et al. 2020). The values for reaction rate constant k and the available reactive surface area s were taken from the literature (Dreybrodt, 1996;Noiriel et al., 2009). For this base case, the nonlinear algebraic Eq. (30) solved with the standard Newton-Raphson method gives m = 0.07483 and λ = 5.471 × 10 -4 ms −1∕2 . Figure 3 compares the radial position of the dissolution front and the front velocity over time for the transient and the quasi-static case. For the typical values of acidizing in Table 1, the difference in the results between the transient quasi-static solutions are less 2u 0 r 0 than 0.7% with the transient front being slower than the quasi-static. The front position is parabolic with time, while the front velocity decreases with the inverse square root of time as one expects due to the ever-increasing area to be dissolved as the distance from the wellbore center increases. In the quasi-static limit, the front position is calculated explicitly as and thus, the corresponding quasi-static coefficient λ qs for the small acid capacity limit is Substituting in Eq.(34) the values of parameters from Table 1, we calculate λ qs = 5.505 × 10 -4 ms −1∕2 which compares well with the transient solution. Figure 4 plots radial profiles of the acid concentration c c 0 and the fluid flux u u 0 at t = 24 h for the transient and the quasi-static cases. The quasi-static solution can be found in Grodzki and Szymczak (2019). The location of the reaction front is also indicated. At the front, there is a step in the fluid flux due to the storage term in the continuity equation for the fluid flux. This step is here relatively small to be visible in the plot. There is as a change in the slope of the flux before and after the front. Initial acid concentration, c 0 (mol/m 3 ) 4.5 × 10 3 Fluid flux, u 0 (m/s) 1 × 10 -4 Stoichiometric coefficient, ν (−) 2 Chalk molar density, ρ s (mol/m 3 ) 2.7 × 10 5 Initial porosity, φ 0 (−) 0.43 Maximum porosity, φ f (−) 0.98 (≈ 1) Diffusion-dispersion coefficient, D 0 (m 2 /s) 1 × 10 -6 Reaction rate constant, k (m/s) 2 × 10 -7 Reaction surface per unit volume, s (m −1 ) 5 × 10 3 Acid capacity number, a c (−) 8.333 × 10 -3 Peclet number, Pe (−) 5.102 To further investigate the range of validity of the quasi-static approximation, the acid concentration is increased to c 0 = 36,180 mol/m 3 which corresponds to the theoretical maximum limit of 100% HCl acid. Figure 6 plots radial profiles of c c 0 at t = 12, 24, and 48 h for u 0 = 10 -5 and 10 -4 m/s. For the higher acid concentration and the higher fluid flux (Fig. 6b), the transient solution front is noticeably lagging the quasi-static.
Another way to look at this is to compare the difference in the values of the parameter λ which controls the position and velocity of the reaction front. In fact, both are  Figure 7 plots the difference λ qsλ between the quasistatic λ qs and the transient λ as a function of fluid flux u 0 for three values of c 0 = 4500, 9000, and 36,180 mol/m 3 which correspond to acid capacity number a c = 0.00833, 0.0167, and 0.067, or HCl solutions 15, 27.5, and 100%. By normalizing the results with a c λ, they collapse into a single curve which now defines the range of applicability of the quasi-static approximation. The difference in the results tends to zero with decreasing fluid flux and for fluid flux u 0 ≤ 10 -6 or 10 -5 m/s the difference is not significant. The difference also tends to an asymptote less than one with increasing fluid flux which means that the percentage error from the quasi-static limit is bounded by approximate 100 × a c %. This corresponds to a maximum error of ca. 0.8% and 7% for the bases case a c = 0.00833 or the maximum a c = 0.067, respectively. Figure 7 shows also that the For Peclet number Pe ≤ 1 , the sum of the three terms in the brackets is negative. Considering that the exponentials, the incomplete lower gamma function, and the remaining terms are always positive, the derivative of the function is always negative. This monotonicity guarantees the uniqueness of the root of f(m) = 0. The assumption that Pe ≤ 1 is quite safe and reasonable at least for applications of carbonate acidizing as demonstrated for the typical acidizing problem in Table 1 where Pe = 0.22. The existence and uniqueness of the solution for a root m is shown graphically in Fig. 8 where the function f(m) is plotted as a function of m for the input values in Table 1.

Solution for Velocity-Dependent Diffusion-Dispersion Coefficient
The diffusion-dispersion coefficient D in the convection-diffusion reaction Eq. (6) is a physical parameter that describes both the effects of molecular diffusion and mechanical dispersion in the flow of a species in a porous medium. As already mentioned, Grodzki and Szymczak (2019) assumed for simplicity that this coefficient is constant. However, as comprehensively explained by Bear (1972), this coefficient depends on the interstitial fluid velocity υ which relates to the fluid flux as  Table 1 demonstrating the existence of a unique, positive root for equation f(m) = 0 According to Phillip (1994), D is related to the fluid velocity υ through a power law where D 0 is the diffusion-dispersion coefficient at υ = υ 0 , and p is an exponent between 1 ≤ p ≤ 2 . Craig and Heidlauf (2009) stated that D is the product of dispersivity and fluid velocity, essentially assuming a linear relationship between D and υ or u, i.e., p = 1, which we also adopt. In linear acid flow, the velocity dependency of D does not play a role because the fluid velocity remains constant in space and time. However, in radial flow, the fluid velocity decreases as we move away from the wellbore due to the diverging flow geometry. The effect may be negligible at low fluid flow velocities where diffusive fluxes dominate over the dispersive. However, in general, the assumption of constant diffusion-dispersion coefficient is not valid.
The solution for a velocity-dependent diffusion-dispersion coefficient can help illuminate the effect of this important assumption made by previous researchers who in addition assumed point acid injection. The analytical solution for velocity-dependent diffusion-dispersion coefficient can be obtained in the quasi-static limit where the assumption of small acid capacity number removes the transient nature of the acidizing problem and solves the concentration field in a quasi-static manner. To reach the quasi-static problem and to reveal important dimensionless control parameters, we nondimensionalize the model by introducing the characteristic dimensions of the system under consideration. We assume a characteristic fluid flux that relates to the characteristic diffusion-dispersion coefficient D 0 (Eq. (38)). This flux is taken as the flux u 0 at wellbore radius r 0 . At this location and for t > 0, the porosity φ = φ f , and thus from Eqs. (37) and (38), the diffusion-dispersion coefficient D can be written as We use the following scaling parameters to scale the time t and the radial coordinate r, respectively With the scaling parameters in Eq. (40) and c 0 and u 0 , we obtain the non-dimensional variables In terms of these non-dimensional variables, and by using Eq. (39) for D, the convection-diffusion equation for the concentration in Eq. (7) upstream where φ = φ f becomes The dimensionless parameter H in the second of Eq. (42) is defined as   (54) shows that the drop of the fluid flux across the interface is a decreasing function of front position or time. The limit value for large time (or front position) is zero.
A difference with the case of point injection and constant diffusivity is that the front position R(t) cannot be solved explicitly (cf. Equation (33)) but the nonlinear algebraic equations (52) must be solved at each point in time. However, Eq. (52) are explicit in time for a given front position. Also, the drop in the flux across the reaction front is not constant but depends on the position of the front. But the most important differences are in the front position and the concentration profiles. Figure 9 plots the front position R r 0 as a function of time for constant and velocitydependent D for the parameters in Table 1, and for injection fluid fluxes, u 0 = 0.001, 0.005, 0.02, and 0.1 mm/s, which correspond to Peclet numbers Pe = 0.051, 0.255, 1.02, and 5.10, respectively. The results show that the difference in the two cases increases as the fluid flux, and thus the convection component of acid transport decreases. The fluid flux decreases with distance from the well and thus in the case of velocity-dependent D, the front position moves slower than for constant D. With increasing fluid flux, convective transport dominates, and the importance of diffusivity diminishes. This is reflected, e.g., to Fig. 9 for u 0 = 0.1 mm/s where there is practically no difference in the results for the front position.
The difference in front position is controlled by the Peclet number Pe, such that for Pe > 3, there is practically little effect of diffusivity on the results. Figure 10 compares the front position for the cases where the acid concentration c 0 is applied at the wellbore radius r = r 0 versus at the wellbore center r = 0, as in the solution by Grodzki and Szymczak (2019). The results show a considerable delay in the front position when injection takes place at the center. Notice that in both cases, the fluid flux u 0 is prescribed at r = r 0 . Figure 11 plots the radial profiles of the concentration c c 0 and the fluid flux u u 0 for constant and velocity-dependent D, u 0 = 0.001 mm/s at time t = 48 h. The plot shows the position of the front in the two cases where the concentration c becomes zero and also shows the discontinuity of the fluid flux across the front which is expressed by Eq. (54). The fluid flux profile is independent of D upstream but downstream depends on D due to the difference in the discontinuity in the two cases across the front. This difference is not necessarily large and decreases as the front moves away from the wellbore. Figure 12 shows the propagation of the concentration u u 0 profile in time, i.e., for t = 48, 240, 480 h for constant and velocity-dependent D and injection fluxes u 0 = 0.001, 0.02, and 0.1 mm/s. The results show that besides the position of the front, there are differences in the concentration profiles as well, namely the velocity-dependent D gives sharper concentration fronts. This is because the concentration tends to propagate faster at smaller r r 0 where the flux and diffusivity D is larger. Even for u 0 = 0.1 mm/s where the front positions practically coincide, the concentration profile is markedly different (Fig. 12c). It can be shown that the profiles are always convex for velocity-dependent D, while for constant D, the profiles are concave for Pe < 0.5 and convex otherwise. Figure 13 shows the propagation of the concentration u u 0 profile in time, i.e., for t = 48, 240, 480 h for constant and velocity-dependent D and injection fluxes u 0 = 0.001, and 0.02 mm/s. For low fluid flux, the discontinuity drop of flux across the front is substantial and results in different fluid flux profiles. With increasing fluid flux, the drop decreases, and for Pe > 1, the effect appears to be minimal.

Conclusions
Analytical solutions to reactive infiltration problems under radial flow conditions were studied. A novel closed-form analytical solution was extracted for the initial transient convection-diffusion equation. Such an analytical solution is possible under the assumption of point acid injection at the wellbore center. This is because of the similarity condition necessary for obtaining an analytical solution. This similarity condition breaks down when an internal radius is introduced for the boundary condition for the concentration. On the other hand, for the fluid flow, there is no such limitation and the boundary condition for the fluid flux can be applied at the wellbore radius. The transient solution requires the solution of an explicit algebraic equation that needs to be solved numerically. This solution, however, is solved only once for a given problem and provides a material and geometry-based parameter that controls the location and velocity of the moving front. In fact, the location and velocity are proportional to this parameter. The comparison of the transient with the quasi-static solution in the literature shows that the transient front is always lagging the quasi-static front. For slow flow, i.e., u 0 ≤ 10 -6 or 10 -5 m/s, the quasi-static limit approximation is valid, and the problem can be solved without the transient term. The percent difference between the two solutions increases with increasing fluid flux but it is bounded by approximate 100 × a c %. This corresponds to a maximum error of ca. 0.8% and 7% for the base case of a 15% HCl solution and the theoretical maximum of 100% HCl, respectively. Despite the assumption that the acid concentration is prescribed at the wellbore center, this comparison shows the range of applicability of the quasi-static limit approximation for any reactive infiltration system. Within the quasi-static limit approximation, the dependence of the diffusion-dispersion coefficient D on the fluid velocity was investigated by deriving the analytical solutions for both constant and velocity-dependent D with acid and fluid flux prescribed at the wellbore radius. Previous solutions had assumed a constant diffusion-dispersion coefficient D and point injection. The constant D assumption can be justified in the case of linear systems where the fluid flux is constant except for a small step in the flux at the reaction front. However, for radial flow, the divergent nature of the flow due to the geometry gives rise to a parabolic decrease in fluid flux in space. The closed-form solution that was developed assumes a linear relation between the diffusivity D and the fluid flux. Such a relation is commonly considered although relations up to the second power of the flux are also reported. The closed-form solution for the dissolution front movement results in a nonlinear algebraic equation for the front position that needs to be solved at each point in time. This equation however is explicit in time at a given front position.
The results show that point injection at the wellbore center gives a significant delay in the position of the reaction front and it is not necessary for a obtaining an analytical solution in the quasi-static limit. In terms of front position, the velocity dependency of D becomes important at low fluid fluxes, i.e., Peclet number Pe < 3. The results showed a significant decrease in front propagation velocity that needs to be considered for realistic predictions. This is a reasonable result since the diffusive and dispersive fluxes that enhance rock dissolution decrease as the front propagates away from the wellbore center. This is because for a velocity-dependent D, the velocity and thus D decrease as the front propagates. At higher Pe, convective transport dominates, and the effect of diffusivity becomes less important for the propagation of the front. However, even then, the concentration profile is markedly different and sharper for the velocity-dependent D. The results for the diffusion-dispersion coefficient are expected to have significant ramifications with respect to the linear stability of the reactive infiltration problem since it gives a wavenumber selection that scales with this coefficient.
where the storage term Q storage expresses the rate of the fluid volume [m 3 /s] stored in the porous rock due to the increase in porosity during the front's movement from R to R + dR. The flow rates Q − and Q + for a unit height of wellbore are given as while the storage term Q storage is equal to the flow rate stored in the newly created voids of the dissolved rock in the green-shaded area in Fig. 14. The dissolved unit height volume can be calculated as where the second-order term of dR has been neglected. The rate of dissolved volume is equal to the rate of stored fluid and thus from Eq. (57), we obtain Substituting Eqs. (56) and (58) in Eq. (55), dividing both sides by 2πR and taking dR∕R ≪ 1 yields which is the interface condition for the fluid flux.

Appendix 2: Stefan condition-Molar mass conservation of acid across the reaction front
This appendix derives the molar mass conservation condition across the reaction front interface, the so-called Stefan condition, under the assumption of the thin front approximation. Figure 14 shows that during an infinitesimal time interval dt, the front moves from position R to position R + dR. The conservation of moles of acid during that time interval can be written as (56) where N − acid and N + acid are the fluxes of moles of acid [mol/(m 2 s)] entering and exiting the interface, respectively, and N consumed acid is the flux of moles of acid consumed from the chemical reaction. The thin front approximation assumes all moles of acid entering the interface are consumed to dissolve the rock to its final porosity and thus no moles of acid exit the interface to the unreacted domain, i.e., N + acid = 0. The molar flux N − acid entering the interface is calculated from Fick's law for a nonstationary fluid (Bear 1972) which states that the molar flux of a species in advective-dispersive transport is equal to − D∇c + cu . Using the continuity condition Eq. (12) for the concentrations across the interface and the axisymmetry of the problem, the rate of moles ̇n − acid [mol/s] entering the interface R(t) can be calculated as On the other hand, the rate of moles consumed to dissolve the rock from an initial porosity φ 0 to its final porosity φ f is according to the stoichiometry of the chemical reaction ν times the rate of moles of dissolved rock. Using the molar density ρ s of the rock and the dissolved volume from Eq. (57), the rate of rock moles dissolved and the rate of acid moles consumed per unit wellbore height are written as Equating the rate of acid moles entering the interface Eq. (61) to the rate of moles being consumed at the interface Eq. (62) and dividing both sides by 2πRφ f yields which is the Stefan condition at the interface.
where the definition Eqs. (26) is used for m. The left-hand side of Eq. (66) is a constant. Thus, the right-hand side must also be a constant. The Pe number is constant, so m must also be a constant and concurrently the ratio R(t) � √ t = λ is also a constant. Solving Eq.
(66) for the integration constant B and substituting it in Eq. (64) together with A from Eq. (65), the solution becomes The constant m in Eq. (67) is obtained from the Stefan condition Eq. (14). The partial derivative of c(r,t) with respect to r is given through the derivative of the lower gamma function (Abramowitz and Stegun 1965). Evaluating this derivative at r = R(t) gives Substituting into the Stefan condition, the front velocity dR∕dt from Eq. (28) and the derivative from Eq. (68) yields a nonlinear algebraic equation for the constant m