Remote Sensing of Soil Moisture in Vineyards Using Airborne and Ground-Based Thermal Inertia Data

Abstract

Thermal remote sensing of soil moisture in vineyards is a challenge. The grass-covered soil, in addition to a standing grape canopy, create complex patterns of heating and cooling and increase the surface temperature variability between vine rows. In this study, we evaluate the strength of relationships between soil moisture, mechanical resistance and thermal inertia calculated from the drop of surface temperature during a clear sky night over a vineyard in the Niagara region. We utilized data from two sensors, an airborne thermal camera (height ≈ 500 m a.g.l.) and a handheld thermal gun (height ≈ 1 m a.g.l.), to explore the effects of different field of views and the high inter-row temperature variability. Spatial patterns of soil moisture correlated more with estimated thermal inertia than with surface temperature recorded at sunrise or sunset. Despite the coarse resolution of airborne thermal inertia images, it performed better than estimates from the handheld thermal gun. Between-row variation was further analyzed using a linear mixed-effects model. Despite the limited spatial variability of soil properties within a single vineyard, the magnitudes of the model coefficients for soil moisture and mechanical resistance are encouraging indicators of the utility of thermal inertia in vineyard management.

1. Introduction

The quality and quantity of grapevine production is controlled by many factors, such as soil characteristics, climate, management system and the frequency of exposure to pests and diseases. Recent studies [1,2] show that productivity within a single vineyard could vary as much as eight-fold. Precision viticulture takes advantage of remote sensing and geomatics to model this variation and estimate yield quality and quantity at the vineyard level [3].

Soil particularly is an important factor in determining the productivity of vineyards. Observations show that high and low production regions within a vineyard tend to be stable over a longer time [4], and these patterns relate to soil spatial distribution, micro-climate patterns and topography variations [5]. Identifying zones with similar soil type helps in the planning of a vineyard, by selecting the suitable grape varieties to soil type and allocating vineyards with homogenous soil to allow easy management [6]. In addition, soil information explains the interplay between year-to-year rainfall and production. Therefore, “considerable effort in precision viticulture research aims at measuring and mapping spatial variability in soils at the single vineyard scale” [7].

Remote sensing provides high quality spatial data for vineyard management. However, it is not applied widely in viticulture [7]. Optical remote sensing is used to sense changes in properties of the few millimeters of the soil surface [8]. Alternatively, researchers apply non-contact electromagnetic survey to map soil variability within a vineyard [4,9]. Measured apparent electric conductivity is used as a proxy for soil moisture content, soil texture and salinity of the soil solution [2,10].

The utility of thermal remote sensing in detecting energy and moisture fluxes at the land surface is well documented [11–15]. For the purpose of monitoring soil moisture content, the common scheme is to decouple the surface thermal properties from ambient temperature (daily temperature cycle) by calculating the thermal inertia (TI), which is a physical property that characterizes the surface resistance to ambient temperature change [16–19]. Various studies report a strong relation between soil moisture content and TI [17,20,21]. However, the thermal inertia method is mostly conducted over bare and dry ground, to avoid complexity added by variations in evapotranspiration patterns [22]. Nevertheless, recent studies [23,24] showed that soil moisture could be estimated over partially vegetated soil if a linear relation between ground flux and surface temperature is maintained.

Verhoef [17] calculated TI using the surface temperature drop, during nights with clear sky and still conditions, to avoid the complex surface energy exchange that occurs during the day. The author found a significant relation between TI calculated over bare soil and volumetric soil moisture content. However, remote thermal inertia techniques were not applied to vineyards. The previous method [17] has a potential in vineyard application, because it avoids the complex heating and evapotranspiration during the day time. However, a careful test of the method is needed to establish the validity of this method over vegetated surfaces [25,26].

In this study, we evaluate a technique for estimating thermal inertia using airborne thermal images acquired over a grass covered soil in a vineyard in the Niagara Region, Ontario, Canada. The technique is based on the drop of surface temperature during the night and has not been tested over grass covered soil. We further explore the functional relationships between estimated thermal inertia in the presence of grass sod (we will refer to it subsequently as TI_c) and subsurface soil properties (moisture and mechanical resistance). Finally, we provide suggestions for improving soil moisture retrieval using the nocturnal thermal inertia method.

2. Background

2.1. Thermal Inertia: Theoretical Background

TI [J·m⁻²·K⁻¹·s^−1/2] of a bare soil is a physical property that describes the response of soil to an ambient temperature change:

$$TI=\sqrt{\rho ck}$$

(1)

