Predicting Soil Water Content on Rainfed Maize through Aerial Thermal Imaging

Restrictions on soil water supply can dramatically reduce crop yields by affecting the growth and development of plants. For this reason, screening tools that can detect crop water stress early have been long investigated, with canopy temperature (CT) being widely used for this purpose. In this study, we investigated the relationship between canopy temperature retrieved from unmanned aerial vehicles (UAV) based thermal imagery with soil and plant attributes, using a rainfed maize field as the area of study. The flight mission was conducted during the late vegetative stage and at solar noon, when a considerable soil water deficit was detected according to the soil water balance model used. While the images were being taken, soil sampling was conducted to determine the soil water content across the field. The sampling results demonstrated the spatial variability of soil water status, with soil volumetric water content (SVWC) presenting 10.4% of variation and values close to the permanent wilting point (PWP), reflecting CT readings that ranged from 32.8 to 40.6 °C among the sampling locations. Although CT correlated well with many of the physical attributes of soil that are related to water dynamics, the simple linear regression between CT and soil water content variables yielded coefficients of determination (R2) = 0.42, indicating that CT alone might not be sufficient to predict soil water status. Nonetheless, when CT was combined with some soil physical attributes in a multiple linear regression, the prediction capacity was significantly increased, achieving an R2 value = 0.88. This result indicates the potential use of CT along with certain soil physical variables to predict crop water status, making it a useful tool for studies exploring the spatial variability of in-season drought stress.


Introduction
Drought stress is considered the most damaging of all of the abiotic stresses, limiting the growth and development of plants and ultimately their yield [1,2]. Climate change is expected to intensify the incidence of extreme weather events, including drought spells [3], which will demand the optimization of agricultural practices that can maximize water use efficiency to sustain actual crop yields, even under rainfed conditions. Precision agriculture (PA) consists of a site-specific form of agriculture that focuses on optimizing farm inputs by conducting the right management practice at the right place, at the right time, and at the appropriate intensity [4][5][6]. To implement this concept, a proper characterization of within-field spatial variability is essential to understand the potential limiting factors.
Because soil physical attributes often present moderate to high spatial variability [7], the amount of water available for plants tends to vary across the field as function of the soil water holding capacity (SWHC) [8][9][10]. When soil water becomes limiting, the plants respond promptly by regulating the stomatal conductance to reduce transpiration as a water saving strategy, which leads to a lower evaporative cooling while increasing leaf temperature [11][12][13]. For this reason, canopy temperature has long been used as a proxy to assess crop water status [14], with remotely sensed temperature being constantly used AgriEngineering 2021, 3 943 for this purpose. Therefore, many studies have successfully used thermal cameras that are able to retrieve canopy temperature to assess plant water stress in crops such as maize [15], soybeans [16], cotton [17], pinto beans [18], citrus orchards [19], and grapevines [20,21].
Most of these studies reflect technology advances, especially when considering the development of cost-effective miniaturized thermal cameras in the last 20 years [22][23][24], which are lightweight and low power enough to fit onto unmanned aerial vehicle (UAV) platforms [25]. From aerial thermal imagery, crop canopy temperature can be monitored with high temporal and spatial resolution at field scale, representing a rapid response tool to detect and quantify plant water stress [13,26]. Although plant water stress expresses restrictions in soil water storage, most studies focus on investigating the relationship between canopy temperature with physiological water stress indicators, such as leaf water potential and stomatal conductance, aiming to validate prediction models for the onset plant water stress. Hence, very few studies have explored how soil water content influences crop canopy temperature measurements. In addition, thermal data have been widely employed as an irrigation resource management tool, with rare studies assessing its use as a screening tool for drought stress on rainfed crops, particularly rainfed maize under field scale, where no studies testing UAV-based thermal imagery as a plant water stress indicator have been reported as of yet.
Therefore, the objective of this study was to explore the use of UAV thermal imagery for monitoring the spatial variability of maize water stress on a rainfed field. The relationship between remotely sensed canopy temperature and maize water stress was assessed based on the plant and soil attributes sampled across the area of study, focusing on direct measurements of soil water content along with crop production metrics and soil physical attributes.

Study Area
The study was conducted in a 40.1 ha field located in Itararé, São Paulo state, southern Brazil (24 • 01 36.93 S-49 • 25 33.46 W, 640 m of altitude) ( Figure 1). The local climate is defined as a subtropical humid (Cfa) climate according to the Köppen classification [27]. Most rainfall occurs between October and March, with a relatively dry period between April and August. The annual average rainfall is 1532 mm, with mean temperature of 17 • C [28]. The main soil type in the region is a Dystrophic Red-Yellow Oxisol [29], with a sandy clay loam texture. The terrain slope average is 6.28 degrees, which necessitates its cultivation in contours and the use of terraces.

