Adjoint-based optimization of the
Rayleigh-Bénard instability with melting
boundary

Tomas Fullana¹,²*, Alejandro Quirós Rodríguez¹,
Vincent Le Chenadec³, Taraneh Sayadi¹,⁴*

¹Institut Jean le Rond d'Alembert, Sorbonne Université, CNRS, Paris,
F-75005, France.

²Laboratory of Fluid Mechanics and Instabilities, EPFL, Laussane,
CH-1015, Switzerland.

³MSME, Université Gustave Eiffel, UPEC, CNRS, Marne-la-Vallée,
F-77454, France.

⁴Mathematical and Numerical Modelling Laboratory, Conservatoire
National des Arts et Métiers, Paris, F-75003, France.

*Corresponding author(s). E-mail(s): tomas.fullana@gmail.com;
taraneh.sayadi@lecnam.net;

Abstract

In this work, we propose an adjoint-based optimization procedure to control the onset of the Rayleigh-Bénard instability with a melting front. A novel cut cell method is used to solve the Navier-Stokes equations in the Boussinesq approximation and the convection-diffusion equation in the fluid layer, as well as the heat equation in the solid phase. To track the interface we use the level set method where its evolution is simply governed by an advection equation. An incomplete continuous adjoint problem is then derived by treating the velocity field obtained from the forward problem as a known variable in the adjoint convection-diffusion equation, thereby avoiding the need to solve a Navier-Stokes adjoint in the fluid phase. To the best of our knowledge, this provides the first adjoint-based optimization framework for Rayleigh-Bénard instability with a melting boundary. Two optimization problems, together with a comparison against a derivative-free particle-swarm method, demonstrate that the proposed incomplete adjoint yields gradients accurate enough to control the front shape while reducing the number of expensive function evaluations by about an order of magnitude.
# 1 Introduction

Stefan problem [1] models transport and transfer phenomena, in particular solid-liquid phase change in evaporating or chemically reacting flows. Such phenomena govern the interface motion in many engineering related problems such as dendritic solidification [2, 3], phase transformation in metallic alloys [4], and solid fuel combustion [5]. At the interface, the Stefan condition arises from the interaction of liquid and solid phases (both considered incompressible), resulting in a moving liquid-solid interface (freezing or melting front). The speed of the front is directly related to the jump in the conductive heat flux across the interface. In one dimension, this problem has been studied in depth [6–8]; in higher dimensions, however, due to the unstable pattern formation [9–11] specific numerical methods, such as the level set method [2, 12–21], have been used to predict the motion and shape of the interface. In crystal growth, for example, under-cooling triggers an instability mechanism, causing the solid phase of the material to grow into the liquid phase in a finger-like or dendritic fashion, resulting in complex interfacial shapes and possible topography changes. Moreover, while convection with a fixed topography is usually studied in the classical Rayleigh-Bénard setup [22], many natural phenomena such as erosion or melting involve a coupling between the flow and the moving boundary [23, 24].

The shape of the resulting interface may strongly affect the outcome and time-
frame of the production processes in many industrial applications involving phase
change [25, 26]. As a result, it is desirable to extract efficient control strategies to
manipulate the motion of the interface, for instance, by tracking a prescribed trajec-
tory. Two major types of optimization methods in use today are (i) gradient-based, and
(ii) derivative free methods. While an efficient class of generic algorithms (belonging
to the class of derivative free methods) based on the surrogate management frame-
work [27] and artificial neural networks [28] have been used for optimization in fluid
mechanics, mainly in the area of aerodynamic shape optimization, they could require
many function evaluations, for training purposes for example. When detailed simu-
lations of interfacial flows are concerned, each function evaluation commands a full
(potentially unsteady) CFD computation, causing gradient-based methods to be at
an advantage. Adjoint-based algorithms present a suitable method, as they allow the
determination of the gradient at a cost comparable to a single function evaluation [29–
32]. The use of adjoint methods for design and optimization has been an active area of
research which started with the pioneering work of Pironneau [33] with applications
in fluid mechanics, and has been extensively used in aeronautical shape optimiza-
tion by Jameson and co-workers [34, 35]. Ever since these groundbreaking studies,
adjoint-based methods have been widely used in fluid mechanics particularly in the
areas of aero- and thermo-acoustics [36, 37]. More recently flow regimes dominated
by nonlinear dynamics, such as separation and mixing have also been analysed using
adjoint-based techniques [38–41]. Adjoint-based methods have also been employed for the purpose of sensitivity analysis or control in flows in the presence of large gradients (flames or interfaces) [42–45], showing great promise, and therefore are adopted here to carry out the optimization procedure.

In the context of Stefan problems various control strategies have been employed to track the location of the interface. In a one-dimensional setting, for example, set-valued fixed point equations [46] or linear-quadratic defect minimization [47] have been used to control the location of the front. Adjoint-based algorithms have also been applied to a Stefan problem caused by heterogeneous reactions on a surface of a one-dimensional solid particle [5] to extract sensitives with respect to various kinetic parameters. Alternatively, in a two-dimensional setting, adjoint-based algorithms have been utilized previously together with finite element and finite difference approaches to track and control the location of interface by imposing heat flux (or temperature) at the boundary in order to realize the desired interface motion [48–50]. In particular, Bernauer & Herzog [51], making use of shape calculus tools, derived the set of adjoint equations to extract control strategies for Stefan problems with a sharp representation of the interface.

Similar to the approach of [51] shape calculus tools have been employed to extract the corresponding adjoint equations in our previous work [17]. However, contrary to the previous studies, control strategies are extracted here to control the interface shape resulting from the melting boundary. For the configuration considered here, deriving the fully coupled analytical adjoint of the Navier–Stokes–Stefan system on a moving domain is not tractable, which motivates the use of an incomplete continuous adjoint tailored to this problem. The incomplete adjoint is derived – by considering the velocity field in the fluid layer as a known variable in the adjoint convection-diffusion equation – and used to minimize a tracking-type cost functional through a gradient-based algorithm. We will show that the derived incomplete adjoint is sufficient to retrieve meaningful information on the system and to drive the cost functional towards its minimum at a substantially reduced computational cost.

The main contributions of this work can be summarized as follows. (i) A sharp-interface optimal-control problem for the Rayleigh–Bénard instability with a melting boundary, based on a level-set/cut-cell discretization of the two-phase Stefan problem under the Boussinesq approximation is formulated. (ii) We then derive an incomplete continuous adjoint that neglects the Navier–Stokes adjoint by treating the velocity field from the forward simulation as a known variable in the adjoint convection-diffusion equation. (iii) We demonstrate, on two non-trivial optimization problems and in comparison with a derivative-free method, that the resulting gradients are sufficiently accurate to control the front shape while significantly reducing the number of function evaluations.

The present study is restricted to a two-dimensional configuration with a sharp interface representation and without surface-tension effects. The control is exerted only through the temperature distribution prescribed at the upper boundary of a periodic domain. As such, the proposed framework is intended as a first step towards adjoint-based control of melting fronts in convection-driven problems.
<figure><img src="image_1.png" /><figcaption>Fig. 1: Schematic of the melting boundary problem.</figcaption></figure>

The paper is organized as follows; in the first section we will define the forward problem and present validations cases, in the second section we will derive the continuous adjoint equations and in the last section we present the optimization results resulting from the solution of the minimization of a tracking-type cost functional.

## 2 Governing equations

Following the recent studies on the Rayleigh-Bénard instability in the presence of a melting boundary [12, 52], we consider the evolution of an initially flat liquid layer heated from below and comprised between a wall and a solid layer. As the solid melts, the liquid layer will grow vertically until a critical height where convection kicks in. A schematic of this configuration is shown in Fig. 1. The gravity is pointing downwards $g = -ge_y$ and the horizontal size of the domain is $bH$, where $b$ is the aspect ratio and where $H = 1$ is the height of the system. The temperature at the bottom wall is $T_b = 0.7$, and the one at the top wall is our control variable $w(x) < 0$, such that $w(x) < T_M < T_b$ where $T_M = 0$ is the melting temperature at the front. The choice of $T_b$ is such that $\Delta T = 1$ when $w(x) = -0.3$. The initial flat interface is located at $h_0 = 0.05$ and the initial temperature field is set to zero.

In the fluid phase, we solve the Navier-Stokes equations in the Boussinesq approximation coupled to the convection-diffusion equation for the temperature field. With this approximation, the variation of the liquid density with the temperature is taken into account in the buoyancy force only. This means that the density and thermal conductivities are assumed to be constant and equal. The equations read

$$ \mathrm{Pr}^{-1} \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = - \nabla p + \mathrm{Ra} T_L \mathbf{e}_y + \nabla^2 \mathbf{u}, \quad (1) $$

