One-dimensional analytical heat transport equations based on temperature time series data have become popular tools to quantify groundwater-surface water interactions. The influence of nonideal field conditions on the use of these equations has been assessed for nonsinusoidal stream temperature signals, uncertainty in thermal parameters, sensor accuracy and multidimensional flow. Given that streambeds are often highly heterogeneous, the influence of streambed heterogeneity on flux estimates from temperature time series requires further investigation. Synthetic streambed temperatures were generated using two-dimensional numerical models with heterogeneous hydraulic conductivity distributions. Streambed temperatures were used to calculate fluxes using methods based on amplitude ratios (Ar), phase shifts (Δφ) and both (ArΔφ). Calculated fluxes were compared to known fluxes from the numerical models for flow fields analogous to losing streams. The influence of streambed structure, degree of heterogeneity, depth of the sensor pair, and location along a flow path were assessed. Errors in calculated fluxes increased with sensor pair depth, position along a flow path, and with the degree of heterogeneity. These errors were larger for streambeds with isotropic structures compared with anisotropic structures, and of the three methods tested; the Δφ method produced the largest errors. The simultaneous estimation of strong fluxes using Δφ, and an inability to obtain a flux estimate from Ar can suggest the presence of low hydraulic conductivity zones. Given the large errors and inability to determine flow direction from the Δφ method, the Ar and ArΔφ methods are recommended for downwelling fluxes.