Signal processing for in situ detection of effective heat pulse probe spacing radius as the basis of a self-calibrating heat pulse probe

A sensor comprised of an electronic circuit and a hybrid single and dual heat pulse probe was constructed and tested along with a novel signal processing procedure to determine changes in the effective dual-probe spacing radius over the time of measurement. The circuit utilized a proportional–integral–derivative (PID) controller to control heat inputs into the soil medium in lieu of a variable resistor. The system was designed for onboard signal processing and implemented USB, RS-232, and SDI-12 interfaces for machine-to-machine (M2M) exchange of data, thereby enabling heat inputs to be adjusted to soil conditions and data availability shortly after the time of experiment. Signal processing was introduced to provide a simplified single-probe model to determine thermal conductivity instead of reliance on late-time logarithmic curve fitting. Homomorphic and derivative filters were used with a dual-probe model to detect changes in the effective probe spacing radius over the time of experiment to compensate for physical changes in radius as well as model and experimental error. Theoretical constraints were developed for an efficient inverse of the exponential integral on an embedded system. Application of the signal processing to experiments on sand and peat improved the estimates of soil water content and bulk density compared to methods of curve fitting nominally used for heat pulse probe experiments. Applications of the technology may be especially useful for soil and environmental conditions under which effective changes in probe spacing radius need to be detected and compensated for over the time of experiment.

1 Introduction

The heat pulse probe (HPP) is widely used to determine the thermal conductivity (Abu-Hamdeh, 2001; Abu-Hamdeh and Reeder, 2000; Jin et al., 2017; Li et al., 2016; Liu et al., 2007; Ochsner and Baker, 2008; Penner, 1970; Yun and Santamarina, 2008), thermal diffusivity, and heat capacity (Bristow, 1998; Ham and Benson, 2004; Kluitenberg et al., 1993; Liu et al., 2007; Ochsner et al., 2001; Zhang et al., 2014) of soil. HPPs have been used to measure the thermal conductivity (Morin et al., 2010; Sturm and Johnson, 1992) and density of snow (Liu and Si, 2008); a review is presented by Kinar and Pomeroy (2015). For soils, HPP measurements provide inputs for mathematical models used to determine volumetric water content (Basinger et al., 2003; Bristow, 1998; Bristow et al., 1993; Ham and Benson, 2004; Heitman et al., 2003; Li et al., 2016; Song et al., 1998) and water flux (Hopmans et al., 2002; Kamai et al., 2008; Mori et al., 2003; Wang et al., 2002). A comprehensive review of HPP sensors used to measure water flux is given by He et al. (2018). Installed into a tree trunk (Green et al., 2003) or plant stem (Miner et al., 2017), HPPs can measure sap flow rates. Multifunctional HPPs can simultaneously measure soil thermal and electrical properties to determine soil water retention and hydraulic conductivity (Bristow et al., 2001; Mori et al., 2003; Valente et al., 2006). Data from HPPs can be used to drive predictive mathematical models for water transport (Liu and Si, 2008; Saito et al., 2006; Trautz et al., 2014) and snowpack evolution (Ochsner and Baker, 2008). These models are also useful for civil and geological engineering applications (Ochsner et al., 2001), as well as for prediction of runoff (Yang and Jones, 2009), agricultural productivity (Pearsall et al., 2014; Sturm and Johnson, 1992), climate change, and avalanche hazards (Morin et al., 2010).

HPPs can be broadly classified into two different types: single-probe (SP) (De Vries, 1952; Li et al., 2016) and dual-probe (DP) (Bristow et al., 1993; Campbell et al., 1991; Ham and Benson, 2004) devices. The SP consists of a single heater needle that is inserted into the geomaterial. A temperature measurement sensor (i.e., a thermistor) placed inside the heater needle is used to determine the change in temperature as the needle releases thermal energy. The DP consists of two needles that are inserted into the geomaterial: one of the needles functions as a heater, whereas the other needle measures the change in temperature of the geomaterial at an offset distance to the heated needle. The DP has an advantage over the SP since the SP can only be used to determine thermal conductivity, whereas the DP can be used to determine the thermal conductivity and diffusivity of the geomaterial (Bristow et al., 1994).

An assumption nominally made in conjunction with DP sensors is that the radius is constant during each measurement. If the DP radius changes after the probe is inserted into the geomaterial or over the time of measurement due to heating and cooling, HPP determination of thermal properties will be inaccurate (Kluitenberg et al., 1993; Mori et al., 2003). The measured thermal conductivity is not sensitive to changes in DP probe spacing radius, whereas the HPP-determined heat capacity and thermal diffusivity exhibit high sensitivity to radius changes (Kluitenberg et al., 2010; Liu et al., 2007). This creates challenges in estimating the moisture content of frozen soils wherein thawing and freezing occur and has required recalibration of individual probes (Zhang et al., 2011).

Most HPP researchers utilized commercially available off-the-shelf (COTS) hardware (i.e., a datalogger) to collect data (Bristow, 1998; Bristow et al., 1994, 2001; Kamai et al., 2008; Li et al., 2016), although recently custom electronic circuits have been proposed. Valente et al. (2006) interfaced a multifunctional soil probe to a processing circuit. Dias et al. (2013) used an NPN transistor as a heat source for an SP device. The temperature of the transistor was determined using a circuit and transistor circuit theory. Sherfy et al. (2016) used an NE555 timer circuit to control the duration of the heating pulse. Miner et al. (2017) and Ravazzani (2017) developed Arduino-based HPP sensors utilizing currently established DP theory. Liu et al. (2013), Wen et al. (2015), and Liu et al. (2016) showed that two or more thermistors placed inside the temperature measurement needles of a DP device can be used to determine probe deflection. Multiple thermistors are required to determine probe deflection, and the method cannot be used to calculate a time series of small changes in the probe spacing radius that occur during the time of measurement when the heater needle increases in temperature.

A self-calibrating heat pulse probe (SCHEPP) system is described that consists of a custom electronic circuit and novel inverse models for the SP and DP. The HPP used for SCHEPP is a hybrid of the SP and DP designs. SP and DP forward models are combined and used to determine changes in the effective probe spacing radius during the time of measurement. This effective radius compensates for model error and, similar to a calibrated probe spacing radius, does not directly coincide with the actual probe spacing radius. Another inverse model is also introduced that allows for the determination of thermal conductivity without the need for an SP model late-time approximation (see Sect. 2.2.1 for the rationale).

2 Materials and methods

Diagrams of the SCHEPP system HPP are shown as Figs. 1a and 2a. A loop of nichrome wire is placed inside a heater needle, along with a measurement thermistor. Another measurement thermistor is placed inside a temperature-sensing needle situated at an offset distance to the heater needle. Figure 2a shows that SCHEPP uses a hybrid SP and DP device. A proportional–integral–derivative (PID) controller is used to precisely control and maintain heat inputs in lieu of a variable resistor. Circuit theory is used to determine the resistance of the nichrome heating wire inside the heater needle during a measurement. This eliminates the need to use a previously measured estimate of the heater wire resistance. The heater wire resistance is directly measured over the time of an experiment. Figure 1b is a conceptual diagram that shows relationships between the models and measurement methods.