$$ \frac{\partial T_L}{\partial t} + \mathbf{u} \cdot \nabla T_L = \nabla^2 T_L, $$
where $T_L$ is the dimensionless temperature field in the liquid phase. The dimensionless numbers governing this equation are the Prandtl number

$$ Pr = \frac{\nu}{k}, \qquad (2) $$

defined by the ratio of the liquid kinematic viscosity $\nu$ and thermal diffusivity $k$, and the Rayleigh number

$$ Ra = \frac{g \alpha_t \Delta T H^3}{\nu k}, \qquad (3) $$

where $g$ is the constant gravitational acceleration, $\alpha_t$ the coefficient of thermal expansion, $H$ a characteristic length of the system and $\Delta T = T_b - T_M$. We apply no slip boundary conditions on the lower boundary and on the interface between the two phases. In the solid phase, the heat equation applies

$$ \frac{\partial T_S}{\partial t} = \nabla^2 T_S, \qquad (4) $$

where $T_S$ is the dimensionless temperature field in the solid phase. At the interface, both temperature fields are equal to the melting temperature

$$ T_L(\mathbf{x}, t) = T_S(\mathbf{x}, t) = T_M \quad \mathbf{x} \in \Gamma, \qquad (5) $$

where $\Gamma$ denotes the front position. In the present example, we disregard surface tension effects in the thermodynamic boundary condition at the front (ie. the Gibbs-Thomson relation).

Finally, as phase change occurs, there will appear a latent heat which is either absorbed or released. The condition of heat conservation at a given point on the moving interface corresponds to the rate at which heat is generated at the boundary, balanced by the rate at which this heat flows in either phase. Along the interface, the Stefan condition [53] states that

$$ v = \mathrm{St} [\nabla T]^S_L \cdot \mathbf{n}, \quad \mathbf{x} \in \Gamma, \qquad (6) $$

where $v$ is the interface velocity, $\mathbf{n}$ is the outward normal unit vector at the interface and $[\nabla T]^S_L = \partial T_S / \partial \mathbf{n} - \partial T_L / \partial \mathbf{n}$ is the jump in normal components of the temperature gradients from the solid phase to the liquid phase. The Stefan number $St$ is defined as the ratio between available heat in the system and the latent heat

$$ St = \frac{k \Delta T}{L_H}, \qquad (7) $$

where $L_H$ is the latent heat. The moving interface is captured using the level set method. A level set function $\phi$ [54] is defined on the computational domain $\Omega$ to map the locus of one of its iso-levels ($\{\mathbf{x} \in \Omega | \phi(\mathbf{x}, t) = \phi_0\}$) to an interface $\Gamma(t)$ that
separates two non-overlapping domains, $\Omega_S(t)$ and $\Omega_L(t)$, each occupied by a different phase. The value $\phi$ is defined as the signed distance to the interface,

$$ \phi(\mathbf{x}, t) = \begin{cases} -d(\mathbf{x}, \Gamma(t)), & \mathbf{x} \in \Omega_S(t) \\ 0, & \mathbf{x} \in \Gamma(t) \\ d(\mathbf{x}, \Gamma(t)), & \mathbf{x} \in \Omega_L(t) \end{cases}, \qquad (8) $$

where $d(\mathbf{x}, \Gamma(t))$ denotes the minimal distance between the point $\mathbf{x}$ and the interface $\Gamma(t)$,

$$ d(\mathbf{x}, \Gamma(t)) = \operatorname*{argmin}_{\mathbf{y} \in \Gamma(t)} \|\mathbf{x} - \mathbf{y}\|, \qquad (9) $$

with $\|\cdot\|$ denoting the Euclidean distance. The main advantage of the level set method is that its evolution is governed by a simple advection equation

$$ \frac{\partial \phi}{\partial t} + v \cdot \nabla \phi = 0. \qquad (10) $$

provided that the velocity field can be smoothly extended from $\Gamma(t)$. The most natural algorithm is to let $v$ be a constant along the lines normal to $\Gamma$. To achieve this, the method described in [55] is adopted here. Using this approach, the velocity is extended in the normal direction by solving the following hyperbolic partial differential equation

