If we want to satisfy the following constraint with the probability of 1-$\epsilon$:
$$ \mathcal{P}(a \leq x) \ge 1 - \epsilon $$ where $x$ is the decision variable, $a$ is a random variable with a normal distribution $ N(\mu,\sigma^2)$, it is enough to have $$ x \geq \mu + \sigma \Phi^{-1}(1-\epsilon) $$
My question is what we can say about the same problem in a higher dimension, for example 2, where $a_1$ and $a_2$ are correlated.
$$ \mathcal{P}(a_1 \leq x_1 , a_2 \leq x_2 ) \ge 1 - \epsilon $$
Partial answer (too long for comment). You can evaluate the bivariate normal CDF numerically (using e.g. Python and Scipy, which is what I did). We see that, even in the correlated case, the feasible region is convex, which is a hopeful sign, computationally (basically the CDF of the bivariate normal distribution is quasi-convex).
For these numerical examples, I took $\mu_1=\mu_2=0$, $\sigma_1=0.2$, $\sigma_2=0.4$ and covariance matrix $$ \Sigma=\begin{bmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2\\\rho\sigma_1\sigma_2 & \sigma_2^2 \end{bmatrix} $$ where $\rho\in[0,1]$ is the correlation parameter. Here's $\rho=0.5$, with $\varepsilon=0.4$.
The simplest naive approximation to this region would be $$ x_i\geq\mu_i+\Phi^{-1}(1-\varepsilon)\sigma_i. $$ In general, this approximation isn't great:
But, counter-intuitively (to me), this approximation seems to be exact in the limit $\rho\to1$. Here's $\rho=0.99$: