Dissolution Phase Diagram in Radial Geometry

Fractures play an important role as flow paths in porous media. When an undersaturated reactive fluid is flowing in a fracture, the dissolution process will alter the flow paths locally. Dissolution patterns grow slowly, but they may lead to a dramatic reorganization of the flow. Here, we study dissolution in a radial geometry, which is relevant for a number of practical applications, e.g., the acidization of oil reservoirs. Different dissolution patterns are presented in a phase diagram with the Péclet and Damköhler numbers as parameters. We also analyze quantitatively the density of wormholes and the relation between the aperture roughness and the ramified patterns.


INTRODUCTION
When a fluid flows in a fractured rock, the permeability of the fractures is usually much larger than the permeability of the porous matrix. If the fluid is reactive, the dissolution process will lead to structural alteration of the fracture surfaces. A variety of dissolution patterns have been observed to form on the fracture surfaces ranging from compact to ramified and fractal [1][2][3]. The detailed geometry of these patterns is of central importance for the flow and dispersion properties of the fractures [4]. Some of the structures resemble the forms seen in other phenomena in non-linear science such as dielectric breakdown [5], diffusion limited aggregation DLA [6], and structures observed in two phase flow in porous media [7,8]. Examples of geological processes where reactive dissolution of fracture surfaces are important is weathering and diagenesis of rocks [9,10], and melt extraction from the mantle [11,12]. Reactive dissolution of fractures has also many engineering applications e.g., stability of dams [13] and CO 2 storage [14]. An important application for the oil industry is acidization: the injection of a reactive liquid to increase the permeability of the oil reservoir near a production well [15,16]. Another example is found in the geothermal industry where fluid dissolution is used to increase the surface available for the fluid flow in order to improve heat exchange in the fractures [17][18][19].
In the presence of flow, the dissolution front may become unstable, due to the so called reactive-infiltration instability [10,[20][21][22]. Studies of the reactive-infiltration instability have been performed both theoretically and numerically [3,15,[23][24][25][26][27][28]. At the onset of instability, a linear stability theory can be applied [20][21][22]29]. The theoretical and numerical predictions of the initial instability and the dissolution phase diagram have been compared to experimental results [3,30,31]. However, there are relatively few experimental studies on the onset of the instability [32], especially on the dissolution in the quasi-2D radial geometry.
In the lab experiments, two different setups are usually used: rock core acidization in Hassler cell [33][34][35] and quasi-2D systems in Hele-Shaw cell which are aimed to study the dissolution in quasi-2D porous media and fractures [1,3,32,36]. In the 3D systems, by using the confining pressure at the perimeter of the sample, the core-flow is forced through the bulk of the rock matrix. However, in the quasi-2d systems it is much harder to avoid the focusing of the flow at the boundary. The injected fluid has a tendency to flow along the difficult-to-detect aperture between the medium and the confining cell plate. This problem can be turned into advantage, if we promote it in a controlled manner instead of avoiding the wall flow. Such a controlled-aperture system can then be considered as an analog of a fracture, and the study is then directed at the investigation of how the fracture aperture evolves in time as a result of dissolution. Another experimental advantage of the quasi-2D systems is that it makes visualization easy [4,32,36,37]. Further, a large number of numerical studies have been performed on these systems [28,[38][39][40][41], which provide the insights into the pattern-formation mechanisms. We have chosen a radial Hele-Shaw cell as our experimental system due to the relevance of this geometry to the practical applications such as reservoir stimulation. Flow patterns and dissolution finger growth in quasi-2d radial geometry have been studied in [1,4,25,37]. In 3D, the dissolution patterns in radial geometry have been analyzed e.g., in [34,42,43].
In the present study we have experimentally analyzed a phase diagram for reactive dissolution inside a single fracture in a radial geometry at different flow rates and aperture widths. We have considered two cases: dissolution by central fluid injection and dissolution at the external rim by withdrawal. In section 3, two dimensionless numbers are introduced for the governing equations of the systems. In section 4, we present our experimental results with discussions of the dissolution phase diagram and the statistical properties of dissolution wormholes. The conclusions are presented in section 5.

DESCRIPTION OF EXPERIMENTS
As an experimental setup we use a Hele Shaw cell which consists of two circular glass plates separated by aluminum spacers (spacing thickness b = 1mm). The bottom glass plate (diameter d 1 = 36.0 cm) is larger than the top plate (diameter d 2 = 25.0 cm) and has an external rim to keep the water inside. An outlet hole at the center of the lower glass plate is connected with a tubing to a syringe pump. With this setup we have the possibility either to inject fresh water from the center or withdraw saturated water with fresh water invading from the external boundary. The experimental setup is illustrated in Figure 1.
Before the experiment, the system was first filled with a gypsum saturated water solution and then we injected a plaster paste from the central hole. The plaster paste was prepared by mixing two parts by weight of water and three parts of plaster powder. After injection, the solidifying paste formed a circular plate of radius R 0 = 8.0cm and was left to rest for about 1 h. During this hydration process the plaster shrinked in a process called bleeding (a segregation process) [44] and saturated water is expelled to form a layer between the upper gypsum surface and the upper glass plate. The distance between the plaster and glass plate was measured to be h s = 50µm. To increase the spacing h between the gypsum sample and the upper glass plate further, a number of plastic film spacers (each of thickness h 1 = 100µm) have been inserted between the upper glass plate and the aluminum spacers. In that way we can create an initial aperture of h 0 = 50, 150 and 250µm.
The permeability of the aperture (κ = h 2 0 /12 ≈ 10 −10 m 2 ) is at least 4 orders of magnitude larger than the permeability of the porous matrix of the plaster (κ p = 6.0 · 10 −14 m 2 ) and therefore almost all the water flows through the aperture. During the experiments the temperature was kept at T = 22 • C and we used distilled water with a pH = 7.17 and a viscosity of µ = 1.0 · 10 −3 Pa · s. The molecular diffusion coefficient of plaster (gypsum) in water is D = 1.0 · 10 −9 m 2 /s [45]. The chemical reaction kinetic constant of gypsum in water is k = 4.6 · 10 −6 m/s [46]. We changed the flow rate within the range Q = 0.6ml/h to Q = 60ml/h.The flow rate is controlled by the syringe pump. In the case of injection experiments, the fresh water is injected by the syringe pump from the center of the bottom glass plate. The water saturated with calcium ions is then collected by a water tank connected to the edge of the Hele-Shaw cell. In the withdrawing experiments, the saturated water is withdrawn by a syringe pump from the center of the system while fresh water is supplied at the outer perimeter. When the syringe is completely emptied in the injection experiment or filled for the withdrawing case, we use the 3-way valve connected to a water reservoir to recharge the syringe and ensure continuity of the flow. The water level at outer boundary of the cell is kept fixed by a spill out system above an open rim for the central injection case, or by manual supply of distilled water at the edge of system for the central extraction case, see Figure 1. In total we performed 7 experiments with injection of water and 6 experiments with withdrawing of water, as detailed in Table 1 (7 injection experiments marked with I and  6 withdrawing experiments marked with W). Each experiment was repeated at least twice. Two dimensionless numbers and a penetration length quoted in the table will be introduced in the next section.
The experiment was illuminated from below by a light box with a homogeneous illumination. Pictures were taken from above with a Nikon D7100 camera every 5 min. The whole process was monitored from the beginning of the water injection to the breakthrough of the dissolution structures toward the opposite end of the system.

DIMENSIONLESS NUMBERS FOR DISSOLUTION IN FRACTURES
The aperture of the fracture is several orders of magnitude smaller than its lateral dimensions, thus the mathematical description of the flow and transport in the fracture can be obtained by depth-averaged equations [22,32,38]. Fluid flow is then described by the Reynolds equation for the local volume flux (per unit length across the fracture), q(x, where h is the aperture of the fracture, p is the pressure and µ is the fluid viscosity. Next, the transport of calcium ions  in the aqueous phase is described in terms of flow-averaged concentration field,c(x, y) = 1 where ∇ is a two-dimensional (x, y) gradient operator, D is the diffusion coefficient and R is the sink term related to the dissolution of the gypsum layer. The kinetics of gypsum dissolution expressed by R is, to a good approximation, linearly related to the difference between the saturated concentration c sat and the concentration of the calcium ions at the mineral surface (interface of solid/fluid), c w [45]: with the intrinsic kinetic constant k. To express the reactive current, R, in terms of average concentration, let us note that R needs to be balanced by the diffusive flux at the surface, where n is the vector normal to the surface. The diffusive current can be expressed in terms of the geometry-dependent Sherwood number, Sh [47].
Sherwood number itself depends on kh/D, but the variation is relatively small [48], bounded by two asymptotic limits: high reaction rates (transport limit) and low reaction rates (reaction limit). For our geometry (slot space with one dissolving wall) these limits correspond to Sh = 4.861 and Sh = 5.385 respectively [49]. We approximate the value of Sh as Sh = 5 in the theoretical calculations. Expressing c w is terms ofc using Equation (3) and Equation (5) leads finally to with the effective reactive constant: which accounts for the diffusive slowdown of reaction as the aperture increases. Two different parameters control the behavior of the system. The Péclet number Pe = Qh 0 /DR 2 0 characterizes the relative rate between convective and diffusive transport calculated from characteristic timescales [37]. Next, the effective Damköhler number Da eff = k eff R 2 0 /Q relates the reaction rate to the rate of convective transport [28]. Importantly, in the experiment we can control both the total flow rate Q and the initial aperture h 0 and thus vary Pe and Da values independently [32]. Note that we have chosen the system size R 0 , as the length scale in the definition of the Damköhler number, since what we will use is mainly for the assessment whether the reactant penetration length is comparable with the system size.
At the initial moment, where the flow is isotropic |q(r)| = Q 2π r and the aperture field -homogeneous h(r, φ) = h 0 , the transport equation (2) can be solved analytically. For large Péclet numbers the undersaturation is of a Gaussian form [50]: where r f is the position of the reaction front, defining the reactant penetration length to be The penetration lengths corresponding to the experiments are quoted in Table 1.

EXPERIMENTAL RESULTS AND DISCUSSIONS Experimental Dissolution Patterns
With freshwater injected into the Hele-Shaw cell, the flow dissolves the plaster sample up to the equilibrium state where the gypsum saturates the water. The dissolution patterns evolution with time are obtained by a sequence of images, see Figures 2, 3. In Figure 2 the aperture is fixed to h 0 = 50µm and the injection rate increases from 0.6 to 60ml/h corresponding to the experiments I1, I2, I3, and I4 in Table 1.
In Figure 3, the aperture thickness is fixed with h 0 = 150µm and the injection rate increases from 1.2 to 60ml/h corresponding to the experiments W2, W3, W4, and W5 in Table 1. In the analogy to the previous studies [3,15] we can classify the pattern in Figure 2A as compact, the pattern in Figure 2B as a dominant wormhole, the pattern in Figure 2C as ramified wormholes and the pattern in Figure 2D as uniform.
For very low injection rates (Q < 2ml/h), the chemical reaction reaches the equilibrium right at the interface. In such case, the front between the undissolved and dissolved phase is very sharp. In the low-velocity regime, the dissolution front is stable to perturbations [50]. As a result, the emerging dissolution pattern is compact, see Figure 2A. The breakthrough time (corresponding to the dissolution front reaching the end of the system) is very long in such case, so it is not an efficient way of augmenting the permeability of the rock.
For medium injection rates (2ml/h < Q < 30ml/h), the fresh water can flow further before it becomes saturated. The dissolution front is then more diffuse with the partially dissolved regions observed at the boundary between undissolved and dissolved regions. The dissolution front is now unstable and soon the first perturbations appear, which break the initial homogeneity. These channels are initially very short, but the longer channels with higher pressure gradient will concentrate more fresh flow which in turn enhances the dissolution. With this positive feedback loop, the longest channels soon transform into pronounced fingers. These dissolution fingers (sometimes called wormholes) compete with each other, surviving ones screening off the shorter ones until the longest one finally breaks through [1-3, 15, 33]. At higher injection rates, the wormhole pattern becomes increasingly ramified. There is now more flow to divide between alternate pathways and the competition between them weakens, allowing many branches to co-exist.
For high injection rates (Q > 30ml/h), the fresh water flows very fast in the fracture, the reactant penetration length l p becomes comparable to the system size and the concentration becomes almost uniform along the flow direction. Based on Equation (9), this takes place when Q ≈ k eff R 2 0 /π, which (for R 0 = 8cm) corresponds to Q ≈ 30ml/h. The fresh water can dissolve the entire domain of the plaster sample with almost the same reaction rate which leads to a uniform dissolution pattern. In this case the increase of permeability comes at a price of a large amount of injected water, thus it is not an optimal injection rate to increase the permeability of fracture. Interestingly, we observe that the initially uniform dissolution patterns eventually become unstable and develop fingers, see Figure 2D. This is a manifestation of the inhomogeneous nature of the flow field in a radial system. The magnitude of the flow decreases with the radius as 1/r and thus with time-as the dissolution front moves away from the center-the patterns transition toward those observed at lower flow rates.    For the withdrawing experiments shown in Figure 3, it is hard to obtain a sharp and stable dissolution front. This is because of much lower flow velocities at the outer rim than at the inner rim at the same volumetric flow rate Q. However, the general features of the flow dependence of the wormhole shapes are similar to those observed in injection experiments: at lower flow rates the wormholes are thicker and less dense, whereas at the higher flow rates their density is higher, with many independent flow channels sharing the flow. The flow rate affects also the instability development time -the growth rate of instabilities is the slowest at low flow rates (Figure 3A) with first perturbations appearing only after about 10 days. As the flow rate increases, the instability development time is reduced to a few hours (for Q = 7.2 ml/h) or tens of minutes (for Q = 10.8 ml/h).

Experimental Phase Diagram
The experimental results are summarized in a phase diagram, Figure 4, where the dissolution patterns are presented as function of the Péclet and the effective Damköhler number defined in section 3. In these experiments, the points corresponding to the patterns shown in Figure 4 are marked by the star dots with the experiment number in Table 1. The dissolution patterns are displayed as binarized images with a threshold applied on the gray-scaled pictures analogously to Figures 2, 3.
From Figure 4, the phase transition boundaries among compact, wormhole and uniform dissolution pattern are drawn based on a qualitative observation of dissolution patterns. For the outward flow, we have observed that the boundary between compact and wormhole phases locates around Pe = 5 · 10 −3 and the boundary between wormhole and uniform phases locates around Da = 4. For the inward flow, we have not observed the boundary between compact and wormhole phases and the boundary between wormhole and uniform phases very clearly. Thus, we classify all these dissolution patterns of the withdrawing experiments as wormhole patterns. The fact that the boundaries between the compact (stable) dissolution and wormholing is different for injection and withdrawal reflects the fact that the local flow rates (q in Equation 1) are much lower at the outer rim of the system than at its inner rim (for the same total volumetric flux). Larger flow velocities more easily trigger the front instability [22]. The screening effect between fingers is another important difference between the injection and withdrawing experiments. For the withdrawing case, the fingers grow as the flow converges to the central outlet and the fingers get closer to each other. This increases the competition between the fingers and the pattern becomes strongly hierarchical, with relatively large number of shorter fingers and a few longer ones, focusing most of the flow. Conversely, the competition is much weaker for the outward flow case with injection at the center. There, the spacing between the fingers increases which makes the screening less intense and promotes splitting of the fingers, with new branches appearing in between the old ones, as clearly visible in Figure 4. Similar phenomena were reported for viscous fingering, where also the radial outward flow is associated with splitting [51]. For the inward flow, the splitting is less intense and associated just with the longest fingers where the flow is the highest [52].
The theoretical determination of the position of this boundary line will require a linear stability analysis of the dissolution problem in radial geometry. So far, such an analysis was only carried out in the limit of infinitely thin reaction front [50], which corresponds to significantly lower velocities than the ones used in the present study. We hope that present results will inspire further theoretical work in this direction.

Density of Wormholes
Dissolution patterns have been quantitatively described by fractal dimension in different works [2,53]. However, particularly in the wormhole regime, our dissolution fingers have cone-like shapes tapering toward the tips and the dissolution pattern is not selfsimilar [4]. Here we will use the concept of wormhole density [54,55] to describe our dissolution patterns. The wormhole density ρ N is defined as the number of wormholes per centimeter along a perimeter of circle with a radius r. A wormhole in our counting procedure is defined as a continuous dissolved part with a minimum width 2 mm. A minimum width is set to avoid miscounting pixel noise as wormholes at the dissolution front. The density analyses are performed on the final states of each dissolution process. In Figure 5, a schematic diagram displays an example of how we calculate the number of wormholes intersecting a circle N wh at a specific radius r.
From Figure 5, we can count the number of peaks of the distribution with a minimum peak width 2 mm as the number of wormholes. We plot the number of wormholes N wh as a function of radius r for the different Péclet numbers Pe, see Figure 6. For the compact dissolution pattern at low Péclet numbers (Pe = 1 · 10 −3 ), as expected, the number of wormholes remains zero except for some variations at the front since the compact pattern is not a perfect circle. For the other 3 experiments carried out at different Pe, the number of wormholes first increases to a maximum, then gradually decreases. The wormhole density shows a similar tendency: as the radius r increases, wormholes grow and split, side fingers start to grow from the roots, which results in an increased number of branches. For the case of dominant wormholes (Pe = 7 · 10 −3 ), some of the branches eventually screen the other, which leads to the appearance of the maximum in the N wh (r) dependence. On the other hand, for the ramified pattern (Pe = 0.01) the side branches continue to grow until nearly the end of the system. Finally, for the uniform dissolution pattern (Pe = 0.1), the situation is somewhat similar to the compact case, except at the front, where a number of relatively wide fingers appear as described in section 4.1.

Roughness and Ramified Fingering
Besides the wormhole density, other significant characteristic features of the wormholes are their widths and the extent of ramification. In some classical instabilities, as Saffman-Taylor instability, the finger width as function of distance to tip is one of the most finely studied quantities [56]. The ramification and thickness of the viscous fingering patterns have been characterized in water and miscible non-newtonian fluid system [57] and in drainage in a porous medium [58]. In the fracture dissolution system, the simulation results show that the rougher the aperture, the more ramified dissolution fingers will be    7 | Ramified or smooth fingers with different initial roughness γ 0 and injection rate Q. From top to bottom, the roughness increases from γ 0 = 0.07 to γ 0 = 0.20 for the initial average aperture h 0 = 50.0µm (top) and 150.0µm (bottom), respectively. From left to right, the flow rate is Q = 3ml/h, 6ml/h, 7ml/h, respectively. Notice that the rightmost two dissolution patterns are obtained in the withdrawing experiments. We observe that as the roughness is increased, the patterns become thinner and more ramified.
Frontiers in Physics | www.frontiersin.org 9 September 2020 | Volume 8 | Article 369 observed [38,59,60]. This has also been qualitatively verified by the observation in both inwards and outwards flow in the radial geometry experiments. The roughness is defined as: where h = h(x, y, t) is the aperture, with a initial mean aperture h 0 = < h(x, y, t = 0) >. The initial aperture field is measured before dissolution using digital optical microscope Keyence VHX5000. The results give the standard deviation of the initial aperture σ 0 = 10.0µm where σ = √ < h 2 > − < h > 2 . Taking into account that h 0 varies from 50.0 to 150.0 µm, we get the initial roughness in the range 0.07 to 0.2. For the higher roughness, experimental images show more ramified and thinner fingers than those observed at lower roughness, see Figure 7.

CONCLUSIONS
The problem of flow-induced fracture dissolution has received a lot of interest in recent years due to its importance in geological phenomena as well as industrial applications. However, to our knowledge, the experimental evidence of radial dissolution in quasi-2d geometry remains scarce. In this paper, we used a plaster disk sample in a Hele-Shaw cell to study the flow-induced dissolution patterns formation in a circular geometry.
We classified the dissolution patterns and analyzed their dependence on the fundamental dimensionless groups characterizing convection, diffusion and reaction, i.e., the Péclet and Damköhler number. By changing the aperture of our Hele-Shaw cell as well as the flow rate, we can vary Péclet and Damköhler number independently and thus obtain a phase diagram of dissolution patterns in Pe-Da variables. The general features of this phase diagram are not unlike those of an analogous diagram for rectangular geometry [3], but there are also important differences. These are connected with the inhomogeneous nature of the flow in the system, with a flow rate decreasing as we move away from the inlet. Because of this, the dissolution features tend to transition with time to the forms characteristic of lower flow rates: uniform patterns become unstable and develop fingers, ramified wormholes transition toward the dominant ones. We hope that these experimental results will constitute a benchmark for the further theoretical and numerical studies of the dissolution process in a radial geometry.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
LX and KM designed the experiment. LX performed the experiments, analyzed the data, and authored the paper. PS, RT, EF, and KM assisted with the interpretation and data analysis and editing of the manuscript. All authors contributed to the article and approved the submitted version.