11institutetext:  Center for Combustion Energy, Key laboratory for Thermal
Science and Power Engineering of Ministry of Education, Department of Energy
and Power Engineering, Tsinghua University, Beijing 100084, China
Univ. Lille, Unité de Mécanique de Lille J. Boussinesq - UML - ULR 7512,
F-59000 Lille, France
Department of Engineering Mechanics, School of Aerospace Engineering, Tsinghua
University, Beijing 100084, China

Rayleigh-Bénard convection, Multiphase flows, Freezing, Melting

# Equilibrium states of the ice-water front
in a differentially heated rectangular cell

Ziqi Wang 11 Enrico Calzavarini<EMAIL_ADDRESS>22 Chao
<EMAIL_ADDRESS>1133112233

###### Abstract

We study the conductive and convective states of phase-change of pure water in
a rectangular container where two opposite walls are kept respectively at
temperatures below and above the freezing point and all the other boundaries
are thermally insulating. The global ice content at the equilibrium and the
corresponding shape of the ice-water interface are examined, extending the
available experimental measurements and numerical simulations. We first
address the effect of the initial condition, either fully liquid or fully
frozen, on the system evolution. Secondly, we explore the influence of the
aspect ratio of the cell, both in the configurations where the background
thermal-gradient is antiparallel to the gravity, namely the Rayleigh-Bénard
(RB) setting, and when they are perpendicular, i.e., vertical convection (VC).
We find that for a set of well-identified conditions the system in the RB
configuration displays multiple equilibrium states, either conductive rather
than convective, or convective but with different ice front patterns. The
shape of the ice front appears to be always determined by the large scale
circulation in the system. In RB, the precise shape depends on the degree of
lateral confinement. In the VC case the ice front morphology is more robust,
due to the presence of two vertically stacked counter-rotating convective
rolls for all the studied cell aspect-ratios.

###### Abstract

Supplementary Information for the main article

###### pacs:

nn.mm.xx

Introduction. – Convective liquids undergoing melting or freezing give rise to
a rich phenomenology of flow patterns and solid-phase morphologies. This bears
great relevance in the geophysical domain, e.g., in volcanology for the
understanding of the solidification of magma chambers [1, 2, 3], in
planetology for the study of magma oceans [4], or in geomorphology and
glaciology for glacier dynamics and their induced erosion [5], and in marine
sciences for the prediction of arctic sea-ice annual cycles [6]. In the
technological context, convection driven phase-change has a key role in
metallurgy [7], in the purification of substances [8] and in the storage of
thermal energy through phase-change materials [9]. Depending on the
application, the fluid can be either very complex in composition and rheology
like a heterogenous molten rock magma, or simpler like muddy-water or sea-
water, or nearly ideal as for purified materials employed in chemical
industrial processes. However, probably the most common instance is the one
where the fluid is relatively clean water. This is a rather intriguing case
due to the peculiar non-monotonous buoyancy force intensity in water above the
freezing point. A large variety of situations of thermal or mechanical driving
mechanisms are encountered. The fluid motion can be the result of natural
convection due to a cooling process, as in the case of volcanic magma, or it
can be steadily driven by localized heat sources, such as a hot/cold boundary,
or by distributed ones, as in the case of internal heating by radioactive
decay (common in rock formation) or by absorption of solar radiation as it
happens for water in glaciers or in the oceans [10, 11]. Finally, the thermal
convection can also be maintained by a mechanical driving of the fluid, as it
happens below the sea-ice and around icebergs [12]. From a physicist’s
perspective, a relevant question is how the convecting fluid flow, which can
be either laminar or turbulent depending on the driving intensity, can
determine the overall shape of ice. Is there just a single general mechanism?
Can one envisage a phenomenological model, i.e. a predictive argument for the
shape of the ice-water front that does not need to take into account the
exceptional complexity of the thermodynamics and fluid dynamics equations
involved in its description?

A way to reply to these questions is to envisage a sufficiently simplified and
well-controllable system that retains part of the physical complexity and rich
phenomenology observed in nature. A well studied and commonly used convective
setup in fluid-dynamics is the Rayleigh-Bénrad (RB) system, a fluid filled
container heated at the bottom surface (of temperature $T_{b}$) and cooled at
the top (of temperature $T_{t}$) and thermally insulated on the lateral sides.
The thermal instability, the onset of convection, the flow bifurcations, the
self-organization of the system in distinct parts, the turbulent regime and
its scaling laws, have been analyzed in great detail [13, 14, 15, 16, 17, 18].
This makes RB an ideal system on top of which the complexity of phase-change
phenomena and its coupling with fluid convection can be added to. Previous
studies of melting in RB have highlighted the complexity of the solid-liquid
interface topography, and its dependencies on the strength of convection
(parametrized by the Rayleigh number), on material properties (Stefan and
Prandtl numbers), on the system dimensionality and geometry[4, 20, 21, 12,
23]. Furthermore, an intriguing phenomenon of bistability depending on whether
the system has been initialized as fully liquid or fully frozen has been
identified [24, 12]. However, most of the above mentioned studies have focused
on the case of pure materials and small thermal gaps so that the resulting
buoyancy is linear with respect to temperature (Boussinesq approximation). For
water, the density anomaly, i.e. the local mass density maximum at around
$T_{c}\sim 4^{\circ}$C at atmospheric pressure, leads to a non-linear buoyancy
force, which plays a great role in the icing dynamics. Recently, [3] combined
experiments and numerical simulations to understand the hydrodynamical
mechanisms that control the global extent and the shape of the ice in a
freezing RB rectangular quasi-two-dimensional cell. They identified four
distinct regimes of coupling between thermal stratification, buoyancy forces,
and phase-change [3] (see Table. 1 for a brief summary). Other research
efforts have explored the effect of the RB system inclination on the overall
icing process [2].

In this work, we aim at extending further the available experimental
measurements and numerical simulations of water melting and icing in an RB
cell, in particular by testing the existence of the bistability in a water-
filled RB cell and at exploring the container aspect ratio dependence and its
role with respect to the system inclination, either straight or rotated by
90∘.

Regime | 1 | 2 | 3 | 4
---|---|---|---|---
Stratification type | SSD | SSD(top) + USD(bottom) | SSD(top) + USC(bottom) | SSD(top) +USC(bottom)
Equilibrium state | Diffusive | Diffusive | Convective | Convective
Ice front shape | Flat | Flat | Flat | Deformed
Table 1: Four regimes of water-freezing RB system for increasing temperature
differences (left to right) as identified in [2]: in the second line, the
first two letters of the acronyms specify the stratification type, either
stably stratified (SS, with temperature ranging from $T_{0}$ to $T_{c}$) or
unstably stratified (US, with temperature ranging from $T_{c}$ to $T_{b}$);
the third letter specifies the mode of heat transport, either diffusion (D) or
convection (C); top and bottom denote the positions of the layers in the RB
system.  Figure 1: Bistability phenomenon. (a) Temporal evolution of the ice
thickness, $H_{i}$, of the freezing process (thick line) and the melting
process (dashed line). (b) Comparison of the simulation (symbols) and the
theoretical model (dashed and solid lines): $H_{i}$ as a function of $T_{b}$.
The blue, yellow, and red shaded areas correspond to the conductive, bisable
and convective equilibrium regimes, respectively. The dashed line denotes one
of the two branches in the bistable regime which predicts the convective
equilibrium state. (c), (d) and (e): panels I-III display the evolution of the
instantaneous temperature field of the freezing process, with the bottom plate
temperature $T_{b}=4.5,5.5,10^{\circ}$C, respectively; panels IV-VI display
the evolution of the instantaneous temperature field of the melting process,
for the same bottom plate temperatures. The visualizations of the temperature
field in the water region are with colors and $T_{0}$ (black line) and $T_{c}$
(gray line) isotherms, along with the velocity vectors when fluid convection
occurs.  Figure 2:  Mechanism of bistability: (a) Sketch of the one-
dimensional model: three plane layers are connected in series. Note that if
$T_{b}<T_{c}$ the unstably stratified layer is absent; The computed heat-
fluxes across the cell as a function of the ice thickness, $h_{i}$, in the
purely conductive regime (b), the bistable regime (c) and the convective one
(d).  Figure 3: Aspect ratio dependence. (a) The temperature field at the
equilibrium state of the RB case. Panels I-I, I-II, I-III, and I-IV display
the experimental results with $\Gamma=0.5,1,2,4$, respectively. Panels II-I,
II-II, II-III, II-IV, and II-V display the simulation results, with
$\Gamma=0.25,0.5,1,2,$ and 4, respectively. Panels III-I and III-II displays
the simulation results for $\Gamma=2$ and 4, same conditions as II-IV and II-V
but with different ice front shapes. (b) The temperature field at the
equilibrium state of the VC case. Panels I-I, II-I, and III-I display the
experimental results with $\Gamma=0.5,1,2$, respectively. Panels I-II, II-II,
III-II, IV-II, and I-III display the simulation results, with
$\Gamma=0.25,0.5,1,2,$ and 4, respectively. The visualizations of the
temperature field in (a-b) display $T_{0}$ (black line) and $T_{c}$ (red line)
isotherms. (c) $H_{i}$ as a function of $\Gamma$ for RB and VC cases from
simulations (S) and experiments (E). The shaded area shows the spatial
variation of the ice front (from minima to maxima of $h_{i}(x,t\to\infty)/h$)
in the RB case for experiments (green) and simulations (red).  Figure 4: Ice
front shape boundary layer based model. (a) Ice front morphology at the
statistical equilibrium in VC as extracted from the numerical simulation with
aspect ratio, $\Gamma=$ 0.25 (light blue line), 0.5 (red line), 1 (green
line), 2 (dark blue line), 4 (orange line). The system is cooled at the left
side and heated at the right side. (b) Sketch of the phenomenological BL
model: the ice front at $T_{0}=0^{\circ}$C (black line), the $T_{m}$ isotherm
which is the outer boundary of the thermal boundary (red dashed line), the
$T_{c}$ isotherm (green dashed line). The angle between the tangential
direction of the ice front and the x-direction is $\gamma$; the thickness
normal to the ice front (the green thick line in the ice) is
$\frac{h_{i}(x)}{cos\gamma}$ (with $cos\gamma=\frac{dS}{dx}$). (c)-(g)
Comparison of ice front morphology between experiments (blue shaded area),
simulations (green line) and the model (red line) for $\Gamma=0.25$ (c); 0.5
(d); 1 (e); 2 (f); 4(g).