where ρ is the soil density [kg·m⁻³], c is soil specific heat capacity [J·kg⁻¹·K⁻¹] and k is soil thermal conductivity [W·m⁻¹·K⁻¹]. TI can be calculated from the night cooling of land surface assuming a constant rate of surface cooling [17,27,28]:

$$TI=\frac{2\left|\overline{R_n}\right|\sqrt{\mathrm{\Delta }t}}{\mathrm{\Delta }T\sqrt{\pi }}$$

(2)

where |̅R_n| [Wm⁻²] is the average net radiation during the night, ΔT [K] is the night temperature drop and Δt [s] is the cooling period in seconds. The common method for calculating thermal inertia depends on the periodic daily heating [29]; in contrast, Equation (2) depends on the non-periodic cooling of the surface under still and clear sky conditions. Theoretically, if one estimated thermal inertia over the same area using both methods, the results should be similar. However, the absence of turbulent heat fluxes (i.e., sensible heat flux and latent heat) during the night simplify the relation between surface temperature and ground heat flux, which cannot be guaranteed during the day [30].

Murray and Verhoef [25] proposed that increasing soil saturation will result in a logistic increase of TI. The authors based their theoretical relation on a model of thermal conductivity as a function of soil saturation by Johansen [31]:

$$TI=Ke\left({TI}_s-{TI}_d\right)+{TI}_d$$

(3)

where the subscripts, s and d, denote the saturated and air-dry conditions, respectively, and Ke is a modified Kersten number, given by:

$$Ke=exp\left\{\gamma \left[1-{\left[\frac{\theta }{{\theta }_s}\right]}^{\gamma -\delta }\right]\right\}$$

(4)

where γ is a soil texture-dependent parameter, δ is a shape parameter and θ/θ_s [-] is the soil saturation ratio. Estimating soil moisture content can be done by inverting Equation (3) with the Kersten number, approximated by [21]:

$$Ke=\frac{TI-{TI}_d}{{TI}_s-{TI}_d}$$

(5)

2.2. Thermal Inertia Modification by a Vegetation Cover

The brightness temperature, measured using a thermal infrared sensor over a grass-covered soil, is modeled as the summation of (a) the energy of the soil surface emission, which passes through the plant canopy, (b) the energy of the plant canopy emission and (c) the reflected energy of plant canopy emission by the soil surface below it, which passes through the canopy [32]:

$$T_b=\varepsilon .\zeta .T_s+\left(1-\omega \right).\left(1-\zeta \right).T_c+\left(1-\varepsilon \right).\left(1-\omega \right).\left(1-\zeta \right).T_c$$

(6)

where Tb [k] is the brightness temperature measured by the thermal infrared (TIR) sensor, T_s [K] is soil surface temperature, T_c [K] is the plant canopy temperature, ε [-] is soil surface emissivity, ω [-] is the single scattering albedo and ζ [-] is the transmissivity of the vegetation canopy. The grass canopy (leaves) temperature differs from ambient air temperature by the net radiation at both the surface of the leaf and by the temperature diffusive resistance, which is a function of leaf size and wind speed [33]. The amount of heat storage, due to photosynthesis, is negligible over a day period.

If remote sensing measurements are taken on a still clear night over a grass-covered soil, it can be assumed that the grass temperature is coupled to the ambient temperature. This will result in a linear reduction of the soil surface temperature, as determined by the transmissivity of the grass canopy and the sensor viewing angle (Equation (6)). Therefore, we postulate that using Equation (2) and surface temperature measured over a grass covered soil will result in an estimated TI_c, which is proportional to the true TI of the soil below the grass. Although Equation (2) has not been applied to a grass-covered surface before, a previous field study by Kim and England [34] reports a significant relation between TI calculated using passive microwave and soil moisture content over a grass covered area.

3. Methods

3.1. Study Area

An aerial survey was conducted over a block of grapes (Vitis vinifera var. Merlot) in a commercial vineyard near Niagara-on-the-Lake, Ontario, Canada. The block was located near the shoreline of Lake Ontario, between latitudes 43°14′50.18″N and 43°14′43.34″N and longitudes 79°8′42.54″W and 79°8′35.07″ W. The vineyard was established on the Winona soil series, which is characterized by fine sandy loam texture originating from deltaic sand deposits over clay [35]. The spacing between vineyard rows was 2.44 m, and the rows were oriented in a north-south direction. The vineyard floor (inter-row) between rows was vegetated with a mixed grass sod, and the vines within the rows were spaced at 1.22 m intervals.

3.2. Fieldwork and Data Collection

Twenty-five plots (4 m² each) were selected, based on the variations seen in the grape canopy during the previous growing season (Figure 1). The inter-row grass cover was cut to less than 1 cm above the ground, using a commercial grass trimmer, to measure soil moisture using the theta probe and soil electrical conductivity using an electromagnetic probe (EM-38); refer to

