Solution pipes and focused vertical water flow: Geomorphology and modelling

Focused flow forms distinctive features termed solution pipes (vertical cylindrical voids) in the epikarst zone of porous calcareous rocks with medium to high matrix porosity. Solution pipes vary in size, but are generally less than 1 m in diameter with variable depths, the deepest reaching 100 m. Most occur in eogenetic rocks, particularly in Quaternary calcareous sandstones (calcarenites), but some pipes are also reported in Miocene-Pliocene and Cretaceous chalk and calcarenites. Most of the published research favours their dissolutional origin, but offers different theories regarding the initial trigger that focused the water flow in the host rock matrix, including vegetation patterns and rock heterogeneities. Here we review the current state of knowledge regarding geomorphology of solution pipes and discuss in greater detail the potential role of dynamic instabilities and self-organisation in the emergence of focusing patterns. We review theoretical and numerical studies which allow estimation of interpipe distances, distribution of pipe lengths, extent of pipe merging and the influence of rock porosity on their shapes. In theory, linking the numerical predictions with field measurements might allow the use of pipe characteristics as potential paleo-climatic indicators, provided the problems of dating can be overcome.


Introduction
Focused flow of water, also often referred to as finger(ed) flow and preferential flow, is a physical process in which water flow becomes organised in preferred pathways.In a wider sense, this includes river systems and ocean currents, but the term is mostly used to describe the effects of non-uniform flow through porous media such as soil, sand and rocks.
In karst geomorphology, focused flow causes inhomogeneous corrosion, which results in the formation of landforms such as cave systems and dolines.Understanding of the initial trigger of focused flow, its characteristics, parameters and timing, can give us insight into the origin of landforms, and consequently, how those landforms interact to form the overall landscape.Differentiating between the Earth processes that are driven by topography and geology, and those driven by climate is fundamental for predictions of the landscape response to climatic changes.
In physics, focusing of flow is one of the prime examples of a self-organising phenomenon, in which the interplay between different processes leads to the formation of qualitatively new patterns of organisation.Flow focusing emerges in a variety of different hydrodynamic instabilities, such as viscous fingering, convective instabilities or reactive-infiltration instability.The initial stages of these processes, with the formation of sinusoidal perturbations in an initially uniform system, are now well understood (Chandrasekhar, 1961;Ortoleva, 1994;Pelcé, 2004).Much less is known, however, about the nonlinear regime, when the initial perturbations of the interface are transformed into finger-like structures that advance into the system.Interestingly, over time these emergent structures often attain well-defined, invariant shapes, like the Saffman-Taylor finger in viscous fingering or the Ivantsov paraboloid in dendritic growth (Ivantsov, 1947;Saffman and Taylor, 1958;Horvay and Cahn, 1961).It can therefore be expected that such invariant forms might also be present in geological/geomorphological systems.Self-focusing of water flow into evolving linear conduits is one of the distinctive features of epigene karst and is strongly dependent on the diagenetic stage of the rock and its porosity.In hard telogenetic limestones the focusing is initially influenced by cracks, joints or bedding partings.However, in porous eogenetic limestones the focusing of flow into vertical linear paths initially originates from diffuse flow through the matrix/inter-granular porosity of the sediment.This results in one of the most distinctive features of such limestones -solution pipes, which are vertical or near vertical finger-like dissolution structures of various sizes, occurring predominantly within the epikarst zone of limestones with high permeability and matrix porosity (mostly Quaternary and Late Neogene calcarenites).
Whilst there is a consensus that focused water flow is responsible for solution pipe formation, there is no agreement as to the initial triggers of the focused flow in limestones with matrix porosity.What is clear, however, is thatsince the medium is solublefocusing of the flow is followed by enhanced dissolution, which eventually transforms initial inhomogeneities into mature pipes (Figs.1-2).Solution pipes are thus prime examples of the positive feedback between focused flow and dissolution.
The objective of this paper is to review the current knowledge of solution pipes and discuss in greater detail their origin, with most attention on focused flow and the role of dynamic instabilities and selforganisation.We attempt to bring together the numerical models from physical sciences and a geomorphological perspective on solution pipe formation in the hope of providing a more comprehensive insight into the phenomenonso far, these two elements have been discussed separately by different research communities.In particular, quantification of the link between shapes and growth rates of the pipes and external conditions such as precipitation or temperature should ultimately allow one to use the pipes as paleoclimate markers.

Solution pipes analysis methodology
The different methods commonly used in solution pipe analysis are reviewed with emphasis on their strengths and limitations, allowing for a more in-depth, critical interpretation of the results.

Morphometrics
The morphometry of solution pipes has received a great deal of attention and is typically based on manual measurements of pipe size and shape in the field, which are then used to calculate areal density and percentage cover (Herwitz, 1993).Large studies that measured over pipes in Poland (Walsh and Morawiecka-Zacharz, 2001) and over pipes in Italy (De Waele et al., 2011) had sufficient data for the calculation of statistics and plotting of logarithmical distributions.Grimes (2004) provided 2-D statistical data based on 5 × 5 m area in Australia which included inside and outside diameter of the pipes, their elongation, the thickness of cemented rim, and distance and direction to nearest neighbour.Advanced morphometrical data were reported from pipes in Turkey ( Öztürk et al., 2018) and Ukraine (including some Polish sites) (Dobrowolski and Mroczek, 2015).In most cases only partial information could be obtained on the pipe characteristics, because the outcrops provide a 2-D section of a threedimensional system.Accurate pipe lengths may be nearly impossible to measure, since they often extend below the outcrop or out of the observed plane.Similar problems arise while quantifying the pipe distribution, since they are often covered by soil.However, it is possible to observe their 3-D shapes in coastal calcarenites where the rimmed pipes often form an inverted landscape, with the rock in between them fully eroded (Fig. 1C).

