[ beforeskip=.2em plus 1pt,]toclinesection

# Two-phase flows with bulk-surface interaction:
thermodynamically consistent Navier–Stokes–Cahn–Hilliard models
with dynamic boundary conditions

Andrea Giorgini111Research supported by MUR grant Dipartimento di Eccellenza
2023-2027.
Politecnico di Milano
Dipartimento di Matematica
Via E. Bonardi 9, 20133 Milano, Italy
<EMAIL_ADDRESS>Patrik Knopf222Research supported by the DFG
(Deutsche Forschungsgemeinschaft): Project 52469428 and RTG 2339.
Universität Regensburg
Fakultät für Mathematik
Universitätsstr. 31, D-93053 Regensburg
<EMAIL_ADDRESS>
###### Abstract

We derive a novel thermodynamically consistent Navier–Stokes–Cahn–Hilliard
system with dynamic boundary conditions. This model describes the motion of
viscous incompressible binary fluids with different densities. In contrast to
previous models in the literature, our new model allows for surface diffusion,
a variable contact angle between the diffuse interface and the boundary, and
mass transfer between bulk and surface. In particular, this transfer of
material is subject to a mass conservation law including both a bulk and a
surface contribution. The derivation is carried out by means of local energy
dissipation laws and the Lagrange multiplier approach. Next, in the case of
fluids with matched densities, we show the existence of global weak solutions
in two and three dimensions as well as the uniqueness of weak solutions in two
dimensions.
Keywords: Two-phase flows, Navier–Stokes–Cahn–Hilliard system, bulk-surface
interaction, dynamic boundary conditions
Mathematics Subject Classification: 35D30, 35Q35, 76D03, 76D05, 76T06, 76T99.

This is a preprint version of the paper. Please cite as:
A. Giorgini and P. Knopf, J. Math. Fluid Mech. 25:65 (2023),
https://doi.org/10.1007/s00021-023-00811-w

###### Contents

  1. 1 Introduction
  2. 2 Model derivation
    1. 2.1 Considerations based on local mass balance laws
    2. 2.2 Local energy dissipation laws
    3. 2.3 Completion of the model derivation by the Lagrange multiplier approach
  3. 3 Mathematical analysis: the case of matched densities
    1. 3.1 Notation and preliminaries
    2. 3.2 Existence and uniqueness of weak solutions
    3. 3.3 Proofs
  4. A A chain rule formula.

##  1 Introduction

The description of multi-phase flows is a primal topic in modern continuum
fluid dynamics with enormous applications in biology, chemistry and
engineering. Two major approaches have been developed according to the
representation of the interface separating the different components (see [19,
5, 30, 46] and the references therein): Sharp Interface (SI) versus Diffuse
Interface (DI) methods (see Figure 1). In the former, the interface is
represented by a hypersurface in the surrounding domain, which leads to the
formulation of free boundary problems. In the latter, also called phase-field
method, the interface is represented by a layer with finite thickness,
corresponding to the level sets of the concentration function whose evolution
is governed by a macroscopic equation. The main advantage of this approach is
the Eulerian formulation of the phase-field equation, which avoids tracking
the interface as in free boundary problems. Although conceptually opposites,
these approaches are strictly related since a proper scaling of DI systems
approximates SI models in the so-called sharp interface limit (see, e.g.,
[5]).

Figure 1: Difference between the Sharp Interface (left) and Diffuse Interface
(right) methods.

In the context of the DI theory, the cornerstone system for the motion of two
viscous and incompressible fluids with matched (constant) densities is the
model H, which was rigorously derived in [35]. The model H consists of the
following Navier–Stokes–Cahn–Hilliard system

$\displaystyle\rho\partial_{t}{\boldsymbol{v}}+\rho\textnormal{div}\big{(}{\boldsymbol{v}}\otimes{\boldsymbol{v}}\big{)}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla
p=-\varepsilon\,\textnormal{div}\big{(}\nabla\phi\otimes\nabla\phi\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.1a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0$
$\displaystyle\quad\text{in $Q$},$ (1.1b)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.1c)
$\displaystyle\mu=-\varepsilon\Delta\phi+\varepsilon^{-1}F^{\prime}(\phi)$
$\displaystyle\quad\text{in $Q$}.$ (1.1d)

Here, $Q:=\Omega\times(0,T)$, where $\Omega\subset\mathbb{R}^{d}$ with $d=2,3$
denotes a bounded, smooth domain with boundary $\Gamma$ and $T>0$ is a given
final time. In (1.1), $\rho$ is the constant density of the mixture,
${\boldsymbol{v}}$ is the velocity of the mixture, $p$ is the pressure, $\phi$
is the difference of the fluid concentrations, $\mu$ is the chemical
potential. In addition, $\nu$ is the concentration depending viscosity,
$\mathbf{D}$ is the symmetric gradient operator, $\varepsilon$ is a positive
parameter related to the thickness of the interface, $m_{\Omega}$ is the
(bulk) mobility function and $F$ is a double-well potential.

The fundamental assumption in deriving the model H is that the density of the
mixture is constant, meaning that the two fluids have the same constant
density $\rho$ (matched densities). In the past two decades, several works
investigated suitable Navier–Stokes–Cahn–Hilliard generalizations of the model
H aiming to describe both incompressible mixtures with unmatched densities and
compressible two-phase flows. Without claiming completeness, we refer the
reader to [6, 43, 12, 18, 49, 50, 30, 36, 21]. Among them is the
thermodynamically consistent model for viscous incompressible mixtures with
unmatched densities proposed by Abels, Garcke and Grün (AGG model) in the
seminal work [6], which reads as follows:

$\displaystyle\partial_{t}\big{(}\rho(\phi){\boldsymbol{v}}\big{)}+\textnormal{div}\big{(}{\boldsymbol{v}}\otimes(\rho(\phi){\boldsymbol{v}}+\mathbf{J})\big{)}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla
p=-\varepsilon\,\textnormal{div}\big{(}\nabla\phi\otimes\nabla\phi\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.2a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0$
$\displaystyle\quad\text{in $Q$},$ (1.2b)
$\displaystyle\rho(\phi)=\tilde{\rho}_{1}\frac{1+\phi}{2}+\tilde{\rho}_{2}\frac{1-\phi}{2},\quad\mathbf{J}=-\tfrac{1}{2}(\tilde{\rho}_{2}-\tilde{\rho}_{1})m_{\Omega}(\phi)\nabla\mu$
$\displaystyle\quad\text{in $Q$},$ (1.2c)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.2d)
$\displaystyle\mu=-\varepsilon\Delta\phi+\varepsilon^{-1}F^{\prime}(\phi)$
$\displaystyle\quad\text{in $Q$}.$ (1.2e)

In (1.2), the novel (coupling) terms are the concentration depending density
$\rho(\phi)$, where $\tilde{\rho}_{1}$ and $\tilde{\rho}_{2}$ are the
different (constant) densities of the two fluids, and the additional mass flux
$\mathbf{J}$.

In the literature, both the model H and the AGG model have mostly been
supplemented with the classical no-slip/homogeneous Neumann boundary
conditions

$\displaystyle{\boldsymbol{v}}$ $\displaystyle=\mathbf{0}$
$\displaystyle\quad\text{on $\Sigma$},$ (1.3a)
$\displaystyle\partial_{\mathbf{n}}\phi$ $\displaystyle=0$
$\displaystyle\quad\text{on $\Sigma$},$ (1.3b)
$\displaystyle\partial_{\mathbf{n}}\mu$ $\displaystyle=0$
$\displaystyle\quad\text{on $\Sigma$}.$ (1.3c)

Here, $\Sigma:=\Gamma\times(0,T)$ and $\mathbf{n}$ denotes the outward normal
vector field on $\Gamma$. The conditions (1.3a) and (1.3c) ensure the mass
conservation in the bulk, namely $\int_{\Omega}\phi(t)\;\mathrm{d}x$ is a
conserved quantity. Sufficiently regular solutions of the system ((1.2),(1.3))
satisfy the energy dissipation law

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E({\boldsymbol{v}},\phi)=-\int_{\Omega}2\nu(\phi)\,\left|\mathbf{D}{\boldsymbol{v}}\right|^{2}\;\mathrm{d}x-\int_{\Omega}m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\;\mathrm{d}x,$
(1.4)

where the total energy $E$ is given by

$\displaystyle E({\boldsymbol{v}},\phi)$
$\displaystyle:=E_{\text{kin}}(\phi,{\boldsymbol{v}})+E_{\text{bulk}}(\phi),$
(1.5)

where

$\displaystyle
E_{\text{kin}}({\boldsymbol{v}},\phi):=\int_{\Omega}\frac{\rho(\phi)}{2}\left|{\boldsymbol{v}}\right|^{2}\;\mathrm{d}x,\quad
E_{\text{bulk}}(\phi):=\int_{\Omega}\frac{\varepsilon}{2}\left|\nabla\phi\right|^{2}+\frac{1}{\varepsilon}F(\phi)\;\mathrm{d}x$
(1.6)

denote the kinetic energy and the free energy in the bulk, respectively. We
point out that in this model, the fluid concentrations on the surface do not
affect the dynamics in the bulk.

Concerning the mathematical analysis, the Cauchy problems corresponding to
(1.1) and (1.2), each endowed with (1.3), have been studied in [1, 11, 22, 33]
and in [3, 4, 7, 31, 32], respectively. We also refer the interested reader to
[5] (see also [6]) for the connection between the above DI systems and the
two-phase Navier–Stokes free boundary problem.

Even though the AGG model ((1.2),(1.3)) is already very advanced and is
capable of describing the complex case where the two fluids have different
specific densities, it still inherits some limitations from the underlying
(convective) Cahn–Hilliard equation with homogeneous Neumann boundary
conditions. The main limitations are:

  1. (L1)

The boundary condition (1.3b) enforces the diffuse interface separating both
fluids to always intersect the boundary $\Gamma$ at a perfect angle of ninety
degrees. Of course, this restrictive condition will not always be satisfied in
concrete applications as the contact angle of the interface might deviate from
ninety degrees and even change dynamically over the course of time. Moreover,
as discussed in [47], the no-slip boundary condition (1.3a) is not well-suited
for describing general moving contact line phenomena. This is because
situations where the velocity field actually contributes to the motion of the
contact line of the diffuse interface in the bulk at the boundary of the
domain cannot be described.

  2. (L2)

The boundary condition (1.3c) can be regarded as a no-flux boundary condition
as it already implies $\mathbf{J}\cdot\mathbf{n}=0$ on $\Sigma$. Therefore,
the system ((1.2),(1.3)) can merely describe the situation where the mass of
both fluids in the bulk $\Omega$ is conserved but it is not capable of
describing a transfer of material between bulk and surface which could, for
instance, be caused by absorption or adsorption processes as well as chemical
reactions taking place on the boundary (see, e.g., [38]).

Due to these well-known limitations of the Cahn–Hilliard model, and in order
to describe short-range interactions between bulk and surface more precisely,
physicists proposed that the total free energy should contain an additional
contribution on the surface being also of Ginzburg–Landau type (see, e.g.,
[10]):

$\displaystyle
E_{\text{free}}(\phi,\psi):=E_{\text{bulk}}(\phi)+E_{\text{surf}}(\psi)\quad\text{where}\quad
E_{\text{surf}}(\psi):=\int_{\Gamma}\frac{\delta\kappa}{2}\left|\nabla_{\Gamma}\psi\right|^{2}+\frac{1}{\delta}G(\psi)\;\mathrm{d}S,$
(1.7)

and $E_{\text{bulk}}$ is defined as in (1.6). Here, $\psi$ is an additional
phase-field variable describing the material distribution on the surface, the
parameter $\delta>0$ corresponds to the width of the diffuse interface on the
surface and the constant $\kappa\geq 0$ acts as a weight for surface diffusion
effects. In most cases, $\psi$ is just assumed to be the trace of the phase-
field $\phi$ on the boundary, i.e.,

$\displaystyle\phi|_{\Sigma}=\psi\quad\text{on $\Sigma$}.$ (1.8)

However, also more general relations between $\phi$ and $\psi$ (so-called
transmission conditions) have been investigated in the literature (see, e.g.,
[16, 39]). The function $G$ stands for an additional potential on the surface.
If phase separating processes on the boundary are to be described, $G$ is also
chosen to be double-well shaped.

Based on this free energy, several dynamic boundary conditions for the
Cahn–Hilliard (CH) equation

$\displaystyle\partial_{t}\phi=\textnormal{div}\big{(}m_{\Omega}(\psi)\nabla\mu\big{)}$
$\displaystyle\;\;\text{in $Q$},$ (1.9a)
$\displaystyle\mu=-\varepsilon\Delta\phi+\delta^{-1}F^{\prime}(\phi)$ in $Q$
(1.9b)

have been introduced in the literature. For instance, the Allen–Cahn type
dynamic boundary condition

$\displaystyle\delta\partial_{t}\psi=\kappa\delta\Delta_{\Gamma}\psi-\varepsilon\partial_{\mathbf{n}}\phi-\frac{1}{\delta}G^{\prime}(\psi)\quad\text{on
$\Sigma$},$ (1.10)

to replace (1.3b) was proposed in [37] and further analyzed in [15]. Here, the
symbol $\Delta_{\Gamma}$ denotes the Laplace–Beltrami operator (see, e.g.,
[20] for differential operators on surfaces).

In recent times, dynamic boundary conditions which also exhibit a
Cahn–Hilliard type structure have become very popular, especially since they
allow for a better description of the transfer of material between bulk and
surface (see, e.g., [38]). In these models, the boundary condition (1.3b) is
replaced by a Cahn–Hilliard type equation on the surface, that is

$\displaystyle\partial_{t}\psi=\textnormal{div}_{\Gamma}\big{(}m_{\Gamma}(\psi)\nabla_{\Gamma}\theta\big{)}-\beta
m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu$ $\displaystyle\quad\text{on
$\Sigma$},$ (1.11a)
$\displaystyle\theta=-\delta\kappa\Delta_{\Gamma}\psi+\delta^{-1}G^{\prime}(\psi)+\varepsilon\partial_{\mathbf{n}}\phi$
$\displaystyle\quad\text{on $\Sigma$}.$ (1.11b) Here, $\theta$ denotes the
chemical potential on the boundary, the function $m_{\Gamma}$ describes the
mobility and the parameter $\kappa\geq 0$ acts as a weight for surface
diffusion effects. Furthermore, the boundary condition (1.3c) needs to be
replaced by a suitable coupling condition for the bulk chemical potential
$\mu$ and the surface chemical potential $\theta$. Recently, for a parameter
$\beta\neq 0$, the following choices have been considered:
$\displaystyle\begin{cases}Lm_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=\beta\theta-\mu&\text{if
$L\in[0,\infty)$},\\\ m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=0&\text{if
$L=\infty$}.\end{cases}$ (1.11c)

  * •

If $L=0$, (1.11c) is to be interpreted as the Dirichlet type coupling
condition $\beta\theta=\mu$ on $\Sigma$. This means that the chemical
potentials $\mu$ and $\theta$ are assumed to always remain in chemical
equilibrium. In this case, the dynamic boundary conditions (1.11) were
proposed in [24] (for $\kappa=0$) and in [34] (for general $\kappa\geq 0$).
Sufficiently regular solutions of the system ((1.8),(1.9),(1.11)) satisfy the
mass conservation law

$\displaystyle\beta\int_{\Omega}\phi(t)\;\mathrm{d}x+\int_{\Gamma}\psi(t)\;\mathrm{d}S=\beta\int_{\Omega}\phi_{0}\;\mathrm{d}x+\int_{\Gamma}\psi_{0}\;\mathrm{d}S$
(1.12)

and the energy dissipation law

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E_{\text{free}}\big{(}\phi(t),\psi(t)\big{)}=-\int_{\Omega}m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\;\mathrm{d}x-\int_{\Gamma}m_{\Gamma}(\psi)\left|\nabla_{\Gamma}\theta\right|^{2}\;\mathrm{d}S$
(1.13)

for all $t\in[0,T]$. In the case $\beta>0$, the well-posedness of the system
((1.8),(1.9),(1.11)) was studied in [34] and its long-time behavior was
investigated in [34, 29].

  * •

For $L=\infty$, the dynamic boundary conditions (1.11) were derived in [42] by
means of an energetic variational approach based on Onsager’s law and the
least action principle. Here, due to (1.11c), the chemical potentials are not
directly coupled. However, mechanical interactions of the materials between
bulk and surface are still taken into account by the trace relation (1.8).
Since (1.11c) with $L=\infty$ implies that the mass flux between bulk and
surface is zero, the bulk mass and the surface mass are conserved separately.
To be precise, sufficiently regular solutions of the system
((1.8),(1.9),(1.11)) satisfy the mass conservation law

$\displaystyle\int_{\Omega}\phi(t)\;\mathrm{d}x=\int_{\Omega}\phi_{0}\;\mathrm{d}x\quad\text{and}\quad\int_{\Gamma}\psi(t)\;\mathrm{d}S=\int_{\Gamma}\psi_{0}\;\mathrm{d}S$
(1.14)

and the energy dissipation law (1.13) for all $t\in[0,T]$. The well-posedness
of the system ((1.8),(1.9),(1.11)) was established in [42, 28] and its long-
time behavior was studied in [42, 45].

  * •

In the case $L\in(0,\infty)$, the boundary conditions (1.11) were proposed and
analyzed in [38]. Here, the chemical potentials $\mu$ and $\theta$ are coupled
by a Robin type boundary condition. The constant $1/L$ is related to the
kinetic rate associated with adsorption/desorption processes or chemical
reactions on the boundary. Sufficiently regular solutions of the system
((1.8),(1.9),(1.11)) satisfy the mass conservation law (1.12) and the energy
dissipation law

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}E_{\text{free}}\big{(}\phi(t),\psi(t)\big{)}&=-\int_{\Omega}m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\;\mathrm{d}x-\int_{\Gamma}m_{\Gamma}(\psi)\left|\nabla_{\Gamma}\theta\right|^{2}\;\mathrm{d}S\\\
&\qquad-\frac{1}{L}\int_{\Gamma}(\beta\theta-\mu)^{2}\;\mathrm{d}S\end{split}$
(1.15)

for all $t\in[0,T]$. The weak and strong well-posedness of the system
((1.8),(1.9),(1.11)) was established in [38]. It was further shown that the
case $L\in(0,\infty)$ can be understood as an interpolation between $L=0$ and
$L=\infty$ as these cases are obtained as asymptotic limits on the level of
strong solutions to the system ((1.8),(1.9),(1.11)) as the parameter $L$ is
sent to zero or to infinity, respectively. For the investigation of long-time
dynamics, we refer to [29].

We further point out that a nonlocal variant of the system
((1.8),(1.9),(1.11)) was proposed and investigated in [41]. Recent reviews of
results concerning the Cahn–Hilliard equation with dynamic boundary conditions
can be found in [44, 54].

In the context of Navier–Stokes–Cahn–Hilliard models for two-phase flows, in
order to overcome the aforementioned limitation (L1), the authors in [47]
proposed a variational derivation through Onsager’s principle of maximal
energy dissipation of a new class of generalized Navier slip boundary
conditions for two-phase flows

$\displaystyle{\boldsymbol{v}}\cdot\mathbf{n}=0,\quad
2\nu(\psi)\left(\mathbf{D}{\boldsymbol{v}}\cdot\mathbf{n}\right)_{\tau}+\gamma{\boldsymbol{v}}_{\tau}=L(\psi)\nabla_{\Gamma}\psi$
$\displaystyle\quad\text{on $\Sigma$},$ (1.16a)
$\displaystyle\partial_{t}\psi+{\boldsymbol{v}}_{\tau}\cdot\nabla_{\Gamma}\psi=-\beta
L(\psi),\quad$ $\displaystyle\quad\text{on $\Sigma$},$ (1.16b)
$\displaystyle\partial_{\mathbf{n}}\mu=0$ $\displaystyle\quad\text{on
$\Sigma$},$ (1.16c)

where $\psi=\phi|_{\Sigma}$ as in (1.8), and

$L(\psi)=-\kappa\Delta_{\Gamma}\psi+\varepsilon\partial_{\mathbf{n}}\psi+G^{\prime}(\psi).$