Table 1. Field measurements.

Measurement	Replicates	Proxy for
Bulk electrical conductivity (by EM38)	Horizontal mode perpendicular to vine rows Horizontal mode parallel to vine rows Vertical mode perpendicular to vine rows Vertical mode parallel to vine rows	Clay mineralogy and soil solution salinity given dry conditions
Mechanical resistance (by soil penetrometer)	Right slice of the vine row Middle of vine row Left slice of the vine row	Soil bulk density, soil compaction and stoniness
Soil dielectric constant (by Theta Probe)	Same as mechanical resistance sampling	Soil volumetric moisture content
TIR emissions (by thermal infrared gun)	Same as mechanical resistance sampling	Surface temperature

for a summary of the field sampling scheme. The average apparent electrical conductivity at the center of each trimmed plot was measured using an EM-38 probe (Geonics Ltd.). Four EM-38 measurements were taken at the vertical and horizontal orientations of the probe coils, as well as parallel and perpendicular to the vineyard rows around sunset. The measurements were corrected to 25 °C [36]. Volumetric soil moisture content of the upper 15 cm was measured using a handheld Theta probe ML2x (Delta-T Devices, Ltd.). The measurements were taken after the sunset of 17 September 2007, with three replicates for each plot, as indicated in Figure 1. Soil mechanical resistance was measured using a Rimik CP20 cone penetrometer (Rimik Agricultural Electronics, Ltd.). The penetrometer measurements were taken after acquiring the sunrise aerial images following the same three replicates per plot scheme. Soil surface temperature was measured using a Raytek Raynger ST2 thermal infrared gun (Fluke, Inc.) with ±1% accuracy. The IR instrument was pointed with a 45° angle at the trimmed surface at a distance of one meter, which results in a ground footprint of approximately 0.25 m². Three replicates were measured at the far left, the right side and middle portion of the each plot (Figure 1). Soil emissivity was set at a value of 1 (perfect blackbody); subsequently, an emissivity correction was applied; refer to Section 3.3. Surface radiation balance were measured using a CNR1 four-component net-radiometer (Kipp & Zonen Ltd.), which was connected to a 21× data logger (Campbell Scientific Ltd.).

A plastic mast was fixed between row seven and eight from the eastern edge of the plot, to mount the net-radiometer (Figure 1). The surface grass around the mast was trimmed to obtain similar conditions to sampling plots. Hourly measurements of air temperature and relative humidity were obtained from the National Climate Data Center (NCDC) archive for the Niagara district weather station (Figure 2, WMO ID USAF:712625).

3.3. Aerial Mapping of Thermal Inertia

Thermal images were recorded using a thermal camera FLIR model ThermaCAM SC2000 (FLIR Ltd.). The thermal camera was pointed by one of the aircraft crew to the nadir view, and images were recorded through an opening in a modified aircraft door. The image recording was set to the automatic collection mode with a 1-s interval. Two flight passes were conducted over the same merlot plot around sunset and sunrise, to measure the drop in temperature during the clear sky night of 17 September 2007 (Figure 2). Organizing a flight around sunset was less limited than the dawn flight, because of Canadian regulations restricting night flights of light aircrafts. Therefore, the sunrise flight was delayed, for half an hour, to benefit from the increased visibility during the civil twilight period. We assumed that this delay did not influence our results, because sun rays at extremely low angles were restricted by the vertical vine canopy ([37,38]; refer to

Table 2. Flight schedules and weather conditions.

Flight Time	Estimated Surface Emissivity	Sky Temperature [°C]	Relative Humidity [%]	Air Temperature [°C]	Flight Height [m AGL]
6:50 p.m.–7:03 p.m. 17 September 2007	0.94	−1.4	40.5	20.6	≈ 490
7:11 a.m.–7:33 a.m. September 18, 2007	0.94	−3.4	78.8	13.1	≈ 490

for details of the flights).

3.4. Calculating Thermal Inertia

