Article Applicability of the Thermal Infrared Spectral Region for the Prediction of Soil Properties Across Semi-Arid Agricultural Landscapes

In this study we tested the feasibility of the thermal infrared (TIR) wavelength region (within the atmospheric window between 8 and 11.5 μm) together with the traditional solar reflective wavelengths for quantifying soil properties for coarse-textured soils from the Australian wheat belt region. These soils have very narrow ranges of texture and organic carbon contents. Soil surface spectral signatures were acquired in the laboratory, using a directional emissivity spectrometer (μFTIR) in the TIR, as well as a bidirectional reflectance spectrometer (ASD FieldSpec) for the solar reflective wavelengths (0.4–2.5 μm). Soil properties were predicted using multivariate analysis techniques (partial least square regression). The spectra were resampled to operational imaging spectroscopy sensor characteristics (HyMAP and TASI-600). To assess the relevance of specific wavelength regions in the prediction, the drivers of the PLS models were interpreted with respect to the spectral characteristics of the soils’ chemical and physical composition. The study revealed the potential of the TIR (for clay: R2 = 0.93, RMSEP = 0.66% and for sand: R2 = 0.93, RMSEP = 0.82%) and its combination with the solar reflective region (for organic carbon: R2 = 0.95, RMSEP = 0.04%) for retrieving soil properties in typical soils of semi-arid regions. The models’ drivers confirmed the opto-physical base of most of the soils’ constituents (clay minerals, silicates, iron oxides), and emphasizes the TIR’s advantage for soils with compositions dominated by quartz and kaolinite.


