Soil organic carbon (SOC) plays a crucial role in global climate change, the carbon cycle, and agricultural productivity, making accurate predictions of SOC changes in a region highly significant. However, due to the complex process of SOC changes, there are many confounding variables and it is not easy to derive robust predictions. The key to the solution is to remove or control these confounding factors. In response to this challenge, this study proposed a method combining causal inference with machine learning to get robust predictions of SOC storage changes. The method first identifies direct and indirect causal variables affecting temporal changes in SOC storage using structural equation modeling (SEM). It then directly predicts the temporal changes with those causal variables based on a newly developed method called two-point machine learning (TPML), rather than comparing spatial interpolation results across different times. In this way, the confounding variables can be removed and it is abbreviated as SemTPML. The SemTPML method was used in a case study of surface SOC (0-10 cm) of Shaanxi Province. The results show that it produces more robust predictions and the highest accuracy. NDVI and average annual precipitation (APre) were identified as the main controlling factors of surface SOC changes in Shaanxi Province. The results also revealed that changes in surface SOC from 1980 to 2020 in Shaanxi Province exhibit a trend of "increasing in the south and decreasing in the north", with the total changes amounting to a reduction of approximately 1.59 x 107 kg.