Thermal infrared gun measurements were corrected for land surface emissivity and reflection of sky temperature. The initial emissivity value for TIR gun measurements was one. After, we corrected the measurements by assigning an emissivity value of 0.94 [39] to represent a mix of soil and vegetation tissues. Sky temperature was derived from the incoming long-wave radiation (5 to 50 μm) measured by the net-radiometer, which had a wider wavelength window than the TIR sensors (8 to 14 μm). We assumed that the incoming long wave radiation was homogeneous over the entire field. Atmospheric correction was not considered for the TIR gun, because of the lack of climatic variables at the site. Using variables from another station could generate a significant bias, because of the different surface cover type. An atmospheric correction was applied to the thermal infrared images using the built-in radio-transfer model in the FLIR analysis software, Quick Report v2.1 (FLIR, Ltd.); details of the atmospheric model were given in [40]. The relative humidity and temperature of an air column with a height of 493 m a.g.l. (flight height) were approximated from the weather station data (Table 2). The approximation was based on the relative stability of the air column during still nights in contrary to a thick daytime mixing layer. However, one should notice that these conditions could create an inversion, which will offset the weather station measurements (colder) from the atmospheric profile average.

A value of 0.94 was assigned to surface emissivity similar to that used to correct the TIR measurements. The vineyard plot was captured by two thermal images at both flight times (sunset and sunrise). A linear regression model was used to match the temperature of one scene to the other. The model was based on 100 randomly sampled clusters with 5 replicates (total of 500) in the overlap region between the two scenes. The two images were then merged into a single mosaic. Subsequently, geographic coordinates were assigned to the mosaic image using the coordinates of the plot corners.

The locations of the 25 plots were identified on the sunset and sunrise thermal image mosaics using the center geographic coordinates of each plot. The average temperatures of matched pixels to field plots locations and the four connected neighbors for each pixel were recorded (i.e., upper, lower, right and left neighbor pixels). We decided to consider the average of five pixels (≈2 m² on ground) instead of using the center value, only because of the susceptibility of single pixels (0.6 × 0.6 m) to errors in image registration and plot coordinates. In addition, averaging land surface temperature over a number of pixels limited the variability of surface temperature caused by local changes in the structure of grass canopy. In a few cases, the center of the plot was found to be located on top of the grape canopy, as indicated by the low temperature associated with the dense plant leaves. In these cases, the sample was relocated to the nearest between-row pixel.

TI_c was calculated from both the TIR gun measurements (field TI_c) and homogenized aerial thermal images (remotely-sensed TI_c) using Equation (2) at each of the 25 plot locations. ΔT was calculated as the difference between land surface temperature (LST) at sunset and sunrise, while Δt was calculated as the time difference (in s) between sunset and sunrise. The net radiation during night (constant cooling in Wm⁻²) was calculated as the average of net radiation at sunset and sunrise. The net radiation at sunrise or sunset was given by an arithmetic sum of incoming and outgoing shortwave and incoming and outgoing longwave radiation. In order to account for local variation in net radiation, radiation balance components, except outgoing long wave, were estimated from the net radiometer measurements. The outgoing longwave radiation was estimated locally from the measured LST using the thermal gun or thermal images.

3.5. Statistical Analysis

The aim of the statistical analyses was to explore empirical relationships between TI_c variables as our response and soil properties (i.e., moisture content, mechanical resistance and electrical conductivity) as predictors. Theta Probe measurements were used as a proxy for soil moisture content, while the horizontal parallel mode of EM38 measurements were interpreted as proxies of soil salinity and clay content, given the dry conditions. We selected the horizontal mode of the EM38 measurements to increase the sensitivity to the uppermost layer of the soil and limited our analysis to the parallel mode to minimize interference with the vines’ metal wires [41]. The averages of vertical profiles of soil mechanical resistance were calculated at each field plot location. In addition, a binary variable was calculated to indicate the existence of missing mechanical resistance measurements at each profile. The lack of such measurements at a location was interpreted as the presence of consolidated substrate or the abundance of rock fragments.

Pearson’s correlation coefficients were calculated between all available quantitative variables, to assess the strength of linear associations between soil characteristics and TI_c. In the case of the mechanical resistance, the area under the receiver operating characteristics curve (AUROC) was used instead, which was a value between 0.5 (no separation) and 1 (perfect separation), which measured a quantitative variable’s ability to discriminate the two classes represented by a binary variable [42]. In addition, we tested the statistical significance of correlations between TI_c calculated from the thermal gun and from aerial imagery at the α = 0.10 level of significance. This elevated significance level was chosen to reduce the Type II error rate in this exploratory analysis with a small sample size [43].

Multiple linear regression models were used to analyze the relationship between the response and Theta Probe, EM38, mechanical resistance (average) and mechanical-resistance (binary) as TI_c candidate predictor variables. We developed separate models for TI_c estimated over grass canopy from aerial imagery and for TI_c using thermal infrared gun over trimmed field plots. We started with all the explanatory variables and, then, performed a backward elimination based on the Akaike Information Criterion (AIC). The AIC was based on a model’s log-likelihood. It penalized larger models in order to give preference to smaller models that fit the data almost as well [44].

