<p>Geothermal heat flow (GHF) comprises comprehensive data on geothermal temperature gradients, rock thermal conductivity, and crustal/mantle heat flow, which is crucial for evaluating regional geothermal resources and conducting other studies. However, the traditional hyperparameter tuning methods such as grid search and Bayesian optimization (BO) are computationally intensive and susceptible to local optima, thereby complicating the accurate prediction of GHF. Moreover, the risk of overfitting poses a significant challenge in developing robust predictive models of GHF. To address these challenges, this study employed a genetic algorithm (GA) to tune the hyperparameters of the gradient boosted regression tree (GBRT) model. This research integrated the GBRT model with GA to establish a GA–GBRT hybrid model for predicting GHF more accurately. The GA–GBRT hybrid model incorporates 17 geological and geophysical features from Henan Province and its surroundings for sample data training. It is demonstrated that the GA–GBRT model improves the generalization performance of the test dataset and increases the <i>R</i>² by 11.75% compared with the GBRT model after BO. Through prediction performance analysis, the GA–GBRT model outperformed the existing optimized random forest (RF), deep neural network (DNN), and traditional interpolation methods in all of the studied cases. Based on the GA–GBRT model, a new GHF map was created and it demonstrated a more rational representation of the GHF distribution within the area of interest compared to traditional interpolation outcomes. The superiority of the GA–GBRT model to predict GHF was validated by the geological background, geodynamics, geophysical information, and high-temperature hot springs data from Henan Province.</p>