Geological analyses
The stratigraphy, mineralogy, fabric and geochemistry of the surrounding rock and pipe rims and fills provide valuable insights into primary depositional conditions of the rock and the processes during pipe formation.Basic mineralogy, fabric and porosity are examined using geological thin sections under petrographic or SEM microscopes, with point counting for grain abundance (Lundberg and Taggart, 1995).X-ray diffraction (XRD) has been often used to analyse the mineralogy of the rims, host rock and fill of the pipes (Lipar et al., 2015), with X-ray fluorescence (XRF) used for determining major-element concentrations (Herwitz, 1993).Stable carbon and oxygen isotopes provide useful data (Lundberg and Taggart, 1995), which can partly be correlated with the age of the host rock (Lipar et al., 2015).
To constrain the age of pipe formation, the host rock can be dated by thermoluminescence (TL) (White, 2000) and optically stimulated luminescence (OSL) (Lipar et al., 2015), providing a maximum age.Pipe infills and the cement in solution pipe rims can also be dated by TL and 14 C respectively (Dobrowolski and Mroczek, 2015;Marsico et al., 2003), providing a minimum age, but the origin and deposition of the fill may be unclear and the presence of allochems in the cemented rims makes the dating challenging.

Mathematical modelling and analogue experiments
From a physical point of view a solution pipe is a relatively simple system (reactive flow through a nearly uniform rock matrix), but mathematical and numerical models of pipe formation are relatively rare.Statistical analyses of field measurements of solution pipe morphology were performed by De Waele et al. (2011), andUpadhyay et al. (2015).On the other hand Petrus and Szymczak (2016) and Kondratiuk et al. (2017) used numerical models to investigate solution pipe growth.
Dissolution structures in rocks have been studied by analogue experiments, e.g. in a mixture of glass balls and salt (Kelemen et al., 1995), quartzitic carbonate sandy sediment (Cha and Santamarina, 2016) or gypsum (Osselin et al., 2016) (Fig. 4).Even though experiments were small-scale, the resulting dissolution structures share a lot of similarities with solution pipes, particularly shorter ones (cf.Fig. 1A, B and D).

Terminology
The terminology of karst features in young porous carbonates with high matrix porosity has been developed independently by different researchers and consequently several different terms for solution pipes exist.Day (1928) originally referred to solution pipes simply as "pipes".The term "solution pipe" was later established by Jennings (1968) and Jennings (1985).Lundberg and Taggart (1995) raised concerns that the term "solution pipe" has been used for dissolution features of differing scales and forms, and therefore suggested "dissolution pipe".Solution pipes have also been named "palmetto stumps", "chimneys", "vertical pipes", "soil roots", "solution pits", "makondos" and "organ pipes" (reviewed in Lundberg and Taggart (1995)).On the Bahamian Archipelago, they have been generally described as "vadose fast flow routes" and termed "pit caves" when they eventually reached large sizes (Harris et al., 1995).Both "solution pipe" and "dissolution pipe" have been widely used, but the term "solution pipe" is used throughout this article, as it has precedence.
Solution pipes occur in a variety of present-day climatic conditions, including tropical, subtropical, Mediterranean, temperate and continental.However, the present climatic conditions are not necessarily those in existence at the time of pipe formation, and rather reflect the localities of carbonates with matrix porosity which are prone to solution pipe formation.
Pipes can occur either on the present surface of young calcareous deposits (e.g., calcarenites), or they may be preserved on older (now buried) surfaces and descend from palaeosol horizons in stratigraphically older calcarenites (palaeokarst), as seen in Australia (Grimes, 2004;Lipar et al., 2015).Most of the Miocene and Cretaceous rocks with observable solution pipes are overlain by younger, generally fluvial and glacial deposits (Walsh and Morawiecka-Zacharz, 2001;Dobrowolski and Mroczek, 2015;Willems and Rodet, 2018).

