Solvation-related interactions strongly influence a wide range of biomolecular processes. However, both our models and our information for parameterizing those models are imperfect. This talk will describe strategies for quantifying the uncertainty in biomolecular solvation models and optimizing those models to provide the best possible accuracy and performance. The first part of the talk will describe generalized polynomial chaos methods for quantifying solvation energy uncertainty due to conformational noise and errors in atomic charge and radius parameters. The second part of the talk will outline Bayesian methods for addressing model uncertainty through statistical aggregation of predictions from multiple models.