thermal diffusivity and gravity terms would also have to be treated implicitly for the time discretization methods to be efficient. Straightforward application of any monolithic solver
would then yield large algebraic systems with multiple couplings between discrete DOF associated to different physical variables. To avoid this, we resort to an operator splitting approach, see e.g. [1]},
as commonly done in numerical models for atmospheric physics, see e.g. the reviews in [2]}, [3]}. More specifically, after spatial discretization, all diffusive terms on the right hand side of (REF ) are split from the hyperbolic part on the left hand side.
The hyperbolic part is treated in a similar fashion to what outlined in [4]}, while the diffusive terms are treated implicitly. For simplicity, the gravity terms will be treated explictly in this first attempt and only a basic, first order splitting will be described, which can be easily improved to second order accuracy by the Strang splitting approach [1]}, [6]}.
