Longitudinal metabolomics holds great promise in understanding the causes, characteristics, and potential interventions for diseases because of the specificity of small-molecule snapshots, temporal ordering of events, and repeated sampling from individuals. There are many challenges with this data, including non-independent observations, large noise components, small sample sizes, abundant structural zeros, nonlinearity, and high dimensionality. In this work we investigate the use of Gaussian processes, a Bayesian nonparametric method to uncover and categorize complex relationships, including temporal, in an automated fashion through kernel search. Simulations show that with sufficient sample size, models with complex temporal dynamics can be recovered, while reverting to uncertainty in higher noise settings thus guarding against overfitting. This approach is then used to fit a two-stage model of metabolomics data of individuals diagnosed with Crohn’s disease from the Integrative Human Microbiome Project to investigate associations between metabolite values and clinical variables of interest.