\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRPHYS,Unicode,Thematic,biblatex]{cedram}
\addbibresource{crphys20240562.bib}
\newenvironment{Table}{\begin{table}}{\end{table}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.98\linewidth}{#1}}
\def\tbody{\noalign{\relax}}
\def\hyphen{\text{-}}
\def\ndash{\text{--}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\mathbi#1{\text{\textbf{\textit{#1}}}}
\usepackage{amssymb}
\RequirePackage{etoolbox}
\def\jobid{crphys20240562}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\back#1{}
\def\rmpi{\uppi}
\def\Lbreak{\newline}
\def\ubreak{\unskip\break}
\def\botline{\\\hline}
\def\0{\phantom{0}}
\def\xmorerows#1#2{#2}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\csdef{Seqnsplit}{\\}
\DOI{10.5802/crphys.230}
\datereceived{2024-07-09}
\daterevised{2024-12-18}
\dateaccepted{2024-12-20}
\ItHasTeXPublished
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\makeatother
}
{\documentclass[crphys]{article}
\def\CDRsupplementaryTwotypes#1#2{}
\usepackage[T1]{fontenc}
\let\printbibliography\relax
\def\CDRdoi{10.5802/crphys.230}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\let\ubreak\relax
\makeatletter
}
\makeatother

\usepackage{upgreek}

\begin{document}

\begin{noXML}

\makeatletter
\def\TITREspecial{\relax}
\def\cdr@specialtitle@english{Geophysical and astrophysical fluid dynamics in the laboratory}
\def\cdr@specialtitle@french{Dynamique des fluides g\'eophysiques et astrophysiques au laboratoire}
\makeatother

\CDRsetmeta{articletype}{review}

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

\alttitle{De l'\'echelle du laboratoire aux paysages plan\'etaires,
quelles sont les empreintes laiss\'ees par des \'ecoulements de surface
sur des substrats solides en ablation soumis \`a des changements de
phase ?}

\author{\firstname{Sabrina} \lastname{Carpy}\CDRorcid{0000-0002-2705-6974}\IsCorresp}
\address{LPG (Laboratoire de Plan\'etologie et G\'eosciences), Nantes
Universit\'e, CNRS (UMR 6112), 44322 Nantes, France}
\email[S. Carpy]{sabrina.carpy@univ-nantes.fr}

\author{\firstname{Michael} \lastname{Berhanu}\CDRorcid{0000-0001-9099-2135}}
\address{MSC (Mati\`ere et Syst\`emes Complexes), Universit\'e Paris
Cit\'e, CNRS (UMR 7057), 75013 Paris, France}

\author{\firstname{Martin} \lastname{Chaigne}\CDRorcid{0009-0001-6937-9931}}
\addressSameAs{2}{MSC (Mati\`ere et Syst\`emes Complexes), Universit\'e
Paris Cit\'e, CNRS (UMR 7057), 75013 Paris, France}

\author{\firstname{Sylvain} \lastname{Courrech~du Pont}\CDRorcid{0000-0002-4133-3536}}
\addressSameAs{2}{MSC (Mati\`ere et Syst\`emes Complexes), Universit\'e
Paris Cit\'e, CNRS (UMR 7057), 75013 Paris, France}

\shortrunauthors

\thanks{French Research Council (SHERPAS ANR-23-CE49-0006 and
PhysErosion ANR-22-CE30-0017)}

\keywords{\kwd{Sublimation}\kwd{Melting}\kwd{Dissolution}\kwd{Natural
solid bedforms}\kwd{Experimental analogues}\kwd{Pattern
\mbox{formation}}\kwd{Classification}}

\altkeywords{\kwd{Sublimation}\kwd{Fusion}\kwd{Dissolution}\kwd{Motifs 
solides naturels}\kwd{Analogues exp\'erimentaux}\kwd{Formation de
motifs}\kwd{Classification}}

\begin{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
\mbox{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,
\mbox{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 meltwater/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 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.
\end{abstract}

\begin{altabstract}
Les processus en ablation, tels que les changements de phase, physiques
(sublimation/fusion) ou chimiques (dissolution), sont \'etudi\'es dans
de nombreuses applications techniques en physique. Pour les interfaces
solide/fluide, l'interaction entre un changement de phase et un
\'ecoulement peut conduire \`a la formation de motifs topographiques.
Dans ce cas, la m\'ecanique des fluides associ\'ee \`a de tels
changements de phase joue un r\^ole primordial dans l'\'evolution des
paysages terrestres et plan\'etaires, observ\'es par les sondes en
orbite autour des plan\`etes et des lunes. Sur Terre, la glace de mer,
les glaciers et les plateaux karstiques s'\'etendent sur des m\`etres
ou des kilom\`etres. L'\'echelle de ces paysages contraste avec
l'\'echelle des m\'ecanismes physiques qui r\'egissent leur dynamique
\'evolutive. C'est en effet la taille typique des couches limites
atmosph\'eriques ou des films d'eau de fonte/solut\'e qui limite le
transfert de chaleur/concentration \`a l'interface de changement de
phase/dissolution, et donc le taux d'ablation du solide. Dans de
nombreuses situations, ces couches sont contr\^ol\'ees par
l'\'ecoulement du fluide, de type convection naturelle ou forc\'ee.
Dans le premier cas, la flottabilit\'e de la fonte/dissolution
elle-m\^eme, entra{\^\i}ne une stratification en densit\'e caus\'ee par
des gradients de temp\'erature/concentration, qui peut \^etre stable ou
instable. Dans le second cas, l'\'ecoulement forc\'e par des vents ou
des pentes peut \^etre consid\'er\'e comme un \'ecoulement de hauteur
infinie ou finie, tels que les \'ecoulements \`a surface libre. Dans
tous les cas, le flux de masse modifie la topographie, qui en retour
impacte les \'ecoulements de couches limites (fines ou \'epaisses) et
donc le taux d'ablation de mani\`ere r\'etroactive. Dans la nature, la
r\'etroaction positive entre la g\'eom\'etrie et le transfert de masse
entra{\^\i}ne la formation spontan\'ee de formes caract\'eristiques \`a
diff\'erentes \'echelles. Ces formes ne sont pas seulement des
curiosit\'es g\'eologiques, comme le sont par exemple les \og Zen stones\fg, ces pierres sur un
pi\'edestal de glace ou encore les c\^ones de poussi\`ere, mais des
marqueurs des processus hydrodynamiques en jeu. De nombreux paysages
sont fa\c{c}onn\'es selon des motifs r\'eguliers et r\'ep\'et\'es,
qu'ils soient pointus, en forme de cupules, \`a ligne de cr\^etes
parall\`eles entre elles ou encore en marches d'escaliers. Dans ce
papier de synth\`ese, nous recensons les \'etudes exp\'erimentales
d\'edi\'ees \`a l'\'etude des motifs g\'en\'er\'es par diff\'erents
modes de transport des flux sur des substrats solides soumis \`a des
changements de phase physiques ou chimiques. Nous souhaitons ainsi
mettre en avant le r\^ole du mode de transport des flux dans la
diversit\'e des motifs observ\'es sur des mat\'eriaux analogues.
Comprendre la diversit\'e de ces motifs est essentiel pour \'evaluer
les conditions environnementales dans lesquelles ils se forment,
notamment sur des plan\`etes comme Mars ou Pluton sur lesquelles les
changements de phase jouent un r\^ole g\'eomorphologique majeur.
\end{altabstract}

\maketitle

{\vspace*{.6pc}}

\end{noXML}

\hyphenpenalty 10000
\exhyphenpenalty 10000
\sloppy

\section{Ablation by mass transfer on solid substrate submitted to
phase change} \label{Intro}

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~\cite{Pascheke2008} or ablation processes in
space capsules~\cite{Combs2014}. Ablation patterns may appear during
the descent phase of an atmospheric space capsules~\cite{Hochrein1976,
Schneider2006}, similar to meteorites' regmaglypts~\cite{Lin1987}, 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~\cite{Thomas1979}, 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~\cite{Smith2013} and the bladed terrain of Pluto, interpreted as
\textit{penitentes}~\cite{Moores2017}. Sublimation of ice combined with
surface flows like wind could be a major geomorphic agent in explaining
these planetary landscapes~\cite{Mangold2011}. On Earth water ice on
its solid phase represents 2.4\% of the global volume of water, but
covers 14\% of the Earth' surface~\cite{Ohmura2004}. 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. $\mathrm{N}_2$,
$\mathrm{CH}_4$, $\mathrm{CO}_2$, $\mathrm{NH}_3$ ices~\cite{Carpy2023}, are
subjected to surface sublimation that prevails in the solar system. 
Soluble rocks on Earth also present solid bedforms such as scallops or
\textit{rillenkarrens} that could be found on limestone, gypsum and
salt~\cite{ford2007karst, kozary1968}. 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~\cite{Gallaire2017}. 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 \textit{rillenkarren}~\cite{Glew1980, Slabe2016, Guerin2020},
scallops and \textit{flutes}, rely on the use of
plaster~\cite{Allen1970, Goodchild1972, Blumberg1974} or even
caramel~\cite{Cohen2020}, for example. In many situations, fluid
dynamic instabilities~\cite{Gallaire2017} (e.g.\ the \textit{Rayleigh
Plateau}, \textit{Kelvin Helmholtz}, \textit{Rayleigh Taylor} or
\textit{Rayleigh B\'enard} 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~\cite{Claudin2017}, or the ripples and cyclic steps on ice due
to run-off flows~\cite{Camporeale2012, Yokokawa2016}, 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
\textit{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~\cite{Guerin2020, bertagni2021hydrodynamic,
Djambov_Gallaire_2023} 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~\cite{Blumberg1974,
Thorsness1979, Chan2007, Camporeale2012a, Claudin2017, Guerin2020}, the
melting~\cite{AshtonKennedy1972, Ashton1972} and the
sublimation~\cite{Bordiec2020} 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~\cite{Thomas1979, Allen1982} in the
case of ice sublimation, ice melting and rock dissolution will be
described in Section~\ref{natural} for: sharp patterns; polygonal
depression; crested bedforms parallel or perpendicular to the flow;
stepped patterns; and debris-covered ice structures. In
Section~\ref{equations}, 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~\ref{experiments} 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~\ref{discussion} focuses on the genetic link between the type
of flow and the bedfoms 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.

{\vspace*{-2pt}}

\section{Diversity of natural solid bedforms by phase change on solid
substrate} \label{natural}

{\vspace*{-2pt}}

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~\ref{experiments}.

{\vspace*{-2pt}}

\subsection{Sharp patterns} \label{Penitents}

{\vspace*{-2pt}}

First described by Charles Darwin in Chile in 1835, \textit{penitentes}
are shapes formed by sublimation in old, compact ice/snow, either as
spikes or as blades pointing skywards
(Figure~\ref{fig:naturalsharp}(a)). On Earth, \textit{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~\cite{Amstutz1958}; on the Quelccaya Ice Cap
in Peru~\cite{Sinclair2016}; and on the Tapado Glacier in northern
Chile~\cite{Sinclair2016}. But similar forms have also been observed at
very low altitudes such as Greenland~\cite{Odell1941}. Given the
climatic conditions that favour their formation, \textit{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 \textit{penitentes}, such as the Bladed
Terrain Deposit of Pluto's \textit{Tartarus Dorsa}
region~\cite{Moores2017, Moore2018}. It has been suggested that
\textit{penitentes} as high as 15~m may exist on Europa's
surface~\cite{Hobley2018}, although the resolution of currently
available images of Europa's surface is insufficient to verify this,
but this hypothesis is debated because of the low atmospheric pressure
on Europa~\cite{Hand2020, Hobley2020}. From base to
tip, these spikes can reach heights of a few
centimetres~\cite{Hand2020} to several metres~\cite{Lliboitry1954},
depending on environmental conditions. The cross-section of the
\textit{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~\cite{Betterton_2001}. The orientation of the
\textit{penitentes} in terms of azimuth is linked to the local geometry
of the insulation~\cite{Hastenrath1981}. In particular, Matthes
{et~al.}~\cite{Matthes1934} observes an inclination of the
\textit{penitentes} following the elevation of the midday sun. For
geometric reasons, sublimation is more effective in hollows. 
Illumination plays a key role~\cite{Betterton2000}, as do heat
transport in the ice and the transport of vapor from the
ice~\cite{Claudin2015}. 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 \textit{penitentes} and other bedforms, such as
scallops (Figure~\ref{fig:naturalscallops}(a)) or transverse linear
waves (Figure~\ref{fig:naturalripples}(a)).

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig:naturalsharp}\textit{Penitentes} and
\textit{Tsingy}. (a)~\textit{Penitentes} (Rio Blanco, Central Andes of
Argentina). Cr\'edits: P.~Dubuc, Arvaki, CC BY-SA 3.0 license.
(b)~{\textit{Tsingy} of Bemaraha~\cite{Veress2013, Huang2020}. Credits:
A.~S.~Maloney, CC BY-SA 3.0 license}.}
\end{figure} 

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:naturalscallops}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\`eche, France).
Credits: M.~Bordiec/S.~Carpy/O.~Bourgeois, reprinted with permission
from~\cite{bordiec_phd2020}. (d)~Limestone scallops by dissolution
(Moulin de Vog\"u\'e, France). Credits:
M.~Bordiec/S.~Carpy/O.~Bourgeois, reprinted with permission
from~\cite{bordiec_phd2020}.}
\end{figure} 

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:naturalripples}Ripple-like bedforms over different
substrates submitted to mass transfer: (a)~by sublimation (blue ice
ripples in \textit{Svea}, Antarctica adapeted from~\cite{Bintanja1999}.
Copyright 1999 by the American Geophysical Union). (b)~By melting,
Rivi\`ere St-Croix ice ripples, Wisconsin, USA~\cite{Carey1966}. (c)~By
sublimation, martian North polar Cap, reprinted with permission
from~\cite{bordiec_phd2020}. \textit{HiRISE} PSP\_009689\_2645. (d)~By
dissolution, \textit{Coups de Gouge}, Grotte de Saint-Marcel
d'Ard\`eche, France. Credits: M.~Bordiec/S.~Carpy/O.~Bourgeois,
reprinted with permission from~\cite{bordiec_phd2020}.}
\end{figure} 

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
\textit{Tsingy}~\cite{Veress2013, salomon2009tsingy}, in
Malaysia~\cite{day2009pinnacle} or in south
China~\cite{knez2002lithologic}. Spectacular sharp and vertical
pinnacles, tens of meters high, rise from the surface of a horizontal
\textit{Karst} plateau (Figure~\ref{fig:naturalsharp}(b)). The general
shape is needle-shaped, with a slope of 15{\textdegree} to
20{\textdegree}, which remains even when the material has been strongly
dissolved. These sharp-edged \textit{karren} forms, which develop in
groups, have several degrees of evolution: from dragon teeth to
sharp-edged pinnacles before evolving to
\textit{Tsingy}~\cite{Gines2009}. The smallest in size belong to the
\textit{lappies}' family (2--3~m) and the largest, by one order of
magnitude, belong to the \textit{mega-lappies}' family (20--30~m). On
their slopes, centimeter scale \textit{rillenkarren} (see
Section~\ref{rillenkarrens}) 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~\cite{Frumkin2013}.

\subsection{Polygonal depression} \label{Polygonal}

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 \textit{penitentes} on Earth, they are located
at lower altitudes, where the ambient temperature is generally higher
(e.g.\ in the Alps~\cite{Betterton2000}), and in areas more exposed to
winds than those where \textit{penitentes} develop. With amplitudes
ranging from a few centimetres to half a metre, they do not reach
heights comparable to those of the \textit{penitentes}, although they
are also controlled by solar radiation and influenced by albedo and the
amount of dust in the ice/snow~\cite{Rhodes1987, Betterton2000}; the
role of wind is not clearly excluded. Suncups could constitute a
surface texture found particularly on the North Polar Cap of
Mars~\cite{Malin2001, Milkovich2006, Banks2010, Masse2010, Mangold2011}
and on the Sputnik plain of Pluto~\cite{Telfer2018}. Similar shapes
appear on the walls of ice caves~\cite{Bella2007, Persoiu2018,
Obleitner2011} 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~\cite{ford2007karst}.
Conceptual models of the climate prevailing in these ice caves have
been built from measurement campaigns conducted over several
decades~\cite{Persoiu2011}, in particular in the \textit{Eisriesenwelt}
cave (Salzburg, Austria)~\cite{Obleitner2011}, 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~{\textdegree}C and the air is dry~\cite{Bordiec2020}, which favors
episodes of sublimation, enhanced by cold, dry air that is blown and
circulates through the cave, at speeds not exceeding
$0.2~\mathrm{m}{\cdot}\mathrm{s}^{-1}$~\cite{Curl1966, Bella2007,
Obleitner2011}. 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 effective, the annual ablation rate in ice caves ranges from
$3~\mathrm{mm}{\cdot}\mathrm{year}^{-1}$~\cite{Marshall1974} to
$35~\mathrm{mm}{\cdot}\mathrm{year}^{-1}$~\cite{Persoiu2018}, resulting in
sublimation rates of between $3.8\times
10^{-10}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$ and $4.5\times
10^{-9}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$.

Scallops are widespread polygonal depression patterns, both on
ice~\cite{Bella2007, Persoiu2018}
(Figure~\ref{fig:naturalscallops}(a--b)) and on soluble rock
walls~\cite{ford2007karst, Gines2009}
(Figure~\ref{fig:naturalscallops}(c--d)), and on many other
substrates~\cite{Thomas1979}. By analogy, Meakin and
Jamtveit~\cite{Meakin2010} suggest that these shapes can be likened to
those observed on certain meteorites~\cite{Lin1987}, on corroding metal
pipe surfaces~\cite{Villien2001} and cohesive surfaces in fluvial
systems~\cite{Allen1971}. 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~\ref{fig:naturalscallops}).
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~\ref{fig:naturalscallops}(d)). Similarly,
scallops have sometimes been assimilated to shapes transverse to the
flow direction~\cite{Thomas1979, Claudin2017}, although they differ
from ripples and flutes, which are transverse waves
(Figure~\ref{fig:naturalripples}), which are considered to be special
cases of scallops~\cite{Blumberg1974} (see Section~\ref{wavy}). 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~\cite{Meakin2010}.

\subsection{Crested patterns perpendicular to the flow direction}
\label{wavy}

Transverse linear waves (Figure~\ref{fig:naturalripples}) 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~\cite{ford2007karst,
Gines2009}. The waves are likely to migrate in the direction of
flow over time. Their cross-sectional profile in relation to the
\mbox{direction} of flow is slightly asymmetrical
(Figure~\ref{fig:naturalripples}(b)), with a downstream slope greater
than the upstream slope \cite{Blumberg1974,Hammer2011}. They exist on
all types of solid substrates, such as salt~\cite{Karcz1987} and ice,
immersed in liquid water~\cite{Ashton1972, Camporeale2012}, in ice
caves~\cite{Curl1966} or on Antarctic ice surfaces exposed to strong
winds~\cite{Bintanja2001}. If dissolution forms are naturally very
common in limestone caves~\cite{Curl1966, Allen1971, Goodchild1972}
with a range of wavelengths from a few centimetres to several metres
(Figure~\ref{fig:naturalripples}(d)), measurements of the flow
characteristics responsible for their formation are rare in these
natural environments. Sublimation waves
(Figure~\ref{fig:naturalripples}(a)) are periodic, linear and
transverse bedforms of different scales (from centimetres to metres) 
depending on their dry windy environment~\cite{Bintanja2001,
Weller1969}. On Earth, Blue Ice Areas that cover 1\% of Antarctica
reveal the presence of sublimation ice ripples~\cite{Schneider2006},
about 10~cm in wavelength that emerge in ${\sim}$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~\cite{Bordiec2020}. Shapes on the North Polar Cap
of Mars (Figure~\ref{fig:naturalripples}(c)) and certain shapes on the
surface of Pluto (plain \textit{Sputnik Planitia}) can also be
interpreted as sublimation waves~\cite{Carpy2023}.

\subsection{Crested patterns parallel to the flow direction}
\label{rillenkarrens}

\looseness=-1
The term \textit{Karren} designates the dissolution patterns, which are
created by the run-off flows on \mbox{denuded} rock
walls~\cite{ford2007karst, Gines2009, Veress2019}. 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 \textit{Karren} patterns occur
mainly by dissolution of limestone and dolomite, but they can also
appear on gypsum and salt (Halite)~\cite{Frumkin2013}. In view of the
wide variety of shapes and front, a specific terminology has been
defined to classify the different types~\cite{Veress2019}. Among the
most important examples, \textit{Rillenkarren} are adjacent grooves or
small channels directed along the main slope~\cite{Lundberg2009,
Veress2019} (Figure~\ref{fig:ImageRillenkarrenCras}(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 \textit{Rillenkarren} found on Madre de
Dios Island in Chile (Figure~\ref{fig:ImageRillenkarrenCras}(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~\cite{Glew1980}. Moreover, \textit{Rillenkarren} grooves are
preferentially observed on the sloped wall (typically 45{\textdegree})
above a smooth area, where the water film is thicker~\cite{Glew1980}.
The \textit{Wandkarren} are a special case of the \textit{Rillenkarren}
and occurs on nearly vertical walls created by large rock fractures.
The \textit{Rinnenkarren} or Runnels are dissolution channels of larger
size, which are formed once the water has been channelized in a
rivulet. \textit{Rundkarren} are similar to \textit{Rinnenkarren}, but
the ridges between the channels are rounded by subsoil
dissolution~\cite{Veress2019} or maybe by snow cover.
\textit{Meanderkarren} are \textit{Rinnenkarren} with meanders. The
flow dynamics of water channelized in \textit{Rinnenkarren} systems
have been studied on the field and in laboratory reproduction with non
soluble materials~\cite{Veress2013}.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig:ImageRillenkarrenCras}(a)~Nearly vertical
\textit{Rillenkarren} on limestone. The photo was taken in Velebit
(Croatia). The typical groove width is about 1~cm. Credits S.~Courrech
du Pont. (b)~\textit{Rillenkarren} of the order of one meter on Madre
de Dios Island in Chile, Ultima Patagonia expedition, \textit{Centre
Terre}. Note the human at the bottom of the picture, to get an idea of
the characteristic scale. From Tourte {et~al.}~\cite{Tourte2017}.
(c)~Inclined limestone pavement or \textit{Lapiaz} at the site
\textit{Combe Oursi\`ere} near Lans en Vercors, France. Rounded
\textit{Rillenkarren} are directed along the main slope of about
45{\textdegree} inclination. Credits M.~Berhanu, M.~Chaigne and
J.~Derr. (d)~3D reconstruction of \textit{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.}
\end{figure} 

\subsection{Stepped bedforms} \label{steps}

Cyclic steps are bedforms that are shaped like a succession of steps,
moving against the flow~\cite{Kennedy1969}. These steps are associated
with hydraulic jumps at the free surface. Cyclic steps occur at the
bottom of \textit{bedi\`eres} (Figure~\ref{fig:naturalsteps}(a)). They
have been studied experimentally on sloping ice over which liquid water
flows~\cite{Yokokawa2016}. Cyclic steps have also been described on
solid rock substrates over which rivers flow~\cite{Izumi2017}. 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~\cite{Smith2013}.

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig:naturalsteps}(a)~Cyclic steps by melting in
\textit{B\'edi\`eres}. Cr\'edit: L.~Karlstrom, adpated
from~\cite{Izumi2017} with permission of John Wiley and Sons.
(b)~\textit{Trittkarrens} steps perpendicular to the crests of
\textit{Rillenkarrens}. Cr\'edit: L.~Perrella, CC BY-SA 3.0 license.}
\end{figure}

There are also stepped, regular undulations known as
\textit{trittkarrens}, \textit{stepkarrens} or \textit{heelprint
karrens} (Figure~\ref{fig:naturalsteps}(b)), which are observed in
dissolving rocks with a gentle slope. They are perpendicular to the
sharp ridge lines of the \textit{rillenkarrens} described in
Section~\ref{rillenkarrens} and therefore perpendicular to the flow.
They are small-scale \textit{karrens} (on the order of centimetres)
that can form on various dissolving rocks, such as limestone and
gypsum~\cite{Mottershead2000}. Morphologically, they are described by
Ford \& Williams~\cite{ford2007karst} as arched-front forms with a flat
floor, 10 to 30~cm in diameter, which are open at the end of the
downslope. Mottershead~\cite{Mottershead2000} characterizes them more
by their ``stepped shape'' in cross-section and the presence of a
slight gradient.

\subsection{Debris-covered ice structures} \label{NaturalCovered}

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~\ref{fig:CappedPattern}~\cite{Evatt_2017, Fowler_2017,
Henot_2021, Henot_2022, Henot_2023}. Glacier tables are large rocks
standing on an ice foot (Figure~\ref{fig:CappedPattern} top), dirt
cones are ice cones covered by a layer of debris
(Figure~\ref{fig:CappedPattern} bottom left), and ice sails (or ice
pyramids) are clean ice structures protruding from the debris-covered
surface (Figure~\ref{fig:CappedPattern} bottom right).

\begin{figure}
{\vspace*{2pt}}
\includegraphics{fig06}
{\vspace*{2pt}}
\caption{\label{fig:CappedPattern}Melting and sublimation structures
due to debris or rock cover. Top: Glacier table (rock supported by a
pillar of ice) forming on the \textit{Mer de Glace} glacier in the
French Alps from H\'enot {et~al.}~\cite{Henot_2022}. 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~\ref{SectionCapExp}), figure adapted under the CC BY 4.0
license. Bottom left: Dirt cones (ice cones covered by a layer of
grains) on S\'olheimaj\"okull 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~\cite{Evatt_2017}. Credits: C.~Mayer, 2011, under the CC BY 4.0
license.}
{\vspace*{4pt}}
\end{figure}

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~\cite{Ostrem_1959, Mattson_1993, Evatt_2015}.
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~\cite{Evatt_2015}. 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~\cite{Evatt_2015}.

\section{Modelization of phase change dynamics in presence of a flow}
\label{equations}

{\vspace*{2pt}}

The natural bedforms reviewed in Section~\ref{natural} 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~\cite{Meakin2010}. The set of
equations which will be found for melting, sublimation or dissolution
will also be valid for solidification, condensation or precipitation.
\looseness=1
 
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 \mbox{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. \mbox{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 \mbox{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. 

\subsection{Boundary conditions at an evolving interface} \label{BC}

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.

\subsubsection{Kinetics of the phase change} \label{sec3.1.1}

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~\cite{Worster2000, Brantley2008}.

\paragraph*{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~\cite{AlexiadesSolomonBook}. 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_{\mathrm{m}}$, which can be determined by the Clapeyron
relation~\cite{Callenbook}. For example, when ice melts at atmospheric
pressure on Earth, the melting temperature during the phase change is
equal to 0~{\textdegree}C. The melting temperature also decreases with
the local curvature of the interface due to Gibbs--Thomson
effect~\cite{Worster2000}, 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~\ref{sec3.1.3}).

\paragraph*{Dissolution}

In contrast, in the case of dissolution, kinetic effects are often
considered relevant~\cite{Brantley2008}. 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_{\mathrm{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_{\mathrm{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~\cite{Aaron1970}. 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.

{\vspace*{2pt}}

\paragraph*{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_{\mathrm{sat}}$. This
quantity can vary strongly as a function of the temperature according
to the Clapeyron relation~\cite{Callenbook}.

{\vspace*{2pt}}

\subsubsection{Relationship between phase change rate and kinetic
effects} \label{sec3.1.2}

{\vspace*{2pt}}

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_{\mathrm{i}} -
c_{\mathrm{sat}}$, where $c_{\mathrm{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~\cite{Colombani2007} or gas kinetics~\cite{Eames1997,
Mas2003, Kossacki2014}.

If $c_{\mathrm{i}} < c_{\mathrm{sat}}$, the solid recedes by dissolution
(sublimation) and vice versa if $c_{\mathrm{i}} > c_{\mathrm{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:
{\vspace*{2pt}}
{\begin{equation} \label{fluxchim}
\mathcal{R}(\mathbf{x},t)=\alpha(c_{\mathrm{sat}} - c_{\mathrm{i}})=-\rho_{\mathrm{s}}
\,\mathbf{v}_{\mathrm{i}} (\mathbf{x},t) \cdot \mathbf{n},
\end{equation}}\unskip
where $\rho_{\mathrm{s}}$ is the solid density and $\mathbf{v}_{\mathrm{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 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~\cite{Eames1997, Mas2003,
Kossacki2014}:
{\begin{equation}
\alpha = \varepsilon \sqrt{\frac{{R}T_{\mathrm{s}}}{2\uppi M_{\mathrm{s}}}}.
\end{equation}}\unskip
$\alpha$ is thus a function of the surface temperature of the solid
$T_{\mathrm{s}}$ and depends on the molar mass of sublimated molecules~$M_{\mathrm{s}}$.

\subsubsection{Moving interface} \label{sec3.1.3} \label{moving}

In Equation~(\ref{fluxchim}), the concentration at the interface,
$c_{\mathrm{i}}$, must be determined to solve the dissolution or sublimation
problems. Therefore, to determine the phase change rate
$\mathbf{v}_{\mathrm{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_{\mathrm{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~\cite{truesdell1960classical,
AlexiadesSolomonBook} 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 
{\begin{equation} \label{localconservation}
\frac{\partial f(\mathbf{x},t)}{\partial t} +\nabla \cdot
(\mathbf{J}_{\mathbf{f}}(\mathbf{x},t))=0, 
\end{equation}}\unskip
where $\mathbf{J}_{\mathbf{f}}(\mathbf{x},t)$ is the incoming flux of
the scalar quantity~$f$, which can be mass, concentration or
temperature and $\nabla \cdot (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~\ref{SchemaJump}, where the two phases considered are the solid
(s) and the liquid (l), the Kotchine's theorem gives:
{\begin{equation} \label{Jumprelation}
(f^{\mathrm{l}}-f^{\mathrm{s}})(\mathbf{v}_{\mathrm{i}} \cdot \mathbf{n})=(\mathbf{J}^{\mathrm{l}}-
\mathbf{J}^{\mathrm{s}})\cdot \mathbf{n}.
\end{equation}}\unskip
This relation gives the boundary condition that relates the interface
velocity $\mathbf{v}_{\mathrm{i}}$ to the quantities on each side of the
interface: $f^{\mathrm{l}}$ and $f^{\mathrm{s}}$, and the flux of~$f$, i.e.\ 
$\mathbf{J}^{\mathrm{l}}$ and $\mathbf{J}^{\mathrm{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~(\ref{localconservation}).

\begin{figure}
\includegraphics{fig07}
\caption{\label{SchemaJump}Scheme for determining the jump conditions
at the evolving interface $S(t)$ separating two phases, for example a
solid $\mathcal{R}^{\mathrm{s}}$ and a liquid $\mathcal{R}^{\mathrm{l}}$. We consider an
infinitesimal spatial domain $\Omega(t)$ on either side of the
interface $S(t)$. $\mathbf{v}_{\mathrm{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^{\mathrm{s}}$ and $\partial \Omega^{\mathrm{l}}$ for the part in contact with
$\mathcal{R}^{\mathrm{s}}$ and $\mathcal{R}^{\mathrm{l}}$, respectively.}
\end{figure} 

\paragraph*{Mass balance}

The mass flux corresponds to an advection term
$\mathbf{J}_{\mathbf{f}}(\mathbf{x},t) = \rho \mathbf{u}$ of the
density field $f(\mathbf{x},t)=\rho$ by the velocity
field~$\mathbf{u}$. Substituting this into the
Equation~(\ref{localconservation}) gives the conservation of mass, in
three dimensions:
{\begin{equation}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\mathbf{u})=0.
\end{equation}}\unskip
Considering an evolving interface $S(t)$ of local velocity
$\mathbf{v}_{\mathrm{i}}$ between a liquid phase and a solid phase (subscripts
$l$ and $s$, respectively), Figure~\ref{SchemaJump}, the jump relation,
Equation~(\ref{Jumprelation}) gives:
{\begin{equation} \label{Jumprelation_mass}
(\rho_{\mathrm{i}}-\rho_{\mathrm{s}})(\mathbf{v}_{\mathrm{i}} \cdot \mathbf{n})=(\rho_{\mathrm{i}}
\mathbf{u}_{\mathrm{i}}\cdot\mathbf{n}-\rho_{\mathrm{s}}\mathbf{u}_{\mathrm{s}}\cdot\mathbf{n})
\end{equation}}\unskip
$\mathbf{u}_{\mathrm{i}}$ and $\rho_{\mathrm{i}}$ are the velocity and the density in the
liquid phase at the position of the interface. $\mathbf{u}_{\mathrm{s}}$ is the
velocity in the solid phase, which is zero (no advection in the solid)
and $\rho_{\mathrm{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:
{\begin{equation} \label{u_interface3D}
\mathbf{u}_{\mathrm{i}} \cdot \mathbf{n}=(\mathbf{v}_{\mathrm{i}} \cdot \mathbf{n})
\left(1-\frac{\rho_{\mathrm{s}}}{\rho_{\mathrm{i}}}\right).
\end{equation}}\unskip
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}$.

\paragraph*{Balance of solute or vapor concentration}

The local conservation relation (Equation~(\ref{localconservation}))
applied to the mass solute concentration $c(\mathbf{x},t)$ reads:
{\begin{equation}
\frac{\partial c}{\partial t}+\nabla\cdot(\mathbf{J}_{\mathbf{c}})=0
\quad\mathrm{with}\ \mathbf{J}_{\mathbf{c}}=c\mathbf{u}-D\nabla c,
\end{equation}}\unskip
where $\nabla$ is the gradient operator. As an example of the incoming
flux $\mathbf{J}_{\mathbf{c}}$, we consider a volume element including
the solid/liquid interface (see Figure~\ref{SchemaJump}) 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: $\mathbf{J}_{\mathbf{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):
{\begin{equation} \label{AdvDiffconc}
\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c=D \nabla^2 c.
\end{equation}}\unskip

When considering a pure solid phase (i.e., concentration and density in
the solid are equal) of the same chemical specie as the solute,
specifying the fluxes in Equation~(\ref{Jumprelation}) for the
concentration gives:
{\begin{equation} \label{Jumprelation_solute}
(c_{\mathrm{i}}-\rho_{\mathrm{s}})(\mathbf{v}_{\mathrm{i}} \cdot \mathbf{n})=c_{\mathrm{i}}\mathbf{u}_{\mathrm{i}}
\cdot \mathbf{n}-D\nabla c \vert_{\mathrm{i}} \cdot \mathbf{n}-\rho_{\mathrm{s}}
\mathbf{u}_{\mathrm{s}}\cdot \mathbf{n}.
\end{equation}}\unskip
The diffusive flux $-D\nabla c\vert_{\mathrm{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
($\mathbf{u}_{\mathrm{s}}=\mathbf{0}$), one finds:
{\begin{equation} \label{cv_solute}
\rho_{\mathrm{s}} \mathbf{v}_{\mathrm{i}}\cdot \mathbf{n}=c_{\mathrm{i}}(\mathbf{v}_{\mathrm{i}}\cdot
\mathbf{n}-\mathbf{u}_{\mathrm{i}}\cdot \mathbf{n})+D(\nabla c\vert_{\mathrm{i}} \cdot
\mathbf{n}). 
\end{equation}}\unskip
Substituting the value of $\mathbf{u}_{\mathrm{i}} \cdot \mathbf{n}$ from
Equation~(\ref{u_interface3D}) into this expression gives a general
boundary condition for the interface velocity $\mathbf{v}_{\mathrm{i}}$:
{\begin{equation} \label{cvd3D}
(\mathbf{v}_{\mathrm{i}}\cdot \mathbf{n})\rho_{\mathrm{s}} \left(1-\frac{c_{\mathrm{i}}}
{\rho_{\mathrm{i}}}\right)=D(\nabla c\vert_{\mathrm{i}} \cdot \mathbf{n}). 
\end{equation}}\unskip
Moreover, the interface velocity is also related to the mass solute
flux $\mathcal{R}$ given by Equation~(\ref{fluxchim}).

\paragraph*{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_{\mathrm{p}}\,T$ where $c_{\mathrm{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~\cite{Chandrasekhar}:
{\vspace*{3pt}}
{\begin{equation} \label{AdvTemp}
\frac{\partial(\rho c_{\mathrm{p}}T)}{\partial t}+ \nabla \cdot (\rho c_{\mathrm{p}}
T\mathbf{u})= \nabla \cdot (\lambda \nabla T),
\end{equation}}\unskip
where $\lambda$ is the thermal conductivity. The heat transport by
diffusion is modeled using the Fourier law, $\mathbf{J}_{\mathbf{T}}
=-\lambda \nabla T$, where $\mathbf{J}_{\mathbf{T}}$ is the heat flux.
It has the dimension of a power per surface unit. If the properties of
the fluid ($\rho$, $c_{\mathrm{p}}$ and $\lambda$) are homogeneous, this last
equation resumes to the classic advection--diffusion equation for the
temperature field:
{\vspace*{3pt}}
{\begin{equation} \label{AdvDiffTemp}
\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla T=\kappa\nabla^2T,
\end{equation}}\unskip
where $\kappa=\lambda/(\rho c_{\mathrm{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_{\mathrm{r}}$, where $h_{\mathrm{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_{\mathrm{p}}$: 
{\vspace*{3pt}}
{\begin{equation} \label{eqcp}
h_{\mathrm{r}}=c_{\mathrm{p},2} T_{\mathrm{i}}-c_{\mathrm{p},1} T_{\mathrm{i}}, 
\end{equation}}\unskip
where $T_{\mathrm{i}}$ is the interface temperature. For the melting of a pure
substance at constant temperature, $T_{\mathrm{i}}$ and $h_{\mathrm{r}}$ identify to the
melting temperature $T_{\mathrm{m}}$ and to the latent heat $L_{\mathrm{m}}$. The melting
is an endothermic process and $L_{\mathrm{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)~\cite{Handbook}.

Applying Kotchine's theorem to energy conservation
(Equation~(\ref{AdvTemp})) gives: 
{\vspace*{3pt}}
{\begin{equation}
[\rho_{\mathrm{i}}c_{\mathrm{p}, \mathrm{l}}T_{\mathrm{i}}-\rho_{\mathrm{s}}c_{\mathrm{p}, \mathrm{s}}T_{\mathrm{m}}] \mathbf{v}_{\mathrm{i}} \cdot
\mathbf{n}=[\rho_{\mathrm{l}}c_{\mathrm{p}, \mathrm{l}}T_{\mathrm{i}}\mathbf{u}_{\mathrm{i}}-\lambda_{\mathrm{l}}\nabla T
\vert_{\mathrm{l},\mathrm{i}} +\lambda_{\mathrm{s}}\nabla T \vert_{\mathrm{s},\mathrm{i}}] \cdot \mathbf{n}.
\end{equation}}\unskip
where subscript $\mathrm{l}$ and $\mathrm{s}$, stand for liquid and solid phase,
respectively. Using Equation~(\ref{u_interface3D}) for the expression
of $\mathbf{u}_{\mathrm{i}} \cdot \mathbf{n}$ gives:
{\vspace*{3pt}}
{\begin{equation}
\left[\rho_{\mathrm{i}}c_{\mathrm{p}, \mathrm{l}}T_{\mathrm{m}} \left[1-\left(1-\frac{\rho_{\mathrm{s}}}
{\rho_{\mathrm{i}}}\right)\right]-\rho_{\mathrm{s}}c_{\mathrm{p}, \mathrm{s}}T_{\mathrm{i}}\right]\mathbf{v}_{\mathrm{i}}
\cdot \mathbf{n}=[-\lambda_{\mathrm{l}}\nabla T \vert_{\mathrm{l}, \mathrm{i}} +\lambda_{\mathrm{s}}
\nabla T \vert_{\mathrm{s},\mathrm{i}}] \cdot \mathbf{n}.
\end{equation}}\unskip
Using the definition of the reaction enthalpy Equation~(\ref{eqcp}),
where~2 and~1, stands for liquid and solid, subscripts~$\mathrm{l}$ and~$\mathrm{s}$,
respectively we obtain the following relation, which identifies to the
classic Stefan condition, when $T_{\mathrm{i}}=T_{\mathrm{m}}$ and $h_{\mathrm{r}}=L_{\mathrm{m}}$:
{\vspace*{3pt}}
{\begin{equation} \label{Stefancond}
\rho_{\mathrm{s}}h_{\mathrm{r}}\mathbf{v}_{\mathrm{i}} \cdot \mathbf{n}=[-\lambda_{\mathrm{l}}\nabla T
\vert_{\mathrm{l},\mathrm{i}} +\lambda_{\mathrm{s}}\nabla T \vert_{\mathrm{s},\mathrm{i}}] \cdot \mathbf{n}.
\end{equation}}\unskip
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.

\paragraph*{Mass transfer versus heat transfer processes}

The same formalism apply to dissolution/{\ubreak}sublimation and
to the melting of multi-component systems. In both cases, the
combination of Equations~(\ref{fluxchim}), (\ref{cvd3D})
and~(\ref{Stefancond}) are coupled because the fluid physical
properties (densities, diffusion coefficient, the saturation
concentration, the latent heat, the heat capacity,\,\ldots) are
\mbox{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
\mbox{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~\cite{Carpy2019}.

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~\cite{Alkattan, Colombani2007,
Brantley2008}. 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~\cite{AlexiadesSolomonBook, Worster2000}. In some conditions,
the phase change can be controlled by the salt concentration and the
ice melting problem is closer to a dissolution
problem~\cite{woods1992melting}. 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~\cite{Malyarenko2020, McCutchan2022, cenedese2023}.

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~\cite{Pegler2020}. In this case, the rate of phase change is
given by Equations~(\ref{fluxchim}) and~(\ref{cvd3D}). Combining these
two equations leads to the boundary condition for the dissolving\break
interface:
{\begin{equation} \label{EqBCdiss}
\rho_{\mathrm{s}}\mathbf{v}_{\mathrm{i}} = -\alpha(c_{\mathrm{sat}} - c_{\mathrm{i}})\mathbf{n}
=\frac{1}{1-c_{\mathrm{i}}/\rho_{\mathrm{i}}} D(\nabla c \vert_{\mathrm{i}} \cdot \mathbf{n})
\mathbf{n}.
\end{equation}}\unskip
This boundary condition combines the concentration value $c_{\mathrm{i}}$ and
the concentration gradient $\nabla c \vert_{\mathrm{i}}$ at the interface. By
introducing the typical length scale $\delta_{c}$ of the concentration
gradient, i.e.\ $\nabla c \sim (c_b-c_{\mathrm{i}})/\delta_c$ where $c_b$ is the
far field concentration, we can define the Damk\"ohler\break number
{\begin{equation}
\mathit{Da}=(\alpha \delta_{c})/D,
\end{equation}}\unskip
which compares the relative importance of the chemical reaction
(middle-hand side of Equation~(\ref{EqBCdiss})) to the diffusive
transport (right-hand side of Equation~(\ref{EqBCdiss})). In the case
of natural rocks, for which $c_{\mathrm{i}} \ll \rho_{\mathrm{i}}$, the interfacial
concentration is given by $c_{\mathrm{i}} = c_{b}/(1+\mathit{Da}) + c_{\mathrm{sat}}
\mathit{Da}/
(1+\mathit{Da})$~\cite{Philippi2019}. If the chemical reaction rate is the
limiting process at low values of \textit{Da} compared to one, then $c_{\mathrm{i}}
\sim c_{b} + \mathit{Da}\, c_{\mathrm{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_{\mathrm{sat}}$~\cite{nicklasson1984relationship}. In contrast,
large values of \textit{Da} give $c_{\mathrm{i}} \approx c_{\mathrm{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~(\ref{AdvDiffconc})) is thus
solved with the boundary conditions $c_{\mathrm{i}}=c_{\mathrm{sat}}$ and
Equation~(\ref{cvd3D}), 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}_{\mathrm{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~\cite{Chaigne2023}.

\subsection{Models for melting, dissolution and sublimation}
\label{sec3.2}

For an evolving interface undergoing a phase change, the interface
velocity $\mathbf{v}_{\mathrm{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. 

\subsubsection{Melting/solidification of a pure substance}
\label{melting}

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:
\advance\jot2pt
{\begin{eqnarray}
\frac{\partial T}{\partial t} &=& \kappa_{\mathrm{s}}\nabla^2 T \quad
\text{(solid phase)}\label{eqAdvDiffTsol}\Seqnsplit
\frac{\partial T}{\partial t}+\mathbf{u} \cdot \nabla T &=& \kappa_{\mathrm{l}}
\nabla^2 T \quad \text{(liquid phase)} \label{eqAdvDiffTliq}\Seqnsplit
\nabla \cdot \mathbf{u} &=& 0 \quad \text{(liquid phase)}
\label{eqIncomp2}\Seqnsplit
\rho_{\mathrm{l},0}\frac{\partial \mathbf{u}}{\partial t} + \rho_{\mathrm{l},0}
(\mathbf{u} \cdot \nabla) \mathbf{u} &=& -\nabla P + \mu \nabla^2
\mathbf{u} +\rho_{\mathrm{l},0} (1-\beta(T-T_0))\mathbf{g} \quad \text{(liquid
phase)}. \label{eqNS2}
\end{eqnarray}}\unskip
Balance in total energy, mass, and momentum are described by the heat
transport equation (Equations~(\ref{eqAdvDiffTsol})
and~(\ref{eqAdvDiffTliq}) in solid and liquid phases, respectively), by
the continuity equation for an incompressible flow
(Equation~(\ref{eqIncomp2})) and the Navier--Stokes equations
(Equation~(\ref{eqNS2})) for a viscous flow. $\beta$ is the thermal
expansion coefficient (for example $\beta=2.06\times
10^{-4}~\mathrm{K}^{-1}$ for liquid water at 20~{\textdegree}C), $P$ is
the pressure field, $\mathbi{g}$ is the gravity acceleration and $\mu$
the dynamic viscosity of the liquid phase (the kinematic viscosity
$\nu$ is equal to $\mu/\rho$). $\kappa_{\mathrm{s}}$ and $\kappa_{\mathrm{l}}$ are the
thermal diffusivities in the solid and liquid phase, respectively. The
liquid density $\rho_{\mathrm{l},0}$ is assumed to be homogeneous and constant
in the Navier--Stokes equation (Equation~(\ref{eqNS2})) according to
the Boussinesq approximation, where the density is only slightly
modified in the gravitational term. $T_0$ is the temperature for which
$\rho_{\mathrm{l}}(T_0)=\rho_{\mathrm{l},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 $\mathbf{v}_{\mathrm{m}}$:
{\begin{eqnarray}
\mathbf{u} \cdot \mathbf{n} &=& v_{\mathrm{m}}\left(1-\frac{\rho_{\mathrm{s}}}
{\rho_{\mathrm{i}}}\right)\Seqnsplit
\mathbf{u} \wedge \mathbf{n} &=& \mathbf{0}\Seqnsplit
T_{\mathrm{i}} &=& T_{\mathrm{m}}\Seqnsplit
\rho_{\mathrm{s}} L_{\mathrm{m}} \mathbf{v}_{\mathrm{m}} \cdot \mathbf{n} &=& -\lambda_{\mathrm{l}}
\nabla T \vert_{\mathrm{l},\mathrm{i}} \cdot \mathbf{n} +\lambda_{\mathrm{s}}\nabla T
\vert_{\mathrm{s},\mathrm{i}} \cdot \mathbf{n}. \label{StefanCond3}
\end{eqnarray}}\unskip
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. 

\begin{table*}[t!] %tab1
\caption{\label{TabDimensionlessnumbersfusion}List of relevant
dimensionless numbers for melting phenomena}
\begin{tabular}{ccc}
\thead
Melting & Definition & Interpretation\\ 
\endthead
\xmorerows{1}{Stefan number} & 
\xmorerows{1}{$\mathit{St}=L_{\mathrm{m}}/(c_{\mathrm{p}}\,\Delta T)$} &
\parbox[t]{5.8cm}{\raggedright Competition of the latent heat with the
stored thermal energy}\\
&& \parbox[t]{5.8cm}{\raggedright One Stefan number for each
phase}{\vspace*{5pt}}\\ 
\morerows{1}{Prandtl number} & $\mathit{Pr}=\nu/\kappa$ &
\morerows{1}{\parbox[t]{5.8cm}{\raggedright Momentum transport by
viscosity versus heat diffusion}}\\ 
& $\mathit{Pr}\sim 10$ (water) & {\vspace*{7pt}}\\ 
\morerows{1}{Density ratio} & $r_\rho={\rho_{\mathrm{l}}}/{\rho_{\mathrm{s}}}$ &
\morerows{1}{\parbox[t]{5.8cm}{\raggedright Ratio liquid phase density to
solid phase density at the interface}}\\
& $r_\rho \sim 1.09$ (water) & {\vspace*{7pt}}\\
Reynolds number & $\mathit{Re}=({\delta_U\,U})/{\nu}$ &
\parbox[t]{5.8cm}{\raggedright Competition between flow advection and the
viscous diffusion}{\vspace*{10pt}}\\
Thermal P\'eclet number & $\mathit{Pe}_{\mathrm{T}}=({\delta_U\,U})/{\kappa}=\mathit{Pr}\,\mathit{Re}$ &
\parbox[t]{5.8cm}{\raggedright Competition between flow advection and the
temperature diffusion}{\vspace*{10pt}}\\
Thermal Rayleigh number & $\mathit{Ra}=\dfrac{\beta\Delta Tg\delta_{T}^3}
{\kappa\nu}$ & \parbox{5.8cm}{\raggedright Competition between
buoyancy driven advection and diffusion phenomena}{\vspace*{10pt}}\\
Interface motion number & $r_{\mathrm{m},U}=\dfrac{\kappa T_{\mathrm{m}}}{\rho_{\mathrm{s}} L_{\mathrm{m}}
\delta_T U}$ & \parbox{5.8cm}{\raggedright Ratio of the melting
interface velocity to the typical flow velocity}{\vspace*{3pt}}
\botline
\end{tabular} 
\tabnote{$\delta_U$ is the
characteristic length scale of the flow, either the typical scale of
velocity gradient or the thickness of a viscous boundary layer.
$\delta_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.}
\end{table*}

Relevant dimensionless numbers can be introduced to define flow and
melting regimes, see Table~\ref{TabDimensionlessnumbersfusion}. The
density ratio $\rho_{\mathrm{i}}/\rho_{\mathrm{s}}$, 
where $\rho_{\mathrm{i}}$ is taken at the melting
temperature, determines the normal velocity at the melting interface.
For each phase, we can define a Stefan number $\mathit{St}=L_{\mathrm{m}}/ (c_{\mathrm{p}}\Delta
T)$, with $\Delta T$ the typical temperature variation in the problem.
This number compares the latent heat to the stored thermal energy, the
lower the $\mathit{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, $\mathit{Pr}=\nu/\kappa_{\mathrm{l}}=
\mu/(\rho\kappa_{\mathrm{l}})$ is the ratio of kinematics viscosity to thermal
diffusivity. \textit{Pr} for water is about 13.5 at 0~{\textdegree}C and 6.99
(smaller) at 20~{\textdegree}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 $\mathit{Re}=(\delta_U U)/\nu$, where $\delta_U$ is a
characteristic length and the thermal Rayleigh number $\mathit{Ra}_{\mathrm{T}}=(\Delta
\rho g \delta_T^3)/(\mu \kappa_{\mathrm{l}})$ with $\Delta \rho=\rho_0 \beta
\Delta T$, where $\delta_T$ is a characteristic scale of the
temperature field. $\delta_T$ and $\delta_U$ depend on the considered
system and flow conditions. The thermal P\'eclet number $\mathit{Pe}_{\mathrm{T}}=(\delta_U
U)/\kappa_{\mathrm{l}}=\mathit{Pr}\,\mathit{Re}$ compares in Equation~(\ref{eqAdvDiffTliq}) 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 \mbox{convection} (Rayleigh--B\'enard
instability) between two horizontal solid plates maintained at constant
temperatures, is predicted for a Rayleigh number above a critical value
equal to 1707.8~\cite{Chandrasekhar}.

\subsubsection{Dissolution/precipitation of a solid in a solvent}
\label{sec3.2.2}

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:
{\begin{eqnarray}
\nabla \cdot \mathbf{u} &=& 0 \label{eqIncomp}\Seqnsplit
\frac{\partial c}{\partial t}+\mathbf{u} \cdot \nabla c &=& D \nabla^2
c \label{eqAdvDiff}\Seqnsplit
\rho_{\mathrm{l},0}\dfrac{\partial \mathbf{u}}{\partial t}+ \rho_{\mathrm{l},0}
(\mathbf{u} \cdot \nabla) \mathbf{u} &=& -\nabla P + \mu\nabla^2
\mathbf{u} +\rho_{\mathrm{l}}\mathbf{g}. \label{eqNS}
\end{eqnarray}}\unskip
Equations~(\ref{eqIncomp}) and~(\ref{eqNS}) are again the
incompressible Navier--Stokes equation in the Boussinesq approximation,
the fluid density is assumed homogeneous and equal to the mean value
$\rho_{\mathrm{l},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_{\mathrm{l}}=c({\rho_{\mathrm{sat}}-
{\rho_{\mathrm{l},0}}})/({c_{\mathrm{sat}}})+\rho_{\mathrm{l},0}$, where
$\rho_{\mathrm{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_{\mathrm{l}}\mathbf{g}$. The boundary conditions at the solid/liquid
interface are:
{\begin{eqnarray}
\mathbf{u} \cdot \mathbf{n} &=& v_{\mathrm{i}}\left(1-\frac{\rho_{\mathrm{s}}}
{\rho_{\mathrm{i}}}\right)\Seqnsplit
\mathbf{u} \wedge \mathbf{n} &=& \mathbf{0}\Seqnsplit
\rho_{\mathrm{s}}\mathbf{v}_{\mathrm{i}} &=& -\alpha (c_{\mathrm{sat}}-c_{\mathrm{i}})
\mathbf{n} \label{cvd4}\Seqnsplit
\rho_{\mathrm{s}}\mathbf{v}_{\mathrm{i}}\left(1-\frac{c_{\mathrm{i}}}{\rho_{\mathrm{i}}}\right) &=&
D (\nabla c\vert_{\mathrm{i}} \cdot \mathbf{n}) \mathbf{n}, \label{cvd5}
\end{eqnarray}}\unskip
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_{\mathrm{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~\ref{TabDimensionlessnumbersdissolution}. We have already
introduced the Damk\"ohler number $\mathit{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
$\mathit{Sc}=\nu/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 (D T)^{1/2}$, $\delta_v \sim (\nu
T)^{1/2}$ and $D\ll \nu$. In addition, the ratios
$c_{\mathrm{sat}}/\rho_{\mathrm{sat}} \approx c_{\mathrm{sat}}/
\rho_{\mathrm{i}}$ and $\rho_{\mathrm{i}}/\rho_{\mathrm{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
$\mathit{Re}=\delta_U\,U/\nu$ 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 $\mathit{Ra}=[(\rho_{\mathrm{i}}-\rho_b)g
\delta_c^3]/(\mu\,D)$, where $\rho_b$ is the liquid density of the
``bath'' i.e., far from the dissolving interface. Like the thermal
Rayleigh number we introduce in Section~\ref{melting}, the solutal
Rayleigh number determines the condition of appearance of convective
motion by a hydrodynamic instability mechanism~\cite{Sullivan1996,
Philippi2019, Cohen2020}. Finally, the P\'eclet number
$\mathit{Pe}=(\delta_U\,U)/D$ compares advection with diffusion. The P\'eclet
number is directly related to the Reynolds and Schmidt numbers
$\mathit{Pe}=\mathit{Sc}\,\mathit{Re}$.

\begin{table*}[t!] %tab2
\caption{\label{TabDimensionlessnumbersdissolution}List of relevant
dimensionless numbers for dissolution phenomena}
\tabcolsep 4.5pt
\begin{tabular}{ccc}
\thead
Dissolution & Definition & Interpretation\\ 
\endthead
Damkh\"oler number & $\mathit{Da}=\alpha \delta_{c}/D$ &
\parbox[t]{5.3cm}{\raggedright Competition between chemical reaction rate
and diffusive transport}{\vspace*{5pt}}\\

\morerows{2}{Schmidt number} & $\mathit{Sc}=\nu/D$ & \morerows{1}{\parbox[t!]{5.3cm}{\raggedright
Momentum transport by viscosity versus solute diffusion}}\\
& $\mathit{Sc}\sim 1000$ &\\
& (ions in water) & {\vspace*{7pt}}\\

\morerows{2}{Density ratio} & $r_\rho={\rho_{\mathrm{sat}}}/{\rho_{\mathrm{s}}}$ &
\morerows{2}{\parbox[t]{5.3cm}{\raggedright Ratio liquid phase density
(solute $+$ solvant) to solid phase density of a soluble body}}\\
& $r_\rho \sim 0.55$ &\\
& (solid and saturated water) & {\vspace*{8pt}}\\

\morerows{1}{Concentration number} & $r_c={c_{\mathrm{sat}}}/{\rho_{\mathrm{sat}}}$ &
\morerows{1}{\parbox[t]{5.3cm}{\raggedright Ratio of the saturation
concentration to the saturation liquid density. Determine the condition
of weak concentration}}\\
& $r_c \sim 0.2$ (salt) & {\vspace*{31pt}}\\

Reynolds number & $\mathit{Re}=({\delta_U\,U})/{\nu}$ &
\parbox[t]{5.3cm}{\raggedright Competition between flow advection and the
viscous diffusion}{\vspace*{10pt}}\\
P\'eclet number & $\mathit{Pe}=({\delta_U\,U})/{D}=\mathit{Sc}\,\mathit{Re}$ &
\parbox[t]{5.3cm}{\raggedright Competition between transport due to flow
advection and the concentration diffusion.}{\vspace*{10pt}}\\
Solutal Rayleigh number & $\mathit{Ra}=\dfrac{(\rho_{\mathrm{sat}}-\rho_b) g
\delta_{c}^3}{D\mu}$ & \parbox{5.3cm}{\raggedright Competition between
buoyancy driven advection and diffusion phenomena}{\vspace*{10pt}}\\
Interface motion number & $r_{d,U}=\dfrac{D c_{\mathrm{sat}}} {\rho_{\mathrm{s}}
\delta_{c}U}$ & \parbox{5.3cm}{\raggedright Ratio of the dissolving
interface velocity to the typical flow velocity}{\vspace*{3pt}}
\botline
\end{tabular} 
\tabnote{$\delta_U$ is a
characteristic scale of the flow, either the typical scale of the
velocity gradient or the thickness of a viscous boundary layer, and
$\delta_{c}$ is the characteristic scale of the concentration field.}
\end{table*}

\subsubsection{Sublimation/condensation between a solid and a gas
phase} \label{sec3.2.3}

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
exemple) are identical to Equations~(\ref{eqIncomp}) and~(\ref{eqNS}).
The fluid is subject to small temperature variations (constant
$\rho_{\mathrm{l},0}$). The transport of the ablated species within the fluid is
treated as a passive scalar, satisfying an advection/diffusion equation
for the concentration (Equation~(\ref{eqAdvDiff})). The boundary
conditions (\ref{cvd4})--(\ref{cvd5}) are identical to those for
dissolution, with $\alpha$ 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 models~used. As with dissolution, dimensionless numbers for
phase change in sublimation are defined in
Table~\ref{TabDimensionlessnumberssublimation}. The Damkh\"oler number
for sublimation, $\mathit{Da}_{\mathrm{s}} = \alpha \delta_{\mathrm{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_\rho={\rho_{\mathrm{sat}}}/{\rho_{\mathrm{s}}}$ ratio is very small,
indicating a concentration jump at the interface. Temperature and
pressure of the environment control the value of $p_{\mathrm{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_{\mathrm{vap}}}/{p_{\mathrm{sat}}}$,
where $p_{\mathrm{vap}}$ is the vapor pressure of the specie and
$p_{\mathrm{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 \textit{Sc} and 
\textit{Re} gives the \textit{Pe} number,
which compares concentration transport by flow advection and diffusion.
If buoyancy forces are driving the flow, the same Rayleigh number \textit{Ra}
as in the dissolution case can be used to estimate the relative
strength of buoyancy due to concentration/density variations and\break
diffusion.

\begin{table*}[t!] %tab3
\caption{\label{TabDimensionlessnumberssublimation}List of relevant
dimensionless numbers for sublimation phenomena}
\begin{tabular}{ccc}
\thead
Sublimation & Definition & Interpretation\\ 
\endthead
\parbox[t]{3cm}{\centering Damkh\"oler number for sublimation} &
$\mathit{Da}_{\mathrm{s}}=\alpha\delta_{\mathrm{s}}/D$ & 
\parbox[t]{6cm}{\raggedright Competition
between sublimation reaction rate and diffusive
transport}{\vspace*{6pt}}\\

\morerows{2}{Schmidt number} & $\mathit{Sc}=\nu/D$ &
\morerows{2}{\parbox[t]{6cm}{\raggedright Momentum transport by
viscosity versus molecular diffusion}}\\
& $\mathit{Sc}\sim 1$ &\\
& (for gas) & {\vspace*{6pt}}\\

Density ratio & $r_{\rho} = \rho_{\mathrm{sat}}/\rho_{\mathrm{s}}$ &
\parbox[t]{6cm}{\raggedright Ratio gaseous phase density (gas $+$
sublimated specie) to solid phase density of a soluble
body}{\vspace*{6pt}}\\

\morerows{1}{Relative humidity} & $RH = \dfrac{p_{\mathrm{vap}}} {p_{\mathrm{sat}}}$
& \morerows{1}{\parbox[t]{6cm}{\raggedright Determine the conditions
favorable to sublimation of the specie}}\\
& $RH<1$ sublimation & {\vspace*{8pt}}\\

Reynolds number & $\mathit{Re}=({\delta_UU})/{\nu}$ &
\parbox[t]{6cm}{\raggedright Competition between flow advection and the
viscous diffusion}{\vspace*{10pt}}\\
P\'eclet number & $\mathit{Pe}=({\delta_UU})/{D}=\mathit{Sc}\,\mathit{Re}$ &
\parbox[t]{6cm}{\raggedright Competition between flow advection and the
concentration diffusion}{\vspace*{10pt}}\\

Mass Rayleigh number & $\mathit{Ra}=\dfrac{(\rho_{\mathrm{sat}}-\rho_b) g
\delta_{\mathrm{s}}^3}{D\mu}$ & \parbox{6cm}{\raggedright Competition between
buoyancy driven advection and and diffusion phenomena}{\vspace*{10pt}}\\

Interface motion number & $r_{\mathrm{s},U}=\dfrac{D c_{\mathrm{sat}}} {\rho_{\mathrm{s}}
\delta_{\mathrm{s}} U}$ & \parbox{6cm}{\raggedright Ratio of the sublimating
interface velocity to the typical flow velocity}{\vspace*{3pt}}
\botline
\end{tabular} 
\tabnote{$\delta_U$ is a
characteristic scale of the flow, either the typical scale of the
velocity gradient or the thickness of a viscous boundary layer.
$\delta_{\mathrm{s}}$ is the characteristic scale of the concentration field.
\looseness=1}
\end{table*}

\subsubsection{Melting, dissolution and sublimation. Analogies and
differences} \label{sec3.2.4}

{\advance\baselineskip by -.07pt

The dissolution/sublimation velocity is set by the normal to the
surface component of the concentration gradient
(Equation~(\ref{cvd3D})) 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\"ohler
number \textit{Da}, the concentration at interface is close to the saturation
concentration ($c_{\mathrm{i}}=c_{\mathrm{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
$\mathit{St}=L_{\mathrm{m}}/ (Cp \Delta T)$), where the advection--diffusion equation of
the temperature is solved with the boundary conditions $T_{\mathrm{i}}=T_{\mathrm{m}}$.
In that case, Equation~(\ref{Stefancond}) 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~\cite{Meakin2010}, as they arise qualitatively by the same
mechanisms. Nevertheless, some of dimensionless numbers have always
very different orders of magnitude. The Schmidt number $\mathit{Sc}$ is about
1000 for dissolution, whereas the Prandtl number $\mathit{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~\cite{Brantley2008}, 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~\cite{ladd2021reactive} 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.}~\cite{DUTKA2020} find an intermediate
situation ($\mathit{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~\cite{Sullivan1996,
Philippi2019, Cohen2020}. Therefore, most macroscopic dissolution
experiments belong to the diffusion limited regime, i.e.\ large \textit{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~\textit{Da}) or in alkaline conditions
(small~\textit{Da})~\cite{covington2014calcite}.

}

\looseness=-1
Then, according to Equation~(\ref{fluxchim}), the kinetics define a
maximal dissolution velocity $v_{d,\mathrm{max}}=\alpha
c_{\mathrm{sat}}/\rho_{\mathrm{s}}$. Although, the values of $\alpha$ are not
always well determined experimentally, the order of magnitude of the
maximal dissolution velocity, $v_{d,\mathrm{max}}$, is of the order of
$10^{-5}$~m/s, for fast dissolving materials like salt~\cite{Alkattan}.
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_{\mathrm{sat}}/(\rho_{\mathrm{s}} \delta_{c} U)$ for dissolution
and $r_{\mathrm{m},U}=\kappa T_{\mathrm{m}}/(\rho_{\mathrm{s}} L_{\mathrm{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 ${}_{\mathrm{l}}$ (liquid phase) by ${}_{\mathrm{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_{\mathrm{p}}$, $\nu$, $\kappa$, $D$,\,\ldots 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~\cite{chittenden2008}. 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 ($\mathit{Sc}=0.63$ for water vapor in air and $\mathit{Pr}=0.7$ for
air)~\cite{Handbook}. Compressibility effects can also become
significant, modifying the hydrodynamic description.

\subsection{Choice of analog material for an experiment} \label{sec3.3}

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~\ref{TabDimensionlessnumbersfusion}--\ref{TabDimensionlessnumberssublimation}.

For melting experiments by a fluid flow, most of experiments are
performed with ice in contact with water. The melting/solidification
temperature of 0~{\textdegree}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~\cite{Lecoanet2021}. Organic
phase change materials are also used to investigate the interaction
between phase change and melting, like cyclohexane~\cite{Davis1984}
($T_{\mathrm{m}}=6.38$~{\textdegree}C), hexadecane~\cite{noel2022natural}
($T_{\mathrm{m}}=17.7$~{\textdegree}C) or polyethylene glycol
wax~\cite{Leitch2003} ($T_{\mathrm{m}}=19$~{\textdegree}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 $\mathit{Pr}=\nu/\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~\cite{Gau1986}
($T_{\mathrm{m}}=29.7$~{\textdegree}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$_2$O) 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_{\mathrm{sat}}$ and $\alpha$ will have similar value for components
belonging to the same family.

The chemical dissolution/precipitation reactions are:
{$$
(\mathrm{NaCl})_{\mathrm{s}} \leftrightarrows \mathrm{Na}^+ + \mathrm{Cl}^-
$$}\unskip
for salt (sodium chloride) and:
{$$
(\mathrm{CaSO}_4{\cdot}2\mathrm{H}_2\mathrm{O})_{\mathrm{s}}\leftrightarrows
\mathrm{Ca}^{2+} + \mathrm{SO}_4^{2-}+2 \mathrm{H}_2\mathrm{O}
$$}\unskip
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~\cite{Kaufmann2009}.
The chemical reactions for removal or deposition of the calcite are
quite complex with several elementary steps~\cite{Brantley2008,
Vesipa2015}. However, they can be summarized in a single chemical
equation:
{$$
\mathrm{CaCO}_3 + \mathrm{CO}_2 + \mathrm{H}_2\mathrm{O}
\leftrightarrows \mathrm{Ca}^{2+} + 2 \mathrm{HCO}_3^-
$$}\unskip

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~\cite{Brantley2008}. The values of the
relevant experimental parameters for the most common soluble materials
are given in Table~\ref{Tabmatprop} for a temperature of about
20~{\textdegree}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~\cite{Colombani2008}. 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~\cite{Colombani2007}. As mentioned
previously, dissolution of a gypsum sample in a microfluidic cell
allows to measure $\alpha$~\cite{DUTKA2020}. 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~\cite{nicklasson1984relationship}. Sugar (sucrose) and
similar organic compound like glucose or benzoic acid are convenient to
perform dissolution experiments~\cite{Nakouzi, Huang2015, Cohen2020},
because they can be molded to the prescribed shape after thermal
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~\cite{Blumberg1974, Glew1980,
Guerin2020}. 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.

\begin{table*}[t!] %tab4
\caption{\label{Tabmatprop}Characteristic parameters for dissolution in
fresh water at a temperature of about 20~{\textdegree}C}
\fontsize{9.7}{11.5}\selectfont
\begin{tabular}{ccccccc}
\thead
& \parbox[t]{1.54cm}{\centering Solid density $\rho_{\mathrm{s}}$ ($\mathrm{kg}{\cdot}\mathrm{m}^{-3}$)} &
\parbox[t]{2cm}{\centering Saturation concentration $c_{\mathrm{sat}}$ ($\mathrm{kg}{\cdot}\mathrm{m}^{-3}$)}{\vspace*{2.5pt}} &
\parbox[t]{1.7cm}{\centering Saturation density $\rho_{\mathrm{sat}}$ ($\mathrm{kg}{\cdot}\mathrm{m}^{-3}$)} &
\parbox[t]{1.56cm}{\centering Dissolution rate $\alpha$ ($\mathrm{m}{\cdot}\mathrm{s}^{-1}$)} &
\parbox[t]{1.59cm}{\centering Diffusion coefficient D~(min|max) ($\mathrm{m}^2{\cdot}\mathrm{s}^{-1}$)} &
\parbox[t]{1.57cm}{\centering Dynamic viscosity $\mu$~(min|max) ($\mathrm{Pa}{\cdot}\mathrm{s}$)}\\
\endthead
\parbox[t]{1.3cm}{\raggedright Limestone \Lbreak CaCO$_3$}\raisebox{10pt}{} & 2711 &
\0\00.312 & 997.432 & \0\,$3\times 10^{-6}$ & $1.4 \times 10^{-9}$ &
$1.002\times 10^{-3}${\vspace*{5pt}}\\

\parbox[t]{1.3cm}{\raggedright Gypsum \Lbreak CaSO$_4{\cdot}$2H$_2$O} &
2320 & \0\02.04\0 & 999.16\0 & $2.7\times 10^{-6}$ & $1.0\times
10^{-9}$ & $1.002\times 10^{-3}${\vspace*{5pt}}\\

\parbox[t]{1.3cm}{\raggedright Salt \Lbreak
NaCl} & 2170 & 317.5\0\0 & 1200 & $5.0\times 10^{-4}$ & 
\parbox[t]{1.5cm}{\centering $1.5\times 10^{-9}$\Lbreak $1.6\times 10^{-9}$} & 
\parbox[t]{1.5cm}{\centering $1.00\times 10^{-3}$ \Lbreak $1.99\times 10^{-3}$}{\vspace*{8pt}}\\

\parbox[b]{1.3cm}{\raggedright Sugar \Lbreak
C$_{12}$H$_{22}$O$_{11}$} & 1540 & 887 & 1327 & ${\gtrsim}5\times 10^{-6}$\, & \0$5.2\times 10^{-10}$ & $1.00\times 10^{-3}$\\
Sucrose &&&&& \0$7.1\times 10^{-11}$ & 0.212

\botline
\end{tabular} 
\tabnote{The density of
fresh water is $\rho_0=997.12~\mathrm{kg}{\cdot}\mathrm{m}^{-3}$ and its
dynamic viscosity $\mu_0=1.002\times
10^{-3}~\mathrm{Pa}{\cdot}\mathrm{s}$~\cite{Handbook}. Limestone, data from
Kaufmann {et~al.}~\cite{Kaufmann2007} at high carbon dioxide pressure.
For limestone, $c_{\mathrm{sat}}$ is strongly dependent on the carbon
dioxide pressure and on the pH. This value is particularly high,
$c_{\mathrm{sat}} \approx 0.015~\mathrm{kg}{\cdot}\mathrm{m}^{-3}$ under
standard conditions. Gypsum, data from Colombani
{et~al.}~\cite{Colombani2007}. Salt, data from Alkattan
{et~al.}~\cite{Alkattan} and from the Handbook~\cite{Handbook}. Sucrose,
data from Cohen {et~al.}~\cite{Cohen2020}, and
others~\cite{english1950diffusion, gladden1953diffusion} and from the
Handbook~\cite{Handbook}.\vspace*{-2pt}}
\end{table*}

Finally, we now consider the case of sublimation of a solid.
Sublimation technique~\cite{Mendes1991} 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~\cite{sherwood1962},
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~\cite{Oneil89}. 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_{\mathrm{sat}} - p_{\mathrm{vap}}$ determines the sublimation
velocity (see Table~\ref{H2OCO2TAB}). 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
environmental conditions, water ice will always sublime about three
orders of magnitude slower than CO$_2$ ice.

\begin{table*}[t!] %tab5
\fontsize{9.3}{11.3}\selectfont
\caption{\label{H2OCO2TAB}Characteristic parameters for sublimation in
N$_2$ atmosphere}
\tabcolsep=4.5pt
\begin{tabular}{ccccccc}
\thead
& \parbox[t]{1.38cm}{\centering Solid density $\rho_{\mathrm{s}}$ ($\mathrm{kg}{\cdot}\mathrm{m}^{-3}$)} &
\parbox[t]{1.38cm}{\centering Saturation pressure $p_{\mathrm{sat}}$ (Pa)} &
\parbox[t]{2cm}{\centering Vapor pressure $p_{\mathrm{vap}}$ Pa~(min|max)}{\vspace*{2.5pt}} &
\parbox[t]{1.58cm}{\centering Sublimation velocity $v_{\mathrm{s}}$ ($\mathrm{m}{\cdot}\mathrm{s}^{-1}$)} &
\parbox[t]{1.73cm}{\centering Diffusion coefficient~D ($\mathrm{m}^2{\cdot}\mathrm{s}^{-1}$)} &
\parbox[t]{1.57cm}{\centering Kinematic viscosity~$\nu$ ($\mathrm{m}^2{\cdot}\mathrm{s}^{-1}$)}\\
\endthead
\parbox[t]{2cm}{\raggedright Sv\'ea, 
Antartica H$_2$O}\raisebox{10pt}{}
& \0865 & 336 & 100--166 & \0\,$3 \times 10^{-9}$ & $2.9\times 10^{-5}$ & $1.5\times 10^{-5}${\vspace*{5pt}}\\

\parbox[t]{2cm}{\raggedright Nantes, 
wind tunnel CO$_2$} & 1560 & $1.4\times 10^{5}$ & 40 & $1.4\times 10^{-6}$ & $1.5\times 10^{-5}$ & $1.5\times 10^{-5}${\vspace*{3pt}}
\botline
\end{tabular} 
\tabnote{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~\cite{Bintanja1999} with $RH \sim 0.3$, CO$_2$ ice data with
$RH \ll 1$ from Bordiec {et~al.}~\cite{Bordiec2020}.}
\end{table*}

\section{Dissolution, melting and sublimating patterns in experiments}
\label{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~\ref{equations} 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. 

\subsection{Buoyancy driven flow and associated bedforms} 
\label{flow1}

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)~\cite{Cohen2016,
DaviesWykes2018, Cohen2020, Pegler2020, Huang2020}, salt (NaCl), or
plaster~\cite{Cohen2020}, or the melting of an icy block in a tank of
water~\cite{Weady2022}.

\subsubsection{Stable density stratification}

The dissolution of axisymmetric sugar bodies in quiescent water has
been studied experimentally by Pegler and Davies
Wykes~\cite{Pegler2020} and by Huang {et~al.}~\cite{Huang2020} for the
case of stable density stratification (Figure~\ref{fig:ExpStable}). 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~\cite{Huang2020} whereas an initial cone blunts with
time~\cite{Pegler2020}; the erosion dynamics depends on the initial
shape~\cite{Pegler2020}. The erosion dynamics and shapes are well
described by a boundary layer model derived by
Acrivos~\cite{Acrivos1960}, in the limit of high Schmidt number, i.e.,
the concentration boundary layer is within the viscous
layer~\cite{Pegler2020, Pegler2021}. 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~\cite{Pegler2021}. For simple shapes, self similar solutions exist
but no asymptotic, stationary shape~\cite{Pegler2020, Pegler2021}.
Huang {et~al.}\ also interestingly ``speculate that this mechanism
contributes to the formation of pinnacles in nature'', such as tsingys
or stone forests~\cite{Huang2020}, the \textit{sharp patterns}
described in Section~\ref{Penitents}, Figure~\ref{fig:naturalsharp}(b).
Weady {et~al.}\ studied the melting of initial vertical cylinder in a
tank of water at different temperature $T_{\infty}$~\cite{Weady2022}.
These cylinders sharpen to form upward pinnacles when $T_{\infty}
\gtrsim 7$~{\textdegree}C and downward pinnacles when when $T_{\infty}
\lesssim 5$~{\textdegree}C, because water density is maximum around
4~{\textdegree}C~\cite{Weady2022}. 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~\cite{Weady2022}.

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig:ExpStable}Left: dissolution of a candy cone (from
Pegler and Davies Wykes~\cite{Pegler2020}). 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.}~\cite{Huang2020}), under the CC BY 4.0 license.}
\end{figure}

\subsubsection{Unstable density stratification} \label{Unstable}

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\"urr described the
appearance of ``convection stripes'' on the surface of immersed soluble
crystals in 1904~\cite{Schurr1905, SchurrThesis}. The appearance of
patterns on the surface of solids under the effect of convection flow
generated by their dissolution~\cite{Grijseels, dacheng1988} or
melting~\cite{Kerr1994, yen1968onset} 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\'enard 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~\cite{Sullivan1995, Liu1996, Sullivan1996}.

This can be explained with the schematic in Figure~\ref{fig:panaches}.
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\'enard instability happens,
analogous to the thermal Rayleigh--B\'enard instability. The
instability threshold is given by a criterion involving the solutal
Rayleigh number:
{\begin{equation} \label{eq:rayleigh}
\frac{\Delta\rho g\delta_{c}^3}{\rho_{b} \nu D}=\mathit{Ra}_{c},
\end{equation}}\unskip
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
$\mathit{Ra}_{c}$ a critical value of the solutal Rayleigh number. Using the
values in Table~\ref{Tabmatprop} and taking $\mathit{Ra}_{c}=1101$, which is the
theoretical value for the onset of the Rayleigh--B\'enard instability
with mixed slip and no-slip boundary conditions~\cite{Chandrasekhar},
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
$\mathit{Ra}_{c}\simeq100$~\cite{Sullivan1995, Sullivan1996}, 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~\cite{Sullivan1995, Cohen2020}. In
Figure~\ref{fig:panaches}, they are visualised using the shadowgraph
technique under a block of plaster suspended in water at an angle of
20{\textdegree} 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.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig:panaches}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\'enard instability and its imprint on the block
(from Cohen {et~al.}~\cite{Cohen2020}, with the permission of the
American Physical Society).}
\end{figure}

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~\cite{Cohen2020}. Under a horizontal surface, this leads to
the appearance of a pattern of cavities surrounded by sharp crests
reminiscent of scallops~\cite{Sullivan1996, Liu1996, Chaigne2023}, 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~\ref{fig:panaches}. 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~\cite{Cohen2016, Cohen2020}. It can be seen in
Figure~\ref{fig:platreMartin}, 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~\ref{fig:caramelCaroline}. 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~\cite{Cohen2016, Cohen2020}. The
same succession of patterns can be seen on salt blocks.
Figure~\ref{fig:selCaroline} 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~\ref{fig:selCaroline}). Thus, the dissolution rate is greater
immediately downstream of a crest than upstream, causing it to
propagate upstream~\cite{Cohen2020}. 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~\cite{Cohen2016, Cohen2020}, 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~\ref{fig:naturalscallops}(b)). Similar pattern
formation is indeed expected when melting induces thermal or solutal
convection. Cohen {et~al.}~\cite{Cohen2016} 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.

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig:platreMartin}Patterning of the bottom surface of a
plaster block immersed in tap water at an angle of 20{\textdegree} from
the horizontal. The water was at room temperature, between
17~{\textdegree}C and 21~{\textdegree}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.}
\end{figure}

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig:caramelCaroline}Patterning of the bottom surface of
a caramel block immersed in tap water (20~{\textdegree}C) at an angle
of 60{\textdegree} from the horizontal (from Cohen
{et~al.}~\cite{Cohen2020}, with the permission of the American Physical
Society).}
\end{figure}

\begin{figure}
\includegraphics{fig12}
\caption{\label{fig:selCaroline}Left: Bottom surface of a salt block
after immersion for 1~h in tap water (20~{\textdegree}C) at an angle of
60{\textdegree} 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.}~\cite{Cohen2020}, with the permission of
the American Physical Society.}
\end{figure}

On a much larger scale, the average thickness of the boundary layer is
fixed by the balance between diffusion and buoyancy
(Equation~(\ref{eq:rayleigh})). 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~\cite{DaviesWykes2018}. 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~\cite{Cohen2020}. This dependency
of angle only becomes significant 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~\cite{Gechter2008}. 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~\cite{sharma2022}.

\subsection{Forced surface flows and associated bedforms} 
\label{flow2}

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 \uppi$, then deep
flow condition should be considered, and if not, shallow water.

\subsubsection{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~\cite{Allen1970, Goodchild1972,
Villien2001, villien2005, wang2016time} report after about a hundred
hours the formation of a scallop pattern, similar to \textit{the
polygonal depressions} observed in the nature and previously described
in Section~\ref{Polygonal}, Figure~\ref{fig:naturalscallops}(c--d). In
experiments, the typical scallop size and the surface roughness
increase in time. In particular Allen~\cite{Allen1970, Allen1971}
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~\ref{fig:ripple_expe}(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 the \textit{crested patterns parallel to the flow
direction} described in Section~\ref{rillenkarrens}, 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.

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig:ripple_expe}Deep flow conditions experiments.
(a)~Dissolution waves from Blumberg and Curl~\cite{Blumberg1974},
copyright Cambridge University Press 1974 and (b)~sublimation waves
from Bordiec {et~al.}~\cite{Bordiec2020}, with permission from Elsevier.
(c)~Scallops from Allen~\cite{Allen1970} 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~\cite{Blumberg1974}, copyright Cambridge
University Press 1974.}
\end{figure} 

Regularly spaced crests, parallel to each other and oriented
perpendicular to the turbulent overlaying flow (Section~\ref{wavy},
Figure~\ref{fig:naturalripples}) 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~\cite{Blumberg1974} experiment, dissolution waves are reproduced
on a plaster block as an analogous material of soluble rock
(Figure~\ref{fig:ripple_expe}b). Their morphodynamic study revealed
wavelengths, amplitudes, propagation angles and dissolution rates
defined in Figure~\ref{fig:ripple_expe}(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~\ref{TabDimensionlessnumbersdissolution}). By comparing the flow
velocity with the ablation rate{\unskip\break} (typically
$1~\mathrm{m}{\cdot}\mathrm{s}^{-1}$ and
$0.5~\upmu\mathrm{mm}{\cdot}\mathrm{s}^{-1}$), these experiments are well
placed in a quasi-static interface evolution, i.e.\ $r_{d,U} \ll 1$ (see
Table~\ref{TabDimensionlessnumbersdissolution}). The boundary layer of
this flow $\delta_U \simeq7.6$~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}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$) in the direction of flow.

\looseness=-1
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~\ref{TabDimensionlessnumbersdissolution}), based on the
wavelength over the viscous length (i.e.\ the kinematic viscosity of the
fluid $\nu$ over the friction velocity~$u_*$). The friction velocity
used for the velocity scale~\cite{schlichting1979} is based on the wall
shear stress such as $u_* = (\tau_{\mathrm{w}}/\rho_{\mathrm{f}})^{1/2}$ with~$\rho_{\mathrm{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~\ref{fig:ripple_expe}(b)) and sublimation waves observed in the
field (Figure~\ref{fig:naturalripples}(a)) and the analog sublimation
experiment on CO$_2$ ice performed in an atmospheric wind
tunnel~\cite{Bordiec2020} (Figure~\ref{fig:ripple_expe}(c)). This
theoretical study~\cite{Carpy2023} is based on a 2D linear stability
analysis of the coupled system between turbulent flow and mass
transfer, following on from pioneer work~\cite{Hanratty1981,
Claudin2017, Bordiec2020}. 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~$\lambda$, growth rate~$\sigma$, angular
frequency~$\omega$). The growth rate and angular frequency are
functions of the wavelength-based Reynolds number ($\mathit{Re}_{\lambda} =
{\lambda u_*}/{2 \uppi \nu} = ({k\nu}/{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\break laws:
{\advance\jot by 5pt\begin{eqnarray}
&\displaystyle \lambda_c = \frac{\nu}{u_*} \frac{2\uppi}{k_{c,
\mathrm{ref}}^+}\mathit{Sc}^{1/5}\Seqnsplit
&\displaystyle v_c = \frac{\omega_{c,\mathrm{ref}}^+}{k_{c,
\mathrm{ref}}^+} \frac{q_0}{\rho_{\mathrm{s}}}\Seqnsplit
&\displaystyle t_c = \frac{1}{\sigma_{c,\mathrm{ref}}^+} \frac{\rho_{\mathrm{s}}}
{q_0} \frac{\nu}{u_*}\mathit{Sc}^{1/6}
\end{eqnarray}}\unskip
with $k_{c,\mathrm{ref}}^+ = 6.7 \times 10^{-3}$, $\omega_{c,
\mathrm{ref}}^+ = 0.013 \times 10^{-3}$, $\sigma_{c, \mathrm{ref}}^+=
8.9 \times 10^{-3}$ where reference values are given for $\mathit{Sc} =
1$~\cite{Carpy2023}. 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~\cite{AshtonKennedy1972, Ashton1972, Hsu1979,
Gilpin1980}. 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.}~\cite{Gilpin1980}, recent experimental work form Bushuk
{et~al.}~\cite{Bushuk2019} 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 $\mathit{Re}_{\lambda}$. 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~\ref{fig:ripple_expe}(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.

\subsubsection{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~\cite{ladd2021reactive} for more information on the subject.
In transport controlled regime, the fracture aperture in gypsum
decreases with the distance~\cite{li2021transport}. However, especially
in experiments with NaCl the general shape and small scale pattern
formation occurs in presence of buoyancy
instabilities~\cite{Oltean2013, hu2023role, guo2023dissolution}, which
have been already discussed in Section~\ref{flow1}. 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 \textit{Karren} patterns (see
Section~\ref{rillenkarrens}) only the generation of
\textit{Rillenkarren} has been reported in quantitative experiments.
Around 1980, Glew and Ford~\cite{Glew1980} applied artificial rain with
an intensity about 35~mm/h upon molded plaster blocks of various
inclination, leading to film flows about $100~\upmu$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~\upmu$m.
The obtained morphology resembles thus well to natural
\textit{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~{\textdegree}C and
45~{\textdegree}C)~\cite{Slabe2016}. More recently, Guerin
{et~al.}~\cite{Guerin2020} 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~\upmu$m, nearly parallel rills of similar width appear after few
minutes on pink salt plates and after few hours on plaster plates (see
Figure~\ref{ImageRillenkarrenCras2b}). 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~\cite{bertagni2021hydrodynamic, Djambov_Gallaire_2023}
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 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 \textit{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 \textit{Rillenkarren} and their
long term evolution. 

\begin{figure}
\includegraphics{fig14}
\caption{\label{ImageRillenkarrenCras2b}(a)~Side view of the
experimental setup of Guerin {et~al.}~\cite{Guerin2020} 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{\textdegree}, 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.}
\end{figure} 

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~\cite{Pitcher2019}. The flow in
these channels may increase the heat exchange and thus their depths by
enhancing melting. Mantelli {et~al.}~\cite{Mantelli2015} proposed an
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 \textit{Rinnenkarren} for dissolution. An analog
of \textit{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.

{\vspace*{2pt}}

\subsection{Pattern formation stimulated by radiation}
\label{SectionCapExp}

{\vspace*{2pt}}

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~\ref{fig:naturalsharp}(a)), glacier tables or dirt cones
(Figure~\ref{fig:CappedPattern}).

Penitentes have been reproduced in a laboratory scale experiment by
Bergeron {et~al.}~\cite{Bergeron_2006}. 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$~{\textdegree}C, penitentes appeared after two hours of
irradiation with a characteristic of a few millimeters
(Figure~\ref{fig:PenitentesExp}). In melting conditions, above
${-}4$~{\textdegree}C, no structures were observed. In sublimation
condition, under ${-}35$~{\textdegree}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~{\textdegree}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)~\cite{Betterton_2001, Bergeron_2006}. 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~\cite{Bergeron_2006}.

\begin{figure}
{\vspace*{4pt}}
\includegraphics{fig15}
{\vspace*{4pt}}
\caption{\label{fig:PenitentesExp}\textit{Penitentes} in a laboratory
scale experiment (from Bergeron {et~al.}~\cite{Bergeron_2006},
reproduced by permission of the American Physical Society).}
{\vspace*{4pt}}
\end{figure}

Pattern formation is generally ascribed to a geometric coupling between
light absorption and topography: light absorption is greater at troughs
than at peaks~\cite{Lliboitry1954, Betterton_2001, Bartels_2012,
Cathles_2014}. 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~\cite{Claudin2015, Mullins1964}.
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~\cite{Claudin2015}. 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.}~\cite{Bergeron_2006}. Cohen {et~al.}\ 
proposed that this length scale could be set by a solutal
Rayleigh--B\'enard instability, because density of water vapor is less
than density of air~\cite{Cohen2020}. This instability in the
buoyancy-driven vapor flow would prevail over other mechanisms in
initiating the pattern and selecting the wavelength{\unskip\break} 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~\ref{fig:CappedPattern}). Here we review recent laboratory
experimental works that have addressed the formation of glacier tables,
dirt cones, and zen stone described in Section~\ref{NaturalCovered},
Figure~\ref{fig:CappedPattern}.

H\'enot {et~al.}\ have studied the formation of glacier tables and, in
particular, the conditions for their formation~\cite{Henot_2021}. 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~{\textdegree}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~\ref{fig:ExpTable}). 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~\cite{Henot_2022}.

\begin{figure}
{\vspace*{4pt}}
\includegraphics{fig16}
{\vspace*{4pt}}
\caption{\label{fig:ExpTable}Ice underlying cylindrical caps melts
faster or slower than the surrounding ice sheet depending on the cap
material (from H\'enot {et~al.}~\cite{Henot_2021}, reproduced by
permission of the American Physical Society).}
\end{figure}

H\'enot {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~\cite{Henot_2023}.
Differential melting resulted in the formation of an ice cone covered
by a granular layer. The authors explain the formation of dirt cones
(Section~\ref{NaturalCovered}, Figure~\ref{fig:CappedPattern}), not
directly by the dependence of the melting rate on the debris
thickness~\cite{Evatt_2015}, 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~\cite{Drewry_1972, Henot_2023}.

Taberlet and Plihon have studied the formation of zen stones by ice
(differential) sublimation in a lyophilizer~\cite{Taberlet_2021}. 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~\cite{Taberlet_2021}. \looseness=1

\section{Discussion and conclusions} \label{discussion}
\subsection{Solid bedforms as imprint of flow structures}
\label{sec5.1}

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~\cite{Gilpin1980}, 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~\cite{Thomas1979}. 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.

\subsection{Effect of non-linearities on the resulting bedform}
\label{nonlin}

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~\cite{Guerin2020} and
interface evolution models~\cite{Chaigne2023}. 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~\cite{Lange1959} has
been recently revisited by Chaigne {et~al.}~\cite{Chaigne2023}. 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~\cite{Claudin2017} in turbulent flows or for Rillenkarren by
run-off~\cite{Guerin2020, bertagni2021hydrodynamic,
Djambov_Gallaire_2023}, 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~\cite{Fauve2005}. 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\'enard convection~\cite{Sullivan1996,
Cohen2016, Cohen2020}, 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~\cite{Cohen2020} or melting
patterns~\cite{Favier2019}. 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 \textit{Rinnenkarren}~\cite{Gines2009} on
limestone pavements or by melting like the ice streams on
Glacier~\cite{Mantelli2015}.

The nonlinear study of ablation pattern in nature could constitute the
next step of future researches.

\subsection{On the use of scaling laws for planetary applications}
\label{sec5.3}

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~\cite{Yokokawa2016}. For cyclic
steps (Section~\ref{steps}, Figure~\ref{fig:naturalsteps}), 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~\cite{Gilpin1980, Bushuk2019}. 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.}~\cite{Carpy2023}, 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 (\ref{nonlin}). In the case of \textit{penitentes}
(Section~\ref{Penitents}, Figure~\ref{fig:naturalsharp}(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.}~\cite{Cohen2020}.

\subsection{Concluding remarks} \label{sec5.4}

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 \mbox{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.

\section*{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.

\section*{Funding}

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

\CDRGrant[FRC]{SHERPAS ANR-23-CE49-0006}
\CDRGrant[FRC]{PhysErosion ANR-22-CE30-0017}

\back{}

\hyphenpenalty 40
\exhyphenpenalty 50
\fussy

\printbibliography
\refinput{crphys20240562-reference.tex}

\end{document}
