A simple Nutrient-Phytoplankton-Zooplankton (NPZ) model was coupled with the non-hydrostatic Finite-Volume Coastal Ocean Model (FVCOM-NH) to study the impact of high-frequency nonlinear internal waves on plankton dynamics in Massachusetts Bay (MB) during the stratified summer season. The temporal and spatial variability of phytoplankton concentration follows the vertical isopycnal displacement to the lowest order as the waves are generated by the semidiurnal tidal flow over Stellwagen Bank (SB) and propagate westward across MB. The tidally-averaged distribution of phytoplankton is characterized by three distinct zones of low subsurface concentration: (I) the western flank of Stellwagen Bank; (II) the center of Stellwagen Basin; and (III) the upper western flank of Stellwagen Basin. The result of a model dye experiment suggests that these zones are created by the following physical processes which are dominant in each zone: (I) hydraulic jump; (II) strong internal wave-tidal current nonlinear interaction; and (III) energetic internal wave dissipation and subsequent mixing processes. The nonlinear interaction of the internal waves and offshore tidal currents significantly enhances the vertical velocity, and increases wave dissipation, thus causing an onshore transport of phytoplankton in zone II. Although the phytoplankton patchy structure can be produced using the hydrostatic FVCOM, the resulting phytoplankton concentration is overestimated due to the unrealistic intensification of vertical velocity and thus vertical nutrient flux from the deep water. It suggests that non-hydrostatic dynamics should be considered for certain small-scale biological processes that are driven primarily by the physics.