Experimental setup and Numerical Methods. – The experiments are conducted in a
battery of standard RB cells with four different geometrical proportions of
the container. The details of the setup are documented in the Supplementary
Information (SI). Here we only give a brief description. The cell aspect
ratio, $\Gamma=l_{x}/h$, is here defined as the quotient between the largest
width of the isothermal plates ($l_{x}$) and the distance between the
isothermal plates (also called height) $h$. We choose $l_{x}=24$ cm and
$h=12,24,48,96$ cm leading to $\Gamma=0.5,1,2,4$ respectively. The third
spatial dimension, $l_{z}=6$ cm, is kept fixed and it is smaller than the
others, a fact that allows to designating the system as quasi-two-dimensional.
We note that $h$ stands for the vertical direction, parallel to the
gravitational acceleration, in the RB configuration, while it is horizontal in
the vertical convection (VC) setting, when the same experimental cells are
tilted by $90$ degrees. In the present study, we only perform solidification
experiments, i.e. experiments whose initial condition is liquid water with
uniform initial temperature same as the hot plate temperature $T_{b}$. Given
the difficulty of manufacturing a cell-sized homogeneous ice block, we could
not perform (well controlled) melting experiments. For this reason, in the
present study we also rely on the results of numerical simulations. We employ
the CH4-project code [1] which is based on a Lattice-Boltzmann equation
algorithm for the description of fluid and temperature dynamics, and on the
enthalpy method for phase-change. Both the numerical methods and the specific
algorithm implementation have been extensively validated, particularly against
RB water freezing experiments in [2, 3]. We use the non-monotonic relationship
of the water density as a function of the temperature,
$\rho=\rho_{c}(1-\alpha^{*}|T_{b}-T_{c}|^{q})$, with
$\rho_{c}=999.972~{}kg/m^{3}$ the maximum density corresponding to
$T_{c}\approx 4^{\circ}$C, $\alpha^{*}=9.30\times 10^{-6}(K^{-q})$, and
$q=1.895$ [5]. We conduct two-dimensional simulations, which have been proved
to match well the experiments for the present quasi-2D cells [2, 3]. The
initial conditions encompass both the fully liquid and fully frozen cases,
i.e., freezing and melting numerical experiments. Throughout the
investigation, the cold plate is fixed both in the experiments and the
simulations at temperature $T_{t}=-10^{\circ}$C. The effects of changing
$T_{t}$ are qualitatively predictable (more details on $T_{t}$ effects can be
found in SI). Both in the experiments and in the simulations we monitor the
evolution of the ice-water interface $h_{i}(x,t)$ ($i$ stand for ice), and its
mean instantaneous value profile
$\overline{h_{i}}(t)=l_{x}^{-1}\int_{0}^{l_{x}}h_{i}(x,t)dx$. Such quantities
are more conveniently analyzed in dimensionless terms, hence
$H_{i}=\overline{h_{i}}/h$ represents the fraction of the cell that is
occupied by the ice ($H_{i}=1$ means complete solidification).

Results: bistability. – We first study whether the melting and freezing system
configurations evolve into the same equilibrium state. Given the fact that the
system in the RB arrangement can display different regimes as the intensity of
the thermal driving increases (Table. 1), we choose four typical cases of the
thermal driving condition from each of the existing four different regimes,
i.e., $T_{b}=3,4.5,5.5,10^{\circ}$C, and perform long-term numerical
simulations for melting (initial condition is $H_{i}=1$) and freezing (initial
condition is $H_{i}=0$). Note that the duration of the simulations is of the
order of 10 days in physical units. Fig. 1(a) shows the normalized global ice
thickness, $H_{i}$, as a function of time under different conditions of
thermal driving in the RB setup. The thick (dashed) lines represent the
freezing (melting) case. It is observed that melting and freezing reach the
same equilibrium state for the cases of $T_{b}=3,4.5,10^{\circ}$C; while for
$T_{b}=5.5^{\circ}$C, melting and freezing reach a conductive and a convective
equilibrium respectively. Hence, in the freezing experiments the global amount
of ice at the equilibrium is smaller. To better understand this phenomenon, we
show the evolution of the temperature field (see Fig. 1(c)-(e), we do not show
the $T_{b}=3^{\circ}$C case because when $T_{b}<T_{c}$ (Regime-1), the fluid
region is trivially stably-stratified throughout the entire simulation). At
$T_{b}=4.5^{\circ}$C (Regime-2), though the stably and unstably-stratified
layers coexist, the effective Rayleigh number (which is based on the unstably-
stratified layer and only defined when $T_{b}>T_{c}$),
$Ra_{e}=\frac{(\Delta\rho/\rho_{\text{c}})g(h_{c})^{3}}{\nu\kappa}=\frac{g\alpha^{*}|T_{\text{b}}-T_{\text{c}}|^{q}(h_{c})^{3}}{\nu\kappa}$
(with $h_{c}$ being the averaged height of the $T_{c}$ isotherm, $g$ the
gravitational acceleration, $\nu$ the kinematic viscosity, and $\kappa$ the
thermal diffusivity), decreases from $Ra_{e}\sim 10^{3}$ as the ice thickness
increases in the freezing case (Fig. 1(c)Panels I-III), while increases to
$Ra_{e}\sim 10$ as the ice thickness decreases in the melting case (Fig.
1(c)Panels IV-VI). Throughout the process of melting and freezing, the thermal
driving remains below the critical Rayleigh number ($Ra_{cr}$ is $\sim 2585$
in a laterally closed unit aspect ratio cell [29], $\sim 1708$ in the
laterally infinite RB system [30], and $\sim 1493$ in a laterally infinite RB
with linear buoyancy and melting [12]), and thus conduction holds [31].
Similarly, when the thermal driving is strong , i.e., $T_{b}=10^{\circ}$C,
(Regime-4, see Fig. 1(e)), the estimated effective Rayleigh number is
$Ra_{e}\sim 10^{7}\gg Ra_{cr}$ at the final equilibrium state, the turbulent
developed convection state prevails both for the freezing and melting
numerical experiments. However, freezing and melting behave differently in
Fig. 1(d) (Regime-3). For the freezing case, as the ice grows, the effective
aspect ratio (i.e., $l_{x}/(h-h_{i})$) of the fluid region increases, and thus
the convective rolls self-organize into several smaller rolls and the $Ra_{e}$
decreases accordingly, but the $Ra_{e}$ is still $\sim 4\times
10^{3}>Ra_{cr}$, leading to retaining convection. While for the melting case,
the system reaches an equilibrium when it is still in the conductive state.

Additionally, we explore the VC configuration in the same conditions.
Differently from the RB case, the global ice thickness $H_{i}$ always reaches
the same final equilibrium independently of the initial conditions or
independently of the history of the system (see SI for the corresponding
visualizations). In VC, the temperature gradient is perpendicular to the
gravity, making the system thermally unstable for any temperature gap or $Ra$.
This leads to the immediate breakdown of the stable stratification that is
observed in the RB.

The bistability observed in the phase-change RB can be quantitatively
predicted by means of a one-dimensional heat flux model. This is as follows:
We consider the system as the juxtaposition in series of compartments each one
characterized by specific thermal properties (we neglect the curvature of the
ice-water interface and the interface between the stably- and unstably-
stratified layers). The system reaches an equilibrium state when there is heat
flux balance among the conductive heat flux through the ice layer (denoted
$q_{i}$), the conductive heat flux through the stably stratified liquid layer
($q_{s}$), and the either conductive or convective heat flux (depending on
$Ra_{e}$) of the unstably stratified layer ($q_{u}$). Note that $q_{u}$ exists
only when $T_{b}>T_{c}$. The corresponding thicknesses for each layer are
denoted as $h_{i}$, $h_{s}$, and $h_{u}$ (see Fig. 12(a)). Depending on
whether $T_{b}$ is above $T_{c}$, two situations need to be considered: 1)
$T_{b}\leq T_{c}$, the fluid layer is gravitationally stable, so the heat
transfers diffusively both in the ice and water, so the definitions of the
heat fluxes are $q_{i}=k_{i}\frac{T_{0}-T_{t}}{h_{i}}$ and
$q_{s}=k_{w}\frac{T_{b}-T_{0}}{h_{\text{s}}}$, where $k_{i}$ and $k_{w}$ are
the thermal conductivity of ice and water, respectively; 2) $T_{b}>T_{c}$, the
density gradient takes opposite signs in the layers of temperature ranging
from $T_{0}$ to $T_{c}$ (gravitationally stably-stratified) and that ranging
from $T_{c}$ to $T_{b}$ (gravitationally unstably-stratified). The heat flux
in the gravitationally unstably-stratified layer can be predicted based on the
relation of the effective Nusselt number, $Nu_{e}$ (the dimensionless heat
flux normalized by the diffusive heat flux based on the unstably-stratified
layer), and the $Ra_{e}$. So the definitions of the heat fluxes are
$q_{i}=k_{i}\frac{T_{0}-T_{t}}{h_{i}}$,
$q_{s}=k_{w}\frac{T_{c}-T_{0}}{h_{\text{s}}}$, and
$q_{u}=\text{Nu}_{e}k_{w}\frac{T_{\text{b}}-T_{\text{c}}}{h_{\text{u}}}$,
where $Nu_{e}$ as a function of the $Ra_{e}$ behaves as in the classical RB
convection (as we already studied in [2]). In summary, depending on the value
of $T_{b}$ the equilibrium condition reads:

Case-1: $T_{b}\leq T_{\text{c}}$. The system is in a conductive state and is
independent of the water layer thickness, the heat flux balance reads

$q_{i}=q_{s}.$ (1)

From Eq. (1) we obtain $h_{i}=\frac{-k_{i}T_{t}}{k_{w}T_{b}-k_{i}T_{t}}h$.

Case-2: $T_{b}>T_{\text{c}}$.

$q_{i}=q_{s}=q_{u}.$ (2)

From Eq. (1) and Eq. (2), $h_{i}$ can be calculated which is reported in Fig.
1(b). The theoretical prediction (dashed and thick lines) agrees well with the
numerical simulation results (symbols). As $T_{b}$ increases, the system
experiences three regimes: 1) the equilibrium state is conduction
independently of the freezing or melting initial system configuration (blue
shaded area in Fig. 1(b)); 2) the freezing case reaches the convective
equilibrium state (dashed line in Fig. 1(b)), while the melting case reaches
the conductive equilibrium state, and this phenomenon is denoted by the
bistability (yellow shaded area in Fig. 1(b)); 3) the convective equilibrium
is reached independently of freezing or melting (red shaded area in Fig.
1(b)). This bistable regime which has been recently observed in numerical
simulation of fluids with a linear buoyancy force [12], is here verified in
the case of freezing pure water.

The one-dimensional model also allows us to understand which equilibrium state
is reached depending on the initial conditions. This is illustrated in Fig.
12, where we plot the values of the heat fluxes across the layers, $q_{i}$,
$q_{s}$, $q_{u}$ as a function of the ice thickness $h_{i}$. We note that in a
freezing experiment $h_{i}$ is initially null and shall grow till a matching
of the fluxes is reached, and vice-versa, i.e., in melting experiment
$h_{i}=h$ and decreases over time till the flux matching occurs. While the
matching is unique in the Case-1 (Fig. 12(b)) and for strongest convection
(Fig. 12(d)), the freezing front (bottom blue arrow) encounters first an
equilibrium position with smaller $h_{i}$ than the melting case (bottom red
arrow) (Fig. 12(c)) in the Case-2 with moderate thermal forcing.

