Essentially, it is a free parameter that you can tune arbitrarily.
In the theoretical limit $\eta \to 0^+$, the spectral function associated to the Green function, i.e. $A(\omega) = - \frac{1}{\pi} \mathrm{Im}G(\omega + i\eta)$ is a Dirac delta peaked at $\varepsilon(k)$. With a finite $\eta$, the peaks are actually broadened into Lorentzian functions, with $\eta$ being proportional to the width of the function. Roughly speaking, the spectral function is non-vanishing in an interval of values $\varepsilon(k) - \eta < \omega < \varepsilon(k) + \eta$.
So numerically you can choose an arbitrary $\eta$, and the effect is simply that you will observe a broader spectral line, which is a numerical artifact that approximates a Dirac delta peak.
To prove this you can just compute $A(\omega)$ for a finite $\eta$:
$$
A(\omega) = - \frac{1}{\pi} \mathrm{Im} \frac{1}{\omega - \varepsilon(k) + i\eta} = - \frac{1}{\pi} \mathrm{Im} \frac{\omega-\varepsilon(k)-i\eta}{[\omega - \varepsilon(k)]^2 + \eta^2} =\frac{1}{\pi} \frac{\eta}{[\omega - \varepsilon(k)]^2 + \eta^2},
$$
which is in fact a Lorentzian curve with $2\eta$ being the full width at half maximum.
If instead we consider the limit $\eta \to 0^+$, we need the Sokhotski–Plemelj theorem: $\lim_{\eta \to 0^+}\mathrm{Im} \frac{1}{\omega-\varepsilon+i\eta} = -\pi \delta(\omega-\varepsilon)$, which leads to
$$
A(\omega) = \delta(\omega - \varepsilon(k)).
$$