A linear mixed-effects model was built to reproduce the regression analysis of field measurements at the level of the individual observations, which were grouped at the plot level. Compared to ordinary linear regression models, linear mixed-effects models account for grouping (e.g., repeated measurements or spatial clustering) by introducing a so-called random-effect term that captures group-specific random variation (e.g., [44]). TI_c (field) calculated from TIR gun measurements was selected as the response variable, and Theta Probe and binary mechanical resistance were chosen as the fixed effects, while the grouping based on the replicates’ locations between vine rows (east, middle, west) was included as a random-effect term. Fixed-effects variables were selected manually based on the AIC. All statistical analyses were performed using the statistical software, R [45], and its package ‘nlme’ for the mixed-effects model [46].

4. Results

4.1. Aerial Thermal Imagery

In the radiometrically corrected and mosaicked thermal images root-mean-square-error (RMSE of georeferencing: 0.9 pixel, 0.4 m), the difference in surface temperature patterns between low temperature grape rows and grass covered ground is clearly distinct as vertical strips (Figure 3). In addition, thermal anomalies, related to soil and vegetation, are observed across the two scenes. For example, an anomaly is visible in the lower half of both scenes; this anomaly is characterized by being warmer (colder) during sunset (sunrise), which indicates a region with lower TI_c.

Remarkably, there is a tonal difference at the edges of the individual images, particularly at sunrise. The residuals of a linear regression model that relates the two images in each mosaic (Figure 4) are small in absolute terms (standard deviation of 0.137 K at sunset, 0.142 K at sunrise), but relatively large compared to the small variation in surface temperatures. The difference between the two images is partly attributed to geo-referencing uncertainties and to the difference in LST, based on the viewing angle. Nevertheless, the two mosaics capture the drop of surface temperature compared to the in situ net radiometer measurements (Figure 5).

4.2. Soil Physical Properties

Remote sensing and field measurements are summarized in

Table 3. Descriptive summary of remotely-sensed and in situ measurements (TI, thermal inertia.)

Variable	Minimum	Lower Quartile	Median	Mean	Upper Quartile	Maximum
Apparent Electrical Conductivity Horizontal-Parallel Mode [mS/m]	18.6	30.7	37.2	40.7	47.8	68.6
Volumetric soil moisture (Theta Probe) [%]	9.7	12.7	15	14.5	16.3	19.4
Mechanical Resistance (average) [kPa]	1,480	2,056	2,368	2,399	2,655	3,306
LST-sunset (field) [°C]	14.2	14.5	14.6	14.6	14.8	15
LST-sunrise (field) [°C]	9.6	10.1	10.2	10.2	10.4	10.7
Lst-sunset (remote) [°C]	15.2	15.5	15.5	15.5	15.6	15.9
Lst-sunrise (remote) [°C]	10.9	11.1	11.2	11.2	11.3	11.5
TIc (field) [J·m−2·K−1·s−1/2]	2,894	2,894	3,357	3,410	3,609	3,872
TIc (remote) [J·m−2·K−1·s−1/2]	3,068	3,287	3,348	3,361	3,430	3,581

. We applied a linear correction term to the TIR gun measurements at sunset to account for a cooling drift observed during data collection. The correction term was estimated by repeating the measurement of the first plot at the end of the data collection session and assuming a linear cooling of the soil during the measurement period. There was no observed trend in TIR gun measurements at sunrise, and therefore, no correction was applied. Field measurements shows a narrow range of surface temperatures at sunrise and sunset with an approximate 4 °C thermal offset between the two.

The average remotely-sensed TI_c is 3,361 J·m⁻²·K⁻¹·s^−1/2 and 3,410 J·m⁻²K⁻¹·s^−1/2 for field TI_c, which is relatively high compared to bare soil values (700–1,200 J·m⁻²·K⁻¹·s^−1/2) [17]. Theta Probe measurements of soil moisture indicate a dry soil with volumetric moisture content ranging from 10–19%. The mechanical penetration resistance is highly variable, likely due to the abundance of subsurface rock fragments and compacted soil layers in the studied field.

4.3. Relationships between Field Measurements and Remote Sensing Imagery

