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 Y^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 representation of the process $$ \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_j \sigma^2_{t-l} + \sum_{m=1}^s \gamma_k Y^2_{t-m}\, \forall t \end{align} $$