Figure 1 (a) Conceptual block diagram of the HPP system. The set point q and the PID controller ensure a constant input of heat into the soil over the length of the heater needle. The PID controller modulates heat inputs by changing the output of a variable DC voltage source, and the feedback path is shown. A Kelvin-connection sense resistor measures the current through the nichrome heater wire. A reduction in voltage over the sense resistor element is ΔE, and the ground-referenced output voltage through the nichrome wire is E_kn. The heater needle temperature sensor and offset needle temperature sensor circuits are also shown. A constant voltage source of 2.5 V is connected to a half-bridge. One element of the bridge is a thermistor with a nominal resistance of 10 kΩ specified at a temperature of 20 ^∘C, and the other is a precision resistor with a fixed resistance. (b) Conceptual diagram showing relationships between models and heat pulse probe (HPP) measurements for the hybrid single and dual probe. Each model in the diagram is described in the associated text. Model text in bold indicates signal processing introduced in this paper. The r(t) indicates a time-variable effective radius.

Figure 2 Diagrams showing the system and microcontroller. (a) The circuit operation is controlled by a 32 bit microcontroller clocked at 300 MHz by a PLL. The microcontroller and associated transceiver circuity implement M2M communications whereby data and system operation can be exchanged between machines using commands sent over USB, RS-232, or SDI-12. The microcontroller acts as a state machine that samples data from the AFE and stores the data in SDRAM where signal processing is conducted. The circuit is powered by a nominal 12 V DC supply that is reduced to 3.3 V by a DC–DC switcher. The HPP is a hybrid SP and DP design comprised of a heater needle and a sense needle. The effective distance between the needles as a function of time is r(t). (b) Picture of the circuit. (c) Soil column used in the experiment. (d) Experimental setup. (e) Image of needle probe prototype.

2.1 Forward models

Assuming that the heater needle is an infinite line source in an infinite medium, the “late-time” change in temperature ΔΓ₁(t) of the heater probe SP device is given by Eq. (22a) of Blackwell (1954):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{1}}\left(t\right)={\displaystyle \frac{q}{\mathrm{4}\mathit{\pi }k}}\log \left(t\right)+B+{\displaystyle \frac{\mathrm{1}}{t}}\left(\left(C\right)\log \left(t\right)+D\right),t>\mathrm{0}.\end{array}$$

In Eq. (1) above, q is the rate of energy transferred per unit length of the probe, k is the thermal conductivity of the geomaterial, B, C, and D are constants, and the natural logarithm is utilized. The assumption of an infinite line source in an infinite medium is valid if the heater needle has a small diameter and the geomaterial is of sufficiently large dimension to be isotropic and homogeneous throughout so that the heat pulse does not interact with dissimilar boundaries (i.e., a container in which the soil is placed) over the time of the measurement (Kluitenberg et al., 1993, 1995; Liu et al., 2007). For $t\gg r_{\mathrm{n}}^{\mathrm{2}}/\mathit{\alpha }$, where r_n is the radius of the needle and α is the thermal diffusivity of the medium, the last term in Eq. (1) can be neglected (Bristow et al., 1994; Li et al., 2016). In this paper, Eq. (1) as a forward model is taken subject to the constraint that ΔΓ₁(t)>0 since negative values are not physically reasonable within the context of the model.

Assuming an infinite line source in an infinite medium for a DP device, the change in temperature γ₁(r,t) sensed at a radial distance r from the heater needle is as follows (Kluitenberg et al., 1993).

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}(r,t)=\left\{\begin{array}{cc}\mathrm{\Delta }{T}_{\mathrm{1}}(r,t),& t_{\mathrm{0}}<t\le t_{\mathrm{h}}\\ \mathrm{\Delta }{T}_{\mathrm{2}}(r,t),& t>t_{\mathrm{h}}\end{array}\right.\text{(3)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{1}}(r,t)={\displaystyle \frac{-q}{\mathrm{4}\mathit{\pi }k}}{E}_{\mathrm{i}}\left({\displaystyle \frac{-r^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }t}}\right),t>\mathrm{0}\text{(4)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{2}}(r,t)={\displaystyle \frac{q}{\mathrm{4}\mathit{\pi }k}}\left[E_{\mathrm{i}}\left({\displaystyle \frac{-r^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }\left(t-t_{\mathrm{h}}\right)}}\right)-E_{\mathrm{i}}\left({\displaystyle \frac{-r^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }t}}\right)\right],t>\mathrm{0}\text{(5)}& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{k}{\mathit{\rho }c}}\end{array}$$

The thermal diffusivity of the medium is α, the density is ρ, and the specific heat capacity is c. The exponential integral function is E_i. The current through the nichrome wire is turned on at time t₀ and turned off at time t_h. Therefore, $t_{\mathrm{0}}<t\le t_{\mathrm{h}}$ is referred to as the heating period and t>t_h as the cooling period. Calibration to determine a radius using least-squares curve fitting will yield an effective radius that is representative of differences between the sensing system and the ideal model described above. This initial radius is referred to as r_initial and is taken as a constant.

2.2 Inverse models

2.2.1 Thermal conductivity

Curve fitting using Eq. (1) can be conducted for the section of the heating curve where $t\gg r_{\mathrm{n}}^{\mathrm{2}}/\mathit{\alpha }$. However, the late-time approximation with $t\gg r_{\mathrm{n}}^{\mathrm{2}}/\mathit{\alpha }$ increases the time of measurement and necessitates that the 1∕t term is negligible. When $t<r_{\mathrm{n}}^{\mathrm{2}}/\mathit{\alpha }$ and for suitable t such that the model of Eq. (1) is valid, the constants $\mathit{\{}k,B,C,D\mathit{\}}$ are difficult to directly determine using curve fitting by optimization whereby multiple values can be found to appropriately fit the same model of Eq. (1).

Given these disadvantages, this section introduces a method, the signal processing SP model, which uses signal processing to reduce the time series associated with Eq. (1) to a simpler model. Least-squares curve fitting is used with a modified version of Eq. (1) and the total SP dataset during heating. Errors introduced during the earlier time of heating are acceptably small, particularly when signal processing has modified the SP dataset as a time-domain signal and k is obtained by a least-squares curve fit.

Equation (1) is subjected to Hadamard (point-by-point) multiplication by t to obtain

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{2}}\left(t\right)=t\left[\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{1}}\left(t\right)\right]={\displaystyle \frac{qt}{\mathrm{4}\mathit{\pi }k}}\log \left(t\right)+Bt+\left(C\right)\log \left(t\right)+D.\end{array}$$

The Hadamard multiplication to produce Eq. (6) is a type of homodyning process whereby later-time values of ΔΓ₁(t) are assigned greater-magnitude weights than earlier-time values. Taking the numerical time derivative is similar to the application of a high-pass filter (Hamming, 1983; pg. 118). The resulting equation is

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{3}}\left(t\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{2}}\left(t\right)={\displaystyle \frac{q}{\mathrm{4}\mathit{\pi }k}}\left[\log \left(t\right)+\mathrm{1}\right]+B+{\displaystyle \frac{C}{t}}.\end{array}$$

Homodyning again by t yields

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{4}}\left(t\right)=t\left[\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{3}}\left(t\right)\right]={\displaystyle \frac{qt}{\mathrm{4}\mathit{\pi }k}}\left[\log \left(t\right)+\mathrm{1}\right]+Bt+C.\end{array}$$

Taking the numerical time derivative again,

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{5}}\left(t\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\mathrm{\Delta }{\mathrm{\Gamma }}_{\mathrm{4}}\left(t\right)\right]={\displaystyle \frac{q}{\mathrm{4}\mathit{\pi }k}}\left[\log \left(t\right)+\mathrm{2}\right]+B.\end{array}$$

