# On the well-posedness of a Hele-Shaw-like system resulting from an inverse
geometry problem formulated through a shape optimization setting

Julius Fergy Tiongson Rabago  Masato Kimura

###### Abstract

The purpose of this study is twofold. First, we revisit a shape optimization
reformulation of a prototypical shape inverse problem and briefly propose a
simple yet efficient numerical approach for solving the minimization problem.
Second, we examine the existence, uniqueness, and continuous dependence of a
classical solution to a Hele-Shaw-like system, derived from the continuous
setting of a numerical discretization of the shape optimization reformulation
for the shape inverse problem. The analysis is based on the methods developed
by G. I. Bizhanova and V. A. Solonnikov in “On Free Boundary Problems for
Second Order Parabolic Equations” (Algebra Anal. 12 (6) (2000) 98-139) and
also by V. A. Solonnikov in “Lectures on Evolution Free Boundary Problems:
Classical Solutions” (Lect. Notes Math., Springer, 2003, pp. 123-175).

###### keywords:

Hele-Shaw problem, quasi-stationary Stefan-type problem, shape inverse
problem, shape optimization

††journal: a journal

[label1]organization=Faculty of Mathematics and Physics, Institute of Science
and Engineering, Kanazawa University, addressline=Kakumamachi, city=Kanazawa,
postcode=920-1192, state=Ishikawa, country=Japan

###### Contents

  1. 1 Introduction
    1. 1.1 A Hele-Shaw-like system
    2. 1.2 Derivation of the system
    3. 1.3 Extended-regularized normal flows for shape identifications
  2. 2 Auxiliary Results
  3. 3 Main Results
  4. 4 Problem Transformation onto a Fixed Domain
  5. 5 Regularity of Solutions on a Fixed Domain
  6. 6 The Nonlinear Problem
  7. 7 The Linear Problem
  8. 8 Proof of the Main Result
  9. A Comparison of normal derivatives
    1. A.1 The case of concentric circles
    2. A.2 The case of concentric spheres
  10. B Mullins-Sekerka analysis for 2D axisymmetric domains
  11. C Change of variables
  12. D Computations of $\delta\textsf{A}^{(0)}$ and $\delta B_{0}$
  13. E Lemmata Proofs
    1. E.1 Proof of Lemma 2.5
    2. E.2 Proof of Lemma 4.2
    3. E.3 Proof of Lemma 7.2
    4. E.4 Proof of Lemma 8.1
    5. E.5 Interpolation inequalities

##  1 Introduction

###  1.1 A Hele-Shaw-like system

In this study, we examine a Hele-Shaw-like system that arises from a shape
optimization reformulation of an inverse geometry problem. We prove the
existence, uniqueness, and continuous dependence of a classical solution to
the quasi-stationary moving boundary problem derived from the discretized
version of the optimization algorithm used to solve the shape optimization
problem. The inverse problem occurs in contexts such as non-destructive
testing to identify inclusions or voids in solids. The available information
or data are collected from a part of the solid’s surface (geometrical
boundary), with the goal of determining the inclusion within the interior of
the solid.

Let $T$ be a positive and number and $t\in[0,T]$. Let
$\Omega(t)\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, be a bounded annular domain
with sufficiently smooth boundary $\partial\Omega(t)=\Gamma(t)\cup\Sigma$
consisting of two non-intersecting boundaries/surfaces, where $\Sigma$ is
fixed while $\Gamma(t)$ is a free boundary. We assume that $\Sigma$ lies
outside the surface $\Gamma(t)$, and denote by $\nu{(t)}$ the unit normal
vector to $\Gamma(t)$ directed inward the domain $\Omega(t)$; that is, it goes
in the same direction with the outward unit normal vector to $\Sigma$. The
main problem we consider here is the following: given a pair of positive-
valued functions $f=f(x,t)$ and $g=g(x,t)$ defined on $\Sigma$, for $t>0$, and
an initial geometry $\Gamma_{\circ}$ of the moving boundary, we seek to find a
surface $\Gamma(t)$ and two functions $u_{\text{D}}=u_{\text{D}}(x,t)$ and
$u_{\text{N}}=u_{\text{N}}(x,t)$ satisfying the following system of partial
differential equations (PDEs):