Physical characteristics
The radii and depths of solution pipes vary enormously; they can be wider than 2 m and deeper than 100 m (Willems and Rodet, 2018).However, the degree of denudation that has occurred since pipe formation can influence what is observed today.For solution pipes formed in the vadose zone, their maximum depths are determined by either the water-table or sea-level (Coetzee, 1975).Nevertheless, because pipes have formed at times of different sea-level, some pipes now reach below sea-level, e.g. in southwestern Western Australia (Fairbridge, 1950).
The radius of a single pipe within a homogeneous rock is usually either constant or tapers towards the base with a cigar-shaped or coneshaped termination (De Waele et al., 2011).There is a correlation between the shape and lengthlonger pipes are usually more cylindrical in shape, whereas shorter ones are more conical or paraboloidal (see Figs. 3-5).The shapes and sizes of pipes may differ considerably between locations.For example, pipes in the Cretaceous chalk in Normandy are usually longer, wider and less cylindrical than pipes in calcarenites (Figs. 1 & 2).
Dissolutional enlargement and cracking of pipes influence their morphology at the ground surface, resulting in wide funnels or wide elongated karren, yet their cylindrical shape is usually preserved deeper in the host rock.From the speleological perspective, solution pipes large enough to be entered by a human are considered to be vertical caves (Fig. 7C).Coalescence of pipes is common (Marsico et al., 2003;Lipar et al., 2015).
Although individual pipes are not uncommon, they tend to occur in clusters (De Waele et al., 2011;Lipar et al., 2015).Fairbridge (1950) observed that the maximum concentration of pipes coincides with depressions in old dune surfaces, explained by additional surface water input from surrounding slopes, whilst for pipes in Apulia, Italy, Marsico et al. (2003) explained clustering by permeability variations of the overlying colluvium.Pipes in Cretaceous chalk have been observed to occur mainly under a sediment cover, and are absent where the chalk does not have a cover.Similar observations were reported for Cretaceous calcarenites along the Belgian-Dutch border where most of the pipes are present under a gravel terrace of the Meuse river, whereas calcarenites covered by Cenozoic clayey sands do not contain any pipes.These observations indicate that permeability of the overlying cover is an important factor in pipe formation (Juvigné, 1992).
rims are usually smoother and more regular on the inside, becoming progressively more uneven outside towards the surrounding rock.Secondary calcium carbonate encrustations may in some places completely fill the pipes (Coetzee, 1975).
The rim can consist of multiple concentric layers (Coetzee, 1975), generally formed by case-hardening around the pipes, as shown by the presence of allochems from the host rock within the rims; the absence of allochems indicates that laminae grew into an open space between a solution pipe rim and the host rock (Lundberg and Taggart, 1995).The rims are usually less porous and more resistant to weathering than the host rock, and consequently stand out after erosion of the surrounding material so that they may resemble tree trunks (Curran et al., 2008), which sparked the mis-identification of solution pipes as petrified trees (Bretz, 1960;Boutakoff, 1963).
Rims form as a result of carbonate precipitation in the fringes of the pipes from the saturated water from the central part of the pipe (Grimes, 2006); this is favoured by a relatively dry and/or seasonal climate (Lipar et al., 2015).Marsico et al. (2003) attempted 14 C dating of the rims of pipes in Plio-Pleistocene calcarenites, but the result was at the limit of the method (~33 ka), or possibly the samples were contaminated with rock-derived dead carbon.
As subsurface features, solution pipes are generally filled with soil, clay or other loose sediment cover, which subsides into them concurrently with their formation.There are no grooves or other turbulent flow related features on the sides of the pipes, which suggests that there was no free flow along the pipe walls, i.e., the pipe was completely filled with sediment during its formation.Besides an organic component, the sediment may contain solutional residue of the host rock and/or aeolian material (Herwitz, 1993;Herwitz and Muhs, 1995;Grimes, 2009;De Waele et al., 2011;Lipar et al., 2015).Solution pipes in formerly glaciated areas usually contain glacial fill (Dobrowolski and Mroczek (2015).Willems et al. (2007) reported a complex fill containing Oligocene sands, sandy clays, fluvial gravel and/or loam from a loess cover.The presence of sediment within solution pipes favours root growth, sometimes visible when solution pipes reach the ceiling of an underlying cave.
An intriguing hypothesis on the genesis of terra-rossa filled pipes has been put forward by Merino and Banerjee (2008).They proposed that terra-rossa is formed by authigenic replacement of the underlying limestone at a narrow reaction front.The hydrogen ions produced in such reaction further dissolve the rock matrix, triggering a reactiveinfiltration instability which results in intense piping.
The fill may become cemented, occasionally containing concentric calcareous laminae (Grimes, 2004;Lipar et al., 2015), and may protrude out from the ground surface (De Waele et al., 2011;Lipar et al., 2015).Uncemented fill can be washed out when pipes are exposed on the edge of cliffs or in the ceiling of a cave (Figs. 5 & 7).Dobrowolski and Mroczek (2015) dated glaciogenic pipe fills in Cretaceous rocks using the TL method as 350 ka to 190 ka, indicating a Saalian age for the deposits.

Solution pipe formation
Numerous hypotheses have been put forward for the formation of solution pipes.The earliest is that of the French orientalist Renan (1864), who observed solution pipes in Lebanon, north of Beirut, near ancient tombs, and hence viewed them as anthropogenically drilled by auger.
Solution pipes have been described as petrified tree trunks.First discussed by Darwin (1845) based on examples from southwestern Australia, the idea was followed by Livingston (1944), Bretz (1960),

Table 1
Selected location of solution pipes, on which detailed research has focused.The numbers in the first column refer to the numbers in Fig. 3. Diameter: 10 cm to 100 cm (some can exceed 100 cm) Depth: up to 2000 cm (usually observed 100 cm to 500 cm) Quaternary aeolianite Yes (Grimes, 2004;Lipar et al., 2015) Boutakoff ( 1963), Neumann and Hearty (1996), Curran et al. (2008) and Hearty and Olson (2011).Decay of a tree trunk rapidly buried by aeolian sand may create a conduit for groundwater, which can cause precipitation of calcrete within or around the pipe, forming a solid cast of an original tree or a hollow pipe, respectively.This process has been convincingly demonstrated for some pipes in Bermuda and The Bahamas, but is not a general mechanism for solution pipe formation.
Although pipes often contain roots, they are frequently surrounded by undistorted aeolianite cross-bedding (i.e.there is no evidence of an original tree), have a higher density than forest trees of that size, and lack palaeosol horizons around their bases with branching root structures (Grimes, 2004).Coetzee (1975) observed nearly cylindrical pipes on intertidal platforms along the south-eastern African coast and postulated that they could have formed by pothole erosion.However, no grinding tools (e.g., resistant pebbles) were found in any of the solution pipes, and pothole erosion decreases with an increase in depth, so at the depths reported by Coetzee (1975) the erosional effects would become insignificant.Day (1928) and Fairbridge (1950) attributed solution pipe formation to dissolution, and Coetzee (1975) used experiments with HCl to conclude that vertical downward dissolution can eventually form a cylindrical pipe in a calcareous porous rock.De Waele et al. (2011) used the statistical analysis of morphological data of pipes to verify that downward corrosion of infiltrating water can cause pipe formation.Morawiecka and Walsh (1997), Walsh and Morawiecka-Zacharz (2001) and Dobrowolski and Mroczek (2015) studied pipes in Miocene calcarenite in Poland that could have formed subglacially under the cover of a continental glacier till.Piping there is restricted to areas beneath a till cover, and pipes did not develop where till is absent.This links the formation of the pipes to the deglaciation at the end of the Elsterian period, when large volumes of cold water were suddenly released into the fully saturated subglacial aquifer (Flowers, 2015;    8).The soft sediments and rock layers below a glacier are strongly over-pressurised, with a significant downward head gradient that induces vertical groundwater flow (Piotrowski, 2006;Siegert et al., 2018).While discussing subglacial groundwater flow during the last glaciation in a similar location in northwestern Germany, Piotrowski (1997) noted that steep hydraulic gradients at the glacier periphery can result in vertical groundwater flows reaching up to 200 m below the glacier bed with velocities up to 30 times faster than at present.Such flows could result in pipe formation within the soluble bedrock under fully phreatic conditions (Fig. 8).Solution pipe formation would have been enhanced by the relatively high solubility of calcite at low temperatures.
The phreatic solution pipes do not seem to be significantly different morphologically from the more common vadose solution pipes described from elsewhere in the world, in terms of aspect ratio, length distribution, average width or average distance between the pipes.Mathematical modelling of solution pipe formation under a pressure gradient in phreatic conditions (Upadhyay et al., 2015;Petrus and Szymczak, 2016;Kondratiuk et al., 2017) showed pipe morphologies similar to those formed under vadose conditions.Marsico et al. (2003) noticed that pipe formation was most effective when there is little secondary porosity (e.g.joints) to guide the flow paths, and Lipar et al. (2015) showed that solution pipes can form rapidly before consolidation of the rock, so pipe development often predates post-depositional rock lithification and deformation.

Heterogeneities of the media
The traditional view is that focused flow in media with matrix porosity (in which solution pipes form) requires pre-existing heterogeneities to focus the flow in a particular spot.The heterogeneities can be at or near the surface where the water enters the media (e.g., stem flow, micro-topography, hydrophobic areas; Dekker and Ritsema, 1994;Lundberg and Taggart, 1995;Grimes, 2009;Lipar et al., 2015) or within the media (macro-or meso-pores, localised water repellency, variable horizontal permeability or bedding structures, porous patches in developing calcrete; Doerr et al., 2000;Grimes, 2006).
A prominent trigger factor may be stemflow down the trunk of a tree; this represents water intercepted by the leaves and branches (Jennings, 1968;Herwitz, 1993;Crockford and Richardson, 2000;Huang et al., 2005;Johnson and Lehmann, 2006;Lipar et al., 2015).The aggressiveness of stemflow is enhanced by wash-off and leaching of leaves and bark (Johnson and Lehmann, 2006), so stemflow may cause intensive rock dissolution beneath trees (Lundberg and Taggart, 1995;Grimes, 2009;Lipar et al., 2015).Buried tree trunk casts can create conduits for groundwater (Neumann and Hearty, 1996;Hearty and Olson, 2011), as can open channels left by decayed tree roots (Mitchell et al., 1995;Devitt and Smith, 2002).
In addition, due to denudation, the surface where the pipes are observed today is often different from the surface where the pipes originated.Hence, the exact primary surficial features of the host rockwhich could potentially trigger pipingmight be lost.Once initiated, a pipe seems to ignore the primary features of the rock enclosing it, but this observation does not mean that the vanished rock, now denuded, did not have features that focused or additionally directed surficial flow.This suggests that attention should be focused on the areas where solution pipes are just starting to form, but to the authors knowledge, none of the published literature has an answer to "are the pipes forming today and where?".Measurements of the water pH inside the pipes, performed by the authors in several piping outcrops around the world, show that the solution is slightly basic, suggesting that the pipes are no longer active.

Dynamic instabilities and self-organisation
However, for many solution pipes no clear trigger is identifiable.For example, Walsh and Morawiecka-Zacharz (2001), discussing the origin of solution pipes in Smerdyna, Poland, noted that neither "structural fabric of the host-rock (nor) irregularities of the interface between host and cover" determined the locations and forms of individual pipes.
Furthermore, in the rare cases where the pipes are visible in plan view (Fig. 9), their spatial distribution is rather uniform, characterised usually by a well-defined interpipe distance, which calls for a physical theory capable of predicting their diameters and distances between them from first principles.
As a result, theories emphasizing the role of dynamic instabilities and self-organisation have been developed, based on nonlinear systems in earth sciences whereby structures form spontaneously with intrinsic length scales independent of initial conditions (Ortoleva, 1994;Turcotte, 1997;Meakin, 1998;Jamtveit and Meakin, 1999).One instability that can focus flow in porous rock or sediment is that of the wetting front infiltrating the vadose zone.Such a front was shown to break into fingers (Hill and Parlange, 1972;Glass et al., 1989;Selker et al., 1992;Dekker and Ritsema, 1994;Babel et al., 1995;Hendrickx and Flury, 2001;Rezanezhad et al., 2006), which then advance into the waterunsaturated sediment.Fingered flow most often occurs in homogeneous sands with grainsize above 0.5-1.0mm (Diment and Watson, 1985;Glass et al., 1989;Rezanezhad et al., 2006), which may be correlated with the extensive occurrence of solution pipes in calcarenites (0.06-2 mm).Fingered flow is particularly favoured if the medium is water repellent, causing uneven infiltration (Roper, 2004).The width of a finger increases until the water saturation within its core becomes quasi-stable, so the diameter is generally greater with increasing water content (De Rooij and Cho, 1999;Rezanezhad et al., 2006).De Waele et al. ( 2011) and Lipar et al. (2015) already suggested that an unstable wetting front (i.e., fingered flow) could be responsible for solution pipe formation without triggers related to vegetation or other heterogeneities in the rock.However, the occurrence and type of fingered flow in sediments is sensitive to many factors such as initial water content, size and Fig. 8.The groundwater drainage pattern through a subglacial aquifer under a glacier resting on unfrozen (A) and partly frozen (B) bed [after (Piotrowski, 2006)].Note that the potentiometric surface is located above the ground, so that the subglacial aquifer remains fully saturated.distribution of sediment particles, rainfall intensity, and water repellence (De Rooij, 2000;Kawamoto et al., 2004).Two points of caution need to be noted here: first, the most rapidly growing wavelength of the wetting front infiltration instability is usually on the centimetre scale, i. e., significantly smaller than observed distances between the pipes.Second, in subtropical climates a generally impermeable calcrete crust develops relatively quickly on dune calcarenites, and this will impact the infiltration patterns.
The second process which has been invoked in relation to the focused flow is reactive-infiltration instability, which develops when a porous matrix is dissolved by a flowing reactive fluid (Chadam et al., 1986;Ortoleva et al., 1987;Glass and Nicholl, 1996;De Rooij, 2000;Wangen, 2013;Szymczak and Ladd, 2014).The underlying mechanism of this instability is the positive feedback between spatial variations in porosity in the initial matrix and the local dissolution rate.A small enhancement in porosity at some point in the reaction front increases the fluid flow in that region, which brings the unsaturated solution deeper into the undissolved rock matrix.By this means, any local variation in porosity is amplified as the reaction front moves downward, eventually developing into fingers.Petroleum engineers were first to notice that the dissolution around an acidized well is often nonuniform and the flow becomes spontaneously localised in pronounced dissolution channels (Rowan, 1959;Hoefner and Fogler, 1988;Panga et al., 2005).They have learned how to use these channels to their advantagethe dissolution channels transport flow in an effective manner while not requiring a lot of reactant for their creation, hence providing an efficient way of increasing the permeability of the rock.These instabilities are not only relevant to acidized wells, but can take place in almost any system in which a surface reaction is coupled with the fluid flow, hence the name reactiveinfiltration instability (Chadam et al., 1986).Ortoleva (1994) in his book on geochemical self-organisation pointed to reactive-infiltration instability as one of the major pattern-forming processes in earth systems.It has been associated with formation of karst conduits (Hanna and Rajaram, 1998;Szymczak and Ladd, 2011), replacement fingers (Kondratiuk et al., 2017;Beaudoin et al., 2018) and conduits of upwelling melt in the mantle (Aharonov et al., 1995;Spiegelman et al., 2001).Reactive-infiltration instability has been demonstrated in laboratory experiments, both in porous cores (Daccord and Lenormand, 1987;Hoefner and Fogler, 1988;Kelemen et al., 1995;McDuff et al., 2010) and in fractures (Detwiler et al., 2003;Garcia-Rios et al., 2015;Osselin et al., 2016).
Interestingly, a very similar instability will be present if dissolution is replaced by mechanical erosion, provided that the increase of porosity (due to erosion) leads to greater permeability and thus higher flow.This process is probably responsible for the formation of pipes in soil (one centimetre to several metres in diameter and depth), with examples also known in quartz beach and dune sands, e.g.Löffler (1974), Dekker and Ritsema (1994) and Thompson (1992).Dams and levees are vulnerable to failure due to internal erosion of this type (Bonelli, 2013).Importantly, there must be an efficient means of transporting away the eroded material, which otherwise would simply clog the pores, halting the process.This is why soil piping can never take place from above, as is the case for solution pipes.Instead, the material is removed at the outlet of the system and erosion continues upstream, resulting in the formation of pipes in the soil matrix (Clibborn and Beresford, 1902;Robbins and Van Beek, 2015;Van Beek et al., 2015;Robbins et al., 2018).This can be regarded as backward piping, and is in many ways a mirror of reactiveinfiltration instability (forward piping; Fig. 10); in the latter the flow is focused in the upstream part and becomes uniform downstream, whereas in backward erosion the focusing takes place downstream and the upstream part remains uniform.Both situations can lead to the formation of pipe-like structures of similar shapes, albeit advancing in different directions.
Forward (chemical) erosion has been suggested as the main mechanism underlying the formation of karst conduit networks (Hanna and Rajaram, 1998;Szymczak and Ladd, 2011), whereas backward (mechanical) erosion is responsible for the growth of seepage channel networks (Devauchelle et al., 2012;Cohen et al., 2015).Interestingly, flank margin caves seem to be an exception to this rule, with an erosion model suggesting upstream development (Mylroie et al., 1990).
The different instabilities discussed in this section are not mutually exclusive and may operate concurrently.For example, wetting front fingers can trigger mechanical erosion in the soil and chemical dissolution in the underlying limestone.Alternatively, chemical erosion within pipes can be accompanied by mechanical erosion as flow rates increase.

Numerical simulations versus natural examples
Modelling of solution pipe formation requires solving coupled equations for groundwater flow, chemical transport and porosity evolution.Numerical modelling of solution pipe formation in subglacial phreatic conditions using Darcy's law (governing the meltwater flow; Ravier and Buoncristiani, 2018) has been conducted by Kondratiuk and Szymczak (2015) and Petrus and Szymczak (2016).The details of this model are presented below.On the other hand, pipes from in subtropical and Mediterranean climates are usually associated with vadose conditions (Lipar et al., 2015), but numerical modelling of these conditions has not yet been undertaken, and remains a future challenge (see Section 5.1).In any case, the morphological similarities between phreatic and vadose solution pipes, as previously discussed, suggest that those conditions might not be the main factor controlling solution pipe formation and evolution.

The governing equations
According to Darcy's law, the rate of groundwater flow (u) through a water-saturated porous medium is proportional to the gradient of the hydraulic head (∇h) where K is the permeability, μ the viscosity of the fluid and φ the porosity.The head itself (h) is expressed as where p is the fluid pressure, g the gravitational constant, ρ the fluid density and z the elevation.The Carman-Kozeny relation is usually adopted for the permeability dependence on the porosity, K(φ) ~ φ 3 .It is also assumed that the Darcy velocity field is incompressible: neglecting contributions to the fluid density from reactants or dissolved products.The chemical erosion equation, which gives the increase of the porosity (φ) as a result of limestone dissolution, is where c is the concentration of calcium ions, which controls the dissolution kinetics of limestone at a pH above 5, typical for nearly all epigene karst groundwater (Plummer and Wigley, 1976;Palmer, 1991).Next, c sat is the saturation concentration at which the reaction stops, k is the kinetic rate and s is the specific reactive surface area per unit volume.In general, the specific reactive surface area changes in the course of the dissolution.It can both increase (as pores become larger) and decrease (as grains become smaller).In Petrus and Szymczak (2016) these changes are neglected and the specific surface area is approximated as a constant (s = s 0 ).The concentration of the calcium ions is calculated through the convection-diffusion-reaction equation where D is the diffusion coefficient.Under typical geophysical conditions, dissolution is slow in comparison to flow and transport processes;,which warrants the assumption of a steady state in both the flow and transport equations.In the numerical implementation, Eqs. ( 1)-( 4) are solved sequentially using finite-difference methods.The flow Eqs. ( 1)-( 2) are solved first, followed by the transport Eq. ( 4).Then, the porosity field is updated in proportion to the local dissolution rate (Eq.( 3)) and the whole process is repeated.

Instability wavelength and the average distance between solution pipes
The theory of reactive-infiltration instability allows us to estimate the initial wavelength of the instability, which should be comparable to the average distance between solution pipes.Since the groundwater velocities in the porous matrix are relatively small, the reactant penetration lengths (defined as the distance over which the infiltrating solution gets saturated) are much smaller than the diffusive length l = D/u, characterising the decay of the concentration in the upstream region.In this case, the so-called thin-front approximation can be adopted (Ladd and Szymczak, 2017).The instability wavelength is then dependent primarily on the flow rate and, to a lesser extent, on the permeability contrast between the rock matrix and the overlaying layers of soil or clay.Typical groundwater velocities are relatively small (10 − 6 -10 − 5 cm/s).For the reaction rate corresponding to the limestone dissolution by karstic waters (Dreybrodt, 1990) and typical reactive surface areas of calcareous sand (Sterianos, 1988) the theory predicts an instability wavelength between tens of centimetres and tens of metres (Szymczak and Ladd, 2013), which is consistent with the average distance between the pipes in different parts of the world.However, these findings need to be interpreted with caution, sinceas we discuss belowthe pipes can merge, thus increasing the characteristic scale of the observed features; both the diameters and the distances between the active (growing) pipes will tend to increase during merging.
For pipes which formed sub-or peri-glacially under the cover of a continental glacier till, the situation is somewhat different, since the groundwater velocities can then be relatively high, up to 4 m per day (Piotrowski, 1997).Taking pipes in Smerdyna as an example, one can assess the total meltwater fluxes which the area had to receive during the deglaciation for such an extensive pipe network to be formed (Kondratiuk, 2017).The pipes in Smerdyna cover about 6% of the surface area, with an average volume of the pipe of the order of 1 m 3 (Morawiecka and Walsh, 1997) This gives a host rock removal volume (per square meter) of V w =0.16 m 3 /m 2 .Assuming an open system with an atmospheric pCO 2 and T = 0 • C, we get the saturation concentration of dissolved calcite of 80 mg/l (Palmer, 1991).Taking into account the density of calcite, ρ = 2.71 g/cm 3 , and the porosity of the bedrock, φ ≈ 50%, the amount of water (per square meter) needed to form the pipe can be estimated to be: which corresponds to a water column height of 2.7 km.This is more than the estimated thickness of the Elsterian ice cover (~500 m), which suggests that the area was infiltrated by meltwater originating form a local glacial drainage system or subglacial lake.Thus, local hydrology was important for the pipe creation, which partially explains why they appear only in certain places, even though the host matrix rock lithology remains the same over a much larger area.The most uncertain factor in pipe genesis is the total time, T, over which the water was infiltrating the rocks, which can be anything from 10 to about 1000 years based on the typical icesheet retreat times reported for recent deglaciations.This gives groundwater fluxes of the order of V w /T ≈ 0.07-0.7 m/day, well within the realistic range of subglacial groundwater velocities (Piotrowski, 1997).
One important prediction of the model is that fingers do not form if the flow rates are too low, since diffusion is then very effective in smoothing out the solute concentration gradients and the front remains stableonly surface dissolution is observed.Thus if the climate is too dry, solution pipes will not form, confirming the absence or minimal formation of solution pipes during the relatively dry periods of marine isotope stages (MIS) 9-8 in southern Australia (Lipar et al., 2015;Lipar et al., 2017).

Hierarchical distribution of pipe length
The fingered pattern evolves over time, driven by two processes: competition between the fingers for the flow and their merging.Competition is driven by a similar instability to that which lies behind finger formationif one finger is slightly longer than another, the density of the flow lines at its tip increases (Fig. 12) and thus the total flow in the finger increases.This is similar to the phenomenon which makes a lightning rod worka long and thin object concentrates the field lines at its tip and screens the surrounding regions from the high voltage.Higher flow drives the undersaturated solution deeper into the rock matrix, and the pipe grows more quickly at the expense of the shorter ones in its immediate neighbourhood.The process then repeats itself, leading to the appearance of a scale-invariant, hierarchical growth pattern, with a large number of short pipes and just a few longer ones, as indeed seen in nature, particularly in places where pipes do not have cemented rims (e.g., Fig. 12).
Numerical simulations (Petrus and Szymczak, 2016) of piping seem to suggest that the distribution of pipe lengths should follow a power law, N(L)~L − α where N denotes a number of pipes longer than L; the exponent in this power law is ~1 for a homogeneous system.Exponents close to unity have also been reported for formation of dissolution fingers in fractures and fractal wormhole formation in acidization of porous rocks (Budek and Szymczak, 2012;Upadhyay et al., 2015).
Analysis of length distribution data collected by De Waele et al. ( 2011) for pipes in the coastal areas of the Mediterranean, if it is restricted to fully visible pipes for which the full length can be measured (Upadhyay et al. (2015)), showed that the distribution fits a power-law model with an exponent of 1.05, close to the value suggested by the modelling (Petrus and Szymczak, 2016).
In Fig. 12 we present exposures of clusters of solution pipes (in Canunda National Park, South Australia and in Smerdyna, Poland), which show hierarchical distribution of lengths with a large number of short pipes and just a few longer ones.We supplement this figure with a panel illustrating the growth of the pipes in a two-dimensional numerical simulation, which shows how the scale-invariant structure is formed.Initially, small perturbations of the dissolution front appear, with the wavelength controlled by the reactive-infiltration instability.Subsequently, however, longer pipes screen the shorter ones off and a hierarchical structure emerges.Fig. 11.A -Solution pipes in Pleistocene aeolianite at Burns Beach, Perth, Australia; Bsolution pipes in Pleistocene aeolianite at Cape Bridgewater, Victoria, Australiaboth formed in Pleistocene aeolianites; Cnumerical simulations of the growth of solution pipes.Note that inclusion of cementation (i.e., cemented rims on both photographs) weakens the interaction between the pipes, which leads to larger density of relatively long pipes of similar length).
In addition, dissolution in the longer pipes may be enhanced by the additional soil and consequently plant activity within these pipes, resulting in additional CO 2 and organic acids increasing the aggressiveness of the water (Ford and Williams, 2007).
Nevertheless, some exposures appear to show solution pipes with relatively uniform lengths (e.g., Fig. 11A & B, which seems different from the field examples and modelling results in Fig. 12).Interestingly, in many of these cases the pipes have low porosity cemented rims.There are two different concepts in the literature regarding rim formation (Grimes, 2004;Lipar et al., 2015).The rims may form in drier periods in geologic time, which follow the wetter periods when the pipes themselves are formed.Alernatively, the formation of the rims might be syngenetic with the pipe formation.In this concept, low pH water first infiltrates the pipes, becoming saturated with calcium ions as it dissolvesthe walls of the pipes.Then, water moves into the limestone matrix, evaporating at the same time, which leads to supersaturation with respect to the calcium ions, precipitation of calcite and rim formation.Importantly, the numerical models show that the presence of the rims has a hindering effect on the interaction between the pipes, leading to a more uniform length distribution.