UAV Thermal Imagery
Since the main goal of this study was to assess the ability of thermal imagery in detecting crop water stress, one of the key points was detecting a suitable time during the growing season to collect thermal images. For this reason, we constantly monitored the  The crop evaluated in the study was maize sown on 18 February 2019 as a second crop with the variety MS30A37PW, with rows spaced in 0.45 m and a final plant population of 63,000 plants ha −1 . The area has been under a no-tillage system for more than 30 years and is cultivated under rainfed conditions, with a soybean-maize rotation system during the summer cropping seasons and with wheat (Triticum spp.) or black oat (Avena strigosa Schreb.) during the winter season.

UAV Thermal Imagery
Since the main goal of this study was to assess the ability of thermal imagery in detecting crop water stress, one of the key points was detecting a suitable time during the growing season to collect thermal images. For this reason, we constantly monitored the weather conditions and implemented a soil water balance model to identify periods with water deficiency that would prompt the aerial mission. Other restrictive factors that were imposed were to conduct aerial missions only after canopy closure of the interrow, avoiding soil background effect on thermal infrared (TIR) images, and before senescence, during which time, the water consumption is not significant. Based on these assumptions, the aerial mission was conducted on 2 April 2019, when the maize was at a V12 growth stage and when the canopy was completely closed, and with a soil water deficit of 6.3 mm according to a 10-day serial water balance [30] performed using meteorological data measured by a nearby automated weather station and regional soil parameters ( Figure 2).

UAV Thermal Imagery
Since the main goal of this study was to assess the ability of thermal imagery in detecting crop water stress, one of the key points was detecting a suitable time during the growing season to collect thermal images. For this reason, we constantly monitored the weather conditions and implemented a soil water balance model to identify periods with water deficiency that would prompt the aerial mission. Other restrictive factors that were imposed were to conduct aerial missions only after canopy closure of the interrow, avoiding soil background effect on thermal infrared (TIR) images, and before senescence, during which time, the water consumption is not significant. Based on these assumptions, the aerial mission was conducted on 2 April 2019, when the maize was at a V12 growth stage and when the canopy was completely closed, and with a soil water deficit of 6.3 mm according to a 10-day serial water balance [30] performed using meteorological data measured by a nearby automated weather station and regional soil parameters ( Figure 2).  The flight mission was conducted around two hours past solar noon (between 13:38 and 14:31 local time), which is considered to be the optimal time to evaluate plant water stress [31,32]. There were no clouds during the flight, with a mean air temperature of 30.5 • C, relative humidity of 55%, vapour pressure deficit (VPD) of 1.96 kPa [33], and wind speeds up to 2.5 ms −1 .
A low-cost radiometric calibrated FLIR Lepton 3.5 thermal camera (FLIR Systems, Inc., Wilsonville, OR, USA) was used to obtain the aerial images. The camera is composed of an uncooled VOx microbolometer FPA, with a spectral range of 8-14 µm and a resolution of 160 × 120 pixels. Temperature readings are extracted from 14-bit images with a resolution of 0.05 • C and absolute measurement accuracy of ±5 • C or 5% over a range between 10 • C and 140 • C. The camera was attached to a DJI Phantom 4 quadcopter (SZ DJI Technology Co., Shenzhen, China) and was adjusted to collect close to nadir images, being powered up 15 min before flight to ensure proper stabilization time [25] (Kelly et al., 2019). During the flight planning, the parameters were adjusted to deliver images with a forward and side overlap of ≥75% and ≥60%, respectively, at 100 m above ground level, which provided a spatial resolution of 0.63 cm pixel −1 .
The image processing necessary to generate the thermal orthomosaic was implemented following the methodology described by Acorsi, Martello, and Gimenez (2020) [34] using the Agisoft Photoscan Professional (v. 1.2.6, Agisoft LLC, St. Petersburg, Russia) to build the orthomosaic. For georeferencing purposes, a total of six ground control points (GCPs) were previously located in the field using distinguishable targets, with their coordinates measured following a post-processing kinematic (PPK) method, using a pair of single-frequency global navigation satellite system (GNSS) receivers (GTR-ABT TechGeo, Brazil).
Once the orthomosaic was finished, it was imported into ArcGIS software (v. 10.2.2, ESRI Ltd., Redlands, CA, USA) along with 18 waypoints corresponding to sampling locations distributed across the field (Figure 1b). The locations were randomly defined within a polygon excluding headland and terrace areas, respecting a minimum and maximum distance between samples of 40 and 200 m, respectively. Circular buffers of a 3 m radius were generated using the sampling location coordinates, which were then used to extract average temperature values from the thermal orthomosaic for each point.

