Influence of Layering on the Formation and Growth of Solution Pipes

In karst systems, hydraulic conduits called solution pipes (or wormholes) are formed as a result of the dissolution of limestone rocks by the water surcharged with CO₂. The solution pipes are the end result of a positive feedback between spatial variations in porosity in the rock matrix and the local dissolution rate. Here, we investigate numerically the effect of rock stratification on the solution pipe growth, using a simple model system with a number of horizontal layers, which are less porous than the rest of the matrix. Stratification is shown to affect the resulting piping patterns in a variety of ways. First of all, it enhances the competition between the pipes, impeding the growth of the shorter ones and enhancing the flow in the longer ones, which therefore grow longer. This is reflected in the change of the pipe length distribution, which becomes steeper as the porosity contrast between the layers is increased. Additionally, stratification affects the shapes of individual solution pipes, with characteristic widening of the profiles in between the layers and narrowing within the layers. These results are in qualitative agreement with the piping morphologies observed in nature.

1. Introduction

The processes leading to the dissolution of a porous or fractured rock by a reactive fluid depend on a subtle interplay between the chemical reactions at mineral surfaces and fluid motion in the pores [1–3]. Strong, non-linear couplings between flow and dissolution may lead to the appearance of instabilities resulting in the formation of a variety of different geological patterns. Particularly relevant for the present study is the so-called reactive-infiltration instability [4, 5]: a relatively small change in pore geometry due to the dissolution may result in a variation of permeability, focusing the flow and leading to further changes in erosion. As a result, small inhomogeneities have a tendency to grow and get transformed into highly permeable, finger-like flow channels, which go by several different names depending on the field. In the petroleum industry, they are called “wormholes” [6] due to the resemblance with the tunnels dug by worms. Geologists call them “solution pipes,” “karst funnels,” or “geological organs” [1, 7, 8], whereas pedologists call them “soil tongues” [9].

Once formed, the channels continue to grow, focusing an ever-increasing portion of the flow, but also competing with each other: since the larger flow in the channel leads to increased dissolution—this again generates the positive feedback loop resulting in fast growth of the longer pipes and starvation of the shorter ones.

The flow-focusing instability has been known for a long time in the petroleum industry [10], where it was used to increase the effectiveness of reservoir stimulation techniques. In that context, it has also been studied experimentally [6, 11] and numerically [12–15]. In the physics community, the interest toward these systems has been sparked by classical papers of Daccord and Lenormand, which analyzed the patterns formed in a dissolving plaster and related them to other pattern-forming systems, such as diffusion-limited aggregation [16, 17]. In earth sciences, the understanding of the profound role that flow-focusing plays in karst development came with the pioneering works of Ewers [18], who has analyzed a number of analog models of karst formation made of sand or plaster of Paris. This has provided him with a number of important insights on the process of cave formation, and the role played in these processes by the competition between evolving flowpaths. Further insight into these effects has been provided by the numerical modeling of early stages of karst development [19–25].

Recent years have seen a growing number of experimental studies of dissolution-induced morphology evolution [26–29], many of them in the context of CO₂ sequestration [30–38]. Since many of the potential reservoirs for CO₂ storage are in carbonate strata, it is crucial to understand how the flow of CO₂-acidified brine impacts the long-term changes in porosity and permeability of the reservoirs and how it affects the caprock properties. Most of the systems analyzed in these studies are relatively small (1 − 10 cm scale), however even at that scale a strongly non-uniform character of the dissolution manifests itself in a profound way, with the appearance of wormholing patterns. Due to the self-similar character of the flow-focusing process [24, 39], one expects to find the dissolution channels at all scales, up to the system size. And this is indeed the case, with the karst solution pipes appearing in many places throughout the world [1, 7, 8, 40, 41], most of them on 1–10 m scale, although much longer structures, up to kilometer length, have also been reported [42].

In this paper, we focus on solution pipes in karst systems, analogous to the ones depicted in Figure 1. In particular, we investigate the influence of rock stratification on the solution pipe growth. As it turns out, the presence of low-porosity layers impacts the shapes of the solution pipes, which tend to narrow as they cut through the layers and then widen up to form bulbous caverns. These features are visible in the pipes in Figure 1 and their detailed investigation is the goal of this study.

FIGURE 1

Figure 1. Examples of the solution pipe morphologies in stratified systems. (A) Solution pipes in gravel with limestone clasts: gravel-pit, west of Rauceby Station, UK. Photo by J. Rhodes courtesy British Geological Survey (catalog no. P206264). (B) Solution pipe in Ramallah, Israel, cutting through stratified dolomites with marl beds. The photograph is courtesy of Sebastian Schmidt of University of Göttingen. (C) Solution pipe cutting through layered limestone rock at Cape Bridgewater, Victoria, Australia (D) three pipes with cemented rims in limestone formations at Burns Beach, Perth, Australia. Pipe no. 3 cuts through the more resistant bed in the lower portion of the photo.

So far, the shapes of wormholes have been analyzed mainly in the context of acidization [11–13, 43, 44], where it has been established that they strongly depend on the flow rate and reaction rate. At low injection rates the growth of the dissolution channels is significantly influenced by diffusion and the resulting wormholes are conical-shaped. At higher flow rates the channels become more elongated—this regime is called in the acidization literature a “dominant wormhole” regime. In the geomorphological context, a systematic analysis of the lengths of solution pipes has been carried out by De Waele et al. [8], who have measured a collection of more than 200 solution pipes in the coastal areas of the Mediterranean and analyzed their length distribution, suggesting that it follows an exponential law. This data has been subsequently re-analyzed in Ref. [45], which led to the conclusion that—if one discards the partially buried pipes the full length of which remain unknown—the length distribution is found to be consistent with a power-law. Such scale-invariant distributions are a characteristic feature of a large number of hierarchical growth processes [46], observed experimentally e.g., in some of the viscous fingering systems [47, 48], dendritic side-branches growth in crystallization [49] or crack propagation in brittle solids [50].

