The rocky interior of the Earth flows viscoelastically over timescales on the order of 1000 years in response to sustained stresses. Such flow is still occurring today as a result of the growth and collapse over the last ice age of massive ice sheets and is evident in changes of the Earth's surface and gravity, a process called glacial isostatic adjustment (GIA). This thesis presents a new technique for computing this viscoelastic deformation and statistical methods for more efficiently inferring properties of the Earth's mantle and the deglaciation from geophysical observations. The first chapter introduces an updated time-domain method for computing the viscoelastic Love numbers --- normalized spherical harmonic responses of an Earth with radially symmetric properties. The method employs a novel normalization and coordinate transformation that, when used in combination with the relaxation method for two-point boundary value problems, yields a very effective method of computation that is applicable to a wide range of possible rheological models. The second chapter describes a geometric perspective of GIA modeling using a heuristic example of the sea level response of a single ice cap melting, a prototype of a full inversion of global rheology and deglaciation. By considering the locus of all possible model predictions, a surface called the model manifold, we demonstrate universal features of nonlinear models, such as edges where parameters unphysically go to infinity, and how these can interfere when inferring parameters from data. Applying geometric corrections to the Levenberg-Marquardt least-squares algorithm facilitate finding the best-fit on the model manifold without getting stuck on an edge, even when started from far away. The final chapter employs a different aspect of this perspective, optimal experiment design, to evaluate the geophysical constraints on the configuration and volume of the Barents Sea Ice Sheet over the last glacial cycle and propose maximally constraining observations. Available observations of GIA in the Barents Sea cannot distinguish between a single, large dome and a more moderate amount of ice in the north. Experimental design identifies an area in the central Barents Sea within which a single observation of uplift would be very constraining.