On Spatial Interpolation of Soil Temperature from Temperature Logs of Monitoring Wells


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 [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.
Soil Temperature Interpolation Based on Temperature Measurement Points

Figure 1: Soil Temperature Interpolation in Frost 3D Universal

We use scattered-data interpolation [1] and our method provides a series of useful features, including the following:

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

  2. 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 [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 depth-proportional intensity) elevations and altitudes. Then, if needed, the smoothed depth field is locally corrected to become strictly decreasing with respect to altitude.

  3. 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).

  4. Along the monitoring wells, the interpolation is shape-preserving, i.e. preserves monotonicity and convexity [3].

Statement of Problem

Let \Pi=[x_0,x_1]\times[y_0,y_1]\neq\emptyset and let h:\Pi\rightarrow\mathbb{R} be the elevation map determining the part \left\{\left.\left(A,h(A)\right)\right|A\in\Pi\right\} of the surface \partial S\subset\mathbb{R}^3 of the infinitely extended and infinitely deep soil S\subset\mathbb{R}^3.

We investigate the problem of interpolating the unknown temperature field T:S\rightarrow\mathrm{Conv}\left\{\left.T_{ij}\right|j=1,2,\ldots,m_i;i=1,2,\ldots,n\right\} (\mathrm{Conv}\,P stands for the convex hull of the set P of points) from the temperature logs

    \[\left[\begin{matrix}T\left(A_i,\mathrm{z}_{i1}\right) \\ T\left(A_i,\mathrm{z}_{i2}\right) \\ \vdots \\ T\left(A_i,\mathrm{z}_{im_i}\right) \end{matrix}\right]=\left[\begin{matrix}T_{i1} \\ T_{i2} \\ \vdots \\ T_{im_i} \end{matrix}\right]\!\!,~i=1,2,\ldots,n,\]

where \mathrm{z}_{i1},\mathrm{z}_{i2},\ldots,\mathrm{z}_{im_i} are altitudes of the measurement points \left(A_i,\mathrm{z}_{i1}\right),\left(A_i,\mathrm{z}_{i2}\right),\ldots,\left(A_i,\mathrm{z}_{im_i}\right) in the ith monitoring well W_i=\left\{A_i\right\}\times\left(-\infty,\hbar(A_i)\right], A_i\in\mathbb{R}^2, i=1,2,\ldots,n. Here \hbar:\mathbb{R}^2\rightarrow\mathrm{Conv}\left(h(\Pi)\right) is the "reasonable" continuation (see stage 3, below) of the function h:\Pi\rightarrow\mathbb{R} to the entire plane \mathbb{R}^2. The exact definitions of the infinitely extended and deep soil S and its surface \partial S are the following:

    \[S=\bigcup\limits_{A\in\mathbb{R}^2}\{A\}\times\left(-\infty,\hbar(A)\right],~\partial S=\left\{\left.\left(A,h(A)\right)\right|A\in\mathbb{R}^2\right\}.\]

Note that \bigcup\limits_{i=1}^nW_i\subset S.


