Normative modeling is a promising statistical framework in clinical neuroscience that characterizes individual deviations from population-based reference distributions. While traditional approaches focus on univariate modeling of individual brain measures, multivariate normative modeling using deep generative models, particularly Variational Autoencoders (VAE), has recently emerged as a powerful alternative. These models capture high-dimensional dependencies across brain features and enable the detection of subtle deviations that are difficult to observe with univariate methods. However, current multivariate approaches lack systematic evaluation of covariate modeling methods and remain underexplored in handling batch effects and clinical applicability. For this work, an experimental platform is developed to train and evaluate VAE-based multivariate normative models. The models are trained on structural MRI data from the Generation R Study and the Healthy Brain Network, incorporating key covariates such as age, sex, and acquisition site. A wide set of covariate modeling methods is systematically evaluated in terms of reconstruction quality, latent space covariate invariance, and alignment with normative priors. This work also investigates whether VAE-based multivariate normative models can accommodate batch effects, such as site variation, and compares them to traditional data harmonization techniques, like ComBat. Finally, the model is applied to the clinically relevant task of brain age estimation. The results show that incorporating covariate modeling into the VAE architecture can significantly improve covariate invariance. When examining the influence of batch effects, covariate modeling methods and ComBat data harmonization both reduce site-related information in low-dimensional latent spaces. However, when the latent dimensionality increases, ComBat data harmonization outperforms all covariate modeling methods. In a proof-of-concept application, the model was successfully extended for brain age estimation, capturing age-related deviations while preserving an age-invariant latent space. Altogether, these findings show the potential of VAE-based multivariate normative models for clinical neuroscience applications.