We numerically investigate turbulent Rayleigh–Bénard convection through two immiscible fluid layers, aiming to understand how the layer thickness and fluid properties affect the heat transfer (characterized by the Nusselt number
$\mbox {Nu}$
) in two-layer systems. Both two- and three-dimensional simulations are performed at fixed global Rayleigh number
$\mbox {Ra}=10^8$
, Prandtl number
$\mbox {Pr}=4.38$
and Weber number
$\mbox {We}=5$
. We vary the relative thickness of the upper layer between
$0.01 \le \alpha \le 0.99$
and the thermal conductivity coefficient ratio of the two liquids between
$0.1 \le \lambda _k \le 10$
. Two flow regimes are observed. In the first regime at
$0.04\le \alpha \le 0.96$
, convective flows appear in both layers and
$\mbox {Nu}$
is not sensitive to
$\alpha$
. In the second regime at
$\alpha \le 0.02$
or
$\alpha \ge 0.98$
, convective flow only exists in the thicker layer, while the thinner one is dominated by pure conduction. In this regime,
$\mbox {Nu}$
is sensitive to
$\alpha$
. To predict
$\mbox {Nu}$
in the system in which the two layers are separated by a unique interface, we apply the Grossmann–Lohse theory for both individual layers and impose heat flux conservation at the interface. Without introducing any free parameter, the predictions for
$\mbox {Nu}$
and for the temperature at the interface agree well with our numerical results and previous experimental data.