# On the solution of the coupled steady-state dual-porosity-Navier-Stokes
fluid flow model with the Beavers-Joseph-Saffman interface
condition11footnotemark: 1

Di Yang<EMAIL_ADDRESS>Yinnian He<EMAIL_ADDRESS>Luling Cao
<EMAIL_ADDRESS>School of Mathematics and Statistics, Xi’an Jiaotong
University, Xi’an, Shaanxi 710049, P.R. China

###### Abstract

In this work, we propose a new analysis strategy to establish an a priori
estimate of the weak solutions to the coupled steady-state dual-porosity-
Navier-Stokes fluid flow model with the Beavers-Joseph-Saffman interface
condition. The most advantage of our proposed method is that the a priori
estimate and the existence result are independent of small data and the large
viscosity restriction. Therefore the global uniqueness of the weak solution is
naturally obtained.

###### keywords:

weak solution, dual-porosity-Navier-Stokes, Beavers-Joseph-Saffman interface
condition, a priori estimate, existence, global uniqueness

††journal:

##  1 Introduction

Coupled free flow and porous medium flow systems play an important role in
many practical engineering fields, e.g., the flood simulation of arid areas in
geological science [20], filtration treatment in industrial production [28,
37], petroleum exploitation in mining, and blood penetration between vessels
and organs in life science [19]. Specifically, the systems are usually
described by Navier-Stokes equations (or Stokes equations) coupled with
Darcy’s equation, and there are amounts of achievements such as [12, 13, 10,
16, 26, 36, 48]. However, the standard Darcy’s equation describes fluids
flowing through only a single porosity medium, which is not accurate to deal
with the complicated multiple porous media similar to naturally fractured
reservoir. Actually, the naturally fractured reservoir is comprised of low
permeable rock matrix blocks surrounded by an irregular network of natural
microfractures, and further they have different fluid storage and conductivity
properties [43, 47]. In 2016, Hou et al. [31] proposed and numerically solved
a coupled dual-porosity-Stokes fluid flow model with four multi-physics
interface conditions. The authors used the dual-porosity equations over
Darcy’s region to describe fluid flowing through the multiple porous medium.
Recently, several related research on the above model can be found in the
literatures [2, 3, 24, 29, 44]. In particular, Gao and Li [24] proposed a
decoupled stabilized finite element method to solve the coupled dual-porosity-
Navier-Stokes fluid flow model in the numerical field.

The steady-state dual-porosity-Navier-Stokes fluid flow model has distinct
features and difficulties in mathematical analysis. Many numerical methods
have been studied for the well-known stationary or time-dependent Navier-
Stokes/Darcy model with Beavers-Joseph or Beaver-Joseph-Saffman interface
condition, including coupled finite element methods [6, 11, 21, 22, 49],
discontinuous Galerkin methods [17, 26, 27], domain decomposition methods [12,
20, 21, 30, 38], and decoupled methods based on two-grid finite element [10,
23, 48, 49]. In spite of the above great contributions to numerical
simulation, the existence of a weak solution to the coupled dual-porosity-
Navier-Stokes fluid flow model with Beavers-Joseph-Saffman interface condition
for general data keeps unresolved. In many literatures [6, 10, 17, 20, 21, 22,
26, 27, 30, 49], a priori estimates and existence of a weak solution need
suitable small data and/or large viscosity restrictions, and therefore only
local uniqueness can be established when the data satisfy additional
restrictions. In [26], the authors pointed out that the difficulty for a
priori estimates and existence with general data is stemmed from the
transmission interface condition, which does not completely compensate the
nonlinear convection term from the Navier-Stokes equations in the energy
balance.

Therefore in this paper, stemming from resolving steady-state Navier-Stokes
equations with mixed boundary conditions in [32], we shall establish a new a
priori estimate of the weak solutions by coupling the model problem with a
designed auxiliary problem in order to completely compensate the nonlinear
convection term from the Navier-Stokes equations. In addition, we shall also
prove existence of a weak solution without small data or large viscosity
restriction. As a result, the global uniqueness of the weak solution is
naturally obtained.

The rest of this paper is organized as follows. In Section 2, we specify the
steady-state dual-porosity-Navier-Stokes fluid flow model with Beavers-Joseph-
Saffman interface condition and provide its Galerkin variational formulation.
In Section 3, we establish a new a priori estimate of the weak solutions by
coupling the model problem with an auxiliary problem, which is designed
subject to the model problem. Finally, in Section 4, we prove existence of the
weak solution without small data and/or large viscosity restriction by the
Galerkin method and Brouwer’s fixed-point theorem, and global uniqueness of
all variables by the inf-sup condition and Babus̆ka–Brezzi’s theory.

##  2 Model specification

###  2.1 Setting of the problem

In this section, we consider the steady-state dual-porosity-Navier-Stokes
fluid flow model in a bounded open polygonal domain
$\Omega\subset\mathbb{R}^{N}(N=2,3)$ with four physically valid interface
conditions.

Figure 1: Schematic representation of fluid flow region $\Omega_{s}$ and dual-
porous medium region $\Omega_{d}$, separated by an interface $\Gamma$.

The domain $\Omega$ consists of a fluid flow region $\Omega_{s}$ and a dual-
porous medium region $\Omega_{d}$, with interface
$\Gamma=\overline{\Omega}_{s}\cap\overline{\Omega}_{d}$ (see Figure 1). Both
$\Omega_{s}$ and $\Omega_{d}$ are open, regular, simply connected, and bounded
with Lipschitz continuous boundaries
$\Gamma_{i}=\partial\Omega_{i}\setminus\Gamma~{}(i=s,d)$, respectively. Here,
$\Omega_{s}\cap\Omega_{d}=\emptyset$,
$\overline{\Omega}_{s}\cup\overline{\Omega}_{d}=\Omega$. The unit normal
vector of the interface $\Gamma$ pointing from $\Omega_{s}$ to $\Omega_{d}$
(from $\Omega_{d}$ to $\Omega_{s}$) is denoted by $\bm{n}_{s}$ (resp.
$\bm{n}_{d}$), and the corresponding unit tangential vectors are denoted by
$\bm{\tau}_{i}$, $i=1,\cdots,N-1$. In two-dimensional case, if we write
$\bm{n}_{s}=(n_{1},n_{2})\in\mathbb{R}^{2}$, then
$\bm{\tau}=\bm{n}_{s}^{\top}=(-n_{2},n_{1})$.

In $\Omega_{s}$, the fluid flow is governed by the Navier-Stokes equation [11,
12, 17, 20, 21, 26, 27]:

$\begin{cases}\left(\bm{u}_{s}\cdot\nabla\right)\bm{u}_{s}-\nabla\cdot\mathbb{T}(\bm{u}_{s},p_{s})=\bm{f}_{s}&\text{in}~{}\Omega_{s},\\\
\nabla\cdot\bm{u}_{s}=0&\text{in}~{}\Omega_{s},\\\
\bm{u}_{s}=\bm{0}&\text{on}~{}\Gamma_{s},\end{cases}$ (2.1)

where
$\mathbb{T}(\bm{u}_{s},p_{s})=-p_{s}\mathbb{I}+2\nu\mathbb{D}(\bm{u}_{s})$ is
the stress tensor, $\nu$ is the kinetic viscosity,
$\mathbb{D}(\bm{u}_{s})=\frac{1}{2}\left[\nabla\bm{u}_{s}+(\nabla\bm{u}_{s})^{\top}\right]$
is the velocity deformation tensor, $\bm{u}_{s}(\bm{x})$ denotes the fluid
velocity in $\Omega_{s}$, $p_{s}(\bm{x})$ denotes the kinematic pressure in
$\Omega_{s}$, and $\bm{f}_{s}(\bm{x})$ denotes a general body force term that
includes gravitational acceleration. As usual, we write formally:

$(\bm{v}\cdot\nabla)\bm{w}=\sum_{i=1}^{N}v_{i}\frac{\partial\bm{w}}{\partial
x_{i}},\quad\nabla\cdot\bm{v}=\sum_{i=1}^{N}\frac{\partial v_{i}}{\partial
x_{i}}.$

The filtration of an incompressible fluid through porous media is often
described using Darcy’s law. So in dual-porous medium domain $\Omega_{d}$, the
flow is governed by a traditional dual-porosity model, which is composed of
matrix and microfracture equations as follows:

$\begin{cases}-\nabla\cdot\big{(}\frac{k_{m}}{\mu}\nabla\phi_{m}\big{)}+\frac{\sigma
k_{m}}{\mu}(\phi_{m}-\phi_{f})=0&\text{in}~{}\Omega_{d},\\\
-\nabla\cdot\big{(}\frac{k_{f}}{\mu}\nabla\phi_{f}\big{)}+\frac{\sigma
k_{m}}{\mu}(\phi_{f}-\phi_{m})=f_{d}&\text{in}~{}\Omega_{d},\\\
\phi_{m}=0&\text{on}~{}\Gamma_{d},\\\
\phi_{f}=0&\text{on}~{}\Gamma_{d},\end{cases}$ (2.2)

where $\phi_{m}(\bm{x})$, $\phi_{f}(\bm{x})$ denote the matrix and
microfracture flow pressure, $\mu$ is the dynamic viscosity, $k_{m}$, $k_{f}$
are the intrinsic permeability of the matrix and microfracture regions,
$\sigma$ is the shape factor characterizing the morphology and dimension of
the microfractures, $f_{d}$ is the sink/source term for the microfractures,
and the term $\frac{\sigma k_{m}}{\mu}(\phi_{m}-\phi_{f})$ describes the mass
exchange between matrix and microfractures.

Based on the fundamental properties of the dual-porosity fluid flow model and
traditional Stokes–Darcy flow model, Hou et al. [31] introduced four
physically valid interface conditions as follows to couple appropriately the
dual-porosity-Stokes model, which are also adapted to our problem.

  * 1.

No-exchange condition between the matrix and the conduits/macrofractures:

$-\frac{k_{m}}{\mu}\nabla\phi_{m}\cdot\bm{n}_{d}=0\quad\text{on}~{}\Gamma.$
(2.3)

  * 2.

Mass conservation:

$\bm{u}_{s}\cdot\bm{n}_{s}=\frac{k_{f}}{\mu}\nabla\phi_{f}\cdot\bm{n}_{d}\quad\text{on}~{}\Gamma.$
(2.4)

  * 3.

Balance of normal forces:

$-\bm{n}_{s}\cdot(\mathbb{T}(\bm{u}_{s},p_{s})\bm{n}_{s})=\frac{\phi_{f}}{\rho}\quad\text{on}~{}\Gamma.$
(2.5)

  * 4.

The Beavers-Joseph-Saffman interface condition [7, 33, 35, 41]: for
$i=1,\cdots,N-1$,

$-\bm{\tau}_{i}\cdot\left(\mathbb{T}(\bm{u}_{s},p_{s})\bm{n}_{s}\right)=\frac{\alpha\nu\sqrt{N}}{\sqrt{\text{trace}(\bm{\Pi})}}\bm{u}_{s}\cdot\bm{\tau}_{i}\quad\text{on}~{}\Gamma.$
(2.6)

In (2.3)–(2.6), $\alpha$ is the Beavers-Joseph constant depending on the
properties of the dual-porous medium, $\bm{\Pi}$ represents the intrinsic
permeability that satisfies the relation $\bm{\Pi}=k_{f}\mathbb{I}$,
$\mathbb{I}$ is the $N\times N$ identity matrix, and $\rho$ is the fluid
density.

###  2.2 Galerkin variational formulation

Throughout this paper we use the following standard function spaces. For a
Lipschitz domain $\mathcal{D}\subset\mathbb{R}^{N}$, $N\geqslant 1$, we denote
by $W^{k,p}(\mathcal{D})$ the Sobolev space with indexes $k\geqslant 0$,
$1\leqslant p\leqslant\infty$ of real-valued functions defined on
$\mathcal{D}$, endowed with the seminorm $|\cdot|_{W^{k,p}(\mathcal{D})}$
denoted by $|\cdot|_{k,p,\mathcal{D}}$ and norm
$\|\cdot\|_{W^{k,p}(\mathcal{D})}$ denoted by $\|\cdot\|_{k,p,\mathcal{D}}$
[1]. When $p=2$, $W^{s,2}(\mathcal{D})$ is denoted as $H^{k}(\mathcal{D})$ and
the corresponding seminorm and norm are written as $|\cdot|_{k,\mathcal{D}}$
and $\|\cdot\|_{k,\mathcal{D}}$, respectively. In addition, with
$|\mathcal{D}|$ we denote the $N$-dimensional Hausdorff measure of
$\mathcal{D}$.

To perform the variational formulation, we define some necessary Hilbert
spaces given by

$\displaystyle\bm{X}_{s}$
$\displaystyle:=\left\\{\bm{v}\in\bm{H}^{1}(\Omega_{s}):\bm{v}=\bm{0}~{}\text{on}~{}\Gamma_{s}\right\\},$
$\displaystyle X_{d}$ $\displaystyle:=\left\\{\psi\in
H^{1}(\Omega_{d}):\psi=0~{}\text{on}~{}\Gamma_{d}\right\\},$ $\displaystyle
Q_{s}$ $\displaystyle:=L^{2}(\Omega_{s}).$

We also need the trace space $\bm{H}_{00}^{1/2}(\Gamma):=\bm{X}_{s}|_{\Gamma}$
(resp. $H_{00}^{1/2}(\Gamma):=X_{d}|_{\Gamma}$), which is a nonclosed subspace
of $\bm{H}^{1/2}(\Gamma)$ (resp. $H^{1/2}(\Gamma)$) and has a continuous zero
extension to $\bm{H}^{1/2}(\partial\Omega_{s})$ (resp.
$H^{1/2}(\partial\Omega_{d})$) [14, 15, 20]. For the trace space
$\bm{H}_{00}^{1/2}(\Gamma)$ and its dual space
$(\bm{H}_{00}^{1/2}(\Gamma))^{\prime}$, we have the following continuous
imbedding result [15]:

$\bm{H}_{00}^{1/2}(\Gamma)\subsetneqq\bm{H}^{1/2}(\Gamma)\subsetneqq\bm{L}^{2}(\Gamma)\subsetneqq\bm{H}^{-1/2}(\Gamma)\subsetneqq(\bm{H}_{00}^{1/2}(\Gamma))^{\prime}.$
(2.7)

One can see more details in [15, 34, 45] and the references therein. For any
bounded domain $\mathcal{D}\subset\mathbb{R}^{N}$,
$(\cdot,\cdot)_{\mathcal{D}}$ denotes the $L^{2}$ inner product on
$\mathcal{D}$, and $\langle\cdot,\cdot\rangle_{\partial\mathcal{D}}$ denotes
the $L^{2}$ inner product (or duality pairing) on the boundary
$\partial\mathcal{D}$. We also consider the following product Hilbert space
$\underline{\bm{Y}}:=\bm{X}_{s}\times X_{d}\times X_{d}$ with norm

$\|\underline{\bm{w}}\|_{\underline{\bm{Y}}}=\Big{(}\|\bm{w}_{s}\|_{1,\Omega_{s}}^{2}+\|\psi_{m}\|_{1,\Omega_{d}}^{2}+\|\psi_{f}\|_{1,\Omega_{d}}^{2}\Big{)}^{1/2},\quad\forall\,\underline{\bm{w}}=(\bm{w}_{s},\psi_{m},\psi_{f})\in\underline{\bm{Y}}.$

In addition, based on the following formula:

$\left((\bm{u}\cdot\nabla)\bm{v},\bm{w}\right)_{\mathcal{D}}=\left\langle\bm{u}\cdot\bm{n},\bm{v}\cdot\bm{w}\right\rangle_{\partial\mathcal{D}}-\left((\bm{u}\cdot\nabla)\bm{w},\bm{v}\right)_{\mathcal{D}}-\left((\nabla\cdot\bm{u})\bm{v},\bm{w}\right)_{\mathcal{D}},\quad\forall\,\bm{u},\bm{v},\bm{w}\in\bm{H}^{1}(\mathcal{D}),$

we introduce the trilinear form $b(\cdot;\cdot,\cdot)$ given by
$\forall\,\bm{u}_{s},\bm{v}_{s},\bm{w}_{s}\in\bm{X}_{s}$,

$\begin{split}b(\bm{u}_{s};\bm{v}_{s},\bm{w}_{s})&=\left((\bm{u}_{s}\cdot\nabla)\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}+\frac{1}{2}\left((\nabla\cdot\bm{u}_{s})\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}\\\
&=\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},\bm{v}_{s}\cdot\bm{w}_{s}\right\rangle_{\Gamma}+\frac{1}{2}\left((\bm{u}_{s}\cdot\nabla)\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}-\frac{1}{2}\left((\bm{u}_{s}\cdot\nabla)\bm{w}_{s},\bm{v}_{s}\right)_{\Omega_{s}}.\end{split}$
(2.8)

Hence, the Galerkin variational formulation of the coupled problem (2.1)–(2.6)
is proposed that: to find
$(\underline{\bm{u}},p_{s})\in\underline{\bm{Y}}\times Q_{s}$ such that

$a(\underline{\bm{u}},\underline{\bm{v}})+d(\bm{v}_{s},p_{s})-d(\bm{u}_{s},q_{s})+b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s})=\rho(\bm{f}_{s},\bm{v}_{s})_{\Omega_{s}}+(f_{d},\psi_{f})_{\Omega_{d}},\quad\forall\,(\underline{\bm{v}},q_{s})\in\underline{\bm{Y}}\times
Q_{s},$ (2.9)

where $\underline{\bm{u}}=(\bm{u}_{s},\phi_{m},\phi_{f})$,
$\underline{\bm{v}}=(\bm{v}_{s},\psi_{m},\psi_{f})$, the bilinear forms
$a(\cdot,\cdot)$ and $d(\cdot,\cdot)$ are defined as

$\displaystyle a_{s}(\underline{\bm{u}},\underline{\bm{v}})$
$\displaystyle=2\rho\nu\left(\mathbb{D}(\bm{u}_{s}),\mathbb{D}(\bm{v}_{s})\right)_{\Omega_{s}},$
$\displaystyle a_{d}(\underline{\bm{u}},\underline{\bm{v}})$
$\displaystyle=\frac{k_{m}}{\mu}\left(\nabla\phi_{m},\nabla\psi_{m}\right)_{\Omega_{d}}+\frac{k_{f}}{\mu}\left(\nabla\phi_{f},\nabla\psi_{f}\right)_{\Omega_{d}}+\frac{\sigma
k_{m}}{\mu}\left(\phi_{m}-\phi_{f},\psi_{m}\right)_{\Omega_{d}}+\frac{\sigma
k_{m}}{\mu}\left(\phi_{f}-\phi_{m},\psi_{f}\right)_{\Omega_{d}},$
$\displaystyle a_{\Gamma}(\underline{\bm{u}},\underline{\bm{v}})$
$\displaystyle=\left\langle\phi_{f},\bm{v}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}-\left\langle\psi_{f},\bm{u}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}+\sum_{i=1}^{N-1}\frac{\alpha\rho\nu}{\sqrt{k_{f}}}\left\langle\bm{u}_{s}\cdot\bm{\tau}_{i},\bm{v}_{s}\cdot\bm{\tau}_{i}\right\rangle_{\Gamma},$
$\displaystyle a(\underline{\bm{u}},\underline{\bm{v}})$
$\displaystyle=a_{s}(\underline{\bm{u}},\underline{\bm{v}})+a_{d}(\underline{\bm{u}},\underline{\bm{v}})+a_{\Gamma}(\underline{\bm{u}},\underline{\bm{v}}),$
$\displaystyle d(\bm{v}_{s},q_{s})$
$\displaystyle=-\rho(\nabla\cdot\bm{v}_{s},q_{s})_{\Omega_{s}}.$

######  Remark 2.1.

We note that the term
$\frac{1}{2}\left((\nabla\cdot\bm{u}_{s})\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}$
vanishes in (2.8) if $\nabla\cdot\bm{u}_{s}=0$.

######  Remark 2.2.

For the sake of clarity, in our analysis we shall adopt homogenous boundary
conditions. In addition, for the general Dirichlet boundary conditions:

$\bm{u}|_{\Gamma_{s}}=\bm{u}^{\text{dir}},\quad\phi_{m}|_{\Gamma_{d}}=\phi_{m}^{\text{dir}},\quad\phi_{f}|_{\Gamma_{d}}=\phi_{f}^{\text{dir}},$

the standard homogenization technique and the lifting operators employed in
[20] can be used to obtain an equivalent system with the homogeneous Dirichlet
boundary conditions.

Thanks to [17, 26], we have the following important result:

######  Lemma 2.1.

Assume that $\bm{f}_{s}\in\bm{L}^{2}(\Omega_{s})$ and $f_{d}\in
L^{2}(\Omega_{d})$. Then if
$(\bm{u}_{s},\phi_{f},\phi_{m},p_{s})\in\bm{X}_{s}\times X_{d}\times
X_{d}\times Q_{s}$ satisfies (2.1)–(2.6), it is also a solution to problem
(2.9). Conversely, any solution of problem (2.9) satisfies (2.1)–(2.6).

##  3 An a priori estimate

In this section, we shall propose an a priori estimate for possible solutions
of (2.9).

###  3.1 Some technique inequalities

Firstly, throughout this paper we use $C$ to denote a generic positive
constant independent of discretization parameters, which may take different
values in different occasions. Then, based on general Sobolev inequalities,
the trace theorem and the Sobolev embedding theorem [1], we have that for any
bounded open set $\mathcal{D}\subset\mathbb{R}^{N}$ with Lipschitz continuous
boundary $\partial\mathcal{D}$ and for all $v\in H^{1}(\mathcal{D})$,

$\displaystyle\|v\|_{0,q,\mathcal{D}}$ $\displaystyle\leqslant
C\|v\|_{1,\mathcal{D}},\quad 1\leqslant q\leqslant 6,$ (3.1)
$\displaystyle\|v\|_{0,\partial\mathcal{D}}$ $\displaystyle\leqslant
C\|v\|_{0,\mathcal{D}}^{1/2}\|v\|_{1,\mathcal{D}}^{1/2}\leqslant
C\|v\|_{1,\mathcal{D}},$ (3.2) $\displaystyle\|v\|_{0,4,\partial\mathcal{D}}$
$\displaystyle\leqslant C\|v\|_{1,\mathcal{D}},$ (3.3)

which can be also found in [17, 26]. We also need Poincaré inequality and
Korn’s inequality [15] that for all $\bm{v}_{s}\in\bm{X}_{s}$ and for all
$\psi\in X_{d}$,

$\displaystyle\|\bm{v}_{s}\|_{0,\Omega_{s}}$ $\displaystyle\leqslant
C|\bm{v}_{s}|_{1,\Omega_{s}},$ (3.4) $\displaystyle\|\psi\|_{0,\Omega_{d}}$
$\displaystyle\leqslant C|\psi|_{1,\Omega_{d}},$ (3.5)
$\displaystyle|\bm{v}_{s}|_{1,\Omega_{s}}$ $\displaystyle\leqslant
C\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}.$ (3.6)

Moreover, for the trilinear form $b(\cdot;\cdot,\cdot)$, the following lemma
holds with the help of [4, 8], (3.1), (3.4) and (3.6).

######  Lemma 3.1.

For any functions $\bm{u}_{s}$, $\bm{v}_{s}$, $\bm{w}_{s}\in\bm{X}_{s}$, we
have

$\displaystyle\left|\left((\bm{u}_{s}\cdot\nabla)\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}\right|$
$\displaystyle\leqslant
C\|\bm{u}_{s}\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}\|\mathbb{D}(\bm{w}_{s})\|_{0,\Omega_{s}},$
(3.7)
$\displaystyle\left|\left((\nabla\cdot\bm{u}_{s})\bm{v}_{s},\bm{w}_{s}\right)_{\Omega_{s}}\right|$
$\displaystyle\leqslant
C\|\bm{v}_{s}\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}\|\mathbb{D}(\bm{w}_{s})\|_{0,\Omega_{s}}.$
(3.8)

######  Proof 1.

First, let us recall the standard interpolation inequality [4, 8]: for any
bounded open set $\mathcal{D}\subset\mathbb{R}^{N}$ with Lipschitz continuous
boundary $\partial\mathcal{D}$, $1\leqslant p<q<\infty$ and
$\theta=N(1/p-1/q)\leqslant 1$,

$\|v\|_{0,q,\mathcal{D}}\leqslant
C\|v\|_{0,p,\mathcal{D}}^{1-\theta}\|v\|_{1,p,\mathcal{D}}^{\theta},\quad\forall\,v\in
W^{1,p}(\mathcal{D}).$ (3.9)

In (3.9), we set $\mathcal{D}=\Omega_{s}$, and

  * 1.

if $N=2$, $p=2$, $q=4$, we then have

$\|\bm{v}_{s}\|_{0,4,\Omega_{s}}\leqslant
C\|\bm{v}_{s}\|_{0,\Omega_{s}}^{1/2}\|\bm{v}_{s}\|_{1,\Omega_{s}}^{1/2},\quad\forall\,\bm{v}_{s}\in\bm{X}_{s};$
(3.10)

  * 2.

if $N=3$, $p=2$, $q=3$, we then have

$\|\bm{v}_{s}\|_{0,3,\Omega_{s}}\leqslant
C\|\bm{v}_{s}\|_{0,\Omega_{s}}^{1/2}\|\bm{v}_{s}\|_{1,\Omega_{s}}^{1/2},\quad\forall\,\bm{v}_{s}\in\bm{X}_{s}.$
(3.11)

Hence, the results (3.7) and (3.8) then follow from Hölder inequality, (3.1),
(3.4), (3.6), (3.10) and (3.11) when $N=2$ or $3$.

###  3.2 The equivalent problem

Let us introduce the following divergence-free space:

$\bm{V}_{s}:=\left\\{\bm{v}_{s}\in\bm{X}_{s}:\,\nabla\cdot\bm{v}_{s}=0\right\\},$

and we consider the new product Hilbert space:

$\underline{\bm{W}}:=\bm{V}_{s}\times X_{d}\times X_{d},$

which equipped with the same norm as $\underline{\bm{Y}}$.

Thanks to [26], it provides a lower constant bound $\beta_{0}>0$ depending
only on $\Omega$ to guarantee the following inf-sup condition:

$\inf_{q_{s}\in
Q_{s}}\sup_{\underline{\bm{v}}=(\bm{v}_{s},\psi_{m},\psi_{f})\in\underline{\bm{Y}}}\frac{d(\bm{v}_{s},q_{s})}{\|q_{s}\|_{0,\Omega_{s}}\|\underline{\bm{v}}\|_{\underline{\bm{Y}}}}\geqslant\beta_{0}.$
(3.12)

Hence, by (3.12) and the same argument in [25], the Galerkin variational
formulation (2.9) is equivalent to the following problem: to find
$\underline{\bm{u}}\in\underline{\bm{W}}$ such that

$a(\underline{\bm{u}},\underline{\bm{v}})+b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s})=\rho(\bm{f}_{s},\bm{v}_{s})_{\Omega_{s}}+(f_{d},\psi_{f})_{\Omega_{d}},\quad\forall\,\underline{\bm{v}}\in\underline{\bm{W}}.$
(3.13)

###  3.3 Discretization and an auxiliary problem