Here, the subscript $\tau$ denotes the tangential component of a vector
${\boldsymbol{w}}$, i.e.,
${\boldsymbol{w}}_{\tau}={\boldsymbol{w}}-({\boldsymbol{w}}\cdot\mathbf{n})\mathbf{n}$,
$\nabla_{\Gamma}$ is the tangential gradient on $\Gamma$ and $\Delta_{\Gamma}$
denotes the Laplace–Beltrami operator on $\Gamma$ (see, e.g., [20] for
differential operators on surfaces). The phenomenological parameters $\beta$,
$\kappa$ and $\gamma$ are positive, and $G$ represents a surface double-well
potential. The surface dynamic of the concentration $\psi$ is described by
(1.16b) that is a convective Allen–Cahn type equation on the surface. It can
be understood as an adaptation of the dynamic boundary condition (1.10) (with
$\delta=1)$ to the situation of an additional volume-averaged velocity field.
Notice that the boundary condition (1.16c) (together with either (1.1) or
(1.2)) entails that $\int_{\Omega}\phi(t)\;\mathrm{d}x$ is invariant (as for
(1.3c)), whereas $\int_{\Gamma}\psi(t)\;\mathrm{d}\Gamma$ is not a conserved
quantity. The model H (1.1) endowed with (1.16) was first studied in [25]. The
authors proved the existence of global weak solutions assuming that both $F$
and $G$ are polynomial-like functions or $F$ is a singular potential and $G$
is polynomial function. In the former case, they also showed the convergence
of any weak solution to a stationary state. Later on, the global existence of
weak solutions for (1.2) endowed with (1.16) has been achieved in [26]. More
recently, the compressible Navier–Stokes–Cahn–Hilliard system supplemented
with (1.16) has been studied in [14], and the Navier–Stokes–Allen–Cahn and the
Navier–Stokes–Voigt–Allen–Cahn systems with (1.16) have been analyzed in [23].
For the stochastic model H with dynamic boundary conditions (1.16), we mention
the recent results in [27]. It is worth pointing out that none of the
aforementioned works established any uniqueness result for weak solutions in
either two or three dimensions. In fact, the only uniqueness result
established so far concerns quasi-strong solutions for the
Navier–Stokes–Voigt–Allen–Cahn systems with (1.16) in [23].

In this paper, we derive and study a new class of thermodynamically consistent
Navier–Stokes–Cahn–Hilliard systems with dynamic boundary conditions. Our
model derivation is based on local mass balance laws in the bulk and on the
surface in which the mass fluxes are still to be determined. Arguing similarly
as in [6], after considering local energy dissipation laws, we apply the
Lagrange multiplier approach to complete our model derivation by identifying
the unknown mass fluxes on the macroscopic level. The resulting
Navier–Stokes–Cahn–Hilliard (NSCH) model reads as follows:

$\displaystyle\partial_{t}\big{(}\rho(\phi){\boldsymbol{v}}\big{)}+\textnormal{div}\big{(}{\boldsymbol{v}}\otimes(\rho(\phi){\boldsymbol{v}}+\mathbf{J})\big{)}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla
p=-\varepsilon\,\textnormal{div}\big{(}\nabla\phi\otimes\nabla\phi\big{)}$
$\displaystyle\;\;\text{in $Q$},$ (1.17a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0,$
$\displaystyle\;\;\text{in $Q$},$ (1.17b)
$\displaystyle\rho(\phi)=\tilde{\rho}_{1}\frac{1+\phi}{2}+\tilde{\rho}_{2}\frac{1-\phi}{2},\quad\mathbf{J}=-\tfrac{1}{2}(\tilde{\rho}_{2}-\tilde{\rho}_{1})m_{\Omega}(\phi)\nabla\mu$
$\displaystyle\;\;\text{in $Q$},$ (1.17c)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\;\;\text{in $Q$},$ (1.17d)
$\displaystyle\mu=-\varepsilon\Delta\phi+\varepsilon^{-1}F^{\prime}(\phi)$
$\displaystyle\;\;\text{in $Q$},$ (1.17e)
$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}}_{\tau})=\textnormal{div}_{\Gamma}\big{(}m_{\Gamma}(\psi)\nabla_{\Gamma}\theta\big{)}-\beta
m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu$ $\displaystyle\;\;\text{on
$\Sigma$},$ (1.17f)
$\displaystyle\theta=-\delta\kappa\Delta_{\Gamma}\psi+\delta^{-1}G^{\prime}(\psi)+\varepsilon\partial_{\mathbf{n}}\phi$
$\displaystyle\;\;\text{on $\Sigma$},$ (1.17g)
$\displaystyle{\boldsymbol{v}}\cdot\mathbf{n}=0,\qquad\phi|_{\Gamma}=\psi,\qquad\begin{cases}Lm_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=\beta\theta-\mu&\text{if
$L\in[0,\infty)$},\\\ m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=0&\text{if
$L=\infty$}\end{cases}$ $\displaystyle\;\;\text{on $\Sigma$},$ (1.17h)
$\displaystyle\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\gamma(\psi){\boldsymbol{v}}\big{]}_{\tau}=\big{[}-\psi\nabla_{\Gamma}\theta+\tfrac{1}{2}\,(\mathbf{J}\cdot\mathbf{n})\,{\boldsymbol{v}}\big{]}_{\tau}$
$\displaystyle\;\;\text{on $\Sigma$},$ (1.17i)
$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}$
$\displaystyle\;\;\text{in $\Omega$},$ (1.17j)
$\displaystyle\psi|_{t=0}=\psi_{0}$ $\displaystyle\;\;\text{on $\Gamma$},$
(1.17k)

In this model, $m_{\Omega}$ and $m_{\Gamma}$ are nonnegative functions
respresenting the mobilities. The positive parameters $\varepsilon$ and
$\delta$ are related to the thickness of the diffuse interfaces in the bulk
and on the boundary, respectively. The constant $\kappa\geq 0$ acts as a
weight for surface diffusion effects. Moreover, $\beta\in\mathbb{R}$ satisfies
$\beta\left|\Omega\right|+\left|\Gamma\right|\neq 0$, and $L\in[0,\infty]$. In
contrast to the dynamic boundary conditions (1.16) which are of Allen–Cahn
type, the dynamic boundary conditions in system (1.17) are of Cahn–Hilliard
type. More precisely, they can be understood as a convective variant of the
dynamic boundary conditions $\big{(}\eqref{BC:TRACE},\eqref{DBC:CH}\big{)}$
presented above. The sharp interface limit corresponding to the system (1.17)
is to be studied in a subsequent work.

The system (1.17) is associated with the total energy

$\displaystyle E({\boldsymbol{v}},\phi,\psi)$
$\displaystyle:=E_{\text{kin}}(\phi,{\boldsymbol{v}})+E_{\text{free}}(\phi,\psi)$
(1.18)

where $E_{\text{kin}}$ is the kinetic energy introduced in (1.6) and
$E_{\text{free}}=E_{\text{bulk}}+E_{\text{surf}}$ is the free energy defined
in (1.7). Sufficiently regular solutions of system (1.17) satisfy the mass
conservation law

$\displaystyle\left\\{\begin{aligned}
&\beta\int_{\Omega}\phi(t)\;\mathrm{d}x+\int_{\Gamma}\psi(t)\;\mathrm{d}S=\beta\int_{\Omega}\phi_{0}\;\mathrm{d}x+\int_{\Gamma}\psi_{0}\;\mathrm{d}S&&\quad\text{if
$L\in[0,\infty)$,}\\\
&\int_{\Omega}\phi(t)\;\mathrm{d}x=\int_{\Omega}\phi_{0}\;\mathrm{d}x\quad\text{and}\quad\int_{\Gamma}\psi(t)\;\mathrm{d}S=\int_{\Gamma}\psi_{0}\;\mathrm{d}S&&\quad\text{if
$L=\infty$,}\end{aligned}\right.$ (1.19)

and the energy dissipation law

$\displaystyle\begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}t}E({\boldsymbol{v}},\phi,\psi)&=-\int_{\Omega}2\nu(\phi)\,\left|\mathbf{D}{\boldsymbol{v}}\right|^{2}\;\mathrm{d}x-\int_{\Gamma}\gamma(\psi)\left|{\boldsymbol{v}}\right|^{2}\;\mathrm{d}S-\int_{\Omega}m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\;\mathrm{d}x\\\
&\qquad-\int_{\Gamma}m_{\Gamma}(\psi)\left|\nabla_{\Gamma}\theta\right|^{2}\;\mathrm{d}S-\mathds{h}(L)\int_{\Gamma}(\beta\theta-\mu)^{2}\;\mathrm{d}S,\end{aligned}$
(1.20)

where

$\displaystyle\mathds{h}(L)=\begin{cases}0&\text{if $L=0$ or $L=\infty$},\\\
L^{-1}&\text{if $L\in(0,\infty)$},\end{cases}$ (1.21)

as long as the solution exists.

After presenting the model derivation of the system (1.17) in Section 2, we
will study the case of fluids having matched densities, namely
$\rho\equiv\tilde{\rho}_{1}=\tilde{\rho}_{2}$. This means that the flux
$\mathbf{J}$ vanishes. We rewrite the system (1.17) as follows

$\displaystyle\rho\,\partial_{t}{\boldsymbol{v}}+\rho\,\textnormal{div}\big{(}{\boldsymbol{v}}\otimes{\boldsymbol{v}}\big{)}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla
p=-\varepsilon\,\textnormal{div}\big{(}\nabla\phi\otimes\nabla\phi\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.22a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0$
$\displaystyle\quad\text{in $Q$},$ (1.22b)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\quad\text{in $Q$},$ (1.22c)
$\displaystyle\mu=-\varepsilon\Delta\phi+\varepsilon^{-1}F^{\prime}(\phi)$
$\displaystyle\quad\text{in $Q$},$ (1.22d)
$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}}_{\tau})=\textnormal{div}_{\Gamma}\big{(}m_{\Gamma}(\psi)\nabla_{\Gamma}\theta\big{)}-\beta
m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu$ $\displaystyle\quad\text{on
$\Sigma$},$ (1.22e)
$\displaystyle\theta=-\delta\kappa\Delta_{\Gamma}\psi+\delta^{-1}G^{\prime}(\psi)+\varepsilon\partial_{\mathbf{n}}\phi$
$\displaystyle\quad\text{on $\Sigma$},$ (1.22f)
$\displaystyle{\boldsymbol{v}}\cdot\mathbf{n}=0,\qquad\phi|_{\Gamma}=\psi,\qquad\begin{cases}Lm_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=\beta\theta-\mu&\text{if
$L\in[0,\infty)$},\\\ m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu=0&\text{if
$L=\infty$}\end{cases}$ $\displaystyle\quad\text{on $\Sigma$},$ (1.22g)
$\displaystyle\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\gamma(\psi){\boldsymbol{v}}\big{]}_{\tau}=\big{[}-\psi\nabla_{\Gamma}\theta\big{]}_{\tau}$
$\displaystyle\quad\text{on $\Sigma$},$ (1.22h)
$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}$
$\displaystyle\quad\text{in $\Omega$},$ (1.22i)
$\displaystyle\psi|_{t=0}=\psi_{0}$ $\displaystyle\quad\text{on $\Gamma$}.$
(1.22j)

The model (1.22) can be regarded as a generalization of the model H with
dynamic boundary conditions of Cahn–Hilliard type. We will focus on the
mathematical analysis of system (1.22) in the case $L=0$ and $\beta>0$ (see
system (3.12) in Section 3). In Theorem 3.3, we first prove the existence of
global weak solutions to system (3.12) with regular (polynomial-like)
potentials $F$ and $G$ in bounded Lipschitz domains (cf.
($\mathrm{A2}$)–($\mathrm{A5}$) below for the specific assumptions). Weak
solutions satisfy the variational formulation of the problem (3.12) and a weak
energy dissipation law (see (3.18)). A weak solution is obtained as the limit
of smooth solutions to suitably regularized problems based on a Faedo–Galerkin
scheme, where appropriate compactness properties are derived on the basis of
uniform and global estimates gained from the energy dissipation law (1.20).
However, an ad hoc construction of the solution ansatz is required since both
$(\phi,\psi)$ and $(\mu,\theta)$ need to be approximated in the same finite-
dimensional subspace in order to show the validity of (1.20) for such
regularized solutions (cf. the change of variables (3.22)). In comparison with
the works in literature [25, 26, 23], our argument is rather simple and relies
only on a one-level Faedo–Galerkin approximation. Furthermore, it is shown
that $(\mu,\theta)\in L^{4}(0,T;\mathcal{L}^{2})$, which is a novel regularity
for NSCH systems (and even CH equations) with regular potentials (cf., for
instance, [11, 22]). Next, in Theorem 3.4, we then establish the uniqueness of
the weak solutions in the two dimensional case. This is the first uniqueness
result in the literature for a NSCH system with dynamic boundary conditions.
It is important to underline that the nonlinear coupling in the NSCH system
with dynamic boundary conditions (1.22) is much stronger than in the model H
(1.1). Roughly speaking, this is due to the elastic (surface) term
$[\psi\nabla\theta]_{\tau}$ in (1.22h), in addition to the classical
capillarity (bulk) term
$-\varepsilon\textnormal{div}(\nabla\phi\otimes\nabla\phi)$ (also referred to
as Korteweg force). More precisely, gradient terms arising from equivalent
formulations of the Korteweg force (cf. (3.11)) disappear in bulk integrals in
view of the incompressibility constraint
$\textnormal{div}\,{\boldsymbol{v}}=0$. Instead, surface integrals involving
the product between a gradient function and ${\boldsymbol{v}}$ does not vanish
since $\textnormal{div}_{\tau}\,{\boldsymbol{v}}\neq 0$ on $\Sigma$. For this
reason, the term $[\psi\nabla\theta]_{\tau}$ cannot be handled as in [33]
through dual estimates. Therefore, the appropriate functional to control the
difference of two solutions corresponds to the total energy. This functional
is capable of reproducing the same cancellation of the nonlinear terms as in
the derivation of (1.20). In order to rigorously justify this fact, two
additional assumptions need to be imposed: the domain $\Omega$ has a
$C^{3}$-boundary, which is needed to gain further spatial regularity of
$(\phi,\psi)$ and to recover the relations of the bulk/surface chemical
potentials $(\mu,\theta)$ almost everywhere; the relation between bulk and
surface potentials $F(s)=\beta G(s)$ for all $s\in\mathbb{R}$, which enables
us to apply a particular chain formula (see Appendix A).

##  2 Model derivation

In this section we will derive the novel Navier–Stokes–Cahn–Hilliard system
with dynamic boundary conditions (1.17).

###  2.1 Considerations based on local mass balance laws

We consider the time evolution of two fluids (indexed with $i=1,2$) in a
(sufficiently smooth) domain $\Omega\subset\mathbb{R}^{d}$ with
$d\in\\{2,3\\}$ on a time interval $[0,T]$ with $T>0$. We write
$Q=\Omega\times(0,T)$ and $\Sigma=\Gamma\times(0,T)$, where
$\Gamma=\partial\Omega$. Let $\rho_{i}:Q\cup\Sigma\to\mathbb{R}$, $i=1,2$,
denote the mass densities of these two fluids. The local mass balance
equations in the bulk read as

$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}\,\widehat{\mathbf{J}}_{i}=0,\quad
i=1,2,\quad\text{in $Q$.}$ (2.1)

Here, $\widehat{\mathbf{J}}_{i}$ are the mass fluxes. As we want to allow for
a transfer of material between bulk and surface, the local mass balance
equation on the boundary is given by

$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}_{\Gamma}\,\widehat{\mathbf{K}}_{i}=\beta\widehat{\mathbf{J}}_{i}\cdot\mathbf{n},\quad
i=1,2,\quad\text{on $\Sigma$.}$ (2.2)

In this relation, $\widehat{\mathbf{K}}_{i}$ stand for the mass fluxes on the
boundary. The right-hand side $\beta\widehat{\mathbf{J}}_{i}$ describes the
transfer of mass between bulk and surface. Here, the constant
$\beta\in\mathbb{R}$ acts as a weight of this mass transfer. Introducing the
individual velocities ${\boldsymbol{v}}_{i}:Q\cup\Sigma\to\mathbb{R}^{d}$,
$i=1,2$, which can be expressed as

$\displaystyle{\boldsymbol{v}}_{i}=\frac{\widehat{\mathbf{J}}_{i}}{\rho_{i}}\quad\text{in
$Q$},\qquad{\boldsymbol{v}}_{i}=\frac{\widehat{\mathbf{K}}_{i}}{\rho_{i}}\quad\text{on
$\Sigma$},\quad i=1,2,$

the mass balances (2.1) and (2.2) can be rewritten as

$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}(\rho_{i}{\boldsymbol{v}}_{i})$
$\displaystyle=0,$ $\displaystyle\quad i=1,2,\quad\text{in $Q$,}$ (2.3)
$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}_{\Gamma}(\rho_{i}{\boldsymbol{v}}_{i})$
$\displaystyle=\beta\widehat{\mathbf{J}}_{i}\cdot\mathbf{n},$
$\displaystyle\quad i=1,2,\quad\text{on $\Sigma$.}$ (2.4)

Let now $\phi_{i}:Q\to\mathbb{R}$ and $\psi_{i}:\Sigma\to\mathbb{R}$, $i=1,2$,
denote the volume fractions of the two fluids in the bulk and on the boundary,
respectively. Provided that each fluid has a constant density
$\tilde{\rho}_{i}$, $i=1,2$, the volume fractions can be identified as

$\displaystyle\phi_{i}=\frac{\rho_{i}}{\tilde{\rho}_{i}}\quad\text{in
$Q$,}\qquad\psi_{i}=\frac{\rho_{i}}{\tilde{\rho}_{i}}\quad\text{on
$\Sigma$,}\quad i=1,2.$ (2.5)

In particular, this naturally entails the trace relation

$\displaystyle\phi|_{\Sigma}=\psi\quad\text{on $\Sigma$.}$ (2.6)

Under the assumption that the excess volume is zero, we have

$\displaystyle\phi_{1}+\phi_{2}=1\quad\text{in
$Q$,}\qquad\psi_{1}+\psi_{2}=1\quad\text{on $\Sigma$.}$ (2.7)

We further set

$\displaystyle\phi=\phi_{2}-\phi_{1}\quad\text{in
$Q$,}\qquad\psi=\psi_{2}-\psi_{1}\quad\text{on $\Sigma$.}$ (2.8)

Furthermore, let ${\boldsymbol{v}}:Q\cup\Sigma\to\mathbb{R}^{d}$ denote the
volume averaged velocity field associated with the two fluids. It can be
expressed as

$\displaystyle{\boldsymbol{v}}$
$\displaystyle=\phi_{1}{\boldsymbol{v}}_{1}+\phi_{2}{\boldsymbol{v}}_{2},$
$\displaystyle\quad i=1,2,\quad\text{in $Q$,}$ (2.9)
$\displaystyle{\boldsymbol{v}}$
$\displaystyle=\psi_{1}{\boldsymbol{v}}_{1}+\psi_{2}{\boldsymbol{v}}_{2},$
$\displaystyle\quad i=1,2,\quad\text{on $\Sigma$.}$ (2.10)

We assume that both fluids cannot permeate the boundary $\Gamma$, which leads
to the condition

$\displaystyle{\boldsymbol{v}}_{\mathbf{n}}={\boldsymbol{v}}\cdot\mathbf{n}=0\quad\text{on
$\Sigma$.}$ (2.11)

Hence, ${\boldsymbol{v}}$ is identical to its tangential component, that is,
${\boldsymbol{v}}_{\tau}={\boldsymbol{v}}-({\boldsymbol{v}}\cdot\mathbf{n})\mathbf{n}={\boldsymbol{v}}$
on $\Sigma$. We further write

$\displaystyle\mathbf{J}_{i}$
$\displaystyle=\widehat{\mathbf{J}}_{i}-\rho_{i}{\boldsymbol{v}},$
$\displaystyle\quad i=1,2,\quad\text{in $Q$,}$ (2.12)
$\displaystyle\mathbf{K}_{i}$
$\displaystyle=\widehat{\mathbf{K}}_{i}-\rho_{i}{\boldsymbol{v}},$
$\displaystyle\quad i=1,2,\quad\text{on $\Sigma$,}$ (2.13)

to denote the mass fluxes relative to the volume averaged velocity. The mass
balance equations (2.3) and (2.4) can thus be expressed as

