Nowadays, multivariate functional data are frequently observed in many scientific fields, and the estimation of quantiles of these data is essential in data analysis. Unlike in the univariate setting, quantiles are more challenging to estimate for multivariate data, let alone multivariate functional data. This article proposes a new method to estimate the quantiles for multivariate functional data with application to air pollution data. The proposed multivariate functional quantile model is a nonparametric, time-varying coefficient model, and basis functions are used for the estimation and prediction. The estimated quantile contours can account for non-Gaussian and even nonconvex features of the multivariate distributions marginally, and the estimated multivariate quantile function is a continuous function of time for a fixed quantile level. Computationally, the proposed method is shown to be efficient for both bivariate and trivariate functional data. The monotonicity, uniqueness, and consistency of the estimated multivariate quantile function have been established. The proposed method was demonstrated on bivariate and trivariate functional data in the simulation studies, and was applied to study the joint distribution of PM 2.5 and geopotential height over time in the Northeastern United States; the estimated contours highlight the nonconvex features of the joint distribution, and the functional quantile curves capture the dynamic change across time.