Results: aspect ratio dependence. – We now focus on the case freezing of an
initially liquid-filled cavity with large thermal driving (high $Ra$). We
perform both experiments and numerical simulations of the aspect ratio
dependence in the RB and VC arrangements. The overall maximal variations of
$H_{i}$ for different $\Gamma$ are modest ($<20\%$) for the RB case, while
they can be large (up to $200\%$) for the VC case. For the RB, it is observed
that as $\Gamma\geq 1$, the global ice content nearly saturates (Fig. 9(c)),
and the only effect induced by larger $\Gamma$ is that the number of the
corresponding large-circulation convective rolls increases, with each roll of
about unit width-to-height ratio (see the one, two, and four convective rolls
for $\Gamma=1,2,4$, respectively, in the visualizations of the temperature
field from simulations in Fig. 9(a) panels II-III, II-IV, and II-V). We also
remark that the ice is locally flat in between the neighboring convective
rolls, in correspondence to the position from which cold plumes detach. The
characteristics of the ice front morphology are due to the shield of the
stably-stratified layer (with temperature ranging from $T_{0}$ to $T_{c}$).
Firstly, acting as a buffer layer, the stably stratified layer is able to
alleviate the effect of the hot plume impingement, creating a locally cold
environment; secondly, at the edge of the two neighboring large scale rolls,
the penetration from the unstably stratified layer does not occur, and this
creates a locally quiescent region, contributing to the local flattening of
the ice (however evident only in the simulations). Another interesting feature
is that multiple forms of stable ice profiles can be observed corresponding to
different convective states under the same external conditions: different
convective roll directions give rise to different ice front shapes (Fig. 9(a)
panels II-IV vs III-I, and panels II-V vs III-II). This multistability (or
bifurcation) of course already exists in the standard RB case where the flow
can rotate either in one direction or in the opposite depending on initial
tiny perturbations. When coupled with phase change, as the ice front reaches
one of the shape forms, the flow is locked to fit in with this specific shape,
leading to a preferred direction of convective rolls. We have not observed
large-scale circulation reversals in the system. In the current experiments,
we have no control over the convective pattern, e.g., we do not have a system
to heat up the liquid locally. While in the simulations we can achieve the
multiple equilibrium states by properly preparing the initial condition of the
simulation. This interesting feature might be exploited in future experiments
for flow control applications. On the other hand, when $\Gamma<1$ the system
equilibrates at smaller $H_{i}$. In these geometrically confined cases
($\Gamma=0.25$ and 0.5), several convective rolls stack upon each other
penetrate the stably-stratified layer and finally affect the ice front (Fig.
9(a) panels II-I and II-II). Note that the $H_{i}$ from the simulations and
the experiments agree well with each other except for very small $\Gamma$,
which is presumably because of the stronger influence from the external
environment in the experiments due to the larger area of the sidewalls.

For the VC cases (Fig. 9(b)), independently of $\Gamma$, it is observed the
robust ice front morphology due to the ubiquitous presence of two counter-
rotating convective rolls, which is unique to the system with density anomaly
working fluids [3]. One roll originates from the cold upward convective
current along the ice front (visualized by the blue arrow in Fig. 9(b)); the
other results from the hot plumes detaching from the hot plate (visualized by
the red arrow in Fig. 9(b)). The competition between these two rolls creates
the ice font morphology under different aspect ratios. As $\Gamma$ increases,
the interplay between these two convective rolls intensifies, leading to
highly mixing of the colder and warmer fluids, as corroborated by the highly
uneven $T_{c}$ isotherm in Fig. 9(b). At $\Gamma=4$, along the wide extent of
the ice front, the penetration of the hot plumes impingement can even modify
the local shape of the ice front. The global ice thickness first increases a
bit with increasing $\Gamma$, and then at large $\Gamma$, $H_{i}$ has a much
weaker dependence of $\Gamma$ (Fig. 9(c)). Note that we do not perform
simulations for the aspect ratio ranging of $\Gamma<0.25$ or $\Gamma>4$, the
reasons are 1) for even smaller $\Gamma$, multiple modes of convective roll
configuration are observed in the classical RB as a consequence of the
elliptical instability [32], and 2) for much larger $\Gamma$, the simulations
are quite expensive to carry out. Future work is required to map out a more
comprehensive picture of the aspect ratio dependence, with possible extensions
to three-dimensional geometries.

As observed above, the VC cases present robust forms of the ice front
morphology due to the presence of two counter-rotating convective rolls. We
measure the ice front shape for different aspect ratio cases we superpose them
in a single figure (Fig. 10(a)). The ice front morphology is similar
particularly at the bottom of the cell. To understand this feature, we use the
theoretical model proposed in [3] based on the concept of developing boundary
layer (BL). Here we only briefly sketch the model assumptions and its
construction. Although phenomenologically, it puts forward the physical
mechanism that dominates the form of the ice front. It can be observed that
there is a recirculating region induced by the upward convective currents
adjacent to the ice (with temperature ranging from $T_{0}$ to $T_{c}$, blue
arrow in Fig. 9(b)). This recirculating region can be separated into the
thermal BL along the ice front (orange area in Fig. 10(b)) and the main flow
region (light red area in Fig. 10(b)). When the system reaches an equilibrium,
there exists a balance between the heat flux through the BL and that through
the ice. A curvilinear coordinate, $S(x)$, is introduced to measure the length
of the ice front starting from its boundary point at $x=0$, and
$S(x)=\int_{0}^{x}\sqrt{1+(\tfrac{d(h_{i}(\xi))}{d\xi})^{2}}\,d\xi$. The BL
thickness is $\delta_{T}(S)$, which is a function of $S$, measured normal to
the ice front (see Fig. 10(b)). To analyze the heat flux across the ice and
the thermal boundary layer, we make the following approximations. For the ice
layer, the temperature gradient normal to the ice front in the ice layer can
be approximately estimated by $h_{i}(x)/\frac{dS(x)}{x}$; for the thermal
boundary layer, which is still nested inside the viscous boundary layer
($Pr\approx 11$), the heat transfer normal to the ice front in the BL is
mostly conduction [33]. The heat flux balance in the direction normal to the
ice surface can be expressed as,
$k_{i}\frac{(T_{0}-T_{t})}{h_{i}(x)\frac{dS(x)}{dx}}=k_{w}\frac{(T_{m}-T_{0})}{\delta_{T}(S(x))}$.
$\delta_{T}(S)$ is assumed to be $\delta_{T}(S)=C_{1}(S+C_{2})^{1/4}$ in which
there are two parameters
$C_{1}=c\cdot(g[1-\rho(T_{m})/\rho_{c}]/(\nu\kappa))^{1/4}$ with $c$ the
proportional constant ($c\approx 5$), and the offset $C_{2}$ because of non-
zero boundary layer thickness at $x=0$, with
$C_{2}=(h_{i}(0)[k_{w}(T_{m}-T_{0})]/[k_{i}(T_{0}-T_{t})]/C_{1})^{1/4}$ [34,
35, 36]. Fig. 10(c)-(g) show the comparison of the ice front morphology from
the experiments (blue shaded area), numerical simulations (green line), and
the theoretical model (red line). A good qualitative agreement is reached at
the bottom of the ice front, indicating that the exact form of the ice front
is indeed the thermal boundary layer dominated. The ice front morphology from
the numerical simulation and from the experiment is slightly different when
the aspect ratio is large (Fig. 10(f) and (g)), which is presumably in that,
near the large extent of the ice front, the mixing behavior between the colder
(with the temperature from $T_{0}$ to $T_{c}$) and the warmer (with the
temperature from $T_{c}$ to $T_{b}$) convective rolls intensifies, which is
attested by the highly irregular $T_{c}$ isotherm in Fig. 9(b)IV-II, and
I-III.

Conclusion and perspectives. – We have explored the equilibria of icing water
in a differentially heated cavity by means of laboratory experiments,
numerical simulations, and (phenomenologically) theoretical modeling. We
showed through simulations that the effect of the initial conditions, either
fully liquid or fully solidified, can lead to different equilibria, a
phenomenon only observed for the RB arrangement of the system in the moderate
thermal stratification case (Regime-3) and not in the VC case. We have also
introduced a simple one-dimensional model that in part relies on the known
heat-flux scaling relation in the RB system without melting, which is capable
of predicting the multiplicity of equilibria in the system and their
occurrence depending on the initial condition of the numerical experiment. In
the VC configuration, the system reaches the same equilibrium state
independently of the system’s initial condition, and this is due to the
intrinsic thermal instability of this configuration, which displays convection
for every Rayleigh number different from zero.
Secondly, we have studied experimentally and numerically the dependence of the
ice front morphology on the $\Gamma$ aspect-ratio in the freezing case for the
RB and VC. It was shown that for $\Gamma\geq 1$ in the RB case, the global ice
thickness is nearly independent of $\Gamma$. The ice front displays large
wavelength modulations and presents a local flatness in between neighboring
convective rolls. Multiple forms of the ice front morphology can be observed
under the same external conditions: different convective roll configurations
lead to different ice front shapes. Once the system reaches a specific ice
front shape, the convective roll appears locked in a certain direction. This
multistability behavior provides possibilities for flow control to maintain
desirable flow structures or solid-fluid interfacial morphology. For the VC
case, despite the major differences in the global ice extend as compared to
RB, a robust form of the ice front morphology can be observed at all $\Gamma$.
Such interfacial form is due to the presence of two counter-rotating
convective rolls and is dominated by a developing boundary layer in the upward
direction.

We hope that both the methodology of the present study on water
melting/freezing in idealized conditions and its findings will stimulate
further investigations. More complex natural and industrial phenomena are
still challenging because many crucial factors should be considered, e.g., the
effect of persistent shear flows [37, 38, 23] or induced by rotation [39], the
concentration gradients of solute components, and the extremely high-pressure
in deep water body environments, which are of great relevance for the
appropriate modeling of geophysical and climatological processes.

###### Acknowledgements.

This work was supported by the Natural Science Foundation of China under
Grants 11988102, 91852202, and 11861131005.

## References

  * [1] Brandeis G. Jaupart C. Earth Planet. Sci. Lett.771986345.
  * [2] Brandeis G. Marsh B. D. Nature3391989613.
  * [3] Sparks R., Annen C., Blundy J., Cashman K., Rust A. Jackson M. Philos. Trans. R. Soc. A377201920180019.
  * [4] Ulvrová M., Labrosse S., Coltice N., Råback P. Tackley P. Phys. Earth Planet. Inter.206201251.
  * [5] Meakin P. Jamtveit B. Proc. R. Soc. A4662010659.
  * [6] Deser C., Walsh J. E. Timlin M. S. J. Clim.132000617.
  * [7] Worster M. G. Annu. Rev. Fluid Mech.29199791.
  * [8] Glicksman M. E. Principles of solidification: an introduction to modern casting and crystal growth concepts (Springer Science & Business Media) 2010.
  * [9] Nazir H., Batool M., Osorio F. J. B., Isaza-Ruiz M., Xu X., Vignarooban K., Phelan P., Kannan A. M. et al. Int. J. Heat Mass Transf.1292019491.
  * [10] Skyllingstad E. D. Paulson C. A. J. Geophys. Res.1122007.
  * [11] Scagliarini A., Calzavarini E., Mansutti D. Toschi F. Modelling sea ice and melt ponds evolution: sensitivity to microscale heat transfer mechanisms in Mathematical Approach to Climate Change and its Impacts (Springer) 2020 pp. 179–198.
  * [12] Hester E. W., McConnochie C. D., Cenedese C., Couston L.-A. Vasil G. Phys. Rev. Fluids62021023802.
  * [13] Siggia E. D. Annu. Rev. Fluid Mech.261994137.
  * [14] Bodenschatz E., Pesch W. Ahlers G. Ann. Rev. Fluid Mech.322000709.
  * [15] Ahlers G., Grossmann S. Lohse D. Rev. Mod. Phys.812009503.
  * [16] Lohse D. Xia K.-Q. Annu. Rev. Fluid Mech.422010335.
  * [17] Chillà F. Schumacher J. Eur. Phys. J. E35201258.
  * [18] Kadanoff L. P. Phys. Today54200134.
  * [19] Rabbanipour Esfahani B., Hirata S. C., Berti S. Calzavarini E. Phys. Rev. Fluids32018053501.
  * [20] Favier B., Purseed J. Duchemin L. J. Fluid Mech.8582019437.
  * [21] Satbhai O., Roy S., Ghosh S., Chakraborty S. Lakkaraju R. Phys. Fluids312019096605.
  * [22] Purseed J., Favier B., Duchemin L. Hester E. W. Phys. Rev. Fluids52020023501.
  * [23] Toppaladoddi S. J. Fluid Mech.9192021A28.
  * [24] Vasil G. Proctor M. J. Fluid Mech.686201177.
  * [25] Wang Z., Jiang L., Du Y., Sun C. Calzavarini E. Phys. Rev. Fluids2020in press.
  * [26] Wang Z., Calzavarini E., Sun C. Toschi F. Proc. Natl. Acad. Sci.118(10)2021e2012870118.
  * [27] Calzavarini E. Software Impacts12019100002.
  * [28] Gebhart B. Mollendorf J. C. Deep Sea Research241977831.
  * [29] Luijkx J. Platten J. J. Non-Equilibrium Thermodynam.61981.
  * [30] Chandrasekhar S. Hydrodynamic and hydromagnetic stability (Courier Corporation) 2013.
  * [31] Satbhai O. Roy S. Phys. Fluids322020064107.
  * [32] Zwirner L., Tilgner A. Shishkina O. Physical Review Letters1252020054502.
  * [33] Sun C., Cheung Y.-H. Xia K.-Q. J. Fluid Mech.605200879.
  * [34] Bejan A. Convection heat transfer (John wiley & sons) 2013.
  * [35] White F. M. Corfield I. Viscous fluid flow Vol. 3 (McGraw-Hill New York) 2006.
  * [36] Shishkina O. Phys. Rev. E932016051102.
  * [37] Toppaladoddi S. Wettlaufer J. S. J. Fluid Mech.8682019648–665.
  * [38] Couston L.-A., Hester E., Favier B., Taylor J. R., Holland P. R. Jenkins A. J. Fluid Mech.9112021A44.
  * [39] Ravichandran S. Wettlaufer J. J. Fluid Mech.9162021A28.

