Multivariate functional data present theoretical and practical complications which are not found in univariate functional data. One of these is a situation where the component functions of multivariate functional data are positive and are subject to mutual time warping. That is, the component processes exhibit a similar shape but are subject to systematic phase variation across their time domains. We introduce a novel model for multivariate functional data that incorporates such mutual time warping via nonlinear transport functions. This model allows for meaningful interpretation and is well suited to represent functional vector data. The proposed approach combines a random amplitude factor for each component with population based registration across the components of a multivariate functional data vector and also includes a latent population function, which corresponds to a common underlying trajectory as well as subject-specific warping component. We also propose estimators for all components of the model. The proposed approach not only leads to a novel representation for multivariate functional data, but is also useful for downstream analyses such as Fr'echet regression. Rates of convergence are established when curves are fully observed or observed with measurement error. The usefulness of the model, interpretations and practical aspects are illustrated in simulations and with application to multivariate human growth curves as well as multivariate environmental pollution data.