$a_0$ will scale quadratically with the scaling of $y_t$, while $a_1$ and $b_1$ will stay the same. Here is why. Take the original model
\begin{aligned}
y_t &= \sqrt h_t \epsilon_t, \\
h_t &= a_0 + a_1 y^2_{t-1} + b_1 h_{t-1}, \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
and scale $y_t$ with a positive constant $c$. This turns the model into
\begin{aligned}
cy_t &= \sqrt{c^2 h_t} \epsilon_t, \\
c^2 h_t &= c^2 a_0 + a_1 (cy_{t-1})^2 + b_1 (c^2 h_{t-1}), \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
which can be expressed as
\begin{aligned}
\tilde y_t &= \sqrt{\tilde h_t} \epsilon_t, \\
\tilde h_t &= \tilde a_0 + a_1\tilde y^2_{t-1} + b_1 (\tilde h_{t-1}), \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
where $\tilde y_t = c y_t$ and $\tilde a_0 = c^2 a_0$.
However, numeric stability issues may create an impression that this is not the case. Try the following code (illustrated with two different R packages for GARCH models) and note how the estimated parameters do not scale all that nicely. (This might be due to the fact that the randomly generated data is i.i.d., not GARCH, making estimated the GARCH coefficients unstable.)
n=1000
set.seed(1); x=rnorm(n) # generate i.i.d N(0,1) data
# With ARMA(1,1) conditional mean
library(rugarch)
spec=ugarchspec()
fit2 =ugarchfit(data=x*100,spec=spec); coef(fit2 ) # scaled by 100
fit1 =ugarchfit(data=x*10 ,spec=spec); coef(fit1 ) # scaled by 10
fit0 =ugarchfit(data=x ,spec=spec); coef(fit0 ) # raw data
fit_1=ugarchfit(data=x/10 ,spec=spec); coef(fit_1) # scaled by 0.1
fit_2=ugarchfit(data=x/100,spec=spec); coef(fit_2) # scaled by 0.01
# With constant conditional mean
library(garchx)
fit2 =garchx(y=x*100); coef(fit2 )
fit1 =garchx(y=x*10 ); coef(fit1 )
fit0 =garchx(y=x ); coef(fit0 )
fit_1=garchx(y=x/10 ); coef(fit_1)
fit_2=garchx(y=x/100); coef(fit_2)
```