Field Data
To investigate the relationship between the remotely sensed canopy temperature (CT) and maize water status, we focused on attributes that could explain the water dynamics in the soil. In this sense, we mainly explored the physical attributes from soil that are linked to its capacity to retain water. On the other hand, we also analyzed the plant attributes that are usually affected by long-term drought stress.

Soil Attributes
On the same day that the aerial mission was carried out, soil samples were collected in the layer 0-0.2 m depth at every pre-defined sampling location. The samples were composed of five sub-samples that were distributed within a three-meter radius. The samples were sealed in plastic bags and were taken to a laboratory for further analysis.
To determine the gravimetric water content (SGWC), the samples were first weighed and then oven dried at 105 • C for 48 h and were re-weighted to measure the water content that was lost. Moreover, granulometry (Sand; Clay and Silt) and organic matter content were determined according to Teixeira et al. (2017) [35].
After the crop was harvested, undisturbed soil samples were taken at each sampling point from trenches that were 0.12 and 0.35 m in depth with two replications, using metal rings with a dimension of 0.05 × 0.05 m (diameter × height) In the laboratory, the samples were saturated and then weighed. Moreover, the samples were transferred and kept on a tension table with a −6 kPa water potential until mass stabilization, and the samples were then re-weighed. Lastly, the samples were oven dried at 105 • C for 48 h and were then weighed once more. Thus, the macroporosity (MaP) was determined by the mass difference between the saturated samples and samples that were submitted to the −6 kPa water potential, whereas the microporosity (MiP) was obtained from the mass difference between the samples that were obtained at a water potential of at −6 kPa and the oven dried samples [35]. In addition, the bulk density (BD) was also determined, dividing the dry mass of the soil by the volume of the metal ring. The BD information was then used to calculate the soil volumetric water content (SVWC), multiplying the gravimetric water content by BD [36].
Along with the undisturbed soil sampling, soil resistance to penetration (SRP) was also measured in the field within a depth of 0-0.30 m using a handheld digital penetrometer (PLG1020, Falker, Porto Alegre, Brazil). The average SRP values were extracted from 10 sub-samples that had been measured within a three-meter radius of every sampling location and that were positioned between the previous crop rows, with soil moisture content close to field capacity.
Moreover, the soil water holding capacity (SWHC) was calculated according to the model proposed by van Genutchen (1980) [37] using water retention parameters derived from a pedotransfer function (level 3) proposed by Tomasella, Hodnett, and Rossato (2000) [38]. The input variables that were used were the soil texture, organic carbon, and bulk density, with the soil volumetric water content difference between the −10 kPa and −1500 kPa water potentials that were used to represent the SWHC at a 0-0.4 m soil depth.

Plant Attributes
To quantify plant biomass across the field, destructive biomass sampling was carried out at each sampling point right after aerial thermal images were obtained. The process consisted of determining the plant population at each sampling location by counting the number of plants in four 20 m length rows. Among the counted plants, 10 individuals were randomly taken and cut at the soil surface to be immediately weighed using a portable digital balance. The fresh biomass (FBM) was obtained by multiplying the plant average weight measured in the field by the plant population.
When the crop reached physiological maturity (7 July 2019), two rows with a length of 10 m each were manually harvested to measure grain yield at the sampling locations. The ears were threshed, and the grains were oven dried at 65 • C for 48 h to calculate and correct the grain yield (YLD) for a 13% grain moisture.

Statistical Analysis
Datasets from soil and plant variables including CT were analyzed using the R software (R Development Core Team, 2020), starting with descriptive statistics, calculating the mean, minimum (Min) and maximum (Max) values, standard deviation (SD), coefficient of variation (CV) as well as the skewness and kurtosis values. Furthermore, the normality of the data was checked using the Shapiro-Wilk's test at 5% significance and when necessary, a Box-Cox transformation was applied to achieve homoscedasticity.
The relationship between the CT and sampled variables was first assessed using a Pearson correlation matrix. Moreover, simple linear regression models were built by relating the CT with the variables that related to the amount of water present in the soil (SVWC). To better understand how CT can contribute to explaining soil water status, other variables that are linked to water dynamics in soil were combined with CT as independent variables in a multiple linear regression at a significance level of 5% (F-test), with predictors being classified according to their importance based on CAR score [39].

