You are conflating the features of a process with its representation. Consider the (return) process $(Y_t)_{t=0}^\infty$. 

* An ARMA(p,q) model specifies the _conditional mean_ of the process as

$$
\begin{align}
\mathbb{E}(Y_t \mid \mathcal{I}_t) &= \alpha_0 + \sum_{j=1}^p \alpha_j Y_{t-j}+ \sum_{k=1}^q \beta_k\epsilon_{t-k}\\
\end{align}
$$
Here, $\mathcal{I}_t$ is the information set at time $t$, which is the $\sigma$-algebra generated by the lagged values of the outcome process $(Y_t)$. 

* The GARCH(r,s) model specifies the _conditional variance_ of the process
$$
\begin{alignat}{2}
& \mathbb{V}(Y_t \mid \mathcal{I}_t) &{}={}& \mathbb{V}(\epsilon_t \mid \mathcal{I}_t) \\
\equiv \,& \sigma^2_t&{}={}& \delta_0 + \sum_{l=1}^r \delta_j \sigma^2_{t-l} + \sum_{m=1}^s \gamma_k \epsilon^2_{t-m}
\end{alignat}
$$

Note in particular the first equivalence $ \mathbb{V}(Y_t \mid \mathcal{I}_t)= \mathbb{V}(\epsilon_t \mid \mathcal{I}_t)$.

_Aside_: Based on this representation, you can write
$$
\epsilon_t \equiv \sigma_t Z_t
$$
where $Z_t$ is a strong white noise process, but this follows from the way the process is defined.

* The two models (for the conditional mean and the variance) are perfectly compatible with each other, in that the mean of the process can be modeled as ARMA, and the variances as GARCH. This leads to the complete specification of an ARMA(p,q)-GARCH(r,s) model for the  process as in the following representation
$$
\begin{align}
Y_t &= \alpha_0 + \sum_{j=1}^p \alpha_j Y_{t-j} + \sum_{k=1}^q \beta_k\epsilon_{t-k} +\epsilon_t\\
 \mathbb{E}(\epsilon_t\mid \mathcal{I}_t) &=0,\, \forall t \\
\mathbb{V}(\epsilon_t \mid \mathcal{I}_t) &= \delta_0 + \sum_{l=1}^r \delta_l \sigma^2_{t-l} + \sum_{m=1}^s \gamma_m \epsilon^2_{t-m}\, \forall t
\end{align}
$$