  Pointwise divergence free velocity field approximations of the Stokes system
are gaining popularity due to their necessity in precise modelling of physical
flow phenomena. Several methods have been designed to satisfy this requirement;
however, these typically come at a greater cost when compared with standard
conforming methods, for example, because of the complex implementation and
development of specialized finite element bases. Motivated by the desire to
mitigate these issues for 2D simulations, we present a $C^0$-interior penalty
Galerkin (IPG) discretization of the Stokes system in the stream function
formulation. In order to preserve a spatially varying viscosity this approach
does not yield the standard and well known biharmonic problem. We further
employ the so-called robust interior penalty Galerkin (RIPG) method; stability
and convergence analysis of the proposed scheme is undertaken. The former,
which involves deriving a bound on the interior penalty parameter is
particularly useful to address the $\mathcal{O}(h^{-4})$ growth in the
condition number of the discretized operator. Numerical experiments confirming
the optimal convergence of the proposed method are undertaken. Comparisons with
thermally driven buoyancy mantle convection model benchmarks are presented.

Numerical Analysis;Computational Engineering, Finance, and Science;Numerical Analysis;Geophysics