A new implicit method for obtaining equilibrium solutions and their sensitivity to changes in parameters is described and applied to an OCMIP-2 type ocean-biogeochemistry model. The method is used to optimize model parameters by minimizing the difference between the observed and simulated PO4 distribution. The optimized parameters include (1) the exponent α in the power law vertical profile for particulate organic matter (POM) fluxes, (2) the fraction σ of biological production allocated to dissolved organic matter (DOM) and (3) the rate constant $kappa$ for the remineralization of DOM. Global PO4 observations constrain σ and $kappa$ but not independently because their sensitivity patterns are highly correlated. In contrast, the sensitivity pattern for α is uncorrelated to those of the other parameters, allowing it to be independently constrained. We show that export production from POC is well constrained by the distribution of PO4 in an OCMIP-2 type model, but that new production and export production from DOC are not well constrained. With the optimal parameter set (α = -1.0, σ = 0.74, and $kappa$ = 1.0 yrs-1) the fraction of the spatial PO4 variance captured by our model increases from 60% with the reference OCMIP-2 parameters to 70%. Combined changes in σ and $kappa$ account for most of the improvements by reducing but not completely eliminating the nutrient trapping effect in the Eastern Equatorial Pacific and northern Indian Ocean that causes the model to over-predict PO4 concentrations. Important remaining model-data misfits in the deep North Atlantic where PO4 is over predicted and in the North Pacific where the model does not produce the observed sharp nutricline are likely attributable to deficiencies in ocean transport. The fact that the fraction of unexplained variance is large at the optimal parameter values highlights the importance of properly simulating physical transport for ocean biogeochemical modeling.