Let the portfolio be given by: $$X=X_1+X_2$$ $(X_1,X_2)$ are dependent through a Copula function $C(u_1,u_2)$, such that the joint distribution is given by: $$F(x_1,x_2)=C(F(x_1),F(x_2))$$
What is the VaR of this portfolio?
Usually VaR is the inverse quantile: $VaR_\alpha=F^{-1}(\alpha)$.
I am not sure how to determine it in this multivariate case?
It may not be possible to determine VaR analytically and you will need to run a Monte Carlo simulation.
As you indicate, for a univariate distribution $F$, VaR is specified as the $\alpha$-percentile:
$$P(X \leq \text{VaR}) = \alpha \implies \text{VaR}=F^{-1}(\alpha).$$
In the multivariate case, with $X = X_1 + X_2$, you need the distribution $F_{X_1+X_2}(x).$ This may not be possible in closed form for a general copula. However, in the case of a Gaussian copula with Gaussian marginals, you can calculate VaR analytically because the sum of normally distributed random variables is also normally distributed.
As an example, consider the determination of VaR using a Student's t copula and a correlation matrix $\Sigma$. You can simulate and sample from a joint distribution, specified by the copula, as follows:
(1) Using a standard random number procedure, generate a vector with a correlated multivariate normal distribution: $\mathbf{Z} \sim N(0,\Sigma)$, and an independent random variable with a chi-squared distribution with $\nu$ degrees of freedom $\chi^2(\nu).$
(2) Let $S = \sqrt{\chi^2(\nu)/\nu}$. Then $\mathbf{X} = S^{-1}\mathbf{Z}$ will have a multivariate Student's t distribution.
(3) For each component $x_i$ of $\mathbf{X}$, calculate $y_i = N^{-1}[F_{t_\nu}(x_i)]$ where $N^{-1}$ is the inverse standard normal distribution function and $F_{t_\nu}$ is the distribution function for a $\chi$-squared distribution with $\nu$ degrees of freedom.
At this point, the vector $\mathbf{Y}$ is a sample from a joint distribution specified by a Student's t copula. The marginal distributions for the components are standard normal. The degree of dependence and kurtosis is controlled through the parameters $\Sigma$ and $\nu$.
To estimate VaR for a portfolio $P = y_1 + \ldots + y_n$, generate many sample vectors $\mathbf{Y}^{(1)}, \mathbf{Y}^{(2)},...$ and find the corresponding portfolio values $P^{(1)}, P^{(2)}, ...$ Sort the values and find the desired quantile estimate.
The actual copula is
$$C_{\nu,\Sigma} (x_1,\ldots,x_n) = F_{\nu,\Sigma}[F_{t_\nu}^{-1}(x_1),\ldots,F_{t_\nu}^{-1}(x_n)],$$
where $F_{\nu, \Sigma}$ is the multivariate Student's t distribution function.
You can use the copula to impose other marginal distributions with the same dependence structure by altering step (3):
$$y_i = \hat{F}^{-1}[F_{t_\nu}(x_i)], $$
where $\hat{F}$ is the desired marginal distribution function.