To reduce noise associated with the derivative operation when working with actual data, a Butterworth low-pass filter with zero-phase filtering and a cutoff frequency of 0.3 Hz is applied to the numerical sequence associated with Eq. (9). The cutoff frequency was chosen to ensure stability of the inverse model within the context of the data used for the experiments reported in this paper. Given a known q, curve fitting is applied to the filtered sequence associated with Eq. (9) to determine k without the need to also determine {C,D} at time $t<r_{\mathrm{n}}^{\mathrm{2}}/\mathit{\alpha }$, where the model of Eq. (1) is still valid. Once the thermal conductivity k is determined by the signal processing, curve fitting using Eq. (1) with a starting value of the determined thermal conductivity is used to estimate parameters for the application of the SP forward model.

The late-time SP model is used for comparison with the signal processing SP model described above. A linear section of the ΔΓ₁(t)−log (t) curve is identified as the time t_d that is 1 s after the time when the second numerical derivative of the ΔΓ₁(t)−log (t) curve is approximately equal to zero. The time t_d is identified using a signal processing zero-crossing detector, and the 1 s delay is introduced to ensure that the curve is approximately linear. Curve fitting using Eq. (1) with $C=D=\mathrm{0}$ is then used to determine k from the linear section of the ΔΓ₁(t)−log (t) curve (Li et al., 2016).

2.2.2 Dual probe and variable radius

The inverse model described in this section uses signal processing to determine r(t) as an effective radius that changes over the time of heating and cooling to compensate for model error and physical changes in the probe spacing. A numerical value for r_initial is required as an estimate of the initial probe spacing radius. Changes in the effective probe spacing r(t) are determined over time and used to obtain an α that is representative of these changes. This is the signal processing DP model.

Thermal conductivity, k, is determined using a measured q and the inverse signal processing model presented in Sect. 2.2.1 for the SP. The thermal conductivity is not directly determined from the DP model using curve fitting since the effective radius r(t) can change over time.

The k is used with a known q to algebraically remove the q∕4πk term from Eq. (2) to obtain γ₂(r(t),t) as an expression written only in terms of E_i. The inverse of the exponential integral is determined using the procedure in Appendix A suitable for an embedded system.

From the inverse, we determine

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{3}}\left(r\left(t\right),t\right)={\displaystyle \frac{-{\left(r\left(t\right)\right)}^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }}},\text{(11)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{4}}\left(r\left(t\right),t\right)=-{\mathit{\gamma }}_{\mathrm{3}}\left(r\left(t\right),t\right)={\displaystyle \frac{{\left(r\left(t\right)\right)}^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }}}.\end{array}$$

The square root function transforms Eqs. (11) and (12) below. To ensure the application of the square root function with real numbers, the numerical implementation must ensure ${\mathit{\gamma }}_{\mathrm{4}}\left(r\left(t\right),t\right)>\mathrm{0}$.

$$\begin{array}{}\text{(12)}& {\mathit{\gamma }}_{\mathrm{5}}\left(r\left(t\right),t\right)={\left[{\mathit{\gamma }}_{\mathrm{4}}\left(r\left(t\right),t\right)\right]}^{\mathrm{1}/\mathrm{2}}={\displaystyle \frac{r\left(t\right)}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{1}/\mathrm{2}}}}\end{array}$$

Taking the logarithm of Eq. (12) results in Eq. (13) below. Applied to each element of the corresponding sequence, this operation is analogous to homomorphic filtering (Oppenheim et al., 1976).

$$\begin{array}{}\text{(13)}& {\mathit{\gamma }}_{\mathrm{6}}\left(r\left(t\right),t\right)=\log \left({\mathit{\gamma }}_{\mathrm{5}}\left(r\left(t\right),t\right)\right)=\log \left(r\left(t\right)\right)-\log \left(\mathrm{2}{\mathit{\alpha }}^{\mathrm{1}/\mathrm{2}}\right)\end{array}$$

Taking the time derivative of Eq. (13) is once again similar to the application of a high-pass filter that suppresses the constant log (2α^1∕2) term. When Eq. (14) is expressed as a discrete sequence sampled at a frequency of f_s the derivative is approximated using a backward-difference method (Eq. 15).

$$\begin{array}{}\text{(14)}& {\mathit{\gamma }}_{\mathrm{7}}\left(r\left(t\right),t\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathit{\gamma }}_{\mathrm{6}}\left(r\left(t\right),t\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\log \left(r\left(t\right)\right)\end{array}$$

The derivative is computed using a backward difference:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathit{\gamma }\left(r\left(t\right),t\right)\approx {\displaystyle \frac{a_{i-\mathrm{1}}-a_i}{\mathrm{\Delta }t}}.\end{array}$$

In Eq. (15) above, a_i=log (r(t_i)), where the index i denotes the element of a discrete sequence and Δt is the time step calculated by $\mathrm{\Delta }t=\mathrm{1}/f_{\mathrm{s}}$. To reduce numerical error at small time t values, the inverse model described in this section is applied with a time step Δt such that the numerical inverse of E_i(x) can be successfully computed (Sect. 2.5).

The solution of Eq. (15) requires a boundary condition a_b=log (r(t_initial)), where t_initial is the time at which the probe spacing radius is r_initial. For application to actual data, we assume that the probe spacing is the calibrated r_initial at $t_{\text{initial}}=t_{\mathrm{p}}+t_{\mathrm{a}}$, where t_p is the time at which the curve associated with γ₁(r,t) is at a maximum and t_a is an additional time delay that compensates for a nonideal system. For this system, the additional time delay was chosen such that $\mathrm{0}\mathrm{s}\le \left|t_{\mathrm{p}}\right|\le \mathrm{2}\mathrm{s}$. Selection of a peak time is a similar idea to the temperature maximum method (Bristow et al., 2001) whereby the calibrated r_initial is used at the time of peak temperature change. The additional time delay t_a is chosen to approximately coincide with the integer-valued time durations of moving-average windows (Sect. 2.5). The input data are trimmed appropriately. This selection of boundary condition is supported by tests on actual soil performed in this paper and a sensitivity analysis that justifies the selection of the additional time t_a (Sect. 3.3).

After r(t) is determined, the α is determined by taking the average of γ₈(r(t),t) over time. Evaluation of Eq. (16) will thereby yield a curve that is a straight line with a slope that is approximately zero when α is approximately constant over the time of heating and cooling. The r(t) is an effective radius that is also affected by temperature drift and deviation of the physical system from an ideal model. Since r(t) is an effective radius, it will not directly coincide with an actual probe spacing radius.

$$\begin{array}{}\text{(16)}& {\mathit{\gamma }}_{\mathrm{8}}\left(r\left(t\right),t\right)={\displaystyle \frac{{\left(r\left(t\right)\right)}^{\mathrm{2}}}{\mathrm{4}{\mathit{\gamma }}_{\mathrm{4}}\left(r\left(t\right),t\right)}}\approx \mathit{\alpha }\end{array}$$

2.3 Measurement of soil water content and density

The heat capacity k and thermal diffusivity α are determined using an inverse model as described in the previous sections of this paper. Neglecting the contribution of air, the volumetric heat capacity of soil C_h is calculated by Kluitenberg (2002):

$$\begin{array}{}\text{(17)}& C_{\mathrm{h}}={\mathit{\theta }}_{\mathrm{m}}{C}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}{C}_{\mathrm{o}}+{\mathit{\theta }}_{\mathrm{w}}{C}_{\mathrm{w}}.\end{array}$$

