In the gradient method, soil heat flux density at a known depth G is determined as the product of soil thermal conductivity λ and temperature T gradient. While measuring λ in situ is difficult, many field studies readily support continuous, long-term monitoring of soil T and water content θ in the vadose zone. In this study, the performance of the gradient method is evaluated for estimating near-surface G using modeled λ and measured T. Hourly λ was estimated using a model that related λ to θ, soil bulk density ρb, and texture at 2-, 6-, and 10-cm depths. Soil heat flux Gm was estimated from modeled λ and measured T gradient (from thermocouples). The Gm results were evaluated with heat flux data GHP determined using independent measured λ and T gradient from heat-pulse probes. The λ model performed well at the three depths with 3.3%–7.4% errors. The Gm estimates were similar to GHP (agreed to within 15.1%), with the poorest agreement at the 2-cm soil depth, which was caused mainly by the relatively greater variability in ρb. Accounting for temporal variations in ρb (with core method) improved the accuracies of λ and Gm at the 2-cm depth. Automated θ monitoring approaches (e.g., time domain reflectometry), rather than gravimetric sampling, captured the temporal dynamics of near-surface λ and G well. It is concluded that with continuous θ and T measurements, the λ model–based gradient method can provide reliable near-surface G. Under conditions of soil disturbance or deformation, including temporally variable ρb, data improves the accuracy of G data.