Let $\mathcal{R}_{i,h}$ be a quasi-uniform regular triangulation of the domain
$\Omega_{i}$, $i=s,d$, respectively. If assuming that the two meshes
$\mathcal{R}_{s,h}$ and $\mathcal{R}_{d,h}$ coincide along the interface
$\Gamma$, then we can define
$\mathcal{R}_{h}:=\mathcal{R}_{s,h}\cup\Gamma\cup\mathcal{R}_{d,h}$, which is
also a quasi-uniform regular triangulation of $\Omega$. The diameter of
element $K\in\mathcal{R}_{h}$ is denoted as $h_{K}$, and we set the mesh
parameter $h=\max_{K\in\mathcal{R}_{h}}h_{K}$.

Then we denote by $\bm{X}_{h}\subset\bm{H}_{0}^{1}(\Omega)$ a finite element
space defined on $\Omega$. The discrete space $\bm{X}_{h}$ can be naturally
restricted to $\Omega_{s}$, so we set
$\bm{X}_{s,h}=\bm{X}_{h}|_{\Omega_{s}}\subset\bm{X}_{s}$. Following the same
technique, we establish other finite element spaces $Q_{s,h}\subset Q_{s}$ and
$X_{d,h}\subset X_{d}$. We assume that $(\bm{X}_{s,h},Q_{s,h})$ is a stable
finite element space pair. Further we define the following vector-valued
Hilbert space on $\Omega_{d}$:

$\bm{X}_{d}:=\left\\{\bm{v}\in\bm{H}^{1}(\Omega_{d}):\,\bm{v}=\bm{0}~{}\mathrm{on}~{}\Gamma_{d}\right\\}.$

In addition, we need a finite element subspace $\bm{V}_{h}\subset\bm{X}_{h}$
defined on $\Omega$ given by

$\bm{V}_{h}:=\left\\{\bm{v}_{h}\in\bm{X}_{h}:\,d(\bm{v}_{h},q_{s,h})=0,\,\forall\,q_{s,h}\in
Q_{s,h}\right\\}.$

Similarly, we have

$\bm{V}_{s,h}=\bm{V}_{h}|_{\Omega_{s}}\subset\bm{X}_{s},\quad\bm{X}_{d,h}=\bm{V}_{h}|_{\Omega_{d}}\subset\bm{X}_{d},$

where $\bm{V}_{s,h}$ is a weakly divergence-free finite element space on
$\Omega_{s}$, and any function $\bm{v}_{d,h}\in\bm{X}_{d,h}$ has an implicit
restriction that $\int_{\Gamma}\bm{v}_{d,h}\cdot\bm{n}_{d}\mathrm{d}s=0$
because of the continuity of $\bm{V}_{h}$ across the interface $\Gamma$.

According to [26], the difficulty of obtaining an a priori estimate of the
Navier-Stokes/Darcy fluid flow model comes from the energy unbalance due to
the nonlinear convection term $(\bm{u}_{s}\cdot\nabla)\bm{u}_{s}$ in (2.1). It
motivates us to construct an auxiliary discrete Galerkin variational problem
defined on $\Omega_{d}$ with compatible boundary conditions, so that the
auxiliary problem can compensate the nonlinear convection term in the energy
balance of the Navier-Stokes equations.

To fix ideas, we first define a lifting operator
$\mathcal{L}:\,\bm{H}_{00}^{1/2}(\Gamma)\rightarrow\bm{X}_{d}$ as follows: for
any $\bm{\eta}\in\bm{H}_{00}^{1/2}(\Gamma)$ with
$\int_{\Gamma}\bm{\eta}\mathrm{d}s=0$ such that

$\mathcal{L}\bm{\eta}\in\bm{X}_{d},\quad(\mathcal{L}\bm{\eta})|_{\Gamma}=\bm{\eta},\quad\nabla\cdot(\mathcal{L}\bm{\eta})|_{\Omega_{d}}=0.$

Then, we introduce the Scott-Zhang interpolator
$\Pi_{s,h}:\,\bm{X}_{s}\rightarrow\bm{X}_{s,h}$ satisfying the following
properties [42]:

$\displaystyle\|\bm{v}_{s}-\Pi_{s,h}\bm{v}_{s}\|_{0,\Omega_{s}}$
$\displaystyle\leqslant
Ch\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}},\quad\forall\,\bm{v}_{s}\in\bm{X}_{s},$
(3.14) $\displaystyle\|\Pi_{s,h}\bm{v}_{s}\|_{1,\Omega_{s}}$
$\displaystyle\leqslant C\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}},\quad\
\;\forall\,\bm{v}_{s}\in\bm{X}_{s}.$ (3.15)

Now, with the mesh parameter $h$, the interpolator $\Pi_{s,h}$ and any given
$\bm{u}_{s}\in\bm{V}_{s}$, we consider the following auxiliary discrete
Galerkin variational problem: to find $\bm{u}_{d,h}\in\bm{X}_{d,h}$ with
$\bm{u}_{d,h}|_{\Gamma}=(\Pi_{s,h}\bm{u}_{s})|_{\Gamma}$ such that for all
$\bm{v}_{d,h}\in\bm{X}_{d,h}$,

$2\kappa\left(\mathbb{D}(\bm{u}_{d,h}),\mathbb{D}(\bm{v}_{d,h})\right)_{\Omega_{d}}+\left((\bm{u}_{d}^{0}\cdot\nabla)\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}-\kappa\left\langle\frac{\partial\bm{u}_{d,h}}{\partial\bm{n}_{d}},\bm{v}_{d,h}\right\rangle_{\Gamma}=0,$
(3.16)

where $\bm{u}_{d}^{0}=\mathcal{L}(\bm{u}_{s}|_{\Gamma})\in\bm{X}_{d}$, and
$\kappa>0$ is a certain positive constant specified later. Furthermore, the
discrete problem (3.16) is uniquely solvable in $\bm{X}_{d,h}$ from Lemma 3.1
in [32].

######  Remark 3.1.

The discrete problem (3.16) is the conforming Galerkin approximation of the
following convection-diffusion equation defined on $\Omega_{d}$: for any given
$\bm{u}_{s}\in\bm{V}_{s}$, to find $\bm{u}_{d}\in\bm{X}_{d}$ with
$\bm{u}_{d}|_{\Gamma}=\bm{u}_{s}|_{\Gamma}$ such that

$\begin{cases}-2\kappa\nabla\cdot\mathbb{D}(\bm{u}_{d})+(\bm{u}_{d}^{0}\cdot\nabla)\bm{u}_{d}=\bm{0}&\mathrm{in}~{}\Omega_{d},\\\
\bm{u}_{d}=\bm{0}&\mathrm{on}~{}\Gamma_{d}.\end{cases}$ (3.17)

Clearly, for any constant $\kappa>0$ and $\bm{u}_{s}\in\bm{V}_{s}$, the
solution of (3.17) is well-posed [40].

###  3.4 An a priori estimate of weak solutions

We realize that (3.13) and (3.16) can form a larger coupled system because of
the convection term $\bm{u}_{d}^{0}$ and the constraint
$\bm{u}_{d,h}|_{\Gamma}=(\Pi_{s,h}\bm{u}_{s})|_{\Gamma}$ for the unknown
variable $\bm{u}_{d,h}$ over the interface $\Gamma$, and we also stress that
the new coupled system has no energy exchange via the interface $\Gamma$. It
implies that (3.16) is subjected to (3.13) but any possible solution of (3.13)
has nothing to do with the auxiliary problem (3.16).

######  Theorem 3.2.

Assume that the data in the auxiliary discrete problem (3.16) satisfy $h$
small enough and $0<\kappa\leqslant C\rho\nu h$. If problem (3.13) exists a
possible solution $(\bm{u}_{s},\phi_{m},\phi_{f})\in\underline{\bm{W}}$, we
then have the following a priori estimate that

$\rho\nu\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\phi_{m}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2},$
(3.18)

where

$\mathcal{C}^{2}=\rho\nu^{-1}\|\bm{f}_{s}\|_{-1,\Omega_{s}}^{2}+\mu\sigma^{-1}k_{f}^{-1}\|f_{d}\|_{-1,\Omega_{d}}^{2}.$

Here we stress that there are no assumptions for the data and physical
parameters of the model problem.

######  Proof 2.

We denote by
$\underline{\bm{u}}=(\bm{u}_{s},\phi_{m},\phi_{f})\in\underline{\bm{W}}$ a
possible solution of problem (3.13), and denote by
$\bm{u}_{d,h}\in\bm{X}_{d,h}$ the solution of problem (3.16). Then, we assume
that there is a positive finite constant $M_{s}<\infty$ such that
$\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}\leqslant M_{s}$. Taking
$\underline{\bm{v}}=\underline{\bm{u}}$ in (3.13), and noting that the terms
$a_{\Gamma}(\underline{\bm{u}},\underline{\bm{u}})$ and $\frac{\sigma
k_{m}}{\mu}\left(\phi_{m}-\phi_{f},\phi_{m}\right)_{\Omega_{d}}+\frac{\sigma
k_{m}}{\mu}\left(\phi_{f}-\phi_{m},\phi_{f}\right)_{\Omega_{d}}$ are non-
negative, we have

$2\rho\nu\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}+\frac{\sigma
k_{m}}{\mu}\|\nabla\phi_{m}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f}\|_{0,\Omega_{d}}^{2}+\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},|\bm{u}_{s}|^{2}\right\rangle_{\Gamma}\leqslant\rho(\bm{f}_{s},\bm{u}_{s})_{\Omega_{s}}+(f_{d},\phi_{f})_{\Omega_{d}}.$
(3.19)

In addition, we note that $\nabla\cdot\bm{u}_{d}^{0}=0$ in $\Omega_{d}$, and
the $\bm{u}_{d}^{0}|_{\Gamma}=\bm{u}_{s}|_{\Gamma}$ by the definition of
$\mathcal{L}$. Hence, the identity (2.8) can also be applied to the second
term in the left hand of (3.16), that is for all
$\bm{v}_{d,h}\in\bm{X}_{d,h}$,

$\begin{split}&\left((\bm{u}_{d}^{0}\cdot\nabla)\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}=\left((\bm{u}_{d}^{0}\cdot\nabla)\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}+\frac{1}{2}\left((\nabla\cdot\bm{u}_{d}^{0})\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}\\\
&=-\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},|\Pi_{s,h}\bm{u}_{s}|^{2}\right\rangle_{\Gamma}+\frac{1}{2}\left((\bm{u}_{d}^{0}\cdot\nabla)\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}-\frac{1}{2}\left((\bm{u}_{d}^{0}\cdot\nabla)\bm{v}_{d,h},\bm{u}_{d,h}\right)_{\Omega_{d}}.\end{split}$
(3.20)

So, taking $\bm{v}_{d,h}=\bm{u}_{d,h}$ in (3.16) and using (3.20), we obtain

$2\kappa\|\mathbb{D}(\bm{u}_{d,h})\|_{0,\Omega_{d}}^{2}-\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},|\Pi_{s,h}\bm{u}_{s}|^{2}\right\rangle_{\Gamma}=\kappa\left\langle\frac{\partial\bm{u}_{d,h}}{\partial\bm{n}_{d}},\Pi_{s,h}\bm{u}_{s}\right\rangle_{\Gamma}.$
(3.21)

It follows from (3.2), (3.3), (3.4), (3.6), (3.14), (3.15), Hölder inequality
and the triangle inequality that

$\begin{split}&\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},|\bm{u}_{s}|^{2}\right\rangle_{\Gamma}-\frac{1}{2}\left\langle\bm{u}_{s}\cdot\bm{n}_{s},|\Pi_{s,h}\bm{u}_{s}|^{2}\right\rangle_{\Gamma}\\\
\leqslant&\frac{1}{2}\left\|\bm{u}_{s}\right\|_{0,4,\Gamma}\|\bm{u}_{s}+\Pi_{s,h}\bm{u}_{s}\|_{0,4,\Gamma}\left\|\bm{u}_{s}-\Pi_{s,h}\bm{u}_{s}\right\|_{0,\Gamma}\\\
\leqslant&Ch^{1/2}M_{s}\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}.\end{split}$
(3.22)

As follows, we now define the dual norms of $\bm{X}_{s}$ and $X_{d}$, which
are denoted as $\|\cdot\|_{-1,\Omega_{s}}$ and $\|\cdot\|_{-1,\Omega_{d}}$,
respectively.

$\|\bm{f}_{s}\|_{-1,\Omega_{s}}:=\inf_{\bm{v}_{s}\in\bm{X}_{s}}\frac{(\bm{f}_{s},\bm{v}_{s})_{\Omega_{s}}}{\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}},\quad\|f_{d}\|_{-1,\Omega_{d}}:=\inf_{\psi\in
X_{d}}\frac{(f_{d},\psi)_{\Omega_{d}}}{\|\nabla\psi\|_{0,\Omega_{d}}}.$

Hence, by using Hölder inequality and Young’s inequality, we have

$\rho(\bm{f}_{s},\bm{u}_{s})_{\Omega_{s}}+(f_{d},\phi_{f})_{\Omega_{d}}\leqslant\frac{\rho\nu}{2}\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}+\frac{\sigma
k_{f}}{2\mu}\|\nabla\phi_{f}\|_{0,\Omega_{d}}^{2}+\frac{\rho}{2\nu}\|\bm{f}_{s}\|_{-1,\Omega_{s}}^{2}+\frac{\mu}{2\sigma
k_{f}}\|f_{d}\|_{-1,\Omega_{d}}^{2}.$ (3.23)

In addition, based on the imbedding result (2.7), the standard trace theorem
[1], Hölder inequality, Korn’s inequality, Young’s inequality, (3.15) and the
following inverse inequality [39]: for any polynomial $\bm{v}$ on $K$,

$\|(\nabla\bm{v})\bm{n}\|_{0,e}\leqslant
Ch^{-1/2}|\bm{v}|_{1,K},\quad\forall\,e\subset\partial K,\
\forall\,K\in\mathcal{R}_{h},$