Supplementary Information for:
Equilibrium states of the ice-water front in a differentially heated
rectangular cell

Experimental Methods. – The experiments for phase-change of pure water are
performed in a standard Rayleigh-Bénard (RB) convection cell (see Fig.
5(a)-(c)). Fig. 5(a) shows the sketch of the experimental cell, which consists
of a top cooling plate (with temperature $T_{t}$), a heating bottom plate
(with temperature $T_{b}$), and plexiglass adiabatic sidewalls (four values of
aspect ratio, $\Gamma=l_{x}/h$, are used, i.e., 0.5, 1, 2, 4 with the length,
$l_{x}=24$ cm, and width, $l_{z}=6$ cm, which are both fixed, and the height
$h=12$ cm, $24$ cm, $48$ cm and $96$ cm, respectively, here the lowercase
letter denotes physical variables with unit, while the capital letter to be
mentioned later is nondimensionalized variables normalized by $h$). Gaskets
are embedded in between the top plate and the sidewalls, and in between the
bottom plate and sidewalls, to create an optimal seal. The top and bottom
plates are well controlled at constant temperatures ($T_{t}<T_{0}$, and
$T_{b}>T_{0}$, with $T_{0}$ being the phase change temperature, which is
$T_{0}\sim 0^{\circ}$C at atmospheric pressure) by the circulating bath
(PolyScience PP15R-40), with the typical temperature fluctuations less than
$\pm 0.2$K. The temperature of the top and bottom plates are monitored by the
thermistors (44000 series thermistor element, see the sketch in Fig. 5(c)).
The working fluid, which is confined between the top and bottom plates, is
deionized and ultrapure water, of which the physical parameters (i.e.,
density, $\rho$, and the thermal expansion, $\alpha$) around the density peak
temperature, $T_{c}$ ($\sim$ 4∘C at atmospheric pressure) are shown in Fig.
5(d). The working fluid is carefully degassed by boiling twice before any
experiments are performed. During the phase-change process, to balance the
pressure induced by the volume change, an expansion vessel is connected to the
experimental cell through a rubber tube, and the expansion vessel is directly
open to the atmosphere so that the pressure of the experimental cell remains
unchanged. To reduce the heat exchange between the experimental cell and the
environment, we adopt two approaches: 1) the cell is wrapped in a sandwich
structure: insulation foam, aluminum plate (aiming to support and is the
temperature measuring spot of the sensor, PT 100, for the Proportional-
Integral-Derivative (PID) system), and insulation foam; 2) a PID controller is
utilized to control the temperature of the surrounding environment outside the
cell (see Fig. 5(b)). Given the difficulty of the manufacture of full-cell-
sized, bubble-free, and transparent ice, we did not perform the melting
experiment. Instead, we mainly perform the solidification experiments of the
aspect ratio dependence and investigate the bistability based on the numerical
simulation and the theoretical modeling.

Numerical Methods. – The simulations are conducted by the CH4-PROJECT code
[1], which adopts a Lattice-Boltzmann algorithm for the description of fluid
and temperature dynamics, and an enthalpy method for the phase change process.
The method has been validated against experiments [2, 3] and the code has been
intensively tested [4]. Here, only a brief description of the method is
presented and a more detailed description can be found in previous works [2,
3]. The simulation is based on the Boussinesq approximation except that for
the buoyancy term we use the non-monotonic relationship of the water density
as a function of the temperature, $\rho=\rho_{c}(1-\alpha^{*}|T_{b}-4|^{q})$,
with $\rho_{c}=999.972~{}kg/m^{3}$ the maximum density corresponding to
$T_{c}\approx 4^{\circ}$C, $\alpha^{*}=9.30\times 10^{-6}(K^{-q})$, and
$q=1.895$ [5]. The physical parameters for all simulations described in this
paper are $\nu_{w}=1.36\times 10^{-7}{m^{2}/s},k_{w}=0.57{W/(m\cdot
K)},k_{i}=2.37{W/(m\cdot K)},C_{pw}=4.21\cdot 10^{3}{J/(kg\cdot
K)},C_{pi}=1.96\cdot 10^{3}{J/(kg\cdot K)}$. We neglect the microscopic
physics leading to kinetic undercooling, Gibbs-Thomson effect and the
anisotropic growth/melting [6] in the simulations. We implement two source
terms to the energy conservation equation: 1) note that the key to accurately
solve such problems is to recover the diffusion term in the energy equation
exactly and, similarly to Ref. [7], the correction source term is implemented
which is induced by the heterogeneous media in the investigated domain; and 2)
the source term resulting from the latent heat contribution at the ice-water
interface. The governing equations read

$\begin{split}&\vec{\nabla}\cdot\vec{u}=0,\\\ &\frac{\partial\vec{u}}{\partial
t}+\vec{u}\cdot\vec{\nabla}\vec{u}=-\frac{\vec{\nabla}p}{\rho_{\text{c}}}+\nu_{w}\nabla^{2}\vec{u}+\alpha^{*}g|T-4|^{q}\hat{e}_{z},\\\
&\sigma(\rho\text{C}_{\text{p}})_{0}\frac{\partial T}{\partial
t}+\vec{\nabla}\cdot(\sigma(\rho\text{C}_{\text{p}})_{0}T\vec{u})=\vec{\nabla}\cdot(k\vec{\nabla}T)+S_{1}+S_{2},\end{split}$
(3)

where $\vec{u}(x,y,t)$, $p(x,y,t)$ and $T(x,y,t)$ are fluid velocity,
pressure, and temperature fields (all temperatures are measured in Celsius),
respectively and we denote with $x$ as the horizontal direction and with $y$
the vertical direction; $\nu_{w}$, $\rho$, $g$ and $k$ are the kinematic
viscosity of water, the density, the acceleration of gravity, and the thermal
conductivity (when it is in water phase $k=k_{w}$,
$\text{C}_{\text{p}}=\text{C}_{\text{pw}}$, and when it is in ice phase
$k=k_{i}$, $\text{C}_{\text{p}}=\text{C}_{\text{pi}}$) respectively.
$\hat{e}_{z}$ is the unit vector pointing in the direction opposite to that of
gravity. The first source term is
$S_{1}=-L\rho\frac{\partial\phi_{w}}{\partial t}$ and the second source term
$S_{2}=-\sigma
k\vec{\nabla}T\vec{\nabla}\frac{1}{\sigma}-\frac{\rho\text{C}_{\text{p}}}{\sigma}T\vec{u}\vec{\nabla}\sigma$.
Here $\sigma=\frac{\rho\text{C}_{\text{p}}}{(\rho\text{C}_{\text{p}})_{0}}$ is
the ratio of heat capacitance (which is variable and depends on the type of
phase, i.e., ice or water) and $(\rho\text{C}_{\text{p}})_{0}$ is reference
heat capacitance which is taken as constant [7].

The boundary conditions corresponding to the governing equations above are
isothermal at the top and bottom plates, no-slip at the bottom plate and at
the ice-water interface, adiabatic at the lateral boundaries, and no-slip and
freezing (namely, Stefan condition [8, 9]) at the phase-changing interface
(see also Figure. 5(e)). The boundary conditions read:

$\begin{split}&T(x,0,t)=T_{\text{b}},\\\ &T(x,1,t)=T_{\text{t}},\\\
&\vec{u}(x,0,t)=0,\\\ &\vec{u}(x,1-h_{\text{i}}(x,t),t)=0,\\\
&\frac{\partial{T(x,y,t)}}{\partial x}|_{x=0}=0,\\\
&\frac{\partial{T(x,y,t)}}{\partial x}|_{x=L_{x}/H}=0,\\\
&L\rho_{\text{i}}V_{n}=\vec{n}\cdot\vec{q}_{\text{w}}-\vec{n}\cdot\vec{q}_{\text{i}},\end{split}$
(4)

where $L$ is the latent heat, $V_{n}$ is the normal speed of the ice-water
interface, and $h_{i}(x,t)$ the dimensionless ice layer thickness at the
position $x$, $\vec{q}$ is the heat flux vector, $\vec{n}$ is a unit normal at
the ice-water interface pointing into the liquid. The subscripts I and w refer
to the ice and the water, respectively. The heat flux reads
$q_{\text{i}}=-k_{\text{i}}\vec{\nabla}T_{\text{i}}$ and
$q_{\text{w}}=-k_{\text{w}}\vec{\nabla}T_{\text{w}}$.

We monitor the global ice thickness and the ice front morphology, which are
expressed as $H_{i}(t)$ and $H_{i}(x,t)$, respectively ($H_{i}=1$ means full
solidification). We consider that the equilibrium is reached when the standard
deviation of $H_{i}(t)$ over a time-window of about 8 minutes is less than
0.5%.

Throughout the investigation, we limit the cold plate to a fixed temperature
of $T_{t}=-10^{\circ}$C. The effects of changing $T_{t}$ are qualitatively
predictable and are not expected to change the overall trends and mechanisms
discussed in the current study. The basic effects that $T_{t}$ may bring are
that, 1) the ice thickness at the equilibrium state of the system is thicker
(thinner) with decreasing (increasing) $T_{t}$; and 2) the corresponding ice
front morphology has a similar shape but with different extension and local
curvatures.

Comparison of one-dimensional model results based on different critical
Rayleigh number –

In the section of Results: Bistability of the main paper, we have mentioned
that the model is based on the critical Rayleigh number, $Ra_{cr}\sim 1708$,
which is estimated by the linear instability analysis and has been intensively
validated in [10, 11, 9]. When the system has phase change, the $Ra_{cr}$ may
change in the range of 1493 to 1708 [12], but nevertheless it does not change
the overall trend of the results in the current study.

Here, we show the comparison of one-dimensional (1D) model results based on
different choice of the critical Rayleigh number of $Ra_{cr}=1493$ and
$Ra_{cr}=1708$ respectively (see Fig. 6). It can be observed that the
equilibrium state is not sensitive to the specific value of the onset of
convection: with different $Ra_{cr}$ the model prediction only represents a
tiny shift of the temperature range of the bistability regime. Based on this,
we choose to use $Ra_{cr}=1708$ in the calculation of the 1D model throughout
the main paper.