Results
The thermal orthomosaic generated from aerial imagery presented a spatial resolution of 0.63 cm pixel −1 , with CT ranging from 23.3 to 59.2 • C across the field (Figure 3). Among the 18 sampling locations, CT varied from 32.8 to 40.6 • C.
Soil and plant attributes sampled in the study generated a total of 15 variables, which are listed in Table 1 The thermal orthomosaic generated from aerial imagery presented a spatial resolution of 0.63 cm pixel −1 , with CT ranging from 23.3 to 59.2 °C across the field (Figure 3). Among the 18 sampling locations, CT varied from 32.8 to 40.6 °C. Soil and plant attributes sampled in the study generated a total of 15 variables, which are listed in Table 1, including the results of the descriptive statistics. The results showed that silt content, SRP at 0-0.1 m, and macroporosity presented greater variation among soil physical variables, with a CV of 38.1, 28.7, and 23.3%, respectively. This variability reflects on SWHC, which ranged from 41.1 to 72.4 mm and on SVWC measurements, indicating that soil moisture content varied between 0.15 and 0.22 m 3 m −3 during the flight campaign.   The Pearson's correlations among the studied variables are listed in Table 2. As expected, CT was negatively correlated with most variables (Sand, Clay, OM, MaP, SVWC, SWHC, FBM, and YLD), indicating that higher temperature readings can be linked to lower values on these attributes. Significant correlation values (p-value ≤ 0.05) were obtained when relating CT with SVWC, SRP at 0-0.1 m, and FBM, with Pearson correlation values of −0.65, 0.58, and −0.56, respectively. Correlation was significant at 0.001 (***); 0.01 (**); and 0.05 (*) probability levels.
In order to test the feasibility of predicting the soil water content using CT information, we first implemented simple linear regression with SVWC as an independent variable ( Figure 4). The results showed that the linear regression model was statistically significant, with a p-value lower than 0.01 and an R 2 of 0.42. Moreover, a multiple linear regression was performed by combining CT and all of the soil physical parameters sampled as independent variables with SVWC as the dependent variable. The regression results are shown in Table 3, demonstrating that by adding soil physical variables, the multiple linear regression model achieved an R 2 of 0.88 (p-value ≤ 0.05). The importance of each predictor in the multiple regression model was determined based on a CAR score [39] in which the clay content, CT, and SRP at 0.2-0.3 m were the most important variables for the prediction SVWC, with CAR values of 0.289, 0.286, and 0.240, respectively.

Discussion
Since the main goal of the study was to assess the suitability of remotely sensed CT on predicting maize water status, the data acquisition and image processing that were employed focused on delivering precise and accurate temperature measurements. According to Mesas-Carrascosa et al. (2018) [40], accuracies around 1 °C are ideal when an- Moreover, a multiple linear regression was performed by combining CT and all of the soil physical parameters sampled as independent variables with SVWC as the dependent variable. The regression results are shown in Table 3, demonstrating that by adding soil physical variables, the multiple linear regression model achieved an R 2 of 0.88 (p-value ≤ 0.05). The importance of each predictor in the multiple regression model was determined based on a CAR score [39] in which the clay content, CT, and SRP at 0.2-0.3 m were the most important variables for the prediction SVWC, with CAR values of 0.289, 0.286, and 0.240, respectively.

