Geothermal heat flow (GHF) is a critical parameter for understanding the thermal structure and dynamics of the lithosphere, offering key insights into geophysical processes and geothermal energy potential. This study investigates the spatial variability of GHF in Germany by applying a Bayesian Markov Chain Monte Carlo method to estimate key thermal parameters, including crustal and mantle thermal conductivities, crustal heat production, and mantle heat flow. The analysis integrates data on surface heat flow, surface temperatures, and the lithosphere-asthenosphere boundary (LAB) depth.   To address the limitations posed by the sparse and uneven distribution of direct borehole measurements, comprising only 595 GHF records, we incorporated a wide range of geophysical and geological datasets, such as gravity, magnetics, seismic velocity, topography, and proximity to faults and volcanic regions. These datasets were analyzed using a Quantile Regression Forest (QRF) approach, which enabled robust GHF estimations while accounting for uncertainties and providing reliable prediction intervals. This methodology significantly improves upon traditional Curie depth-based methods, providing a more accurate and comprehensive GHF model for Germany.   The probabilistic multi-geoobservable modelling presented here enhances our understanding of the geothermal regime in Germany, contributing to a more precise assessment of its geothermal resources and the thermal state of its lithosphere.