To the best of our knowledge, the shapes of the solution pipes in the layered porous medium have not been analyzed in the literature. However, some of the insight into these questions can be gained from studying of the related problem of the growth of karst cutters [51, 52], i.e., solutionaly enlarged joints. The dissolution of cutters is somewhat different from that of solution pipes. In the former case the flow is focused from the very beginning within the joints, whereas in the latter the dominant flow pathways are created spontaneously due to the reactive-infiltration instability and the subsequent competition between the emerging wormholes. Nevertheless, some of the morphological features of the cutters and pipes seem to be similar. In particular, Howard [51] observes undulations on the sides of the cutters, which correlate with the strata and notes that cavernous widenings have the tendency to form predominantly beneath the layers. Analogous morphological characteristics are observed in the case of solution pipes (cf. Figure 1 and the numerical results of Section 3). Understanding of the mechanism underlying the formation of these features is one of the goals of this study. Additionally, we will also go beyond the scale of a single pipe and analyze the impact of the layering on the competition between different pipes for the available flow.

The paper is organized as follows. In Section 2 the Darcy-scale mathematical model of a dissolving porous rock is briefly recalled. In Section 3 we introduce the model of a stratified medium and analyze numerically its impact on the dissolution patterns. Finally, the conclusions are drawn in Section 4.

2. Governing Equations

