Here we briefly outline the analysis used in proving that our model has a solution.

Theorem

There exists a unique $(\mathbf{u},\theta,q_v,q_c)$ such that:

$$\mathbf{u}\in L^2(0,\tau;V)\cap L^{\infty}(0,\tau;H),$$

$$\theta \in L^2(0,\tau;H^1_{q}(\Omega))\cap L^{\infty}(\Omega \times (0,\tau)),$$

and

$$q_v,q_c \in L^{\infty}(0,\tau;L^2(\Omega)),$$

$$\mathbf{u}_t \in L^2(0,\tau;V^*)$$

and

$$\mathbf{u}(\cdot,0)=\mathbf{u}_0 \text{ a.e. in } \Omega \times (0,\tau) \text{ and } \forall w \in V,$$

$$\langle\mathbf{u}_t,w\rangle +a(\mathbf{u},w)+b(\mathbf{u},\mathbf{u},w)=\langle \mathbf{B},w\rangle$$

$$ \theta_t\in L^2(0,\tau;H^{-1}(\Omega)) \text{ and } \forall \xi \in L^2(0,\tau;H^1_q(\Omega))$$

$$\langle\theta,\xi \rangle +a(\theta,\xi)+b(\mathbf{u},\theta,\xi)=\langle QC_c+u_n,\xi\rangle$$

$$q_v(0)=q_{v0}, \text{ }q_c(0)=q_{c0} \text{ a.e. in } {\Omega} \text{ and } \forall \eta \in L^2(\Omega)$$

$$\langle q_{ct},\eta\rangle=-b(\mathbf{u},q_c,\eta)+\langle C_c,\eta \rangle $$

$$\langle q_{vt},\eta\rangle=-b(\mathbf{u},q_c,\eta)-\langle C_c,\eta\rangle.$$

We begin by proving some a priori estimates on the quantities in the model, assuming sufficient regularity of the initial conditions. After this we decouple the system, associating to each $n \in \mathbb{N}$ the functions $(\mathbf{u}_n,T_n, q_{vn}, q_{cn})$ and show that each solves its own equation for $n$, given that the others have a solution for $n-1$.

To achieve this we use a Galerkin method to show existence of an approximate solution on a finite dimensional subspace of our Hilbert space, we then combine this with our a priori estimates to ensure that these approximations converge to a solution of the full equation. Once we have this we let $n \to \infty$ and show we have a weak solution of the full model. One of the main difficulties in this was in the term $q_v$. After considering the characteristic $\mathbf{X}$ where
\[
\dot{\mathbf{X}} = \mathbf{u}(\mathbf{X},t)
\]
and representing $q_v$ in terms of $\mathbf{X}$ we must show convergence of the term $\mathbf{u}_n(\mathbf{X}_n,t)$ to $\mathbf{u}(\mathbf{X},t)$.

For uniqueness, let us define $$
(\mathbf{u},\theta,q_v,q_c) = (\mathbf{u}_1-\mathbf{u_2},\theta_1-\theta_2,q_{v1}-q_{v2},q_{c1}-q_{c2}).
$$
as difference of two solutions, and $\mathcal{L}: L^2(0,\tau;V)\rightarrow L^2([0,\tau]\times \Omega) $ as
$$
\mathbf{u} \cdot (\mathcal{L} \mathbf{w}:\mathbf{u}_2):=\mathbf{u}_2\cdot(\mathbf{u}\cdot \nabla)\mathbf{w}= \mathbf{u}\cdot (\frac{\partial \mathbf{w}}{\partial x} \cdot \mathbf{u}_2, \frac{\partial \mathbf{w}}{\partial y} \cdot \mathbf{u}_2).
$$
Clearly $\mathcal{L}$ is linear and continuous.
Also define
$$
h_m :=
\begin{cases}
h & h\le m,\\
m & h <m,
\end{cases}
$$
with a constant $m>0$ large enough, and $h:[0,\tau]\times \Omega\rightarrow [0,\infty]$ given by
$$
h := \chi_{[\theta \ne 0]}.
$$

Let us take smooth test functions $\mathbf{w}(\cdot, t)$ and $\theta(\cdot, t)$ such that $\mathbf{w}(\cdot, \tau) = 0$ and $\xi(\cdot,\tau)=0$ in $\Omega$. Integrating by part we have
\begin{equation}\label{theta2}
\begin{aligned}
\int_{\Omega} \theta (\tau)\xi(\tau)
& = \int_0^\tau \int_{\Omega} \mathbf{u}\cdot(\mathbf{w}_t + (\mathbf{u}_1\cdot \nabla )\xi + \mathcal{L} \mathbf{w}:\mathbf{u}_2 -\xi \nabla \theta_2 +\Delta \mathbf{w})\\
&+ \int_0^\tau \int_\Omega \theta(\xi_t+(\mathbf{u}_1 \cdot \nabla) \xi +\mathbf{f}\cdot \mathbf{w}+ h_m \Delta \xi+C\xi)\\
&+ \int_0^\tau \int_{\Omega}(h-h_m)\theta \Delta \xi.
\end{aligned}
\end{equation}

So we seek test functions $\xi$ and $\mathbf{w}$ satisfying
\begin{equation}
\begin{aligned}
\mathbf{w}_t + (\mathbf{u}_1\cdot \nabla )\xi + \mathcal{L} \mathbf{w}:\mathbf{u}_2 -\xi \nabla \theta_2 +\Delta \mathbf{w} = 0 & \quad [0,\tau]\times \Omega\\
\nabla \cdot \mathbf{w} = 0 &\quad [0,\tau]\times \Omega\\
\xi_t+(\mathbf{u}_1 \cdot \nabla) \xi +\mathbf{f}\cdot \mathbf{w}+ h_m \Delta \xi+C\xi =0 &\quad [0,\tau]\times \Omega \\
\mathbf{w} =\mathbf{0}, \,\xi =0 &\quad [0,\tau]\times \partial \Omega\\
\mathbf{w}(\cdot, \tau)=\mathbf{0},\, \xi(\cdot, \tau)= 0 & \quad \text{in } \Omega.
\end{aligned}
\end{equation}

Lemma
There is a unique weak solution with regularity of the test functions and
\begin{equation}
\begin{aligned}
\mathbf{w}\in H^1(0,\tau;H)\cap L^\infty(0,\tau; V) \cap L^2(0,\tau;V^2)\\
\xi \in H^1(0,\tau;L^2)\cap L^\infty(0,\tau;H_0^1) \cap L^2(0,\tau;H^2)\\
\end{aligned}
\end{equation}
and
$$
\int_0^\tau \int_{\Omega}|\Delta \xi |^2 \le C.
$$

The lemma is proved in [Diaz]*. We have $$\int_0^\tau \int_\Omega \theta^2 = \int_0^\tau \int_{\Omega}(h-h_m)\theta \Delta \xi.$$ The dominated convergence theorem gives
$$
\|(h-h_m)\theta\|_{L^2([0,\tau]\times\Omega)} \rightarrow 0 \quad \text{as }m\rightarrow \infty.
$$
Therefore we have
$$
\int_0^\tau\int_{\Omega} \theta^2 =0.
$$

So $\theta_1 = \theta_2$. Finally standard arguments for the Navier-Stokes equations in 2D imply $\mathbf{u}_1 =\mathbf{u}_2$.

* J. I. Diaz & G. Galiano, Existence and uniqueness of solutions of the boussinesq system with nonlinear thermal diffusion, Top. Meth. in Nonlin. Anal. Vol. 11, 1998, 59-82