SUMMARY Significant compositional differences may exist in the lithospheric mantle and above the core–mantle boundary (CMB) relative to the ambient mantle. The intrinsic density differences may affect the development of thermal boundary layer (TBL) instabilities associated with lithospheric delamination and formation of thermochemical plumes. In this study, we explored the instability of two-layer thermochemical fluid using two different techniques: marginal stability analysis with a propagator-matrix method and finite element modelling. We investigated both the instabilities in lithospheric mantle (i.e. lithospheric instability) and the mantle above the CMB (i.e. plume-forming instability) using a background temperature Tbg(z) with the TBL. For lithospheric instability, we found that two-layer fluid with free-slip boundary conditions mainly undergoes the same three different convective modes (i.e. two oscillatory convection modes and one layered convection regime) as that with no-slip boundary conditions reported in Jaupart et al. However, with free-slip boundary conditions, the transitions between these convection modes occur at larger values of buoyancy number B. Free-slip boundary conditions lead to smaller critical Rayleigh number Rac, but larger convective wavelength and oscillation frequency ωc, compared with those with no-slip boundary conditions. Our numerical modelling results demonstrate that Rac and ωc predicted from the classical marginal stability analyses using Tbg(z) with TBL temperature may have significant errors when the oscillatory period is comparable with or larger than the timescale of lithospheric thermal diffusion that causes Tbg(z) to vary with time significantly. In this case, using a more gently sloped background temperature profile ignoring the TBL temperature, the stability analysis predicts more accurate stability conditions, thus presenting an effective remedy to the stability analysis. For plume-forming instability, because of the reduced viscosity in the hot and compositionally dense bottom layer, the transition to the layered convection occurs at significantly smaller B values, and in the oscillatory convection regime, Rac is larger but ωc is smaller, compared with those for lithospheric instability. Finally, our study provides a successful benchmark of numerical models of thermochemical convection by comparing Rac and ωc from numerical models with those from the marginal stability analysis.