$$ \left\{ \begin{array}{l} \frac{\partial F}{\partial t^*} + S(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla F = 0 \quad \text{in } \Omega \\ F(x, 0) = v \quad \text{on } \Gamma \end{array} \right. \qquad (11) $$

where $F$ is the extended velocity field equal to $v$ at the front, $t^*$ denotes a pseudo-time and $S(\phi)$ is the signature function

$$ S(\phi) = \begin{cases} -1 & \text{if } \phi < 0 \\ 0 & \text{if } \phi = 0 \\ +1 & \text{if } \phi > 0 \end{cases} \qquad (12) $$

Equation 11 is then discretized using a first order upwind scheme and integrated in time by a forward Euler method until steady state. Taking $\mathbf{n}$ as the normal vector defined as

$$ n = (n_x, n_y) = \left( \frac{\phi_x}{\sqrt{(\phi_x^2 + \phi_y^2)}}, \frac{\phi_y}{\sqrt{(\phi_x^2 + \phi_y^2)}} \right), \qquad (13) $$

the discretisation leads to

$$ F_{ij}^{n+1} = F_{ij}^{n} - \tau^{*} \left( \left(S_{ij}n_{ij}^{x}\right)^{+} \frac{F_{ij} - F_{i-1j}}{\Delta} + \left(S_{ij}n_{ij}^{x}\right)^{-} \frac{F_{i+1j} - F_{ij}}{\Delta} + \left(S_{ij}n_{ij}^{y}\right)^{+} \frac{F_{ij} - F_{ij-1}}{\Delta} + \left(S_{ij}n_{ij}^{y}\right)^{-} \frac{F_{ij+1} - F_{ij}}{\Delta} \right) \qquad (14) $$

where $\Delta$ is the uniform grid spacing, $(x)^{+} = \max(0, x)$ and $(x)^{-} = \min(0, x)$, and the time step $\tau^*$ is chosen so that $\tau^*/\Delta^2 = 0.45$. The pseudo-time spawn in the velocity
extension algorithm is purely fictitious and the number of iterations in Equation 14 corresponds to the width of the narrow-band (NB) around the 0-level set where the velocities are initialized. A semi-implicit scheme described in [56, 57] is used here to solve the level set advection equation. Using this method, the equation is written in an alternative form

$$ \frac{\partial \phi}{\partial t} + F \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla \phi = 0, \qquad (15) $$

which is then divided into conservative and non-conservative terms

$$ \underbrace{\frac{\partial \phi}{\partial t} + \nabla \cdot \left( F \phi \frac{\nabla \phi}{|\nabla \phi|} \right)}_{a} - \underbrace{\phi \nabla \cdot \left( F \frac{\nabla \phi}{|\nabla \phi|} \right)}_{b} = 0, \qquad (16) $$

resulting in a second order partial differential equation akin to a weighted diffusion equation. The first term (a) has a diffusion coefficient $F\phi$ that depends on the solution and represents a nonlinear curvature flow whereas in the second term (b) the solution is multiplied by the curvature of its level-sets. The main idea behind this scheme is to distinguish two cases: if the product $F\phi$ is negative (positive, respectively) then a represents a forward (backward, respectively) diffusion and b represents a backward (forward, respectively) diffusion. The forward diffusion is treated implicitly while the backward diffusion is treated explicitly leading to a semi-implicit scheme with a dif- fusive CFL number (fixed to 0.5 in the test cases presented here).

In order to discretize Equation 16, we use the same notation as in [56]. Consider p to be a finite volume of a cell and $e_{pq}$ the edge between p and q, $q \in \mathcal{N}(p)$ where $\mathcal{N}(p)$ is the set of neighboring finite volumes. The length of $e_{pq}$ is normalized to 1. Let $\mathbf{n}_{pq}$ be the outer normal vector to $e_{pq}$ with respect to p. Finally, let us denote $\bar{\phi}_p$ the constant reconstruction of $\phi_p$ in the finite volume p and $\bar{\phi}_{pq}$ the constant reconstruction of $\phi$ on $e_{pq}$. Integrating volume p, Equation 16 yields

$$ \int_{p} \frac{\partial \phi}{\partial t} dx + \int_{p} \nabla \cdot \left( F \phi \frac{\nabla \phi}{|\nabla \phi|} \right) dx - \int_{p} \phi \nabla \cdot \left( F \frac{\nabla \phi}{|\nabla \phi|} \right) dx = 0. \qquad (17) $$

Applying Stokes' theorem and using the constant reconstructions of $\phi$, we obtain

$$ \int_{p} \frac{\partial \phi}{\partial t} dx + \sum_{q \in \mathcal{N}(p)} \bar{\phi}_{pq} \int_{e_{pq}} F \frac{1}{|\nabla \phi|} \nabla \phi \cdot \mathbf{n}_{pq} ds \\ - \sum_{q \in \mathcal{N}(p)} \bar{\phi}_{p} \int_{e_{pq}} F \phi \frac{1}{|\nabla \phi|} \nabla \phi \cdot \mathbf{n}_{pq} ds = 0. $$

Let us denote $|\nabla \phi_{pq}|$ the reconstructed Hamiltonian $|\nabla \phi|$ on the edge $e_{pq}$ and $(\phi_q - \phi_p)/1$ the normal derivative $\nabla \phi \cdot \mathbf{n}_{pq}$ on the same edge. We obtain the following expression

$$ \int_{p} \frac{\partial \phi}{\partial t} dx + \sum_{q \in \mathcal{N}(p)} \frac{F \bar{\phi}_{pq}}{|\nabla \phi_{pq}|} (\phi_q - \phi_p) $$
$$ - \sum_{q \in \mathcal{N}(\mathbf{p})} \frac{F \bar{\phi}_{\mathbf{p}}}{|\nabla \phi_{\mathbf{p}q}|} (\phi_{\mathbf{q}} - \phi_{\mathbf{p}}) = 0, $$

leading to

$$ \int_{\mathbf{p}} \frac{\partial \phi}{\partial t} dx + \sum_{q \in \mathcal{N}(\mathbf{p})} \frac{F(\bar{\phi}_{\mathbf{p}} - \bar{\phi}_{\mathbf{pq}})}{|\nabla \phi_{\mathbf{pq}}|} (\phi_{\mathbf{p}} - \phi_{\mathbf{q}}) = 0. \quad (18) $$

Looking at the term $F(\bar{\phi}_{\mathbf{p}} - \bar{\phi}_{\mathbf{pq}})$, we can distinguish two cases:

* If the term is positive, we have a 'forward diffusion' or inflow towards the cell.

* If the term is negative, we have a 'backward diffusion' or outflow from the cell.

We therefore define the diffusion coefficient $a_{\mathbf{pq}}$ as

$$ a_{\mathbf{pq}} = \frac{F(\bar{\phi}_{\mathbf{p}} - \bar{\phi}_{\mathbf{pq}})}{|\nabla \phi_{\mathbf{pq}}|}, \quad (19) $$

and the related dominant forward and backward diffusion parts as

$$ a_{\mathbf{pq}}^{f} = \max(a_{\mathbf{pq}}, 0), \quad a_{\mathbf{pq}}^{b} = \min(a_{\mathbf{pq}}, 0). \quad (20) $$

Using a backward Euler time discretization, taking the forward contribution explicitly and the backward contribution implicitly, Equation 18 gives the following linear system

$$ \underbrace{\phi_{\mathbf{p}}^{n} + \frac{\tau}{\Delta^2} \sum_{q \in \mathcal{N}(\mathbf{p})} a_{\mathbf{pq}}^{f} (\phi_{\mathbf{p}}^{n} - \phi_{\mathbf{q}}^{n})}_{\text{implicit}} = \phi_{\mathbf{p}}^{n-1} + \frac{\tau}{\Delta^2} \underbrace{\sum_{q \in \mathcal{N}(\mathbf{p})} a_{\mathbf{pq}}^{b} (\phi_{\mathbf{p}}^{n-1} - \phi_{\mathbf{q}}^{n-1})}_{\text{explicit}} \quad (21) $$

where $\tau$ is the time step, $\Delta$ the uniform grid spacing and $n$ a given time step. To reconstruct of $\bar{\phi}_{\mathbf{p}}$, $\bar{\phi}_{\mathbf{pq}}$ and $|\nabla \phi_{\mathbf{pq}}|$, we use the diamond-cell strategy described in [57]. This evolution equation however does not preserve the signed distance property, which can be restored by periodically employing a reinitialization procedure by solving an Eikonal equation. There exist many numerical methods for solving this equation and, here we have adopted that of Min [58], also recently used in the work of Limare [12]. For further details on the numerical implementation of the level set coupled to a cut cell method, we refer the reader to [17, 19, 20].

The forward problem, setting $St = 1$, can be recast in the level set framework as follows.
Find a function $T: \Omega \times [0, t_f] \to \mathbb{R}$ and a function $\phi: \Omega \times [0, t_f] \to \mathbb{R}$ such that

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial T_S}{\partial t} = \nabla^2 T_S & \text{in } \Omega_S(t) \quad (\text{FP.a}) \\[2ex]
\displaystyle \text{Pr}^{-1} \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \text{Ra} T_L \mathbf{e}_y + \nabla^2 \mathbf{u} & \text{in } \Omega_L(t) \quad (\text{FP.b}) \\[2ex]
\displaystyle \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega_L(t) \quad (\text{FP.c}) \\[2ex]
\displaystyle \frac{\partial T_L}{\partial t} + \mathbf{u} \cdot \nabla T_L = \nabla^2 T_L & \text{in } \Omega_L(t) \quad (\text{FP.d}) \\[2ex]
\displaystyle \frac{\partial T_L}{\partial t} = \nabla^2 T_L & \text{in } \Omega_L(t) \quad (\text{FP.e}) \\[2ex]
\displaystyle T(x, 0) = T_0(x) & \text{in } \Omega_L(t) \quad (\text{FP.f}) \quad (\text{FP}) \\[2ex]
\displaystyle T_S = w(x) & \text{on } \partial \Omega_S \quad (\text{FP.g}) \\[2ex]
\displaystyle T_L = T_b & \text{on } \partial \Omega_L \quad (\text{FP.h}) \\[2ex]
\displaystyle T(x, t) = T_M & \text{on } \Gamma(t) \quad (\text{FP.i}) \\[2ex]
\displaystyle \frac{\partial \phi}{\partial t} = -[\nabla T_i]_L^S \cdot \nabla \phi & \text{on } \Gamma(t) \quad (\text{FP.j}) \\[2ex]
\displaystyle \phi(x, 0) = \phi_0(x) & \text{in } \Omega_L(t) \quad (\text{FP.k})
\end{array}
\right.
$$

The global Rayleigh number in Eq. (3) will control the onset of the Rayleigh-Bénard instability. The effective Rayleigh number is defined as

$$ \mathrm{Ra}_e = \mathrm{Ra} (1 - T_M) \bar{h}^3, \qquad (22) $$

where $\bar{h}$ is the average fluid height

$$ \bar{h}(t) = \frac{1}{b} \int_{0}^{b} h(x, t) dx. \qquad (23) $$

When the effective Rayleigh number reaches the critical value $Ra_c = 1707.76$, the initial diffusion-driven motion is transformed into a convection-driven one. We consider a domain with aspect ratio $b = 4$, a fixed temperature $w = -0.3$ at the top boundary and vary the global Rayleigh number. Figure 2a shows the average height as function of time for the different cases. In the $Ra = 10^4$ case, the critical Rayleigh number is not reached and the motion of the fluid layer is not affected, thus remaining a diffusion-driven one, similarly to the planar motion. In the rest of the cases, the effective Rayleigh number $Ra_e$ (Fig. 2b) crosses the threshold indicating the onset of
the instability, characterized by a sharp increase in interface speed. Figure 3 shows a time series of the temperature field and interface position for the Ra = 10⁵ case and Figure 4 shows the vorticity field for the same case. When the critical Rayleigh number is reached, the first bifurcation appears, creating the convection cells. The size of the convection cells will then vary with the secondary bifurcations mechanism. When the averaged height $\bar{h}$ matches the characteristic wavelength of the convection rolls, the convection cells have sufficient time to merge and then stabilize. We also note that the interface is deformed according to the shape of the cells.

<figure><img src="image_2.png" /><figcaption>Fig. 2: Average height and effective Rayleigh number as a function of time for Ra = 10⁵, 8 × 10⁴, 4 × 10⁴, 10⁴.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>(a) t = 0.1</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) t = 0.2</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) t = 0.3</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) t = 0.4</figcaption></figure>

**Fig. 3:** Times series of the temperature field and interface position for the Ra = 10⁵ case. The color map is the temperature field in both phases and the interface corresponding to the 0-level set is denoted in red.
<figure><img src="image_7.png" /><figcaption>Fig. 4: Times series of the vorticity and interface position for the Ra = 10⁵ case. The interface corresponding to the 0-level set is denoted in red.</figcaption></figure>
### 3 Adjoint problem

We now aim at solving a minimization problem by using the adjoint of (FP). We assume that the desired temperature field $T^d$ and the desired position of the interface $\phi^d$ are known a priori. The control variable $w$ is then iteratively calculated through an optimization procedure to reach these desired quantities and will drive our cost functional towards its minimum. The following tracking-type cost functional provides a mathematical description of the control goals stated above

