AbstractSteady state 1‐D analytical solutions to estimate groundwater fluxes from temperature profiles are an attractive option because they are simple to apply, with no complex boundary or initial conditions. Steady state solutions have been applied to estimate both aquifer scale fluxes as well as to estimate groundwater discharge to streams. This study explores the sources of uncertainty in flux estimates from regional scale aquifers caused by sensor precision, aquifer heterogeneity, multidimensional flow and variations in surface temperature due to climate change. Synthetic temperature profiles were generated using 2‐D groundwater flow and heat transport models with homogeneous and heterogeneous hydraulic and thermal properties. Temperature profiles were analyzed assuming temperature can be determined with a precision between 0.1°C and 0.001°C. Analysis of synthetic temperature profiles show that the Bredehoeft and Papadopulos (1965) method can provide good estimates of the mean vertical Darcy flux over the length of the temperature profile. Reliable flux estimates were obtained when the ratio of vertical to horizontal flux was as low as 0.1, and in heterogeneous media, providing that temperature at the upper boundary was constant in time. However, temporal increases in surface temperature led to over‐estimation of fluxes. Overestimates increased with time since the onset of, and with the rate of surface warming. Overall, the Bredehoeft and Papadopulos (1965) method may be more robust for the conditions with constant temperature distributions than previously thought, but that transient methods that account for surface warming should be used to determine fluxes in shallow aquifers.