Introduction
Large scale monitoring of soil properties is widely demanded due to a range of global issues such as food production, climate change, and environmental degradation (http://www.globalsoilmap.net/).Remote sensing in the solar reflective spectral range (visible, VIS: 0.4-0.7 µm; near infrared, NIR: 0.7-1.1 µm; shortwave infrared, SWIR: 1.1-2.5 µm) has already been demonstrated to be beneficial for digital mapping of essential pedogenic surface properties.The use of hyperspectral techniques has also demonstrated the quantification, which is attractive for soil scientists due to the reduction in the need for time-and cost-intensive soil laboratory analyses and field campaigns [1].Several studies have demonstrated that such soil characteristics can be quantified and predicted statistically via their spectral signatures in the commonly used/accessible VNIR-SWIR wavelength region [2][3][4][5].
Although the VNIR-SWIR region has provided suitable results [6], the thermal infrared (TIR) wavelength within the atmospheric window between 8 and 14 µm has the potential to provide extended capabilities to soil scientists: In particular, mapping texture from sandy soils with very low clay contents is reaching the limits of the capabilities of the VNIR-SWIR regions.The grain size distribution of these soils is crucially driven by its major constituent mineral, quartz, but the VNIR-SWIR region is missing any distinctive Si−O related spectral features [7].These coarse textured soils should benefit from the use of the TIR wavelength region, as it can detect the strong fundamental molecular vibrations of the Si−O-stretching, which is characterised by intense spectral contrast in the Reststrahlen bands between 8 and 12 µm [8].Furthermore, to support this aim, very sensitive spectral features are traced within these Reststrahlen bands, in the presence of clay minerals (i.e., kaolinite).
Organic carbon content is typically low in semi-arid soils, whereas these soils are often dominated optically by iron bearing oxides, which are known for obscuring the spectral information used to determine the organic carbon content, as they occur in the same region of the VIS [6].Although most of the distinct spectral features of organics occur in the midwave infrared (MWIR, 3-5 µm) [9], the TIR region can provide indirect information through intercorrelations with other soil parameters, which have features in this region.Such correlations are often spread throughout the continuous bands of the spectra and are difficult to specify, but can be profitable input for prediction models.
Moreover, the prediction of clay content in the VNIR-SWIR region mainly depends on its common spectral features in the SWIR (e.g., kaolinite, illite, or smectite all possess absorption features near 2.20 µm), which can be influenced by the proximity of the cellulose feature at 2.08 µm [10].For agricultural applications, operating on conservation tillage systems such as no-till farming (often employed on semi-arid paddocks to protect the soil surface from wind erosion processes during the summer fallow-stubble retention) the cellulose influence from the vegetation residues can be significant.The use of the TIR promises more flexibility in this regard, even if it also contains a distinct cellulose spectral feature at 11.1 µm, but that is more separated from most of the clay features of interest.
Thermal infrared hyperspectral remote sensing is gaining more attention, as can be seen from increasing availability of TIR airborne imagery systems: Spatially Enhanced Broadband Array Spectrograph System (SEBASS), Hypercam (Telops), AisaOWL (Speciem), and TASI-600 (Itres).These TIR imaging sensors detect energy emitted by the target area itself.Laboratory studies aiming to gain spectral signatures, which are quantitatively comparable with remote sensing TIR imagery data, are limited to either emission FTIR spectroscopy (passive sensory mode) or using an active system with an integrating sphere for obtaining directional-hemispherical reflectance measurements (DHR) [11,12].The direct measurement of the samples emitted energy simulates the remote sensing data more closely, which can be of importance when thin coatings, such as desert varnishes, are present [11].
High-dimensional datasets, typically acquired in all disciplines of spectroscopy, require adequate statistical analyses techniques, which are able to deal with many highly collinear spectral bands (predictor variables) from relatively few observations.Therefore, ordinary multiple regression is no longer practical, due to such multicolinearities.Multivariate analysis techniques, such as Partial Least Squares Regression (PLSR) analysis, are capable of managing such associated overlapping features.PLSR is an adequate technique to handle the difficulties inherited from interpret overtones by extracting response variable relevant information from the spectra, while ignoring redundant information [13,14].
The aim of this study is to assess the potential of the TIR for modelling soil surface properties from airborne/satellite platforms.In this context we investigated the possibilities of this wavelength region to quantitatively derive soil's clay-, sand-, and organic carbon content.Based on resampled HyMap and TASI-600 spectral bands, we compared the efficacy of the TIR and VNIR-SWIR methods.For the local farmers of Mullewa and the Department of Agriculture and Food of Western Australia (DAFWA), who are dealing with the issue of soils prone to wind erosion, the outcomes of this study represent an important step towards a monitoring system of such soil properties, crucially influencing the soils' surface stability [15].To the knowledge of the authors, this is the first time emission FTIR spectroscopy has been applied for investigating the quantification of soil properties of semi-arid areas in Western Australia vulnerable to wind erosion.The few previous studies, in which such soil properties have been predicted with PLSR from the TIR are all based on standard DRIFT optic measurements [16][17][18], which do not comply with the requirements necessary for a quantitative comparison of TIR remote sensing imagery data.

Soil Sampling and Sample Analyses
Subsurface (upper ∼2 cm) soil samples (n = 55) were collected during a field campaign in February 2010 from paddocks within the investigation area to the west of the township of Mullewa (28 • 32 15 S, 115 • 30 25 E), which is located in the semi-arid wheat belt region of Western Australia.The area is dominated by an undulating sand plain system, which includes yellow and red sands and an alluvial valley system with relict red loams over a red-brown hardpan, which are classified as Tenosols and Kandosols, respectively, in the Australian Soil Classification [19].A few additional samples were collected further south of the area, where grey sands (Tenosols) and alkaline clays (Vertisols) were sampled to cover a wider range of soil characteristics.
Two datasets were selected from the soil samples collected at Mullewa.An extended dataset, which includes the extreme clay rich samples (n = 55) and the original dataset (n = 45), which includes only samples that are representative for the mineralogy of the Mullewa paddocks.The original dataset was divided into a training or calibration dataset (n = 30) and an external validation dataset (n = 15) for the multivariate statistical analysis.
The soil samples were air dried and passed through a 1 mm sieve.The grain size distribution was determined by using the standard pipette method following Stokes law [20], to determine the fractions: clay (<2 µm), silt (2-20 µm), and subsumed sand (20-1,000 µm).Soil organic carbon content (%) was measured using the Walkley Black method [21].The inherent mineralogy and chemistry of the soils were determined from XRD analyses.Note that an assumption in this study was that the clay mineral content, which is the spectral driver in the PLSR models, and the clay texture, derived by Stokes Method, are essentially the same and that non-clay minerals such as very fine quartz, iron oxide particles, or organic material less than 2 µm are insignificant.In this study area the XRD results produced only traces of iron oxide and the lab analysis for organic carbon gave only 0.09% to 1.11%, which makes the assumption fair in this situation.For other areas with increased iron oxide and organic carbon a pretreatment to Stokes Method is recommended to improve PLSR relationships.

TIR-Emission FTIR Spectroscopy
Soil radiance spectra (passive mode emission measurements) in the TIR were recorded at Commonwealth Scientific and Industrial Research Organisation (CSIRO) laboratories (Perth, Western Australia) using a micro Fourier Transform Interferometer (µFTIR) (model 102, Design&Prototypes instruments, http://www.dpinstruments.com)with a spectral resolution of 6 cm −1 via a fore optic with a 4.5 • field of view from ∼250 mm height, resulting in an GFOV of ∼20 mm.The instrument was equipped with a liquid-nitrogen-cooled two element detector (dual InSb/MCT) and measured radiance within the spectral range from 2 to 16 µm.However the range used in this study was limited to the wavelength within the atmospheric window from 7 to 14 µm, restricted by the range of the MCT detector.
(1) and (2): The µFTIR was radiometrically calibrated using a blackbody (model 41P, D&P instruments) attached to the fore optics, adjusted for two different temperatures.For the first reading it was set just above the sample temperature (T warm : 65 • C) and for the second one just below the ambient temperature (T cold : 25 • C).The temperature settings were held within this region to minimize potential nonlinearities from the instrument [25].
Offsets, incorporating electrical factors and the instrument's self-emission, could be retrieved from either blackbody by Equation (2): (3) and (4): During the measurement with the µFTIR the background downwelling radiance, which is emitted from the ceiling, the walls and the atmosphere (in particular above the sample) in the laboratory, is reflected from the surface of the sample and recorded by the instrument on top of the sample radiance measurement.Downwelling radiance was measured via a diffuse reflective brass plate just before the soil sample reading.
The brass plate is spectrally flat and has a emissivity of ε∼0.25.Following Kirchhoff's law (ε = 1 -R) [26], the complementary reflectance is R∼0.75, and accordingly can be assumed that the at sensor radiance recorded by the instrument, when the brass plate is measured, is a good estimate for the downwelling radiance (L DW (λ)), which was calculated from Equation (3): where L BP (λ) is the spectral radiance, and ε BP (λ) the emissivity of the brass plate.B(T BP , λ) is the Planck blackbody spectral radiance of the brass plate at its temperature T BP .Note, that for these measurements, highly defuse gold plates (ε∼0.03-0.06)are also common.L DW (λ) measured off the brass plate will be more accurately calibrated to radiance units than off a gold plate when outdoor measurements are done under extreme cold background temperature conditions [27].For laboratory measurements, however, where the plate's temperature is close to the ambient temperature, the absolute emissivity/reflectance of the plate is not an issue anymore, as long as the background temperature is estimated well enough.The background temperature was measured by a thermocouple thermometer physically attached to the back of the brass plate and provided ±0.1 • C accuracy.Finally, the sample surface was measured, which was preheated (60 • C) before in an oven.The measurement of relatively hot samples above the background room temperature ensured an improved signal to noise ratio (SNR).The instruments integration time (16 scan cycles) was adjusted to keep its temperature decrease to a minimum during recording time.Having calculated the soil surface radiance (L S (λ)), its emissivity (ε S (λ)), corrected for downwelling radiance, could be derived by Equation (4).
However, solving Equation ( 4) requires the soil surface temperature to be known (which is especially difficult for remote sensing data).For that purpose an inverse Planck curve fitting algorithm was applied to estimate the temperature.Assuming that for a specific wavelength range the emissivity is constant and known, identical temperatures for these wavelengths could be derived by using the inverse Planck function as a function of radiance.For the soils in the Mullewa region, which are dominated by silicate minerals, this assumption is true, where typically the silicate Christiansen frequency features occur.In this region the emissivity values are at their highest (close to ε = 1, see Figure 1), and thus the highest apparent radiant temperatures of the sample.Considering increasing instrument noise at shorter wavelength the range between 7.7 and 7.8 µm was selected for the temperature emissivity separation algorithm.Within this range the algorithm identified the maximum emissivity, approaching 1.0.Software developed by Andy Green (pers.comm.) of CSIRO was used to apply the above processing steps.Note that the path atmospheric spectra transmissivity between the target and the instrument has not been considered in the above equations, as it becomes negligible for short distances, as has been used in our setup [23].Bands between 7 and 7.63 µm were excluded from the spectral range in further analyses due to extreme noise levels, resulting in 208 spectral bands.

VNIR-SWIR-Bidirectional Measurements
Soil surface bidirectional reflectance in the VNIR-SWIR spectral region was measured in the laboratory on the samples using an Analytical Spectral Devices (ASD) Inc FieldSpec-Pro (http://www.asdi.com)spectrometer with a contact probe (GFOV∼10 mm).The Spectral resolution is 3-4 nm in the 0.35-to 1.0-µm region (spectral sampling 1.4 nm), and 10-12 nm in the 1.0-to 2.5-µm region (spectral sampling 2 nm).The entire spectrum is resampled at 1 nm for display purposes, which resulted in 2,151 spectral bands.The measurements were performed straight after the TIR measurements, without any interference to ensure that the soil surfaces were similar as possible for both TIR and VNIR-SWIR measurements.The ASD spectra were subsequently corrected to absolute reflectance using the white reference measurement of a spectralon panel (http://www.labsphere.com/)taken before soil radiance measurements.

Spectral Resampling to Operational Hyperspectral Remote Sensors
The spectra were resampled to operational imagery sensor spectral specifications for both wavelength regions, the VNIR-SWIR and the thermal infrared range.ASD spectra were resampled to HyMap spectral characteristics (2010 HyVISTA Co., http://www.hyvista.com),with 125 bands and a full width half maximum (FWHM) between 15-20 nm, covering the spectral range from 450-to 2,480-nm.The TIR spectra were resampled to TASI-600 specifications (2009 ITRES Research Limited), which possess 32 bands with a FWHM of 125 nm and a spectral range from 8 to 11.4 µm.Figure 1 shows the defined sensors' spectral response functions (SRF), which follows their defined spectral characteristics.
HyMAP and TASI-600 were chosen as representative hyperspectral sensors that possessed specifications most suited to the aim of this study.Data from the HyMAP instrument has been used extensively in the past 10-15 years in all areas of the world and is nowadays one of the most used and known airborne hyperspectral imager for the VNIR-SWIR region.The TASI-600 is one of the very few airborne TIR sensors currently available.The spatially enhanced broadband array spectrograph system (SEBASS) covers a broader spectral range from 7.5 to 13.5 µm (128 contiguous spectral bands).We have decided to test the TIR potential by resampling to the TASI-600 sensor characteristics for two reasons.The TASI-600 is more accessible than the SEBASS sensor.We also wanted to see whether the slightly reduced spectral coverage and resolution of TASI-600 is sufficient for the aim of this study.
The TIR spectra were converted to reflectance using Kirchhoff's law.The resampled TIR and the VNIR-SWIR laboratory soil spectra were transformed from reflectance (R) to absorbance (log 10 (1/R)), as input to the statistical analyses, to be consistent with the Beer-Lambert law [28].

Multivariate Regression Modelling
Prediction models were established, using PLSR analysis based on the NIPALS algorithm [13].The dataset for a PLSR basically consists of two matrices, X and Y , where X represents the measured soil spectra (independent or predictor variables) and Y the known quantities of the soil properties (dependent or response variables).The PLSR comprises two independent calculation stages, the calibration and the validation (internal or external).

PLSR-Calibration and Validation
The 45 samples from the original dataset were divided into a training or calibration dataset of 30 samples, and 15 samples were used as a test dataset for an independent external validation.The extended dataset contains all of the 55 samples, including the extreme clay rich samples, thus did not allow a representative individual test dataset for an independent external validation.For this case, an internal validation method (full cross-validation) was used to assess the prediction accuracy of the PLSR models.
In the calibration, PLSR establishes a multivariate regression model, based on the training or calibration dataset.This involved the reducing of the high dimensional X data to a few factors (latent variables), which concentrate the spectral information that is relevant to predict Y .This essential information is stored within the first factors, while higher factors contain decreasing Y -relevance and increasing noise data.
To assess the models performance with respect to its dimensionality, a full cross validation (leave-one-out validation) was employed.The models were evaluated based on the following properties: lowest number of factors as a measure of model robustness, the lowest root mean square error of prediction (RMSEP) as a measure of expected error in future predictions, and the coefficient of determination (R 2 ).
Predicted versus observed plots are used to display the Y -values, which are predicted from the regression models.The closer the trend line of the calibration stage follows the trend line of the validation stage, the more robust the model can be rated in terms of the samples used for the validation.

Spectral Interpretation of the PLSR Models
The influence of each of the spectral bands in terms of their informative contribution for predicting the soil parameter of interest was assessed by interpreting the pattern of the model's loading weights and its regression coefficients.
Loading weights (w) describe the relationships between the Xand the Y -variables.High loading weights (in the positive and in the negative dimension) indicate spectral bands being important in this relationship, and hence have a positive or negative link with the predicted soil parameter in that specific factor.Loading weights from different factors can be judged by their proportion of the total explained variance [14].If the loading weights as a function of wavelength reveal spectral characteristics related to known physical and chemical sample constituents, they can be defined as the spectral driver for the model within that factor.
For the final PLSR model, the regression coefficients are generated, which can be interpreted spectrally in the same way as the loading weights.However, they combine the loading weights information from each of the selected factors.Thus, the regression coefficients reveal the most influential spectral bands of the selected PLSR model for the soil parameter of interest.

Soil Analyses
Analysis of the 45 samples within the original dataset shows that the clay content ranges were between 3.8% and 14.8%, and sand content between 77.3% and 95.3%.Organic carbon content ranges were between 0.09% and 1.11%.The extended dataset (n = 55), including the extreme clay rich samples, has a much broader textural range.Figure 2 displays the descriptive statistics of the determined soil properties.The asymmetric distribution towards coarse textured soil samples is apparent.Figure 2(c,d) shows that the textural properties are highly correlated and are both independent of the organic carbon content, which appears from their cross correlation.Based on the XRD and grain size analyses, the majority of the samples are coarse in texture, dominated by quartz, with minor amounts of clay-sized minerals of kaolinite, illite, feldspar and smectite.Some traces of iron oxides and calcite were also detected.Altogether the soil sample set represents a typical distribution of the Western Australian wheat belt region, with sandy soils as the dominant composition.

Soil Spectral Characteristics
The measured spectra show absorption features typical for the ancient, deeply weathered soils from the semi-arid regions in Australia.Figure 3 displays their spectral variations for a selection of the soil samples in the VNIR-SWIR (as reflectance) and the TIR (as emissivity) regions.The soil spectra show iron-oxide (hematite) and -hydroxide (goethite) features in the VNIR wavelength region (between 0.5 and 1.0 µm), as well as clay features in the SWIR.The doublet feature close to 2.2 microns is distinctive for kaolinite [1], which is just as abundant as iron oxides in most Australian soils.They are among the few remaining minerals from predominant highly weathered and leached Australian mantle [29].Another residual mineral in this context is quartz, particularly for the sandy soils of Western Australia.But quartz does not provide any usable feature in the VNIR-SWIR wavelength region, as can be seen from the dashed spectrum (99.9% SiO 2 content) from Figure 3(a).
However, quartz displays a broad emissivity doublet in the TIR between 8 and 10 µm (Figure 3(b)).This doublet is the assertive feature within the Reststrahlen bands (>7.5 and <12.3 µm), which are associated with the strong fundamental molecular vibrations of the Si−O stretching and exhibit strong spectral contrast [8].Furthermore, located just before the start of the Reststrahlen bands, there appears another very distinct spectral feature, typical for quartz, the Christiansen Frequency Feature (CFF).The first CCF shows up at around 7.5 µm (a second, smaller one, appears at around 12.3 µm), where the refractive index of quartz approaches that of the surrounding medium.For these CCF wavelengths scattering is minimum, penetration into the material and thus emissivity in contrast is maximum [30].The majority of the TIR spectra from the Mullewa soils display the features prominently.
The quartz doublet, can appear to display an additional lobe at 9 µm (e.g., a triplet) with the presence of clay minerals, such as kaolinite, montmorillonite, or illite, in the soils.Kaolinite has the strongest effect in reducing the quartz Reststrahlen signature at that specific wavelength (Figure 4), as its absorption contrast is most pronounced.An additional clay mineral feature is found at 8.3 µm (with a small spectral shift among the three of them).Several additional kaolinite features can be identified, which helps to discriminate it from other clay minerals at 9.8 µm, at 10.2 µm and a broader one from 10.6-10.8µm.The kaolinite dominance in the Mullewa soils' clay minerals is spectrally distinctive, as all of its features are represented in most of the Mullewa soil spectra, and verified throughout the XRD analyses.Iron-bearing oxides also have diagnostic features within the TIR as shown in Figure 4. Hematite and goethite, which are both present among most of the Mullewa soils, do show their abundance in the quartz Reststrahlen bands.The iron influence appears to reduce the reflectance of the first lobe of the quartz, with increasing effect towards the longer wavelength side.This can result in a peak shift from the short wavelength side (pure quartz) to the long wavelength side (iron influence) within the first lobe, as can be seen for the Mullewa soil spectra in Figure 3(b).(Note that the figure shows emissivity values, so the peaks are displayed as trough in this context.)The effect of this result can be associated with desert varnish coating on the quartz sand grains, as described by Lyon [32], Rivard et al. [33].
However, this feature does not assist in distinguishing between these two iron minerals.Goethite has two diagnostic reflectance peaks at 11 and 12.3 µm, which are absent for the hematite.Whereas the feature at 11 µm comes clearly through in many of the soil spectra, the 12.3 µm feature probably gets masked by the more dominant second CFF, which appears at the very same spot (note that the TASI-resampled spectra will not cover the wavelength of this feature anymore).
Although the spectra were resampled to imaging sensor spectral characteristics and some of the spectral information was lost, the relevant spectral information of the described features remained in the spectra.Figure 5 shows the spectral resampled (HyMap and TASI-600, respectively) and transformed (to log 10 (1/R)) soil absorption spectra for the VNIR-SWIR and the TIR as they were entered into the prediction modelling.

Prediction of Soil Properties
Calibration and validation were performed for the selected soil properties (% clay (<2 µm), % sand (20-2,000 µm), and % organic carbon) as the response variables using the resampled and transformed absorption spectra of the VNIR-SWIR, the TIR, and their combination (VNIR-SWIR-TIR) as predictors for the modelling.The PLSR results are summarised in Table 1 and discussed in the next sections.
Table 1.Summary of PLSR modelling validation results (internal cross-validation and external validation using a test dataset) for each of the soil properties in the specific spectral region for the original and the extended dataset: number of samples (n), number of used model factors (f ), coefficient of determination (R 2 ), root mean square error of prediction (RM SEP ), not possible due to unstable model (n.p.).

Clay Content
All of the three spectral ranges (VNIR-SWIR, TIR, and VNIR-SWIR-TIR) were able to predict the soil's very small range of clay content (<15%) from the original dataset.The TIR performed best amongst them, as it was able to predict the clay content with a 3 factor model (Figure 6).Applying this model to the test dataset, the clay content could be predicted with a R 2 of 0.93 and a RM SEP of 0.66%.The combination of the TIR with the VNIR-SWIR spectral range did not improve the TIR performance, since the VNIR-SWIR-TIR model needed 5 factors to predict the clay content with a R 2 of 0.84 and a RM SEP of 0.98%.Using only the VNIR-SWIR spectral range did result in a slightly better performance than its combination with the TIR (5 factor model with a R 2 of 0.89 and a RM SEP of 0.81%).

Sand Content
The models for predicting the sand content of the soils in the original dataset revealed good results for the TIR and the VNIR-SWIR-TIR combination (Figure 7).The PLSR analysis was not able to create an adequate model from the VNIR-SWIR wavelength alone.Similar to the clay content, the TIR model on its own performed best, as 3 factors were sufficient to built a stable model, which predicted the sand content with a R 2 of 0.93 and a RM SEP of 0.82%.The VNIR-SWIR-TIR produced the best results with 5 factor model (R 2 = 0.9 and RM SEP = 0.97%).

Organic Carbon Content
The soil's content of organic carbon could be predicted from all of the three spectral ranges, despite of the very low content of <1.5%.The TIR, which does not contain diagnostic organic spectral features, produced good results in the cross-validation stage (5 factors, R 2 of 0.91 and a RM SEP of 0.08%).However, the model resulted in moderate prediction abilities for the test dataset (R 2 of 0.62 and a RM SEP of 0.12%).The best result was achieved with the combination of the VNIR-SWIR-TIR spectral range (Figure 8).A PLSR model with 5 factors was able to predict the organic carbon content with a R 2 of 0.95 and a RM SEP of 0.04%.()(*+,$<(=(>+&'4"*(((   !   Figure 9 provides an insight to the drivers (also referred to as "inner relations") for the VNIR-SWIR-TIR models.The regression coefficients for the combined VNIR-SWIR-TIR models identify the significance of a spectral region for the prediction of the soil parameter of interest.Similar to the interpretation of the loading weights within a specific factor, the correlation coefficients relate to the spectral bands of importance for the final all factor combined PLSR model.As expected, the variables clay and sand reveal their inverse relationship in their regression coefficient patterns, at least within the most important regions of the TIR and the SWIR. Figure 9 shows that the SWIR and the TIR are the most important spectral regions for the prediction of clay and sand content, which fit with the occurrence of the most distinctive textural spectral features (clay minerals, such as kaolinite and silicates, and quartz) for these wavelengths.In spite of the regression coefficients being composed of the combined information from all of the used loading weights, they still are consistent with the diagnostic spectral features.For example, the kaolinite features at 2.2 µm, at 9.0 µm, at 9.8 µm, and between 10.6-10.8µm (red arrows).There is a subtle suggestion of the 2.2 µm kaolinite doublet feature too.
Furthermore the first lobe of the quartz Reststrahlen bands starting at 8 µm displays a steep stair-like descent towards longer wavelength (blue arrows), indicating its prominent role as the driver for grain size variations.This effect is likely to be associated with the reduction in spectral contrast of the Reststrahlen bands with decreasing grain size [30,34].However, it should be noted that more complicated interactions of radiation with mixtures of fine quartz and/or clay particles coating quartz sand grains can also reduce the Reststrahlen spectral contrast [35].
The organic carbon content PLSR model indicates that the TIR wavelengths are valuable input for its prediction.However, the most important spectral bands are found in the VNIR region, which is consistent with the literature, where spectral features in the VIS region are identified as being related to the presence of organic compounds (green arrow).The change in slope in the VIS is a function of the soil organic matter content [6,36], which is represented as increasing regression coefficient values from 588 nm to 710 nm (Figure 9).Subsequently they decrease towards a trough with a minimum at 900 nm, which is associated with features of iron bearing minerals such as haematite (880 nm) and goethite (920 nm).The feature located at around 500 nm is also due to the iron oxide charge transfer feature in the ultraviolet/visible wavelength region [37].Note, since the soils' iron oxide content has not been measured to this point of the research, we cannot ascertain if the relationship between soil organic carbon and the reflectance spectra is due to an association with this or any other unmeasured parameter.
The interpretation of the PLSR models' loading weights for a particular factor reveals a more detailed view of the models drivers relative to the shape of the soils' spectral features.With the focus on the TIR spectral region, the texture's inverse relationship of the two response variables sand and clay is clearly evident again.Figure 10 shows the loading weights for the 3 factors (f1, f2, and f3), involved with both models.The loading weights for sand (start at negative values) and clay (start at positive values) are following very similar patterns in each of the model's factors, however the same spectral information is used inversely in the prediction.Factor 1 (dotted lines) represent the overall trend of the spectral contrast, mainly driven by the quartz Reststrahlen bands, since there are no other spectral features evident in the loading weights pattern (note that the distinct 8.65 µm feature of quartz is not considered in f1, f2, and f3).Factor 2 is already dominated by the spectral features of kaolinite (9, 9.8, 10.2, and 10.8 µm).Most of these kaolinite feature are being confirmed in factor 3, which also contains the clay feature at 8.3 µm.This also reveals the influence of the goethite feature just before 11 µm, which might be interfering with the clay minerals information in that region.4. Discussion

Extended Data Set
The extended dataset included an additional few clay rich soil samples, which were more common south of the investigation area (clay content 60%-70%) and a few rarely found among the investigated paddocks (clay content 30%-40%).As they are also underrepresented among the collected soil samples, they have been analysed separately within the extended dataset.The PLSR results show that for both, the sand and clay content, stable prediction models could be established with a similar number of model factors as in the original dataset.However, due to the unbalanced distribution, the RM SEP is higher (in the TIR: RM SEP of 2.55% for clay-and RM SEP of 2.73% for sand content).Still the TIR models (3 factors) are able to predict with a smaller number of model factors than the combined VNIR-SWIR-TIR models (5 factors).

Silt Fraction
Since the silt fraction (2-20 µm) represents the complementary part to the sum of the clay and the sand content, we decided to focus this paper on clay and sand prediction.As can be seen from the modelling results of the two other texture properties, their prediction performance always depends on both the mineralogical and the grain size information.The silt however does not contain such distinct mineralogical characteristics (such as clay minerals, or quartz) but is rather a mixture of them.For the Mullewa soils, it was shown that the silt fraction is actually more likely to be related to the mineralogy of the sand fraction [29].Thus, we would suggest that the indirect prediction of % silt, via the sand and clay content, would be a feasible method, as their spectral drivers are more distinct and thus more powerful for the modelling of a soil's texture characteristics.

Conclusions
The study reveals that the prediction of all three soil properties (% Corg, % clay, and % sand) is promising for the TIR based on resampled TASI-600 spectral bands, in spite of the small textural range for the sand-dominated Mullewa soils.Good results were achieved for the properties of clay and the sand content based on the TIR spectral signatures with a R 2 = 0.93 for % clay (RM SEP = 0.66%), and 0.93 for % sand (RM SEP = 0.82%).The high spectral contrast of the TIR spectral region, as well as the distinct spectral features for both, clay minerals and quartz, were responsible for good quantification performances using the PLSR analyses.
This paper shows that predicting clay and sand content can be achieved more accurately and reliably using the TIR alone than in combination with the VNIR-SWIR region, as the 3 factor TIR models were more robust than the 5 factor VNIR-SWIR-TIR models.Modelling with the TIR could be realised with a reduced number of model factors, although the TIR provides much less numbers of spectral bands (1/3 of the VNIR-SWIR).
The PLSR analyses were not able to create an adequate model from the VNIR-SWIR wavelength alone to predict the sand content, which as we interpret is related to the very sandy character of the soil samples (>75% sand).In this case, the quartz dominance does not allow enough spectral features within the VNIR-SWIR wavelengths.
Even if the TIR model for the organic carbon content showed promising results in the internal cross-validation (R 2 = 0.91 and RM SEP = 0.08%), the external validation based on the test dataset did not give adequate results (R 2 = 0.62 and RM SEP = 0.12%).However, the VNIR-SWIR spectral range was able to predict it with a R 2 of 0.81 and RM SEP of 0.06%).The VNIR-SWIR-TIR combination revealed the best prediction result for organic carbon content, as the PLSR model (5 factors) performed with a R 2 of 0.95 and a RM SEP of 0.04%.
Overall this paper shows that the TIR region, using 32 spectral bands, has a high potential for the prediction of the clay-and sand-content across semi-arid sandy agricultural landscapes, which possess an extremely narrow textural range.Based on emission IR spectroscopy measurements, the results meet the requirements for a quantitative comparison with TIR remote sensing imagery data.In combination with the VNIR-SWIR wavelength, the VNIR-SWIR-TIR models also displayed great propose in the prediction of the low ranges in organic carbon content, typical for such semi-arid soils.Prediction from remote sensing data has many more limitations and has numerous perturbing effects (e.g., atmospheric influences, lower signal to noise ratio, coarser spatial resolutions, spectral intermixtures within pixels, and the influence of other soil variables such as soil roughness and moisture).Here we prospect that the additional use of TIR sensing will lead to a more flexible approach for quantifying soil properties.For predicting texture characteristics from very sandy soils, the TIR bands hold the promise of greater prediction accuracy through higher spectral contrast.Regarding this, the availability of airborne and future satellite TIR systems with an increased spectral range/resolution, such as SEBASS, or NASA's planned Hyperspectral Infrared Imager (HyspIRI) satellite system, will be a big asset for global digital soil mapping.

Figure 2 .
Figure 2. Descriptive statistics for the original (a) and the extended dataset (b): Min, Max, Mean, Median, and Std Deviation of the used soil properties % organic carbon (Corg), % clay, and % sand, and the cross correlations for the the original (c) and the extended dataset (d).

Figure 3 .
Figure 3. Spectral variations for a selection of the Mullewa soil samples and their corresponding parameters (% clay, % sand, and % Corg), in the (a): VNIR-SWIR (reflectance) and (b): the TIR (emissivity) region.The arrow indicates the position of the Christiansen Frequency Feature (CFF).

Figure 4 .
Figure 4. MWIR and TIR spectral reflectance signatures (DHR-measurements) of typical arid and semi-arid soil minerals, from the ASTER spectral library [31].Grey and grey-pink bars indicate the spectral characteristics of the TIR bands from the satellites ASTER and HyspIRI (TIR module on a planned NASA mission).

Figure 5 .
Figure 5. Soil absorption (transformed to log 10 (1/R)) spectra for the VNIR-SWIR and the TIR, resampled to HyMap and TASI-600 spectral characteristics, as they are entered into the prediction modelling.Locations of the spectral soil features (RB: Reststrahlen bands).

Figure 6 .
Figure 6.Observed (reference) clay content versus clay content predicted from the PLSR models for the TIR and the VNIR-SWIR-TIR spectral regions.(a) calibration (red) and internal cross-validation (blue); (b) external validation using a test dataset (blue) for the original dataset; (c) the results for the external dataset, calibration (red) and internal cross-validation (blue), where the dashed bordered area displays the data range from the original dataset.

Figure 7 .
Figure 7. Observed (reference) sand content versus sand content predicted from the PLSR models for the TIR and the VNIR-SWIR-TIR spectral regions.(a) calibration (red) and internal cross-validation (blue); (b) external validation using a test dataset (blue) for the original dataset; (c) the results for the external dataset, calibration (red) and internal cross-validation (blue), where the dashed bordered area displays the data range from the original dataset.

Figure 8 . 4 .
Figure 8. Observed (reference) organic carbon content versus organic carbon content predicted from the PLSR models for the TIR and the VNIR-SWIR-TIR spectral regions.(a) calibration (red) and internal cross-validation (blue); (b) external validation using a test dataset (blue) for the original dataset.

Figure 9 .
Figure 9. PLSR regression coefficients for the combined VNIR-SWIR-TIR models for the prediction of clay, sand and organic carbon.The arrows indicate the corresponding spectral features (red: kaolinite; blue: quartz; green: organic matter).

Figure 10 .
Figure 10.PLSR loading weights (factors 1-3) for the prediction of clay and sand in the TIR spectral region, and spectra of the relevant minerals from the ASTER spectral library[31], resampled to TASI-600 spectral characteristics.