The effect of changing the cold plate temperature $T_{t}$ – It has been
mentioned in the main paper that throughout the investigation, the cold plate
is fixed at temperature $T_{t}=-10^{\circ}$C. The effects of changing $T_{t}$
are qualitatively predictable and are not expected to change the overall
trends and mechanisms discussed in the current study. The basic effects that
$T_{t}$ may bring are that, 1) the ice thickness at the equilibrium state of
the system is thicker (thinner) with decreasing (increasing) $T_{t}$; and 2)
the corresponding ice front morphology has a similar shape but with different
extension and local curvatures. About the bistability, different $T_{t}$ only
changes the temperature range when the bistability occurs but the overall
trend remains similar.

Here, by means of theoretical modeling, we explain better the effects brought
by changing $T_{t}$.

Fig. 7 shows the theoretical prediction of the $H_{i}$ as a function of
$T_{b}$ for three different cold plate temperatures, $T_{t}=-5^{\circ}$C (blue
curves), $T_{t}=-10^{\circ}$C (red curves), and $T_{t}=-20^{\circ}$C (green
curves). The blue, red, and green shaded areas correspond to the bisability
regime for the three values of $T_{t}$ respectively.

It can be observed that a smaller amplitude of $T_{t}$ leads to thinner global
ice thickness, and the bistability regime occurs earlier as $T_{b}$ increases.
The temperature range of the bistability regime (shaded areas in Fig. 7)
becomes wider as $T_{t}$ increases. When the $T_{b}$ is normalized by the
temperature difference, $(T_{b}-T_{t})$, the bistability regime has a similar
range (see the inset of Fig. 7).

But nevertheless, the overall trend remains similar on different $T_{t}$
conditions. So we do not expect the mechanisms discussed in the current study
to change when varying $T_{t}$.

Melting and freezing in the vertical convection. – Here, we show detailed
results of the melting and freezing cases in the vertical convection.

It has been mentioned in the main article that four typical cases are selected
from the four distinct regimes as thermal driving increases (each regime has
different levels of coupling among ice front, stably- and unstably-stratified
layers [2]). To compare the Rayleigh-Bénard convection (RB) with the vertical
convection (VC), we conduct six sets of simulations for the VC case (we focus
on the cases where $T_{b}>T_{c}$), with all the other external conditions
remaining the same as those in the RB case, i.e., $4.5^{\circ}$C,
$5.5^{\circ}$C, and $10^{\circ}$C, and the cold plate temperature is fixed at
$T_{t}=-10^{\circ}$C. We perform long-term numerical simulations for melting
(initial condition is $H_{i}=1$) and freezing (initial condition is
$H_{i}=0$).

Fig. 8 shows the global ice thickness, $H_{i}$, as a function of time under
different conditions of thermal driving in the VC case. The thick (dashed)
lines represent the freezing (melting) case (denoted by “F” (“M”) in the
legend). It is observed that the $H_{i}$ reaches the same final equilibrium
state independent of the initial conditions or the history of the system.

We also show the evolution of the temperature field for the four thermal
driving conditions (see Figs. 9-11 for $4.5^{\circ}$C, $5.5^{\circ}$C, and
$10^{\circ}$C, respectively). In the VC configuration, there are always two
main convection rolls: 1) the colder one that is near the ice front and
originates from the cold plumes detaching from the ice-water interface; 2) the
warmer one originating from the hot plumes detaching from the hot plate.
Different hot plate temperatures are able to adjust the relative intensity of
the two convection rolls, which gives rise to different forms of ice front
shapes. If the colder roll wins (Fig. 9-10), the ice front tends to be
protected by this colder environment from deforming, and thus the ice front is
flat but only tilted. On the other hand, if the warmer roll wins (Fig. 11),
the ice front loses the shield of the cold roll and tends to display a local
maximum of the ice thickness in the middle, the thinnest at the top due to the
penetration of the hot plumes, and an increasing ice thickness from the
bottom.

Model of mechanism of bistability: – The flow chart in Fig. 12 summarizes the
iterative process used to compute the system equilibria. Firstly, $h_{i}$
varies in the range of $(0,~{}h)$, and the corresponding $q_{i}$ can be
calculated and plot as a function of $h_{i}$. Then an initial value of $h_{s}$
is set to start the iteration process. Once $h_{s}$ is set, $q_{s}$ and
$q_{u}$ are obtained which are also a function of $h_{i}$. Plotting $q_{s}$
and $q_{u}$ on top of the $q_{i}$ vs $h_{i}$, if one intersection is observed,
then the chosen value for $h_{s}$ is the solution for the equilibrium state,
otherwise we update $h_{s}$ and begin a new round of iteration until an
intersection is observed. By means of this method, it is observed two
intersections when the system is in the bistability regime, while only one
intersection in the convective equilibrium regime.

Figure 5: System configuration. (a) Sketch of the experimental system for
Rayleigh-Bénard convection coupled with phase-change of pure water. (b)
Temperature control facilities: the PID (Proportional-Integral-Derivative)
controller. (c) Sketch of resistance thermistor: 44000 series thermistor
element, which is used to measure the temperature of the top and bottom plates
. (d) The nonmonotonic relationship of the water density (black line) and
thermal expansion coefficient (red line) with temperature for cold water near
the density peak temperature, $T_{c}$ ($\sim 4^{\circ}$C at the atmospheric
pressure), from Ref. [5]. (e) Schematic description of the problem considered
in the numerical simulation (which is the same as the problem considered in
the experiment). The blue shaded region corresponds to the solid phase (ice)
and the red shaded region to the liquid phase (water). The top plate is cooled
at constant temperature $T_{t}$, and the bottom plate is heated at the
constant temperature $T_{b}$. The dashed line indicates the ice-water
interface with temperature $T_{0}$ and is of no slip condition. The two side
walls are of adiabatic conditions. The formulation for the boundary conditions
read,
$T(x,0,t)=T_{\text{b}},~{}T(x,1,t)=T_{\text{t}},~{}\vec{u}(x,0,t)=0,~{}\vec{u}(x,1-H_{\text{i}}(x,t),t)=0,~{}\frac{\partial{T(x,y,t)}}{\partial
x}|_{x=0}=0,~{}\frac{\partial{T(x,y,t)}}{\partial
x}|_{x=L_{x}/H}=0,~{}L\rho_{\text{i}}V_{n}=\vec{n}\cdot\vec{q}_{\text{w}}-\vec{n}\cdot\vec{q}_{\text{i}}$,
where $L$ is the latent heat, $V_{n}$ is the normal speed of the ice-water
interface, and $h_{i}(x,t)$ the dimensionless ice layer thickness at the
position $x$, $\vec{q}$ is the heat flux vector, $\vec{n}$ is a unit normal at
the interface pointing into the liquid. The subscripts I and w refer to the
ice and the water, respectively. The heat flux reads
$q_{\text{i}}=-k_{\text{i}}\vec{\nabla}T_{\text{i}}$ and
$q_{\text{w}}=-k_{\text{w}}\vec{\nabla}T_{\text{w}}$.  Figure 6: Theoretical
prediction of the bistability phenomena for RB case: $H_{i}$ as a function of
$T_{b}$ The red and blue curves correspond to the model prediction using
$Ra_{cr}=1493$ and $Ra_{cr}=1708$ respectively, for the condition of
$T_{t}=-10^{\circ}$. The symbols are results from the simulation.  Figure 7:
Theoretical prediction of the bistability phenomena for RB case: $H_{i}$ as a
function of $T_{b}$ for three different cold plate temperatures,
$T_{t}=-5^{\circ}$C (blue curves), $T_{t}=-10^{\circ}$C (red curves), and
$T_{t}=-20^{\circ}$C (green curves). The blue, red, and green shaded areas
correspond to the bisability regime for the three values of $T_{t}$
respectively. The inset is $H_{i}$ as a function of the normalized temperature
$T_{b}/(T_{b}-T_{t})$.  Figure 8: Evolution of the ice thickness, $H_{i}$, of
the freezing process (thick line) and the melting process (dashed line) for
the three thermal driving conditions (i.e., $4.5^{\circ}$C, $5.5^{\circ}$C,
and $10^{\circ}$C, respectively)  Figure 9: The evolution of the instantaneous
temperature field of the freezing process (top row) and melting process
(bottom row), with the bottom plate temperature $T_{b}=$4.5∘C. The black and
blue thick lines denote the ice-water interface with temperature $0^{\circ}$C
and the density peak temperature (i.e., $4^{\circ}$C) isotherm, respectively.
The arrows denotes the velocity vectors.  Figure 10: The evolution of the
instantaneous temperature field of the freezing process (top row) and melting
process (bottom row), with the bottom plate temperature $T_{b}=$5.5∘C. The
black and blue thick lines denote the ice-water interface with temperature
$0^{\circ}$C and the density peak temperature (i.e., $4^{\circ}$C) isotherm,
respectively. The arrows denotes the velocity vectors.  Figure 11: The
evolution of the instantaneous temperature field of the freezing process (top
row) and melting process (bottom row), with the bottom plate temperature
$T_{b}=$10∘C. The black and blue thick lines denote the ice-water interface
with temperature $0^{\circ}$C and the density peak temperature (i.e.,
$4^{\circ}$C) isotherm, respectively. The arrows denotes the velocity vectors.
Figure 12: The flow chat that summarizes the overall iterative process to
identify the solution of the thickness of the ice ($h_{i}$, with the
temperature ranging from $T_{t}$ to $T_{0}$), of the stably stratified layer
($h_{s}$, with the temperature ranging from $T_{0}$ to $T_{c}$), and the
unsteably stratified layer ($h_{u}$, with the temperature ranging from $T_{c}$
to $T_{b}$), and the corresponding heat flux through each layer is expressed
as $q_{i}$, $q_{s}$, and $q_{u}$, respectively. This flow chart applies to the
case of $T_{b}>T_{c}$ where three unknows exist (i.e., $h_{i}$, $h_{s}$, and
$h_{u}$). When $T_{b}\leq T_{c}$, $h_{c}\equiv 0$, so there are only two
unknowns that the heat flux $q_{s}$ and $q_{i}$ can be directly calculated
without using iterations. Inset: sketch of the one dimensional model: there
are three layers connecting in series and we assume the interface of the
neighboring layers is flat. (b) The model prediction of the conductive
equilibrium that has only one intersection.

## References

  * [1] Calzavarini E. Software Impacts12019100002.
  * [2] Wang Z., Calzavarini E., Sun C. Toschi F. Proceedings of the National Academy of Sciences118(10)2021e2012870118.
  * [3] Wang Z., Jiang L., Du Y., Sun C. Calzavarini E. Phys. Rev. Fluids2020in press.
  * [4] Rabbanipour Esfahani B., Hirata S. C., Berti S. Calzavarini E. Physical Review Fluids32018053501.
  * [5] Gebhart B. Mollendorf J. C. Deep Sea Research241977831.
  * [6] Dash J., Rempel A. Wettlaufer J. Reviews of modern physics782006695.
  * [7] Chen S., Yan Y. Y. Gong W. International Journal Heat and Mass Transfer1072017862.
  * [8] Alexiades V. Solomon A. D. Journal of Solar Energy Engineering1151993121.
  * [9] Bodenschatz E., Pesch W. Ahlers G. Annual review of fluid mechanics322000709.
  * [10] Pellew A. Southwell R. V. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences1761940312.
  * [11] Dominguez-Lerma M., Ahlers G. Cannell D. S. Physics of fluids271984856.
  * [12] Purseed J., Favier B., Duchemin L. Hester E. W. Physical Review Fluids52020023501.


# Cost-Based Budget Active Learning for Deep Learning