$$ \mathcal{J}(T, \phi, w) = \frac{\beta_1}{2} \int_{\Omega} |T^f - T^d|^2 dx + \frac{\beta_2}{2} \int_{\Gamma^f} |\phi^f - \phi^d|^2 ds + \frac{\beta_3}{2} \int_{0}^{t_f} \int_{\partial\Omega_S} |w|^2 dsdt \quad (24) $$

with constants $\beta_1$ to $\beta_3$ acting as weights, $T^f$ and $\phi^f$ the final temperature field and level set function of (FP) and $\Gamma^f$ the final position of the interface (0-level set of $\phi^f$). The first term monitors the final temperature distribution and is mostly used as an initializer for the adjoint temperature field for the adjoint problem. The second term controls the relative position of the level set function with respect to the desired one. Taking advantage of the signed distance function property of both $\phi^d$ and $\phi^f$, the discrete form is simply the square of the difference between both functions computed in the points belonging to the final interface $\Gamma^f$. Finally, the last term penalizes the control cost and leads to the gradient equation in the adjoint problem.

We assume that each choice of the control variable $w$ leads to unique states $T(w)$ and $\phi(w)$. Therefore, the minimization problem (MP) reads

$$ \min_{w} \mathcal{J}(T(w), \phi(w), w) \qquad (MP) \\ \text{subject to (FP).} $$

We will now proceed with the derivation of the corresponding continuous adjoint problem. The main novelty behind our derivation is that the velocity field $u$, known at each time step from the forward solution, is treated as a known variable in the adjoint equations. This allows us to discard the adjoint of the Navier-Stokes equations in the fluid phase and leads to an incomplete adjoint formulation. The minimization procedure will serve as a validation of this hypothesis.

Let $\Theta$ be the adjoint temperature and $\psi$ the adjoint level set function. In order to compute the gradient of $\mathcal{J}(w)$, we introduce the following Lagrange functional $\mathcal{L}(dx, dt)$ and $ds$ are omitted for brevity)

$$ \begin{aligned} \mathcal{L}(T, \phi, w, \Theta, \Theta_S^D, \Theta_L^D, \Theta^I, \psi) = & \mathcal{J}(T, \phi, w) \\ & - \int_0^{t_f} \int_{\Omega_S} \left( \frac{\partial T_S}{\partial t} - \nabla^2 T_S \right) \Theta_S - \int_0^{t_f} \int_{\Omega_L} \left( \frac{\partial T_L}{\partial t} + \mathbf{u} \cdot \nabla T_L - \nabla^2 T_L \right) \Theta_L \\ & - \int_0^{t_f} \int_{\partial \Omega_S} (T - w) \Theta_S^D - \int_0^{t_f} \int_{\partial \Omega_L} (T - T_b) \Theta_L^D - \int_0^{t_f} \int_{\Gamma} (T - T_M) \Theta^I \\ & - \int_0^{t_f} \int_{\Gamma} \left( \frac{\partial \phi}{\partial t} + [\nabla T_i]_L^S \cdot \nabla \phi \right) \psi. \end{aligned} \quad (25) $$
Similarly to the temperature field, $\Theta_{S,L}: (x,t) \to \mathbb{R}^+$ denotes the adjoint temperature field in either phase. The Lagrange multipliers $\Theta_{S,L}^D$ and $\Theta^I$ are used for the boundary conditions on the domain and the interface respectively.

The adjoint system is obtained by setting to zero the derivatives of the Lagrange functional with respect to $T$ and $\phi: \mathcal{L}_T(\cdot) = \mathcal{L}_\phi(\cdot) = 0$. Setting the initial conditions to

$$ 
\begin{aligned}
T(x, 0) &= T_0(x), \\
\phi(x, 0) &= \phi_0(x),
\end{aligned} 
\qquad (26) $$

leads to $\delta T(x, 0) = \delta \phi(x, 0) = 0$ when calculating the derivatives in the direction $\delta T$ and $\delta \phi$.

We start by deriving the adjoint temperature equations by setting $\mathcal{L}_T(\cdot) = 0$. The terms in the Lagrangian that do not depend on the temperature vanish and we move the spatial and temporal derivatives towards the adjoint state $\Theta$. By applying integration by parts, once with respect to time and twice with respect to space, we obtain

$$ 
\begin{aligned}
\mathcal{L}_T \delta T = & \int_0^{t_f} \int_{\Omega_S} \left( \frac{\partial \Theta_S}{\partial t} + \nabla^2 \Theta_S \right) \delta T + \int_0^{t_f} \int_{\Omega_L} \left( \frac{\partial \Theta_L}{\partial t} - \mathbf{u} \nabla \Theta_L + \nabla^2 \Theta_L \right) \delta T \\
& - \int_{\Omega_S} \left( \Theta_S^f - \beta_1 (T^f - T^d) \right) \delta T^f - \int_{\Omega_L} \left( \Theta_L^f - \beta_1 (T^f - T^d) \right) \delta T^f \\
& + \int_0^{t_f} \int_{\partial \Omega_S} \left( -\frac{\partial \Theta_S}{\partial n} \delta T + \frac{\partial \delta T}{\partial n} (\Theta_S - \Theta_S^D) \right) + \int_0^{t_f} \int_{\partial \Omega_L} \left( -\frac{\partial \Theta_L}{\partial n} \delta T + \frac{\partial \delta T}{\partial n} (\Theta_L - \Theta_L^D) \right) \\
& + \int_0^{t_f} \int_{\Gamma} \left( \Theta_S \vec{v}_S \cdot n - \Theta_L \vec{v}_L \cdot n - \frac{\partial \Theta_S}{\partial n} + \frac{\partial \Theta_L}{\partial n} - \Theta^I \right) \delta T \\
& + \int_0^{t_f} \int_{\Gamma} (\Theta - \psi |\nabla \phi|) [\nabla \delta T]_L^S \cdot n \quad \forall \delta T = 0.
\end{aligned}
 $$

