a) $|f_{\epsilon}(z) - f(z)| = |\epsilon g(z)| < |f(z)|$ on $\partial D(0,1)$ for $\epsilon$ sufficiently small. So we apply Rouché's theorem and we conclude.
b) We use the generalized argument principle: with the usual hypotheses and $g$ holomorphic $$\frac{1}{2\pi i}\int_{\gamma} \frac{f'(z)}{f(z)}g(z) dz = \sum_{k}g(z_{k}) - g(p_{k})$$ where $z_{k}$ are the zeros of $f$ inside $\gamma$ and $p_{k}$ are the poles.
So in our case we have $$\frac{1}{2\pi i}\int_{\gamma} \frac{f'(z)+\epsilon g'(z)}{f(z) + \epsilon g(z)}z \ dz = z_{\epsilon}$$ where $\gamma = \partial D(0, 1)$
If a holomorphic function $h$ has a zero of multiplicity $k$ at $z_0$, then we have $h(z) = (z-z_0)^k\cdot \tilde{h}(z)$ with $\tilde{h}$ holomorphic in a neighbourhood of $z_0$ and $\tilde{h}(z_0) \neq 0$. Thus we have
$$\frac{h'(z)}{h(z)} = \frac{k(z-z_0)^{k-1}\tilde{h}(z) + (z-z_0)^k\tilde{h}'(z)}{(z-z_0)^k\tilde{h}(z)} = \frac{k}{z-z_0} + \frac{\tilde{h}'(z)}{\tilde{h}(z)},$$
and the second summand is holomorphic in a neighbourhood of $z_0$. The residue of $\frac{h'}{h}$ in $z_0$ is therefore $k$. Then
$$z\frac{h'(z)}{h(z)} = \frac{k\cdot z}{z-z_0} + z\frac{\tilde{h}'(z)}{\tilde{h}(z)},$$
so
$$\operatorname{Res} \biggl(z\frac{h'}{h}; z_0\biggr) = \operatorname{Res} \biggl(\frac{k\cdot z}{z-z_0}; z_0\biggr) = k\cdot z_0.$$
Here, $f_\epsilon$ has only one zero in the unit disk (counting multiplicities), and therefore
$$z_\epsilon = \frac{1}{2\pi i} \int_{\lvert z\rvert = 1} z\cdot\frac{f_\epsilon'(z)}{f_\epsilon(z)}\,dz = \frac{1}{2\pi i} \int_{\lvert z\rvert = 1} z\cdot\frac{f'(z) + \epsilon g'(z)}{f(z) + \epsilon g(z)}\,dz,$$
which shows that the function $\epsilon \mapsto z_\epsilon$ is continuous [even holomorphic], since the integrand depends continuously [holomorphically] on $\epsilon$.
Best Answer
A reference is IV.7.4 in J.B. Conway's Functions of one complex variable.
$_{\text{This was copied from a comment in an attempt to get the question off of the Unanswered list.}}$