Rearranging Eq. (17) and solving for volumetric water content yields the following.

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\mathit{\theta }}_{\mathrm{w}}& ={\displaystyle \frac{\mathrm{1}}{C_W}}\left[C_{\mathrm{h}}-\left({\mathit{\theta }}_{\mathrm{m}}{C}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}{C}_{\mathrm{o}}\right)\right]\\ C_{\mathrm{h}}& =k/\mathit{\alpha }\\ C_{\mathrm{m}}& =c_{\mathrm{m}}{\mathit{\rho }}_{\mathrm{m}}\\ C_{\mathrm{o}}& =c_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}\end{array}\end{array}$$

In the equations above, C_m, C_o, and C_w are the volumetric heat capacities of mineral content, organic content, and water; θ_m, θ_o, and θ_w are the associated volume fractions; and C_h is the total volumetric heat capacity of the soil. The c_m and c_o are the specific heat capacities of the mineral and organic content, and the ρ_m and ρ_o are the associated densities. Experiments in which numerical values of θ_w<0 or θ_w>1 are not valid and indicate improper contact between the probe and the soil medium. The θ_w is numerically constrained to be within the range $\mathrm{0}\le {\mathit{\theta }}_w\le \mathrm{1}$. The mineral content of the soil θ_m is known, and the organic content θ_o can be easily determined from laboratory testing or an organic carbon soil map of a geographic area. For implementation using a microcontroller, θ_m and θ_o are stored in flash (non-volatile) memory, and these values change based on the geographic location of the soil. The density of the soil is determined by volume fractions:

$$\begin{array}{}\text{(19)}& \mathit{\rho }={\mathit{\theta }}_{\mathrm{m}}{\mathit{\rho }}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}{\mathit{\rho }}_{\mathrm{o}}+{\mathit{\theta }}_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}.\end{array}$$

In Eq. (19) above, the constituent densities and heat capacities are known. The θ_w is determined using Eq. (18).

2.4 Circuit theory and PID control

Figure 1 shows a conceptual block diagram of the system. Thermistors in half-bridge configurations are used to determine the temperatures of HPP needles.

The heat input into the soil by the heater probe is

$$\begin{array}{}\text{(20)}& {\displaystyle }P=I^{\mathrm{2}}{R}_{\mathrm{w}},\text{(21)}& {\displaystyle }q={\displaystyle \frac{P}{\mathrm{\ell }}}.\end{array}$$

In the above Eqs. (20) and (21), the electrical power is P, the total resistance of the nichrome heater wire is R_w, and ℓ is the length of the heater needle. Given a measured voltage drop ΔE over a four-terminal Kelvin sense resistor with known resistance R_s, the current though the heater wire is calculated using Ohm's law:

$$\begin{array}{}\text{(22)}& I={\displaystyle \frac{\mathrm{\Delta }E}{R_{\mathrm{s}}}}.\end{array}$$

For a current I through the heater wire and sense resistor, the output voltage is measured as E_kn by an analog-to-digital converter (Fig. 1). Using Kirchhoff's voltage law for this circuit, the resistance R_w of the nichrome wire is determined at each sampling time step by

$$\begin{array}{}\text{(23)}& R_{\mathrm{w}}={\displaystyle \frac{E_{\text{kn}}-\mathrm{\Delta }E}{I}}.\end{array}$$

To set a constant q, a proportional–integral–derivative (PID) controller (Ang et al., 2005) is utilized. The variable voltage source is adjusted at each discrete time step by a digital-to-analog-converter (DAC). Since ΔE and E_kn are measured at each discrete time step at a sampling rate of f_s, Eqs. (20) to (23) are used with the feedback loop shown in Fig. 1 to ensure that the q remains close to a set-point value during the time of experiment. The use of the PID controller requires a higher sampling rate f_s than a nominal HPP experiment to adjust the output q. The PID controller thereby ensures that the soil can heat up in a controlled fashion and considers resistance changes in the nichrome wire in lieu of using an assumed resistance. Figure 2a is a block diagram indicating how the system incorporates a microcontroller and communication interfaces.

2.5 Determination of temperature change curves

The sampled temperature inside the heater needle is denoted as Γ(t), and the sampled temperature inside the second needle at an offset distance from the heater needle is denoted as T(t). The sampled temperatures Γ(t) and T(t) are low-pass filtered using a fifth-order Butterworth filter with a cutoff frequency of 10 Hz applied as a zero-phase filter to reduce noise. The Butterworth filter was chosen since it is maximally flat in the passband, and the zero-phase filtering ensures that time shifts are minimized to ensure accurate application of the inverse models described in this paper using the collected data.

For the DP model calibration to find an initial r₀ using curve fitting, the sampled temperatures are processed by a moving-average filter over 1 s windows to further reduce noise before curve fitting. Alternately, for the DP inverse model (Sect. 2.2.2), a moving-average filter is used to obtain an equivalent sampling rate of 12 Hz to ensure that the E_i(x) inverse can be accurately computed using floating-point number representations.

The respective averaged initial temperatures of the needles before heating are determined as Γ_av and T_av. Therefore, the temperature changes are calculated as $\mathrm{\Delta }\mathrm{\Gamma }\left(t\right)=\mathrm{\Gamma }\left(t\right)-{\mathrm{\Gamma }}_{\text{av}}$ and $\mathrm{\Delta }T\left(t\right)=T\left(t\right)-T_{\text{av}}$ after the application of any initial filtering. The temperature changes are used for the application of inverse models related to Eqs. (1) and (2).

2.6 Determination of heat inputs

Heat inputs into the soil are determined during the time of experiment and calculated using Eq. (21). When the current is applied and travels through the heater wire, there is a short time delay (≪1 s) before the set point q is attained when the heater needle increases in temperature. The interval of the time series for a constant calculated q during the time of experiment is determined using a step detection algorithm (Carter et al., 2008) based on the Student's t test (Ebdon, 1991; pp. 61–64) with a null hypothesis at a significance level of 1 % and a window size of 31 elements. The significance level and window size are dependent on the implementation of the sampling system and are thereby chosen to detect the step within the context of this experiment. To find an estimate of a constant value of q, the time series is averaged over the plateau of the step. The time series location of the step associated with a constant q is found by application of a sliding mean filter with a window size of 31 elements applied to a time series of electrical power P used for the computation of q using Eq. (21). The mean filter is applied to a binary sequence created by mapping non-rejection of the null hypothesis to binary 0 and rejection of the null hypothesis to binary 1. The plateau is coincident with a sequence of zeros away from the edges of the step. The edges of the step are indicated by nonzero elements in this sequence surrounded by zeros. The window size is appropriate for the sampling system described in the context of this paper.

2.7 Apparatus

A custom electronic circuit board was designed and constructed for the SCHEPP system (Fig. 2a and b). The circuit board was placed into an enclosure box and connected to the HPP by a cable and mating circular connectors. The HPP body was epoxied into a circular hole cut in the bottom of a cylindrical polyvinyl chloride (PVC) container of 10.0 cm diameter and 10.5 cm height that held 825 cm³ of soil. The HPP needles thereby protruded into the soil column formed by the PVC container (Fig. 2c).

