Non-linear biphasic mixture model: Existence and uniqueness results

Abstract

This paper is concerned with the development and analysis of a mathematical model that is motivated by interstitial hydrodynamics and tissue deformation mechanics (poro-elasto-hydrodynamics) within an in-vitro solid tumour. The classical mixture theory is adopted for mass and momentum balance equations for a two-phase system. A main contribution of this study is we treat the physiological transport parameter (i.e., hydraulic resistivity) as anisotropic and heterogeneous, thus the governing system is strongly coupled and non-linear. We derived a weak formulation and then formulated the equivalent fixed-point problem. This enabled us to use the Galerkin method, and the classical results on monotone operators combined with the well-known Schauder and Banach fixed-point theorems to prove the existence and uniqueness of results.

1. Introduction

The study of fluid flow through porous media has gained the attention of many researchers over the years. Examples of natural porous materials are living tissues, rocks, soils, etc. On the other hand, there are manmade porous materials, e.g., concrete, foam, ceramic, etc. Because of their wide presence and hard-to-estimate effective properties, flow through porous material is studied by engineers and scientists. This leads to a coupled phenomenon where fluid flow and solid deformation in porous materials interplay. This is a classical problem in geomechanics and biomechanics. Recently, one of the most studied topics in the field of fluid mechanics is flow through biological tissues such as tumours, gylcocalyx layers and articular cartilage. This paper deals with the mathematical modelling and analysis of the coupled phenomena of fluid flow and solid deformation (in short poro-elasto-hydrodynamics) within an in-vitro tumour model. Typically, a tumour is assumed to be a deformable porous medium that consists of multiple phases, e.g., one fluid phase and many solid ones. A tumour may exist in isolation (in-vitro) or may be present in normal tissue (in-vivo) [1]. Although the internal geometrical structure of tumours is complicated, developing mathematical models for approximate situations is very useful. Theoretical predictions generated from such approximate models may help to reduce the number of animal experiments that need to be carried out and also suggest new experimental programmes that identify optimal tumour therapy schedules [2].

Mathematical models of tumours in general can be divided into three categories: discrete, continuum and hybrid. Here, we focus on continuum models that treat cells as averaged populations and are based on the continuum mechanics approach to porous media combined with mixture theory [3]. The early mathematical models on the avascular tumour growth [4] assumed that tumours are made of single type of cells having a constant density. However, various experimental, and theoretical evidence have shown that such a description is not sufficient to study the tumour dynamics [4, 5]. Hence, multiphase models came into play. In this case, one can consider density variations within mixture components to evaluate the evolution of partial stresses. Biot’s theory of poroelasticity and the theory of mixture are commonly adapted models to explore poro-elasto-hydrodynamics [1]. Alike continuum-level approaches involve the development of a set of equations to represent the mechanical behaviour of a soft tissue such as a tumour (which is assumed as a deformable porous material) at the macro scale, using a porous media approach.

The first multiphase model for tumours is proposed by Please et al. [6]. The authors proposed a diffusion equation for cell concentration and generalised Darcy equation for cells and water motion inside the tumour. Further, such multiphase models have been studied analytically and numerically by several authors, one can refer [5, 7, 8]. Sumets et al. [9] described a new boundary-integral representation for biphasic mixture theory, where they solved elastohydrodynamic-mobility problems using boundary element methods. Dey and Raja Sekhar [1] used a biphasic mixture model for poro-elasto-hydrodynamics and nutrient transport inside an in-vitro solid tumour. The authors assumed the presence of an unknown fluid source/sink in the model and biphasic mixture equations have been solved explicitly in the case of one-dimensional spherical geometry. Slvia and Wheeler [10] presented a coupled geomechanics and unsteady reservoir flow model using the theory of poroelasticity. They established the existence and uniqueness of a weak solution and computed a priori error estimates for the numerical solution with stress-dependent permeability. In [11], a non-linear model for a poroelastic medium (described by quasi-static Biot-equations) coupled to transport equations of substances was considered. They have modelled time and space-dependent processes in deformable cellular tissues by the method of homogenisation, starting from a reactive flow system coupling mechanics at the pore scale. The model was analysed and the global-in-time existence and uniqueness of the solution were shown. Cao et al. [12] have considered a non-linear steady flow model in a deformable biological domain based on the theory of poroelasticity (non-linearity is due to the assumption of dilation-dependent interstitial permeability of the solid matrix). They established the existence and uniqueness of a weak solution. Looking through mentioned literature, we observe that there is a gap in dealing with the existence and uniqueness of corresponding biphasic mixture equations that describe the coupled phenomena of fluid flow and solid deformation within biomaterials such as tumours. Attempting to fill such a gap, Alam et al. [13, 14] developed a well-posedness theory and certain regularity results in $2d$ and $3d$ for poro-elasto-hydrodynamics model inside an in-vitro solid tumour. Further, in the case of an in-vivo solid tumour, Alam et al. [15] developed existence and uniqueness resulting in a weak sense for poro-elasto-hydrodynamics while assuming the hydraulic resistivity is heterogeneous and deformation dependent.

We note that, in general, poro-elasto-hydrodynamics models within a tumour may not lead directly to linear biphasic mixture equations. In practice, due to the non-uniform blood vessel distribution, the supply of fluids and macromolecules within a tumour is heterogeneous. As a consequence, physiological transport parameters (e.g., hydraulic resistivity or permeability) depend on space and deformation. For instance, in the case of soft permeable tissue and gel, Barry and Aldis [16] and Holmes and Mow [17] considered permeability depending exponentially on the strain. Also, some of the biological tissues and cells display anisotropic permeability [18]. In particular, articular cartilage typically exhibits anisotropic behaviour [19]. Further, in the case of multicellular tumour tissue, Giverso and Preziosi [20] followed Holmes and Mow [17] and considered that the permeability depends on the volumetric deformation (or strain). Dey and Sekhar [1] assumed that the hydraulic conductivity of soft tumour depends on the radial distance. These non-linear effects inserted in the physiological parameter yield non-linear biphasic mixture equations. As far as we know, for non-linear biphasic mixture models, there is a lack of literature regarding the existence, uniqueness and regularity of the solution. In this paper, we present a non-linear biphasic mixture model that represents poro-elasto-hydrodynamics, which do not account for any new growth of tumour cells. The physiological transport parameter (hydraulic resistivity) is assumed to be deformation dependent, which yields the non-linearity in the model. We develop a local weak solvability theory.

1.1 Biphasic mixture theory

In this subsection, we introduce the generic governing equations. We use biphasic mixture theory to represent the fluid and solid phases of the tumour. Following [1, 8, 21], we apply the conservation of mass and momentum to the fluid and solid phases, viewing the fluid as viscous Newtonian and the solid as deformable, and accounting for momentum exchange between the two phases. Let $\mathbf{V}_f$ and $\mathbf{V}_s$ denote the velocities of the fluid and solid phases, respectively. The apparent densities of the fluid and solid phases are denoted by $\tilde{\rho }_f$ and $\tilde{\rho }_s,$ respectively, and their corresponding volume fractions by $\phi _f$ and $\phi _s.$ The true densities of the fluid and solid phases are then $\rho _f=\phi _f\tilde{\rho }_f$ and $\rho _s=\phi _s\tilde{\rho }_s.$ Accordingly, in $\Omega$ , the mass and linear momentum balance equations for the fluid phase are given by

(1.1) \begin{equation} \frac{\partial (\tilde{\rho }_f\phi _f)}{\partial t} +\nabla \cdot [(\tilde{\rho }_f\phi _f)\mathbf{V}_f]= \tilde{\rho }_f S_f, \end{equation}

(1.2) \begin{equation} \rho _f\left (\frac{\partial \mathbf{V}_f}{\partial t}+(\mathbf{V}_f\cdot \nabla )\mathbf{V}_f \right )= \nabla \cdot \mathbf{T}_f+\boldsymbol{\Pi }_f+\mathbf{b}_f, \end{equation}

where $\textbf{T}_{f}$ denotes the stress tensor for the fluid phase

(1.3) \begin{eqnarray} \mathbf{T}_{f}=-[\phi _{f}P-\lambda _{f}\nabla \cdot \mathbf{V}_{f}]\mathbf{I} +\mu _{f}[\nabla \mathbf{V}_{f}+(\nabla \mathbf{V}_{f})^{T}]. \end{eqnarray}

The corresponding mass and linear momentum equations for the solid phase are

(1.4) \begin{equation} \frac{\partial (\tilde{\rho }_s\phi _s)}{\partial t} +\nabla \cdot [(\tilde{\rho }_s\phi _s)\mathbf{V}_s]= \tilde{\rho }_f S_s, \end{equation}

(1.5) \begin{equation} \rho _s\left (\frac{\partial \mathbf{V}_s}{\partial t}+(\mathbf{V}_s\cdot \nabla )\mathbf{V}_s \right )= \nabla \cdot \mathbf{T}_s+\boldsymbol{\Pi }_s+\mathbf{b}_s, \end{equation}

where $\mathbf{T}_{s}$ denotes the stress tensor for the solid phase

(1.6) \begin{eqnarray} \mathbf{T}_{s}=-[\phi _{s}P-\chi _s (\nabla \cdot \textbf{U}_{s})]\mathbf{I}+\mu _{s}[\nabla \mathbf{U}_{s}+(\nabla \mathbf{U}_{s})^{T}]. \end{eqnarray}

In Equations (1.1) and (1.4), $S_f$ and $S_s$ are fluid and solid source terms, respectively. $\mathbf{U}_s$ denotes the displacement of the solid phase. Hence, $\mathbf{V}_{s}=\frac{\partial \mathbf{U}_{s}}{\partial t}.$ The average interstitial fluid pressure is $P$ and $\mathbf{b}_j$ $j=\{1,2\}$ denotes the body force. Further, the volume fractions $\phi _f$ and $\phi _s$ are assumed to satisfy the following saturation assumption

(1.7) \begin{equation} \phi _f+\phi _s= 1. \end{equation}

We suppose that the two phases interact together via the drag forces $\boldsymbol{\Pi }_{s}$ and $\boldsymbol{\Pi }_{f},$ , which by Newton’s third law, are equal and opposite. Following [7, 8], we define

(1.8) \begin{eqnarray} -\boldsymbol{\Pi }_{s}=\boldsymbol{\Pi }_{f}={K}(\textbf{V}_{s}-\textbf{V}_{f})- (\nabla \phi _{s})P, \end{eqnarray}

where $K=\frac{\mu _f}{k}$ is the hydraulic resistivity or drag coefficient, where $k$ is the permeability of the porous matrix (the precise nature of $K$ will be defined in the next section). Furthermore, $\mu _f$ ( $\mu _s$ ) is the dynamic viscosity of the fluid phase (solid phase), while $\lambda _f$ and $\chi _s$ denote the Lame coefficient and shear modulus of the fluid and solid phases. The elastic modulli $\chi _s$ and $\mu _s$ are related to the Young’s modulus $(\mathcal{Y})$ and Poisson’s ratio $(\nu _p)$ via the relationship

\begin{equation*}\chi _s=\frac {\nu _p\mathcal {Y}}{(1+\nu _p)(1-2\nu _p)}\,\mbox {and}\,\mu _s =\frac {\mathcal {Y}}{2(1+\nu _p)}.\end{equation*}

1.2 Main modelling assumptions

Having presented a more generic set of mixture theory equations in the absence of any assumptions or boundary conditions, here we list some biologically suitable assumptions restricted to a specific tumour model. Our choice of model is motivated by the study of an isolated (in-vitro) tumour that behaves as a heterogeneous deformable porous medium. Usually, tumour tissues are considered incompressible fluids with no voids present. It is assumed that each phase has an equal constant density. Suppose $\Omega \subset \mathbb{R}^d,$ $d=\{2,3\}$ is a bounded Lipschitz domain that is filled by the tumour. Let $\partial \Omega$ be its boundary (Figure 1). One may note that the solid tumour is essentially a multicellular spheroid. When nutrients perfuse the interstitial space, a large number of cells receive adequate food for survival and proliferation. As a consequence, the tumour grows in size. For a growing tumour, the permeability and the effective mechanics parameters ( $\chi _s$ and $\mu _s$ ) may depend on the volume fraction of the cell population [8]. Moreover, the volume fractions depend on space. Hence, it is extremely difficult to analyse mathematically the growth model and fluid transport model simultaneously. To simplify this, we assume that the tumour tissue is not growing and all elastic parameters ( $\chi _s$ and $\mu _s$ etc.) are independent of volume fractions. Further, we make the following modelling assumptions (A1)–(A4):

1. (A1) Nutrient perfusion and transport occur on much shorter timescales than the timescale for tumour cell growth. Accordingly, we view the tumour as a static, perfused biological domain. On the short timescale associated with nutrient transport to (and within) the tumour, cell death and proliferation are assumed to be negligible. Therefore, we fix $S_s = 0$ in the tumour and normal tissue regions. Further, on the timescale of interest, the solid volume fraction $\phi _{s}$ remains constant, and for simplicity, we assume it to be independent of spatial position and time as (see [1]).

2. (A2) A fluid source is attached to the mixture, hence $S_f\neq 0,$ [1, 22].

3. (A3) Motion of the cells and interstitial fluid flow is so slow that the inertial terms can be neglected in both phases, see e.g. [8, 23].