Patrick K. Gikunda
Computer Science Dept.
Paris8 University, France
<EMAIL_ADDRESS>Nicolas Jouandeau
Computer Science Dept.
Paris8 University, France
<EMAIL_ADDRESS>
###### Abstract

Majorly classical _Active Learning_ (AL) approach usually uses statistical
theory such as entropy and margin to measure instance utility, however it
fails to capture the data distribution information contained in the unlabeled
data. This can eventually cause the classifier to select outlier instances to
label. Meanwhile, the loss associated with mislabeling an instance in a
typical classification task is much higher than the loss associated with the
opposite error. To address these challenges, we propose a _Cost-Based Bugdet
Active Learning_ (CBAL) which considers the classification uncertainty as well
as instance diversity in a population constrained by a _budget_. A principled
approach based on the min-max is considered to minimize both the labeling and
decision cost of the selected instances, this ensures a near-optimal results
with significantly less computational effort. Extensive experimental results
show that the proposed approach outperforms several state-of -the-art active
learning approaches.

##  1 Introduction

_Active Learning_ (AL) considers that if a Learning Algorithm can choose an
instance to label for training a model, then the instances choosen should
maximize the learning performance under a fixed budget [16]. Typically this
process involves randomly sampling large amount of data from underlying
distribution for training a model. In _Deep Networks_ (DN) acquiring labels
for training set is very costly and time consuming even when using state-of-
the-art computing resources. For example, its expensive to hire dermatologists
to annotate 129,450 skin cancer images [5]. However, _Machine Learning_ (ML)
algorithms does not require all the training data to be labeled [12].

In AL there are three scenarios in which the ML algorithm will query the label
of an instance, they include: a) _Membership Query Synthesis_ that generates
constructs of an instance from underlying distribution [2]. The queries of an
instance are generated then labels are requested. In this scenario the quality
of the randomly generated instances is not guaranteed; b) _Stream-Based
Selective Sampling_ that uses query strategy to determine whether to query the
label of an instance or reject it based on some utility metrics [20]. The
model sequencially picks data and checks the data one by one to determine
whether to label the picked sample or not. In this scenario, selecting
unlabeled instances comes at no or minimal cost; c) _Pool-Based Sampling_
large pool of unlabled instances gathered at once and then examples are ranked
in order of informativenes. The labels of most informative instances are
queried [13]. The remainder of the paper is structured as follows: Section 2
discusses recent relevant literature on AL approaches. Section 3 describes the
instance selection method used. Section 4 describes the datasets used, other
AL methods and experimental results. Finally, Section 5 concludes the work by
presenting some insights for further research.

##  2 Related Work

AL methods in deep networks use different strategies or a combination of the
strategies to query for labels [7]. The strategies range from density
estimation to multi-factor methods. The strategies can be categorized into two
groups: _population-based_ strategies and _pool-based_ strategies [18]. In
population-based AL, training and test sets are drawn from the same
distribution with an assumption that training and test data both follow the
same conditional distribution $p(y|x)$. In this type, the objective is to find
the optimal training input density to generate the training input instances.
In pool-based AL, the objective is to optimally select some unlabeled
instances from a pool so that a model trained from them can best label the
remaining samples. Regardless of whether it is population-based or pool-based,
AL is an iterative process [14]. It first builds a base model from a small
number of labeled training instances, and then using different or a
combination of utility metrics it selects unlabeled instances and queries for
their labels. The newly labeled instances are added to the training labeled
set and the model is updated. This process iterates until a termination
criterion is met, for example, when the labeling budget is exhausted or the
maximum number of iterations is met. Based on the number of unlabeled
instances to query at each iteration, AL methods can be grouped as either
sequence-mode AL, where one instance is queried each time or batch-mode AL,
where multiple instances are queried at each iteration [14].

This paper focuses on pool-based batch-mode AL for DN. Although numerous AL
approaches have been proposed in the literature [16], a number of them are for
standard pool-based AL problems e.g. Donmez _et al._ presents a dynamic
approach that updated selection parameters based on estimated future residual
error reduction after each actively sampled instance [3], Settles and Craven
formulates the implementation of AL using information density [17] and Krempl
_et al._ implements cost-sensitive probabilistic approach for binary
classification [11]. Among those approaches limited to AL for DN include the
work presented by Ducoffe and Precioso that uses margin theory to compute
instances along the decision boundary, [4], [1], [8] and [9].

The emphasis in AL is to evaluate the informativeness of an instance, with an
assumption that an instance with higher classification uncertainty is more
crucial to label. This classical approach usually uses statistical theory such
as entropy and margin to measure instance utility, however it fails to capture
the data distribution information contained in the unlabeled data. This can
eventually cause the classifier to select outlier instances to label.
Therefore, its important to consider the classification uncertainty as well as
instances diversity in a population while developing an AL solution. In our
approach we consider both the uncertainty and correlation measure to calculate
the most informative and representative instance, which we refer as a high
confidence instance.

In AL where there is a pre-determined budget on labeling, it is important to
estimate the objective function for data selection. This guarantees near-
optimal results with significantly less computational effort. The aim is to
maximize the objective function while minimizing data acquisition costs (or to
remaining within a budget). To deal with this problem we formulate the most
informative budget selection task as a continuous optimization problem. The
aim is to determine possible queries that maximize the improvement to the
classifier’s strategy, without overspending the budget. The proposed model
addresses the cost-sensitive learning problem based on learning algorithms
that construct models for class probability estimation $p(y|x)$. The
probability estimates provide an easy means for factoring in the
misclassification losses in the classification decision making step. To
address the labeling cost problem, we employ the same probability estimation
techniques over the selected high confidence instances. The two methods are
then combined into an algorithm that minimizes the combined cost. At each
iteration (while the total sum of annotation cost is under a given budget) a
high confidence instance 111High confidence instances are the most informative
and representative instances selected from the unlabeled set to label is
selected.

##  3 Selecting Most Informative Instance

We consider the problem of actively selecting a batch of instances to label,
where the contents of the batch must be constrained by some budget. We will
use the following notation in this paper. Let $x_{i}$ represents an instance
and $y_{i}$ where $y_{i}\in\\{-1,+1\\}$ represents the class label for
$x_{i}$, $D=D^{L}\cup D^{U}$, $D^{L}$ denotes labeled instances where
$D^{L}=\\{(x_{1},y_{1}),(x_{2},y_{2}),....,(x_{n},y_{n})\\}$, $D^{U}$ denotes
unlabeled instances where $D^{U}=\\{(x_{1},?),(x_{2},?),....,(x_{n},?)\\}$,
$D^{H}$ denotes high confidence instances and $\Theta$ denotes the model
defined by model parameters. For label space $Y$ with $K$ classes in $D$ we
use the class probability estimator $P(y|x)$ to compute the estimate of a
label. In order to avoid the problem of generalization of unseen instances and
to learn an accurate model, we present a robust approach that uses different
utility metrics and a cost function. The utility metrics considered in this
work are uncertainty, correlation and informativeness measure, thus we present
three main components of our approach: a) uncertainty measure, b) correlation
measure and the cost-based labeling.

###  3.1 Uncertainty Measure

Given a label space $Y$ the uncertainty measure $f_{u}$ of an instance
considering both the features and the label can be defined as:

$f_{u}(x):\begin{cases}L^{S}\rightarrow R,&\mbox{(i) }\mbox{features view}\\\
(D^{U}\times L^{S})\rightarrow R,&\mbox{(ii) }\mbox{features-label
view}\end{cases}$ (1)

to a real number space $R$. From Equation 1 above: (i) the uncertainty measure
is computed from instance-features only while (ii) the uncertainty measure is
computed from both the instance-features and instance-label. In our method we
consider the uncertainty measure computed from instance features and label
view which is considered the most effective [10]. Out of the three common
uncertainty measure criteria namely least confidence, sample margin and
entropy, sample margin is prefered since it integrates the second most
probable class label in the uncertainty metric hence able to reduce the error
rate by defining the decision boundary. We therefore define uncertainty
measure as:

$f_{u}(x)=P(y_{i}=l_{1}|x_{i};Y)-P(y_{i}=l_{2}|x_{i};Y)$ (2)

With $l_{1}$ and $l_{2}$ being the most likely and second most likely labels.
High uncertainty value $f_{u}$ implies current model have little knowledge of
the instance, and including it into the training set can help improve to the
prediction performance of the model.

###  3.2 Correlation Measure

When developing efficient AL methods, it is critical to consider samples
distribution information [19]. The instance diversity information aids in
selecting most representative instances. In order to have more information
about the unlabeled instances it is appropriate to select a candidate instance
in a more dense region. In addition, selecting an instance to label only based
on uncertainty measure may lead to selecting an outlier instance, therefore
exploiting sample instance diversity will provide the most informative
instance to label. Our method is based on the fact that the trade off between
instance uncertainty and correlation is an essential AL problem to address.
Given a label space $Y$, we can define different groups of correlation of an
instance $x$ in a set of unlabeled set as;

$f_{c}(x):\begin{cases}D^{U}\times D^{U}\rightarrow R,&\mbox{feature view}\\\
Y\times Y\rightarrow R,&\mbox{label view}\\\
(D^{U},y)\times(D^{U},y)\rightarrow R,&\mbox{combined view}\end{cases}$ (3)

to a real number space $R$. In Equation 3, the combination of feature and
label correlation is called combined view. Different algorithms exist for
exploiting this type of combination [10]. Majorly these algorithms are used in
a multi-label learning tasks when an instance has more than one label. This
setting is ideal for mining tasks on instances with complex structure. In this
work we focus on exploiting the pairwise similarities of instances, therefore
the informativeness of an instances is weighed by average similarity to its
neighbours. Let $x_{i}$ and $x_{j}$ be a pair of instances. To cope with the
drawback of uncertainty based selection, we then consider the diversity by
evaluating the correlation of the instances. Given a label space $Y$ the
correlation measure $f_{c}(x_{i},x_{j})$ between a pair of instances $x_{i}$
and $x_{j}$ can be defined as:

$f_{c}(x)=\frac{1}{D^{U}}\sum_{x_{j}\in D^{U}}(x_{i},x_{j})$ (4)

The value of $f_{c}(x_{i})$ represents the instance density of $x_{i}$ in the
unlabeled set. The larger the value, the more densely an instance is
correlated with others. A low value of the correlation measure indicates an
outlier instance which should not be considered for labeling. Our motivation
is that the most representative instances of a distribution can be very
informative for improving the generalization performance. Therefore, given
correlation measure $f_{c}(x_{i})$ and uncertainty measure $f_{u}(x_{i})$ the
informativeness of an instance can be defined as:

$f_{i}(x)=f_{u}(x_{i})\cdot f_{c}(x_{i})$ (5)

It can be rewritten as:

$x_{i}=\underset{i}{\mathrm{argmax}}(u_{i}\cdot c_{i})$ (6)

###  3.3 Cost-Based labeling

In our approach the high confidence instance evaluation is based on the
instance informativeness which is computed from both uncertainty and
correlation metrics. The model is trained on labeled instances: feature and
label views. After querying for an high confidence unlabeled instance, a model
prediction result is generated based on output probability distribution. Each
instance $x_{i}=\\{f_{1}^{i},f_{2}^{i},...f_{q}^{i},y^{i}\\}$ in labeled set
$D^{L}=\\{x_{1},x_{2},....x_{s}\\}$ is represented in a feature space $F$
consisting of a feature space and its class label $y^{i}$. The size of $D^{L}$
is denoted by $s$ and $x_{i}$ denoted the $i$th instance in $D^{L}$. The
prediction can be denoted as a mapping function from the feature space $F$ to
the class label space $Y$ which can be expressed as;

$P(x):F\mapsto Y$ (7)