Figures 1 and 2a graphically show the different subsystems of the PCB (Fig. 2b). In the analog front end (AFE), a two-channel 24 bit ADC with a precision 10 k (0.01 %, ± 5 ppm ^∘C) resistor half-bridge for each channel and a 2.5 V voltage reference (2 ppm ^∘C, ± 0.02 % voltage error) was used to determine the resistance of the thermistors inside each needle. The resistance of a thermistor was related to temperature by the Steinhart–Hart equation (Steinhart and Hart, 1968).

The variable voltage source was constructed from a DC–DC switcher for efficiency and thereby reduced power consumption in lieu of a linear regulator. The DC–DC switcher could be turned off for an output voltage of 0 V or turned on and adjusted from ∼ 0.49 to ∼ 8.965 V using a DAC that injected current into the feedback loop of the switcher. With design criteria including a heater needle length of ℓ=3.0 cm and a nichrome wire resistance of R_w≈34 Ω, the heater output was limited by software to be within the range of q≈0.24 to $q\approx \mathrm{79}\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}$. The precision sense resistor was chosen as R_s=0.01 Ω (0.1 %, ± 15 ppm ^∘C) to reduce the voltage drop over this circuit element and to ensure that the output voltage could be accurately adjusted. This resistance is smaller than the 1 Ω current sense resistor nominally used in other HPP experiments (Bristow et al., 1994; Li et al., 2016; Liu and Si, 2011, 2008; Valente et al., 2006; Zhang et al., 2011). Moreover, the precision sense resistor had a Kelvin terminal connection for precision and was physically large to reduce self-heating by current flow. The voltage drop over a sense resistor was determined by a precision difference amplifier and a 16 bit ADC, allowing for a 1 LSB step size of 2.5 µV. The output voltage E_kn (Fig. 1) was also measured by a 16 bit ADC and amplifier, resulting in a 1 LSB step size of 1.25 mV.

A 32 bit microcontroller with a system clock of 300 MHz was used to control the HPP experiment and perform floating-point calculations (Fig. 2a). The system clock had to be set at 300 MHz to allow the microcontroller to sample all ADCs in the system at f_s=120 Hz and also perform floating-point calculations associated with this application. The system clock speed is provided here to provide a starting point for engineering of similar designs. The DAC used to control the output voltage was also updated at the same sampling rate with the PID controller output. The 120 Hz sampling rate enabled functioning of the PID controller feedback loop and allowed for digital filtering for signal processing.

The microcontroller had an integrated USB transceiver for communication with a computer. RS-232 and SDI-12 interfaces were also integrated into the system for communication with a computer or datalogger as machine-to-machine (M2M) interfaces. SDRAM stored data from the experiment and provided temporary memory for heap allocation of arrays and data structures. Code for the microcontroller was written in the C programming language.

A command-line serial port interface permitted changing the duration of the experiment, the set-point q value, and the time of heating. For each experiment, the microcontroller monitored the maximum temperature rise at the heater needle and terminated the experiment if the temperature rise exceeded the maximum operating temperature of the thermistors.

The mechanical construction and design of the HPP used for this paper have been reported and rationalized in other papers (Li et al., 2016; Liu and Si, 2008, 2010). The needles of length ℓ=3.0 cm were constructed from stainless-steel tubing (1.28 mm o.d. and 0.84 mm i.d.) and filled with thermally conductive epoxy. The sense thermistor was placed in the geometric middle of each needle to prevent edge effects associated with heat conduction, and the needles were filled with thermal epoxy (Saito et al., 2007). The nominal spacing between the heater needle and the sense needle was 6 mm. During laboratory testing of SCHEPP (Fig. 2c and d), the experiment was initiated by a laptop computer connected to the circuit's USB port, and communication was conducted over the USB interface.

Table 1 Quantities utilized for sand and peat HPP experiments.

2.8 Data collection

Following Campbell et al. (1991), calibration was conducted to find r_initial using a 5 g L⁻¹ agar gel solution. The thermal conductivity k of the agar gel was taken to be the same as the thermal conductivity of water (Saito et al., 2007). Reported values for thermal conductivity (Ramires et al., 1995) and heat capacity (Wagner and Pruß, 2002) of water were used.

For all experiments, the temperature of the probe needles was measured for 1 s at a sampling rate of f_s=120 Hz before electrical current was applied to the nichrome wire. This initial temperature measurement for each trial was averaged over the 1 s period.

The agar gel was washed out using distilled water, and the cylindrical container was packed with soil. Two types of soil were used for the HPP tests: sand and peat. These soils are indicative of the physical extent of soil thermal properties. The soils were collected from field sites near Fort McMurray, Alberta. The sand contained small amounts of bitumen as representative of the Alberta Oil Sands area. Two independent laboratory analyses with incineration at 1100 ^∘C were conducted on the sand, finding the total carbon content to range between a mean of 0.44 % and 0.75 % by mass. The soil properties are summarized in Table 1.

The water content for the sand was chosen so that the sand was saturated, whereas the water content for the peat was chosen so that the soil would remain as wet as possible (Table 1). Due to the absorbent characteristics of the peat soil, it was not possible during the time of the laboratory experiment to completely saturate the pore spaces of the soil column. However, the volumetric water content θ for both soils was chosen to ensure adequate contact between the probe and the soil and also to reduce air gaps that can increase thermal contact resistance and decrease the accuracy of the measurement (Liu and Si, 2010). These air gaps can occur in drier soils with lower water contents. Since comparisons are required to be made between heat pulse probe and gravimetric measurements for testing the in situ calibration procedure described in this paper, the presence of air gaps represents an additional source of error that was controlled.

Table 2 Heat pulse strengths, time of heating and cooling, and total time for each experiment conducted on sand and peat. The experiment identifier is an alphabetical letter that identifies the experiment set. “No.” indicates the total number of experiments conducted per set.

Since the soil dried out over the time of multiple experiments, some additional water was added between successive days to ensure that the volumetric water content θ was close to the target value. Between trials, the top of the container was covered with a cap to reduce evaporation of water from the soil. Changes in water content occurred over the time of the experiment due to evaporation since the cap did not create a hermetic seal between the top of the container and the soil column. Table 1 shows quantities used for the application of HPP forward and inverse models to peat and sand.

Experiment sampling durations, q heat inputs, and heat durations are summarized in Table 2. Trial numbers of each experiment refer to groups of experiments conducted temporally close together.

The heat pulse strength and time of heating were chosen to minimize interaction of the heat pulse with the container boundaries. Due to the short time span over which each experiment was conducted in a laboratory setting, explicit correction was not applied for changes in ambient temperature (Young et al., 2008; Zhang et al., 2014). Between each experiment, the temperature of the soil column returned to a level that approximated the initial temperature before the probe was heated again for the next trial. All experiments were conducted at room temperature (∼ 20 ^∘C).

Numerical comparisons were made using the root mean squared difference (RMSD) and mean bias (MB). The RMSD indicates the overall differences between two datasets. The MB indicates whether the model underpredicts or overpredicts relative to the observations.

3 Results

3.1 Synthetic experiments

Synthetic heating curves were constructed using Eqs. (1) and (2) to serve as a forward model and provide a test of the signal processing. The SP and DP curves are shown as Figs. 3 and 4 and were generated using the model inputs given in Table 3. For the DP, the assumed change in probe spacing radius is shown for a linear increase (Fig. 4b), decrease (Fig. 4e), and Brownian random-walk-scaled so that the numerical values are between a starting and ending radius (Fig. 4h).

