We predict present-day rates of change of the long-wavelength components of the Earth's geopotential due to the late Pleistocene glacial cycles (hereinafter referred to as glacial isostatic adjustment, or GIA). These predictions are generated using spherically symmetric, self-gravitating, (Maxwell) viscoelastic Earth models. Previous studies have generally considered only zonal (i.e., azimuthally independent) harmonics, or the so-called coefficients. We extend these efforts to focus on the nonzonal harmonics and explore the sensitivity of our predictions to changes in a variety of model inputs. As an example, we examine the influence on the GIA predictions of ignoring rotational perturbations and assuming an elastically incompressible model and find that these assumptions can have a significant impact on the degree-2, order-1 coefficients, while the other harmonics are generally affected by <10%. Though the GIA predictions are generally insensitive to variations in lithospheric thickness and upper mantle viscosity, reasonable variations in either lower mantle viscosity or the late Pleistocene ice history will have a more significant effect. Mindful of upcoming dedicated gravity missions, we show that predictions of the GIA-induced changes of the geoid can vary by at least 0.5 mm yr⁻¹ over a plausible range of Earth and ice models.