Summary: A two-dimensional quadrilateral spectral multidomain penalty method (SMPM) model has been developed for the simulation of high Reynolds number incompressible stratified flows. The implementation of higher-order quadrilateral subdomains renders this model a nontrivial extension of a one-dimensional subdomain SMPM model built for the simulation of the same type of flows in vertically nonperiodic domains [P. J. Diamessis et al., J. Comput. Phys. 202, No. 1, 298–322 (2005; Zbl 1061.76054)]. The nontrivial aspects of this extension consist of the implementation of subdomain corners, the penalty formulation of the pressure Poisson equation (PPE), and, most importantly, the treatment of specific challenges that arise in the iterative solution of the SMPM-discretized PPE. The two primary challenges within the framework of the iterative solution of the PPE are its regularization to ensure the consistency of the associated linear system of equations and the design of an appropriate two-level preconditioner. A qualitative and quantitative assessment of the accuracy, efficiency, and stability of the quadrilateral SMPM solver is provided through its application to the standard benchmarks of the Taylor vortex, lid-driven cavity, and double shear layer. The capacity of the flow solver for the study of environmental stratified flow processes is shown through the simulation of long-distance propagation of an internal solitary wave of depression in a manner that is free of numerical dispersion and dissipation. The methods and results presented in this paper make it a point of reference for future studies oriented toward the reliable application of the quadrilateral SMPM model to more complex environmental stratified flow process studies.