Merging of pipes
In coastal areas with good exposures of pipe clusters (Figs. 13 & 14), coalescing pipes (doublets or triplets) are common.Very similar features are observed in the numerical simulations (Fig. 13B).An analysis of the interactions between fingers in unstable growth processes (Budek et al., 2017) shows that coalescence of fingers (pipes) takes place if the permeability contrast between weathered and unweathered rock is not too strong.The shorter finger is then usually attracted to the longer one and eventually merges with it.Note that when pipes merge at some distance from the inlet,coalescence can only be observed after considerable denudation, which is indeed the case in Figs. 13 and 14.
Competition and coalescence are two main mechanisms of pattern coarsening, leading to an increase in the diameters of the largest pipes, as they intercept the flow directed previously to their neighbours.The extent to which one mechanism is dominant over the other again depends on the permeability contrast between the phases (Budek et al., 2017).If it is large, then the pattern is largely hierarchical and the pipes are relatively straight, rarely merging.For smaller permeability ratios, competition plays a lesser role than coalescence and the differences in lengths between the pipes are less pronounced.

Pipes in layered systems
In most sediments the layers have different porosity or chemical composition, and this affects the shapes of solution pipes, which are narrower in some layers and wider in others (Petrus and Szymczak, 2016).As the pipe emerges from a low-porosity layer into a layer with higher porosity, large lateral currents around its tip result in a widening of the pipe, as observed both in natural examples and in numerical simulations (Fig. 15).As the porosity contrast between the layers increases, the competition between the pipes is enhanced, impeding the growth of the shorter ones and enhancing the flow in the longer ones, so that the pipe length distribution becomes more differentiated (Petrus and Szymczak, 2016).Similar effects are evident in dissolutionally enlarged joints in epikarst; Howard (1963) observed undulations on the sides of joints which correlated with the strata, noting that the cavernous widenings have the tendency to form predominantly beneath the more resistant layers.
We note that a study on the effect of differential solubility of the layers on pipe morphology has not been conducted and remains as a future challenge.

