Mapping stressed wheat plants by soil aluminum effect using C-band SAR images: implications for plant growth and grain quality

Under toxic aluminum (Al) levels in the soil, wheat (Triticum aestivum L.) suffers stress and plant growth is affected. A method for diagnosis of plants is proposed that includes the following as a strategy: to analyze total Al in the soil, employ satellite radar imagery and calculate a vegetation index. The objective of this research, conducted at the field scale, was to explore how radar backscattering coefficients from a winter wheat canopy, combined with the normalized difference vegetation index (NDVI) and geographic information system (GIS) technology, can be used as a mapping tool for the variability of Al-stressed canopies. As a result, an analysis of covariance showed significant differences, and the lowest plant height was obtained at a high level of soil Al, as well as the minimum grain weight and magnesium content. It was found that a simple model could be used to estimate plant height from the backscattering coefficient of vertical transmit-vertical receive polarization (σ°VV), with a strong correlation (r − 0.84). In turn, a third-order polynomial regression model (R2 0.70) was proposed to estimate the NDVI from σ°VV. This model provided a good estimate of the NDVI at the stem elongation stage of growth (50 days after sowing). Detected NDVI patterns were associated with variation in canopy stress depending on polarimetric information, which, in turn, was related to soil Al levels. Thus, the maps derived from the model can monitor spatial variability, where NDVI values < 0.68 indicate stressed areas. This study provides guidance for in-season stress spatial variability caused by Al.


Introduction
Cereals are by far the world's main food source, both for direct and indirect human consumption. Wheat (Triticum aestivum L.) is the most widely consumed cereal in the world, accounting for approximately 31% of the total. In the European Union, according to area and production, Spain ranks sixth (after France, Germany, Poland, the United Kingdom and Italy) producing 6.3 Mt from around 1.9 Mha, with approximately 6% of this extension land classified as having acidic soils. Acidic soils are widespread worldwide; they cover approximately 50% of potential arable land and contain aluminum (Al).
The presence of high Al levels in the soil is associated with mechanisms that determine root growth. Thus, Al induces alterations in the molecules of the genetic material, which has consequences in mitosis, causing modifications in radicle elongation and changes in the external appearance of roots (they become thicker) and considerably reducing the contact area of the epidermis with the soil solution (Xue et al., 2003). Plant roots are particularly susceptible to Al stress.
A relationship has also been found between the soil Al concentration and the phenological stage at which stress occurs. It is reported to hinder root growth in the first 10 days following the coleoptile emergence phase, varying in intensity according to the degree of tolerance of the cultivar, which is a dominant hereditary characteristic (Aniol, 1990). It also interferes with the absorption, transport and use of essential elements such as phosphorus (P), potassium (K), calcium (Ca) and magnesium (Mg), as well as in the absorption of water (Lofton et al., 2010). During leaf development, it has been observed that the Al concentration increases in plant tissues, with damage occurring in the mitochondrial membrane of the chloroplast, and the chlorophyll content in leaves is reduced (Meriño-Gergichevich et al., 2010). At maturity, there is a reduction in yield components (number of grains per spike and thousand-grain weight).
To monitor large areas of soil and identify their limitations, it is particularly convenient to use data acquisition techniques by remote sensing (RS), in which a sensor is used for the electromagnetic radiation measurement emitted, which is reflected by objects. Sensors are classified into two types of technologies: passive (optical images) and active (radar images). Using this approach, satellite imagery in RS studies is limited by various factors: optical image acquisition is only possible during the day (source of emission, the sun), but clouds and fog can affect the reflectance captured by the sensor. However, radar images have the advantage of operating at night (own emission source) and under conditions of cloudiness, fog and rain.
A plant stressed modifies in multiple ways the energy flow through the leaf, with the result that the absorption, reflectance and transmittance properties of leaves are changed (Lichtenthaler, 1996). Hence, vegetation indices can be estimated to monitor crop stress, calculated from reflectivity values at different wavelengths. The normalized difference vegetation index (NDVI) is the most commonly used index and is calculated from the difference of two spectral bands (red and near infrared), with values between − 1 and 1. Furthermore, Veloso et al. (2017) compared C-band synthetic aperture radar (SAR) image data correlated to NDVI estimations from optical data and demonstrated that useful information on crop development was obtained. This insight has also been corroborated by recent research (e.g., Filgueiras et al., 2019).
Optical RS satellites have been widely used to map the spatial distribution of some soil properties, such as pH, electrical conductivity or cation exchange capacity. However, this technology has not been used to study spatial distribution of soil Al. Recently, a model was proposed on the basis that the soil Al content could be estimated and mapped using spectral bands from satellite imagery (Hernández & Francisco-Bethencourt, 2017). Additionally, optical RS techniques have been used to study changes in spectral reflectance patterns as a result of water deficit, nutrient deficiency and crop diseases (Campos et al., 2018). Compared to optical RS, radar RS techniques focus principally on monitoring soil moisture, crop classification, crop growth and its water content. To date, the radar RS technique has not been applied in research on plant stressed by Al toxicity.
Furthermore, the success of RS in agriculture depends strongly upon an efficient and accurate method for determining within-field variation in soil properties. It has been proven that RS is a tool necessary for multiple studies and has the advantages of low cost, rapidity and relatively high spatial resolution; e.g., soil maps that can enable lime application in specific sites according to the spatial variability of acidity are needed (Adamchuk et al., 2004;Aggelopoulou et al., 2011). In fact, these soils will require applications of lime periodically to neutralize acidity. Additionally, an appropriate use of calcium will displace Al 3+ , which causes a decline in crop yields. The literature review shows that a variable lime application rate can lead to improved soil acidity at a low cost (e.g., Holland et al., 2019).
The current study focuses on the potential of radar images for monitoring wheat cultivation areas with high cloud coverage. The within-field variability of soil Al and its effects were evaluated. For this, SAR backscatter and NDVI were used, as also employed by other researchers (Bousbih et al., 2017;Filgueiras et al., 2019). The overall goal was to develop an alternative method to estimate and map the canopy of wheat plants stressed by soil Al. The specific objectives were to (i) evaluate the relationship between phenological wheat development and NDVI and how it is affected by soil Al levels, (ii) investigate the correlation of reflectance and backscatter values with NDVI and plant height in these soil conditions, (iii) develop a model for the estimation of NDVI through polarimetric data, (iv) determine whether spatial variability in canopy backscatter follows a pattern of stressassociated changes by soil Al, and (v) analyze the effect of three levels of soil Al on plant height, yield components and grain quality.