$\left\\{\begin{array}[]{rcll}-\Delta u_{\text{D}}(x,t)&=&0,&\quad
x\in\Omega(t),\quad t\in[0,T],\\\ u_{\text{D}}(x,t)&=&f(x,t),&\quad
x\in\Sigma,\quad t\in[0,T],\\\ u_{\text{D}}(x,t)&=&0,&\quad
x\in\Gamma(t),\quad t\in[0,T],\\\ -\Delta u_{\text{N}}(x,t)&=&0,&\quad
x\in\Omega(t),\quad t\in[0,T],\\\
\dfrac{\partial{}}{\partial\nu}u_{\text{N}}(x,t)&=&g(x,t),&\quad
x\in\Sigma,\quad t\in[0,T],\\\ u_{\text{N}}(x,t)&=&0,&\quad
x\in\Gamma(t),\quad t\in[0,T],\\\
V_{n}(x,t)&=&-\dfrac{\partial{}}{\partial\nu}\left(u_{\text{D}}(x,t)-u_{\text{N}}(x,t)\right)&\quad
x\in\Gamma(t),\quad t\in[0,T],\\\
\Gamma(0)&=&\Gamma_{\circ},&\end{array}\right.$ (1)

where $\partial/\partial{\nu}=\partial/\partial{\nu(t)}$ is the inward normal
derivative operator on $\Gamma(t)$ while $V_{n}(x,t)$ represents the velocity
of movement of $\Gamma(t)$ in the direction of the normal $\nu{(t)}$ for all
$t>0$. As alluded above, the problem originates from a reformulation of a
shape inverse problem into a shape optimization setting. Shape optimization is
a well-established technique extensively applied in engineering sciences for
addressing shape identification problems. The modern mathematical theory of
shape optimization was laid in monographs [1, 2, 3].

Motivation and review of related work. As far as we are concerned, shape
inverse problems are typically solved numerically from an optimization
perspective through shape optimization settings; see, for example, [4, 5] and
the references therein. The model problem from which the Hele-Shaw-like system
(1) was derived is a specific case of the more general conductivity
reconstruction problem. This type of problem is known to be severely ill-posed
in the sense of Hadamard [4]. Nevertheless, it has been widely studied in
literature, both theoretically and numerically; see, e.g., [4, 6, 7, 8, 9, 10,
11, 12, 13] and references therein. For instance, the issue of the existence
and uniqueness of the solution to the problem from boundary measurement data
has been studied by several authors; see, e.g., [8, 9, 10, 13]. Aside from
these aforementioned investigations, shape optimization reformulations of
shape inverse problems, in general, are hardly examined from a different
theoretical and numerical point of view. The main purpose of this
investigation, therefore, is to carry out a rigorous analysis of the
existence, uniqueness, and continuous dependence on the data of the classical
solution of (1) in short-time horizon. We assert that, to our knowledge,
little to no work has been done on the well-posedness of the shape
optimization problem from which (1) is derived, especially in the direction of
our current study. We believe that the system (1) itself, originating from a
shape optimization context, which in turn was derived from a shape inverse
problem, is novel. Therefore, the analysis conducted in this work, inspired by
the work of Bizhanova and Solonnikov [14] and Solonnikov [15], offers a new
perspective. We emphasize, however, that Escher and Simonett demonstrated the
existence and uniqueness of classical solutions for the multidimensional
expanding Hele-Shaw problem in [16]. In their study, the problem is formulated
over the exterior boundary component, while the interior component remains
fixed. It is also worth mentioning that geometric flows related to shape
optimization problems similar to Bernoulli-type free boundary problems are
explored by Cardaliaguet and Ley in [17]. The evolution process examined in
their work involves a combination of curvature and Hele-Shaw type nonlocal
terms. They introduce the concept of generalized set solutions, drawing
heavily from viscosity solutions (see, e.g., [18] or [19, Chap. 2, p. 69]),
which diverges from our current focus. The main result in [17] is a principle
of inclusion preservation for generalized solutions, ensuring the existence,
uniqueness, and stability of solutions. Furthermore, they prove that the
solutions converge to a generalized Bernoulli exterior free-boundary problem.
The study also investigates the asymptotic behavior of the flow, demonstrating
the convergence of solutions towards a generalized Bernoulli exterior free-
boundary problem. Meanwhile, from geometric analysis perspective, Plotnikov
and Sokołowski [20] recently provided a comprehensive review of existing
results in shape optimization emphasizing the theory of gradient flows for
objective functions and their regularizations. For the sake of simplicity,
their mathematical exposition is confined to two spatial dimensions, and they
illustrate their findings through a model problem, leading to an initial
result concerning convergence in shape optimization. In an earlier work [21],
Plotnikov and Sokołowski rigorously established the well-posedness of the
Cauchy problem arising from a penalized gradient flow and conducted an in-
depth examination of the case involving the Kohn-Vogelius cost functional.

######  Remark 1.1.

The analysis conducted in this investigation can be applied to a class of
shape optimization problems. Specifically, the same methodology presented in
the subsequent sections can be utilized to prove the well-posedness of (1),
where the normal velocity of the boundary of the evolving domain is given by
$V_{n}(x,t)=-G(x,u(x))$ for $x\in\Gamma(t)$ and $t>0$. Here, $G$ can be the
exact shape gradient of a cost functional such as (3) and is assumed to be
positive for all $x\in\Omega(t)$, for all $t>0$. Additionally, this approach
is applicable to other quasi-stationary moving boundary problems that arise
from the discretization of numerical shape optimization schemes used to solve
the shape optimal control reformulation of the inverse geometry problem, such
as $L^{2}$ tracking of boundary data, under appropriate assumptions on the
kernel of the shape gradient of the associated cost functional.

In the next subsection, we will explore the derivation of system (1) and its
connection to a shape optimization problem to provide further context and
motivation.

###  1.2 Derivation of the system

Consider a bounded (simply connected) domain $D$ with boundary
$\partial{D}=\Sigma$ and assume the existence of a simply connected open set
${{\omega}}$ such that $\overline{{\omega}}\subset D$, composed of a material
with a constant conductivity that is essentially different from the constant
conductivity of the material in the annular subregion
$\Omega:=D\setminus\overline{{\omega}}$. We are interested in determining the
shape and location of the inclusion from the knowledge of the Cauchy data of
the electrical potential $u$ on the boundary $\Sigma$. In other words, through
the concept of non-destructive testing and evaluation, we want to identify the
inclusion given that the pair of boundary data $f=u\big{|}_{\Sigma}$ and
$g=(\nabla u\cdot\nu)\big{|}_{\Sigma}$ is known, where $\nu$ is the outward
unit normal to $\Sigma$. As alluded in the previous section, the problem under
consideration is a special case of the general conductivity reconstruction
problem and has been extensively studied in the literature, particularly in
the context of inverse problems.

Various numerical techniques can solve the inverse problem, including applying
shape optimization methods as demonstrated by Roche and Sokołowski in [5], and
by Eppler and Harbrecht in [4]. The former study utilized first-order shape
optimization algorithms to analyze and numerically address the proposed shape
optimization formulation. On the other hand, [4] investigated second-order
methods related to the previous study, developed and applied by the same
authors in their prior work. Here, we introduce an alternative numerical
approach closely aligned with the shape optimization techniques employed in
these references. Notably, our method does not directly require knowledge of
the exact shape gradient typically used in shape optimization procedures.
Nonetheless, our approach can be seen as stemming from the specific shape
optimization formulation in [4, 5]. In essence, our method simplifies the
gradient-based optimization procedure outlined in [5]. In fact, we will choose
a suitable descent direction that improves descent speed for the relevant
shape optimization problem. It is important to highlight in advance that our
method operates as a Lagrangian-type numerical scheme, utilizing finite
element methods. This is in contrast to the approaches used in [4, 5] which
rely on the concepts of boundary integral equations.

$\Sigma$$\omega$$\boldsymbol{D}$$\Omega$$\Gamma$$\nu$ Figure 1: Conceptual
model

To further explain the idea about our method, we first review the shape
optimization formulation proposed in [5] and [4] for the model problem
introduced previously; see Figure 1 for a conceptual model. To determine the
inclusion ${{\omega}}$ bounded by $\Gamma$, we apply the concept of
nondestructive testing and evaluation, a well-known technique used in
engineering and related sciences to evaluate the properties of a material
without causing damage. That is, we identify ${{\omega}}$ by measuring, for a
given current distribution $g\in H^{-1/2}(\Sigma)$, the voltage distribution
$f\in H^{1/2}(\Sigma)$ at the boundary $\Sigma$. This means that we are
seeking to find a connected domain $\Omega:=D\setminus\overline{{\omega}}$ and
an associated harmonic function $u$ which satisfies the overdetermined
boundary value problem:

$-\Delta u=0\ \text{in $\Omega$},\qquad u=0\ \text{on $\Gamma$},\qquad
u=f\quad\text{and}\quad\nabla u\cdot{\nu}=g\ \text{on $\Sigma$}.$ (2)

Equation (2) admits a solution only when the true inclusion (cavity or
void)111While the terms inclusion, cavity, and void are related and commonly
used in various contexts, such as materials science, geology, and engineering,
we have used them interchangeably in this discussion. ${{\omega}}$ is
considered.

We recall in the following theorem a key result about identifiability for this
inverse problem. It guarantees the uniqueness of the inclusion ${{\omega}}$
and consequently, the potential $u$ from a pair of boundary data
$(f,g)\neq(0,0)$.

######  Theorem 1.2 ([10]).

The Cauchy pair $(f,g)\neq(0,0)$ uniquely determine $\Gamma$ and $u$
satisfying (2).

In the context of thermal imaging, the Dirichlet data $f$ represent a
prescribed surface heat, and correspondingly, $g$ stands for the surface heat
flux of the material. This perspective reinterprets the state variable $u$ as
the internal thermal distribution within the material. It is important to note
that the homogeneous Dirichlet boundary condition applied to the unknown shape
$\Sigma$ models a perfectly conducting inclusion. Note that when the Dirichlet
data on the surface $\Sigma$ is positive, the solution $u$ to (2) is also
positive due to the maximum principle and unique continuation property [22].
This implies that the outgoing flux $\nabla u\cdot\nu$ on $\Sigma$ is
positive, while on $\Gamma$, the outgoing flux is negative. Moreover, because
$u$ takes homogeneous Dirichlet data on the surface $\Sigma$ of the unknown
inclusion, we have the identity $\left|{\nabla u}\right|\equiv-\nabla
u\cdot\nu$ on $\Sigma$. For a non-conducting inclusion, a Neumann boundary
condition governs the problem. Additionally, if the conductivity of the
inclusion differs from that of the surrounding material or object $\Omega$, an
inverse transmission problem arises.

Roche and Sokołowski in [5] proposed a shape optimization reformulation of (2)
as follows:

$J(\Omega)=\int_{\Omega}{\left|{\nabla(u_{\text{D}}-u_{\text{N}})}\right|^{2}}\,{{d}}x=\int_{\Sigma}{\left(g-\frac{\partial{u_{\text{D}}}}{\partial{\nu}}\right)\left(u_{\text{N}}-f\right)}\,{{d}}s\to\inf,$
(3)

subject to the underlying well-posed state problems

$\left\\{\begin{aligned} -\Delta u_{\text{D}}=0\quad\text{in $\Omega$},\qquad
u_{\text{D}}=f\quad\text{on $\Sigma$},\qquad u_{\text{D}}=0\quad\text{on
$\Gamma$},\\\ -\Delta u_{\text{N}}=0\quad\text{in $\Omega$},\qquad\nabla
u_{\text{N}}\cdot\nu=g\quad\text{on $\Sigma$},\qquad
u_{\text{N}}=0\quad\text{on $\Gamma$}.\end{aligned}\right.$ (4)

In (3), the infimum has to be taken over all domains containing an inclusion
with sufficiently regular boundary. The existence of optimal solutions with
respect to the shape optimization problem (3) was established in [5].

######  Remark 1.3.

  1. 1.

The proof of the existence of an optimal solution to (3) over some set of
admissible domains is based on properties of harmonic functions in
$\mathbb{R}^{2}$ and uses the same arguments as those given by Šverák [23] for
a slightly different shape optimization problem. In particular, the method
cannot be directly applied in $\mathbb{R}^{d}$ for $d\geqslant 3$.

  2. 2.

To establish an existence result in dimension $d=3$, we need to assume that
the family of admissible domains, say $U_{\text{ad}}$, meets the following
compactness criterion (see, e.g., [2, Chap. 2]): given any sequence
$\\{\Omega_{n}\\}_{n\in\mathbb{N}}\subset U_{\text{ad}}$, there exists a
subsequence
$\\{\Omega_{n_{k}}\\}_{k\in\mathbb{N}}=\\{D\setminus\overline{\omega}_{n_{k}}\\}_{k\in\mathbb{N}}$,
such that the sequence of compacts $\\{\omega_{n_{k}}\\}_{k\in\mathbb{N}}$
converges in the Hausdorff metric to the compact $\omega$. Moreover, the
associated sequence of metric projections $P_{n}:H_{0}^{1}(D)\to
H_{0}^{1}(\Omega_{n})$ converges strongly to the metric projection
$P:H_{0}^{1}(D)\to H_{0}^{1}(D\setminus\overline{\omega})$. For such a family
of admissible domains, there is a solution to the shape identification problem
in $\mathbb{R}^{3}$. The hypothesis generally requires some form of uniform
regularity for the boundaries $\Gamma$, such as the cone condition or more
intricate conditions related to capacity (refer to [1, 2]). For further
details, see Henrot, Horn, and Sokołowski [24] for an overview of stability
results for the Dirichlet problem, and Bucur and Zolesio [25] for an approach
to existence problems. Additionally, [2] provides more discussion on the
existence of optimal shapes in shape optimization problems. For a recent
related result on the existence issue related to Problem (3), see [26].

  3. 3.

For technical purposes, in order to establish the existence optimal solution,
we require that any admissible inclusion or obstacle $\omega\Subset D$
satisfies the condition that there exists a constant $d_{0}>0$ such that
$\operatorname{d}(x,\partial D)>d_{0}$ for all $x\in\omega$. This condition is
assumed without further notice throughout the paper. We emphasize that such a
requirement is also needed for computing the shape derivatives of the states
$u$ and the cost function $J$ with respect to $\omega$ (see, e.g., [26]).

The equivalence between (2) and (3) issues from the first-order necessary
condition of a minimizer of the shape functional $J(\Omega)$; that is,222Here,
the representation of the shape gradient is non-conventional. We must note
that $\nu$ denotes the inward unit normal to $\Gamma$ and
$\partial/\partial\nu$ is the inward normal derivative on $\Gamma$.

$\displaystyle{d}J(\Omega)[\boldsymbol{\mathbf{V}}]=\lim_{\varepsilon\searrow
0}\frac{J(\Omega_{\varepsilon})-J(\Omega)}{\varepsilon}=\left.\frac{{d}}{{d}\varepsilon}J(\Omega_{\varepsilon})\right|_{\varepsilon=0}=\int_{\Gamma}{G\nu\cdot\boldsymbol{\mathbf{V}}}\,{{d}}s,$

where

$G=\left(\frac{\partial
u_{\text{D}}}{\partial\nu}\right)^{2}-\left(\frac{\partial
u_{\text{N}}}{\partial\nu}\right)^{2}=\left(\frac{\partial
u_{\text{D}}}{\partial\nu}+\frac{\partial
u_{\text{N}}}{\partial\nu}\right)\left(\frac{\partial
u_{\text{D}}}{\partial\nu}-\frac{\partial
u_{\text{N}}}{\partial\nu}\right)=:G^{+}G^{-},$

has to hold for all sufficiently smooth perturbation fields
$\boldsymbol{\mathbf{V}}$. Here, $\Omega_{\varepsilon}$ stands for a
perturbation of $\Omega$ along $\boldsymbol{\mathbf{V}}$ that vanishes on the
known surface $\Sigma$. Obviously, by the calculus of variations, the above
statement implies the necessary condition $\nabla
u_{\text{N}}\cdot\nu\equiv\nabla u_{\text{D}}\cdot\nu$ on $\Gamma$. For
details on how to compute ${d}J(\Omega)[\boldsymbol{\mathbf{V}}]$, and for
more discussion on shape optimization methods in general, we refer readers to
[1] and [3].

To resolve the minimization problem (3), a typical approach is to utilize the
kernel of the shape derivative ${d}J(\Omega)[\boldsymbol{\mathbf{V}}]$, the
so-called shape gradient (see, e.g., [1, Thm. 3.6, p. 479–480]), in a
gradient-based descent algorithm. For example, for a smooth boundary
$\partial\Omega$ and sufficiently regular states $u_{\text{D}}$ and
$u_{\text{N}}$, one can take
$\boldsymbol{\mathbf{0}}\not\equiv\boldsymbol{\mathbf{V}}=-G{\nu}\in
L^{2}(\Gamma)^{d}$. This implies that, by formal expansion, and $t>0$
sufficiently small, the following inequality holds

$\displaystyle
J(\Omega_{t})=J(\Omega)+t\left.\frac{{d}}{{d}\varepsilon}J(\Omega_{\varepsilon})\right|_{\varepsilon=0}+O(t^{2})$
$\displaystyle=J(\Omega)+t\int_{\Gamma}{G\nu\cdot\boldsymbol{\mathbf{V}}}\,{{d}}s+O(t^{2})$
$\displaystyle=J(\Omega)-t\int_{\Gamma}{\left|{\boldsymbol{\mathbf{V}}}\right|^{2}}\,{{d}}s+O(t^{2})<J(\Omega).$

Hence, in practice, opting for the descent direction $\mathbf{V}=-G\nu$ is
straightforward. However, there are cases where one can deviate from this
choice. For instance, if $f>0$, we can choose $\mathbf{V}=-G^{-}\nu$, which
still provides a descent direction for $J$. Indeed, in such a case, $G^{+}<0$
on $\Gamma$. Hence, it is evident that when $\mathbf{V}=-G^{-}\nu$, the
following holds:

$\displaystyle J(\Omega_{t})$
$\displaystyle=J(\Omega)+t\int_{\Gamma}{{G^{+}}{G^{-}}\nu\cdot\boldsymbol{\mathbf{V}}}\,{{d}}s+O(t^{2})$
$\displaystyle=J(\Omega)-t\int_{\Gamma}\underbrace{{G^{+}}}_{>\
0}\left|{\boldsymbol{\mathbf{V}}}\right|^{2}\,{d}s+O(t^{2})<J(\Omega).$ (5)

A naive idea to numerically solve (2) then involves an iterative procedure
that minimizes the cost functional $J$. For simplicity, let ${\Delta{t}}>0$ be
a fixed number. Given an initial shape $\Gamma_{\circ}$, for all
$k=1,2,\ldots$, denote the free boundary at iteration $k$ as $\Gamma_{k}$.
Then, the update at iteration $k+1$ is given by
$\Gamma_{k+1}=\\{x+\Delta{t}\boldsymbol{\mathbf{V}}_{k}(x)\mid
x\in\Gamma_{k}\\}$, for $k=0,1,\ldots$. Accordingly, the main steps for
computing the $k$th free boundary $\Gamma_{k}$ using $J$ via a naive gradient-
based descent method is given in Algorithm 1.

######  Algorithm 1.

Boundary variation algorithm:

– Initialization

Fix a number $\Delta{t}>0$ and choose an initial shape $\Gamma_{\circ}$.

– Iteration

For $k=0,1,\ldots$

  1. 1.

Solve the corresponding variational formulations of (4) on $\Omega_{k}$.

  2. 2.

Define $V_{n,k}:=\nabla u_{\text{D},k}\cdot\nu_{k}-\nabla
u_{\text{N},k}\cdot\nu_{k}$ on $\Gamma^{k}$.

  3. 3.

Set $\Gamma_{k+1}=\\{x+\Delta{t}\boldsymbol{\mathbf{V}}_{k}(x)\mid
x\in\Gamma_{k}\\}$.

– Stop Test

Repeat Iteration until convergence.

The identification of the unknown boundary $\Gamma$ in (2) according to
Algorithm 1 describes a similar evolutionary equation for the Hele-Shaw-like
system (1) (see, e.g., [27, 28, 29] for the classical Hele-Shaw flow problem)
in discrete setting. Indeed, let $T>0$ be a given final time of interest (this
can be interpreted as the condition for Stop Test in Algorithm 1), $N_{T}$ be
a fixed positive integer, and set the time discretization step-size as
$\Delta{t}:=T/N_{T}$. For each time-step index $k=0,1,\cdots,N_{T}$, we let
$t_{k}=k\Delta{t}$ and denote the time-discretized domain and free boundary by
$\Omega_{k}\approx\Omega(k\Delta{t})$ and
$\Gamma_{k}\approx\Gamma(k\Delta{t})$, respectively, and the associated time-
discretized functions as $u_{\text{D},k}\approx
u_{\text{D}}(\cdot,k\Delta{t})$ and $u_{\text{N},k}\approx
u_{\text{N}}(\cdot,k\Delta{t})$. Then, given an initial free boundary
$\Gamma_{\circ}$, Algorithm 1 reduces to solving the moving boundary value
problem

$\left\\{\begin{array}[]{rcll}-\Delta u_{\text{D},k}&=&0,&\quad\text{in}\
\Omega_{k}\\\ u_{\text{D},k}&=&f,&\quad\text{on}\ \Sigma,\\\
u_{\text{D},k}&=&0,&\quad\text{on}\ \Gamma_{k}\\\ -\Delta
u_{\text{N},k}&=&0,&\quad\text{in}\ \Omega_{k}\\\ \nabla
u_{\text{N},k}\cdot\nu_{k}&=&g,&\quad\text{on}\ \Sigma,\\\
u_{\text{N},k}&=&0,&\quad\text{on}\ \Gamma_{k}\\\
V_{n,k}&=&-\left(\nabla{u_{\text{D},k}}\cdot\nu_{k}-\nabla{u_{\text{N},k}}\cdot\nu_{k}\right)&\quad\text{on}\
\Gamma_{k}\\\ \Gamma(0)&=&\Gamma_{\circ},&\end{array}\right.$ (6)

Notice that the discrete system (6) is clearly the time-discretized version of
system (1) with $f=f(x)$ and $g=g(x)$, $x\in\Sigma$.

###  1.3 Extended-regularized normal flows for shape identifications

In various studies, it has been observed that directly using
$\boldsymbol{\mathbf{V}}=-G\nu$ can lead to unwanted oscillations on the
boundary, which may cause numerical instabilities, especially in finite
element methods. To address this issue, one typically employs a Sobolev
gradient method in numerical optimization; see [30, Chap. 11]. A typical
approach includes the traction method [31] or the $H^{1}$ gradient method
[32], both widely used smoothing techniques in shape optimization problems
[33]. On the other hand, it is worth noting that treating inverse geometry
problems as shape optimization problems presents significant challenges due to
their inherent ill-posed nature [4]. These problems are inherently unstable
and are often stabilized numerically using regularization techniques [34]. For
example, in a related work [35], the authors stabilized their shape
optimization algorithm, similar to (2), by employing perimeter regularization.

Following the discussed idea, Algorithm 1 can be modified to formulate a more
stable numerical shape optimization algorithm to solve (2). Given an initial
shape $\Omega_{0}$ and denoting the domain shape at iteration $k$ as
$\Omega_{k}$, the update at iteration $k+1$ is
$\Omega_{k+1}=\\{x+\Delta{t}\boldsymbol{\mathbf{V}}_{k}(x)\mid
x\in\Omega_{k}\\}$, where ${\Delta{t}}>0$ is a sufficiently small step size
and $\boldsymbol{\mathbf{V}}_{k}$ is the descent deformation field at
iteration $k$. The main steps for computing $\Omega_{k}$ through formulation
(3) via a gradient-based descent method in Hilbert space can then be
summarized in the following optimization algorithm.

######  Algorithm 2.

Domain variation algorithm:

– Initialization

Fix a number $\Delta{t}>0$ and choose an initial shape $\Omega_{\circ}$.

– Iteration

For $k=0,1,\ldots$

  1. 1.

Solve the corresponding variational formulations of (4) on $\Omega_{k}$.

  2. 2.

Compute $\boldsymbol{\mathbf{V}}_{k}\in V(\Omega_{k})^{d}$ by solving the
variational equation

$a(\boldsymbol{\mathbf{V}}_{k},\boldsymbol{\mathbf{\varphi}})=\int_{\Gamma_{k}}{\tilde{G}_{k}\nu_{k}\cdot\boldsymbol{\mathbf{\varphi}}\,ds},\quad\forall\boldsymbol{\mathbf{\varphi}}\in
V(\Omega_{k})^{d},$

where $V(\Omega_{k}):=\\{\varphi\in H^{1}(\Omega_{k})\mid\varphi=0\ \text{on
$\Sigma$}\\}$ and $a$ is a bounded and coercive bilinear form on
$V(\Omega_{k})^{d}$.

  3. 3.

Set $\Omega_{k+1}=\\{x+\Delta{t}\boldsymbol{\mathbf{V}}_{k}(x)\mid
x\in\Omega_{k}\\}$.

– Stop Test

Repeat Iteration until convergence.

In Step 2 of Iteration, we can take $\tilde{G}$ as either exactly $G$ or
${G^{-}}$, depending on the sign of $f$ and $g$. Meanwhile, for the Stop Test,
we simply terminate the algorithm after a finite number of iterations.

As emphasized in previous subsections, our main intent in this exposition is
to prove the existence and uniqueness of classical solution to (1). Before
delving with this problem, however, we first exhibit a numerical example in
two spatial dimensions illustrating the feasibility of the proposed
optimization algorithm with the choice of preconditioned descent vector field
$\boldsymbol{\mathbf{V}}$ obtained through the vector ${G^{-}}{\nu}$. Here, we
define
$a(\boldsymbol{\mathbf{\varphi}},\boldsymbol{\mathbf{\psi}}):=(\nabla{\boldsymbol{\mathbf{\varphi}}},\nabla{\boldsymbol{\mathbf{\psi}}})_{\Omega}$,
for $\boldsymbol{\mathbf{\varphi}},\boldsymbol{\mathbf{\psi}}\in V(\Omega)$.
Hereon, we shall refer to our algorithm as quasi-stationary Stefan type-scheme
or QSSTS as it also provides a kind of a comoving mesh method [36, 37] for
(1).

######  Example 1.4.

Using Algorithm (2)333Here, we used a non-uniform time step size to clearly
highlight the potential of taking $\tilde{G}=G^{-}$. In fact, we determine the
step size using a backtracking line search procedure, calculated as
${\Delta{t}_{k}}=cJ(\Omega_{k})/|\boldsymbol{\mathbf{V}}_{k}|^{2}_{\mathbf{H}^{1}(\Omega_{k})}$
at each time step, where $c>0$ is a scaling factor. This factor is adjusted to
prevent reversed triangles within the mesh after updating (see [38, p. 281]).,
we reconstruct a kite-shaped cavity under the following setup. The specimen
under investigation is a circular disk with a unit radius, centered at the
origin. For the input data, we set $f\equiv 1$ on $\Sigma$, while synthetic
data is used for extra boundary measurements on the outer part of $\Sigma$. To
avoid “inverse crimes” (see [39, p. 179]), the synthetic data is generated
using a different numerical scheme, employing a larger number of
discretization points and ${P2}$ finite element basis functions in the
FreeFem++ code [40]. In the inversion procedure, all variational problems are
solved using ${P1}$ finite elements with coarse meshes.

We contrast the traditional shape optimization approach (utilizing
$\tilde{G}=G$) with the proposed method (employing $\tilde{G}=G^{-}$). The
geometry of the initial free boundary $\Gamma_{\circ}$ is a circle centered at
the origin with a radius of $0.9$. The experimental findings, outlined in
Figure 2, indicate a clear advantage of the proposed method (QSSTS) over the
conventional shape optimization method (SO). The QSSTS method not only
achieves faster identification of the shape and location of the unknown
inclusion but also does so with greater accuracy. Although not illustrated
here, we note that similar advantages of the proposed method over traditional
shape optimization approaches were observed in other instances as well.

Figure 2: Identification of a kite-shape cavity (left plot), histories of
Hausdorff distances (middle plot), and residuums (right plot).

The rest of the paper is organized as follows. In Section 2, we establish
several auxiliary results related to a diffeomorphic mapping between a
reference domain and its perturbation. In Section 3, we present the main
result of the paper (see Theorem 3.1), which addresses the local-in-time
solvability of system (1) by employing techniques developed by Solonnikov in
[15]. The remaining sections sequentially outline the key steps in proving the
main result. Section 4 details the transformation of the original problem onto
a fixed domain. In Section 5, we establish the regularity of solutions on the
fixed domain. Section 6 addresses the nonlinear problem with respect to a
parameter. This is followed by an analysis of the linearized problem in
Section 7. Finally, the proof of the main result, as stated in Theorem 3.1, is
completed in Section 8.

This work also includes several appendices that provide preliminary analyses
of the well-posedness of (1) in the case of axis-symmetric domains (see
Appendices A and B), demonstrations of transformations and computations of key
equations and identities (Appendices C and D) utilized in the paper, and
detailed proofs of auxiliary lemmas used in the proof of the main result
(Appendix E).

##  2 Auxiliary Results

To state clearly our main result, we start by fixing some notations and prove
some preliminary results in this section. Throughout the paper, for
$k=0,1,\ldots$, ${{C}}^{k}(\overline{\Omega})$ denotes the usual space of all
functions having continuous and bounded derivatives in $\overline{\Omega}$ up
to $k$th order while, with $\alpha\in(0,1)$, ${{C}}^{k+\alpha}(\Omega)$ is the
Banach space of all functions $u\in{{C}}^{k}(\Omega)$ for which the norm444See
also the norms issued at the beginning of Section 3 (cf. notations in [41,
42].

$\left|{u}\right|^{(k+\alpha)}_{\Omega}=\left|{u}\right|^{(k)}_{\Omega}+[u]^{(k+\alpha)}_{\Omega},$

is finite, where

$\left|{u}\right|^{(k)}_{\Omega}=\sum_{j=0}^{k}[u]^{(j)}_{\Omega},\qquad[u]^{(j)}_{\Omega}=\max_{\left|{\beta}\right|=j}\left|{D^{\beta}{u}}\right|^{(0)}_{\Omega},\qquad\left|{u}\right|^{(0)}_{\Omega}=[u]^{(0)}_{\Omega}=\sup_{\Omega}\left|{u}\right|,$

and

$[u]^{(k+\alpha)}_{\Omega}=\max_{\left|{\beta}\right|=k}[D^{\beta}{u}]^{(\alpha)}_{\Omega},\qquad[u]^{(\alpha)}_{\Omega}\sup_{\begin{subarray}{c}x,x^{\prime}\in\Omega\\\
x\neq
x^{\prime}\end{subarray}}\frac{\left|{u(x)-u(y)}\right|}{\left|{x-y}\right|^{\alpha}}.$

Given a bounded simply connected domain $D\subset\mathbb{R}^{d}$ with boundary
$\partial{D}=\Sigma$, we define the set of admissible geometries (the boundary
of inclusions) as follows:

$\mathscr{A}^{2+\alpha}:=\left\\{\Gamma=\partial{\omega}\mid\overline{\Omega}\subset
D,\ \text{$\omega$ is a simply connected bounded domain and
$\partial{\omega}\in{{C}}^{2+\alpha}$}\right\\}.$ (7)

For $\Gamma\in\mathscr{A}^{2+\alpha}$, we introduce the notation
$\Omega(\Gamma)$ to denote an annular domain in $\mathbb{R}^{d}$ with boundary
$\partial\Omega(\Gamma)=\Gamma\cup\Sigma$. Accordingly, given
$f\in{{C}}^{2+\alpha}(\Sigma)$, by $u_{\text{D}}(\Gamma)$,
$\Gamma\in\mathscr{A}^{2+\alpha}$, we denote by $u_{\text{D}}(\Gamma)$ a
unique solution of

$\left\\{\begin{array}[]{rcll}u_{\text{D}}&\in&{{C}}^{2+\alpha}(\overline{\Omega(\Gamma)})\\\
-\Delta u_{\text{D}}&=&0,&\quad\text{in}\ \Omega(\Gamma)\\\
u_{\text{D}}&=&f,&\quad\text{on}\ \Sigma,\\\ u_{\text{D}}&=&0,&\quad\text{on}\
\Gamma.\end{array}\right.$ (8)

Similarly, with $g\in{{C}}^{1+\alpha}(\Sigma)$ given, $u_{\text{N}}(\Gamma)$,
$\Gamma\in\mathscr{A}^{2+\alpha}$, denotes the unique solution of

$\left\\{\begin{array}[]{rcll}u_{\text{N}}&\in&{{C}}^{2+\alpha}(\overline{\Omega(\Gamma)})\\\
-\Delta u_{\text{N}}&=&0,&\quad\text{in}\ \Omega(\Gamma)\\\
\dfrac{\partial{u_{\text{N}}}}{\partial\nu}&=&g,&\quad\text{on}\ \Sigma,\\\
u_{\text{N}}&=&0,&\quad\text{on}\ \Gamma.\end{array}\right.$ (9)

Hereinafter, unless otherwise stated, we have
$\Gamma\in\mathscr{A}^{2+\alpha}$.

For technical purposes and for economy of space, we need to introduce a notion
of quasi-normal vector on $\Gamma$.

######  Definition 2.1.

We say that a vector field $\mathbf{N}$ is quasi-normal on
$\Gamma\in\mathscr{A}^{2+\alpha}$, inheriting the regularity of $\Gamma$, if

$\boldsymbol{\mathbf{N}}\in{{C}}^{2+\alpha}(\Gamma;\mathbb{R}^{d})$ and it is
such that $\left|{\boldsymbol{\mathbf{N}}(\xi)}\right|=1$ and
$\boldsymbol{\mathbf{N}}(\xi)\cdot\nu(\xi;\Gamma)>0$, for all $\xi\in\Gamma$.
(10)

We let ${I}_{\varepsilon_{0}}:=[-\varepsilon_{0},\varepsilon_{0}]$ and fix a
constant $\varepsilon_{0}=\varepsilon_{0}(\Gamma,\boldsymbol{\mathbf{N}})>0$
such that the map

$X:\Gamma\times{I}_{\varepsilon_{0}}\ \longrightarrow\
\Gamma^{\varepsilon_{0}}\subset\mathbb{R}^{d},\qquad X(\xi,\rho)\ \longmapsto\
\xi+\rho\boldsymbol{\mathbf{N}}(\xi)\subset{D},$

is a ${{C}}^{2+\alpha}$-diffeomorphism, where

$\Gamma^{{\varepsilon}}:=\\{X(\xi,r):=\xi+r\boldsymbol{\mathbf{N}}(\xi)\mid(\xi,r)\in\Gamma\times{I}_{\varepsilon}\\},$

for ${\varepsilon}>0$. In fact, we have the following proposition.

######  Proposition 2.2.

There exists a constant $\varepsilon_{0}>0$ such that

$X\in\operatorname{Diffeo}^{2+\alpha}(\Gamma\times\overline{I}_{\varepsilon_{0}};{\Gamma^{\varepsilon_{0}}}),\qquad{\Gamma^{\varepsilon_{0}}}:=X(\Gamma\times\overline{I}_{\varepsilon_{0}}).$

###### Proof.

Let $\varrho\in B_{1}:=B_{1}{(\mathbb{R}^{d-1})}$ (i.e., $B_{1}$ denotes the
unit ball in $\mathbb{R}^{d-1}$). For all $\xi_{0}\in\Gamma$, there exist an
open set $B(\xi_{0})$ of $\mathbb{R}^{d}$ containing $\xi_{0}$ and a mapping
$\varphi\in\operatorname{Diffeo}^{2+\alpha}(\overline{B_{1}},\overline{\Gamma\cap
B(\xi_{0})})$ such that

$\Upsilon:\overline{B}_{1}\times{I}_{\varepsilon}\to\mathbb{R}^{d},\qquad\Upsilon(\varrho,r):=X(\varphi(\varrho),r)=\varphi(\varrho)+r\boldsymbol{\mathbf{N}}(\varphi(\varrho)).$

We shall prove first that the determinant of the Jacobian of the map
$\Upsilon$ is non-zero for all
$(\varrho,r)\in\overline{B}_{1}\times{I}_{\varepsilon}$.

We let $\widetilde{\boldsymbol{\mathbf{N}}}$ be a smooth extension of
$\boldsymbol{\mathbf{N}}$ in a tubular neighborhood of $\Gamma$ and define
$\textsf{N}(\varrho,r)=\textsf{I}+r\nabla_{\Gamma}^{\top}\widetilde{\boldsymbol{\mathbf{N}}}(\varphi(\varrho))$.
By straightforward computations, we get

$\displaystyle\mathbb{R}^{d\times
d}\ni\nabla_{(\varrho,r)}^{\top}\Upsilon{(\varrho,r)}$
$\displaystyle=\begin{pmatrix}\nabla_{\varrho}^{\top}\Upsilon{(\varrho,r)}&\partial_{r}\Upsilon{(\varrho,r)}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\nabla_{\varrho}^{\top}\varphi(\varrho)+r\left(\nabla_{x}^{\top}\widetilde{\boldsymbol{\mathbf{N}}}(x)\big{|}_{x=\varphi(\varrho)}\right)\nabla_{\varrho}^{\top}\varphi(\varrho)&\boldsymbol{\mathbf{N}}(\varphi(\varrho))\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\nabla_{\varrho}^{\top}\varphi(\varrho)+r\nabla_{\Gamma}^{\top}\widetilde{\boldsymbol{\mathbf{N}}}(\varphi(\varrho))\nabla_{\varrho}^{\top}\varphi(\varrho)&\boldsymbol{\mathbf{N}}(\varphi(\varrho))\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\left(\textsf{I}+r\nabla_{\Gamma}^{\top}\widetilde{\boldsymbol{\mathbf{N}}}(\varphi(\varrho))\right)\nabla_{\varrho}^{\top}\varphi(\varrho)&\boldsymbol{\mathbf{N}}(\varphi(\varrho))\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\nabla_{\varrho}^{\top}\varphi(\varrho)&\boldsymbol{\mathbf{N}}(\varphi(\varrho))\end{pmatrix}$
$\displaystyle=\textsf{N}(\varrho,r)\begin{pmatrix}\textsf{N}^{-1}(\varrho,r)\nabla_{\varrho}^{\top}\varphi(\varrho)&\textsf{N}^{-1}(\varrho,r)\boldsymbol{\mathbf{N}}(\varphi(\varrho))\end{pmatrix}.$

For $\xi:=\varphi(\varrho)$, we claim that

$\left[\textsf{N}^{-1}(\varrho,r)\boldsymbol{\mathbf{N}}(\xi)\right]\cdot\nu(\xi)>0.$

Indeed, since $\boldsymbol{\mathbf{N}}$ is a quasi-normal vector on $\Gamma$,
there exists a constant $c>0$ such that
$\boldsymbol{\mathbf{N}}(\xi)\cdot\nu(\xi)\geqslant{c}>0$, for all
$\xi\in\Gamma$. Noting that

$\textsf{N}^{-1}-\textsf{I}=\textsf{N}^{-1}(\textsf{I}-\textsf{N})=r\textsf{N}^{-1}\nabla_{\Gamma}^{\top}\boldsymbol{\mathbf{N}},$

we can choose $\varepsilon_{0}\leqslant\varepsilon_{1}$ where

$\varepsilon_{1}={c}\min\left\\{1,\frac{1}{2\left\|{\textsf{N}^{-1}\nabla_{\Gamma}^{\top}\boldsymbol{\mathbf{N}}}\right\|}\right\\}$

so that

$\displaystyle\left[\textsf{N}^{-1}(\varrho,r)\boldsymbol{\mathbf{N}}(\xi)\right]\cdot\nu(\xi)$
$\displaystyle=\left[\textsf{N}^{-1}(\textsf{I}-\textsf{N})(\varrho,r)\boldsymbol{\mathbf{N}}(\xi)\right]\cdot\nu(\xi)+\boldsymbol{\mathbf{N}}(\xi)\cdot\nu(\xi)$
$\displaystyle\geqslant{c}-\left\|{\textsf{N}^{-1}-\textsf{I}}\right\|\geqslant\frac{c}{2},$

where $\left\|{\cdot}\right\|$ is the usual entry-wise matrix norm.

We next prove that the map $X$ is injective. For this purpose, we show that

if ${\xi},{\tilde{\xi}}\in\Gamma$ and
$X({\xi},{\rho})=X({\tilde{\xi}},{\tilde{\rho}})$, then
$({\xi},{\rho})=({\tilde{\xi}},{\tilde{\rho}})$.

Again, we note that for each $\xi\in\Gamma$, there exist an open set
$B({\xi})$ of $\mathbb{R}^{d}$ containing $\xi$. Meanwhile, since $\Gamma$ is
compact, there exists a finite collection of open covers
$\\{O_{j}\\}_{j=1}^{M}$, $M\in\mathbb{N}$, such that
$\Gamma=\bigcup_{j=1}^{M}O_{j}$. For each $j=1,\ldots,M$, we let
$\xi_{j}\in\Gamma$ and define, for each open set $B(\xi_{j})\ni\xi_{j}$, the
set $O_{j}=\Gamma\cap B(\xi_{j})\ni\xi_{j}$. Clearly, for
$\xi_{i}\neq\xi_{j}$, there exists $r_{0}>0$ such that

$\inf_{\begin{subarray}{c}\\{\xi_{i},\xi_{j}\\}\not\subset O_{j}\\\
\xi_{i},\xi_{j}\in\Gamma\end{subarray}}\left|{\xi_{i}-\xi_{j}}\right|=r_{0}>0.$
(11)

By the previous argument, it suffices to prove that

if ${\xi},{\tilde{\xi}}\in\Gamma$ and
$X({\xi},{\rho})=X({\tilde{\xi}},{\tilde{\rho}})$, then ${\xi}={\tilde{\xi}}$.

Let us assume that ${\xi},{\tilde{\xi}}\in\Gamma\times{I}_{\varepsilon}$,
$X({\xi},{\rho})=X({\tilde{\xi}},{\tilde{\rho}})$, and
${\xi}\neq{\tilde{\xi}}$. The latter condition implies that there exist
$i,j\in\\{1,\ldots,M\\}$, $i\neq j$, such that
${\xi}\in{O}_{i}\setminus{O}_{j}$ and $\tilde{\xi}\in O_{j}\setminus O_{i}$.
Now, we let

${r}_{0}:=\inf\\{\left|{{\xi}-{\tilde{\xi}}}\right|\mid{\xi},{\tilde{\xi}}\in\Gamma\
\text{and}\ \\{\xi_{i},\xi_{j}\\}\not\subset O_{j},\ \text{for all
$j=1,\ldots,M$}\\}.$

Meanwhile, we have

$0\neq{\xi}-{\tilde{\xi}}={\rho}\boldsymbol{\mathbf{N}}({\xi})-{\tilde{\rho}}\boldsymbol{\mathbf{N}}({\tilde{\xi}})=({\rho}-{\tilde{\rho}})\boldsymbol{\mathbf{N}}({\xi})+{\tilde{\rho}}({\boldsymbol{\mathbf{N}}}({\xi})-{\boldsymbol{\mathbf{N}}}({\tilde{\xi}})).$

Hence,

$\left|{{\xi}-{\tilde{\xi}}}\right|\leqslant\left|{{\rho}-{\tilde{\rho}}}\right|+L_{0}\left|{{\tilde{\rho}}}\right|\left|{{\xi}-{\tilde{\xi}}}\right|,\qquad\text{where}\
L_{0}:=\sup_{\begin{subarray}{c}x,y\in\Gamma\\\ x\neq
y\end{subarray}}\frac{\left|{\boldsymbol{\mathbf{N}}(x)-\boldsymbol{\mathbf{N}}(y)}\right|}{\left|{x-y}\right|},$

Since $\left|{{\tilde{\rho}}}\right|\leqslant\varepsilon_{0}$, we get

$(1-\varepsilon_{0}L_{0})\left|{{\xi}-{\tilde{\xi}}}\right|\leqslant\left|{{\rho}-{\tilde{\rho}}}\right|\leqslant
2\varepsilon_{0}.$

From (11), we deduce that

$r_{0}=\inf_{\begin{subarray}{c}\\{\xi_{i},\xi_{j}\\}\not\subset O_{j}\\\
\xi_{i},\xi_{j}\in\Gamma\end{subarray}}\left|{\xi_{i}-\xi_{j}}\right|\leqslant\left|{{\xi}-{\tilde{\xi}}}\right|\leqslant\frac{2\varepsilon_{0}}{1-\varepsilon_{0}L_{0}}.$

Taking $\varepsilon_{0}<\dfrac{r_{0}}{2+r_{0}L_{0}}$ we arrive at a
contradiction unless otherwise $r_{0}=0$. Thus, ${\xi}={\tilde{\xi}}$.

For the sake of being specific, we choose

$\varepsilon_{0}=c\min\left\\{\dfrac{r_{0}}{3+r_{0}L_{0}},\varepsilon_{1}\right\\},\qquad\varepsilon_{1}={c}\min\left\\{1,\frac{1}{2\left\|{\textsf{N}^{-1}\nabla_{\Gamma}^{\top}\boldsymbol{\mathbf{N}}}\right\|}\right\\}$

to conclude the proof. ∎

######  Remark 2.3 (Extensions of the normal vector).

In the previous proof, it suffices to assume that
$\widetilde{\boldsymbol{\mathbf{N}}}$ is ${C}^{1+\alpha}$. This is guaranteed
because $\Gamma\in\mathscr{A}^{2+\alpha}$ implies that $\nu\in C^{1+\alpha}$,
and therefore we can create a ${C}^{1+\alpha}$ extension of $\nu$ in
$D\supset\Gamma$. More specifically, there exists an
$\varepsilon_{\star}$-neighborhood of $\Gamma$
($\varepsilon_{\star}\geqslant\varepsilon_{0}$), say
$\mathscr{N}^{\varepsilon_{\star}}(\Gamma)\supset
X(\Gamma\times{I}_{\varepsilon_{0}})$, such that
$\widetilde{\boldsymbol{\mathbf{N}}}\in{C}^{1+\alpha}(\mathscr{N}^{\varepsilon_{\star}}(\Gamma);\mathbb{R}^{d})$.
Additionally, we see that the distance function
$\operatorname{d}(\cdot,\Gamma)\in{C}^{2+\alpha}(\mathscr{N}^{\varepsilon_{\star}}(\Gamma))$
where $\operatorname{d}(x,\Gamma):=\pm\operatorname{dist}(x,\Gamma)$.
Consequently, for any $x\in\mathscr{N}^{\varepsilon_{\star}}(\Gamma)$, we can
define
$\widetilde{\boldsymbol{\mathbf{N}}}(x):=\boldsymbol{\mathbf{N}}(\xi(x))$
where $\xi(x)=x-\operatorname{d}(x,\Gamma)\nabla\operatorname{d}(x,\Gamma)$.
Furthermore, for later use, we emphasize that we can create another extension
of $\boldsymbol{\mathbf{N}}$ which is ${C}^{2+\alpha}$. That is, for
$(\xi,\rho)\in\Gamma\times{I}_{\varepsilon_{0}}$, we can define
${\mathchoice{\widetilde{\hbox{\set@color$\displaystyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\textstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\scriptstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\scriptscriptstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}}(X(\Gamma\times{I}_{\varepsilon_{0}})):=\boldsymbol{\mathbf{N}}(\xi)$
such that
${\mathchoice{\widetilde{\hbox{\set@color$\displaystyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\textstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\scriptstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}{\widetilde{\hbox{\set@color$\scriptscriptstyle\widetilde{\boldsymbol{\mathbf{N}}}$}}}}\big{|}_{\Gamma}=\boldsymbol{\mathbf{N}}$.

For fixed real numbers $a$, $b$ where $b>a$, the scalar valued $\rho$ is such
that it belongs to the Banach space

${{R}}_{[a,b]}(\Gamma,\boldsymbol{\mathbf{N}}):=\left\\{\rho:\Gamma\times[a,b]\to{I}_{\varepsilon_{0}(\Gamma,\boldsymbol{\mathbf{N}})}\mid\rho\in{{C}}^{0}([a,b];{{C}}^{2+\alpha}(\Gamma))\cap{{C}}^{1}([a,b];{{C}}^{1+\alpha}(\Gamma))\right\\}.$
(12)

We also introduce the set

$\displaystyle{{R}}_{0}(\Gamma,\boldsymbol{\mathbf{N}})$
$\displaystyle:=\left\\{\rho\in{{C}}^{2+\alpha}(\Gamma)\mid\left|{\rho(\xi)}\right|\leqslant\varepsilon_{0}(\Gamma,\boldsymbol{\mathbf{N}}),\
\text{$\forall\xi\in\Gamma$}\right\\}.$

For $\rho\in{{R}}_{0}(\Gamma,\boldsymbol{\mathbf{N}})$, it can be shown that

${\mathscr{S}}(\rho):=\\{X(\xi,\rho)\mid\xi\in\Gamma\\}$ is a
${{C}}^{2+\alpha}$ boundary.

In fact, ${\mathscr{S}}(\rho)\in\mathscr{A}^{2+\alpha}$ as we claim in the
following proposition.

######  Proposition 2.4.

There exists
$\varepsilon_{1}:=\varepsilon_{1}(\Gamma,\boldsymbol{\mathbf{N}})\in(0,\varepsilon_{0}(\Gamma,\boldsymbol{\mathbf{N}})]$
such that ${\mathscr{S}}(\rho)\in\mathscr{A}^{2+\alpha}$ holds for
$\rho\in{{R}}_{1}(\Gamma,\boldsymbol{\mathbf{N}})$, where

${{R}}_{1}(\Gamma,\boldsymbol{\mathbf{N}}):=\left\\{\rho\in{{R}}_{0}(\Gamma,\boldsymbol{\mathbf{N}})\mid\text{$\left|{\rho(\xi)}\right|\leqslant\varepsilon_{1}$,
$\left|{\nabla_{\Gamma}\rho(\xi)}\right|\leqslant\varepsilon_{1},\forall\xi\in\Gamma$}\right\\}.$

To prove the proposition, we need the next lemma.

######  Lemma 2.5.

Let $k\in\mathbb{N}$, $\alpha\in[0,1)$, and $\Omega\subset\mathbb{R}^{d}$ be
an open bounded set with ${{C}}^{k+\alpha}$ boundary. Let
$\phi\in{{C}}_{0}^{k+\alpha}(\Omega)$ and consider $\varphi(x)=x+\phi(x)$,
$x\in\Omega$. Assume that
$\max_{x\in\overline{\Omega}}\left\|{\nabla^{\top}\phi(x)}\right\|<1$. Then,
$\operatorname{{det}}(\nabla^{\top}\varphi)>0$ and
$\varphi\in\operatorname{Diffeo}^{k+\alpha}(\Omega,\Omega)$; i.e., the map
$\varphi:\Omega\to\Omega$ is a ${{C}}^{k+\alpha}$-diffeomorphism.

See Appendix Appendix E.1 for the proof.

###### Proof of Proposition 2.4.

Let $\varepsilon_{0}:=\varepsilon_{0}(\Gamma,\boldsymbol{\mathbf{N}})$. The
map $X:\Gamma\times{I}_{\varepsilon_{0}}\to\Gamma^{\varepsilon_{0}}$ is a
${{C}}^{2+\alpha}$-diffeomorphism. We choose and fix
$\varepsilon>\varepsilon_{0}$, with
$\left|{\varepsilon-\varepsilon_{0}}\right|$ sufficiently small, such that the
map
$X:\Gamma\times{I}_{\varepsilon}\to\Gamma^{\varepsilon}\subset{D}\subset\mathbb{R}^{d}$
is ${{C}}^{2+\alpha}$-diffeomorphism. Consequently,
$X^{-1}:\Gamma^{\varepsilon}\to\Gamma\times{I}_{\varepsilon}$ is also a
${{C}}^{2+\alpha}$-diffeomorphism; i.e., $X^{-1}=(\hat{\xi}(x),\hat{r}(x))$
where $\hat{\xi}\in{{C}}^{2+\alpha}(\Gamma^{\varepsilon},\Gamma)$ and
$\hat{r}\in{{C}}^{2+\alpha}(\Gamma^{\varepsilon},{I}_{\varepsilon})$. We
consider $\eta\in{{C}}_{0}^{\infty}(\mathbb{R})$ ($0\leqslant\eta(s)\leqslant
1$) such that $\eta(s)=1$ for $\left|{s}\right|\leqslant\varepsilon_{0}$ and
$\operatorname{supp}\eta\subset(-\varepsilon,\varepsilon)$.

Let us define $\Phi$ as follows:

$\Phi(x)=\begin{cases}x&\text{if
}x\in\overline{D}\setminus\Gamma^{\varepsilon},\\\
x+\eta(\hat{r}(x))\rho(\hat{\xi}(x))\boldsymbol{\mathbf{N}}(\hat{\xi}(x))&\text{if
}x\in\Gamma^{\varepsilon}.\end{cases}$

Clearly, $\Phi\in{{C}}^{2+\alpha}(\overline{D};\mathbb{R}^{d})$. We claim that
$\Phi\in\operatorname{Diffeo}^{2+\alpha}(\Omega,\Omega)$. Indeed, this follows
from Lemma 2.5. We only need to verify that the norm of the Jacobian of its
perturbation – given that $\left|{\rho(\xi)}\right|\leqslant\varepsilon_{1}$,
$\left|{\nabla_{\Gamma}\rho(\xi)}\right|\leqslant\varepsilon_{1}$, for all
$\xi\in\Gamma$ – has magnitude less than one. To facilitate the proof, let us
write
$\boldsymbol{\mathbf{G}}(x):=\eta(\hat{r}(x))\boldsymbol{\mathbf{N}}(\hat{\xi}(x))$
and define $\Phi(x)=x+\phi(x)$ where

$\phi(x)=\begin{cases}\boldsymbol{\mathbf{0}}&\text{if
}x\in\overline{D}\setminus\Gamma^{\varepsilon},\\\
\boldsymbol{\mathbf{G}}(x)\rho(\hat{\xi}(x))&\text{if
}x\in\Gamma^{\varepsilon}.\end{cases}$

From the definition of $\phi$, we have

$\nabla^{\top}\phi(x)=\begin{cases}\boldsymbol{\mathbf{0}}&\text{if
}x\in\overline{D}\setminus\Gamma^{\varepsilon},\\\
\nabla^{\top}\boldsymbol{\mathbf{G}}(x)\rho(\hat{\xi}(x))+\boldsymbol{\mathbf{G}}(x)\nabla^{\top}_{\Gamma}\rho(\hat{\xi}(x))\nabla^{\top}_{x}\hat{\xi}(x)&\text{if
}x\in\Gamma^{\varepsilon}.\end{cases}$

Computing its norm, while noting that
$\left|{\rho(\xi)}\right|\leqslant\varepsilon_{1}$,
$\left|{\nabla_{\Gamma}\rho(\xi)}\right|\leqslant\varepsilon_{1}$, for all
$\xi\in\Gamma$, we get

$\displaystyle\left\|{\nabla^{\top}\phi(x)}\right\|$
$\displaystyle\leqslant\left\|{\rho}\right\|_{\infty}\left\|{\nabla^{\top}\boldsymbol{\mathbf{G}}}\right\|+\left|{\boldsymbol{\mathbf{G}}(x)}\right|\left|{\nabla^{\top}_{\Gamma}\rho(\hat{\xi}(x))}\right|\left\|{\nabla^{\top}_{x}\hat{\xi}(x)}\right\|$
$\displaystyle\leqslant\varepsilon_{1}\left[\max_{x\in\Gamma^{\varepsilon}}\left\|{\nabla^{\top}\boldsymbol{\mathbf{G}}(x)}\right\|+\max_{x\in\Gamma^{\varepsilon}}\left(\left|{\boldsymbol{\mathbf{G}}(x)}\right|\left\|{\nabla^{\top}_{x}\hat{\xi}(x)}\right\|\right)\right]<1.$

This proves the proposition. ∎

For $\rho\in{{R}}_{[a,b]}(\Gamma,\boldsymbol{\mathbf{N}})$, we define the
moving boundary

$\mathscr{M}(\rho,[a,b]):=\bigcup_{t\in[a,b]}{\mathscr{S}}(\rho(t))\times\\{t\\},$
(13)

with normal velocity
$V_{n}(t)=V_{n}(\cdot,t)\in{{C}}^{0}({\mathscr{S}}(\rho(t)))$ where

$V_{n}(x,t):=\rho_{t}(\xi,t)\boldsymbol{\mathbf{N}}(\xi)\cdot\nu(x;{\mathscr{S}}(\rho(t))),\qquad
x=\xi+\rho(\xi,t)\boldsymbol{\mathbf{N}}(\xi)\in{\mathscr{S}}(\rho(t)),\quad\xi\in\Gamma.$

We define the set of moving boundaries

${{M}}_{[a,b]}(\Gamma,\boldsymbol{\mathbf{N}}):=\left\\{\mathscr{M}(\rho,[a,b])\subset\mathbb{R}^{d}\times\mathbb{R}\mid\text{$\exists\rho\in{{R}}_{[a,b]}(\Gamma,\boldsymbol{\mathbf{N}})$
such that \eqref{eq:M_rho_ab} is satisfied}\right\\}.$

Now, we consider the quasi-stationary moving boundary problem stated as
follows:

######  Problem 2.6.

For given $\Gamma_{\circ}\in\mathscr{A}^{2+\alpha}$,
$f\in{{C}}^{2+\alpha}(\Sigma)$, and $g\in{{C}}^{1+\alpha}(\Sigma)$, find $T>0$
and $\mathscr{M}=\bigcup_{0\leqslant t\leqslant T}\Gamma(t)\times\\{t\\}$ such
that

$\left\\{\begin{aligned}
V_{n}(t)&=-\dfrac{\partial{}}{\partial\nu}\left[u_{\text{D}}(\Gamma(t))-u_{\text{N}}(\Gamma(t))\right]\quad\text{on
$\Gamma(t)$, \quad$(0\leqslant t\leqslant T)$},\\\
\Gamma(0)&=\Gamma_{\circ},\end{aligned}\right.$ (14)

where $u_{\text{D}}(\Gamma(t))$ and $u_{\text{N}}(\Gamma(t))$ are defined by
(8) and (9), respectively.

We formally define the solution of Problem 2.6 as follows.

######  Definition 2.7.

We say $\mathscr{M}=\bigcup_{0\leqslant t\leqslant
T}\Gamma(t)\times\\{t\\}\subset\mathbb{R}^{d}\times\mathbb{R}$ a solution of
Problem 2.6, if for $\Gamma(0)=\Gamma_{\circ}$, there exists a collection of
closed intervals $\\{I_{k}\\}_{k=1}^{n}$ such that
$\bigcup_{k=1}^{n}{I_{k}}=[0,T]$, and for each $k$, there exists $t_{k}\in
I_{k}$, $\Gamma_{k}\in\mathscr{A}^{2+\alpha}$, and quasi-normal
$\boldsymbol{\mathbf{N}}_{k}$ on $\Gamma_{k}$ such that

$\mathscr{M}\big{|}_{I_{k}}\in{{M}}_{I_{k}}(\Gamma_{k},\boldsymbol{\mathbf{N}}_{k})\quad\text{where}\quad\mathscr{M}\big{|}_{I_{k}}=\bigcup_{t\in
I_{k}}\Gamma(t)\times\\{t\\},$

is a solution of

$\text{$V_{n}(t)=-\dfrac{\partial{}}{\partial\nu}\left[u_{\text{D}}(\Gamma(t))-u_{\text{N}}(\Gamma(t))\right]$
on $\Gamma(t)$, \quad for $t\in I_{k}$, for each $k=1,\ldots,n$},$

where $u_{\text{D}}(\Gamma(t))$ and $u_{\text{N}}(\Gamma(t))$ are defined by
(8) and (9), respectively.

We note that the definition of $V_{n}$ does not depends on the choice of
$\Gamma_{k}$ and $\boldsymbol{\mathbf{N}}_{k}$.

Before we proceed further, we ask the following question in view of Definition
2.7: suppose $\mathscr{M}_{[0,T]}$ is a solution to (1), then is it true that
$\mathscr{M}_{[{t}_{\ast},T]}$ is also a solution to (1) for
${t}_{\ast}\in[0,T]$? The next lemma answers this question affirmatively.

######  Lemma 2.8.

Let $\Gamma\in\mathscr{A}^{2+\alpha}$, $\boldsymbol{\mathbf{N}}$ be a quasi-
normal vector on $\Gamma$,

$\mathscr{M}=\bigcup_{a\leqslant t\leqslant
b}\Gamma(t)\times\\{t\\}\in{{M}}_{[a,b]}(\Gamma,\boldsymbol{\mathbf{N}}),$

${t}_{\ast}\in[a,b]$, and ${\boldsymbol{\mathbf{N}}}_{\ast}$ be quasi-normal
on $\Gamma({t}_{\ast})$. Then, there exists a ${\tau}>0$ such that, setting
$I_{\ast}:=[a,b]\cap[{t}_{\ast}-{\tau},{t}_{\ast}+{\tau}]$, we have

$\mathscr{M}\big{|}_{I_{\ast}}\in{{M}}_{I_{\ast}}(\Gamma({t}_{\ast}),{\boldsymbol{\mathbf{N}}}_{\ast}).$
(15)

###### Proof.

We start by noting that there exist $\varepsilon_{\ast}>0$ such that

${{X}_{\ast}}\in\operatorname{Diffeo}^{2+\alpha}\left(\Gamma({t}_{\ast})\times{I}_{\ast},{X}_{\ast}(\Gamma({t}_{\ast})\times{I}_{\ast})\right),\qquad{{X}_{\ast}}(x,r):=x+r{\boldsymbol{\mathbf{N}}}_{\ast}(x),$

where

$(x,r)\in\Gamma({t}_{\ast})\times{I}_{\ast},\qquad
I_{\ast}:={I}_{\varepsilon_{\ast}}\cap[a,b],\qquad{I}_{\varepsilon_{\ast}}=[-\varepsilon_{\ast},\varepsilon_{\ast}].$

We claim the following:

There exists a ${\tau}>0$ such that
$\Gamma(t)\subset{X}_{\ast}(\Gamma({t}_{\ast})\times{I}_{\ast\ast})$ for
$t\in{I}_{\ast\ast}:={I}_{\ast}\cap[{t}_{\ast}-{\tau},{t}_{\ast}+{\tau}]$.

Now, for $(y,t)\in\Gamma\times{I}_{\ast\ast}$ we define

$({x}_{\ast}(y,t),{\rho}_{\ast}(y,t)):={X}_{\ast}^{-1}(X(y,\rho(y,t)))$

and for $(x,t)\in\Gamma({t}_{\ast})\times{I}_{\ast\ast}$, we set

$\widetilde{\rho}(x,t):={\rho}_{\ast}({x}_{\ast}(\cdot,t)^{-1}(x),t),$

where ${x}_{\ast}(\cdot,t):\Gamma\to\Gamma({t}_{\ast})$. Hence, by
construction

$X(y,\rho(y,t))={X}_{\ast}(x,\tilde{\rho}(x,t))\Big{|}_{x={x}_{\ast}(y,t)}.$

Now, to validate (15), we need to show that

$\tilde{\rho}\in{{C}}^{0}({I}_{\ast\ast};{C}^{2+\alpha}(\Gamma({t}_{\ast}))\cap{C}^{1}({I}_{\ast\ast};{C}^{1+\alpha}(\Gamma({t}_{\ast})).$
(16)

This follows directly from the following assumptions

$\displaystyle(y,t)$
$\displaystyle\in\Gamma\times{I}_{\ast\ast}\subset\Gamma\times{I}_{\ast}\subset\Gamma\times[a,b],$
$\displaystyle\rho$
$\displaystyle\in{{C}}^{0}({I}_{\ast};{{C}}^{2+\alpha}(\Gamma))\cap{{C}}^{1}({I}_{\ast};{{C}}^{1+\alpha}(\Gamma)),$
$\displaystyle X$ $\displaystyle\in{C}^{2+\alpha}(\Gamma\times
I_{\varepsilon_{0}}),$ $\displaystyle(x,t)$
$\displaystyle\in\Gamma({t}_{\ast})\times{I}_{\ast\ast}\subset\Gamma({t}_{\ast})\times{I}_{\ast},$

as well as from the regularities of the following mappings:

$\displaystyle[(y,t)\mapsto X(y,\rho(y,t))]$
$\displaystyle\in{{C}}^{0}({I}_{\ast\ast};{C}^{2+\alpha}(\Gamma;\mathbb{R}^{d}))\cap{C}^{1}({I}_{\ast\ast};{C}^{1+\alpha}(\Gamma;\mathbb{R}^{d})),$
$\displaystyle{X}_{\ast}$
$\displaystyle\in\operatorname{Diffeo}^{2+\alpha}(\Gamma({t}_{\ast})\times{I}_{\ast\ast},{X}_{\ast}(\Gamma({t}_{\ast})\times{I}_{\ast\ast})),$
$\displaystyle{\rho}_{\ast}$
$\displaystyle\in{{C}}^{0}({I}_{\ast\ast};{{C}}^{2+\alpha}(\Gamma))\cap{{C}}^{1}({I}_{\ast\ast};{{C}}^{1+\alpha}(\Gamma)),$
$\displaystyle{x}_{\ast}$
$\displaystyle\in{{C}}^{2+\alpha}(\Gamma;\Gamma({t}_{\ast})),$
$\displaystyle{x}_{\ast}^{-1}$
$\displaystyle\in{{C}}^{2+\alpha}(\Gamma({t}_{\ast});\Gamma).$

In view of Definition 2.7, the composition of these mappings leads to (16). ∎

By the previous lemma, we can state an equivalent definition of the solution
of Problem 2.6 given as follows:

######  Definition 2.9.

We say $\mathscr{M}=\bigcup_{0\leqslant t\leqslant
T}\Gamma(t)\times\\{t\\}\subset\mathbb{R}^{d}\times\mathbb{R}$ a solution of
Problem 2.6, if for all ${t}_{\ast}\in[0,T]$, there exist $a<b$ and ${\tau}>0$
such that
$[{t}_{\ast}-{\tau},{t}_{\ast}+{\tau}]\cap[0,T]\subset[a,b]\subset[0,T]$ and
there exists a quasi-normal vector ${\boldsymbol{\mathbf{N}}}_{\ast}$ on
$\Gamma({t}_{\ast})$ such that

$\mathscr{M}_{[a,b]}\in{{M}}_{[a,b]}(\Gamma({t}_{\ast}),{\boldsymbol{\mathbf{N}}}_{\ast})\quad\text{where}\quad\mathscr{M}_{[a,b]}:=\bigcup_{a\leqslant
t\leqslant b}\Gamma(t)\times\\{t\\},$

and $u_{\text{D}}(\Gamma(t))$ and $u_{\text{N}}(\Gamma(t))$ defined by (8) and
(9) solve (14).

######  Remark 2.10.

The definition of the solution of Problem 2.6 given in Definition 2.9 is
indeed equivalent to Definition 2.7. To see this, we note that for all
$t\in[0,T]$, there exists a ${\tau}(t)>0$ such that
$[t-{\tau}(t),t+{\tau}(t)]\cap[0,T]\subset[a,b]\subset[0,T]$, and there exists
a quasi-normal vector $\boldsymbol{\mathbf{N}}$ on $\Gamma(t)$ such that

$\mathscr{M}_{[a,b]}\in{{M}}_{[a,b]}(\Gamma(t),\boldsymbol{\mathbf{N}}),$

and $u_{\text{D}}(\Gamma(t))$ and $u_{\text{N}}(\Gamma(t))$ defined by (8) and
(9) solve (14). Now, we observe that

$\emptyset\neq\mathscr{O}(t):=\left\\{\begin{aligned}
&(t-{\tau}(t),t+{\tau}(t))\cap(0,T)&\text{for }t&\in(0,T),\\\
&[0,{\tau}(0))&\text{for }t&=0,\\\ &(t-{\tau}(t),T]&\text{for
}t&=T.\end{aligned}\right.$

Note that $\mathscr{O}$ is open in $[0,T]$ and
$\bigcup_{t\in[0,T]}\mathscr{O}(t)=[0,T]$.

##  3 Main Results

Our main objective is to demonstrate the solvability (local-in-time) of system
(1) by applying the techniques outlined by Solonnikov in [15]. Specifically,
we refer to the proofs of [15, Thm. 1.1 and Thm. 3.1] and also to [43] for a
relevant application of this method.

Notations. Let $l$ be a nonnegative real number. We denote by
${{C}}^{0}([0,T];{{C}}^{l}(\overline{\Omega}))$ the space of continuous
functions with respect to

$(x,t)\in\left\\{(x,t)\mid t\in[0,T],x\in\overline{\Omega}\right\\}$

with the finite norm

$\max_{0\leqslant t\leqslant T}|u(\cdot,t)|^{(l)}_{\overline{\Omega}},$

where

$\displaystyle\left|{u}\right|^{(l)}_{\overline{\Omega}}$
$\displaystyle:=\left\|{u}\right\|_{{{C}}^{[l],l-[l]}(\overline{\Omega})}=\left|{u}\right|_{[l],l-[l];\,\overline{\Omega}}=\sum_{\left|{j}\right|<l}{{\max_{\overline{\Omega}}}}|{{D}}^{j}u(x)|+[u]^{(l)}_{\overline{\Omega}},$
$\displaystyle[u]^{(l)}_{\overline{\Omega}}$
$\displaystyle:=[u]_{[l],l-[l];\,\overline{\Omega}}=\sum_{\left|{j}\right|=[l]}\max_{x,{{\hat{x}}}\in\overline{\Omega}}\frac{|{{D}}^{j}u(x)-{{D}}^{j}u({{\hat{x}}})|}{|x-{{\hat{x}}}|^{l-[l]}}.$

For example, $\left|{u}\right|_{\overline{\Omega}}^{(0)}$ is the maximum norm
of $u(x)$, i.e.,
$\left|{u}\right|_{\overline{\Omega}}^{(0)}=\max_{x\in\overline{\Omega}}|u(x)|$.
Also, in a more familiar notation,
$\left|{u}\right|^{(k+\alpha)}_{\overline{\Omega}}:=\left\|{u}\right\|_{{{C}}^{k,\alpha}(\overline{\Omega})}$
where $l=k+\alpha$, $k=[l]\in\mathbb{N}\cup\\{0\\}$, $\alpha=[l]-l\in(0,1)$;
see [44, Eqs. (4.5)–(4.6), p. 53]. The spaces
${{C}}^{0}([0,T];{{C}}^{l}(\Sigma))$ and ${{C}}^{0}([0,T];{{C}}^{l}(\Gamma))$
are introduced in a similar manner. Here, ${{C}}^{[l],l-[l]}:={{C}}^{k}$ for
integer $l=k\in\mathbb{N}$. Throughout the paper, we use $c$ as a general
positive constant, meaning it can vary in value depending on the context.

Let $a$, $b$ be fixed real numbers such that $b>a$,
$k\in\mathbb{N}\cup\\{0\\}$, and $\alpha\in(0,1)$. For functions
$\varphi=\varphi(\xi,t)\in{{C}}^{0}([a,b];{{C}}^{k+\alpha}(\overline{\Omega}))$,
the norm

${\max_{a\leqslant\tau\leqslant
b}}\left|{\varphi(\cdot,\tau)}\right|^{(k+\alpha)}_{\overline{\Omega}},$

will appear many times in the paper (especially from Section 7 onwards and in
Appendix E.3). So, for economy of space, we introduce the following notation

$\left|{\varphi}\right|^{(k+\alpha)}_{\Xi;\,[a,b]}:={\max_{a\leqslant\tau\leqslant
b}}\left|{\varphi(\cdot,\tau)}\right|^{(k+\alpha)}_{\Xi},\qquad\Xi\in\\{\Omega,\overline{\Omega},\Gamma,\Sigma\\}.$

For nested maximum norms, we write

$\left|{\varphi}\right|^{\infty}_{\Xi;\,[a,b]}:={\max_{a\leqslant\tau\leqslant
b}}\,{\max_{\Xi}}\left|{\varphi(\cdot,\tau)}\right|,\qquad\Xi\in\\{\overline{\Omega},\Gamma,\Sigma\\}.$

Additionally, for functions $u=u(\xi,t)$,
$v=v(\xi,t)\in{{C}}^{0}([a,b];{{C}}^{k+\alpha}(\Xi))$,
$\Xi\in\\{\overline{\Omega},\Gamma,\Sigma\\}$, we introduce the following
specially defined norms

$\left\|{(u,v)}\right\|^{(k+\alpha)}_{\Xi;\,[a,b]}:=\left|{u}\right|^{(k+\alpha)}_{\Xi;\,[a,b]}+\left|{v}\right|^{(k+\alpha)}_{\Xi;\,[a,b]}={\max_{a\leqslant
t\leqslant
b}}\left|{u(\cdot,\tau)}\right|^{(k+\alpha)}_{\Xi}+{\max_{a\leqslant
t\leqslant b}}\left|{v(\cdot,\tau)}\right|^{(k+\alpha)}_{\Xi}.$

For a pair of functions with subscripts ${\cdot}_{\text{D}}$ and
${\cdot}_{\text{N}}$, we simply write

$\left\|{\varphi_{{\text{D},\text{N}}}}\right\|^{(k+\alpha)}_{\Xi;\,[a,b]}:=\left\|{(\varphi_{\text{D}},\varphi_{\text{N}})}\right\|^{(k+\alpha)}_{\Xi;\,[a,b]},$

especially when such a norm appears many times in a sequence of arguments.
Moreover, and specifically for the boundary data
$f\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\Sigma))$ and
$g\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Sigma))$, $[a,b]\subseteq[0,T]$, we
simply write

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\Sigma;\,[a,b]}:={\max_{a\leqslant\tau\leqslant
b}}\left|{f(\cdot,\tau)}\right|^{(2+\alpha)}_{\Sigma}+{\max_{a\leqslant\tau\leqslant
b}}\left|{g(\cdot,\tau)}\right|^{(1+\alpha)}_{\Sigma},\qquad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}:=\left\|{(f,g)}\right\|^{(2+\alpha)}_{\Sigma,\,[0,T]}.$

Additionally, for $k\in\mathbb{N}$ and
$\rho\in{{C}}^{0}([a,b];{{C}}^{k+\alpha}(\Xi))\cap{{C}}^{1}([a,b];{{C}}^{k-1+\alpha}(\Xi))$
and $\Xi\in\\{\overline{\Omega},\Gamma\\}$, we introduce the norm notation

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\rho\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(k+\alpha)}_{\Xi;\,[a,b]}:={\max_{a\leqslant\tau\leqslant
b}}\left|{\rho(\cdot,\tau)}\right|^{(k+\alpha)}_{\Xi}+{\max_{a\leqslant\tau\leqslant
b}}\left|{\frac{d}{d\tau}\rho(\cdot,\tau)}\right|^{(k-1+\alpha)}_{\Xi},$

and another special notation

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\rho\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(k+\alpha)}_{\overline{\Omega},\,\Gamma;\,[a,b]}:={\max_{a\leqslant\tau\leqslant
b}}\left|{\rho(\cdot,\tau)}\right|^{(k+\alpha)}_{\overline{\Omega}}+{\max_{a\leqslant\tau\leqslant
b}}\left|{\frac{d}{d\tau}\rho(\cdot,\tau)}\right|^{(k-1+\alpha)}_{\Gamma}.$

Finally, for any pair of well-defined functions $\varphi_{\text{D}}$ and
$\varphi_{\text{N}}$, we will extensively use the following special notations
for convenience:

${\varphi}_{\text{D}\text{N}}:=\varphi_{\text{D}}-\varphi_{\text{N}}\quad\text{and}\quad{\varphi}_{\text{N}\text{D}}:=\varphi_{\text{N}}-\varphi_{\text{D}}.$

For example, we write ${u}_{\text{D}\text{N}}=u_{\text{D}}-u_{\text{N}}$ and
${u}_{\text{D}\text{N}}(\Gamma)=u_{\text{D}}(\Gamma)-u_{\text{N}}(\Gamma)$.

Prior to stating the primary result we seek to establish in this study, we
observe that, in the original inverse geometry problem setup, $f$ and $g$ are
solely space-dependent. However, for the remainder of the discussion,
otherwise specified, we shall assume that these functions are also
(pseudo-)time-dependent. Specifically, we let $f=f(x,t)$ and $g=g(x,t)$, where
$x\in\Omega(t)$ and $t\in[0,T]$, in (1), and assume that they are both
positive-valued.

Our main theorem is formulated as follows.

######  Theorem 3.1.

Let the following assumption be satisfied:

(A1)

For some $\alpha\in(0,1)$,

$\Sigma,\,\Gamma=\Gamma_{\circ}\in{{C}}^{2+\alpha},\quad
f\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\Sigma)),\ f>0,\quad
g\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Sigma)),\ g>0,$

such that

$\dfrac{\partial{}}{\partial\nu}\left({u}_{\text{D}\text{N}}(\Gamma_{\circ})\right)>0,$
(17)

where $u_{\text{D}}$ and $u_{\text{N}}$ respectively solves (8) and (9) in
$\Omega(\Gamma_{\circ})$.

Then, there exists a unique solution $\Gamma(t)$, $u_{\text{D}}(x,t)$, and
$u_{\text{N}}(x,t)$ to (1) defined on some small time-interval
${{I}^{\star}}=[0,{t}^{\star}]$, where ${{t}^{\star}}<T$. The free surface
$\Gamma(t)$ is described by the equation

$x=\xi+\rho(\xi,t)\boldsymbol{\mathbf{N}}(\xi),\qquad\xi\in\Gamma,$ (18)

where $\xi$ is the local coordinate on the surface $\Gamma$ and
$\boldsymbol{\mathbf{N}}$ is a smooth vector field on $\Gamma$ such that
$\boldsymbol{\mathbf{N}}\cdot\nu_{\circ}\geqslant\nu_{\star}>0$, where
$\nu_{\circ}$ is the unit normal vector to the surface $\Gamma$ directed
inward the domain $\Omega(\Gamma)$. The function
$\rho\in{{C}}^{0}({{I}^{\star}};{{C}}^{2+\alpha}(\Gamma))$ has extra
smoothness with respect to the variable $t$; namely,
$\rho_{t}\in{{C}}^{0}({{I}^{\star}};{{C}}^{1+\alpha}(\Gamma))$. Meanwhile, the
functions $u_{\text{D}}(x,t)$ and $u_{\text{N}}(x,t)$ are defined in
$\Omega(t)$ for $t\in{{I}^{\star}}$ and both belong to the space
${{C}}^{0}({{I}^{\star}};{{C}}^{2+\alpha}(\overline{\Omega(t)}))$. Moreover,
the following estimate hold

$\left\|{{u_{{\text{D},\text{N}}}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\rho\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{{\Gamma;\,[0,t]}}\leqslant
c{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{{\Sigma;\,[0,t]}}\leqslant
c\,{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.$
(19)

for some constant $c>0$, for all $t\in{{I}^{\star}}$.

Given the short-time existence of solution to (1), we can also prove the
uniqueness of short-time solution to the system.

######  Theorem 3.2.

A solution of Problem 2.6 is unique.

###### Proof.

If there exist two solutions
$\mathscr{M}_{i}:=\bigcup\Gamma^{i}(t)\times\\{t\\}$, $i=1,2$, to Problem 2.6,
then we need to show that $\mathscr{M}_{1}=\mathscr{M}_{2}$. We prove the
assertion via a contradiction. Let us assume that
$\mathscr{M}_{1}\neq\mathscr{M}_{2}$. Then, there exists ${t}_{\ast}\in[0,T)$
and a sequence $\\{t_{k}\\}_{k=1}^{\infty}\in({t}_{\ast},T]$ such that

$\left\\{\begin{aligned}
&\mathscr{M}_{1}\big{|}_{[0,{t}_{\ast}]}=\mathscr{M}_{2}\big{|}_{[0,{t}_{\ast}]},\\\
&T\geqslant t_{1}>t_{2}>\cdots>{t}_{\ast},\quad\text{where \
$\lim_{k\to\infty}t_{k}={t}_{\ast}$,\ and}\\\
&\Gamma^{1}(t_{k})\neq\Gamma^{2}(t_{k}),\quad\text{for
$k=1,2,\ldots$}.\end{aligned}\right.$ (20)

Since $\Gamma({t}_{\ast}):=\Gamma^{1}({t}_{\ast})=\Gamma^{2}({t}_{\ast})$
satisfies the condition in Theorem 3.1, there exists
${t}_{\ast\ast}\in({t}_{\ast},T]$ such that there is a unique $\Gamma(t)$ for
$t\in[{t}_{\ast},{t}_{\ast\ast}]$. However, this contradicts the last two
lines in (20). Thus, $\mathscr{M}_{1}$ and $\mathscr{M}_{2}$ have to be the
same solution to Problem 2.6. ∎

Before we move forward, we will further comment on condition (17) (cf. (A3)).
For this purpose, we first consider the following lemma.

######  Lemma 3.3.

Let $\Omega\subset\mathbb{R}^{d}$, of class ${{C}}^{2+\alpha}$, be an open
bounded connected set with non-intersecting boundaries $\Gamma$ and $\Sigma$.
Assume that
$v\in{{C}}^{2+\alpha}(\Omega)\cap{{C}}^{0+\alpha}(\overline{\Omega})$ and

$-\Delta{v}=0\quad\text{in $\Omega$},\qquad{v}=0\quad\text{on
$\Gamma$},\qquad\dfrac{\partial{v}}{\partial\nu}>0\quad\text{on $\Sigma$}.$

Then,

${v}>0\quad\text{in $\Omega$}.$

###### Proof.

The proof proceeds via contradiction and using the maximum principle [45,
Chap. 6.4, p. 344]. We start by observing that $v\not\equiv\text{constant}$ in
$\Omega$; otherwise, $v=0$ in $\Omega$ because $v=0$ on $\Gamma$, and so
$\dfrac{\partial{v}}{\partial\nu}=0$ on $\Sigma$, which contradicts the
assumption.

Now, we suppose that

$\inf_{\Omega}{v}<0.$

Then, by the maximum principle,

$0>\inf_{\Omega}{v}=\min_{\bar{\Omega}}{v}=\min_{\Sigma}{v}.$

That is,

$\text{there exists }x_{0}\in\Sigma\text{ such that
}{v}(x_{0})=\min_{\Sigma}{v}=\min_{\bar{\Omega}}{v}.$

Because $v=0$ on $\Gamma$, we deduce that

$\dfrac{\partial{v}}{\partial\nu}(x_{0})<0,\qquad x_{0}\in\Sigma.$

This is a contradiction to our assumption that
$\dfrac{\partial{v}}{\partial\nu}>0$ on $\Sigma$. Therefore,
$\inf_{\Omega}{v}>0$, and thus $v>0$ in $\Omega$, proving the lemma. ∎

Using Lemma 3.3, we will prove in the next proposition that there exists a
suitable choice of $\Gamma_{\circ}$ such that condition (17) holds. In fact,
we will use the same idea as in Appendices A.1 and A.2, which compares the
normal derivatives on the free boundary of concentric circles and spheres.

######  Proposition 3.4.

Let $\Omega=D\setminus\overline{\omega}\subset\mathbb{R}^{d}$, of class
${{C}}^{2+\alpha}$, be an open bounded connected set with non-intersecting
boundaries $\partial{\omega}=\Gamma\in\mathscr{A}^{2+\alpha}$ and
$\Sigma=\partial{D}$. Assume that
$\partial{\omega^{\star}}=\Gamma^{\star}\in\mathscr{A}^{2+\alpha}$ is the
exact interior boundary that satisfies (2) and $\omega$ strictly contained
$\overline{\omega}^{\star}$ (i.e., $\Gamma$ lies entirely in the interior of
$\overline{\Omega}^{\star}=\overline{D\setminus\overline{\omega}^{\star}}$).
Let $f\in{{C}}^{2+\alpha}(\Sigma)$ and $g\in{{C}}^{1+\alpha}(\Sigma)$. Then,
the functions $u_{\text{D}}(\Gamma)$ and $u_{\text{D}}(\Gamma)$ satisfying (8)
and (9), respectively, satisfy the following condition

$u_{\text{D}}>u_{\text{N}}\quad\text{in $\Omega$}.$

Consequently,

$\dfrac{\partial{}}{\partial\nu}\left(u_{\text{D}}-u_{\text{N}}\right)>0\quad\text{on
$\Gamma$}.$

###### Proof.

Let the assumptions of the proposition be satisified. Let us denote by
$(\Omega^{\star},u^{\star}(\Omega^{\star}))$,
$\Omega^{\star}=D\setminus\overline{\omega}^{\star}$, $\Omega^{\star}$ is of
class ${{C}}^{2+\alpha}$,
$u^{\star}\in{{C}}^{2+\alpha}(\Omega)\cap{{C}}^{0+\alpha}(\overline{\Omega})$,
the exact solution pair of the free boundary problem (2) with the
corresponding exact interior boundary
$\Gamma^{\star}=\partial{\omega}^{\star}\in\mathscr{A}^{2+\alpha}$. We note
the following observation

$\left\\{\begin{array}[]{rcll}\Delta u^{\star}=\Delta u_{\text{D}}=\Delta
u_{\text{N}}&=&0,&\quad\text{in}\ \Omega\\\
u^{\star}=u_{\text{D}}&>&0,&\quad\text{on}\ \Sigma,\\\
u_{\text{D}}=u_{\text{N}}&=&0,&\quad\text{on}\ \Gamma,\\\
u^{\star}&>&0,&\quad\text{on}\ \Gamma\end{array}\right.$ (21)

Let us define

${v}_{\text{D}}:=-u_{\text{D}}+u^{\star}\qquad\text{and}\qquad{v}_{\text{N}}:=-u_{\text{N}}+u^{\star}.$

Then, from (21), we have the following

$\left\\{\begin{array}[]{rcll}\Delta{v}_{\text{D}}=\Delta{v}_{\text{D}}&=&0,&\quad\text{in}\
\Omega\\\ {v}_{\text{D}}&=&0,&\quad\text{on}\ \Sigma,\\\
\dfrac{\partial{{v}_{\text{N}}}}{\partial\nu}&=&0,&\quad\text{on}\ \Sigma,\\\
{v}_{\text{D}}={v}_{\text{N}}=u^{\star}&>&0,&\quad\text{on}\
\Gamma\end{array}\right.$ (22)

We deduce from (22) that

${v}_{\text{D}}>0\quad\text{in
$\Omega$}\qquad\text{and}\qquad\dfrac{\partial{v_{\text{D}}}}{\partial\nu}<0\quad\text{on
$\Sigma$}.$ (23)

Now, we define

${v}={v}_{\text{N}}-{v}_{\text{D}}=(-u_{\text{N}}+u^{\star})-(-u_{\text{D}}+u^{\star})=u_{\text{D}}-u_{\text{N}}.$

Hence, we have

$-\Delta{v}=0\quad\text{in $\Omega$},\qquad{v}=0\quad\text{on
$\Gamma$},\qquad\dfrac{\partial{v}}{\partial\nu}>0\quad\text{on $\Sigma$}.$

Therefore, by Lemma 3.3, we get

$u_{\text{D}}-u_{\text{N}}={v}>0\quad\text{in $\Omega$}.$

Consequently, when $\Gamma\in\mathscr{A}^{2+\alpha}$ lies entirely in the
interior of
$\overline{\Omega}^{\star}=\overline{D\setminus\overline{\omega}^{\star}}$, we
get

$\dfrac{\partial{}}{\partial\nu}\left(u_{\text{D}}-u_{\text{N}}\right)>0\quad\text{on
$\Gamma$},$

which concludes the proof. ∎

Theorem 3.1 follows from an auxiliary result (see Theorem 4.4) which we issue
in the next section.

##  4 Problem Transformation onto a Fixed Domain

To carry out our analysis, we first need to transform system (1) into a
problem over a fixed domain via the change of variables to be described below.
Such technique has been used in many studies (see, e.g., [15]) and can be
achieved by a special mapping – a modification of the well-known Hanzawa
transform – of the pseudo-time dependent domain $\Omega(t)$ onto the fixed
domain $\Omega$. The target equation is given by (30), and the main result we
want to state here is emphasized in Theorem 4.4.

From this point forward, we assume that Assumption (A1) holds true.
Additionally, we consider $\boldsymbol{\mathbf{N}}$ to be a quasi-normal
vector on $\Gamma$. Also, for some technical purposes, we assume there is a
constant $\nu_{\star}>0$ such that

$\boldsymbol{\mathbf{N}}(\xi)\cdot\nu_{\circ}(\xi)\geqslant\nu_{\star}>0,\qquad\xi\in\Gamma,$
(A2)

where $\nu_{\circ}$ denotes the unit inward normal vector to $\Gamma$. This
requirement will be made clear in Section 6.

Given a constant $\varepsilon_{\circ}>0$, it was shown in [15, Sec. 3, p.134]
that the tubular strip

${\mathscr{S}}_{\varepsilon_{\circ}}=\\{\xi+\boldsymbol{\mathbf{N}}(\xi)\lambda\mid\left|{\lambda}\right|<\varepsilon_{\circ},\xi\in\Gamma\\}$

contains some neighborhood of $\Gamma$ where the equation
$x=\xi+\boldsymbol{\mathbf{N}}(\xi)\lambda$ determines some functions $\xi(x)$
and $\lambda(x)$ of class ${{C}}^{2+\alpha}$ (see the proof of Proposition
2.4). For sufficiently small $t>0$, the free boundary $\Gamma(t)$ is contained
in ${\mathscr{S}}_{\varepsilon_{\circ}}$ (cf. Proposition 2.4) and can be
described by the equation $\lambda=\rho(\xi,t)$. Accordingly, for

$\rho\in
R_{\circ}(\Gamma,\boldsymbol{\mathbf{N}}):=\\{\rho\in{{C}}^{2+\alpha}(\Gamma;\mathbb{R})\mid\left\|{\rho}\right\|_{\infty}<\varepsilon_{\circ}\\}$

we define

${\mathscr{S}}(\rho):=\\{x\in\mathbb{R}^{d}\mid
x=\xi+\boldsymbol{\mathbf{N}}(\xi)\rho(\xi),\ \xi\in\Gamma\\}.$ (24)

We introduce the admissible set of hypersurfaces consisting of $d-1$
dimensional ${{C}}^{2+\alpha}$ manifold embedded in $\mathbb{R}^{d}$ as
follows

${{H}}_{\varepsilon_{\circ}}(\Gamma,\boldsymbol{\mathbf{N}}):=\\{{\mathscr{S}}(\rho)\mid\rho\in
R_{\circ}(\Gamma,\boldsymbol{\mathbf{N}})\\}.$

We define a family of admissible moving surfaces (or moving boundaries) in the
following definition.

######  Definition 4.1.

We say that a moving surface (or moving boundary) $\mathscr{M}$ is admissible,
and we write $\mathscr{M}\in{{M}}_{[0,T]}(\Gamma,\boldsymbol{\mathbf{N}})$ if
and only if

$\left\\{\begin{aligned} \ \mathscr{M}&=\bigcup_{0\leqslant t\leqslant
T}\Gamma(t)\times\\{t\\}\subset\mathbb{R}^{d}\times\mathbb{R},\\\ \
\Gamma(t)&={\mathscr{S}}(\rho(t))\in{{H}}_{\varepsilon_{\circ}}(\Gamma,\boldsymbol{\mathbf{N}}),\
\text{for all }t\in[0,T],\\\ \
\rho&\in{{R}}_{[0,T]}(\Gamma,\boldsymbol{\mathbf{N}}),\end{aligned}\right.$

where ${{R}}_{[0,T]}(\Gamma,\boldsymbol{\mathbf{N}})$ is the set given in (12)
with $\varepsilon_{0}=\varepsilon_{\circ}$.

For every function $\rho\in{{R}}_{[0,T]}(\Gamma,\boldsymbol{\mathbf{N}})$, we
put into correspondence its extension

$\tilde{\rho}(y,t):={E}\rho(y,t),$ (25)

where $E$ is a linear and bounded map

${E}:{{R}}_{[0,T]}(\Gamma,\boldsymbol{\mathbf{N}})\to{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))\cap{{C}}^{1}([0,T];{{C}}^{1+\alpha}(\overline{\Omega}))$
(26)

satisfying

$\left\|{{E}{\rho}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant
c{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\rho\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{{\Gamma;\,[0,t]}},$

for some constant $c>0$. Such an extension of $\rho$ from $\Gamma\times[0,T]$
to $\overline{\Omega}\times[0,T]$, can be constructed in different ways. In
this investigation, we assume that
$\tilde{\rho}={\tilde{\rho}}(y,t):=E\rho(y,t)$ satisfies the following
boundary value problem.

$-\Delta{\tilde{\rho}}=0,\ y\in\Omega,\
t>0,\qquad\left.{\tilde{\rho}}\right|_{y\in\Gamma}=\rho,\qquad\left.{\tilde{\rho}}\right|_{y\in\Sigma}=0,$
(27)

where the latter condition is essential because $\Sigma$ is fixed. Such a
choice of extension will streamline many of our succeeding arguments. In the
rest of our discussion, we will always assume that ${\tilde{\rho}}$ satisfies
(27) without further notice.

To meet certain technical requirements, we will extend the vector field
$\boldsymbol{\mathbf{N}}$ (retaining the same notation) to the domain
${\mathscr{S}}_{\varepsilon_{\circ}}$ by defining
$\boldsymbol{\mathbf{N}}(x)=\boldsymbol{\mathbf{N}}(\xi(x))$. Then, we will
further extend it to the domain
$\Omega\cup{\mathscr{S}}_{\varepsilon_{\circ}}$, where $\Omega=\Omega(0)$,
ensuring that the smoothness class
$\boldsymbol{\mathbf{N}}\in{{C}}^{2+\alpha}(\Omega\cup{\mathscr{S}}_{\varepsilon_{\circ}})$
is maintained (cf. Remark 2.3).

Let us introduce the change of variables $y=Y(x,t)$ under which $\Omega(t)$,
$t>0$, is transformed to $\Omega$ and define the inverse transform
$Y^{-1}:\Omega\to\Omega(t)$ as follows:

$Y^{-1}(y,t)=x\in\Omega(t),\quad
x=y+\boldsymbol{\mathbf{N}}(y){\tilde{\rho}}(y,t),$ (28)

for $y\in\Omega$ and $t>0$. The given map is bijective, as stated in the
following lemma.

######  Lemma 4.2.

Let $Z(y):=y+\boldsymbol{\mathbf{N}}(y)\tilde{\rho}(y)$,
$y\in\overline{\Omega}$, such that $Z\left.{=}\right|_{y\in\Sigma}y$ and
$Z\big{|}_{y\in\Gamma}=y+\boldsymbol{\mathbf{N}}(y)\rho(y)$. Then, the map
$Z:\overline{\Omega}\to\overline{\Omega}(\rho)$, where
$\overline{\Omega}(\rho)$ is the annular domain bounded by $\Sigma$ and
${\mathscr{S}}(\rho)$ is bijective.

###### Proof.

See Appendix E.2. ∎

To transform the Hele-Shaw-like system (1) onto a fixed domain, we make a few
additional preparations. The Jacobi matrix $\textsf{J}=(\textsf{J}_{km})_{km}$
($k,m=1,\ldots,d$) of the transform $Y^{-1}$ has entries

$\textsf{J}_{km}=(D_{y}Y^{-1}(y,t))_{km}=\delta_{km}+N_{k}\frac{\partial{\tilde{\rho}}}{\partial
y_{m}}+\frac{\partial N_{k}}{\partial y_{m}}{\tilde{\rho}}.$

Correspondingly, we denote by $\textsf{J}^{km}$ the entries of the inverse
matrix $\textsf{J}^{-1}$, i.e.,

$\textsf{J}^{-1}=\left(\textsf{J}^{km}\right)_{km}=\left((D_{x}Y(x,t))_{km}\right)_{km}=\left((D_{y}Y^{-1}(y,t))_{km}\right)_{km}^{-1},$

which is the Jacobi matrix of the transform $Y(x,t)$. Denoting
$\textsf{J}^{-\top}=(\textsf{J}^{-1})^{\top}$, the operator
$\nabla_{x}:=\left(\partial_{x_{1}},\ldots,\partial_{x_{d}}\right)$,
$\partial_{x_{k}}:=\partial/\partial{x_{k}}$, $k=1,\ldots,d$, takes the form

$\textsf{J}^{-\top}\nabla_{y}^{\top}:=\left(\sum_{k=1}^{d}\textsf{J}^{km}\frac{\partial}{\partial
y_{k}}\right),\quad
m=1,\ldots,d,\qquad\nabla_{y}:=\left(\frac{\partial}{\partial
y_{1}},\ldots,\frac{\partial}{\partial y_{d}}\right).$ (29)

The following relation between the inward unit normal vectors $\nu_{\circ}$ to
$\Gamma$ and $\nu(t)$ to $\Gamma(t)$ (cf. [15, p. 135] or see [1, Thm. 4.4, p.
488]) holds:

$\nu(t)\circ
Y^{-1}=\frac{\textsf{J}^{-\top}\nu_{\circ}}{|\textsf{J}^{-\top}\nu_{\circ}|}\qquad\Longleftrightarrow\qquad\nu(t)=\frac{\textsf{J}^{-\top}\nu_{\circ}}{|\textsf{J}^{-\top}\nu_{\circ}|}\circ
Y.$

Applying the above change of variables and identities, we can pass from the
$t$-dependent system (1) to the problem in the given fixed domain $\Omega$
with respect to the three unknown functions
$U_{\text{D}}(y,t)=u_{\text{D}}(Y^{-1}(y,t),t)$,
$U_{\text{N}}(y,t)=u_{\text{N}}(Y^{-1}(y,t),t)$, and ${\tilde{\rho}}(y,t)$,
$y\in\Omega$. More precisely, given the harmonic function $\tilde{\rho}$
satisfying (27), we have the following transformation of (1) using the map
$Y:\Omega(t)\to\Omega$ where $Y(x,t)=y\in\Omega$, for $t>0$:

$\left\\{\begin{aligned}
\displaystyle\sum_{m,p=1}^{d}{{\textsf{A}}}_{mp}\frac{\partial^{2}U_{\text{D}}}{\partial
y_{m}\partial
y_{p}}+\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{D}}}{\partial y_{p}}&=0,\quad y\in\Omega,\
t>0,\\\
\left.{U_{\text{D}}}\right|_{y\in\Sigma}=f(y,t),\qquad\left.{U_{\text{D}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\displaystyle\sum_{m,p=1}^{d}{{\textsf{A}}}_{mp}\frac{\partial^{2}U_{\text{N}}}{\partial
y_{m}\partial
y_{p}}+\sum_{k,m,p=1}^{d}\textsf{J}^{mk}\frac{\partial\textsf{J}^{pk}}{\partial
y_{m}}\frac{\partial U_{\text{N}}}{\partial y_{p}}&=0,\quad y\in\Omega,\
t>0,\\\
\left.{\dfrac{\partial{U_{\text{N}}}}{\partial\nu}}\right|_{y\in\Sigma}=g(y,t),\qquad\left.{U_{\text{N}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\left.{\left(\dfrac{\partial{\tilde{\rho}}}{\partial{t}}+B_{{\tilde{\rho}}}\dfrac{\partial{{{{U}_{\text{D}\text{N}}}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}&=0,\\\
{\tilde{\rho}}\big{|}_{y\in\Gamma,\,t=0}&=0.\end{aligned}\right.$ (30)

where

${{\textsf{A}}}_{mp}:=\sum_{k=1}^{d}\textsf{J}^{mk}\textsf{J}^{pk}$

are the entries of the matrix
${{\textsf{A}}}:=\textsf{J}^{-1}\textsf{J}^{-\top}$, and

$B_{\tilde{\rho}}:=\left(\boldsymbol{\mathbf{N}}\cdot\textsf{J}^{-\top}\nu_{\circ}\right)^{-1}\left(\nu_{\circ}\cdot{{\textsf{A}}}\nu_{\circ}\right).$
(31)

We note here that
$\dfrac{\partial{{{{U}_{\text{D}\text{N}}}}}}{\partial{\nu_{\circ}}}>0$ in
accordance with Assumption (A1).

For convenience of later use, we will write problem (30) in compact form. For
this purpose, for $\rho\in{{C}}^{2+\alpha}(\Gamma)$, we define the following
linear and bounded operators:

$\left\\{\begin{aligned}
\pazocal{{L}}_{\tilde{\rho}}&:{{C}}^{2+\alpha}(\overline{\Omega})\to{{C}}^{\alpha}(\overline{\Omega}),\quad\qquad\pazocal{{L}}_{\tilde{\rho}}=\sum_{m,p=1}^{d}{{\textsf{A}}}_{mp}\frac{\partial^{2}}{\partial
y_{m}\partial y_{p}},\\\
\pazocal{{K}}_{\tilde{\rho}}&:{{C}}^{2+\alpha}(\overline{\Omega})\to{{C}}^{1+\alpha}(\overline{\Omega}),\qquad\pazocal{{K}}_{\tilde{\rho}}=-\sum_{j,k=1}^{d}{N}_{k}\textsf{J}^{jk}\frac{\partial}{\partial
y_{j}},\\\
\pazocal{{M}}_{\tilde{\rho}}&:{{C}}^{2+\alpha}(\overline{\Omega})\to{{C}}^{1+\alpha}(\overline{\Omega}),\\\
\quad\pazocal{{M}}_{\tilde{\rho}}&=-\sum_{p=1}^{d}\left[\sum_{j,k,m,q=1}^{d}\textsf{J}^{km}\textsf{J}^{pj}\textsf{J}^{qm}\left(\frac{\partial{N}_{j}}{\partial
y_{q}}\frac{\partial{\tilde{\rho}}}{\partial
y_{k}}+\frac{\partial{N}_{j}}{\partial
y_{k}}\frac{\partial{\tilde{\rho}}}{\partial
y_{q}}+\frac{\partial^{2}{N}_{j}}{\partial y_{k}\partial
y_{q}}{\tilde{\rho}}\right)\right]\frac{\partial}{\partial
y_{p}}.\end{aligned}\right.$ (32)

The main equations in system (30) can then be written as follows

$\pazocal{{L}}_{\tilde{\rho}}{U_{\text{D}}}+\pazocal{{K}}_{\tilde{\rho}}{U_{\text{D}}}\pazocal{{L}}_{\tilde{\rho}}{{\tilde{\rho}}}+\pazocal{{M}}_{\tilde{\rho}}{U_{\text{D}}}=0\quad\text{and}\quad\pazocal{{L}}_{\tilde{\rho}}{U_{\text{N}}}+\pazocal{{K}}_{\tilde{\rho}}{U_{\text{N}}}\pazocal{{L}}_{\tilde{\rho}}{{\tilde{\rho}}}+\pazocal{{M}}_{\tilde{\rho}}{U_{\text{N}}}=0.$
(33)

We note here that $\pazocal{{M}}_{\tilde{\rho}}\equiv 0$ when
${\tilde{\rho}}\equiv 0$.

######  Lemma 4.3.

The following regularities hold:

$\left\\{\begin{aligned}
&\in{{C}}^{\omega}({{C}}^{2+\alpha}(\Gamma),\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{\alpha}(\overline{\Omega})))\\\
[\rho\mapsto\pazocal{{K}}_{\tilde{\rho}}],\
[\rho\mapsto\pazocal{{M}}_{\tilde{\rho}}]&\in{{C}}^{\omega}({{C}}^{2+\alpha}(\Gamma),\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{1+\alpha}(\overline{\Omega}))).\end{aligned}\right.$
(34)

###### Proof.

We will only our claim that

$[\rho\mapsto\pazocal{{L}}_{\tilde{\rho}}]\in{{C}}^{\omega}({{C}}^{2+\alpha}(\Gamma),\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{\alpha}(\overline{\Omega})))$
(35)

The other two results can be proven in a similar fashion.

We let $\rho\in{{C}}^{2+\alpha}(\Gamma)$. In view of (25) and (27), we define

$\tilde{\rho}(y):={E}\rho(y),\quad y\in\Gamma,\qquad\text{where}\
E\in\mathscr{B}({{C}}^{2+\alpha}(\Gamma);{{C}}^{2+\alpha}(\overline{\Omega})),$

and such that

$-\Delta{\tilde{\rho}}=0,\
y\in\Omega,\qquad\left.{\tilde{\rho}}\right|_{y\in\Gamma}=\rho,\qquad\left.{\tilde{\rho}}\right|_{y\in\Sigma}=0.$

We claim that

$\pazocal{{L}}_{\tilde{\rho}}\in\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{\alpha}(\overline{\Omega})).$
(36)

Let us note that the following regularities hold:

  * •

$\left[\tilde{\rho}\mapsto\dfrac{\partial{\tilde{\rho}}}{\partial{y_{m}}}\right]\in\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{1+\alpha}(\overline{\Omega}))$,
for all $m=1,\ldots,d$;

  * •

$D_{mp}:=\dfrac{\partial^{2}}{\partial_{y_{m}}\partial_{y_{p}}}\in\mathscr{B}({{C}}^{2+\alpha}(\overline{\Omega}),{{C}}^{\alpha}(\overline{\Omega}))$,
for all $m,p=1,\ldots,d$;

  * •

$A_{mp}\in{{C}}^{\omega}({{C}}^{2+\alpha}(\overline{\Omega}));{{C}}^{1+\alpha}(\overline{\Omega})))$,
for all $m,p=1,\ldots,d$;

  * •

for $u\in{{C}}^{2+\alpha}(\overline{\Omega}))$, we have
$\Lambda_{mp}\in{{C}}^{\omega}({{C}}^{\alpha}(\overline{\Omega}));{{C}}^{\alpha}(\overline{\Omega})))$,
for all $m,p=1,\ldots,d$, where
$\Lambda_{mp}^{\rho}(x):=A_{mp}(\tilde{\rho})(x)u(x)$,
$x\in\overline{\Omega}$.

The latter regularity result follows from the fact that the following map

$\mathscr{F}:{{C}}^{\alpha}(\overline{\Omega})\times{{C}}^{\alpha}(\overline{\Omega})\longrightarrow{{C}}^{\alpha}(\overline{\Omega}),\qquad(\varphi,\psi)\longmapsto\varphi\psi,$

is analytic (cf. [44, p. 53]); i.e.
$\mathscr{F}\in{{C}}^{\omega}({{C}}^{\alpha}(\overline{\Omega})\times{{C}}^{\alpha}(\overline{\Omega});{{C}}^{\alpha}(\overline{\Omega})))$.
Hence, using $\Lambda_{mp}^{\rho}(x)=\mathscr{F}(A_{mp}(\tilde{\rho}),u)$, we
establish the final regularity result mentioned above. Consequently, by
composing the given maps, we validate (36).

To complete the proof, we use the composition of the map
$\rho\mapsto\tilde{\rho}$ (defined through the extension operator $E$) and the
operator $\pazocal{{L}}_{\tilde{\rho}}$, which allows us to confirm our claim
(35). ∎

Finally, our main result given by Theorem 3.1 follows from the result given
next asserting the local-in-time solvability of (30).

######  Theorem 4.4.

Let Assumption (A1) hold true. Specifically,

$\dfrac{\partial{}}{\partial{\nu_{\circ}}}\left({{{U}_{\text{D}\text{N}}(y,0)}}\right)>0,\qquad
y\in\Gamma.$ (A3)

Then, problem (30) has a unique solution
$({\tilde{\rho}}(y,t),U_{\text{D}}(y,t),U_{\text{N}}(y,t))$, $y\in\Omega$,
that is defined for $t\in{{I}^{\star}}$ with ${{t}^{\star}}<T$, and such that
the following regularities hold:

${\tilde{\rho}}\in{{C}}^{0}({{I}^{\star}};{{C}}^{2+\alpha}(\overline{\Omega})),\quad{\tilde{\rho}}_{t}\big{|}_{\Gamma}\in{{C}}^{0}({{I}^{\star}};{{C}}^{1+\alpha}(\Gamma)),\quad
U_{\text{D}}(y,t),\
U_{\text{N}}(y,t)\in{{C}}^{0}({{I}^{\star}};{{C}}^{2+\alpha}(\overline{\Omega})).$

Moreover, the following estimate holds

$\left\|{{U_{{\text{D},\text{N}}}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\tilde{\rho}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{{\Sigma;\,[0,t]}}\leqslant
c\,{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(f,g)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},$
(37)

for some constant $c>0$, for all $t\in{{I}^{\star}}$.

##  5 Regularity of Solutions on a Fixed Domain

Our immediate aim is to differentiate between the linear parts of the main
conditions and the dynamic boundary conditions in the nonlinear problem stated
in Equation (47). Part of the analysis requires understanding the unique
solvability of a pure Dirichlet problem and a mixed Dirichlet-Neumann problem,
as well as obtaining estimates for the solutions of such equations. In this
section, we lay the groundwork for these preparations.

Let $u_{\text{D}}^{0}:=u_{\text{D}}(x,0)$,
$f_{0}:=f(\cdot,0)\in{{C}}^{2+\alpha}(\Sigma)$, and consider the pure
Dirichlet problem

$\Delta u_{\text{D}}^{0}=0\quad\text{in}\ \Omega,\qquad
u_{\text{N}}^{0}=f_{0}\quad\text{on}\ \Sigma,\qquad
u_{\text{D}}^{0}=0\quad\text{on}\ \Gamma.$ (38)

The solvability of (38) in the space ${{C}}^{k+\alpha}(\overline{\Omega})$
($\Omega$ a bounded region) is well-known and is given, for instance, in [42,
Thm. 1.1, p. 107] or [46, Chap. V.36, I, p. 166] (see also Kellogg’s Theorem
[47, 48]). In particular, we have the following result.

######  Lemma 5.1.

For any $f_{0}\in{{C}}^{2+\alpha}(\Sigma)$, there exists a unique solution
$u_{\text{D}}^{0}\in{{C}}^{2+\alpha}(\overline{\Omega})$ to (38).

###### Proof.

See, for example, [44, Thm. 6.14, p. 107]. ∎

######  Remark 5.2.

It is worth to emphasize here that, for $k=0,1,\ldots$ and $\alpha\in(0,1)$,
the space ${{C}}^{k+\alpha}(\Omega)$ essentially coincides with
${{C}}^{k+\alpha}(\overline{\Omega})$ for any open (bounded) set
$\Omega\subset\mathbb{R}^{d}$ [49, Prop. 1.1.7, p. 6].

Now, let $u_{\text{N}}^{0}:=u_{\text{N}}(x,0)$,
$g_{0}:=g(\cdot,0)\in{{C}}^{1+\alpha}(\Sigma)$, and consider the following
mixed Dirichlet-Neumann boundary value problem

$\Delta u_{\text{N}}^{0}=0\quad\text{in}\
\Omega,\qquad\dfrac{\partial{u_{\text{N}}^{0}}}{\partial{\nu_{\circ}}}=g_{0}\quad\text{on}\
\Sigma,\qquad u_{\text{N}}^{0}=0\quad\text{on}\ \Gamma.$ (39)

In the next proposition we want to prove that
$u_{\text{N}}^{0}\in{{C}}^{2+\alpha}(\overline{\Omega})$ provided that
$g_{0}\in{{C}}^{1+\alpha}(\Sigma)$.

######  Lemma 5.3.

Let $\Omega\subset\mathbb{R}^{d}$ be an annular open bounded set with
${{C}}^{2+\alpha}$ boundary $\partial\Omega=\Sigma\cup\Gamma$ and
$g_{0}\in{{C}}^{1+\alpha}(\Sigma)$. Then Equation (39) has the unique solution
$u_{\text{N}}^{0}\in{{C}}^{2+\alpha}(\overline{\Omega})$.

###### Proof.

Let us choose an annular connected open bounded set $\Omega_{1}$ with exterior
boundary $\Gamma$ and interior boundary $S$ of class
${{C}}^{0,\alpha}\cap{{C}}^{\infty}$. We denote by $\Omega_{2}$ the annular
open bounded set contained in $\Omega$ whose boundary are $S$ and $\Sigma$.
Obviously, $\Omega_{1}\cap\Omega_{2}=\emptyset$ and
$\Omega_{1}\cup\Omega_{2}=\Omega$. We then consider the following pure
Dirichlet boundary value problem

$\Delta{{V}}=0\quad\text{in}\
\Omega_{1},\qquad{{V}}=u_{\text{N}}^{0}\quad\text{on}\
S,\qquad{{V}}=0\quad\text{on}\ \Gamma.$ (40)

By [45, Thm. 3, p. 316], we know that
$u_{\text{N}}^{0}\in{{C}}^{\infty}(\Omega)$, in particular,
$u_{\text{N}}^{0}\in{{C}}^{2}(\Omega)$. Hence,
$u_{\text{N}}^{0}\in{{C}}^{2+\alpha}(\overline{\Omega}_{1})$ by [44, Thm. 4.6,
p. 60]. This implies that the above problem admits a unique solution
$w\in{{C}}^{2+\alpha}(\overline{\Omega}_{1})$ because of [44, Thm. 6.8, p.
100] (see also [42, Thm. 1.3, p. 107]). Since
$u_{\text{N}}^{0}\big{|}_{\overline{\Omega}_{1}}$ also solves (40), then by
uniqueness, we have $w\equiv u_{\text{N}}^{0}\big{|}_{\overline{\Omega}_{1}}$.

Let us next consider the following boundary value problem:

$\Delta v-v=-u_{\text{N}}^{0}\quad\text{in}\
\Omega_{2},\qquad\dfrac{\partial{v}}{\partial{\nu_{\circ}}}=g_{0}\quad\text{on}\
\Sigma,\qquad\dfrac{\partial{v}}{\partial{\nu_{\circ}}}=\dfrac{\partial{u_{\text{N}}^{0}}}{\partial{\nu_{\circ}}}\quad\text{on}\
S.$ (41)

Since $u_{\text{N}}^{0}\in{{C}}^{2+\alpha}(\overline{\Omega}_{2})$, $\nabla
u_{\text{N}}^{0}\cdot{\nu}\in{{C}}^{1+\alpha}(S)$, and
$g_{0}\in{{C}}^{1+\alpha}(\Sigma)$, then by [50, Thm. 5.2], we see that (41)
admits a unique solution $v\in{{C}}^{2+\alpha}(\overline{\Omega}_{2})$.
Because $u_{\text{N}}^{0}\big{|}_{\overline{\Omega}_{2}}$ also solves (41),
uniqueness then implies that $v\equiv
u_{\text{N}}^{0}\big{|}_{\overline{\Omega}_{2}}$. This concludes that Equation
(39) admits a unique solution in ${{C}}^{2+\alpha}(\overline{\Omega})$. ∎

Let $U^{0}_{\text{D}}:=U^{0}_{\text{D}}(y,t)$, $y\in\Omega$, $t>0$, and
$f(\cdot,t)\in{{C}}^{2+\alpha}(\Sigma)$ for all $t\geqslant 0$. Then, it can
be checked that the boundary value problem

$\Delta U^{0}_{\text{D}}=0,\quad
y\in\Omega,\qquad\left.{U^{0}_{\text{D}}}\right|_{y\in\Sigma}=f(y,t),\qquad\left.{U^{0}_{\text{D}}}\right|_{y\in\Gamma}=0,$
(42)

admits a unique solution
$U^{0}_{\text{D}}\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$, for
all $t\geqslant 0$, where $t$ is a parameter appearing in the boundary
condition on $\Sigma$. The solution to (42) satisfies the inequality condition

${\max_{0\leqslant\tau\leqslant
t}}\left|{U^{0}_{\text{D}}(\cdot,\tau)}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c(d,\Omega,\alpha){\max_{0\leqslant\tau\leqslant
t}}\left|{f(\cdot,\tau)}\right|_{\Sigma}^{(2+\alpha)},$ (43)

Similarly, for $U^{0}_{\text{N}}:=U^{0}_{\text{N}}(y,t)$, $y\in\Omega$, $t>0$,
and $g(\cdot,t)\in{{C}}^{1+\alpha}(\Sigma)$ for all $t\geqslant 0$, the
boundary value problem

$\Delta U^{0}_{\text{N}}=0,\quad
y\in\Omega,\qquad\left.{\dfrac{\partial{U^{0}_{\text{N}}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=g(y,t),\qquad\left.{U^{0}_{\text{N}}}\right|_{y\in\Gamma}=0,$
(44)

admits a unique solution
$U^{0}_{\text{N}}\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$.
Again, $t$ is a parameter appearing in the boundary condition on $\Sigma$.

We want to prove next an estimate of
$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)}$ in terms of
$\left|{g}\right|_{\Sigma}^{(1+\alpha)}$. More precisely, we want to prove the
validity of the following proposition.

######  Proposition 5.4.

Let $\alpha\in(0,1)$ and $\Omega\subset\mathbb{R}^{d}$ be an open, bounded,
and connected set with ${{C}}^{2+\alpha}$ regularity and
$g\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Sigma))$ be given. If
$U^{0}_{\text{N}}\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$ is a
solution to (44), then we have the following estimate

${\max_{0\leqslant\tau\leqslant
t}}\left|{U^{0}_{\text{N}}(\cdot,\tau)}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c(d,\Omega,\alpha){\max_{0\leqslant\tau\leqslant
t}}\left|{g(\cdot,\tau)}\right|_{\Sigma}^{(1+\alpha)},$ (45)

where $c(d,\Omega,\alpha)>0$ is a constant that depends only on the dimension
$d$, the set $\Omega$, and the number $\alpha\in(0,1)$.

###### Proof.

Let $g\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Sigma))$ for some $\alpha\in(0,1)$,
$t\in[0,T]$, and $\tau<t$ be fixed. Note that it is enough to prove that
$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c\left|{g}\right|_{\Sigma}^{(1+\alpha)}$, for some constant $c>0$. In the
proof, we abuse some notations (i.e., we write $\varphi$ instead of
$\varphi(\xi,\tau)=:\varphi(\xi)$ as $\tau$ is fixed for the functions
involved here) and assume that
$U^{0}_{\text{N}}\in{{C}}^{2+\alpha}(\overline{\Omega})$ is a solution to
(44). We consider a cutoff function $\phi\in{{C}}^{\infty}(\mathbb{R}^{d})$
such that $\phi=1$ near $\Sigma$ and $\phi=0$ near $\Gamma$. Next, we define
the functions ${v}_{\text{N}}:=\phi U^{0}_{\text{N}}$ and
${v}_{\text{D}}:=(1-\phi)U^{0}_{\text{N}}$ such that

$\displaystyle\Delta{v}_{\text{N}}={h}_{\phi}\quad\text{in}\ \Omega,$
$\displaystyle\dfrac{\partial{{v}_{\text{N}}}}{\partial{\nu_{\circ}}}=g\quad\text{on}\in\Sigma,$
$\displaystyle\dfrac{\partial{{v}_{\text{N}}}}{\partial{\nu_{\circ}}}=0\quad\text{on}\
\Gamma;$ $\displaystyle\ \Delta{v}_{\text{D}}=-{h}_{\phi}\quad\text{in}\
\Omega,$ $\displaystyle{v}_{\text{D}}=0\quad\text{on}\in\Sigma,$
$\displaystyle{v}_{\text{D}}=0\quad\text{on}\ \Gamma,$

where ${h}_{\phi}:=(\Delta\phi)U^{0}_{\text{N}}+\nabla\phi\cdot\nabla
U^{0}_{\text{N}}$. By these constructions, notice that
$U^{0}_{\text{N}}={v}_{\text{N}}+{v}_{\text{D}}$. Now, by [44, Thm. 6.30, equ.
(6.77), p. 127], we immediately get the following estimates

$\displaystyle\left|{{v}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)}$
$\displaystyle\leqslant
c_{\text{0N}}(d,\Omega,\alpha)\left(\left|{{v}_{\text{N}}}\right|_{\overline{\Omega}}^{(0)}+\left|{g}\right|_{\Sigma}^{(1+\alpha)}+\left|{{h}_{\phi}}\right|_{\overline{\Omega}}^{(0+\alpha)}\right)\leqslant
c_{\text{1N}}(d,\Omega,\alpha)\left(\left|{g}\right|_{\Sigma}^{(1+\alpha)}+\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(1+\alpha)}\right),$
$\displaystyle\left|{{v}_{\text{D}}}\right|_{\overline{\Omega}}^{(2+\alpha)}$
$\displaystyle\leqslant
c_{\text{0D}}(d,\Omega,\alpha)\left|{{h}_{\phi}}\right|_{\overline{\Omega}}^{(0+\alpha)}\leqslant
c_{\text{1D}}(d,\Omega,\alpha)\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(0+\alpha)},$

where $c_{\text{0N}},c_{\text{1N}},c_{\text{0D}},c_{\text{1D}}>0$. These
inequalities clearly implies that

$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c_{0}(d,\Omega,\alpha)\left(\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(1+\alpha)}+\left|{g}\right|_{\Sigma}^{(1+\alpha)}\right).$

The quantity
$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(1+\alpha)}$ can be
estimated as follows (see, e.g., [44, Lem. 6.35, p. 135]):

$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(1+\alpha)}\leqslant
c_{1}(\varepsilon,\Omega)\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(0)}+\varepsilon\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)},$

