Propensity scores are commonly employed in observational study settings where the goal is to estimate average treatment effects. This paper introduces a flexible propensity score modeling approach, where the probability of treatment is modeled through a Gaussian process framework. To evaluate the effectiveness of the estimated propensity score, a metric of covariate imbalance is developed that quantifies the discrepancy between the distributions of covariates in the treated and control groups. It is demonstrated that this metric is ultimately a function of the hyperparameters of the covariance matrix of the Gaussian process and therefore it is possible to select the hyperparameters to optimize the metric and minimize overall covariate imbalance. The effectiveness of the GP method is compared in a simulation against other methods of estimating the propensity score and the method is applied to data from Dehejia and Wahba (1999) to demonstrate benchmark performance within a relevant policy application.