This dissertation encompasses a breadth of topics in the area of functional data analysis where each function is modeled as a Gaussian process within the framework of a Bayesian hierarchical model. As Gaussian processes cannot be worked with directly in this context, a foundational aspect of this work illustrates that using a finite approximation to each process is sufficient to provide good estimates throughout the entire process. More importantly, it is established that using a finite approximation of a bivariate random process within the estimation procedure also results in providing good estimates throughout the entire bivariate process. With this result, the mean and covariance functions associated with a Gaussian process can be considered as random effects within a Bayesian hierarchical model. Inference for both parameters is based upon their posterior distributions which provide not only estimates of these parameters, but also quantifies variation in these parameters. Here we also propose Bayesian hierarchical models for smoothing, functional linear regression, and functional registration. The registration model introduced here is shown to favorably compare with the best registration methods currently available as measured by the Sobolev Least Squares criterion. Within this registration framework, an Adapted Variational Bayes algorithm is introduced to address the computational costs associated with inference in high-dimensional Bayesian models. With multiple examples, both simulated and using real data, it is shown that this algorithm results in registered function estimates that closely agree with corresponding estimates obtained from an MCMC sampling scheme. With this algorithm, functional prediction is considered for the first time in a registration context. The final area of inference for functional data that is proposed for the first time here is a combined registration and factor analysis model. This model is shown to outperform currently available registration methods for data in which the registered functions vary in more than one functional direction. The models presented here are applied to several simulated data sets as well as data from the Berkeley Growth Study, functional sea-surface temperature data, and a juggling data set.