Multivariate quantile modeling and estimation for functional data is challenging. Unlike in the univariate setting, multivariate quantile is not uniquely defined, which makes the extension to the multivariate functional setting even more challenging. In this paper, we propose a flexible quantile model for multivariate functional data and develop the estimation and prediction procedures. Marginally, the multivariate contour estimation can characterize nonGaussian and even nonconvex multivariate distributions, and cubic spline models are adopted for estimating the multivariate quantile curves. We establish the consistency of the quantile contour estimation using the uniform convergence property. For numerical studies, we simulate bivariate functional data with nonconvex marginal distribution and illustrate the estimation of quantile contours and functional quantile curves. For a bivariate normal marginal distribution, we show that the estimated contours coincide with the density contours. Finally, we use the proposed method to analyze air pollution data, where we study the joint distribution of particulate matter (PM2.5) and geopotential height at pressure level 850 Pa along time.