Operational modal analysis is the primary tool for modal parameter identification in civil engineering. Bayesian statistics offers an ideal framework for analyzing uncertainties associated with the identified modal parameters. However, the exact Bayesian formulation is usually intractable due to the high computational demand in obtaining the posterior distributions of modal parameters. In this paper, the variational Bayes method is employed to provide an approximate solution. Unlike the Laplace approximation and Monte Carlo sampling, the variational Bayes approach provides a gradient-free algorithm to analytically approximate the posterior distributions. Working with the state-space representation of a dynamical system, the variational Bayes approach for identification of modal parameters is derived by ignoring statistical correlation between latent variables and the model parameters. In this approach, the joint distribution of the state-transition and observation matrices as well as the joint distribution of the process noise and measurement error are firstly calculated analytically using conjugate priors. The distribution of modal parameters is extracted from these obtained joint distributions using a first-order Taylor series expansion. A robust implementation of the method is discussed by using square-root filtering and Cholesky decomposition. The proposed approach is illustrated by its application to an example mass-spring system and the One Rincon Hill Tower in San Francisco.