One main geodetic objective of the European Space Agency’s satellite mission GOCE (gravity field and steady-state ocean circulation explorer) is the contribution to global height system unification. This can be achieved by applying the Geodetic Boundary Value Problem (GBVP) approach. Thereby one estimates the unknown datum offsets between different height networks (datum zones) by comparing the physical (e.g. orthometric) height values H of benchmarks in different datum zones to the corresponding values derived from the difference between ellipsoidal heights h (e.g. determined by means of global navigation satellite systems) and geoid heights N . In the ideal case, i.e. neglecting data errors, the misfit between H and ( h − N ) is constant inside one datum zone and represents the datum offset. In practise, the accuracy of the offset estimation depends on the accuracy of the three quantities H, h and N , where the latter can be computed from the combination of a GOCE-derived Global Potential Model (GPM) for the long to medium wavelength and terrestrial data for the short wavelength content. Providing an optimum estimation of the datum offsets along with realistic error estimates, theoretically, requires propagation of the full error variance and covariance information of the GOCE spherical harmonic coefficients to geoid heights, respectively geoid height differences. From a numerical point of view, this is a very demanding task which cannot simply be run on a single PC. Therefore it is worthwhile to investigate on different levels of approximation of the full variance-covariance matrix (VCM) with the aim of minimizing the numerical effort. In this paper, we compare the estimation error based on three levels of approximation, namely (1) using the full VCM, (2) using only elements of the dominant m-block structure of the VCM and (3) using only the main diagonal of the VCM, i.e. neglecting all error covariances between the spherical harmonic coefficients. We show that the m-block approximation gives almost the same result as provided by the full VCM. The diagonal approximation however over- or underestimates the geoid height error, depending on the geographic location and therefore is not regarded to be a suitable approximation.