We consider the problem of inferring the basal sliding coefficient field for an uncertain Stokes ice sheet forward model from surface velocity measurements. The uncertainty in the forward model stems from uncertain auxiliary parameters (e.g., rheology parameters). This inverse problem is posed within the Bayesian framework. To account for the associated model uncertainty, we employ the Bayesian Approximation Error (BAE) approach to approximately premarginalize simultaneously over both the noise in measurements and uncertainty in the forward model. We also carry out approximative posterior uncertainty quantification based on a linearization of the parameter-to-observable map centered at the maximum a posteriori (MAP) basal sliding coefficient estimate. The MAP estimate is found by minimizing the negative log posterior using an inexact Newton conjugate gradient method. The gradient and Hessian actions to vectors are efficiently computed using adjoints. Sampling from the approximate covariance is made tractable by invoking a low-rank approximation of the data misfit component of the Hessian. We study the performance of the BAE approach in the context of three numerical examples in two and three dimensions. For each example the basal sliding coefficient field is the parameter of primary interest, which we seek to infer, and the rheology parameters (e.g., the flow rate factor, or the Glen’s flow law exponent coefficient field) represent so-called nuisance parameters. Our results indicate that accounting for model uncertainty stemming from the presence of nuisance parameters is crucial. Namely our findings suggest that using nominal values for these parameters, as is often done in practice, without taking into account the resulting modeling error, can lead to overconfident and heavily biased results.