The query strategy used in this work is based on the value of $f_{i}$
discussed in equation 6. Instances are ranked based on the value $f_{i}$ with
top ranked instances being the most appropriate to label. Under the current
distribution $P(y_{i}|x_{i};Y)$ each possible instance $(x_{i},?)$ from the
selected instances $D^{H}$ will be labeled with label $y_{i}$. When $y_{i}=1$,
$x_{i}$ is regarded as a high confidence sample. The model update strategy is
to train a model based on the information provided by model weights computed
from model validation of the performance. We employ the probability estimators
in our approach to both minimize the labeling costs and the misclassification
decisions. We make the assumption that the loss function associated with the
decisions is represented as a static $K$ and a loss matrix $L$ available at
learning time. The contents of $L(i,j)$ specify the cost incurred when an
example is predicted to be in class $i$ when in fact it belongs to class $j$.
Therefore, high confidence instance selection criteria in this study will be
based on probability of $x_{i}$ belonging to $K^{th}$ class which can be
expressed as:

$y_{i}={argmin_{\xi\in
D^{H}}(C(\xi)}{\mathscr{x+}{\sum_{k=1}^{K}}P(k|\xi){\sum_{i}}{\sum_{j=1}^{K}}P_{\xi,k}(j|x_{i})L(y_{i},j))}$
(8)

The Algorithm 1 describes the _Cost-Based Budget Active Learning_ (CBAL) with
budget labeling.

1 Input: labeled instance set $D^{L}$, unlabeled instance set $D^{U}$, loss
matrix $L$, labeling cost $C$, empty set $D^{H}$, a budget $m$;

2 Output: model $\Theta$;

3 $\Theta\leftarrow$ getModel $(D^{L})$;

4 while  _$|D^{L}| <m$_ do

5 for  _each $x_{i}$ in $D^{U}$_ do

6 $u_{i}\leftarrow f_{u}(x_{i})$;

7 $c_{i}\leftarrow f_{c}(x_{i})$;

8 $x^{*}\leftarrow{\underset{i}{\mathrm{argmax}}(u\cdot c)}$;

9 $D^{H}\leftarrow D^{H}\cup\\{x\\}$;

10

11 for each $j$ learn $P(y|D^{H})$;

12 $y_{i}\leftarrow$ getLabel using Eq.8;

13 $D^{H}\leftarrow D^{H}\setminus\\{y_{i}\\}$;

14 $D^{L}\leftarrow D^{L}\cup\\{y_{i}\\}$;

15 $\Theta\leftarrow$ updateModel $(D^{L})$;

16

17return $\Theta$;

Algorithm 1 Cost-Based Budget Active Learning (CBAL).

In Algorithm 1, the labeling is defined by the budget $m$ with model updates
after each iteration (lines 4-14). At first the base model is trained using
the initial set of labeled data $D^{L}$. Instance evaluation is done to
identify the most informative and representative instance to label (lines
5-9). This evaluation returns the high confidence instances $D^{H}$ selected
from the unlabeled population (line 9). For each of the selected instance, its
label is queried and consequently the labeled set is updated. The model
selection strategy is updated with the learned parameters after every
iteration. CBAL is designed to train a classification model using a small
labeled population sample proportion.

##  4 Experiments

We conduct experiments with 12 real-world data sets (_wine
222http://archive.ics.uci.edu/ml/datasets/Wine,
seeds333https://archive.ics.uci.edu/ml/datasets/seeds, v2-plant
seedling444https://www.kaggle.com/vbookshelf/v2-plant-seedlings-dataset,
liver555https://archive.ics.uci.edu/ml/datasets/Liver+Disorders,
sonar666https://archive.ics.uci.edu/ml/datasets/Connectionist+Bench+(Sonar,+Mines+vs.+Rocks,
vehicle777https://archive.ics.uci.edu/ml/datasets/Statlog+Vehicle+Silhouettes,
breast888https://www.kaggle.com/paultimothymooney/breast-histopathology-
images, diabetic999https://www.kaggle.com/sovitrath/diabetic-
retinopathy-224x224-gaussian-filtered,
heart101010https://archive.ics.uci.edu/ml/datasets/Heart+Disease,
isolet111111https://archive.ics.uci.edu/ml/datasets/isolet,
plant121212https://www.kaggle.com/vipoooool/new-plant-diseases-dataset,
svhn131313https://www.kaggle.com/stanfordu/street-view-house-numbers_)
previously considered by other authors in similar domain. For each data set,
we split 80% of the instances as the training set, and the balance 20% as the
test set to evaluate the prediction accuracy of the models. We select a subset
of instances from the training data to query (100 instances per query) for
labels and then construct a base classification model according to these
labeled data. The goal is to pick out high confidence instances such that the
constructed model maintains effective classificaion ability. The training is
implemented in a batch-mode AL. We compare our method with other state-of-the-
art methods: a). Random Sampling (RS) which selects a certain number of
samples from a given set and quire labels; b). Core-Set AL (CSAL) [15] which
defines the AL problem as a competitive sample core-set selection which is
then applied to a CNN in a batch setting; c). Deep Bayesian Active Learning
(DBAL): a Bayesian framework proposed in [6] for high dimensional data which
considers Deep Learning problem of dependence on big amount of data; d).
Adversarial AL for deep networks (AAL) a margin based approach proposed in [4]
for deep networks with intention of reducing the number of queries to the
oracle during training. The both the budget $m$ and the initial labeled set is
specified before start of iterations. A batch size of 64 was considered for
all iterations for both training and testing selection. 100 queries were
considered for each iteration.

###  4.1 Results and Discussion

Figure 1 shows the classification accuracy of different active learning
approaches with varied number of queries. From the observation RS tends to
yield better performance when the number of queries is small but as the number
of queries increases it starts to slow its effectivenes in prediction. This
observation might be as a result of sampling bias induced by an intelligent
selection strategy. CSAL that define AL as a core-set problem, is not
performing well at the start of training. As the number of queries increases,
there is improvement and yields better performance. This is because with few
training instances, the learned decision boundary tends to be inaccurate, and
as a result, the unlabeled instances near the decision boundary may not be the
most high confidence instances to label. The performance of DBAL is better on
some datasets but performs poorly on others. This inconsistency might be as a
result of identified cluster structure of unlabeled data that is not always
consistent with the target classification model. The behavior of AAL is
similar to that of DBAL. Finally, we observe that for most cases, CBAL is able
to outperform the baseline methods significantly and we attribute this success
to the principle selecting high confident samples at each iteration, and the
specially aspect of minimizing the labeling and decision cost after each
subsequent iteration.

##  5 Conclusion and Future Work

We propose a new near optimal AL approach called CBAL, that measure both the
informative and representative of an instance using instance utility to get a
high confidence instance to lable while minimizing the labeling and decision
cost. The proposed approach of minimizing cost is based on the minmax
principled view. Our current work is based on budget constrained learning with
pairwise similarities of instances. In the future, we plan to extend this work
to multi-label learning tasks by considering instances with more than one
label. In addition we plan to consider the expert knowledge in the training by
allowing the user to control tradeoff between selection and labeling, this
will lead to incoporating domain knowledge into AL

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$30$$40$$50$$60$$70$$80$# of queries
(100)CBALRSCSALAALDBAL (a) wine

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$30$$40$$50$$60$$70$# of queries
(100)CBALRSCSALAALDBAL (b) seeds

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$30$$40$$50$$60$$70$$80$$90$# of queries
(100)CBALRSCSALAALDBAL (c) v2-plant seedling

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$# of queries
(100)Prediction AccuracyCBALRSCSALAALDBAL (d) liver

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$# of queries
(100)CBALRSCSALAALDBAL (e) sonar

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$# of queries
(100)CBALRSCSALAALDBAL (f) vehicle

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$30$$40$$50$$60$$70$$80$$90$# of queries
(100)CBALRSCSALAALDBAL (g) breast

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$$80$$90$# of
queries (100)CBALRSCSALAALDBAL (h) diabetic

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$# of queries
(100)CBALRSCSALAALDBAL (i) heart

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$20$$30$$40$$50$$60$$70$$80$$90$# of queries
(100)Prediction AccuracyCBALRSCSALAALDBAL (j) isolet

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10$$20$$30$$40$$50$$60$$70$$80$# of queries
(100)CBALRSCSALAALDBAL (k) plant

$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$# of
queries (100)Prediction AccuracyCBALRSCSALAALDBAL (l) svhn

Figure 1: Classfication Accuracy on different datasets

## References

  * [1] Hamed H Aghdam, Abel Gonzalez-Garcia, Joost van de Weijer, and Antonio M López, ‘Active learning for deep detection neural networks’, in Proceedings of the IEEE International Conference on Computer Vision, pp. 3672–3680, (2019).
  * [2] Dana Angluin, ‘Queries and concept learning’, Machine learning, 2(4), 319–342, (1988).
  * [3] Pinar Donmez, Jaime G Carbonell, and Paul N Bennett, ‘Dual strategy active learning’, in European Conference on Machine Learning, pp. 116–127. Springer, (2007).
  * [4] Melanie Ducoffe and Frederic Precioso, ‘Adversarial active learning for deep networks: a margin based approach’, arXiv preprint arXiv:1802.09841, (2018).
  * [5] Andre Esteva and Brett Kuprel  _et al_ , ‘Dermatologist-level classification of skin cancer with deep neural networks’, Nature, 542(7639), 115, (2017).
  * [6] Yarin Gal and Riashat Islam  _et al_ , ‘Deep bayesian active learning with image data’, in ICML, pp. 1183–1192, (2017).
  * [7] M Hasenjäger and H Ritter, ‘Active learning in neural networks’, in New learning paradigms in soft computing, 137–169, Springer, (2002).
  * [8] Manuel Haussmann, Fred A Hamprecht, and Melih Kandemir, ‘Deep active learning with adaptive acquisition’, arXiv preprint arXiv:1906.11471, (2019).
  * [9] Steven CH Hoi, Rong Jin, Jianke Zhu, and Michael R Lyu, ‘Batch mode active learning and its application to medical image classification’, in ICML, pp. 417–424, (2006).
  * [10] Sheng-Jun Huang, Nengneng Gao, and Songcan Chen, ‘Multi-instance multi-label active learning.’, in IJCAI, pp. 1886–1892, (2017).
  * [11] Georg Krempl, Daniel Kottke, and Vincent Lemaire, ‘Optimised probabilistic active learning (opal)’, Machine Learning, 100(2-3), 449–476, (2015).
  * [12] Mingsheng Long and Han Zhu  _et al_ , ‘Deep transfer learning with joint adaptation networks’, in ICML, pp. 2208–2217, (2017).
  * [13] Kamal Nigam and Andrew McCallum, ‘Pool-based active learning for text classification’, in Conference on Automated Learning and Discovery (CONALD), (1998).
  * [14] Hiranmayi Ranganathan, Hemanth Venkateswara, Shayok Chakraborty, and Sethuraman Panchanathan, ‘Deep active learning for image classification’, in 2017 IEEE International Conference on Image Processing (ICIP), pp. 3934–3938, (2017).
  * [15] Ozan Sener and Silvio Savarese, ‘Active learning for convolutional neural networks: A core-set approach’, arXiv:1708.00489, (2017).
  * [16] Burr Settles, ‘Active learning literature survey’, Technical report, University of Wisconsin-Madison, Dept of Computer Sciences, (2009).
  * [17] Burr Settles and Mark Craven, ‘An analysis of active learning strategies for sequence labeling tasks’, in Proceedings of the 2008 Conference on Empirical Methods in Natural Language Processing, pp. 1070–1079, (2008).
  * [18] Masashi Sugiyama and Shinichi Nakajima, ‘Pool-based active learning in approximate linear regression’, Machine Learning, 75(3), 249–274, (2009).
  * [19] Martin Szummer and Tommi S. Jaakkola, ‘Information regularization with partially labeled data’, in Advances in Neural Information processing systems, pp. 1049–1056, (2003).
  * [20] Xingquan Zhu and Peng Zhang  _et al_ , ‘Active learning from data streams’, in IEEE ICDM, pp. 757–762, (2007).


