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 G_m was estimated from modeled λ and measured T gradient (from thermocouples). The G_m results were evaluated with heat flux data G_HP 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 G_m estimates were similar to G_HP (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 G_m 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.