with $\vec{v}_S$ and $\vec{v}_L$ the velocities of the control volumes $\Omega_S$ and $\Omega_L$ respectively, only non-zero on $\Gamma$. The second to last term, corresponds to the normal jump in adjoint temperature, which gives $\Theta^I = -[\nabla \Theta]_L^S \cdot n$ as $\Theta_S = \Theta_L$ on $\Gamma$. By altering the directions of variations and eliminating certain terms, we obtain the adjoint temperature problem
$$
\left\{
\begin{array}{ll}
-\frac{\partial \Theta_S}{\partial t} = \nabla^2 \Theta_S & \text{in } \Omega_S(t) \\
-\frac{\partial \Theta_L}{\partial t} - \mathbf{u} \nabla \Theta_L = \nabla^2 \Theta_L & \text{in } \Omega_L(t) \\
\Theta(x, t_f) = \beta_1 (T^f - T^d) & \text{in } \Omega \\
\frac{\partial \Theta_S}{\partial n} = 0 & \text{on } \partial \Omega_S \\
\frac{\partial \Theta_L}{\partial n} = 0 & \text{on } \partial \Omega_L \\
\Theta(x, t) = \psi |\nabla \phi| & \text{on } \Gamma(t) \\
\Theta^I = -[\nabla \Theta]^S_L \cdot n & \text{on } \Gamma(t)
\end{array}
\right.
\qquad (\text{AT})
$$

The first two equations of the adjoint temperature problem (AT) are the heat and convection-diffusion equations in reverse time. The third equation is the initial condition for the adjoint temperature field, that depends on the desired temperature distribution $T^d$. The fourth and fifth equations are the homogeneous Neumann boundary conditions for the adjoint field mapped from the Dirichlet boundary conditions of the forward problem. The second to last equation is the Dirichlet boundary condition at the interface that now depends on the value of the adjoint level set function $\psi$. We can note that the adjoint level set function will not behave as a signed distance function but as an auxiliary variable that acts on the temperature field through the boundary condition at the interface $\Gamma$. Finally, the last equation states that the multiplier $\Theta^I$, defined on $\Gamma$, must be equal to the jump in normal gradient of $\Theta$.

Turning now to the adjoint level set problem; due to the geometric non-linearity induced by $\phi$ on the domains of integration a careful treatment of each term is required as described in [51]. Denoting $D[[H(\phi); \delta\phi]]$ the variation of $H(\cdot)$ in the direction $\delta\phi$, the adjoint level set equations are derived by setting $\mathcal{L}_\phi(\cdot) = 0$

$$
\mathcal{L}_{\phi} \delta \phi = D[[\mathcal{J}(T, \phi, w); \delta \phi]] \\
+D\left[\int_{0}^{t_f} \int_{\Omega_S} \left(\frac{\partial T_S}{\partial t} - \nabla^2 T_S\right) \Theta_S; \delta \phi\right] + D\left[\int_{0}^{t_f} \int_{\Omega_L} \left(\frac{\partial T_L}{\partial t} + \mathbf{u} \cdot \nabla T_L - \nabla^2 T_L\right) \Theta_L; \delta \phi\right] \\
-D\left[\int_{0}^{t_f} \int_{\partial \Omega} (T_S - w) \Theta_S^D; \delta \phi\right] - D\left[\int_{0}^{t_f} \int_{\partial \Omega} (T_L - T_b) \Theta_L^D; \delta \phi\right] \\
-D\left[\int_{0}^{t_f} \int_{\Gamma} (T - T_M) \Theta^I; \delta \phi\right] - D\left[\int_{0}^{t_f} \int_{\Gamma} \left(\frac{\partial \phi}{\partial t} + [\nabla T_i]^S_L \cdot \nabla \phi\right) \psi; \delta \phi\right].
$$
We divide the contributions of the adjoint level set term-by-term. The contribution from the cost functional, removing the terms that do not depend on $\phi$, simplifies to

$$D[\mathcal{J}(T, \phi, w); \delta\phi] = D\left[\frac{\beta_2}{2} \int_{\Gamma^f} |\phi^f|^2; \delta\phi\right]$$

By using the theorem on the derivative of a boundary integral (Theorem A), we obtain

$$D[\mathcal{J}(T, \phi, w); \delta\phi] = -\frac{\beta_2}{2} \int_{\Gamma^f} \frac{\delta\phi^f}{|\nabla\phi|} \left( \frac{\partial|\phi^f|^2}{\partial n} + \kappa|\phi^f|^2 \right).$$

By specifying that $\delta\phi = 0$ on $\delta\Omega$, the second and third terms are then equal to zero. For the fourth term, we need to assume that the adjoint level set $\phi$ and the multiplier $\Theta^I$ are defined on all of $\Omega$. Under such assumptions, we can apply the same boundary integral theorem (Theorem A). Moreover, using the fact that $T_M$ is constant along the interface, this term simplifies to

$$-D\left[\int_{0}^{t_f} \int_{\Gamma} (T - T_M) \Theta^I; \delta\phi\right] = -\int_{0}^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \frac{\partial T}{\partial n} \Theta^I.$$

Finally, by applying the same boundary integral theorem together with the chain rule, the last term simplifies to

$$\begin{aligned} & -D\left[\int_{0}^{t_f} \int_{\Gamma} \left(\frac{\partial\phi}{\partial t} + [\nabla T_i]^S_L \cdot \nabla\phi\right) \psi; \delta\phi\right] \\ &= -\int_{0}^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \left(\frac{\partial}{\partial n}\left(\left(\frac{\partial\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\phi\right)\psi\right) + \kappa\left(\frac{\partial\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\phi\right)\psi\right) \\ &\quad + \int_{0}^{t_f} \int_{\Gamma} \left(\frac{\partial\delta\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\delta\phi\right) \psi \\ &= -\int_{0}^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \frac{\partial}{\partial n}\left(\frac{\partial\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\phi\right) \psi + \int_{0}^{t_f} \int_{\Gamma} \left(\frac{\partial\delta\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\delta\phi\right) \psi. \end{aligned}$$

By using the assumption of constant velocity $v$ extended in the normal direction on all of $\Omega$, we obtain that

$$-D\left[\int_{0}^{t_f} \int_{\Gamma} \left(\frac{\partial\phi}{\partial t} + [\nabla T_i]^S_L \cdot \nabla\phi\right) \psi; \delta\phi\right] = -\int_{0}^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \left(\frac{\partial\delta\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\delta\phi\right) \psi.$$

In order to obtain the adjoint Stefan condition, we need to move the derivatives from $\delta\phi$ to $\psi$. Assuming that $\psi$ is defined on all of $\Omega$, we can use the corollary B on integration by parts in time on a moving surface

$$-\int_{0}^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \left(\frac{\partial\phi}{\partial t} + [\nabla T]^S_L \cdot \nabla\phi\right) \psi = -\int_{\Gamma^f} \delta\phi^f \psi^f + \int_{\Gamma^0} \delta\phi^0 \psi^0$$
$$ 
\begin{aligned}
& -\int_0^{t_f} \int_{\Gamma} \left( \delta\phi \frac{\partial\psi}{\partial t} + \nabla(\delta\phi\psi) \cdot v + \delta\phi\psi \operatorname{div}_{\Gamma} v + \psi[\nabla T]_L^S \cdot \nabla\delta\phi \right) \\
& = -\int_{\Gamma_f} \delta\phi^f \psi^f + \int_{\Gamma^0} \delta\phi^0 \psi^0 - \int_0^{t_f} \int_{\Gamma} \delta\phi \frac{\partial\psi}{\partial t} \\
& -\int_0^{t_f} \int_{\Gamma} \delta\phi (\nabla(\psi) \cdot v + \psi \operatorname{div}_{\Gamma} v) + \int_0^{t_f} \int_{\Gamma} (\psi[\nabla T]_L^S \cdot \nabla\delta\phi - \psi\nabla\delta\phi \cdot v).
\end{aligned}
$$

The last term cancels due to the fact that $[\nabla T]_L^S = v$ on $\Gamma$. The resulting terms of the adjoint level set, using $\delta\phi^0 = 0$, give

$$ 
\begin{aligned}
\mathcal{L}_{\phi}\delta\phi = & -\frac{\beta_2}{2} \int_{\Gamma_f} \frac{\delta\phi^f}{|\nabla\phi|} \frac{\partial|\phi^f|^2}{\partial n} + \kappa|\phi^f|^2 \\
& -\int_0^{t_f} \int_{\Gamma} \frac{\delta\phi}{|\nabla\phi|} \frac{\partial T}{\partial n} \Theta^I \\
& -\int_{\Gamma_f} \delta\phi^f \psi^f - \int_0^{t_f} \int_{\Gamma} \delta\phi \frac{\partial\psi}{\partial t} \\
& -\int_0^{t_f} \int_{\Gamma} \delta\phi (\nabla(\psi) \cdot v + \psi \operatorname{div}_{\Gamma} v).
\end{aligned}
$$

By regrouping the terms by their domains of integration and requiring $\mathcal{L}_{\phi}\delta\phi = 0 \forall \delta\phi$, we obtain that

$$ 
\begin{aligned}
\psi^f = & -\frac{1}{|\nabla\phi|} \frac{\beta_2}{2} \left( \frac{\partial|\phi^f|^2}{\partial n} + \kappa|\phi^f|^2 \right) \\
-\frac{\partial\psi}{\partial t} = & \nabla\psi \cdot v + \psi \operatorname{div}_{\Gamma} v + \frac{1}{|\nabla\phi|} \frac{\partial T}{\partial n} \Theta^I.
\end{aligned}
$$

This equation can be interpreted as a first-order conservation law on $\Gamma$ and the source term on the right-hand-side can be extended in the neighborhood of $\Gamma$. For details on the derivation and simplifications of $\mathcal{L}_{\phi}(\cdot) = 0$, we refer the reader to [17].

Finally, the last term to complete the adjoint Stefan problem is the gradient equation. By setting $\mathcal{L}_w\delta w = 0$ in Equation 25, we have

$$ \mathcal{L}_w \delta w = \left[ \left[ \frac{\beta_3}{2} \int_0^{t_f} \int_{\partial\Omega_S} |w|^2; \delta w \right] - \left[ \int_0^{t_f} \int_{\partial\Omega_S} (T_S - w) \Theta_S^D; \delta w \right] \right]. $$

Using the previously defined multiplier $\Theta_S^D$ that is identically equal to $\Theta_S$ on $\partial\Omega_S$ we obtain the gradient equation $0 = \beta_3 w + \Theta$ on $\partial\Omega_S$. Putting it all together, we can now formulate the adjoint problem as follows
Find a function $\Theta : \Omega \times [t_f, 0] \to \mathbb{R}$ and a function $\psi : \Omega \times [t_f, 0] \to \mathbb{R}$ such that

$$
\left\{
\begin{array}{ll}
-\frac{\partial \Theta_S}{\partial t} = \nabla^2 \Theta_S & \text{in } \Omega_S(t) \quad (\text{AP.a}) \\
-\frac{\partial \Theta_L}{\partial t} - \mathbf{u} \nabla \Theta_L = \nabla^2 \Theta_L & \text{in } \Omega_L(t) \quad (\text{AP.b}) \\
\Theta(x, t_f) = \beta_1(T(t_f) - T_{t_f}) & \text{in } \Omega \quad (\text{AP.c}) \\
\frac{\partial \Theta_S}{\partial n} = 0 & \text{on } \partial \Omega_S \quad (\text{AP.d}) \\
\frac{\partial \Theta_L}{\partial n} = 0 & \text{on } \partial \Omega_L \quad (\text{AP.e}) \quad (\text{AP}) \\
\Theta(x, t) = \psi |\nabla \phi| & \text{on } \Gamma(t) \quad (\text{AP.f}) \\
\frac{\partial \psi}{\partial t} + \text{div}(\psi v) = \frac{1}{|\nabla \phi|} \frac{\partial T}{\partial n} [\nabla \Theta^I]_L^S \cdot n & \text{on } \Gamma(t) \quad (\text{AP.g}) \\
\psi(x, t_f) = -\frac{\beta_2}{2} \left( \frac{\partial}{\partial n} |\phi^f|^2 + \kappa |\phi^f|^2 \right) & \text{in } \Omega \quad (\text{AP.h}) \\
0 = \beta_4 w + \Theta & \text{on } \partial \Omega \quad (\text{AP.i})
\end{array}
\right.
$$

We again emphasis on the fact that, here, the velocity field **u** is considered as a known variable, eliminating the need to derive the adjoint of the Navier-Stokes equations n the fluid phase and leading to an incomplete adjoint formulation.

4 Optimization results

Adjoint-based methods—being efficient gradient-based optimization approaches—are often the methods of choice in PDE-constrained optimal control. By solving an adjoint equation once, one obtains the gradient of the objective with respect to all control variables at minimal cost, making the optimal-control formulation both tractable and efficient even for large-scale parameter spaces. We choose to solve the minimization problem (MP) by using the limited memory BFGS (L-BFGS) method, a quasi-Newton method originally described in [59]. The main characteristic of this method is that it determines the descent direction by preconditioning the gradient with an approximation of the Hessian matrix. This information is obtained using past approximations – the number of approximations is determined by the memory length parameter which is set to *m* = 10 – as well as the gradient. As an initial guess for the initial Hessian, we use the scaled identity matrix as described in [60]. Algorithm 1 summarizes the L-BFGS method used in our numerical example. The algorithm is stopped at a given iteration *n* if one of the following criteria is fulfilled

• The relative difference in control variable $\left| \frac{w^{n+1} - w^n}{w^n} \right| < 10^{-8}$.
* The relative difference in cost functional $\left| \frac{\mathcal{J}^{n-1} - \mathcal{J}^n}{\mathcal{J}^n} \right| < 10^{-8}$. This criterion can be relaxed to allow temporary increase of the cost functional, for example to 'escape' a local minimum.

* The norm of the gradient $|\nabla \mathcal{J}^n| < 10^{-6}$.

<caption>**Algorithm 1:** Optimization procedure using the L-BFGS method</caption>

input : $w^0$, $m = 10$
output: $w$, $T$, $\phi$, $\Theta$, $\psi$
$\mathcal{J}^k = \beta_3 w^k + \Theta^k$

k ← 0, l ← 0

**while** not converged **do**

Solve the forward Stefan problem (FP) for $T^k$ and $\phi^k$

Solve the adjoint Stefan problem (AP) for $\Theta^k$ and $\psi^k$

Compute the gradient:

$$ \nabla \mathcal{J}^k = \beta_3 w^k + \Theta^k $$

**if** $k \ge 1$ **then**

$s^{k-1} = \sigma^{k-1} d^{k-1}$, $g^{k-1} = \nabla \mathcal{J}^k - \nabla \mathcal{J}^{k-1}$

**if** $(s^{k-1})^T g^{k-1} \le 0$ **then**

| $l \leftarrow 0$

**else if** $(s^{k-1})^T g^{k-1} > 0$ **then**

$l \leftarrow l+1$

**if** $l > m$ **then**

| Remove $\{s^{l-m}, g^{l-m}\}$

**end**

Add $\{s^{l-m}, s^{l-m}\}$

**end**

**end**

Choose an initial approximation to the inverse of the Hessian $H_0^k$

Construct the direction $d^k = -H^k \nabla \mathcal{J}^k$

Determine $\sigma^k$ using a Line Search algorithm with backtracking where
$$ \sigma^k = \operatorname*{argmin}_{\mathcal{J}} (w^k + \sigma^k d^k) $$

Update $w^{k+1} = w^k + \sigma^k d^k$

$k \leftarrow k+1$

**end**
To test our gradient-based minimization procedure, we now consider the numerical setup defined in 1. The goal of this optimization procedure is to actuate on the top boundary condition to diminish the fluid layer growth, therefore reducing the onset of the instability. The actuator w is parametrized with the following basis

$$w = -|a_1| - |a_2|(1 - \tanh(2x)^2), \quad (27)$$

where x corresponds to the bounds of the domain and $a_1$ and $a_2$ to the basis coefficients. Through the optimization process, the amplitude of each coefficient is determined using the gradient equation (AP).i. Practically, we fit the temperature values on the boundary with a given basis, and update the coefficient at each iteration of the optimization procedure. By opting for a parameterized distribution we ensure the smoothness of the actuation function. Indeed, due to the high sensitivity of the cost functional with respect to the basis considered – too many parameters will create multiple local minima – the number of parameters are kept at a low enough value to ensure the convexity of the problem while allowing spatial variation of the actuation function.

The global Rayleigh number is set to $Ra = 10^5$ and an initial guess $w = 0$. The coefficients in the cost functional (24) are set to $\beta_1 = 1$, $\beta_2 = 1$ and $\beta_3 = 10^{-3}$. The desired solution is computed beforehand with coefficients $a_1 = 0.3$ and $a_2 = 2$. Figure 6 shows the interface position and temperature field at different iterations of the optimization procedure. The initial guess with $w = 0$ allows the interface to grow unperturbed, creating self-similar convection cells that affect the interface shape. The localized heat flux at the top boundary will create a heat gradient at the center of the domain leading to a destabilization of the interface and the appearance of the convection-driven regime before the unperturbed case (Fig. 5).

As the optimization procedure advances, we observe that the final interface position converges to the desired one, with a localized heat flux that prevents the melting interface growth at the center of the domain. In Figure 7 we show the evolution of the normalized cost functional $J/J_0$ as a function of the iterations. The procedure stops at iteration 15, with a relative value three orders of magnitude lower, as the norm of the gradient goes below the prescribed threshold and the basis coefficients have reached the desired ones (Fig. 8).

We now show a second optimization example, with a more complex basis

$$w = \sum_{n=1}^{4} a_n \sin^n(2\pi x) + \sum_{n=1}^{4} a_{n+4} \cos^n(2\pi x). \quad (28)$$

The global Rayleigh number is identical to the previous test case. The initial guess is 0 and the desired solution is such that the interface remains flat, with a constant heat flux $w = -1$. As shown in Figs. 9 and 10, the optimization procedure is able to recover the desired solution even with a higher number of fitting coefficients.

We now compare the results to a derivative-free method (where no information on the gradient is required), the Particle Swarm optimization algorithm [61]. This method attempts to improve global convergence by switching between four evolutionary states:
<figure><img src="image_8.png" /><figcaption>**Fig. 5:** Average height and effective Rayleigh number as a function of time for the desired solution and initial guess.</figcaption></figure>

exploration, exploitation, convergence, and jumping out. In the jumping out state it tries to take the best particle and move it away from its local optimum, to improve the ability to find a global one. We test the derivative-free optimization procedure on the cases presented using the same parameters for the cost functional. Table 1 summarizes the comparison between both methods. As expected, the derivative-free one attains a lower minimum at a cost of a much higher number of function evaluations (one order of magnitude higher). The L-BFGS algorithm driven by the incomplete adjoint thus attains objective values close to those of the particle-swarm optimizer at a fraction of the computational cost.

<table><caption>**Table 1:** Comparison between the L-BFGS and the Particle swarm method for the two considered cases.</caption><thead><tr><th rowspan="2">Case</th><th colspan="3">L-BFGS</th><th colspan="2">Particle Swarm</th></tr><tr><th>J calls</th><th>∇J calls</th><th>J<sub>final</sub>/J<sub>0</sub></th><th>J calls</th><th>J<sub>final</sub>/J<sub>0</sub></th></tr></thead><tbody><tr><td>Case 1</td><td>31</td><td>16</td><td>1.02 × 10<sup>-3</sup></td><td>928</td><td>2.04 × 10<sup>-4</sup></td></tr><tr><td>Case 2</td><td>42</td><td>16</td><td>1.01 × 10<sup>-3</sup></td><td>2321</td><td>2.14 × 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 6: Iterations 0, 1, 2 and 15 of the optimization procedure. The red curve represents the desired shape and the black one the final position of the interface at a given iteration. The color map corresponds to the final temperature field of the (FP).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7:** Normalized cost functional as a function of the iterations of the L-BFGS optimization procedure in the first test case.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8:** Basis coefficients as a function of the iterations of the L-BFGS optimization procedure in the first test case.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9:** Normalized cost functional as a function of the iterations of the L-BFGS optimization procedure in the second test case.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10:** Basis coefficients as a function of the iterations of the L-BFGS optimization procedure in the second test case.</figcaption></figure>
# 5 Conclusion

This paper demonstrates the feasibility of extending the optimization of Stefan problems to scenarios where the flow in the fluid phase is significant and cannot be ignored. By treating the velocity field as a known variable in the adjoint equations, we were able to derive an incomplete continuous adjoint formulation for the two-phase problem. Our results indicate that the adjoint-based gradient descent method effectively controls the shape of melting fronts in this highly nonlinear system, while greatly reducing the number of function evaluations compared with a derivative-free particle-swarm algorithm. The present formulation is limited to two-dimensional, laminar flows with a prescribed top boundary temperature in a periodic domain and an incomplete adjoint in which the Navier-Stokes adjoint problem is not solved. While the numerical examples indicate that the resulting gradients are accurate enough for practical optimization, the method does not capture sensitivities associated with changes in the flow field itself. This approach nevertheless shows great potential for applications in more complex and dynamic melting processes. Future work could explore more complex actuator basis, non-periodic domains, and three-dimensional systems to further enhance the optimization framework and broaden its applicability.

# Funding

There is no source of funding.

# Author contributions

T.F. and T.S. conceived the study. T.F. performed the numerical simulations. T.F. wrote the original draft with feedback from all authors. T.S. supervised the study.

# Data availability

The data used in this study are available from the corresponding authors upon request.

# Replication of results

The results are obtained with an in-house Julia code Flower.jl
https://github.com/flnt/Flower.jl/.

# Ethics approval and Consent to participate

Not applicable.

# Conflict of interest

The authors declare no competing interests.

# A Shape calculus theorems
**Theorem A** (Derivative of boundary integral). Let $J(\Omega) = \int_{\Gamma} f ds$ be a boundary integral, the derivative is given by

$$dJ(\Omega; V) = \left. \left( \frac{d}{d\lambda} \int_{\Gamma_{\lambda}} f ds_{\lambda} \right) \right|_{\lambda=0} \\ = - \int_{D} \frac{\delta\phi}{|\nabla\phi|} \left( \frac{\partial f}{\partial n} + f\kappa \right) ds$$

where $V$ is the velocity field and $\kappa$ is the mean curvature of $\Gamma$.

**Theorem B** (Surface transport theorem). Let $f(\cdot, t) : S_t \to \mathbb{R}$ be a scalar field defined on the moving surface $S_t$. Then

$$\frac{d}{dt} \int_{S_t} f(x, t) dS_t = \int_{S_t} \dot{f}(x, t) + f(x, t) \operatorname{div} S_t w(x, t) dS_t$$

where $w$ is the normal velocity of the moving surface $S_t$ and $\dot{f}$ is the parameter-time derivative of $f$. If $f(\cdot, t)$ is the restriction of a function $\hat{f}(\cdot, t)$ to $S_t$, then

$$\dot{f}(x, t) = w(x, t) \cdot \nabla \hat{f}(x, t) + \frac{\partial}{\partial t} \hat{f}(x, t)$$

**Corollary B** (Integration by Parts in Time on a Moving Surface).

$$\int_{0}^{T} \int_{S_t} \hat{g}(x, t) \hat{h}_t(x, t) dS_t dt = \int_{S_T} g(x, T) h(x, T) dS_T - \int_{S_0} g(x, 0) h(x, 0) dS_0 \\ - \int_{0}^{T} \int_{S_t} \hat{g}_t(x, t) \hat{h}(x, t) + w(x, t) \cdot \nabla (\hat{g}(x, t) \hat{h}(x, t)) \\ + g(x, t) h(x, t) \operatorname{div}_{S_t} w dS_t dt$$

where $g$ and $h$ are restrictions of $\hat{g}$ and $\hat{h}$ to $S_t$.

## References

[1] Šarler, B.: Stefan's work on solid-liquid phase changes. Engineering Analysis with Boundary Elements 16(2), 83–92 (1995) https://doi.org/10.1016/0955-7997(95) 00047-X

[2] Osher, S., Sethian, J.A.: Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J Comput. Phys. 79(1), 12–49 (1988) https://doi.org/10.1016/0021-9991(88)90002-2

[3] Juric, D., Tryggvason, G.: A front-tracking method for dendritic solidification. J Comput. Phys. 123, 127–148 (1996)
[4] Segal, G., Vuik, C., Vermolen, F.: A conserving discretization for the free boundary in a two-dimensional stefan problem. J Comput. Phys. **141**, 1–21 (1998)

[5] Hassan, A., Sayadi, T., Le Chenadec, V., Pitsch, H., Attili, A.: Adjoint-based sensitivity analysis of steady char burn out. Com. Theo. Modeling **25**(1), 96–120 (2021) https://doi.org/10.1080/13647830.2020.1838614

[6] Javierre, E., Vuik, C., Vermolen, F.J., Zwaag, S.: A comparison of numerical models for one-dimensional Stefan problems. J Comp. Applied Mathematics **192**(2), 445–459 (2006) https://doi.org/10.1016/j.cam.2005.04.062

[7] Rose, M.E.: An enthalpy scheme for stefan problems in several dimensions. App. Numerical Math. **12**, 229 (1993)

[8] Brattkus, K., Meiron, D.I.: Numerical simulations of unsteady crystal growth. SIAM J. Appl. Math. **52**, 1303 (1992)

[9] Langer, J.S.: Instabilities and pattern formation in crystal growth. Rev. Mod. Phys. **52**(1), 1–28 (1980) https://doi.org/10.1103/RevModPhys.52.1

[10] Mullins, W.W., Sekerka, R.F.: Stability of a planar interface during solidification of a dilute binary alloy. Journal of Applied Physics **35**(2), 444–451 (1964) https://doi.org/10.1063/1.1713333. Publisher: American Institute of Physics

[11] Woods, A.W.: Melting and dissolving. J. Fluid Mech. **239**, 429 (1992) https://doi.org/10.1017/S0022112092004476

[12] Limare, A., Popinet, S., Josserand, C., Xue, Z., Ghigo, A.: A hybrid level-set / embedded boundary method applied to solidification-melt problems. arXiv (2022) arXiv:2202.08300 [math.NA]

[13] Gibou, F., Fedkiw, R., Caflisch, R., Osher, S.: A Level Set Approach for the Numerical Simulation of Dendritic Growth. Journal of Scientific Computing **19**(1), 183–199 (2003) https://doi.org/10.1023/A:1025399807998

[14] Chen, H., Min, C., Gibou, F.: A numerical scheme for the Stefan problem on adaptive Cartesian grids with supralinear convergence rate. Journal of Computational Physics **228**(16), 5803–5818 (2009) https://doi.org/10.1016/j.jcp.2009.04.044

[15] Jenkins, N., Maute, K.: Level set topology optimization of stationary fluid-structure interaction problems. Struct Multidisc Optim **52**(1), 179–195 (2015) https://doi.org/10.1007/s00158-015-1229-9. Accessed 2025-07-19

[16] Bayat, E., Egan, R., Bochkov, D., Sauret, A., Gibou, F.: A sharp numerical method for the simulation of Stefan problems with convective effects. Journal of Computational Physics **471**, 111627 (2022) https://doi.org/10.1016/j.jcp.2022.111627
[17] Fullana, T.: Simulation and optimization of complex phenomena in multi-phase flows. phdthesis, Sorbonne Université (December 2022). https://theses.hal.science/tel-04053531

[18] Silva, K.E.S., Sivapuram, R., Ranjbarzadeh, S., Gioria, R.S., Silva, E.C.N., Picelli, R.: Topology optimization of stationary fluid-structure interaction problems including large displacements via the TOBS-GT method. Struct Multidisc Optim 65(11), 1-18 (2022) https://doi.org/10.1007/s00158-022-03442-3. Company: Springer Distributor: Springer Institution: Springer Label: Springer Number: 11 Publisher: Springer Berlin Heidelberg. Accessed 2025-07-19

[19] Fullana, T., Le Chenadec, V., Sayadi, T.: Adjoint-based optimization of two-dimensional Stefan problems. Journal of Computational Physics 475, 111875 (2023) https://doi.org/10.1016/j.jcp.2022.111875

[20] Quirós Rodríguez, A., Fullana, T., Sayadi, T., Le Chenadec, V.: A level-set based cut cell method for interfacial flows. Part 2: Free-surface flows and dynamic contact angle treatment. Acta Mechanica (2024)

[21] Castañar, I., Baiges, J., Codina, R.: Topology optimization of incompressible structures subject to fluid-structure interaction. Struct Multidisc Optim 67(6), 90 (2024) https://doi.org/10.1007/s00158-024-03770-6. Accessed 2025-07-19

[22] Bodenschatz, E., Pesch, W., Ahlers, G.: Recent Developments in Rayleigh-Bénard Convection. Annual Review of Fluid Mechanics 32 (Volume 32, 2000), 709-778 (2000) https://doi.org/10.1146/annurev.fluid.32.1.709. Publisher: Annual Reviews Type: Journal Article

[23] Yang, R., Howland, C.J., Liu, H.-R., Verzicco, R., Lohse, D.: Morphology evolution of a melting solid layer above its melt heated from below. Journal of Fluid Mechanics 956, 23 (2023) https://doi.org/10.1017/jfm.2023.15

[24] Perissutti, D., Marchioli, C., Soldati, A.: Morphodynamics of melting ice over turbulent warm water streams (2024). https://arxiv.org/abs/2406.12116

[25] Picelli, R., Ranjbarzadeh, S., Sivapuram, R., Gioria, R.S., Silva, E.C.N.: Topology optimization of binary structures under design-dependent fluid-structure interaction loads. Struct Multidisc Optim 62(4), 2101-2116 (2020) https://doi.org/10.1007/s00158-020-02598-0. Accessed 2025-07-19

[26] Siqueira, L.O., Cortez, R.L., Sivapuram, R., Ranjbarzadeh, S., Gioria, R.d.S., Silva, E.C.N., Picelli, R.: Topology optimization for stationary fluid-structure interaction problems with turbulent flow via sequential integer programming and smooth explicit boundaries. Advances in Engineering Software 190, 103599 (2024) https://doi.org/10.1016/j.advengsoft.2024.103599. Accessed 2025-07-19
[27] Marsden, A.L., Feinstein, J.A., Taylor, C.A.: A computational framework for derivative-free optimization of cardiovascular geometries. Computer Methods in Applied Mechanics and Engineering 197(21-24), 1890-1905 (2008)

[28] Pierret, S., Coelho, R.F., Kato, H.: Multidisciplinary and multiple operating points shape optimization of three-dimensional compressor blades. Structural and Multidisciplinary Optimization 33(1), 61-70 (2007)

[29] Giles, M.B., Pierce, N.A.: An Introduction to the Adjoint Approach to Design. Flow Turb. Combus. 65(3), 393-415 (2000)

[30] Griewank, A., Walther, A.: Algorithm 799: revolve: an implementation of check-pointing for the reverse or adjoint mode of computational differentiation. ACM Transactions on Mathematical Software (TOMS) 26(1), 19-45 (2000). Publisher: ACM New York, NY, USA

[31] Rodrigues, S.S., Marta, A.C.: Adjoint-based shape sensitivity of multi-row turbomachinery. Struct Multidisc Optim 61(2), 837-853 (2020) https://doi.org/10.1007/s00158-019-02386-5. Accessed 2025-07-19

[32] Burghardt, O., Gomes, P., Kattmann, T., Economon, T.D., Gauger, N.R., Palacios, R.: Discrete adjoint methodology for general multiphysics problems. Struct Multidisc Optim 65(1), 28 (2022) https://doi.org/10.1007/s00158-021-03117-5. Accessed 2025-07-19

[33] Pironneau, O.: On optimum design in fluid mechanics. J Fluid Mech. 64(1), 97-110 (1974)

[34] Jameson, A.: Aerodynamic design via control theory. J Scientific Computing 3(3), 233-260 (1988)

[35] Jameson, A., Martinelli, L., Pierce, N.A.: Optimum Aerodynamic Design Using the Navier-Stokes Equations. Theor. Comp. Fluid Dyn. 10(1), 213-237 (1998)

[36] Juniper, M.: Triggering in the horizontal Rijke tube: non-normality, transient growth and bypass transition. J. Fluid Mech. 667, 272-308 (2010)

[37] Lemke, M., Reiss, J., Sesterhenn, J.: Adjoint-based analysis of thermoacoustic coupling. ICNAAM, 2163-2166 (2013)

[38] Schmidt, S., Ilic, C., Schulz, V., Gauger, N.R.: Three-dimensional large-scale aerodynamic shape optimization based on shape calculus. AIAA J 51(11), 2615-2627 (2013)

[39] Rabin, S., Caulfield, C., Kerswell, R.: Designing a more nonlinearly stable laminar flow via boundary manipulation. J Fluid Mech. 738, 1-112 (2014)

[40] Foures, D., Caulfield, C., Schmid, P.: Optimal mixing in two-dimensional plane
Poiseuille flow at finite Péclet number. J Fluid Mech. **748**, 241–277 (2014)

[41] Duraisamy, K., Alonso, J.: Adjoint-based techniques for uncertainty quantification in turbulent flows with combustion. 42nd AIAA Fluid Dynamics Conference and Exhibit, 25–28 (2012)

[42] Fikl, A., Le Chenadec, V., Sayadi, T.: Control and optimization of interfacial flows using adjoint-based techniques. Fluids **5**(3) (2020) https://doi.org/10.3390/fluids5030156

[43] Ou, K., Jameson, A.: Unsteady adjoint method for the optimal control of advective and burger's equations using high-order spectral difference method. American Institute of Aeronautics and Astronautics, 1–18 (2011) https://doi.org/10.2514/6.2011-24

[44] Braman, K., Oliver, T.A., Raman, V.: Adjoint-based sensitivity analysis of flames. Comb. Theo. Modelling **19**(1), 29–56 (2015)

[45] Lemke, M., Cai, L., Reiss, J., Pitsch, H., Sesterhehn, J.: Adjoint-based sensitivity analysis of quantities of interest of complex combustion models. Comb. Theo. Modelling **23**(1), 180–196 (2019)

[46] Hoffmann, K.H., Sprekels, J.: Real-time control of the free boundary in a two-phase stefan problem. Numerical Functional Analysis and Optimization **5**(1), 47–76 (1982)

[47] Knabner, P.: Control of stefan problems by means of linear-quadratic defect minimization. Numerische Mathematik **46**(3), 429–442 (1985)

[48] Kang, S., Zabaras, N.: Control of the freezing interface motion in two-dimensional solidification processes using the adjoint method. International J Numerical Methods in Engineering **38**, 63–80 (1995)

[49] Yang, Z.: The adjoint method for the inverse design of solidification processes with convection. PhD thesis, Cornell University (1997)

[50] Hinze, M., Ziegenbalg, S.: Optimal control of the free boundary in a two-phase stefan problem. J Comput. Phys. **223**(2), 657–684 (2007)

[51] Bernauer, M.K., Herzog, R.: Optimal Control of the Classical Two-Phase Stefan Problem in Level Set Formulation. SIAM J. Sci. Comput. **33**(1), 342–363 (2011) https://doi.org/10.1137/100783327

[52] Favier, B., Purseed, J., Duchemin, L.: Rayleigh–Bénard convection with a melting boundary. J. Fluid Mech. **858**, 437–473 (2019) https://doi.org/10.1017/jfm.2018.773

[53] Stefan, J.: Über die theorie der eisbildung, insbesondere über die eisbil- dung im
polarmeere. Ann. Physik Chemie **42**, 269–286 (1891)

[54] Sethian, J.A.: Level Set Methods and Fast Marching. Cambridge Monographs on Applied and Computational Mathematics, vol. 1. Cambridge University Press, University of California, Berkeley (1999)

[55] Peng, D., Merriman, B., Osher, S., Zhao, H., Kang, M.: A PDE-Based Fast Local Level Set Method. J Comp. Phys. **155**(2), 410–438 (1999) https://doi.org/10.1006/jcph.1999.6345

[56] Mikula, K., Ohlberger, M., Urbán, J.: Inflow-implicit/outflow-explicit finite volume methods for solving advection equations. Applied Numerical Mathematics **85**, 16–37 (2014) https://doi.org/10.1016/j.apnum.2014.06.002

[57] Mikula, K., Ohlberger, M.: A New Level Set Method for Motion in Normal Direction Based on a Semi-Implicit Forward-Backward Diffusion Approach. SIAM J. Sci. Comput. **32**(3), 1527–1544 (2010) https://doi.org/10.1137/09075946X

[58] Min, C.: On reinitializing level set functions. J Comp. Phys. **229**(8), 2764–2772 (2010) https://doi.org/10.1016/j.jcp.2009.12.032

[59] Liu., D.C., Nocedal, J.: On the limited memory method for large scale optimization. Mathematical Programming **3**(45), 503–528 (1989)

[60] Wright, S.J., Nocedal, J.: Numerical Optimization. Springer, ??? (2006)

[61] Zhan, Z.-H., Zhang, J., Li, Y., Chung, H.S.-H.: Adaptive particle swarm optimization. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) **39**(6), 1362–1381 (2009) https://doi.org/10.1109/TSMCB.2009.2015956