$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}(\rho_{i}{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{i}$
$\displaystyle=0,$ $\displaystyle\quad i=1,2,\quad\text{in $Q$,}$ (2.14)
$\displaystyle\partial_{t}\rho_{i}+\textnormal{div}_{\Gamma}(\rho_{i}{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{i}$
$\displaystyle=\beta\mathbf{J}_{i}\cdot\mathbf{n},$ $\displaystyle\quad
i=1,2,\quad\text{on $\Sigma$.}$ (2.15)

In this context, $\mathbf{J}_{i}$ and $\mathbf{K}_{i}$, $i=1,2$, can be
regarded as diffusive flow rates. We now define the total mass density as

$\displaystyle\rho$ $\displaystyle=\rho_{1}+\rho_{2},$ $\displaystyle\quad
i=1,2,\quad\text{in $Q\cup\Sigma$.}$ (2.16)

This leads to the relations

$\displaystyle\partial_{t}\rho+\textnormal{div}(\rho{\boldsymbol{v}})+\textnormal{div}(\mathbf{J}_{1}+\mathbf{J}_{2})$
$\displaystyle=0$ in $Q$, (2.17)
$\displaystyle\partial_{t}\rho+\textnormal{div}_{\Gamma}(\rho{\boldsymbol{v}})+\textnormal{div}_{\Gamma}(\mathbf{K}_{1}+\mathbf{K}_{2})$
$\displaystyle=\beta(\mathbf{J}_{1}+\mathbf{J}_{2})\cdot\mathbf{n}$ on
$\Sigma$. (2.18)

Next, assuming a conservation of linear momentum with respect to the velocity
field ${\boldsymbol{v}}$, we have

$\displaystyle\partial_{t}(\rho{\boldsymbol{v}})+\textnormal{div}(\rho{\boldsymbol{v}}\otimes{\boldsymbol{v}})=\textnormal{div}\,\tilde{\mathbf{T}}\quad\text{in
$Q$.}$ (2.19)

Here, $\tilde{\mathbf{T}}$ denotes the stress tensor that needs to be
specified by constitutive assumptions. In the following, we write

$\displaystyle\mathbf{J}=\mathbf{J}_{1}+\mathbf{J}_{2}\quad\text{in
$Q$,}\qquad\mathbf{K}=\mathbf{K}_{1}+\mathbf{K}_{2}\quad\text{on $\Sigma$.}$
(2.20)

Using (2.17) to rewrite (2.19), we obtain

$\displaystyle\begin{aligned}
\rho(\partial_{t}{\boldsymbol{v}}+{\boldsymbol{v}}\cdot\nabla{\boldsymbol{v}})&=\textnormal{div}\,\tilde{\mathbf{T}}+(\textnormal{div}\,\mathbf{J}){\boldsymbol{v}}\\\
&=\textnormal{div}(\tilde{\mathbf{T}}+{\boldsymbol{v}}\otimes\mathbf{J})-\nabla{\boldsymbol{v}}\,\mathbf{J}\quad\text{in
$Q$.}\end{aligned}$ (2.21)

This allows us to define the objective tensor

$\displaystyle\mathbf{T}=\tilde{\mathbf{T}}+{\boldsymbol{v}}\otimes\mathbf{J}\quad\text{in
$Q$,}$ (2.22)

which is frame indifferent. Hence, (2.21) can be rewritten as

$\displaystyle\rho\partial_{t}{\boldsymbol{v}}+\nabla{\boldsymbol{v}}\,(\rho{\boldsymbol{v}}+\mathbf{J})=\textnormal{div}\,\mathbf{T}\quad\text{in
$Q$.}$ (2.23)

We now multiply (2.3) by ${1}/{\tilde{\rho}_{i}}$, $i=1,2$, and we add the
resulting equations. Recalling (2.7) and (2.9), this yields

$\displaystyle\begin{aligned}
\textnormal{div}\,{\boldsymbol{v}}&=\partial_{t}(\phi_{1}+\phi_{2})+\textnormal{div}\,{\boldsymbol{v}}\\\
&=\partial_{t}\left(\frac{\rho_{1}}{\tilde{\rho}_{1}}+\frac{\rho_{2}}{\tilde{\rho}_{2}}\right)+\textnormal{div}\left(\frac{\rho_{1}}{\tilde{\rho}_{1}}{\boldsymbol{v}}_{1}+\frac{\rho_{2}}{\tilde{\rho}_{2}}{\boldsymbol{v}}_{2}\right)=0&&\quad\text{in
$Q$,}\end{aligned}$ (2.24)

which means that the volume averaged velocity field is divergence-free.
Multiplying (2.4) by ${1}/{\tilde{\rho}_{i}}$, $i=1,2$, and adding the
resulting equations, we infer

$\displaystyle\begin{aligned}
\textnormal{div}_{\Gamma}\,{\boldsymbol{v}}&=\partial_{t}(\psi_{1}+\psi_{2})+\textnormal{div}_{\Gamma}\,{\boldsymbol{v}}\\\
&=\partial_{t}\left(\frac{\rho_{1}}{\tilde{\rho}_{1}}+\frac{\rho_{2}}{\tilde{\rho}_{2}}\right)+\textnormal{div}_{\Gamma}\left(\frac{\rho_{1}}{\tilde{\rho}_{1}}{\boldsymbol{v}}_{1}+\frac{\rho_{2}}{\tilde{\rho}_{2}}{\boldsymbol{v}}_{2}\right)\\\
&=\beta\left(\frac{\mathbf{J}_{1}}{\tilde{\rho}_{1}}+\frac{\mathbf{J}_{2}}{\tilde{\rho}_{2}}\right)\cdot\mathbf{n}&&\quad\text{on
$\Sigma$.}\end{aligned}$ (2.25)

We now define

$\displaystyle\tilde{\mathbf{J}}_{i}$
$\displaystyle=\frac{\mathbf{J}_{i}}{\tilde{\rho}_{i}},$ $\displaystyle\quad
i=1,2,$ $\displaystyle\qquad\mathbf{J}_{\phi}$
$\displaystyle=\tilde{\mathbf{J}}_{2}-\tilde{\mathbf{J}}_{1}$ in $Q$, (2.26)
$\displaystyle\tilde{\mathbf{K}}_{i}$
$\displaystyle=\frac{\mathbf{K}_{i}}{\tilde{\rho}_{i}},$ $\displaystyle\quad
i=1,2,$ $\displaystyle\qquad\mathbf{K}_{\psi}$
$\displaystyle=\tilde{\mathbf{K}}_{2}-\tilde{\mathbf{K}}_{1}$ on $\Sigma$.
(2.27)

We next multiply (2.14) by ${1}/{\tilde{\rho}_{i}}$, $i=1,2$. Subtracting the
resulting equations, we deduce that

$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{\phi}=0\quad\text{in
$Q$.}$ (2.28)

Moreover, multiplying (2.14) by ${1}/{\tilde{\rho}_{i}}$, adding the resulting
equations, and recalling (2.7) and (2.24), we conclude that

$\displaystyle\textnormal{div}(\tilde{\mathbf{J}}_{1}+\tilde{\mathbf{J}}_{2})=\partial_{t}(\phi_{1}+\phi_{2})+\textnormal{div}\big{(}(\phi_{1}+\phi_{2}){\boldsymbol{v}}\big{)}+\textnormal{div}(\tilde{\mathbf{J}}_{1}+\tilde{\mathbf{J}}_{2})=0\quad\text{in
$Q$.}$ (2.29)

Proceeding similarly with (2.15), and recalling (2.7) and (2.25), we obtain

$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{\psi}=\beta\mathbf{J}_{\phi}\cdot\mathbf{n}\quad\text{on
$\Sigma$,}$ (2.30)

and

$\displaystyle\begin{aligned}
\textnormal{div}_{\Gamma}(\tilde{\mathbf{K}}_{1}+\tilde{\mathbf{K}}_{2})&=\partial_{t}(\psi_{1}+\psi_{2})+\textnormal{div}_{\Gamma}\big{(}(\psi_{1}+\psi_{2}){\boldsymbol{v}}\big{)}\\\
&\quad+\textnormal{div}_{\Gamma}(\tilde{\mathbf{K}}_{1}+\tilde{\mathbf{K}}_{2})-\beta(\tilde{\mathbf{J}}_{1}+\tilde{\mathbf{J}}_{2})\cdot\mathbf{n}=0\quad\text{on
$\Sigma$.}\end{aligned}$ (2.31)

Using (2.26) and (2.29), we infer that

$\displaystyle\begin{aligned}
\textnormal{div}\,\mathbf{J}&=\frac{1}{2}\bigg{[}\tilde{\rho}_{2}\textnormal{div}\,\tilde{\mathbf{J}}_{2}+\tilde{\rho}_{1}\textnormal{div}\,\tilde{\mathbf{J}}_{1}\bigg{]}+\frac{1}{2}\bigg{[}\tilde{\rho}_{2}\textnormal{div}\,\tilde{\mathbf{J}}_{2}+\tilde{\rho}_{1}\textnormal{div}\,\tilde{\mathbf{J}}_{1}\bigg{]}\\\\[4.30554pt]
&=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\textnormal{div}\,\tilde{\mathbf{J}}_{2}-\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\textnormal{div}\,\tilde{\mathbf{J}}_{1}=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\textnormal{div}\,\mathbf{J}_{\phi}&&\quad\text{in
$Q$.}\end{aligned}$ (2.32)

This means that $\mathbf{J}$ and
$\tfrac{1}{2}(\tilde{\rho}_{2}-\tilde{\rho}_{1})\mathbf{J}_{\phi}$ can merely
differ by an additive divergence-free function. As in [6], we thus assume that

$\displaystyle\mathbf{J}=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\mathbf{J}_{\phi}\quad\text{in
$Q$.}$ (2.33)

In the same fashion, we derive the identity

$\displaystyle\textnormal{div}_{\Gamma}\,\mathbf{K}=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\textnormal{div}_{\Gamma}\,\mathbf{K}_{\psi}\quad\text{on
$\Sigma$.}$ (2.34)

From (2.7), (2.8) and (2.16), assuming that $\phi$ and $\psi$ only attain
values in $[-1,1]$, we further infer that the density $\rho=\rho(\phi)$ is
given as

$\displaystyle\rho(\phi)=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\phi+\frac{\tilde{\rho}_{2}+\tilde{\rho}_{1}}{2}\quad\text{in
$Q$,}$ (2.35)

and due to (2.6), the trace $\rho(\phi)|_{\Sigma}$ can be expressed as

$\displaystyle\rho(\phi)|_{\Sigma}=\rho(\psi)=\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}\psi+\frac{\tilde{\rho}_{2}+\tilde{\rho}_{1}}{2}\quad\text{on
$\Sigma$.}$ (2.36)

Plugging (2.35) into (2.17), and (2.36) into (2.18), and using the identities
(2.29) and (2.31), we eventually derive the equations

$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{\phi}$
$\displaystyle=0$ $\displaystyle\quad i=1,2,\quad\text{in $Q$,}$ (2.37)
$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{\psi}$
$\displaystyle=\beta\mathbf{J}_{\phi}\cdot\mathbf{n}$ $\displaystyle\quad
i=1,2,\quad\text{on $\Sigma$.}$ (2.38)

###  2.2 Local energy dissipation laws

Proceeding as in [6], we introduce the following energy density in the bulk:

$\displaystyle
e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)=\frac{\rho(\phi)}{2}\left|{\boldsymbol{v}}\right|^{2}+f(\phi,\nabla\phi)\quad\text{in
$Q$.}$ (2.39)

The first summand on the right-hand side stands for the kinetic energy density
whereas the second summand denotes the free energy density in the bulk. On the
boundary, we introduce the additional energy density

$\displaystyle
e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)=g(\psi,\nabla_{\Gamma}\psi)\quad\text{on
$\Sigma$,}$ (2.40)

where $g$ stands for the free energy density on the surface.

We now consider an arbitrary test volume $V(t)\subset\Omega$, $t\in[0,T]$,
that is transported by the flow associated with ${\boldsymbol{v}}$. Let $\nu$
denote the outer unit normal vector field of $V(t)$. In an isothermal
situation, the second law of thermodynamics leads to the dissipation
inequality

$\displaystyle\begin{aligned}
0&\geq\frac{\mathrm{d}}{\mathrm{d}t}\left[\int_{V(t)}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)\;\mathrm{d}x+\int_{\partial
V(t)\cap\Gamma}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)\;\mathrm{d}S\right]\\\
&\qquad+\int_{\partial
V(t)\cap\Omega}\mathbf{J}_{e}\cdot\nu\;\mathrm{d}S+\int_{\partial_{\Gamma}(\partial
V(t)\cap\Gamma)}\mathbf{K}_{e}\cdot\nu_{\Gamma}\;\mathrm{d}S_{\Gamma}.\end{aligned}$
(2.41)

Here $\mathbf{J}_{e}$ and $\mathbf{K}_{e}$ are energy fluxes that will be
specified later. As we consider a closed system, there is no transfer of
energy over the boundary $\Gamma$ and thus, the domain of the first integral
in the second line is just $\partial V(t)\cap\Omega$ instead of $\partial
V(t)$. We further point out that $\partial_{\Gamma}(\partial
V(t)\cap\Gamma)\subseteq\Gamma$ is to be understood as the relative boundary
of the set $\partial V(t)\cap\Gamma$ within the submanifold $\Gamma$, and
$\nu_{\Gamma}$ stands for the corresponding outer unit normal vector field.
Applying Gauß’s divergence theorem on both integrals in the second line, and
recalling that $\nu=\mathbf{n}$ on $\partial V(t)\cap\Gamma$, we reformulate
(2.41) as

$\displaystyle\begin{aligned}
0&\geq\frac{\mathrm{d}}{\mathrm{d}t}\left[\int_{V(t)}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)\;\mathrm{d}x+\int_{\partial
V(t)\cap\Gamma}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)\;\mathrm{d}S\right]\\\
&\qquad+\int_{V(t)}\textnormal{div}\,\mathbf{J}_{e}\;\mathrm{d}x+\int_{\partial
V(t)\cap\Gamma}\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}\;\mathrm{d}S.\end{aligned}$
(2.42)

Applying the Reynolds transport theorem, we find that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{V(t)}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)\;\mathrm{d}x=\int_{V(t)}\partial_{t}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)+\textnormal{div}\big{(}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi){\boldsymbol{v}}\big{)}\;\mathrm{d}x.$
(2.43)

Similarly, using the transport theorem for evolving hypersurfaces (see, e.g.,
[9, Theorem 32]), we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{\partial
V(t)\cap\Gamma}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)\;\mathrm{d}S=\int_{V(t)\cap\Gamma}\partial_{t}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)+\textnormal{div}_{\Gamma}\big{(}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi){\boldsymbol{v}}\big{)}\;\mathrm{d}S.$
(2.44)

Combining (2.42), (2.43) and (2.44), we thus get

$\displaystyle\begin{aligned}
0&\geq\int_{V(t)}\partial_{t}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)+\textnormal{div}\big{(}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi){\boldsymbol{v}}\big{)}+\textnormal{div}\,\mathbf{J}_{e}\;\mathrm{d}x\\\
&\quad+\int_{\partial
V(t)\cap\Gamma}\partial_{t}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)+\textnormal{div}_{\Gamma}\big{(}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi){\boldsymbol{v}}\big{)}+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}\;\mathrm{d}S.\end{aligned}$
(2.45)

In particular, the above inequality holds true for all test volumes
$V(t)\subset\Omega$ with $\partial V(t)\cap\Gamma=\emptyset$. We thus infer
the local dissipation law in the bulk, which reads as

$\displaystyle
0\geq-\mathcal{D}_{\Omega}:=\partial_{t}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)+\textnormal{div}\big{(}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi){\boldsymbol{v}}\big{)}+\textnormal{div}\,\mathbf{J}_{e}\quad\text{in
$Q$.}$ (2.46)

Let now $\alpha>0$ be arbitrary and suppose that $V(t)\subset\Omega$ is a test
volume with $\left|V(t)\right|$ being sufficiently small such that

$\displaystyle\int_{V(t)}\mathcal{D}_{\Omega}\;\mathrm{d}x<\alpha.$

Then, we use (2.45) to infer that

$\displaystyle\int_{\partial
V(t)\cap\Gamma}\partial_{t}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)+\textnormal{div}_{\Gamma}\big{(}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi){\boldsymbol{v}}\big{)}+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}\;\mathrm{d}S<\alpha.$

Since $\alpha>0$ and the test volume $V(t)$ were arbitrary (except for the
above restriction on $V(t)$), we conclude the following local dissipation law
on the boundary:

$\displaystyle
0\geq-\mathcal{D}_{\Gamma}:=\partial_{t}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)+\textnormal{div}_{\Gamma}\big{(}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi){\boldsymbol{v}}\big{)}+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}\quad\text{on
$\Sigma$.}$ (2.47)

###  2.3 Completion of the model derivation by the Lagrange multiplier
approach

We now complete the model derivation by means of the Lagrange multiplier
approach. We introduce the functions $\mu$ and $\theta$ that will be fixed in
the subsequent approach. In view of the identities (2.37) and (2.38), the
local energy dissipation laws (2.46) and (2.47) can be written as

$\displaystyle 0\geq-\mathcal{D}_{\Omega}$
$\displaystyle=\partial_{t}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi)+\textnormal{div}\big{(}e_{\Omega}({\boldsymbol{v}},\phi,\nabla\phi){\boldsymbol{v}}\big{)}+\textnormal{div}\,\mathbf{J}_{e}$
$\displaystyle\quad-\mu\big{(}\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{\phi}\big{)}$
in $Q$, (2.48) $\displaystyle 0\geq-\mathcal{D}_{\Gamma}$
$\displaystyle=\partial_{t}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi)+\textnormal{div}_{\Gamma}\big{(}e_{\Gamma}(\psi,\nabla_{\Gamma}\psi){\boldsymbol{v}}\big{)}+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}$
$\displaystyle\quad-\theta\big{(}\partial_{t}\psi+\textnormal{div}_{\Gamma}({\boldsymbol{v}}\psi)+\textnormal{div}_{\Gamma}\,\mathbf{K}_{\psi}-\beta\mathbf{J}_{\phi}\cdot\mathbf{n}\big{)}$
on $\Sigma$. (2.49)

Here, the functions $\mu$ and $\theta$ can be understood as Lagrange
multipliers.

In the following, for brevity, we will just write $\rho$, $f$ and $g$ instead
of $\rho(\phi)$, $f(\phi,\nabla\phi)$ and $g(\psi,\nabla_{\Gamma}\psi)$. By
the definition of the energy densities $e_{\Omega}$ and $e_{\Gamma}$ (see
(2.39) and (2.40)), we rewrite (2.3) and (2.3) as

$\displaystyle 0\geq-\mathcal{D}_{\Omega}$
$\displaystyle=\partial_{t}\big{(}\tfrac{1}{2}\rho\left|{\boldsymbol{v}}\right|^{2}\big{)}+\textnormal{div}\big{(}\tfrac{1}{2}\rho\left|{\boldsymbol{v}}\right|^{2}{\boldsymbol{v}}\big{)}+\partial_{t}f+\textnormal{div}(f{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{e}$
$\displaystyle\quad-\mu\big{(}\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})+\textnormal{div}\,\mathbf{J}_{\phi}\big{)}$
in $Q$, (2.50) $\displaystyle 0\geq-\mathcal{D}_{\Gamma}$
$\displaystyle=\partial_{t}g+\textnormal{div}_{\Gamma}(g{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}-\mathbf{J}_{e}\cdot\mathbf{n}$
$\displaystyle\quad-\theta\big{(}\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{\psi}-\beta\mathbf{J}_{\phi}\cdot\mathbf{n}\big{)}$
on $\Sigma$. (2.51)

By the chain rule, the derivatives $\partial_{t}f$ and $\partial_{t}g$ can be
expressed as

$\displaystyle\partial_{t}f$
$\displaystyle=\partial_{\phi}f\,\partial_{t}\phi+\partial_{\nabla\phi}f\,\partial_{t}\nabla\phi$
in $Q$, (2.52) $\displaystyle\partial_{t}g$
$\displaystyle=\partial_{\psi}g\,\partial_{t}\psi+\partial_{\nabla_{\Gamma}\psi}g\,\partial_{t}\nabla_{\Gamma}\psi$
on $\Sigma$. (2.53)

In the following, for any function $h:Q\to\mathbb{R}$, we use the notation

$\displaystyle D_{t}h=\partial_{t}h+(v\cdot\nabla)h$

to denote its material derivative. In the bulk, we proceed exactly as in [6]
to reformulate (2.3) as

$\displaystyle\begin{aligned}
0&\geq\textnormal{div}\big{[}\mathbf{J}_{e}-\tfrac{1}{2}\left|{\boldsymbol{v}}\right|^{2}\mathbf{J}+\mathbf{T}^{\top}{\boldsymbol{v}}-\mu\mathbf{J}_{\phi}+\partial_{\nabla\phi}fD_{t}\phi\big{]}\\\
&\quad+\big{[}\partial_{\phi}f-\textnormal{div}(\partial_{\nabla\phi}f)-\mu\big{]}D_{t}\phi\\\
&\quad-\big{[}\mathbf{T}+\nabla\phi\otimes\partial_{\nabla\phi}f\big{]}:\nabla{\boldsymbol{v}}+\nabla\mu\cdot\mathbf{J}_{\phi}&&\quad\text{in
$Q$.}\end{aligned}$ (2.54)

To ensure (2.54), we now choose the chemical potential $\mu$ and the energy
flux $\mathbf{J}_{e}$ as

$\displaystyle\mu$
$\displaystyle=\partial_{\phi}f-\textnormal{div}\big{(}\partial_{\nabla\phi}f\big{)}$
in $Q$, (2.55) $\displaystyle\mathbf{J}_{e}$
$\displaystyle=\tfrac{1}{2}\left|{\boldsymbol{v}}\right|^{2}\mathbf{J}-\mathbf{T}^{\top}{\boldsymbol{v}}+\mu\mathbf{J}_{\phi}-\partial_{\nabla\phi}fD_{t}\phi$
in $Q$. (2.56)

By these choices, the first two lines of the right-hand side in (2.54) vanish.
We further assume the mass flux $\mathbf{J}_{\phi}$ to be of Fick’s type, that
is

$\displaystyle\mathbf{J}_{\phi}=-m_{\Omega}(\phi)\nabla\mu\quad\text{in $Q$,}$
(2.57)

where $m_{\Omega}=m_{\Omega}(\phi)$ is a nonnegative function representing the
mobility. Hence, (2.54) reduces to

$\displaystyle
0\geq-\big{[}\mathbf{T}+\nabla\phi\otimes\partial_{\nabla\phi}f\big{]}:\nabla{\boldsymbol{v}}-m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\quad\text{in
$Q$.}$ (2.58)

We now define the tensor

$\displaystyle\mathbf{S}=\mathbf{T}+p\mathbf{I}+\nabla\phi\otimes\partial_{\nabla\phi}f\quad\text{in
$Q$.}$

Here, the variable $p$ denotes the pressure, and $\mathbf{I}$ denotes the
identity matrix. $\mathbf{S}$ is the viscous stress tensor that corresponds to
irreversible changes of the energy due to friction. For Newtonian fluids,
$\mathbf{S}$ is usually chosen as

$\displaystyle\mathbf{S}=2\nu(\phi)\mathbf{D}{\boldsymbol{v}},$ (2.59)

where $\nu=\nu(\phi)$ is a nonnegative function representing the viscosity of
the fluids, and $\mathbf{D}{\boldsymbol{v}}$ is the symmetric gradient of
${\boldsymbol{v}}$. By this choice, (2.58) is satisfied since we obtain

$\displaystyle
0\geq-\nu(\phi)\left|\mathbf{D}{\boldsymbol{v}}\right|^{2}-m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\quad\text{in
$Q$,}$ (2.60)

by means of the identity
$p\mathbf{I}:\nabla{\boldsymbol{v}}=p\,\textnormal{div}\,{\boldsymbol{v}}=0$
in $Q$.

We now assume that the energy density $f$ is of Ginzburg–Landau type, that is

$\displaystyle
f(\phi,\nabla\phi)=\frac{\varepsilon}{2}\left|\nabla\phi\right|^{2}+\frac{1}{\varepsilon}F(\phi)\quad\text{in
$Q$.}$ (2.61)

Here, the parameter $\varepsilon>0$ is related to the thickness of the diffuse
interface that separates the two fluids, and $F$ is a potential that usually
exhibits a double-well structure. Hence, the chemical potential $\mu$ reads as

$\displaystyle\mu=-\varepsilon\Delta\phi+\frac{1}{\varepsilon}F^{\prime}(\phi)\quad\text{in
$Q$,}$ (2.62)

and the total stress tensor $\mathbf{T}$ is given by

$\displaystyle\begin{aligned}
\mathbf{T}&=\mathbf{S}-p\mathbf{I}-\nabla\phi\otimes\partial_{\nabla\phi}f\\\
&=2\nu(\phi)\mathbf{D}{\boldsymbol{v}}-p\mathbf{I}-\varepsilon\nabla\phi\otimes\nabla\phi\quad\text{in
$Q$.}\end{aligned}$ (2.63)

Plugging (2.57) into (2.37), we further obtain the equation

$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}\quad\text{in
$Q$.}$ (2.64)

