Abstract: Numerical ocean models have become powerful tools for providing a realistic view of the ocean state and for describing ocean processes that are difficult to observe. Recent improvements in model performance focus on simulating realistic ocean interior mixing rates, as ocean mixing is the main physical process that creates water masses and maintains their properties. Below the mixed layer, diapycnal mixing primarily arises from the breaking of internal waves, whose energy is largely supplied by winds and tides. This is particularly true in abyssal regions, where the barotropic tide interacts with rough topography and where high levels of diapycnal mixing have been recorded (e.g., the Hawaiian Archipelago). Many studies have discussed the representation of internal wave generation, propagation, and evolution in ocean numerical models. Expanding on these studies, this work seeks to better understand the representation of internal wave dynamics, energetics, and their associated mixing in several different classes of widely used ocean models (e.g., the HYbrid Coordinate Ocean Model, HYCOM; the Regional Ocean Modeling System, ROMS; and the MIT general circulation model, MITgcm). First, a multi-model study investigates the representation of internal waves for a wide spectrum of numerical choices, such as the horizontal and vertical resolution, the vertical coordinate, and the choice of the numerical advection scheme. Idealized configurations are compared to their corresponding analytical solutions. Some preliminary results of realistic baroclinic tidal simulations are shown for the Gulf of Mexico. Second, the spurious diapycnal mixing that exists in models with fixed vertical coordinates (i.e., geopotential and terrain following) is documented and quantified. This purely numerical error arises because, in fixed-coordinate models, the numerical framework cannot properly maintain the adiabatic properties of an advected water parcel. This unrealistic mixing of water masses can be a source of major error in both regional and global ocean models. We use the tracer flux method to compute the spurious diapycnal diffusivities for both a lockexchange scenario and a propagating internal wave field using all three models. Results for the lock exchange experiments are compared to the results of a recent study. Our results, obtained by using three different model classes, provide a comprehensive analysis of the impact of model resolution choice and numerical framework on the magnitude of the spurious diapycnal mixing and the representation of internal waves.