Figure 3 Synthetic example of SP signal processing. (a) Forward model and reconstruction of the forward model by the inverse model. (b) Numerical difference between the forward and inverse models. (c) Forward model ΔΓ₅(t) from the signal processing compared to the theoretical model of ΔΓ₅(t). (d) Numerical difference between the forward model ΔΓ₅(t) and the theoretical model.

Figure 4 Synthetic example of DP forward and inverse models. (a, d, g) DP temperature change ΔT with a fixed and variable radius. (b, e, h) Known change in r(t) as a variable radius. The forward and inverse model (reconstruction by signal processing) is shown. (c, f, i) Numerical difference between the forward and inverse models. Each row of Fig. 4 corresponds to a numerical experiment with an associated change in r(t) given in the second column (b, e, h).

Table 3 Synthetic SP and DP model inputs.

The time-variable radius is r(t), and the associated curve is shown on the plots as a DP variable radius. The DP fixed radius curve is calculated using the first element of the r(t) used for a particular DP variable radius curve. The fixed radius is taken to be constant over the time of heating and cooling.

For the SP, Fig. 3a shows the forward model and the reconstruction of the forward model by the inverse model proposed in Sect. 2.2.1. The numerical difference between the forward and inverse models is shown by Fig. 3b and is on the order of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$. This difference occurs due to discretization of the numerical derivatives and floating-point round-off error from the homodyning process. The simplified model given by ΔΓ₅(t) is shown by Fig. 3c and demonstrates the reduction of terms from the original model (Eq. 1). Figure 3d is the numerical difference between the forward and inverse models associated with ΔΓ₅(t). The numerical difference remains small over the time of heating.

Figure 5 Example signal processing for sand showing (a) the measured and modelled SP heating curves; (b) the measured and modelled DP heating curves along with the detected peak time; (c) the detected change in effective radius r(t) from the signal processing. The effective radius r(t) as obtained from signal processing compensates for model error, experimental error, and physical changes in the spacing of the heater needles.

Table 4 Comparisons for sand example.

Figure 4 demonstrates the heating and cooling curves for a DP model with fixed radius and variable radius. The assumed time-variable radius r(t) is given along with the application of the inverse model proposed in Sect. 2.2.2. The first row of Fig. 4 (a–c) is for a linear increase; the second row of Fig. 4 (d–f) is for a linear decrease; and the third row of Fig. 4 (g–i) is for the Brownian random walk. The difference between the forward and inverse models is on the order of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{15}}$ for all changes in probe spacing radius (Fig. 4c–i), indicating that for a synthetic model the error is mostly associated with floating-point calculations and that the inverse model is accurate.

Figure 6 Example signal processing for peat showing (a) the measured and modelled SP heating curves; (b) the measured and modelled DP heating curves along with the detected peak time; and (c) the detected change in effective radius r(t)from the signal processing. The effective radius r(t) as obtained from signal processing compensates for model error, experimental error, and physical changes in the spacing of the heater needles.

Table 5 Comparisons for peat example.

3.2 Soil data

Since the signal processing SP and signal processing DP models are applied together, hereafter these will be referred to as the signal processing SP and DP model. The late-time SP model is described in Sect. 2.2.1. The heating and cooling DP model refers to the nominal curve fitting using Eq. (2).

Figure 5 provides an example of the models applied to sand (Table 4). The peak time corresponding to t_p+t_a is indicated as a vertical line in Fig. 5b. Changes in effective radius r(t) determined by signal processing are shown as Fig. 5c. The rapid fluctuations in the effective radius r(t) occur due to temperature drift, model error, and experimental error. The signal processing thereby compensates for these effects using r(t) as an effective radius. Table 4 shows that the quantities found using all models have the same orders of magnitude and indicates that the signal processing SP and DP model is more accurate than the heating and cooling DP model compared to the gravimetric values used for the laboratory experiment.

Figure 6 shows the signal processing SP and DP model applied to peat (Table 5). Compared to the sand example given above, the thermal conductivity k and diffusivity α are lower for the peat, demonstrating that the peat takes longer to warm up during the time of experiment. The heating and cooling curves are thereby distinctively different between sand and peat. However, in the same manner as the sand example, the quantities found using all models are of similar orders of magnitude. Figure 6c shows that there are fewer rapid fluctuations in the effective radius r(t) determined for peat compared to sand. Moreover, the change in effective radius is less pronounced for peat compared to sand due to smaller temperature drift associated with lower k and α. Also, in a similar fashion to sand, the percentage difference (PD) and numerical difference demonstrate that the signal processing models introduced in this paper are more accurate than the nominal models.

Figure 7 shows the thermal conductivity k and diffusivity α for the sand and peat determined for the experimental trials. For the sand, the nominal heating and cooling DP model provides estimates of k that are mostly higher than the signal processing SP and late-time SP models. The late-time SP model provides estimates of k that are intermediate between the other two models. The signal processing SP model produces the lowest estimates of k. However, the estimates of k made by all three models are the same orders of magnitude and remain relatively constant over experiments conducted on each soil type. The thermal diffusivity estimates provided by the heating and cooling DP model are slightly higher than the estimates provided by the signal processing SP and DP model for sand (Fig. 7c), whereas for peat (Fig. 7d) the thermal diffusivity for both models is approximately similar.

Figure 7 Thermal conductivity k and thermal diffusivity α for each of the experiments. (a) Thermal conductivity for sand. (b) Thermal conductivity for peat. (c) Thermal diffusivity for sand. (d) Thermal diffusivity for peat.

Results corresponding to Figs. 7 an 8 are shown by Tables 6 and 7. These tables show the RMSD, MB, and PD, and the results of these calculations are discussed below.

Table 6 RMSD, MB, and PD comparisons for θ_w.

Table 7 RMSD, MB, and PD comparisons for ρ.

Figure 8 shows the determined values of water content θ_w and density ρ. Overall, for the sand and peat experiments, the nominal heating and cooling DP model has a higher-magnitude RMSD, MB, and PD compared to the signal processing SP and DP model introduced in this paper, indicating that the determination of the effective r(t) by signal processing improves estimates of θ_w and ρ.

Figure 8 Water content θ_w and density ρ for each of the experiments. (a) Water content for sand. (b) Water content for peat. (c) Density for sand. (d) Density for peat.

For all sand experiments, the signal processing SP and DP model reduces the θ_w RMSD by 0.10 (10 %), the MB by 0.33 (33 %), and the PD by 4.2 % compared to the nominal heating and cooling DP model. For all peat experiments, the corresponding reduction for θ_w is 0.07 (7 %) for the RMSD, 0.41 (41 %) for the MB, and 12.3 % for the PD.

For the density ρ of sand, the signal processing SP and DP model reduces the RMSD by 102 kg m⁻³, the MB by 338 kg m⁻³, and the PD by 0.85 %. For the density ρ of peat, the signal processing SP and DP model reduces the RMSD by 94 kg m⁻³, the MB by 407 kg m⁻³, and the PD by 5.1 %.

Figure 9 Sensitivity analysis for water content θ_w for sand and peat with respect to the different models. (a) RMSD for changes in the calibrated initial radius r. (b) MB for changes in the calibrated initial radius. (c) RMSD for changes in the organic volume fraction θ_o. (d) MB for changes in the organic volume fraction θ_o. (e) RMSD for changes in the mineral volume fraction θ_m. (f) MB for changes in the mineral volume fraction θ_m. (g) RMSD for changes in the time delay t_a. (h) MB for changes in the time delay t_a.