Eventually, recalling (2.22) and (2.63), we use (2.19) to derive the equation

$\displaystyle\partial_{t}(\rho{\boldsymbol{v}})+\textnormal{div}\big{(}{\boldsymbol{v}}\otimes(\rho(\phi){\boldsymbol{v}}+\mathbf{J})\big{)}-\textnormal{div}\big{(}2\nu(\phi)\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla
p=-\varepsilon\textnormal{div}(\nabla\phi\otimes\nabla\phi)\quad\text{in
$Q$,}$ (2.65)

where

$\displaystyle\mathbf{J}=-\frac{\tilde{\rho}_{2}-\tilde{\rho}_{1}}{2}m_{\Omega}(\phi)\nabla\mu\quad\text{in
$Q$.}$ (2.66)

We now consider the local energy dissipation law (2.3) on the boundary.
Recalling the formulas for $\mathbf{J}_{e}$ (see (2.56)) and $\mathbf{T}$ (see
(2.63)), we infer that

$\displaystyle\begin{aligned}
0&\geq\partial_{\psi}g\,\partial_{t}\psi+\textnormal{div}_{\Gamma}\big{(}\partial_{\nabla_{\Gamma}\psi}g\partial_{t}\psi\big{)}-\textnormal{div}_{\Gamma}\big{(}\partial_{\nabla_{\Gamma}\psi}g\big{)}\partial_{t}\psi+\textnormal{div}_{\Gamma}(g{\boldsymbol{v}})+\textnormal{div}_{\Gamma}\,\mathbf{K}_{e}\\\
&\quad-\tfrac{1}{2}(\mathbf{J}\cdot\mathbf{n})\,({\boldsymbol{v}}\cdot{\boldsymbol{v}})+2\nu(\phi)\big{(}\mathbf{D}{\boldsymbol{v}}\,\mathbf{n}\big{)}\cdot{\boldsymbol{v}}-\varepsilon(\nabla\phi\otimes\nabla\phi)\mathbf{n}\cdot{\boldsymbol{v}}-\mu\mathbf{J}_{\phi}\cdot\mathbf{n}\\\
&\quad+(\partial_{\nabla\phi}f\cdot\mathbf{n})\partial_{t}\psi+(\partial_{\nabla\phi}f\cdot\mathbf{n})({\boldsymbol{v}}\cdot\nabla_{\Gamma}\psi)-\theta\partial_{t}\psi-\textnormal{div}_{\Gamma}(\theta\psi{\boldsymbol{v}})+\psi\nabla_{\Gamma}\theta\cdot{\boldsymbol{v}}\\\
&\quad-\textnormal{div}_{\Gamma}(\mathbf{K}_{\psi}\theta)+\nabla_{\Gamma}\theta\,\mathbf{K}_{\psi}+\beta\theta\mathbf{J}_{\phi}\cdot\mathbf{n}\\\\[4.30554pt]
&=\textnormal{div}_{\Gamma}\big{[}\mathbf{K}_{e}+\partial_{\nabla_{\Gamma}\psi}g\partial_{t}\psi+g{\boldsymbol{v}}-\theta\psi{\boldsymbol{v}}-\theta\mathbf{K}_{\psi}\big{]}\\\
&\quad+\big{[}\partial_{\psi}g-\textnormal{div}_{\Gamma}\big{(}\partial_{\nabla_{\Gamma}\psi}g\big{)}+\partial_{\nabla\phi}f\cdot\mathbf{n}-\theta\big{]}\partial_{t}\psi\\\
&\quad+\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\psi\nabla_{\Gamma}\theta-\tfrac{1}{2}(\mathbf{J}\cdot\mathbf{n}){\boldsymbol{v}}\big{]}\cdot{\boldsymbol{v}}\\\
&\quad+(\partial_{\nabla\phi}f\cdot\mathbf{n})({\boldsymbol{v}}\cdot\nabla_{\Gamma}\psi)-\varepsilon(\nabla\phi\otimes\nabla\phi)\mathbf{n}\cdot{\boldsymbol{v}}\\\
&\quad+\nabla_{\Gamma}\theta\,\mathbf{K}_{\psi}+(\beta\theta-\mu)\mathbf{J}_{\phi}\cdot\mathbf{n}\end{aligned}$
(2.67)

on $\Sigma$. Here, we have used that

$\displaystyle{\boldsymbol{v}}\cdot\nabla_{\Gamma}\psi={\boldsymbol{v}}\cdot\nabla_{\Gamma}\phi={\boldsymbol{v}}\cdot\big{[}\nabla\phi-\mathbf{n}(\nabla\phi\cdot\mathbf{n})\big{]}={\boldsymbol{v}}\cdot\nabla\phi\quad\text{on
$\Sigma$}$ (2.68)

due to the definition of the surface gradient and the boundary condition
(2.11). In order to ensure that the inequality (2.67) is fulfilled, we choose
the chemical potential $\theta$, the mass flux $\mathbf{K}_{\psi}$, and the
energy flux $\mathbf{K}_{e}$ as follows:

$\displaystyle\theta$
$\displaystyle=\partial_{\psi}g-\textnormal{div}_{\Gamma}\big{(}\partial_{\nabla_{\Gamma}\psi}g\big{)}+\partial_{\nabla\phi}f\cdot\mathbf{n}$
on $\Sigma$, (2.69) $\displaystyle\mathbf{K}_{\psi}$
$\displaystyle=-m_{\Gamma}(\psi)\nabla_{\Gamma}\theta$ on $\Sigma$, (2.70)
$\displaystyle\mathbf{K}_{e}$
$\displaystyle=-\partial_{\nabla_{\Gamma}\psi}gD_{t}\psi+\theta\psi{\boldsymbol{v}}+\theta\mathbf{K}_{\psi}-g{\boldsymbol{v}}$
on $\Sigma$. (2.71)

In (2.70), $m_{\Gamma}=m_{\Gamma}(\psi)$ is a nonnegative function
representing the mobility. This means that the first two lines of the right-
hand side in (2.67) vanish. Moreover, since ${\boldsymbol{v}}_{\mathbf{n}}=0$
(see (2.11)), we have ${\boldsymbol{v}}={\boldsymbol{v}}_{\tau}$ and thus,

$\displaystyle\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\psi\nabla_{\Gamma}\theta-\tfrac{1}{2}(\mathbf{J}\cdot\mathbf{n}){\boldsymbol{v}}\big{]}\cdot{\boldsymbol{v}}$
$\displaystyle=\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\psi\nabla_{\Gamma}\theta-\tfrac{1}{2}(\mathbf{J}\cdot\mathbf{n}){\boldsymbol{v}}\big{]}_{\tau}\cdot{\boldsymbol{v}}_{\tau}$

Therefore, in order to ensure that (2.67) is satisfied, we make the
constitutive assumption

$\displaystyle\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\psi\nabla_{\Gamma}\theta-\tfrac{1}{2}(\mathbf{J}\cdot\mathbf{n}){\boldsymbol{v}}\big{]}_{\tau}=-\gamma(\psi){\boldsymbol{v}}_{\tau}.$
(2.72)

where $\gamma=\gamma(\psi)$ is a nonnegative function. This equation can be
regarded as an inhomogeneous Navier slip boundary condition. We now assume
that the energy density $g$ is also of Ginzburg–Landau type, that is

$\displaystyle
g(\psi,\nabla_{\Gamma}\psi)=\frac{\kappa\delta}{2}\left|\nabla_{\Gamma}\psi\right|^{2}+\frac{1}{\delta}G(\psi)\quad\text{on
$\Sigma$.}$ (2.73)

Here, $\delta>0$ is related to the thickness of the diffuse interface on the
boundary, $\kappa\geq 0$ acts as a weight for surface diffusion effects, and
$G$ is a potential that usually exhibits a double-well structure. Hence,
recalling the definition of the energy density $f$ (see (2.61)), the chemical
potential $\theta$ is given as

$\displaystyle\theta=-\kappa\delta\Delta_{\Gamma}\psi+\frac{1}{\delta}G^{\prime}(\psi)-\varepsilon\partial_{\mathbf{n}}\phi\quad\text{on
$\Sigma$.}$ (2.74)

Thanks to (2.68), we further have

$\displaystyle(\partial_{\nabla\phi}f\cdot\mathbf{n})({\boldsymbol{v}}\cdot\nabla_{\Gamma}\psi)=\varepsilon(\nabla\phi\cdot\mathbf{n})({\boldsymbol{v}}\cdot\nabla\phi)=\varepsilon(\nabla\phi\otimes\nabla\phi)\mathbf{n}\cdot{\boldsymbol{v}}\quad\text{on
$\Sigma$.}$ (2.75)

Hence, (2.67) reduces to

$\displaystyle 0$
$\displaystyle\geq-\gamma(\psi)\left|{\boldsymbol{v}}\right|^{2}-m_{\Gamma}(\psi)\left|\nabla_{\Gamma}\theta\right|^{2}-(\beta\theta-\mu)\,m_{\Omega}(\phi)\nabla\mu\cdot\mathbf{n}\quad\text{on
$\Sigma$.}$ (2.76)

The first two terms on the right-hand side are clearly nonpositive as
$\gamma(\psi)$ and $m_{\Gamma}(\psi)$ are nonnegative. This means that (2.76)
is fulfilled if one of the following boundary conditions holds:

$\displaystyle\beta\theta-\mu$ $\displaystyle=0$ on $\Sigma$, (2.77a)
$\displaystyle m_{\Omega}(\phi)\nabla\mu\cdot\mathbf{n}$
$\displaystyle=\tfrac{1}{L}(\beta\theta-\mu)$ on $\Sigma$ for a constant
$L\in(0,\infty)$, (2.77b) $\displaystyle
m_{\Omega}(\phi)\nabla\mu\cdot\mathbf{n}$ $\displaystyle=0$ on $\Sigma$.
(2.77c)

Finally, by substituting (2.57) and (2.70) into (2.38), we obtain

$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}})=\textnormal{div}_{\Gamma}(m_{\Gamma}(\psi)\nabla_{\Gamma}\theta)-\beta
m_{\Omega}(\psi)\nabla\mu\cdot\mathbf{n}\quad\text{on $\Sigma$.}$ (2.78)

Collecting (2.65), (2.24), (2.66), (2.64), (2.62), (2.78), (2.74), (2.11),
(2.6), (2.77) and (2.72), and imposing the initial conditions

$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}\quad\text{in
$\Omega$},\qquad\psi|_{t=0}=\psi_{0}\quad\text{on $\Sigma$,}$ (2.79)

for prescribed initial data ${\boldsymbol{v}}_{0}$, $\phi_{0}$ and $\psi_{0}$,
we have thus derived the system (1.17).

##  3 Mathematical analysis: the case of matched densities

###  3.1 Notation and preliminaries

We fix some notation and assumptions that are supposed to hold throughout the
remainder of this paper.

#### Notation.

  1. $(\mathrm{N1})$

The set of natural numbers excluding zero is denoted by $\mathbb{N}$, and we
write $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$.

  2. $(\mathrm{N2})$

Let $\Omega$ be a bounded Lipschitz domain in $\mathbb{R}^{d}$ with
$d\in\left\\{2,3\right\\}$. We write $\Gamma:=\partial\Omega$. For any real
numbers $k\geq 0$ and $p\in[1,\infty]$, the Lebesgue and Sobolev spaces for
functions defined on $\Omega$ with values in $\mathbb{R}$ are denoted as
$L^{p}(\Omega)$ and $W^{k,p}(\Omega)$. We write $\|\cdot\|_{L^{p}(\Omega)}$
and $\|\cdot\|_{W^{k,p}(\Omega)}$ to denote the standard norms on these
spaces. If $p=2$, we use the notation $H^{k}(\Omega)=W^{k,2}(\Omega)$. We
point out that $H^{0}(\Omega)$ can be identified with $L^{2}(\Omega)$.
Analogously, the Lebesgue and Sobolev spaces on $\Gamma$ are denoted by
$L^{p}(\Gamma)$ and $W^{k,p}(\Gamma)$ with corresponding norms
$\|\cdot\|_{L^{p}(\Gamma)}$ and $\|\cdot\|_{W^{k,p}(\Gamma)}$, respectively.
In the case of vector fields on $\Omega$ with values in $\mathbb{R}^{n}$ for
some $n\in\mathbb{N}$ with $n>1$, we use the notation
$\mathbf{L}^{p}(\Omega)$, $\mathbf{W}^{k,p}(\Omega)$ and
$\mathbf{H}^{k}(\Omega)$. For simplicity, their norms are denoted as in the
scalar case by $\|\cdot\|_{L^{p}(\Omega)}$, $\|\cdot\|_{W^{k,p}(\Omega)}$ and
$\|\cdot\|_{H^{k}(\Omega)}$, respectively.

Let $I$ be a closed interval in $\mathbb{R}$ and $X$ be a Banach space. The
space $C(I;X)$ denotes the set of continuous functions from $I$ to $X$ and,
for $k\in\mathbb{N}$, $C^{k}(I;X)$ denotes the space of $k$-times continuously
differentiable functions from $I$ to $X$. In particular, we simply write
$C(I)$ and $C^{k}(I)$ if $X=\mathbb{R}$. Moreover, $C_{w}(I;X)$ denotes the
space of functions mapping from $I$ to $X$ which are continuous on $I$ with
respect to the weak topology on $X$. This means that for any function $f\in
C_{w}(I;X)$ and every sequence $(t_{k})_{k\in\mathbb{N}}$ in $I$ with
$t_{k}\to t$ in $I$ as $k\to\infty$, it holds $f(t_{k})\to f(t)$ weakly in $I$
as $k\to\infty$. Furthermore, for any real numbers $k\geq 0$ and
$p\in[1,\infty]$, the Bochner spaces of functions defined on an interval $I$
in $\mathbb{R}$ with values in $X$ are denoted by $L^{p}(I;X)$ and
$W^{k,p}(I;X)$.

  3. $(\mathrm{N3})$

For any Banach space $X$, we write $X^{\prime}$ to denote its dual space. The
associated duality pairing of elements $\phi\in X^{\prime}$ and $\zeta\in X$
is denoted as $\langle\phi,\zeta\rangle_{X}$. If $X$ is a Hilbert space, we
write $(\cdot,\cdot)_{X}$ to denote its inner product.

  4. $(\mathrm{N4})$

For any bounded Lipschitz domain $\Omega$, we define

$\displaystyle\langle
u\rangle_{\Omega}:=\begin{cases}\frac{1}{\left|\Omega\right|}\langle
u,1\rangle_{H^{1}(\Omega)}&\text{ if }u\in H^{1}(\Omega)^{\prime},\\\
\frac{1}{\left|\Omega\right|}\int_{\Omega}u\;\mathrm{d}x&\text{ if }u\in
L^{1}(\Omega)\end{cases}$

to denote the (generalized) spatial mean of $u$. Here, $\left|\Omega\right|$
denotes the $d$-dimensional Lebesgue measure of $\Omega$. The spatial mean of
a function $v\in H^{1}(\Gamma)^{\prime}$ (or $v\in L^{1}(\Gamma)$,
respectively) is denoted as $\langle v\rangle_{\Gamma}$ and defined
analogously.

  5. $(\mathrm{N5})$

For any bounded Lipschitz domain $\Omega$ and $d\in\\{2,3\\}$, we introduce
the space

$\mathbf{L}^{2}_{\textnormal{div}}(\Omega):=\left\\{{\boldsymbol{v}}\in
L^{2}(\Omega;\mathbb{R}^{d}):\ \textnormal{div}\,{\boldsymbol{v}}=0\
\text{in}\ \Omega,\ {\boldsymbol{v}}\cdot\mathbf{n}=0\ \text{on}\
\Gamma\right\\}=\overline{C_{0,\sigma}^{\infty}(\Omega;\mathbb{R}^{d})}^{L^{2}(\Omega)},$

where
$C_{0,\sigma}^{\infty}(\Omega;\mathbb{R}^{d}):=\left\\{{\boldsymbol{v}}\in
C_{c}^{\infty}(\Omega;\mathbb{R}^{d}):\ \textnormal{div}\,{\boldsymbol{v}}=0\
\text{in}\ \Omega\right\\}$. Moreover, we define

$\mathbf{H}^{1}_{\textnormal{div}}(\Omega)=\mathbf{H}^{1}(\Omega)\cap\mathbf{L}^{2}_{\textnormal{div}}(\Omega).$

#### Assumptions.

To prove the existence of weak solutions in the case of matched densities, we
make the following assumptions:

  1. $(\mathrm{A1})$

We consider a bounded domain $\Omega\subset\mathbb{R}^{d}$ with
$d\in\\{2,3\\}$ with Lipschitz boundary $\Gamma=\partial\Omega$ and a final
time $T>0$. We further use the notation

$\displaystyle Q:=\Omega\times(0,T),\quad\Sigma:=\Gamma\times(0,T).$

  2. $(\mathrm{A2})$

The constants occurring in the system (1.22) satisfy
$\varepsilon,\delta,\beta,\kappa>0$ and $L=0$. Since the choice of $\delta$,
$\varepsilon$ and $\kappa$ has no impact on the mathematical analysis, we will
simply set $\delta=\varepsilon=\kappa=1$ without loss of generality.

  3. $(\mathrm{A3})$

The potentials $F:\mathbb{R}\to[0,\infty)$ and $G:\mathbb{R}\to[0,\infty)$ are
twice continuously differentiable and there exist two exponents $p$ and $q$
satisfying

$\displaystyle p\in\begin{cases}[2,\infty),&\text{if $d=2$,}\\\
[2,6],&\text{if $d=3$,}\end{cases}\quad\text{and}\quad q\in[2,\infty)$

as well as constants $c_{F^{\prime\prime}},c_{G^{\prime\prime}}\geq 0$ such
that the second-order derivatives satisfy the growth conditions

