Accurate modeling of physical and biogeochemical dynamics in coastal ocean regions is required for multiple scientific and societal applications, covering a wide range of time and space scales. However, in light of the strong nonlinearities observed in coastal regions and in biological processes, such modeling is challenging. An important subject that has been largely overlooked is the numerical requirements for regional ocean simulation studies. Major objectives of this thesis are to address such computational questions for non-hydrostatic multiscale flows and for biogeochemical interactions, and to derive and develop numerical schemes that meet these requirements, utilizing the latest advances in computational fluid dynamics. We are interested in studying nonlinear, transient, and multiscale ocean dynamics over complex geometries with steep bathymetry and intricate coastlines, from sub-mesoscales to basin-scales. These dynamical interests, when combined with our requirements for accurate, efficient and flexible ocean modeling, led us to develop new variable resolution, higher-order and non-hydrostatic ocean modeling schemes. Specifically, we derived, developed and applied new numerical schemes based on the novel hybrid discontinuous Galerkin (HDG) method in combination with projection methods. The new numerical schemes are first derived for the Navier-Stokes equations. To ensure mass conservation, we define numerical fluxes that are consistent with the discrete divergence equation. To improve stability and accuracy, we derive a consistent HDG stability parameter for the pressure-correction equation. We also apply a new boundary condition for the pressure-corrector, and show the form and origin of the projection method’s time-splitting error for a case with implicit diffusion and explicit advection. Our scheme is implemented for arbitrary, mixed-element unstructured grids using a novel quadrature-free integration method for a nodal basis, which is consistent with the HDG method. To prevent numerical oscillations, we design a selective high-order nodal limiter. We demonstrate the correctness of our new schemes using a tracer advection benchmark, a manufactured solution for the steady diffusion and stokes equations, and the 2D lock-exchange problem. These numerical schemes are then extended for non-hydrostatic, free-surface, variable-density regional ocean dynamics. The time-splitting procedure using projection methods is derived for non-hydrostatic or hydrostatic, and nonlinear free-surface or rigid-lid, versions of the model. We also derive consistent HDG stability parameters for the free-surface and non-hydrostatic pressure-corrector equations to ensure stability and accuracy. New boundary conditions for the free-surface-corrector and pressure-corrector are also introduced. We prove that these conditions lead to consistent boundary conditions for the free-surface and pressure proper. To ensure discrete mass conservation with a moving free-surface, we use an arbitrary LagrangianEulerian (ALE) moving mesh algorithm. These schemes are again verified, this time using a tidal flow problem with analytical solutions and a 3D lock-exchange benchmark. We apply our new numerical schemes to evaluate the numerical requirements of the coupled biological-physical dynamics. We find that higher-order schemes are more accurate at the same efficiency compared to lower-order (e.g. second-order) accurate schemes when modeling a biological patch. Due to decreased numerical dissipation, the higher-order schemes are capable of modeling biological patchiness over a sustained duration, while the lower-order schemes can lose significant biomass after a few non-dimensional times and can thus solve erroneous nonlinear dynamics. Finally, inspired by Stellwagen Bank in Massachusetts Bay, we study the effect of non-hydrostatic physics on biological productivity and phytoplankton fields for tidally-driven flows over an idealized bank. We find that the non-hydrostatic pressure and flows are important for biological dynamics, especially when flows are supercritical. That is, when the slope of the topography is larger than the slope of internal wave rays at the tidal frequency. The non-hydrostatic effects increase with increasing nonlinearity, both when the internal Froude number and criticality parameter increase. Even in cases where the instantaneous biological productivity is not largely modified, we find that the total biomass, spatial variability and patchiness of phytoplankton can be significantly altered by non-hydrostatic processes. Our ultimate dynamics motivation is to allow quantitative simulation studies of fundamental nonlinear biological-physical dynamics in coastal regions with complex bathymetric features such as straits, sills, ridges and shelfbreaks. This thesis develops the necessary numerical schemes that meet the stringent accuracy requirements for these types of flows and dynamics.