we have

$\begin{split}&\kappa\left\langle\frac{\partial\bm{u}_{d,h}}{\partial\bm{n}_{d}},\Pi_{s,h}\bm{u}_{s}\right\rangle_{\Gamma}\leqslant\kappa\left\|(\nabla\bm{u}_{d,h})\bm{n}_{d}\right\|_{\left(\bm{H}_{00}^{1/2}(\Gamma)\right)^{\prime}}\left\|\Pi_{s,h}\bm{u}_{s}\right\|_{\bm{H}_{00}^{1/2}(\Gamma)}\\\
&\quad\leqslant
C\kappa\left\|(\nabla\bm{u}_{d,h})\bm{n}_{d}\right\|_{0,\Gamma}\left\|\Pi_{s,h}\bm{u}_{s}\right\|_{1,\Omega_{s}}\leqslant
C\kappa\left(\sum_{e\in\Gamma}\|(\nabla\bm{u}_{d,h})\bm{n}\|_{0,e}^{2}\right)^{1/2}\left\|\mathbb{D}(\bm{u}_{s})\right\|_{0,\Omega_{s}}\\\
&\quad\leqslant C\kappa
h^{-1/2}\left\|\mathbb{D}(\bm{u}_{d,h})\right\|_{0,\Omega_{d}}\left\|\mathbb{D}(\bm{u}_{s})\right\|_{0,\Omega_{s}}\leqslant\frac{C\kappa^{2}}{\rho\nu
h}\left\|\mathbb{D}(\bm{u}_{d,h})\right\|_{0,\Omega_{d}}^{2}+\frac{\rho\nu}{2}\left\|\mathbb{D}(\bm{u}_{s})\right\|_{0,\Omega_{s}}^{2}.\end{split}$
(3.24)

Finally, with the assumptions that $h$ is small enough such that
$Ch^{1/2}M_{s}<\rho\nu/2$, and $\kappa$ satisfies $0<\kappa\leqslant C\rho\nu
h$, gathering (3.19), (3.21), (3.22), (3.23) and (3.24) yields (3.18), where

$\mathcal{C}^{2}=\rho\nu^{-1}\|\bm{f}_{s}\|_{-1,\Omega_{s}}^{2}+\mu\sigma^{-1}k_{f}^{-1}\|f_{d}\|_{-1,\Omega_{d}}^{2}.$

##  4 Existence and global uniqueness of the solution

In this section, we shall use the technique of the Galerkin method to verify
that problem (3.13) has at least one weak solution, and then we can prove the
global uniqueness of the weak solution due to the a priori estimate (3.18)
obtained in Section 3.

###  4.1 The solvability of the conforming Galerkin approximation problem

We denote by $\underline{\bm{W}}_{h}$ the product finite element space
$\bm{V}_{s,h}\times X_{d,h}\times X_{d,h}$. Then, we consider the following
conforming Galerkin approximation problem of (3.13): to find
$\underline{\bm{u}}_{h}=(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})\in\underline{\bm{W}}_{h}$
such that
$\forall\,\underline{\bm{v}}_{h}=(\bm{v}_{s,h},\psi_{m,h},\psi_{f,h})\in\underline{\bm{W}}_{h}$,

$a(\underline{\bm{u}}_{h},\underline{\bm{v}}_{h})+b(\bm{u}_{s,h};\bm{u}_{s,h},\bm{v}_{s,h})=\rho(\bm{f}_{s},\bm{v}_{s,h})_{\Omega_{s}}+(f_{d},\psi_{f,h})_{\Omega_{d}}.$
(4.1)

To show the solvability of (4.1), stemming from Theorem 3.2, we shall consider
the following constructed coupled discrete system: to find
$(\underline{\bm{u}}_{h},\bm{u}_{d,h})\in\underline{\bm{W}}_{h}\times\bm{X}_{d,h}$
such that

$\displaystyle
a(\underline{\bm{u}}_{h},\underline{\bm{v}}_{h})+b(\bm{u}_{s,h};\bm{u}_{s,h},\bm{v}_{s,h})=\rho(\bm{f}_{s},\bm{v}_{s,h})_{\Omega_{s}}+(f_{d},\psi_{f,h})_{\Omega_{d}}$
$\displaystyle\forall\,\underline{\bm{v}}_{h}\in\underline{\bm{W}}_{h},$ (4.2)
$\displaystyle
2\xi\left(\mathbb{D}(\bm{u}_{d,h}),\mathbb{D}(\bm{v}_{d,h})\right)_{\Omega_{d}}+\left((\bm{\beta}\cdot\nabla)\bm{u}_{d,h},\bm{v}_{d,h}\right)_{\Omega_{d}}-\xi\left\langle\frac{\partial\bm{u}_{d,h}}{\partial\bm{n}_{d}},\bm{v}_{d,h}\right\rangle_{\Gamma}=0$
$\displaystyle\forall\,\bm{v}_{d,h}\in\bm{X}_{d,h},$ (4.3)
$\displaystyle\bm{u}_{d,h}|_{\Gamma}=\bm{u}_{s,h}|_{\Gamma},$ (4.4)

where $\xi>0$ is a positive constant specified later,
$\bm{\beta}:=\mathcal{L}(\bm{u}_{s,h}|_{\Gamma})$, and
$\nabla\cdot\bm{\beta}=0$ in $\Omega_{d}$ by the definition of the lifting
operator $\mathcal{L}$. Furthermore, it follows the head statements of Section
3.4 that if
$(\underline{\bm{u}}_{h},\bm{u}_{d,h})\in\underline{\bm{W}}_{h}\times\bm{X}_{d,h}$
is a solution of problem (4.2)–(4.4), then
$\underline{\bm{u}}_{h}\in\underline{\bm{W}}_{h}$ will solve (4.1).

######  Lemma 4.1.

If problem (4.1) exists a possible solution
$\underline{\bm{u}}_{h}=(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})\in\underline{\bm{W}}_{h}$,
we have the following a priori estimate that

$\rho\nu\|\mathbb{D}(\bm{u}_{s,h})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\phi_{m,h}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f,h}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2},$
(4.5)

where $\mathcal{C}^{2}$ is defined in Theorem 3.2.

######  Proof 3.

The proof is quite close to the proof of Theorem 3.2. To avoid repeating, we
just present the differences. Since
$\bm{\beta}=\mathcal{L}(\bm{u}_{s,h}|_{\Gamma})$, we have

$\left((\bm{\beta}\cdot\nabla)\bm{u}_{d,h},\bm{u}_{d,h}\right)_{\Omega_{d}}=-\frac{1}{2}\left\langle\bm{u}_{s,h}\cdot\bm{n}_{s},\left|\bm{u}_{s,h}\right|^{2}\right\rangle_{\Gamma},$

and thus the term (3.22) vanishes here. As a result, some subtle differences
occur to the following estimate:

$\xi\left\langle\frac{\partial\bm{u}_{d,h}}{\partial\bm{n}_{d}},\bm{u}_{s,h}\right\rangle_{\Gamma}\leqslant\frac{C\xi^{2}}{2\rho\nu
h}\|\mathbb{D}(\bm{u}_{d,h})\|_{0,\Omega_{d}}^{2}+\rho\nu\|\mathbb{D}(\bm{u}_{s,h})\|_{0,\Omega_{s}}^{2}.$

Thus, for any given mesh parameter $h$ , the result (4.5) follows the
assumption that $0<\xi\leqslant C\rho\nu h$.

Now we start to verify the solvability of problem (4.2)–(4.4). For any
$\widehat{\bm{v}}_{h}\in\bm{V}_{h}$, following the similar technique proposed
in [32], we denote
$\bm{v}_{h}^{s}=\widehat{\bm{v}}_{h}|_{\Omega_{s}}\in\bm{V}_{s,h}$,
$\bm{v}_{h}^{d}=\widehat{\bm{v}}_{h}|_{\Omega_{d}}\in\bm{X}_{d,h}$, and

$\widehat{\widetilde{\bm{v}}}_{h}:=\left\\{\begin{array}[]{cl}\bm{v}_{h}^{s}&\mathrm{in}~{}\Omega_{s},\\\
\bm{v}_{d}^{0}&\mathrm{in}~{}\Omega_{d},\end{array}\right.$

where $\bm{v}_{d}^{0}=\mathcal{L}(\bm{v}_{h}^{s}|_{\Gamma})\in\bm{X}_{d}$ with
$\nabla\cdot\bm{v}_{d}^{0}=0$ in the domain $\Omega_{d}$. Then, we denote by
$\underline{\bm{Z}}_{h}$ the product space $\bm{V}_{h}\times X_{d,h}\times
X_{d,h}$, and define a mapping
$\mathcal{F}_{h}:\underline{\bm{Z}}_{h}\rightarrow\underline{\bm{Z}}_{h}$ as:
for each
$\widehat{\underline{\bm{v}}}_{h}=(\widehat{\bm{v}}_{h},\psi_{m,h},\psi_{f,h})\in\underline{\bm{Z}}_{h}$
such that
$\forall\,\widehat{\underline{\bm{w}}}_{h}=(\widehat{\bm{w}}_{h},\varphi_{m,h},\varphi_{f,h})\in\underline{\bm{Z}}_{h}$,

$\begin{split}&\left(\mathcal{F}_{h}\widehat{\underline{\bm{v}}}_{h},\widehat{\underline{\bm{w}}}_{h}\right)_{\underline{\bm{Z}}_{h}}:=2\rho\nu\left(\mathbb{D}(\bm{v}_{h}^{s}),\mathbb{D}(\bm{w}_{h}^{s})\right)_{\Omega_{s}}+\frac{k_{m}}{\mu}\left(\nabla\psi_{m,h},\nabla\varphi_{m,h}\right)_{\Omega_{d}}+\frac{k_{f}}{\mu}\left(\nabla\psi_{f,h},\nabla\varphi_{f,h}\right)_{\Omega_{d}}\\\
&\quad+\frac{\sigma
k_{m}}{\mu}\left(\psi_{m,h}-\psi_{f,h},\varphi_{m,h}\right)_{\Omega_{d}}+\frac{\sigma
k_{m}}{\mu}\left(\psi_{f,h}-\psi_{m,h},\varphi_{f,h}\right)_{\Omega_{d}}+2\xi\left(\mathbb{D}(\bm{v}_{h}^{d}),\mathbb{D}(\bm{w}_{h}^{d})\right)_{\Omega_{d}}\\\
&\quad+\left\langle\psi_{f,h},\bm{w}_{h}^{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}-\left\langle\varphi_{f,h},\bm{v}_{h}^{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}+\sum_{i=1}^{N-1}\frac{\alpha\rho\nu}{\sqrt{k_{f}}}\left\langle\bm{v}_{h}^{s}\cdot\bm{\tau}_{i},\bm{w}_{h}^{s}\cdot\bm{\tau}_{i}\right\rangle_{\Gamma}-\xi\left\langle\frac{\partial\bm{v}_{h}^{d}}{\partial\bm{n}_{d}},\bm{w}_{h}^{d}\right\rangle_{\Gamma}\\\
&\quad+\left((\widehat{\widetilde{\bm{v}}}_{h}\cdot\nabla)\widehat{\bm{v}}_{h},\widehat{\bm{w}}_{h}\right)_{\Omega}+\frac{1}{2}\left((\nabla\cdot\widehat{\widetilde{\bm{v}}}_{h})\widehat{\bm{v}}_{h},\widehat{\bm{w}}_{h}\right)_{\Omega}-\rho(\bm{f}_{s},\bm{w}_{h}^{s})_{\Omega_{s}}-(f_{d},\varphi_{f,h})_{\Omega_{d}}.\end{split}$
(4.6)

Clearly, $\mathcal{F}_{h}$ define a mapping from $\underline{\bm{Z}}_{h}$ into
itself, and a zero of $\mathcal{F}_{h}$ is a solution of the coupled system
(4.2)–(4.4). Further we introduce the Brouwer’s fixed-point theorem:

######  Lemma 4.2.

[18] Let $U$ be a nonempty, convex, and compact subset of a normed vector
space and let $F$ be a continuous mapping from $U$ into $U$. Then $F$ has at
least one fixed point.

Based on Lemma 4.1, we define $\underline{\bm{U}}_{h}$ a subset of
$\underline{\bm{Z}}_{h}$ as:

$\underline{\bm{U}}_{h}:=\left\\{\widehat{\underline{\bm{v}}}_{h}\in\underline{\bm{Z}}_{h}:\,\rho\nu\|\mathbb{D}(\bm{v}_{h}^{s})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\psi_{m,h}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\psi_{f,h}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2}\right\\},$

where $\mathcal{C}^{2}$ is defined in Theorem 3.2. Then, taking
$\widehat{\underline{\bm{w}}}_{h}=\widehat{\underline{\bm{v}}}_{h}\in\underline{\bm{U}}_{h}$
in (4.6), and following the steps of proving Theorem 3.2 and Lemma 4.1, we
obtain that

$\left(\mathcal{F}_{h}\widehat{\underline{\bm{v}}}_{h},\widehat{\underline{\bm{v}}}_{h}\right)_{\underline{\bm{Z}}_{h}}\geqslant
0,\quad\forall\,\widehat{\underline{\bm{v}}}_{h}\in\underline{\bm{U}}_{h}.$

Hence, it follows Lemma 4.2 that there is at least one zero of
$\mathcal{F}_{h}$ in the ball $\underline{\bm{U}}_{h}$ centered at the origin.

Gathering all the above results, we conclude the following theorem:

######  Theorem 4.3.

For any given mesh parameter $h>0$, the conforming Galerkin approximation
problem (4.1) has at least one solution
$(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})\in\bm{V}_{s,h}\times X_{d,h}\times
X_{d,h}$ with the following estimate:

$\rho\nu\|\mathbb{D}(\bm{u}_{s,h})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\phi_{m,h}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f,h}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2},$
(4.7)

