In this paper we investigate the Bayesian approach to inverse Robin problems. These are problems for certain elliptic boundary value problems of determining a Robin coefficient on a hidden part of the boundary from Cauchy data on the observable part. Such a nonlinear inverse problem arises naturally in the initialisation of large-scale ice sheet models that are crucial in climate and sea-level predictions. We motivate the Bayesian approach for a prototypical Robin inverse problem by showing that the posterior mean converges in probability to the data-generating ground truth as the number of observations increases. Related to the stability theory for inverse Robin problems, we establish a logarithmic convergence rate for Sobolev-regular Robin coefficients, whereas for analytic coefficients we can attain an algebraic rate. The use of rescaled analytic Gaussian priors in posterior consistency for nonlinear inverse problems is new and may be of separate interest in other inverse problems. Our numerical results illustrate the convergence property in two observation settings.