Retrieving the subsurface and deeper ocean (SDO) dynamic parameters from satellite observations is crucial for effectively understanding ocean interior anomalies and dynamic processes, but it is challenging to accurately estimate the subsurface thermal structure over the global scale from sea surface parameters. This study proposes a new approach based on Random Forest (RF) machine learning to retrieve subsurface temperature anomaly (STA) in the global ocean from multisource satellite observations including sea surface height anomaly (SSHA), sea surface temperature anomaly (SSTA), sea surface salinity anomaly (SSSA), and sea surface wind anomaly (SSWA) via in situ Argo data for RF training and testing. RF machine‐learning approach can accurately retrieve the STA in the global ocean from satellite observations of sea surface parameters (SSHA, SSTA, SSSA, SSWA). The Argo STA data were used to validate the accuracy and reliability of the results from the RF model. The results indicated that SSHA, SSTA, SSSA, and SSWA together are useful parameters for detecting SDO thermal information and obtaining accurate STA estimations. The proposed method also outperformed support vector regression (SVR) in global STA estimation. It will be a useful technique for studying SDO thermal variability and its role in global climate system from global‐scale satellite observations.

Accurately retrieving and describing subsurface temperature on a large scale can provide valuable information that can be used for subsurface dynamic and variability studies. This study develops a new satellite‐based geographically weighted regression (GWR) model to estimate a subsurface temperature anomaly (STA) in the upper 2,000 m of the Indian Ocean by combining satellite observations (sea surface height, sea surface temperature, sea surface salinity, and sea surface wind) and Argo in situ data (STA). This model improves the estimation accuracy by considering the significant spatial nonstationarity feature between the surface and subsurface parameters in the ocean. The performance of the GWR model is measured by using Akaike Information Criterion combined with root‐mean‐square error and R2. The results showed that the proposed GWR model can easily retrieve the STA and outperform the ordinary least squares model. The GWR model can also explain the contribution from each variable via a local regression coefficient distribution. The sea surface height from altimetry is the most significant variable for GWR estimation. This study demonstrates the great potential and advantage of the GWR model for large‐scale subsurface modeling and information retrieving. Thus, we have developed a novel approach for investigating subsurface thermal anomaly and variability from satellite observations.