Time-varying graph signals are alternative representation of multivariate (or multichannel) signals in which a single time-series is associated with each of the nodes or vertex of a graph. Aided by the graph-theoretic tools, time-varying graph models have the ability to capture the underlying structure of the data associated with multiple nodes of a graph -- a feat that is hard to accomplish using standard signal processing approaches. The aim of this contribution is to propose a method for the decomposition of time-varying graph signals into a set of graph modes. The graph modes can be interpreted in terms of their temporal, spectral and topological characteristics. From the temporal (spectral) viewpoint, the graph modes represent the finite number of oscillatory signal components (output of multiple band-pass filters whose center frequencies and bandwidths are learned in a fully data-driven manner), similar in properties to those obtained from the empirical mode decomposition and related approaches. From the topological perspective, the graph modes quantify the functional connectivity of the graph vertices at multiple scales based on their signal content. In order to estimate the graph modes, a variational optimization formulation is designed that includes necessary temporal, spectral and topological requirements relevant to the graph modes. An efficient method to solve that problem is developed which is based on the alternating direction method of multipliers (ADMM) and the primal-dual optimization approach. Finally, the ability of the method to enable a joint analysis of the temporal and topological characteristics of time-varying graph signals, at multiple frequency bands/scales, is demonstrated on a series of synthetic and real time-varying graph data sets.