for some constant $c_{1}(\varepsilon,\Omega)>0$. Using the above interpolation
inequality, we get

$\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c_{2}(n,\Omega,\alpha,\varepsilon)\left(\left|{U^{0}_{\text{N}}}\right|_{\overline{\Omega}}^{(0)}+\left|{g}\right|_{\Sigma}^{(1+\alpha)}\right),$
(46)

for some constant $c_{2}(n,\Omega,\alpha,\varepsilon)>0$, and we want to prove
that we actually have (45). To do this, we mimic an argument in showing
inequalities between two equivalent norms which is also similar to a proof of
Poincaré’s inequality via a contradiction (see, e.g., [45, Proof of Theorem 1,
pp. 275–276] or [51, Lem. 49.30, p. 1037]). Analogous proof can also be found
in [52, Proof of Theorem 3.28, pp. 194–195].

The argument to prove (45) will proceed by contradiction. So, let us suppose
that inequality (45) does not hold. This means that we can find a sequence
$\\{U^{0}_{\text{N}k}\\}_{k\in\mathbb{N}}\subset{{C}}^{2+\alpha}(\Omega)$ and
$\\{g_{k}\\}\subset{{C}}^{1+\alpha}(\Sigma)$ such that

$\Delta U^{0}_{\text{N}k}=0\quad\text{in}\
\Omega,\qquad\dfrac{\partial{U^{0}_{\text{N}k}}}{\partial{\nu_{\circ}}}={g_{k}}\quad\text{on}\in\Sigma,\qquad
U^{0}_{\text{N}k}=0\quad\text{on}\ \Gamma,$