$\displaystyle\left|F^{\prime\prime}(s)\right|$ $\displaystyle\leq
c_{F^{\prime\prime}}(1+\left|s\right|^{p-2}),$ (3.1)
$\displaystyle\left|G^{\prime\prime}(s)\right|$ $\displaystyle\leq
c_{G^{\prime\prime}}(1+\left|s\right|^{q-2}),$ (3.2)

for all $s\in\mathbb{R}$.

These assumptions already entail that there exist constants
$c_{F},c_{G},c_{F^{\prime}},c_{G^{\prime}}\geq 0$ such that $F^{\prime}$,
$G^{\prime}$, $F$ and $G$ satisfy the growth conditions

$\displaystyle\left|F^{\prime}(s)\right|$ $\displaystyle\leq
c_{F^{\prime}}(1+\left|s\right|^{p-1}),$ (3.3)
$\displaystyle\left|G^{\prime}(s)\right|$ $\displaystyle\leq
c_{G^{\prime}}(1+\left|s\right|^{q-1}),$ (3.4) $\displaystyle F(s)$
$\displaystyle\leq c_{F}(1+\left|s\right|^{p}),$ (3.5) $\displaystyle G(s)$
$\displaystyle\leq c_{G}(1+\left|s\right|^{q}),$ (3.6)

for all $s\in\mathbb{R}$.

  4. $(\mathrm{A4})$

The mobility functions $m_{\Omega}:\mathbb{R}\to\mathbb{R}$ and
$m_{\Gamma}:\mathbb{R}\to\mathbb{R}$ are continuous, bounded and uniformly
positive. This means that there exist positive constants $m_{\Omega}^{*}$,
$M_{\Omega}^{*}$, $m_{\Gamma}^{*}$, $M_{\Gamma}^{*}$ such that for all
$s\in\mathbb{R}$,

$\displaystyle 0<m_{\Omega}^{*}\leq m_{\Omega}(s)\leq
M_{\Omega}^{*}\quad\text{and}\quad 0<m_{\Gamma}^{*}\leq m_{\Gamma}(s)\leq
M_{\Gamma}^{*}.$ (3.7)

  5. $(\mathrm{A5})$

The viscosity of the mixture $\nu:\mathbb{R}\to\mathbb{R}$ and the friction
parameter $\gamma:\mathbb{R}\to\mathbb{R}$ are continuous, bounded and
uniformly positive. This means that there exist positive constants $\nu_{0}$,
$\nu_{1}$, $\gamma_{0}$, $\gamma_{1}$ such that for all $s\in\mathbb{R}$,

$\displaystyle 0<\nu_{0}\leq\nu(s)\leq\nu_{1},\quad\text{and}\quad
0<\gamma_{0}\leq\gamma(s)\leq\gamma_{1}.$ (3.8)

######  Remark 3.1.

We point out that the polynomial double-well potential

$\displaystyle W_{\text{dw}}(s)=\tfrac{1}{4}(s^{2}-1)^{2},\quad
s\in\mathbb{R},$

is a suitable choice for $F$ and $G$ as it satisfies ($\mathrm{A3}$) with
$p=4$ and $q=4$. However, singular potentials like the logarithmic
Flory–Huggins potential or the double-obstacle potential are not admissible in
this setting as they do not even satisfy (3.1)–(3.2).

#### Preliminaries.

We next introduce several function spaces, inner products, norms and operators
that will be used throughout this paper.

  1. $(\mathrm{P1})$

For any real numbers $k\geq 0$ and $p\in[1,\infty]$, we set

$\displaystyle\mathcal{L}^{p}:=L^{p}(\Omega)\times
L^{p}(\Gamma),\quad\text{and}\quad\mathcal{H}^{k}:=H^{k}(\Omega)\times
H^{k}(\Gamma),$

and we identify $\mathcal{L}^{2}$ with $\mathcal{H}^{0}$. Note that
$\mathcal{H}^{k}$ is a Hilbert space with respect to the inner product

$\displaystyle\big{(}(\phi,\psi),(\zeta,\xi)\big{)}_{\mathcal{H}^{k}}:=\big{(}\phi,\zeta\big{)}_{H^{k}(\Omega)}+\big{(}\psi,\xi\big{)}_{H^{k}(\Gamma)}\quad\text{for
all $(\phi,\psi),(\zeta,\xi)\in\mathcal{H}^{k}$,}$

and its induced norm
$\|\cdot\|_{\mathcal{H}^{k}}:=\left(\cdot,\cdot\right)_{\mathcal{H}^{k}}^{1/2}$.

  2. $(\mathrm{P2})$

For $\beta>0$, we introduce the subspace

$\displaystyle{\mathcal{D}_{\beta}}:=\left\\{(\phi,\psi)\in\mathcal{H}^{1}\;\big{|}\;\phi|_{\Gamma}=\beta\psi\;\;\text{a.e.
on}\;\Gamma\right\\}\subset\mathcal{H}^{1},$

endowed with the inner product
$\left(\cdot,\cdot\right)_{\mathcal{D}_{\beta}}:=\left(\cdot,\cdot\right)_{\mathcal{H}^{1}}$
and its induced norm. The space ${\mathcal{D}_{\beta}}$ is a Hilbert space.
Moreover, we define the product

$\displaystyle\big{<}(\phi,\psi),(\zeta,\xi)\big{>}_{{\mathcal{D}_{\beta}}}:=\left(\phi,\zeta\right)_{L^{2}(\Omega)}+\left(\psi,\xi\right)_{L^{2}(\Gamma)}$

for all $(\phi,\psi),(\zeta,\xi)\in\mathcal{L}^{2}$. By means of the Riesz
representation theorem, this product can be extended to a duality pairing on
${\mathcal{D}_{\beta}^{\prime}}\times{\mathcal{D}_{\beta}}$, which will also
be denoted as $\langle\cdot,\cdot\rangle_{{\mathcal{D}_{\beta}}}$.

In particular, the spaces
$\big{(}{\mathcal{D}_{\beta}},\mathcal{L}^{2},{\mathcal{D}_{\beta}^{\prime}}\big{)}$
form a Gelfand triplet, and the operator norm on
${\mathcal{D}_{\beta}^{\prime}}$ is given by

$\displaystyle\|(\phi,\psi)\|_{{\mathcal{D}_{\beta}^{\prime}}}:=\sup\Big{\\{}\,\big{|}\langle(\phi,\psi),(\zeta,\xi)\rangle_{{\mathcal{D}_{\beta}}}\big{|}\;\Big{|}\;(\zeta,\xi)\in{\mathcal{D}_{\beta}}\text{
with }\|(\zeta,\xi)\|_{{\mathcal{D}_{\beta}}}=1\Big{\\}},$

for all $(\phi,\psi)\in{\mathcal{D}_{\beta}^{\prime}}$.

  3. $(\mathrm{P3})$

We further recall the following bulk-surface Poincaré inequalities:

    1. $(\mathrm{P3.1})$

For any $\alpha>0$, there exists a constant $C_{P}$ depending only on $\alpha$
and $\Omega$ such that

$\|(\phi,\psi)\|_{\mathcal{L}^{2}}\leq
C_{P}\|(\nabla\phi,\nabla_{\Gamma}\psi)\|_{\mathcal{L}^{2}}$ (3.9)

for all $(\phi,\psi)\in{\mathcal{D}_{\alpha}}$ with
$\alpha\left|\Omega\right|\langle\phi\rangle_{\Omega}+\left|\Gamma\right|\langle\psi\rangle_{\Gamma}=0$.
This bulk-surface Poincaré inequality is a special case of the one established
in [40, Lemma A.1] (with the parameters therein being chosen as $K=0$ and
$\alpha,\beta>0$).

    2. $(\mathrm{P3.2})$

There exists a constant $C_{P}^{\prime}$ depending only on $\Omega$ such that

$\|u\|_{L^{2}(\Omega)}\leq C_{P}^{\prime}\left(\|\nabla
u\|_{L^{2}(\Omega)}+\|u\|_{L^{2}(\Gamma)}\right)$ (3.10)

for all $u\in H^{1}(\Omega)$.
This bulk-surface Poincaré inequality is established in [53, Chapter II,
Section 1.4].

###  3.2 Existence and uniqueness of weak solutions

In this section, we consider the case of matched densities, that is
$\rho\equiv\tilde{\rho}_{1}=\tilde{\rho}_{2}$. This means we consider the
system (1.22) supplemented with the boundary condition (1.22g) with $L=0$ and
$\beta>0$. For the analysis, without loss of generality, we set
$\rho=\kappa=\varepsilon=\delta=1$. We further employ the usual decomposition

$\displaystyle-\,\textnormal{div}\big{(}\nabla\phi\otimes\nabla\phi\big{)}=-\nabla\left(\frac{1}{2}\left|\nabla\phi\right|^{2}+F(\phi)\right)+\mu\nabla\phi.$
(3.11)

Hence, replacing the pressure $p$ by

$\displaystyle\overline{p}:=p+\left(\frac{1}{2}\left|\nabla\phi\right|^{2}+F(\phi)\right),$

the system (1.22) can be restated as follows

$\displaystyle\partial_{t}{\boldsymbol{v}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla\overline{p}=\mu\nabla\phi$
$\displaystyle\quad\text{in $Q$},$ (3.12a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0,$
$\displaystyle\quad\text{in $Q$},$ (3.12b)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\quad\text{in $Q$},$ (3.12c)
$\displaystyle\mu=-\Delta\phi+F^{\prime}(\phi)$ $\displaystyle\quad\text{in
$Q$},$ (3.12d)
$\displaystyle\partial_{t}\psi+\textnormal{div}_{\Gamma}(\psi{\boldsymbol{v}}_{\tau})=\textnormal{div}_{\Gamma}\big{(}m_{\Gamma}(\psi)\nabla_{\Gamma}\theta\big{)}-\beta
m_{\Omega}(\psi)\,\partial_{\mathbf{n}}\mu$ $\displaystyle\quad\text{on
$\Sigma$},$ (3.12e)
$\displaystyle\theta=-\Delta_{\Gamma}\psi+G^{\prime}(\psi)+\partial_{\mathbf{n}}\phi$
$\displaystyle\quad\text{on $\Sigma$},$ (3.12f)
$\displaystyle{\boldsymbol{v}}\cdot\mathbf{n}=0,\qquad\phi|_{\Gamma}=\psi,\qquad\mu|_{\Gamma}=\beta\theta$
$\displaystyle\quad\text{on $\Sigma$},$ (3.12g)
$\displaystyle\big{[}2\nu(\psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\gamma(\psi){\boldsymbol{v}}\big{]}_{\tau}=\big{[}-\psi\nabla_{\Gamma}\theta\big{]}_{\tau}$
$\displaystyle\quad\text{on $\Sigma$},$ (3.12h)
$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}$
$\displaystyle\quad\text{in $\Omega$},$ (3.12i)
$\displaystyle\psi|_{t=0}=\psi_{0}$ $\displaystyle\quad\text{on $\Gamma$}.$
(3.12j)

Here, we recall that in the case of matched densities, the term $\mathbf{J}$
vanishes. Therefore, in the Navier–Stokes equation (3.12a) and in the
corresponding boundary condition (3.12h), the terms related to $\mathbf{J}$ do
not appear any more. The total energy associated with this system is

$\displaystyle
E({\boldsymbol{v}},\phi,\psi)=\int_{\Omega}\frac{1}{2}\left|{\boldsymbol{v}}\right|^{2}\;\mathrm{d}x+\int_{\Omega}\frac{1}{2}\left|\nabla\phi\right|^{2}+F(\phi)\;\mathrm{d}x+\int_{\Gamma}\frac{1}{2}\left|\nabla_{\Gamma}\psi\right|^{2}+G(\psi)\;\mathrm{d}S.$
(3.13)

The notion of weak solutions to system (3.12) is defined as follows:

######  Definition 3.2 (Weak solutions of system (3.12)).

Suppose that the assumptions ($\mathrm{A1}$)–($\mathrm{A5}$) hold. Let
${\boldsymbol{v}}_{0}\in\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$ and
$(\phi_{0},\psi_{0})\in\mathcal{D}_{1}$ be arbitrary. The quintuplet
$({\boldsymbol{v}},\phi,\psi,\mu,\theta)$ is called a weak solution of the
system (3.12) on $[0,T]$ if the following properties hold:

  1. (i)

The functions ${\boldsymbol{v}}$, $\phi$, $\psi$, $\mu$ and $\theta$ have the
following regularity:

$\displaystyle{\boldsymbol{v}}\in
L^{\infty}\big{(}0,T;\mathbf{L}_{\textnormal{div}}^{2}(\Omega)\big{)}\cap
L^{2}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)\big{)},$ (3.14a)
$\displaystyle\left\\{\begin{aligned} &{\boldsymbol{v}}\in
W^{1,2}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime}\big{)}\cap
C\big{(}[0,T];\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\big{)}\ &&\text{if }\
d=2,\\\ &{\boldsymbol{v}}\in
W^{1,\frac{4}{3}}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime}\big{)}\cap
C_{w}\big{(}[0,T];\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\big{)}\ &&\text{if
}\ d=3,\end{aligned}\right.$ (3.14b) $\displaystyle(\phi,\psi)\in
C([0,T];\mathcal{L}^{2})\cap H^{1}(0,T;{\mathcal{D}_{\beta}^{\prime}})\cap
L^{\infty}(0,T;\mathcal{D}_{1}),$ (3.14c) $\displaystyle(\mu,\theta)\in
L^{4}(0,T;\mathcal{L}^{2})\cap L^{2}(0,T;\mathcal{D}_{\beta}).$ (3.14d)

  2. (ii)

The functions ${\boldsymbol{v}}$, $\phi$ and $\psi$ satisfy the initial
conditions

$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}\quad\text{a.e.\leavevmode\nobreak\
in
$\Omega$},\quad\text{and}\quad\psi|_{t=0}=\psi_{0}\quad\text{a.e.\leavevmode\nobreak\
on $\Gamma$}.$ (3.15)

  3. (iii)

The functions ${\boldsymbol{v}}$, $\phi$, $\psi$, $\mu$ and $\theta$ satisfy
the weak formulation

$\displaystyle\begin{aligned}
&\langle\partial_{t}{\boldsymbol{v}},{\boldsymbol{w}}\rangle_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}-\int_{\Omega}({\boldsymbol{v}}\otimes{\boldsymbol{v}}):\nabla{\boldsymbol{w}}\;\mathrm{d}x+\int_{\Omega}2\nu(\phi)\mathbf{D}{\boldsymbol{v}}:\mathbf{D}{\boldsymbol{w}}\;\mathrm{d}x\\\
&\quad=\int_{\Omega}\mu\nabla\phi\cdot{\boldsymbol{w}}\;\mathrm{d}x-\int_{\Gamma}\gamma(\psi){\boldsymbol{v}}\cdot{\boldsymbol{w}}+\psi\nabla_{\Gamma}\theta\cdot{\boldsymbol{w}}\;\mathrm{d}S,\end{aligned}$
(3.16a) $\displaystyle\begin{aligned}
&\big{<}(\partial_{t}\phi,\partial_{t}\psi),(\zeta,\xi)\big{>}_{{\mathcal{D}_{\beta}}}-\int_{\Omega}\phi{\boldsymbol{v}}\cdot\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}\psi{\boldsymbol{v}}\cdot\nabla_{\Gamma}\xi\;\mathrm{d}S\\\
&\quad=-\int_{\Omega}m_{\Omega}(\phi)\nabla\mu\cdot\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}m_{\Gamma}(\psi)\nabla_{\Gamma}\theta\cdot\nabla_{\Gamma}\xi\;\mathrm{d}S,\end{aligned}$
(3.16b) $\displaystyle\begin{aligned}
&\int_{\Omega}\mu\,\eta\;\mathrm{d}x+\int_{\Gamma}\theta\,\vartheta\;\mathrm{d}S\\\
&\quad=\int_{\Omega}\nabla\phi\cdot\nabla\eta+F^{\prime}(\phi)\eta\;\mathrm{d}x+\int_{\Gamma}\nabla_{\Gamma}\psi\cdot\nabla_{\Gamma}\vartheta+G^{\prime}(\psi)\vartheta\;\mathrm{d}S\end{aligned}$
(3.16c)

for all test functions
${\boldsymbol{w}}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$,
$(\zeta,\xi)\in{\mathcal{D}_{\beta}}$, $(\eta,\vartheta)\in\mathcal{D}_{1}$.

  4. (iv)

The functions $\phi$ and $\psi$ satisfy the mass conservation law

$\displaystyle\beta\int_{\Omega}\phi(t)\;\mathrm{d}x+\int_{\Gamma}\psi(t)\;\mathrm{d}S=\beta\int_{\Omega}\phi_{0}\;\mathrm{d}x+\int_{\Gamma}\psi_{0}\;\mathrm{d}S=:m\quad\text{for
all $t\in[0,T]$}.$ (3.17)

  5. (v)

The functions ${\boldsymbol{v}}$, $\phi$, $\psi$, $\mu$ and $\theta$ satisfy
the weak energy dissipation law

$\displaystyle\begin{aligned}
&E\big{(}{\boldsymbol{v}}(t),\phi(t),\psi(t)\big{)}+\int_{0}^{t}\int_{\Omega}2\nu(\phi)\,\left|\mathbf{D}{\boldsymbol{v}}\right|^{2}\;\mathrm{d}x\;\mathrm{d}t+\int_{0}^{t}\int_{\Gamma}\gamma(\psi)\left|{\boldsymbol{v}}\right|^{2}\;\mathrm{d}S\;\mathrm{d}t\\\
&\quad+\int_{0}^{t}\int_{\Gamma}m_{\Gamma}(\psi)\left|\nabla_{\Gamma}\theta\right|^{2}\;\mathrm{d}S\;\mathrm{d}t+\int_{0}^{t}\int_{\Omega}m_{\Omega}(\phi)\left|\nabla\mu\right|^{2}\;\mathrm{d}x\;\mathrm{d}t\leq
E({\boldsymbol{v}}_{0},\phi_{0},\psi_{0})\end{aligned}$ (3.18)

for all $t\in[0,T]$.

The existence of such a weak solution is ensured by the following theorem.

######  Theorem 3.3.

Suppose that the assumptions ($\mathrm{A1}$)–($\mathrm{A5}$) hold and let
${\boldsymbol{v}}_{0}\in\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$ and
$(\phi_{0},\psi_{0})\in\mathcal{D}_{1}$ be arbitrary. Then, there exists a
weak solution of system (3.12) in the sense of Definition 3.2. If we
additionally suppose that the domain $\Omega$ is of class $C^{k}$ with
$k\in\\{2,3\\}$, and in the case $d=3$, we further assume that assumption
($\mathrm{A3}$) holds with $p<6$, then $(\phi,\psi)\in
L^{2}(0,T;\mathcal{H}^{k})$ as well as

$\displaystyle\mu$ $\displaystyle=-\Delta\phi+F^{\prime}(\phi)$
$\displaystyle\quad\text{a.e.\leavevmode\nobreak\ in $Q$},$ (3.19a)
$\displaystyle\theta$
$\displaystyle=-\Delta_{\Gamma}\psi+G^{\prime}(\psi)+\partial_{\mathbf{n}}\phi$
$\displaystyle\quad\text{a.e.\leavevmode\nobreak\ on $\Sigma$}.$ (3.19b)

In two dimensions, we show the uniqueness of the weak solutions under
additional assumptions on the domain $\Omega$, the viscosity $\nu$, the
mobilities $m_{\Omega}$ and $m_{\Gamma}$ and the potentials $F$ and $G$.

######  Theorem 3.4.

Suppose that the assumptions ($\mathrm{A1}$)–($\mathrm{A5}$) hold with $d=2$.
We additionally assume that $\Omega$ is of class $C^{3}$, $\gamma$ is
differentiable and $\gamma^{\prime}$ is locally bounded, the functions $\nu$,
$m_{\Omega}$ and $m_{\Gamma}$ are constant, and the potentials $F$ and $G$
satisfy the relation

$\displaystyle F(s)=\beta G(s)\quad\text{for all $s\in\mathbb{R}$.}$ (3.20)

Then, the weak solution of system (1.22) given by Theorem 3.3 is unique. In
addition, given two weak solutions
$({\boldsymbol{v}}_{1},\phi_{1},\psi_{1},\mu_{1},\theta_{1})$ and
$({\boldsymbol{v}}_{2},\phi_{2},\psi_{2},\mu_{2},\theta_{2})$ corresponding to
the initial data $({\boldsymbol{v}}_{0}^{1},(\phi_{0}^{1},\psi_{0}^{1}))$,
$({\boldsymbol{v}}_{0}^{2},(\phi_{0}^{2},\psi_{0}^{2}))\in\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\times\mathcal{D}_{1}$,
the following continuous dependence estimate

