<p>The paleo-geothermal gradient is a crucial parameter for converting the thermal history to the exhumation history. However, the precise estimation of this parameter has been a challenge. This paper presents a simple two-step method to model the paleo-geothermal gradient using low-temperature thermochronology. (1) It uses the Monte Carlo approach to generate thermal histories in a vertical section randomly and calculates the entire thermal history within the goodness-of-fit thresholds based on different paleo-geothermal gradients. (2) It selects the optimum paleo-geothermal gradient by comparing the entire thermal history within different goodness-of-fit thresholds. We validated the method with apatite (U-Th)/He and fission track data collected from two drill cores in the Haiyuan-Liupanshan region. The result revealed that the best-fit paleo-geothermal gradient was ∼42 °C/km during the Early Cretaceous–Miocene and has decreased rapidly to 20 °C/km since ∼10 Ma. The crust thickening in the study area may explain the rapid reduction in the paleo-geothermal gradient since ∼10 Ma. Our results are consistent with earlier studies in the region, suggesting that our simple and more intuitive approach provides an alternative method for paleo-geothermal gradient modeling.</p>