and, for each $k\in\mathbb{N}$, $U^{0}_{\text{N}k}$ satisfies

$\left|{U^{0}_{\text{N}k}}\right|_{\overline{\Omega}}^{(2+\alpha)}=1\qquad\text{and}\qquad\left|{U^{0}_{\text{N}k}}\right|_{\overline{\Omega}}^{(2+\alpha)}\geqslant
k\left|{{g_{k}}}\right|_{\Sigma}^{(1+\alpha)}.$

Letting $k\to\infty$, we get ${g_{k}}\to 0$ in ${{C}}^{1+\alpha}(\Sigma)$
because
$\frac{1}{k}\left|{U^{0}_{\text{N}k}}\right|_{\overline{\Omega}}^{(2+\alpha)}\to
0$. By the first condition above, we see that for every multi-index $\beta$,
with $\left|{\beta}\right|=0,1,2$, the sequence
$\\{D^{\beta}U^{0}_{\text{N}k}\\}$ is uniformly bounded in
${{C}}^{0}(\overline{\Omega})$. Consequently, the sequence is equicontinuous
because the inequality

$\left|{D^{\beta}U^{0}_{\text{N}k}(\xi_{1})-D^{\beta}U^{0}_{\text{N}k}(\xi_{2})}\right|\leqslant\left|{\xi_{1}-\xi_{2}}\right|^{\alpha},\qquad\text{for
all}\ \xi_{1},\xi_{2}\in\Omega,\quad\text{for all}\ \left|{\beta}\right|=2,$

