Surface heat flux is highly sensitive to small temperature fluctuations at varying timescales. The magnitude of these variations depend on the duration of climatic fluctuations and add significant uncertainty to present-day surface heat flow estimates. We assimilate heat flow data with multiple geophysical observations that resolve the crust in its present-day state to improve the constraints on subsurface thermal models and quantify the uncertainty in surface heat flux. The Curie depth isotherm is computed from magnetic data to constrain temperatures in the middle-lower crust, while P and S-wave velocities that we extract from tomographic models are sensitive to temperature to varying degrees throughout the lithosphere. We integrate these within an adjoint inversion framework that we apply to Southeastern Australia to invert the structure of thermal conductivity and heat sources within the crust. Based on previous inversions solely constrained by surface heat flow points and seismic velocity, we found that relatively high rates of heat production in Proterozoic crust control the variation of heat flux at the surface. This Proterozoic crust shares tectonic provenance with Antarctica and may have significant implications for its heat flow regime. Here, we will quantify the uncertainty reduction of thermal structure from assimilating Curie depth in addition to seismic velocity and heat flow observations. Our inversion framework can be easily adapted to integrate additional data types available for Antarctica to improve the precision of geothermal heat flux estimates.