Palaeoclimatic implications
Solution pipes have been used as palaeoclimatic indicators (Lipar, 2018), with wetter climates favourable to solution pipe formation whereas more arid climates result in calcrete rim development (Marsico et al., 2003;De Waele et al., 2009;Lipar et al., 2015;Öztürk et al., 2018).In southern Italy the deeper pipes in Sardinia compared to those in Apulia could indicate the exposure of Sardinia to Atlantic cyclones at the time of pipe formation, and the different lengths of the pipes in Apulia could indicate wetter and drier climates (De Waele et al., 2011).In southern Australia, pipes developed during times of higher precipitation at the end of interglacial episodes; the most extensive pipe development occurred during the last interglacial, indicating that this was one of the wettest interglacials, and the wider pipes in Western Australia compared to Victoria may indicate a wetter climate in Western Australia at that time (Lipar et al., 2015;Lipar et al., 2017).Solution pipes developed in regions related to Pleistocene glaciations may be indicators of substantial icemelt (Walsh and Morawiecka-Zacharz, 2001).
Nevertheless, to use the pipes as robust palaeoclimate indicators, the timing of their formation must be well understood.Comparing pipe lengths from two different regions with relatively similar bedrock (e.g, Apulia and Sardinia) is only valid if the pipes formed over the same time interval.The modelling described in Section 4 may give estimates of the time involved based on the width and length of the solution pipes, provided data such as flow rates and the permeability contrast between inside and outside the pipes are known or can be estimated; it may also be important to measure other variables like pCO 2 , saturation of the infiltrating water and temperature.
As stated by De Waele et al. ( 2011), "there is still an unavoidable bias in the pipe measurements due to the surface erosion that often cuts their upper parts".Denudation after pipe formation is clearly evident at many locations by the protruding cemented pipe rims or fillings, and has rarely been taken into account in analyses of pipe morphology, especially their depths.In some cases there has been no denudatione.g., in Poland the rock was covered by glacial till before pipe formation (Walsh and Morawiecka-Zacharz, 2001).
Therefore, before the lengths of pipes can be used as precise palaeoclimatic indicators, denudation rates need to be known, along with the time frame that the pipe was exposed to subaerial conditions (including possible burial by younger sediments).Although the problem has been discussed (Mylroie and Mylroie, 2018), there are few studies that have calculated the individual denudation rates required for particular regions due to local climatic and other variability, e.g.Miklavič (2011) studied calcarenites from Guam, reporting a denudation rate of ~64 mm/ka.

