Past changes in the Earth's surface energy balance propagate into the subsurface and appear as perturbations of the subsurface thermal regime. This paper presents a singular value decomposition inversion method used to reconstruct surface heat flux histories (SHFH) from the heat flux anomalies detected in the shallow subsurface. Synthetic tests were used to assess the robustness of the inversion procedure. It was found that data noise can have a significant effect on the stability of the SHFH inferred from inversion. This translates in SHFHs having lower temporal resolution than ground surface temperature histories (GSTHs) obtained from the same data, but the long-term trends are robust. Results are encouraging for temperature data noise levels typically encountered in field measurements. Synthetic data tests yield results in agreement with analytical expressions derived from GSTHs for the same parameterization. Temperature-depth profiles from Canada's geothermal database were used to illustrate the inversion procedure. Individual temperature profile inversions are shown as examples. All 112 temperature logs in the database were used to obtain a mean heat flux history for the region. Results indicate that the ground heat flux has increased an average of 24 mW m-2 over the last 200 years in Canada. Application of this method to the existing global geothermal data base should allow for a quantification of the global energy balance at the Earth's surface for the past few centuries and may be useful for land surface models. ¿ 2001 American Geophysical Union