A statistical turbulence model is proposed for ensemble calculations with two fluids coupled across a flat interface, motivated by atmosphere-ocean interaction. For applications, like climate research, the response of an equilibrium climate state to variations in forcings is important to interrogate predictive capabilities of simulations. The method proposed here focuses on the computation of the ensemble mean-flow fluid velocities. In particular, a closure model is used for the Reynolds stresses that accounts for the fluid behavior at the interface. The model is shown to converge at long times to statistical equilibrium and an analogous, discrete result is shown for two numerical methods. Some matrix assembly costs are reduced with this approach. Computations are performed with monolithic (implicit) and partitioned coupling of the fluid velocities; the former being too expensive for practical computing, but providing a point of comparison to see the effect of partitioning on the ensemble statistics. It is observed that the partitioned methods reproduce the mean-flow behavior well, but may introduce some long-time statistical bias.