Extension of the modelling to vadose conditions and the impact of climate change
The morphological similarities between phreatic and vadose solution pipes documented in this paper suggest that it might not be decisive whether the conditions are phreatic or vadose.It is possible that these similarities reflect the fact that most of the growth in vadose pipes takes place under episodic large precipitation events, when the pipes are flooded and the conditions are close to phreatic.Alternatively, the differences between the Richards equation (describing the water flow at vadose conditions) and the Darcy law (valid for phreatic conditions) may not be crucial for the morphologies of the pipes.Nevertheless, an indepth assessment of this problem requires an extension of the modelling to account for vadose flow by replacing the Darcy flow solver with that of Richards equation.This would also allow study of wetting front instability in the vadose zone.It may be possible to therefore assess the relative importance of reactive infiltration instability and wetting front instability in the formation of the solution pipes.Both may play a role in piping dynamics, depending on the prevailing conditions.The numerical modelling of vadose and phreatic pipe formation should be supplemented by an in-depth comparison of their morphological properties, such as length distributions, shapes, aspect ratios etc., to address the question as to how and to what extent a different origin of the pipes is reflected in their morphology.
The importance of climate change on solution pipe formation can be assessed by the modelling, in particular the effect of extreme events and changes in the periodicity of the driving forces, in order to address issues such as "will pipe length change if annual rainfall is distributed more evenly or unevenly?"and "what is the impact of extreme rains and flooding events, where the water table rises?"

