Correlation of UAV and satellite-derived vegetation indices with cotton physiological parameters and their use as a tool for scheduling variable rate irrigation in cotton

Current irrigation management zones (IMZs) for variable rate irrigation (VRI) systems are static. They are delineated in the beginning of the season and used thereafter. However, recent research has shown that IMZ boundaries are transient and change with time during the growing season. The primary goal of this study was to explore the potential of using vegetation indices (VIs) developed from unmanned aerial vehicle (UAV) and satellite images to predict cotton physiological parameters that can be used to delineate in-season boundaries of IMZs. A 2 year study was conducted in a 38 ha commercial cotton field in southwestern Georgia, USA. Throughout the two growing seasons, VIs were calculated from UAV, PlanetScope, and Sentinel-2 images. Predawn leaf water potential (LWPPD) and plant height were measured at 37 locations in the field on the same day as the flights and correlated with UAV and satellite based-VIs. GNDVI (Green normalized difference vegetation index) was the best predictor of plant height with correlation values of 0.72 (p < .0001) and 0.84 (p < .0001) for 2019 and 2020, respectively. A secondary goal was to compare the performance of dynamic VRI (DVRI) to conventional irrigation. The field was divided into alternating parallel conventional, and DVRI strips to compare the two scheduling methods. The conventional strips were irrigated using the farmer’s standard method and individual IMZs within the DVRI strips were irrigated based on soil water tension (SWT) measured with a wireless soil moisture sensor network. LWP and SWT measurements correlated well. IMZs were initially delineated using soil texture, apparent soil electrical conductivity (ECa), and yield maps and satellite images from previous years and were modified in-season to reflect patterns observed in the plant height maps. In 2020, the DVRI system prescribed an average irrigation amount of 50.8 mm while conventional irrigation applied an average of 58.4 mm. Average yields for DVRI and conventional were 1248 and 1191 kg ha−1, respectively. The DVRI system resulted in average yield 4.6% higher than conventional irrigation, while applying 14.0% less water. Despite the lower water application by the DRVI system, the performance comparison between the DRVI and the conventional irrigation was not conclusive.