4. (A4) Structure of the hydraulic resistivity: Various experimental and theoretical investigations indicate clearly that for the deformable porous medium (or soft biological tissue such as articular cartilage, arterial tissue and tumour), the permeability also called in this context hydraulic resistivity depends on stress, dilatation, volume fractions (porosity), etc. There are several analytical expressions for permeability that are used in literature such as $k(x)=\exp\![mU'(x)],$ where $m$ is a constant and $U'=\frac{dU}{dx}$ with $U$ as displacement, while $k(x)=k_0/[1 - mU'(x)]$ or $k(x)=k_0[1 + mU'(x)]$ for small $mU'(x)$ , or $k(x)=k_0(\phi/\phi _0)^n$ [16, 17]. Here, $k_0$ is the permeability at reference porosity $\phi _0$ and $n$ is a variable that may be determined by fitting experimental data. In this framework, we propose two different cases. We assume that the hydraulic resistivity $K$ depends explicitly on (a) the solid phase displacement $\mathbf{U}_s$ and (b) the strain/dilatation i.e. on $\nabla \cdot \mathbf{U}_s.$ In both of these cases, we admit anisotropic effects, i.e. $\mathbf{K}$ is a square matrix of order $d =\{2,3\}.$ Note that in case (a), we do not use any specific, explicit expression of $\mathbf{K}(\mathbf{U}_s)$ in our analysis. One may think through the following choice

   \begin{equation*}\mathbf {K}(\mathbf {U}_s)=(\alpha _1|\mathbf {U}_s|+\alpha _3)\mathbf {I} +(\alpha _1-\alpha _2) \frac {\mathbf {U}_s\mathbf {U}_s^T}{|\mathbf {U}_s|},\end{equation*}
   where $\alpha _i$ ( $i=1,\,2,\,3$ ) are real constants such that $\alpha _1\geqslant 0,\alpha _2\geqslant 0$ and $\alpha _3\gt 0.$ One can easily show that $\mathbf{K}$ is Lipschitz and uniformly positive [24]. For case (b), one can think of a form $\mathbf{K}(\nabla \cdot \mathbf{U}_s) = (\gamma _1+\gamma _2|\nabla \cdot \mathbf{U}_s|)\mathbf{I},$ where $\gamma _i\geqslant 0$ are real constants, and $\mathbf{I}$ is the identity matrix. One can observe that $\mathbf{K}$ is Lipschitz and uniformly positive. We will state further assumptions on $\textbf{K}$ in the next Section 3, [16, 17].

The assumptions (A1)–(A4) are taken into account in Equations (1.1)–(1.8). This helps us to define a coupled non-linear system of steady-state mass and momentum balance equations in unknowns $(\mathbf{V}_f,\mathbf{U}_s, P)$ as follows:

(1.9) \begin{equation} -\nabla \cdot (2\mu _fD(\mathbf{V}_f)+\lambda _f(\nabla \cdot \mathbf{V}_f)\mathbf{I} -\phi _fP\mathbf{I}) +\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f=\mathbf{b}_f\,\,\mbox{in}\,\,\Omega, \end{equation}

(1.10) \begin{equation} -\nabla \cdot (2\mu _sD(\mathbf{U}_s)+\chi _s(\nabla \cdot \mathbf{U}_s)\mathbf{I}-\phi _sP\mathbf{I}) -\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f=\mathbf{b}_s\,\,\mbox{in}\,\,\Omega, \end{equation}

(1.11) \begin{equation} \nabla \cdot (\phi _f\mathbf{V}_f)=S_f\,\,\mbox{in}\,\,\Omega, \end{equation}

where $\boldsymbol{\varsigma }$ is either $\nabla \cdot \mathbf{U}_s$ or $\mathbf{U}_s.$ Further, $D({\cdot})$ denotes the deviatoric matrix, which is defined as $D(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u})^t),$ where $(\nabla \mathbf{u})^t$ denotes the transpose of the matrix $\nabla \mathbf{u}.$ We have made the approximation $\textbf{V}_{s}=\frac{\partial \mathbf{U}_{s}}{\partial t}\approx 0$ which is consistent with the infinitesimal strain theory that is used in e.g. [1, 25]. The mass balance equations for the fluid phase include a source term $S_f$ , which models fluid exchange with vasculature and lymph vessels. For a closed mixture (for instance the case of avascular tumours), we consider $S_f=0$ so that $\nabla \cdot (\phi _f\mathbf{V}_f)=0.$ This is the counterpart of the incompressibility constraint. Note that in Equation (1.11), even though the density of each phase is constant, the vector $\mathbf{V}_f$ is not solenoidal. On the other hand, when the external sources/sinks are attached to the mixture, we have $S_f\neq 0$ [5, 22]. Typically, the fluid source $S_f$ is assumed to be driven by the average transmural pressure and (trusting [1, 22]) takes the form

(1.12) \begin{equation} S_f=-{L_p}\left (\frac{A}{V}\right )\{1+L_rA_r\}(P-P_F), \end{equation}

where $L_p$ is the average hydraulic conductivity coefficient of capillary. In (1.12), $A/{V}$ denotes the capillary surface area per unit tissue volume in the tumour tissue and $L_rA_r$ denotes the ratio of the strength of distributed solute source through the vasculature and solute sink through the lymph vessels and $P_F$ is the weighted vascular pressure.

Figure 1. Geometry of the problem.

1.3 Boundary conditions

The model that we have considered here is supposed to mimic an in-vitro tumour. Accordingly, we prescribe

(1.13) \begin{equation} \mathbf{T}_f\cdot \mathbf{n}=\mathbf{T}_\infty \,\,\mbox{and}\,\mathbf{U}_s=0\,\mbox{on}\,\,\partial \Omega, \end{equation}

where $\mathbf{n}$ is the outward normal unit vector to the boundary $\partial \Omega .$

2. Non-dimensional equations

Using the transformations $\hat{\textbf{x}}=\frac{\textbf{x}}{R},$ $\nabla '=R\nabla,$ $\hat{P}=\frac{P}{P_F},$ $\hat{\mathbf{V}}_f =\frac{\mathbf{V}_f}{\frac{RP_F}{\mu _f}},$ $\hat{\mathbf{U}}_s=\frac{\mathbf{U}_s}{\frac{R^3P_F}{\mu _f\nu }},$ $\mathbf{K}=\hat{\mathbf{K}}K_d,$ where $K_{d}$ is the hydraulic resistivity (drag coefficient) of the tumour tissue in the absence of deformation, $R$ is the length of the edge of the $d$ -cube in which $\Omega$ is contained. The following dimensionless form of the governing Equations (1.9)–(1.11) (‘hat’ is dropped for convenience) is available in $\Omega$

(2.1) \begin{equation} -\nabla \cdot \left (2D(\mathbf{V}_f)+\frac{\lambda _f}{\mu _f}(\nabla \cdot \mathbf{V}_f)\mathbf{I} -\phi _fP\mathbf{I}\right ) +\frac{1}{Da}\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f=\mathbf{b}_f, \end{equation}

(2.2) \begin{equation} -\nabla \cdot \left (\frac{\varrho _t}{(1+\nu _p)}D(\mathbf{U}_s)+\frac{\nu _p\varrho _t}{(1+\nu _p)(1-2\nu _p)} (\nabla \cdot \mathbf{U}_s)\mathbf{I}-\phi _sP\mathbf{I}\right ) -\frac{1}{Da}\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f=\mathbf{b}_s, \end{equation}

(2.3) \begin{equation} \nabla \cdot (\phi _f\mathbf{V}_f)=-\alpha ^2_t(1+L_rA_r)(P-1). \end{equation}

In (2.1) and (2.2), $\mathbf{b}_f$ and $\mathbf{b}_s$ are modified non-dimensional body forces, and $\varrho _{t}={\mathcal{Y}R^{2}\rho _{f}}/{\mu ^{2}_{f}}$ is the dimensionless Young’s modulus $Y$ associated with the solid phase. It contains the response of the solid phase (cellular phase + extracellular matrix) towards viscous drag due to interstitial fluid movement. $\alpha ^2_{t}=L_P(A/V)\mu _f$ is the strength of solute source, and $Da=\frac{K_d\mu _f}{R^2}$ is the Darcy number (permeability parameter). The corresponding boundary conditions are

(2.4) \begin{equation} \left (2D(\mathbf{V}_f)+\frac{\lambda _f}{\mu _f} (\nabla \cdot \mathbf{V}_f)\mathbf{I}-\phi _fP\mathbf{I} \right )\cdot \mathbf{n}=\mathbf{T}_\infty \,\mbox{and}\,\mathbf{U}_s =0\,\,\mbox{on}\,\,\partial \Omega . \end{equation}

For the sake of writing convenience, set $\lambda =\frac{\lambda _f}{\mu _f},$ $\alpha _1=\frac{\varrho _t}{2(1+\nu _p)},$ $\alpha _2=\frac{\nu _p\varrho _t}{(1+\nu _p)(1-2\nu _p)},$ and $a_0=\alpha ^2_t(1+L_rA_r).$ Observe that the system of Equations (2.1)–(2.3) is non-linear and fully coupled whenever $\boldsymbol{\varsigma }$ is equals to either $\mathbf{U}_s$ or $\nabla \cdot \mathbf{U}_s$ , which is our primary interest for now. The main aim is to study the well-posedness of the non-linear system (2.1)–(2.3) subject to the data (2.4).

3. Well-posedness of the auxiliary sub-problems

$\mathbf{V}_f,$ $P,$ $\mathbf{U}_s$ are the unknown functions in the system of Equations (2.1)–(2.3). We assume the following:

1. (A) The parameters $\phi _f\gt 0,$ $\phi _s\gt 0,$ $\lambda \geqslant 0,$ $\alpha _1\gt 0,$ $\alpha _2\gt 0,$ $a_0\gt 0,$ $Da\gt 0$ are known real constants, and the functions $\mathbf{b}_j\in L^2(\Omega )^d$ where $j=f,\,s,$ $\mathbf{T}_\infty \in L^{2}(\partial \Omega )^d$ are also known. $c_k\gt 0,\,c_p\,{\gt}\,0, c_t\gt 0,\,c_s\gt 0$ are some real constants that appear in Korn’s, Poincare’s, trace and Sobolev’s inequalities, respectively.1

2. (B) Let $\mathbf{K}\,{:}\,\mathbb{R}^d\rightarrow \mathbb{R}^{d^2}$ be a symmetric, uniformly bounded and positive definite matrix. This ensures that there exist positive constants $k_{1}$ and $k_2$ such that for all $\boldsymbol{\xi },\,\mathbf{x}\in \mathbb{R}^d$ , we have:

   (3.1) \begin{equation} \mbox{(i)}\,k_{1}\boldsymbol{\xi }\cdot \boldsymbol{\xi }\leqslant \mathbf{K}(\mathbf{x})\boldsymbol{\xi }\cdot \boldsymbol{\xi }\,\, \mbox{and}\,\,\mbox{(ii)}\,\,||\mathbf{K}(\mathbf{x})||\leqslant k_2, \end{equation}
   $||\cdot ||$ denotes the Euclidean norm.

3. (C) We assume that the hydraulic resistivity $\mathbf{K}$ is Lipschitz continuous. To this extent, let us assume that there exists a constant $k_L\gt 0$ such that

   (3.2) \begin{equation} ||\mathbf{K}(\mathbf{x})-\mathbf{K}(\mathbf{y})||\leqslant k_L||\mathbf{x}-\mathbf{y}||\,\mbox{for all}\,\mathbf{x},\,\mathbf{y}\in \mathbb{R}^d. \end{equation}

3.1 Concept of weak formulation

Choose the triplet of test functions $(\mathbf{W},\mathbf{Z},q)\in H^1(\Omega )^d \times H^1_{0}(\Omega )^d \times L^2(\Omega )$ . Taking the scalar product of (2.1) with $\mathbf{W}$ , (2.2) with $\mathbf{Z}$ and (2.3) with $q,$ and using the boundary conditions (2.4), we get the following non-linear weak formulation:

Find the triplet $(\mathbf{V}_f, \mathbf{U}_s,P)\in H^1(\Omega )^d\times H^1_{0}(\Omega )^d\times L^2(\Omega )$ such that

(3.3) \begin{align} 2(D(\mathbf{V}_f)\,{:}\,D(\mathbf{W}))_{\Omega }+\lambda (\nabla \cdot \mathbf{V}_f,\nabla \cdot \mathbf{W})_{\Omega } -\phi _f(P,\nabla \cdot \mathbf{W})_{\Omega }\\ +\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f,\mathbf{W})_{\Omega } =(\mathbf{b}_f,\mathbf{W})_{\Omega } +(\mathbf{T}_\infty,\mathbf{W})_{\partial \Omega }\nonumber \end{align}

(3.4) \begin{align} \phi _f(\nabla \cdot \mathbf{V}_f,q)_{\Omega }+a_0(P,q)_{\Omega }=(a_0,q)_{\Omega } \end{align}

(3.5) \begin{align} 2\alpha _1(D(\mathbf{U}_s)\,{:}\,D(\mathbf{Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf{U}_s,\nabla \cdot \mathbf{Z})_{\Omega } -\phi _s(P,\nabla \cdot \mathbf{Z})_{\Omega }\nonumber \\ -\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f,\mathbf{Z})_{\Omega }= (\mathbf{b}_s,\mathbf{Z})_{\Omega } \end{align}

holds for all $(\mathbf{W},\mathbf{Z},q)\in H^1(\Omega )^d\times H^1_{0}(\Omega )^d\times L^2(\Omega ).$ Here in $\boldsymbol{\varsigma }$ is $\mathbf{U}_s$ and $\nabla \cdot \mathbf{U}_s.$

Lemma 1. (Equivalence of weak formulations) Suppose that parameters and data satisfy assumptions (A) and (B). Then any solution (in the sense of distributions) $(\mathbf{V}_f, \mathbf{U}_s,P)\in H^1(\Omega )^d\times H^1_{0}(\Omega )^d\times L^2(\Omega )$ of the coupled problem (2.1)–(2.4) is also a solution to the variational problem (3.3)–(3.5). Conversely, any solution to the weak problem (3.3)–(3.5) satisfies (2.1)–(2.4) in the sense of distributions.

The proof follows using standard arguments. We omit to show it.

3.2 Decoupled problem corresponding to (3.3)–(3.5), Case (a): $\mathbf{K}\mathrm{(}\boldsymbol{\varsigma }\mathrm{)}=\mathbf{K}\mathrm{(}\mathbf{U}_s\mathrm{)}$

We note that the weak formulation (3.3)–(3.5) can be decoupled concerning the unknowns $(\mathbf{V}_f, P)$ and $\mathbf{U}_s.$ In case (a), we are dealing with $\mathbf{K}(\boldsymbol{\varsigma })=\mathbf{K}(\mathbf{U}_s)$ , which is a non-linear function of $\mathbf{U}_s$ and satisfies assumptions (3.1) and (3.2). In this case, we can solve (3.3)–(3.5) sequentially, that is, given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ find $(\mathbf{V}_f,P)\in H^1{(\Omega )}^d \times L^2(\Omega )$ such that

\begin{align*} (Q_{w_1}(\boldsymbol{\varsigma })){\left \{\begin{array}{ll} 2(D(\mathbf{V}_f)\,{:}\,D(\mathbf{W}))_{\Omega }+\lambda (\nabla \cdot \mathbf{V}_f,\nabla \cdot \mathbf{W})_{\Omega } -\phi _f(P,\nabla \cdot \mathbf{W})_{\Omega }\\[6pt] +\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f,\mathbf{W})_{\Omega } +\phi _f(\nabla \cdot \mathbf{V}_f,q)_{\Omega } +a_0(P,q)_{\Omega }\\[6pt] =(\mathbf{b}_f,\mathbf{W})_{\Omega } +(\mathbf{T}_\infty,\mathbf{W})_{\partial \Omega }+(a_0,q)_{\Omega } \end{array} \right .} \end{align*}

holds for all $(\mathbf{W},q)\in H^1(\Omega )^d \times L^2(\Omega ),$ and then for a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega ),$ find $\mathbf{U}_s\in H^1_0(\Omega )^d$ such that