Conclusions
Solution pipes, vertical or near vertical tubular karst formations within calcareous media with matrix porosity, are a global phenomenon, but their formation has often been understood on a local/regional scale.Through a review of current knowledge of solution pipes, we conclude that the unifying process responsible for their formation is a focused vertical flow of water.Consequently, the nature of focussing of the water flow and its triggers are fundamental for understanding the formation of solution pipes across a diverse range of significantly different climatic settings.In most cases, the focused water flow occurred in the vadose environment, but solution pipes could have also been formed under phreatic conditions beneath the melting glaciers.In either case, the general morphology and distribution of the pipes are strikingly similar, and the numerical model based on phreatic conditions shows promising compatibility with general geomorphic characteristics of solution pipes.
In order to develop solution pipes as a palaeoclimatic indicator, there are still several challenges to consider.Amongst them are precise dating of solution pipe formation and denudation rates, and the development of a numerical model based on vadose conditions.

Declaration of Competing Interest
The authors declare that they have no conflict of interest.

Fig. 4 .
Fig. 4. Analogue experiments showing dissolution channel/solution pipe formation: Athe experiments of of Cha and Santamarina (2016) on the dissolution of quartzitic-carbonate sandy sediment by an acidic solution, Bthe experiments of Kelemen et al. (1995) on dissolution of salt and glass ball mixture, Cthe experiments of Osselin et al. (2016) on the dissolution of gypsum.