Let us briefly recall the equations for the dissolution of a porous matrix [53]. The rate of groundwater flow through the porous medium is taken to be proportional to the pressure gradient (Darcy's law)

$\begin{array}{cc}u=-\frac{K(\varphi )}{\mu }\nabla p& (1)\end{array}$

where u is the Darcy velocity, K-the permeability, μ-the viscosity of the fluid, ϕ-the porosity, and p-the fluid pressure. Here we will adopt the Carman-Kozeny relation for the permeability dependence on the porosity, K(ϕ) ~ ϕ³. We also assume that the Darcy velocity field is incompressible,

$\begin{array}{cc}\nabla ·u=0,& (2)\end{array}$

neglecting contributions to the fluid volume from reactants or dissolved products. Under typical geophysical conditions, dissolution is slow in comparison to flow and transport processes; we can therefore assume a steady state in both the flow and transport equations.

The concentration of the solute, c is calculated through the convection-diffusion-reaction equation

$\begin{array}{cc}\nabla ·(uc)-\nabla ·(D\varphi ·\nabla c)=R(c),& (3)\end{array}$

where D is the diffusion coefficient and R(c) is the reactive flux into the matrix. Finally, the erosion equation, which gives the increase of the porosity as a result of the dissolution, is

$\begin{array}{cc}{c}_{sol}{\partial }_t\varphi =R(c)& (4)\end{array}$

where c_sol is the mineral concentration of the solid species and R(c) is the reaction rate. Here, we assume a linear kinetic equation for the reaction rate

$\begin{array}{cc}R(c)=ks(c_{sat}-c)\theta ({\varphi }_{max}-\varphi )& (5)\end{array}$

where k is the kinetic rate and s is the specific reactive surface area per unit volume that we assume to be constant, s = s₀. Next, c_sat is the saturation concentration at which adding further solute results in precipitation. In the context of karst formation, c corresponds to the concentration of calcium ions, which controls the dissolution at moderate to large pH [54]. The Heaviside step function, θ(ϕ_max − ϕ) guarantees that the reaction term vanishes in the region where all of the soluble material has dissolved (which corresponds to the porosity ϕ = ϕ_max). In this region Equation (3) reduces to a convection-diffusion equation.

The above equations are supplemented by boundary conditions on the velocity and concentration field. Along the x direction (parallel to the flow):

$\begin{array}{cc}p(x=0,y)=p_{0}c(x=0,y)=0,& (6)\end{array}$

and

$\begin{array}{cc}\begin{array}{ll}p(x=L_x,y)=0\hfill & c(x=L_x,y)=c_{sat},\hfill \end{array}& (7)\end{array}$

which accounts for the applied pressure drop p₀ and the inflow of unsaturated solution imposed at the inlet. In the case of solution pipes, the flow is mostly due to gravity. In the simulations, we keep the total volumetric flow rate through the system constant, and p₀ is changed dynamically as the permeability of the system evolves. Along the y direction periodic boundary conditions are imposed. The above equations are solved numerically using a finite difference method, described in more detail in Ref. [23]. The initial porosity field is sampled from a lognormal distribution characterized by the standard deviation σ and the correlation length λ. In most of the simulations presented here, we choose a relatively small standard deviation of the initial porosity field, σ = 10⁻³.

The inverse of the specific surface area, $s_0^{-1}$ is a characteristic length in our system, which will be chosen as a unit of length in all the reported data, with the dimensionless coordinates $\widehat{x}$ and $\widehat{y}$ defined as

$\begin{array}{cc}\widehat{x}=x{s}_0\widehat{y}=y{s}_0& (8)\end{array}$

On the other hand, the time is scaled by a characteristic dissolution time:

$\begin{array}{cc}\widehat{t}=t/\tau & (9)\end{array}$

$\begin{array}{cc}\tau =\frac{c_{sol}}{k{s}_0{c}_{sat}}.& (10)\end{array}$

Finally, the dimensionless concentration is defined as

$\begin{array}{cc}\widehat{c}=\frac{c_{sat}-c}{c_{sat}}& (11)\end{array}$

which is constrained to be in the interval (0,1). We also scale porosity by its maximum value, ϕ_max,

$\begin{array}{cc}\widehat{\varphi }=\frac{\varphi }{{\varphi }_{max}}.& (12)\end{array}$

The relative importance of convective, diffusive and reaction processes in the system are measured by the Péclet number and Damkölher number. The former is defined as [23]

$\begin{array}{cc}\text{Pe}=\frac{u_0}{D{s}_0{\varphi }_0}& (13)\end{array}$

and measures the ratio of the convective fluxes to the diffusive fluxes. Here u₀ is the average Darcy flux in the initial system whereas ϕ₀ is an initial porosity of the system. On the other hand, the Damkölher number

$\begin{array}{cc}\text{Da}=\frac{k}{u_0}& (14)\end{array}$

relates the surface reaction rate to the fluid flux.

As shown in a number of studies [11–13, 43, 55] the wormhole patterns depend strongly on the values of these two numbers. In particular, it has been shown that the Péclet number needs to be large enough for the instability to occur. This is why piping is often observed in unconsolidated rock and gravel, where the permeabilities (and thus flows) are large [56].

3. Dissolution of a Stratified Medium

To study the influence of stratification on pipe shapes and dynamics, we have adopted a simple conceptual model of a stratified medium (see Figure 2). The porous medium is assumed to consist of alternate layers of higher and lower porosity, ϕ_a (of depth L_a) and ϕ_b (of depth L_b), respectively. The two parameters characterizing such a setup are the porosity ratio C = ϕ_a∕ϕ_b and the length ratio κ = L_a∕L_b. In the following, we limit ourselves to two dimensions. However, the same methodology can be applied to the analysis of the 3D system, which will be the subject of a further study.

FIGURE 2

Figure 2. Schematic of a model of stratified medium: alternating layers of porosities ϕ_a and ϕ_b.

Figure 3 presents the results of the dissolution simulations for such a system for Pe = 100, Da = 0.5 with a porosity contrast of C = 2 and length ratio κ = 5. The characteristic undulations in the shape of the dissolution pipes are clearly visible here—the pipes become narrower in the less porous layers and widen up in between the layers. Another effect, clearly observed in Figure 3 is that the less porous layer acts as a dissolution barrier, which stops the growth of shorter pipes, as there is not enough reactant transported through their bodies to be able to break through the layer. As a result, the tips of the shorter pipes often coincide with the layer boundary. This effect can also be observed in Figure 1D, where only the longest pipe (pipe no. 3) has succeeded in breaking through the packed layer.

FIGURE 3

Figure 3. Solution patterns at Pe = 100, Da = 0.5, the porosity contrast C = 2 and length ratio κ = 5 for a domain of the size $L_x=L_y=1000s_0^{-1}$. The colors represent rescaled porosity, $\widehat{\varphi }$, as indicated in the colorbar on the right.

3.1. Flow and Concentration Fields Near the Tips of the Pipes

Some insights into the mechanisms underlying this behavior can be obtained by the analysis of the flow and concentration fields in the vicinity of the tip of the pipe. These are shown in Figures 4, 5 for three different positions of the tip with respect to the layers: the pipe terminating in a porous region between the layers (left panels), the pipe terminating within the packed layer (center panels), and the pipe emerging from the packed strata (right panels).

FIGURE 4

Figure 4. The Darcy flow field around the tip of the pipe before hitting one of the layers (left), with the tip located within the layer (center), and after emerging from the layer (right). The colors represent the magnitude of the flow scaled by v₀. The lengths of the vectors in the overlaying vector plot has been scaled to v(x, y)^1∕5 to aid visualization.

FIGURE 5

Figure 5. The concentration field, ĉ(x, y), around the tip of the pipe at the same times as in Figure 4.

There are several things to be noted here. First of all, since the permeability of the pipe is much larger than that of the surrounding medium, the tip of the pipe is a local pressure maximum, distributing the flow to the porous matrix, both in the axial and in lateral direction (cf. Figure 4). However, due to the small permeability of the matrix in the packed layers, the lateral flows are diminished when the pipe enters a layer (center panel of Figure 4) with respect to the situation when the tip find itself in a more porous layer (left panel of this Figure 4). Smaller lateral flows give rise to a narrower region of elevated concentration (cf. Figure 5), which decreases the diameter of the pipe as it enters the packed strata.

Next, as the pipe emerges from the layer, the flow is distributed almost uniformly over a semi-circle in the downstream direction with the tip as its center, analogously to the flow created by the point source at the wall (right panel of Figure 4). The unsaturated solution is carried deeper into the matrix (right panel of Figure 5) which gives rise to a significant widening of the pipes as they emerge from the layer. However, as the pipe tip moves away from the layer, this effects weakens and the pipe comes back to its original width.

As observed in Figure 5, the undersaturated solution is concentrated along the length of the pipe with a halo around the tip. However, since the inside of the pipe is fully dissolved (with ϕ = ϕ_max), the reaction term vanishes there and the concentration is simply transported along the length of the pipe. The situation changes in the vicinity of the tip, where strong lateral currents push the undersaturated solution into the undissolved matrix. This is amply illustrated in Figure 6, where we plot the time derivative of the porosity field, $d\widehat{\varphi }∕d\widehat{t}$, which marks the regions of the active reaction. Comparing Figure 6 with the corresponding concentration field (right panel of Figure 5), we see that the dissolution is mainly taking place around the tip and in a thin layer along the perimeter of the pipe.

FIGURE 6

Figure 6. The porosity change, $d\widehat{\varphi }∕d\widehat{t}$, corresponding to the concentration field in the right panel of Figure 5.

3.2. Comparison with Natural Systems

The characteristic undulations of the solution pipes observed in the simulations are consistent with the features observed in natural systems. Figure 1A presents a number of spontaneously formed solution pipes in gravel with limestone clasts. Undulations in different pipes are spatially correlated, which suggests the presence of at least three layers in this system. The dark rims along the perimeter of the pipes are most likely the result of clay and iron oxide accumulation due to the illuviation processes [57]. The clay particles are transferred by water from the upper parts of the soil and then flocculated at the clay-limestone boundary, where the pH changes from mildly acidic to alkaline. Such rims are commonly seen in solution pipes in which the infill is still present [7, 58].

In Figure 1B we see a solitary solution pipe, which might have been triggered by the presence of the tree. In fact, stemflow is usually significantly more acidic than incident precipitation [59] and thus it is often considered as one of the factors speeding up solution pipe formation [40, 41]. The pipe in Figure 1B cuts through several dolomite layers, with characteristic widenings, very similar to those observed in the numerical simulations. Finally, Figures 1C,D show the pipes formed in beach calcarenites in southern and western Australia. Again, the narrowings and widenings of the pipes as they cut through the layers are clearly visible here.

3.3. Hierarchical Growth of the Pipes and Their Length Distribution

The value of porosity contrast has a profound impact on the patterns, as it affects the intensity of the competition between the growing pipes. As observed in Figure 7, as C gets larger the number of long pipes in the system diminishes and the number of short ones increases. This means, however, that the longest pipes are getting significantly more flow in these systems (since the flow is then distributed between a smaller number of pipes). As a result, the longest pipe grows faster in the stratified systems with a large porosity contrast.

FIGURE 7

Figure 7. Piping patterns for C = 1 (no layers, left), C = 2 (center), and C = 4 (right) for Pe = 100, Da = 0.1, and κ = 5.

These observation can be quantified by the analysis of the distribution of the lengths of the pipes. As shown in a number of studies [24, 39, 45] a strong screening between the dissolution channels usually leads to the power-law distribution of their lengths,

$\begin{array}{cc}N(L_p)~L_p^{-\alpha },& (15)\end{array}$

where N denotes a number of pipes longer than L_p. In a homogeneous system, the exponent α is close to unity [24, 39, 45]. Figure 8 shows the length distributions corresponding to the patterns presented in Figure 7. In the homogenous case (C = 1) the results are in full agreement with the earlier studies, indicating $L_p^{-1}$ distribution over the entire length range. However, the situation is different for the layered system. The region characterized by the power-law distribution is now much narrower followed by a clearly visible shoulder in the distribution extending up to the location of the first layer (at $100s_0^{-1}$), which corresponds to the surplus of pipes which have been stopped in their growth by the presence of the layer. Additionally, the power-law exponents for the layered case are significantly larger than 1, which confirms that the stratification indeed increases the screening between the pipes.

FIGURE 8

Figure 8. The cumulative distribution N(L) of solution pipe lengths for the patterns presented in Figure 7 with C = 1 (red points), C = 2 (blue), and C = 4 (green). The solid lines have slopes of 1.02 (red), 1.95 (blue), and 2.45 (green).

A quantity directly linked to the length distribution is the mean spacing between adjacent active pipes (i.e., pipes which continue to grow due to the presence of undersaturated solution). The initial spacing between the spontaneously formed pipes is determined by the spectral characteristics of the reactive-infiltration instability [4, 23, 60], which triggers the pipe formation process. Soon, however, the competition between the pipes for the available flow leads to the starvation of the shorter ones and the spacing between active pipes increases. Note that for the power-law length distribution (15) the characteristic separation between pipes of length L_p is equal to $L_y∕N(L_p)~L_p^{\alpha }$, where L_y is the width of the domain. Hence, the increase of α exponent in the stratified case results in a larger spacing between the active pipes.

Before we move on, let us analyze the impact of the amplitude of the initial noise on the dissolution patterns. Figure 9 presents the same systems as in Figure 7 but with a much higher standard deviation, σ = 0.15. As observed, the pipes are now much more irregular and winding, but the overall patterns are very similar to those in Figure 7. This conclusion agrees with a more detailed analysis of Upadhyay et al. [45] which shows that the statistical properties of the dissolution patterns are only weakly sensitive to the parameters governing the randomness in the initial porosity field. In the remainder of the paper, we focus on the low-noise regime, since a more regular geometry of the pipes in this case allows us to better discern the effects of the stratification of the porous matrix.

FIGURE 9

Figure 9. Same as in Figure 7 but for a larger value of the initial noise.

3.4. Two-Pipe Case

To get more insight into the nature of the competition process, we have considered a case in which two pipes grow symmetrically at equal distances from the centerline of the system. To trigger the growth of the solution pipes at these places, small regions of increased porosity (“pipe seeds”) have been placed in an initial porosity field. These seeds have a length of 3% of the system length and are placed at L_y∕4 and 3L_y∕4, respectively. Although the seeds are equal at the beginning, small noise in the initial porosity field triggers the instability which eventually leads to one of them outcompeting the other, as shown in Figure 10.

FIGURE 10

Figure 10. The positions of the pipe tip vs. time in a two-pipe system in non-stratified (left) and stratified case. Red color marks the faster, winning pipe, and blue - the slower one. The 2d contour-plots show the corresponding dissolution patterns.

An intriguing observation inferred from Figure 10 is that in the low porosity layers the speed of the pipe growth increases. This is connected with the narrowing of the pipe in the packed layers and the related increase of the axial fluid velocities. This is a strong effect, which allows the pipe tip to advance faster, even though locally, at each specific point of the packed strata, the total amount of rock to dissolve is larger than at a corresponding point in the less porous layers.

Another observation is that—in the stratified case—just after the growth of the shorter pipe is stopped, the longer pipe starts to grow faster, since it is getting a larger fluid flux. Its speed is higher than that in the case of non-stratified medium (for example, at $\widehat{t}=200$ it reaches $\widehat{x}=540$ whereas the tip of pipe in a non-stratified medium is located at $\widehat{x}=480$ at the same time). With time, however, the growth of the longer pipe in a stratified medium becomes impeded.

To understand the origin of such a slowdown of the dissolution in the stratified case, we have analysed the concentration profiles along the centreline of the winning solution pipe. Figure 11 presents such profiles for non-stratified medium at Pe = 8 and Pe = 100 and for stratified medium (C = 10, κ = 5) at Pe = 100. As observed, the undersaturation gradually decreases along the centreline. Since there is no reaction inside the pipe (the material there is fully dissolved), such a decrease must be connected with the diffusion of the calcium ions from the sides of the wormhole (where the reaction takes place) toward its center. This is further confirmed by the fact that the centreline concentration profile is strongly dependent on the Péclet number, with a characteristic lengthscale which is a function of the flow rate and diffusion within the pipe (cf. the left and center panel of Figure 11). In fact, the influence of Péclet number on the flow within the pipe is twofold. First, the total fluid flux through the cross-section of the pipe is smaller as Pe decreases. Second, the characteristic value of the local Darcy velocity within the pipe is further diminished, since the pipe diameter gets larger with the decreasing Pe.

FIGURE 11

Figure 11. Concentration (red) and porosity (blue) profiles along the line y = 3∕4L_y, corresponding to the centreline of the winning pipe. Left: non-stratified case at Pe = 8, center: non-stratified case at Pe = 100, right: stratified case at Pe = 100.

Interestingly, in the stratified case, such a “diffusive” effect within the solution pipes is significantly stronger, and the profile at Pe = 100 (right panel of Figure 11) has now a decay rate similar to than in the Pe = 8 case in non-stratified medium. Although the quantitative model of this behavior still needs to be developed, the most plausible explanation is connected to the presence of widening in the pipe shape due to the stratification. The undulated shape has a much larger perimeter than its smooth counterpart, which increases the total diffusive flux of the calcium ions from the pipe walls.

3.5. Effects of the Layer Thickness, Porosity Contrast, Péclet, and Damkölher Number on the Dissolution Patterns

Let us now analyze the influence of the thickness of the packed layer on the wormhole shape. Figure 12 compares the shapes of the pipes which emerge in the system with a single packed layer of two different widths. There are several things to be noted here. First, the shape of a bulbous widening carved by the pipe as it emerges from the layer is independent of the layer thickness. Similarly, the width of the pipe at the end of the packed layer, just before the emergence, again seems to be thickness-independent. The main differences between the shapes are visible in the upstream region: whereas in the thick-layer case the pipe is seen to have a necking at the entrance of the packed strata, its counterpart in the thin-layer case bulges out slightly, first widening and then gently narrowing as it approaches the layer.

FIGURE 12

Figure 12. Influence of thickness of the layer on the wormhole shape (Pe = 100, Da = 0.1, C = 10). The solution pipe in red corresponds to the case where the layer extends from the red vertical line to the black vertical line. The solution pipe in blue corresponds to the case where the layer extends from the blue line to the black line. Note that the aspect ratio of the plot has been altered to better visualize details.

Next, let us investigate how the shapes of the pipes are affected by the values of porosity contrast C and the Damkölher number, Da. To this end we have carried out a number of simulations on the 1000 × 400 system with Pe = 100 and Da = 0.1 with different porosity contrasts (1 < C < 4) followed by a series of simulations with fixed C = 10 and different Da. The resulting wormhole shapes are presented in Figures 13, 14, respectively.

FIGURE 13

Figure 13. Piping patterns for different contrasts of porosity (as marked) for Pe = 100, Da = 0.1. The aspect ratio of the plot has been altered to better visualize details.

FIGURE 14

Figure 14. Piping patterns for different Damkölher numbers for a stratification with a porosity contrast of C = 10. The aspect ratio of the plot has been altered to better visualize details.

A porosity contrast is seen to affect the patterns in a significant way. Whereas, in the non-stratified case the pipe has a very regular, cigar-like shape, introduction of even a small porosity contrast (C = 1.25) results in the appearance of pronounced undulations. Their amplitude increases with C-the pipes narrow down progressively in the packed layers and then widen up in the next porous zone, following the rule that the greater the narrowing, the more pronounced the subsequent widening. This can be rationalized by noting that—due to mass conservation—a small diameter of a pipe entails large flow velocities, which then get diverted in the lateral direction as the pipe emerges from a packed layer, resulting in larger cavities. Interestingly, after the bulge, the pipes quickly come back to their original widths, which do not seem to depend on the porosity contrast.

Furthermore, a close inspection of Figure 13 reveals that the pipe begins narrowing already at some distance upstream from the low porosity layer and the higher the contrast is, the earlier the narrowing is taking place. It is particularly pronounced in the space between the first and the second layer—e.g., the C = 10 pipe (marked by an orange curve) begins to narrow down just after it has emerged from the bulge, whereas C = 2 and C = 1.25 pipes continue to grow for a while with a constant width, before eventually narrowing.

Let us now analyze Figure 14, which shows the influence of the Damkölher number on the pipe shapes. Strikingly, in this case the wormholes seem to have the same basic shape and the effect of Da comes down to the appropriate rescaling: the smaller the Da number, the larger the widenings of the pipe. In order to quantify these effects, we have measured the width and length of the bulbous widening (W and L in Figure 14) starting from the point where all of the pipe outlines overlap. Both W and L are plotted as a function of the Damkölher number in Figure 15. Interestingly, both seem to scale like $1∕\sqrt{\text{Da}}$, which can be understood in terms of a simple model in which we assume the tip of the pipe emerging from the packed layer distributes the flow uniformly in the semi-circle (−π∕2 < θ < π∕2). The Darcy flux near the tip is then

$\begin{array}{cc}u=\frac{Q}{\pi r}{e}_r& (16)\end{array}$

where r is the position vector in the polar coordinate system centered at the tip and Q = v₀L_y is the total fluid flux in the system (which we assume to be focused in the pipe). Since the flows within the pipe are high, we can neglect diffusion and obtain from the transport equation (near the tip).

$\begin{array}{cc}\frac{Q}{\pi r}\frac{dc}{dr}=k(c_{sat}-c),& (17)\end{array}$

the solution of which reads simply

$\begin{array}{cc}c(r)=c_{sat}\left(1-e^{-\frac{\pi k{r}^2}{2Q}}\right)& (18)\end{array}$

with a characteristic decay length given by

$\begin{array}{cc}l=\sqrt{\frac{Q}{\pi k}}=\sqrt{\frac{L_y}{\pi \text{Da}}}& (19)\end{array}$

which indeed gives the Da^−1∕2 scaling, as observed in the simulations. Using L_y = 400 (which corresponds to the width of the system in Figure 13), we get the coefficient of 11.3 in this power-law. The difference between this value and the values reported in the fits in Figure 15 (i.e., 6.5 and 41.5 respectively) is connected with the highly simplified character of the model. In reality, the outflow from the tip is expected to be anisotropic, with a larger flux directed along the axial direction [which is in fact reflected in a larger value of the coefficient of the L scaling (41.5) with respect to that of the W scaling (6.5)]. Next, the scaling Equation (19) has been derived for the concentration field. Porosity increase, however, is proportional to a time integral of the concentration [cf. Equation (4)]. Hence, in principle, one would need to know the entire time evolution of the c field to be able to estimate the respective porosity growth. However, as illustrated in Figure 6, the dissolution is active almost exclusively in a small region at the tip of the pipe. Hence, the cavern grows over a relatively short time only, T = L∕V, where V is the growth velocity of the pipe. Our approximation then corresponds to neglecting the flow and concentration changes over that time.

FIGURE 15

Figure 15. The width (W) and length (L) of the widenings, as marked in Figure 13, as a function of the Damkölher number. The points represent the simulation data, whereas the curves are fit to $1∕\sqrt{\text{Da}}$ dependencies: $L(\text{Da})=41.5∕\sqrt{\text{Da}}$ (blue curve) and $W(\text{Da})=6.5∕\sqrt{\text{Da}}$ (red curve).

Moreover, the assumption that the flow from a pipe breaking through the layer corresponds to a point source at the wall is expected to hold only at large values of the porosity contrast, C. There are two reasons for this. The first is that at large C the narrowing of the pipes in the packed layers is more pronounced. Second, the smaller C ratios mean that a larger portion of the flow is transported through the porous matrix in the layer. By performing the simulations at different porosity contrast we have checked that the scaling (19) indeed deteriorates progressively as C is decreased, with both L and W decreasing less pronouncedly with Da.

It is worth noting that in three dimensions, an analysis analogous to Equations (16–19) would yield Da^−1∕3 scaling, thus the dependence on the dissolution rate is weaker in 3D.

Finally, let us discuss the influence of the Péclet number on the patterns. Figure 16 presents the counterpart of Figure 7 for a relatively small value of the Péclet number (Pe = 1). Arguably the most pronounced effect that the decrease of the flow rate has on a dissolution patterns is the change of the overall shape of the pipes. Instead of the linear structures, the pipes are now more conical, with the diameter tapering toward the tip. However, the effects of the stratification remain qualitatively similar to those at higher Pe, with characteristic undulations superimposed on the conical shapes of the pipes. Getting back to the natural examples of Figure 1 we can infer that the pipe in Figure 1C has been formed at a relatively small value of Pe, since the tapering is clearly visible here, whereas more linear forms of Figures 1A,B,D are probably associated with higher Péclet numbers.

FIGURE 16

Figure 16. Piping patterns for C = 1 (no layers, left), C = 2 (center), and C = 4 (right) for Pe = 1, Da = 0.1 and κ = 5.

4. Conclusions

In this paper, we have been analyzing the formation of solution pipes in stratified soluble rocks. We have shown that stratification influences the dissolution patterns in a significant way. By serving as a flow and dissolution barrier, the packed layers impede the growth of the shorter pipes and direct the entire fluid flux into the winning pipes. This affects the distribution of pipe lengths, which becomes steeper as the porosity contrast between the layers is increased.

Interestingly, the shapes of the pipes itself are also affected by the stratification, with characteristic narrowings within the packed layers and large bulbous widenings, carved by the pipes as they emerge from the layers. We have connected the appearance of these forms to the interplay of flow and dissolution in the pipes. In particular, we have argued that as the pipe emerges from a packed layer, the flow rate around its tip resembles that produced by a point source near a wall. Large lateral currents result then in a dramatic widening of the pipe. In two dimensions, the radius of a cavity thus formed is, to a first approximation, inversely proportional to the square root of the Damkölher number.

The results presented here can help in the interpretation of the natural patterns, such as those in Figure 1. For example, if one can measure the porosity and permeability contrast between the layers, as well as estimate the reaction rate based on the petrography of the rock, then the diameters of the pipes in the widest point can be used to estimate the flow conditions under which these structures were formed. Conversely, if the flow conditions can be guessed then these data can be used to calculate the permeability contrast in the layers.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This manuscript was edited by Bjornar Sandnes; a Topic Editor of the Research Topic Flow and Transformations in porous media, in the absence of bias and to the high standards of ethical conduct expected by Frontiers.

Acknowledgments

The authors benefited from discussions with Tony Ladd, Andrew Farrant, Sebastian Schmidt, Elena Abarca, Agnieszka Budek, Paweł Kondratiuk and Michał Pecelerowicz. Tony Ladd is also acknowledged for his finite-difference code dissol, which has been used in the simulations. Sebastian Schmidt is gratefully acknowledged for the photo and helpful discussion. We also thank Ken Grimes and Matej Lipar for their invaluable help with the exploration of solution pipes in southern and western Australia. KP acknowledges funding from the European Union Seventh Framework Programme for research, technological development and demonstration under grant agreement no 316889. PS acknowledges the support from National Science Centre (Poland) under research Grant No. 2012/07/E/ST3/01734.

References

1. Jennings JN. Karst Geomorphology. Oxford: Blackwell (1985).

Google Scholar

2. Palmer, Arthur, N. Origin and morphology of limestone caves. Geol Soc Am Bull. (1991) 103:1–21.

3. Jamtveit B, Meakin P. Growth, Dissolution and Pattern Formation in Geosystems. Dordrecht: Springer (1999).

Google Scholar

4. Chadam D, Hoff D, Merino E, Ortoleva P, Sen A. Reactive infiltration instabilities. IMA J Appl Math. (1986) 36:207–21.

Google Scholar

5. Ortoleva P, Chadam J, Merino E, Sen A. Geochemical self-organization II: the reactive-infiltration instability. Am J Sci. (1987) 287:1008–40.

Google Scholar

6. Hoefner ML, Fogler HS. Pore evolution and channel formation during flow and reaction in porous media. AIChE J. (1988) 34:45–54.

Google Scholar

7. Walsh P, Morawiecka-Zacharz I. A dissolution pipe palaeokarst of mid-Pleistocene age preserved in Miocene limestones near Staszów, Poland. Palaeogeogr Palaeoclimatol Palaeoecol. (2001) 174:327–50. doi: 10.1016/S0031-0182(01)00317-0

CrossRef Full Text | Google Scholar

8. De Waele J, Lauritzen SE, Parise M. On the formation of dissolution pipes in quaternary coastal calcareous arenites in mediterranean settings. Earth Surf Proc Land. (2011) 36:143–57. doi: 10.1002/esp.2022

CrossRef Full Text | Google Scholar

9. Yehle LA. Soil tongues and their confusion with certain indicators of periglacial climate [Wisconsin]. Am J Sci. (1954) 252:532–46.

Google Scholar

10. Rowan G. Theory of acid treatment of limestone formations. J Inst Pet. (1959) 45:321.

11. Fredd CN, Fogler HS. Influence of transport and reaction on wormhole formation in porous media. AIChE J. (1998) 44:1933–49.

Google Scholar

12. Panga M, Ziauddin M, Balakotaiah V. Two-scale continuum model for simulation of wormhole formation in carbonate acidization. AIChE J. (2005) 51:3231–48. doi: 10.1002/aic.10574

CrossRef Full Text | Google Scholar

13. Cohen C, Ding D, Quintard M, Bazin B. From pore scale to wellbore scale: impact of geometry on wormhole growth in carbonate acidization. Chem Eng Sci. (2008) 63:3088–99. doi: 10.1016/j.ces.2008.03.021

CrossRef Full Text | Google Scholar

14. Maheshwari P, Ratnakar RR, Kalia N, Balakotaiah V. 3-D simulation and analysis of reactive dissolution and wormhole formation in carbonate rocks. Chem Eng Sci. (2013) 90:258–74. doi: 10.1016/j.ces.2012.12.032

CrossRef Full Text | Google Scholar

15. Ishibashi T, McGuire T, Watanabe N, Tsuchiya N, Elsworth D. Permeability evolution in carbonate fractures: competing roles of confining stress and fluid pH. Water Resour Res. (2013) 49:2828–42. doi: 10.1002/wrcr.20253

CrossRef Full Text | Google Scholar

16. Daccord G. Chemical dissolution of a porous medium by a reactive fluid. Phys Rev Lett. (1987) 58:479–82.

PubMed Abstract | Google Scholar

17. Daccord G, Lenormand R. Fractal patterns from chemical dissolution. Nature (1987) 325:41–3.

Google Scholar

18. Ewers RO. Cavern Development in the Dimensions of Length and Breadth. Ph.D. thesis, McMaster University (1982).

19. Groves CG, Howard AD. Early development of karst systems. I. Preferential flow path enlargement under laminar flow. Water Resour Res. (1994) 30:2837–46.

Google Scholar

20. Hanna RB, Rajaram H. Influence of aperture variability on dissolutional growth of fissures in karst formations. Water Resour Res. (1998) 34:2843–53.

Google Scholar

21. Cheung W, Rajaram H. Dissolution finger growth in variable aperture fractures: role of the tip-region flow field. Geophys Res Lett. (2002) 29:2075. doi: 10.1029/2002GL015196

CrossRef Full Text | Google Scholar

22. Worthington S R., Ford D. C. Self-organized permeability in carbonate aquifers. Groundwater (2009) 47:326–36. doi: 10.1111/j.1745-6584.2009.00551.x

CrossRef Full Text

23. Szymczak P, Ladd AJC. The initial stages of cave formation: beyond the one-dimensional paradigm. Earth Planet Sci Lett. (2011) 301:424–32. doi: 10.1016/j.epsl.2010.10.026

CrossRef Full Text | Google Scholar

24. Budek A, Szymczak P. Network models of dissolution of porous media. Phys Rev E (2012) 86:056318. doi: 10.1103/PhysRevE.86.056318

CrossRef Full Text | Google Scholar

25. Chaudhuri A, Rajaram H, Viswanathan H. Early-stage hypogene karstification in a mountain hydrologic system: a coupled thermohydrochemical model incorporating buoyant convection. Water Resour Res. (2013) 49:5880–99. doi: 10.1002/wrcr.20427

CrossRef Full Text | Google Scholar

26. Polak A, Elsworth D, Liu J, Grader AS. Spontaneous switching of permeability changes in a limestone fracture with net dissolution. Water Resour Res. (2004) 40:W03502. doi: 10.1029/2003WR002717

CrossRef Full Text | Google Scholar

27. Noiriel C, Gouze P, Made B. 3D analysis of geometry and flow changes in a limestone fracture during dissolution. J Hydrol. (2013) 486:211–23. doi: 10.1016/j.jhydrol.2013.01.035

CrossRef Full Text | Google Scholar

28. Luquot L, Rodriguez O, Gouze P. Experimental characterization of porosity structure and transport property changes in limestone undergoing different dissolution regimes. Transport Porous Med. (2014) 101:507–32. doi: 10.1007/s11242-013-0257-4

CrossRef Full Text | Google Scholar

29. Garcia-Rios M, Luquot L, Soler JM, Cama J. Influence of the flow rate on dissolution and precipitation features during percolation of CO₂-rich sulfate solutions through fractured limestone samples. Chem Geol. (2015) 414:95–108. doi: 10.1016/j.chemgeo.2015.09.005

CrossRef Full Text | Google Scholar

30. Luquot L, Gouze P. Experimental determination of porosity and permeability changes induced by injection of CO₂ into carbonate rocks. Chem Geol. (2009) 265:148–59. doi: 10.1016/j.chemgeo.2009.03.028

CrossRef Full Text | Google Scholar

31. Ellis BR, Peters CA, Fitts JP, Bromhal GS, McIntyre D, Warzinski R, et al. Deterioration of a fractured carbonate caprock exposed to CO₂-acidified brine flow. Greenh Gases Sci Technol. (2011) 1:248–60. doi: 10.1002/ghg.25

CrossRef Full Text | Google Scholar

32. Elkhoury JE, Ameli P, Detwiler RL. Dissolution and deformation in fractured carbonates caused by flow of co₂-rich brine under reservoir conditions. Int J Greenh Gas Control (2013) 16:S203–15. doi: 10.1016/j.ijggc.2013.02.023

CrossRef Full Text | Google Scholar

33. Hao Y, Smith M, Sholokhova Y, Carroll S. CO₂ -induced dissolution of low permeability carbonates. Part II : numerical modeling of experiments. Adv Water Res. (2013) 62:388–408. doi: 10.1016/j.advwatres.2013.09.009

CrossRef Full Text | Google Scholar

34. Deng H, Ellis BR, Peters CA, Fitts JP, Crandall D, Bromhal GS. Modifications of carbonate fracture hydrodynamic properties by CO₂-acidified brine flow. Energy Fuels (2013) 27:4221–31. doi: 10.1021/ef302041s

CrossRef Full Text | Google Scholar

35. Carroll S, Hao Y, Smith M, Sholokhova Y. Development of scaling parameters to describe CO₂-rock interactions within Weyburn-Midale carbonate flow units. Int J Greenh Gas Control (2013) 16:S185–93. doi: 10.1016/j.ijggc.2012.12.026

CrossRef Full Text | Google Scholar

36. Smith M, Sholokhova Y, Hao Y, Carroll S. CO₂-induced dissolution of low permeability carbonates. Part I: characterization and experiments. Adv Water Resour. (2013) 62:370–87. doi: 10.1016/j.advwatres.2013.09.008

CrossRef Full Text | Google Scholar

37. Elkhoury J, Detwiler R, Ameli P. Can a fractured caprock self-heal? Earth Planet Sci Lett. (2015) 417:99–106. doi: 10.1016/j.epsl.2015.02.010

CrossRef Full Text | Google Scholar

38. Ott H, Oedai S. Wormhole formation and compact dissolution in single-and two-phase CO₂-brine injections. Geophys Res Lett. (2015) 42:2270–6. doi: 10.1002/2015GL063582

CrossRef Full Text | Google Scholar

39. Szymczak P, Ladd AJC. A network model of channel competition in fracture dissolution. Geophys Res Lett. (2006) 33:L05401. doi: 10.1029/2005GL025334

CrossRef Full Text | Google Scholar

40. Herwitz SR. Stemflow influences on the formation of solution pipes in Bermuda eolianite. Geomorphology (1993) 6:253–71.

Google Scholar

41. Lipar M, Webb JA, White SQ, Grimes KG. The genesis of solution pipes: evidence from the middle–late pleistocene bridgewater formation calcarenite, Southeastern Australia. Geomorphology (2015) 246:90–103. doi: 10.1016/j.geomorph.2015.06.013

CrossRef Full Text | Google Scholar

42. Sun Q, Cartwright J, Wu S, Chen D. 3D seismic interpretation of dissolution pipes in the South China Sea: genesis by subsurface, fluid induced collapse. Mar Geol. (2013) 337:171–81. doi: 10.1016/j.margeo.2013.03.002

CrossRef Full Text | Google Scholar

43. Golfier F, Zarcone C, Bazin B, Lenormand R, Lasseux D, Quintard M. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium. J Fluid Mech. (2002) 457:213–54. doi: 10.1017/S0022112002007735

CrossRef Full Text | Google Scholar

44. Golfier F, Bazin B, Lenormand R, Quintard M. Core-scale description of porous media dissolution during acid injection - Part I: theoretical development. Comput Appl Math (2004) 23:21. doi: 10.1590/S0101-82052004000200005

CrossRef Full Text | Google Scholar

45. Upadhyay VK, Szymczak P, Ladd AJC. Initial conditions or emergence: what determines dissolution patterns in rough fractures? J Geophys Res Solid Earth (2015) 120:6102–21. doi: 10.1002/2015jb012233

CrossRef Full Text | Google Scholar

46. Krug J. Origin of scale invariance in growth processes. Adv Phys. (1997) 46:139–282.

Google Scholar

47. Budek A, Garstecki P, Samborski A, Szymczak P. Thin-finger growth and droplet pinch-off in miscible and immiscible displacements in a periodic network of microfluidic channels. Phys Fluids (2015) 27:112109. doi: 10.1063/1.4935225

CrossRef Full Text | Google Scholar

48. Roy A, Roy S, Bhattacharyya AJ, Banerjee S, Tarafdar S. Discrete scale invariance in viscous fingering patterns. Eur Phys J B (1999) 12:1–3. doi: 10.1007/s100510050966

CrossRef Full Text | Google Scholar

49. Couder Y, Maurer J, González-Cinca R, Hernández-Machado A. Side-branch growth in two-dimensional dendrites. I. experiments. Phys Rev E (2005) 71:31602. doi: 10.1103/PhysRevE.71.031602

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Huang Y, Ouillon G, Saleur H, Sornette D. Spontaneous generation of discrete scale invariance. Phys Rev E (1997) 55:6433–47.

Google Scholar

51. Howard AD. The development of karst features. Bull Natl Spel Soc. (1963) 25:45–65.

Google Scholar

52. Palmer AN. Cutters and pinnacles in the Salem limestone of Indiana. In: Ginés A, Knez M, Slabe T, Dreybrodt W, editors. Karst Rock Features: Karren Sculpturing. Ljubljana: ZRC Publishing (2009). pp. 349–58.

53. Steefel CI. Geochemical kinetics and transport. In: Brantley SL, Kubicki JD, White AF, editors. Kinetics of Water-Rock Interaction. New York, NY: Springer (2008). pp. 545–89.

54. Dreybrodt W. The role of dissolution kinetics in the development of karst aquifers in limestone: a model simulation of karst evolution. J Geol. (1990) 98:639–55.

Google Scholar

55. Kalia N, Balakotaiah V. Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks. Chem Eng Sci. (2007) 62:919–28. doi: 10.1016/j.ces.2006.10.021

CrossRef Full Text | Google Scholar

56. Grimes KG. Solution pipes and pinnacles in syngenetic karst. In: Ginés A, Knez M, Slabe T, Dreybrodt W, editors. Karst Rock Features. Karren Sculpturing. Ljubljana: ZRC Publishing (2009). pp. 513–26.

57. Schaetzl R, Anderson S. Soils: Genesis and Geomorphology. New York, NY: Cambridge University Press (2005).

Google Scholar

58. Kondratiuk P, Szymczak P. Steadily translating parabolic dissolution fingers. SIAM J Appl Math. (2015) 75:2193–213. doi: 10.1137/151003751

CrossRef Full Text | Google Scholar

59. Adriano DC, Havas M, editors. Acidic Precipitation. New York, NY: Springer (1989).

Google Scholar

60. Szymczak P, Ladd AJC. Reactive infiltration instabilities in rocks. Part 2. Dissolution of a porous matrix. J Fluid Mech. (2014) 738:591–630. doi: 10.1017/jfm.2013.586

CrossRef Full Text | Google Scholar