Solution (source code)

= Solution

Assume the original ice sheet had time to reach local isostatic compensation, melted rapidly relative to subsequent rebound, and left one exponentially relaxing depression mode. Neglect elastic layering and assign the whole present <gravity> deficit to that unrelaxed depression. Load balance gives
$$
h_0\simeq\frac{\rho_i}{\rho_m}H_i=900\,\mathrm m.
$$
The supplied sheet-gravity approximation gives the residual amplitude
$$
h(t)=\frac{|\Delta g|}{2\pi G\rho_m}
=\frac{2.0\times10^{-4}}{2\pi(6.67\times10^{-11})(3000)}
\simeq159\,\mathrm m.
$$
From $h(t)=h_0e^{-t/\tau}$,
$$
\boxed{\tau=\frac{10^4\,\mathrm{yr}}{\ln(900/159)}
\simeq5.77\times10^3\,\mathrm{yr}.}
$$
For an order-of-magnitude single-mode estimate, choose the ice-sheet diameter as <wavelength>, $\lambda\simeq2R=6000\,\mathrm{km}$. Taking $g=9.81\,\mathrm{m\,s^{-2}}$, the more accurate specified relation then gives
$$
\boxed{\eta=\frac{\rho_mg\lambda\tau}{4\pi}
\simeq2.6\times10^{21}\,\mathrm{Pa\,s}.}
$$
A circular load has a spectrum, not one uniquely specified <wavelength>. Choosing $\lambda\sim R$ instead gives about $1.3\times10^{21}\,\mathrm{Pa\,s}$; another <wavelength> convention rescales the answer linearly. \b[The robust inference is <viscosity> of order $10^{21}\,\mathrm{Pa\,s}$ within this simple model.] A realistic estimate would also model gradual melting, other <gravity> anomalies, <density> interfaces, <lithospheric flexure> and depth-dependent <mantle viscosity>.