Discussion
Since the main goal of the study was to assess the suitability of remotely sensed CT on predicting maize water status, the data acquisition and image processing that were employed focused on delivering precise and accurate temperature measurements. According to Mesas-Carrascosa et al. (2018) [40], accuracies around 1 • C are ideal when analyzing crop water status with thermal cameras, which is similar to reported values obtained with the same instrument and methodology employed in our study, where an accuracy of 1.32 • C (RMSE) and R 2 of 0.99 were observed under equivalent conditions [34].
Another important aspect regarding plant water status assessment using CT is the timing for data acquisition. In general, CT is a function of the evapotranspiration rate, which depends on the atmospheric evaporative demand and soil water storage [41]. Considering the VPD of 1.96 kPa measured during image acquisition, the meteorological condition was not restricting the evapotranspiration rate of the plants [42][43][44], and therefore, CT would be mainly influenced by the soil water status. Results from the 10-day soil water balance calculated with SWHC adjusted for each sampling location revealed a water deficiency between 3.0 and 5.8 mm during this period, which was confirmed by the soil sampling that was carried out while aerial images were taken. The SVWC results varied from 0.15 to 0.22 m 3 m −3 and compared well with the SVWC values at the permanent wilting point (1500 kPa), which, according to the pedotransfer function employed to estimate SWHC, varies between 0.13 and 0.20 m 3 m −3 among the sampling locations.
In addition to soil water deficit and VPD condition, the timing for thermal image acquisition tends to be key when assessing crop water stress. For this reason, the aerial imagery was carried out close to noontime, which is the time of the day when the leaf water potential and stomatal conductance tend to be the lowest, resulting in maximum expression of stress [20,31], which is also when the CT readings between 32.81 and 40.58  [15] reported similar results from a maize field thermal orthomosaic, with CT values starting around 20 • C ranging up to 60 • C. According to the authors, portions of the field with higher CT values occurred due to severe water stress, in which plants with curled leaves were more frequent. When the leaves are curled, the leaf area index is reduced, and that might increase the proportion of soil pixels that are able to be detected by the sensor, resulting in higher temperature values being extracted from the thermal images. The presence of terraces was found to be the main cause of the higher temperature values observed across the field, in which the accumulation of water from rainfall that occurred during early growth stages reduced the plant population dramatically on the upper portion of some terraces, with temperature readings mainly being influenced the by soil. On the other hand, lower CT readings were more prominent on flatter portions of the field and lower altitudes, where the water accumulation tends to be higher.
The Pearson's correlation analysis relating CT with plant variables was based on fresh biomass (FBM) and grain yield (YLD). Negative correlation values were obtained for both variables, with −0.56 for FBM and −0.45 for YLD. Similar values were reported by Bhandari (2016) [46] in a study assessing different maize hybrids under dryland conditions in north Texas, USA, with Person's correlation values between −0.32 and −0.64 when comparing CT and FBM. Moreover, comparable results regarding the relationship between CT and YLD were observed by Zia et al. (2013) [47] during experiments with different maize genotypes in Mexico, achieving a correlation value of −0.55.
Among the physical attributes of the soil, the clay content (r = −0.44), macroporosity (r = −0.45), and SRP (r = 0.21 to 0.58) yielded the highest correlation values. These variables were well correlated with SVWC, presenting Pearson's correlation values of 0.61, 0.37, and −0.40, respectively. The results corroborate the simple linear regression model obtained by plotting CT against SVWC, where the linear relationship was statistically significant with an R 2 of 0.42, demonstrating that CT tends to increase when the soil water content is reduced. The correlation was somewhat weaker in comparison to the values that were reported by Zhang et al. (2019) [15], in which the linear regression models between SVWC at 0-0.2 m soil depth and CT yielded R 2 values between 0.40 and 0.53. The higher R 2 values can be explained by the treatments that were tested in their study, where five levels of irrigation were applied, ranging from fully irrigated to severe deficit irrigation, resulting in a wider range of soil water contents during the assessment, which contributed to obtaining regression models with well-distributed points favoring higher correlation values. Even though the correlation value obtained with SVWC indicates that predicting soil water content using CT may lead to substantial errors, thermal orthomosaics can be used to understand the spatial variability of soil water content, helping to identify regions within the field with lower and higher SVWC.
When the CT was associated with soil physical attributes in a multiple linear regression model with SVWC set as the dependent variable, the prediction performance increased significantly, achieving an R 2 value of 0.88. According to the CAR score [39], which classifies the importance of independent variables in multiple regressions, CT was the second most important variable (CAR = 0.286) among the 11 attributes that were included in the model, along with clay content (CAR = 0.289) and SRP at the 0.2-0.3 m soil depth (CAR = 0.240). This result indicates the benefit of including CT on crop water status assessments, even when the soil physical attributes that are directly related to soil water dynamics are available. While the physical attributes of soil tend to be more stable over time, CT presents high temporal variability, matching the needs of in-season analysis.

Conclusions
From UAV-based thermal imagery, we were able to generate a thermal orthomosaic representing the spatial variability of a rainfed maize canopy temperature along a field of 40.1 ha. Among the sampling locations, temperature readings varied from 32.81 to 40.58 • C, indicating different levels of plant drought stress across the field.
Exploring the hypothesis that canopy temperature is primarily a function of soil water storage, soil sampling was carried out to determine the soil water content and other key variables that could be directly related to soil water content. The results indicated that canopy temperature alone is not sufficient to predict soil water status (R 2 = 0.42). Alternatively, the performance of soil water status prediction was substantially increased when the thermal data were associated with soil physical attributes (R 2 = 0.88), demonstrating a combination that can be used as an in-season tool for onset drought stress assessments.