Under the assumption that the Earth's thermal field is one-dimensional and purely conductive, the temperature w is related to the Earth model m through a partial differential equation (PDE), where m is the set of model parameters consisting of the ground surface temperature, the background heat flow density, the thermal conductivity, the specific heat capacity, and the rate of heat production; by setting the time origin sufficiently far in the past, the initial temperature field may be taken as the steady state temperature field. Given data (d0, Cd) on w and a priori information (m0, Cm) on m, where Cd and Cm are covariance operators describing uncertainties in d0 and m0, respectively, the aim of the least squares inversion is to determine the most probable model that minimizes the misfit function S=1/2⟨C-1d(d-d0), d-d0⟩ +1/2⟨C-1m(m-m0), m-m0⟩. We formulate this problem in the functional space as opposed to the conventional discrete formulation and solve it using iterative gradient methods. The formulation reduces the computation in each iteration to essentially two forward solutions of the PDE, the first for the primal problem: given m, solve for the actual field w, and the second for the dual problem: using the weighted data residuals as heat source, solve for the residual temperature field in the same medium, but with homogeneous boundary conditions and with time reversed. The correlation of the residual and the actual fields, gives the gradient and also the Hessian of S, the latter of which evaluated at the most probable model is the approximate a posteriori covariance operator. Because discretization is required only when solving the forward problems, we avoid the computing and storing of partial derivatives of d with respect to discretized m, which can be a prohibitive task when the number of data and the number of discretized m are large. ¿American Geophysical Union 1991