$\begin{split}&\|{\boldsymbol{v}}_{1}(t)-{\boldsymbol{v}}_{2}(t)\|_{L^{2}(\Omega)}^{2}+\|\phi_{1}(t)-\phi_{2}(t)\|_{H^{1}(\Omega)}^{2}+\|\psi_{1}(t)-\psi_{2}(t)\|_{H^{1}(\Gamma)}^{2}\\\
&\leq
C\mathrm{e}^{C\int_{0}^{T}\mathcal{F}(t)\;\mathrm{d}t}\left(\|{\boldsymbol{v}}_{0}^{1}-{\boldsymbol{v}}_{0}^{2}\|_{L^{2}(\Omega)}^{2}+\|\phi_{0}^{1}-\phi_{0}^{2}\|_{H^{1}(\Omega)}^{2}+\|\psi_{0}^{1}-\psi_{0}^{2}\|_{H^{1}(\Gamma)}^{2}\right),\end{split}$
(3.21)

holds for all $t\in[0,T]$, where

$\mathcal{F}:=1+\|{\boldsymbol{v}}_{1}\|_{L^{4}(\Omega)}^{4}+\|{\boldsymbol{v}}_{2}\|_{H^{1}(\Omega)}^{2}+\|\mu_{2}\|_{H^{1}(\Omega)}^{2}+\|\theta_{2}\|_{H^{1}(\Gamma)}^{2}+\|\phi_{1}\|_{W^{1,4}(\Omega)}^{2},$

and the constant $C>0$ depends only on $\Omega$, the parameters of the system
and the norms of the initial data.

###  3.3 Proofs

  * Proof of Theorem 3.3.

We construct a weak solution to system (3.12) by discretizing the weak
formulation (3.16) through a Faedo–Galerkin scheme.

In order to obtain suitable uniform estimates for the Galerkin approximation,
it is necessary that the phase-fields $(\phi,\psi)$ and the chemical
potentials $(\mu,\theta)$ are approximated in the same finite-dimensional
subspace. In turn this is only possible if the phase-fields and the chemical
potentials are coupled by the same Dirichlet type boundary condition (cf.
(3.12g)), which is not the case for general $\beta>0$. However, by a change of
variables, we can reformulate system (3.12) in such a way that the Dirichlet
type boundary conditions for the redefined variables are exactly the same. To
this end, we set $\alpha:=\sqrt{\beta}>0$, and we introduce the functions

$\displaystyle\left\\{\;\begin{aligned}
\Psi&:=\alpha^{-1}\psi,&&\quad\text{i.e.,}\quad\psi=\alpha\Psi,\\\
\Psi_{0}&:=\alpha^{-1}\psi_{0},&&\quad\text{i.e.,}\quad\psi_{0}=\alpha\Psi_{0},\\\
\Theta&:=\alpha\theta,&&\quad\text{i.e.,}\quad\theta=\alpha^{-1}\Theta,\\\
n_{\Gamma}(\Psi)&:=\alpha^{-2}m_{\Gamma}(\alpha\Psi)=\alpha^{-2}m_{\Gamma}(\psi).\end{aligned}\right.$
(3.22)

By this change of variables, system (3.12) is equivalent to

$\displaystyle\partial_{t}{\boldsymbol{v}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}-\textnormal{div}\big{(}2\nu(\phi)\,\mathbf{D}{\boldsymbol{v}}\big{)}+\nabla\overline{p}=\mu\nabla\phi$
$\displaystyle\quad\text{in $Q$},$ (3.23a)
$\displaystyle\textnormal{div}\,{\boldsymbol{v}}=0,$
$\displaystyle\quad\text{in $Q$},$ (3.23b)
$\displaystyle\partial_{t}\phi+\textnormal{div}(\phi{\boldsymbol{v}})=\textnormal{div}\big{(}m_{\Omega}(\phi)\nabla\mu\big{)}$
$\displaystyle\quad\text{in $Q$},$ (3.23c)
$\displaystyle\mu=-\Delta\phi+F^{\prime}(\phi)$ $\displaystyle\quad\text{in
$Q$},$ (3.23d)
$\displaystyle\partial_{t}\Psi+\textnormal{div}_{\Gamma}(\Psi{\boldsymbol{v}}_{\tau})=\textnormal{div}_{\Gamma}\big{(}n_{\Gamma}(\Psi)\nabla_{\Gamma}\Theta\big{)}-\alpha\,m_{\Omega}(\alpha\Psi)\,\partial_{\mathbf{n}}\mu$
$\displaystyle\quad\text{on $\Sigma$},$ (3.23e)
$\displaystyle\Theta=-\alpha^{2}\,\Delta_{\Gamma}\Psi+\alpha
G^{\prime}(\alpha\Psi)+\alpha\partial_{\mathbf{n}}\phi$
$\displaystyle\quad\text{on $\Sigma$},$ (3.23f)
$\displaystyle{\boldsymbol{v}}\cdot\mathbf{n}=0\qquad\phi|_{\Gamma}=\alpha\Psi,\qquad\mu|_{\Gamma}=\alpha\Theta$
$\displaystyle\quad\text{on $\Sigma$},$ (3.23g)
$\displaystyle\big{[}2\nu(\alpha\Psi)(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+\gamma(\alpha\Psi){\boldsymbol{v}}\big{]}_{\tau}=\big{[}-\Psi\nabla_{\Gamma}\Theta\big{]}_{\tau}$
$\displaystyle\quad\text{on $\Sigma$},$ (3.23h)
$\displaystyle{\boldsymbol{v}}|_{t=0}={\boldsymbol{v}}_{0},\quad\phi|_{t=0}=\phi_{0}$
$\displaystyle\quad\text{in $\Omega$},$ (3.23i)
$\displaystyle\Psi|_{t=0}=\Psi_{0}$ $\displaystyle\quad\text{on $\Gamma$}.$
(3.23j)

In view of Definition 3.2, a weak solution
$({\boldsymbol{v}},\phi,\Psi,\mu,\Theta)$ of system (3.22) needs to have the
regularity

$\displaystyle{\boldsymbol{v}}\in
L^{\infty}\big{(}0,T;\mathbf{L}_{\textnormal{div}}^{2}(\Omega)\big{)}\cap
L^{2}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)\big{)},$ (3.24a)
$\displaystyle\left\\{\begin{aligned} &{\boldsymbol{v}}\in
W^{1,2}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime}\big{)}\cap
C\big{(}[0,T];\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\big{)}\ &&\text{if }\
d=2,\\\ &{\boldsymbol{v}}\in
W^{1,\frac{4}{3}}\big{(}0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime}\big{)}\cap
C_{w}\big{(}[0,T];\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\big{)}\ &&\text{if
}\ d=3,\end{aligned}\right.$ (3.24b) $\displaystyle(\phi,\Psi)\in
C([0,T];\mathcal{L}^{2})\cap H^{1}(0,T;\mathcal{D}_{\alpha})\cap
L^{\infty}(0,T;\mathcal{D}_{\alpha}),$ (3.24c) $\displaystyle(\mu,\Theta)\in
L^{4}(0,T;\mathcal{L}^{2})\cap L^{2}(0,T;\mathcal{D}_{\alpha}),$ (3.24d)

and satisfy the weak formulation which reads as

$\displaystyle\begin{aligned}
&\langle\partial_{t}{\boldsymbol{v}},{\boldsymbol{w}}\rangle_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}-\int_{\Omega}({\boldsymbol{v}}\otimes{\boldsymbol{v}}):\nabla{\boldsymbol{w}}\;\mathrm{d}x+\int_{\Omega}2\nu(\phi)\mathbf{D}{\boldsymbol{v}}:\mathbf{D}{\boldsymbol{w}}\;\mathrm{d}x\\\
&\quad=\int_{\Omega}\mu\nabla\phi\cdot{\boldsymbol{w}}\;\mathrm{d}x-\int_{\Gamma}\gamma(\alpha\Psi){\boldsymbol{v}}\cdot{\boldsymbol{w}}+\Psi\nabla_{\Gamma}\Theta\cdot{\boldsymbol{w}}\;\mathrm{d}S,\end{aligned}$
(3.25a) $\displaystyle\begin{aligned}
&\big{<}(\partial_{t}\phi,\partial_{t}\Psi),(\zeta,\xi)\big{>}_{\mathcal{D}_{\alpha}}-\int_{\Omega}\phi{\boldsymbol{v}}\cdot\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}\Psi{\boldsymbol{v}}\cdot\nabla_{\Gamma}\xi\;\mathrm{d}S\\\
&\quad=-\int_{\Omega}m_{\Omega}(\phi)\nabla\mu\cdot\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}\eta_{\Gamma}(\Psi)\nabla_{\Gamma}\Theta\cdot\nabla_{\Gamma}\xi\;\mathrm{d}S,\end{aligned}$
(3.25b) $\displaystyle\begin{aligned}
&\int_{\Omega}\mu\,\eta\;\mathrm{d}x+\int_{\Gamma}\Theta\,\vartheta\;\mathrm{d}S\\\
&\quad=\int_{\Omega}\nabla\phi\cdot\nabla\eta+F^{\prime}(\phi)\eta\;\mathrm{d}x+\int_{\Gamma}\alpha^{2}\nabla_{\Gamma}\Psi\cdot\nabla_{\Gamma}\vartheta+\alpha
G^{\prime}(\alpha\Psi)\vartheta\;\mathrm{d}S\end{aligned}$ (3.25c)

for all test functions
${\boldsymbol{w}}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ and
$(\zeta,\xi),\,(\eta,\vartheta)\in\mathcal{D}_{\alpha}$.

We now prove the existence of a weak solution of system (3.23) (in the
aforementioned sense) and afterwards, we use the change of variables (3.22) to
infer the existence of a weak solution to the original model (3.12). The weak
formulation (3.25) can indeed be discretized by a suitable Faedo–Galerkin
scheme.

Step 1: Discretization via a bulk-surface Faedo–Galerkin scheme. It has
already been established in [40, Theorem 4.4] that the second-order elliptic
eigenvalue problem

$\displaystyle\left\\{\begin{aligned} -\Delta\zeta&=\lambda\zeta&&\text{in
$\Omega$},\\\
-\Delta_{\Gamma}\xi+\alpha\partial_{\mathbf{n}}\zeta&=\lambda\xi&&\text{on
$\Gamma$},\\\ \zeta|_{\Gamma}&=\alpha\xi&&\text{on
$\Gamma$}\end{aligned}\right.$ (3.26)

has countably many eigenvalues $\\{\lambda_{i}\\}_{i\in\mathbb{N}}$ with
$\lambda_{i}\geq 0$ for all $i\in\mathbb{N}$, and the corresponding
eigenfunctions
$\\{(\zeta_{k},\xi_{k})\\}_{k\in\mathbb{N}}\subset\mathcal{D}_{\alpha}$ can be
chosen such that they form an orthonormal basis of $\mathcal{L}^{2}$. Here, we
fix the eigenfunction associated with the first eigenvalue $\lambda_{1}=0$ as

$\displaystyle(\zeta_{1},\xi_{1})=\frac{1}{\sqrt{\alpha^{2}|\Omega|+|\Gamma|}}(\alpha,1).$
(3.27)

Moreover, in view of their construction (see [40]), the eigenfunctions
$(\zeta_{k},\xi_{k})$ with $k\geq 2$ satisfy the mean value constraint

$\displaystyle\alpha|\Omega|{\langle\zeta_{k}\rangle}_{\Omega}+|\gamma|{\langle\xi_{k}\rangle}_{\Gamma}=0.$
(3.28)

For any $k\in\mathbb{N}$, we introduce the finite dimensional subspace

$\displaystyle\mathcal{Z}_{k}:=\text{span}\,\\{(\zeta_{1},\xi_{1}),...,(\zeta_{k},\xi_{k})\\}\subset{\mathcal{D}_{\alpha}},$
(3.29)

and we write

$\displaystyle\mathbb{P}_{\mathcal{Z}_{k}}(\zeta,\xi)=\big{(}\mathbb{P}^{\Omega}_{\mathcal{Z}_{k}}(\zeta,\xi),\mathbb{P}^{\Gamma}_{\mathcal{Z}_{k}}(\zeta,\xi)\big{)}\quad\text{for
any $(\zeta,\xi)\in\mathcal{D}_{\alpha}$},$

to denote the $\mathcal{L}^{2}$-orthogonal projection onto $\mathcal{Z}_{k}$.

Next, we consider the Stokes operator with Navier boundary conditions (see,
e.g., [2, Appendix A])

$\mathbf{A}:\mathcal{D}(\mathbf{A})\subset\mathbf{L}^{2}_{\textnormal{div}}(\Omega)\rightarrow\mathbf{L}^{2}_{\textnormal{div}}(\Omega):\quad{\boldsymbol{v}}\mapsto\mathbf{A}{\boldsymbol{v}}:=-\mathbb{P}\textnormal{div}(2D{\boldsymbol{v}}),$
(3.30)

where $\mathbb{P}$ is the Leray projection operator, with domain

$\mathcal{D}(\mathbf{A})=\Big{\\{}{\boldsymbol{v}}\in\mathbf{H}^{2}(\Omega)\cap\mathbf{L}^{2}_{\textnormal{div}}(\Omega):\big{[}2(\mathbf{D}{\boldsymbol{v}}\,\mathbf{n})+{\boldsymbol{v}}\big{]}_{\tau}=0\;\text{on
$\Gamma$}\Big{\\}}.$ (3.31)

Since $\mathbf{A}$ is a positive and self-adjoint operator on
$\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$ with compact inverse, and

$\left(\mathbf{L}^{2}_{\textnormal{div}}(\Omega),\mathcal{D}(\mathbf{A})\right)_{\frac{1}{2},2}=\mathbf{H}^{1}_{\textnormal{div}}(\Omega),$

there exists a sequence of countably many positive eigenvalues
$\big{\\{}\widetilde{\lambda}_{i}\big{\\}}_{i\in\mathbb{N}}$ and the
corresponding eigenfunctions
$\\{{\boldsymbol{w}}_{i}\\}_{i\in\mathbb{N}}\subset\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$,
which are determined by

$\int_{\Omega}2D{\boldsymbol{w}}_{i}:\nabla{\boldsymbol{v}}\;\mathrm{d}x+\int_{\Gamma}{\boldsymbol{w}}_{i}\cdot{\boldsymbol{v}}\;\mathrm{d}S=\widetilde{\lambda}_{i}\int_{\Omega}{\boldsymbol{w}}_{i}\cdot{\boldsymbol{v}}\;\mathrm{d}x\quad\text{for
all ${\boldsymbol{v}}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$},$ (3.32)

can be chosen in such a way that they form an orthonormal basis of
$\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$. For any $k\in\mathbb{N}$, we
define the finite dimensional subspace

$\mathcal{U}_{k}=\text{span}\,\\{{\boldsymbol{w}}_{1},...,{\boldsymbol{w}}_{k}\\}\subset\mathbf{H}^{1}_{\textnormal{div}}(\Omega),$

and the $\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$-orthogonal projection onto
$\mathcal{U}_{k}$ is denoted by $\mathbb{P}_{\mathcal{U}_{k}}$.

Next, for any $k\in\mathbb{N}$ and $t\in[0,T]$, we make the ansatz

$\displaystyle{\boldsymbol{v}}_{k}(t)$
$\displaystyle:=\sum_{i=1}^{k}a^{k}_{i}(t)\,{\boldsymbol{w}}_{i},$
$\displaystyle\quad\text{a.e. in $\Omega$},$ (3.33a)
$\displaystyle\big{(}\phi_{k}(t),\Psi_{k}(t)\big{)}$
$\displaystyle:=\sum_{i=1}^{k}b^{k}_{i}(t)\,(\zeta_{i},\xi_{i}),$
$\displaystyle\quad\text{a.e. in $\Omega\times\Gamma$},$ (3.33b)
$\displaystyle\big{(}\mu_{k}(t),\Theta_{k}(t)\big{)}$
$\displaystyle:=\sum_{i=1}^{k}c^{k}_{i}(t)\,(\zeta_{i},\xi_{i}),$
$\displaystyle\quad\text{a.e. in $\Omega\times\Gamma$}.$ (3.33c)

Here, the scalar, time-dependent coefficients $a^{k}_{i}$, $b^{k}_{i}$ and
$c^{k}_{i}$, $i=1,...,k$ are assumed to be continuously differentiable
functions that are still to be determined. They need to be designed in such a
way that the discretized weak formulation

$\displaystyle\begin{aligned}
&\langle\partial_{t}{\boldsymbol{v}}_{k},{\boldsymbol{w}}\rangle_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}-\int_{\Omega}({\boldsymbol{v}}_{k}\otimes{\boldsymbol{v}}_{k})\nabla{\boldsymbol{w}}\;\mathrm{d}x+\int_{\Omega}2\nu(\phi_{k})\mathbf{D}{\boldsymbol{v}}_{k}:\mathbf{D}{\boldsymbol{w}}\;\mathrm{d}x\\\
&\quad=\int_{\Omega}\mu_{k}\nabla\phi_{k}\cdot{\boldsymbol{w}}\;\mathrm{d}x-\int_{\Gamma}\gamma(\alpha\Psi_{k}){\boldsymbol{v}}_{k}\cdot{\boldsymbol{w}}+\Psi_{k}\nabla_{\Gamma}\Theta_{k}\cdot{\boldsymbol{w}}\;\mathrm{d}S,\end{aligned}$
(3.34a) $\displaystyle\begin{aligned}
&\big{<}(\partial_{t}\phi_{k},\partial_{t}\Psi_{k}),(\zeta,\xi)\big{>}_{{\mathcal{D}_{\alpha}}}-\int_{\Omega}\phi_{k}{\boldsymbol{v}}_{k}\cdot\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}\Psi_{k}{\boldsymbol{v}}_{k}\cdot\nabla_{\Gamma}\xi\;\mathrm{d}S\\\
&\quad=-\int_{\Omega}m_{\Omega}(\phi_{k})\nabla\mu_{k}:\nabla\zeta\;\mathrm{d}x-\int_{\Gamma}n_{\Gamma}(\Psi_{k})\nabla_{\Gamma}\Theta_{k}:\nabla_{\Gamma}\xi\;\mathrm{d}S,\end{aligned}$
(3.34b) $\displaystyle\begin{aligned}
&\int_{\Omega}\mu_{k}\cdot\zeta\;\mathrm{d}x+\int_{\Gamma}\Theta_{k}\cdot\xi\;\mathrm{d}S\\\
&\quad=\int_{\Omega}\nabla\phi_{k}\cdot\nabla\zeta+F^{\prime}(\phi_{k})\zeta\;\mathrm{d}x+\int_{\Gamma}\alpha^{2}\nabla_{\Gamma}\Psi_{k}\cdot\nabla_{\Gamma}\xi+\alpha
G^{\prime}(\alpha\Psi_{k})\xi\;\mathrm{d}S\end{aligned}$ (3.34c)

holds for all ${\boldsymbol{w}}\in\mathcal{U}_{k}$ and all
$(\zeta,\xi)\in\mathcal{Z}_{k}$, supplemented with the initial conditions

$\displaystyle{\boldsymbol{v}}_{k}|_{t=0}$
$\displaystyle=\mathbb{P}_{\mathcal{U}_{k}}({\boldsymbol{v}}_{0})\in\mathcal{U}_{k}$
in $\Omega$, (3.35) $\displaystyle(\phi_{k},\Psi_{k})|_{t=0}$
$\displaystyle=\mathbb{P}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\in\mathcal{Z}_{k}$
$\displaystyle\quad\text{in $\Omega\times\Gamma$},$ (3.36)

where $\Psi_{0}:=\alpha^{-1}\psi_{0}$.

Let now $\mathbf{a}^{k}:=(a_{1}^{k},...,a_{k}^{k})$,
$\mathbf{b}^{k}:=(b_{1}^{k},...,b_{k}^{k})$ and
$\mathbf{c}^{k}:=(c_{1}^{k},...,c_{k}^{k})$ denote the coefficient vectors.
Testing (3.34a) with ${\boldsymbol{w}}_{1},...,{\boldsymbol{w}}_{k}$ and
(3.34b) with $(\zeta_{1},\xi_{1}),...,(\zeta_{k},\xi_{k})$, we conclude that
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ is determined by a system of $2k$
ordinary differential equations whose right-hand side depends continuously on
$\mathbf{a}^{k}$, $\mathbf{b}^{k}$ and $\mathbf{c}^{k}$. Due to (3.35) and
(3.36), this ODE system is subject to the initial condition

$\displaystyle[\mathbf{a}^{k}]_{i}(0)$
$\displaystyle=a_{i}^{k}(0)=\big{(}{\boldsymbol{v}}_{0},{\boldsymbol{w}}_{i}\big{)}_{\mathbf{L}^{2}(\Omega)},$
$\displaystyle\quad i\in\\{1,...,k\\},$ $\displaystyle[\mathbf{b}^{k}]_{i}(0)$
$\displaystyle=b_{i}^{k}(0)=\big{(}(\phi_{0},\Psi_{0}),(\zeta_{i},\xi_{i})\big{)}_{\mathcal{L}^{2}},$
$\displaystyle\quad i\in\\{1,...,k\\}.$