actually implies that

$\left|{D^{\beta}U^{0}_{\text{N}k}(\xi_{1})-D^{\beta}U^{0}_{\text{N}k}(\xi_{2})}\right|\leqslant
m_{0}\left|{\xi_{1}-\xi_{2}}\right|^{\alpha},\qquad\text{for all}\
\xi_{1},\xi_{2}\in\Omega,\quad\text{for all}\ \left|{\beta}\right|=0,\ 1,$

for some constant $m_{0}=m_{0}(\Omega)>0$ according to [49, Prop. 1.5.2, p.
34, Thm. 4.4.1, p. 109, or Thm. 4.4.2, p. 120]. By iteratively applying the
Arzelà-Ascoli theorem (see, e.g., [45, Sec. C.7, pp. 634–635] or [53, Thm.
2.86, p. 48]), we obtain a subsequence (which we denote with the same
notation) such that

$U^{0}_{\text{N}k}\to U^{0}_{\text{N}}\quad\text{in
${{C}}^{0}(\overline{\Omega})$},\qquad\text{and}\qquad
D^{\beta}U^{0}_{\text{N}k}\to D^{\beta}U^{0}_{\text{N}}\quad\text{in
${{C}}^{0}(\overline{\Omega})$},\quad\text{for all}\ \left|{\beta}\right|=1,\
2,$