We will first explore the correlations between field/remotely-sensed surface temperature and thermal inertia before we present empirical models of TI_c, as explained by measured soil properties. In this context, we expect that the effect of grass cover is more pronounced in the estimates of remotely-sensed LST and TI_c compared to field LST and TI_c. Consequently, remotely-sensed variables are based on an average of five pixels between the rows (at each plot location), and the grass covered pixels are dominant. Field and remotely-sensed TI_c have a moderate positive correlation of 0.37 (p-value of 0.065;

Table 4. Pearson’s correlation between field and remote measurements of land surface temperature and thermal inertia (TI_c).

Variables	Correlation Coefficient (p-value)
LST (remote), LST (field) at local sunrise	0.075 (p = 0.721)
LST (remote) vs. LST (field) at local sunset	0.392 (p = 0.053)
TIc (remote) vs. TIc (field)	0.374 (p = 0.065)

), similar to the correlation between field and remotely-sensed LST at sunset (0.39, p-value of 0.053). At sunrise, there is no correlation of remote and field LST (0.08, p-value of 0.721). The low correlation indicates non-proportional variations between remotely-sensed and field measurements over the 25 locations, which is likely due to collecting the measurements over a longer time (not in an instantaneous fashion), the limited field of view for ground measurements and the nonstandard view angle of the non-contact thermometer.

The correlation analyses show a weak to near-zero correlation of TI_c and surface temperature variables with the soil-related predictor variables (

Table 5. Pearson’s correlations of thermal inertia and surface temperatures with the predictor variables. The area under the receiver operating characteristics curve (AUROC) is reported for the binary mechanical resistance.

	EM38 Horizontal Parallel	Theta Probe	Mechanical Resistance Average	Mechanical Resistance (0/1)
TIc (field)	−0.057	0.392	0.108	0.699
TIc (remote)	−0.351	0.500	0.167	0.706
LST-sunset (field)	0.199	−0.274	−0.326	0.688
LST-sunrise (field)	0.085	0.316	−0.066	0.614
Lst-sunset (remote)	0.128	−0.331	−0.195	0.772
Lst-sunrise (remote)	−0.267	0.249	−0.001	0.522

). Correlation is strongest between soil moisture and remotely-sensed TI_c, as well as field TI_c (0.50 and 0.39, respectively). Soil moisture content also correlates positively (negatively) and non-significantly with surface temperature at sunrise (sunset). These relations clearly indicate the sensitivity of TI_c to moisture content, even under a grass canopy, compared to LST at sunset or sunrise. There is also some evidence of a negative correlation of TI_c with EM38 and possibly a positive, but weak and non-significant, correlation with mechanical resistance.

4.4. Regression Relationships between Thermal Inertia and Soil Properties

All three linear models relating TI_c measurements to soil properties show remarkably similar structure and comparable coefficient estimates (

Table 6. Summary of statistical models, including the coefficient estimates and their standard errors in square brackets; in parentheses, p-values for two-sided tests. Average mechanical resistance is not included in any of the models.

Response Variable	Intercept	Theta Probe	EM38	Mechanical Resistance (0/1)	R2 (adjusted R2)	Residual Standard Deviation
TIc (field)	2,863.72	34.57	-	138.17	0.224	227.2
		[17.59]		[97.73]	(0.153)
		(0.062)		(0.171)
TIc (remote)	3,118.101	21.365	−2.313	81.965	0.416	105.6
		[8.265]*	[1.589]	[45.721]	(0.333)
		(0.01))	(0.160)	(0.087)
TIc (field) In the mixed-effects model	2,950.321	31.0470	-	139.143	0.167*	287.9***
		[11.385]		[86.944]	0.624**
		(0.009)		(0.116)

^*Pseudo-R² without random effects;

^**Pseudo-R² with random effects;

^***without random effects.

). Locations with five percentage point higher soil moisture have an about 100 unit higher remotely-sensed TI_c or a 150 unit higher field TI_c (all p-values < 0.10) when accounting for the effects of the other variables included in the models. The mechanical resistance indicator variable is also included in all models; locations with the presence of consolidated or rock substrate are associated with an estimated ∼80 unit higher remotely-sensed TI_c or ∼140 unit higher field TI_c, subject to substantial uncertainty (p-value < 0.10, only for the former). EM38 measurements are only selected into the model for remotely-sensed TI_c, suggesting that locations with 20 m·S/m higher apparent electrical conductivity had an estimated ∼40 unit lower TI_c, although this effect is not significant.

The models account for less than half of the variation in the response variables (R² between 17% and 42% of total TI_c variation). The largest coefficient of determination is achieved for remotely-sensed TI_c, which can be attributed to the variance reduction by spatial aggregation (Table 3). This variance reduction may also explain the consistently smaller coefficient estimates in the model for remotely-sensed TI_c. The results of the mixed-effects models (where the inter-row sampling position is selected as a random effect) suggest that there is a significant inter-row variability, as indicated by the pseudo R² values.

5. Discussion

5.1. Thermal Inertia and Soil Properties

