Description: Glacial isostatic adjustment is largely governed by the rheological properties of the Earth's mantle. Large mass redistributions in the ocean–cryosphere system and the subsequent response of the viscoelastic Earth have led to dramatic sea level changes in the past. This process is ongoing, and in order to understand and predict current and future sea level changes, the knowledge of mantle properties such as viscosity is essential. In this study, we present a method to obtain estimates of mantle viscosities by the assimilation of relative sea level rates of change into a viscoelastic model of the lithosphere and mantle. We set up a particle filter with probabilistic resampling. In an identical twin experiment, we show that mantle viscosities can be recovered in a glacial isostatic adjustment model of a simple three-layer Earth structure consisting of an elastic lithosphere and two mantle layers of different viscosity. We investigate the ensemble behaviour on different parameters in the following three set-ups: (1) global observations data set since last glacial maximum with different ensemble initialisations and observation uncertainties, (2) regional observations from Fennoscandia or Laurentide/Greenland only, and (3) limiting the observation period to 10 ka until the present. We show that the recovery is successful in all cases if the target parameter values are properly sampled by the initial ensemble probability distribution. This even includes cases in which the target viscosity values are located far in the tail of the initial ensemble probability distribution. Experiments show that the method is successful if enough near-field observations are available. This makes it work best for a period after substantial deglaciation until the present when the number of sea level indicators is relatively high.