\begin{align*} (Q_{w_2}(\mathbf{U}_s))\left \{ \begin{array}{ll} 2\alpha _1(D(\mathbf{U}_s)\,{:}\,D(\mathbf{Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf{U}_s,\nabla \cdot \mathbf{Z})_{\Omega } = \phi _s(P,\nabla \cdot \mathbf{Z})_{\Omega }\\[6pt] +\frac{1}{Da}(\mathbf{K}(\mathbf{U}_s)\mathbf{V}_f,\mathbf{Z})_{\Omega } +(\mathbf{b}_s,\mathbf{Z})_{\Omega } \end{array} \right . \end{align*}

holds for all $\mathbf{Z}\in H^1_0(\Omega )^d.$ For notational convenience, we denote combined problem as $(Q_{w_1}(\boldsymbol{\varsigma }))-(Q_{w_2}(\mathbf{U}_s)).$

3.3 Decoupled problem corresponding to (3.3)–(3.5), Case (b): $\mathbf{K}$ ( $\nabla \cdot \boldsymbol{\varsigma }$ ) = $\mathbf{K}$ ( $\nabla \cdot \mathbf{U}_s$ )

We can solve (3.3)–(3.5) sequentially, that is, for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ find $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega )$ such that

\begin{equation*}(Q_{w_1}(\nabla \cdot \boldsymbol {\varsigma })){\left \{ \begin {array}{ll} 2(D(\mathbf {V}_f)\,{:}\,D(\mathbf {W}))_{\Omega }+\lambda (\nabla \cdot \mathbf {V}_f,\nabla \cdot \mathbf {W})_{\Omega } -\phi _f(P,\nabla \cdot \mathbf {W})_{\Omega } \\[6pt] +\frac {1}{Da}(\mathbf {K}(\nabla \cdot \boldsymbol {\varsigma })\mathbf {V}_f,\mathbf {W})_{\Omega } +\phi _f(\nabla \cdot \mathbf {V}_f,q)_{\Omega } +a_0(P,q)_{\Omega }\\[6pt] ={(\mathbf {b}_f,\mathbf {W})}_{\Omega } +{(\mathbf {T}_\infty,\mathbf {W})}_{\partial \Omega }+(a_0,q)_{\Omega } \end {array} \right .}\end{equation*}

holds for all $(\mathbf{W},q)\in H^1(\Omega )^d \times L^2(\Omega ),$ and then for a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega ),$ find $\mathbf{U}_s\in H^1_0(\Omega )^d$ such that

\begin{equation*}(Q_{w_2}(\nabla \cdot \mathbf {U}_s))\left \{ \begin {array}{ll} 2\alpha _1(D(\mathbf {U}_s)\,{:}\,D(\mathbf {Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf {U}_s,\nabla \cdot \mathbf {Z})_{\Omega } = \phi _s(P,\nabla \cdot \mathbf {Z})_{\Omega }\\[6pt] +\frac {1}{Da}(\mathbf {K}(\nabla \cdot \mathbf {U}_s)\mathbf {V}_f,\mathbf {Z})_{\Omega } +(\mathbf {b}_s,\mathbf {Z})_{\Omega } \end {array} \right .\end{equation*}

holds for all $\mathbf{Z}\in H^1_0(\Omega )^d.$ For notational convenience, we denote combined problem as $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s)).$

3.4 Case (a): existence and uniqueness results for $(Q_{w_1}(\boldsymbol{\varsigma }))$

In order to solve weak formulation $(Q_{w_1}(\boldsymbol{\varsigma }))$ , we use the following method. We rephrase the weak formulation $(Q_{w_1}(\boldsymbol{\varsigma }))$ into an abstract setting. Set $\mathbb{Y}=H^1(\Omega )^d \times L^2(\Omega ).$ To do so, define a mapping $\mathcal{H}_{\boldsymbol{\varsigma }}$ from $\mathbb{Y}$ to $\mathbb{Y}$ by

(3.6) \begin{align} \nonumber \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P),(\mathbf{W},q)\rangle = 2(D(\mathbf{V}_f)\,{:}\,D(\mathbf{W}))_{\Omega }+\lambda (\nabla \cdot \mathbf{V}_f,\nabla \cdot \mathbf{W})_{\Omega } -\phi _f(P,\nabla \cdot \mathbf{W})_{\Omega }\nonumber \\ +\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f,\mathbf{W})_{\Omega } +\phi _f(\nabla \cdot \mathbf{V}_f,q)_{\Omega }+a_0(P,q)_{\Omega }-[(\mathbf{b}_f,\mathbf{W})_{\Omega } +(\mathbf{T}_\infty,\mathbf{W})_{\partial \Omega } +(a_0,q)_{\Omega }]. \end{align}

Using the mapping $\mathcal{H}_{\boldsymbol{\varsigma }},$ the variational formulation $(Q_{w_1}(\boldsymbol{\varsigma }))$ can equivalently be written as: for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ , find $(\mathbf{V}_f,P)\in \mathbb{Y}$ such that

(3.7) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P), (\mathbf{W},q)\rangle =0\,\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}. \end{eqnarray}

Conversely, if (3.7) holds, then (2.1) and (2.3) with the first boundary condition in the equation (2.4) satisfy in the sense of distributions. Hence, our immediate task is to find a pair $(\mathbf{V}_f, P)\in \mathbb{Y}$ that satisfies (3.7). In order to do so, we proceed as follows. The mapping $\mathcal{H}_{\boldsymbol{\varsigma }}$ satisfies the following lemma:

Lemma 3.1. Suppose that parameters and data satisfy assumptions (A) and (B). If $\mathcal{H}_{\boldsymbol{\varsigma }}$ is a mapping from $\mathbb{Y}$ into itself defined by (3.6), then the following statements hold:

1. (i) $\mathcal{H}_{\boldsymbol{\varsigma }}$ is continuous.

2. (ii) There exists a real number $r\gt 0$ such that

   \begin{equation*}\langle \mathcal {H}_{\boldsymbol {\varsigma }}(\mathbb {V}),\mathbb {V}\rangle \gt 0,\,\,\mbox {for all}\,\,\mathbb {V}\in \mathbb {Y}\,\, \mbox {with}\,\,||\mathbb {V}||_{\mathbb {Y}}=r,\end{equation*}

i.e., $\mathcal{H}_{\boldsymbol{\varsigma }}$ is coercive on a ball of radius $r$ in $\mathbb{Y}.$ Here, for any $\mathbb{V}=(\mathbf{V}_f,P)\in \mathbb{Y}= H^1(\Omega )^d \times L^2(\Omega ),$ $||\cdot ||_{\mathbb{Y}}$ is defined as

\begin{equation*}||\mathbb {V}||^2_{\mathbb {Y}}=||(\mathbf {V}_f,P)||^2_{\mathbb {Y}} =||\mathbf {V}_f||^2_{1,\Omega }+||P||^2_{0,\Omega }.\end{equation*}

Proof. (i) The continuity of the mapping $\mathcal{H}_{\boldsymbol{\varsigma }}$ can be shown using the continuity of scalar product. Indeed, let $\{\mathbb{V}^m\}_{m\geqslant 1}= \{(\mathbf{V}^{m}_f,P^{m})\}_{m\geqslant 1}$ be any sequence in $\mathbb{Y}$ that converges strongly to $\mathbb{V}=(\mathbf{V}_f,P)\in \mathbb{Y}$ as $m\rightarrow \infty,$ i.e.

(3.8) \begin{align} ||\mathbf{V}^{m}_f-\mathbf{V}_f||_{1,\Omega }\rightarrow 0,\,\, ||P^{m}-P||_{0,\Omega }\rightarrow 0\,\mbox{as }\,m\rightarrow \infty . \end{align}

Relying on the definition of $\mathcal{H}_{\boldsymbol{\varsigma }}$ and on Cauchy-Schwarz inequality, we get

\begin{align*} |\langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}^m) -\mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}),(\mathbf{W},q)\rangle |\leqslant 2||\mathbb{D}^{\,f}(\mathbf{V}^m_f-\mathbf{V}_f)||_{0,\Omega }||\mathbb{D}^{\,f}(\mathbf{W})||_{0,\Omega } \\ +\lambda ||\nabla \cdot (\mathbf{V}^m_f-\mathbf{V}_f)||_{0,\Omega }||\nabla \cdot \mathbf{W}||_{0,\Omega } +\phi _f||P^m-P||_{0,\Omega }||\nabla \cdot \mathbf{W}||_{0,\Omega }\\ +\frac{k_2}{\mbox{Da}}||\mathbf{V}^m_f-\mathbf{V}_f||_{0,\Omega }||\mathbf{W}||_{0,\Omega } +\phi _f||\nabla \cdot (\mathbf{V}^m_f-\mathbf{V}_f)||_{0,\Omega }||q||_{0,\Omega } \\ +a_0||P^m-P||_{0,\Omega }||q||_{0,\Omega }. \end{align*}

Using (3.8), we obtain

\begin{eqnarray*} |\langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}^{m}_f,P^{m}) -\mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P), (\mathbf{W},q)\rangle |\rightarrow 0\,\, \,\forall \,(\mathbf{W},q)\,\mbox{as}\, m\rightarrow \infty . \end{eqnarray*}

This argument establishes the continuity of $\mathcal{H}_{\boldsymbol{\varsigma }}.$

(ii) For any $\mathbb{V}=(\mathbf{V}_f,P)\in \mathbb{Y},$ we have

(3.9) \begin{align} \nonumber \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}),\mathbb{V}\rangle = 2||D(\mathbf{V}_f)||^2_{0,\Omega }+\lambda ||\nabla \cdot \mathbf{V}_f||^2_{0,\Omega } -\phi _f(P,\nabla \cdot \mathbf{V}_f)_{\Omega }+\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })\mathbf{V}_f,\mathbf{V}_f)_{\Omega } \\ +\phi _f(\nabla \cdot \mathbf{V}_f,P)_{\Omega }+a_0||P||^2_{0,\Omega } -[(\mathbf{b}_f,\mathbf{V}_f)_{\Omega } +(\mathbf{T}_\infty,\mathbf{V}_f)_{\partial \Omega }+(a_0,P)_{\Omega }]. \end{align}

Using Cauchy-Schwarz, Korn’s and trace inequalities, we obtain

(3.10) \begin{align} \nonumber \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}),\mathbb{V}\rangle \geqslant \alpha ||\mathbf{V}_f||^2_{1,\Omega }+\lambda ||\nabla \cdot \mathbf{V}_f||^2_{0,\Omega } +{a_0}||P||^2_{0,\Omega }\\ -(||\mathbf{b}_f||_{0,\Omega }+\sqrt{c_t}||\mathbf{T}_{\infty }||_{0,\partial \Omega })||\mathbf{V}||_{1,\Omega } -||a_0||_{0,\Omega }||P||_{0,\Omega }, \end{align}

where $\alpha =\frac{1}{c_k}\min \{2,\frac{k_1}{Da}\}.$ Further, (3.10) can be rewritten as

(3.11) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}),\mathbb{V}\rangle \geqslant \alpha _3(||\mathbf{V}_f||^2_{1,\Omega } +||P||^2_{0,\Omega })-\alpha _4||\mathbb{V}||_{1,0,\Omega }, \end{eqnarray}

where

(3.12) \begin{align} \alpha _3&=\min \left \{\alpha,{a_0}\right \}, \end{align}

(3.13) \begin{align} \alpha _4&=[(||\mathbf{b}_f||_{0,\Omega }+\sqrt{c_t}||\mathbf{T}_{\infty }||_{0,\partial \Omega })^2 +||a_0||^2_{0,\Omega }]^{1/2}. \end{align}

If $||\mathbb{V}||_{\mathbb{Y}}=r_0$ for some $r_0\gt 0,$ then we have

(3.14) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbb{V}),\mathbb{V}\rangle \gt 0\,\,\,\forall \,\mathbb{V}\in \mathbb{Y}, \,\,\mbox{when}\,\,r_0\gt \frac{\alpha _4}{\alpha _3}. \end{eqnarray}

This completes the proof of Lemma 3.1.

Remark 3.2. Note that $\lambda =\frac{\lambda _f}{\mu _f}$ (which is the ratio of the two viscosity coefficients) plays a significant role in the coercivity proof. The literature [26, 27] suggests that researchers have debated on the sign (or value) of $\lambda$ . According to the well-known Stokes-hypothesis, $\lambda =-2/3$ [27]. On the other hand, the existing literature also suggests that this ratio can be non-negative [26]. Thus, we consider both of these possibilities. If $\lambda \geqslant 0$ , then the coercivity of $\mathcal{H}_{\boldsymbol{\varsigma }},$ as shown by us in (3.14), holds. However, when $\lambda \lt 0$ (i.e. a typical Stokes hypothesis), the coercivity of $\mathcal{H}_{\boldsymbol{\varsigma }}$ holds with relevant restrictions on the constants. For instance, when $\lambda \lt 0,$ the mapping $\mathcal{H}_{\boldsymbol{\varsigma }}$ is coercive if $\alpha =\frac{1}{c_k}\min \{2,\frac{k_1}{Da}\}\gt 2/3.$ For convenience from here onward, we assume $\lambda$ to be a non-negative constant.

Based on Lemma A.1 (see Appendix A), we now present the following existence results.

Theorem 3.3. Suppose that (A) and (B) hold. Then for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ , the problem (3.7) has at least one solution $(\mathbf{V}_f,P)\in \mathbb{Y}=H^1(\Omega )^d\times L^2(\Omega )$ satisfying the problems (2.1) and (2.3) with the first boundary condition in the equation (2.4) in the sense of distributions. Moreover, the solution $(\mathbf{V}_f,P)$ satisfies the following a priori bound

(3.15) \begin{eqnarray} ||\mathbf{V}_f||^2_{1,\Omega }+||P||^2_{0,\Omega }\leqslant \left (\frac{\alpha _4}{\alpha _3}\right )^{2}. \end{eqnarray}

Proof. To prove this result, we use the Galerkin method. The spaces $H^1(\Omega )^d,$ $L^2(\Omega )$ are separable Hilbert spaces. Hence, there exist corresponding bases $\{\mathbf{W}_{i}\}_{i=1}^{\infty }$ and $\{q_i\}_{i=1}^{\infty }$ of smooth functions. Let $\mathbb{Y}_m$ be the space spanned by $\{(\mathbf{W}_{i},q_i)\}_{i=1}^{m}.$ The scalar product on $\mathbb{Y}_m$ is induced by the scalar product of $\mathbb{Y}.$ We define the approximate solution $(\mathbf{V}^m_f,P^m)$ as follows:

(3.16) \begin{equation} \mathbf{V}^m_f=\sum _{i=1}^{m}a_i\mathbf{W}_{i},\,\,\, P^m=\sum _{i=1}^{m}c_iq_i, \end{equation}

with

\begin{equation*}(Q_{m}(\boldsymbol {\varsigma }))\left \{ \begin {array}{ll} 2(D(\mathbf {V}^{m}_f)\,{:}\,D(\mathbf {W}))_{\Omega }+\lambda (\nabla \cdot \mathbf {V}^{m}_f,\nabla \cdot \mathbf {W})_{\Omega } -\phi _f(P^{m},\nabla \cdot \mathbf {W})_{\Omega }\\[3pt] +\frac {1}{Da}(\mathbf {K}(\boldsymbol {\varsigma })\mathbf {V}^{m}_f,\mathbf {W})_{\Omega } +\phi _f(\nabla \cdot \mathbf {V}^{m}_f,q)_{\Omega }+a_0(P^{m},q)_{\Omega } \\[3pt] =(\mathbf {b}_f,\mathbf {W})_{\Omega } +(\mathbf {T}_\infty,\mathbf {W})_{\partial \Omega }+(a_0,q)_{\Omega }, \end {array} \right .\end{equation*}

holding for all $(\mathbf{W},q)\in \mathbb{Y}_m$ with $a_i,\,b_i,\,c_i\in \mathbb{R},$ for $i\,=\,1,2,\ldots,m.$ The task now is to ensure the existence of solutions to $(Q_{m}(\boldsymbol{\varsigma }))$ and show that $(Q_{m}(\boldsymbol{\varsigma }))$ recovers $(Q_{w_1}(\boldsymbol{\varsigma }))$ as $m\rightarrow \infty .$ The linear structure of $(Q_{m}(\boldsymbol{\varsigma }))$ suggests that weak convergence is enough to pass the limit. Hence, in order to do so, we define a mapping $\mathcal{H}_m$ inspired by the structure of mapping $\mathcal{H}$ as

(3.17) \begin{equation} \langle \mathcal{H}^m_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P),(\mathbf{W},q)\rangle = \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P),(\mathbf{W},q)\rangle,\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}_m \end{equation}

