Water botton temperature variations (BTV) introduce many uncertainties into heat flow density (HFD) values determined from the temperature-depth (T-z) plot from one-time probe insertion into the top few meters of sediments. The BTV have to be carefully monitored for a number of years up to the date of the probe measurements; this makes the method time consuming and economically unattractive. In this paper we present the theory for a promising new technique. The method does not require monitoring of the BTV but does require at least two T-z profiles separated in time by only a few months. A generalized inversion of the BTV history using the separated T-z data, and any other a priori information that is available, is then performed in the frequency domain to give the most probable Bayesian estimates of various parameters and their error bounds. Numerical examples show that the inversion can resolve not only on the HFD but also the mean sediment surface temperature and gives information on the longer-term trend as well as details of the BTV over the previous year; these, together with thermal properties structure, are the factors that contribute most to the shape of the time-dependent T-z profiles used. In a typical situation where a uniform half-space with a background HFD=60 mW m-2 has been perturbed by complex BTV with a maximum variation of ¿1.8 K, using a misguessed HFD of 90¿100 mW m-2 in the inversion of two T-z profiles taken 3 months apart over a 2- to 5-m section results in a posteriori HFD of 58¿11 mW m-2. Although the technique gives an optimally determied mean HFD, it principal value lies in the ability to place reliable error bounds on the mean. ¿ American Geophysical Union 1987