Moreover, testing (3.34c) with $(\zeta_{1},\xi_{1}),...,(\zeta_{k},\xi_{k})$,
we infer that $\mathbf{c}^{k}$ is explicitly given by a system of $k$
algebraic equations whose right-hand side depends continuously on
$\mathbf{b}^{k}$. We now replace the variable $\mathbf{c}^{k}$ appearing in
the right-hand side of the aforementioned ODE system to obtain a closed
$2k$-dimensional ODE system describing
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ whose right-hand side depends
continuously on $(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$. We can thus apply
the Cauchy–Peano theorem to infer the existence of at least one local solution
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}:[0,T_{k}^{*})\cap[0,T]\to\mathbb{R}^{2k}$
with $T_{k}^{*}>0$ to the corresponding initial value problem. Without loss of
generality, we assume that $T_{k}^{*}\leq T$ and that
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ is right-maximal, meaning that
$T_{k}^{*}$ is chosen as large as possible. Now, we can reconstruct
$\mathbf{c}^{k}:[0,T_{k}^{*})\to\mathbb{R}^{k}$ by the aforementioned system
of $k$ algebraic equations. In view of (3.33), we thus obtain functions

$\displaystyle{\boldsymbol{v}}_{k}\in
C^{1}\big{(}[0,T_{k}^{*});\mathbf{H}^{1}_{\textnormal{div}}(\Omega)\big{)},\quad(\phi_{k},\Psi_{k})\in
C^{1}\big{(}[0,T_{k}^{*});{\mathcal{D}_{\alpha}}\big{)},\quad(\mu_{k},\Theta_{k})\in
C^{1}\big{(}[0,T_{k}^{*});{\mathcal{D}_{\alpha}}\big{)},$

which fulfill the discretized weak formulation (3.34) on the time interval
$[0,T_{k}^{*})$.

Step 2: Uniform estimates. We establish suitable estimates for each
approximate solution
$({\boldsymbol{v}}_{k},\phi_{k},\psi_{k},\mu_{k},\theta_{k})$ which are
uniform with respect to the index $k$. In particular, let $T_{k}<T_{k}^{*}$ be
arbitrary. In the following, let $C$ denote generic positive constants
depending only on the initial data and the constants introduced in
($\mathrm{A2}$)–($\mathrm{A5}$) including the final time $T$, but is
independent of $k$ and $T_{k}$.

Testing (3.34a) with ${\boldsymbol{v}}_{k}$, (3.34b) with
$(\mu_{k},\theta_{k})$, and (3.34c) with
$-(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})$, adding the resulting
equations, and integrating with respect to time from $0$ to $t$, we derive the
discretized energy identity

$\displaystyle\begin{aligned}
&E\big{(}{\boldsymbol{v}}_{k}(t),\phi_{k}(t),\alpha\Psi_{k}(t))+\int_{0}^{t}\int_{\Omega}2\nu(\phi_{k})\,\left|\mathbf{D}{\boldsymbol{v}}_{k}\right|^{2}\;\mathrm{d}x\;\mathrm{d}t+\int_{0}^{t}\int_{\Gamma}\gamma(\alpha\Psi_{k})\left|{\boldsymbol{v}}_{k}\right|^{2}\;\mathrm{d}S\;\mathrm{d}t\\\
&\qquad+\int_{0}^{t}\int_{\Omega}m_{\Omega}(\phi_{k})\left|\nabla\mu_{k}\right|^{2}\;\mathrm{d}x\;\mathrm{d}t+\int_{0}^{t}\int_{\Gamma}n_{\Gamma}(\psi_{k})\left|\nabla_{\Gamma}\Theta_{k}\right|^{2}\;\mathrm{d}S\;\mathrm{d}t\\\\[4.30554pt]
&=E\big{(}{\boldsymbol{v}}_{k}(0),\phi_{k}(0),\alpha\Psi_{k}(0))\end{aligned}$
(3.37)

for all $t\in[0,T_{k}]$. Here, $E$ stands for the total energy introduced in
(3.13). Due to the growth conditions on $F$ and $G$ (see ($\mathrm{A3}$)), the
Sobolev embeddings $H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$ and
$H^{1}(\Gamma)\hookrightarrow L^{p}(\Gamma)$ for any $1\leq p<\infty$, and the
properties of the projections $\mathbb{P}_{\mathcal{U}_{k}}$ and
$\mathbb{P}_{\mathcal{Z}_{k}}$, we use the initial conditions (3.35) and
(3.36) to infer that

$\displaystyle E\big{(}{\boldsymbol{v}}_{k}(0),\phi_{k}(0),\alpha\Psi_{k}(0))$
$\displaystyle\leq\frac{1}{2}\|\mathbb{P}_{\mathcal{U}_{k}}({\boldsymbol{v}}_{0})\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\|\nabla\mathbb{P}^{\Omega}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{L^{2}(\Omega)}^{2}+\frac{\alpha^{2}}{2}\|\nabla\mathbb{P}^{\Gamma}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle\quad+c_{F}\left(|\Omega|+\|\mathbb{P}^{\Omega}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{L^{p}(\Omega)}^{p}\right)+c_{G}\left(|\Gamma|+\|\mathbb{P}^{\Gamma}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{L^{q}(\Gamma)}^{q}\right)$
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{v}}_{0}\|_{L^{2}(\Omega)}^{2}+C\|(\phi_{0},\Psi_{0})\|_{\mathcal{H}^{1}}^{2}$
$\displaystyle\quad+C\left(1+\|\mathbb{P}^{\Omega}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{H^{1}(\Omega)}^{p}\right)+C\left(1+\|\mathbb{P}^{\Gamma}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\|_{H^{1}(\Gamma)}^{q}\right)$
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{v}}_{0}\|_{L^{2}(\Omega)}^{2}+C\|(\phi_{0},\Psi_{0})\|_{\mathcal{H}^{1}}^{2}$
$\displaystyle\quad+C\left(1+\|(\phi_{0},\Psi)\|_{\mathcal{H}^{1}}^{p}\right)+C\left(1+\|(\phi_{0},\Psi_{0})\|_{\mathcal{H}^{1}}^{q}\right),$

which, in turn, entails that

$\displaystyle
E\big{(}{\boldsymbol{v}}_{k}(0),\phi_{k}(0),\alpha\Psi_{k}(0))\leq C.$ (3.38)

Recalling that $F,G\geq 0$ and that the functions $\nu$, $\gamma$,
$m_{\Omega}$ and $m_{\Gamma}$ (and thus also $n_{\Gamma}$) are uniformly
positive (see ($\mathrm{A4}$)–($\mathrm{A5}$)), we conclude the uniform bounds

$\displaystyle\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Omega))}+\|\nabla{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};L^{2}(\Omega))}+\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};L^{2}(\Gamma))}$
$\displaystyle\leq C,$ (3.39)
$\displaystyle\|\nabla\phi_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Omega))}+\|\nabla_{\Gamma}\Psi_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Gamma))}$
$\displaystyle\leq C,$ (3.40)
$\displaystyle\|\nabla\mu_{k}\|_{L^{2}(0,T_{k};L^{2}(\Omega))}+\|\nabla_{\Gamma}\Theta_{k}\|_{L^{2}(0,T_{k};L^{2}(\Gamma))}$
$\displaystyle\leq C.$ (3.41)

Testing (3.34a) with the first eigenfunction $(\zeta_{1},\xi_{1})$ (see
(3.27)), we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left(\alpha\int_{\Omega}\phi_{k}(t)\;\mathrm{d}x+\int_{\Gamma}\Psi_{k}(t)\;\mathrm{d}S\right)=0,$

which implies

$\alpha\int_{\Omega}\phi_{k}(t)\;\mathrm{d}x+\int_{\Gamma}\Psi_{k}(t)\;\mathrm{d}S=\alpha\int_{\Omega}\phi_{0}\;\mathrm{d}x+\int_{\Gamma}\Psi_{0}\;\mathrm{d}S,$
(3.42)

for all $t\in[0,T_{k}]$. Here, we have used that

$\displaystyle\alpha\left|\Omega\right|\langle\mathbb{P}^{\Omega}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\rangle_{\Omega}+\left|\Gamma\right|\langle\mathbb{P}^{\Gamma}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0})\rangle_{\Gamma}=\big{(}\mathbb{P}_{\mathcal{Z}_{k}}(\phi_{0},\Psi_{0}),(\alpha,1)\big{)}_{\mathcal{L}^{2}}$
$\displaystyle\quad=\big{(}(\phi_{0},\Psi_{0}),(\alpha,1)\big{)}_{\mathcal{L}^{2}}=\alpha\left|\Omega\right|\langle\phi_{0}\rangle_{\Omega}+\left|\Gamma\right|\langle\Psi_{0}\rangle_{\Gamma},$

which follows by means of (3.27) and (3.28). Defining the constant

$c_{m}:=\frac{\alpha\left|\Omega\right|\langle\phi_{0}\rangle_{\Omega}+\left|\Gamma\right|\langle\Psi_{0}\rangle_{\Gamma}}{\alpha^{2}\left|\Omega\right|+\left|\Gamma\right|},$

we have $(\alpha c_{m},c_{m})\in{\mathcal{D}_{\alpha}}$ and invoking (3.42), a
straightforward computation yields

$\displaystyle\big{(}\phi_{k}(t)-\alpha
c_{m},\Psi_{k}(t)-c_{m}\big{)}\in{\mathcal{D}_{\alpha}}\quad\text{and}\quad\alpha\left|\Omega\right|\langle\phi_{k}(t)-\alpha
c_{m}\rangle_{\Omega}+\left|\Gamma\right|\langle\Psi_{k}(t)-c_{m}\rangle_{\Gamma}=0$

for all $t\in[0,T_{k}]$. Hence, applying the bulk-surface variant of the
Poincaré inequality (3.9) (with $\beta=\alpha^{2}$), we deduce that

$\displaystyle\|(\phi_{k},\Psi_{k})\|_{\mathcal{L}^{2}}$
$\displaystyle\leq\|(\phi_{k}-\alpha
c_{m},\Psi_{k}-c_{m})\|_{\mathcal{L}^{2}}+\|(\alpha
c_{m},c_{m})\|_{\mathcal{L}^{2}}$ $\displaystyle\leq
C_{P}\|(\nabla\phi_{k},\nabla_{\Gamma}\Psi_{k})\|_{\mathcal{L}^{2}}+C$

for all $t\in[0,T_{k}]$. In view of (3.40), we thus conclude

$\displaystyle\|\phi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Omega))}+\|\Psi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Gamma))}$
$\displaystyle\leq C.$ (3.43)

Next, we derive a uniform estimate for $(\mu_{k},\Theta_{k})$ in the full
$\mathcal{H}^{1}$-norm. For this aim, we test (3.34c) by
$(\zeta,\xi):=\mathbb{P}_{\mathcal{Z}_{k}}(\zeta^{*},\xi^{*})$, where
$(\zeta^{*},\xi^{*})\in{\mathcal{D}_{\alpha}}$ is an arbitrary test function.
By the growth conditions (3.3)–(3.4) from ($\mathrm{A3}$) and the Sobolev
embeddings $H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$ and
$H^{1}(\Gamma)\hookrightarrow L^{p}(\Gamma)$ for any $1\leq p<\infty$, we have

$\displaystyle\begin{aligned}
\left|\big{<}(\mu_{k},\Theta_{k}),(\zeta^{*},\xi^{*})\big{>}_{{\mathcal{D}_{\alpha}}}\right|&=\left|\big{<}(\mu_{k},\Theta_{k}),(\zeta,\xi)\big{>}_{{\mathcal{D}_{\alpha}}}\right|\\\\[4.30554pt]
&\leq\|\nabla\phi_{k}\|_{L^{2}(\Omega)}\|\nabla\zeta\|_{L^{2}(\Omega)}+\|F^{\prime}(\phi_{k})\|_{L^{\frac{6}{5}}(\Omega)}\|\zeta\|_{L^{6}(\Omega)}\\\
&\quad+\alpha^{2}\|\nabla_{\Gamma}\Psi_{k}\|_{L^{2}(\Gamma)}\|\nabla_{\Gamma}\xi\|_{L^{2}(\Gamma)}+\alpha\|G^{\prime}(\alpha\Psi_{k})\|_{L^{2}(\Gamma)}\|\xi\|_{L^{2}(\Gamma)}\\\
&\leq
C\left(1+\|\nabla\phi_{k}\|_{L^{2}(\Omega)}+\|\phi_{k}\|_{L^{6}(\Omega)}^{5}\right)\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}\\\
&\quad+C\left(1+\|\nabla_{\Gamma}\Psi_{k}\|_{L^{2}(\Gamma)}+\|\Psi_{k}\|_{L^{2q}(\Gamma)}^{q}\right)\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}\\\
&\leq
C\left(1+\|\phi_{k}\|_{H^{1}(\Omega)}^{5}+\|\Psi_{k}\|_{H^{1}(\Gamma)}^{q}\right)\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}},\end{aligned}$
(3.44)

for all $t\in[0,T_{k}]$. Taking the supremum over all
$(\zeta^{*},\xi^{*})\in{\mathcal{D}_{\alpha}}$ with
$\|(\zeta^{*},\xi^{*})\|_{{\mathcal{D}_{\alpha}}}\leq 1$, and exploiting
(3.43), we find that

$\|(\mu_{k},\Theta_{k})\|_{L^{\infty}(0,T_{k};\mathcal{D}_{\alpha}^{\prime})}\leq
C.$ (3.45)

By a duality argument, we infer the existence of a positive constant
$C_{{\mathcal{D}_{\alpha}}}$ such that

$\|(\bar{\zeta},\bar{\xi})\|_{\mathcal{L}^{2}}^{2}=\big{<}(\bar{\zeta},\bar{\xi}),(\bar{\zeta},\bar{\xi})\big{>}_{{\mathcal{D}_{\alpha}}}\leq
C_{{\mathcal{D}_{\alpha}}}\|(\bar{\zeta},\bar{\xi})\|_{\mathcal{D}_{\alpha}^{\prime}}\|(\nabla\bar{\zeta},\nabla_{\Gamma}\bar{\xi})\|_{\mathcal{L}^{2}}+C_{{\mathcal{D}_{\alpha}}}\|(\bar{\zeta},\bar{\xi})\|_{\mathcal{D}_{\alpha}^{\prime}}^{2}$
(3.46)

for all $(\bar{\zeta},\bar{\xi})\in{\mathcal{D}_{\alpha}}$. Hence, applying
(3.46) with $(\bar{\zeta},\bar{\xi})=(\mu_{k},\Theta_{k})$ and using (3.41)
and (3.45), we arrive at

$\|(\mu_{k},\Theta_{k})\|_{L^{4}(0,T_{k};\mathcal{L}^{2})}\leq C.$ (3.47)

In addition, we infer from (3.41) and (3.47) that

$\displaystyle\|(\mu_{k},\Theta_{k})\|_{L^{2}(0,T_{k};\mathcal{H}^{1})}$
$\displaystyle\leq C.$ (3.48)

Lastly, we derive uniform estimates on the time derivatives. Therefore, let
${\boldsymbol{w}}^{*}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ and
$(\zeta^{*},\xi^{*})\in{\mathcal{D}_{\alpha}}$ be arbitrary test functions. We
now test (3.34a) with
${\boldsymbol{w}}:=\mathbb{P}_{\mathcal{U}_{k}}({\boldsymbol{w}}^{*})$.

In three dimensions, recalling that the functions $\nu$ and $\gamma$ are
bounded, we use Hölder’s inequality, the continuous embeddings
$H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$ and $H^{1}(\Omega)\hookrightarrow
L^{4}(\Gamma)$ as well as the Gagliardo-Nirenberg interpolation inequality to
obtain the estimate

$\displaystyle\begin{aligned}
\left|\big{<}\partial_{t}{\boldsymbol{v}}_{k},{\boldsymbol{w}}^{*}\big{>}_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}\right|&=\left|\big{<}\partial_{t}{\boldsymbol{v}}_{k},{\boldsymbol{w}}\big{>}_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}\right|\\\\[4.30554pt]
&\leq\|{\boldsymbol{v}}_{k}\|_{L^{4}(\Omega)}^{2}\|\nabla{\boldsymbol{w}}\|_{L^{2}(\Omega)}+C\|\nabla{\boldsymbol{v}}_{k}\|_{L^{2}(\Omega)}\|\nabla{\boldsymbol{w}}\|_{L^{2}(\Omega)}\\\
&\qquad+\|\mu_{k}\|_{L^{6}(\Omega)}\|\nabla\phi_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{w}}\|_{L^{3}(\Omega)}+C\|{\boldsymbol{v}}_{k}\|_{L^{2}(\Gamma)}\|{\boldsymbol{w}}\|_{L^{2}(\Gamma)}\\\
&\qquad+\|\Psi_{k}\|_{L^{4}(\Gamma)}\|\nabla\Theta_{k}\|_{L^{2}(\Gamma)}\|{\boldsymbol{w}}\|_{L^{4}(\Gamma)}\\\\[4.30554pt]
&\leq
C\|{\boldsymbol{v}}_{k}\|_{L^{2}(\Omega)}^{\frac{1}{2}}\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}^{\frac{3}{2}}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}+C\|\nabla{\boldsymbol{v}}_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}\\\
&\qquad+C\|\mu_{k}\|_{H^{1}(\Omega)}\|\nabla\phi_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}+C\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}\\\
&\qquad+\|\Psi_{k}\|_{H^{1}(\Gamma)}\|\Theta_{k}\|_{H^{1}(\Gamma)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}\end{aligned}$

for all $t\in[0,T_{k}]$. Taking the supremum over all
${\boldsymbol{w}}^{*}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ with
$\|{\boldsymbol{w}}^{*}\|_{\mathbf{H}^{1}_{\textnormal{div}}(\Omega)}\leq 1$,
taking the power $\frac{4}{3}$ and integrating over $[0,T_{k}]$, we now deduce
that

$\displaystyle\begin{aligned}
\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{\frac{4}{3}}(0,T_{k};\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}^{\frac{4}{3}}&\leq
C\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Omega))}^{\frac{2}{3}}\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{2}\\\
&\quad+C\,T_{k}^{\frac{1}{3}}\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{\frac{4}{3}}\\\
&\quad+C\,T_{k}^{\frac{1}{3}}\|\phi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Omega))}^{\frac{4}{3}}\|\mu_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{\frac{4}{3}}\\\
&\quad+C\,T_{k}^{\frac{1}{3}}\|\Psi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Gamma))}^{\frac{4}{3}}\|\Theta_{k}\|_{L^{2}(0,T_{k};H^{1}(\Gamma))}^{\frac{4}{3}}.\end{aligned}$

By exploiting (3.39), (3.43) and (3.48), and recalling that $T_{k}\leq T$, we
simply reach

$\displaystyle\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{\frac{4}{3}}(0,T_{k};\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}\leq
C,\quad\text{if $d=3$}.$ (3.49)

In two dimensions, arguing as above and using the Ladyzhenskaya inequality for
the convective term, we first obtain

$\displaystyle\begin{aligned}
\left|\big{<}\partial_{t}{\boldsymbol{v}}_{k},{\boldsymbol{w}}^{*}\big{>}_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}\right|&=\left|\big{<}\partial_{t}{\boldsymbol{v}}_{k},{\boldsymbol{w}}\big{>}_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}\right|\\\\[4.30554pt]
&\leq
C\|{\boldsymbol{v}}_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}+C\|\nabla{\boldsymbol{v}}_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}\\\
&\qquad+C\|\mu_{k}\|_{H^{1}(\Omega)}\|\nabla\phi_{k}\|_{L^{2}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}+C\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}\\\
&\qquad+\|\Psi_{k}\|_{H^{1}(\Gamma)}\|\Theta_{k}\|_{H^{1}(\Gamma)}\|{\boldsymbol{w}}^{*}\|_{H^{1}(\Omega)}.\end{aligned}$

Thus, after taking the supremum over all
${\boldsymbol{w}}^{*}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ with
$\|{\boldsymbol{w}}^{*}\|_{\mathbf{H}^{1}_{\textnormal{div}}(\Omega)}\leq 1$,
taking the square and integrating over $[0,T_{k}]$, we find

$\displaystyle\begin{aligned}
\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}^{2}&\leq
C\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Omega))}^{2}\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{2}+C\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{2}\\\
&\quad+C\|\phi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Omega))}^{2}\|\mu_{k}\|_{L^{2}(0,T_{k};H^{1}(\Omega))}^{2}\\\
&\quad+C\|\Psi_{k}\|_{L^{\infty}(0,T_{k};H^{1}(\Gamma))}^{2}\|\Theta_{k}\|_{L^{2}(0,T_{k};H^{1}(\Gamma))}^{2},\end{aligned}$

which immediately entails that

$\displaystyle\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{2}(0,T_{k};\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}\leq
C,\quad\text{if $d=2$}.$ (3.50)

