\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 110, pp. 1--15.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2013/110\hfil 
 Existence and uniqueness of stationary solutions]
{Existence and uniqueness of stationary solutions
 to bioconvective flow equations}

\author[J. L. Boldrini, M. A. Rojas-Medar, M. D. Rojas-Medar 
\hfil EJDE-2013/110\hfilneg]
{Jos\'e Luiz Boldrini, Marko Antonio Rojas-Medar,\\
 Maria Drina Rojas-Medar}  % in alphabetical order

\address{Jos\'e Luiz Boldrini \newline
Unicamp-Imecc, Rua S\'ergio Buarque de Holanda,
651: 13083-859 Campinas, SP, Brazil}
\email{josephbold@gmail.com}

 \address{Marko Antonio Rojas-Medar \newline
Dpto. de Ciencias B\'asicas, Facultad de Ciencias,
Universidad del B\'{\i}o-B\'{\i}o, Campus Fernando May,
Casilla 447, Chill\'an, Chile}
\email{marko@ubiobio.cl}

\address{Maria Drina Rojas-Medar \newline
Dpto. de Matem\'aticas, Facultad de Ciencias B\'asicas, 
Universidad de Antofagasta, Casilla 170, Antofagasta, Chile}
\email{mrojas@uantof.cl}

\thanks{Submitted January 17, 2013. Published April 29, 2013.}
\subjclass[2000]{35Q80, 76Z10, 92B05}
\keywords{Bioconvective flow; stationary solutions}

\begin{abstract}
 We analyze a system of nonlinear partial  differential equations
 modeling the stationary flow induced by the upward swimming of certain
 microorganisms in a fluid. We consider the realistic  case in which
 the effective viscosity of the fluid depends on the concentration of
 such microorganisms. Under certain conditions, we prove the existence
 and uniqueness  of solutions for such generalized bioconvective flow equations
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks

\section{Introduction}\label{Intro}

In this work we perform a mathematical analysis of the 
system of partial differential equations
\begin{equation}\label{sistema1}
\begin{gathered}
-2 \operatorname{div} (\nu(m)D(\mathbf{u})) + \mathbf{u}\cdot \nabla \mathbf{u} + \nabla q
= -m \cdot \chi + \mathbf{f} \quad \text{in } \Omega ,\\
 \operatorname{div} \mathbf{u} = 0 \quad \text{in } \Omega ,\\
-\theta \Delta m + \mathbf{u}\cdot \nabla m + U \frac{\partial m}{\partial x_3} = 0 \quad 
\text{in } \Omega ,
\end{gathered}
\end{equation}
subject to the boundary and total amount conditions
\begin{equation} \label{sistema2}
\begin{gathered}
\mathbf{u} = \mathbf{0} \quad \text{on } S,\\
\mathbf{u}\cdot\mathbf{n} = 0 \quad \text{on } \Gamma,\\
\nu(m) [D(\mathbf{u})\mathbf{n} - \mathbf{n}\cdot (D(\mathbf{u})\mathbf{n}) \mathbf{n}] = \mathbf{b}_1 \quad \text{on } \Gamma,\\
\theta \frac{\partial m}{\partial \mathbf{n}}- Un_3 m = 0 \quad 
\text{on } \partial \Omega,\\
\int_\Omega m dx = \alpha .
\end{gathered}
\end{equation}

This system is a mathematical model for the stationary flow induced by the 
upward swimming of certain microorganisms in a fluid
in the realistic situation that the concentration of such microorganisms
 may affect the effective viscosity of the fluid.
The last condition in \eqref{sistema2} fix the total mass of microorganism
as  $\alpha$, which is a given positive constant.

The flow occurs in a set $\Omega \subset \mathbb{R}^3$, which 
is assumed to be a bounded domain with smooth boundary $\partial \Omega $;
 $S$ and $\Gamma$ are disjoint open subsets of $\partial \Omega$ such 
that $\partial\Omega = S \cup \overline{\Gamma}$, and the superficial 
measure of $S$ is strictly positive.

The unknowns in the problem are the fluid velocity $\mathbf{u}$, its associated 
pressure $p$ and the concentration of microorganisms $m$.
In the previous equations, the following are given:
the fluid viscosity function $\nu(\cdot) > 0$;
the constant rate of diffusion of microorganisms $\theta >0$;
the constant vector $\chi = (0, 0, 1)^t$, meaning that the coordinate 
system is placed such that the gravitational force acts along the vertical.
Other given data are the following: $\mathbf{f}$, the external force field;
$U$, the average upward speed of swimming of the microorganisms;
$\alpha$, a positive constant given total mass of microorganisms.

As usual, the symbols $\nabla, \Delta$ and $\operatorname{div}$ denote respectively
the gradient, Laplacian and divergence operators;
$\mathbf{u}\cdot\nabla \mathbf{u}$ denotes the convection operator, whose component are
$i$-th  in cartesian coordinates is given by $(\mathbf{u}\cdot \nabla \mathbf{u})_i
= \sum^3_{j=1} u_j u_{i, x_j}$.
Also, $D(\mathbf{u}) = (\nabla \mathbf{u} + (\nabla \mathbf{u})^t)/2$ is the symmetric part 
of the deformation  rate tensor.


Bioconvective flows have been studied by many autors along the years; 
here we just mention the book by
Levandowsky, Childress, Hunter and Spiegel \cite{levandowsky}, 
and the articles by Moribe \cite{moribe} and
Kan-On, Narukawa and Teramoto \cite{kan}; the interested reader can 
consult also the references mentioned in these works.
Next, we briefly comment on previously published articles that are 
directly related to the present one.

Kan-on, Narukawa and Teramoto in \cite{kan} analyzed the 
classical bioconvective equations  (i.e., the case of constant fluid viscosity)
with Dirichlet boundary conditions for the fluid velocity. 
By using  fixed point arguments,  they proved the existence of generalized 
solutions; with the help of  classical regularity results, they also prove 
the existence of strong solutions in the case that of small enough upward 
swimming speed $U$.


On the other hand, Lorca and Boldrini in \cite{lorca} obtained  results 
on existence and uniqueness of weak solutions for the generalized Boussinesq 
equations, which are equations governing thermally driven flows in the case 
that the viscosity and thermal conductivity coefficients may depend on the 
temperature. They also considered Dirichlet boundary conditions for the 
velocity and used Galerkin approximations and fixed point arguments, 
together with estimates for the ''pressure" associated to a Helmholtz 
decomposition of a $L^2$-field, to prove their results.

The present work generalizes the results of Kan-on, Narukawa and Teramoto 
in \cite{kan} to  the generalized bioconvective equations,
that is,  to the case where the fluid viscosity may depend on the 
concentration of microorganisms.
As in \cite{kan}, we prove results on existence of weak and strong solutions 
of problem when $U$ is small; we also give a result on uniqueness of weak 
solutions.
For the proofs, we have to adapt some of the techniques presented in Lorca 
and Boldrini \cite{lorca} to the case of our boundary conditions, including 
the  estimates for the pressure associated to a Helmholtz decomposition.
We remark that one major difficulty to obtain more regular solutions of 
our system of equations,
as well as in the case of the generalized Boussinesq equations,
is to estimate the nonlinear terms; we are able to do this with the help 
of the previously mentioned estimates for the pressure associated to 
a Helmholtz decomposition.

We finally observe that we are also analyzing the associated evolution 
problem; the results of such analysis will appear elsewhere;
we also remark that Climent-Ezquerra {\it et al.} recently proved in  
\cite{Climent}  the existence of time reproductive solutions for this 
same evolution problem.

\section{Preliminaries}

Let  $\Omega \subseteq \mathbb{R}^{3}$ be a bounded domain of class $C^{3}$.
We consider the Lebesgue spaces  $L^p(\Omega)$, $1 \leq p \leq \infty$, 
with the usual norms $|u|_p$;
for simplicity, we just denote $|\cdot|_2 = |\cdot|$ and the usual inner 
product in $L^2(\Omega)$ by $(\cdot, \cdot)$.
For $m \geq 0$ and $1 \leq p < \infty$, we consider the usual Sobolev
spaces
$W^{m,p} (\Omega) = \{u \in L^p(\Omega); D^\alpha u \in L^p(\Omega), 
 \forall |\alpha| \leq m\}$,
with the norm $\|u\|_{m,p} = [ \sum_{|\alpha|\leq
m} \|D^\alpha u(t)\|^p_{L^p(\Omega)}]^{1/p}$;
when $p=2$, we denote, as usual, $W^{m,p} (\Omega) = H^m(\Omega)$.
Also, $W^{1-\frac1p, p} (\partial\Omega)$ is the space of traces
on $\partial\Omega$ of functions in $W^{1,p}(\Omega)$, equipped with the norm
$|\gamma|_{W^{1-\frac1p, p_{(\partial\Omega)}}}
= \inf \{\|v\|_{W^{k,p}(\Omega)}; v \in W^{1,p} (\Omega), v = \gamma 
\text{ on } \partial\Omega\} $;
when $p = 2$, we denote $W^{1/2,2} (\partial\Omega) = H^{1/2} (\partial\Omega)$. 
For details and properties of such spaces, see Adams \cite{adams}.


We will also need the following classical results.

\begin{lemma}[Poincar\'e-Friedrichs inequality] \label{Lema1.1}
Let $\Sigma \subseteq \partial\Omega$ a portion of boundary with strictly
 positive superficial measure;
then there exists a positive constant $C_P$ depending only on $\Omega$ 
and $\Sigma$ such that $|u| \leq C_P |\nabla u|$, for all $u \in H^1(\Omega)$ 
such that $u |_\Sigma = 0$.
\end{lemma}

\begin{lemma} \label{Lema1.2} 
There exists a constant $C_\Omega$, depending only on $\Omega$, such that
$|\phi| \leq C_\Omega |\nabla\phi|$, for all $\phi \in B = H^1(\Omega)\cap Y$.
\end{lemma}


To treat the unknown velocity, we will need the following functional spaces:
being $S$ and $\Gamma$ as described in the Introduction, we define the 
following functional spaces
\begin{gather*}
\dot {\mathbf{H}}(\Omega) = \{\mathbf{u} \in (C^\infty (\overline\Omega))^3:
\mathbf{u}|_S = 0,\; \mathbf{u}\cdot \mathbf{n}|_{\Gamma} = 0\}, \\
\mathbf{H}(\Omega) \text{ closure of $\dot{\mathbf{H}}(\Omega)$  with respect
to norm }\|\cdot\|_{\mathbf{H}(\Omega)},
\end{gather*}
where
\begin{equation}\label{eq1}
\|\mathbf{u}\|_{\mathbf{H}(\Omega)} = \Big[\int_\Omega \nabla \mathbf{u} :
\nabla \mathbf{u} dx\Big]^{1/2} = [(\nabla \mathbf{u}, \nabla \mathbf{u})]^{1/2} = |\nabla \mathbf{u}| .
\end{equation}


Given $\mathbf{u}:\Omega \to \mathbb{R}^3$ with suitable regularity, the rate 
of strain tensor is defined as
$D(\mathbf{u}) = \frac12 (\nabla \mathbf{u} + (\nabla\mathbf{u})^t)$.
Let us also define
\[
(D(\mathbf{u}), D(\mathbf{v})) = \int_\Omega D(\mathbf{u}) :
D(\mathbf{v}) dx = \int_{\Omega}\sum_{i, j=1}^{3}
(\frac{\partial u_{i}}{\partial x_{j}} +
\frac{\partial u_{j}}{\partial x_{i}})
(\frac{\partial v_{i}}{\partial x_{j}} +
\frac{\partial v_{j}}{\partial x_{i}})dx,
\]
and thus, $(D(\mathbf{u}), D(\mathbf{u})) = \int_\Omega D(\mathbf{u}): D(\mathbf{u}) dx \equiv |D(\mathbf{u})|^2$.

We will also need the following results:

\begin{lemma}[{Korn inequality \cite[p. 191]{solon}}] \label{Lema1.4} 
There exists a positive constant $\overline c$, such that
$$
\| \mathbf{u}\|_{\mathbf{H}(\Omega)} = |\nabla \mathbf{u}| \leq \overline c
|D(\mathbf{u})|, \quad \forall \mathbf{u} \in \mathbf{H}(\Omega).
$$
\end{lemma}

As a consequence of this lemma, we have the following result.

\begin{lemma} \label{Lema1.5} 
There exists a positive constant $\gamma$ such that
$|\mathbf{u}|^2 \leq \gamma |D(\mathbf{u})|^2$, for all $\mathbf{u} \in \mathbf{H}(\Omega) $.
\end{lemma}

The previous results imply that  the norms $|\nabla \mathbf{u}|$ and $|D(\mathbf{u})|$ 
are equivalent  in $\mathbf{H}(\Omega)$.
Next,  we consider
$\mathbf{C}^\infty_{0,\sigma} (\Omega)= \{ \mathbf{f} \in (C^\infty_0
(\Omega))^3 : \ \operatorname{div}\mathbf{f} = 0 \}$,
and then we take
$$
\mathbf{X}(\Omega)= \text{closure of $\mathbf{C}^\infty_{0,\sigma} (\Omega)$ 
in }  (L^2(\Omega))^3.
$$
It is well known \cite{temam} that
$$
(L^2(\Omega))^3 = \mathbf{X}(\Omega) \oplus \mathbf{G}(\Omega),
$$
with
$\mathbf{G}(\Omega) = \{\varphi \in (L^2(\Omega))^3, 
\varphi = \nabla q, \; q \in H^1 (\Omega)\}$.


We will also need  the following functional spaces:
\begin{gather*}
\dot{\mathbf{J}}(\Omega) = \{\mathbf{u} \in \dot{\mathbf{H}}(\Omega),  
\operatorname{div}
\mathbf{u} = 0\}, \\
\mathbf{J}_0(\Omega)=  \text{ closure of $\dot{\mathbf{J}}(\Omega)$ in the norm  \eqref{eq1}}.
\end{gather*}
Next, we consider the following applications:
\begin{gather*}
B_0 : \mathbf{J}_0(\Omega) \times \mathbf{J}_0(\Omega) 
\times \mathbf{J}_0(\Omega) \to \mathbb{R}, \\
B_1 : \mathbf{J}_0(\Omega) \times (H^1(\Omega))^3 \times (H^1(\Omega))^3 \to
\mathbb{R}
\end{gather*}
given by
\begin{equation} \label{e2}
\begin{gathered}
B_0(\mathbf{u}, \mathbf{v}, \mathbf{w}) = (\mathbf{u}\cdot\nabla \mathbf{v}, \mathbf{w}) 
= \int_{\Omega}\sum_{i,j=1}^{N}u_{j}(x)((\partial v_{i})
/(\partial x_j))(x)w_{i}(x)dx ,
\\
B_1 (\mathbf{u}, \phi, \psi) = (\mathbf{u}\cdot \nabla \phi,
\psi) = \int_{\Omega}\sum_{j=1}^{N}u_{j}(x)((\partial \phi)
/(\partial x_j))(x)\psi(x)dx.
\end{gathered}
\end{equation}
They are well defined trilinear forms with the following properties:
\begin{equation} \label{e3}
\begin{gathered}
B_0 (\mathbf{u}, \mathbf{v}, \mathbf{w}) = -B_0 (\mathbf{u}, \mathbf{w}, \mathbf{v}), \quad 
B_1 (\mathbf{u}, \phi, \psi) = -B_1 (\mathbf{u}, \psi, \phi),\\
B_0 (\mathbf{u}, \mathbf{v}, \mathbf{v}) = 0, \quad B_1 (\mathbf{u}, \phi, \phi) = 0.
\end{gathered}
\end{equation}

Let $P$ be the orthogonal projection from $(L^2 (\Omega))^3$ onto 
$\mathbf{X}(\Omega)$. Then, we define the operator $A$  as the 
Friedrichs extension of the symmetric operator $P \Delta$,
with  $D(A) = \{\mathbf{u} \in \mathbf{J}_0 (\Omega) \cap (H^2 (\Omega))^3; D(\mathbf{u})
\mathbf{n} - \mathbf{n}\cdot(D(\mathbf{u})\mathbf{n})\mathbf{n}|_{\Gamma}=0\}$.

The proofs of the following results concerning this operator $A$ 
can be found in Rionero and Mulone \cite[pp. 478-481]{rionero}.

\begin{lemma} \label{Lema1.6}
The Stokes operator $A: \mathbf{X}(\Omega) \to \mathbf{X} (\Omega)$ 
defined as the Friedrichs extension of $-P \Delta$,  with domain
$D(A) = \{ \mathbf{u} \in \mathbf{J}_0 (\Omega) \cap (H^2 (\Omega))^3; 
D(\mathbf{u})\mathbf{n} - \mathbf{n} \cdot (D(\mathbf{u})\mathbf{n})\mathbf{n}|_{\Gamma}=0\}$, is a selfadjoint,
positive definite operator with compact inverse.
\end{lemma}

Thus (see Brezis \cite{brezis}), $A$ has a sequence $\{\alpha_i\}_{i=1}^\infty$
of eigenvalues satisfying $0 < \alpha_1 \leq \alpha_2 \leq  \dots$ and 
$\lim_{i \to +\infty} \alpha_i = + \infty$,
whose associated eigenfunctions  $\{\overline{\mathbf{w}}^i\}_{i=1}^\infty$ 
form a complete orthogonal system  in
$\mathbf{X}(\Omega)$, $\mathbf{J}_0(\Omega)$ and $D(A)$,
with their natural inner products.

The following result concerning the Helmholtz decomposition is 
analogous to the one in  Lemma 3.4   of Lorca and Boldrini \cite{lorca1};
 its proof can be done similarly as in \cite{lorca1}.

\begin{lemma} \label{Lema1.8}
Let $\mathbf{v} \in J_0(\Omega) \cap (H^2(\Omega))^3$ and consider the Helmholtz 
decomposition of $-\Delta \mathbf{v}$:
$$
-\Delta \mathbf{v} = A\mathbf{v} + \nabla \overline q ,
$$
where $\overline q \in H^1(\Omega)$ and $\int_\Omega
\overline q dx = 0$. Then, there exists a positive constants $C > 0$ and, 
for any $\varepsilon > 0$, a associated positive constant $C_\varepsilon$,
such that
$$
 \|\overline q\|_1 \leq C |A\mathbf{v}|
\quad \text{and} \quad
|\overline q| \leq C_\varepsilon |\nabla \mathbf{v}| + \varepsilon |A\mathbf{v}|, \quad 
\forall \mathbf{v} \in D(A) .
$$
\end{lemma}

To treat the unknown microorganism concentration, we will need the following
functional spaces:

$Y$ is the closed subspace of $L^2(\Omega)$ consisting of functions that 
are orthogonal to the constants; i.e.,
$$
Y=\{f \in L^2(\Omega) : \int_\Omega f(x)dx =0 \}.
$$
We then define
$$
B=H^1 (\Omega) \cap Y.
$$

Next, let $\overline P$ be the orthogonal projection from $L^2(\Omega)$ onto $Y$.
As before,  an operator $A_1$ can be defined as the Friedrichs extension 
of the symmetric operator $\overline P (-\theta \Delta )$,
with domain
$D(A_1) = \{\varphi \in Y \cap H^2(\Omega); \theta
\frac{\partial\varphi}{\partial \mathbf{n}} - Un_3
\varphi = 0$ on $\partial\Omega\}$.
The proofs of the following results, which are similar to the ones for 
the operator $A$, can be found in
Kan-On, Narukawa and Teramoto \cite[pp. 150-152]{kan}.

\begin{lemma}\label{Lema1.7}
The operator $A_1 : Y \to Y $, defined as the Friedrichs extension of
$ \overline P (- \theta \Delta )$, with domain
$D(A_1) = \{\varphi \in Y \cap H^2(\Omega); \theta
 \frac{\partial\varphi}{\partial \mathbf{n}} - Un_3 \varphi = 0$ on $\partial\Omega\}$,
is a selfadjoint, positive definite operator with compact inverse.
\end{lemma}

From the definition of $A_1$, it  follows that $D(A_1^{1/2} ) = B$ and
$(\theta - 2 U C_P)^{1/2} |\nabla \varphi|  
\leq |A_1^{1/2} \varphi| \leq (\theta + 2 U C_P)^{1/2} |\nabla \varphi|   $ 
for all $\varphi \in B$, for this, see again 
Kan-On, Narukawa and Teramoto \cite[pp. 145]{kan}, 

Operator $A_1$ has a sequence $\{\beta_i\}_{i=1}^\infty$ of
eigenvalues satisfying
$0 < \beta_1 \leq \beta_2 \leq  \dots$ and
$\lim_{i \to +\infty} \beta_i = + \infty$,
and whose corresponding eigenfunctions $\{\overline \phi^i\}$ form a 
complete orthogonal system  in $Y$, $B$ and $D(A_1)$, with their natural 
inner products; we assume that it is normalized in $Y$.

We will also use the following  orthogonal projections: for each $n$, define
\begin{equation} \label{ProjectionBar-n}
\begin{gathered}
\overline P_n : Y   \to M_n = \operatorname{span}
\{ \overline \phi^1, \overline \phi^2, \ldots,  \overline \phi^n \},
\\
 f  \to \overline P_n(f) = \sum_{\ell=1}^n (f, \overline \phi^\ell)
 \overline \phi^\ell ,
\end{gathered}
\end{equation}

Standard computations with the fractional powers of $A_1$, using the 
previous results, give us the following:
\begin{gather}\label{Projection-prop1}
|\nabla  \overline P_n \varphi | \leq  |A_1^{1/2} \varphi |
 \leq  (\theta + 2 U C_P)^{1/2} | \nabla \varphi |, \quad \forall \varphi \in B ,
\\
\label{Projection-prop2}
\begin{aligned}
(\theta - 2 U C_P)^{1/2} |\nabla (\varphi - \overline P_n \varphi) |
&\leq |A_ 1^{1/2}  (\varphi - \overline P_n \varphi) |
\leq \frac{1}{\beta_{n+1}^{1/2}} |A_1  \varphi  |\\
&\leq \frac{\theta}{\beta_{n+1}^{1/2}} | \Delta \varphi |, \quad \forall 
\varphi \in D(A_1) .
\end{aligned}
\end{gather}

We will need an inequality similar to the last one, but holding for a 
larger set of functions.
For this, we consider the following extension of the operator $A_1$:
let $\tilde A_1$ be the Friedrichs extension of the symmetric operator
$\overline P (-\theta \Delta )$, but now with domain
$D(A_1) = \{\varphi \in H^2(\Omega); \theta 
\frac{\partial\varphi}{\partial \mathbf{n}} - Un_3
\varphi = 0$ on $\partial\Omega\}$.
By observing that the subspace $\operatorname{span} \{ 1 \}$ is the 
orthogonal complement of $Y$ in $L^2 (\Omega)$,
and in particular, $L^2 (\Omega) = \operatorname{span} \{ 1 \}  \oplus Y $, 
we have that $\tilde A_ 1 = 0 \oplus A_ 1$;
we conclude that $\tilde A_ 1$ is a semipositive selfadjoint operator with
eigenvalues $\{\beta_i\}_{i=0}^\infty$,
where $\beta_0 =0$ and its associated eigenfunction is the constant 
function $1$; the other eigenvalues are exactly the same ones of $A_1$, 
with the same eigenfunctions.
As before, we have that
$ |\tilde A_ 1^{1/2}  (\varphi - \tilde P_n \varphi) | 
\leq \frac{1}{\beta_{n+1}^{1/2}} |\tilde A_1  \varphi  | $
for all $\varphi \in D(\tilde A_ 1)$, where $\tilde P_ n$ is now the 
$L^2(\Omega)$-orthogonal projection
on $\operatorname{span} \{ 1 \} \oplus M_n$.
For $n \geq 1$, by using the definitions of fractional powers in 
terms of the eigenvalues and eigenfunctions,
we see that 
$\tilde A_ 1^{1/2}  (\varphi - \tilde P_n \varphi) 
=  A_ 1^{1/2}  (\varphi - P_n \varphi) $.
Finally, from the previous results, by proceeding as before, 
we finally also obtain
\begin{equation} \label{Projection-prop3}
(\theta - 2 U C_P)^{1/2} |\nabla (\varphi - \overline P_n \varphi) |
\leq
 \frac{\theta}{\beta_{n+1}^{1/2}} | \Delta \varphi |, \quad 
\forall \varphi \in D(\tilde A_1)
\end{equation}

\section{Existence of weak solutions}

To analyze our problem, it is convenient to introduce the following 
change of variables
$$
\overline m = m - E,
$$
where
$$
E(x) = C_\alpha \exp \Big(\frac U\theta x_3 \Big),
$$
and the constant $C_\alpha$ is chosen such that
$\int_\Omega E(x) dx = \alpha$.

The idea behind this change of variable is the following: since
$ -\theta \Delta E + U \frac{\partial E}{\partial x_3} = 0$
and
$  \theta \frac{\partial E}{\partial \mathbf{n}} - Un_3 E = 0$, we have that
$E(\cdot)$ is a particular solution of equation \eqref{sistema1} (iii) in the special case of no fluid motion, that is, when $\mathbf{u} \equiv 0$;
thus, our intension to look for solutions in a neighborhood of  this special microorganism distribution.

By writing Problem \eqref{sistema1} - \eqref{sistema2} in terms of variables
 $\mathbf{u}$ and $\overline m$, we obtain
\begin{equation}  \label{nuevasvariables}
\begin{gathered}
- 2\operatorname{div}(\nu(\overline m + E)D(\mathbf{u})) + \mathbf{u}\cdot \nabla \mathbf{u} +
\nabla (q+\frac{\theta}U E)
= \overline m \cdot \chi + \mathbf{f} ,
 \\
\operatorname{div} \mathbf{u} = 0 ,
 \\
- \theta \Delta \overline m
+ \mathbf{u}\cdot \nabla(\overline m + E)
+  U \frac{\partial\overline m}{\partial x_3} = 0
\quad \text{in }  \Omega.
\\
%&
\mathbf{u} = 0 \quad \text{on } S ,
\\
\mathbf{u}\cdot\mathbf{n} = 0  \quad \text{on } \Gamma ,
\\
\nu(\overline m + E)[D(\mathbf{u})\mathbf{n}-\mathbf{u}\cdot (D(\mathbf{u})\mathbf{n})\mathbf{n}] =
\mathbf{b}_1  \quad \text{on } \Gamma ,
\\
\theta \frac{\partial \overline m}{\partial \mathbf{n}} - Un_3 \overline m = 0
\quad \text{on } \partial\Omega,
\\
\int_\Omega \overline m dx = 0.
\end{gathered}
\end{equation}
Next, we give the definition of a weak solution of our problem.

\noindent {\bf Definition.}
 Let $\mathbf{f} \in X(\Omega)$; a pair of functions
$(\mathbf{u}, \overline m) \in J_0 (\Omega) \times B$ is called a 
\emph{weak solution} of \eqref{nuevasvariables}  when the following 
two equalities  are satisfied
for all $(\mathbf{v}, \phi) \in J_0(\Omega) \times B$:
\begin{gather}
2(\nu(\overline m+E)D(\mathbf{u}), D(\mathbf{v}))
+ B_0 (\mathbf{u},\mathbf{u},\mathbf{v}) + (\overline m \cdot \chi, \mathbf{v}) 
-2 \int_\Gamma \mathbf{b}_1 \mathbf{v} d\sigma = (\mathbf{f},\mathbf{v}),  \label{e7}
 \\
 \theta(\nabla\overline m, \nabla\phi) + B_1 (\mathbf{u}, \overline m + E, \phi) 
- U(\overline m, \frac{\partial\phi}{\partial x_3}) = 0.
\label{e8}
\end{gather}

This variational formulation is obtained, as usual, by working in a 
formal way. To give an idea on how to do that, here we remark that it 
is obtained by using the following computations done, by simplicity,
on the original form of the first equation. 
Assume that $\mathbf{u}, \mathbf{v} \in \dot{\mathbf{H}}(\Omega)$ and $q\in C^{1}$; then we have
\begin{align*}
&\int_\Omega [-2\operatorname{div} (\nu(m)D(\mathbf{u})) + \nabla q]\mathbf{v} dx
\\
&= 2 \int_{\Omega}\nu(m)D(\mathbf{u}):\nabla \mathbf{v} dx
- 2 \int_{\partial \Omega} \nu(m) D(\mathbf{u}) \mathbf{n} \cdot \mathbf{v} dS
+ \int_{\partial \Omega}q \mathbf{n}\cdot\mathbf{v} dS.
\end{align*}
The third term is zero since $\mathbf{v} \in \dot{\mathbf{H}}(\Omega)$ and thus $\mathbf{v}|_S = 0$ 
and $\mathbf{v}\cdot \mathbf{n}|_{\Gamma} = 0$.
Next, we can split $\nabla \mathbf{v}$ in its symmetric and antisymmetric parts, 
and observe that the symmetric part is exactly $D(\mathbf{v})$;
since $D(\mathbf{u})$ is also a symmetric matrix, a standard computation then 
gives that we can rewrite the first term as
\[
\int_{\Omega} \nu (m) D(\mathbf{u}):\nabla \mathbf{v} dx = \int_{\Omega} \nu (m) D(\mathbf{u}):D(\mathbf{v})dx.
\]
On the other hand, by using the tangential and normal components at 
each point of the boundary, the second term becomes
\begin{align*}
& -2 \int_{\partial \Omega} \nu(m)D(\mathbf{u}) \cdot \mathbf{n} \  \mathbf{v} d S
\\
& = -2 \int_{\partial \Omega} ( \nu(m)D(\mathbf{u}) \mathbf{n}) \cdot \mathbf{n} \  (\mathbf{v} \cdot \mathbf{n}) dS
 \\
&\quad  - 2 \int_{\partial \Omega}[ \nu(m)D(\mathbf{u})\mathbf{n} - (\nu(m)D(\mathbf{u})\mathbf{n}) \cdot \mathbf{n}) \mathbf{n}] \cdot(\mathbf{v} - (\mathbf{v}\cdot\mathbf{n})\mathbf{n})dS
\\
&\quad - 2 \int_{\Gamma} \nu(m)[D(\mathbf{u})\mathbf{n} - \mathbf{n} \cdot (D(\mathbf{u})\mathbf{n})  \mathbf{n}] \cdot \mathbf{v} dS \\
&= - 2 \int_{\Gamma} \mathbf{b}_1 \cdot \mathbf{v} dS
\end{align*}
were we have used the fact that $\mathbf{v} \in \dot{\mathbf{H}}(\Omega)$ and the boundary 
condition on $\Gamma$.
Thus,
\[
\int_\Omega [-2\operatorname{div} (\nu(m)D(\mathbf{u})) + \nabla q] \mathbf{v} dx 
= 2(\nu(m)D(\mathbf{u}),D(\mathbf{v}))- 2\int_{\Gamma}\mathbf{b}_{1}\cdot\mathbf{v} dS.
\]

Analogously, the second equation in the model (already with the previous
 change of variable) can be formally treated as
\begin{align*}
& \int_{\Omega}-\theta \Delta \overline m  \phi\,dx 
+ \int_{\Omega}U\frac{\partial\overline m}{\partial x_3}  \phi\,dx \\
& = \theta\int_{\Omega}\nabla\overline m  \nabla \phi\,dx- \int_{\partial\Omega}\theta\frac{\partial\overline m}{\partial n}\phi\,dS - U  \int_{\Omega}\overline m\frac{\partial \phi}{\partial x_3} dx + \int_{\partial\Omega}U n_{3}\overline m \phi\,dS.\\
& = \theta\int_{\Omega}\nabla\overline m  \nabla \phi  dx - U
\int_{\Omega}\overline m\frac{\partial \phi}{\partial x_3} dx\,.
\end{align*}
Then we have the following result.

\begin{theorem} \label{thm2.1}
Let $\nu$ be a continuous function satisfying
\begin{equation} \label{e9}
\nu_0 = \inf \{\nu(m), m\in \mathbb{R}\} > 0, \quad 
\nu_1 = \sup \{\nu(m), m \in \mathbb{R}\} < +\infty ;
\end{equation}
$\mathbf{f} \in X(\Omega)$ and $\frac U\theta <(C_P)^{-1}$,
where $C_{P}$ is the Sobolev constant appearing in the 
Poincar\`e-Friedrichs inequality.
Then there exists a weak solution of  Problem \eqref{nuevasvariables} 
satisfying
\[
|D(\mathbf{u})|^2 + |\nabla \overline m|^2 \leq C (|\mathbf{f}|^2 +
\|\mathbf{b}_1\|^2_{H^{1/2}(\Gamma)} + |\nabla E|^2),
\]
with a constant $C$ independent of $\mathbf{f}$, $\mathbf{b}_1$ and $E$.
\end{theorem}

\begin{proof} 
We consider the following Schauder bases formed by the eigenfunctions 
described in the previous section:
$(\overline{\mathbf{w}}^j)_1^\infty$ for $\mathbf{J}_0(\Omega)$ and 
$(\overline \phi^j)_1^\infty$ for $B$.
For each $n \in \mathbb{N}$, we define
$\mathbf{W}_n = \operatorname{span} \{\overline{\mathbf{w}}^j, \ 1 \leq j \leq n \}$ and
$M_n = \operatorname{span} \{\phi^\ell, \ 1 \leq \ell \leq n \}$
and consider the Galerkin approximations
\[
\mathbf{u}^n = \sum^n_{j=1} c_{n,j} \overline{\mathbf{w}}^j \in \mathbf{W}_n.
\quad \overline m^n = \sum^n_{\ell=1} d_{n, \ell} \overline \phi^\ell \in M_n,
\]
satisfying the following approximate problem
\begin{gather}
\begin{aligned}
&2(\nu(\overline m^n+E)D(u^n), D(\overline v))
+ B_0(\mathbf{u}^n, \mathbf{u}^n, \overline v)\\
&+ (\overline m^n \cdot \chi, \overline v)
- 2 \int_\Gamma \mathbf{b}_1 \overline v d\sigma
= (\mathbf{f},\overline v),
\end{aligned} \label{e10}
\\
\theta (\nabla \overline m^n, \nabla \overline \phi)
+ B_1 (\mathbf{u}^n, \overline m^n + E, \overline \phi)
- U(\overline m^n, \frac{\partial \overline \phi}{\partial x_3}) = 0,
\label{e11}
\end{gather}
for all $\overline v \in \mathbf{W}_n$ and all $\overline \phi \in M_n$.

Firstly, by assuming the existence of  $(\mathbf{u}^n, \overline m^n)$ for all 
$n \in \mathbb{N}$ (such existence will be proved later on),
we will prove that they indeed converge, along subsequences, to a solution 
of our problem.
To do this, it will be necessary to obtain estimates for the gradients 
of the unknowns.

We set  $\overline v =  u^n$ in \eqref{e10} to obtain
\begin{equation}
\begin{aligned}
& 2(\nu(\overline m^n+E) D(\mathbf{u}^n), D(\mathbf{u}^n)) + B_0
(\mathbf{u}^n, \mathbf{u}^n, \mathbf{u}^n)\\
&+ (\overline m^n\cdot \chi, \mathbf{u}^n)
 - 2 \int_\Omega \mathbf{b}_1 \mathbf{u}^n d\sigma = (\mathbf{f},
\mathbf{u}^n).
\end{aligned}
\end{equation}
By observing that $B_0 (\mathbf{u}^n, \mathbf{u}^n, \mathbf{u}^n) =
0$ (see \eqref{e3}) and using H\"older inequality together
with \eqref{e9},
its follows that
\begin{align*}
2\nu_0 |D(\mathbf{u}^n)|^2
& \leq  |(\overline m^n \cdot \chi, \mathbf{u}^n)+ (\mathbf{f}, \mathbf{u}^n)
 + 2 \int_\Omega \mathbf{b}_1 \mathbf{u}^n d \sigma|\\
& \leq   |\overline m^n| \ |\mathbf{u}^n| + |\mathbf{f}| \ |\mathbf{u}^n|
+ 2c(\Gamma) \Big(\int_\Gamma \mathbf{b}^2_1 d\sigma\Big)^{1/2} |D(\mathbf{u}^n)|.
\end{align*}
By the Sobolev embeddings,  Korn, Young and Poincar\`e-Friedrichs inequalities,
we then obtain
\begin{equation} \label{e14}
\nu_0 |D(\mathbf{u}^n)|^2
\leq \frac{3}{4\nu_0} (C^2_p \gamma |\nabla \overline m^n|^2 + \gamma |\mathbf{f}|^2 + 4c^2
(\Gamma)\|\mathbf{b}_1\|^2_{H^{1/2}(\Gamma)}).
\end{equation}

Next, for each $n$ we consider the orthogonal projection $\overline P_n$ 
defined in \eqref{ProjectionBar-n}, denote by $|\Omega|$ the Lebesgue 
measure of $\Omega$, and take 
$\overline \phi = \overline m^n + \overline P_n(E - \alpha/|\Omega|) $
to get
\begin{equation}\label{estimativaTheta}
\begin{aligned}
&\theta |\nabla \overline m^n|^2
+ \theta (\nabla \overline m^n, \nabla \overline P_n (E- \alpha/|\Omega|) )
+ B_1 (\mathbf{u}^n, \overline m^n + E, \overline m^n \\
&+ \overline P_n (E - \alpha/|\Omega|))
 - U(\overline m^n, \frac{\partial}{\partial x_3} 
(\overline m^n + \overline P_n ( E - \alpha/|\Omega| ) )) = 0.
\end{aligned}
\end{equation}
Now, we observe that
\begin{equation}\label{splitting}
\begin{aligned}
&B_1 (\mathbf{u}^n, \overline m^n + E , \overline m^n 
 + \overline P_n (E - \alpha/|\Omega| ))
\\
&= B_1 (\mathbf{u}^n, \overline m^n  + E -  \alpha/|\Omega|  , \overline m^n 
+ \overline P_n (E - \alpha/|\Omega| ))
\\
&=  \hspace{0.1cm} B_1 (\mathbf{u}^n, \overline m^n 
+ \overline P_n ( E  - \alpha|/\Omega|  ) ,  \overline m^n 
+ \overline P_n ( E  - \alpha/|\Omega|  ) )
\\
&\quad + B_1 (\mathbf{u}^n,  E  - \alpha|/\Omega| 
 - \overline P_n ( E  - \alpha|/\Omega|  ) ,  \overline m^n )
\\
&\quad  + B_1 (\mathbf{u}^n, E  - \alpha|\Omega| 
 - \overline P_n ( E  - \alpha|/\Omega|  ) ,  
\overline P_n (E - \alpha/|\Omega| )),
\end{aligned}
\end{equation}

We also remark that, due to the fact that $(1, \overline \phi_j) = 0$ 
for all $j \geq 1$, that the gradients of a constant is zero,
and that $ E \in D(\tilde A_1)$,  by using \eqref{Projection-prop3}, 
we obtain
\begin{equation}\label{Projection-prop4}
 |\nabla (E- \alpha/|\Omega|  - \overline P_n (E- \alpha/|\Omega|) ) |
 = |\nabla (E  - \overline P_n E ) |
\leq  \frac{\theta}{\beta_{n+1}^{1/2} (\theta - 2 U C_P)^{1/2} } | \Delta E |
 \end{equation}
Using \eqref{splitting} in \eqref{estimativaTheta}, recalling that
$B_1 (\mathbf{u}^n, \overline m^n + \overline P_n ( E  - \alpha|/\Omega|  ) ,  
\overline m^n + \overline P_n ( E  - \alpha/|\Omega|  ) ) = 0$ 
(see \eqref{e3} ),
suitably using H\"older inequality and properties \eqref{Projection-prop1} 
and \eqref{Projection-prop4}
and the fact that $z \leq 1 + z^2$ for all real $z$, together with the
 hypothesis that $\theta - U C_P > 0$, we obtain that
\begin{equation}
\label{e15}
|\nabla\overline m^n|^2
\leq
K_ 1 \Big(\frac{1}{\beta_{n+1}} + \frac{1}{\beta_{n+1}^{1/2}} \Big) 
| \Delta E | |D( u^n )| + K_ 2 |\nabla E |^2 ,
\end{equation}
with positive constants $K_ 1$ and $K_2$ that do not depend on $n$.

By multiplying \eqref{e15} by $ \frac{3}{2\nu_0} C^2_p$,
adding the result to \eqref{e14}, for $n \geq N_0$ such that
\[
  \frac{3}{2\nu_0} C^2_p K_ 1 \Big(\frac{1}{\beta_{N_0+1}} 
+ \frac{1}{\beta_{N_0+1}^{1/2}} \Big) |\Delta E | < \frac{\nu_0}{2},
\]
 we  have
\begin{equation} \label{e17}
|D(\mathbf{u}^n)|^2 + |\nabla\overline m^n|^2 \leq C(|\mathbf{f}|^2
+ \|\mathbf{b}_1\|^2_{H^{1/2}\Gamma)} + |\nabla E|^2).
\end{equation}
Therefore, the sequence $\{(\mathbf{u}^n, \overline m^n)\}$ is
bounded in $\mathbf{J}_0(\Omega) \times B$.

Next, since $\mathbf{J}_0(\Omega)$ is compactly immersed in 
$\mathbf{X}(\Omega)$, and $B$ is compactly immersed in $Y$,
there are elements  $\mathbf{u} \in \mathbf{J}_0(\Omega), \overline m \in B$
and a subsequence, which for simplicity we still denote by 
$\{(\mathbf{u}^n, \overline m^n)\}$), such that
\begin{gather*}
\mathbf{u}^n \to \mathbf{u} \quad \text{weakly in $\mathbf{J}_0(\Omega)$ and strongly in
$\mathbf{X} (\Omega)$},
\\
\overline m^n \to m  \quad \text{weakly in $B$ and strongly in } Y,
\\
D(\mathbf{u}^n) \to D(\mathbf{u}) \quad \text{weakly in } (L^2(\Omega))^9,
\\
\nabla(\overline m^n) \to \nabla (\overline m) \quad \text{weakly in } 
 (L^2(\Omega))^3.
\end{gather*}
These convergences are enough to allow us to take the limit as 
$n \to +\infty$
in \eqref{e10} and \eqref{e11}  to obtain
\begin{equation}\label{e18}
\begin{gathered} 
2(\nu(\overline m+E)D(\mathbf{u}), D(\overline{\mathbf{w}}^j)) + B_0(\mathbf{u}, \mathbf{u}, \overline{\mathbf{w}}^j)
\\
+ (\overline m \cdot \chi, \overline{\mathbf{w}}^j) 
 - 2  \int_\Gamma \mathbf{b}_1 \cdot\overline{\mathbf{w}}^j d\sigma = (\mathbf{f}, \overline{\mathbf{w}}^j),
\\
\theta (\nabla\overline m, \nabla \overline \phi^\ell) + B_1 (\mathbf{u},
\overline m+E, \overline \phi^\ell) - U\Big(\overline m,
\frac{\partial \overline \phi^\ell}{\partial
x_3}\Big) = 0,
\end{gathered}
\end{equation}
for all $j, \ell \in \mathbb{N}$.

In fact, it is well known   \cite{temam} that with the previous
 convergences we obtain
\begin{gather*}
 B_0 (\mathbf{u}^n, \mathbf{u}^n, \mathbf{w}) \to B_0 (\mathbf{u}, \mathbf{u}, \mathbf{w}), \quad
  \forall \mathbf{w} \in \mathbf{J}_0(\Omega),\\
 B_1 (\mathbf{u}^n, \overline m^n, \phi) \to
B_1(\mathbf{u}, \overline m, \phi), \quad \forall \phi \in B .
\end{gather*}
We also note that
\[
(\nu(\overline m^n+E)D(\mathbf{u}^n), D(\mathbf{w}^j)) \to (\nu(\overline m+E) D(\mathbf{u}), D(\mathbf{w}^j))
\]
since 
\begin{gather*}
(\nu(\overline m^n+E)D(\mathbf{u}^n), D(\mathbf{w}^j)) 
= (D(\mathbf{u}^n),\nu(\overline m^n+E)D(\mathbf{w}^j)),\\
(D(\mathbf{u}), \nu(\overline m + E) D(\mathbf{w}^j))= (\nu(\overline m+E) D(\mathbf{u}), D(\mathbf{w}^j)),\\
\nu(\overline m^n + E) D(\mathbf{w}^j) \to \nu (\overline m + E)D(\mathbf{w}^j)
\end{gather*}
strongly in $(L^2(\Omega))^9$ due to
the Lebesgue dominated convergence theorem.

Finally, as the $\{\overline{\mathbf{w}}^j\}$ and
 $\{\overline \phi^\ell\}$ are Schauder bases, respectively in 
$\mathbf{J}_0(\Omega)$ and $B$, by using \eqref{e18},
we conclude that $(\mathbf{u}, \overline m)$ satisfies \eqref{e7}, \eqref{e8},
and thus $(\mathbf{u}, \overline m)$ is the required weak solution.

\subsection{Existence of approximate solutions}

It remains to prove that, for each $n \in \mathbb{N}$, equations 
\eqref{e10}-\eqref{e11} have
solutions. For this, we proceed similarly as in 
Lorca and Boldrini \cite{lorca}.

Let $\mathbf{W}_n$ and $M_n$ be as in the last section.
Given any $(\mathbf{z}, \xi) \in \mathbf{W}_n \times M_n$,
we consider the unique solution $(\mathbf{v},\Psi) \in \mathbf{W}_n \times M_n$ of the
linearized equations
\begin{gather}
\begin{aligned}
&2(\nu(\xi+E)D(\mathbf{v}), D(\overline{\mathbf{w}}^j)) + B_0 (\mathbf{z}, \mathbf{v},\overline{\mathbf{w}}^j)\\
&+ (\Psi \cdot \chi,\overline{\mathbf{w}}^j) 
-2 \int_\Gamma \mathbf{b}_1 \cdot \overline{\mathbf{w}}^j d\sigma 
- (\mathbf{f}, \overline{\mathbf{w}}^j) = 0, 
\end{aligned}\label{e19}
 \\
 \theta(\nabla \Psi, \nabla \overline \phi^\ell) + B_1
(\mathbf{z}, \Psi + E,\overline \phi^\ell) - U
\Big(\Psi, \frac{\partial \overline \phi^\ell}{\partial
x_3}\Big) = 0, \label{e20}
\end{gather}
for $1 \leq j,\ell \leq n$.

To prove that there is only one such solution  
$(\mathbf{v}, \Psi) \in \mathbf{W}_n \times M_n$,
we observe that \eqref{e19}, \eqref{e20} constitute in fact a 
linear system with $2n$ equations for the $2n$ coefficients of the expansions
\[
\mathbf{v} = \sum^n_{j=1} c_j \overline{\mathbf{w}}^j, \quad
\Psi = \sum^n_{\ell=1} d_\ell \overline \phi^\ell .
\]
Thus, to show the existence and uniqueness of solutions of  
system \eqref{e19}, \eqref{e20},
it is enough to prove that the only solution of its associated homogeneous 
linear system, that is, the corresponding equations with $ \mathbf{b}_1 = 0$ 
and $\mathbf{f}=0$, is the trivial null solution.
For this, let $(\mathbf{v}, \Psi)$ be any solution of the such homogeneous system, 
and $c_j$ and $d_\ell$ its corresponding coefficients as in the previous 
expansions.
Then,  by multiplying \eqref{e19} by $c_j$ and \eqref{e20} by $d_\ell$, 
and adding in $j$ and $\ell$ from 1 to $n$, we obtain
\begin{gather}
2 (\nu(\xi+E) D(\mathbf{v}), D(\mathbf{v})) + B_0 (\mathbf{z}, \mathbf{v}, \mathbf{v})
+ (\Psi \cdot \chi, \mathbf{v}) = 0, \label{e21} \\
\theta |\nabla \Psi|^2 + B_1 (\mathbf{z}, \Psi, \Psi) 
- U \Big(\Psi, \frac{\partial\Psi}{\partial x_3}\Big) = 0. \label{e22}
\end{gather}

As $B_1 (\mathbf{z}, \Psi, \Psi) = 0$, by using  H\"older and 
Poincar\`e-Friedrichs inequalities, \eqref{e22} becomes
$(\theta - U C_P) |\nabla\Psi|^2 \leq 0 $.
Thus, $(\theta-UC_P) > 0$ implies  $|\nabla W|^2=0$;
i.e., $\Psi$ is constant; since $\int_\Omega \Psi dx = 0$,
we conclude that $\Psi = 0$.
Since $B_0 (\mathbf{z}, \mathbf{v}, \mathbf{v}) = 0$ and $\Psi = 0$, \eqref{e21} 
gives $2 \nu_0 |D(\mathbf{u})|^2 \leq 0$, which in turn implies that $\mathbf{v} = 0$.


Thus, for each $n$ we have a well defined operator
\begin{equation}\label{T}
T_n  : \mathbf{W}_n \times M_n \to \mathbf{W}_n \times M_n,
\end{equation}
such that to each  $(\mathbf{z}, \xi) \in \mathbf{W}_n \times M_n$ associates
$T(\mathbf{z}, \xi) = (\mathbf{v},\Psi)$, where $(\mathbf{v},\Psi) \in \mathbf{W}_n \times M_n$ 
is the unique solution of \eqref{e19}.
Moreover, since $\mathbf{W}_n$ and $M_n$ are finite dimensional vector spaces,
it is rather standard to prove that $T_n$ is continuous for any chosen norms.

Next, by proceeding exactly as before in the derivation of the estimates 
for $(\mathbf{u}^n,\overline m^n)$, we obtain for system \eqref{e19}, \eqref{e20} 
 the same kind of estimate as the one in \eqref{e17}; i. e.,
$$
|D(\mathbf{v})|^2 + |\nabla \Psi|^2 \leq C (|\mathbf{f}|^2 +
\|\mathbf{b}_1\|^2_{H^{1/2}(\Gamma)} + |\nabla E|^2) =  R_1^2,
$$
with $R_1$ independent of $n$ and $(\mathbf{z}, \xi) \in \mathbf{W}_n \times M_n$.
By denoting
$F_n = \{(\mathbf{z}, \xi) \in \mathbf{W}_n
\times M_n; (|D(\mathbf{z})|^2 + |\nabla\xi|^2) \leq R_1^2\}$,
and restricting the operator $T_n$ to such $F_n$. we have a continuous operator
$T_n : F_n \to F_n$ acting from a finite dimensional closed convex 
set convex $F_n$ into itself.
Brower fixed point theorem then gives us the existence of at least one 
fixed point, $(\mathbf{u}, \overline m)$,
which is a solution of \eqref{e10}-\eqref{e11}. This completes the proof.
\end{proof}

\section{Existence of strong solutions}

Here we prove the existence of solutions that are more regular 
than the ones obtained in the previous section.
The main difficulty for this will be to obtain the necessary higher 
order estimates for the nonlinear terms present in the equations.

\begin{theorem}\label{thm2.2}
Assume that $\mathbf{b}_1 = 0$, $\mathbf{f} \in \mathbf{X}(\Omega)$ and  $\nu$ is a
 $C^1$-function satisfying \eqref{e9}.
Then, for a small enough $U$, there exists a strong solution of 
\eqref{nuevasvariables}; that is,
there exists a pair of functions  
$(\mathbf{u},\overline m) \in (J_0(\Omega) \cap H^2(\Omega)) 
\times (Y \cap H^2(\Omega))$ 
satisfying
\begin{gather*}
 P[-2 \operatorname{div} (\nu(\overline m + E)D(\mathbf{u})) + \mathbf{u}\cdot \nabla \mathbf{u}
+ \overline m \cdot \chi - \mathbf{f}] = 0 \quad \text{in } (L^2 {(\Omega))^3,}
\\
 \overline P \big[-\theta \Delta \overline m + \mathbf{u}\cdot
\nabla (\overline m + E) + U \frac{\partial\overline{m}}{\partial
x_3}\big] = 0 \quad \text{in } L^2 (\Omega).
\end{gather*}
\end{theorem}

\begin{proof}
We start by recalling that the $L^2(\Omega)$-norm of the Stokes operator
 $A$ (respectively the operator $A_1$)
applied to an element and the norm of 
$\mathbf{J}_0(\Omega) \cap (H^2 (\Omega))^3$
(respectively the norm of $Y \cap H^2 (\Omega)$ ) 
of the same element are equivalent.

Next, we repeat the construction used in the proof of Theorem \ref{thm2.1}
to show the existence of approximations $\mathbf{u}^n$ and $\overline m^n$.
In the present situation, however, under the conditions of 
Theorem \ref{thm2.2}, we will show that each of the previous
operators $T_n$ admits fixed points satisfying an estimate independent 
of $n$ in a more regular space.

Since all  the estimates obtained in the previous section hold true
 with the same proofs,  we proceed with the derivation of  
further $H^2(\Omega)$-estimates for $(\mathbf{v}, \Psi)$ solution of 
\eqref{e19}- \eqref{e20}.

For this, we multiply \eqref{e19} by $\alpha_j c_j$, \eqref{e20} by
$\beta_\ell d_\ell$ and add in $j$ and $\ell$ from 1 to $n$ to obtain
\begin{gather}
 -2 (\operatorname{div}(\nu(\xi+E)D(\mathbf{v})), A\mathbf{v}) + B_0 (\mathbf{z}, \mathbf{v}, A\mathbf{v})
 + (\overline m \cdot \chi, A\mathbf{v}) - (\mathbf{f}, A\mathbf{v}) = 0, \label{e23}
\\
 \theta (A_1 \Psi, A_1 \Psi) + B_1(\mathbf{z}, \Psi + E, A_1 \Psi)  
+ U \left( \frac{\partial}{\partial x_3} \Psi, A_1 \Psi\right) = 0, 
\label{e24}
\end{gather}
again with $B_0$ and $B_1$ given by \eqref{e2}.

Using the identity $2\operatorname{div} (\nu(\xi+E) D(\mathbf{v}))
= \nu(\xi+E) \Delta \mathbf{v} + 2\nu'(\xi+E) \nabla (\xi+E) D(\mathbf{v})$
in \eqref{e23}, we have
\begin{equation}
\begin{aligned}
-(\nu(\xi+E)\Delta \mathbf{v}, A\mathbf{v})
& =  -B_0 (z, \mathbf{v}, A\mathbf{v}) - (\Psi \cdot \chi, A\mathbf{v}) +(\mathbf{f}, A\mathbf{v}) \\
&\quad +(2\nu'(\xi+E) \nabla (\xi+E) D(\mathbf{v}), A\mathbf{v}).
\end{aligned}\label{e25}
\end{equation}
Next, from Helmholtz decomposition, we know that there exists
$\overline q \in H^1(\Omega)$ such that
 $-\Delta \mathbf{v} = A\mathbf{v} + \nabla \overline q$,
and we have the estimate
\begin{equation}\label{e26}
\|\overline q\|_1 \leq c |\Delta \mathbf{v}| .
\end{equation}
Therefore, \eqref{e25} becomes
\begin{align*}
(\nu(\xi+E) A\mathbf{v}, A\mathbf{v})
& =  - B_0 (\mathbf{z}, \mathbf{v}, A\mathbf{v})- (\Psi \cdot \chi, A\mathbf{v}) + (\mathbf{f}, A\mathbf{v})
\\
&\quad + 2(\nu'(\xi+E)\nabla(\xi+E)D(\mathbf{v}), A\mathbf{v})
 - (\nu(\xi+E)\nabla \overline q, A\mathbf{v}).
\end{align*}
Using Korn, H\"older and Poincar\`e-Friedrichs inequalities and
 Sobolev embeddings, together with \eqref{e9}, we obtain
\begin{align*}
\nu_0|A\mathbf{v}|^2 & \leq  C_0 |D(\mathbf{z})| |A\mathbf{v}|^2
+ \overline C_P |A_1 \Psi||A\mathbf{v}| + |\mathbf{f}||A\mathbf{v}|  \\
&\quad + 2\nu_1' (|A_1\xi| + |\nabla E|_4)|A\mathbf{v}|^2) +
|(\nu(\xi+E)\nabla \overline q, A \mathbf{v})|
\end{align*}
where $\nu_1'= \sup \{|\nu'(r)|, r \in \mathbb{R}\}$.
Since $A\mathbf{u} \in \mathbf{W}_n$, we observe that
$$
(\nu(\xi+E) \nabla \overline q, A\mathbf{v}) = -(\overline q, \ \operatorname{div}
(\nu(\xi+E)A\mathbf{v})) = - (\overline q, (\nu '(\xi+E) \nabla(\xi+E) A\mathbf{v})),
$$
and thus
\begin{equation} \label{e28}
|(\nu(\xi+E)\nabla \overline q, A\mathbf{v})|
\leq \nu_1'|\overline q|_4 |\nabla (\xi+E)|_4 |A\mathbf{v}|
\leq  c \nu_1'(|A_1 \xi| + |\nabla E|_4) \|\overline q\|_1 A\mathbf{v}| .
\end{equation}
Combining \eqref{e26}-\eqref{e28}, using Young  and Poincar\`e-Friedrichs
inequalities, its follows that
\begin{equation} \label{e29}
\frac{\nu_0}{2} |A\mathbf{v}|^2
\leq C_1 (|A\mathbf{z}|^2 + |A_1\xi|^2 + |\nabla E|^2_4) \ |A\mathbf{v}|^2
+ \frac{2\overline C^2_P}{\nu_0} |A_1 \Psi|^2 + \frac{2}{\nu_0} |\mathbf{f}|^2.
\end{equation}
Similarly, by using H\"older, Korn and Friedrichs inequalities
and Sobolev embeddings, from \eqref{e24} we obtain
\begin{equation} \label{e30}
\theta |A_1 \Psi|^2
\leq C_2 |A\mathbf{z}|^2 (|A_1 \Psi|^2 + |\nabla E|^2_4)
+ U\overline C_P |A_1 \Psi|^2 .
\end{equation}
Multiplying \eqref{e29} by
$\delta = \frac{\nu_0}{2\overline C^2_P} (\theta - U\overline
C_P)$ and adding the result to \eqref{e30}, we obtain
\begin{align*}
&\frac{\nu^2_0}{4\overline C^2_P} (\theta - U\overline C_P) |A\mathbf{v}|^2
+ \frac12 (\theta - U\overline C_P) |A_1 \Psi|^2
\\
&\leq \delta (C_1(|A\mathbf{z}|^2 + |A_1\xi|^2 + |\nabla E|^2_4) |A\mathbf{v}|^2
+ \frac{2}{\nu_0} |\mathbf{f}|^2) + C_2
|A\mathbf{z}|^2 (|A_1 \Psi|^2 + |\nabla E|^2_4) .
\end{align*}
Since
$\frac U\theta < (\overline C_P)^{-1}$,  there exists a positive
constant $C_4$ such that
\[
|A\mathbf{v}|^2 + |A_1 \Psi|^2
\leq C_4([|A\mathbf{z}|^2 + |A_1 \xi|^2 + |\nabla E|^2_4] |A\mathbf{v}|^2 + |\mathbf{f}|^2
+ |A\mathbf{z}|^2 (|A_1 \Psi|^2 + |\nabla E|^2_4)),
\]
which implies
\[
|A\mathbf{v}|^2 + |A_1 \Psi|^2
\leq
C_4(|A\mathbf{z}|^2  +|A_1\xi|^2 + |\nabla E|^2_4)(|A\mathbf{v}|^2 + |A_1\Psi|^2
+ |\nabla E|^2_4)  + C_4|\mathbf{f}|^2 .
\]
This expression can be rewritten for a positive constant $C$ as
\begin{equation} \label{e31}
[1-C(|A\mathbf{z}|^2  +|A_1\xi|^2](|A\mathbf{v}|^2 + |A_1\Psi|^2)
\leq
C(|A\mathbf{z}|^2  +|A_1\xi|^2 ) |\nabla E|^2_4 + C |\nabla E|^4_4+ C|\mathbf{f}|^2 ,
\end{equation}

Next, we will use \eqref{e31} to find $R_2 >0$ and suitable conditions 
that will guarantee  that when
$|A\mathbf{z}|^2 + |A_1\xi|^2 \leq R_2^2$ we also have 
$|A\mathbf{v}|^2 + |A_1 \Psi|^2 \leq R_2^2$.
For this, let $R_2^2= 1/(4C)$ and then take $U$ small enough such 
that $|\nabla E|_4^2 \leq 1/(4C)$.
Under these choices, when $|A\mathbf{z}|^2 + |A_1\xi|^2 \leq R_2^2$, 
inequality \eqref{e31} implies 
\[
(1/2)(|A\mathbf{v}|^2 + |A_1\Psi|^2)
\leq
(1/4) R_2^2 + (1/64C)+ C|\mathbf{f}|^2  = (1/4) R_2^2 + (1/16)R_2^2+ C|\mathbf{f}|^2,
\]
which is rewritten as $(|A\mathbf{v}|^2 + |A_1\Psi|^2) \leq  (5/8)R_2^2+ C|\mathbf{f}|^2$.
Hence, when $|\mathbf{f}| \leq (3/8) R_2^2 = 3/(32C)$, we obtain that
$(|A\mathbf{v}|^2 + |A_1\Psi|^2 \leq R_2^2$, as required, and such $R_2$ 
is independent of $n$ and $(\mathbf{z}, \xi)$.

Under these conditions, we can consider again the operators $T_n$, 
but now as continuous operators
$T_n: G_n \to G_n$ acting from a finite dimensional closed convex set
$G_n = \{(\mathbf{z},\xi) \in W_n \times M_n;  |A\mathbf{z}|^2 + |A_1 \xi|^2 \leq R_2^2\}$ 
in itself.

Now, proceeding again similarly as before, we can apply Brower 
fixed point theorem for each of such operator $T_n$;
this gives a sequence of approximate solutions  $(\mathbf{u}^n, \overline m^n)$ 
satisfying \eqref{e10}, \eqref{e11},
which are now uniformly bounded in $H^2(\Omega)$ .
By taking the limit as $n \to + \infty$ and proceeding exactly as before,
we obtain a solution of our original problem. This proves the theorem.
\end{proof}



\subsection*{Remarks}
To each given solution of Problem \eqref{nuevasvariables}, 
it is trivially associated a solution of  \eqref{sistema1}-\eqref{sistema2}.
Moreover, the regularities obtained for $(\mathbf{u}, \overline m)$  
also hold for $(\mathbf{u}, m)$.

As in Temam \cite{temam},  there exists a unique function $q$ 
(the hydrostatic pressure) in
$H^1(\Omega) \cap L^2_0(\Omega)$, where $L^2_0(\Omega) = \{h \in
L^2(\Omega); \ (h,1)=0\}$, which satisfies
$$
-2 \operatorname{div} (\nu(m) D(\mathbf{u})) + \mathbf{u}\cdot  \nabla
\mathbf{u} + m \cdot \chi - \mathbf{f} = - \nabla q .
$$

Finally, the positivity of the concentration established by
 Kan-On, Narukawa and Teramoto  in \cite{kan}
also holds true for our problem, with exactly the same proof. 
That is, we have
\begin{lemma}[Positivity of the concentration] \label{Lema2.1}
Let $(\mathbf{u}, q, m)$ be a solution of  \eqref{sistema1}-\eqref{sistema2} 
given by Theorem \ref{thm2.2}.
Then $m(x)> 0$ for all $x \in \Omega$.
\end{lemma}


\section{Uniqueness of the solution}

\begin{theorem} \label{thm2.3}
Assume that $\nu(\cdot)$ is  a Lipschitz continuous function and that 
$(\mathbf{u}, \overline m)$ is a weak solution of the Problem \eqref{nuevasvariables}
such that $(|D(\mathbf{u})| + |A_1 \overline m| + |\nabla E|_4))$ 
is small enough; then such solution is unique.
\end{theorem}

Note that as before, this last requirement on $|\nabla E|_4$ can be attained 
for small enough $U$.


\begin{proof} 
Let $(\mathbf{u}_1, \overline m_1), (\mathbf{u}_2, \overline m_2)$ be two weak solutions 
of \eqref{nuevasvariables}, with
$\overline m_1$ and $\overline m_2$ in $(H^2(\Omega) \cap Y)$ and both
 satisfying the stated smallness condition.
By denoting
$\mathbf{z} = \mathbf{u}_1-\mathbf{u}_2$ and $\varphi= \overline m_1 - \overline m_2$,
then we have that $\mathbf{z} \in \mathbf{J}_0$,  $\varphi \in (H^2 (\Omega) \cap Y)$
and the pair $(\mathbf{z}, \varphi)$ satisfies, for all  $\mathbf{w} \in \mathbf{J}_0(\Omega)$ 
and  $\phi \in L^2(\Omega)$,
the following equations:
\begin{equation}\label{e32}
\begin{gathered}
2(\nu(\overline m_1) D(\mathbf{z}), D(\mathbf{w}))
+ 2((\nu(\overline m_1) - \nu(\overline m_2)) D(\mathbf{u}_2), D(\mathbf{w}))
\\
+ B_0 (\mathbf{z}, \mathbf{u}_1, \mathbf{w})  + B_0 (\mathbf{u}_2, \mathbf{z}, \mathbf{w}) + (\varphi \cdot \chi, \mathbf{w}) = 0,
\\
\theta (\nabla\varphi, \nabla\phi) + B_1 (\mathbf{z}, \overline m_1 + E, \phi)
+ B_1 (\mathbf{u}_2, \varphi, \phi) - U 
\Big(\varphi, \frac{\partial\phi}{\partial x_3}\Big) = 0 .
\end{gathered}
\end{equation}
By setting $\mathbf{w} = \mathbf{z}$ and $\phi = - A_1 \varphi$ in
\eqref{e32}, we obtain
\begin{equation}
\begin{aligned}
2 \nu_0 |D(\mathbf{z})|^2 
& \leq  |2 ((\nu(\overline m_1) - \nu(\overline m_2)) D(\mathbf{u}_2), D(\mathbf{z}))
+ B_0 (\mathbf{z},\mathbf{u}_1, \mathbf{z}) + (\varphi \cdot \chi, \mathbf{z})|
\\
& \leq  2C_1 |\varphi|_\infty |D(\mathbf{u}_2)| \ |D(\mathbf{z})| 
 + C_0 |D(\mathbf{z})| \ |D(\mathbf{u}_1)| \ |D(\mathbf{z})| \\
& \quad + \overline C_P C_P |A_1 \varphi| \ |D(\mathbf{z})|,
\end{aligned}\label{e33}
\end{equation}
\begin{align*}
\theta |A_1 \varphi|^2  
& \leq  \Big| B_1(\mathbf{z}, \overline m_1 + E,
A_1 \varphi) + B_1(\mathbf{u}_2, \varphi, A_1 \varphi)
 + U \Big( \frac{\partial\varphi}{\partial x_3} , A_1 \varphi \Big) \Big|
\\
&  \leq  C_0 |D(\mathbf{z})| (|A_1 \overline m_1| + |\nabla E|_4) |A_1\varphi|
 + C_0 |D(\mathbf{u}_2)| \ |A_1\varphi|^2
+ U\overline C_P |A_1\varphi|^2,
\end{align*}
where $C_1$ is the Lipschitz constant of $\nu(\cdot)$; that is,
$|\nu(r) - \nu(s)| \leq C_1 |r-s|$, for all $r, s \in \mathbb{R}$.
Thus, by taking $U$ small enough so that $\theta - U \overline C_P > 0$, 
we obtain
$$
|A_1\varphi|
\leq [C_0 /  (\theta-U\overline C_P)]
(|D(\mathbf{z})|(|A_1 \overline m_1| + |\nabla E|_4) + |D(\mathbf{u}_2)| \ |A_1 \varphi|),
$$
which, under the condition that $(C_0/(\theta-U\overline C_P)) |D(\mathbf{u}_2)| < 1/2$, 
gives
\begin{equation} \label{e35}
|A_1 \varphi|
\leq [C_0 / (\theta-U\overline C_P)] (|A_1 \overline m_1| + |\nabla E|_4) |D(\mathbf{z})|.
\end{equation}

Since $|\varphi|_\infty \leq C |A_1 \varphi|$, by combining \eqref{e33} and
\eqref{e35}, we obtain
$|D(\mathbf{z})|   \leq  D_3|D(\mathbf{z})|$,
where
\begin{align*}
D_3 & =    \frac{C_0}{\nu_0(\theta-U\overline C_P)} 
\Big(  C_1C (|A_1\overline m_1| +  |\nabla E|_4) |D(\mathbf{u}_2)|
 +  \frac{\overline C_P C_P}{2} (|A_1 \overline m_1| + |\nabla E|_4)\Big)
\\
&\quad  + \frac{C_0}{2 \nu_0} |D(\mathbf{u}_1)|.
\end{align*}
When $D_3 < 1$, we conclude that $|D(\mathbf{z})| = 0$ and, from \eqref{e35}, 
that $|A_1 \varphi| = 0$.
Since $\mathbf{z} \in \mathbf{J}_0(\Omega)$ and $\varphi \in H^2(\Omega) \cap Y$,
 its follows that $\mathbf{z} = 0$, $\varphi = 0$ in $\Omega$; consequently, 
$\mathbf{u}_1 = \mathbf{u}_2$ and $\overline m_1 = \overline m_2$.
\end{proof}

\subsection*{Acknowledgments}
J. L. Boldrini was partially supported by grant 305467/2011-5 from
CNPq (Brazil),  grants 2007/06638-6 and 2009/15098-0
from FAPESP (Brazil), and grant 1120260 from
Fondecyt (Chile). 

M. A. Rojas-Medar was partially supported by grant MTM2012-32325
from DGI-MEC (Spain),
 and grants 1120260, GI/C-UBB 121909 from Fondecyt-Chile.


\begin{thebibliography}{99}

\bibitem{adams} R. A. Adams;
\emph{Sobolev spaces}, Academic Press, New York, 1975.

\bibitem{brezis} H. Brezis;
\emph{Analyse fonctionnelle-th\'eorie et
Applications}, 2 tirage, Masson, 1987.

\bibitem{Climent}  B. Climent-Ezquerra, L. Friz, M. A. Rojas-Medar;
\emph{Time-reproductive solutions for a bioconvective flow}, 
Annali di Matematica Pura ed Applicata, In Press. DOI 10.1007/s10231-011-0245-7

\bibitem{kan} Y. Kan-On, K. Narukawa, Y. Teramoto;
\emph{On the equations of bioconvective flow},
J. Math. Kyoto Univ. (JMKYAZ), Vol.32, N$^{\mathbf{b}{0}}$ 1, 1992.

\bibitem{levandowsky} M. Levandowsky, W. S. Childress, S. H. Hunter,
 E. A. Spiegel;
\emph{A mathematical model of pattern formation by swimming
microorganisms}, J. Protozoology, Vol.22, 1975.

\bibitem{lorca} S. Lorca, J. L. Boldrini;
\emph{Stationary solutions for generalized Boussinesq models},
J. Differential Equations, Vol. 124, N. 2, 1996.

\bibitem{lorca1} S. Lorca, J. L. Boldrini;
\emph{The initial value problem for a generalized Boussinesq model},
Nonlinear Analysis: Theory, Methods and Applications,
Volume 36, Issue 4, 1999

\bibitem{moribe} Y. Moribe;
\emph{On the bioconvection of tetrahymena pyriformis},
Master's Thesis, Osaka University, 1973.

\bibitem{rionero} S. Rionero, G. Mulone;
\emph{Existence and Uniqueness Theorems for a Steady
Thermo-Diffusive Mixture in a Mixed Problem}, Nonlinear Analysis,
TMA, col.12, n 5, 1988.

\bibitem{solon} V. Solonnikov, V. Scadilov;
\emph{On a boundary value problem for a stationary system
of Navier-Stokes equations}, Proc. Steklov Inst. Math., Vol.125,
1973.

\bibitem{temam} R. Temam;
\emph{Navier-Stokes equations-Theory and numerical analysis}
(revised edition), North-Holland Publ. Comp., Amsterdam, 1979.

\end{thebibliography}

\end{document}