where $\mathcal{H}_{\boldsymbol{\varsigma }}$ is defined in (3.6). From Lemma 3.1, we deduce that the mapping $\mathcal{H}^m_{\boldsymbol{\varsigma }}$ satisfies the conditions needed for Lemma A.1 (see Appendix A) and hence, there exists a solution $(\mathbf{V}^m_f,P^m)\in \mathbb{Y}_m$ for each $m$ such that

(3.18) \begin{equation} \langle \mathcal{H}^m_{\boldsymbol{\varsigma }}(\mathbf{V}^m_f,P^m),(\mathbf{W},q)\rangle =0, \,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}_m. \end{equation}

It follows that $(\mathbf{V}^m_f,P^m)$ satisfy $(Q_{m}(\boldsymbol{\varsigma }))$ and $a_i,$ $c_i$ can be determined.

Energy Estimates: Let $\mathbf{W}=\mathbf{V}^m_f,$ and $q=P^m$ in $(Q_{m}(\boldsymbol{\varsigma }))$ , then by performing calculations similar to those leading to (3.11), we obtain

(3.19) \begin{eqnarray} \alpha _3(||\mathbf{V}^m_f||^2_{1,\Omega } +||P^m||^2_{0,\Omega })-\alpha _4||(\mathbf{V}^m_f,P^m)||_{1,0,\Omega }\leqslant 0, \end{eqnarray}

where $\alpha _3$ and $\alpha _4$ are defined in (3.12)–(3.13). Consequently,

(3.20) \begin{eqnarray} ||\mathbf{V}^m_f||^2_{1,\Omega } +||P^m||^2_{0,\Omega }\leqslant \left (\frac{\alpha _4}{\alpha _3}\right )^{2}. \end{eqnarray}

Inequality (3.20) implies that the sequence $\{(\mathbf{V}^m_f,P^m)\}_{m\geqslant 1}$ is uniformly bounded in $\mathbb{Y}.$ Hence, it has a sub-sequence $\{(\mathbf{V}^m_f, P^m)\}_{m\geqslant 1}$ (for convenience, we denote it by the same symbol) and a pair $(\mathbf{V}_f, P)\in \mathbb{Y}$ such that

(3.21) \begin{eqnarray} (\mathbf{V}^m_f,P^m)\rightharpoonup (\mathbf{V}_f,P)\,\,\,\mbox{as}\, m\rightarrow \infty \,\,\mbox{weakly}\,\textrm{in}\,\mathbb{Y}. \end{eqnarray}

By taking the limit in (3.18) and using the weak convergence (3.21), we get

(3.22) \begin{equation} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P),(\mathbf{W},q)\rangle =0,\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}_m. \end{equation}

A continuity argument shows that (3.22) holds for any $(\mathbf{W},q)\in \mathbb{Y}.$ Hence, $(\mathbf{V}_f, P)$ is a solution of (3.7) and equivalently, of the weak formulation $(Q_{w_1}(\boldsymbol{\varsigma })).$ Using the lower semi-continuity property of norm in (3.20), we can achieve the following a priori bound on solution $(\mathbf{V}_f, P)$ given by

(3.23) \begin{eqnarray} ||\mathbf{V}_f||^2_{1,\Omega }+||P||^2_{0,\Omega }\leqslant \left (\frac{\alpha _4}{\alpha _3}\right )^{2}. \end{eqnarray}

Proposition 3.4. Suppose the hypotheses of Theorem 3.3 hold. Then, the weak formulation $(Q_{w_1}(\boldsymbol{\varsigma }))$ has a unique solution that depends continuously on the given data.

Proof. Uniqueness: Let $(\mathbf{V}^{1}_f,P^{1})$ and $(\mathbf{V}^{2}_f,P^{2})$ be two solutions that satisfy Equation (3.7) or equivalently the weak formulation $(Q_{w_1}(\boldsymbol{\varsigma })).$ Define $(\mathbf{V}_f,P)=(\mathbf{V}^{1}_f-\mathbf{V}^{2}_f,P^{1}-P^{2}).$ Then, from (3.7), we have

(3.24) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}^{1}_f,P^{1})-\mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}^{2}_f,P^{2}),(\mathbf{W},q) \rangle =0 \end{eqnarray}

for all $(\mathbf{W},q)\in \mathbb{Y}.$ Replace $(\mathbf{W},q)$ by $(\mathbf{V}_f,P)$ in (3.24) and using the definition of $\mathcal{H},$ we find

\begin{eqnarray*} \alpha ||\mathbf{V}_f||^2_{1,\Omega }+\lambda ||\nabla \cdot \mathbf{V}_f||^2_{0,\Omega } +{a_0}||P||^2_{0,\Omega }\leqslant 0. \end{eqnarray*}

The above implies $\mathbf{V}_f=0,$ $\mathbf{U}_s=0$ and $P=0$ almost everywhere in $\Omega .$ Hence, the weak formulation $(Q_{w_1}(\boldsymbol{\varsigma }))$ has a unique solution.

Continuous dependence: Let $(\mathbf{V}^1_f,P^1)$ and $(\mathbf{V}^2_f,P^2)$ be two solutions of $(Q_{w_1}(\boldsymbol{\varsigma }))$ corresponding to the two sets of data $(\mathbf{T}_{\infty,1},\mathbf{b}_{f,1},a_{0,1})$ and $(\mathbf{T}_{\infty,2},\mathbf{b}_{f,2},a_{0,2}),$ then the difference $(\mathbf{V}_f,P)=(\mathbf{V}^{1}_f-\mathbf{V}^{2}_f,P^{1}-P^{2})$ satisfies

\begin{eqnarray*} \alpha ||\mathbf{V}_f||^2_{1,\Omega }+\lambda ||\nabla \cdot \mathbf{V}_f||^2_{0,\Omega } +{a_0}||P||^2_{0,\Omega }\\ -(||\mathbf{b}_f||_{0,\Omega }+\sqrt{c_t}||\mathbf{T}_{\infty }||_{0,\partial \Omega })||\mathbf{V}||_{1,\Omega } -||a_0||_{0,\Omega }||P||_{0,\Omega }\leqslant 0, \end{eqnarray*}

or,

(3.25) \begin{align} \nonumber ||\mathbf{V}^1_f-\mathbf{V}^2_f||^2_{1,\Omega } +||P^1-P^2||^2_{0,\Omega } \\\leqslant \frac{1}{\alpha ^2_3}\left [(||\mathbf{b}_{f,1}-\mathbf{b}_{f,2}||_{0,\Omega } +\sqrt{c_t}||\mathbf{T}_{\infty,1}-\mathbf{T}_{\infty,2}||_{0,\partial \Omega })^2 +||a_{0,1}-a_{0,2}||^2_{0,\Omega }\right ]. \end{align}

Thus, if $(\mathbf{T}_{\infty,1},\mathbf{b}_{f,1}, a_{0,1})$ is close to $(\mathbf{T}_{\infty,2},\mathbf{b}_{f,2}, a_{0,2}),$ then the left-hand side of (3.25) (the difference of solutions) must be small. This establishes the well-posedness of the auxiliary linear problem $(Q_{w_1}(\boldsymbol{\varsigma })).$ Next, we would like to consider the sub-problem, $(Q_{w_2}(\nabla \cdot \mathbf{U}_s)).$

3.5 Case (a): existence and uniqueness results for $(Q_{w_2}(\mathbf{U}_s))$

We note that $\mathbf{K}(\boldsymbol{\varsigma })=\mathbf{K}(\mathbf{U}_s)$ is a non-linear function of $\mathbf{U}_s$ (see assumption ( $A4$ ) in Subsection 1.2) which makes $(Q_{w_2}(\mathbf{U}_s))$ a semilinear problem. By introducing a semilinear form $B(\cdot,\cdot ) \,{:}\, H^1_0(\Omega )^d\times H^1_0(\Omega )^d\rightarrow \mathbb{R}$ that is given by

(3.26) \begin{equation} B(\mathbf{U}_s,\mathbf{Z})=2\alpha _1(D(\mathbf{U}_s)\,{:}\,D(\mathbf{Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf{U}_s,\nabla \cdot \mathbf{Z})_{\Omega } -\frac{1}{Da}(\mathbf{K}(\mathbf{U}_s)\mathbf{V}_f,\mathbf{Z})_{\Omega } \end{equation}

and a linear form $L\,{:}\,H^1_0(\Omega )^d\rightarrow \mathbb{R}$ defined as

(3.27) \begin{equation} L(\mathbf{Z})=\phi _s(P,\nabla \cdot \mathbf{Z})_{\Omega }+ (\mathbf{b}_s,\mathbf{Z})_{\Omega } , \end{equation}

weak problem $(Q_{w_2}(\mathbf{U}_s))$ can be rewritten as an abstract formulation. For a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega ),$ find $\mathbf{U}_s\in H^1_0(\Omega )^d$ such that

(3.28) \begin{equation} B(\mathbf{U}_s,\mathbf{Z})=L(\mathbf{Z})\, \mbox{for all}\, \mathbf{Z}\in H^1_0(\Omega )^d. \end{equation}

In order to show existence and uniqueness results for problem (3.28), we will use the Browder-Minty theorem (see Thereom A.3, Appendix A), which is based on the monotone operator approach. To justify the hypotheses of the Browder-Minty theorem, we prove the following results in the form of lemmas.

Lemma 2. The correspondence $\mathbf{Z}\mapsto B(\mathbf{U}_s,\mathbf{Z})$ is a bounded linear operator and ${L}\in (H^1_0(\Omega )^d)^*.$

Proof. Clearly, the mapping $\mathbf{Z}\mapsto B(\mathbf{U}_s,\mathbf{Z})$ is linear (obvious) and bounded. Indeed, using Cauchy-Schwarz, Hölder’s and Sobolev inequalities, we find

(3.29) \begin{align} \nonumber |B(\mathbf{U}_s,\mathbf{Z})|\leqslant 2\alpha _1||D(\mathbf{U}_s)||_{0,\Omega }||D(\mathbf{Z})||_{0,\Omega } +\alpha _2||\nabla \cdot \mathbf{U}_s||_{0,\Omega }||\nabla \cdot \mathbf{Z}||_{0,\Omega } \\ \nonumber + \frac{1}{Da}||\mathbf{K}(\mathbf{U}_s)||_{L^\infty (\Omega )}||\mathbf{V}_f||_{0,\Omega }||\mathbf{Z}||_{0,\Omega } \leqslant \Big{(}(2\alpha _1+\alpha _2)||\nabla \mathbf{U}_s||_{0,\Omega } \\ +\frac{k_2\sqrt{c_p}}{Da}||\mathbf{V}_f||_{0,\Omega }\Big{)}||\nabla \mathbf{Z}||_{0,\Omega }. \end{align}

Now, $L$ is linear (obvious) and bounded. Indeed, we have

(3.30) \begin{eqnarray} |L(\mathbf{Z})| \leqslant (\phi _s||P||_{0,\Omega }+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega })||\nabla \mathbf{Z}||_{0,\Omega }. \end{eqnarray}

This implies $L\in (H^1_0(\Omega )^d)^*.$

The Lemma 2 implies that there exists an operator (non-linear) $\mathcal{A}\,{:}\, H^1_0(\Omega )^d\rightarrow (H^1_0(\Omega )^d)^*=H^{-1}(\Omega )^d$ with

(3.31) \begin{equation} (\mathcal{A}\mathbf{U}_s,\mathbf{Z})=B(\mathbf{U}_s,\mathbf{Z}). \end{equation}

Thus, the variational problem (3.28) equivalently reduces to the operator equation: find $\mathbf{U}_s\in H^1_0(\Omega )^d$ such that

(3.32) \begin{equation} \mathcal{A}\mathbf{U}_s=L \end{equation}

in the sense

(3.33) \begin{equation} (\mathcal{A}\mathbf{U}_s,\mathbf{Z})=L(\mathbf{Z}),\,\, \mbox{for all}\,\mathbf{Z}\in H^1_0(\Omega )^d. \end{equation}

Further, estimate (3.29) implies the non-linear operator $\mathcal{A}$ is bounded. Indeed,

(3.34) \begin{equation} ||\mathcal{A}\mathbf{U}_s||_{H^{-1}(\Omega )^d} \leqslant \left [(2\alpha _1+\alpha _2)||\nabla \mathbf{U}_s||_{0,\Omega }+\frac{k_2\sqrt{c_p}}{Da}||\mathbf{V}_f||_{0,\Omega }\right ]. \end{equation}

Lemma 3. If $\frac{2\alpha _1}{c_k}\gt \frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}$ , then the semilinear form $B(\cdot,\cdot )$ is elliptic that is there exists a constant $c\gt 0$ such that

(3.35) \begin{align} B(\mathbf{U}^1_s,\mathbf{U}^1_s-\mathbf{U}^2_s)- B(\mathbf{U}^2_s,\mathbf{U}^1_s-\mathbf{U}^2_s) \geqslant c||\nabla \mathbf{U}^1_s-\nabla \mathbf{U}^2_s||^2_{0,\Omega }\,\,\mbox{for all}\,\mathbf{U}^1_s,\,\mathbf{U}^2_s \in H^1_0(\Omega )^d. \end{align}

Note: This lemma implies that $\mathcal{A}$ is strongly monotone.

Proof: Indeed, consider

(3.36) \begin{align} \nonumber (\mathcal{A}\mathbf{U}^1_s-\mathcal{A}\mathbf{U}^2_s, \mathbf{U}^1_s-\mathbf{U}^2_s)= B(\mathbf{U}^1_s,\mathbf{U}^1_s-\mathbf{U}^2_s)- B(\mathbf{U}^2_s,\mathbf{U}^1_s-\mathbf{U}^2_s)\\ \nonumber = 2\alpha _1||D(\mathbf{U}^1_s)-D(\mathbf{U}^2_s)||^2_{0,\Omega }+ \alpha _2||\nabla \cdot \mathbf{U}^1_s-\nabla \cdot \mathbf{U}^2_s||^2_{0,\Omega } \\ \nonumber -\frac{1}{Da}((\mathbf{K}(\mathbf{U}^1_s)-\mathbf{K}(\mathbf{U}^2_s))\mathbf{V}_f,\mathbf{U}^1_s-\mathbf{U}^2_s)_{\Omega } \geqslant 2\alpha _1||D(\mathbf{U}^1_s)-D(\mathbf{U}^2_s)||^2_{0,\Omega } \\ \nonumber + \alpha _2||\nabla \cdot \mathbf{U}^1_s-\nabla \cdot \mathbf{U}^2_s||^2_{0,\Omega } -\frac{1}{Da}||\mathbf{K}(\mathbf{U}^1_s)-\mathbf{K}(\mathbf{U}^2_s)||_{0,\Omega }||\mathbf{V}_f||_{L^4(\Omega )}||\mathbf{U}^1_s-\mathbf{U}^2_s||_{L^4(\Omega )} \\ \nonumber \geqslant \frac{2\alpha _1}{c_k}||\nabla \mathbf{U}^1_s-\nabla \mathbf{U}^2_s||^2_{0,\Omega } -\frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}||\nabla \mathbf{U}^1_s-\nabla \mathbf{U}^2_s||^2_{0,\Omega } \\ = \left (\frac{2\alpha _1}{c_k} -\frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )||\nabla \mathbf{U}^1_s-\nabla \mathbf{U}^2_s||^2_{0,\Omega } \end{align}

