A spectral Navier-Stokes equation solver has been developed for the accurate prediction of two-dimensional buoyancy-driven convection in confined two-layer systems. In particular, the analysis is focussed on fluids with larger density disparity (heavier fluid is at the bottom), wherein the corresponding interfacial deformation can be neglected. The coupled evolution of flow in both fluid layers is handled through a domain decomposition method involving the influence matrix approach to incorporate the interfacial boundary conditions. The flow variables in the individual domains are transiently evaluated by expanding them in terms of Lagrangian interpolation polynomials defined over the Gauss-Lobatto-Chebyshev points. Using these approaches, the occurrence of mechanical (viscous) and thermal coupling modes has been examined at various interfacial heights in the confined two-layer system. The presence of oscillatory switching (standing wave mode) between the thermal and mechanical coupling modes has also been identified when the fluid layers are of nearly equal height, for a chosen combination of fluid properties. Additionally, benchmark results of spectral accuracy are listed for the sake of comparison with the results of future studies on two-layer convection. © 2014 Taylor and Francis Group, LLC.