Fig. 5 .
Fig. 5. Pipe cross sections: A & C -Plio-Pleistocene calcarenite in Apulia, Italy; B, D, E & F -Miocene calcarenite in Smerdyna, Poland; G -Pleistocene calcarenite at Cape Perron, Perth, Australia.Note that some pipes are rim-free (e.g., D, E) while some have well-developed rims (e.g., C, G); in the case of G, almost filling the pipe, which could be an indicator of pipe senescence.

Fig. 6 .
Fig. 6.Surface exposures of solution pipes around the world: A -Solution pipe with concentric rim in coastal Plio-Pleistocene calcarenite north from Torre dell'Orso in Apulia, Italy; Bsolution pipe with well cemented rim and fill in Pleistocene aeolianites of San Salvador, The Bahamas; C -Solution pipes in Plio-Pleistocene calcarenite on the coast of Marina Serra in Apluia, Italy; Dsolution pipes with thin rims in Pleistocene aeolianites of Victoria, Australia..

Fig. 7 .
Fig. 7. Solution pipes from the cave perspective.A & Bhollow solution pipe formed in the Miocene calcarenites in Poland with the preserved out-washed fill beneath it; Csolution pipe, large enough to represent a vertical pit cave, but still containing cylindrical structure (Tertiary calcarenite, South Australia); Dsolution pipe as a preferential path for root growth (Pleistocene calcarenite, Western Australia).

Fig. 10 .
Fig. 10.Reactive-infiltration instability resulting in a forward erosion piping (left), and backward erosion piping characteristic for seepage erosion in soil (right).The blue arrows mark the direction of pipe growth, whereas the black arrows follow the direction of fluid flow.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. (A) Soil-filled solution pipes in Pleistocene aeolianite beneath a younger aeolian deposit in Canunda NP, South Australia.(B) Solution pipes in the Miocene calcarenite quarry at Smerdyna (Poland).Photo courtesy of Dr. Łukasz Użarowicz (Warsaw University of Life Sciences -SGGW, Poland).(C) Computer simulations of solution pipes formation showing the appearance of hierarchical structure.Successive panels correspond to different moments of time.

Fig. 13 .
Fig. 13.(A) A cluster of solution pipes in calcarenite exposed due to the denudation of the overlying layers -Canunda NP (Australia); the arrow showing knapsack for scale.(B) A cross-section through the group of pipes obtained by numerical simulations.Note several coalesced pipes both in nature (e.g., near the knapsack) and in the simulation.

Fig. 15 .
Fig. 15.(A) Solution pipe cutting through aeolianite interbedded with better cemented and less porous protosol layers at Cape Bridgewater, Victoria, Australia.(B) Solution pipe cutting through aeolianite (upper part) and calcarenite (lower part) layers at Burns Beach, Perth, Australia.(C) Solution pipe in Ramallah, Israel, cutting through stratified dolomites with marl beds.The photograph is courtesy of Sebastian Schmidt of University of Göttingen.(D) Numerical simulations of solution pipe (blue) cutting through the layered rock (Petrus and Szymczak, 2016) (dark red denoted layers of a lower porosity than the rest of the rock matrix).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)