TI_c values were found to be higher than the bare soil TI range. In Equation (2), small changes in surface temperature or net radiation resulted in substantial changes in calculated TI values. It is suggested that the observed TI_c values were due to reduction in ΔT, which are more susceptible to error than net radiometer measurements. Lamb et al. [2] had shown that unless the pixel size was much smaller than the inter-row distance pixel, digital number values would always be a mix of grape canopy and the inter-row surface. In our study, the inter-row (2.4 m) was captured in four pixels, thanks to the sub-meter resolution (0.6 m). Hall et al. [7] indicated three categories of pixels, when pixel length was smaller than the inter-row distance: (a) vine only; (b) non-vine and (c) mixed pixels. Given that we selected the corresponding pixels to the sampling plots and their four neighbors, which translates to three pixels-wide (center pixel and the two neighbors), it is likely that the averaged remotely-sensed LST contained both non-vines and mixed pixels. The mixed pixels were influenced by sticking out branches of the vine rows. Unlike the grass leaves, grape canopy temperature was significantly decoupled from the soil surface temperature and could cause further reduction of estimated ΔT. Adding to that, the grass sod also modified the surface temperature, increased its buffer capacity and resisted change in surface temperature.

Nevertheless, soil properties explained a substantial proportion of the variation of thermal inertia within a vineyard, both in the remotely-sensed and the field analysis (Figure 6). While the limited overall spatial variability of surface and substrate properties within a single vineyard had to be taken into consideration when interpreting coefficients of determination and hypothesis tests in this context, the estimated effect sizes for soil moisture content and mechanical resistance were encouraging indicators of the influence of these two variables on thermal inertia. Apparent electrical conductivity possibly showed too little variation in this study to provide results of practical relevance.

A significant linear relationship was found between TI_c (field and remotely-sensed) and averaged soil moisture content and mechanical resistance. The narrow range of soil moisture content allowed us to approximate the relation using a linear relation, even though the theoretical model (Equations (3) and (4)) was characterized by a logarithmic increase over a wider soil moisture range [25]. Furthermore, the theoretical model was sensitive to soil texture (variable γ in Equation (4)). The sandy texture of the vineyard had increased the sensitivity of TI_c at lower moisture content, which was consistent with the steep thermal inertia response to low soil moisture increase for coarse textured soils; Figure 1b in [25].

5.2. Advantages of Remotely-Sensed Thermal Inertia

While remotely-sensed data evidently had a reduced spatial resolution compared to field measurements, this spatial aggregation also reduced the effects of small-scale variation on statistical results, while still providing data at a level of detail that was relevant to viticulture practice. Remotely-sensed TI_c was more indicative of a spatial average of soil moisture and less susceptible to local variation, as indicated by the high statistical significance of Theta Probe regression coefficients compared to the case of field TI_c. Similar observation of enhancing the soil moisture estimation by using aggregated thermal remote sensing data was noticed in a recent study by Minacapilli et al. [20].

Although the grass cover had been trimmed at each sampling spot, the average Theta Probe measurements at the trimmed plots were representative of the spatial average of soil moisture content over a larger radius than the trimmed plot area, because removing the grass was done on the day of aerial imaging. The short-term exposure was unlikely to influence the soil moisture content of the subsurface (≈15 cm) under a dry condition. The observation that remote measurements (500 m) were more correlated to soil moisture measurements than on ground electromagnetic probe (r values) was remarkable. A positive relation was also found with the binary indicator of soil mechanical resistance, although it was not significant.

The vegetative canopy was found to have strong effects on remotely-sensed land surface temperature, as found by Kim and England [34]. For example, there was a noticeable difference between the two scenes in each image mosaic, as indicated by the random sample from the overlap region (Figure 4). This systematic difference in LST was likely induced by coregistration errors and LST variability based on the structure of grass canopy and the instantaneous view angle (Equation (6)). Adding to that, our estimated values of field TI_c were approximately double the thermal inertia values reported for bare soil surface. This increase was attributed to a partial plant cover, which is characterized by a high TI_c compared to the TI_c of a bare soil surface.

5.3. Limitations of Field Estimated Thermal Inertia

The use of TIR guns to measure surface temperature and estimate thermal inertia was found to impose several practical challenges. For example, there was a time lag between measuring the surface temperature from one plot to another. This time lag resulted in further cooling of the plots at the end of the measuring route. Adding to this, it was difficult to maintain a constant view angle of the unit under manual operation.

Remote and field estimates of TI_c and LST were poorly correlated, especially LST around sunrise. We speculated that the accuracy and representativeness of field measurements given a limited field of view and a high inter-row variability might have caused this discrepancy. There was no clear explanation as to why low correlations were obtained around sunrise, but the isothermal conditions and lack of temperature contrast made it susceptible to noise.