To reach (3.36), we have applied Hölder’s, Poincare’s and Sobolev’s inequalities and Lipschitz continuous property of $\mathbf{K}.$ Thus, if $\frac{2\alpha _1}{c_k}\gt \frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}$ , then $B$ is elliptic with $c=\left (\frac{2\alpha _1}{c_k} -\frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}\right ).$

Lemma 4. The non-linear operator $\mathcal{A}$ as in (3.31) is continuous from $H^1_0(\Omega )^d$ to $H^{-1}(\Omega )^d.$

Proof. Let $\mathbf{U}^n_s \rightarrow \mathbf{U}_s$ in $H^1_0(\Omega )^d$ as $n\rightarrow \infty .$ Consider

\begin{align*} |(\mathcal{A}\mathbf{U}^n_s-\mathcal{A}\mathbf{U}_s,\mathbf{Z})| \leqslant 2\alpha _1||D(\mathbf{U}^n_s)-D(\mathbf{U}_s)||_{0,\Omega }||D(\mathbf{Z})||_{0,\Omega } \\ + \alpha _2||\nabla \cdot \mathbf{U}^n_s-\nabla \cdot \mathbf{U}_s||_{0,\Omega } ||\nabla \cdot \mathbf{Z}||_{0,\Omega } +\frac{1}{Da}||\mathbf{K}(\mathbf{U}^n_s)-\mathbf{K}(\mathbf{U}_s)||_{0,\Omega }||\mathbf{V}_f||_{L^4(\Omega )}||\mathbf{Z}||_{L^4(\Omega )} \end{align*}

On applying Lipschitz continuity of $\mathbf{K}$ and Sobolev inequality, we find

\begin{align*} ||\mathcal{A}\mathbf{U}^n_s-\mathcal{A}\mathbf{U}_s||_{H^{-1}(\Omega )^d} \leqslant 2\alpha _1||D(\mathbf{U}^n_s)-D(\mathbf{U}_s)||_{0,\Omega }+ \alpha _2||\nabla \cdot \mathbf{U}^n_s-\nabla \cdot \mathbf{U}_s||_{0,\Omega } \\ +\frac{k_Lc_s}{Da}||\mathbf{U}^n_s-\mathbf{U}_s||_{0,\Omega } \end{align*}

and $||\mathcal{A}\mathbf{U}^n_s-\mathcal{A}\mathbf{U}_s||_{H^{-1}(\Omega )^d}\rightarrow 0$ as $n\rightarrow \infty .$ That is, $\mathcal{A}$ is continuous.

Theorem 3.5. Suppose that assumptions (A), (B) and (C) hold. Further, if the given data and non-dimensional parameters satisfy the following assumption

(3.37) \begin{equation} \frac{2\alpha _1}{c_k}\gt \frac{k_Lc_s\sqrt{c_p}\alpha _4}{\alpha _3Da}, \end{equation}

then for a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega )^d$ , the problem $(Q_{w_2}(\mathbf{U}_s))$ has a unique solution $\mathbf{U}_s\in H^1_0(\Omega )^d.$ Moreover, if

(3.38) \begin{equation} \frac{2\alpha _1}{c_k}\gt \frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da} \end{equation}

holds, then $\mathbf{U}_s$ satisfies the following a priori estimate

(3.39) \begin{eqnarray} ||\nabla \mathbf{U}_s||_{0,\Omega } \leqslant \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )} \left (\frac{\phi _s\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right ). \end{eqnarray}

Proof. The analysis shown in Lemmas 3 and 4 justifies that $\mathcal{A}$ satisfies the hypothesis of the Browder-Minty theorem (see Theorem A.3, Appendix A). Consequently, the operator Equation (3.32) or the problem (3.28) has a unique solution $\mathbf{U}_s\in H^1_0(\Omega )^d$ for any given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega ).$ Further, if (3.38) holds, then $\mathbf{U}_s$ satisfies the apriory estimate (3.39). Indeed, from (3.28) replacing $\mathbf{Z}=\mathbf{U}_s$ and using Poincare’s inequality, we have

\begin{eqnarray*} B(\mathbf{U}_s,\mathbf{U}_s)=L(\mathbf{U}_s)\leqslant (\phi _s||P||_{0,\Omega }+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega })||\nabla \mathbf{U}_s||_{0,\Omega } \end{eqnarray*}

Making use of the definition of $B(\mathbf{U}_s,\mathbf{U}_s)$ and Korn’s, Hölder’s and Sobolev inequalities and boundedness property of $\mathbf{K}$ , we obtain

\begin{eqnarray*} ||\nabla \mathbf{U}_s||_{0,\Omega } \leqslant \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )} \left (\frac{\phi _s\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right ). \end{eqnarray*}

The analysis in Subsections 3.4–3.5 describes the existence and uniqueness of problems $(Q_{w_1}(\boldsymbol{\varsigma }))-(Q_{w_2}(\mathbf{U}_s)).$ In the next subsections, we focus on developing existence and uniqueness results corresponding to $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ that is the case (b).

3.6 Case (b): existence and uniqueness of a solution to $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))$

Analysis in this subsection is analogous to the arguments in the Subsection 3.4. Thus, we only state the main theorem and omit the proof.

Theorem 3.6. Suppose that the assumptions (A) and (B) hold. Then for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ , the problem $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))$ has a unique solution $(\mathbf{V}_f,P)\in \mathbb{Y}=H^1(\Omega )^d\times L^2(\Omega )$ that depends continuously on the given data. Moreover, the solution $(\mathbf{V}_f,P)$ satisfies the following a priori bound

(3.40) \begin{eqnarray} ||\mathbf{V}_f||^2_{1,\Omega }+||P||^2_{0,\Omega }\leqslant \left (\frac{\alpha _4}{\alpha _3}\right )^{2}. \end{eqnarray}

Proof. Proof of this theorem follows from Theorem 3.3 and Proposition 3.4.

3.7 Case (b): existence and uniqueness of a solution to $(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$

We note that $\mathbf{K}(\nabla \cdot \mathbf{U}_s)$ is a non-linear function of $\nabla \cdot \mathbf{U}_s$ (see assumption ( $A4$ ) in Subsection 1.2) which makes $(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ a semilinear problem.

By introducing a semilinear form $B(\cdot,\cdot ) \,{:}\, H^1_0(\Omega )^d\times H^1_0(\Omega )^d\rightarrow \mathbb{R}$ that is given by

(3.41) \begin{equation} B(\mathbf{U}_s,\mathbf{Z})=2\alpha _1(D(\mathbf{U}_s)\,{:}\,D(\mathbf{Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf{U}_s,\nabla \cdot \mathbf{Z})_{\Omega } -\frac{1}{Da}(\mathbf{K}(\nabla \cdot \mathbf{U}_s)\mathbf{V}_f,\mathbf{Z})_{\Omega } \end{equation}

and a linear form $L\,{:}\,H^1_0(\Omega )^d\rightarrow \mathbb{R}$ defined as

(3.42) \begin{equation} L(\mathbf{Z})=\phi _s(P,\nabla \cdot \mathbf{Z})_{\Omega }+ (\mathbf{b}_s,\mathbf{Z})_{\Omega }, \end{equation}

weak problem $(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ can be rewritten as an abstract formulation. For a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega ),$ find $\mathbf{U}_s\in H^1_0(\Omega )^d$ such that

(3.43) \begin{equation} B(\mathbf{U}_s,\mathbf{Z})=L(\mathbf{Z})\, \mbox{for all}\, \mathbf{Z}\in H^1_0(\Omega )^d. \end{equation}

Similar to Subsection 3.5 to show existence and uniqueness results for problem (3.43), we will use the Browder-Minty theorem (see Theorem A.3, Appendix A), which is based on the monotone operator approach. We state the main theorem without proof, which can be done similarly to the proof of Theorem 3.5.

Theorem 3.7. Suppose that assumptions (A), (B) and (C) hold. Further, if the given data and non-dimensional parameters satisfy the following assumption

(3.44) \begin{equation} \frac{2\alpha _1}{c_k}\gt \frac{k_Lc_s\alpha _4}{\alpha _3Da}, \end{equation}

then for a given pair $(\mathbf{V}_f,P)\in H^1(\Omega )^d \times L^2(\Omega )^d$ , the problem $(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ has a unique solution $\mathbf{U}_s\in H^1_0(\Omega )^d.$ Moreover, if

(3.45) \begin{equation} \frac{2\alpha _1}{c_k}\gt \frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da} \end{equation}

holds, then $\mathbf{U}_s$ satisfies the following a priori estimate

(3.46) \begin{eqnarray} ||\nabla \mathbf{U}_s||_{0,\Omega } \leqslant \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )} \left (\frac{\phi _s\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right ). \end{eqnarray}

Proof. The proof of this theorem can be done similarly to the proof of Theorem 3.5. The analysis in Subsections 3.6–3.7 completes the existence and uniqueness of solution to problem $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s)).$ In the next section, we focus on the development of existence and uniqueness results corresponding to coupled non-linear problems $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ and $(Q_{w_1}(\nabla \cdot \mathbf{U}_s))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s)),$ respectively, by converting them into a fixed-point problem.

4. Case (a): reduction to a fixed-point problem for $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$

We note that for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d,$ the problem $(Q_{w_1}(\boldsymbol{\varsigma }))$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ (see Subsection 3.4). Consequently, we can define a mapping $T_1\,{:}\, H^1_0(\Omega )^d\rightarrow H^1(\Omega )^d\times L^2(\Omega )$ such that $T_1(\boldsymbol{\varsigma })=(\mathbf{V}_f, P).$ Further, for a given pair $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega ),$ the problem $(Q_{w_2}(\mathbf{U}_s))$ has a unique solution $\mathbf{U}_s\in H^1_0(\Omega )^d$ (see Subsection 3.5). Therefore, we can define a mapping $T_2\,{:}\, H^1(\Omega )^d\times L^2(\Omega )\rightarrow H^1_0(\Omega )^d$ such that $T_2(\mathbf{V}_f, P)=\mathbf{U}_s.$ Now, in order to get the fixed-point problem corresponding to $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s)),$ we define a composition map $T=T_2\circ T_1\,{:}\, H^1_0(\Omega )^d \rightarrow H^1_0(\Omega )^d$ such that

(4.1) \begin{equation} T(\boldsymbol{\varsigma })=(T_2\circ T_1)(\boldsymbol{\varsigma })=T_2(T_1(\boldsymbol{\varsigma }))=T_2(\mathbf{V}_f, P)=\mathbf{U}_s. \end{equation}

Thus, a fixed-point of mapping $T$ solves the coupled non-linear problem $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ or equivalently, (3.3)–(3.5) when $\boldsymbol{\varsigma }=\mathbf{U}_s.$ In order to show that the mapping $T$ has a fixed point, we use Schauder’s fixed-point theorem (see Thereom A.2, Appendix A).2 Thus, in the following analysis, we prove some results in the form of lemmas to verify the hypotheses of Schauder’s fixed-point theorem.

4.1 Analysis of the fixed-point problem

Throughout this subsection, we assume that hypotheses of Theorem 3.3, Proposition 3.4 and Theorem 3.5 hold. Then, we have

Lemma 5. Given $r\gt 0,$ let $\mathbf{W}_r$ be a closed and convex subset of $H^1(\Omega )^d$ defined by

\begin{equation*}\mathbf {W}_r=\{\mathbf {w}\in H^1_0(\Omega )^d \,{:}\, ||\mathbf {w}||_{H^1_0(\Omega )^d}=||\nabla \mathbf {w}||_{0,\Omega }\leqslant r\},\end{equation*}

and assume that the data satisfies

(4.2) \begin{equation} \frac{1}{\left (\frac{2\alpha _1}{c_k}- \frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )}\left (\phi _s\frac{\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right )\leqslant r. \end{equation}

Then, ${T}(\mathbf{W}_r)\subseteq \mathbf{W}_r.$

Proof. For any $\mathbf{w}\in \mathbf{W}_r,$ using estimate (3.39), we find

(4.3) \begin{align} \nonumber ||T(\mathbf{w})||_{H^1_0(\Omega )^d}=||T_2(T_1(\mathbf{w}))||_{H^1_0(\Omega )^d}=||T_2(\mathbf{V}_f,P)||_{H^1_0(\Omega )^d}=||\mathbf{U}_s||_{H^1_0(\Omega )^d}\\ = ||\nabla \mathbf{U}_s||_{0,\Omega }\leqslant \frac{1}{\left (\frac{2\alpha _1}{c_k}- \frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )}\left (\phi _s\frac{\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right ). \end{align}

The above estimate (4.3) together with the assumption (4.2) imply $T(\mathbf{w}) \in \mathbf{W}_r,$ which proves $T(\mathbf{W}_r)\subseteq \mathbf{W}_r.$

Lemma 6. The map ${T}_1 \,{:}\, H^1_0(\Omega )^d \rightarrow H^1(\Omega )^d\times L^2(\Omega )$ satisfies

(4.4) \begin{eqnarray} ||{T}_1(\boldsymbol{\varsigma })-{T}_1(\tilde{\boldsymbol{\varsigma }})||_{\mathbb{Y}} \leqslant \frac{\sqrt{2}k_L\alpha _4c_s}{\alpha ^2_3Da}||\boldsymbol{\varsigma } -\tilde{\boldsymbol{\varsigma }}||_{0,\Omega }. \end{eqnarray}

Proof. Given $\boldsymbol{\varsigma }, \tilde{\boldsymbol{\varsigma }}\in H^1_0(\Omega )^d,$ let $(\mathbf{V}_f, P),\,(\tilde{\mathbf{V}}_f, \tilde{P})\in H^1(\Omega )^d\times L^2(\Omega )$ be the corresponding solutions of $(Q_{w_1}(\boldsymbol{\varsigma })).$ Then, Equation (3.7) implies

(4.5) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P),(\mathbf{W},q)\rangle =0\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}, \end{eqnarray}

(4.6) \begin{eqnarray} \langle \mathcal{H}_{\tilde{\boldsymbol{\varsigma }}}(\tilde{\mathbf{V}}_f, \tilde{P}),(\mathbf{W},q)\rangle =0\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}. \end{eqnarray}

Taking the difference between the above equations, we get

(4.7) \begin{eqnarray} \langle \mathcal{H}_{\boldsymbol{\varsigma }}(\mathbf{V}_f,P)-\mathcal{H}_{\tilde{\boldsymbol{\varsigma }}}(\tilde{\mathbf{V}}_f, \tilde{P}),(\mathbf{W},q)\rangle =0\,\mbox{for all}\, (\mathbf{W},q)\in \mathbb{Y}. \end{eqnarray}