Introduction
Cotton is an economically important fiber crop in the United States and is considered relatively drought tolerant (Chastain et al., 2014). Nevertheless, drought negatively affects seed germination and seedling establishment (Bradow & Bauer, 2010), root system growth (McMichael et al., 2010), plant height and node development (Marani et al., 1985;Wells et al., 2010), and lint yield (Wanjura et al., 2002). These negative effects are observed even in regions where average precipitation during the growing season meets crop demands (Bednarz et al., 2003) because episodic drought caused by irregular distribution of precipitation can lead to slower crop development and lower yield and water use efficiency (WUE) (Bednarz et al., 2003;Chastain et al., 2014). Additionally, crop water requirements can vary within a field due to the spatial variability of soil texture, moisture, and slope (Vellidis et al., 2016a(Vellidis et al., , 2016b. One approach to addressing the spatially variable water needs of crops is dynamic variable rate irrigation (DVRI). DVRI consists of site-specific application of irrigation water to crops based on real-time measurements or estimates of soil moisture conditions within individual IMZs (Liakos et al., 2017). In contrast to static prescriptions which remain constant over time, DVRI prescriptions change during the season. In this context, DVRI is an important practice that aims to optimize WUE in crops by applying the needed amount of water at the appropriate time to individual IMZs. DVRI can help mitigate the effects of the spatial variability in soils and drought stress, increasing yield and maximizing WUE. Liakos et al. (2017) developed and demonstrated a soil moisture sensor-based DVRI system. Their system installed a soil moisture sensor node with sensors at three depths within each IMZ and used the data from the sensors to drive an irrigation scheduling decision support tool (DST) (Liakos et al., 2015;Liang et al., 2016). Thus, each IMZ received a unique prescription based on the soil moisture in that zone. Their system resulted in WUE gains ranging from 15 to 40% and yield gains of up to 4% in peanut when compared to farmer standard methods (Liakos et al., 2017). Currently, the boundaries of IMZs used for VRI systems are static. They are usually delineated when a VRI system is installed and used thereafter, ignoring temporal changes in wetting and drying patterns during the growing season caused by the interaction between soil, plant, and environment (O'Shaughnessy et al., 2015).
Recent research has shown that it may be possible to adjust in-season IMZ boundaries using crop water status data (Cohen et al., 2017;O'Shaughnessy et al., 2015). Leaf water potential (LWP) has been identified as a reliable indicator of plant water status (Argyrokastritis et al., 2015;Bellvert et al., 2016;Paço et al., 2013) and Cohen et al. (2017) used measured LWP in cotton to develop temporally-variable IMZ boundaries. LWP has also been reported to directly affect other cotton physiological parameters that affect yield such as photosynthetic assimilation rate (Turner et al., 1986), stomatal conductance (Lacerda et al., 2022), and plant growth (Chastain et al., 2016a(Chastain et al., , 2016b, and these parameters can in turn be used to indicate field variability in water stress. However, collecting LWP data in large fields at a scale that can characterize spatial variability and be used to implement DVRI is time consuming and has high labor costs. Many authors have identified leaf temperature as a good indicator of plant water status González-Dugo et al., 2006;Conaty et al., 2015). The premise that canopy temperature can indicate the level of water stress in plants assumes that in well-watered plants under specific environmental conditions, transpiration is at its potential rate, thereby cooling the crop. As water becomes limiting, transpiration is decreased due to stomatal closure, causing canopy temperature to increase (Jackson et al., 1988). The crop water stress index (CWSI) can be used to calculate crop water stress from canopy temperature Jackson et al., 1981). Recent advances in thermal imagery allow the collection of canopy temperature for large areas (Cohen et al., 2005). However, the thermal resolution must be within ± 1 ºC for accurate calculation of CWSI (Cohen et al., 2017) which requires expensive thermal infrared (IR) cameras.
The overall goal of this study was to evaluate whether multispectral, remotely sensed images can be used to predict spatial variability of cotton physiological parameters that are indicative of plant water status as an alternative to in-field measurements or thermal imaging. Specific objectives were to 1) explore the potential of using multispectral VIs developed from UAV and satellite images to predict LWP and plant height, 2) adjust in-season boundaries of IMZs based on plant parameters that addresses in-season temporal variability, and 3) compare performance of DVRI to conventional irrigation.

Study area
The study was conducted during the 2019 and 2020 growing seasons in a 38 ha commercial cotton field located in Miller County in southwestern, Georgia, USA (31°11′21″ N, 84°45′42″ W). The region has a humid subtropical climate (Kottek et al., 2006) with an average temperature ranging from 7.8 °C in January to 26.5 °C in July (Knox & Mogil, 2020). Precipitation is well distributed year-round with the driest months (April and October) receiving 85.3 mm and 86.6 mm, respectively, and the wettest months (July and March) receiving 133.9 mm and 122.2 mm, respectively. In Miller County, the 20 year minimum and maximum average temperatures between April and November were 17.3 °C and 29.7 °C, respectively (NOAA, 2021).
Cotton was planted on May 31st in 2019 and on June 2nd in 2020. Cotton rows were planted parallel to the western edge of the field (Fig. 1). The field was irrigated by a center pivot system equipped with FarmScan (Advanced Ag Systems, Dothan, AL, USA) VRI controls. In both years, the field was divided into four approximately 200 m wide parallel Fig. 1 Conventional and VRI strip layout, soil moisture sensor location, delineated IMZs, and 10 m radius circle plots around sensor locations in 2019 and 2020. In 2019, 6 sensor locations and 2 zones (indicated by red numbers) were excluded from data collection due to time and personnel limitations 1 3 strips with two strips irrigated uniformly (conventional) and two irrigated using the DVRI system described earlier (Fig. 1). The conventional strips were irrigated using the grower's standard irrigation method which was to apply 15 mm at each irrigation event on a weekly basis early in the season and twice weekly during the reproductive stage unless precipitation of that amount or greater was received. The strips that were irrigated using DVRI in 2019 were irrigated uniformly in 2020, and the strips irrigated uniformly in 2019 were irrigated using DVRI in 2020 (Fig. 1).

IMZ delineation
Prior to the growing season, DVRI strips were delineated into 25 IMZs in 2019 and 29 in 2020 based on five data layers that included apparent soil electrical conductivity (EC a ), elevation, NDVI and yield maps from previous years, and soil texture data (Fig. 2). EC a was measured at 30 cm (shallow) and 90 cm (deep) using the Veris 3100 instrument (Veris, Salinas, KS, USA). Because deep EC a measurements integrate the entire root zone, deep EC a data were used to create an EC a map of the whole field and used as a layer for IMZ delineation in both years. A digital elevation model (DEM) of the field was developed from elevation data collected during ECa mapping with RTK GPS with approximately 20 mm x, y, and z accuracy. In 2019, 75 soil core samples were collected from the field and analyzed. Each core sample was divided into one surficial layer of 7.6 cm and three additional layers of approximately 15 cm reaching a total depth of 53 cm. Percent clay content from 38 to 53 cm layer was used to create a percent clay content map of the field that was used in both years to delineate IMZs. The deepest layer was selected because the majority of cotton roots can be found between 30 and 90 cm (Ritchie et al., 2007). Kriging was used to interpolate between sampling locations. Yield and NDVI maps from 2017 were used as data layers in 2019, while in 2020, yield and NDVI maps from 2019 were used. PlanetScope CubeSat (Planet Labs, San Francisco, CA, USA) satellite images from July 2017 and July 2019 were used to develop the NDVI maps. The IMZ delineation process was done manually using ArcGIS (ESRI, Redlands, CA, USA) after visual analysis of all data layers to identify areas of variability. Each data layer was classified in 5 different classes using quantile classification, in exception to the clay % map that was classified in 7 groups. The zones resulted from the classification in each map together with the researcher's expertise and knowledge of the field were used to guide the IMZs boundary delineation. Boundaries were manually delineated around areas that displayed uniform yield, clay content, ECa, elevation and NDVI as much as possible. The goal of this process was to maximize the uniformity of the measured parameters within each IMZ. In 2020, the boundaries of the initial IMZ map were modified mid-season based on plant height maps generated before the flowering stage.

Soil moisture measurement and irrigation scheduling
UGA Smart Sensor Array (UGA SSA) nodes were installed in 37 different locations in the field (Fig. 3). The UGA SSA is a smart, wireless soil moisture sensing system that measures soil moisture at three depths that vary depending on the crop (Vellidis et al., 2013). A UGA SSA node consists of the probe used to measure soil moisture and the electronics package used to transmit the data to a cloud server. In this study, sensors on the soil moisture probe were at depths of 0.15, 0.30 and 0.45 m when installed. Soil moisture was measured in terms of soil water tension (SWT) measured in units of kPa, which is an indicator of the pressure/suction a plant must exert to extract water from the soil matrix. The UGA SSA is capable of sensing soil moisture from 0 to 200 kPa. Each IMZ had at least one UGA SSA node, with a total of 31 nodes installed in the DVRI strips. For monitoring purposes, three UGA SSA nodes were installed in each of the two conventional strips. Nodes were installed within the row between plants about two weeks after emergence. Nodes were removed from the field a few days before harvest. SWT data were collected hourly for the duration of the growing season.
Irrigation thresholds of 50 kPa prior to first flower and 40 kPa after first flower were used to initiate irrigation (Chastain et al., 2016a(Chastain et al., , 2016bMeeks et al., 2017). The daily weighted average SWT of the three sensors at each node at 07:00 AM was used to determine if an IMZ required irrigation. The weighting factors used to calculate the weighted average SWT characterize the estimated distribution of roots at the three sensor depths and changed as the rooting depth increased over time. Equation 1 shows the weighting factors used when cotton is flowering. Table 1 shows how the weighting factors changed during the growing season. The DVRI DST was used to determine the amount of irrigation water needed to bring the soil profile in that IMZ back to the desired soil moisture condition using the weighted average SWT. In this study, the desired soil moisture condition was 80% of field capacity or approximately 12-15 kPa SWT which leaves "room" for the soil profile to absorb precipitation events.
Because the grower typically irrigated once or twice per week, irrigation applications in the DVRI strips mostly occurred when the grower ran the pivot to apply water to the conventional strips. In these instances, water was applied to IMZs where SWT was within 10 kPa below the threshold to optimize the logistics of operating the pivot daily. On one or two occasions per season, the pivot was run because SWT exceeded the set threshold in one or more IMZs and irrigation was applied to the needed zones and to the conventional strips.

Field measurements
Physiological measurements of cotton plants were carried out weekly from July until the last week of irrigation in late September in both years at the location of the UGA SSA nodes. There were 9 sampling days in 2019 and 10 in 2020. Predawn LWP (LWP PD ) was collected using a model 615 Scholander pressure chamber (PMS Instruments, Albany, OR) from two plants at 31 node locations in 2019, and 37 node locations in 2020. Measurements were taken between 3:00 and 6:00 AM. Specifically, leaves from the fourth node below the terminal of the plant were severed near the base of the petiole from two plants at each sample location. The petiole was then placed in a compression gasket and the leaf blade was sealed in the chamber and immediately pressurized at a rate of 0.05 MPa s −1 . The pressure required for xylem sap to appear at the cut surface of the petiole was converted to LWP PD in negative MPa. Selected plants were adjacent to the UGA SSA node. The dates of field measurements and the number of each measurement made are shown in Table 2.
Plant height measured from the soil surface to the apical meristem was collected from five plants adjacent to the UGA SSA nodes. The same five plants were used throughout the growing season. A LI-6800 portable photosynthesis system (LI-COR, Lincoln, NE, USA) was used to collect data on midday photosynthetic assimilation rate (A n ) and stomatal (1) conductance to H 2 O (g s ). Leaf chamber settings included a flow rate of 600 µmol s −1 , reference CO 2 = 400 µmol mol −1 , air temperature = ambient temperature, relative humidity = 60 ± 15%, and chamber light intensity = 1500 µmol m −2 s −1 photosynthetically active radiation (PAR). The leaf was clamped into the chamber until steady-state conditions were reached (60 to 120 s per sample). Midday measurements were performed from 11:30 AM to 2:30 PM on the fourth mainstem node below the terminal for two plants adjacent to the sensor node. Because of the time needed to make these measurements, they were limited to the two southern strips (Fig. 1). Measurements were made at 18 sensor node locations in 2019 and at 15 sensor node locations in 2020. Yield data were collected at the time of harvest in the whole field with a cotton picker equipped with a yield monitor. Yield data was cleaned and processed using ArcGIS (ESRI, Redlands, CA, USA) to create yield maps for both years.

UAV and satellite imagery acquisition and processing
Remote sensing data were collected from a satellite platform and with a UAV. Satellite images acquired from up to 3 days before and after field data were collected were downloaded from two different satellite system platforms. The additional days were included to ensure the availability of a cloud-free image. Images were used to calculate a variety of VIs that are described below. The PlanetScope satellite constellation is one of the three constellations operated by Planet Labs (San Francisco, CA, USA). It currently has approximately 130 satellites imaging the earth's land surface with a spatial resolution of 3 m and a temporal resolution of one day (Planet labs, 2020). The second satellite platform used was Sentinel-2 operated by the European Space Agency (ESA). Sentinel-2 is a two-satellite platform operating in the same orbit 180° apart with a revisit time of 5 days and spatial resolution of 10, 20, and 60 m depending on spectral band. PlanetScope images provide data in the blue (455-515 nm), green (500-590 nm), red (590-670 nm) and NIR (780-860 nm) regions of the spectrum. Surface reflectance data were downloaded directly from the Planet Labs website (https:// www. planet. com) and a factor (10,000) was applied to obtain correct reflectance data (Planet Labs, 2019). Sentinel-2 images were downloaded from the United States Geological Survey (USGS) Earth Explorer website (https:// earth explo rer. usgs. gov/) and bands 2 (490 nm), 3 (560 nm), 4 (665 nm) and 8 (842 nm) at 10 m spatial resolution, and 11 (1610 nm), and 12 (2190 nm) at 20 m resolution were used to calculate VIs. The Semi-Automatic Classification Plugin (SCP) in QGIS 3.14.0 (QGIS Development Team, 2020) was used to perform atmospheric correction on all Sentinel-2 images.
The UAV was a 3DR Solo quadcopter (3D Robotics, Berkeley, CA, USA) equipped with a Parrot Sequoia (Parrot Drone SAS, Paris, France) camera that acquires images in four spectral bands: Green (530-570 nm), Red (640-680 nm), NIR (770-810 nm) and Red Edge (730-740 nm). UAV flights were performed at an altitude of 120 m with 75% front and side overlap resulting in a resolution of approximately 125 mm per pixel. UAV image processing was performed using Pix4Dmapper software (Pix4D SA, Lausanne, Switzerland) version 4.4.12. Single images were stitched to create reflectance map mosaics of all four bands. During the stitching process, geographic correction was performed using six ground control points with known coordinates distributed around the perimeter of the field. Calibration panel images taken after each flight were used to perform radiometric calibration of the final reflectance maps.
ArcGIS (ESRI, Redlands, CA, USA) software was used for remote sensing data extraction. Reflectance maps created from the UAV and satellite images were used to calculate the VIs shown in Table 3 that were selected based on their ability to quantify specific plant growth parameters and the spectral bands available. Average VI values were extracted from the 10 m radius area around each soil moisture sensor node and from the delineated IMZs. A 5 m buffer was created between zones, and images from all 3 platforms were overlapped to minimize potential data extraction errors due to geographic accuracy.

Data analysis
Pearson correlation at 95% confidence interval was used to assess season-long and weekly correlations between VIs and plant height, and season-long correlations between VIs and LWP PD . The VI with the highest correlation with plant height was selected. An exponential function with 2 parameters was used to model the selected VI response to changes in plant height during the entire growing season (Eq. 2).
where, y represents the plant height, a is growth rate, b is the scale parameter, and x represents the GNDVI value. A linear regression model was then developed to model the relationship between VI and plant height at a specific period during the season and used to create a predicted plant height map. Curve fitting and statistical analysis was performed using JMP® PRO 16.0.0 (SAS, Cary, NC).

Meteorological data
The two growing seasons presented sharply differing weather conditions (Fig. 4). Precipitation during the 2019 growing season was 485 mm and skewed towards the beginning of the growing season. Precipitation was higher (699 mm) and more evenly distributed during the 2020 growing season. In 2020, the field received 180 mm of precipitation from planting to Table 3 Vegetation indices calculated from UAV, PlanetScope and Sentinel-2 images SWIR 11: Sentinel-2 band 11; SWIR 12: Sentinel-2 band 12 CIgreen green chlorophyll index, GNDVI green normalized difference vegetation index, GRVI green-red vegetation index, NDVI normalized difference vegetation index, OSAVI optimized soil-adjusted vegetation index, NLI non-linear index, GBNI: green-blue normalized index, NDWI normalized difference water index, NMDI normalized multi-band drought index Cotton plant sensitivity to water stress starts to increase during squaring (Wrona et al., 1999), reaching its maximum during the flowering stage (Snowden et al., 2014). Although the cotton was primarily in the boll-filling stage during the 2 months without precipitation during the 2019 growing season, the plants were still highly sensitive to water stress. Minimum and maximum temperatures were similar during the first 3 months of the season in both years. Maximum temperatures from June to September in 2019 were 33.6 °C with the highest average temperatures in September (34.4 °C). During 2020, maximum temperature for the same period averaged 32.3 °C with the warmest temperatures observed in August. October was the coolest month in both years with minimum and maximum temperatures of 16.9 and 29.4 °C in 2019 and 12.1 and 24.7 °C in 2020. Figure 5 shows the season-long SWT, LWP PD, A n , and g s distributions among the IMZs delineated in 2019 and 2020. Data analysis was divided into before canopy closure (BCC) and after canopy closure (ACC) due to saturation of the NIR-based indices at full canopy (Hunsaker et al., 2003) and the differences in the cotton plant physiological responses between the vegetative and reproductive stages (Cohen et al., 2017). In 2019, data for zones 2 and 9 were not collected and are not presented. In 2020, data at the ACC stage refers to IMZ map 2, when IMZ boundaries was adjusted based on predicted crop height map

growing season
In 2019, five field-measurement days took place during the BCC period and four during the ACC period (Fig. 5a). Canopy closure was achieved approximately 100 days after planting (DAP) between the last week of August and first week of September. All four measurements exhibited within-season and within-zone temporal variability as well as between-zone spatial variability. SWT was significantly higher ACC indicating drier soil conditions than BCC. During early plant development, the soil profile in IMZs 1,4,6,10,12,13,15,16,17,20,21 and 25 maintained adequate soil moisture from precipitation alone with SWT values averaging 14 kPa. As indicated earlier, the pre-bloom irrigation scheduling threshold was 50 kPa. Zone 23 was driest with an average SWT of approximately 122 kPa despite irrigation. The highest within-zone variability in soil moisture was observed in zones 3, 5, 8, 18, 23 and 24. The SWT in these zones varied by an average of 90 kPa suggesting that for one or more days during pre-bloom, rainfall, and irrigation together were not sufficient to keep the soil profile within the desirable SWT range. This is indicative of sandier soils where SWT can increase sharply because of the soil's low water holding capacity. SWT greater than 70 kPa at the pre-bloom stage may result in slower plant growth and significant yield loss (Collins et al., 2011;Meeks et al., 2017). The remaining IMZs exhibited moderate variability with SWT values all below 75 kPa, indicating that in these zones soil moisture was adequate to ensure good growth.
With the onset of the reproductive stage, cotton slows its vegetative growth but increases its water requirements. As indicated earlier, the irrigation threshold during the reproductive stage was 40 kPa. Increased water demand coupled with decreased rainfall from September onwards should have resulted in more irrigation. Instead, the grower who managed this field decreased the frequency of both conventional and DVRI irrigation events because of disease pressure caused by the frequent rains earlier in the season. As a result, during the latter half of the reproductive period, soils were drier and higher SWT was recorded in most of the field. Of the 25 IMZs, only zones 4 and 25 had average SWT below or close to 40 kPa ACC. There was high soil moisture variability between zones at this stage, and relatively high within-zone temporal variability in at least 14 IMZs with SWT values varying from 40 to 122 kPa. Spatial and temporal variability was also seen in the plant's physiological indicators. LWP PD showed negative correlation to SWT, where higher SWT values were typically associated with lower LWP PD . Prior to bloom, LWP PD in all zones varied between 0 and -0.6 MPa which is above the threshold (− 0.8 MPa) for mainstem growth inhibition (Chastain et al., 2016a(Chastain et al., , 2016b. ACC, LWP PD values decreased significantly below − 0.6 MPa, displaying varied levels of stress. The same trends were observed in stomatal conductance and photosynthetic assimilation rate. The water stress observed after bloom caused a decrease in g s that in turn affected A n . The largest difference in average assimilation rate was observed between zones 20 and 23 with values of 30.8 and 14.5 µmol m −2 s −1 , respectively. Turner et al. (1986) observed a decrease in photosynthetic rates of 20 µmol m −2 s −1 when LWP decreased from − 0.4 to− 1.6 MPa. Furthermore, it is generally accepted that reductions in stomatal conductance are the primary drivers of drought-induced photosynthetic decline in cotton (Chastain et al., 2014(Chastain et al., , 2016a(Chastain et al., , 2016b.

growing season
In 2020, new IMZ boundaries were delineated at the end of the vegetative stage based on a predicted plant height map (Fig. 6). The quantile classification method was used to classify predicted plant height in different classes. Boundary of IMZs delineated at the beginning of the season were then adjusted to match plant growth patterns observed when peak growth was achieved before the onset of reproductive stage. The new map was used for irrigation management for the remaining of the season. Hence, the data distribution at the ACC stage shown in Fig. 5b refers to IMZs in map 2.
Soil moisture was consistently higher during the 2020 growing season than the 2019 growing season (Fig. 5b). This was especially true for ACC. In contrast to 2019, SWT between-and within-zone variability was higher BCC with zones 1, 3, 6, 7, 9, 14, 17, 22, and 24 showing values at field capacity or below with little change over time, and zones 23, 27 and 28 showing a greater range of values from wet soil (< 33 kPa) to very dry soil (> 100 kPa). During September and October 2020, there were frequent precipitation events associated with tropical storms which led to very low SWT values in all IMZs (Fig. 5b). This was reflected in the physiological measurements. LWP, g s and A n were similar to 2019 BCC but much higher ACC. Although the measurements were higher ACC in 2020 than in 2019, in some areas of the field, the plants were damaged from the high winds associated with the tropical storms and this likely depressed the photosynthetic rate.

Relationships among yield, plant height, SWT and physiological parameters
Regression analyses were used to evaluate the relationship between soil moisture and plant water status and how the water status related to leaf stomatal conductance (g s ) and assimilation rate (A n ). Linear regression and second-order polynomial regressions were Fig. 6 IMZs boundary mid-season adjustment based on predicted height map before the onset of reproductive stage. Maps 1 and 2 represent initial IMZs and adjusted IMZs, respectively developed between SWT, LWP PD , A n , and g s (Fig. 7). Average values of the measurements taken at each sensor node at each sampling date were used in all regressions. Measurements were taken at 31 nodes on 9 dates in 2019 and 37 nodes on 10 dates in 2020. A strong linear relationship (R 2 = 0.7, p < 0.0001) between LWP PD and SWT was observed in 2019 (Fig. 7a). A strong non-linear relationship was observed between LWP PD and A n (R 2 = 0.62, p < 0.0001) (Fig. 7b) and g s (R 2 = 0.60, p < 0.0001) (Fig. 7c). In 2020, relationships between SWT and LWP PD (Fig. 7d), LWP PD and A n (Fig. 7e), and LWP PD and g s (Fig. 7f) were all weak with coefficients of determination of 0.05, 0.18 and 0.12, respectively. This is primarily because LWP PD measurements were consistently high in 2020 because the soil moisture was also consistently high and did not vary as much as in 2019. The data point cluster observed in the regression between SWT and LWP PD when SWT values were low (Fig. 7d) may be related to that the soil moisture sensors used in the UGA SSA probes are less accurate at near-saturation conditions.
In 2019, the relationship between LWP PD and the other two physiological parameters indicated a limiting effect of varying degrees of water stress in the cotton photosynthesis, in which both g s and A n tended to decrease as LWP PD become more negative. Most data points furthest from the curve in this year were from plots 22 and 23 that were located along the eastern edge of the field where there is extensive sand deposition from erosion ( Fig. 1) and as a result plant growth was limited leading to lower g s and A n values regardless of water status. On the other hand, the absence of water stress in 2020 indicated that the variability in g s and A n were caused by factors other than water availability.
Higher average plant growth was generally observed in IMZs that received a higher total amount of water during the growing season except for IMZs located in eroded areas of the field (Fig. 1). In these areas, plant growth was limited irrespective of the amount of water received. In 2019, IMZs that received less than 700 mm of precipitation and irrigation presented an average final height of 0.9 m, while IMZs that received more than 700 mm of water had an average final height of 1.12 m (Table 4). In 2020, all IMZs received a total water amount higher than 700 mm. IMZs with average total water between 700 and 750 had an average height of 1.06 m, and IMZs that received more than 750 mm had an average height of 1.19 m.
Final plant height averages from each IMZ were plotted against final yield (Fig. 8). In 2019, a strong linear relationship between height and yield was observed (R 2 = 0.64, p < 0.0001) with taller plants achieving higher yields. A weaker relationship was observed in 2020 (R 2 = 0.54, p < 0.0001) due to a decreased yield observed in plants taller than 1 m. The decrease in yield can be attributed to the damage caused during a series of tropical storms that passed over the field. Taller plants suffered from strong winds and were blown over and lodging was observed in these areas of the field. As a result, these plants were likely not harvested as efficiently. High soil moisture may also have resulted in additional vegetative growth at the expense of reproductive growth, as has been observed previously under conditions of water excess (Ermanis et al., 2021). In spite of weaker relationship in 2020, results from this study show similarities with those of Sui et al. (2013) in which they observed higher plant heights and yield in plots with more available water a (irrigated) when compared to less available water (rainfed).

Plant height prediction
Pearson's correlation was used to compare the response of VIs developed from UAV imagery to plant height for 2019 and 2020 (Table 5). The overall correlation of the entire season indicated that all VIs showed high positive correlation with plant height during both years. GNDVI had the strongest overall correlation with r values of 0.72 (p < 0.0001) and 0.84 (p < 0.0001) for 2019 and 2020, respectively, followed by CIgreen with r values of 0.72 (p < 0.0001) for 2019, and 0.82 (p < 0.0001) for 2020. Within-season correlations from 60 to 77 days after planting (DAP) in 2019, and from 64 to 71 DAP were the highest with all values above 0.7 (p < 0.0001). VIs developed from PlanetScope images showed similar correlations to plant height as VIs from the UAV. However, due to their higher spatial resolution, UAV results were used for further analysis. An exponential model was used to further quantify this relationship between GNDVI and plant height (Fig. 9). Plant height response curve to GNDVI had higher accuracy in 2020 with an R 2 value of 0.71 and RMSE of 15.05. In 2019, the weaker relationship between the two variables was most likely caused by the higher number of dates collected at the end of the season, compared to early season, when increased assimilate demands by the fruit causes vegetative growth to significantly decrease or to cease completely (Wells et al., 2010). Therefore, changes in the VI values at this stage are not associated with plant growth.
In 2020, UAV data collection started at an earlier stage than in the previous year. At 43 to 64 DAP the correlation coefficient between GNDVI and plant height was below 0.7 (Table 5), and points were more spread around the curve (Fig. 9). While in 2019, correlation between the two variables was already high at the first day of data collection at 60 DAP. These results indicate that plant height prediction accuracy using GNDVI is higher after around 60 DAP, and image collection at an earlier stage might not be necessary. In both years, it is possible to observe that when the crop reached a plant height of approximately 120 cm, which occurred between 85 and 90 DAP (Table 5), GNDVI values remained constant between 0.77 and 0.85 (Fig. 9). At this time, the point of saturation was reached and the VI was not able to sense the small changes in height. The exponential The period between the dates in which the correlation coefficient increased (at 60 to 64 DAP) and later started to decrease (at 82 to 85 DAP) represents the period in which UAVbased multispectral flights are the most suitable for plant height prediction. Hence, based on results from both years, it is possible to infer that the optimal period for the development of height prediction models is between 60 and 80 DAP. Similar results were observed by Raper  and Varco (2015). NDVI and GNDVI showed a strong relationship with plant height during the period from bud formation to first flower, which commonly occurs around 60 to 70 days after planting (Robertson et al., 2007). Based on the results from 2019, in 2020 the UAV flight performed at 64 DAP was used to calculate GNDVI, and a linear model was fitted between GNDVI and measured plant height at the same date resulting in a R 2 of 0.72 (p < 0.0001) and a root mean square error (RMSE) of 11.2 cm (Fig. 10). This model was then used to develop a predictive plant height map using the 2020 GNDVI data for in-season dynamic IMZ delineation. The predicted map (Fig. 11a) shows great similarity to the final yield map (Fig. 11b).

Relationship between plant height and yield
A strong linear relationship (R 2 = 0.55, p < 0.0001) was observed between zone plant height and zone average yield (Fig. 12). This relationship was even stronger showing an R 2 of 0.75 when zones 15, 19 and 23, that did not follow the same relationship, were excluded from the regression. Zones 19 and 23 showed higher yield compared to other zones despite their lower plant height average, while zone 15 had one of the lowest recorded yields among IMZs despite its average predicted height ( Table 8 in Appendix).
Based on the strong relationship between plant height and final yield, the use of predicted plant height maps during the season can give great insight on high and low yielding areas of a field. In addition, the effect of water availability in plant growth suggests that predicted plant height maps can serve as an indicator of plant status at the beginning of the season during the vegetative stage and enable the delineation of new IMZs that address temporal changes caused by the interaction of plant, soil, and environment during that period.

Relationship between VIs and LWP PD
Overall VI values calculated from UAV, PlanetScope and Sentinel-2 had good to moderate correlations with LWP PD measurements in both years (Table 6). The strongest correlations were observed in 2019 for CIgreen and GNDVI (-0.63 and -0.61, p < 0.0001, respectively) calculated from Sentinel-2 images. In 2020 only OSAVI (-0.49, p < 0.0001) calculated from PlanetScope and GBNI (0.42, p < 0.0001) calculated from Sentinel-2 showed moderate correlation with LWP PD. In contrast to the results from this study, Beeri et al. (2018) found that GBNI, NDWI and NMDI had a strong relationship with LWP (R 2 = 0.71, 0.71, 0.53 p < 0.0001, respectively) in cotton fields in Israel and Australia. The higher relationships in the Beeri et al. (2018) study can be attributed to the fact that LWP values were collected from the wettest and driest parts of the field, while in the current study, LWP data from the whole field were used independently of the level of stress. The consistently high LWP PD values clustered above -0.6 MPa (Fig. 7d) were likely the cause of low correlations in 2020. Additional results from Beeri et al. (2018) reported strong relationship between LWP and vigor VIs such as NDVI (R 2 = 0.69), which corroborates present results. The correlations observed in this study were not strong enough to develop LWP prediction models with high accuracy. This is likely because the range of measured LWP PD was small compared to studies conducted in semiarid or arid environments.

Comparison between DVRI and conventional irrigation scheduling
The precipitation difference between the two growing seasons can be illustrated by the difference in irrigation applied and number of irrigation events (Table 7). A total of 7 irrigation events occurred in 2019 to compliment the 485 mm of precipitation. Irrigation was triggered mainly during September when precipitation was low. The DVRI strips received an average of 86 mm of water from irrigation and an average of 87 mm was prescribed to the uniform strips. As described earlier, in 2019, the crop experienced water stress in September and October as the grower reduced the frequency of irrigation because of disease pressure (see Fig. 4). On several occasions when DVRI IMZs needed irrigation, the grower was not willing to apply water because of the disease concerns. This resulted in SWT values much higher than the established thresholds.
Although there were no significant differences in the overall water applied across the field, the amount of water applied to individual IMZs differed at every irrigation event and was determined by the dynamic DVRI DST. The average yield for the conventional strips was 1890 kg ha −1 while yield of the DVRI strips averaged 1794 kg ha −1 which is a 5.21% difference in yield. Irrigation water use efficiency (IWUE) for the DVRI strips was 20.7 kg/mm and 21.7 kg/mm for the conventional. Although, the conditions in this study led to very little difference in water applied between conventional and DRVI systems, similar DVRI studies conducted in southern Georgia on peanut resulted in yield increases of up to 4% and IWUE gains of up to 40% (Liakos et al., 2017). The lower DVRI yields may have been affected by the IMZs that included highly eroded areas of the field as well as areas with sand deposition and where yields in 2019 were very low (Figs. 1, 2 and 11). There was little difference in IWUE because overall, the same amount of water was applied to both treatments. Because there were regular precipitation events throughout the 2020 growing season, the field was irrigated only four times, three of which took place in October. The DVRI system prescribed an average irrigation amount of 51 mm while the conventional irrigation applied an average of 58 mm. Average yield for DVRI and conventional was 1,248 and 1,191 kg ha −1 , respectively. In 2020, the DVRI system resulted in an average yield 4.6% higher than conventional irrigation, while applying 14.0% less water. IWUE was 24.6 kg/mm in the DVRI strips and 20.4 kg/mm in the conventional strips. As in 2019, the average yields may be affected by the highly eroded areas in the field.
The physiological parameters comparison between DVRI and conventional irrigation is shown in Fig. 13. Despite the small differences in irrigation between the DVRI and conventional treatments, the average LWP PD in all IMZs tended to be equal or higher than the LWP PD collected in the 6 locations in the conventional strips in almost all date of both years. This difference was more prominent in the ACC (after canopy closure) period of 2019 because of the low precipitation. Although the DVRI irrigation prescriptions were not correctly followed by the grower in that year, the few irrigation events based on the SWT thresholds in September were enough to guarantee a higher LWP PD in the DVRI strips. When comparing both years, the overall season-long values were much lower in 2019 than in 2020, as a result of the difference in total water (irrigation + precipitation) received in the field in both years. The A n and g s measured in both treatment strips tended to follow the same trends as the LWP PD, with higher values for the DVRI treatment when compared to conventional in almost all dates. As a result of the more frequent and higher precipitation in 2020, the differences in A n and g s are smaller than in 2019.

Fig. 13
Physiological parameters comparison between DVRI and conventional treatments, including average LWP PD (a-d), A n (b-e), and g s (c-f) in 2019 and 2020 respectively. DVRI averages considered 23 IMZs in 2019, and all 29 IMZs in 2020, while average for the conventional irrigation treatment were calculated from the 6 measurements collected in each year

Conclusions
Results from this study showed that management zones changed spatially and temporarily (within and between growing seasons). This temporal variation within a growing season caused by the interaction of soil, plant and environment can be addressed by monitoring plant growth patterns and physiological responses. Based on the results of this study, remotely sensed data in the visible and NIR regions of the spectrum can be used in the form of VIs to estimate plant height in cotton with high accuracy between 60 and 80 days after planting when cotton plants are close to the flowering stage. The ability to predict plant height for the whole field at high spatial resolution facilitates the identification of cotton growing patterns in the field during crop development. IMZ boundaries can then be adjusted according to within season plant feedback during the vegetative stage. Results showed that VIs were also significantly correlated with LWP; however, the correlations observed were weak and not feasible for creating any prediction models.
Further research is needed to explore the advantages of the delineation of IMZs that change during the season. Although leaf reflectance in the visible to NIR wavebands are not directly related to water status, additional studies can evaluate the relationship between visible and infrared-based VIs and water status. A larger data set from different fields exposed to different conditions can help test stability of these correlations. The potential of estimating LWP from remotely sensed multispectral data at a later stage in the season combined with plant height estimation during the vegetative stage could make dynamic IMZ delineation more feasible.
The comparison of the DVRI performance to conventional irrigation was not conclusive. This is because the conditions during both years of the study resulted in little difference in water applied to the conventional and DRVI treatments. The technology should be further evaluated on cotton as it has been shown to be effective during peanut production.
material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.