Reconstructing the vertical structures of the ocean from sea surface information is of great importance for ocean and climate studies. In this study, an ensemble machine learning (Ens-ML) model is proposed to retrieve ocean subsurface thermal structure (OSTS) by using satellite-derived sea surface data and Argo data in the South China Sea (SCS). The input data include sea surface height (SSH), sea surface temperature (SST), sea surface salinity (SSS), sea surface wind (SSW), and geographic information (including longitude and latitude). We select three stable machine learning models, namely, extreme gradient boosting (XGBoost), RandomForest and light gradient boosting machine (LightGBM) as our benchmark models, and then use an artificial neural network (ANN) technique to combine outputs from the three individual models. The proposed Ens-ML model using sea surface data only by SSH, SST, SSS, and SSW performs less satisfactorily than that considering the contribution of geographical information, indicating that the geographical information is essential to estimate the OSTS accurately. The estimated OSTS from the Ens-ML model are compared with Argo data. The results show that the proposed Ens-ML model can accurately estimate the OSTS (upper 1000 m) in the SCS, which is relatively more accurate and precise than the individual models. The performance of the Ens-ML model also varies with season, and better estimation is obtained in winter, which is probably due to stronger mixing and weaker stratification. This study shows the great potential and advantage of the multi-model ensemble of machine learning algorithm for the ocean’s interior information retrieving, showing great potential in expanding the scope of ocean observations.