An analytical closed form solution to mean and variance in heat flow is obtained by solving the stochastic heat conduction equation incorporating randomness in the thermal conductivity. The method has a wide range of applications in quantifying the thermal state of the crust and in obtaining closed form expressions for subsurface heat flow structure along with its error bounds. The exact formulae employed can be used to better evaluate the thermal state for related oil and natural gas applications and also in tectonic studies and in studies related to the crystallization of minerals.