where $\mathcal{C}^{2}$ is defined in Theorem 3.2.

###  4.2 Existence and global uniqueness

It stems from Theorem 4.3 and the conforming property
$\bm{V}_{s,h}\times\bm{X}_{d,h}\times\bm{X}_{d,h}\subset\bm{V}_{s}\times\bm{X}_{d}\times\bm{X}_{d}$
that there exists a 3-tuple function $(\bm{u}_{s},\phi_{m},\phi_{f})$ in the
Hilbert space $\underline{\bm{W}}=\bm{V}_{s}\times\bm{X}_{d}\times\bm{X}_{d}$,
and a uniformly bounded subsequence
$\left\\{(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})\right\\}_{h>0}$ such that

$\lim_{h\rightarrow
0}(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})=(\bm{u}_{s},\phi_{m},\phi_{f})\quad\mathrm{weakly~{}in}~{}\underline{\bm{W}}.$
(4.8)

Furthermore, the Sobolev imbedding implies that the above convergence results
are strong in $L^{q}(\Omega)$ for any $1\leqslant q<6$ whenever $N=2$ or $3$.
In particular, by extracting another subsequence, still denoted by $h$, we
obtain

$\lim_{h\rightarrow
0}(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})=(\bm{u}_{s},\phi_{m},\phi_{f})\quad\mathrm{strongly~{}in}~{}\bm{L}^{2}(\Omega_{s})\times
L^{2}(\Omega_{d})\times L^{2}(\Omega_{d}).$ (4.9)

######  Theorem 4.4.

If the data satisfies that

$\mathcal{N}\left(\rho^{-1}\nu^{-2}\|\bm{f}_{s}\|_{-1,\Omega_{s}}+\rho^{-3/2}\nu^{-3/2}\mu^{1/2}\sigma^{-1/2}k_{f}^{-1/2}\|f_{d}\|_{-1,\Omega_{d}}\right)<1,$
(4.10)

the problem (3.13) then admits a unique solution
$(\bm{u}_{s},\phi_{m},\phi_{f})$ in $\bm{V}_{s}\times X_{d}\times X_{d}$ such
that

$\rho\nu\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\phi_{m}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2},$

where $\mathcal{C}^{2}$ is defined in Theorem 3.2.

######  Proof 4.

For $(\bm{u}_{s},\phi_{m},\phi_{f})$ denoted as $\underline{\bm{u}}$ and the
subsequences $\\{(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})\\}_{h>0}$ denoted as
$\\{\underline{\bm{u}}_{h}\\}_{h>0}$ defined in (4.8) and (4.9), we can easily
obtain

$\lim_{h\rightarrow
0}\\{a_{s}(\underline{\bm{u}}_{h},\underline{\bm{v}})+a_{d}(\underline{\bm{u}}_{h},\underline{\bm{v}})\\}=a_{s}(\underline{\bm{u}},\underline{\bm{v}})+a_{d}(\underline{\bm{u}},\underline{\bm{v}}),\quad\forall\,\underline{\bm{v}}=(\bm{v}_{s},\psi_{m},\psi_{f})\in\underline{\bm{W}},$
(4.11)

because of (4.8). Then, for the trace bilinear term
$a_{\Gamma}(\underline{\bm{u}}_{h},\underline{\bm{v}})$, it follows from
Hölder inequality, (3.2), and (3.4)–(3.6) that

$\displaystyle a_{\Gamma}(\underline{\bm{u}}_{h},\underline{\bm{v}})$
$\displaystyle=\left\langle\phi_{f,h}-\phi_{f},\bm{v}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}-\left\langle\psi_{f},(\bm{u}_{s,h}-\bm{u}_{s})\cdot\bm{n}_{s}\right\rangle_{\Gamma}+\sum_{i=1}^{N-1}\frac{\alpha\rho\nu}{\sqrt{k_{f}}}\left\langle(\bm{u}_{s,h}-\bm{u}_{s})\cdot\bm{\tau}_{i},\bm{v}_{s}\cdot\bm{\tau}_{i}\right\rangle_{\Gamma}$
$\displaystyle\quad+\left\langle\phi_{f},\bm{v}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}-\left\langle\psi_{f},\bm{u}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}+\sum_{i=1}^{N-1}\frac{\alpha\rho\nu}{\sqrt{k_{f}}}\left\langle\bm{u}_{s}\cdot\bm{\tau}_{i},\bm{v}_{s}\cdot\bm{\tau}_{i}\right\rangle_{\Gamma}$
$\displaystyle\leqslant
C\|\phi_{f,h}-\phi_{f}\|_{0,\Omega_{d}}^{1/2}\|\nabla(\phi_{f,h}-\phi_{f})\|_{0,\Omega_{d}}^{1/2}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}+C\|\bm{u}_{s,h}-\bm{u}_{s}\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s,h}-\bm{u}_{s})\|_{0,\Omega_{s}}^{1/2}\|\nabla\psi_{f}\|_{0,\Omega_{d}}$
$\displaystyle\quad+C\|\bm{u}_{s,h}-\bm{u}_{s}\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s,h}-\bm{u}_{s})\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}+a_{\Gamma}(\underline{\bm{u}},\underline{\bm{v}}).$

Since the uniform boundedness of $(\bm{u}_{s,h},\phi_{m,h},\phi_{f,h})$ shown
in (4.7) and the strong $L^{2}$-convergence result (4.9), we derive that

$\lim_{h\rightarrow
0}a_{\Gamma}(\underline{\bm{u}}_{h},\underline{\bm{v}})=a_{\Gamma}(\underline{\bm{u}},\underline{\bm{v}}),\quad\forall\,\underline{\bm{v}}\in\underline{\bm{W}}.$
(4.12)

In addition, for the limit of the trilinear form
$b(\bm{u}_{s,h};\bm{u}_{s,h},\bm{v}_{s})$, gathering (2.8) and Lemma 3.1
yields

$\displaystyle\left|b(\bm{u}_{s,h};\bm{u}_{s,h},\bm{v}_{s})-b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s})\right|$
$\displaystyle\leqslant\left|\left((\bm{u}_{s,h}-\bm{u}_{s})\cdot\nabla)\bm{u}_{s,h},\bm{v}_{s}\right)_{\Omega_{s}}\right|+\frac{1}{2}\left|\left(\nabla\cdot\bm{u}_{s,h},(\bm{u}_{s,h}-\bm{u}_{s})\cdot\bm{v}_{s}\right)_{\Omega_{s}}\right|$
$\displaystyle\quad+\left|\left((\bm{u}_{s}\cdot\nabla)(\bm{u}_{s,h}-\bm{u}_{s}),\bm{v}_{s}\right)_{\Omega_{s}}\right|+\frac{1}{2}\left|\left(\nabla\cdot(\bm{u}_{s,h}-\bm{u}_{s}),\bm{u}_{s}\cdot\bm{v}_{s}\right)_{\Omega_{s}}\right|$
$\displaystyle\leqslant\underbrace{C\|\bm{u}_{s,h}-\bm{u}_{s}\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s,h}-\bm{u}_{s})\|_{0,\Omega_{s}}^{1/2}\|\mathbb{D}(\bm{u}_{s,h})\|_{0,\Omega_{s}}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}}_{I_{h}}$
$\displaystyle\quad+\underbrace{\left|\left((\bm{u}_{s}\cdot\nabla)(\bm{u}_{s,h}-\bm{u}_{s}),\bm{v}_{s}\right)_{\Omega_{s}}\right|+\frac{1}{2}\left|\left(\nabla\cdot(\bm{u}_{s,h}-\bm{u}_{s}),\bm{u}_{s}\cdot\bm{v}_{s}\right)_{\Omega_{s}}\right|}_{II_{h}}.$

We can easily obtain that $\lim_{h\rightarrow 0}I_{h}=0$ by (4.7) and (4.9),
and $\lim_{h\rightarrow 0}II_{h}=0$ by (4.8). Therefore the following limit
holds:

$\lim_{h\rightarrow
0}b(\bm{u}_{s,h};\bm{u}_{s,h},\bm{v}_{s})=b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s}),\quad\forall\,\bm{v}_{s}\in\bm{V}_{s}.$
(4.13)

It follows from (4.11)–(4.13) that

$a(\underline{\bm{u}},\underline{\bm{v}})+b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s})=\rho(\bm{f}_{s},\bm{v}_{s})_{\Omega_{s}}+(f_{d},\psi_{f})_{\Omega_{d}},\quad\forall\,\underline{\bm{v}}\in\underline{\bm{W}},$

which implies that
$\underline{\bm{u}}=(\bm{u}_{s},\phi_{m},\phi_{f})\in\underline{\bm{W}}$ is a
solution of (3.13).

Finally, we assume that there are two solutions
$\underline{\bm{u}}^{1},\underline{\bm{u}}^{2}\in\underline{\bm{W}}$ to
(3.13). Then, there differences
$\bm{\mathrm{e}}_{s}=\bm{u}_{s}^{1}-\bm{u}_{s}^{2}$,
$\mathrm{e}_{m}=\phi_{m}^{1}-\phi_{m}^{2}$ and
$\mathrm{e}_{f}=\phi_{f}^{1}-\phi_{f}^{2}$ satisfy that
$\forall\,(\bm{v}_{s},\psi_{m},\psi_{f})\in\bm{V}_{s}\times X_{d}\times
X_{d}$,

$\begin{split}&2\rho\nu\left(\mathbb{D}(\bm{\mathrm{e}}_{s}),\mathbb{D}(\bm{v}_{s})\right)_{\Omega_{s}}+\frac{k_{m}}{\mu}\left(\nabla\mathrm{e}_{m},\nabla\psi_{m}\right)_{\Omega_{d}}+\frac{k_{f}}{\mu}\left(\nabla\mathrm{e}_{f},\nabla\psi_{f}\right)_{\Omega_{d}}+\left((\bm{\mathrm{e}}_{s}\cdot\nabla)\bm{u}_{s}^{1},\bm{v}_{s}\right)_{\Omega_{s}}+\left((\bm{u}_{s}^{2}\cdot\nabla)\bm{\mathrm{e}}_{s},\bm{v}_{s}\right)_{\Omega_{s}}\\\
&+\frac{\sigma
k_{m}}{\mu}\left(\mathrm{e}_{m}-\mathrm{e}_{f},\psi_{m}\right)_{\Omega_{d}}+\frac{\sigma
k_{m}}{\mu}\left(\mathrm{e}_{f}-\mathrm{e}_{m},\psi_{f}\right)_{\Omega_{d}}+\left\langle\mathrm{e}_{f},\bm{v}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}-\left\langle\psi_{f},\bm{\mathrm{e}}_{s}\cdot\bm{n}_{s}\right\rangle_{\Gamma}+\sum_{i=1}^{N-1}\left\langle\bm{\mathrm{e}}_{s}\cdot\bm{\tau}_{i},\bm{v}_{s}\cdot\bm{\tau}_{i}\right\rangle_{\Gamma}=0.\end{split}$
(4.14)

Taking $\bm{v}_{s}=\bm{\mathrm{e}}_{s}$, $\psi_{m}=\mathrm{e}_{m}$ and
$\psi_{f}=\mathrm{e}_{f}$ in (4.14) and using another version of (3.7), which
is

$\left|\left((\bm{v}_{s}\cdot\nabla)\bm{w}_{s},\bm{z}_{s}\right)_{\Omega_{s}}\right|\leqslant\mathcal{N}\|\mathbb{D}(\bm{v}_{s})\|_{0,\Omega_{s}}\|\mathbb{D}(\bm{w}_{s})\|_{0,\Omega_{s}}\|\mathbb{D}(\bm{z}_{s})\|_{0,\Omega_{s}},\quad\forall\,\bm{v}_{s},\bm{w}_{s},\bm{z}_{s}\in\bm{X}_{s},$
(4.15)

we obtain

$\left[2\rho\nu-\mathcal{N}(\|\mathbb{D}(\bm{u}_{s}^{1})\|_{0,\Omega_{s}}+\|\mathbb{D}(\bm{u}_{s}^{2})\|_{0,\Omega_{s}})\right]\|\mathbb{D}(\bm{\mathrm{e}}_{s})\|_{0,\Omega_{s}}^{2}+\frac{k_{m}}{\mu}\|\nabla\mathrm{e}_{m}\|_{0,\Omega_{d}}^{2}+\frac{k_{f}}{\mu}\|\nabla\mathrm{e}_{f}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{m}}{\mu}\|\mathrm{e}_{m}-\mathrm{e}_{f}\|_{0,\Omega_{d}}^{2}\leqslant 0.$
(4.16)

Hence, if we assume (4.10) holds, then (4.16) shows the unique solution to
(3.13) based on a priori estimate (3.18).

In order to prove the existence and uniqueness of the solution
$(\underline{\bm{u}},p_{s})\in\underline{\bm{Y}}\times Q_{s}$ to the model
problem (2.9), we shall use the inf-sup condition (3.12) and the
Babus̆ka–Brezzi’s theory [5, 9, 25, 46].

######  Theorem 4.5.

Under the assumption (4.10) of Theorem 4.4, the model problem (2.9) admits a
unique solution $(\underline{\bm{u}},p_{s})\in\underline{\bm{Y}}\times Q_{s}$
such that

$\displaystyle\rho\nu\|\mathbb{D}(\bm{u}_{s})\|_{0,\Omega_{s}}^{2}+\frac{2\sigma
k_{m}}{\mu}\|\nabla\phi_{m}\|_{0,\Omega_{d}}^{2}+\frac{\sigma
k_{f}}{\mu}\|\nabla\phi_{f}\|_{0,\Omega_{d}}^{2}\leqslant\mathcal{C}^{2},$
(4.17)
$\displaystyle\|p_{s}\|_{0,\Omega_{s}}\leqslant\beta_{0}^{-1}\left(2\rho^{1/2}\nu^{1/2}\mathcal{C}+\rho^{-1}\nu^{-1}\mathcal{N}\mathcal{C}^{2}+\rho\|\bm{f}_{s}\|_{-1,\Omega_{s}}+\|f_{d}\|_{-1,\Omega_{d}}\right),$
(4.18)