Concerning the bulk and surface phase-fields, we test (3.34b) with
$(\zeta,\xi):=\mathbb{P}_{\mathcal{Z}_{k}}(\zeta^{*},\xi^{*})$. Recalling that
the functions $m_{\Omega}$ and $m_{\Gamma}$ (and thus also $n_{\Gamma}$) are
bounded, we use Hölder’s inequality along with the continuous embeddings
$H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$, $H^{1}(\Gamma)\hookrightarrow
L^{4}(\Gamma)$ and $H^{1}(\Omega)\hookrightarrow L^{4}(\Gamma)$ to infer that

$\displaystyle\begin{aligned}
\left|\big{<}(\partial_{t}\phi_{k},\partial_{t}\Psi_{k}),(\zeta^{*},\xi^{*})\big{>}_{{\mathcal{D}_{\alpha}}}\right|&=\left|\big{<}(\partial_{t}\phi_{k},\partial_{t}\Psi_{k}),(\zeta,\xi)\big{>}_{{\mathcal{D}_{\alpha}}}\right|\\\\[4.30554pt]
&\leq\|\phi_{k}\|_{H^{1}(\Omega)}\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}+C\|\mu_{k}\|_{H^{1}(\Omega)}\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}\\\
&\quad+\|\Psi_{k}\|_{H^{1}(\Gamma)}\|{\boldsymbol{v}}_{k}\|_{H^{1}(\Omega)}\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}+C\|\Theta_{k}\|_{H^{1}(\Gamma)}\|(\zeta^{*},\xi^{*})\|_{\mathcal{H}^{1}}.\end{aligned}$

Hence, after taking the supremum over all
$(\zeta^{*},\xi^{*})\in{\mathcal{D}_{\alpha}}$ with
$\|(\zeta^{*},\xi^{*})\|_{{\mathcal{D}_{\alpha}}}\leq 1$, we take the square
of the resulting inequality and integrate over $[0,T_{k}]$. This yields

$\displaystyle\|(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})\|_{L^{2}(0,T_{k};\mathcal{D}_{\alpha}^{\prime})}^{2}$
$\displaystyle\leq\|\phi_{k}\|_{L^{\infty}(0,T;H^{1}(\Omega))}^{2}\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+C\|\mu_{k}\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}$
$\displaystyle\qquad+\|\Psi_{k}\|_{L^{\infty}(0,T;H^{1}(\Gamma))}^{2}\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+C\|\Theta_{k}\|_{L^{2}(0,T;H^{1}(\Gamma))}^{2}.$

In view of the uniform estimates (3.39), (3.43) and (3.48), we thus conclude

$\displaystyle\|(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})\|_{L^{2}(0,T_{k};\mathcal{D}_{\alpha}^{\prime})}\leq
C.$ (3.51)

Step 3: Extension of the approximate solution onto $[0,T]$. We recall from
Step 1 that the coefficient vector $(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ is
determined as a solution of a nonlinear ODE system. Using the definition of
the approximate solutions given in (3.33) as well as the uniform a priori
estimates (3.39) and (3.43), we obtain that for any $T_{k}\in[0,T_{k}^{*})$,
all $t\in[0,T_{k}]$, and all $i\in\\{1,...,k\\}$,

$\displaystyle|a_{i}^{k}(t)|+|b_{j}^{k}(t)|$
$\displaystyle=\big{|}\left({\boldsymbol{v}}_{k}(t),{\boldsymbol{w}}_{i}\right)_{L^{2}(\Omega)}\big{|}+\big{|}\left(\big{(}\phi_{k}(t),\Psi_{k}(t)\big{)},(\zeta_{i},\xi_{i})\right)_{\mathcal{L}^{2}}\big{|}$
$\displaystyle\leq\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T_{k};L^{2}(\Omega))}+\|(\phi_{k},\Psi_{k})\|_{L^{\infty}(0,T_{k};\mathcal{L}^{2})}\leq
C.$

This means that the solution $(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ is
bounded on the time interval $[0,T_{k}^{*})$ by a constant that does neither
depend on $T_{k}$ nor on $k$. Consequently, by the classical ODE theory, the
solution can thus be extended beyond the time $T_{k}^{*}$. However, as
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ was assumed to be a right-maximal
solution, meaning that $T_{k}$ was assumed to be chosen as large as possible,
this is a contradiction. We thus conclude that the solution
$(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ actually exists on the whole time
interval $[0,T]$. As the coefficients $\mathbf{c}^{k}$ can be reconstructed
from $(\mathbf{a}^{k},\mathbf{b}^{k})^{\top}$ via the vector-valued algebraic
equation mentioned in Step 1, we further infer that the vector-valued
coefficient function $(\mathbf{c}^{k})^{\top}$ also exists on $[0,T]$. This
directly entails that the approximate solution
$({\boldsymbol{v}}_{k},\phi_{k},\Psi_{k},\mu_{k},\Theta_{k})$ actually exists
on $[0,T]$ and satisfies the discretized weak formulation (3.34) everywhere in
$[0,T]$. Moreover, as the choice of $T_{k}$ did not play any role in the
derivation of the uniform estimates, we further conclude that the estimates
(3.39), (3.42), (3.43), (3.47), (3.48), (3.49), (3.50) and (3.51) remain true
when $T_{k}$ is replaced by $T$. In summary, we thus know that the approximate
weak solution $({\boldsymbol{v}}_{k},\phi_{k},\Psi_{k},\mu_{k},\Theta_{k})$
satisfies the uniform estimates

$\displaystyle\begin{aligned}
&\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{\frac{4}{3}}(0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}+\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;H^{1}(\Omega))}+\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;L^{2}(\Gamma))}\\\
&\;\;+\|(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})\|_{L^{2}(0,T;\mathcal{D}_{\alpha}^{\prime})}+\|(\phi_{k},\Psi_{k})\|_{L^{\infty}(0,T;\mathcal{D}_{\alpha})}\\\
&\;\;+\|(\mu_{k},\Theta_{k})\|_{L^{4}(0,T;\mathcal{L}^{2})}+\|(\mu_{k},\Theta_{k})\|_{L^{2}(0,T;\mathcal{H}^{1})}\leq
C,\quad\text{if $d=3$},\end{aligned}$ (3.52a) $\displaystyle\begin{aligned}
&\|\partial_{t}{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})}+\|{\boldsymbol{v}}_{k}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;H^{1}(\Omega))}+\|{\boldsymbol{v}}_{k}\|_{L^{2}(0,T;L^{2}(\Gamma))}\\\
&\;\;+\|(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})\|_{L^{2}(0,T;\mathcal{D}_{\alpha}^{\prime})}+\|(\phi_{k},\Psi_{k})\|_{L^{\infty}(0,T;\mathcal{D}_{\alpha})}\\\
&\;\;+\|(\mu_{k},\Theta_{k})\|_{L^{4}(0,T;\mathcal{L}^{2})}+\|(\mu_{k},\Theta_{k})\|_{L^{2}(0,T;\mathcal{H}^{1})}\leq
C,\quad\text{if $d=2$}.\end{aligned}$ (3.52b)

Step 4: Convergence to a weak solution. In view of the uniform estimates
(3.52a) and (3.52b), the Banach–Alaoglu theorem and the Aubin–Lions–Simon
lemma (see, e.g., [13, Theorem II.5.16]) imply the existence of functions
${\boldsymbol{v}}$, $\phi$, $\Psi$, $\mu$ and $\theta$ such that

$\displaystyle\partial_{t}{\boldsymbol{v}}_{k}$
$\displaystyle\to\partial_{t}{\boldsymbol{v}}$ weakly in
$L^{\frac{4}{3}}(0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})$ if
$d=3$, (3.53) weakly in
$L^{2}(0,T;\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})$ if $d=2$,
(3.54) $\displaystyle{\boldsymbol{v}}_{k}$ $\displaystyle\to{\boldsymbol{v}}$
weakly-star in $L^{\infty}(0,T;\mathbf{L}_{\textnormal{div}}^{2}(\Omega))$,
weakly in $L^{2}(0,T;\mathbf{H}^{1}(\Omega))\cap
L^{2}(0,T;\mathbf{L}^{2}(\Gamma))$, strongly in
$L^{2}(0,T;\mathbf{L}_{\textnormal{div}}^{2}(\Omega))\cap
C([0,T];\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime})$, (3.55)
$\displaystyle(\partial_{t}\phi_{k},\partial_{t}\Psi_{k})$
$\displaystyle\to(\partial_{t}\phi,\partial_{t}\Psi)$ weakly in
$L^{2}(0,T;\mathcal{D}_{\alpha}^{\prime})$, (3.56)
$\displaystyle(\phi_{k},\Psi_{k})$ $\displaystyle\to(\phi,\Psi)$ weakly-star
in $L^{\infty}(0,T;{\mathcal{D}_{\alpha}})$, strongly in
$C([0,T];\mathcal{H}^{s})$ for all $s\in[0,1)$, (3.57)
$\displaystyle(\mu_{k},\Theta_{k})$ $\displaystyle\to(\mu,\Theta)$ weakly in
$L^{4}(0,T;\mathcal{L}^{2})\cap L^{2}(0,T;{\mathcal{D}_{\alpha}})$, (3.58)

as $k\to\infty$, along a non-relabeled subsequence. In view of (3.22), we set
$\psi:=\alpha\Psi$ and $\theta:=\alpha^{-1}\Theta$. Using Sobolev’s embedding
theorem, we infer from (Proof of Theorem 3.3.) that

$\displaystyle\phi_{k}$ $\displaystyle\to\phi$ strongly in
$C([0,T];L^{r}(\Omega))$ for all $r\in[2,6)$, and a.e. in $Q$, (3.59)
$\displaystyle\Psi_{k}$ $\displaystyle\to\Psi$ strongly in
$C([0,T];L^{r}(\Gamma))$ for all $r\in[2,\infty)$, and a.e. on $\Sigma$,
(3.60)

as $k\to\infty$, after another subsequence extraction. Recalling the growth
assumptions on $G$ and $G^{\prime}$ imposed in ($\mathrm{A3}$), we use the
above convergences along with Lebesgue’s general convergence theorem (see [8,
Section 3.25]) to deduce that

$\displaystyle G(\alpha\Psi_{k})$ $\displaystyle\to G(\alpha\Psi)$ strongly in
$L^{1}(\Sigma)$, and a.e. on $\Sigma$, (3.61) $\displaystyle
G^{\prime}(\alpha\Psi_{k})$ $\displaystyle\to G^{\prime}(\alpha\Psi)$ strongly
in $L^{2}(\Sigma)$, and a.e. on $\Sigma$, (3.62)

as $k\to\infty$. Moreover, using the growth assumptions on $F^{\prime}$ from
($\mathrm{A3}$), the uniform estimate (3.43) and the continuous embedding
$H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)$, we infer that
$F^{\prime}(\phi_{k})$ is bounded in $L^{6/5}(Q)$ uniformly in
$k\in\mathbb{N}$. Hence, there exists a function $f\in L^{6/5}(Q)$ such that
$F^{\prime}(\phi_{k})\to f$ weakly in $L^{6/5}(Q)$ as $k\to\infty$. In view of
(3.59), we also have $F^{\prime}(\phi_{k})\to F^{\prime}(\phi)$ a.e. in $Q$ as
$k\to\infty$. By a convergence principle based on Egorov’s theorem (see [17,
Proposition 9.2c]), we infer $f=F^{\prime}(\phi)$. We have thus shown that

$\displaystyle F^{\prime}(\phi_{k})$ $\displaystyle\to F^{\prime}(\phi)$
weakly in $L^{6/5}(Q)$, and a.e. in $Q$, (3.63)

as $k\to\infty$. Furthermore, it follows from (3.59) and (3.60) that

$\displaystyle\nu(\phi_{k})$ $\displaystyle\to\nu(\phi)$
$\displaystyle\quad\text{and}\quad m_{\Omega}(\phi_{k})$ $\displaystyle\to
m_{\Omega}(\phi)$ a.e. in $Q$, (3.64) $\displaystyle\gamma(\alpha\Psi_{k})$
$\displaystyle\to\gamma(\alpha\Psi)$ $\displaystyle\quad\text{and}\quad
n_{\Gamma}(\Psi_{k})$ $\displaystyle\to n_{\Gamma}(\Psi)$ a.e. in $\Sigma$,
(3.65)

as $k\to\infty$. Moreover, we also have from ($\mathrm{A4}$) and
($\mathrm{A5}$) that

$\displaystyle\nu(\phi_{k})$ $\displaystyle\to\nu(\phi)$
$\displaystyle\quad\text{and}\quad m_{\Omega}(\phi_{k})$ $\displaystyle\to
m_{\Omega}(\phi)$ in $L^{r}(Q)$, (3.66) $\displaystyle\gamma(\alpha\Psi_{k})$
$\displaystyle\to\gamma(\alpha\Psi)$ $\displaystyle\quad\text{and}\quad
n_{\Gamma}(\Psi_{k})$ $\displaystyle\to n_{\Gamma}(\Psi)$ in $L^{r}(\Sigma)$,
(3.67)

for all $r\in[2,\infty)$, as $k\to\infty$. Lastly, by ($\mathrm{A4}$) and
($\mathrm{A5}$), the weak-strong convergence principle and the Lebesgue
convergence theorem, we infer that

$\displaystyle\nu(\phi_{k})\;\mathbf{D}{\boldsymbol{v}}_{k}$
$\displaystyle\to\nu(\phi)\;\mathbf{D}{\boldsymbol{v}}$ weakly in $L^{2}(Q)$,
(3.68) $\displaystyle m_{\Omega}(\phi_{k})\;\nabla\phi_{k}$ $\displaystyle\to
m_{\Omega}(\phi)\;\nabla\phi$ weakly in $L^{2}(Q)$, (3.69)
$\displaystyle\gamma(\alpha\Psi_{k})\;{\boldsymbol{v}}_{k}$
$\displaystyle\to\gamma(\alpha\Psi)\;{\boldsymbol{v}}$ weakly in
$L^{2}(\Sigma)$, (3.70) $\displaystyle
n_{\Gamma}(\Psi_{k})\;\nabla_{\Gamma}\Psi_{k}$ $\displaystyle\to
n_{\Gamma}(\Psi)\;\nabla_{\Gamma}\Psi$ weakly in $L^{2}(\Sigma)$, (3.71)

as $k\to\infty$.

We now multiply all equations of the discretized weak formulation (3.34) by an
arbitrary test function $\sigma\in C_{0}^{\infty}([0,T])$ and integrate the
resulting equations with respect to time from $0$ to $T$. Using the
convergences (3.53)–(3.71), we can then pass to the limit $k\to\infty$ to
conclude that

$\displaystyle\begin{aligned}
&\int_{0}^{T}\langle\partial_{t}{\boldsymbol{v}},{\boldsymbol{w}}_{i}\rangle_{\mathbf{H}_{\textnormal{div}}^{1}(\Omega)}\;\sigma\;\mathrm{d}t-\int_{Q}({\boldsymbol{v}}\otimes{\boldsymbol{v}})\nabla{\boldsymbol{w}}_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t+\int_{Q}2\nu(\phi)\mathbf{D}{\boldsymbol{v}}:\mathbf{D}{\boldsymbol{w}}_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t\\\
&\quad=\int_{Q}\mu\nabla\phi\cdot{\boldsymbol{w}}_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t-\int_{\Sigma}\gamma(\alpha\Psi){\boldsymbol{v}}\cdot{\boldsymbol{w}}_{i}\;\sigma+\Psi\nabla_{\Gamma}\Theta\cdot{\boldsymbol{w}}_{i}\;\sigma\;\mathrm{d}S\;\mathrm{d}t,\end{aligned}$
(3.72a) $\displaystyle\begin{aligned}
&\int_{0}^{T}\big{<}(\partial_{t}\phi,\partial_{t}\Psi),(\zeta_{i},\xi_{i})\big{>}_{{\mathcal{D}_{\alpha}}}\;\sigma\;\mathrm{d}t-\int_{Q}\phi{\boldsymbol{v}}\cdot\nabla\zeta_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t-\int_{\Sigma}\Psi{\boldsymbol{v}}\cdot\nabla_{\Gamma}\xi_{i}\;\sigma\;\mathrm{d}S\;\mathrm{d}t\\\
&\quad=-\int_{Q}m_{\Omega}(\phi)\nabla\mu\cdot\nabla\zeta_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t-\int_{\Sigma}n_{\Gamma}(\Psi)\nabla_{\Gamma}\Theta\cdot\nabla_{\Gamma}\xi_{i}\;\sigma\;\mathrm{d}S\;\mathrm{d}t,\end{aligned}$
(3.72b) $\displaystyle\begin{aligned}
&\int_{Q}\mu\,\zeta_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t+\int_{\Sigma}\Theta\,\xi_{i}\;\sigma\;\mathrm{d}S\;\mathrm{d}t\\\
&\quad=\int_{Q}\nabla\phi\cdot\nabla\zeta_{i}\;\sigma+F^{\prime}(\phi)\zeta_{i}\;\sigma\;\mathrm{d}x\;\mathrm{d}t+\int_{\Sigma}\alpha^{2}\nabla_{\Gamma}\Psi\cdot\nabla_{\Gamma}\xi_{i}\;\sigma+\alpha
G^{\prime}(\alpha\Psi)\xi_{i}\;\sigma\;\mathrm{d}S\;\mathrm{d}t\end{aligned}$
(3.72c)

hold for all $i\in\mathbb{N}$ and all test functions $\sigma\in
C_{0}^{\infty}([0,T])$. Since
$\mathrm{span}\\{{\boldsymbol{w}}_{i}\,|\,i\in\mathbb{N}\\}$ is dense in
$\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ and
$\mathrm{span}\\{(\zeta_{i},\xi_{i})\,|\,i\in\mathbb{N}\\}$ is dense in
${\mathcal{D}_{\alpha}}$, this proves that the quintuplet
$({\boldsymbol{v}},\phi,\Psi,\mu,\Theta)$ satisfies the weak formulation
(3.16) for all test functions
${\boldsymbol{w}}\in\mathbf{H}^{1}_{\textnormal{div}}(\Omega)$ and
$(\zeta,\xi),(\eta,\vartheta)\in{\mathcal{D}_{\alpha}}$. Reversing the change
of variables (3.22), we conclude that Definition 3.2(iii) is fulfilled.

Integrating the weak formulation (3.16b) with respect to time from $0$ to an
arbitrary $t\in[0,T]$ and choosing $(\zeta,\xi)\equiv(\alpha,1)$, we obtain
(cf. (3.42))

$\displaystyle\alpha\int_{\Omega}\phi(t)\;\mathrm{d}x+\int_{\Gamma}\Psi(t)\;\mathrm{d}S=\alpha\int_{\Omega}\phi_{0}\;\mathrm{d}x+\int_{\Gamma}\Psi_{0}\;\mathrm{d}S$

By the change of variables (3.22), this proves the mass conservation law
(3.17) meaning that the condition in Definition 3.2(iv) is fulfilled.

Furthermore, it follows from Strauss’ lemma (see [52, Corollary 2.1]) that

${\boldsymbol{v}}\in
C_{w}([0,T];\mathbf{L}_{\textnormal{div}}^{2}(\Omega)),\quad\text{and}\quad(\phi,\Psi)\in
C_{w}([0,T];\mathcal{D}_{\alpha}).$ (3.73)

In the two dimensional case, in view of the enhanced time integrability of
$\partial_{t}{\boldsymbol{v}}$, we even obtain ${\boldsymbol{v}}\in
C([0,T];\mathbf{L}_{\textnormal{div}}^{2}(\Omega))$. In summary, we conclude
that the regularity condition (3.24) is satisfied and thus, by the change of
variables (3.22), Definition 3.2(i) is fulfilled.

In view of (3.35) and (3.36), we have

$\displaystyle{\boldsymbol{v}}_{k}(0)$ $\displaystyle\to{\boldsymbol{v}}_{0}$
$\displaystyle\quad\text{strongly in
$\mathbf{L}^{2}_{\textnormal{div}}(\Omega)$},$ (3.74)
$\displaystyle\big{(}\phi_{k}(0),\Psi_{k}(0)\big{)}$
$\displaystyle\to(\phi_{0},\Psi_{0})$ $\displaystyle\quad\text{strongly in
$\mathcal{L}^{2}$},$ (3.75)

as $k\to\infty$, due to the convergence properties of the projections
$\mathbb{P}_{\mathcal{Z}_{k}}$ and $\mathbb{P}_{\mathcal{U}_{k}}$. On the
other hand, we infer from (Proof of Theorem 3.3.) and (Proof of Theorem 3.3.)
that

$\displaystyle{\boldsymbol{v}}_{k}(0)$ $\displaystyle\to{\boldsymbol{v}}(0)$
$\displaystyle\quad\text{strongly in
$\mathbf{H}^{1}_{\textnormal{div}}(\Omega)^{\prime}$},$ (3.76)