Site description and soil sampling criteria
Tenerife (Spain) is an island in the Canary archipelago, situated in the Eastern Central Atlantic Ocean (27° 60′ to 28° 35′ lat. N; 16° 05′ to 16° 55′ long. W), with a volcanic origin. The study was conducted at Icod el Alto (municipality Los Realejos) over a portion of territory at an altitude between 520 and 1000 masl and comprised an area of 556 ha (Fig. 1). This area, typical of agricultural land use in winter for cereal production, is characterized by a rainfed production system and local variety cultivation. It has a very stable climate throughout the year because it is under the direct influence of the trade wind belt, with an annual precipitation of approximately 662 mm and an average annual temperature of 15.1 °C.
In the study area, the soils are classified in the alfisols order (USDA soil taxonomy) and ustic moisture regime; moreover, acidification together with the mineral transformation of allophane and imogolite to crystalline phyllosilicates and gibbsite would lead to the predominance of Al, according to Tejedor et al. (2013).
Previous field work was carried out, consisting of visiting agriculture farms sown with wheat in the study area and identifying each farm by taking a point using a handheld GNSS device (Magellan™ Meridian, San Dimas, USA) with an accuracy < 8 m. The boundary data of each farm were obtained from the geographic information system of agricultural parcels (SIGPAC) of the Spain Ministry of Agriculture, Fisheries and Food (URL: http:// wms. mapama. es/ wms/ wms. aspx). A set of 207 polygons was manually drawn by employing GIS tools and mapped using QGIS software (QGIS Development Team, 2018). These polygons were small, with an average surface area of 1314.27 m 2 and a total of 27.21 ha wheat (Fig. 1b). To determine the spatial distribution of soil Al and its effect on wheat development, sampling plots were selected from these polygons on 30 farms according to the following procedure: First, the information from total Al levels in the soils from the mapping that was obtained in a previous study made by Hernández and Francisco-Bethencourt (2017) was employed. The mapped polygons contained topsoil (200 mm) data categorized into three classes < 3, 3-4 and > 4 g kg −1 . These classes are defined as qualitatively homogeneous soil units (Hengl et al., 2003). It is considered that entities conform to the same order (alfisols), a uniform climatic environment (nearest weather station, REALE to < 2.5 km) and have average altitude variation less than 110 masl (with respect to the weather station). Second, data from soils and farm polygons were overlaid and then reclassified based on total Al using QGIS. Three Al units were delimited (one per class), hereafter referred to as Al levels. Third, a simple stratified sample was used to select farm polygons randomly until 10 farms were obtained from each Al level, for a total of 30 farms. Each polygon contained one sampling plot. One of these polygons was extracted randomly again at each Al level and was named the reference farm (Table 1, Fig. 1c).
The sampling plots were larger than 30 m × 30 m (in accordance with the pixel resolution of Landsat-8 images) with a set of five points in each sampling area and approximately located in the center of the farm polygon (Fig. 1b). In this study, sampling was carried out prior to sowing only in reference farms (5, 19 and 25), taking approximately 1 kg of soil from the composite sample for laboratory analysis after plants and debris covering the soil surface were removed. The composite sample from topsoil (0-200 mm) was collected using a portable sampler (Eijkelkamp™, Morrisville, USA), taking subsamples and mixing collections from five points (corners and center) of each sampling plot. The location of the center point was recorded.

Soil analysis
Soil samples were air-dried and passed through a 2 mm mesh sieve. The pH was measured in a soil-water mixture (ratio 1:2.5), the available cations of Ca, K, Mg and Na were extracted with a 1 M solution of ammonium acetate, while available P was extracted by the Olsen method. The electrical conductivity (EC) and sulfur (S) from a saturated paste extract and the cation exchange capacity (CEC) were determined by extraction with ammonium acetate at pH 7. Other total concentrations of Al, B, Cu, Fe, Mn and Zn were extracted with ethylenediaminetetraacetic acid (EDTA), and total N was determined by the Kjeldahl method. These methods were applied as described by Olsen et al. (1954) and MAPA (1994).

Monitoring crop phenology
The canopy reflectance pattern is a relationship of reflectance/absorption of specific wavelengths, which is influenced by specific plant traits and has been used for the assessment of biomass and environmental stresses using a vegetation index, e.g., NDVI. In fact, Campos et al. (2018) employed a model based on the relationship of NDVI and biomass. They considered that biomass production is conditioned by stress factors. Additionally, they suggested that under stress conditions, the variability affects crop growth at the field scale. To estimate the biomass production of wheat, in their model, knowledge of growing degreedays (GDD) is required. In this study, GDD was used to evaluate the climatic variation in wheat development due to the different sowing dates that farmers in the study area used. For this purpose, it was calculated according to McMaster and Wilhelm (1997). The phenological development of wheat was visually observed in each growth stage based on the BBCH (Biologische Bundesantalt, Bundessortenamt and Chemische Industrie, Germany) scale proposed by Lancashire et al. (1991). Five main stages were covered: stem elongation, booting, flowering, milk grain and ripening (see Supplementary material S8). This was complemented by calculating NDVI using multi-temporal optical images from Landsat-8 and applying: where NIR and R indicate the near-infrared B5 and red B4 spectral bands, respectively (see Supplementary material S1).
The sampling plots were monitored at least once in the main phenological stages using time series of radar images from Sentinel-1. However, the relationship between backscatter signal radar and NDVI values from optical images was assessed only at the reference farms (5, 19 and 25) to obtain a model.

Landsat-8/OLI optical images preprocessing
Georeferenced multispectral images captured by the operational land imager sensor (OLI) of the Landsat-8 (L8) satellite were used, which are freely downloadable (Table 2). They were acquired from the United States Geological Survey (USGS) server through the Earthexplorer (https:// earth explo rer. usgs. gov/). The L8/OLI images have 9 spectral bands, but for this study, only 7 bands were considered, with a spatial resolution of 30 m (see Supplementary material S1). The temporal resolution was 16 days for all of these.
Satellite optical images were acquired during the 2015 and 2017 crop growing seasons. There were 13 images in 2015, where the cloudiness affecting 10 of them was high (> 70%), rendering them unusable, while it was low in the remaining 3 (< 5%). Thus, to complete the key dates of phenological development, the available images with low cloud cover in the next season were considered valid. Crop rotation was practiced by the farmers in the sampling plots the following year. For this reason, the nearest was the 2017 season, and 4 images were acquired with less than 5% cloud cover (Table 2). This decision was made because, in the study area, the sowing dates varied little when comparing wheat cropping seasons, the crop management practices were the same, and the farmer was not in the habit of successfully neutralizing Al. In L8 imagery, radiometric, atmospheric, topographical and emissivity corrections were made. For this purpose, the iCOR plugin was installed in the sentinel application platform (SNAP) software (De Keukelaere et al., 2018). The plugin (https:// blog. vito. be/ remot esens ing/ icor_ avail able) and software (https:// step. esa. int/ main/ downl oad/ snap-downl oad) can be freely downloaded. After using the plugin, corrected images were obtained in UTM/WGS84 projection with bottom of the atmosphere (BOA) reflectance values.

Sentinel-1/SAR-C radar images preprocessing
Sentinel-1 (A/B) provides C-band radar images, single polarization (VV or VH) and interferometric wide-swath (IW) mode. The time series imagery of Sentinel-1A during the 2015 crop growing season was downloaded free from the ESA-Copernicus open access hub (https:// scihub. coper nicus. eu/ dhus/#/ home). These had a spatial resolution of 5 × 20 m and temporal resolution of 12 days, but currently, the repeat cycle is 3 days in the study area. All images were taken at level 1 with the high resolution and ground range detected (GRD) product type (Table 3).
The Sentinel-1Toolbox in SNAP was utilized to preprocess the SAR images, so the process went through the following stages: orbit calibration, thermal noise calibration, radiometric calibration, speckle filtering and geometric correction (Bousbih et al., 2017). The values obtained were the backscattering coefficient (σ°), and these were converted to decibels (dB) using: In the speckle filtering stage, the gamma map 3 × 3 option was used, selecting the range Doppler terrain correction option in the geometric correction based on an external digital elevation model (DEM) of Tenerife Island, which was obtained from the Spain National Geographic Institute (IGN) for free download (http:// centr odede scarg as. cnig. es/ Centr oDesc argas/ index. jsp). It is an image in raster format that presents a resolution of 5 m, which was resampled to 20 m by the bilinear interpolation method using QGIS.

Test design and data recorded
A field campaign was conducted in 2015, and 30 sampling plots were chosen for an equal number of commercial wheat farms (see "Site description and soil sampling criteria" section). The research was performed following a completely random design, so that soil strata (three levels of Al) were considered as treatments, each with 10 replications (sampling plots), to evaluate wheat stress by Al on canopy and crop development. It was conducted in the winter season, with the crop cv Barbilla (a local soft wheat, Triticum aestivum L.) sown from February 7 to 11 under rainfed conditions. Field management followed usual farm procedures, which included variable seeding rates between farms, without soil fertilization (there was a residue of organic manure previously applied to potato in the crop rotation), and no pesticides were used. In addition, uniform lime distribution is the conventional method used by farmers as a soil amendment to correct soil acidity. In this study, wheat phenological stages were visually monitored (see "Monitoring crop phenology" section) by visiting the sampling plots at least once in each stage. Additionally, during the first visit, farmers were surveyed to determine the seeding rate on their farm (in kg), which was then transformed to kg ha −1 by applying the empirical formula to determine the sowing density: where SOD is the sowing density in seeds m −2 , SER is the seeding rate in kg ha −1 and TGW is the thousand-grain weight in g.
Among the observations of crop development to be carried out were days to heading (DTH) when 50% of the spike or panicle of the pod had emerged, days to ripening (DTR) when 80% of the grains were hard (difficult to split with the thumbnail) and plant height at maturity. The time of harvesting was determined when the grain moisture content was approximately 14%, for which periodic measurements were made with a portable moisture tester (Pfeuffer™ HE 50, Kitzingen, Germany). To complete the assessment of changes during plant growth, multitemporal series of satellite images were used (see "Remote sensing data" section). This included the relationship of a single polarization channel and plant height. An overview of this study is given in Fig. 2.
Plant height was measured in each sampling plot at harvest time. This was done in four corners and the center of each sampling plot within a flexible frame of 0.5 m × 0.5 m, and five random tillers per plot were measured for a total of 25. At the same time, 12 spikes were randomly harvested from the reference farms (5, 19 and 25) to finally obtain a total of 60 spikes for each plot. These spikes were threshed by hand and cleaned to obtain the grains.
In situ data were recorded for the following parameters. Plant height (PHT) was measured as the average distance from the soil surface to the tip of the spike on the main tiller, with awns excluded, in cm. Spike length (SL) was measured in the main tiller from its ear base to the tip of the uppermost spikelet (excluding awns), average in cm. Grain per spike (GS) was determined by counting and averaging the grain number among 10 randomly selected spikes. These grains were then weighed in g, and the grain weight per spike (GWS) was calculated. Thousand-grain weight (TGW) was calculated by weighing 100 kernels randomly selected from the total sample. This was converted to TGW (g) by multiplying by 10. Hectoliter weight (HW) was calculated as the quotient between the weight (g) and volume (ml) of a 10 g kernel sample taken randomly and then converted from g ml −1 to kg hl −1 by multiplying by 100.

3
Finally, the form-density factor (FFD) was calculated for a sample of 25 kernels. For this purpose, kernel length (KL) and kernel width or diameter (KD) were measured with a digital caliper. The formula for calculating FFD is (Giura & Saulescu, 1996): where FFD is unitless, TGW is thousand-grain weight in g, and KL and KD are in mm.

Grain analysis
Dry samples of impurity-free grains were ground and then weighed to 10 g in a porcelain capsule using disposable plastic material to avoid possible contamination. The sample was then incinerated in a muffle furnace at 450 ± 25 °C, gradually reaching this temperature by increasing by no more than 50 °C for 30 min, over a 24 h period. Subsequently, the resulting ash was subjected to the established procedure (MAPA, 1994): proteins were determined by the Kjeldahl method, N was determined by titration with standardized HCl, and macroelement concentrations (Na, P, K, Ca, Mg and S) and microelement concentrations (Cu, Fe, Zn, Mn, Mo, B and Al) were quantified with an ICP-EOS PerkinElmer™ spectrometer.

Statistical analysis
The Shapiro-Wilk test was carried out to confirm the assumptions of normality, randomness and independence of the data. In addition, two-way repeated measures analysis of variance (RM-ANOVA) was used to compare dates (phenological stages) of various plots (reference farms) in search of possible differences. However, to remove the variation that could not be effectively controlled for by test design, analysis of covariance (ANCOVA) (4) FFD = TGW∕(KL*KD) Fig. 2 Flow chart of the procedure for evaluating phenological development and the collection of remote sensing data was applied according to Yang and Juskin (2011). Additionally, Tukey´s test was used to find significant differences between treatments (α = 0.05). Pearson correlation analysis was conducted between the NDVI, plant height and the backscattering coefficient. Then, to develop predictions, a regression analysis was conducted, comparing the performance of these models and the best model used for mapping wheat canopy stress caused by Al levels in soil. Specifically, the root mean square error (RMSE) and the coefficient of determination (R 2 ) were calculated to evaluate the prediction. The statistical calculations were made using Minitab software.

Soil aluminum and nutrient balance
In this study, the concentration of total Al in soil was analyzed. Although the exchangeable Al method is the most commonly used by laboratories, total Al was chosen for this type of research because it is an agile and less cumbersome measurement procedure. As a result, total Al at the 0-200 mm depth was influenced by soil pH, with an inverse relationship. Based on this analysis, total Al varied between 2.86 and 6.7 g kg −1 . Additionally, a summary of other results obtained in the reference farms is shown in Table 4. This is consistent with results from a previous study that reported a relationship between pH and total Al (Hernández & Francisco-Bethencourt, 2017). This relationship clearly shows that at pH < 6.5, the Al level increased due to soil acidity.
The total Al concentration in the soil of the reference farms was classified as low (3 g kg −1 ), medium (3-4 g kg −1 ) and high (> 4 g kg −1 ). It was then compared with the available nutrient content in those soils. A critical reference level was then used to make a balance, evaluating the sufficiency of those nutrients for adequate wheat growth (Goulding et al., 2008) and the Al effect. As a result, an increase in soil Al level decreases Ca (11.25 to 1.38 cmol c kg −1 ) and Mg (4.65 to 0.71 cmol c kg −1 ) concentrations, below the critical level, when compared to a low Al level (3 g kg −1 ) see Supplementary material S3 and S4. This result agrees with that reported by Rahman et al. (2018).

Environmental conditions
During the 2015 season, almost 50% of annual rainfall accumulations occurred throughout the wheat development cycle, between February and August (258.5 mm). All sampling plots received 17% less rainfall than the 14 year average of 313.2 mm. Nonetheless, the accumulated rainfall between February and March (vegetative stage) was 123 mm, which

Wheat development and growth
Winter wheat cv Barbilla displayed similar phenological development in the reference farms and was not affected by sowing dates or soil Al levels; it was sown between February 7 (DOY 38) and 11 (DOY 42) and harvested between August 18 (DOY 231) and 26 (DOY 239). The duration of the total growth cycle (from sowing to harvest) was between 193 and 197 days. Seedling emergence took place seven days after the sowing in winter, followed by the phenophases "flag leaf sheath swollen (booting)", "flowering", and "seed ripening" during spring-summer. Different development stages covered by satellite images are shown in Table 5. Days to heading (DTH) did not show differences regarding the soil Al levels of the reference farms, which ranged from 67 to 70 days after sowing. Physiological maturity or ripening was assumed to be the period when the flag leaf and spikes turned yellow; it was not different between the reference farms (soil Al levels). Days to ripening (DTR) were between 181 and 184 after sowing. Table 7 shows these data.
For monitoring wheat growth, height is an important parameter. In this study, SOD management of wheat growth was observed. SOD ranged from 330 seeds m −2 to 724 seeds m −2 , with a mean of 531.1 ± 87.7 seeds m −2 (± standard deviation), while in PHT, measurements, varied between 1090 and 1320 mm, average of 1191 ± 54 mm. A linear regression analysis showed a highly significant result between the two variables. Nevertheless, the analysis of variance did not reveal significant differences when PHT was compared with different soil Al levels. Therefore, an analysis of covariance (ANCOVA) was used to compare PHT means because the covariate (SOD) affected the data. The ANCOVA results showed a significant (p < 0.05) effect of soil Al levels on PHT (Table 6). Tukey´s test revealed three homogeneous groups, resulting in a mean higher PHT of 1212 mm at the lowest soil Al level (Table 7).

Yield components and grain quality under stress conditions
The effect of soil Al levels on yield components such as SL, GS and GWS was not significant, but it was observed that these were less than local reference values reported by Afonso (2012). However, under high Al stress, TGW was most significantly affected, with a decline to 33.9 g, when compared to the reference value of grain weight (37.6 g) cited by the same author, as shown in Table 8. For other parameters associated with yield, such as FFD and HW, it was observed that they decreased significantly when comparing a low Al level and high Al stress. With FFD, the change was from 2.61 to 1.68, while HW varied from 64.7 to 57.9 kg hl −1 (see Supplementary material S6). In addition, it was observed that an increase in the soil Al level reduced TGW but did not affect KL. This was also reported by Zhang et al. (2018), which is associated with poor grain filling under stress conditions. Studies exploring the stress-induced effects on wheat grain composition and quality are scarce. In this regard, this study found that grain protein and N and S concentrations varied  1 3 as a function of soil Al levels. In protein content, the highest values (13.6% and 17.7%) were observed with increasing Al level and were positively associated with increasing N from 23.8 to 31.1 mg g −1 and S from 0.21 to 0.37 mg g −1 . In contrast, elements such as Al, Mo and B did not show significant variations (Table 9). Other elements analyzed in the grain showed a higher content when compared to a reference value, except Mg, which declined at high Al levels (see Supplementary material S7).

Landsat-8 time series
The measurement of canopy reflectance during the growing season provided information using optical data from L8 imagery. Values from DOY 123 (BBCH 61 and DAS 85) showed a decrease with increasing Al level in the soil, and a significant change was observed in the NIR spectral range (750-1300 nm) at the beginning of flowering (anthesis) stage. Al stress manifests as a change in the spectral reflectance values of the wheat canopy from 0.48 to 0.36 with low and high Al levels, respectively. Instead, significant changes in visible (VIS) at 400-700 nm and shortwave infrared (SWIR) at 1300-2500 nm did not occur, as shown in Fig. 3a.
In addition, as the growing season progressed, significant differences were observed on the dates corroborated by repeated measures analysis of variance (RM-ANOVA) see Supplementary material S10. Therefore, the variation in the observed NDVI occurs in response  to different soil Al levels. As is also evident, according to the crop growth stage (DOY dates), the curves in Fig. 3b show differences. At the beginning of wheat growth (from leaf development to stem elongation), DOY 70 to 91, the differences in the NDVI were clearer. This can be attributed to the soil influence because in early crop development stages, plant size does not allow full ground cover, and most surfaces in reference farms are bare (see Supplementary material S8a). Additionally, as for the differences found at the end of the crop cycle (after flowering), as of DOY 123, can be attributed to differences in the leaf architecture, suggesting a variation in their internal structure due to the Al effect. Instead, differences in NDVI values between DOY 91 and 123, corresponding to stages ranging from stem elongation to flowering, are sufficiently relevant to indicate early effects of canopy stress by Al (Fig. 3b). In this regard, it was observed that on DOY 91, the NDVI was 0.76 at a low Al level and 0.62 at a high Al level.

Sentinel-1 time series
Data provided by Sentinel-1 C-band SAR from the backscatter time series correlate significantly (r -0.84) in VV polarization (σ°V V ) with plant height, which can be estimated by a simple model (see Supplementary material S9). For the visual observation of phenological stages, σ°V V values varied significantly during the crop growing season (DOY dates), as shown in Supplementary material S10 from the RM-ANOVA results. Backscatter variation was particularly complex to analyze; nevertheless, it was observed that these values generally fluctuated between − 1.50 and − 11.94 dB, with differences also observed between soil Al levels (Fig. 4).
During the leaf development and stem elongation stages (DOY 68 to 92), σ°V V increased, which could be explained by an increase in soil moisture (bare soil affects backscatter). Instead, from DOY 92 to 124, σ°V V decreased, possibly influenced by doublebounce between wheat stems and soil at the advanced stem elongation stage. This was slight during the booting stage and stronger at the heading and flowering stages, where the flag leaf opened and inflorescences emerged. This is all associated with an increase in canopy biomass leading to attenuation of the backscatter signal. Later, from DOY 124 to 171, σ°V V increased again during the flowering stage until the onset of the milk grain caused by a canopy with a lower moisture content. The σ°V V decreased slightly in the advanced phase of milk grain stage at DOY 176. It was assumed to be due to grains having a low moisture content. These changes were greater during ripening, e.g., at DOY 184 (Fig. 4). The above results, as described, are associated with stress-free conditions for the crop (low Al level, < 3 g kg −1 ) and correspond with those reported by Harfenmeister et al. (2019), with σ°V V values varying between − 4.47 and − 9.93 dB. However, distortion in the backscattering signal was observed at medium and high levels of soil Al at the stem elongation (DOY 92), booting (DOY 104) and milk-grain (DOY 171 to 176) stages, which has not been reported in previous studies and is assumed to be due to stress caused by structural changes in the canopy due to Al.

Estimation models and predicting the NDVI
PHT and σ°V V were evaluated at three growth stages: booting (DOY 104), flowering (DOY 124) and milk-grain (DOY 176). Regression analysis was applied using linear maximum likelihood (ML) according to linear, quadratic and cubic models. At the beginning of the flowering stage, PHT was higher. This result was similar to that obtained by Bousbih et al. (2017) when they studied plant stress due to soil moisture changes. A linear model was highly significant (p < 0.01) and presented a good fit (R 2 0.69), as shown in Table 10.
A multitemporal analysis was conducted for SAR backscattering (VV and VH) over reference farms with corresponding NDVI values calculated from L8 imagery for both non stressed (low Al level) and stressed (medium and high Al level) canopies during the growing season. The results revealed that σ°V V was positively associated with the NDVI. Curve fitting for the regression function of the data was acceptable for low Al levels (R 2 0.84) and a lower accuracy for both medium and high Al levels in the soil (R 2 0.35 and 0.12, respectively), as shown in Fig. 5a.
In this study, a good correlation was not obtained at Al levels > 3 g kg −1 , as shown in Fig. 5a. It was assumed that this occurred because the dielectric constant was affected by the presence of Al in the soil. In this case, σ°V V is sensitive to the contrast between the dielectric constant and the mineralogical composition of the soil, according to Sharif (1995). Therefore, in soils with higher Al content, the dielectric constant increases, which results in attenuation of the backscattering, as reported by Chen et al. (2007) and Nahm (2010). This may explain the almost flat shape of the NDVI curves in the graph for medium and high Al levels over the wheat crop.
When comparing the NDVI data calculated from multi-temporal optical images of L8 (observed) with the data obtained by applying the proposed regression model (predicted) shown in Table 11, the result gives a good estimate with a significant correlation (Pearson r 0.96, p < 0.05) at low Al levels (< 3 g kg −1 ), while at medium and high levels, it is not significant (r − 0.27 and − 0.17, respectively), as shown in Fig. 5b.

Mapping of within-field variability in the NDVI
Understanding the spatial variability of early-stage stressed wheat plants is important for the subsequent implementation of management strategies. Related to this issue, it is highlighted that in-season Sentinel-1 imagery performs better at the stem elongation stage (DOY 92) when assessing crop development, giving a highly significant correlation between the NDVI (from L8 data) and σ°V V (from Sentinel-1). In addition, a third-order polynomial model generated by the regression analysis was highly significant (p < 0.01) and presented a good fit (R 2 0.70), which could be an interesting option since it can be easily applied (Table 11).
The method developed in this study for NDVI mapping at the field scale employs free and open access data. Previously, optical and radar images were resampled to the same resolution (5 m) by the bilinear method in SNAP, as seen in Fig. 6a. To obtain the model, two deterministic dates were compared, and 2015 was labeled as the predicted NDVI because the model was applied to it. In contrast, the label in 2017 was the observed NDVI because it was calculated from L8 data. To examine the ability to predict wheat canopy stress from soil Al effects, a Pearson correlation analysis was performed between the observed and predicted values with significant results. The descriptive statistics obtained are shown in Supplementary Material S5. The regression model tested presented a low RMSE (0.05) and good accuracy (R 2 0.68) for non-stressed canopies, confirming their usefulness in acidic soils with Al problems (see table in Fig. 6).

3
Finally, the resulting maps of the NDVI in three reference farms show the spatial variability of pixels, where higher values (≥ 0.68) can be classified as non-stressed zones (Fig. 6). This indicates the close link between NDVI and stressed wheat plants due to the effect of soil Al. The effectiveness of this method allows liming within-field management by zone, delineating NDVI pixels with low values (< 0.68). This method is based on predictable data over time being the best method of diagnosis and could be implemented to improve management in precision agriculture.

Growth conditions and stress-induced changes
In this study, a favorable mean temperature from the sowing to flowering stages of 11 °C to 11.9 °C was recorded. These climatic conditions are correlated with high-weight grain at harvest when the daily minimum mean temperature is higher than 6.9 °C (Villegas et al., 2016). Their relationship with the accumulation of GDD (mean 2943.1) was positive, being higher than other GDDs reported, e.g., Salazar-Gutierrez et al. (2013). Therefore, unfavorable climatic conditions were ruled out for research development. Instead, it was observed that wheat is sensitive to the soil Al level with increasing toxicity at pH values below 6.5. PHT showed a negative response in our field-scale study, in agreement with the   Baquy et al. (2017), who, using greenhouse experiments, reported that PHT decreased as the Al concentration increased.
In contrast, Al caused a reduction in TWG and an increase in the protein content of grain, which is in agreement with the finding reported by Flagella et al. (2010) when studying drought stress. In addition, grain N and S concentrations were also increased by Al, and this response has likewise been observed under drought stress (Gooding et al., 2003). Few studies have clearly reported that Al can cause changes in wheat grain composition. Indeed, an adverse change was also observed in the Ca and Mg grain contents because they declined with increasing Al levels in the soil. Other changes associated with Al stress affect yield components such as FFD and HW, with a significant reduction reported in this study when the soil Al level was increased, which has not been reported in the references reviewed.
Spectral reflectance of wheat canopies through the NDVI is a good indicator of plant health. In this study, the crop development stage influenced the correlation between the Al stress level and NDVI, which was consistent with the results obtained by Cattani et al. (2017). In contrast, just as SAR signal variations of Sentinel-1 imagery have been analyzed in previous reports, in this study, the backscattering of VV polarization was significantly correlated with NDVI at the three soil Al levels, and a model was proposed to estimate it (R 2 0.70). Recently, Veloso et al. (2017) found the same relationship when studying crop growth, but obtained a weaker correlation (R 2 0.58).

Spatial variability of the estimated NDVI
Al toxicity symptoms in plants, including stress, are not easily identifiable because they are frequently mixed with other symptoms, e.g., deficiencies in nutrients and drought. Depending on the intensity and severity of this stress, the effects may vary. Recognition of plant stress and its spatial distribution is possible with the regression model reported in this study. Additionally, use of this model could help identify areas of high Al stress at the field scale within seasons. The NDVI mapping results revealed good accuracy of the model and showed that NDVI pixel values higher than 0.68 represented an area with unstressed plants if the evaluation of crop growth at the stem elongation stage was carried out. However, to date, no practical studies relating Al stress and changes in canopy properties are available to compare these results. Instead, there are some reports relating a non-linear NDVI model to wheat canopy biomass, similar to this study; e.g., Mirasi et al. (2019) report a model with NDVI data from a hyperspectral sensor located within a field (precision R 2 0.70).
The value of this research is in the use of management zones; for example, the NDVI derived map can be used to delineate zones with stressed plants where subsurface soil Al level could be a problem and implement variable rate lime application, such that rates can be increased in zones with lower NDVI values and/or decreased (or even not applied) in areas with high NDVI values.

Conclusions
The current trend in remote sensing is to use radar imagery for crop monitoring, which allows data to be obtained under both favorable and adverse weather conditions (e.g., rain, fog). In the present study, a model based on time series data from Sentinel-1 C-band satellite imagery is proposed to assess wheat growth in soils with low to high Al levels. Thus, large areas of cultivated land could be remotely monitored, making the model useful for the diagnosis of Al, since its presence and effects are difficult to detect visually. The results showed patterns of spatial variability of predicted NDVI values as a function of soil Al content, allowing the development of a field-scale map, which could be used for the management of acidic soils through precision agriculture practices. In this regard, the model was more accurate for the early detection of crop canopy Al stress at the stem elongation stage (BBCH 39) using backscattering VV polarization. The methodology developed is generic and could be applied to many more crops or relevant soil properties. In the future, a more complex model would consider interactions between Al effects and other factors, such as soil type, previous crops or fertilization, as well as assess liming efficiency through variable rate application.