A spectral multidomain penalty method model has been developed for the simulation of high Reynolds number localized stratified turbulence. This is the first time that a penalty method, with a particular focus on subdomain interface treatment, has been used with a multidomain scheme to simulate incompressible flows. The temporal discretization ensures maximum temporal accuracy by combining third order stiffly stable and backward differentiation schemes with a high-order boundary condition for the pressure. In the non-periodic vertical direction, a spectral multidomain discretization is used and its stability for under-resolved simulations at high Reynolds numbers is ensured through use of penalty techniques, spectral filtering and strong adaptive interfacial averaging. The penalty method is implemented in different formulations for both the explicit non-linear term advancement and the implicit treatment of the viscous terms. The multidomain model is validated by comparing results of simulations of the mid-to-late time stratified turbulent wake with non-zero net momentum to the corresponding laboratory data for a towed sphere. The model replicates correctly the characteristic vorticity and internal wave structure of the stratified wake and exhibits robust agreement with experiments in terms of the temporal power laws in the evolution of mean profile characteristic velocity and lengthscales. (C) 2004 Elsevier Inc. All rights reserved.