The law of total (co)variance writes as
$$
\text{Cov}(\mathbf{x}) = \mathbb{E}[\text{Cov}(\mathbf{x} \, \vert \, y)] +
\text{Cov}[\mathbb{E}(\mathbf{x} \,\vert\, y)].
$$
With $q := 1 -p$, the first part (a.k.a. within group) is
$$
\mathbb{E}[\text{Cov}(\mathbf{x} \, \vert \, y)]
= q \boldsymbol{\Sigma}_0 + p \boldsymbol{\Sigma}_1.
$$
Using $\boldsymbol{\mu} := q \boldsymbol{\mu}_0 + p\boldsymbol{\mu}_1$,
the second part (a.k.a. between groups) is
$$
\text{Cov}[\mathbb{E}(\mathbf{x} \,\vert\, y)] =
q \left[\boldsymbol{\mu}_0 - \boldsymbol{\mu} \right]
\left[\boldsymbol{\mu}_0 - \boldsymbol{\mu} \right]^\top +
p \left[\boldsymbol{\mu}_1 - \boldsymbol{\mu} \right]
\left[\boldsymbol{\mu}_1 - \boldsymbol{\mu} \right]^\top =
qp \left[\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0 \right]
\left[\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0 \right]^\top
$$
where the second equality comes by simple algebra. This shows that
$\text{Cov}[\mathbb{E}(\mathbf{x} \,\vert\, y)]$ is a rank-one matrix.
Problem Given the diagonal conditional covariances
$\boldsymbol{\Sigma}_0$ and $\boldsymbol{\Sigma}_1$, given a
covariance matrix $\boldsymbol{\Sigma}$, we want to find
$\boldsymbol{\mu}_0$ and $\boldsymbol{\mu}_1$ such that
$\text{Cov}(\mathbf{x}) = \boldsymbol{\Sigma}$.
The matrix $\boldsymbol{\Sigma}$ should be such that the difference
$\boldsymbol{\Delta} := \boldsymbol{\Sigma} - q \boldsymbol{\Sigma}_0 -
p \boldsymbol{\Sigma}_1$ is positive (hence has positive diagonal
elements). Moreover, the difference should be a rank-one matrix.
If these two conditions are fulfilled, then $\boldsymbol{\Delta} = \boldsymbol{\delta}\boldsymbol{\delta}^\top$
for some vector $\boldsymbol{\delta}$. By choosing any $\boldsymbol{\mu}$
and then
$$
\boldsymbol{\mu}_0 := \boldsymbol{\mu} + \alpha_0 \, \boldsymbol{\delta},
\qquad
\boldsymbol{\mu}_1 := \boldsymbol{\mu} + \alpha_1 \, \boldsymbol{\delta}
$$
one gets the wanted covariance $\boldsymbol{\Sigma}$ for a suitable choice
of $\alpha_0$ and $\alpha_1$. We have then $\boldsymbol{\mu}_1 -\boldsymbol{\mu}_0 = [\alpha_1 - \alpha_0]\boldsymbol{\delta}$ so the constraints
are $qp [\alpha_1 - \alpha_0]^2 = 1$ and $q\alpha_0 + p \alpha_1 =0$
so we can take $\alpha_0 := \sqrt{p/q}$ and $\alpha_1:= - \sqrt{q/p}$.
The result does not restrict to the bivariate case as in OP. However with a
mixture of two $d$-dimensional Gaussian distributions, the rank-one
condition becomes more difficult to fulfil. We can then consider a
mixture of $m >2$ Gaussian distributions, the constraint then being
that the between-groups covariance $\boldsymbol{\Delta}$ has rank $\leqslant m-1$. Note also that taking the covariance matrices of the component
as zero, we consider a mixture of Dirac distributions. In order to get
an arbitrary mean $\boldsymbol{\mu}$ and an arbitrary covariance
$\boldsymbol{\Sigma}$ for a mixture of $d$ dimensional Gaussians with
given weights we need to use $d + 1$ distributions.
EDIT For a generalisation, consider the case where $\mathbf{x}$ has length
$d$ and a mixture of $d + 1$ Gaussian distributions with covariance matrices
$\boldsymbol{\Sigma}_i$ and with a vector $\mathbf{p}$ of $d+1$ weights $p_i >0$
with $\sum_{i=1}^{d+1}p_i = 1$. We claim that provided that the
matrix $\boldsymbol{\Delta}:= \boldsymbol{\Sigma} - \sum_{i=1}^{d+1}
p_i \boldsymbol{\Sigma}_i$ is positive, we can find $d+1$ mean vectors
$\boldsymbol{\mu}_i$ so that the mixture $\sum_i p_i \,\texttt{Norm}(\boldsymbol{\mu}_i,\,
\boldsymbol{\Sigma}_i)$
has the given mean $\boldsymbol{\mu}$
and the given covariance $\boldsymbol{\Sigma}$.
Since $\boldsymbol{\Delta}$ is positive we can write it as
$\boldsymbol{\Delta} = \mathbf{V}\mathbf{V}^\top$ where $\mathbf{V}$
is a $d \times d$ matrix. For that aim, the eigendecomposition of
$\boldsymbol{\Delta}$ can be used or a Cholesky decomposition. Let us
temporarily admit that we can find a $d \times (d +1)$ matrix
$\mathbf{A}$ such that
$$
\tag{1}
\left\{
\begin{array}{c c}
\mathbf{A} \text{diag}(\mathbf{p}) \mathbf{A}^\top &= \mathbf{I}_d,\\
\mathbf{A} \mathbf{p} &= \mathbf{0}_d,\rule{0pt}{1.2em}
\end{array}
\right.
$$
where $\mathbf{I}_d$ and $\mathbf{0}_d$ are the identity matrix and
the vector of zeros.
Then with $\boldsymbol{\alpha}_i$ being the $i$-th column of the $d \times (d + 1)$ matrix
$\mathbf{A}$, take
$$
\boldsymbol{\mu}_i := \boldsymbol{\mu} + \mathbf{V} \boldsymbol{\alpha}_i
\qquad i=1, \, \dots,\, d+1.
$$
Let us check that we get: the wanted covariance and the wanted mean. Firstly
$$
\sum_{i=1}^{d+1} p_i \left[\boldsymbol{\mu}_i - \boldsymbol{\mu} \right]
\left[\boldsymbol{\mu}_i - \boldsymbol{\mu} \right]^\top =
\sum_{i=1}^{d+1} p_i \mathbf{V} \boldsymbol{\alpha}_i
\boldsymbol{\alpha}_i^\top \mathbf{V}^\top
= \mathbf{V}
\left\{ \sum_{i=1}^{d+1} p_i
\boldsymbol{\alpha}_i
\boldsymbol{\alpha}_i^\top
\right\} \mathbf{V}^\top = \mathbf{V} \mathbf{V}^\top = \boldsymbol{\Delta},
$$
since the matrix between the curly brackets is
$\mathbf{A} \text{diag}(\mathbf{p}) \mathbf{A}^\top = \mathbf{I}_d$ from (1). Secondly
$$
\sum_{i=1}^{d+1} p_i \left[\boldsymbol{\mu}_i - \boldsymbol{\mu} \right] =
\sum_{i=1}^{d+1} p_i \mathbf{V} \boldsymbol{\alpha}_i = \mathbf{V}
\left\{ \sum_{i=1}^{d+1} p_i
\boldsymbol{\alpha}_i \right\} = \mathbf{V} \mathbf{A} \mathbf{p} = \mathbf{0}
$$
because of the second condition in (1). So $\sum_{i=1}^{d+1} p_i \boldsymbol{\mu}_i = \boldsymbol{\mu}$ as wanted.
Now let us show how the matrix $\mathbf{A}$ in (1) can be found. By separating the
first $d$ columns and the last one, let $\mathbf{A} =: [\mathbf{A}_1 \vert \mathbf{A}_2]$.
Let $\mathbf{D}_1 := \text{diag}(\mathbf{p}_{1:d})$ and $\mathbf{D}_2 := p_{d+1}$. We want that
$$
\left\{
\begin{array}{c c}
\mathbf{A}_1 \mathbf{D}_1 \mathbf{A}_1^\top +
\mathbf{A}_2 \mathbf{D}_2 \mathbf{A}_2^\top
&= \mathbf{I}_d,\\
\mathbf{A}_1 \mathbf{p}_{1:d} + \mathbf{A}_2 p_{d+1} &= \mathbf{0}_d.\rule{0pt}{1.2em}
\end{array}
\right.
$$
The second equation gives $\mathbf{A}_2 = -\mathbf{A}_1 \mathbf{p}_{1:d} /p_{d+1}$, and then the first one
$$
\mathbf{A}_1 \left\{ \mathbf{D}_1 + \frac{1}{p_{d+1}^2}
\mathbf{p}_{1:d}\mathbf{p}_{1:d}^\top \right\} \mathbf{A}_1^\top
= \mathbf{I}_d
$$
we can find $\mathbf{A}_1$ from the Cholesky decomposition of the
matrix between the curly brackets.
