We analyse the entrainment and energetics of equal and opposite axisymmetric tur-bulent air plumes in a vertically confined space at a Rayleigh number of1.24×107using theory and direct numerical simulation. On domains of sufficiently large aspectratio, the steady-state consists of turbulent plumes penetrating an interface betweentwo layers of approximately uniform buoyancy. As described by Baines & Turner (J.Fluid Mech.vol. 37, 1969, pp. 51-80), upon penetrating the interface the flow in eachplume becomes forced and behaves like a constant-momentum jet, due to a reduction inits mean buoyancy relative to the local environment. To observe the behaviour of theplumes we partition the domain into sub-domains corresponding to each plume. Domainsof relatively small aspect ratio produce a single primary mean-flow circulation betweenthe sub-domains that is maintained by entrainment into the plumes. At larger aspectratios the mean flow between the sub-domains bifurcates, indicating the existence of asecondary circulation within each layer associated with entrainment into the jets. Thelargest aspect ratios studied here exhibit an additional, tertiary, circulation in the vicinityof the interface. Consistency between independent calculations of an effective entrainmentcoefficient allows us to identify aspect ratios for which the flow can be modelled usingplume theory, under the assumption of a two-layer stratification.To study the flow’s energetics we use a local definition of available potential energy(APE). For plumes with Gaussian velocity and buoyancy profiles, the theory we developsuggests that the kinetic energy dissipation is split equally between the jets and theplumes and, collectively, accounts for almost half of the input of APE at the boundaries.In contrast,1/4of the APE dissipation and background potential energy (BPE) pro-duction occurs in the jets, with the remaining3/4occurring in the plumes. These bulktheoretical predictions agree with observations of BPE production from simulations towithin1%and form the basis of a similarity solution that models the vertical dependenceof APE dissipation and BPE production. Unlike results concerning the dissipation ofbuoyancy variance and the strength of the circulations described above, the model forthe flow’s energetics does not involve an entrainment coefficient.