Mean dynamic topography (MDT) is crucial for research in oceanography and climatology. The optimal interpolation method (OIM) is applied to MDT modeling, where the error variance–covariance information of the observations is established. The global geopotential model (GGM) derived from GOCE (Gravity Field and Steady-State Ocean Circulation Explorer) gravity data and the mean sea surface model derived from satellite altimetry data are combined to construct MDT. Numerical experiments in the Kuroshio over Japan show that the use of recently released GOCE-derived GGM derives a better MDT compared to the previous models. The MDT solution computed based on the sixth-generation model illustrates a lower level of root mean square error (77.0 mm) compared with the ocean reanalysis data, which is 2.4 mm (5.4 mm) smaller than that derived from the fifth-generation (fourth-generation) model. This illustrates that the accumulation of GOCE data and updated data preprocessing methods can be beneficial for MDT recovery. Moreover, the results show that the OIM outperforms the Gaussian filtering approach, where the geostrophic velocity derived from the OIM method has a smaller misfit against the buoy data, by a magnitude of 10 mm/s (17 mm/s) when the zonal (meridional) component is validated. This is mainly due to the error information of input data being used in the optimal interpolation method, which may obtain more reasonable weights of observations than the Gaussian filtering method.