Introduction
Nonasymptotic numerical simulation of transient heat transfer requires knowledge of initial conditions. If the heat transfer medium is soil, they typically use temperature logs (depthtemperature tables) of monitoring wells for the initial time point. Then one may either use scattereddata interpolation [1] or solve a steady heat transfer problem with Dirichlet boundary conditions at temperature measurement points. If you choose the second way and the problem is nonlinear, there might be some convergence issues.
Figure 1: Soil Temperature Interpolation in Frost 3D Universal
We use scattereddata interpolation [1] and our method provides a series of useful features, including the following:

In case of single monitoring well, isothermal surfaces follow the terrain, or more precisely, its smoothed (Gaussian blurred [2]) versions. The blur radius is depthproportional.

In the multiwell case, isothermal surfaces adapt to the terrain. In order to achieve this, we use a straightforward piecewiselinear triangulationbased interpolation [1] along isosurfaces of the smoothed depth field. The smoothed depths of soil points are initially computed as the differences between their smoothed (with Gaussian blur [2] of depthproportional intensity) elevations and altitudes. Then, if needed, the smoothed depth field is locally corrected to become strictly decreasing with respect to altitude.

Globally, our interpolation is rangerestricted. This guarantees that the resulting temperature won't go below the absolute zero. Also, if all measured temperatures are above (below) the freezing point, then the entire interpolation domain will be thawed (frozen).

Along the monitoring wells, the interpolation is shapepreserving, i.e. preserves monotonicity and convexity [3].
Statement of Problem
Let and let be the elevation map determining the part of the surface of the infinitely extended and infinitely deep soil .
We investigate the problem of interpolating the unknown temperature field
where are altitudes of the measurement points in the ^{th} monitoring well , , . Here is the "reasonable" continuation (see stage 3, below) of the function to the entire plane . The exact definitions of the infinitely extended and deep soil and its surface are the following:
Note that .
Algorithm
One may cope with the given interpolation problem in several stages. Let's list them.

Merge the duplicate wells and remove the empty (containing no measurement points) ones. If there are no wells now, stop. Renumber the remaining wells, redefine .

Merge the duplicate measurement points (for multiple temperature values, use their arithmetic mean as the new temperature). Renumber the updated measurement points in and redefine , .

Construct such a function that for all points , one has
where is the closest to point of . In particular, if , then and .

Construct the depth field ,
Obviously, .

For each nonempty 2dimensional soil layer
() , obtain a smoothed (with depthproportional intensity) version ,of the extended elevation field . Here
is Gaussian blur [2] with the standard deviation (radius) . We take
where the physically dimensionless constant doesn't depend on and is derived from empirical considerations.
Note that for all

Construct the smoothed depth field ,
For all
i.e. , just as for . Where needed, correct the field in such a way that the strict decreasing of with respect to would hold for all , but would remain true. From here we obtain that .

Compute the smoothed depths
of the measurement points. These depths are pairwise different within the monitoring well due to the strict monotonicity of with respect to , .
Construct such linear interpolation (for extrapolation, the nearest neighbor method is used) splines that

If , put
and stop.

Construct the convex hull and the Delaunay triangulation of the same set of points .

In each fixed point (let's denote its smoothed depth with
) , compute the temperature using the following algorithm. Let be the closest to point of the convex hull . In particular, if , then . If the point lies in some nondegenerate Delaunay triangle , and therefore there exist such uniquely determined scalars (barycentric coordinates) that and , then putThere's also the case when the Delaunay triangulation is empty (contains no nondegenerate triangles). In this situation (keep in mind that ), the point lies on some nondegenerate segment
() , i.e. , , where the choice of the scalars is unique. Put
It is clear from construction that isosurfaces of the resulting temperature field adapt to the terrain, but, as the depth increases, this adaptation weakens monotonically.
References

I. Amidror. Scattered data interpolation methods for electronic imaging systems: a survey. J. Electron. Imaging, 11(2), 157–176 (April 2002).

Wikipedia contributors. Gaussian blur. Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Gaussian_blur&oldid=672686407 (accessed August 3, 2015).

F. N. Fritsch and R. E. Carlson. Monotone Piecewise Cubic Interpolation. SIAM J. Numer. Anal., 17, 238–246 (1980).