The focus of the study is optimizing the technique for estimating geocenter motion and variations in J2 by combining data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission with output from an Ocean Bottom Pressure model and a Glacial Isostatic Adjustment (GIA) model. First, we conduct an end-to-end numerical simulation study. We generate input time-variable gravity field observations by perturbing a synthetic Earth model with realistically simulated errors. We show that it is important to avoid large errors at short wavelengths and signal leakage from land to ocean, as well as to account for self-attraction and loading effects. Second, the optimal implementation strategy is applied to real GRACE data. We show that the estimates of annual amplitude in geocenter motion are in line with estimates from other techniques, such as satellite laser ranging (SLR) and global GPS inversion. At the same time, annual amplitudes of C10 and C11 are increased by about 50% and 20%, respectively, compared to estimates based on Swenson et al. (2008). Estimates of J2 variations are by about 15% larger than SLR results in terms of annual amplitude. Linear trend estimates are dependent on the adopted GIA model but still comparable to some SLR results.