which implies that $U^{0}_{\text{N}k}\to U^{0}_{\text{N}}$ in
${{C}}^{2}(\overline{\Omega})$. Consequently, we arrive at the following
limits

$\Delta U^{0}_{\text{N}}=\lim_{k\to\infty}\Delta
U^{0}_{\text{N}k}=0,\qquad\dfrac{\partial{}}{\partial{\nu_{\circ}}}U^{0}_{\text{N}}=\lim_{k\to\infty}\dfrac{\partial{}}{\partial{\nu_{\circ}}}U^{0}_{\text{N}k}=\lim_{k\to\infty}{g_{k}}=0,$

from which we obtain

$\Delta U^{0}_{\text{N}}=0\quad\text{in}\
\Omega,\qquad\dfrac{\partial{U^{0}_{\text{N}}}}{\partial{\nu_{\circ}}}=0\quad\text{on}\in\Sigma,\qquad
U^{0}_{\text{N}}=0\quad\text{on}\ \Gamma.$

This implies that $U^{0}_{\text{N}}\equiv 0$. Comparing this with the first
estimate (73) yields a contradiction because the subsequence
$\\{U^{0}_{\text{N}k}\\}_{k\in\mathbb{N}}$ satisfies

$1=\left|{U^{0}_{\text{N}k}}\right|_{\overline{\Omega}}^{(2+\alpha)}\leqslant
c_{0}(d,\Omega,\alpha)\left[\left|{U^{0}_{\text{N}k}}\right|_{\overline{\Omega}}^{(0)}+\left|{{g_{k}}}\right|_{\Sigma}^{(1+\alpha)}\right]\longrightarrow
0.$