Replacing $\mathbf{W}=\mathbf{V}_f-\tilde{\mathbf{V}}_f$ and $q=P-\tilde{P}$ in the above equation and using the definition of $\mathcal{H}_{\boldsymbol{\varsigma }}$ and $\mathcal{H}_{\tilde{\boldsymbol{\varsigma }}},$ we get

(4.8) \begin{eqnarray} 2||D(\mathbf{V}_f)-D(\tilde{\mathbf{V}}_f)||^2_{0,\Omega }+\lambda ||\nabla \cdot \mathbf{V}_f-\nabla \cdot \tilde{\mathbf{V}}_f||^2_{0,\Omega } +a_0||P-\tilde{P}||^2_{0,\Omega }\nonumber\\=-\frac{1}{Da}(\mathbf{K}(\boldsymbol{\varsigma })(\mathbf{V}_f-\tilde{\mathbf{V}}_f),\mathbf{V}_f-\tilde{\mathbf{V}}_f)_{\Omega }-\frac{1}{Da}((\mathbf{K}(\boldsymbol{\varsigma })-\mathbf{K}(\tilde{\boldsymbol{\varsigma }}))\tilde{\mathbf{V}}_f,\mathbf{V}_f-\tilde{\mathbf{V}}_f)_{\Omega }. \end{eqnarray}

Using Korn’s, Hölder’s inequalities and estimate (3.23), we obtain

\begin{eqnarray*} ||\mathbf{V}_f-\tilde{\mathbf{V}}_f||^2_{1,\Omega }+||P-\tilde{P}||^2_{0,\Omega } \leqslant \frac{\alpha _4}{\alpha ^2_3 Da}||\mathbf{K}(\boldsymbol{\varsigma }) -\mathbf{K}(\tilde{\boldsymbol{\varsigma }})||_{0,\Omega }||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{L^4(\Omega )}. \end{eqnarray*}

Further, using Lipschitz continuous property of $\mathbf{K}$ and Sobolev’s inequality, we get

\begin{eqnarray*} \left (||\mathbf{V}_f-\tilde{\mathbf{V}}_f||^2_{1,\Omega }+||P-\tilde{P}||^2_{0,\Omega }\right ) \leqslant \frac{k_L\alpha _4{c_s}}{\alpha ^2_3Da}||\boldsymbol{\varsigma } -\tilde{\boldsymbol{\varsigma }}||_{0,\Omega }||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{1,\Omega }, \end{eqnarray*}

or,

(4.9) \begin{align} ||(\mathbf{V}_f, P)-(\tilde{\mathbf{V}}_f, \tilde{P})||_{\mathbb{Y}}=\left (||\mathbf{V}_f-\tilde{\mathbf{V}}_f||^2_{1,\Omega }+||P-\tilde{P}||^2_{0,\Omega }\right )^{1/2} \leqslant \frac{\sqrt{2}k_L\alpha _4c_s}{\alpha ^2_3Da}||\boldsymbol{\varsigma } -\tilde{\boldsymbol{\varsigma }}||_{0,\Omega }. \end{align}

This establishes (4.4).

Lemma 7. The map ${T}_2 \,{:}\, H^1(\Omega )^d\times L^2(\Omega )\rightarrow H^1_0(\Omega )^d$ satisfies

(4.10) \begin{align} ||{T}_2(\mathbf{V}_f, P)-{T}_2(\tilde{\mathbf{V}}_f, \tilde{P})||_{H^1_0(\Omega )} \leqslant{\beta }\left [||P-\tilde{P}||_{0,\Omega } +||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{1,\Omega }\right ] \end{align}

where ${\beta }= \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_Lc_s\alpha _4\sqrt{c_p}}{\alpha _3Da}\right )}\max \left \{\phi _s, \frac{k_2\sqrt{c_p}}{Da}\right \}.$

Proof. Given $(\mathbf{V}_f, P),\,(\tilde{\mathbf{V}}_f, \tilde{P})\in H^1(\Omega )^d\times L^2(\Omega ),$ let $\mathbf{U}_s,\,\tilde{\mathbf{U}}_s\in H^1_0(\Omega )^d$ be the corresponding solutions of $(Q_{w_2}(\mathbf{U}_s)).$ Then, (3.28) implies

(4.11) \begin{eqnarray} B(\mathbf{U}_s,\mathbf{Z})=L(\mathbf{Z})\, \mbox{for all}\, \mathbf{Z}\in H^1_0(\Omega )^d, \end{eqnarray}

(4.12) \begin{eqnarray} B(\tilde{\mathbf{U}}_s,\mathbf{Z})=L(\mathbf{Z})\, \mbox{for all}\, \mathbf{Z}\in H^1_0(\Omega )^d. \end{eqnarray}

Taking the difference between the above equations, we get

\begin{align*} 2\alpha _1(D(\mathbf{U}_s)-D(\tilde{\mathbf{U}}_s)\,{:}\,D(\mathbf{Z}))_{\Omega }+\alpha _2(\nabla \cdot \mathbf{U}_s-\nabla \cdot \tilde{\mathbf{U}}_s,\nabla \cdot \mathbf{Z})_{\Omega }\\ =\phi _s(P-\tilde{P},\nabla \cdot \mathbf{Z})_{\Omega } +\frac{1}{Da}(\mathbf{K}(\mathbf{U}_s)(\mathbf{V}_f-\tilde{\mathbf{V}}_f),\mathbf{Z})_{\Omega } \\ +\frac{1}{Da}((\mathbf{K}(\mathbf{U}_s)-\mathbf{K}(\tilde{\mathbf{U}}_s))\tilde{\mathbf{V}}_f,\mathbf{Z})_{\Omega }. \end{align*}

Replace $\mathbf{Z}=\mathbf{U}_s -\tilde{\mathbf{U}}_s$ and using Cauchy-Schwarz and Hölder’s inequalities, we get

\begin{align*} 2\alpha _1||D(\mathbf{U}_s)-D(\tilde{\mathbf{U}}_s)||^2_{0,\Omega }+\alpha _2||\nabla \cdot \mathbf{U}_s-\nabla \cdot \tilde{\mathbf{U}}_s||^2_{0,\Omega } \leqslant \phi _s||P-\tilde{P}||_{0,\Omega }||\nabla \cdot (\mathbf{U}_s -\tilde{\mathbf{U}}_s)||_{0,\Omega } \\ +\frac{1}{Da}||\mathbf{K}(\mathbf{U}_s)||_{L^\infty (\Omega )}||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{0,\Omega }||\mathbf{U}_s -\tilde{\mathbf{U}}_s||_{0,\Omega } \\ +\frac{1}{Da}||\mathbf{K}(\mathbf{U}_s)-\mathbf{K}(\tilde{\mathbf{U}}_s)||_{0,\Omega }||\tilde{\mathbf{V}}_f||_{L^4(\Omega )}||\mathbf{U}_s -\tilde{\mathbf{U}}_s||_{L^4(\Omega )}. \end{align*}

Moreover, using Lipschitz property of $\mathbf{K}$ , Korn’s, Poincare’s and Sobolev’s embedding inequalities and estimate (3.23), we get

\begin{align*} ||T_2(\mathbf{V}_f,P)-T_2(\tilde{\mathbf{V}}_f,\tilde{P})||_{H^1_0(\Omega )^d}=||\mathbf{U}_s-\tilde{\mathbf{U}}_s||_{H^1_0(\Omega )^d}= ||\nabla \mathbf{U}_s -\nabla \tilde{\mathbf{U}}_s||_{0,\Omega }\\ \leqslant \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_Lc_s\alpha _4\sqrt{c_p}}{\alpha _3Da}\right )}\left [\phi _s||P-\tilde{P}||_{0,\Omega } +\frac{k_2\sqrt{c_p}}{Da}||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{0,\Omega }\right ], \end{align*}

or,

(4.13) \begin{align} ||T_2(\mathbf{V}_f,P)-T_2(\tilde{\mathbf{V}}_f,\tilde{P})||_{H^1_0(\Omega )^d} \leqslant \beta \left [||P-\tilde{P}||_{0,\Omega } +||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{1,\Omega }\right ] \end{align}

where ${\beta }= \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_Lc_s\alpha _4\sqrt{c_p}}{\alpha _3Da}\right )}\max \left \{\phi _s, \frac{k_2\sqrt{c_p}}{Da}\right \}.$

Lemma 8. The map ${T}={T}_2\circ{T}_1 \,{:}\, H^1_0(\Omega )^d \subset L^2(\Omega )^d\rightarrow H^1_0(\Omega )^d$ satisfies

(4.14) \begin{align} ||{T}(\boldsymbol{\varsigma })-{T}(\tilde{\boldsymbol{\varsigma }})||_{H^1_0(\Omega )^d}\leqslant \frac{2{\beta } k_L\alpha _4{c_s}}{\alpha ^2_3Da}||\boldsymbol{\varsigma } -\tilde{\boldsymbol{\varsigma }}||_{0,\Omega }\leqslant \frac{2{\beta } k_L\alpha _4c_s\sqrt{c_p}}{\alpha ^2_3Da}||\nabla \boldsymbol{\varsigma } -\nabla \tilde{\boldsymbol{\varsigma }}||_{0,\Omega }. \end{align}

Proof. From (4.13), we have

(4.15) \begin{align} \nonumber ||T(\boldsymbol{\varsigma })-T(\tilde{\boldsymbol{\varsigma }})||_{H^1_0(\Omega )^d)}= ||T_2(\mathbf{V}_f,P)-T_2(\tilde{\mathbf{V}}_f,\tilde{P})||_{H^1_0(\Omega )^d} \\ \leqslant \beta \left [||P-\tilde{P}||_{0,\Omega } +||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{1,\Omega }\right ] \end{align}

We achieve (4.14) with the help of (4.9) and (4.15).

Theorem 4.1. The mapping ${T} \,{:}\, \mathbf{W}_r\subset H^1_0(\Omega )^d \rightarrow \mathbf{W}_r\subset H^1_0(\Omega )^d$ is continuous and $\overline{{T}(\mathbf{W}_r)}$ is compact.

Proof. The continuity of $T$ follows in a straightforward manner from (4.14). Now, given a sequence $\{\boldsymbol{\varsigma }_k\}_{k\in \mathbb{N}}$ of $\mathbf{W}_r$ which is clearly bounded, there exists a sub-sequence $\{\boldsymbol{\varsigma }_{k_j}\}\subset \{\boldsymbol{\varsigma }_k\}_{k\in \mathbb{N}}$ and $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d$ such that $\boldsymbol{\varsigma }_{k_j}\rightharpoonup \boldsymbol{\varsigma }$ in $H^1_0(\Omega )^d.$ In this way, thanks to the compact embedding of $H^1_0(\Omega )^d$ in $L^2(\Omega )^d,$ which implies $\boldsymbol{\varsigma }_{k_j} \rightarrow \boldsymbol{\varsigma }$ in $L^2(\Omega )^d,$ which combined with (4.14) implies that $T(\boldsymbol{\varsigma }_{k_j})\rightarrow T(\boldsymbol{\varsigma }).$ This proves that $\overline{T(\mathbf{W}_r)}$ is compact.

Theorem 4.2. Suppose that the hypotheses of Theorem 3.3, Proposition 3.4 and Theorem 3.5 hold. Then, the mapping $T$ defined in 4.1 has a fixed point $\mathbf{U}_s\in H^1_0(\Omega )^d$ , which in turn implies that coupled problem $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ has a solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$ Further, if

(4.16) \begin{equation} \frac{2{\beta } k_L\alpha _4{c_s}\sqrt{c_p}}{\alpha ^2_3Da}\lt 1 \end{equation}

then $T$ has a unique fixed point $ \mathbf{U}_s\in H^1_0(\Omega )^d$ , which in turn implies that coupled problem $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$

Proof. The above Theorem 4.1 implies that $T$ satisfies all hypotheses of Schauder’s fixed-point theorem (see Thereom A.2, Appendix A). Consequently, the mapping $T$ has a fixed-point $\mathbf{U_s}\in H^1_0(\Omega ).$ This implies that coupled problem $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ or equivalently, (3.3)–(3.5) when $\boldsymbol{\varsigma }=\mathbf{U}_s$ has a solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$ Further, if (4.12) is true, then $T\,{:}\, \mathbf{W}_r \rightarrow \mathbf{W}_r$ is a strict contraction mapping. That implies $T$ has a unique fixed point $ \mathbf{U}_s\in H^1_0(\Omega )^d$ due to Banach’s fixed-point theorem (see p. 415 [34]). This implies that coupled problem $(Q_{w_1}(\mathbf{U}_s))-(Q_{w_2}(\mathbf{U}_s))$ or equivalently, (3.3)–(3.5) when $\boldsymbol{\varsigma }=\mathbf{U}_s$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$

Remark 4.3 (Continuous dependence). If the non-dimensional parameters and constants satisfy the following assumption

(4.17) \begin{align} 2\alpha ^* Da\gt \frac{k_L\alpha _4c_s}{\alpha _3},\, \frac{4\alpha _1Da}{c_k}\gt \left (\frac{c_pk^2_2}{k_1}+\frac{k_L\alpha _4c_s(c_p+{2}\sqrt{c_p})}{\alpha _3}\right ),\, 2\alpha _2\geqslant \frac{\phi ^2_s}{a_0}, \end{align}

where $\alpha ^*=\frac{1}{c_k}\min \{2,\frac{k_1}{2Da}\}.$ Then, the continuous dependence (3.25) holds. Indeed, we have

(4.18) \begin{align} \nonumber ||\mathbf{V}^1_f-\mathbf{V}^2_f||^2_{1,\Omega } +||\nabla (\mathbf{U}^1_s- \mathbf{U}^2_s)||^2_{0,\Omega } +||P^1-P^2||^2_{0,\Omega } \leqslant \frac{1}{\alpha ^2_6}[(||\mathbf{b}_{f,1}-\mathbf{b}_{f,2}||_{0,\Omega } \\ +\sqrt{c_t}||\mathbf{T}_{\infty,1}-\mathbf{T}_{\infty,2}||_{0,\partial \Omega })^2 +||a_{0,1}-a_{0,2}||^2_{0,\Omega }+c_p||\mathbf{b}_{s,1}-\mathbf{b}_{s,2}||^2_{0,\Omega }], \end{align}

where $\alpha _5=\min \left \{\alpha ^*-\frac{k_L\alpha _4c_s}{2\alpha _3Da}, \left (\frac{2\alpha _1}{c_k}-\frac{c_pk^2_2}{2k_1Da}-\frac{k_L\alpha _4c_s(c_p+{2}\sqrt{c_p})}{2\alpha _3Da} \right ), \frac{a_0}{2}\right \}.$

5. Case (b): reduction to a fixed-point problem to $(Q_{w_1}(\nabla \cdot \mathbf{U}_s))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$