Figure 10 Sensitivity analysis for density ρ for sand and peat with respect to the different models. (a) RMSD for changes in the calibrated initial radius r. (b) MB for changes in the calibrated initial radius. (c) RMSD for changes in the organic volume fraction θ_o. (d) MB for changes in the organic volume fraction θ_o. (e) RMSD for changes in the mineral volume fraction θ_m. (f) MB for changes in the mineral volume fraction θ_m. (g) RMSD for changes in the time delay t_a. (h) MB for changes in the time delay t_a.

For experiment identifiers A, B, C, and D associated with sand, the RMSD, MB, and PD magnitudes are lowest for θ_w determined from the signal processing SP and DP model compared to the nominal heating and cooling DP model. The RMSD, MB, and PD magnitudes are also lower for most sand estimates of ρ other than experiment identifier D, for which the RMSD and MB magnitudes are larger than magnitudes associated with the nominal heating and cooling DP model. The heat pulse strength and duration were of greater magnitude for experiment identifier D, but since ρ is being calculated as a function of θ_w this indicates a change in soil constituents (see Eq. 19) in the vicinity of the probe for experiment identifier D and demonstrates sensitivity of the signal processing to this change.

An increase in water content is apparent for experiment identifiers E to G associated with peat. The increase in water content occurred since the HPP experiments were initially conducted less than an hour after water was added to the soil column. Since the soil column was opaque, the infiltration of water in the column and the associated wetting front could not be tracked, and thereby localized volumes of water surrounding the HPP caused a rise in water content. Since the rise is consistent and shown by the heating and cooling DP model as well as the signal processing SP and DP model for θ_w and ρ (Fig. 8b and d), this indicates that both models are in physical agreement. Lower RMSD, MB, and PD values for the signal processing SP and DP model indicates that the signal processing introduced in this paper also improves estimates of θ_w and ρ when the water content changes during the time interval of experiment identifier E.

For experiment identifiers F and G, the determined θ_w and ρ remain approximately constant over time for peat. Compared to experiment identifier E, a reduction in water content has occurred due to infiltration over time and some loss of water due to evaporation. For these experiments, the RMSD, MB, and PD are lower for the signal processing SP and DP model compared to the heating and cooling DP model.

3.3 Sensitivity analysis

To identify the effects of model parameters on the RMSD and MB of model outputs θ_w and ρ, a sensitivity analysis (Figs. 9 and 10) was conducted over all data collected. The sensitivity analysis utilized the OAT (one-at-a-time) approach, whereby one variable at a time is changed, whereas the other model inputs are held constant (Hamby, 1994). Overall, for a range of nominal model inputs, Figs. 9 and 10 demonstrate that the signal processing associated with the signal processing SD and DP model reduces the RMSD and MB compared to the nominal heating and cooling DP model. This also indicates that the signal processing method produces more accurate estimates than the curve-fitting models nominally used for heat pulse probe experiments.

For all models and soils used to determine θ_w and ρ, the RMSD and MB are lowest when the initial radius r_initial is close to the calibrated value, indicating the importance of calibration for all models. If the r_initial is underpredicted, the MB indicates an overprediction of θ_w and ρ, whereas an overprediction of r₀ indicates an underprediction of θ_w and ρ.

For sand θ_w, an organic content θ_o<0.1 produces the lowest-magnitude RMSD and MB, whereas for peat θ_w, an organic content close to θ_o≈0.40 produces the lowest-magnitude RMSD and MB. This physically approximates the composition of the sand and peat soils used for these experiments. The MB for sand and peat indicates a model overprediction for θ_o values lower than these thresholds and a model underprediction for θ_o values higher than these thresholds. A similar effect is also shown for the mineral content θ_m, with θ_m≈0.60 resulting in the lowest-magnitude RMSD and MB for sand. For peat, θ_m≈0.15 results in the lowest-magnitude RMSD and MB.

For sand θ_w the time delay t_a≈1.5 s is a good assumption to provide the lowest values of RMSD and MB for the signal processing SD and DP model for both sand and peat. For peat ρ the time delay is t_a≈0 s. The time delay t_a is not a parameter for the nominal heating and cooling DP model, and a sensitivity analysis is not conducted for t_a when using this model. Applied to peat, the signal processing SD and DP model is relatively insensitive to the time delay t_a due to the lower thermal conductivity k and diffusivity α relative to sand that dampens changes in the effective radius r(t) as determined by signal processing.

In the context of the sensitivity analysis, as θ_o increases for sand, the RMSD and MB related to ρ also increase. For peat, a concomitant increase in θ_o is associated with an increase in the RMSD and MB, indicating that for ρ it is not possible to calibrate for θ_o and an approximation of θ_o must be known for model application within the context of these experiments. A mineral content of θ_m≈0.55 for sand produces the lowest RMSD and MB related to ρ. As θ_m increases for peat, the RMSD and MB also increase, indicating that for ρ it is once again not possible to calibrate for θ_m and an approximation must be utilized.

4 Conclusions

• A novel circuit was designed and tested using a hybrid SP and DP heat pulse probe (HPP) design. The circuit utilized a PID controller to precisely control the heat input into the soil. In lieu of a variable resistor, this enabled the heat input q to be changed by a computer or a datalogger. When deployed at a remote or inaccessible field site, the HPP heat input can be set to a given value using the communication interfaces. This enables the heat input to be appropriately selected for soil type.

• Instead of using a 1 Ω sense resistor to infer heat inputs into the soil, the circuit used a resistor with a 0.01 Ω nominal resistance. This reduced the voltage drop across the sense resistor and still allowed the current through the nichrome wire to be adequately determined during the time of experiment, although a differential amplifier was required to detect the voltage difference before digitization by an ADC.

• A sampling rate of 120 Hz was required for application of the PID controller and theory associated with signal processing related to a hybrid SP and DP heat pulse probe. The higher sampling rate allowed digital filtering to be applied.

• Signal processing was used to determine thermal conductivity using an SP model that did not rely on a late-time SP approximation. A DP model was used to determine changes in the effective DP probe spacing radius.

• The DP and SP signal processing models introduced in this paper improved overall estimates of soil water content θ_w and bulk soil density ρ for sand and peat soils, indicating that the detection of effective changes in the probe spacing radius using signal processing is useful to correct for model error and physical changes in the probe spacing. This improvement is associated with standard HP and DP probes that are used together in a novel fashion along with signal processing.

• Further research is required to test the signal processing models introduced in this paper and to compare the estimates of soil water content θ_w and bulk soil density ρ to estimates made using other measurement systems and technologies. The effective radius calibration has an advantage for soil types wherein expansion and contraction of the soil can cause changes in the effective probe spacing radius r(t).

Appendix A

To efficiently obtain the inverse of E_i(x), endpoints for a search interval need to be appropriately selected, particularly when the inverse model runs on an embedded resource-constrained microcontroller. To ensure numerical continuity and accuracy between the forward and inverse models, the same E_i(x) function is used in the inverse of E_i(x) in lieu of alternative numerical approximations.

Take $E_{\mathrm{1}}\left(x\right)=-E_{\mathrm{i}}(-x)$, where E₁(x)=E_n(x) with n=1. Let x>0 and we need to show

