The Bayesian Hierarchical Model (BHM) methodology is increasingly applied in large state-space geophysical fluid systems. Recent success with a BHM to build surface vector wind ensembles for ocean forecast applications prompted an ambitious effort to identify parameters of a lower- trophic level ocean ecosystem model. Early results were far less encouraging than for the surface wind applications. The ocean ecosystem process model equations are adapted from a Nutrient- Phytoplankton-Zooplankton-Detritus (NPZD) model that is a specific formulation of classical predator-prey systems with added source and sink terms. The NPZD ocean ecosystem model is severely undetermined in that the equations involve O(20) parameters that are not constrained by data. Seven parameters were treated as random variables in an NPZD BHM applied to the coastal Gulf of Alaska (CGOA) ecosystem given station data from GLOBEC and surface phytoplankton retrievals from SeaWiFS. The NPZD BHM parameters were not identifiable without the addition of mixing terms in the process model to represent seasonal thermocline formation and shoaling of the upper ocean mixed layer. Given considerable guidance from ensemble forward model calculations, posterior estimations for phytoplankton growth rate (VmNO3) and zooplankton grazing rate (ZooGR) are obtained.