The iteration scheme starts with \(u^0 = 0\) . Note that since \({\rm div} u = 0\) , we also have \(\mathrm {div}(\Psi _\delta (u^{n})) = 0\) ,
for \(t\in \mathbb {R}^+\) . At each step \(n\) , one solves recursively
\(n+1\)  linear equations: first one solves the transport-diffusion
equation (with smooth convective velocity) for the density and
temperature; after \(a^n, \theta ^n\)  is computed, solving the third part
of (REF ) amounts to solving a linear equation
on each strip \(\mathbb {R}^3 \times (m\delta , (m+1)\delta )\) , for \(m= 0, 1, \cdots , n-1\) . So we can get the existence of the mollified
Boussinesq system (REF ) like the Proposition
3.1 in [1]}. Then there is no additional difficulty in the
passage to the limit in the non-linear terms arises in the equation
of the temperature other than those already existing for the
Navier-Stokes equations. Hence, the distributional limit \((a,\theta ,u, \Pi )\)  of a convergent subsequence of \((a^n, \theta ^n, u^n,\Pi ^n)\)  is a weak solution of the Boussinesq system.
