The coexistence of water, vapor, and solid phases and the spatiotemporal changes in thermal conductivities in response to variations in water content pose a challenge to the heat‐tracing method used to calculate water flux in variably saturated zones. In this study, a new approach, which adopted finite difference numerical schemes along with the limited‐memory Broyden–Fletcher–Goldfarb–Shanno bounded (L‐BFGS‐B) parameter optimization algorithm, was developed to inversely trace vertical soil water fluxes in variably saturated zones using temperature‐depth time‐series data. The newly developed method did not utilize hydraulic parameters, and only temperature time‐series data were used at the adjacent depths as boundary conditions. The soil water, vapor flux, and spatiotemporal variations in thermal conductivity at a given location were optimized by fitting the temperature time‐series at the target depth using the L‐BFGS‐B algorithm. Further, the method was verified through laboratory experiments under different saturated conditions that were conducted in two columns with saturated and free drainage bottom boundaries, and five sets of different flow rates with a relatively stable flow condition. The soil water flux calculated using the proposed method was consistent with the observed data, thus, verifying the reliability of the approach in applications with relatively stable flow rates. Furthermore, the new method was applied to a field site to estimate the time‐varying soil water flux using the measured temperature data. The trends of the calculated soil water flux at different depths were consistent with the field hydrologic conditions and observations. Thus, this study provides a novel computationally efficient and easy‐to‐implement tool to determine soil water flux using temperature measurements.