We note that for a given $\boldsymbol{\varsigma }\in H^1_0(\Omega )^d,$ the problem $(Q_{w_1}(\nabla \cdot \boldsymbol{\varsigma }))$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ (see Subsection 3.6). Consequently, we can define a mapping $T_1\,{:}\, H^1_0(\Omega )^d \rightarrow H^1(\Omega )^d\times L^2(\Omega )$ such that $T_1(\boldsymbol{\varsigma })=(\mathbf{V}_f, P).$ Further, for a given pair $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega ),$ the problem $(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ has a unique solution $\mathbf{U}_s\in H^1_0(\Omega )^d$ (see Subsection 3.7). Therefore, we can define a mapping $T_2\,{:}\, H^1(\Omega )^d\times L^2(\Omega )\rightarrow H^1_0(\Omega )^d$ such that $T_2(\mathbf{V}_f, P)=\mathbf{U}_s.$ Now, in order to get the fixed-point problem corresponding to $(Q_{w_1}(\nabla \cdot \mathbf{U}_s))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s)),$ we define a composition map $T=T_2\circ T_1\,{:}\, H^1_0(\Omega )^d \rightarrow H^1_0(\Omega )^d$ such that

(5.1) \begin{equation} T(\boldsymbol{\varsigma })=T_2(T_1(\boldsymbol{\varsigma }))=T_2(\mathbf{V}_f, P)=\mathbf{U}_s. \end{equation}

Thus, a fixed-point of mapping $T$ solves the coupled non-linear problem $(Q_{w_1}((\nabla \cdot \mathbf{U}_s))$ and $(Q_{w_2}((\nabla \cdot \mathbf{U}_s))$ or equivalently, (3.3)–(3.5) when $\boldsymbol{\varsigma }=\nabla \cdot \mathbf{U}_s.$ In order to show that the mapping $T$ has a fixed point, we use Schauder’s fixed-point theorem (see Thereom A.2, Appendix A). Thus, in the following analysis, we state some results in the form of lemmas to verify the hypotheses of Schauder’s fixed-point theorem. The proof of the following results is almost similar to the proof presented in Subsection 4.1, so we omit the details.

5.1 Analysis of the fixed-point problem

Throughout this subsection, we assume that hypotheses of Theorems 3.6 and 3.7 hold. Then, we have

Lemma 9. Given $r\gt 0,$ let $\mathbf{W}_r$ be a closed and convex subset of $H^1(\Omega )^d$ defined by

\begin{equation*}\mathbf {W}_r=\{\mathbf {w}\in H^1_0(\Omega )^d \,{:}\, ||\mathbf {w}||_{H^1_0(\Omega )^d}=||\nabla \mathbf {w}||_{0,\Omega }\leqslant r\}\end{equation*}

and assume that the data satisfy

(5.2) \begin{equation} \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_2\sqrt{c_p}\alpha _4}{\alpha _3Da}\right )} \left (\frac{\phi _s\alpha _4}{\alpha _3}+\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }\right )\leqslant r. \end{equation}

Then, $T(\mathbf{W}_r)\subseteq \mathbf{W}_r.$

Proof. This lemma is proved similar to Lemma 5.

Lemma 10. The map $T_1 \,{:}\, H^1_0(\Omega )^d \rightarrow H^1(\Omega )^d\times L^2(\Omega )$ satisfies

(5.3) \begin{eqnarray} ||T_1(\boldsymbol{\varsigma })-T_1(\tilde{\boldsymbol{\varsigma }})||_{\mathbb{Y}} \leqslant \frac{\sqrt{2}k_L\alpha _4{c_s}}{\alpha ^2_3Da}||\nabla \cdot \boldsymbol{\varsigma } -\nabla \cdot \tilde{\boldsymbol{\varsigma }}||_{0,\Omega }, \end{eqnarray}

(see Lemma 3.1 for the definition of the space $\mathbb{Y}$ ).

Proof. This lemma is proved similar to Lemma 6.

Lemma 11. The map $T_2\,{:}\, H^1(\Omega )^d\times L^2(\Omega )\rightarrow H^1_0(\Omega )^d$ satisfies

(5.4) \begin{align} ||T_2(\mathbf{V}_f, P)-T_2(\tilde{\mathbf{V}}_f, \tilde{P})||_{H^1_0(\Omega )^d} \leqslant \tilde{\beta }\left [||P-\tilde{P}||_{0,\Omega } +||\mathbf{V}_f-\tilde{\mathbf{V}}_f||_{1,\Omega }\right ], \end{align}

where $\tilde{\beta }= \frac{1}{\left (\frac{2\alpha _1}{c_k}-\frac{k_L\alpha _4c_s}{\alpha _3Da}\right )}\max \left \{\phi, \frac{k_2\sqrt{c_p}}{Da} \right \}$ .

Proof. This lemma is proved similar to Lemma 7.

Lemma 12. The map $T=T_2\circ T_1 \,{:}\, H^1_0(\Omega )^d \rightarrow H^1_0(\Omega )^d$ satisfies

(5.5) \begin{align} ||T(\boldsymbol{\varsigma })-T(\tilde{\boldsymbol{\varsigma }})||_{H^1_0(\Omega )^d)} \leqslant \frac{2\tilde{\beta } k_L\alpha _4{c_s}}{\alpha ^2_3Da}||\nabla \cdot \boldsymbol{\varsigma } -\nabla \cdot \tilde{\boldsymbol{\varsigma }}||_{0,\Omega }\leqslant \frac{2 \tilde{\beta } k_L\alpha _4c_s}{\alpha ^2_3Da}||\nabla \boldsymbol{\varsigma } -\nabla \tilde{\boldsymbol{\varsigma }}||_{0,\Omega }. \end{align}

Proof. This lemma is proved similar to Lemma 8.

Theorem 5.1. Suppose that the hypotheses of Theorems 3.6 and 3.7 hold and if the following assumption

(5.6) \begin{equation} \frac{2 \tilde{\beta } k_L\alpha _4c_s}{\alpha ^2_3Da}\lt 1 \end{equation}

holds then $T$ has a unique fixed point $ \mathbf{U}_s\in H^1_0(\Omega )^d$ , which in turn implies that coupled problem $(Q_{w_1}(\nabla \cdot \mathbf{U}_s))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$

Proof. If $ \frac{4 \tilde{\beta } k_L\alpha _4c_s}{\alpha ^2_3Da}\lt 1$ , then Lemma 12 implies $T\,{:}\, \mathbf{W}_r\subset H^1_0(\Omega )^d \rightarrow \mathbf{W}_r\subset H^1_0(\Omega )^d$ is a strict contraction mapping. That implies $T$ has a unique fixed point $ \mathbf{U}_s\in H^1_0(\Omega )^d$ due to Banach’s fixed-point theorem. That means the coupled problem $(Q_{w_1}(\nabla \cdot \mathbf{U}_s))-(Q_{w_2}(\nabla \cdot \mathbf{U}_s))$ or equivalently, the problem (3.3)–(3.5) with $\boldsymbol{\varsigma }=\nabla \cdot \mathbf{U}_s$ has a unique solution $(\mathbf{V}_f, P)\in H^1(\Omega )^d\times L^2(\Omega )$ and $ \mathbf{U}_s\in H^1_0(\Omega )^d.$

Remark 5.2.

• • If the non-dimensional parameters and the constants satisfy the following constraints

  (5.7) \begin{align} 2\alpha ^* Da\gt \frac{k_Lc_s\alpha _4}{\alpha _3}, \frac{4\alpha _1Da}{c_k}\gt \left (\frac{c_pk^2_2}{k_1} +\frac{3k_Lc_s\alpha _4}{\alpha _3}\right ), 2\alpha _2\geqslant \frac{\phi ^2_s}{a_0}. \end{align}
  then the continuous dependence (4.18) holds with the modified constant
  \begin{equation*}\alpha _6=\min \left \{\alpha ^*-\frac {k_Lc_s\alpha _4} {2\alpha _3Da},\left (\frac {2\alpha _1}{c_k} -\frac {c_pk^2_2}{2k_1Da} -\frac {3k_Lc_s\alpha _4}{2\alpha _3Da} \right ), \frac {a_0}{2}\right \}.\end{equation*}

6. Unbounded $\mathbf{K}$

One may note that theorems in Sections 3, 4 and 5 are proven under the boundedness assumption (i) of (3.1) and Lipschitz continuity (3.2) of $\mathbf{K}=\mathbf{K}(\boldsymbol{\varsigma }).$ In this section, we would like to relax such assumptions for case (a) $\boldsymbol{\varsigma }=\mathbf{U}_s$ and case (b) $\boldsymbol{\varsigma }=\nabla \cdot \mathbf{U}_s$ . Instead of boundedness property (ii) of (3.1), we assume there exists a constant $k_0\gt 0,$ such that following sub-linear growth condition holds

(6.1) \begin{eqnarray} {||\mathbf{K}(\mathbf{x})||\leqslant k_0(1+||\mathbf{x}||)\,\mbox{for all}\,\mathbf{x}\in \mathbb{R}^d}. \end{eqnarray}

The above growth condition yields for case (a)

(6.2) \begin{eqnarray} ||\mathbf{K}(\mathbf{U}_s)||_{0,\Omega }\leqslant \sqrt{2}k_0(\sqrt{|\Omega |}+||\mathbf{U}_s||_{0,\Omega }). \end{eqnarray}

We have the following result for case (a).

Theorem 6.1. Assume that the data and parameters satisfy the assumption (A), $(i)$ of (3.1) and (3.2), (3.37), (3.38) and $\mathbf{K}$ satisfy the growth condition (6.2). Further, if the following parameter constraint

(6.3) \begin{align} \frac{2\alpha _1}{c_k}\gt \frac{\sqrt{2c_p}\,c_sk_0\alpha _4}{\alpha _3Da} \end{align}

is satisfied, then (3.3)–(3.5) together with $\boldsymbol{\varsigma }=\mathbf{U}_s$ has a solution $(\mathbf{V}_f, \mathbf{U}_s,P)\in H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega )$ such that

(6.4) \begin{align} ||(\mathbf{V}_f,P)||_{\mathbb{Y}} \leqslant \frac{\alpha _4}{\alpha _3} \end{align}

and

(6.5) \begin{align} ||\nabla \mathbf{U}_s||_{0,\Omega }\leqslant \frac{1}{{\left (\frac{2\alpha _1}{c_k}-\frac{\sqrt{2c_p}\,c_sk_0\alpha _4}{\alpha _3Da}\right )}}\left [ \sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }+\frac{\alpha _4}{\alpha _3}\left (\phi _s+\frac{\sqrt{2|\Omega |}\,c_sk_0}{Da}\right )\right ] \end{align}

where $|\Omega |$ denotes the area or volume of the domain $\Omega .$

Proof. We are inspired by the working techniques from [28]. Here, we are dealing with particular non-linear structures of the hydraulic resistivity. We approximate the operators $\mathbf{K}$ with a sequence of uniformly positive definite, bounded operators, $\{\mathbf{K}^m\}_{m\geqslant 1}$ defined by

(6.6) \begin{eqnarray} \mathbf{K}^m_{ij}(\mathbf{U}_s)\,{:\!=}\,\min (m,\mathbf{K}_{ij}(\mathbf{U}_s)), \,\mbox{for all}\,m\in \mathbb{N}. \end{eqnarray}

Since $\mathbf{K}^m$ is bounded for each $m,$ Theorem 4.2 implies that there exists a triplet

$(\mathbf{V}^m_f, \mathbf{U}^m_s,P^m)$ in $H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega )$ such that

(6.7) \begin{align} \nonumber 2(D(\mathbf{V}^{m}_f)\,{:}\,D(\mathbf{W}))_{\Omega }+\lambda (\nabla \cdot \mathbf{V}^{m}_f,\nabla \cdot \mathbf{W})_{\Omega } -\phi _f(P^{m},\nabla \cdot \mathbf{W})_{\Omega } \\ \nonumber +\frac{1}{Da}(\mathbf{K}^m(\mathbf{U}^{m}_s)\mathbf{V}^{m}_f,\mathbf{W})_{\Omega } +\phi _f(\nabla \cdot \mathbf{V}^{m}_f,q)_{\Omega } +a_0(P^{m},q)_{\Omega }\\ =(\mathbf{b}_f,\mathbf{W})_{\Omega } +(\mathbf{T}_\infty,\mathbf{W})_{\partial \Omega }+(a_0,q)_{\Omega }, \end{align}

(6.8) \begin{align} \nonumber 2\alpha _1(D(\mathbf{U}^{m}_s)\,{:}\,D(\mathbf{Z}))_{\Omega } +\alpha _2(\nabla \cdot \mathbf{U}^{m}_s,\nabla \cdot \mathbf{Z})_{\Omega } = \phi _s(P^{m},\nabla \cdot \mathbf{Z})_{\Omega } \\ +\frac{1}{Da}(\mathbf{K}^m(\mathbf{U}^{m}_s)\mathbf{V}^{m}_f,\mathbf{Z})_{\Omega }+ (\mathbf{b}_s,\mathbf{Z})_{\Omega }. \end{align}

We desire a uniform bound for the sequence $\{(\mathbf{V}^m_f, \mathbf{U}^m_s,P^m)\}$ , which is independent of $m.$ In order to do so, replace $(\mathbf{W},q)$ by $(\mathbf{V}^m_f,P^m)$ in (6.7). Using the positive definiteness of $\mathbf{K}^m$ and the Cauchy-Schwarz, trace, Korn’s inequalities, we get the following estimate

(6.9) \begin{align} ||(\mathbf{V}^m_f,P^m)||_{\mathbb{Y}} \leqslant \frac{\alpha _4}{\alpha _3}. \end{align}

One can observe the right-hand side of (6.9) is independent of $m.$ Next, replacing $\mathbf{Z}$ by $\mathbf{U}^m_s$ in (6.8) and use Hölder’s, Poincare’s, Korn’s, Sobolev inequalities and (6.1), to get

\begin{align*} \frac{2\alpha _1}{c_k}||\nabla \mathbf{U}^m_s||^2_{0,\Omega } \leqslant \sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }||\nabla \mathbf{U}^m_s||_{0,\Omega }+\phi _s||P^m||_{0,\Omega }||\nabla \cdot \mathbf{U}^m_s||_{0,\Omega } \\ +\frac{1}{Da}||\mathbf{K}(\mathbf{U}^{m}_s)||_{0,\Omega }||\mathbf{V}^{m}_f||_{L^4(\Omega )}||\mathbf{U}^{m}_s||_{L^4(\Omega )} \end{align*}

\begin{align*} \frac{2\alpha _1}{c_k}||\nabla \mathbf{U}^m_s||^2_{0,\Omega } \leqslant \sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }||\nabla \mathbf{U}^m_s||_{0,\Omega }+\phi _s||P^m||_{0,\Omega }||\nabla \mathbf{U}^m_s||_{0,\Omega } \\ +\frac{c_s}{Da}[\sqrt{2}k_0(\sqrt{\Omega }+\sqrt{c_p}||\nabla \mathbf{U}^m_s||_{0,\Omega })]||\mathbf{V}^{m}_f||_{1,\Omega }||\nabla \mathbf{U}^{m}_s||_{0,\Omega } \end{align*}

(6.10) \begin{align} ||\nabla \mathbf{U}^m_s||_{0,\Omega }\leqslant \frac{1}{{\left (\frac{2\alpha _1}{c_k}-\frac{\sqrt{2c_p}\,c_sk_0\alpha _4}{\alpha _3Da}\right )}}\left [\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }+\frac{\alpha _4}{\alpha _3}\left (\phi _s+\frac{\sqrt{2|\Omega |}\,c_sk_0}{Da}\right )\right ] \end{align}

