This paper introduces a scheme for the numerical approximation of a model for two turbulent flows with coupling at an interface. We consider the variational formulation of the coupled model, where the turbulent kinetic energy equation is formulated by transposition. We prove the convergence of the approximation to this formulation for 3D flows for large turbulent viscosities and smooth enough flows, whenever bounded in W 1,p Sobolev norms for p large enough. Under the same assumptions, we show that the limit is a solution of the initial problem. Finally, we give some numerical experiments to enlighten the theoretical work.