Non-asymptotic numerical simulation of transient heat transfer requires knowledge of initial conditions. If the heat transfer medium is soil, they typically use temperature logs (depth-temperature tables) of monitoring wells for the initial time point. Then one may either use scattered-data interpolation  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 scattered-data interpolation  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 ) versions. The blur radius is depth-proportional.
In the multi-well case, isothermal surfaces adapt to the terrain. In order to achieve this, we use a straightforward piecewise-linear triangulation-based interpolation  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  of depth-proportional 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 range-restricted. 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 shape-preserving, i.e. preserves monotonicity and convexity .
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 .
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 ,
For each nonempty 2-dimensional soil layer
(), obtain a smoothed (with depth-proportional intensity) version ,
of the extended elevation field . Here
is Gaussian blur  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 ,
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
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 non-degenerate Delaunay triangle , and therefore there exist such uniquely determined scalars (barycentric coordinates) that and , then put
There's also the case when the Delaunay triangulation is empty (contains no non-degenerate triangles). In this situation (keep in mind that ), the point lies on some non-degenerate 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.