Field TI_c values were found to be higher than the bare soil TI range and in the same range of remotely-sensed TI_c. Although the grass was trimmed a few hours prior to taking the observations, the root mass and remaining tissues still constituted a large portion of the exposed soil surface, which could account for the elevated range of field TI_c. Adding to that, it was difficult to calculate the exact time difference (Δt) between sunset and sunrise. We used the net-radiometer incoming shortwave radiation as a guide. However, one should consider that vineyards were designed intentionally to increase the sunlight on grapes’ canopy rather than on the vineyard floor. Previous studies [34,35] showed a delay of the heating cycle at the ground surface for a few hours.

5.4. Improving Moisture Retrieval Using Thermal Inertia

The empirical results indicated that TI_c had a potential for monitoring variation in soil moisture content within a uniformly compacted field (i.e., equal mechanical resistance). Moreover, we showed the possibility to extend TI_c application to moderately vegetated vineyard soils. Despite the challenges encountered in this study (i.e., limited range of soil moisture content, limited sample size, difficulties in scheduling flights), our findings supported the promise of the TI_c technique.

Several issues need to be considered for future development of remotely-sensed TI_c from nocturnal cooling. For example, the complete characterization of atmospheric temperature and water vapor content is a challenge. Measuring atmospheric variables simultaneously at the land surface and on the aircraft could provide data to interpolate profiles of air temperature and water vapor and estimate representative values for atmospheric correction. In our study, we have used incoming long wave radiation to estimate reflected sky temperature. This method could yield a positive bias, because of the different sensitivity range of each instrument. We suggest measuring the sky temperature using the TIR thermometer directed at a low emissivity surface (i.e., aluminum surface). Moreover, plant canopy is a major source of uncertainty. Therefore, including remote sensing products, such as the Normalized Difference Vegetation Index (NDVI) [18,47], full waveform analysis of LiDAR data [48] and leaf area index [49] to characterize the canopy density and leaf area index will succeed in reducing uncertainty. In addition, NDVI can be used to estimate a spatially distributed surface emissivity map, rather than using a constant emissivity value.

Soil moisture content contributed less to the variation in field TI_c, compared to remotely-sensed TI_c. The results of the mixed effect models suggested that field measurements were susceptible to inter-row variability. Therefore, a spatially distributed infrared sensor network could be used, to obtain a temporally consistent measurement of surface temperature, at a specific time. Adding to that, surface temperature had a significant inter-row variability, as indicated by the mixed-effects model (Figure 6). Therefore, for vineyard applications, the height of the sensor should be increased in order to expand the field of view and to avoid between row variability.

Furthermore, the presented technique could be used for continuous monitoring of vineyards by conducting regular flights during the growing season. A thermal inertia image, at a certain time in the growing season, could be normalized using the minimum and maximum TI_c of this season, assuming that these points corresponded to the wilting point and field capacity. Further research must be conducted to establish the use of TI_c to approximate the Kersten function (Equation (5)) in order to retrieve soil moisture content.

6. Conclusions

In the present study, we evaluated the relationships between soil moisture content, mechanical resistance and thermal inertia (TI_c) over a grass covered vineyard in the Niagara region, Canada. TI_c was calculated using surface temperature from two sensors; an airborne thermal camera (height ≈ 500 meters AGL) and a handheld thermal gun (height ≈ 1 meter AGL). Applying an expression to calculate thermal inertia from cooling during the night revealed elevated TI_c values as a consequence of the influence of grape canopy (average TI_c of 3,361 and 3,410 J·m⁻²·K⁻¹·s^−1/2 for remotely-sensed and field TI_c, respectively). However, we found a significant relationship between TI_c and soil moisture and mechanical resistance, despite the limited range of soil moisture. For example, soil with five percentage points higher soil moisture has an approximately 100 unit higher remotely-sensed TI^c or a 150 unit higher field TI_c (p-values < 0.10), considering that the mechanical resistance of that field does not change. On the other hand, elevated mechanical resistance (due to the presence of consolidated rock or substrate) was found to have ∼80 unit higher remotely-sensed TI_c or ∼140 unit higher field TI_c with higher uncertainty (p-value < 0.10, only for the former). In general, both soil variables explained less than half of the TI_c variability given the complex heating and cooling patterns associated with vineyards. Furthermore, remotely-sensed TI_c was more sensitive to in situ moisture variability compared to handheld TI_c, due to the strong small-scale variability of surface temperature between vine rows. Our results were encouraging for future research on extending thermal remote sensing of soil moisture in heterogeneous land cover regions.