# Fingerprints of surface flows on solid substrates ablated by phase change: from laboratory experiments to planetary landscapes

Sabrina Carpy, Michael Berhanu, Martin Chaigne, Sylvain Courrech Du Pont

► To cite this version:

Sabrina Carpy, Michael Berhanu, Martin Chaigne, Sylvain Courrech Du Pont. Fingerprints of surface flows on solid substrates ablated by phase change: from laboratory experiments to planetary landscapes. Comptes Rendus. Physique, 2025, 25, pp.95 - 142. 10.5802/crphys.230. hal-05133122

HAL Id: hal-05133122

https://hal.science/hal-05133122v1

Submitted on 27 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Physique

Sabrina Carpy, Michael Berhanu, Martin Chaigne and Sylvain Courrech du Pont

Fingerprints of surface flows on solid substrates ablated by phase change: from laboratory experiments to planetary landscapes

Volume 25, Special Issue S3 (2024), p. 95-142

Online since: 5 February 2025

**Part of Special Issue:** Geophysical and astrophysical fluid dynamics in the laboratory

**Guest editors:** Stephan Fauve (Laboratoire de Physique de l'ENS, CNRS, PSL Research University, Sorbonne Université, Université Paris Cité, Paris, France) and Michael Le Bars (CNRS, Aix Marseille Univ, Centrale Marseille, IRPHE, Marseille, France)

https://doi.org/10.5802/crphys.230

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Review article / *Article de synthèse*

Geophysical and astrophysical fluid dynamics in the laboratory /
*Dynamique des fluides géophysiques et astrophysiques au*
*laboratoire*

Fingerprints of surface flows on solid
substrates ablated by phase change: from
laboratory experiments to planetary
landscapes

*De l'échelle du laboratoire aux paysages planétaires, quelles sont les empreintes laissées par des écoulements de surface sur des substrats solides en ablation soumis à des changements de phase?*

Sabrina Carpy<sup>®,*,a</sup>, Michael Berhanu<sup>®,b</sup>, Martin Chaigne<sup>®,b</sup> and Sylvain Courrech du Pont<sup>®,b</sup>

<sup>a</sup> LPG (Laboratoire de Planétologie et Géosciences), Nantes Université, CNRS (UMR 6112), 44322 Nantes, France

<sup>b</sup> MSC (Matière et Systèmes Complexes), Université Paris Cité, CNRS (UMR 7057), 75013 Paris, France

*E-mail:* sabrina.carpy@univ-nantes.fr (S. Carpy)

**Abstract.** Physical or chemical phase changes in ablation, such as sublimation, melting or dissolution, are studied in physics for their many engineering applications. At solid/fluid interfaces, the interaction between a phase change and a flow can lead to pattern formation. In this case, the fluid mechanics associated with such phase changes play a key role in the evolution of terrestrial and planetary landscapes, observed by probes orbiting planets and moons. On Earth, sea ice, glaciers and karst plateaus extend over meters or kilometers. The scale of these landscapes contrasts with the scale of the physical mechanisms that govern their evolutionary dynamics. Indeed, it is the typical size of atmospheric boundary layers or melt-water/vapor/solute films that constrain the heat/concentration transfer at the phase change/dissolution interface, and hence the rate of solid ablation. In many situations, these layers are controlled by fluid flow, either natural or forced convection. In the former case, the flow may be buoyancy driven by the melting/dissolution/sublimation itself, resulting in density stratification caused by, for example, temperature/concentration gradients. This stratification may be stable or unstable. In the second case, the flow forced by winds or slopes can be considered as a flow of an infinite height or of a finite height, such as shallow water flow. In all cases, the mass flux modifies topography, which in turn affects the boundary layer flows and

<footnote>* Corresponding author</footnote>
thus the ablation rate in a retroactive way. In nature, the positive feedback between geometry and mass transfer drives the spontaneous formation of characteristic patterns at different scales. These patterns are not just geological curiosities, such as Zen stones or dirt cones but markers of the hydrodynamic processes at work. Many landscapes are shaped by regular, repeated patterns, whether sharp-edged, scalloped, parallel-crested, or stepped. By experimentally investigating different modes of flow transport on solid substrates undergoing physical or chemical phase change, this review aims to highlight the role of the flow transport mode in the diversity of patterns observed on analogous materials. Understanding the diversity of these patterns is key to assessing the environmental conditions under which they form on planets such as Mars or Pluto, where phase changes play a very important geomorphological role.

**Résumé.** Les processus en ablation, tels que les changements de phase, physiques (sublimation/fusion) ou chimiques (dissolution), sont étudiés dans de nombreuses applications techniques en physique. Pour les interfaces solide/fluide, l'interaction entre un changement de phase et un écoulement peut conduire à la formation de motifs topographiques. Dans ce cas, la mécanique des fluides associée à de tels changements de phase joue un rôle primordial dans l'évolution des paysages terrestres et planétaires, observés par les sondes en orbite autour des planètes et des lunes. Sur Terre, la glace de mer, les glaciers et les plateaux karstiques s'étendent sur des mètres ou des kilomètres. L'échelle de ces paysages contraste avec l'échelle des mécanismes physiques qui régissent leur dynamique évolutives. C'est en effet la taille typique des couches limites atmosphériques ou des films d'eau de fonte/soluté qui limite le transfert de chaleur/concentration à l'interface de changement de phase/dissolution, et donc le taux d'ablation du solide. Dans de nombreuses situations, ces couches sont contrôlées par l'écoulement du fluide, de type convection naturelle ou forcée. Dans le premier cas, la flottabilité de la fonte/dissolution elle-même, entraîne une stratification en densité causée par des gradients de température/concentration, qui peut être stable ou instable. Dans le second cas, l'écoulement forcé par des vents ou des pentes peut être considéré comme un écoulement de hauteur infinie ou finie, tels que les écoulements à surface libre. Dans tous les cas, le flux de masse modifie la topographie, qui en retour impacte les écoulements de couches limites (fines ou épaisses) et donc le taux d'ablation de manière rétroactive. Dans la nature, la rétroaction positive entre la géométrie et le transfert de masse entraîne la formation spontanée de formes caractéristiques à différentes échelles. Ces formes ne sont pas seulement des curiosités géologiques, comme le sont par exemple les « Zen stones », ces pierres sur un piédestal de glace ou encore les cônes de poussière, mais des marqueurs des processus hydrodynamiques en jeu. De nombreux paysages sont façonnés selon des motifs réguliers et répétés, qu'ils soient pointus, en forme de cupules, à ligne de crêtes parallèles entre elles ou encore en marches d'escaliers. Dans ce papier de synthèse, nous recensons les études expérimentales dédiées à l'étude des motifs générés par différents modes de transport des flux sur des substrats solides soumis à des changements de phase physiques ou chimiques. Nous souhaitons ainsi mettre en avant le rôle du mode de transport des flux dans la diversité des motifs observés sur des matériaux analogues. Comprendre la diversité de ces motifs est essentiel pour évaluer les conditions environnementales dans lesquelles ils se forment, notamment sur des planètes comme Mars ou Pluton sur lesquelles les changements de phase jouent un rôle géomorphologique majeur.

**Keywords.** Sublimation, Melting, Dissolution, Natural solid bedforms, Experimental analogues, Pattern formation, Classification.

**Mots-clés.** Sublimation, Fusion, Dissolution, Motifs solides naturels, Analogues expérimentaux, Formation de motifs, Classification.

**Funding.** French Research Council (SHERPAS ANR-23-CE49-0006 and PhysErosion ANR-22-CE30-0017).

Manuscript received 9 July 2024, revised 18 December 2024, accepted 20 December 2024.

# 1. Ablation by mass transfer on solid substrate submitted to phase change

In engineering applications, mass transfer techniques such as dissolution or sublimation have been used to constrain fluid/structure interactions. Naphthalene sublimation is a process that has been widely used in aerodynamics in a broad range of topics such as pollutant dispersion in urban areas [1] or ablation processes in space capsules [2]. Ablation patterns may appear during the descent phase of an atmospheric space capsules [3,4], similar to meteorites' regmaglypts [5], due to the laminar-turbulent transition of the flow. Other solid substrates such as bitumen, copper, brass, aluminum or magnetite, when subjected to liquid water flow, develop ablation patterns [6], which are similar to those appearing on molten ice subjected to water flow or
sublimated ice subjected to wind flow. Whether the phase change is a solid–gas phase transition (sublimation), a solid–liquid phase transition (melting) or a dissolution phenomenon, substrate ablation is strongly linked to a feedback between the fluid flow and the underlying topography, which generally leads to the emergence of patterns. Hereafter, we will refer as solid bedforms to these periodic patterns that are induced by mass transfer built up at the interface between a moving fluid and a solid surface. Although previous experimental works have observed the emergence of solid bedforms, their morphology and dynamics were not the main focus of these engineering applications.

Interest in solid bedforms has been renewed by space missions exploration that have revealed many surprising landscapes, such as the spiral troughs on the north polar cap of Mars, interpreted as cyclic steps [7] and the bladed terrain of Pluto, interpreted as *penitentes* [8]. Sublimation of ice combined with surface flows like wind could be a major geomorphic agent in explaining these planetary landscapes [9]. On Earth water ice on its solid phase represents 2.4% of the global volume of water, but covers 14% of the Earth' surface [10]. Water ice on Earth is subjected to melting and sublimation processes at its surface. However, for most planetary bodies, excluding Earth and Titan, environmental conditions (pressure, temperature) at the surface do not allow the existence of three phases of their exotic ice. N₂, CH₄, CO₂, NH₃ices [11], are subjected to surface sublimation that prevails in the solar system. Soluble rocks on Earth also present solid bedforms such as scallops or *rillenkarrens* that could be found on limestone, gypsum and salt [12,13]. Solid bedforms therefore have a wide variety of compositions, environments, shapes and physical phenomena involved.

The apparent regularity of these solid bedforms is surprising, given the complexity of the mechanisms behind their formation [14]. Lab-scale experiments contribute significantly to our understanding of their formation and dynamic processes. It requires constant control of experimental parameters (pressure, temperature, partial pressure/humidity, acidity) to reproduce natural cases. Often, analogues substrates are selected to accelerate process kinetics to experimentally more feasible times. So, experimental approaches that simulate rock dissolution with the aim of creating *rillenkarren* [15–17], scallops and *flutes*, rely on the use of plaster [18–20] or even caramel [21], for example. In many situations, fluid dynamic instabilities [14] (e.g. the *Rayleigh Plateau, Kelvin Helmholtz, Rayleigh Taylor* or *Rayleigh Bénard* instabilities) explain the emergence of natural patterns with a characteristic scale. For example, the growth of transverse patterns (perpendicular to the mean flow), like the dissolution waves on limestone walls in the laminar–turbulent regime [22], or the ripples and cyclic steps on ice due to run-off flows [23,24], have been explained as a positive feedback between the topography, the thickness of the turbulent boundary layer and the ablation rate. Moreover, longitudinal patterns (parallel to the mean flow), like the *Rillenkarren* formed by run-off flows on limestone, could be explained by a linear stability study of the coupled system of equations involving mass transfer and fluid motion, that affect the topography of the bed [17,25,26] under specific boundary conditions.

What environmental conditions are required for the formation of these very large natural solid bedforms? What physical parameters control their size, shape and spatial distribution? How can we explain the fact that the same patterns are found in rock and ice on Earth, but also on other icy planetary surfaces, when the phase changes at work are not always the same? In order to understand the mechanisms responsible for the formation of solid bedforms subjected to a surface flow, and to deduce the links between their shape, dimensions and the characteristics of the environment in which they develop, we describe some experimental results and if any, the scaling laws that describe the development of these solid bedforms during the dissolution [17,20,22,27–29], the melting [30,31] and the sublimation [32] processes. Since the time scales involved for the formation of natural bedforms are too large, experimental models are used to test formation hypotheses over realistic experimental timescales, in order to reconstruct
the formation history of these patterns and the environmental conditions under which they emerged.

The morphological similarities between some solid bedforms resulting from different ablation processes [6,33] in the case of ice sublimation, ice melting and rock dissolution will be described in Section 2 for: sharp patterns; polygonal depression; crested bedforms parallel or perpendicular to the flow; stepped patterns; and debris-covered ice structures. In Section 3, we explore the positive feedback that influence the evolution of the solid/fluid interface by identifying the relevant equations of the problem (hydrodynamic-morphodynamic coupling) with a jump condition at the interface. We will only explore cases where one mechanism is dominant over the other (e.g. double diffusion is out of the scope of this paper). Analogous modeling presented in Section 4 can then be used to constrain the role of fluid mechanics in the evolution of these forms and test hypothesis. The diversity of these patterns is reinforced by the complexity of the flows that carry a passive scalar such as temperature or concentration: a change in density in a natural convection flow would give rise to buoyancy instabilities; a disturbance to the topography in an advective flow can lead to space modulations of concentration. On a large scale, these mechanisms are responsible for the formation of surprising landscapes when hydrodynamic surface flows drain through karst environments, or when katabatic winds descend the slopes of terrestrial and extraterrestrial glaciers. The discussion in Section 5 focuses on the genetic link between the type of flow and the bedforms that could be helpful in the classification process. We propose to discuss the effect of non-linearity on the resulting bedforms and the perspective to use some scaling laws in a planetological context.

## 2. Diversity of natural solid bedforms by phase change on solid substrate

There are a wide variety of natural bedforms carved out of ice or soluble rock. These forms are fascinating in their uniqueness or repetitiveness. Examples include stones that look as if they have been placed on a pedestal of ice, blade-like shapes made of ice or soluble rock, steps of ice or soluble rock in streams of water, gullies that incise the soluble rock with spiky crests, or gentle undulations on expanses of blue ice, which can take the form of scallops on cave walls. All these forms are subject to phase changes, whether liquid/solid or solid/vapor for ice, or dissolution for soluble rock. This section highlights the similarities in shape despite the differences in mass transfer mechanisms, and lists the natural forms that have been observed. Experimental studies have been carried out to explain some of these shapes, which we review in Section 4.

### 2.1. *Sharp patterns*

First described by Charles Darwin in Chile in 1835, *penitentes* are shapes formed by sublimation in old, compact ice/snow, either as spikes or as blades pointing skywards (Figure 1(a)). On Earth, *penitentes* are mainly located in high altitudes: they have been particularly well described in the Santiago Highlands, where all snow and icy surfaces above 4000 m are covered with them; on the summit and flanks of the Tutupaca volcano in Peru [34]; on the Quelccaya Ice Cap in Peru [35]; and on the Tapado Glacier in northern Chile [35]. But similar forms have also been observed at very low altitudes such as Greenland [36]. Given the climatic conditions that favour their formation, *penitentes* are generally found in regions where humidity and temperature are low, on downslopes where the wind is very weak but sunshine is intense. As a result, similar patterns on the icy surfaces of other Solar System bodies have been interpreted as *penitentes*, such as the Bladed Terrain Deposit of Pluto's *Tartarus Dorsa* region [8,37]. It has been suggested that *penitentes* as high as 15 m may exist on Europa's surface [38], although the resolution of currently available images of Europa's surface is insufficient to verify this, but this hypothesis is debated
<figure><img src="image_3.png" /><figcaption>**Figure 1.** *Penitentes* and *Tsingy*. (a) *Penitentes* (Rio Blanco, Central Andes of Argentina). Crédits: P. Dubuc, Arvaki, CC BY-SA 3.0 license. (b) *Tsingy* of Bemaraha [47,48]. Credits: A. S. Maloney, CC BY-SA 3.0 license.</figcaption></figure>

because of the low atmospheric pressure on Europa [39,40]. From base to tip, these spikes can reach heights of a few centimetres [39] to several metres [41], depending on environmental conditions. The cross-section of the *penitentes* is rather lenticular, tapering towards the pointed summit, and their blades are parallel to the sunlight at midday. The horizontal spacing is of the same order of magnitude as the height of the tips [42]. The orientation of the *penitentes* in terms of azimuth is linked to the local geometry of the insulation [43]. In particular, Matthes et al. [44] observes an inclination of the *penitentes* following the elevation of the midday sun. For geometric reasons, sublimation is more effective in hollows. Illumination plays a key role [45], as do heat transport in the ice and the transport of vapor from the ice [46]. The bedforms detected on the surface of icy substrates, even if their formation is induced by sublimation, can be interpreted in different ways. Knowledge of the mechanism of formation of the different sublimation induced shapes can help to classify solid bedforms between *penitentes* and other bedforms, such as scallops (Figure 2(a)) or transverse linear waves (Figure 3(a)).

Pointed shapes can also be found in the dissolution of rocks, in karst. In that case, the limestone substrate is rigid and not very porous. One of the most intriguing sharp-edged landscapes formed by rock dissolution are limestone forests, found at tropical locations such as in the western part of Madagascar in the area called *Tsingy* [47,52], in Malaysia [53] or in south China [54]. Spectacular sharp and vertical pinnacles, tens of meters high, rise from the surface of a horizontal *Karst* plateau (Figure 1(b)). The general shape is needle-shaped, with a slope of 15° to 20°, which remains even when the material has been strongly dissolved. These sharp-edged *karren* forms, which develop in groups, have several degrees of evolution: from dragon teeth to sharp-edged pinnacles before evolving to *Tsingy* [55]. The smallest in size belong to the *lappies'* family (2–3 m) and the largest, by one order of magnitude, belong to the *mega-lappies'* family (20–30 m). On their slopes, centimeter scale *rillenkarren* (see Section 2.4) can be observed. Because the material is pure and not porous, dissolution remains on the surface, not in the bulk. Formation of the limestone forests appears to be complex, but is believed to be caused by a dissolution process under the run-off caused by tropical heavy rain. The rainfalls channelized by large scale fractures are collected into an underground cave network. Subjected to run-off flow, soluble rock are less dissolved as it flows along the slope, such as the layer gain in concentration. This could preserve the homothetic character of these forms (recession of these faces parallel to itself). If the material is porous, such as dolomite, then the shape is not sharp-edged, but resembles a mushroom. In very seldom locations, in very arid area like the Atacama desert in Chile or the banks of the Dead Sea in Israel, salt pinnacles of few ten centimeters can be observed on salt deposit [56].
<figure><img src="image_4.png" /></figure>

**Figure 2.** Scallops on different substrates. (a) Icy scallops formed by sublimation (French Alps). Credits: A. Mounier. (b) Icy scallops formed by melting on an overturned iceberg (Paulet Island, Antarctique). Credits: P. Colla. (c) Limestone scallops formed by dissolution (Grotte de Saint-Marcel d'Ardèche, France). Credits: M. Bordiec/S. Carpy/O. Bourgeois, reprinted with permission from [49]. (d) Limestone scallops by dissolution (Moulin de Vogüé, France). Credits: M. Bordiec/S. Carpy/O. Bourgeois, reprinted with permission from [49].

## 2.2. Polygonal depression

Suncups and ablation hollows are bowl-shaped, circular or oval in plan, and parabolic or tapered in vertical section, and occur on ice and snow. In comparison to *penitentes* on Earth, they are located at lower altitudes, where the ambient temperature is generally higher (e.g. in the Alps [45]), and in areas more exposed to winds than those where *penitentes* develop. With amplitudes ranging from a few centimetres to half a metre, they do not reach heights comparable to those of the *penitentes*, although they are also controlled by solar radiation and influenced by albedo and the amount of dust in the ice/snow [45,57]; the role of wind is not clearly excluded. Suncups could constitute a surface texture found particularly on the North Polar Cap of Mars [9,58–61] and on the Sputnik plain of Pluto [62]. Similar shapes appear on the walls of ice caves [63–65] and are known as “scallops”. Ice caves, although they represent only a tiny fraction of the Earth’s cryosphere, exist on every continent outside the tropics and are defined as rock caves containing ice, in seasonal or perennial form [12]. Conceptual models of the climate prevailing in these ice caves have been built from measurement campaigns conducted over several decades [66], in particular in the *Eisriesenwelt* cave (Salzburg, Austria) [65], which is one of the largest known ice caves in the world, with its more than 40 km of galleries. In winter, the temperature of this cave is below 0 °C and the air is dry [32], which favors episodes of sublimation, enhanced by cold, dry air that is blown and circulates through the cave, at speeds not exceeding 0.2 m·s⁻¹ [63,65,67]. In most ice caves, the temperature of the ice surface remains relatively stable throughout the year, with varying rates of sublimation depending on the geometry of the cave. For example, although ablation is restricted to the winter months when sublimation is most
<figure><img src="image_5.png" /></figure>

**Figure 3.** Ripple-like bedforms over different substrates submitted to mass transfer: (a) by sublimation (blue ice ripples in *Svea*, Antarctica adapted from [50]. Copyright 1999 by the American Geophysical Union). (b) By melting, Rivière St-Croix ice ripples, Wisconsin, USA [51]. (c) By sublimation, martian North polar Cap, reprinted with permission from [49]. *HiRISE* PSP_009689_2645. (d) By dissolution, *Coups de Gouge*, Grotte de Saint-Marcel d'Ardèche, France. Credits: M. Bordiec/S. Carpy/O. Bourgeois, reprinted with permission from [49].

effective, the annual ablation rate in ice caves ranges from 3 mm·year⁻¹ [68] to 35 mm·year⁻¹ [64], resulting in sublimation rates of between 3.8 × 10⁻¹⁰ m·s⁻¹ and 4.5 × 10⁻⁹ m·s⁻¹.

Scallops are widespread polygonal depression patterns, both on ice [63,64] (Figure 2(a–b)) and on soluble rock walls [12,55] (Figure 2(c–d)), and on many other substrates [6]. By analogy, Meakin and Jamtveit [69] suggest that these shapes can be likened to those observed on certain meteorites [5], on corroding metal pipe surfaces [70] and cohesive surfaces in fluvial systems [71]. The characteristic sizes of scallops range from centimetres to metres. They are shell-shaped, ablation-induced depressions. When they cover a surface, they exhibit uniformity in terms of size, shape and spacing (Figure 2). They are formed when solid surfaces are immersed in a fluid different from their composition when subjected to mass transfer induced ablation. The scallops are asymmetrical compared to suncups, which are circular depressions. They are oriented with respect to the flow: the downstream side forms a sharp shape and the upstream side is steeper in longitudinal section (Figure 2(d)). Similarly, scallops have sometimes been assimilated to shapes transverse to the flow direction [6,22], although they differ from ripples and flutes, which are transverse waves (Figure 3), which are considered to be special cases of scallops [20] (see Section 2.3). Due to the complexity of flows and possible variations in the conditions that lead to the formation of these solid bedforms, the understanding of their development is not completely acquired [69].
## 2.3. *Crested patterns perpendicular to the flow direction*

Transverse linear waves (Figure 3) form and develop on the surface of icy substrates (by sublimation or melting) or soluble rocks (by dissolution), and they are known as “sublimation waves”, “ices ripples”, “fluted scallops”, “flutes” or “solution ripples” in karstology. Their crests are regularly spaced and parallel, oriented perpendicular to the flow [12,55]. The waves are likely to migrate in the direction of flow over time. Their cross-sectional profile in relation to the direction of flow is slightly asymmetrical (Figure 3(b)), with a downstream slope greater than the upstream slope [20,72]. They exist on all types of solid substrates, such as salt [73] and ice, immersed in liquid water [23,31], in ice caves [67] or on Antarctic ice surfaces exposed to strong winds [74]. If dissolution forms are naturally very common in limestone caves [19,67,71] with a range of wavelengths from a few centimetres to several metres (Figure 3(d)), measurements of the flow characteristics responsible for their formation are rare in these natural environments. Sublimation waves (Figure 3(a)) are periodic, linear and transverse bedforms of different scales (from centimetres to metres) depending on their dry windy environment [74,75]. On Earth, Blue Ice Areas that cover 1% of Antarctica reveal the presence of sublimation ice ripples [4], about 10 cm in wavelength that emerge in ~2 month and migrate at 2 cm/month. The difficulties of field campaigns in such remote places with uncontrolled environmental conditions explain the very few field measurements collected on sublimation waves. The compilation of documented natural examples of sublimation waves has revealed a real lack of morphological and kinetic data related to the development of these solid bedforms [32]. Shapes on the North Polar Cap of Mars (Figure 3(c)) and certain shapes on the surface of Pluto (plain *Sputnik Planitia*) can also be interpreted as sublimation waves [11].

## 2.4. *Crested patterns parallel to the flow direction*

The term *Karren* designates the dissolution patterns, which are created by the run-off flows on denuded rock walls [12,55,76]. We focus here on the patterns, which are caused by the flows of water driven by gravity along the bare slope and originated from rain and snow falls. Some effects of seepage are possible for slightly porous and fractured rocks and can affect the morphologies. The *Karren* patterns occur mainly by dissolution of limestone and dolomite, but they can also appear on gypsum and salt (Halite) [56]. In view of the wide variety of shapes and front, a specific terminology has been defined to classify the different types [76]. Among the most important examples, *Rillenkarren* are adjacent grooves or small channels directed along the main slope [76,77] (Figure 4(a–d)). In cross-section, their profiles are parabolic, separated by sharp ridge lines. Viewed from above, these lines are remarkably parallel and spatially regular, typically 1.5 to 2.5 cm wide and 2 to 6 mm deep in limestone. However, the typical size can be significantly larger than the centimeter, like for the *Rillenkarren* found on Madre de Dios Island in Chile (Figure 4(c)). In the field, the flow creating these patterns is a thin film of water that runs down the slope. The run-off flow is fed by the impact of droplets and the role of the droplet impact remains an open question [15]. Moreover, *Rillenkarren* grooves are preferentially observed on the sloped wall (typically 45°) above a smooth area, where the water film is thicker [15]. The *Wandkarren* are a special case of the *Rillenkarren* and occurs on nearly vertical walls created by large rock fractures. The *Rinnenkarren* or Runnels are dissolution channels of larger size, which are formed once the water has been channelized in a rivulet. *Rundkarren* are similar to *Rinnenkarren*, but the ridges between the channels are rounded by subsoil dissolution [76] or maybe by snow cover. *Meanderkarren* are *Rinnenkarren* with meanders. The flow dynamics of water channelized in *Rinnenkarren* systems have been studied on the field and in laboratory reproduction with non soluble materials [47].
<figure><img src="image_6.png" /></figure>

**Figure 4.** (a) Nearly vertical *Rillenkarren* on limestone. The photo was taken in Velebit (Croatia). The typical groove width is about 1 cm. Credits S. Courrech du Pont. (b) *Rillenkarren* of the order of one meter on Madre de Dios Island in Chile, Ultima Patagonia expedition, *Centre Terre*. Note the human at the bottom of the picture, to get an idea of the characteristic scale. From Tourte et al. [78]. (c) Inclined limestone pavement or *Lapiaz* at the site *Combe Oursière* near Lans en Vercors, France. Rounded *Rillenkarren* are directed along the main slope of about 45° inclination. Credits M. Berhanu, M. Chaigne and J. Derr. (d) 3D reconstruction of *Rillenkarren* on the previous field site at small scale using a structured light scanner. The slope is aligned along the *x* axis. Sharp ridges can be seen.

## 2.5. Stepped bedforms

Cyclic steps are bedforms that are shaped like a succession of steps, moving against the flow [79]. These steps are associated with hydraulic jumps at the free surface. Cyclic steps occur at the bottom of *bedières* (Figure 5(a)). They have been studied experimentally on sloping ice over which liquid water flows [24]. Cyclic steps have also been described on solid rock substrates over which rivers flow [80]. In the case of sublimation, the spirals of the North Cap of Mars have been interpreted as cyclic steps of several kilometers in wavelength, influenced by jumps of katabatic origin [7].

There are also stepped, regular undulations known as *trittkarrens*, *stepkarrens* or *heelprint karrens* (Figure 5(b)), which are observed in dissolving rocks with a gentle slope. They are perpendicular to the sharp ridge lines of the *rillenkarrens* described in Section 2.4 and therefore
<figure><img src="image_7.png" /></figure>

**Figure 5.** (a) Cyclic steps by melting in *Bédières*. Crédit: L. Karlstrom, adapted from [80] with permission of John Wiley and Sons. (b) *Trittkarrens* steps perpendicular to the crests of *Rillenkarrens*. Crédit: L. Perrella, CC BY-SA 3.0 license.

perpendicular to the flow. They are small-scale *karrens* (on the order of centimetres) that can form on various dissolving rocks, such as limestone and gypsum [81]. Morphologically, they are described by Ford & Williams [12] as arched-front forms with a flat floor, 10 to 30 cm in diameter, which are open at the end of the downslope. Mottershead [81] characterizes them more by their "stepped shape" in cross-section and the presence of a slight gradient.

## 2.6. Debris-covered ice structures

Natural pattern can also form because an object changes boundary conditions. Rock or debris covering the ice can insulate the underlying ice or increase the melting rate. The resulting differential melting forms characteristic structures observed on glaciers such as glacier tables, dirt cones or ice sails as shown in Figure 6 [82–86]. Glacier tables are large rocks standing on an ice foot (Figure 6 top), dirt cones are ice cones covered by a layer of debris (Figure 6 bottom left), and ice sails (or ice pyramids) are clean ice structures protruding from the debris-covered surface (Figure 6 bottom right).

Many glaciers are covered by a layer of debris, which modifies the melting rate of the ice. Several field measurements have shown that the melting rates initially increases with the thickness of debris cover, reaches a maximum for a thickness of about 2 cm, and then decreases with debris thickness [87–89]. Debris insulate the underlying ice, which explains the decrease in melting rate with the thickness of debris. The initial increase with thickness is ascribed to the patchiness of the debris coverage, which gradually changes the albedo of the surface as the coverage proportion increases and to a wind-shielding effect of the cover; the wind velocity at the ice/debris interface rapidly decreases with debris thickness [89]. This decrease in wind velocity reduces the rate of evaporation of the melt water and the consumption of energy in latent heat (for evaporation), which increases the energy available for melting [89].
<figure><img src="image_8.png" /></figure>

**Figure 6.** Melting and sublimation structures due to debris or rock cover. Top: Glacier table (rock supported by a pillar of ice) forming on the *Mer de Glace* glacier in the French Alps from Hénot et al. [85]. Inset: Zen stone (stone supported by a pillar of ice) on Lake Baikal, southern Siberia. Credits: Olga Zima. Glacier tables and zen stones form at different scales and by slightly different mechanisms (Section 4.3), figure adapted under the CC BY 4.0 license. Bottom left: Dirt cones (ice cones covered by a layer of grains) on Sólheimajökull Glacier, Iceland. Credits: Patagonier, CC BY-SA 3.0 license. Bottom right: Ice Sails on the Baltoro Glacier (Karakoram mountain range, Pakistan), from Evatt et al., 2017 [82]. Credits: C. Mayer, 2011, under the CC BY 4.0 license.

**3. Modelization of phase change dynamics in presence of a flow**

The natural bedforms reviewed in Section 2 emerge when there are heterogeneities in the ablation velocity of a bed, which can result from a coupling between the rate of phase change (dissolution, sublimation or melting) and a flow. Here we introduce the equations that allow to model these processes occurring at the solid/fluid interface. Solidification, condensation and precipitation are not discussed in this review, but in some natural examples deposition and ablation effects may be intertwined [69]. The set of equations which will be found for melting, sublimation or dissolution will also be valid for solidification, condensation or precipitation.

During a phase change of a solid, the chemical bonds between the particles constituting the solid are broken to form a liquid phase or gaseous phase depending on the thermodynamic
conditions, usually the temperature and the pressure. A phase is a region of material in which the chemical and physical properties are locally homogeneous at the molecular scale. Thermodynamics determines the equilibrium state, with the occurrence of one or another specific phase, or even the coexistence of several phases. The dynamics of the phase transition, the rate of phase change from one phase to another, is an out of thermodynamic equilibrium phenomenon.

In some cases, the dynamics is set by kinetic effects at the molecular scale associated to chemical bonds breaking and formation. However, when these processes are "fast", the limiting mechanism that determines the dynamics of the phase change is a transport mechanism that either transfers the energy necessary for the phase change, or removes chemical species to maintain the interface in a non-equilibrium state, or even a combination of both phenomena. For example, the melting of ice in salty (NaCl) water could be limited either by the heat transport or by the salt transport to the ice through the melt water, because the melting temperature depends on the salt concentration. In the absence of flow, these transfers occur by thermal diffusion or by chemical species diffusion. In the presence of a flow in the fluid phase, they are generally greatly enhanced by advection phenomena. The modeling of phase change dynamics can be thus often described as an advection-diffusion problem in the fluid phase, where the rate of phase change is expressed as a boundary condition.

## 3.1. Boundary conditions at an evolving interface

We consider the phase transition between a solid and a fluid phase, liquid or gas. For simplicity, we focus on the melting of a pure substance, the dissolution of a homogeneous solid in a fluid phase composed of a solvent and a solute and the sublimation of a solid in a gaseous phase containing a small proportion of the vapor generated by the phase change. In the second case, the solute molecules are produced by the dissolution of the solid and in the third cases, the vapor molecules are produced by the sublimation of the solid. These two cases are formally equivalent. Dissolution patterns can therefore have analogous counterparts for sublimation.

### 3.1.1. Kinetics of the phase change

At the solid/fluid interface, the rate of detachment of the molecules and atoms constituting the solid is modeled by the chemical kinetics of phase change. In a first approximation, this rate can be taken to be proportional to the distance to thermodynamic equilibrium [90,91].

**Melting.** For the melting/solidification, except for the rapid phase change of metastable materials such as supercooled water, the kinetic molecular effects are usually negligible [92]. In this case, the solid/interface can be considered locally at the thermodynamic equilibrium. This means that the pressure and the temperature at the interface belong to the coexistence curve for two phases at thermodynamic equilibrium. If the pressure is imposed externally, then the temperature is the corresponding melting temperature $T_m$, which can be determined by the Clapeyron relation [93]. For example, when ice melts at atmospheric pressure on Earth, the melting temperature during the phase change is equal to 0 °C. The melting temperature also decreases with the local curvature of the interface due to Gibbs-Thomson effect [90], as the local environment differs when a molecule on the solid is located at a tip. This effect, which is important for crystal growth, is often negligible for melting at macroscopic scales. Under these conditions, the melting rate is determined by the thermal energy input (see Section 3.1.3).

**Dissolution.** In contrast, in the case of dissolution, kinetic effects are often considered relevant [91]. We consider the case of a pure solid in contact with a liquid phase constituted of a solvent and a solute, whose molecules or ions are the same as those in the solid phase. The
proportion of solute particles in the solution is quantified by the mass solute concentration field $c(\mathbf{x}, t)$, i.e. the local mass of solute per unit volume. This case corresponds for example to the dissolution of solid salt (NaCl, halite crystal), which decomposes into the ions sodium $Na^+$ and chloride $Cl^-$ in water. Due to electroneutrality, the concentrations of positive and negative ions are equal at the macroscopic scale and the solute concentration is numerically equal to the sum of the concentrations of positive and negative ions. For example the mass solute concentration of NaCl: the mass of solute per unit volume is equal to the mass concentration of ions $Na^+$ plus the mass concentration of ions $Cl^-$ . The saturation concentration $c_{sat}$ is the solute concentration when the solid and the liquid phases are at thermodynamic equilibrium. The attachment and detachment rates of ions are then equal. As a thermodynamic parameter, $c_{sat}$ depends on the temperature, the pressure, the presence of other chemical species and also on the local curvature by a phenomenon analogous to the Gibbs–Thomson effect [94]. This last effect is relevant for scales below the micron and is neglected in the remainder. For small concentrations, the chemical activities can be assimilated to concentrations.

**Sublimation.** In the case of sublimation of a solid (e.g. water ice), the concentration of the volatile substance (e.g. water vapor), is usually expressed in terms of the partial vapor pressure. According to the ideal gas law, which holds for diluted gases, the partial pressure is proportional to the mass concentration, i.e. $p_x = c_x RT/M_x$, where $M_x$ is the molar mass of the specie x and R is the universal gas constant. Similarly, when the interface between the atmosphere and the sublimating solid is at thermodynamic equilibrium, the partial pressure is equal to the saturation vapor pressure $p_{sat}$. This quantity can vary strongly as a function of the temperature according to the Clapeyron relation [93].

### 3.1.2. Relationship between phase change rate and kinetic effects

Dissolution and sublimation rates can be related to the rate of the chemical reaction $\mathcal{R}$, which corresponds to the breaking of the chemical bonds in the solid phase. We define $\mathcal{R}$ as the solute mass flux for dissolution or the vapor mass flux for sublimation. $\mathcal{R}$ is homogeneous to a mass per time and surface unit. The simplest models for dissolution/sublimation processes assume that $\mathcal{R}$ is proportional to the distance to equilibrium expressed in terms of concentrations, $c_i - c_{sat}$, where $c_i$ is the solute or vapor concentration at the solid interface. We call $\alpha$, the proportionality coefficient, which expresses the speed of the chemical reaction. This linear relationship can be justified by more rigorous considerations of chemical kinetics [95] or gas kinetics [96–98].

If $c_i < c_{sat}$, the solid recedes by dissolution (sublimation) and vice versa if $c_i > c_{sat}$, the solid grows by precipitation (condensation). Writing the conservation of solute mass, the mass flux $\mathcal{R}$ is also equal to the product of the density (mass per unit volume) and the velocity of the interface. The following equation then relates the velocity of the interface, i.e. the rate of phase change, to the speed of chemical reaction in the laboratory reference frame for a motionless solid phase, at a point of the interface:

$$ \mathcal{R}(\mathbf{x}, t) = \alpha(c_{sat} - c_i) = -\rho_s \mathbf{v}_i(\mathbf{x}, t) \cdot \mathbf{n}, \quad (1) $$

where $\rho_s$ is the solid density and $\mathbf{v}_i$ is the interface velocity (or dissolution/sublimation velocity) and $\mathbf{n}$ is the vector normal to the interface pointing outwards the solid surface.

To summarize, the mass flux for the two ablation processes, sublimation and dissolution, actually share a common formulation, in the sense that the existence of the transfer and the direction in which it takes place are conditioned by a concentration, either in the atmosphere or in the solution. For dissolution and sublimation, $\alpha$ is a positive coefficient expressing the speed of the chemical reaction and is homogeneous to a velocity. $\alpha$ qualitatively translates the finite time required to break and create chemical bonds in the dissolution and sublimation processes. The
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Scheme for determining the jump conditions at the evolving interface $S(t)$ separating two phases, for example a solid $\mathbb{R}^s$ and a liquid $\mathbb{R}^l$. We consider an infinitesimal spatial domain $\Omega(t)$ on either side of the interface $S(t)$. $v_i$ is the velocity field of displacement of the interface. The outer boundary of the domain is denoted by $\partial\Omega(t)$, which is composed of $\partial\Omega^s$ and $\partial\Omega^l$ for the part in contact with $\mathbb{R}^s$ and $\mathbb{R}^l$, respectively.</figcaption></figure>

interface velocity can vary along the solid surface and its heterogeneity causes pattern formation. For sublimation, using the theory of gas kinetics, the coefficient $\alpha$ for vacuum or very tenuous atmospheres is the product of a thermal velocity of the gas molecules times a dimensionless adsorption coefficient $\varepsilon$ on the surface [96–98]:

$$ \alpha = \varepsilon \sqrt{\frac{RT_s}{2\pi M_s}}. \qquad (2) $$

$\alpha$ is thus a function of the surface temperature of the solid $T_s$ and depends on the molar mass of sublimated molecules $M_s$.

### 3.1.3. Moving interface

In Equation (1), the concentration at the interface, $c_i$, must be determined to solve the dissolution or sublimation problems. Therefore, to determine the phase change rate $v_i$, in addition to these kinetic considerations, it is necessary to take into account the balance equations through the continuity equation, the solute (or vapor) mass equation and the energy equation, in the presence of a moving interface. In melting problems, kinetic effects are usually negligible, so the interface temperature is the melting temperature $T_m$ given by the thermodynamics, and the melting rates are in most cases determined by the balance equations.

The Kotchine's theorem or jump relations [92,99] provides an elegant method to derive the conservation relations in three dimensions. We consider a quantity $f(\mathbf{x}, t)$ obeying to a local conservation equation with no source or sink term

$$ \frac{\partial f(\mathbf{x}, t)}{\partial t} + \nabla \cdot (\mathbf{J}_f(\mathbf{x}, t)) = 0, \qquad (3) $$

where $\mathbf{J}_f(\mathbf{x}, t)$ is the incoming flux of the scalar quantity $f$, which can be mass, concentration or temperature and $\nabla \cdot (\mathbf{s})$ is the divergence operator. By applying the Reynolds transport theorem for the quantity $f(\mathbf{x}, t)$ to the infinitesimal domain $\Omega(t)$ shown in Figure 7, where the two phases considered are the solid (s) and the liquid (l), the Kotchine's theorem gives:

$$ (f^l - f^s)(v_i \cdot \mathbf{n}) = (\mathbf{J}^l - \mathbf{J}^s) \cdot \mathbf{n}. \qquad (4) $$

This relation gives the boundary condition that relates the interface velocity $v_i$ to the quantities on each side of the interface: $f^l$ and $f^s$, and the flux of $f$, i.e. $\mathbf{J}^l$ and $\mathbf{J}^s$, where $\mathbf{n}$ is the vector normal to the interface pointing outwards the solid surface. This theorem can be applied to any physical quantity following a local conservation Equation (3).
**Mass balance.** The mass flux corresponds to an advection term $J_f(x, t) = \rho u$ of the density field $f(x, t) = \rho$ by the velocity field $u$. Substituting this into the Equation (3) gives the conservation of mass, in three dimensions:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0. \qquad (5) $$

Considering an evolving interface $S(t)$ of local velocity $v_i$ between a liquid phase and a solid phase (subscripts $l$ and $s$, respectively), Figure 7, the jump relation, Equation (4) gives:

$$ (\rho_i - \rho_s)(v_i \cdot \mathbf{n}) = (\rho_i u_i \cdot \mathbf{n} - \rho_s u_s \cdot \mathbf{n}) \qquad (6) $$

$u_i$ and $\rho_i$ are the velocity and the density in the liquid phase at the position of the interface. $u_s$ is the velocity in the solid phase, which is zero (no advection in the solid) and $\rho_s$ is the density of the solid, supposed homogeneous.

By generalizing to any type of fluid, this jump relation determines the normal fluid velocity at the evolving interface:

$$ u_i \cdot \mathbf{n} = (v_i \cdot \mathbf{n}) \left( 1 - \frac{\rho_s}{\rho_i} \right). \qquad (7) $$

This condition for the normal velocity in the fluid translates the non penetration of the liquid phase in the solid. It must be completed by an expression for the velocity tangential to the interface in the liquid. Here, the tangential velocity of the fluid is zero at the dissolving boundary, as usual for a viscous fluid in contact with a solid wall, which reads $\mathbf{u} \wedge \mathbf{n} = \mathbf{0}$.

**Balance of solute or vapor concentration.** The local conservation relation (Equation (3)) applied to the mass solute concentration $c(x, t)$ reads:

$$ \frac{\partial c}{\partial t} + \nabla \cdot (J_c) = 0 \quad \text{with} \quad J_c = c\mathbf{u} - D\nabla c, \qquad (8) $$

where $\nabla$ is the gradient operator. As an example of the incoming flux $J_c$, we consider a volume element including the solid/liquid interface (see Figure 7) in dissolution. The chemical species can be transported by advection and diffusion. According to the Fick's law, the diffusive flux is proportional to the solute gradient: $J_c = -D\nabla c$, where $D$ is the solute diffusion coefficient. Applied to an incompressible flow and homogeneous diffusion coefficient, this equation identifies with the classical advection-diffusion equation for the concentration (solute or vapor):

$$ \frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = D\nabla^2 c. \qquad (9) $$

When considering a pure solid phase (i.e., concentration and density in the solid are equal) of the same chemical species as the solute, specifying the fluxes in Equation (4) for the concentration gives:

$$ (c_i - \rho_s)(v_i \cdot \mathbf{n}) = c_i u_i \cdot \mathbf{n} - D \nabla c |_{i} \cdot \mathbf{n} - \rho_s u_s \cdot \mathbf{n}. \qquad (10) $$

The diffusive flux $-D\nabla c|_i$ at the solid interface depends on the concentration gradient at the interface in the fluid. After simplification by assuming that the solid phase is immobile ($u_s = \mathbf{0}$), one finds:

$$ \rho_s v_i \cdot \mathbf{n} = c_i (v_i \cdot \mathbf{n} - u_i \cdot \mathbf{n}) + D (\nabla c|_i \cdot \mathbf{n}). \qquad (11) $$

Substituting the value of $u_i \cdot \mathbf{n}$ from Equation (7) into this expression gives a general boundary condition for the interface velocity $v_i$:

$$ (v_i \cdot \mathbf{n}) \rho_s \left( 1 - \frac{c_i}{\rho_i} \right) = D (\nabla c|_i \cdot \mathbf{n}). \qquad (12) $$

Moreover, the interface velocity is also related to the mass solute flux $\mathcal{R}$ given by Equation (1).
**Heat balance.** Finally, we express the conservation of energy. Similarly, the spatial variation of the temperature in the solid and fluid phases leads to thermal diffusion according to the Fourier's law. The heat flux due to conduction is equal to $-\lambda\nabla T$, where T is the temperature field and $\lambda$ is the thermal conductivity (different for solid and liquid phases). The enthalpy $h$ per unit mass is $h = (1/2)\mathbf{u}^2 + c_p T$ where $c_p$ is the specific heat capacity at constant pressure (different for each phases) and T is the temperature. In absence of source term, negligible fluid dissipation and for an incompressible flow, the heat balance is [100]:

$$ \frac{\partial(\rho c_p T)}{\partial t} + \nabla \cdot (\rho c_p T \mathbf{u}) = \nabla \cdot (\lambda \nabla T), \qquad (13) $$

where $\lambda$ is the thermal conductivity. The heat transport by diffusion is modeled using the Fourier law, $\mathbf{J}_T = -\lambda\nabla T$, where $\mathbf{J}_T$ is the heat flux. It has the dimension of a power per surface unit. If the properties of the fluid ($\rho$, $c_p$ and $\lambda$) are homogeneous, this last equation resumes to the classic advection-diffusion equation for the temperature field:

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \kappa \nabla^2 T, \qquad (14) $$

where $\kappa = \lambda/(\rho c_p)$ is the thermal diffusivity.

In general, a phase change reaction leads to an energy variation due to the breaking of molecular bonds. If a mass $m_x$ has experienced a phase change from a phase 1 to 2 at constant temperature and pressure, the difference of enthalpy reads: $\Delta H = m_x h_r$, where $h_r$ is the reaction enthalpy per mass unit. For condensed phases (liquid or solid), the reaction enthalpy reads as a function of the heat capacity per unit mass $c_p$:

$$ h_r = c_{p,2} T_i - c_{p,1} T_i, \qquad (15) $$

where $T_i$ is the interface temperature. For the melting of a pure substance at constant temperature, $T_i$ and $h_r$ identify to the melting temperature $T_m$ and to the latent heat $L_m$. The melting is an endothermic process and $L_m$ is always positive, as well as in sublimation. In contrast the reaction enthalpy for dissolution can be positive (e.g. salt NaCl or sucrose) or negative (e.g. like caustic soda NaOH) [101].

Applying Kotchine's theorem to energy conservation (Equation (13)) gives:

$$ [\rho_i c_{p,1} T_i - \rho_s c_{p,s} T_m] \mathbf{v}_i \cdot \mathbf{n} = [\rho_l c_{p,l} T_i \mathbf{u}_i - \lambda_l \nabla T]_{l,i} + \lambda_s \nabla T|_{s,i} \cdot \mathbf{n}. \qquad (16) $$

where subscript l and s, stand for liquid and solid phase, respectively. Using Equation (7) for the expression of $\mathbf{u}_i \cdot \mathbf{n}$ gives:

$$ [\rho_i c_{p,l} T_m \left[1 - \left(1 - \frac{\rho_s}{\rho_i}\right)\right] - \rho_s c_{p,s} T_i] \mathbf{v}_i \cdot \mathbf{n} = [-\lambda_l \nabla T]_{l,i} + \lambda_s \nabla T|_{s,i} \cdot \mathbf{n}. \qquad (17) $$

Using the definition of the reaction enthalpy Equation (15), where 2 and 1, stands for liquid and solid, subscripts l and s, respectively we obtain the following relation, which identifies to the classic Stefan condition, when $T_i = T_m$ and $h_r = L_m$:

$$ \rho_s h_r \mathbf{v}_i \cdot \mathbf{n} = [-\lambda_l \nabla T]_{l,i} + \lambda_s \nabla T|_{s,i} \cdot \mathbf{n}. \qquad (18) $$

This last equation determines the rate of phase change of a pure substance, which is controlled by the difference of thermal energy flux on both sides of the interface.
**Mass transfer versus heat transfer processes.** The same formalism apply to dissolution/sublimation and to the melting of multi-component systems. In both cases, the combination of Equations (1), (12) and (18) are coupled because the fluid physical properties (densities, diffusion coefficient, the saturation concentration, the latent heat, the heat capacity, ...) are function of the temperature or of the composition. The combined resolution of the three equations is not addressed in the literature, to our knowledge, for ablation problems. However, with the hypothesis of local thermodynamic equilibrium (fast kinetics), the phase change of binary or multi-components systems can be modeled, using energy conservation and mass conservation of each specie [102].

Nevertheless, in most practical cases, the interface motion is governed by one limiting process, either the temperature or the solute transport. For example, for dissolution, the phase change rate is assumed to be not controlled by the energy input, but by the concentration field at the interface [91,95,103]. Similarly, for the melting of pure ice into salt water, the local salt modifies the melting temperature, which can lead to complex phenomena, where the hypothesis of thermodynamic equilibrium is not valid [90,92]. In some conditions, the phase change can be controlled by the salt concentration and the ice melting problem is closer to a dissolution problem [104]. In this review, we do not address in details the case of melting of pure water ice in salt water, which is a rich topic with a large variety of phenomena. This topic of crucial environmental importance on Earth is covered by recent reviews [105–107].

In many natural dissolution phenomena, the temperature can be considered constant and set by the external conditions, i.e. the heat transport is fast enough to compensate the dissolution reaction enthalpy [108]. In this case, the rate of phase change is given by Equations (1) and (12). Combining these two equations leads to the boundary condition for the dissolving interface:

$$ \rho_s \mathbf{v}_i = -\alpha(c_{\text{sat}} - c_i)\mathbf{n} = \frac{1}{1 - c_i / \rho_i} D(\nabla c|_i \cdot \mathbf{n})\mathbf{n}. \quad (19) $$

This boundary condition combines the concentration value $c_i$ and the concentration gradient $\nabla c|_i$ at the interface. By introducing the typical length scale $\delta_c$ of the concentration gradient, i.e. $\nabla c \sim (c_b - c_i)/\delta_c$ where $c_b$ is the far field concentration, we can define the Damköhler number

$$ Da = (\alpha \delta_c) / D, \quad (20) $$

which compares the relative importance of the chemical reaction (middle-hand side of Equation (19)) to the diffusive transport (right-hand side of Equation (19)). In the case of natural rocks, for which $c_i \ll \rho_i$, the interfacial concentration is given by $c_i = c_b/(1+Da) + c_{\text{sat}}Da/(1+Da)$ [109]. If the chemical reaction rate is the limiting process at low values of Da compared to one, then $c_i \sim c_b+Dac_{\text{sat}}$. This corresponds to the case of slow dissolution reaction, valid for most weakly soluble minerals, as the low values of $\alpha$ correspond often to small values of $c_{\text{sat}}$ [110]. In contrast, large values of Da give $c_i \approx c_{\text{sat}}$ and correspond to fast dissolution regime. This occurs for very soluble substances like salt or sugar in water. In this case, the advection–diffusion equation (Equation (9)) is thus solved with the boundary conditions $c_i = c_{\text{sat}}$ and Equation (12), which relates the normal concentration gradient to the interface velocity.

Finally, we remark, that the previous boundary conditions determine the normal projection of the interface velocity ($\mathbf{v}_i \cdot \mathbf{n}$) and do not constrain the tangential part. We assume thus, that the interfacial velocity is directed along the local normal to the interface $\mathbf{n}$. This condition has important geometric consequences in the evolution of interfaces undergoing phase changes [111].
## 3.2. Models for melting, dissolution and sublimation

For an evolving interface undergoing a phase change, the interface velocity $v_i$ is thus coupled by boundary conditions to advection-diffusion equations, which set the transport of the solute concentration field and/or of the temperature field, by the velocity field $\mathbf{u}(\mathbf{x}, t)$ in the fluid phase. To close the problem, the determination of the velocity field is given by the Navier-Stokes equations, which express the local conservation of momentum. However, the spatial variations in composition due to the solute concentration and temperature affect the fluid properties. The density, the fluid viscosity, the diffusion coefficient, the diffusivity, the latent heat then vary spatially. The heterogeneity of the properties makes the question difficult to address theoretically. For sake of simplicity, we assume hereafter weak solute concentrations and weak variations of fluid properties with temperature. The density variations, which depend on concentration and temperature, are only taken into account when multiplied by the gravitational term (Boussinesq hypothesis). This is true for small variations in concentration and temperature. Finally, we assume that the properties of the fluid are homogeneous.

### 3.2.1. Melting/solidification of a pure substance

In these conditions, the model describing the melting/solidification of a pure substance in presence of a flow in the liquid phase is given by the following set of partial equations:

$$ \frac{\partial T}{\partial t} = \kappa_s \nabla^2 T \quad (\text{solid phase}) \qquad (21) $$

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \kappa_1 \nabla^2 T \quad (\text{liquid phase}) \qquad (22) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad (\text{liquid phase}) \qquad (23) $$

$$ \rho_{1,0} \frac{\partial \mathbf{u}}{\partial t} + \rho_{1,0} (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla P + \mu \nabla^2 \mathbf{u} + \rho_{1,0} (1 - \beta(T - T_0)) \mathbf{g} \quad (\text{liquid phase}). \qquad (24) $$

Balance in total energy, mass, and momentum are described by the heat transport equation (Equations (21) and (22) in solid and liquid phases, respectively), by the continuity equation for an incompressible flow (Equation (23)) and the Navier-Stokes equations (Equation (24)) for a viscous flow. $\beta$ is the thermal expansion coefficient (for example $\beta = 2.06 \times 10^{-4} \text{ K}^{-1}$ for liquid water at 20 °C), $P$ is the pressure field, $\mathbf{g}$ is the gravity acceleration and $\mu$ the dynamic viscosity of the liquid phase (the kinematic viscosity $\nu$ is equal to $\mu/\rho$). $\kappa_s$ and $\kappa_1$ are the thermal diffusivities in the solid and liquid phase, respectively. The liquid density $\rho_{1,0}$ is assumed to be homogeneous and constant in the Navier-Stokes equation (Equation (24)) according to the Boussinesq approximation, where the density is only slightly modified in the gravitational term. $T_0$ is the temperature for which $\rho_1(T_0) = \rho_{1,0}$. Note that the coefficient $\beta$ corresponds to a linearization of the equation of state in the liquid phase, whose knowledge is required to solve the Navier-Stokes equation.

These equations must be associated to the following boundary conditions at the melting interface, whose velocity is $v_m$:

$$ \mathbf{u} \cdot \mathbf{n} = v_m \left( 1 - \frac{\rho_s}{\rho_i} \right) \qquad (25) $$

$$ \mathbf{u} \wedge \mathbf{n} = \mathbf{0} \qquad (26) $$

$$ T_i = T_m \qquad (27) $$

$$ \rho_s L_m v_m \cdot \mathbf{n} = -\lambda_1 \nabla T |_{l,i} \cdot \mathbf{n} + \lambda_s \nabla T |_{s,i} \cdot \mathbf{n}. \qquad (28) $$
<table><caption>**Table 1.** List of relevant dimensionless numbers for melting phenomena</caption><thead><tr><th>Melting</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>Stefan number</td><td>St = L<sub>m</sub>/(c<sub>p</sub>ΔT)</td><td>Competition of the latent heat with the stored thermal energy<br>One Stefan number for each phase</td></tr><tr><td>Prandtl number</td><td>Pr = v/κ<br>Pr ~ 10 (water)</td><td>Momentum transport by viscosity<br>versus heat diffusion</td></tr><tr><td>Density ratio</td><td>r<sub>ρ</sub> = ρ<sub>l</sub>/ρ<sub>s</sub><br>r<sub>ρ</sub> ~ 1.09 (water)</td><td>Ratio liquid phase density to solid<br>phase density at the interface</td></tr><tr><td>Reynolds number</td><td>Re = (δ<sub>U</sub> U)/ν</td><td>Competition between flow advection<br>and the viscous diffusion</td></tr><tr><td>Thermal Péclet number</td><td>Pe<sub>T</sub> = (δ<sub>U</sub> U)/κ = Pr Re</td><td>Competition between flow advection<br>and the temperature diffusion</td></tr><tr><td>Thermal Rayleigh number</td><td>Ra = βΔTgδ<sub>T</sub><sup>3</sup>/κν</td><td>Competition between buoyancy driven<br>advection and diffusion phenomena</td></tr><tr><td>Interface motion number</td><td>r<sub>m,U</sub> = κT<sub>m</sub>/ρ<sub>s</sub>L<sub>m</sub>δ<sub>T</sub>U</td><td>Ratio of the melting interface velocity<br>to the typical flow velocity</td></tr></tbody></table>

$δ_U$ is the characteristic length scale of the flow, either the typical scale of velocity gradient or the thickness of a viscous boundary layer. $δ_T$ is the characteristic scale of the temperature field (such as the thickness of a boundary layer). Note that the Stefan number may be defined differently in the literature.

For the other interfaces, the relevant boundary conditions depend on the considered problem, like for example cancellation of the velocity field on solid walls and of thermal gradient for insulating walls. The pressure field depends also on the external boundary conditions.

Moreover, as the solid/fluid interface evolves the domain of each phases as a function of time, the interface evolution introduces thus geometrical non-linearities by domain variation and induces computational challenges for direct resolutions.

Relevant dimensionless numbers can be introduced to define flow and melting regimes, see Table 1. The density ratio $ρ_i/ρ_s$, where $ρ_i$ is taken at the melting temperature, determines the normal velocity at the melting interface. For each phase, we can define a Stefan number $St = L_m/(c_pΔT)$, with $ΔT$ the typical temperature variation in the problem. This number compares the latent heat to the stored thermal energy, the lower the $St$ value, the higher the melting speed. Note that the inverse definition of the Stefan number is also commonly used in the literature. The Prandtl number, $Pr = v/κ_1 = μ/(ρκ_1)$ is the ratio of kinematics viscosity to thermal diffusivity. $Pr$ for water is about 13.5 at 0 °C and 6.99 (smaller) at 20 °C, because the fluid viscosity decreases with temperature, whereas the thermal conductivity increases. The momentum transport is thus faster than the heat transport, but without a clear separation of scales. For the hydrodynamics, we can introduce the Reynolds number $Re = (δ_U U)/ν$, where $δ_U$ is a characteristic length and the thermal Rayleigh number $Ra_T = (Δρgδ_T^3)/(μκ_1)$ with $Δρ = ρ_0βΔT$, where $δ_T$ is a characteristic scale of the temperature field. $δ_T$ and $δ_U$ depend on the considered system and flow conditions. The thermal Péclet number $Pe_T = (δ_U U)/κ_1 = Pr Re$ compares in Equation (22) the heat advection to the heat diffusion. For sufficiently large Rayleigh numbers and unstable density stratification, the liquid can move by thermal convection: the onset of convection (Rayleigh-Bénard instability)
between two horizontal solid plates maintained at constant temperatures, is predicted for a Rayleigh number above a critical value equal to 1707.8 [100].

### 3.2.2. Dissolution/precipitation of a solid in a solvent

For the dissolution/precipitation of a solid in a liquid phase (solution) consisting of a solvent and a solute of concentration c, the set of equations to solve in the liquid phase is, with the previous simplifying hypotheses:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (29) $$

$$ \frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = D \nabla^2 c \qquad (30) $$

$$ \rho_{1,0} \frac{\partial \mathbf{u}}{\partial t} + \rho_{1,0} (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla P + \mu \nabla^2 \mathbf{u} + \rho_{1g} \mathbf{g}. \qquad (31) $$

Equations (29) and (31) are again the incompressible Navier-Stokes equation in the Boussinesq approximation, the fluid density is assumed homogeneous and equal to the mean value $\rho_{1,0}$, except for the gravity term. The equation of state of the liquid phase relates the local density to the temperature and the chemical composition (here the solute concentration). For negligible temperature variations, the liquid density is related to the concentration by the relation: $\rho_l = c(\rho_{sat} - \rho_{1,0})/(c_{sat}) + \rho_{1,0}$, where $\rho_{sat}$ is the saturation liquid density. This linear dependency of density with concentration is valid for most solutes. The increase in density with c leads to buoyancy effects through the term $\rho_{lg}$. The boundary conditions at the solid/liquid interface are:

$$ \mathbf{u} \cdot \mathbf{n} = v_i \left( 1 - \frac{\rho_s}{\rho_i} \right) \qquad (32) $$

$$ \mathbf{u} \wedge \mathbf{n} = \mathbf{0} \qquad (33) $$

$$ \rho_s \mathbf{v}_i = -\alpha(c_{sat} - c_i)\mathbf{n} \qquad (34) $$

$$ \rho_s \mathbf{v}_i \left( 1 - \frac{c_i}{\rho_i} \right) = D (\nabla c)_i \cdot \mathbf{n} \mathbf{n}, \qquad (35) $$

where $\alpha$ is the chemical dissolution rate coefficient.

To close the set of equation, like for melting, additional boundary conditions are required at non dissolving interfaces of the domain. Then, the dissolution velocity $v_i$ at each point of the interface can be in principle computed to predict the erosion rate and a possible pattern emergence. Several dimensionless numbers are relevant to distinguish the hydrodynamic and dissolution regimes, see Table 2. We have already introduced the Damköhler number $Da = \alpha\delta/D$, which compares the chemical kinetics of the dissolution reaction with the diffusion transport at the solid/liquid interface. The Schmidt number $Sc = v/D = \mu/(\rho D)$ is the ratio of the momentum diffusion due to viscosity to the solute diffusion. Note that for the dissolution of ions or organics substances in liquid water, the Schmidt number is always of the order of 1000, i.e. the particle diffusion is always slow compared to the momentum diffusion. Consequently, the concentration boundary layers are usually significantly thinner than the viscous boundary layers, because on a given time scale momentum diffuses farther than solute $\delta_c \sim (DT)^{1/2}$, $\delta_v \sim (vT)^{1/2}$ and $D \ll v$. In addition, the ratios $c_{sat}/\rho_{sat} \approx c_{sat}/\rho_i$ and $\rho_i/\rho_s$ are also useful to define the weak concentration regime and the density jump at the interface, respectively. To characterize the hydrodynamics, the Reynolds number $Re = \delta_U U/v$ evaluates the importance of the nonlinear term in the Navier-Stokes equation (velocity field advection) in front of the viscosity term, where $\delta_U$ and $U$ are the characteristic length and velocity of the flow. Using $\delta_c$ the characteristic scale for the concentration field, usually small in front of $\delta_U$, we can define the solutal Rayleigh number $Ra = [(\rho_i - \rho_b)g\delta_c^3]/(\mu D)$, where $\rho_b$ is the liquid density of the “bath” i.e., far from
<table><caption>**Table 2.** List of relevant dimensionless numbers for dissolution phenomena</caption><thead><tr><th>Dissolution</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>Damkhöler number</td><td>Da = αδ<sub>c</sub>/D</td><td>Competition between chemical reaction rate and diffusive transport</td></tr><tr><td>Schmidt number</td><td>Sc = v/D<br>Sc ~ 1000<br>(ions in water)</td><td>Momentum transport by viscosity versus solute diffusion</td></tr><tr><td>Density ratio</td><td>r<sub>ρ</sub> = ρ<sub>sat</sub>/ρ<sub>s</sub><br>r<sub>ρ</sub> ~ 0.55<br>(solid and saturated water)</td><td>Ratio liquid phase density (solute + solvant) to solid phase density of a soluble body</td></tr><tr><td>Concentration number</td><td>r<sub>c</sub> = c<sub>sat</sub>/ρ<sub>sat</sub><br>r<sub>c</sub> ~ 0.2 (salt)</td><td>Ratio of the saturation concentration to the saturation liquid density. Determine the condition of weak concentration</td></tr><tr><td>Reynolds number</td><td>Re = (δ<sub>U</sub> U)/ν</td><td>Competition between flow advection and the viscous diffusion</td></tr><tr><td>Péclet number</td><td>Pe = (δ<sub>U</sub> U)/D = ScRe</td><td>Competition between transport due to flow advection and the concentration diffusion.</td></tr><tr><td>Solutal Rayleigh number</td><td>Ra = (ρ<sub>sat</sub> - ρ<sub>b</sub>)gδ<sub>c</sub><sup>3</sup>/Dμ</td><td>Competition between buoyancy driven advection and diffusion phenomena</td></tr><tr><td>Interface motion number</td><td>r<sub>d,U</sub> = Dc<sub>sat</sub>/ρ<sub>s</sub>δ<sub>c</sub>U</td><td>Ratio of the dissolving interface velocity to the typical flow velocity</td></tr></tbody></table>

δ<sub>U</sub> is a characteristic scale of the flow, either the typical scale of the velocity gradient or the thickness of a viscous boundary layer, and δ<sub>c</sub> is the characteristic scale of the concentration field.

the dissolving interface. Like the thermal Rayleigh number we introduce in Section 3.2.1, the solutal Rayleigh number determines the condition of appearance of convective motion by a hydrodynamic instability mechanism [21,109,112]. Finally, the Péclet number Pe = (δ<sub>U</sub> U)/D compares advection with diffusion. The Péclet number is directly related to the Reynolds and Schmidt numbers Pe = ScRe.

### 3.2.3. Sublimation/condensation between a solid and a gas phase

The hydrodynamic and mass transfer equations for the sublimation of a solid in an incompressible fluid whose composition is different from that of the sublimating material (sublimation of water ice in air for example) are identical to Equations (29) and (31). The fluid is subject to small temperature variations (constant ρ<sub>l,0</sub>). The transport of the ablated species within the fluid is treated as a passive scalar, satisfying an advection/diffusion equation for the concentration (Equation (30)). The boundary conditions (34)–(35) are identical to those for dissolution, with α representing the sublimation reaction rate coefficient. The sublimation velocity will be influenced by the type of flow (natural or forced convection) and the nature of the flow (laminar or turbulent). Special attention must therefore be paid to the turbulence
<table><caption>**Table 3.** List of relevant dimensionless numbers for sublimation phenomena</caption><thead><tr><th>Sublimation</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>Damkhöler number for sublimation</td><td>Da<sub>s</sub> = αδ<sub>s</sub>/D</td><td>Competition between sublimation reaction rate and diffusive transport</td></tr><tr><td>Schmidt number</td><td>Sc = v/D<br>Sc ~ 1<br>(for gas)</td><td>Momentum transport by viscosity versus molecular diffusion</td></tr><tr><td>Density ratio</td><td>r<sub>ρ</sub> = ρ<sub>sat</sub>/ρ<sub>s</sub></td><td>Ratio gaseous phase density (gas + sublimated specie) to solid phase density of a soluble body</td></tr><tr><td>Relative humidity</td><td>RH = p<sub>vap</sub>/p<sub>sat</sub><br>RH &lt; 1 sublimation</td><td>Determine the conditions favorable to sublimation of the specie</td></tr><tr><td>Reynolds number</td><td>Re = (δ<sub>U</sub>U)/ν</td><td>Competition between flow advection and the viscous diffusion</td></tr><tr><td>Péclet number</td><td>Pe = (δ<sub>U</sub>U)/D = ScRe</td><td>Competition between flow advection and the concentration diffusion</td></tr><tr><td>Mass Rayleigh number</td><td>Ra = (ρ<sub>sat</sub> - ρ<sub>b</sub>)gδ<sub>s</sub><sup>3</sup>/Dμ</td><td>Competition between buoyancy driven advection and diffusion phenomena</td></tr><tr><td>Interface motion number</td><td>r<sub>s,U</sub> = Dc<sub>sat</sub>/ρ<sub>s</sub>δ<sub>s</sub>U</td><td>Ratio of the sublimating interface velocity to the typical flow velocity</td></tr></tbody></table>

$δ_U$ is a characteristic scale of the flow, either the typical scale of the velocity gradient or the thickness of a viscous boundary layer. $δ_s$ is the characteristic scale of the concentration field.

models used. As with dissolution, dimensionless numbers for phase change in sublimation are defined in Table 3. The Damkhöler number for sublimation, $Da_s = αδ_s/D$, is used to compare sublimation kinetics with diffusion transport at the gas/solid interface. The Schmidt number is the ratio of momentum diffusion to viscosity and molecular diffusion in the gas. For gases, this number is often close to 1, i.e. 1000 times lower than for dissolution, meaning that molecular diffusion is comparable to momentum diffusion. In this case, the velocity and concentration boundary layers are of the same order of magnitude. In the case of sublimation, the $r_ρ = ρ_{sat}/ρ_s$ ratio is very small, indicating a concentration jump at the interface. Temperature and pressure of the environment control the value of $p_{sat}$. Heat can come from radiation, which would add a source term in volume in the heat equation that has to be solved. To ensure a weak concentration regime, in the case of sublimation we are interested in the specific humidity, $RH = p_{vap}/p_{sat}$, where $p_{vap}$ is the vapor pressure of the specie and $p_{sat}$ the saturated pressure of the specie. A low value of this number corresponds to a dry atmosphere (and therefore favorable to sublimation). As in the case of dissolution, the Reynolds number is used to assess the competition between viscous diffusion and advection by the flow and the combination of Sc and Re gives the Pe number, which compares concentration transport by flow advection and diffusion. If buoyancy forces are driving the flow, the same Rayleigh number Ra as in the dissolution case can be used to estimate the relative strength of buoyancy due to concentration/density variations and diffusion.
### 3.2.4. Melting, dissolution and sublimation. Analogies and differences

The dissolution/sublimation velocity is set by the normal to the surface component of the concentration gradient (Equation (12)) and the concentration field in the fluid phase is found by solving an advection diffusion equation for the concentration. As seen previously, for large values of the Damköhler number $Da$, the concentration at interface is close to the saturation concentration ($c_i = c_{sat}$). The problem becomes thus mathematically identical to the melting situation in case of negligible temperature diffusion in the solid (often valid for large Stefan number $St = L_m/(C_p\Delta T)$), where the advection-diffusion equation of the temperature is solved with the boundary conditions $T_i = T_m$. In that case, Equation (18) relates the interface velocity to the normal temperature gradient. In both cases, the kinetics of the molecular mechanisms (attachment and detachment) is assumed fast enough to be neglected and the interface state corresponds to the thermodynamics equilibrium, when the coexistence between two phases occurs. The motion of the interface is yet an out-of-equilibrium process caused by the difference of fluxes on each side of the interface. Therefore, dissolution and melting patterns can often be paralleled [69], as they arise qualitatively by the same mechanisms. Nevertheless, some of dimensionless numbers have always very different orders of magnitude. The Schmidt number $Sc$ is about 1000 for dissolution, whereas the Prandtl number $Pr$ is about 10, meaning that the heat diffuses faster than the concentration and that the diffusion of both quantities is less efficient than momentum diffusion. Typically, the sizes of the concentration structures like plumes, eddies, boundary layers are smaller than those corresponding to the temperature.

Another important difference, is related to the external boundary conditions of the problem. Despite the mathematical analogy, the saturation concentration does not play the same role than the melting temperature in practical cases. Usually in dissolution experiments in a closed domain, the process stops, when the liquid concentration reaches the saturation value, due to accumulation of the solute. Providing fresh water for concentration is the most suitable way to avoid this effect. For melting examples, an energy input is externally supplied to sustain the phase change and if it is maintained constant, the melting continues until the solid has completely disappeared. However, it can also stop if the system is closed without loss of energy, allowing it to reach equilibrium temperature. In both case, for the ablation process to be effective, the system must be out of equilibrium.

For dissolution, the kinetics effects are often considered more relevant that for melting. For most low solubility minerals, such as quartz and other silicates [91], the rate of dissolution is often considered to be limited by the chemical kinetics. For minerals of intermediate solubility, such as limestone and gypsum, the interaction between transport and reaction has been addressed in a recent review [113] of dissolution in porous media, wormhole formation and fracture growth by dissolution. The reaction seems to be confined to a small distance necessary to reach conditions close to saturation. In addition, the role of the reaction in the presence of a flow, and thus of the coefficient $\alpha$ in the dissolution rate, has only been directly shown in very few experiments. Using a microfluidic flow experiment combined to dedicated numerical simulations, Dutka et al. [114] find an intermediate situation ($Da \sim 3$) for dissolution of gypsum samples. In this example, the dissolution rate is mainly limited by diffusion, but is significantly influenced by the reaction rate coefficient $\alpha$. For very soluble materials like salt (NaCl) or sucrose, to our knowledge, in macroscopic experiments the role of the reaction is not perceptible for dissolution dynamics and pattern formation [21,109,112]. Therefore, most macroscopic dissolution experiments belong to the diffusion limited regime, i.e. large $Da$. The case of the generation of meter scale natural conduits in limestone caves is debatable and could depends if the dissolution occurs in acid conditions (large $Da$) or in alkaline conditions (small $Da$) [115].

Then, according to Equation (1), the kinetics define a maximal dissolution velocity $v_{d,max} = \alpha c_{sat}/\rho_s$. Although, the values of $\alpha$ are not always well determined experimentally, the
order of magnitude of the maximal dissolution velocity, $v_{d,max}$, is of the order of $10^{-5}$ m/s, for fast dissolving materials like salt [103]. Consequently in the presence of a flow at a macroscopic scale, the interface velocity is usually very small in front of the hydrodynamic velocity, in the case of dissolution or sublimation. The evolution of the interface and of the hydrodynamic acts thus with two separated time scales. For the melting of a pure substance in extreme situations, the separation of time scales is not always guaranteed. A dimensionless number to quantify this time scale separation can be defined as the ratio between the interface velocity to the typical flow velocity: $r_{d,U} = D c_{sat} / (\rho_s \delta_c U)$ for dissolution and $r_{m,U} = \kappa T_m / (\rho_s L_{sat} \delta_T U)$ for melting respectively. These numbers are usually very small compared to one, but must be evaluated in each situation.

We note also an important analogy between dissolution/sublimation and melting in the context of erosion processes. Once thermodynamic equilibrium is broken, the flow in the fluid phase always affects the rate of phase change. At a melting interface, an increase of the local flow velocity increases the shear rate in the thermal boundary layer or the concentration boundary layer. The thickness of boundary layers decreases and the local phase change rate increases. Contrary to mechanical erosion, for which erosion occurs when the shear stress overcomes a threshold, there is no such threshold for dissolution or melting. An infinitesimally weak velocity can thus participate to the emergence of a dissolution/sublimation or melting pattern.

Finally, as we pointed out before, the sublimation of a solid in a atmosphere, where the sublimated vapor is in small proportion is formally analogous to the dissolution of a solid in a solvent. As the theoretical description is identical, the same dimensionless numbers arise naturally, by considering the concentration field of vapor in the atmosphere (the index 1 (liquid phase) by *g* (gas phase)). However, the physical properties of the liquid and gas phases differ strongly, leading to different orders of magnitude for the physical coefficients $c_p$, *ν*, $\kappa$, *D*,... In addition, because the heat capacity of gases is smaller, the coupling between concentration and temperature fields is stronger. Sublimation is usually modeled by addressing the vapor transport and taking into account the temperature variations, which change in turn the saturation vapor pressure, with the hypothesis of local equilibrium of the interface [116]. As the molecular mechanisms at the origin of diffusion and viscosity are similar in gas (which is not the case for liquids), the Schmidt number and the Prandtl numbers are of order one ($Sc$ = 0.63 for water vapor in air and $Pr$ = 0.7 for air) [101]. Compressibility effects can also become significant, modifying the hydrodynamic description.

### 3.3. *Choice of analog material for an experiment*

As many natural patterns created by a phase change often occur in complex situations with intermittent conditions and entangled physical mechanisms, model laboratory experiments are useful to identify their formation mechanism and dynamics. However, the timescales in which they appear and evolve are often very long compared to human time scales. It is therefore necessary to accelerate the phenomenon in order to carry out experiments within a manageable time, typically one hour. The question is then whether the physics studied for fast evolution is qualitatively similar to the natural case, and whether the laboratory results can be extrapolated to the field, by examining the relevant dimensionless number. A list of these dimensionless numbers are summarized in the Tables 1–3.

For melting experiments by a fluid flow, most of experiments are performed with ice in contact with water. The melting/solidification temperature of 0 °C is relatively close to ambient temperature compared to other common pure materials and is accessible using standard laboratory equipment. The melting rate can be conveniently adjusted by varying the temperature of the water flow. Ice can also be used to mimic the melting of minerals or materials at high
temperatures. This has been done in an experiment where an impinging jet of molten material hits a solid wall [117]. Organic phase change materials are also used to investigate the interaction between phase change and melting, like cyclohexane [118] ($T_m$ = 6.38 °C), hexadecane [119] ($T_m$ = 17.7 °C) or polyethylene glycol wax [120] ($T_m$ = 19 °C). The latent heat and thermal parameters are of a similar order of magnitude, but the organic liquid are much more viscous and less thermally conductive. Consequently the Prandtl number $Pr = v/\kappa$ is large compared to one and of order 1000 instead of about 10 with water. Fluid mechanics experiments with phase change material have been also performed with liquid metals like gallium [121] ($T_m$ = 29.7 °C), which have larger densities and allow to reach small values of the Prandtl numbers (about $2 \times 10^{-2}$).

For phase changes by dissolution, the three most important types of soluble rocks are salt (Halite or sodium chloride, NaCl), gypsum ($CaSO_4 \cdot 2H_2O$) and limestone (Calcite, $CaCO_3$), whose solubility (approximately the saturation concentration) in water are around 360 g/L, 2.0 g/L and 15 mg/L, respectively. The solubility of limestone is strongly dependent on carbon dioxide concentration and the acidity, which is measured with the pH scale. These three rock examples belong respectively to three mineral families with similar chemical and physical properties, the halides, the sulfates and the carbonates, respectively. In general, the dissolution parameters, $c_{sat}$ and $\alpha$ will have similar value for components belonging to the same family.

The chemical dissolution/precipitation reactions are:

$$(NaCl)_s \rightleftharpoons Na^+ + Cl^-$$

for salt (sodium chloride) and:

$$(CaSO_4 \cdot 2H_2O)_s \rightleftharpoons Ca^{2+} + SO_4^{2-} + 2H_2O$$

for gypsum (dihydrated calcium sulfate).

They are described by simple dissolution/precipitation reaction. Limestone is the most common and important case, due to the large distribution all over the continents of regions shaped by the hydrochemical erosion of limestone called Karsts [122]. The chemical reactions for removal or deposition of the calcite are quite complex with several elementary steps [91,123]. However, they can be summarized in a single chemical equation:

$$CaCO_3 + CO_2 + H_2O \rightleftharpoons Ca^{2+} + 2HCO_3^-$$

Note that the reaction involves $CO_2$ as a gas dissolved in the aqueous phase, whose presence in the solution acidifies the water. Therefore, the rate at which calcite dissolves depends on exposure to the atmosphere and pH. To simplify, the limestone dissolves at low pH (acid) with the consumption of $CO_2$ and precipitates at high pH (basic) with the production of $CO_2$ especially for acid enough water, erosion is often assumed to be controlled by solute transport rather than dissolution kinetics [91]. The values of the relevant experimental parameters for the most common soluble materials are given in Table 4 for a temperature of about 20 °C. We note that the dissolution rate $\alpha$ conveying the chemical kinetics, is not well known, especially for fast dissolving materials, because the hydrodynamic methods (like the rotating sample one) depends on a laminar model for the boundary layer and could suffer from biases [124]. In absence of external flow, the use of a holographic method to follow the evolution of the density profile by diffusion provides a consistent method to measure $\alpha$ for gypsum [95]. As mentioned previously, dissolution of a gypsum sample in a microfluidic cell allows to measure $\alpha$ [114]. The obtained values are compatible with the holographic measurement without flow motion. Nevertheless, it appears that in general the largest values of $\alpha$ are associated to large saturation concentration, as it has been tested for the dissolution of some organic substances [110]. Sugar (sucrose) and similar organic compound like glucose or benzoic acid are convenient to perform dissolution experiments [21,125,126], because they can be molded to the prescribed shape after thermal
<table><caption>**Table 4.** Characteristic parameters for dissolution in fresh water at a temperature of about 20 °C</caption><thead><tr><th></th><th>Solid<br>density ρ<sub>s</sub><br>(kg·m<sup>-3</sup>)</th><th>Saturation<br>concentration<br>c<sub>sat</sub> (kg·m<sup>-3</sup>)</th><th>Saturation<br>density ρ<sub>sat</sub><br>(kg·m<sup>-3</sup>)</th><th>Dissolution<br>rate α<br>(m·s<sup>-1</sup>)</th><th>Diffusion<br>coefficient<br>D (min|max)<br>(m<sup>2</sup>·s<sup>-1</sup>)</th><th>Dynamic<br>viscosity<br>μ (min|max)<br>(Pa·s)</th></tr></thead><tbody><tr><td>Limestone<br>CaCO<sub>3</sub></td><td>2711</td><td>0.312</td><td>997.432</td><td>3 × 10<sup>-6</sup></td><td>1.4 × 10<sup>-9</sup></td><td>1.002 × 10<sup>-3</sup></td></tr><tr><td>Gypsum<br>CaSO<sub>4</sub>·2H<sub>2</sub>O</td><td>2320</td><td>2.04</td><td>999.16</td><td>2.7 × 10<sup>-6</sup></td><td>1.0 × 10<sup>-9</sup></td><td>1.002 × 10<sup>-3</sup></td></tr><tr><td>Salt<br>NaCl</td><td>2170</td><td>317.5</td><td>1200</td><td>5.0 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-9</sup><br>1.6 × 10<sup>-9</sup></td><td>1.00 × 10<sup>-3</sup><br>1.99 × 10<sup>-3</sup></td></tr><tr><td>Sugar<br>C<sub>12</sub>H<sub>22</sub>O<sub>11</sub></td><td>1540</td><td>887</td><td>1327</td><td>≥5 × 10<sup>-6</sup></td><td>5.2 × 10<sup>-10</sup><br>7.1 × 10<sup>-11</sup></td><td>1.00 × 10<sup>-3</sup><br>0.212</td></tr><tr><td>Sucrose</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

The density of fresh water is $ρ_0 = 997.12 \text{ kg·m}^{-3}$ and its dynamic viscosity $μ_0 = 1.002 \times 10^{-3} \text{ Pa·s}$ [101]. Limestone, data from Kaufmann et al. [127] at high carbon dioxide pressure. For limestone, $c_{sat}$ is strongly dependent on the carbon dioxide pressure and on the pH. This value is particularly high, $c_{sat} \approx 0.015 \text{ kg·m}^{-3}$ under standard conditions. Gypsum, data from Colombani et al. [95]. Salt, data from Alkattan et al. [103] and from the Handbook [101]. Sucrose, data from Cohen et al. [21], and others [128,129] and from the Handbook [101].

melting. However, viscosity increases by several orders of magnitude with concentration, which is not the case for mineral rocks, and some of the observations may not be relevant for soluble rocks. Plaster has the same chemical composition than gypsum. It has an intermediate solubility between salt and limestone and can be easily molded. This is why plaster has been used in several important studies of dissolution in water flows within a matter of hours [15,17,20]. Yet, the obtained samples are often porous, which may influence the surface patterning and the measured erosion rates. Furthermore, some experiments can be performed with synthetic large salt crystals (like NaCl, KCl or KBr), but the samples are expensive and limited to a few centimeters. Then, mineralogical samples from mines and quarries can be also used in experiments. However, they often show defects in composition if they contain impurities, and in mechanical properties due to grain boundaries or fractures.

Finally, we now consider the case of sublimation of a solid. Sublimation technique [130] has been used to obtain informations about heat transfer in experiments. Analogies have been made between mass and heat transfer. Naphthalene ($C_{10}H_8$), with its high sublimation rate at room temperature [131], has been used as sublimating substrate for its good casting and machining properties but it is considered carcinogenic and to our knowledge has not been used to study bedforms. Camphor has already been used to test the influence of material surface roughness on the triggering of the laminar-turbulent transition when objects coming from space enter the atmosphere. However, it requires high temperatures for sublimation [132]. Water ice under controlled cold and dry atmospheres appears to be the most obvious substrate. The use of water ice as the substrate for sublimation experiment has some advantages over the previous ones: it is not toxic, it requires less extreme temperatures, the material is inexpensive to produce and can be shaped in any desired manner. But here the duration of experiments can be a problem. For a fixed ambient and partial pressure, the saturation margin $p_{sat} - p_{vap}$ determines the sublimation velocity (see Table 5). Thus, water-ice experiments would have to run for weeks to achieve the same depth of ablation that a CO$_2$-ice experiment can achieve in hours. Under similar
<table><caption>**Table 5.** Characteristic parameters for sublimation in N₂ atmosphere</caption><thead><tr><th></th><th>Solid<br>density ρ<sub>s</sub><br>(kg·m<sup>-3</sup>)</th><th>Saturation<br>pressure<br>p<sub>sat</sub> (Pa)</th><th>Vapor<br>pressure p<sub>vap</sub><br>Pa (min|max)</th><th>Sublimation<br>velocity v<sub>s</sub><br>(m·s<sup>-1</sup>)</th><th>Diffusion<br>coefficient D<br>(m<sup>2</sup>·s<sup>-1</sup>)</th><th>Kinematic<br>viscosity ν<br>(m<sup>2</sup>·s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Svéa,<br>Antartica H<sub>2</sub>O</td><td>865</td><td>336</td><td>100–166</td><td>3 × 10<sup>-9</sup></td><td>2.9 × 10<sup>-5</sup></td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>Nantes, wind<br>tunnel CO<sub>2</sub></td><td>1560</td><td>1.4 × 10<sup>5</sup></td><td>40</td><td>1.4 × 10<sup>-6</sup></td><td>1.5 × 10<sup>-5</sup></td><td>1.5 × 10<sup>-5</sup></td></tr></tbody></table>

Kinematic viscosities are computed from Sutherland's law for the viscosity of gases. In both cases, sublimation velocities are estimated in the absence of wind. Water ice data from Bintanja [50] with RH ~ 0.3, CO₂ ice data with RH << 1 from Bordiec et al. [32].

environmental conditions, water ice will always sublimate about three orders of magnitude slower than CO₂ ice.

## 4. Dissolution, melting and sublimating patterns in experiments

Pattern can be shaped by an instability of the flow, be periodic with a characteristic wavelength, but recognizable shapes can also be driven by specific boundary conditions. Here, we have chosen to classify experiments as a function of the external forcing and boundary conditions, which, we believe, modify the coupling processes described in Section 3 that drive pattern formation. Because diffusion tends to homogenize concentration and temperature fields, it cannot drive the spatial variation of the concentration/temperature field that leads to pattern formation (but when associated with specific boundary conditions); a flow is required. This flow can be buoyancy driven by the phase change itself or externally driven (or both). When externally driven, the flow can be deep, as in caves or oceans, or shallow and possibly interacting with the free surface, which is the usual case of a surface exposed to rainfall. Riverbeds can be subject to both deep and shallow flows.

### 4.1. Buoyancy driven flow and associated bedforms

Melting or dissolution of a body in quiescent water leads to density stratification at the eroding surface, either because the solution density depends on temperature or on solute concentration. This density stratification may be stable or unstable with respect to gravity. In either case, this density stratification drives a buoyant flow (but in the singular case of an infinite and perfectly flat surface). Experiments usually consist of simply following the dissolution of a block of sugar (glucose, sucrose, or a mixture) [21,48,108,133,134], salt (NaCl), or plaster [21], or the melting of an icy block in a tank of water [135].

#### 4.1.1. *Stable density stratification*

The dissolution of axisymmetric sugar bodies in quiescent water has been studied experimentally by Pegler and Davies Wykes [108] and by Huang et al. [48] for the case of stable density stratification (Figure 8). In this case, the erosion creates a thin convective boundary layer of concentration that grows while flowing along the body surface. These experiments show that an initial upright cylinder sharpens with time while receding [48] whereas an initial cone blunts with time [108]; the erosion dynamics depends on the initial shape [108]. The erosion dynamics and shapes are well described by a boundary layer model derived by Acrivos [136], in the limit of
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: dissolution of a candy cone (from Pegler and Davies Wykes [108]). Copyright Cambridge University Press 2020. The cone slowly blunts with time and the receding velocity slowly decreases with time. Right: A block seeded with vertical pores dissolves and form an array of sharp pinnacles (from Huang et al. [48]), under the CC BY 4.0 license.</figcaption></figure>

high Schmidt number, i.e., the concentration boundary layer is within the viscous layer [108,137]. They have shown in particular that in this 2D-model convex upward shapes sharpen with time with a receding rate that increases with time, whereas triangles and concave upward shapes blunt with time with a receding rate that decreases with time [137]. For simple shapes, self similar solutions exist but no asymptotic, stationary shape [108,137]. Huang et al. also interestingly “speculate that this mechanism contributes to the formation of pinnacles in nature”, such as tsingys or stone forests [48], the *sharp patterns* described in Section 2.1, Figure 1(b). Weady et al. studied the melting of initial vertical cylinder in a tank of water at different temperature $T_\infty$ [135]. These cylinders sharpen to form upward pinnacles when $T_\infty \gtrsim 7$ °C and downward pinnacles when $T_\infty \lesssim 5$ °C, because water density is maximum around 4 °C [135]. At intermediate temperatures, scallop-like patterns form at the ice surface. Authors explain these patterns as imprints of a Kelvin-Helmholtz instability between the rising fluid near the surface and the sinking outer flow [135].

### 4.1.2. *Unstable density stratification*

If the density stratification induced by the dissolution or melting of an object is unstable rather than stable, i.e. if the density gradient is opposed to the gravity direction, the results are significantly different, both at the global scale and at a smaller scale, where patterns appear. The first observations of this phenomenon date back to at least the beginning of the 20th century, when Schürr described the appearance of “convection stripes” on the surface of immersed soluble crystals in 1904 [138,139]. The appearance of patterns on the surface of solids under the effect of convection flow generated by their dissolution [140,141] or melting [142,143] was subsequently observed, but the first quantitative studies of patterning date from the late 1990s. Drawing an analogy with turbulent thermal convection, Sullivan et al. proposed that a solutal Rayleigh-Bénard instability determines the thickness of the concentration boundary layer, and therefore the dissolution rate, but also the characteristic size of the patterns. It proved to be in good agreement with experiments conducted on horizontal salt crystals [112,144,145].

This can be explained with the schematic in Figure 9. When a soluble solid is brought into contact with water, a boundary layer of concentration is created around it, which grows initially by diffusion. At downward facing surfaces, the boundary layer eventually destabilises when buoyancy overcomes diffusion in concentration transport, i.e., a solutal Rayleigh-Bénard
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left: Shadowgraph of a plaster block in water showing the sinking plumes of the solutal convection. The bottom surface is 20 cm long. Credits: Martin Chaigne. Right: Sketch of the solutal Rayleigh-Bénard instability and its imprint on the block (from Cohen et al. [21], with the permission of the American Physical Society).</figcaption></figure>

instability happens, analogous to the thermal Rayleigh-Bénard instability. The instability threshold is given by a criterion involving the solutal Rayleigh number:

$$ \frac{\Delta\rho g \delta_c^3}{\rho_b v D} = Ra_c, \qquad (36) $$

with $\Delta\rho$ the difference in density between the saturated fluid and the fluid in the bath, $\eta$ the dynamic viscosity of the fluid, $D$ the diffusion coefficient of the solute in the fluid, and $Ra_c$ a critical value of the solutal Rayleigh number. Using the values in Table 4 and taking $Ra_c = 1101$, which is the theoretical value for the onset of the Rayleigh-Bénard instability with mixed slip and no-slip boundary conditions [100], we find that the thickness $\delta_c$ is of the order of 0.1 mm for salt, 0.4 mm for gypsum and between 0.7 and 2 mm for limestone. Experimental values of the dissolution rate rather suggest $Ra_c \approx 100$ [112,144], which modifies the estimated thickness of the boundary layer by a factor of 2. The nonlinear evolution of the destabilized boundary layer then leads to the emission of thin plumes of solute-laden fluid. Their typical thickness is that of the boundary layer, while the typical spacing between two plumes is of the order of a few boundary layer thicknesses [21,144]. In Figure 9, they are visualised using the shadowgraph technique under a block of plaster suspended in water at an angle of 20° to the horizontal. From the time the first plumes are emitted, the flow enters a quasi-permanent regime during which the boundary layer remains, in average, marginally stable.

At a large scale, the boundary layer therefore has a constant thickness that leads to a constant dissolution rate. At a smaller scale, however, the convection flow induces a spatially variable dissolution rate: the boundary layer is thicker (and therefore the dissolution rate lower) vertically across a plume than between two plumes, where it is thinned by the arrival of fresh water rising from the bottom. These local variations in the dissolution rate then imprint on the soluble solid, all the more so as the emergence of a topography locks the position of the plumes [21]. Under a horizontal surface, this leads to the appearance of a pattern of cavities surrounded by
sharp crests reminiscent of scallops [111,112,145], the size of which reflects the wavelength of the hydrodynamic instability.

If the surface is inclined, the patterns obtained are qualitatively different. The dissolution flow tends to follow the inclination of the block as observed on the shadowgraph in Figure 9. This flow causes a spatial coherence of flow structures along the block and leads to the appearance of parallel stripes in the direction of gravity [21,133]. It can be seen in Figure 10, which shows photos at different time intervals of the underside of a block of plaster dissolving in water. The block of plaster was suspended in an aquarium of fresh water at room temperature, frequently replenished so that the concentration of solute remained far from saturation. Stripes appear, becoming deeper and wider, until they form grooves with very sharp edges. These same stripes can be observed on inclined caramel plates immersed in water, as shown in Figure 11. In this case, however, the initial pattern of stripes is itself unstable. The stripes rapidly cross to form chevrons and then asymmetrical concave, open at the bottom, scallop-like patterns, which run up the block in the opposite direction to the flow while widening [21,133]. The same succession of patterns can be seen on salt blocks. Figure 12 shows the (transverse) pattern, made up of cavities surrounded by sharp crests. The propagation of the patterns in the opposite direction to the flow can be explained by the fact that the plumes of solute-laden water detach at the crests (Figure 12). Thus, the dissolution rate is greater immediately downstream of a crest than upstream, causing it to propagate upstream [21]. Cohen et al. suggest that differential upstream velocities could explain transversalization of the initial stripe pattern because solute flow channelize on crests. Yet the precise mechanism explaining why this initial pattern destabilises in caramel and salt to form propagating scallops, which is not observed in plaster, remains to be elucidated. One of the major differences between the two cases is that the large differences in density between pure water and water saturated with salt or sugar ($\Delta\rho/\rho_0 = 0.2$ for salt and $0.33$ for sugar) lead to highly turbulent convection flows, with typical velocities of the order of a few cm/s, much larger than for plaster, for which $\Delta\rho/\rho_0 = 1.7 \times 10^{-3}$. For limestone, the density differences are even smaller ($\Delta\rho/\rho_0 < 3 \times 10^{-4}$), which makes it difficult to extrapolate these experimental results directly to the natural case. Although Cohen et al. point out the morphological similarity between the scallops appearing on blocks of salt or caramel and those adorning the walls of limestone caves [21,133], the possibility that solutal convection plays a role in the formation of the latter remains an open question, all the more so when an externally imposed flow is present. Underground rivers flowing in cave conduits can indeed reach typical velocities of the order of 1 m/s, much larger than the typical velocities induced by solutal convection. How such flows might influence the convection instability and the associated pattern formation is unclear. It is more likely, however, that the natural convection instability is at the origin of the large arrays of scallops observed on the immersed walls of icebergs, which melt in the polar oceans (see Figure 2(b)). Similar pattern formation is indeed expected when melting induces thermal or solutal convection. Cohen et al. [133] reported the formation of scallops at the top face of an ice block of fresh water melting in a bath of (densier) salt water. Between fresh water and seawater, the relative density difference $\Delta\rho/\rho_0$ equal to 0.025.

On a much larger scale, the average thickness of the boundary layer is fixed by the balance between diffusion and buoyancy (Equation (36)). For a flat surface inclined downwards, the dissolution rate is therefore constant over time and globally uniform along the surface. This was verified by Davies Wykes et al., who showed that sugar cones or prisms suspended upside down in an aquarium of water retain their shape [134]. The authors found no significant variation in erosion rate as a function of surface inclination angle, in contrast to Cohen et al. who found that erosion rate decreases with the inclination of the block. They argue that inclination should be taken into account in the gravity term, which predicts a dissolution rate proportional to $\cos(\theta)^{1/3}$, in agreement with their measurements [21]. This dependency of angle only becomes significant
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Patterning of the bottom surface of a plaster block immersed in tap water at an angle of 20° from the horizontal. The water was at room temperature, between 17 °C and 21 °C. The plaster block was regularly removed from the tank to take photographs of the surface and to record the topography. The scale is the same for each photograph. The block on the left is 20 cm long. Gravity points downward. *d* is the mean eroded thickness of the bottom surface. From left to right, photographs were taken after 148, 343, 541, and 928 h of immersion. Credits: Martin Chaigne.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Patterning of the bottom surface of a caramel block immersed in tap water (20 °C) at an angle of 60° from the horizontal (from Cohen et al. [21], with the permission of the American Physical Society).</figcaption></figure>

for large inclination, which may explain why Davies Wykes et al. did not observe it. We can therefore expect that an object of arbitrary shape dissolving from below will not always retain its shape. The large-scale convection can also shape cavities in confined geometries. By injecting water at a low flow rate into a cylindrical hole dug in salt, Gechter et al. observed that, due to the convective instability, the ceiling dissolved much faster than the sides and floor [146].
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left: Bottom surface of a salt block after immersion for 1 h in tap water (20 °C) at an angle of 60° from the horizontal. Right: Detail showing plumes of solute-laden fluid detaching from the crests of a caramel block. Both images are from Cohen et al. [21], with the permission of the American Physical Society.</figcaption></figure>

The initially circular profile deforms to become, over time, a triangular prism with a vertical ceiling. Sharma et al. obtained a similar shape by studying the formation of alcoves from defects on horizontal sugar walls dissolving in water [147].

## 4.2. Forced surface flows and associated bedforms

In the case of wind over sublimating surface or river current under a melting or dissolving surface, the velocity of the flow can be much larger than the velocity of the ablation process. In that case, two boundary conditions for convective flows have to be considered: the deep flow and the shallow water flow. Experiments should therefore be considered with care, ensuring that the wavelength of the patterns is either small in front of the flow thickness or much greater than it, respectively. To this end, the dimensionless number $k\delta_U$ has to be considered, where k is the wave number and $\delta_U$, the thickness of the fluid layer under consideration. This thickness corresponds to the height of the boundary layer thickness in atmospheric wind tunnels or in a closed water channel and to the height of the flow in a free surface flow. If $k\delta_U > 2\pi$, then deep flow condition should be considered, and if not, shallow water.

### 4.2.1. Deep flows

The patterning in response to a unidirectional water flow at high Reynolds number has been mainly investigated in deep flow regime for the dissolution of plaster plates, which have the same chemical composition than gypsum. By imposing a flow velocity of order the meter per second, several studies [18,19,70,148,149] report after about a hundred hours the formation of a scallop pattern, similar to the polygonal depressions observed in the nature and previously described in Section 2.2, Figure 2(c–d). In experiments, the typical scallop size and the surface roughness increase in time. In particular Allen [18,71] studies the temporal evolution of the pattern and the emergence of the scallops and pointed out the crucial role of initial defects, which are unavoidable with plaster due to its slight porosity. Defects larger than a critical size (about the centimeter) depending on the hydrodynamics give rise to “heel-shaped hollows of parabolic plan called flutes”. The subsequent pattern evolution generates an assemblage of scallops (Figure 13(c)). Smaller defects in contrast, grow mainly in length to produce narrow grooves parallel to the flow and of width about the centimeters. These patterns correspond to
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Deep flow conditions experiments. (a) Dissolution waves from Blumberg and Curl [20], copyright Cambridge University Press 1974 and (b) sublimation waves from Bordiec et al. [32], with permission from Elsevier. (c) Scallops from Allen [18] by dissolution of plaster. Copyright Cambridge University Press 1971. (d) Schematic of the parameters used in the scaling laws: wavelength $\lambda$, amplitude $\Delta h$, propagation angle $\phi$ and ablation rates $q_0$, adapted from Blumberg and Curl [20], copyright Cambridge University Press 1974.</figcaption></figure>

the *crested patterns parallel to the flow direction* described in Section 2.4, but are surprisingly generated rather in deep flow regime, whereas the natural examples correspond mostly to thin film flows. Allen proposed that these dissolution structures may be related to structures of the wall turbulence, the streaks. However, these grooves can generate for longer flow exposure in a second step, flutes and finally a scallop assemblage.

Regularly spaced crests, parallel to each other and oriented perpendicular to the turbulent overlaying flow (Section 2.3, Figure 3) are known to be formed under deep flow conditions. We will show in the remainder, that they have been reproduced in dissolution, sublimation and melting experiments. The mass transfer at work, coupled with the flow, induces a modulation of the surface on which transverse linear waves appear. In Blumberg and Curl [20] experiment, dissolution waves are reproduced on a plaster block as an analogous material of soluble rock (Figure 13b). Their morphodynamic study revealed wavelengths, amplitudes, propagation angles and dissolution rates defined in Figure 13(d) as a function of flow velocity and the ratio between fluid viscosity and diffusion coefficient, a dimensionless quantity known as the Schmidt number (see Table 2). By comparing the flow velocity with the ablation rate (typically $1 \text{ m} \cdot \text{s}^{-1}$ and $0.5 \text{ }\mu\text{mm} \cdot \text{s}^{-1}$), these experiments are well placed in a quasi-static interface evolution, i.e. $r_{d,U} \ll 1$ (see Table 2). The boundary layer of this flow $\delta_U \approx 7.6 \text{ cm}$ is imposed by the water tunnel geometry. This height must be greater than the wavelength of the patterns. Depending on the set of parameters, that we can reduce here to the dependence on two dimensionless numbers (Schmidt and Reynolds), dissolution waves have wavelengths between 2–5 cm. They have a millimetric amplitude and migrate at different speeds ($1–3.2 \times 10^{-7} \text{ m} \cdot \text{s}^{-1}$) in the direction of flow.

The aim of this approach is to understand the process of dissolution wave formation in a qualitative way, and help to enrich a database. It highlights a regime of appearance of this undulating surface roughness formed by dissolution in a transition regime, by considering a dimensionless number, the Reynolds number (see Table 2), based on the wavelength over the viscous length (i.e. the kinematic viscosity of the fluid $v$ over the friction velocity $u_*$). The friction velocity used for the velocity scale [150] is based on the wall shear stress such as $u_* = (\tau_w/\rho_f)^{1/2}$ with $\rho_f$, the fluid density. However, this approach does not allow the exact determination of the growth time of the patterns, which in this case do not emerge from an initially perfectly flat surface. In order to promote the appearance of patterns, the authors dug grooves transverse to
the flow at regular intervals to disrupt the flow. Nevertheless, according to the authors, these initial disturbances to the topography did not appear to affect either the wavelengths or the final profiles, as demonstrated by the tests they were able to carry out by varying the distance between intervals for the same set of parameters. The patterns obtained evolved gradually by increasing or decreasing the wavelength initially imposed by the grooves, in response to hydrodynamic instability.

A correlation has been established between dissolution waves (Figure 13(b)) and sublimation waves observed in the field (Figure 3(a)) and the analog sublimation experiment on CO₂ ice performed in an atmospheric wind tunnel [32] (Figure 13(c)). This theoretical study [11] is based on a 2D linear stability analysis of the coupled system between turbulent flow and mass transfer, following on from pioneer work [22,32,151]. It highlights the unstable wavelength range for solid bedforms by varying the Schmidt number, showing the dependence on it of the parameters involved in the dispersion relations (wavenumber λ, growth rate σ, angular frequency ω). The growth rate and angular frequency are functions of the wavelength-based Reynolds number ($Re_λ = λu_*/2πν = (kv/u_*)^{-1} = k_+^{-1}$), and a positive growth rate is found in the transitional near-wall regime, whereas modes in the turbulent or laminar regime decay due to their negative growth rate in the deep flow condition. The critical mode is the mode with the strongest positive growth rate, meaning that with time it will become the dominant wavelength visible on the ice surface. The dimensionless values of this critical growth rate, wave number and angular frequency are found numerically, and can be dimensionalised to give rise to three scaling laws:

$$ \lambda_c = \frac{\nu}{u_*} \frac{2\pi}{k_{c,ref}^+} Sc^{1/5} \quad (37) $$

$$ v_c = \frac{\omega_{c,ref}^+ q_0}{k_{c,ref}^+ \rho_s} \quad (38) $$

$$ t_c = \frac{1}{\sigma_{c,ref}^+} \frac{\rho_s}{q_0 u_*} \frac{\nu}{Sc^{1/6}} \quad (39) $$

with $k_{c,ref}^+ = 6.7 \times 10^{-3}$, $\omega_{c,ref}^+ = 0.013 \times 10^{-3}$, $\sigma_{c,ref}^+ = 8.9 \times 10^{-3}$ where reference values are given for $Sc = 1$ [11]. The first scaling law links the friction velocity to the critical wavelength that correspond to the most amplified mode. The second scaling law links the migration velocity to the ablation rate. The third scaling law combines the ablation rate and friction velocity estimated before into formation timescale, allowing to constrain the age of the surface.

Similar work could be carried out to extend the database in order to include melting bedforms [30,31,152,153]. Experimental works showed that under turbulent flow conditions, ice-water ripple patterns developed in response to an initial perturbation in the ice interface. While the topography was initially disturbed in previous experiments from Gilpin et al. [153], recent experimental work from Bushuk et al. [154] shows that the waves can develop from an initially flat bed with millimeter-scale imperfections but observes the formation of these patterns in the same flow regime known as laminar-turbulent transition based on $Re_λ$. These two-dimensional ice waves, which develop perpendicularly to the direction of flow, migrate in the direction of flow, indicating that the maximum mass flux is not found in the troughs, but slightly upstream from the creation of these patterns (see Figure 13(a)). These patterns can then evolve into three-dimensional structures, known as ice-scallops. They identify 3 stages of the evolution, starting from the flat bed to the scalloped interface. Recirculating 3D effects could induce reinforcing feedback due to turbulence redistribution for the three component of the kinetic energy due to wall effect. They identify a length of the shear production zone that could
constrain the wavelength of the scalloped bedforms. Finally, we note that the understanding of the evolution of an initial flat bed by melting or dissolution remains an open question. The conditions of the flow and of the bed, which will give rise either to wavy bedforms perpendicular to the flow or to straight patterns parallel to the flow, are not known. The initial defects or bed preparation seem to play a crucial role. Yet, further experimental and theoretical studies seem necessary to elucidate the question.

### 4.2.2. Shallow flow

In shallow water, the flow may strongly depend on the thickness of the flow and may couple with channel shape or free surface deformation, which can promote the growth of specific patterns. The case of melting/dissolution in presence of shallow water has been mainly investigated experimentally in creeping flow regime in the context of fracture dissolution and in the context of wormhole formation in a soluble porous medium. We refer to the recent review of Ladd and Szymczak [113] for more information on the subject. In transport controlled regime, the fracture aperture in gypsum decreases with the distance [155]. However, especially in experiments with NaCl the general shape and small scale pattern formation occurs in presence of buoyancy instabilities [156–158], which have been already discussed in Section 4.1. For faster flows, the shallow water situation corresponds often to the case of thin films flows, which occur in the nature as the run-off caused by the rainfalls. In this case the water free-surface dynamics may affect the pattern emergence and its evolution with time.

To our knowledge, among the diversity of *Karren* patterns (see Section 2.4) only the generation of *Rillenkarren* has been reported in quantitative experiments. Around 1980, Glew and Ford [15] applied artificial rain with an intensity about 35 mm/h upon molded plaster blocks of various inclination, leading to film flows about 100 µm depth. After several hundred hours an array of parabolic centimeters wide grooves were obtained. Their cross sections are well approximated by parabola. Sharp ridges separate the rills. The grooves disappear at the bottom of the blocks, when the film thickness exceeds approximately 150 µm. The obtained morphology resembles thus well to natural *Rillenkarren* reported on limestone. The authors estimate that an important condition to generate these pattern consists in the impact of the rain droplets on the water film, which may generate small scale turbulence and efficient solute mixing in the film. With a similar protocol, some morphological difference on the rain flutes development are found as a function of the water temperature (24 °C and 45 °C) [16]. More recently, Guerin et al. [17] demonstrate that dissolution grooves appear without drop impacts, by subjecting inclined soluble plates to runoff water flows. For flowing film of homogeneous thickness of order 200 µm, nearly parallel rills of similar width appear after few minutes on pink salt plates and after few hours on plaster plates (see Figure 14). Their width and depth increase linearly with time by pattern coarsening, until their crests emerge and channelize the flow, leading a slower increase of channel width. The initial typical channel width (of order 0.5 mm) appears to be proportional to the water film thickness, which is set by the hydrodynamic conditions. By rescaling the measured channel width and depth by the film thickness, the experimental points obtained for different flow parameters values and during all the experiment duration (few 10 h for plaster) gather on one master curve. Two formation scenarios have been proposed, a bed instability mechanism and an imprint of flow structures. In the first case, a locally greater fluid depth corresponds to a higher flow velocity, which increases the local dissolution rate. Any perturbation is thus amplified, leading to a linear instability of the bed. In the second scenario, a streamwise modulation of the flow velocity would produce a modulation of the dissolution rate, which could explain the linear growth of pattern amplitude with time. Recent theoretical studies [25,26] investigate more thoroughly the first scenario. The spontaneous emergence of rills is confirmed. The initial pattern wavelength selection arises from a balance between the film flattening at large scale by gravity and capillarity
<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) Side view of the experimental setup of Guerin et al. [17] with a block of pink Himalayan salt. (b) Front view. A homogeneous film of streaming water flows on a plaster block. (c) Evolution of the surface of a plaster block subjected to a runoff water film in the experiment. Flow rate 2.8 L/min, inclination 39°, average velocity 0.84 m/s. A pattern of dissolution grooves develops progressively in the direction of flow. Reproduced by permission of the American Physical Society.</figcaption></figure>

and the viscosity that acts as a low-pass space filter for small scale perturbations. However, the comparison with experiments and the field remains challenging, because once the pattern becomes noticeable with naked eye, its amplitude is likely large enough to perturb the flow hydrodynamics of the thin film. The regime of validity of the linear instability analysis would be limited to short times. Finally, we note that in the field, the age of the *Karren* are not precisely known, as well as the meteorological conditions that accompanied their formation. Further experimental and numerical studies are required to identify the mechanisms explaining the emergence of *Rillenkarren* and their long term evolution.

Although less investigated in the literature, dissolution patterns emerging from runoff flows have also some counterparts generated by melting. In particular, at the surface of glaciers, the meltwater is channelized into liquid water streams [159]. The flow in these channels may increase the heat exchange and thus their depths by enhancing melting. Mantelli et al. [160] proposed an
<figure><img src="image_17.png" /><figcaption>**Figure 15.** *Penitentes* in a laboratory scale experiment (from Bergeron et al. [161], reproduced by permission of the American Physical Society).</figcaption></figure>

explanation of the regular spacing between the ice streams by a linear stability analysis coupling the bed topography, the water flow and the heat transport. To our knowledge, no laboratory experiments have been reported in the literature. Finally, we note that the ice streams are an equivalent of the *Rinnenkarren* for dissolution. An analog of *Rillenkarren* is missing. Arrays of adjacent parallel grooves on steep slopes do not seem indeed commonly observed. The reason of this statement deserves also further studies.

### 4.3. *Pattern formation stimulated by radiation*

The heat required for melting or sublimation can come largely from the absorption of light, which is often the case in glaciers under sublimation conditions or where the ambient temperature is low. Then, modulation in albedo, or exposure to light, plays an important role in the formation of patterns such as penitentes (Figure 1(a)), glacier tables or dirt cones (Figure 6).

Penitentes have been reproduced in a laboratory scale experiment by Bergeron et al. [161]. A block of snow was kept in a chamber with a plexiglas cover at low temperature in a dried air environment (relative humidity was typically between 5 and 10%) and illuminated from above with a 120 W lamp, disposed outside the chamber. In sublimation conditions, for temperature between -10 and -20 °C, penitentes appeared after two hours of irradiation with a characteristic of a few millimeters (Figure 15). In melting conditions, above -4 °C, no structures were observed. In sublimation condition, under -35 °C, no penitentes formed after irradiation time as long as ten hours. Authors suggest that the block may be saturated in water vapor because vapor pressure of water rapidly decreases with temperature below 0 °C. They also report that a breeze of 2.5 m/s above the snow prevents penitentes formation. Authors propose that initial pattern are the smallest possible ablation structures, which are limited by the optical extinction length (characteristic length on which light is absorbed in the snow) [42,161]. The optical extinction length depends on particle radius and ice volume fraction in the snow, and is estimated to be of the order of 1 mm in their experiment [161].

Pattern formation is generally ascribed to a geometric coupling between light absorption and topography: light absorption is greater at troughs than at peaks [41,42,162,163]. Claudin et al. did a linear stability analysis and found two driving mechanisms: the above mentioned geometric effect but also a Mullins-Sekerka type mechanism [46,164]. Since the light is absorbed by the snow (or ice) on a given distance from the interface, and the temperature of the interface is held at
the (cold) ambient air temperature, the temperature increases with depth in the snow. This drives a heat flux in the snow toward the interface, which is greater at troughs than at peaks. In their calculation, the balancing stabilizing effect comes from vapor diffusion, which authors arbitrarily limit to a boundary layer of finite thickness, thereby providing the characteristic wavelength of the instability and pattern formation [46]. Authors suggest that this characteristic thickness could be set by wind, forcing turbulent mixing above a certain height above the snow. However, this seems to be challenged by the observation of Bergeron et al. [161]. Cohen et al. proposed that this length scale could be set by a solutal Rayleigh-Bénard instability, because density of water vapor is less than density of air [21]. This instability in the buoyancy-driven vapor flow would prevail over other mechanisms in initiating the pattern and selecting the wavelength of formation.

Patterns are not exclusively shaped by instabilities. Recognizable structures can also be formed by differential ablation due to specific boundary conditions, such as the interaction with a non-melting (or less eroding) cover. This is the case for fairy chimneys (also named hoodoo), glacier tables, dirt cones, and zen stone, for example (Figure 6). Here we review recent laboratory experimental works that have addressed the formation of glacier tables, dirt cones, and zen stone described in Section 2.6, Figure 6.

Hénot et al. have studied the formation of glacier tables and, in particular, the conditions for their formation [84]. They followed the melting of an ice sheet partially covered by a cylindrical cap made of different materials and with different height-to-radius aspect ratios. The ice melts at room temperature (25 °C) due to infrared radiation from the walls (at room temperature) and convection in the air (because cold air at the ice interface is denser than the air at room temperature). They found that, depending on the cap material, the cap may sink into the ice sheet or a glacier table may form, which corresponds to the ice under the cap melting faster or slower than the surrounding sheet, respectively (Figure 16). They explained the onset of table formation by the competition between two effects. On the one hand, the cap insulates the underlying ice because the cap temperature is higher than the ice temperature, which decreases the infrared radiative influx (and increases the infrared radiative outflux) and the convective flux. On the other hand, the cap exposes a larger area to the incoming heat than the covered (contact) area, which may increase the (conductive) heat flux to the covered ice. The relative weight of these two effects depends on the thermal conductivity of the cap the material. The conductive flux (from the cap to the ice) decreases and the cap temperature (and insulation) increases as the thermal conductivity decreases, which is favorable for table formation. This modeling proved to be predictive in the field [85].

Hénot et al. used the same approach to study the formation of dirt cones: they followed the melting at room temperature of an ice sheet partially covered with a pile of plastic grains [86]. Differential melting resulted in the formation of an ice cone covered by a granular layer. The authors explain the formation of dirt cones (Section 2.6, Figure 6), not directly by the dependence of the melting rate on the debris thickness [89], but by the same mechanisms as for glacier tables with the additional coupling to the granular medium: as the cone forms, the exposed capped area increases and the shape can reach a steady state [86,165].

Taberlet and Plihon have studied the formation of zen stones by ice (differential) sublimation in a lyophilizer [166]. The initial flat ice sheet is partially covered with a 3 cm-diameter 5 mm-thick disk. Using various metals for the covering disk, they showed that the material conductivity is of little importance in the formation of zen stone. They observed the formation of a zen stone for aluminum, grooved aluminium, and copper disks. They conclude that the main mechanism is an umbrella effect, the covering disk shields the ice from external radiations [166].
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Ice underlying cylindrical caps melts faster or slower than the surrounding ice sheet depending on the cap material (from Hénot et al. [84], reproduced by permission of the American Physical Society).</figcaption></figure>

**5. Discussion and conclusions**

5.1. *Solid bedforms as imprint of flow structures*

The objectives of this review paper were: (i) to examine the natural bedforms induced by phase
change that occur at the solid surface of Earth and of planets and (ii) to demonstrate how
laboratory experiments can contribute to a better understanding of the pattern formation and
their evolution. These bedforms are produced largely by the combined action of a flow and a
phase change, which may be melting, dissolution, or sublimation. Although the phase changes
at work are not the same, the natural objects inventoried in this study have similar shapes,
whether the substrate is rock or ice. We proposed a classification of patterns in the field, based
on morphologic indicators, into 6 categories: sharp patterns, polygonal depression, crested
bedforms perpendicular or parallel to the flow, stepped patterns, and structures due to capping.
Then, we presented the theoretical framework for modeling the interaction between a flow and
a phase change, with the objective of elucidating the evolution of the solid/fluid interface and
the emergence of patterns. Subsequently, a review of the previous experimental studies in
the literature reporting bedforms was conducted, and a classification of experiments has been
established based on the type of flow (flow generated by the phase change or imposed externally)
or boundary conditions.

To be able to relate the classification of these shapes to that of the type of flow (forced and
natural convection) and boundary conditions, it is necessary to consider whether some of these
shapes form part of a continuum or not. For example, by experimentally reproducing in a
controlled environment ice scallops in a water-tunnel under unidirectional turbulent flow [153],
2D structures similar to the sublimation ripples observed in blue ice-areas in Antarctica, were
initially obtained, some of which subsequently evolved into 3D structures. There could therefore
be a continuum between these forms, which would explain why they all line up with the curve
proposed by Thomas [6]. In addition, these scallops exist in both forced and natural convection
flows, and one can ask if there is a morphometric parameter that could provide insight into
the competition between inertia and buoyancy? The dynamic nature of these objects, which
propagate downstream (forced convection) or upstream (natural convection), should help to
refine the processes involved and better understand the environments in which they were
formed.

5.2. *Effect of non-linearities on the resulting bedform*

Among the solid bedforms in ablation, some of them are generated by instabilities of the flow or
of the bed that can initiate a characteristic scale related to the parameters of the flow. However,
as the pattern amplitude grows several effects must be taken into account to explain its shape
at a given time, and these often cannot be analyzed linearly. Firstly, geometric coarsening
occurs when the small structures are absorbed by the large ones. By this merging process, the
characteristic wavelength increases with time. A measured pattern size must be put in relation
with the pattern age. This effect has been particularly evidenced in laboratory experiments [17]
and interface evolution models [111]. In addition, other geometric focusing effects can lead
to crest formations by normal ablation. This phenomenon initially proposed by A. Lange in
the fifties to explain the scallop shapes in limestone wall [167] has been recently revisited by
Chaigne et al. [111]. They have shown that any ablation process on an interface displaying sign
changes of curvatures will show in finite time emergence of singularities, i.e. discontinuities of
the gradient. This effect can explain, why the sinusoidal patterns predicted by linear instability
mechanism evolve generically to structures with sharp tips and crests. Although the singularities
are regularized at small scale by material properties or the boundary layer hydrodynamics, this
concept describes well the qualitative shape of a solid surface undergoing an ablation process
like melting, dissolution or sublimation.

Once a pattern has reached a sufficient size, it can retroact on the flow that created it. In
the case of the patterns emerging from a bed instability like for the mechanisms proposed for
scallops [22] in turbulent flows or for Rillenkarren by run-off [17,25,26], this mechanism is a
part of the instability mechanism. The pattern amplitude increases the flow modulation, which
amplifies in turn the pattern growth. However, the saturation of the instability relies often on
nonlinear mechanisms, which may involve higher order terms in amplitude equations [168].
Physically, a nonlinear feedback can decrease the flow modulation, once the pattern amplitude
exceeds a characteristic scale, typically the boundary layer thickness of the relevant field
(temperature, solute concentration, vapor concentration), the free-surface in thin films or the
system size. For patterns emerging by an instability of the flow like those created by from
Rayleigh-Bénard convection [21,112,133], the hydrodynamic feedback is more complex. In linear
regime, the surface is a fingerprint of the flow structures. However, we noted the need of plumes
position locking to explain the emergence of coherent dissolution [21] or melting patterns [169].
These effects are currently not well addressed quantitatively and deserve further study.

Hydrodynamic focusing constitutes a more specific situation of feedback of the pattern
topography on the flow. For example, in a free-surface flow, once separated channels have been
formed, they will collect more water, until a drainage network emerges. This phenomenon is
more described for rivers carved by erosion, but analogous cases occur by dissolution like the
*Rinnenkarren* [55] on limestone pavements or by melting like the ice streams on Glacier [160].

The nonlinear study of ablation pattern in nature could constitute the next step of future
researches.
### 5.3. On the use of scaling laws for planetary applications

Because of the large timescale required to form natural solid bedforms, the choice of analogous materials that enable to sustain experimental timescale is a real challenge. In this context, experiments on the dissolution of plaster, caramel, sugar or salt have provided an opportunity to explore a wide range of natural bedforms under different flow conditions (deep flow, free surface flow, buoyancy driven flow with or without stable stratification). It is then possible to reproduce most of the natural shapes identified, such as sharp-edged, crested bedforms parallel or perpendicular to the flow, and scalloped. To our knowledge, stepped bedforms have not yet been reproduced in dissolution, to date just in melting [24]. For cyclic steps (Section 2.5, Figure 5), therefore, modelling and comparing models with laboratory experiments has only been tackled to a limited extent. Other experimental studies in melting concern crested bedforms parallel to the flow and scallops [153,154]. Experiments in sublimation are difficult to implement, requiring highly constrained environments (temperature, humidity, hygrometry). As a result, they only exist for penitentes and crested bedforms perpendicular to the flow, without any systematic exploration of the parameters involved.

However, sublimation dominates on many icy surfaces of planetary bodies. Therefore, complementary phase change studies would provide additional databases that could be useful for validating scaling laws that are commonly used in planetology. Further work is needed both on modeling and on comparative analogue studies between different phase changes to confirm some of these scaling laws. A first step in this direction has been initiated by Carpy et al. [11], in the case of crested waves perpendicular to the flow direction, considering dissolution and sublimation. They use sublimation and dissolution database to validate scaling laws obtained theoretically on the basis of a linear stability analysis of the coupled system of equations between mass transfer and the turbulent boundary layer in a deep flow case. Further development is required to include melting database. These scaling laws link the friction velocity to the critical wavelength and the critical migration velocity to the sublimation rate. The third scaling law that combines the sublimation rate and friction velocity estimated before into a sublimation wave formation timescale, allowing to constrain the age of the ice surface has yet to be validated. Before applying scaling laws to planetary surfaces, it is necessary to bear in mind the constraints on the use of these laws, independently of the non-linear effects discussed in the previous paragraph (5.2). In the case of penitentes (Section 2.1, Figure 1(a)), new experiments (and models) are needed to understand the role of fluid flow in their formation in order to ascertain if buoyancy-driven vapor flow plays a role in their formation as suggested by Cohen et al. [21].

### 5.4. Concluding remarks

Although geomorphology and planetology have been the main motivation for the experimental studies reported here, they raise fundamental questions that combine several fields, including hydrodynamics, thermodynamics, materials science, nonlinear physics, and morphogenesis. Laboratory-scale experiments are key to understanding the role of flows in the formation processes and dynamics of these solid bedforms undergoing ablation phase change. Since experiments take place over shorter timescales, forcings and boundary conditions can be tested. As we have seen from the various experiments collected, these solid bedforms correspond to imprints of the flow, driven either by boundary conditions, by instabilities intrinsic to the flow, or by instabilities associated with a perturbation of the topography at the interface. This is an important distinction. However, the non-linear pattern dynamics and geometric effects result in a time evolution of the pattern, with a typical slow increase in characteristic length scales through coarsening and the occurrence of sharp structures. In some specific cases, in particular when the
pattern emergence can be explained by a linear instability mechanism of the interface, scaling laws can be derived relating the typical bedform length to hydrodynamic parameters.

Comparing phase changes with each other provides a better understanding of the processes occurring at the interface, and allows to investigate the large dimensionless numbers and different scales involved. Due to the complex nature of the flow or the complicated boundary conditions, as well as the geometric shapes of natural objects, numerical simulations are limited, so experiments play a key role. Many natural pattern shapes have not yet been addressed, which could lead to new experimental projects. In addition, further experimental studies are required to address some pending questions, such as the evolution of a flat, dissolving or melting bed in the presence of a turbulent flow. Similar work could be carried out on accumulation patterns (solidification, precipitation, condensation). In this case, will the change in the sign of the transfer lead to a "photo negative", or will this have an effect on the shape, size and diversity of the patterns obtained? This question opens up new experimental perspectives to be explored. Improving models through the laboratory experiments could help predict landscape evolution or infer past conditions from patterns.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Funding

SC is supported by French Research Council project through grant SHERPAS ANR-23-CE49-0006. MB, MC and SCdP are supported by French Research Council project through grant PhysErosion ANR-22-CE30-0017. We thank M. Bordiec and O. Bourgeois for scientific discussion on natural examples.

## References

[1] F. Pascheke, J. F. Barlow and A. Robins, "Wind-tunnel modelling of dispersion from a scalar area source in urban-like roughness", *Bound.-Layer Meteorol.* **126** (2008), no. 23, pp. 103–124.

[2] C. Combs, N. T. Clemens and P. M. Danehy, "Development of naphthalene PLIF for visualizing ablation products from a space capsule heat shield", in *52nd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics: Reston, 2014. AIAA 2014-1152*. Online at https://arc.aiaa.org/doi/10.2514/6.2014-1152 (accessed January 31, 2025).

[3] G. Hochrein and G. Wright Jr., "Analysis of the TATER nosetip boundary layer transition and ablation experiment", in *14th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics: Reston, 1976. AIAA 1976-167*. Online at https://arc.aiaa.org/doi/abs/10.2514/6.1976-167 (accessed January 31, 2025).

[4] S. P. Schneider, "Laminar-turbulent transition on reentry capsules and planetary probes", *J. Spacecr. Rockets* **43** (2006), no. 6, pp. 1153–1173.

[5] T. C. Lin and P. Qun, "On the formation of regmaglypts on meteorites", *Fluid Dyn. Res.* **1** (1987), no. 3-4, pp. 191–199.

[6] R. M. Thomas, "Size of scallops and ripples formed by flowing water", *Nature* **277** (1979), pp. 281–283.

[7] I. B. Smith, J. W. Holt, A. Spiga, A. D. Howard and G. Parker, "The Spiral troughs of Mars as cyclic steps", *J. Geophys. Res.: Planets* **118** (2013), no. 9, pp. 1835–1857.

[8] J. E. Moores, C. L. Smith, A. D. Toigo and S. D. Guzewich, "Penitentes as the origin of the bladed terrain of Tartarus Dorsa on Pluto", *Nature* **541** (2017), no. 7636, pp. 188–190.

[9] N. Mangold, "Ice sublimation as a geomorphic process: A planetary perspective", *Geomorphology* **126** (2011), no. 1-2, pp. 1–17.

[10] A. Ohmura, "Cryosphere during the twentieth century", *Geophys. Monogr. Ser.* **150** (2004), pp. 239–257.
[11] S. Carpy, M. Bordiec and O. Bourgeois, "Scaling laws for ablation waves formed by ice sublimation and rock dissolution: applications to the Earth, Mars and Pluto", *Front. Astron. Space Sci.* **10** (2023), article no. 1176158.

[12] D. Ford and P. D. Williams, *Karst Hydrogeology and Geomorphology*, John Wiley & Sons, Springer: Netherlands, 2007.

[13] M. T. Kozary, J. C. Dunlap and W. E. Humphrey, "Incidence of saline deposits in geologic time", in *Saline Deposits: A Symposium based on Papers from the International Conference on Saline Deposits, Houston, Texas, 1962* (R. B. Mattox, W. T. Holser, H. Ode, W. L. McIntire, N. M. Short, R. E. Taylor and D. C. Van Siclen, eds.), Geological Society of America: New York, 1968, pp. 43–57.

[14] F. Gallaire and P. T. Brun, "Fluid dynamic instabilities: Theory and application to pattern forming in complex media", *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* **375** (2017), no. 2093, pp. 1–44.

[15] J. R. Glew and D. C. Ford, "A simulation study of the development of Rillenkarren", *Earth Surf. Process.* **5** (1980), no. 1, pp. 25–36.

[16] T. Slabe, A. Hada and M. Knez, "Laboratory modeling of karst phenomena and their rock relief on plaster, subsoil karren, rain flutes karren and caves", *Acta Carsologica* **45** (2016), no. 2, pp. 187–204.

[17] A. Guérin, J. Derr, S. Courrech du Pont and M. Berhanu, "Streamwise dissolution patterns created by a flowing water film", *Phys. Rev. Lett.* **125** (2020), no. 19, article no. 194502.

[18] J. R. L. Allen, "Bed forms due to mass transfer in turbulent flows: a kaleidoscope of phenomena", *J. Fluid Mech.* **49** (1971), no. 1, pp. 49–63.

[19] M. F. Goodchild and D. C. Ford, "Analysis of scallop patterns by simulation under controlled conditions", *J. Geol.* **80** (1972), no. 1, pp. 121–122.

[20] P. N. Blumberg and R. L. Curl, "Experimental and theoretical studies of dissolution roughnesses", *J. Fluid Mech.* **65** (1974), no. 4, pp. 735–751.

[21] C. Cohen, M. Berhanu, J. Derr and S. Courrech du Pont, "Buoyancy driven dissolution of inclined blocks: Erosion rate and pattern formation", *Phys. Rev. Fluids* **5** (2020), article no. 053802.

[22] P. Claudin, O. Durán and B. Andreotti, "Dissolution instability and roughening transition", *J. Fluid Mech.* **832** (2017), no. R2, pp. 1–14.

[23] C. Camporeale and L. Ridolfi, "Ice ripple formation at large Reynolds numbers", *J. Fluid Mech.* **694** (2012), pp. 225–251.

[24] M. Yokokawa, N. Izumi, K. Naito, G. Parker, T. Yamada and R. Greve, "Cyclic steps on ice", *J. Geophys. Res: F. Earth Surf.* **121** (2016), no. 5, pp. 1023–1048.

[25] M. B. Bertagni and C. Camporeale, "The hydrodynamic genesis of linear karren patterns", *J. Fluid Mech.* **913** (2021), article no. A34.

[26] S. Djambov and F. Gallaire, "Boundary-layer approach to the linear karren instability", *J. Fluid Mech.* **977** (2023), article no. R3.

[27] C. B. Thorsness and T. J. Hanratty, "Stability of dissolving or depositing surfaces", *AIChE J.* **25** (1979), no. 4, pp. 697–701.

[28] P. Y. Chan and N. Goldenfeld, "Steady states and linear stability analysis of precipitation pattern formation at geothermal hot springs", *Phys. Rev. E* **76** (2007), no. 4, pp. 1–11.

[29] C. Camporeale and L. Ridolfi, "Hydrodynamic-driven stability analysis of morphological patterns on stalactites and implications for cave paleoflow reconstructions", *Phys. Rev. Lett.* **108** (2012), no. 23, pp. 1–5.

[30] G. D. Ashton and J. F. Kennedy, "Ripples on the underside of river ice covers", *J. Hydraul Eng.* **98** (1972), pp. 1603–1624.

[31] G. D. Ashton, "Turbulent heat transfer to wavy boundaries", in *Proceedings of the Heat Transfer and Fluid Mechanics Institute*, 1972, pp. 200–213.

[32] M. Bordiec, S. Carpy, O. Bourgeois, et al., "Sublimation waves: Geomorphic markers of interactions between icy planetary surfaces and winds", *Earth-Sci. Rev.* **211** (2020), article no. 103350.

[33] J. R. L. Allen (ed.), *Sedimentary Structures, Their Character and Physical Basis*, Developments in Sedimentology, Elsevier: Amsterdam, Oxford, New York, 1982, p. 592.

[34] G. C. Amstutz, "On the formation of snow penitentes", *J. Glaciol.* **3** (1958), no. 24, pp. 304–311.

[35] K. E. Sinclair and S. MacDonell, "Seasonal evolution of penitente glaciochemistry at Tapado Glacier, Northern Chile", *Hydrol. Process.* **30** (2016), no. 2, pp. 176–186.

[36] N. E. Odell, "Ablation at high altitudes and under high solar incidence", *Am. J. Sci.* **239** (1941), no. 5, pp. 379–382.

[37] J. M. Moore, A. D. Howard, O. M. Umurhan, et al., "Bladed Terrain on Pluto: Possible origins and evolution", *Icarus* **300** (2018), pp. 129–144.

[38] D. E. J. Hobley, J. M. Moore, A. D. Howard and O. M. Umurhan, "Formation of metre-scale bladed roughness on Europa's surface by ablation of ice", *Nat. Geosci.* **11** (2018), no. 12, pp. 901–904.

[39] K. P. Hand, D. Berisford, T. Daimaru, J. Foster, A. E. Hofmann and B. Furst, "Penitente formation is unlikely on Europa", *Nat. Geosci.* **13** (2020), no. 1, pp. 17–19.
[40] D. E. J. Hobley, J. M. Moore, A. D. Howard and O. M. Umurhan, "Reply to: Penitente formation is unlikely on Europa", *Nat. Geosci.* **13** (2020), no. 1, pp. 20–21.

[41] L. Lliboity, "The origin of penitentes", *J. Glaciol.* **2** (1954), no. 15, pp. 331–338.

[42] M. D. Betterton, "Theory of structure formation in snowfields motivated by penitentes, suncups, and dirt cones", *Phys. Rev. E* **63** (2001), no. 5, article no. 056129.

[43] S. Hastenrath and B. Koci, "Micro-morphology of the snow surface at the Quelccaya ice cap, Peru", *J. Glaciol.* **27** (1981), no. 97, pp. 423–428.

[44] F. E. Matthes, "Ablation of snow-fields at high altitudes by radiant solar heat", *Trans. Am. Geophys. Union* **15** (1934), article no. 380.

[45] M. D. Betterton, "Theory of structure formation in snowfields motivated by penitentes, suncups, and dirt cones", *Phys. Rev. E* **63** (2001), no. 5, pp. 1–12.

[46] P. Claudin, H. Jarry, G. Vignoles, M. Plapp and B. Andreotti, "Physical processes causing the formation of penitentes", *Phys. Rev. E* **92** (2015), no. 3, article no. 033015.

[47] M. Veress, Z. Zentai, K. Péntek, Z. Mitre, G. Deák and S. Samu, "Flow dynamics and shape of rinnenkarren systems", *Geomorphology* **198** (2013), pp. 115–127.

[48] J. M. Huang, J. Tong, M. Shelley and L. Ristroph, "Ultra-sharp pinnacles sculpted by natural convective dissolution", *Proc. Natl. Acad. Sci. USA* **117** (2020), no. 38, pp. 23339–23344.

[49] M. Bordiec, *Formation d'ondes topographiques par transfert de masse et diffusion dans un écoulement turbulent : analyse d'exemples naturels terrestres et planétaires et modélisations*, PhD thesis, Nantes Université, 2020.

[50] R. Bintanja, "On the glaciological, meteorological, and climatological significance of Antarctic Blue Ice areas", *Rev. Geophys.* **37** (1999), no. 3, pp. 337–359.

[51] K. L. Carey, "Observed configuration and computed roughness of the underside of rive ice, St. Croix River, Wisconsin", *US Geol. Survey Prof. Pap.* **550** (1966), B192–B198.

[52] J. N. Salomon, "The tsingy karrenfields of Madagascar", in *Karst Rock Features: Karren Sculpturing* (A. Ginés, M. Knez, T. Slabe and W. Dreybrodt, eds.), Karst Research Institute, Založba ZRC, 2009, pp. 411–422.

[53] M. Day and T. Waltham, "The pinnacle karrenfields of Mulu", in *Karst Rock Features. Karren Sculpturing* (A. Ginés, M. Knez, T. Slabe and W. Dreybrodt, eds.), Karst Research Institute, Založba ZRC, 2009, pp. 423–432.

[54] M. Knez and T. Slabe, "Lithologic and morphological properties and rock relief of the Lunan stone forests", in *Karst Rock Features. Karren Sculpturing* (A. Ginés, M. Knez, T. Slabe and W. Dreybrodt, eds.), Karst Research Institute, Založba ZRC, 2002, pp. 259–266.

[55] A. Ginés, "Karrenfield Landscapes and Karren Landforms", in *Karst Rock Features, Karren Sculpturing* (A. Ginés, M. Knez, T. Slabe and W. Dreybrodt, eds.), Karst Research Institute, Založba ZRC, 2009, pp. 13–24.

[56] A. Frumkin, "Salt karst", in *Treatise on Geomorphology*, Academic Press: San Diego, CA, 2013, pp. 407–424.

[57] J. J. Rhodes, R. L. Armstrong and S. G. Warren, "Mode of formation of "ablation hollows" controlled by dirt content of snow", *J. Glaciol.* **33** (1987), no. 114, pp. 135–139.

[58] M. C. Malin and K. S. Edgett, "Mars global surveyor Mars orbiter camera: Interplanetary cruise through primary mission", *J. Geophys. Res.* **106** (2001), no. E10, pp. 23429–23570.

[59] S. Milkovich, "Surface textures of Mars' north polar layered deposits: A framework for interpretation and future exploration", *Mars J.* **2** (2006), pp. 21–45.

[60] M. E. Banks, S. Byrne, K. Galla, et al., "Crater population and resurfacing of the Martian north polar layered deposits", *J. Geophys. Res. E: Planets* **115** (2010), no. 8, pp. 1–11.

[61] M. Massé, O. Bourgeois, S. Le Mouélic, C. Verpoorter, L. Le Deit and J. P. Bibring, "Martian polar and circum-polar sulfate-bearing deposits: Sublimation tills derived from the North Polar Cap", *Icarus* **209** (2010), no. 2, pp. 434–451.

[62] M. W. Telfer, E. J. R. Parteli, J. Radebaugh, et al., "Dunes on Pluto", *Science* **360** (2018), no. 6392, pp. 992–997.

[63] P. Bella, "Morphology of ice surfaces in the Dobšiná Ice Cave", in *Proceedings of the 2nd International Workshop on Ice Caves* (S. Turri and J. Zelinka, eds.), Slovak Caves Administration: Liptovský Mikulá, 2006, pp. 15–23. Online at http://www.ssj.sk/en/user_files/IWIC_2_15.pdf (accessed January 31, 2025).

[64] A. Persoiu and S.-E. Lauritzen (eds.), *Ice Caves*, Elsevier Inc.: Amsterdam, Oxford, Cambridge, 2018, p. 752.

[65] F. Obleitner and C. Spötl, "The mass and energy balance of ice within the Eisriesenwelt cave, Austria", *Cryosphere* **5** (2011), no. 1, pp. 245–257.

[66] A. Perşoiu, B. P. Onac and I. Perşoiu, "The interplay between air temperature and ice mass balance changes in Scărişoara Ice Cave, Romania", *Acta Carsologica* **40** (2011), no. 3, pp. 445–456.

[67] R. L. Curl, "Scallops and flutes", *Trans. Cave Res. Group of Great Britain* **7** (1966), no. 2, pp. 121–160.

[68] P. Marshall and M. C. Brown, "Ice in Coulthard Cave, Alberta", *Canad. J. Earth Sci.* **11** (1974), no. 4, pp. 510–518.

[69] P. Meakin and B. Jamtveit, "Geological pattern formation by growth and dissolution in aqueous systems", *Proc. R. Soc. A: Math. Phys. Eng. Sci.* **466** (2010), no. 2115, pp. 659–694.

[70] B.-P. Villien, Y. Zheng and D. H. Lister, "The scalloping phenomenon and its significance in flow-assisted corrosion", in *Proceedings of the Twenty Sixth Annual CNS-CNA Student Conference*, Canadian Nuclear Society:
Toronto, 2001, pp. 1–10. Online at https://www.nuceng.ca/univcomm/studconf2001/S16.PDF (accessed January 31, 2025).

[71] J. R. Allen, "Transverse erosional marks of mud and rock: their physical basis and geological significance", *Sediment. Geol.* **5** (1971), no. 3-4, pp. 167-385.

[72] O. Hammer, S. E. Lauritzen and B. Jamtveit, "Stability of dissolution flutes under turbulent flow", *J. Cave Karst Studies* **73** (2011), no. 3, pp. 181-186.

[73] I. Karcz and I. Zak, "Bedforms in salt deposits of the Dead sea brines", *J. Sediment. Res.* **57** (1987), no. 4, pp. 723-735.

[74] R. Bintanja, C. H. Reijmer and S. J. M. H. Hulscher, "Detailed observations of the rippled surface of Antarctic blue-ice areas", *J. Glaciol.* **47** (2001), no. 158, pp. 387-396.

[75] G. Weller, "The heat and mass balance of snow dunes on the central Antarctic plateau", *J. Glaciol.* **8** (1969), no. 53, pp. 277-284.

[76] M. Veress, "The KARREN and KARREN formation of bare slopes", *Earth-Sci. Rev.* **188** (2019), pp. 272-290.

[77] J. Lundberg and A. Ginés, "Rillenkarren", in *Karst Rock Features, Karren Sculpturing*, Karst Research Institute, Zalozba ZRC, 2009, pp. 185-210.

[78] B. Tourte, F. Chenu, L. Morel, et al., "Ruissellement et dissolution dun bassin versant lapiazé des karsts de Patagonie, île Tarlton, archipel de Madre de Dios (Chili)", *Collection EDYTEM. Cahiers de géographie* **19** (2017), no. 1, pp. 137-146. Included in a thematic issue: Monitoring en milieux naturels. Retours d'expériences en terrains difficiles.

[79] J. F. Kennedy, "The formation of sediment ripples, dunes, and antidunes", *Annu. Rev. Fluid Mech.* **1** (1969), pp. 147-168.

[80] N. Izumi, M. Yokokawa and G. Parker, "Incisional cyclic steps of permanent form in mixed bedrock-alluvial rivers", *J. Geophys. Res.: Earth Surf.* **122** (2017), no. 1, pp. 130-152.

[81] D. Mottershead, "Identification and mapping of rock weathering surface forms and features", *Zeitschrift für Geomorphologie Supplementband* **120** (2000), pp. 5-22.

[82] G. W. Evatt, C. Mayer, A. Mallinson, I. D. Abrahams, M. Heil and L. Nicholson, "The secret life of ice sails", *J. Glaciol.* **63** (2017), no. 242, pp. 1049-1062.

[83] A. C. Fowler and C. Mayer, "The formation of ice sails", *Geophys. Astrophys. Fluid Dyn.* **111** (2017), no. 6, pp. 411-428.

[84] M. Hénot, N. Plihon and N. Taberlet, "Onset of glacier tables", *Phys. Rev. Lett.* **127** (2021), no. 10, article no. 108501.

[85] M. Hénot, V. J. Langlois, J. Vessaire, N. Plihon and N. Taberlet, "Formation of glacier tables caused by differential ice melting: field observation and modelling", *Cryosphere* **16** (2022), no. 6, pp. 2617-2628.

[86] M. Hénot, V. J. Langlois, N. Plihon and N. Taberlet, "How dirt cones form on glaciers: Field observation, laboratory experiments, and modeling", *Phys. Rev. E* **107** (2023), no. 3, article no. 034905.

[87] G. Östrem, "Ice melting under a thin layer of moraine, and the existence of ice cores in moraine ridges", *Geografiska Annaler* **41** (1959), no. 4, pp. 228-230.

[88] L. E. Mattson, "Ablation on debris covered glaciers: an example from the Rakhiot Glacier, Punjab, Himalaya", *Snow Glacier Hydrol.* **218** (1993), pp. 289-296.

[89] G. W. Evatt, I. D. Abrahams, M. Heil, C. Mayer, J. Kingslake, S. L. Mitchell, A. C. Fowler and C. D. Clark, "Glacial melt under a porous debris layer", *J. Glaciol.* **61** (2015), no. 229, pp. 825-836.

[90] M. G. Worster, "Solidification of fluids", in *Perspectives in Fluid Dynamics - A Collective Introduction to Current Research* (G. K. Batchelor, H. K. Moffatt and M. G. Worster, eds.), Cambridge University Press: Cambridge, 2000, pp. 393-446. Online at http://www.itg.cam.ac.uk/people/grae/44.pdf (accessed January 31, 2025).

[91] S. L. Brantley, in *Kinetics of Water-Rock Interaction* (S. L. Brantley, J. D. Kubicki and A. F. White, eds.), Springer: New York, 2008.

[92] V. Alexiade and A. D. Solom, *Mathematical Modeling of Melting and Freezing Processes*, Hemisphere Publishing Corporation, Taylor & Francis: Washington DC, 1993.

[93] H. B. Callen, *Thermodynamics and An Introduction to Thermostatistics*, John Wiley & Sons: New York, Chichester, Brisbane, Toronto, Singapore, 1985.

[94] H. B. Aaron and G. R. Kotler, "The effects of curvature on the dissolution kinetics of spherical precipitates", *Metal Sci. J.* **4** (1970), no. 1, pp. 222-225.

[95] J. Colombani and J. Bert, "Holographic interferometry study of the dissolution and diffusion of gypsum in water", *Geochim. Cosmochim. Acta* **71** (2007), no. 8, pp. 1913-1920.

[96] I. W. Eames, N. J. Marr and H. Sabir, "The evaporation coefficient of water: a review", *Int. J. Heat Mass Transfer* **40** (1997), no. 12, pp. 2963-2973.

[97] C. Mas, *Modélisation physique du procédé de découpe de métaux par laser*, PhD thesis, 2003. Online at https://apisemanticscholar.org/CorpusID:169339509.
[98] K. J. Kossacki and J. Leliwa-Kopystynski, "Temperature dependence of the sublimation rate of water ice: Influence of impurities", *Icarus* **233** (2014), pp. 101–105.

[99] C. Truesdell and R. Toupin, "The classical field theories", in *Principles of Classical Mechanics and Field Theory* (S. Flügge, ed.), *Encyclopedia of Physics*, Springer: New York, 1960.

[100] S. Chandrasekhar, *Hydrodynamic and Hydromagnetic Stability*, Clarendon Press: Oxford, 1961.

[101] D. R. Lide (ed.), *The Handbook of Chemistry and Physics*, CRC Press: Boca Raton, FL, 2004.

[102] S. Carpy and H. Mathis, "Modeling binary alloy solidification by a random projection method", *Numer. Methods Partial Differ. Equ.* **35** (2019), no. 2, pp. 733–760.

[103] M. Alkattan, E. H. Oelkers, J.-L. Dandurand and J. Schott, "Experimental studies of halite dissolution kinetics, 1 The effect of saturation state and the presence of trace metals", *Chem. Geol.* **137** (1997), no. 8, pp. 201–219.

[104] A. W. Woods, "Melting and dissolving", *J. Fluid Mech.* **239** (1992), pp. 429–448.

[105] A. Malyarenko, A. J. Wells, P. J. Langhorne, N. J. Robinson, M. J. Williams and K. W. Nicholls, "A synthesis of thermo-dynamic ablation at ice-ocean interfaces from theory, observations and models", *Ocean Model.* **154** (2020), article no. 101692.

[106] A. L. McCutchan and B. A. Johnson, "Laboratory experiments on ice melting: a need for understanding dynamics at the ice-water interface", *J. Marine Sci. Eng.* **10** (2022), no. 8, article no. 1008.

[107] C. Cenedese and F. Straneo, "Icebergs melting", *Annu. Rev. Fluid Mech.* **55** (2023), pp. 377–402.

[108] S. S. Pegler and M. S. Davies Wykes, "Shaping of melting and dissolving solids under natural convection", *J. Fluid Mech.* **900** (2020), article no. A35.

[109] J. Philippi, M. Berhanu, J. Derr and S. Courrech du Pont, "Solutal convection induced by dissolution", *Phys. Rev. Fluids* **4** (2019), article no. 103801.

[110] M. Nicklasson and A. Brodin, "The relationship between intrinsic dissolution rates and solubilities in the water-ethanol binary solvent system", *Int. J. Pharm.* **18** (1984), no. 1-2, pp. 149–156.

[111] M. Chaigne, S. Carpy, M. Massé, J. Derr, S. Courrech du Pont and M. Berhanu, "Emergence of tip singularities in dissolution patterns", *Proc. Natl. Acad. Sci. USA* **120** (2023), no. 48, article no. e2309379120.

[112] T. S. Sullivan, Y. Liu and R. E. Ecke, "Turbulent solutal convection and surface patterning in solid dissolution", *Phys. Rev. E* **54** (1996), no. 1, article no. 486.

[113] A. J. C. Ladd and P. Szymczak, "Reactive flows in porous media: challenges in theoretical and numerical methods", *Annu. Rev. Chem. Biomol. Eng.* **12** (2021), no. 1, pp. 543–571.

[114] F. Dutka, V. Starchenko, F. Osselin, S. Magni, P. Szymczak and A. J. C. Ladd, "Time-dependent shapes of a dissolving mineral grain: Comparisons of simulations with microfluidic experiments", *Chem. Geol.* **540** (2020), article no. 119459.

[115] M. D. Covington, "Calcite dissolution under turbulent flow conditions: a remaining conundrum", *Acta Carsologica* **43** (2014), no. 1, pp. 195–202.

[116] J. Chittenden, V. Chevrier, L. Roe, K. Bryson, R. Pilgrim and D. Sears, "Experimental study of the effect of wind on the stability of water ice on Mars", *Icarus* **196** (2008), no. 2, pp. 477–487.

[117] A. Lecoanet, F. Payot, C. Journeau, N. Rimbert and M. Gradeck, "Classification of ablation mode during impact of hot liquid jet on a solid", *Int. J. Heat Mass Transfer* **181** (2021), article no. 121883.

[118] S. H. Davis, U. Müller and C. Dietsche, "Pattern selection in single-component systems coupling Bénard convection and solidification", *J. Fluid Mech.* **144** (1984), pp. 133–151.

[119] J. Noel, C. Metivier, S. Becker and S. Leclerc, "Natural convection in phase change material: experimental study", *Int. J. Heat Mass Transfer* **183** (2022), article no. 122047.

[120] A. M. Leitch, "Analog experiments on melting and contamination at the roof and walls of magma chambers", *J. Volcanol. Geotherm. Res.* **129** (2004), no. 1, pp. 173–197.

[121] C. Gau and R. Viskanta, "Melting and solidification of a pure metal on a vertical wall", *J. Heat Transfer* **108** (1986), no. 1, pp. 174–181.

[122] G. Kaufmann, "Modelling karst geomorphology on different time scales", *Geomorphology* **106** (2009), pp. 62–77.

[123] R. Vesipa, C. Camporeale and L. Ridolfi, "Thin-film induced morphological instabilities over calcite surfaces", *Proc. R. Soc. A* **471** (2015), no. 2176, article no. 20150031.

[124] J. Colombani, "Measurement of the pure dissolution rate constant of a mineral in water", *Geochim. Cosmochim. Acta* **72** (2008), no. 23, pp. 5634–5640.

[125] E. Nakouzi, R. E. Goldstein and O. Steinbock, "Do dissolving objects converge to a universal shape?", *Langmuir* **31** (2014), no. 14, pp. 41145–41150.

[126] J. M. Huang, M. N. Moore and L. Ristroph, "Shape dynamics and scaling laws for a body dissolving in fluid flow", *J. Fluid Mech.* **765** (2015), pp. 1–10.

[127] G. Kaufmann and W. Dreybrodt, "Calcite dissolution kinetics in the system CaCO₃–H₂O–CO₂ at high undersaturation", *Geochim. Cosmochim. Acta* **71** (2007), no. 6, pp. 1398–1410.

[128] A. C. English and M. Dole, "Diffusion of sucrose in supersaturated solutions", *J. Am. Chem. Soc.* **72** (1950), no. 7, pp. 3261–3267.
[129] J. K. Gladden and M. Dole, "Diffusion in supersaturated solutions. II. Glucose solutions", *J. Am. Chem. Soc.* **75** (1953), no. 16, pp. 3900-3904.

[130] P. S. Mendes, "The naphthalene sublimation technique", *Exp. Therm. Fluid Sci.* **4** (1991), no. 5, pp. 510-523.

[131] T. K. Sherwood and C. Johannes, "The maximum rate of sublimation of solids", *AIChE J.* **8** (1962), no. 5, pp. 590-593.

[132] M. J. O'Neil, *The Merck Index: An Encyclopedia of Chemicals, Drugs, and Biologicals*, Royal Society of Chemistry: Great Britain, 2013.

[133] C. Cohen, M. Berhanu, J. Derr and S. Courrech du Pont, "Erosion patterns on dissolving and melting bodies (2015 Gallery of Fluid motion)", *Phys. Rev. Fluids* **1** (2016), article no. 050508.

[134] M. S. Davies Wykes, J. M. Huang, G. A. Hajjar and L. Ristroph, "Self-sculpting of a dissolvable body due to gravitational convection", *Phys. Rev. Fluids* **3** (2018), article no. 043801.

[135] S. Weady, J. Tong, A. Zidovska and L. Ristroph, "Anomalous convective flows carve pinnacles and scallops in melting ice", *Phys. Rev. Lett.* **128** (2022), no. 4, article no. 044502.

[136] A. Acrivos, "Mass transfer in laminar-boundary-layer flows with finite interfacial velocities", *AIChE J.* **6** (1960), no. 3, pp. 410-414.

[137] S. S. Pegler and M. S. Davies Wykes, "The convective Stefan problem: Shaping under natural convection", *J. Fluid Mech.* **915** (2021), article no. A86.

[138] J. Schürr, "Sur la vitesse de dissolution des sels dans leurs solutions aqueuses", *J. Phys. Théor. Appl.* **4** (1905), no. 1, pp. 17-26.

[139] J. Schürr, "Recherches sur la vitesse de dissolution des sels dans leurs solutions aqueuses", *J. Chim. Phys. Phys.-Chim. Biol.* **2** (1904), pp. 245-306.

[140] H. Grijseels, D. J. A. Crommelin and C. J. De Blaey, "Hydrodynamic approach to dissolution rate", *Int. J. Clin. Pharm.* **3** (1981), pp. 129-144.

[141] F. Dacheng and A. Hellawell, "The surface morphology of crystals melting under solutions of different densities", *J. Cryst. Growth* **92** (1988), no. 3-4, pp. 364-370.

[142] R. C. Kerr, "Melting driven by vigorous compositional convection", *J. Fluid Mech.* **280** (1994), pp. 255-285.

[143] Y.-C. Yen, "Onset of convection in a layer of water formed by melting ice from below", *Phys. Fluids* **11** (1968), no. 6, pp. 1263-1270.

[144] T. S. Sullivan, W. J. Brown, M. D. Kerkel and R. E. Ecke, "Turbulent convection and surface patterning in solid dissolution", *Phys. Lett. A* **206** (1995), pp. 49-53.

[145] Y. Liu, L. Ning and R. E. Ecke, "Dynamics of surface patterning in salt-crystal dissolution", *Phys. Rev. E* **53** (1996), no. 6, pp. 5572-5575.

[146] D. Gechter, P. Huggenberger, P. Ackerer and H. N. Waber, "Genesis and shape of natural solution cavities within salt deposits", *Water Resour. Res.* **44** (2008), no. 11, article no. W11409.

[147] R. S. Sharma, M. Berhanu and A. Kudrolli, "Alcove formation in dissolving cliffs driven by density inversion instability", *Phys. Fluids* **34** (2022), no. 5, article no. 054118.

[148] B. Villien, Y. Zheng and D. Lister, "Surface dissolution and the development of scallops", *Chem. Eng. Commun.* **192** (2005), no. 1, pp. 125-136.

[149] D. Wang, D. Ewing and C. Y. Ching, "Time evolution of surface roughness in pipes due to mass transfer under different Reynolds numbers", *Int. J. Heat Mass Transfer* **103** (2016), pp. 661-671.

[150] H. Schlichting, *Boundary-Layer Theory*, 7th edition, McGraw Hill: New York, 1979.

[151] T. J. Hanratty, "Stability of surfaces that are dissolving or being formed by convective diffusion", *Annu. Rev. Fluid Mech.* **13** (1981), no. 1, pp. 231-252.

[152] K.-S. Hsu, F. A. Locher and J. F. Kennedy, "Forced-convection heat transfer from irregular melting wavy boundaries", *J. Heat Transfer* **101** (1979), pp. 598-602.

[153] R. R. Gilpin, T. Hirata and K. C. Cheng, "Wave formation and heat transfer at an ice-water interface in the presence of a turbulent flow", *J. Fluid Mech.* **99** (1980), pp. 619-640.

[154] M. Bushuk, D. M. Holland, T. P. Stanton, A. Stern and C. Gray, "Ice scallops: a laboratory investigation of the ice-water interface", *J. Fluid Mech.* **873** (2019), pp. 942-976.

[155] W. Li, J. T. Germaine and H. H. Einstein, "Transport-controlled dissolution in an evolving fracture: The extended Purday solution and fracture flow tests", *Water Resour. Res.* **57** (2021), no. 3, article no. e2020WR029166.

[156] C. Oltean, F. Golfier and M. A. Bués, "Numerical and experimental investigation of buoyancy-driven dissolution in vertical fracture", *J. Geophys. Res.: Solid Earth* **118** (2013), no. 5, pp. 2038-2048.

[157] R. Hu, K. Li, C.-X. Zhou, T. Wang, Z. Yang and Y.-F. Chen, "On the role of gravity in dissolving horizontal fractures", *J. Geophys. Res.: Solid Earth* **128** (2023), no. 3, article no. e2022JB025214.

[158] W. Guo, R. Hu, C.-X. Zhou, Z. Yang and Y.-F. Chen, "Dissolution regimes of a horizontal channel in a gravity field", *Phys. Rev. Fluids* **8** (2023), no. 12, article no. 123902.

[159] L. H. Pitcher and L. C. Smith, "Supraglacial streams and rivers", *Annu. Rev. Earth Planet. Sci.* **47** (2019), pp. 421-452.
[160] E. Mantelli, C. Camporeale and L. Ridolfi, "Supraglacial channel inception: Modeling and processes", *Water Resour. Res.* **51** (2015), no. 9, pp. 7044–7063.

[161] V. Bergeron, C. Berger and M. D. Betterton, "Controlled irradiative formation of penitentes", *Phys. Rev. Lett.* **96** (2006), no. 9, article no. 098502.

[162] T. Bartels-Rausch, V. Bergeron, J. H. E. Cartwright, et al., "Ice structures, patterns, and processes: A view across the icefields", *Rev. Modern Phys.* **84** (2012), no. 2, pp. 885–944.

[163] L. M. Cathles, D. S. Abbot and D. R. MacAyeal, "Intra-surface radiative transfer limits the geographic extent of snow penitents on horizontal snowfields", *J. Glaciol.* **60** (2014), no. 219, pp. 147–154.

[164] W. W. Mullins and R. F. Sekerka, "Stability of a planar interface during solidification of a dilute binary alloy", *J. Appl. Phys.* **35** (1964), no. 2, pp. 444–451.

[165] D. J. Drewry, "A quantitative assessment of dirt-cone dynamics", *J. Glaciol.* **11** (1972), no. 63, pp. 431–446.

[166] N. Taberlet and N. Plihon, "Sublimation-driven morphogenesis of Zen stones on ice surfaces", *Proc. Natl. Acad. Sci. USA* **118** (2021), no. 40, article no. e2109107118.

[167] A. L. Lange, "Introductory notes on the changing geometry of cave structures", *Cave Stud.* **11** (1959), pp. 69–90.

[168] S. Fauve, "Pattern forming instabilities", in *Hydrodynamics and Nonlinear Instabilities* (C. Godrèche and P. Manneville, eds.), Cambridge University Press: Cambridge, 2005.

[169] B. Favier, J. Purseed and L. Duchemin, "Rayleigh-Bénard convection with a melting boundary", *J. Fluid Mech.* **858** (2019), pp. 437–473.