$\newcommand\vpi\varphi\newcommand\de\delta\newcommand\ep\varepsilon\newcommand\R{\Bbb R}$Continuous compactly supported functions are dense in $L^1$. So, the problem reduces to the following: 
>Given any continuous p.d.f. $f$ supported on a compact subset $K$ of $\R^d$, approximate $f$ in $L^1$ by a Gaussian mixture. 

Take any real $\ep>0$. Then for all small enough real $\de>0$ we have 
$$\|f_\de-f\|_1<\ep,$$
where $f_\de:=f*\vpi_\de$, $\vpi_\de$ is the p.d.f. of the centered normal distribution with covariance matrix $\de I_d$, and $I_d$ is the $d\times d$ identity matrix. 

Since $K$ is compact and $f$ is continuous, there is a finite partition $(K_j)_{j\in J}$ of $K$ into nonempty measurable sets $K_j$ such that $\|y-z\|<\de^2$ and $|f(y)-f(z)|<\ep$ for any $j\in J$ and all $y$ and $z$ in $K_j$, where $\|\cdot\|$ is the Euclidean norm. For each $j\in J$, take arbitrarily some $y_j\in K_j$. For $x\in\R^d$, let 
$$g(x):=\sum_{j\in J}f(y_j)|K_j|\vpi_\de(x-y_j), \tag{10}\label{10}$$
where $|\cdot|$ is the Lebesgue measure. Then 
\begin{align}
\|f_\de-g\|_1
\le\sum_{j\in J}\int dx\int_{K_j}dy\,|\vpi_\de(x-y)f(y)-\vpi_\de(x-y_j)f(y_j)|
\le S_1+S_2,
\end{align}
where
\begin{align}
S_1&:=\sum_{j\in J}\int dx\int_{K_j}dy\,|\vpi_\de(x-y)f(y)-\vpi_\de(x-y)f(y_j)| \\
&\le\ep\sum_{j\in J}\int_{K_j}dy\,\int dx\,\vpi_\de(x-y)
=\ep|K|,
\end{align}
\begin{align}
S_2&:=\sum_{j\in J}\int dx\int_{K_j}dy\,|\vpi_\de(x-y)f(y_j)-\vpi_\de(x-y_j)f(y_j)| \\
&=\sum_{j\in J}f(y_j)\int_{K_j}dy\,\int dx\,|\vpi_\de(x-y)-\vpi_\de(x-y_j)| \\ 
&=\sum_{j\in J}f(y_j)\int_{K_j}dy\;O(\de^2/\de)=O(M|K|\de),
\end{align}
where $M$ is a finite upper bound on $f$. 

Thus, letting $\ep$ and $\de$ be small enough, we can make $\|f_\de-g\|_1$ however small. Note that the sum $s:=\sum_{j\in J}f(y_j)|K_j|$ of the "weights" in \eqref{10} is $\int g$, which will be close to $\int f_\de=1$ (if $\|f_\de-g\|_1$ is small). So, by a small "vertical" re-scaling adjustment, we approximate, however closely, $f$ in $L^1$ by the Gaussian mixture 
$h:=g/s$.