  In the present work we propose and analyze a fully coupled virtual element
method of high order for solving the two dimensional nonstationary Boussinesq
system in terms of the stream-function and temperature fields. The
discretization for the spatial variables is based on the coupling $C^1$- and
$C^0$-conforming virtual element approaches, while a backward Euler scheme is
employed for the temporal variable. Well-posedness and unconditional stability
of the fully-discrete problem is provided. Moreover, error estimates in $H^2$-
and $H^1$-norms are derived for the stream-function and temperature,
respectively. Finally, a set of benchmark tests are reported to confirm the
theoretical error bounds and illustrate the behavior of the fully-discrete
scheme.