where $\mathcal{C}^{2}$ is defined in Theorem 3.2.

######  Proof 5.

For the solution
$\underline{\bm{u}}=(\bm{u}_{s},\phi_{m},\phi_{f})\in\underline{\bm{W}}$ to
(3.13), the following mapping:

$\underline{\bm{v}}=(\bm{v}_{s},\psi_{m},\psi_{f})\in\underline{\bm{Y}}\;\mapsto\;a(\underline{\bm{u}},\underline{\bm{v}})+b(\bm{u}_{s};\bm{u}_{s},\bm{v}_{s})-\rho(\bm{f}_{s},\bm{v}_{s})_{\Omega_{s}}-(f_{d},\psi_{f})_{\Omega_{d}}$

defines an element $L(\underline{\bm{v}})$ of the dual space
$\underline{\bm{Y}}^{\prime}$, and furthermore, $L$ vanishes on
$\underline{\bm{W}}$. As a result, the inf-sup condition (3.12) implies that
there exists exactly one $p_{s}\in Q_{s}$ such that

$L(\underline{\bm{v}})=d(\bm{v}_{s},p_{s}),\quad\forall\,\underline{\bm{v}}=(\bm{v}_{s},\psi_{m},\psi_{f})\in\underline{\bm{Y}}.$
(4.19)

Therefore the fact $\bm{u}_{s}\in\bm{X}_{s}$ and (4.19) show that the model
problem (2.9) admits a unique solution
$(\underline{\bm{u}},p_{s})\in\underline{\bm{Y}}\times Q_{s}$. Finally, the
result (4.18) is a straightforward application of the inf-sup condition (3.12)
with the help of (4.15) and (4.17).

## References

  * [1] R. Adams, J. Fournier, Sobolev Spaces, Acadamics Press, New York, 2003.
  * [2] M. A. Al Mahbub, X.-M. He, N. J. Nasu, C. Qiu, H. Zheng, Coupled and decoupled stabilized mixed finite element methods for nonstationary dual–porosity–Stokes fluid flow model, Int. J. Numer. Methods. Eng. 120 (6) (2019) 803–833.  doi:10.1002/nme.6158.
  * [3] M. A. Al Mahbub, F. Shi, N. J. Nasu, Y. Wang, H. Zheng, Mixed stabilized finite element method for the stationary Stokes–dual–permeability fluid flow model, Comput. Methods Appl. Mech. Engrg. 358 (2020) 112616.  doi:10.1016/j.cma.2019.112616.
  * [4] T. Aubin, Nonlinear Analysis on Manifolds. Monge-Ampère Equations, Grundlehren der mathematischen Wissenschaften 252, Springer-Verlag New York, NY, 1982.  doi:10.1007/978-1-4612-5734-9.
  * [5] I. Babus̆ka, The finite element method with Lagrangian multipliers, Numer. Math. 20 (1973) 179–192.  doi:10.1007/BF01436561.
  * [6] L. Badea, M. Discacciati, A. Quarteroni, Numerical analysis of the Navier-Stokes/Darcy coupling, Numer. Math. 115 (2010) 195–227.  doi:10.1007/s00211-009-0279-6.
  * [7] G. Beavers, D. Joseph, Boundary conditions at a naturally permeable wall, Journal of Fluid Mechanics 30 (1) (1967) 197–207.  doi:10.1017/S0022112067001375.
  * [8] J. Bergh, J. Löfström, Interpolation Spaces: An Introduction, Grundlehren der mathematischen Wissenschaften 223, Springer-Verlag Berlin Heidelberg, Berlin, 1976.  doi:10.1007/978-3-642-66451-9.
  * [9] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, RAIRO Anal. Numer. R2 (1974) 129–151.  doi:10.1051/m2an/197408R201291.
  * [10] M. Cai, M. Mu, J. Xu, Numerical solution to a mixed Navier–Stokes/Darcy model by the two-grid approach, SIAM J. Numer. Anal. 47 (5) (2009) 3325–3338.  doi:10.1137/080721868.
  * [11] L. Cao, Y. He, J. Li, D. Yang, Decoupled modified characteristic FEMs for fully evolutionary Navier-Stokes-Darcy model with the Beavers-Joseph interface condition, Journal of Computational and Applied Mathematics 383 (2021) 113128.  doi:10.1016/j.cam.2020.113128.
  * [12] Y. Cao, Y. Chu, X. He, M. Wei, Decoupling the Stationary Navier-Stokes-Darcy System with the Beavers-Joseph-Saffman Interface Condition, Abstract and Applied Analysis 2013 (2013) 136483.  doi:10.1155/2013/136483.
  * [13] Y. Cao, M. Gunzburger, X. He, X. Wang, Robin–Robin domain decomposition methods for the steady-state Stokes–Darcy system with the Beavers–Joseph interface condition, Numer. Math. 117 (4) (2011) 601–629.  doi:10.1007/s00211-011-0361-8.
  * [14] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, W. Zhao, Finite element approximations for Stokes–Darcy flow with Beavers–Joseph interface conditions, SIAM. J. Numer. Anal. 47 (6) (2010) 4239–4256.  doi:10.1137/080731542.
  * [15] Y. Cao, M. Gunzburger, F. Hua, X. Wang, Coupled Stokes–Darcy model with Beavers–Joseph interface boundary condition, Commun. Math. Sci. 8 (1) (2010) 1–25.  doi:10.4310/CMS.2010.v8.n1.a2.
  * [16] W. Chen, M. Gunzburger, S. Dong, X. Wang, Efficient and long-time accurate second-order methods for Stokes-Darcy system, SIAM J. Numer. Anal. 51 (5) (2013) 2563–2584.  doi:doi:10.1137/120897705.
  * [17] P. Chidyagwai, B. Rivière, On the solution of the coupled Navier–Stokes and Darcy equations, Comput. Methods Appl. Mech. Engrg. 198 (47) (2009) 3806–3820.  doi:10.1016/j.cma.2009.08.012.
  * [18] D. Cioranescu, V. Girault, K. R. Rajagopal, Mechanics and Mathematics of Fluids of the Differential Type, Advances in Mechanics and Mathematics, vol. 35. Springer International Publishing, Cham, Switzerland, 2016.  doi:10.1007/978-3-319-39330-8.
  * [19] C. D’Angelo, P. Zunino, Robust numerical approximation of coupled Stokes’ and Darcy’s flows applied to vascular hemodynamics and biochemical transport, ESAIM: Mathematical Modelling and Numerical Analysis 45 (3) (2011) 447–476.  doi:10.1051/m2an/2010062.
  * [20] M. Discacciati, Domain decomposition methods for the coupling of surface and groundwater flows, Ph.D. thesis, Ecole Polytechnique Fédérale de Lausanne, Switzerland, 2004.  doi:10.5075/epfl-thesis-3117.
  * [21] M. Discacciati, A. Quarteroni, Navier-Stokes/Darcy Coupling: Modeling, Analysis, and Numerical Approximation, Rev. Mat. Complut. 22 (2) (2009) 315–426.  doi:10.5209/rev_REMA.2009.v22.n2.16263.
  * [22] M. Discacciati, R. Oyarzúa, A conforming mixed finite element method for the Navier–Stokes/Darcy coupled problem, Numer. Math. 135 (2017) 571–606.  doi:10.1007/s00211-016-0811-4.
  * [23] J. Fang, P. Huang, Y. Qin, A two-level finite element method for the steady-state Navier-Stokes/Darcy model, J. Korean Math. Soc. 57 (4) (2020) 915–933.  doi:10.4134/JKMS.j190449.
  * [24] L. Gao, J. Li, A decoupled stabilized finite element method for the dual-porosity-Navier–Stokes fluid flow model arising in shale oil, Numer. Methods Partial Differential Eq. (2020) 1–18.  doi:10.1002/num.22718.
  * [25] V. Girault, P. A. Raviart, Finite element methods for Navier–Stokes equations, Springer-Verlag, Berlin, 1986.  doi:10.1007/978-3-642-61623-5.
  * [26] V. Girault, B. Rivière, DG Approximation of Coupled Navier–Stokes and Darcy Equations by Beaver–Joseph–Saffman Interface Condition, SIAM J. Numer. Anal. 47 (3) (2009) 2052–2089.  doi:10.1137/070686081.
  * [27] V. Girault, G. Kanschat, B. Rivière, On the Coupling of Incompressible Stokes or Navier-Stokes and Darcy Flows Through Porous Media. In: J. Ferreira, S. Barbeiro, G. Pena, M. Wheeler(eds) Modelling and Simulation in Fluid Dynamics in Porous Media. Springer Proceedings in Mathematics & Statistics, vol 28. Springer, New York, NY, 2013.  doi:10.1007/978-1-4614-5055-9_1.
  * [28] N. Hanspal, A. Waghode, V. Nassehi, R. Wakeman, Numerical analysis of coupled Stokes/Darcy flows in industrial filtrations, Transport in Porous Media 64 (1) (2006) 73.  doi:10.1007/s11242-005-1457-3.
  * [29] X. He, N. Jiang, C. Qiu, An artificial compressibility ensemble algorithm for a stochastic Stokes-Darcy model with random hydraulic conductivity and interface conditions, International Journal for Numerical Methods in Engineering 121 (4) (2020) 712–739.  doi:10.1002/nme.6241.
  * [30] X. He, J. Li, Y. Lin, J. Ming, A domain decomposition method for the steady-state Navier-Stokes-Darcy model with Beavers-Joseph interface condition, SIAM J. Sci. Comput. 37 (5) (2015), 264–290.  doi:10.1137/140965776.
  * [31] J. Hou, M. Qiu, X. He, C. Guo, M. Wei, B. Bai, A dual–porosity–Stokes model and finite element method for coupling dual–porosity flow and free flow, SIAM J. Sci. Comput. 38 (5) (2016) B710–B739.  doi:10.1137/15M1044072.
  * [32] Y. Hou, S. Pei, On the weak solutions to steady Navier-Stokes equations with mixed boundary conditions, Mathematische Zeitschrift 291 (2019) 47–54.  doi:10.1007/s00209-018-2072-7.
  * [33] I. P. Jones, Low Reynolds number flow past a porous spherical shell, Mathematical Proceedings of the Cambridge Philosophical Society 73 (1) (1973) 231–238.  doi:10.1017/S0305004100047642.
  * [34] R. Li, J. Li, Z. Chen, Y. Gao, A stabilized finite element method based on two local Gauss integrations for a coupled Stokes–Darcy problem, Journal of Computational and Applied Mathematics 292 (15) (2016), 92–104.  doi:10.1016/j.cam.2015.06.014.
  * [35] A. Mikelić, W. Jäger, On the interface boundary condition of Beavers, Joseph, and Saffman, SIAM Journal on Applied Mathematics 60 (4) (2000) 1111–1127.  doi:10.1137/S003613999833678X.
  * [36] M. Mu, X. Zhu, Decoupled schemes for a non-stationary mixed Stokes-Darcy model, Mathematics of Computation 79 (270) (2010) 707–731.  doi:10.1090/S0025-5718-09-02302-3.
  * [37] V. Nassehi, Modelling of combined Navier-Stokes and Darcy flows in crossflow membrane filtration, Chemical Engineering Science 53 (6) (1998) 1253–1265.  doi:10.1016/S0009-2509(97)00443-0.
  * [38] C. Qiu, X. He, J. Li, Y. Lin, A domain decomposition method for the time-dependent Navier-Stokes-Darcy model with Beavers-Joseph interface condition and defective boundary condition, Journal of Computational Physics 411 (2020) 109400.  doi:10.1016/j.jcp.2020.109400.
  * [39] B. Rivière, Discontinuous Galerkin methods for solving elliptic and parabolic equations. Theory and implementation, Front. Appl. Math. 35, SIAM, Philadelphia, 2008.  doi:10.1137/1.9780898717440.
  * [40] H. G. Roos, M. Stynes, L. Tobiska, Robust numerical methods for singularly perturbed differential equations. Convection–Diffusion–Reaction and Flow Problems, 2nd ed., Springer Ser. Comput. Math. 24, Springer, Berlin, 2008.  doi:10.1007/978-3-540-34467-4.
  * [41] P. G. Saffman, On the boundary condition at the surface of a porous medium, Studies in Applied Mathematics 50 (2) (1971) 93–101.  doi:10.1002/sapm197150293.
  * [42] L. Scott, S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (190) (1990) 483–493.  doi:10.1090/S0025-5718-1990-1011446-7.
  * [43] K. Serra, A. C. Reynolds, R. Raghavan, New pressure transient analysis methods for naturally fractured reservoirs, Journal of Petroleum Technology 35 (12) (1983) 2271–2283.  doi:10.2118/10780-PA.
  * [44] L. Shan, J. Hou, W. Yan, J. Chen, Partitioned time stepping method for a dual–porosity–Stokes model, J. Sci. Comput. 79 (1) (2019) 389–413.  doi:10.1007/s10915-018-0879-3.
  * [45] L. Shan, H. Zheng, Partitioned time stepping method for fully evolutionary Stokes–Darcy flow with Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 51 (2) (2013) 813–839.  doi:10.1137/110828095.
  * [46] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, 3rd Edition, North-Holland, Amsterdam, 1984.  doi:10.1090/chel/343.
  * [47] J. E. Warren, P. J. Root, The behavior of naturally fractured reservoirs, Society of Petroleum Engineers Journal 3 (3) (1963) 245–255.  doi:10.2118/426-PA.
  * [48] J. Zhao, T. Zhang, Two-grid finite element methods for the steady Navier-Stokes/Darcy model, East Asian Journal on Applied Mathematics 6 (1) (2016) 60–79.  doi:10.4208/eajam.080215.111215a.
  * [49] L. Zuo, Y. Hou, Numerical analysis for the mixed Navier–Stokes and Darcy Problem with the Beavers–Joseph interface condition, Numer. Methods Partial Differential Eq. 31 (2015) 1009–1030.  doi:10.1002/num.21933.