One may cope with the given interpolation problem in several stages. Let's list them.

  1. 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 n.

  2. Merge the duplicate measurement points (for multiple temperature values, use their arithmetic mean as the new temperature). Renumber the updated measurement points in W_i and redefine m_i, i=1,2,\ldots,n.

  3. Construct such a function \hbar:\mathbb{R}^2\rightarrow\mathrm{Conv}\left(h(\Pi)\right) that for all points A\in\mathbb{R}^2, one has


    where A_\Pi is the closest to A\in\mathbb{R}^2 point of \Pi. In particular, if A\in\Pi, then A_\Pi=A and \hbar(A)=h(A).

  4. Construct the depth field d:S\rightarrow\left[0,\infty\right),


    Obviously, d\left(\partial S\right)=\{0\}.

  5. For each nonempty 2-dimensional soil layer S_\mathrm{z}=\left\{\left.A\in\mathbb{R}^2\right|(A,\mathrm{z})\in S\right\} (\left(\mathbb{R}^2\times\{\mathrm{z}\}\right)\cap S\neq\emptyset), obtain a smoothed (with depth-proportional intensity) version \tilde{\hbar}_\mathrm{z}:S_\mathrm{z}\rightarrow\mathrm{Conv}(h(\Pi)),

        \[\tilde{\hbar}_\mathrm{z}=B_{\lambda d(A,\mathrm{z})}\hbar\]

    of the extended elevation field \hbar. Here

        \[\left(B_{\sigma}f\right)(A)\equiv\begin{cases} \frac{1}{2\pi\sigma^2}\displaystyle{\int\limits_{\mathbb{R}^2}}e^{-\frac{\left(|X-A|/\sigma\right)^2}{2}}f(X)dX, & \sigma>0; \\ f(A), & \sigma=0;\end{cases}\]

    B_\sigma is Gaussian blur [2] with the standard deviation (radius) \sigma\geq0. We take

        \[\sigma=\lambda d\left(A,\mathrm{z}\right),~\left(A,\mathrm{z}\right)\in S,\]

    where the physically dimensionless constant \lambda>0 doesn't depend on \left(A,\mathrm{z}\right) and is derived from empirical considerations.

    Note that for all \left(A,\mathrm{z}\right)\in\partial S


  6. Construct the smoothed depth field \tilde{d}:S\rightarrow\mathbb{R},


    For all \left(A,\mathrm{z}\right)\in\partial S


    i.e. \tilde{d}\left(\partial S\right)=\{0\}, just as for d. Where needed, correct the field \tilde{d} in such a way that the strict decreasing of \tilde d\left(A,\mathrm{z}\right) with respect to \mathrm{z}\in\left(-\infty,\hbar(A)\right] would hold for all A\in\mathbb{R}^2, but \tilde{d}\left(\partial S\right)=\{0\} would remain true. From here we obtain that \tilde{d}\geq0.

  7. Compute the smoothed depths


    of the measurement points. These depths are pairwise different within the monitoring well W_i due to the strict monotonicity of \tilde{d}\left(A_i,\mathrm{z}\right) with respect to \mathrm{z}\in\left(-\infty,\hbar(A_i)\right], i=1,2,\ldots,n.

    Construct such linear interpolation (for extrapolation, the nearest neighbor method is used) splines T_i:\left[0,\infty\right)\rightarrow\mathrm{Conv}\left\{T_{i1},T_{i2,},\ldots,T_{im_i}\right\} that

        \[\left[\begin{matrix}T_i\left(\tilde{d}_{i1}) \\ T_i\left(\tilde{d}_{i2}) \\ \vdots \\ T_i\left(\tilde{d}_{im_i}) \end{matrix}\right]=\left[\begin{matrix} T_{i1} \\ T_{i2} \\ \vdots \\ T_{im_i} \end{matrix}\right]\!\!,~i=1,2,\ldots,n.\]

  8. If n=1, put

        \[T\left(A,\mathrm{z}\right)\equiv T_1\left(\tilde{d}\left(A,\mathrm{z}\right)\right)\]

    and stop.

  9. Construct the convex hull C=\mathrm{Conv}\left\{A_1,A_2,\ldots,A_n\right\} and the Delaunay triangulation of the same set of points \left\{A_1,A_2,\ldots,A_n\right\}.

  10. In each fixed point \left(A^*,\mathrm{z}^*\right)\in S (let's denote its smoothed depth \tilde{d}\left(A^*,\mathrm{z}^*\right) with \tilde{d}^*), compute the temperature T^* using the following algorithm. Let A_C^* be the closest to A^* point of the convex hull C. In particular, if A^*\in C, then A_C^*=A^*. If the point A_C^* lies in some non-degenerate Delaunay triangle A_iA_jA_k, i,j,k\in\left\{1,2,\ldots,n\right\} and therefore there exist such uniquely determined scalars (barycentric coordinates) \alpha,\beta,\gamma\geq0 that A_C^*=\alpha A_i+\beta A_j+\gamma A_k and \alpha+\beta+\gamma=1, then put

        \[T^*=\alpha T_i\left(\tilde{d^*}\right)+\beta T_j\left(\tilde{d^*}\right)+\gamma T_k\left(\tilde{d^*}\right).\]

    There's also the case when the Delaunay triangulation is empty (contains no non-degenerate triangles). In this situation (keep in mind that n\geq2), the point A_C^* lies on some non-degenerate segment A_iA_j (i,j\in\{1,2,\ldots,n\}), i.e. A_C^*=\alpha A_i+\beta A_j, \alpha+\beta=1, where the choice of the scalars \alpha,\beta\geq0 is unique. Put

        \[T^*=\alpha T_i\left(\tilde{d^*}\right)+\beta T_j\left(\tilde{d^*}\right).\]

It is clear from construction that isosurfaces of the resulting temperature field T adapt to the terrain, but, as the depth d increases, this adaptation weakens monotonically.


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

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

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

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>