Your result is true -- and you have tremendous freedom in selecting the underlying distribution of $(\eta_1,\eta_2).$
Your regions consist of three (infinite) angular sectors located at $(-V_1,-V_2).$ Let $F$ be any continuous distribution in the plane. According to a result which I state and prove rigorously below, it is possible to shift and (uniformly) scale $F$ to assign any specified probabilities $p_i$ to their respective regions $\mathcal{R}_i,$ at least with arbitrarily small error (and with no error when all the $p_i$ are positive), assuming only that the $p_i$ are consistent with the axioms of probability: they must be non-negative and sum to unity.
Given $F$ and $(p_1,p_2,p_3),$ assume (with no loss of generality) that $F$ has been shifted and scaled as just described. Let $(\epsilon_0,\eta_1,\eta_2)$ be any 3D random vector for which $(\eta_1,\eta_2)$ has distribution $F.$ By construction, the random vector $(\epsilon_0, \epsilon_0+\eta_1, \epsilon_0+\eta_2)$ satisfies all your requirements.
This figure illustrates the ideas using the notation introduced below.
The boundaries of the three regions of the question are shown at left. For this construction I have chosen $F$ to be the standard Normal distribution (in $\mathbb{R}^2$) and kept $\sigma$ fixed at $1,$ thereby varying only its origin $\mu.$ When $\mu$ follows one of the colored circles at the left, the corresponding probabilities assigned to the three regions are plotted by the same colored curve at the right, which shows a map of the 2-simplex $S_2$ representing all possible such probability triples. It is evident that any probability triple $(p_1,p_2,p_3)$ you might care to assign to the three regions can be realized in this way.
Note that in the illustration we can take $\epsilon_0$ to be standard Normal and independent of the other variables, whence *every one of these solutions is the standard Normal distribution (in $\mathbb{R}^3$) shifted by $(0,\mu_1,\mu_2)$ for a suitable vector $\mu=(mu_1,mu_2)$ determined uniquely by $(p_1,p_2,p_3).$
I recommend finding $\mu$ numerically. I used the pmvnorm in the mvtnorm package for R to compute the sector probabilities -- the rest is then easy (just minimize the squared Euclidean norm $|\Phi(\mu,1)-(p_1,p_2,p_3)|^2:$ see the function f in the code at the end.) This required some work because this package only computes probabilities of rectangles: you first have to transform an angular sector into a rectangle (either a quadrant, half plane, or complement of a quadrant) using an area-preserving transformation and then apply pmvnorm. The details are given by the function psector at the end of this post.
The regions $\mathcal{R}_i$ don't need to be infinite angular sectors in the plane. In order to carry out the following construction, we need only assume the plane has been partitioned into three disjoint measurable regions $\mathcal{R}_i,$ each with nonempty interior, in such a way that between any two regions indexed by $i$ and $j$ there is a continuous path $\gamma_{ij}$ lying wholly in the interior of $\mathcal{R}_i\cup \mathcal{R}_j$ connecting some point in the interior of $\mathcal{R}_i$ to some point in the interior of $\mathcal{R}_j.$ Let's call this a "good" partition. (There exist partitions of the plane that do not have this property!) For future reference, designate three "base points" $x_i\in\mathcal{R}_i$ lying in the interiors.
Let $X$ be any continuous random variable in the plane associated with a distribution function $F(\mathcal A) = \Pr(X\in\mathcal A)$ to any measurable plane set $\mathcal A.$ To find a solution, we are going to vary $F$ by (a) shifting it and (b) uniformly rescaling it. As a matter of notation, then, write the shifted, rescaled version of $F$ as
$$F(\mathcal A;\mu,\sigma) = \Pr\left(\sigma X + \mu\in\mathcal A\right)$$
for $\mu\in\mathbb{R}^2$ and $\sigma\gt 0.$
The reason for the continuity assumption is that for any $\mathcal A,$ it implies the function
$$(\mu,\sigma)\to F(\mathcal{A};\mu,\sigma)$$
is continuous. (Proof: a small change in $\mu$ or $\sigma$ amounts to applying $F$ to a slightly modified version of $\mathcal A.$ When $\mu$ and $\sigma$ change by sufficiently small amounts, the change in the region is so small that--since $F$ is continuous--the change in its probability is small, too: that means the map is continuous.)
Consequently the map
$$\Phi: \mathbb{R}^2\times\mathbb{R}^{+}\to \mathbb{R}^3$$
given by the three probabilities
$$\Phi(\mu,\sigma) = (F(\mathcal{R}_1;\mu,\sigma),F(\mathcal{R}_1;\mu,\sigma),F(\mathcal{R}_1;\mu,\sigma))$$
is continuous, too. Because the $\mathcal{R}_i$ are a partition, the sum of these three probabilities is always $1,$ showing the image of $\Phi$ is a subset of the simplex
$$S_2 = \{(x,y,z)\in\mathbb{R}^3\mid x+y+z+1;\ x\ge 0, y\ge 0, z\ge 0\}.$$
From these minimal assumptions we may still deduce a great deal about $\Phi:$
By setting $\mu$ to one of the basepoints $x_i$ and shrinking $\sigma,$ we can focus almost all the probability within $\mathcal{R}_i.$ That is, $$\lim_{\sigma\to 0^+} \Phi(x_1,\sigma) = (1,0,0)$$ with a similar statement for $x_2$ and $x_3.$ Thus, the image of $\Phi$ is arbitrarily close to the vertices of the simplex. In fact, the image will include the vertices when the support of $F$ is bounded (simply choose $\sigma$ smaller than the distance from $x_i$ to the boundary of $\mathcal{R}_i$ divided by the diameter of the support of $F$).
By following a path $\gamma_{12}$ from $\mathcal{R}_1$ to $\mathcal{R}_2$ that avoids $\mathcal{R}_3$ and shrinking $\sigma$ sufficiently, we obtain a continuous map $$t \to \Phi(\gamma_{12}(t),\sigma)$$ whose endpoints are arbitrarily close to $(1,0,0)$ and $(0,1,0)$ (by virtue of the preceding observation) and for which the third component of the image is as small as we may like. Thus, the image of this path is arbitrarily close to the edge of $S_2$ from vertex $(1,0,0)$ to vertex $(0,1,0).$ Since the same construction holds for vertices $2,3$ and vertices $3,1,$ we conclude the image of $\Phi$ is arbitrarily close to the edges of the simplex. (Again, when the support of $F$ is bounded, the image of $\Phi$ includes the edges.)
It immediately follows from topological considerations that the image of $\Phi$ is the entire simplex. (A rigorous proof requires homology theory or homotopy theory, but intuitively it amounts to the idea that if the image of $\Phi$ omitted any point in the interior of $S_2,$ then $\Phi$ would have to be discontinuous--it would have to "tear" its domain--in order to create that hole.)
In particular, any triple of prescribed probabilities $(p_1,p_2,p_3),$ being a point of $S_2,$ must lie in the closure of the image of $\Phi$ (and if none of these probabilities is zero, it must actually be in the image of $\Phi$).
This has proven that
For any good partition $(\mathcal{R}_1,\mathcal{R}_2,\mathcal{R}_3)$ of the plane, any probability triple $(p_1,p_2,p_3),$ and any continuous 2D distribution $F,$ there exists a location $\mu\in\mathbb{R}^2$ and a scale factor $\sigma\gt 0$ for which $F(\mathcal{R}_i;\mu,\sigma) \approx p_i$ with to desirable degree of accuracy (and perfect accuracy when all the $p_i$ are nonzero or $F$ has bounded support).
#
# Find the probability of the angular sector at `origin` with nonzero oriented
# direction vectors x1 and x2 for a Binormal(mu, Sigma) distribution.
#
psector <- function(origin, x1, x2, mu=c(0,0), Sigma=diag(1,2)) {
require(mvtnorm)
# Are x1 and x2 collinear?
z <- zapsmall(c(sum(x1 * (rev(x2)*c(-1,1))), sqrt(sum(x1^2)*sum(x2^2))), digits=8)[1]
if (z == 0 && sum(x1 * x2) > 0) {
#
# Degenerate sector: interpret as the whole thing.
#
q <- 1
} else {
#
# Shift the mean.
#
mu <- mu - origin
#
# Rotate x1 to (1,0).
#
x1 <- x1 / sqrt(sum(x1^2))
O <- matrix(c(x1[1], -x1[2], x1[2], x1[1]), 2)
Sigma <- O %*% Sigma %*% t(O)
mu <- O %*% mu
x2 <- O %*% x2
#
# Transform x2 to (0,+) if possible, using an area-preserving transformation
#
complement <- isTRUE(sign(x2[2]) == -1) # Sector angle exceeds pi
if (zapsmall(x2)[2] == 0) { # Sector angle equals pi
lower <- c(-Inf, 0) # Upper half plane
upper <- c(Inf, Inf)
} else {
A <- matrix(c(1, 0, -x2[1]/x2[2], sign(x2[2])), 2)
Sigma <- A %*% Sigma %*% t(A) # Transforms x2 to (0,+) and x1 to (1,0)
mu <- A %*% mu
lower <- c(0, 0) # First quadrant
upper <- c(Inf, Inf)
}
q <- pmvnorm(lower=lower, upper=upper, mean=c(mu), sigma=Sigma)
if (complement) q <- 1 - q
}
return(q)
}
#
# Given three sectors originating at `origin` bounded by directions `x1`, `x2`,
# and `x3` (oriented positively), along with target probabilities (p1, p2, p3)
# (summing to unity), find `mu` for which a standard Normal distribution
# centered at `mu` has the given probabilities in the three sectors.
#
f <- function(p, origin, x1, x2, x3, ...) {
p <- p / sum(p)
mu <- function(x) origin + exp(x[1]) * c(cos(x[2]), sin(x[2]))
obj <- function(x) {
m <- mu(x)
(p[1] - psector(origin, x1, x2, m))^2 +
(p[2] - psector(origin, x2, x3, m))^2 +
(p[3] - psector(origin, x3, x1, m))^2
}
x.hat <- nlm(obj, c(0,0), fscale=1e-14, ...)$estimate
mu(x.hat)
}
#
# Example
#
x1 <- c(1,1)
x2 <- c(-1,0)
x3 <- c(0,-1)
p <- c(0.2, 0.1, 0.7)
origin <- c(0,0)
mu <- f(p, origin, x1, x2, x3)