Analysis of Complex Survival Data: a tutorial using the Shiny MSM.app
application

Gustavo Soutinho1 and Luís Meira-Machado2 1EPIUnit - University of Porto
Rua das Taipas 135, 4050-600 Porto - Portugal<EMAIL_ADDRESS>

2Centre of Mathematics, University of Minho<EMAIL_ADDRESS>
Keywords: R language, Shiny package, Survival analysis, Multi-state models.

The development of applications for obtaining interpretable results in a
simple and summarized manner in multi-state models is a research field with
great potential, namely in terms of using open source tools that can be easily
implemented in biomedical applications. In this tutorial, we introduce
MSM.app, an interactive web application using the Shiny package for the R
language. In the following sections, we present the main functionalities of
the MSM.app and an explanation of the outputs obtained for better
understanding, independent of the statistical knowledge of users.

1\. What is the Shiny MSM.app?

The appearance of the shiny R package allowed to automatically share results
obtained from the R language to be analyzed for users without any prior
knowledge in terms of informatics via the internet [1, 2, 3].

In this context, the MSM.app web application arises with the goal of carrying
out an analysis of multi-state survival data sets. Among its functionalities,
we can highlight the possibility to conduct a traditional survival analysis
regarding the following topics: (i) estimation of survival functions (using
the classical Kaplan-Meier estimator); (ii) comparison of survival functions
between groups; and (iii) use of semi-parametric and parametric regression
models to study the relationship between explanatory variables and survival
time.

The MSM.app can also be used for the analysis of multi-state survival data
that can be seen as a generalization of survival analysis in which survival is
the ultimate outcome of interest but where information is available about
intermediate events that individuals may experience during the study period
[4, 5, 6].

To this statistical analysis, the Shiny MSM.app application combines shiny
with some other packages such as survival [7], mstate [8] and survidm [9].

2\. Features of the Shiny MSM.app

The MSM.app is a web application that was developed to perform dynamic
analysis through a set of dynamic web forms, tables, and graphics.

It is built upon two components: the user-interface scripts for the layout of
the application where the outputs are displayed (ui.R); and the other given by
the server scripts with the instructions of the application (server.R) [10].

Among the main aspects of this Shiny web tool that improve the flexibility and
productivity, we could highlight the easy rendering of the contents without
multiple reloads, the feature to add computed (or processed) outputs from R
scripts, or interactively add reports and visualizations.

The MSM.app also provides integration with other R packages, javascript
libraries or CSS customization, being under the GPL-2 open source license. The
communication between the client and server is done over the normal TCP
connection [11].

3\. How is the Shiny MSM.app different from other applications?

There exists a set of available web tools specifically aimed at carrying out
some parts of multi-state analysis.

Among them, we can stand out the MSDshiny which provides a useful and
streamlined way to plan and power clinical trials with multi-state outcomes
such as a view of the multi-state structure, treatment effects, or the results
of different types of simulations. Another example is MSM-shiny application.
This tool uses a CSV file containing multi-state data and provides the
modeling and comparison of transition hazard models and the prediction of
occupation probabilities [12]. Recently, MSMplus provides a flexible
visualization of the transition probabilities, transition intensities, or
probability of visiting a particular state [13].

After analyzing existing web applications, we have concluded that their use by
non-statisticians has been limited. A possible reason for this is the lack of
friendly software that covers the main goals involving survival analysis and
multistate models on the same platform.

The MSM.app allows users to explore various types of multi-state models and
perform regression inference as well as obtain several predictive measures of
interest, such as the occupation probabilities, the transition probabilities,
and the cumulative incidence functions. Recent methods for checking the Markov
assumption are also implemented.

4\. How to get the Shiny MSM.app

The Shiny MSM.app is available for free open access at the Shiny Apps
repository https://gsoutinho.shinyapps.io/appmsm/.

5\. Structure of the Shiny MSM.app

The web application consists of three parts representing different aspects of
the survival analysis and its extension to complex multi-state models.

The first one allows to perform the survival analysis from mainly of most
common functions of the survival R packages.

The second enables one to obtain some of the main goals of a multi-state
analysis, such as the inference of regression models and the estimation of
transition probabilities, through the survidm and mstate R packages.

Finally, MSM.app also includes local and global statistical tests to check the
Markov assumption for multi-state using the markovMSM package.

6\. Type of input data files and requirements

The MSM.app only requires CSV files as input. By default, in terms of
structure, the values of the data set should be separated by a comma. However,
files that use a semicolon to separate the values are also accepted.

Three examples of data sets are available for consultation in the MSM.app from
which it is possible to check the requirements that the files must have to
pursue the data analysis. Each one corresponds to a specific type of structure
of data: survival, illness-death model, or more general multi-state models.

The first data set corresponds to survival data in patients with acute
myelogenous leukemia [14]. Figure 1 shows some registers and the three columns
that the file must have. “Time1” and “status” correspond to the time to the
event and the status of the censoring, respectively. “x” is the only covariate
in this data that can be used for regression or obtaining the survival
estimates for groups.

Figure 1: Structure of the CSV files for survival analysis (given by the aml
data set).

The second corresponds to a data set from a clinical trial on colon cancer,
which can be modeled using the progressive illness-death model [15]. Figure 2
shows the schematic diagram of transitions involved in the model. Among the
variables of the colonIDM data set, the first four correspond to the times or
the status indicator for the illness-death model. “time1” represents the
sojourn time in the initial state and “Stime” and “event1” and “event” the
corresponding censoring indicates. The other variables could be used for
transition probabilities or to check the effect on the transition intensities
(Figure 3).

Figure 2: Illness-death model for the colon cancer study. Figure 3: Structure
of the CSV files for Illness-death models analysis (given by the colonIDM data
set).

Finally, extensions to progressive processes beyond the three-state illness-
death model are discussed using data from the European Group for Blood and
Marrow Transplantation (EBMT) [4]. The movement of the patients among the six
states can be modelled through the multi-state model with the following six
states: ‘Alive and in remission, no recovery or adverse event’ (State 0);
‘Alive in remission, recovered from the treatment’ (state 1); ‘Alive in
remission, occurrence of the adverse event’ (state 2); ‘Alive, both recovered
and adverse event’ (state 3); ‘Alive, in relapse’ (treatment failure) (state
4) and ‘Dead (treatment failure)’ (state 5). In total there are 12
transitions, three intermediate events given by recovery (Rec), adverse event
(AE) and a combination of the two (AE and Rec), and two absorbing states:
Relapse and Death (Figure 4).

Figure 4: A six-states model for leukemia patients after bone marrow
transplantation.

In terms of format, the CSV file should be in wide format. In the case of the
ebmt4 data (Figure 5), the transition times from the initial to the ultimate
state (‘srv’) or to the intermediate states (‘rec’, ‘ae’, ‘recae’ and ‘rel’),
and the corresponding censoring variables are rec.s’, ‘ae.s’, ‘recae.s’,
‘rel.s’ and ‘srv.s’. Other variables correspond to the covariates for
regression models.

Figure 5: Structure of the CSV files for multi-state models analysis (given by
the ebmt4 data set).

For all these data sets, the names of the variables could be different. In
these cases, the user must correctly indicate the corresponding name in the
web forms of the application. The data files can be accessed at
https://w3.math.uminho.pt/~lmachado/shiny/.

7\. How to select the input files?

On the “input files” page we can find an interactive form. For each type of
data set (“survival data”, “illness-death model”, or “multi-state model”), a
new web form appears below the radio buttons, in which we indicate the times
to the events or the status indicator (Figure 6, left hand side and center).
In case of the more complex multi-state models, it is also necessary to
indicate the transition schema as well as the number of states, for instance
(Figure 6, right hand side).

Figure 6: The event time and status variables for the survival analysis from
the veteran data set (left); an indication of the two event times and their
corresponding status for the illness-death model given by the colonIDM data
(center); and a description of the MSM model through the number and the state
names, the transition schema, and the event times and corresponding status for
the ebmt4 data (right).

After submitting a data set, a table appears on the right hand side of the
page which can be dynamically. This table can be changed using filters or by
searching for specific words in the table. Figure 7 shows a partial view the
veteran data set that can be found in the survival package.

Figure 7: The input file page with the data table with some results for the
veteran data set and three radio buttons for each type of model.

8\. How to carry out a Survival analysis?

From the “survival analysis” button, we can carry out the classical methods
for survival analysis. In this section, we analyze the main aspects of the
outputs displayed on the right sides of the pages: Kaplan-Meier estimator,
Compare survival curves, Cox PH models and parametric models, by using the
veteran data set.

Kaplan-Meier estimator

Non-parametric estimation of the survival function is traditionally performed
using the Kaplan-Meier estimator, which can be desegregated for different
groups of categorical variables. By default, a summary of the estimates is
presented on the right side. Figure 8 shows, respectively, times with events,
the number of individuals at risk at that time, and the survival estimates
with corresponding standard errors and confidence intervals. It is also
possible to display an interactive plot of the survivals with (or without) the
confidence interval intervals (Figure 9).

Figure 8: Summary of the survival estimation for the categorical covariate
“celltype” of the veteran data set using the Kaplan-Meier estimator. Figure 9:
The output of the survival estimation for the categorical covariate “celltype”
of the veteran data set using the Kaplan-Meier estimator. Survival curves for
each group with confidence intervals could be shown at the bottom by choosing
“Yes”.

Compare survival curves

Statistical tests can be used to compare survival rates between groups. The
null hypothesis states that there is no difference in survival between groups.
The log-rank test and the Gehan-Wilcoxon test are the most commonly used. Both
tests are available in the MSM.App. The output of the tests comprises the chi-
squared statistics, degrees of freedom, and the corresponding p-value. The
number of individuals for each group and the number of observed and expected
values are also presented. Results of Figure 10 show significant differences
in the survival curves among groups given by the cell types ($p<$0.0001).

Figure 10: The output of the Log-rank test for the categorical covariate
“celltype” of the veteran data set for testing the null hypothesis of no
difference in survival between the two groups.

Cox PH models The semi-parametric Cox proportional hazards model [16] is
usually used to evaluate the effects of several factors on survival. The
output for the Cox model is updated as the user selects the variables to be
included in the model in the web form. Figure 11 shows, respectively, the
estimate of the coefficients for the covariates “celltype”, Karnofsky
performance score (“karno”) and “age”; the hazard ratio given by $exp(coef)$
and the respective standard error; and the statistical significance of the
model. Result of the global test for the proportional hazards assumption
confirms that this requirement is no fulfilled and consequently a parametric
accelerated failure time (AFT) model must be preferable to be fitted to this
survival data.

Figure 11: Results of the Cox model for the variables “celltype”, “karno” and
“age” of the veteran data set. Result of the test for the proportional hazard
assumption shows that a Cox PH model

Parametric AFT models The idea behind the outputs for parametric survival
models is quite similar to those provided by the Cox PH models. In this case,
six possibilities of distributions are available to model the baseline hazard
function of the models: exponential, weibull, gaussian, logistic, lognormal or
loglogistic. Besides the summary of the models, the Akaike Information
Criterion (AIC) values are also presented to make it easier to compare which
model fits better to assess the effect of several risk factors on survival
time. Figures 12, 13 and 14 show, respectively, the outputs of the AFT models
given by the exponential, weibull and loglogistic distributions. Through the
AUC values we can assume that loglogistic parametric model is preferable. From
this model, we can observe that only “age” and “small cells” do not have a
significant effect on survival.

Figure 12: The outputs of the parametric models page with the results of the
fitted model for the variables “celltype”, “karno” and “age” of the veteran
data set using the exponential distribution. Figure 13: The outputs of the
parametric models page with the results of the fitted model for the variables
“celltype”, “karno” and “age” of the veteran data set using the weibull
distribution. Figure 14: The outputs of the parametric models page with the
results of the fitted model for the variables “celltype”, “karno” and “age” of
the veteran data set using the loglogistic distribution.

9\. How to carry out a illness-death model Analysis?

From the “ IDM-Analysis” button, we can carry out a data analysis of a
progressive illness-death (IDM) model given by two events, three states, and
three transitions. Number of events First, we are interested in having an idea
of the movement of individuals among the three states. It is also possible to
see the proportions of transitions (Figure 15).

Figure 15: The number of transitions among the three states of the colonIDM
data set.

Regression models In the multi-state models, there are so many transition
intensities as there are transitions. For each one, we can then check the
effect of the individual characteristics of the individuals by fitting
separate intensities using semi-parametric Cox proportional hazard regression
models [16]. In the inference of the regression models, we must take into
account if the Markov condition, in which the past and future are independent
given the present state, is verified. Regarding the dependence of the
transition intensities and time, in the case of the failure of the
markovianity, we can use a semi-Markov model in which the future of the
process does not depend on the current time but rather on the duration of the
current state. These two methods are available in MSM.app. In terms of
interpretation, the outputs of the regression, for each transition, are very
similar to those presented for the survival analysis. Results of Figures 16
indicate that save the treatment “Lev(amisole)+5-FU” all the other covariates
have no effect for recurrence transition and only age could be considered
important on mortality without recurrence. At least, only the covariates
“perfor”, “sex” and the treatment “Obs” have no association for the transition
from recurrence to mortality.

Figure 16: Results of the application of the Cox PH model with the following
covariates: “rx”, “sex”, “age”, “obstruct” and “perfor”. Results for each of
the three transition intensities of the colonIDM. A Markovian process is
assumed.

It is also possible to obtain the outputs of ANOVA tests and the p-values of
the tests for nonlinearity. For both outputs, the summaries show the chi-
squared statistics and the p-values. For ANOVA tests, the log-likelihood
values for each parameter are also presented (Figures 17 and 18).

Figure 17: ANOVA results for transitions $0\longrightarrow 1$ and
$0\longrightarrow 2$ with the following covariates: “rx”, “sex”, “age”,
“obstruct” and “perfor” of the colonIDM. A Markovian process is assumed.
Figure 18: The proportional hazards assumption was tested for the transitions
$0\rightarrow 1$ and $0\rightarrow 2$ with the following covariates: “rx”,
“sex”, “age”, “obstruct” and “perfor” of the colonIDM. A Markovian process is
assumed.

Transition probabilities The transition probabilities are quantities of
particularly interest since they allow for long-term predictions of the multi-
state process. The MSM.app allows estimating these quantities using the Aalen-
Johansen estimator [17], the landmark methods (LM) proposed by [18], as well
as its presmoothed version (PLM) proposed by [19], and the landmark Aalen-
Johansen [20]. Related references can also be seen in the papers by [21] and
[22]. Categorical covariates can be included using all four of these methods
by splitting the sample for each level of the covariate and repeating the
described procedures for each subsample. Through the IPCW estimator proposed
by [23] is also possible to estimate transition probabilities conditional on
one single continuous covariate. Finally, the MSM.app also provides the
estimation of transition probabilities conditional on several covariates
through Breslow’s method for estimating the baseline hazard function of the
Cox models fitted marginally to each transition. The outputs for all these
methods are identical. As an example, Figure 19 shows the estimates of the
transition probabilities, using the landmark approach, from the initial single
time $s=365$ days to the next four years (730, 1095, 1460, and 1825 days).
Results are presented combining the values of the corresponding times and
transitions which are labeled from “00” to “12”. For instance, “01”
corresponds to the transition from the initial state (State 0) to intermediate
state (Recurrence, State 1) and, in similar way, “12” the transition from the
intermediate state (Recurrence, State 1) to death (State 2). Plots for each
five transition probabilities are shown in Figure 20.
$\widehat{p}_{00}(s=365,t)$ corresponds to the probability of a individual to
occupy the initial state at time $t$ conditional to be in the same state at
time $s=365$. In a similar way, $\widehat{p}_{11}(s,t)$ represents the
conditional probability of those individuals observed in State 1 at time
$s=365$ to remain in the intermediate state at a later time $t$. Plots for
these transition probabilities report, respectively, survival fractions along
times among the individuals that belong to initial state and the intermediate
state at time $s=365$, being represented by monotone non-increasing functions.
Plots for $\widehat{p}_{02}(s,t)$, report one minus the survival fraction
along time, among the individuals in the initial state at time $s$. In this
case, the plot is given by a monotone non-decreasing function. Finally, plots
for $\widehat{p}_{01}(s,t)$ allows for an inspection along time of the
probability of being in State 1 for the individuals who belong to State 0 at
time $s$. As expected the confidence bands become wider with greater lags
times $t-s$. Finally, Figure 21 depict the estimates of the transition
probabilities for individuals with covariate values (“rx”,“sex”, and “age”) as
Obs,1,48, respectively.

Figure 19: Results of the estimates of the transition probabilities and the
corresponding confidence intervals for s = 365 and times 730, 1090, 1460, and
1825 days for the colonIDM data set using the landmark estimators. Figure 20:
Transition probability estimates with confidence intervals for each transition
for s = 365 using the landmark estimators. Figure 21: Results of the estimates
of the transition probabilities and the corresponding confidence intervals for
s = 365 and times 730, 1090, 1460, and 1825 days for the colonIDM data set
using the Breslow estimator.

Cumulative Incidence Function (CIF) The cumulative incidence of the illness
(intermediate state) is another quantity of interest in IDM models [24]. This
quantity denotes the probability of the individual or item being or having
been in the intermediate ‘diseased’ state at some particular time $t$. It can
be estimated conditional on a covariate, continuous or categorical. Figure 22
shows the estimates, and respective bounds of confidence, of the CIF
conditional to “age” at 50 years for three specific times.

Figure 22: Cumulative recurrence incidence with 95% bootstrap confidence
intervals. Data from a colon cancer study.

10\. How to carry out a multi-state analysis for other models?

By clicking on the “MSM-analysis”, we can extend some of the methods addressed
for the IDM models to more complex multi-state models (MSM) with more than
three states and possible reversible transitions. Number of events Figure 23
shows the movement of the individuals among the six states of the multi-state
model represented by the data set of the European Group for Blood and Marrow
Transplantation (ebmt4 data set).

Figure 23: The number of transitions among the states of the ebmt4 data set.

Regression models Figure 24 shows the results of the Cox regression model for
the transition $1\rightarrow 2$ which includes the covariates “year”, “age”
and “proph”. As we can see, in terms of interpretation, the output is quite
similar to those of IDM models but also shows the global significance of the
model through different tests.

Figure 24: The output of the Cox regression model for the transition
$1\rightarrow 2$ that include the ‘year”, “age” and “proph” covariates. ebmt4
data set.

Transition probabilities The steps for obtaining the transition probabilities
are identical to those used for the illness-death model. In Figure 25, the
output shows all the estimates from the indicated start and last states of the
transition probabilities as well as the corresponding confidence intervals. A
plot with the transition probability can also be presented for each
transition.

Figure 25: Estimates of all possible transition probabilities from the state 1
to 5, for s = 365 and times equal to 730, 1095, 1460, and 1825 using the AJ
estimators. ebmt4 data set.

11\. How to check the Markov condition?

Traditionally the Markov assumption is checked by including covariates
depending on the history through a proportional hazards model. Since the
landmark methods of the transition probabilities are free of the Markov
assumption, they can also be used to introduce such tests by measuring their
discrepancy to Markovian (AJ) estimators.

The MSM.app web application offers two types of tests for checking this
assumption using recent literature methods: (i) local tests, which are
obtained by fixing a specific time value, $s$ and are particularly useful for
estimating transition probabilities; and (ii) global tests, which may be
preferable for regression purposes.

Local tests

Two types of methods are available for checking the local tests of the Markov
condition: (i) the AUC method, which is based on measuring the discrepancy
between the AJ estimator of the transition probabilities (which provides
consistent estimates when the process is Markovian) and the landmark
estimators (which are free of the Markov condition). In this case, the web
tool uses the LM estimator for the progressive illness-death models and LMAJ
in the case of more complex MSM models; (ii) the Log-rank method, which
considers summaries from families of log-rank statistics where patients are
grouped by the state occupied at different times. For both types of models
(IDM and MSM), the output using the AUC is the same. As an example, in Figure
26, we obtained the p-values for the local tests based on the $s$ times 365,
730, 1095, 1460, and 1825 for the transition from state 2 to state 3. Results
were obtained for an IDM model based on the colon cancer data using 100
replicas. Even though the web form is quite similar to the AUC local test, the
output of the log-rank test only provides the results for a specific
transition and the times chosen (Figure 27).

Figure 26: Results of the local test for the illness-death model using the
colon cancer data set, for s = 365, 730, 1095, 1460, and 1825 days, using the
AUC test, from state 1 to state 3. Figure 27: Results of the local test for
the illness-death model using the colon cancer data set, for s = 365, 730 and
1460 days, using the Log-rank test, for state 2 to state 3.

Global tests

In the MSM.app application, three global tests are available for both IMD and
more complex MSM models: (i) the first one is based on Cox models, from which
it is possible to evaluate the effect of history on the process. In this case,
this can be done by checking the significance of the covariate time until
entering the first state of a particular transition. As illustrated in Figure
28, we can conclude that there is no effect of the time spent in the initial
state on the transition $2\rightarrow 3$ (p-value = 0.1543195), which does not
induce the failure of the Markovianity. (ii) the recent proposed global test
propose [25], based on the area under the curves (AUC), can be used. This test
is based on the (AUC) local test results for specific percentiles. The outputs
on the right hand of Figure 29 shows the proportion of rejections of the test
for all possible transitions between state 1 and state 5. (iii) it is also
possible to use the global test based on the log-rank statistics [26]
throughout the similar steps of the previous methods, after selecting Log-rank
in the radio button HTML element. The outputs provide only the results of the
tests for each transition (Figure 30).

Figure 28: Results for this global test given by the Cox PH model to our data
indicated that the effect of the time spent in State 1 is not significant
(p-value of 0.154), revealing no evidence against the Markov model for the
colon data set. Figure 29: Outputs of the global test for the illness-death
model based on the ebmt data set using the AUC test, from state 1 to state 5.
Results for the AUC local test are also shown. Figure 30: Results of the
global test for the illness-death model using the colon cancer data set using
the Log-rank test.

Acknowledgements

This research was financed within the research grants PTDC/MAT-STA /28248/2017
and PD/BD/142887/2018.

## References

  * [1] Wojciechowski, J., Hopkins, A. and Upton, R. Interactive pharmacometric applications using R and the Shiny package. CPT: pharmacometrics & systems pharmacology, 4(3), 146–159; 2015
  * [2] Chang, W. Web Application Framework for R, CRAN; 2017
  * [3] Kaushik, S. Creating Interactive data visualization using Shiny App in R (with examples), Analytics Vidhya; 2016
  * [4] Putter, H., Fiocco, M. and Geskus, R.B. Tutorial in biostatistics: Competing risks and multi-state models, Statistics in Medicine, 26(11), 2389–2430.
  * [5] Meira-Machado, L., de Uña-Álvarez, J. and Cadarso-Suárez, C. and Andersen, P.K. Multi-state models for the analysis of time to event data. Statistical Methods in Medical Research, 18, 195–222; 2009.
  * [6] Meira-Machado, L. and Sestelo, M. Estimation in the progressive illness-death model: A nonexhaustive review. Biometrical Journal, 61 (2), 245–263; 2019.
  * [7] Therneau, T.M. A Package for Survival Analysis in R. https://CRAN.R-project.org/ package=survival; 2021.
  * [8] Putter, H., de Wreede, L.C., Fiocco, M., Geskus, R.B., Bonneville, E.F., Manevski, D. mstate: Data Preparation, Estimation and Prediction in Multi-State Models. The R Journal; 2020.
  * [9] Soutinho, G., Sestelo, M. and Meira-Machado, L. survidm: An R package for Inference and Prediction in an Illness-Death Model. The R Journal, Accept for publication; 2021.
  * [10] Govan, P. eAnalytics: Dynamic Web-based Analytics for the Energy Industry. Journal of Open Research Software, 4: e45; 2016. DOI: http://dx.doi.org/10.5334/jors.144
  * [11] Seal, A., Wild, D.J. Netpredictor: R and Shiny package to perform Drug-Target Bipartite network analysis and prediction of missing links. Cold Spring Harbor Laboratory; 2016. doi: https://doi.org/10.1101/080036
  * [12] Lacy, S. msm-shiny, https://stulacy.shinyapps.io/msm-shiny/; 2021
  * [13] Skourlis, N., Crowther, M.J., Andersson , T., Lambert, P.C. MSMplus, https://nskbiostatistics.shinyapps.io/MSMplus/; 2021
  * [14] Miller, R.G. Survival Analysis. John Wiley & Sons; 1997
  * [15] Moertel, G., Fleming, T.R., Macdonald, J.S., Haller, D.G., Laurie, J.A., Goodman, P.J., Ungerleider, J.S., Emerson, W.A., Tormey, D.C., Glick, J.H., Veeder, M.H. and Mailliard, J.A., Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma. New England Journal of Medicine, 322(6):352–358; 1990
  * [16] Cox, D.R. Regression models and life tables. Journal of the Royal Statistical Society Series B, 34, 187–220; 1972.
  * [17] Aalen, O., Johansen, S. An empirical transition matrix for non homogeneous Markov and chains based on censored observations. Scandinavian Journal of Statistics, 5: 141–150; 1978.
  * [18] de Uña-Álvarez, J. and Meira-Machado, L. Nonparametric estimation of transition probabilities in the non-markov illness-death model: A comparative study. Biometrics, 71(2), 364–375; 2015.
  * [19] Meira-Machado, L. Smoothed landmark estimators of the transition probabilities. SORT-Statistics and Operations Research Transactions, 40: 375–398; 2016.
  * [20] Putter, H. and Spitoni, C. Non-parametric estimation of transition probabilities in non-markov multistate models: The landmark aalen-johansen estimator. Statistical Methods in Medical Research, 27, 2081–2092; 2018.
  * [21] Moreira, A., de Uña-Álvarez, J. and Meira-Machado, L. Presmoothing the aalen-johansen estimator inthe illness-death model. Electronical Journal of Statistics; 2013, 7:1491–1516. DOI: https://doi.org/10.1214/13-EJS816.
  * [22] Araújo A, Meira-Machado L, Roca-Pardiñas J (2014). “TPmsm: Estimation of the Transition Probabilities in 3-State Models.” Journal of Statistical Software, 62(4), 1–29. doi: 10.18637/jss.v062.i04.
  * [23] Meira-Machado, L., de Uña-Álvarez, J. and Datta, S. Nonparametric estimation of conditional transition probabilities in a non-Markov illness-death model. Computational Statistics, 2015, 30(2), 377–397. doi: https://doi.org/10.1007/s00180-014-0538-6
  * [24] Kalbfleisch, J. D. and Prentice R. L. The statistical analysis of failure time data. John Wiley & Sons; 1980.
  * [25] Soutinho, G. and Meira-Machado, L. Methods for checking the Markov condition in multi-state survival data. Comput Stat; 2021. https://doi.org/10.1007/s00180-021-01139-7
  * [26] Titman, A. and Putter, H. General tests of the Markov property in multi-state models. Biostatistics; 2020. doi:10.1093/biostatistics/kxaa030.