The estimates (6.9) and (6.10) imply that $(\mathbf{V}^m_f, \mathbf{U}^m_s,P^m)\in H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega )$ uniformly bounded for all $m\geqslant 1.$ Hence, there exists a triplet $(\mathbf{V}_f, \mathbf{U}_s,P)\in H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega )$ and a sub-sequence of $(\mathbf{V}^m_f, \mathbf{U}^m_s, P^m)$ (we denote by the same symbol) such that

(6.11) \begin{eqnarray} (\mathbf{V}^m_f, \mathbf{U}^m_s,P^m)\rightharpoonup (\mathbf{V}_f, \mathbf{U}_s,P)\,\mbox{weakly in}\,H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega ), \end{eqnarray}

and the compact embedding $H^1(\Omega )^d\hookrightarrow L^4(\Omega )^d$ yields

(6.12) \begin{eqnarray} (\mathbf{V}^m_f, \mathbf{U}^m_s)\rightarrow (\mathbf{V}_f, \mathbf{U}_s)\,\mbox{strongly in}\,L^4(\Omega )^d\times L^4(\Omega )^d. \end{eqnarray}

We seek to pass the limit in (6.7)–(6.8) as $m\rightarrow \infty .$ We observe that the weak convergence (6.11) is sufficient to pass the limit in the linear terms of (6.7)–(6.8); however, the non-linear terms demand strong convergence (6.12). Subsequently, the non-linear terms of the problem (6.7)–(6.8), which can be rewritten as

(i) for (6.7)

(6.13) \begin{align} \underbrace{((\mathbf{K}^m(\mathbf{U}^{m}_s)-\mathbf{K}(\mathbf{U}_s))\mathbf{V}^{m}_f,\mathbf{W})_{\Omega }} +\underbrace{(\mathbf{K}(\mathbf{U}_s)(\mathbf{V}^{m}_f-\mathbf{V}_f),\mathbf{W})_{\Omega }}+ (\mathbf{K}(\mathbf{U}_s)\mathbf{V}_f,\mathbf{W})_{\Omega } \end{align}

(ii) for (6.8)

(6.14) \begin{align} \underbrace{((\mathbf{K}^m(\mathbf{U}^m_s)-\mathbf{K}(\mathbf{U}_s))\mathbf{V}^{m}_f,\mathbf{Z})_{\Omega }} +\underbrace{(\mathbf{K}(\mathbf{U}_s)(\mathbf{V}^{m}_f-\mathbf{V}_f),\mathbf{Z})_{\Omega }} +(\mathbf{K}(\mathbf{U}_s)\mathbf{V}_f, \mathbf{Z})_{\Omega }, \end{align}

Since $\mathbf{U}^m_s$ converges to $\mathbf{U}_s$ strongly in $L^2(\Omega )^d,$ it implies $\mathbf{U}^{m}_s-\mathbf{U}_s\rightarrow 0$ a.e. in $\Omega$ up to a sub-sequence. The fact that $\mathbf{K}$ is Lipschitz guarantees that

(6.15) \begin{eqnarray} {\mathbf{K}^m(\mathbf{U}^{m}_s)-\mathbf{K}(\mathbf{U}_s)\rightarrow 0\,\mbox{a.e. in}\,\,\Omega .} \end{eqnarray}

Indeed, we have

\begin{align*}{|\mathbf{K}^m(\mathbf{U}^{m}_s)(\mathbf{x})-\mathbf{K}(\mathbf{U}_s)(\mathbf{x})|}\,{=}\,|\mathbf{K}^m(\mathbf{U}^{m}_s)(\mathbf{x})-\mathbf{K}^m(\mathbf{U}_s)(\mathbf{x})+ \mathbf{K}^m(\mathbf{U}_s)(\mathbf{x})+\mathbf{K}(\mathbf{U}_s)(\mathbf{x})|\\{\leqslant k_L|\mathbf{U}^{m}_s(\mathbf{x})-\mathbf{U}_s(\mathbf{x})| +|\mathbf{K}^m(\mathbf{U}_s)(\mathbf{x})-\mathbf{K}(\mathbf{U}_s)(\mathbf{x})|.} \end{align*}

As $\mathbf{U}^{m}_s-\mathbf{U}_s\rightarrow 0$ a.e. in $\Omega$ together with the definition of the truncation function implies $\mathbf{K}^m(\mathbf{U}_s)-\mathbf{K}(\mathbf{U}_s)\rightarrow 0$ a.e. in $\Omega .$ This establishes (6.15). Using the bound (6.9) and the convergence result (6.15), we get

(6.16) \begin{align} \lim _{m\rightarrow \infty }((\mathbf{K}^m(\mathbf{U}^{m}_s)-\mathbf{K}(\mathbf{U}_s)) \mathbf{V}^{m}_f,\mathbf{W})_{\Omega }=0,\,\,\lim _{m\rightarrow \infty }((\mathbf{K}^m(\mathbf{U}^{m}_s)-\mathbf{K}(\mathbf{U}_s)) \mathbf{V}^{m}_f,\mathbf{Z})_{\Omega }=0. \end{align}

The bound (6.10) together with estimate (6.2), and the strong convergence $\mathbf{V}^m_f\rightarrow \mathbf{V}_f$ in $L^4(\Omega )^d$ leads to

(6.17) \begin{align} \lim _{m\rightarrow \infty }(\mathbf{K}(\mathbf{U}_s) (\mathbf{V}^{m}_f-\mathbf{V}_f),\mathbf{W})_{\Omega }=0, \,\, \lim _{m\rightarrow \infty }(\mathbf{K}(\mathbf{U}_s) (\mathbf{V}^{m}_f-\mathbf{V}_f),\mathbf{Z})_{\Omega }=0 \end{align}

Thus, the terms with under braces in (6.13) and (6.14) tend to zero as $m\rightarrow \infty .$ Hence, we can say that (6.7)–(6.8) recovers (3.3)–(3.5) as $m\rightarrow \infty .$ This completes the proof of the Theorem 6.1.

Remark 6.2.

• • We note that in the above case i.e. when $\mathbf{K}(\mathbf{U}_s)$ is not bounded, however, satisfies the sub-linear growth condition (6.1), the uniqueness of solutions holds whenever the non-dimensional parameters and constants satisfy the following inequalities:

  (6.18) \begin{align} \frac{2\alpha Da}{c_s}\gt \left [\frac{k_L\sqrt{c_p}\alpha _4}{\alpha _3}+{\sqrt{2}k_0} (\sqrt{|\Omega |}+\sqrt{c_p}\alpha _7)\right ],\,\,2\alpha _2\geqslant \frac{\phi ^2_s}{a_0}, \end{align}
  (6.19) \begin{align} \frac{4\alpha _1Da}{c_kc_s}\gt \left [\frac{3k_L\sqrt{c_p}\alpha _4}{\alpha _3} +{\sqrt{2}k_0}(\sqrt{|\Omega |}+\sqrt{c_p}\alpha _7)\right ], \end{align}

  where $\alpha _7= \frac{1}{{\left (\frac{2\alpha _1}{c_k}-\frac{\sqrt{2c_p}\,c_sk_0\alpha _4}{\alpha _3Da}\right )}}\left [\sqrt{c_p}||\mathbf{b}_s||_{0,\Omega }+\frac{\alpha _4}{\alpha _3}\left (\phi _s+\frac{\sqrt{2|\Omega |}\,c_sk_0}{Da}\right )\right ].$

• • In case (b) that is when $\boldsymbol{\varsigma }=\nabla \cdot \mathbf{U}_s$ , the sub-linear growth condition (6.1) becomes

  (6.20) \begin{equation} ||\mathbf{K}(\nabla \cdot \mathbf{U}_s)||_{0,\Omega }\leqslant \sqrt{2}k_0(\sqrt{|\Omega |}+||\nabla \cdot \mathbf{U}_s||_{0,\Omega }). \end{equation}
  The specific structure $\mathbf{K}(\nabla \cdot \mathbf{U}_s) = [\gamma _1 + \gamma _2 |\nabla \cdot \mathbf{U}_s|]\mathbf{I}$ falls under this case, and we can clearly see that $\mathbf{K}$ satisfies (6.20). In this case, existence and uniqueness analysis can be developed based on similar lines as in Theorem 6.1. Indeed, we have the following theorem.

Theorem 6.3. Assume that the data and parameters satisfy the assumption (A), $(i)$ of (3.1) and (3.2), (3.44), (3.45) and $\mathbf{K}$ satisfy the growth condition (6.20). Further, if the following parameter constraint

(6.21) \begin{align} \frac{2\alpha _1}{c_k}\gt \frac{\sqrt{2}\,c_sk_0\alpha _4}{\alpha _3Da} \end{align}

is satisfied, then (3.3)–(3.5) together with $\boldsymbol{\varsigma }=\nabla \cdot \mathbf{U}_s$ has a solution $(\mathbf{V}_f, \mathbf{U}_s,P)\in H^1(\Omega )^d \times H^1_0(\Omega )^d \times L^2(\Omega )$ such that

(6.22) \begin{align} ||(\mathbf{V}_f,P)||_{\mathbb{Y}} \leqslant \frac{\alpha _4}{\alpha _3},\,\,\,\, ||\nabla \mathbf{U}_s||_{0,\Omega }\leqslant \alpha _8 \end{align}

where $|\Omega |$ denotes the area or volume of the domain $\Omega$ and

\begin{equation*}\alpha _8=\frac {1}{{\left (\frac {2\alpha _1}{c_k}-\frac {\sqrt {2}\,c_sk_0\alpha _4}{\alpha _3Da}\right )}}\left [ \sqrt {c_p}||\mathbf {b}_s||_{0,\Omega }+\frac {\alpha _4}{\alpha _3}\left (\phi _s+\frac {\sqrt {2|\Omega |}\,c_sk_0}{Da}\right )\right ].\end{equation*}

Further, the solution is unique and subject to the following constraints

(6.23) \begin{align} \frac{2\alpha Da}{c_s}\gt \left [\frac{k_L\alpha _4}{\alpha _3}+{\sqrt{2}k_0} (\sqrt{|\Omega |}+\alpha _8)\right ],\,\,2\alpha _2\geqslant \frac{\phi ^2_s}{a_0}, \end{align}

(6.24) \begin{align} \frac{4\alpha _1Da}{c_kc_s}\gt \left [\frac{3k_L\alpha _4}{\alpha _3} +{\sqrt{2}k_0}(\sqrt{|\Omega |}+\alpha _8)\right ]. \end{align}

Proof. The proof of this theorem is similar to Theorem 6.1, we omit the details.

6.1 Comments on parameter restrictions

It may be noted that the existence and uniqueness of results (e.g. see Theorems 3.3, 3.5, 3.6, 3.7, 4.2, 5.1, 6.1 and 6.3, etc.) that are established in this work hold subject to certain parameter restrictions see e.g. (3.37), (3.38), (3.44), (3.45), (4.16), (4.17), (5.6), (5.7), (6.3), (6.18), (6.19), (6.21), (6.23) and (6.24) etc. Such a situation is typical in the case of multiphase mixture models where extra care needs to be paid to the physical admissibility of the parameters. See also Vromans et al. [29]. Further, some of the parameters in these inequalities do show a dependency on material properties of the tissue, like Poisson ratio, Young’s modulus etc. Thus, one has to ensure that these inequalities are satisfied simultaneously. This certainly depends on the choice of relevant parameters. We have collected data from the existing literature on parameters like $Da$ , $\alpha _t,$ $\varrho _t$ etc. which are relevant to biological tissues. We have then verified that there do exist parameter combinations within the given ranges in Table 1 that obey all the inequalities. This ensures that these assumptions (3.37), (3.38), (3.44), (3.45), (4.16), (4.17), (5.6), (5.7), (6.3), (6.18), (6.19), (6.21), (6.23) and (6.24) etc. can be interpreted from the point of view of various applications. We list the dimensionless poro-elasto-hydrodynamics parameters in Table 1. Further, we choose $\Omega$ as a d-dimensional ( $d=3$ ) sphere of unit radius (in dimension 1 mm) then $|\Omega |=\frac{4\pi }{3},$ $|\partial \Omega |=4\pi .$ We set $\mathbf{b}_i=0,$ $i\in \{f,s\},$ $\mathbf{T}_\infty =(1,0,0),$ and choose the numerical values for constants $c_k(\Omega )\gt 0,\,c_p(\Omega )\gt 0,\,c_s(d)\gt 0,\,c_t(\Omega,d)\gt 0$ as follows: $c_k=2 \,(or\,3)$ for $H^1_0(\Omega )^d$ (or $H^1(\Omega )^d$ ) [33, 34] and $c_p=1/2,$ $c_s=1/2$ [35], $c_t=2$ [34]. To verify the inequalities (3.37), (3.38), (3.44), (3.45), (4.16), (4.17), (5.6), (5.7), (6.3) and (6.21), we use the following combination of parameters from Table 1: $Da=2\times 10^{-2},$ $\alpha _t=1,$ $\varrho _t=10^4,$ $\nu _p=0.45,$ $k_1=0.5,$ $k_2=1.4,$ $\phi _s=0.4,$ $L_rA_r=1,$ $k_L=2\times 10^{-3},$ and $0\lt k_0\leqslant 1.$ Further, the restrictions (6.18), (6.19) and (6.23), (6.24), which ensure the uniqueness for arbitrary $\mathbf{K}$ , hold when we choose $Da=10^{-1},$ and $k_0=10^{-2}$ with the remaining parameters as chosen above. Note that the choice of parameters mentioned above may not be unique; there could be other parameter combinations for which these restrictions hold true.

Table 1. Dimensionless poro-elasto-hydrodynamics parameters corresponding to tumour tissue with their value range.

7. Summary

In this work, we have modelled the poro-elasto-hydrodynamics that mimic an in-vitro solid tumour using biphasic mixture theory. We simplified the generic biphasic mixture equations using certain assumptions based on the biological context and treated hydraulic resistivity as anisotropic, which depends on the deformation. This made our model non-linear and coupled. We derive an equivalent variational (or weak) formulation and developed existence and uniqueness results using the Galerkin method, monotone operator theory and fixed-point theory. The detailed analysis is done by considering two cases: (a) $\mathbf{K}(\boldsymbol{\varsigma })=\mathbf{K}(\mathbf{U}_s)$ (b) $\mathbf{K}(\nabla \cdot \boldsymbol{\varsigma })=\mathbf{K}(\nabla \cdot \mathbf{U}_s)$ . In both cases, we first developed existence and uniqueness analysis for auxiliary linear and semilinear sub-problems using the Galerkin method and monotone operator theory. Then, we convert the corresponding coupled non-linear problem to the fixed-point problem in both cases. Further, to develop the existence of solutions for the corresponding fixed-point problems, we used the Schauder fixed-point theorem. Uniqueness is assured via the Banach contraction theorem. For the case where $\mathbf{K}$ is not bounded but satisfies the sub-linear growth condition, we have developed the existence and uniqueness results. Moreover, we have collected certain realistic ranges of parameters involved in the model and ensured that the theoretical restrictions derived by us are compatible with these parameter ranges.