should be inseparable either and hence, so should be the Frobenius map

$\hat{\phi}:\left\\{\begin{array}[]{ccc}\mathscr{E}^{(p)}&\longrightarrow&\mathscr{E}^{(p^{2})}\\\
(x,y)&\longmapsto&(x^{p},y^{p}).\end{array}\right.$

Thus, $\mathscr{E}^{(p^{2})}=\mathscr{E}$ and hence $\mathscr{E}$ is defined
over $\mathbb{F}_{p^{2}}$. ∎

######  Theorem 81.

The number of $\overline{\mathbb{F}}_{p}$–isomorphism classes of supersingular
elliptic curves over $\overline{\mathbb{F}}_{p}$ equals

$\left\lfloor\frac{p}{12}\right\rfloor+\left\\{\begin{array}[]{ccccc}0&\text{if}&p&\equiv&1\mod
12\\\ 1&\text{if}&p&\equiv&5\mod 12\\\ 1&\text{if}&p&\equiv&7\mod 12\\\
2&\text{if}&p&\equiv&11\mod 12.\end{array}\right.$

###### Proof.

See [Sil09, Thm. V.4.1]. ∎

###  6.3 Proof of the main theorem

With these two last statements at hand we can finally provide the proof of
Theorem 41. We restrict the proof to the case $p\geqslant 5$. Note that
Tsfasman–Vlăduţ–Zink theorem remains true when $p=2,3$ but the coding
theoretic interest is rather limited.

###### Proof of Theorem 41.

Consider the sequence of curves $X_{0}(\ell)$ for $\ell\equiv 11\mod 12$. From
Theorem 74 it has genus $g_{\ell}=\frac{\ell+1}{12}$. From Theorem 81, the
curve $X_{0}(1)$ has at least $\frac{p-1}{12}$ $\mathbb{F}_{p^{2}}$–rational
points corresponding to isomorphism classes of supersingular elliptic curves.
Such an elliptic curve with no nontrivial automorphism has $\ell+1$ pre-images
in $X_{0}(\ell)$ which also correspond to supersingular elliptic curves and
hence are $\mathbb{F}_{p^{2}}$–rational points. Depending on the class of $p$
modulo $12$ the curves with $j$–invariant $0$ and $1728$ may be supersingular.
More precisely, from [Sil09, Ex. V.4.4 & V.4.5],

  * •

for $p\equiv 1\mod 12$, both curves are ordinary (i.e. non supersingular) and
then any supersingular curve has no nontrivial automorphism and hence has
$\ell+1$ pre-images in $X_{0}(\ell)(\mathbb{F}_{p^{2}})$. Therefore, from
Theorem 81,

$\sharp{X_{0}(\ell)(\mathbb{F}_{p^{2}})}\geqslant(\ell+1)\frac{p-1}{12}\cdot$

  * •

for $p\equiv 5\mod 12$, the curve with $j=0$ is supersingular and the one with
$j=1728$ is ordinary. Therefore, there are $\frac{p-5}{12}+1$ supersingular
curves and all of them but one have $\ell+1$ distinct pre-images. The
remaining curve is the one with $j=0$ and an automorphism group of order $6$.
Its treatment is very similar to the proof of Theorem 74. Consider the action
of the automorphism of order $6$ on the $\ell$–torsion. From Lemma 77, this
induces an automorphism $\eta$ or order $6$ of
$\mathscr{E}[\ell]\simeq\mathbb{F}_{\ell}^{2}$. Since $\ell\equiv 11\mod 12$,
then $\ell\equiv 2\mod 3$ and hence $\mathbb{F}_{\ell}$ does not contains the
sixth roots of unity. Thus, $\eta$ is not diagonalisable and hence cannot fix
a line. Consequently, one proves that the orbit of any line is the union of
$3$ distinct lines ($\eta^{3}=-\text{Id}$, which fixes the lines). Therefore,
the curve with $j$–invariant $0$ has $\frac{\ell+1}{3}$ pre-images and
Consequently, using Theorem 81:

$\sharp{X_{0}(\ell)(\mathbb{F}_{p^{2}})}\geqslant(\ell+1)\frac{p-5}{12}+\frac{\ell+1}{3}=(\ell+1)\frac{p-1}{12}\cdot$

  * •

for $p\equiv 7\mod 12$, the curve with $j=0$ is ordinary and the one with
$j=1728$ is supersingular. A similar reasoning yields

$\sharp{X_{0}(\ell)(\mathbb{F}_{p^{2}})}\geqslant(\ell+1)\frac{p-7}{12}+\frac{\ell+1}{2}=(\ell+1)\frac{p-1}{12}\cdot$

  * •

for $p\equiv 11\mod 12$, both curves are supersingular and we get:

$\sharp{X_{0}(\ell)(\mathbb{F}_{p^{2}})}\geqslant(\ell+1)\frac{p-11}{12}+\frac{\ell+1}{2}+\frac{\ell+1}{3}=(\ell+1)\frac{p-1}{12}\cdot$

In summary, we always have a lower bound $(\ell+1)\frac{p-1}{12}$ on the
number of rational points. Then

$\frac{\sharp{X_{0}(\ell)(\mathbb{F}_{p^{2}})}}{g_{\ell}}\geqslant\frac{\frac{p-1}{12}(\ell+1)}{\left(\frac{\ell+1}{12}\right)}=p-1.$

Thus, over $\mathbb{F}_{q}$ for $q=p^{2}$, we identified a family of curves
whose number of $\mathbb{F}_{q}$–rational points goes to infinity and whose
ratio, number of $\mathbb{F}_{q}$–points divided by the genus goes to
$\sqrt{q}-1$. Which turns out to be optimal. ∎

## References

  * [BH08] Peter Beelen and Tom Høholdt.  The decoding of algebraic geometry codes.  In Advances in algebraic geometry codes, volume 5 of Ser. Coding Theory Cryptol., pages 49–98. World Sci. Publ., Hackensack, NJ, 2008\.
  * [Cou16] Alain Couvreur.  Introduction to coding theory, 2016.  Personal lecture notes.
  * [CR21] Alain Couvreur and Hugues Randriambololona.  Algebraic geometry codes and some applications.  In W. Cary Huffman, Jon-Lark Kim, and Patrick Solé, editors, Concise Encyclopedia of Coding Theory. Chapman and Hall/CRC, 2021.
  * [DR73] P. Deligne and M. Rapoport.  Les schémas de modules de courbes elliptiques.  In Modular functions of one variable, II (Proc. Internat. Summer School, Univ. Antwerp, Antwerp, 1972), Lecture Notes in Math., Vol. 349, pages 143–316. Springer, Berlin, 1973.
  * [Duu08] Iwan M. Duursma.  Algebraic geometry codes: general theory.  In Advances in algebraic geometry codes, volume 5 of Ser. Coding Theory Cryptol., pages 1–48. World Sci. Publ., Hackensack, NJ, 2008.
  * [Ful89] William Fulton.  Algebraic curves.  Advanced Book Classics. Addison-Wesley Publishing Company Advanced Book Program, Redwood City, CA, 1989.  An introduction to algebraic geometry, Notes written with the collaboration of Richard Weiss, Reprint of 1969 original.
  * [Gop81] Valerii D. Goppa.  Codes on algebraic curves.  Dokl. Akad. Nauk SSSR, 259(6):1289–1290, 1981.  In Russian.
  * [GS95] Arnaldo Garcia and Henning Stichtenoth.  A tower of Artin-Schreier extensions of function fields attaining the Drinfeld-Vladut bound.  Inventiones Mathematicae, 121:211–222, 1995.
  * [Har77] Robin Hartshorne.  Algebraic geometry, volume 52 of Graduate Texts in Mathematics.  Springer-Verlag, New York, 1977.
  * [HP95] Tom Høholdt and Ruud Pellikaan.  On the decoding of algebraic–geometric codes.  IEEE Trans. Inform. Theory, 41(6):1589–1614, Nov 1995.
  * [HvLP98] Tom Høholdt, Jacobus Hendricus van Lint, and Ruud Pellikaan.  Algebraic geometry of codes.  In Handbook of coding theory, Vol. I, II, pages 871–961. North-Holland, Amsterdam, 1998.
  * [Iha81] Y. Ihara.  Some remarks on the number of rational points of algebraic curves over finite fields.  J. Fac. Sci. Univ. Tokyo Sect. IA Math., 28:721–724, 1981.
  * [KM85] Nicholas M. Katz and Barry Mazur.  Arithmetic moduli of elliptic curves, volume 108 of Annals of Mathematics Studies.  Princeton University Press, Princeton, NJ, 1985.
  * [Lor96] Dino Lorenzini.  An invitation to arithmetic geometry, volume 9 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 1996.
  * [LPS88] A. Lubotzky, R. Phillips, and P. Sarnak.  Ramanujan graphs.  Combinatorica, 8(3):261–277, 1988.
  * [Mar88] G. A. Margulis.  Explicit group-theoretic constructions of combinatorial schemes and their applications in the construction of expanders and concentrators.  Problemy Peredachi Informatsii, 24(1):51–60, 1988.
  * [Mil17] James S. Milne.  Modular functions and modular forms (v1.31), 2017.  Available at www.jmilne.org/math/.
  * [Mor90] Carlos J. Moreno.  Algebraic curves over finite fields.  Cambridge tracts in mathematics. Cambridge University Press, Cambridge, 1990.
  * [Sha94] Igor R. Shafarevich.  Basic algebraic geometry. 1.  Springer-Verlag, Berlin, second edition, 1994.
  * [Sil09] Joseph Silverman.  The arithmetic of elliptic curves, volume 106.  Springer-Verlag New York, second edition, 2009.
  * [Ste99] Serguei A. Stepanov.  Codes on algebraic curves.  Kluwer Academic/Plenum Publishers, New York, 1999.
  * [Sti09] Henning Stichtenoth.  Algebraic function fields and codes, volume 254 of Graduate Texts in Mathematics.  Springer-Verlag, Berlin, second edition, 2009.
  * [TVN07] Michael A. Tsfasman, Serge Vlăduţ, and Dmitry Nogin.  Algebraic geometric codes: basic notions, volume 139 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2007.
  * [TVZ82] Michael A. Tsfasman, Serge G. Vlăduţ, and Th. Zink.  Modular curves, Shimura curves, and Goppa codes, better than Varshamov-Gilbert bound.  Math. Nachr., 109:21–28, 1982.
  * [VD83] Sergei G. Vlăduţ and Vladimir G. Drinfeld.  Number of points of an algebraic curve.  Funct. Anal. Appl., 17:53–54, 1983.
  * [VM84] S. G. Vlăduţ and Yu. I. Manin.  Linear codes and modular curves.  In Current problems in mathematics, Vol. 25, Itogi Nauki i Tekhniki, pages 209–257. Akad. Nauk SSSR Vsesoyuz. Inst. Nauchn. i Tekhn. Inform., Moscow, 1984.
  * [Wal00] Judy L. Walker.  Codes and curves, volume 7 of Student Mathematical Library.  American Mathematical Society, Providence, RI; Institute for Advanced Study (IAS), Princeton, NJ, 2000.  IAS/Park City Mathematical Subseries.