$\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 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(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$.