The mixed layer depth (MLD) is generally estimated using in situ or model data. However, MLD analyses have limitations due to the sparse resolution of the observed data. Therefore, this study reconstructs three-dimensional (3D) ocean thermal structures using only satellite sea surface measurements for a higher spatial and longer temporal resolution than that of Argo and diagnoses the decadal variation of global MLD variability. To simulate the ocean thermal structures, the relationship between the ocean subsurface temperature and the sea surface fields was computed based on gridded Argo data. Based on this relationship, high spatial resolution and extended periods of satellite-derived altimeter, sea surface temperature (SST), and wind stress data were used to estimate the 3D ocean thermal structures with 0.25° spatial resolution and 26 standard depth levels (5–2000 m) for 24 years (1993–2016). Then, the MLD was calculated using a temperature threshold method (∆T = 0.2 °C) and correlated reasonably well (>0.9) with other MLD datasets. The extended 24-year data enabled us to analyze the decadal variability of the MLD. The global linear trend of the 24-year MLD is −0.110 m yr−1; however, from 1998 to 2012, the linear trend is −0.003 m yr−1 which is an order of magnitude smaller than that of other periods and corresponds to a global warming hiatus period. Via comparisons between the trends of the SST anomalies and the MLD anomalies, we tracked how the MLD trend changes in response to the global warming hiatus.