$$\begin{array}{}\text{(A1)}& -{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}<E_{\mathrm{1}}\left(x\right)<{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

Using 5.1.20 of Abramowitz and Stegun (1964), we need to show that with x>0

$$\begin{array}{rl}-{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}& <{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\exp (-x)\log \left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}}{x}}\right)\\ & <E_{\mathrm{1}}\left(x\right)<\exp (-x)\log \left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}\right)<{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

Since $\exp (-x)>\mathrm{0}$ and $\log \left(\mathrm{1}+\frac{\mathrm{2}}{x}\right)>\mathrm{0}$ as x→∞, then $\frac{\mathrm{1}}{\mathrm{2}}\exp (-x)\log \left(\mathrm{1}+\frac{\mathrm{2}}{x}\right)>\mathrm{0}$ as x→∞. Therefore, $-\frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}<\frac{\mathrm{1}}{\mathrm{2}}\exp (-x)\log \left(\mathrm{1}+\frac{\mathrm{2}}{x}\right)$. We next show

$$\exp (-x)\log \left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}\right)<{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}.$$

Algebraic rearrangement yields

$${\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\log \left(x\right)+\log \left(\log \left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}\right)\right)<x.$$

Since x>0, we need to show

$${\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\log \left(x\right)+\log \left(\log \left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}\right)\right)<\mathrm{0}.$$

By algebraic rearrangement, we find

$$\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}<\exp \left({\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}\right).$$

The right-hand side can be replaced by a power series representation:

$$\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}<\sum _{i=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{x^{-i/\mathrm{2}}}{i\mathrm{!}}}.$$

The proof proceeds by contradiction and shows that the inequality holds. We assume

$$\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x}}\ge \mathrm{1}+{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}+\mathrm{\cdots }$$

However, this is a contradiction, so inequality Eq. (A1) holds.

From inequality Eq. (A1), the search interval for E_i(x) is $\left[-\mathrm{1}/{\left(E_{\mathrm{i}}\left(x\right)\right)}^{\mathrm{2}},\mathrm{0}\right)$ for x<0. The inverse of E_i(x) is computed using Golden-section search (Kiefer, 1953) on the bounded interval.

The following shows how endpoints of the search interval for the inverse are selected for the cooling section of the DP curve using the result given above. Without the q∕(4πk) term, Eq. (4) is of the form

$$\begin{array}{}\text{(A2)}& f(x,y)=f\left(h\right)=E_{\mathrm{i}}(-x)-E_{\mathrm{i}}(-y).\end{array}$$

Let $x=h/\left(t-t_{\mathrm{h}}\right)$ and $y=h/t$, where $h={\left(r\left(t\right)\right)}^{\mathrm{2}}/\left(\mathrm{4}\mathit{\alpha }\right)$ and x>0 and y>0. Equation (A2) can be rewritten as

$$f\left(h\right)=f(x,y)=-E_{\mathrm{1}}\left(x\right)+E_{\mathrm{1}}\left(y\right).$$

From the above and inequality (A1),

$${\displaystyle \frac{-\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{y^{\mathrm{1}/\mathrm{2}}}}<-E_{\mathrm{1}}\left(x\right)+E_{\mathrm{1}}\left(y\right)<{\displaystyle \frac{\mathrm{1}}{x^{\mathrm{1}/\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{y^{\mathrm{1}/\mathrm{2}}}}.$$

Algebraic manipulation yields the following search interval for t>t_h:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& \mathrm{0}<h<\left(t-t_{\mathrm{h}}\right){\left[{\displaystyle \frac{\mathrm{1}+z}{z}}{\displaystyle \frac{\mathrm{1}}{-E_{\mathrm{1}}\left(x\right)+E_{\mathrm{1}}\left(y\right)}}\right]}^{\mathrm{2}}\\ & z={\left({\displaystyle \frac{t-t_{\mathrm{h}}}{t}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

Golden-section search with the inequality Eq. (A3) can cause numerical underflow when computing f(h) near the right endpoint of the search interval. In lieu of Golden-section search, the inverse for the cooling section is computed using Nelder–Mead optimization with the inequality Eq. (A3) used as a Box constraint (Box, 1965).

Appendix B

Acronyms
ADC	Analog-to-digital converter
AFE	Analog front end
COTS	Commercial off-the-shelf
DAC	Digital-to-analog converter
DC	Direct current
DP	Dual probe
HPP	Heat pulse probe
MB	Mean bias
M2M	Machine-to-machine
NA	Not applicable
OAT	One-at-a-time
PCB	Printed circuit board
PD	Percentage difference
PID	Proportional–integral–derivative controller
PLL	Phase-locked loop
PVC	Polyvinyl chloride
RMSD	Root mean squared difference
RS-232	Recommended standard 232 serial port
SDI-12	Serial digital interface at 1200 baud
SDRAM	Synchronous dynamic random-access memory
SP	Single probe
USB	Universal serial bus

Symbols and SI units
ai	log (r(ti))
α	Thermal diffusivity of soil (m2 s−1)
$\mathit{\{}B,C,D\mathit{\}}$	Coefficients used for the SP model
Ch	Volumetric heat capacity ($\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$)
$\left\{C_{\mathrm{m}},C_{\mathrm{o}},C_{\mathrm{w}}\right\}$	Volumetric heat capacity (mineral, organic, and water)
c	Specific heat capacity ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
{cm,co}	Specific heat capacity (mineral, organic) ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
ΔΓ(t)	Change in temperature (K) as a function of time for SP
γ(r,t) or γ(r(t),t)	Signal processing computation step
d∕dt	Time derivative
ΔE	Voltage drop over resistive element (V)
Ei(x)	Exponential integral function
Ekn	Known (measured) output voltage (V)
f(x,y)	Function of x,y
fs	Sampling rate (Hz)
H	r2∕4α
$\left\{{\mathit{\theta }}_{\mathrm{m}},{\mathit{\theta }}_{\mathrm{o}},{\mathit{\theta }}_{\mathrm{w}}\right\}$	Volume fractions (mineral, organic, and water)
I	Current through nichrome heater wire (A)
i	Integer index
k	Thermal conductivity ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
ℓ	Length of heater needle (m)
log ( )	Natural logarithm function
N	Index number as integer
π	≈3.14159
P	Electrical power (W)
ρ	Density (kg m−3)
{ρm,ρo}	Densities (mineral, organic) (kg m−3)
q	Energy transfer per time per length of heater needle (W m−1)
Rs	Resistance of sensor resistor (Ω)
Rw	Resistance of nichrome heater wire (Ω)
r(t)	Effective DP radius as a function of time (m)
r	DP radius (m)
rinitial	Initial DP radius (m)
rn	Radius of heater needle (m)
t	Time (s)
t0	Start time of heating nichrome wire (s)
ta	Additional time delay (s)
tc	Time of cooling (s)
td	Time at which curve is assumed to be linear (s)
th	Stop time of heating nichrome wire (s)
ti	Time used to determine initial temperature (s)
tp	Time at which curve is at a maximum peak (s)
$t_{\mathrm{T}}=t_{\mathrm{i}}+t_{\mathrm{h}}+t_{\mathrm{c}}$	Total time of experiment (s)
ΔT(t)	Change in temperature (K) as a function of time for DP
$\mathrm{\Delta }t=\mathrm{1}/f_{\mathrm{s}}$	As the time step (s)
$x,y,z$	Real numbers
||	Absolute value function