We present a statistical method for reconstructing continuous background density profiles that embeds incomplete measurements and a physically intuitive density stratification model within a Bayesian hierarchal framework. A double hyperbolic tangent function is used as a parametric density stratification model that captures various pycnocline structures in the upper ocean and offers insight into several density profile characteristics (e.g., pycnocline depth). The posterior distribution is used to quantify uncertainty and is estimated using recent advances in Markov chain Monte Carlo sampling. Temporally evolving posterior distributions of density profile characteristics, isopycnal heights, and nonlinear ocean process models for internal gravity waves are presented as examples of how uncertainty propagates through models dependent on the density stratification. The results show 0.95 posterior interval widths that ranged from 2.5% to 4% of the expected values for the linear internal wave phase speed and 15%-40% for the nonlinear internal wave steepening parameter. The data, collected over a year from a through-the-column mooring, and code, implemented in the software package Stan, accompany the article.