This work presents the details behind each step in the development of a framework for two-dimensional quadrilateral discontinuous Spectral Multidomain Penalty Method (SMPM) solvers for environmental flow processes: a shallow water equation (SWE) solver and an incompressible Navier-Stokes equations (NSE) (under the Boussinesq approximation) solver, with additional emphasis given to the associated pressure solver. The potential for environmental flow simulations through spectral methods is very strong since these methods are exponentially accurate, non-dissipative and nondispersive. These characteristics translate into capturing the smallest resolved scales of the flow and the propagation of ocean/lake waves with minimum numerical error. In addition, the element-based capability of the method enables the appropriate resolution of the important scales of the processes being modeled, the localization of specific events, and the treatment of complex boundary conditions and geometries. Finally, the discontinuous character of the method add enhanced stability to the method for highlynonlinear under-resolved simulations, an intrinsic characteristic of environmental flow simulations. In the SWE solver, the SMPM is compared with a nodal discontinuous Galerkin method (DGM), where the equations are solved with an explicit SSP-RK34 method. The comparison is done by applying both methods to a suite of six commonly considered geophysical flow test cases; we also include results for a classical continuous Galerkin (i.e., spectral element) method for comparison. Both the analysis and numerical experiments show that the SMPM and DGM are essentially identical; both methods can be shown to be equivalent for very special choices of quadrature rules and Riemann solvers in the DGM along with special choices in the type of penalty term in the SMPM. In the NSE solver time is discretized with a high-order fractional step projection method, where the non-linear advection and forcing terms are advanced explicitly via a stiffly stable scheme. After that, an implicit solution of a Poisson pressure equation (PPE) is solved in order to introduce the incompressibility constraint. In the final fractional time-step linear viscosity forces are also solved implicitly by means of a modified Helmholtz equation. Stability of the numerical scheme for under-resolved simulations at high Reynolds numbers is ensured through use of penalty techniques, spectral filtering, dealiasing, and strong adaptive interfacial averaging. Special attention is given to the solution of the PPE linear system of equations, where the fundamental building blocks of the PPE solver presented here are a Kronecker (tensor) product-based computation of the left null singular value of the non-symmetric SMPM-discretized Laplacian matrix and a custom-designed two-level preconditioner. Both of these tools are essential towards ensuring existence and uniqueness of the solution of the discrete linear system of equations and enabling its efficient iterative calculation. Accuracy, efficiency, and stability of the multidomain model are assessed through the solution of the Taylor vortex, lid-driven cavity flow and double shear layer. The propagation of a non-linear internal wave of depression type is also presented to assess the potential of the solver for the study of environmental stratified flows. The availability of the quadrilateral SMPM solver allows the numerical investigation of a much broader range of environmental processes, namely those in streamwise,vertical non periodic domains with both horizontal and vertical localization.