Thus, inequality (45) holds true. Taking the supremum of $\tau<t\in[0,T]$
yields the desired estimate. ∎

##  6 The Nonlinear Problem

In this section, we distinguish the linear part of the main and dynamic
boundary conditions of the nonlinear problem (47). To do this, we compute the
variations of the operators given in (32) with respect to $\tilde{\rho}$. We
introduce the new unknown functions

${V}_{\text{D}}=U_{\text{D}}-U^{0}_{\text{D}}$ and
${V}_{\text{N}}=U_{\text{N}}-U^{0}_{\text{N}}$.

Hence, using Equations (33), we may write (30) in the following form

$\left\\{\begin{aligned}
\pazocal{{L}}_{\tilde{\rho}}{{V}_{\text{D}}}+(\pazocal{{L}}_{\tilde{\rho}}-\Delta)U^{0}_{\text{D}}+(\pazocal{{K}}_{\tilde{\rho}}{{V}_{\text{D}}}+\pazocal{{K}}_{\tilde{\rho}}{U^{0}_{\text{D}}})\pazocal{{L}}_{\tilde{\rho}}{{\tilde{\rho}}}+(\pazocal{{M}}_{\tilde{\rho}}{{V}_{\text{D}}}+\pazocal{{M}}_{\tilde{\rho}}{U^{0}_{\text{D}}})&=0,\quad
y\in\Omega,\ t>0,\\\
\left.{{V}_{\text{D}}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{D}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\pazocal{{L}}_{\tilde{\rho}}{{V}_{\text{N}}}+(\pazocal{{L}}_{\tilde{\rho}}-\Delta)U^{0}_{\text{N}}+(\pazocal{{K}}_{\tilde{\rho}}{{V}_{\text{N}}}+\pazocal{{K}}_{\tilde{\rho}}{U^{0}_{\text{N}}})\pazocal{{L}}_{\tilde{\rho}}{{\tilde{\rho}}}+(\pazocal{{M}}_{\tilde{\rho}}{{V}_{\text{N}}}+\pazocal{{M}}_{\tilde{\rho}}{U^{0}_{\text{N}}})&=0,\quad
y\in\Omega,\ t>0,\\\
\left.{\dfrac{\partial{{V}_{\text{N}}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{N}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\left.{\left(\dfrac{\partial{\tilde{\rho}}}{\partial{t}}+B_{{\tilde{\rho}}}\dfrac{\partial{{{V}_{\text{D}\text{N}}}}}{\partial{\nu_{\circ}}}+B_{{\tilde{\rho}}}\dfrac{\partial{{{U}_{\text{D}\text{N}}^{0}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}&=0,\\\
{\tilde{\rho}}\big{|}_{y\in\Gamma,\,t=0}&=0,\end{aligned}\right.$

(47)

where $\tilde{\rho}=\tilde{\rho}(y,t)=E\rho(y,t)$, $E$ being the operator
(26), satisfies (27).

In order to identify the linear part of the main equations and of the dynamic
boundary condition in (47), we introduce the variation

$\delta\pazocal{{F}}_{0}=\left.\frac{d}{d\lambda}\pazocal{{F}}_{\lambda{\tilde{\rho}}}\right|_{\lambda=0},$

of the operator
$\pazocal{{F}}\in\\{\pazocal{{L}},\pazocal{{K}},\pazocal{{M}}\\}$ which
depends on ${\tilde{\rho}}$. Note that the map
$[\tilde{\rho}\mapsto\pazocal{{F}}_{\tilde{\rho}}]$, where
$\pazocal{{F}}_{\tilde{\rho}}\in\\{\pazocal{{L}}_{\tilde{\rho}},\pazocal{{K}}_{\tilde{\rho}},\pazocal{{M}}_{\tilde{\rho}}\\}$,
are analytic; see (34). Also, we note of the following identities:

$\displaystyle\pazocal{{L}}_{0}{V}_{i}$
$\displaystyle=\Delta{V}_{i},\qquad\pazocal{{M}}_{0}{V}_{i}=0,\qquad\pazocal{{K}}_{0}{V}_{i}=\boldsymbol{\mathbf{N}}\cdot\nabla{V}_{i},\qquad
B_{0}=\frac{1}{\boldsymbol{\mathbf{N}}\cdot\nu_{\circ}},\qquad
i=\text{D},\text{N}.$

Because $\pazocal{{M}}_{0}\equiv 0$ and
$\pazocal{{L}}_{0}{{\tilde{\rho}}}\equiv 0$, then for $i=\text{D},\text{N}$,
we can write the main equations in (47) posed over $\Omega$ in the following
way:

$\displaystyle\pazocal{{L}}_{0}{{V}_{i}}+\delta\pazocal{{L}}_{0}{U_{i}}+\delta\pazocal{{M}}_{0}{U_{i}}+\pazocal{{K}}_{0}{U_{i}}\pazocal{{L}}_{0}{{\tilde{\rho}}}$
$\displaystyle\qquad\qquad=-(\pazocal{{L}}_{{\tilde{\rho}}}-\pazocal{{L}}_{0}){V}_{i}-(\pazocal{{L}}_{{\tilde{\rho}}}-\pazocal{{L}}_{0}-\delta\pazocal{{L}}_{0})U_{i}$
$\displaystyle\qquad\qquad\qquad-(\pazocal{{M}}_{{\tilde{\rho}}}-\pazocal{{M}}_{0}){V}_{i}-(\pazocal{{M}}_{{\tilde{\rho}}}-\pazocal{{M}}_{0}-\delta\pazocal{{M}}_{0})U_{i},$
$\displaystyle\qquad\qquad\qquad-(\pazocal{{K}}_{\tilde{\rho}}{{V}_{i}}+\pazocal{{K}}_{{\tilde{\rho}}}{U_{i}})(\pazocal{{L}}_{{\tilde{\rho}}}{{\tilde{\rho}}}-\pazocal{{L}}_{0}{{\tilde{\rho}}})$
$\displaystyle\qquad\qquad=:\pazocal{{N}}[{V}_{i},{\tilde{\rho}}].$

On the other hand, the dynamic boundary condition on $\Gamma$ in (47) can be
written as follows:

$\displaystyle\dfrac{\partial{\tilde{\rho}}}{\partial{t}}+B_{0}\dfrac{\partial{{{{V}_{\text{D}\text{N}}}}}}{\partial{\nu_{\circ}}}+\left.\delta
B_{0}\dfrac{\partial{{{U}_{\text{D}\text{N}}^{0}}}}{\partial{\nu_{\circ}}}\right|_{\Gamma}$
(48) $\displaystyle\qquad=-\left[(B_{{\tilde{\rho}}}-B_{0}-\delta
B_{0})\dfrac{\partial{{{U}_{\text{D}\text{N}}^{0}}}}{\partial{\nu_{\circ}}}+(B_{{\tilde{\rho}}}-B_{0})\dfrac{\partial{{{{V}_{\text{D}\text{N}}}}}}{\partial{\nu_{\circ}}}+B_{0}\dfrac{\partial{{{U}_{\text{D}\text{N}}^{0}}}}{\partial{\nu_{\circ}}}\right]$
$\displaystyle\qquad=:\pazocal{{B}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}]-B_{0}\dfrac{\partial{{{U}_{\text{D}\text{N}}^{0}}}}{\partial{\nu_{\circ}}}.$

The computations of the variations $\delta\pazocal{{L}}_{0}$,
$\delta\pazocal{{K}}_{0}$, $\delta\pazocal{{M}}_{0}$, and
$\delta\pazocal{{B}}_{0}$, which consists of the variations of $\delta B_{0}$
and/or $\delta\textsf{A}^{(0)}$, can be done without difficulty (see Appendix
D). In fact, it can be checked, for instance, that

$\displaystyle\delta\pazocal{{L}}_{0}U_{i}$
$\displaystyle=\sum_{m,p=1}^{d}\delta{{\textsf{A}}}_{mp}^{(0)}\frac{\partial^{2}U_{i}}{\partial
y_{m}\partial
y_{p}},\qquad(i=\text{D},\text{N}),\qquad\delta{{\textsf{A}}}_{mp}^{(0)}=-\sum_{m=1}^{d}\left(\frac{\partial(N_{m}{\tilde{\rho}})}{\partial
y_{p}}+\frac{\partial(N_{p}{\tilde{\rho}})}{\partial y_{m}}\right),$
$\displaystyle\delta B_{0}$
$\displaystyle=-2\nu_{\circ}\cdot\nabla{\tilde{\rho}}-\frac{\boldsymbol{\mathbf{N}}\cdot\nabla{\tilde{\rho}}}{\boldsymbol{\mathbf{N}}\cdot\nu_{\circ}}+h(\boldsymbol{\mathbf{N}},\partial_{y}\boldsymbol{\mathbf{N}}){\tilde{\rho}},$

where

$h(\boldsymbol{\mathbf{N}},\partial_{y}\boldsymbol{\mathbf{N}}):=-\frac{2\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\nu_{\circ}}{\boldsymbol{\mathbf{N}}\cdot\nu_{\circ}}+\frac{1}{\left(\boldsymbol{\mathbf{N}}\cdot\nu_{\circ}\right)^{2}}\left(\nu_{\circ}\cdot(\nabla\otimes\boldsymbol{\mathbf{N}})^{\top}\boldsymbol{\mathbf{N}}\right).$

Note that $\delta B_{0}$ is well-defined because of Assumption (A2).

Observe from the above expressions that the variations
$\delta\pazocal{{L}}_{0}$ and $\delta B_{0}$ only consist of first-order
derivatives of ${\tilde{\rho}}$. Because ${\tilde{\rho}}$ is satisfies the
Laplace equation (27), the linear part of the first and fourth equation in
problem (47) do not contain the second-order partial derivatives of the
function ${\tilde{\rho}}$. Therefore, we are able to reduce the problem as
follows:

${\left\\{\begin{aligned}
\Delta{V}_{\text{D}}+\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla{\tilde{\rho}}+{Q}_{\text{D}}{\tilde{\rho}}=\pazocal{{N}}[{V}_{\text{D}},{\tilde{\rho}}],&\quad
y\in\Omega,\ t>0,\\\
\left.{{V}_{\text{D}}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{D}}}\right|_{y\in\Gamma}=0,&\\\\[5.0pt]
\Delta{V}_{\text{N}}+\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla{\tilde{\rho}}+{Q}_{\text{N}}{\tilde{\rho}}=\pazocal{{N}}[{V}_{\text{N}},{\tilde{\rho}}],&\quad
y\in\Omega,\ t>0,\\\
\left.{\dfrac{\partial{{V}_{\text{N}}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=0,\qquad\left.{{V}_{\text{N}}}\right|_{y\in\Gamma}=0,&\\\\[5.0pt]
\left.{\left(\dfrac{\partial{\tilde{\rho}}}{\partial{t}}+b_{\perp}\dfrac{\partial{{\tilde{\rho}}}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}{\tilde{\rho}}+h{\tilde{\rho}}+B_{0}\dfrac{\partial{{V}_{\text{D}\text{N}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}&=\pazocal{{B}}[{V}_{\text{D}},{V}_{\text{N}},{\tilde{\rho}}]-B_{0}\dfrac{\partial{{{{U}_{\text{D}\text{N}}^{0}}}}}{\partial{\nu_{\circ}}},\\\
{\tilde{\rho}}\big{|}_{y\in\Gamma,\,t=0}=0,&\end{aligned}\right.}$ (49)

where $h:=h(\boldsymbol{\mathbf{N}},\partial_{y}\boldsymbol{\mathbf{N}})$.
Here, $\pazocal{{B}}$ and $\pazocal{{N}}$ are the nonlinear terms, the
functions $\boldsymbol{\mathbf{q}}_{\text{D}}$ and ${Q}_{\text{D}}$ are
dependent on $U^{0}_{\text{D}}$ and $\boldsymbol{\mathbf{N}}$ and their
derivatives, while the functions $\boldsymbol{\mathbf{q}}_{\text{N}}$ and
${Q}_{\text{N}}$ are dependent on $U^{0}_{\text{N}}$ and
$\boldsymbol{\mathbf{N}}$ and their derivatives. The functions $b_{\perp}$ and
$\boldsymbol{\mathbf{b}}_{\parallel}$ are dependent on
$\boldsymbol{\mathbf{N}}$ and its derivatives.

At this point, it is important to know the sign of the coefficient $b_{\perp}$
of $\dfrac{\partial{{\tilde{\rho}}}}{\partial{\nu_{\circ}}}$ appearing in the
boundary equation on $\Gamma$ (cf. (48)). By examining the explicit form of
$\delta B_{0}$, (A2), and (A3), we deduce that $b_{\perp}<0$.

##  7 The Linear Problem

Our goal here is to prove the existence of classical solution to a linear
problem corresponding to (49). That is, we study the following system of
partial differential equations:

$\left\\{\begin{aligned} \Delta
w_{\text{D}}+\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla\theta+{Q}_{\text{D}}\theta&=F_{\text{D}}(y,t),\quad
y\in\Omega,\ t>0,\\\
\left.{w_{\text{D}}}\right|_{y\in\Sigma}=0,\qquad\left.{w_{\text{D}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\Delta
w_{\text{N}}+\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla\theta+{Q}_{\text{N}}\theta&=F_{\text{N}}(y,t),\quad
y\in\Omega,\ t>0,\\\
\left.{\dfrac{\partial{w_{\text{N}}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=0,\qquad\left.{w_{\text{N}}}\right|_{y\in\Gamma}&=0,\\\\[5.0pt]
\left.{\left(\dfrac{\partial{\theta}}{\partial{t}}+b_{\perp}\dfrac{\partial{\theta}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}\theta+h\theta+B_{0}\dfrac{\partial{{{w}_{\text{D}\text{N}}}}}{\partial{\nu_{\circ}}}\right)}\right|_{y\in\Gamma}&=\psi(y,t),\\\
{\theta}\big{|}_{y\in\Gamma,\,t=0}&=0,\end{aligned}\right.$ (50)

where $\theta=\theta(y,t)$, for $t>0$, is harmonic in $y\in\Omega$ and
vanishes on $\Sigma$ (i.e., $\theta\big{|}_{\Sigma}=0$).

With regards to the above problem, we will prove the following result.

######  Theorem 7.1.

Let $\Sigma,\Gamma\in{{C}}^{2+\alpha}$ for some $\alpha\in(0,1)$, and suppose
that the coefficients in (50) satisfy the following conditions

$\left\\{\begin{aligned} b_{\perp},\
h&\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma)),\\\
\boldsymbol{\mathbf{b}}_{\parallel}&\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma)^{d}),\\\
{Q}_{\text{D}},\
{Q}_{\text{N}}&\in{{C}}^{0}([0,T];{{C}}^{0+\alpha}(\overline{\Omega})),\\\
\boldsymbol{\mathbf{q}}_{\text{D}},\
\boldsymbol{\mathbf{q}}_{\text{N}}&\in{{C}}^{0}([0,T];{{C}}^{0+\alpha}(\overline{\Omega})^{d}),\\\
b_{\perp}&<0.\end{aligned}\right.$

For any
$F_{\text{D}},F_{\text{N}}\in{{C}}^{0}([0,T];{{C}}^{0+\alpha}(\overline{\Omega}))$
and $\psi\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma))$, system (50) has a
unique solution

$w_{\text{D}},\ w_{\text{N}},\
\theta\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$

with $\theta$ having additional smoothness with respect to $t$ on the surface
$\Gamma$. That is, $\theta_{t}\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma))$.
In addition, the following estimate hold:

$\left\|{{w_{{\text{D},\text{N}}}}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{\theta}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant
c\left(\left\|{{F_{{\text{D},\text{N}}}}}\right\|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\psi}\right|^{(1+\alpha)}_{\Gamma;\,[0,t]}\right),$
(51)

for some constant $c>0$.

To solve system (50), we will apply the method of successive approximations
(see, e.g., [44, p. 74] or [55, Sec. 1.1.1, p. 124]). At each step, three
problems are solved: the first two are elliptic equations whose coefficients
and unknown functions depend on time like a parameter (this corresponds to the
pure Dirichlet problem and mixed Dirichlet-Neumann problem in (50)), and the
problem for an elliptic equation with the time-derivative in the boundary
condition (this corresponds to the last two equations in (50)). In the next
several lines, we focus on the existence of solution to the last mentioned
problem.

Let $\Omega\subset\mathbb{R}^{d}$ be an open, bounded, connected set with
boundary $\partial\Omega=\Gamma\cup\Sigma$, where
$\Gamma,\Sigma\in{{C}}^{2+\alpha}$ are disjoint surfaces, and $\Gamma$ is
interior to $\Sigma$. Let us consider and examine the following system:

$\left\\{\begin{aligned} -\Delta\Theta=0,\quad y\in\Omega,\
t>0,\qquad\left.{\Theta}\right|_{y\in\Sigma}=\psi_{2},\qquad{\Theta}\big{|}_{y\in\Gamma,\,t=0}&=0,\\\
\left.{\left(\dfrac{\partial\Theta}{\partial
t}+\boldsymbol{\mathbf{b}}\cdot\nabla\Theta+{\kappa}\Theta\right)}\right|_{y\in\Gamma}&=\psi_{1}.\end{aligned}\right.$
(52)

######  Lemma 7.2.

Let the coefficients in (52) satisfy the following conditions

$b_{i}\ (i=1,\ldots,d),\
{\kappa}\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma)),\qquad\boldsymbol{\mathbf{b}}\cdot\nu<-b_{0}<0,\quad(\boldsymbol{\mathbf{b}}=(b_{1},\ldots,b_{d})^{\top}),$

for some constant $b_{0}>0$. For any given boundary data

$\psi_{1}\in{{C}}^{0}([0,T];{{C}}^{1+\alpha}(\Gamma))$ and
$\psi_{2}\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\Sigma))$,

there exists a unique solution
$\Theta\in{{C}}^{0}([0,T];{{C}}^{2+\alpha}(\overline{\Omega}))$ to problem
(52), and such that the following estimate hold

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left(\left|{\psi_{1}}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\psi_{2}}\right|^{(2+\alpha)}_{{\Sigma;\,[0,t]}}\right),\qquad
t\leqslant T,$ (53)

where the constant $c>0$ depends on the coefficients in the boundary
condition.

###### Proof.

The proof is given in Appendix E.3. ∎

We are now in the position to prove Theorem 7.1.

###### Proof of Theorem 7.1.

We first confirm estimate (51). To this end we write the main equations in
(50) in the form

$\left\\{\begin{aligned} \Delta w_{\text{D}}&={H}_{\text{D}}\equiv
F_{\text{D}}-\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla\theta-{Q}_{\text{D}}\theta,\\\
\Delta w_{\text{N}}&={H}_{\text{N}}\equiv
F_{\text{N}}-\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla\theta-{Q}_{\text{N}}\theta.\end{aligned}\right.$
(54)

We consider the Poisson equation above with the homogeneous Dirichlet
condition on the whole boundary $\partial\Omega$. By the well-known result
about the pure Dirichlet problem (see, e.g., [42, Eq. (1.11), p. 110]) we have
the estimate
$\left|{w_{\text{D}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant
c\left(\left|{{H}_{\text{D}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+{{\max_{\overline{\Omega}}}}\left|{w_{\text{D}}}\right|\right)$,
for some constant $c>0$. We can disregard the term
${{\max_{\overline{\Omega}}}}\left|{w_{\text{D}}}\right|$ because, as we will
verify later in the proof, the Poisson problem has a unique solution (cf. [42,
p. 110]). Hence, we get the bound

$\left|{w_{\text{D}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant
c_{\text{D}}\left(\left|{\theta}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{F_{\text{D}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\right),$
(55)

where $c_{\text{D}}>0$ is a constant that depends only on ${Q}_{\text{D}}$ and
$\boldsymbol{\mathbf{q}}_{\text{D}}$.

For the second equation in (54), we consider the boundary problem for the
Poisson equation with homogenous Neumann condition on $\Sigma$ and a
homogenous Dirichlet condition on $\Gamma$. Using the result from [52, Thm.
3.28(ii), p. 194], we have the estimate

$\left|{w_{\text{N}}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}\leqslant
c\left|{{H}_{\text{N}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\\\
\leqslant
c_{\text{N}}\left(\left|{\theta}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{F_{\text{N}}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\right),$
(56)

for some constant $c,c_{\text{N}}>0$ that depends only on ${Q}_{\text{N}}$ and
$\boldsymbol{\mathbf{q}}_{\text{N}}$.

Next, we consider the last two equations in system (50) and apply estimate
(53) with

$\psi_{1}=\psi-
B_{0}\dfrac{\partial{{{w}_{\text{D}\text{N}}}}}{\partial{\nu_{\circ}}}$

to obtain

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left|{\psi_{1}}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}\\\ \leqslant
c\left(\left|{\psi}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left\|{{w_{{\text{D},\text{N}}}}}\right\|^{(2+\alpha)}_{{\Gamma;\,[0,t]}}\right),$

for some constant $c>0$. Combining this estimate with (55) and (56), we get a
new estimate

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\theta\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}\leqslant
c\left(\left|{\psi}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left|{\theta}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{{F_{{\text{D},\text{N}}}}}\right\|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\right).$

It only remains to estimate the sum
$\left|{\theta}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}$.
To get rid of these terms on the right hand side of the above inequality, we
simply apply the same argument used in the latter part of the proof of Lemma
7.2. That is, we utilize the interpolation inequalities (E.4), apply the
maximum principle (since $\theta$ is harmonic), and then use Grönwall lemma,
noting that $\theta(y,t)=0$ for $y\in\Gamma$ at $t=0$, to eventually get the
desired estimate (51). This ends the verification of estimate (51).

We next establish the solvability of system (50). Our approach is to apply the
method of successive approximation (see, e.g., [56]). To this end, the initial
approximation $(w_{\text{D}}^{(0)},w_{\text{N}}^{(0)},\theta^{(0)})$ is found
by solving the problem

$\mathscr{Q}(w_{\text{D}}^{(0)},w_{\text{N}}^{(0)},\theta^{(0)})(y,t)=(F_{\text{D}},F_{\text{N}})(y,t),\qquad
y\in\overline{\Omega},\quad t>0,$ (57)

where, for $n=0,1,2,\ldots$,

$\begin{gathered}\begin{aligned}
&\mathscr{Q}(w_{\text{D}}^{(n)},w_{\text{N}}^{(n)},\theta^{(n)})(y,t)=(\Psi_{\text{D}},\Psi_{\text{N}})(y,t),\qquad
y\in\overline{\Omega},\quad t>0,\\\\[10.00002pt]
&\qquad\Updownarrow\\\\[10.00002pt] &\left\\{\begin{aligned} &\Delta
w_{\text{D}}^{(n)}=\Psi_{\text{D}}(y,t),\quad y\in\Omega,\
t>0,\quad\quad\left.{w_{\text{D}}^{(n)}}\right|_{y\in\Sigma}=0,\quad\left.{w_{\text{D}}^{(n)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta w_{\text{N}}^{(n)}=\Psi_{\text{D}}(y,t),\quad y\in\Omega,\ t>0,\quad\
\left.{\dfrac{\partial{w_{\text{N}}^{(n)}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=0,\quad\left.{w_{\text{N}}^{(n)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta\theta^{(n)}=0,\quad y\in\Omega,\
t>0,\quad\left.{\theta^{(n)}}\right|_{y\in\Sigma}=0,\qquad\theta^{(n)}\big{|}_{y\in\Gamma,\,t=0}=0,\\\
&\left.{\left(\dfrac{\partial{{\theta}^{(n)}}}{\partial{t}}+b_{\perp}\dfrac{\partial{\theta^{(n)}}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}\theta^{(n)}+h\theta^{(n)}\right)}\right|_{y\in\Gamma}=\psi(y,t)-B_{0}\displaystyle\dfrac{\partial{{w}^{(n)}_{\text{D}\text{N}}}}{\partial{\nu_{\circ}}}.\end{aligned}\right.\end{aligned}\end{gathered}$
(58)

for $n=0,1,\ldots$, where
$(w_{\text{D}}^{(n)},w_{\text{N}}^{(n)},\theta^{(n)})$ are the unknown
triplet, and $(\Psi_{\text{D}},\Psi_{\text{N}})$ are a pair of given
functions.

Then, we define the approximants
$(w_{\text{D}}^{(m+1)},w_{\text{N}}^{(m+1)},\theta^{(m+1)})$, for
$m=0,1,\ldots$, as solution to (58) with $n=m+1$; that is,
$(w_{\text{D}}^{(m+1)},w_{\text{N}}^{(m+1)},\theta^{(m+1)})$ solves the
equation

$\mathscr{Q}(w_{\text{D}}^{(m+1)},w_{\text{N}}^{(m+1)},\theta^{(m+1)})(y,t)=({H}_{\text{D}}^{(m)},{H}_{\text{N}}^{(m)})(y,t),\qquad
y\in\overline{\Omega},\quad t>0,$ (59)

where

${H}_{i}^{(m+1)}\equiv{F}_{i}-\boldsymbol{\mathbf{q}}_{i}\cdot\nabla\theta^{(m)}-{Q}_{i}\theta^{(m)},\qquad
i=\text{D},\text{N}.$

The solvability of the Poisson problem with pure Dirichlet boundary condition
is well-known (see, e.g., [44, Thm. 4.3, p. 56]).666Because $\Omega$ is
bounded, its closure is compact. Meanwhile, because
$w_{\text{N}}^{(m+1)}(y,t)$ vanishes for $y\in\Gamma$, $t>0$, for all
$m=0,1,\ldots$, then the Poisson problem with mixed Dirichlet-Neumann boundary
condition is also guaranteed to be solvable (see [52, Chap. 3]). These
together with Theorem 7.1 and Lemma 7.2 we infer the unique solvability of
(59) for all $t\in[0,T]$, for all $m=0,1,\ldots$. In addition, the solutions
satisfy the following estimates777Note here that
$\left|{\Theta^{(m)}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}={{\max_{\overline{\Omega}}}}\left|{\Theta^{(m)}}\right|+[\Theta^{(m)}]^{(\alpha)}_{{{{\overline{\Omega}}}}}$.

$\displaystyle\left|{w_{i}^{(m+1)}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}$
$\displaystyle\leqslant
c_{i}^{m}\left(\left|{F_{i}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta^{(m)}}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta^{(m)}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}\right),\qquad
i=\text{D},\text{N},$
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\theta^{(m+1)}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}$
$\displaystyle\leqslant
c\left(\left|{\psi}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}+\left\|{{w_{{\text{D},\text{N}}}^{(m+1)}}}\right\|^{(2+\alpha)}_{{\Gamma;\,[0,t]}}\right),$

for some constants
$c_{i}^{m}:=c_{\text{D}}^{m}(Q_{i},\boldsymbol{\mathbf{q}}_{i})>0$,
$i=\text{D},\text{N}$, and $c>0$.

We first estimate the term
$\left|{\Theta^{(m)}}\right|^{(\alpha)}_{{{{\overline{\Omega}}}}}$ via the
interpolation inequality

$[{\Theta^{(m)}}]_{{{{\overline{\Omega}}}}}^{(\alpha)}\leqslant\varepsilon_{4,m}\left|{\Theta^{(m)}}\right|_{{{{\overline{\Omega}}}}}^{(1)}+c_{\varepsilon_{4,m}}\left|{\Theta^{(m)}}\right|_{{{{\overline{\Omega}}}}}^{(0)},$

and choose $\varepsilon_{4,m}>0$ small enough so that we can bound the terms
above only by ${{\max_{\overline{\Omega}}}}\left|{\Theta^{(m)}}\right|$. Then,
we argue as in the latter part of the proof of Lemma 7.2, to write the
estimates above as follows:

$\displaystyle\left|{w_{i}^{(m+1)}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}$
$\displaystyle\leqslant\tilde{c}_{i}^{m}\left(\left|{{F}_{i}}\right|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\theta^{(m)}}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right),\qquad
i=\text{D},\text{N},$
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\theta^{(m+1)}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{(2+\alpha)}_{\overline{\Omega},\,\Gamma;\,[0,t]}$
$\displaystyle\leqslant\tilde{c}\left(K(t)+\left|{\theta^{(m)}}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}}\right),$

for some constants
$\tilde{c}_{i}^{m}:=\tilde{c}_{i}^{m}(Q_{i},\boldsymbol{\mathbf{q}}_{i})>0$,
$i=\text{D},\text{N}$, and $\tilde{c}>0$, where

$K(t):=\left\|{{F_{{\text{D},\text{N}}}}}\right\|^{(\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left|{\psi}\right|^{(1+\alpha)}_{{\Gamma;\,[0,t]}}.$

Moreover, the initial approximant
$(w_{\text{D}}^{(0)},w_{\text{N}}^{(0)},\theta^{(0)})$ satisfies the estimate

$\left\|{w_{{\text{D},\text{N}}}^{(0)}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{\theta^{(0)}}\right\|^{(2+\alpha)}_{\Gamma;\,[0,t]}\leqslant
c_{0}K(t),$ (60)

for some constant $c_{0}>0$.

Let us consider the differences

$\varpi_{i}^{(m+1)}=w_{i}^{(m+1)}-w_{i}^{(m)},\quad
i=\text{D},\text{N},\qquad\text{and}\qquad\vartheta^{(m+1)}=\theta^{(m+1)}-\theta^{(m)},$

where $m=0,1,\ldots$. Clearly, for $m=0,1,\ldots$,
$(\varpi_{\text{D}}^{(m+1)},\varpi_{\text{N}}^{(m+1)},\vartheta^{(m+1)})$
satisfies

$\left\\{\begin{aligned}
&\Delta\varpi_{\text{D}}^{(m+1)}=-\boldsymbol{\mathbf{q}}_{\text{D}}\cdot\nabla\vartheta^{(m)}-{Q}_{\text{D}}\vartheta^{(m)},\quad
y\in\Omega,\
t>0,\quad\left.{\varpi_{\text{D}}^{(m+1)}}\right|_{y\in\Sigma}=0,\quad\left.{\varpi_{\text{D}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta\varpi_{\text{N}}^{(m+1)}=-\boldsymbol{\mathbf{q}}_{\text{N}}\cdot\nabla\vartheta^{(m)}-{Q}_{\text{N}}\vartheta^{(m)},\quad
y\in\Omega,\
t>0,\quad\left.{\dfrac{\partial{\varpi_{\text{N}}^{(m+1)}}}{\partial{\nu_{\circ}}}}\right|_{y\in\Sigma}=0,\quad\left.{\varpi_{\text{N}}^{(m+1)}}\right|_{y\in\Gamma}=0,\\\\[5.0pt]
&\Delta\vartheta^{(m+1)}=0,\quad y\in\Omega,\
t>0,\qquad\left.{\vartheta^{(m+1)}}\right|_{y\in\Sigma}=0,\qquad{\vartheta^{(m+1)}}\Big{|}_{y\in\Gamma,\,t=0}=0,\\\
&\qquad\qquad\qquad\left.{\left(\dfrac{\partial{\vartheta}}{\partial{t}}^{(m+1)}+b_{\perp}\dfrac{\partial{\vartheta^{(m+1)}}}{\partial{\nu_{\circ}}}+\boldsymbol{\mathbf{b}}_{\parallel}\cdot\nabla_{\Gamma}\vartheta^{(m+1)}+h\vartheta^{(m+1)}\right)}\right|_{y\in\Gamma}=-B_{0}\dfrac{\partial{{{\varpi}_{\text{N}\text{D}}^{(m+1)}}}}{\partial{\nu_{\circ}}}.\end{aligned}\right.$

(61)

Following the previous estimations, it can be deduced that
$(\varpi_{\text{D}}^{(m+1)},\varpi_{\text{N}}^{(m+1)},\vartheta^{(m+1)})$
satisfy the estimates

$\left\|{\varpi_{\text{D},\text{N}}^{(m+1)}}\right\|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+\left\|{\vartheta^{(m+1)}}\right\|^{(2+\alpha)}_{{\Gamma;\,[0,t]}}\leqslant
c\left|{\vartheta^{(m)}}\right|^{(1+\alpha)}_{{\overline{\Omega};\,[0,t]}},$
(62)

for some constants $c>0$.

We estimate in the next few lines the right hand side of (62). For this
purpose, we use the following property of norms on ${{C}}^{k+\alpha}$ (see,
e.g., [42, eq. (5.7), p. 403]):

$\left|{u}\right|_{{\overline{\Omega};\,[0,t]}}^{(k-1+\alpha)}\leqslant\varepsilon_{6}\left|{u}\right|_{{\overline{\Omega};\,[0,t]}}^{(k+\alpha)}+c_{\varepsilon_{6}}{{\max_{\overline{\Omega}}}}\left|{u}\right|,\quad(k\geqslant
2).$

Here, ${\varepsilon_{6}}>0$ is an arbitrary small number and
$c_{\varepsilon_{6}}\to\infty$ as $\varepsilon_{6}\to 0$. This, together with
the fact that $\vartheta^{(m)}(y,0)=0$, for $y\in\Gamma$, for each
$m=0,1,\ldots$, leads us to the following estimate

$\displaystyle\left|{\vartheta^{(m)}}\right|_{{\overline{\Omega};\,[0,t]}}^{(1+\alpha)}$
$\displaystyle\leqslant\left(\varepsilon_{6}\left|{\vartheta^{(m)}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+c_{\varepsilon_{6}}{{\max_{\overline{\Omega}}}}\left|{\vartheta^{(m)}}\right|\right)$
(63)
$\displaystyle\leqslant\varepsilon_{6}\left|{\vartheta^{(m)}}\right|^{(2+\alpha)}_{{\overline{\Omega};\,[0,t]}}+c_{\varepsilon_{6}}\int_{0}^{t}{{\max_{\Gamma}}}\left|{\frac{\partial}{\partial
s}}\vartheta^{(m)}(\cdot,s)\right|{{{d}}s}.$

For $m=0$, it is easy to see that