APPLICATIONS OF DRONE-BASED REMOTE SENSING IN CARROT AND TOMATO CROPPING SYSTEMS By Michael Abraham Metiva A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Horticulture—Master of Science 2021 ABSTRACT APPLICATIONS OF DRONE-BASED REMOTE SENSING IN CARROT AND TOMATO CROPPING SYSTEMS By Michael Abraham Metiva Using models of canopy height and vegetation indices (VIs), drone-based remote sensing (RS) can allow for large-scale assessments of plant growth and nutrient status. The goal of this thesis was to assess drone-based RS in addressing research and production challenges in processing carrots and fresh-market tomatoes. In processing carrots, a two-year trial was conducted to investigate the effects of topdress N rate, number of topdresses, and timing of split applications on processing carrot production, as well as the potential for VIs to guide N topdress decisions. Yield and shoot biomass were found to increase with higher N rates. Splitting applications did not affect yield but increased shoot biomass and N uptake in a wet year. Both early and late split applications showed potential to increase N loss. VI-based sufficiency indices explained at most 66% and 29% of the variation in carrot root yield in 2019 and 2020, respectively, but explained greater variation on average than petiole sap nitrate (6%) and recommended N applications 26% less often. In fresh-market tomatoes, RS was integrated into a cover crop by N fertilizer rate experiment to compare RS measurements to manual measurements of plant height, leaf tissue N, and leaf chlorophyll meter (SPAD) readings. Crop surface model plant height estimates were good estimators of measured heights (R2 = 0.89-0.96), with comparable correlations to final yields and ability to resolve significant treatment differences. Foliar N identified more significant differences between treatments than SPAD or VIs. In both experiments, drone-based RS demonstrated the potential to detect relevant in-season plant treatment responses comparably to manual measurements, with possible advantages in scalability, cost, and resolution. Copyright by MICHAEL ABRAHAM METIVA 2021 This thesis is dedicated to my family, who got me this far by their excellent example, to my advisor, Zack Hayden, who demonstrates the best ideals of academia, and to me, who worked very hard on it. iv ACKNOWLEDGEMENTS No one accomplishes anything alone, and that is especially true of this thesis. First and foremost, I’d like to thank Dr. Zack Hayden for his constant support and quality feedback on everything I’ve done in graduate school. By his nature, he creates an ideal environment for learning, working, and most importantly enjoying doing so. My other committee members, Dr. Erin Bunting and Dr. Kurt Steinke, have similarly shaped this thesis from beginning to end and I’ve learned more than I hoped from each of them. I’d also like to thank the friends and colleagues who shaped my experience into such a positive one. Olivia Davidson and Dr. Lin Liu helped me the most early on through their positive examples and friendship. Colin Phillippo deserves the credit for helping me into the Hayden Lab as well as being the ideal road trip duet partner. He also managed our research assistants, including Austin Green, Genny Feister, Patrick Squire, Sam Callow, and Laura Hudacek, without whom this work would have taken twice as long and been half as fun. Rounding out the vegetable crew were Alyssa Tarrant, Daniel Priddy, Monique Hemker, Nicole Soldan, and Dr. Dan Brainard, all of whom provided valuable aid and friendship during my time in graduate school. Celia Hallan, though not an official member of the vegetable crew, deserves an honorable mention for all her love and support. Last, but certainly not least, this project was supported by the Agriculture and Food Research Initiative Competitive Grant no. 2019-67013-29203 from the USDA National Institute of Food and Agriculture, the USDA AMS/MDARD Specialty Crop Block Grant #190000000309, the Michigan Carrot Committee, and the North Central SARE Graduate Student Grant GNC19-283. v TABLE OF CONTENTS LIST OF TABLES ....................................................................................................................... viii LIST OF FIGURES ....................................................................................................................... ix KEY TO ABBREVIATIONS ........................................................................................................ xi CHAPTER I: Introduction ...............................................................................................................1 Drone-Based Remote Sensing in Agriculture ......................................................................1 Nitrogen Management in Processing Carrot Production .....................................................2 Tracking Tomato Plant Health .............................................................................................3 Research Needs and Objectives ...........................................................................................4 REFERENCES ..........................................................................................................................6 CHAPTER II: Topdress Strategies and Remote Sensing for Nitrogen Management in Processing Carrots .......................................................................10 ABSTRACT .............................................................................................................................10 INTRODUCTION ...................................................................................................................11 MATERIALS AND METHODS .............................................................................................14 Field Trials .........................................................................................................................14 Soil and Plant Data.............................................................................................................16 Remote Sensing Data Collection .......................................................................................17 Image Processing ...............................................................................................................18 Vegetation and Sufficiency Indices ...................................................................................19 Statistical Analysis .............................................................................................................21 RESULTS .................................................................................................................................24 Weather ..............................................................................................................................24 Carrot Yield and Quality ....................................................................................................24 Nitrogen Uptake and Removal ...........................................................................................27 Sufficiency Indices .............................................................................................................30 Topdress Application Decisions ........................................................................................31 DISCUSSION ..........................................................................................................................35 Tradeoffs in N Topdress Strategies....................................................................................35 Topdress Decision Making Tools ......................................................................................38 CONCLUSIONS......................................................................................................................45 REFERENCES ........................................................................................................................46 CHAPTER III: In-Season Remote Sensing of Growth and N Status in Fresh-Market Tomato ...........................................................................................51 ABSTRACT .............................................................................................................................51 INTRODUCTION ...................................................................................................................52 MATERIALS AND METHODS .............................................................................................55 Field Trial...........................................................................................................................55 vi Data Collection ..................................................................................................................56 Image Processing ...............................................................................................................57 Crop Surface Models .........................................................................................................59 Statistical Analysis .............................................................................................................61 RESULTS .................................................................................................................................64 Plant Height .......................................................................................................................64 Plant N Status .....................................................................................................................66 Relationship to Total Yield ................................................................................................70 DISCUSSION ..........................................................................................................................73 Accuracy of Tomato Plant Height Estimations .................................................................73 Comparison of Methods for Assessing Plant N Status ......................................................75 Utility of Remote Sensing Metrics in Research and Management ....................................78 CONCLUSIONS......................................................................................................................82 REFERENCES ........................................................................................................................84 vii LIST OF TABLES Table 1.1. Summary of initial soil chemical characteristics at experimental sites ........................14 Table 1.2. Timeline of topdress events for different treatment groups ..........................................16 Table 1.3. Selected vegetation indices. R indicates measured reflectance with a subscript denoting the wavelength in nm. Subscripts r, g, and b represent the normalized red, green, and blue channels of the RGB camera, which does not have well-defined spectral bands....................................................................20 Table 1.4. Mean temperature and monthly rainfall for Hart, MI in 2019 and 2020, with NOAA 30-year (1981-2010) monthly averages for the same area .......................24 Table 1.5. Effects of N rate and topdress division on carrot root and shoot weight at harvest .......................................................................................................................25 Table 1.6. Effects of split topdress timing on carrot root and shoot weight at harvest ..................26 Table 1.7. Effects of N rate and topdress division on carrot root cull categories at harvest ...........................................................................................................................27 Table 1.8. Effects of N rate and division on plant N uptake, net N uptake, net N removal, and residual soil N .........................................................................................28 Table 1.9. Split topdress timing effects on plant N uptake, net N, and residual soil N .................29 Table 1.10. Predicted SI values corresponding to the mean total yield of the reference treatment in 2019 and 2020 ........................................................................31 Table 1.11. Percentages of plots receiving topdress N during experiment or indicated to receive topdress N through decision support tools. Data aggregated across sampling dates. Only dates with all three data types included: 2 and 16 July, 1, 19, and 29 August, 11 and 25 September 2019; 20 and 29 July, 12 and 26 August, 22/23 September 2020 ..............................34 Table 2.1. Selected vegetation indices related to plant N status. R indicates measured reflectance with a subscript denoting the wavelength in nm. Subscripts r, g, and b represent the normalized channels of the RGB camera, which does not have well-defined spectral bands ...........................................59 viii LIST OF FIGURES Figure 1.1. Lines of best fit for total yield and shoot fresh weight at harvest. Different fits were used for subsets of the data (front + none, split + none) where topdress division effects were found. Fit equations and coefficients of determination (R2) values are shown for corresponding regression fit lines ..........................................................................26 Figure 1.2. 2019 effects of interaction effects between N rate and topdress division on N uptake and net N uptake and removal calculations. ** and *** indicate significant differences in marginal means of p < 0.01 and p < 0.001, respectively, at a given N rate....................................................................29 Figure 1.3. Coefficients of determination between total yield and visible VIs, multispectral VIs, and manual sampling data (petiole nitrate, shoot N percentage, root fresh weight, shoot dry weight) over the 2019 and 2020 seasons ................................................................................................................30 Figure 1.4. GLI- and NDVI-based sufficiency indices of treatments over the late July through early October sampling period, 2019 and 2020. Error bars show standard error of the treatment mean for each sampling date. Mean minimum thresholds for each year are shown as solid lines.............................32 Figure 1.5. Petiole sap nitrate concentrations for different treatments. The minimum thresholds published by Michigan State University are shown as a solid line in both years for all treatments .............................................................................33 Figure 1.6. Contingency table summarizing hypothetical topdress N decisions (Y “yes”; N “no”) indicated by three testing methods (petiole sap nitrate test, GLI-based SI, and NDVI-based SI) for all plots over sampling dates with all three data types (n=432) ........................................................34 Figure 2.1. Mean measured plant heights compared to CSM1 (single initial DTM) and CSM2 (unique DTM per drone flight) estimates of plant height with all data from each year combined (2018 n=316; 2019 n=448). The dotted line represents a 1:1 ratio, and the solid line represents the ordinary least squares fit of the data ............................................................................64 Figure 2.2. R2 values for correlations between measured plant heights and estimates of plant height derived from CSM1 (single initial DTM) and CSM2 (unique DTM per drone flight) for each individual sampling date in 2018 and 2019 .......................................................................................................................65 ix Figure 2.3. Contingency tables summarizing the significance of pairwise comparisons of marginal means calculated from mixed-effects models with different response variables (measured plant heights vs. CSM1, CSM2 height estimates). Pairwise comparisons were performed between marginal means computed on each date for levels of cover crop (n=72), levels of cash crop N rate (n=72), and all individual treatments (n=1440) .................66 Figure 2.4. R2 values for linear correlations between foliar N concentration and SPAD readings, visible indices, and multispectral indices for individual sampling dates. Both foliar N data and any other data type (in this case, visible indices) were collected on only one date in 2018; three dates with shared data are available and shown for all data types in 2019 ...................................67 Figure 2.5. Contingency tables summarizing the significance of Tukey-adjusted pairwise comparisons of marginal means calculated from mixed-effects models with different response variables (foliar N concentration vs. SPAD, visible indices, and multispectral indices). Pairwise comparisons were performed over levels of the cover crop and cash crop N rate factors (SPAD n=18; visible indices n=24; multispectral indices n=18) as well as the interactions (SPAD n=360; visible indices n=480; multispectral indices n=360) .............................................................................................................69 Figure 2.6. R2 values for linear correlations between total tomato yield and plant heights (measured, CSM1 estimated, CSM2 estimated) in 2018 and 2019 ................70 Figure 2.7. Coefficients of determination between total tomato yield and foliar N % (measured leaf N content), SPAD measurements (based on leaf spectral transmittance) and visible and multispectral indices from remote sensing (canopy-level averages) for individual sampling dates. Multispectral data available in 2019 only .................................................................................................71 x KEY TO ABBREVIATIONS Above ground level Blue-green index Crop surface model Digital surface model Digital terrain model Ground control point Green leaf index Green normalized difference vegetation index Normalized difference red-edge index Normalized difference vegetation index Normalized pigment chlorophyll index Red green blue Sufficiency index Vegetation index AGL BGI CSM DSM DTM GCP GLI GNDVI NDRE NDVI NPCI RGB SI VI xi CHAPTER I: Introduction Drone-Based Remote Sensing in Agriculture As drone technology has become more accessible over the past decade, its use in agriculture has become an intense area of interest for researchers and growers alike. Platforms for remote sensing had previously been quite limited. Satellite imagery, while often accessible for free through government sources, suffers from low spatial resolution and essentially fixed temporal resolution along with significant atmospheric interference. Imagery collection using airplanes offers increased spatial resolution and more flexible timing, improving its relevance for agricultural research and management, but it is also limited by higher prices. By contrast, remote sensing using drones allows for very high-resolution (~0.1 m or less) imagery to be collected essentially on-demand, with a relatively low financial and temporal cost of acquiring and utilizing the required equipment (Mulla, 2013). Information about plant stress can be inferred based on the spectral reflectance of plant canopies, and imagery measuring this reflectance across various spectral bands can allow for large-scale quantitative assessments of nutrient deficiencies, canopy cover, and other biophysical parameters of interest (Jackson, 1986). Many vegetation indices have been developed to use canopy reflectance to assess nutrient deficiencies, with N deficiencies in particular linked to increases in red reflectance (i.e. a yellowing of the canopy); the popular normalized difference vegetation index (NDVI) for example takes advantage of this change as well as smaller changes in the infrared reflectance to assess plant N status and growth for a wide variety of applications (Pinter et al., 2003; Rouse et al., 1973). In addition to recording reflectance data for the calculation of vegetation indices, the high-resolution images taken using drones can be used as inputs in structure-from-motion techniques to generate 3-dimensional digital surface models. 1 Estimates of plant height and biomass have been generated using DSMs in a variety of field crops, particularly for high-throughput phenotyping applications (Bendig et al., 2013, 2015; Brocks & Bareth, 2018; De Souza et al., 2017; Holman et al., 2016), but considerably less research has been conducted in high value, lower acreage specialty crops like vegetables. Nitrogen Management in Processing Carrot Production Plant nitrogen (N) status is a major determining factor for yield, quality, and mechanical harvest efficiency in processing carrots, so effective N management is a key aspect of the crop’s production (Batra & Kalloo, 1990). Processing carrots have relatively high N demands and a long growing season, typically April through October; these factors, combined with a high risk of nitrate leaching due to irrigation and sandy soils, present many important decision points for growers to ensure the timely delivery of N fertilizer to the plants (Warncke et al., 2004). A combination of starter N and topdress N applications are typically used to ensure a steady supply of N throughout the season while minimizing the risk of leaching. When exactly to apply these topdresses, how many topdresses are needed, and the tradeoffs associated with different application strategies are still areas of concern, however. In order to make more informed topdressing decisions, some carrot farmers use petiole sap nitrate testing to monitor the N status of the carrot crop. This test involves manually sampling carrot shoots, pressing the leaf petioles to express their sap, and testing the sap with a handheld nitrate meter. Growers can then determine the crop’s need for topdress N by comparing the sampled nitrate values to minimum thresholds. These thresholds are not well-established, however, and use of the petiole sap nitrate test is limited by the amount of labor required for sampling and processing (Selk et al., 2004; Westerveld et al., 2003, 2007). 2 Tracking Tomato Plant Health Methods for non-destructive assessment of plant health are a key component of data collection regimes in both agricultural research and production systems, allowing data collection without killing limited or valuable plants. Such measurements are often used as proxies or estimators for biophysical characteristics which would be too costly or labor-intensive to sample directly. For example, plant height can be used as a metric for evaluating plant health and nutrient status for research applications, often as an indicator of biomass. In grain crops, for example, plant height is a phenotype that is regularly evaluated in plant breeding and genetics research as it is related to certain yield parameters (Jimenez-Berni et al., 2018; Pérez- Harguindeguy et al., 2013; Salas Fernandez et al., 2009). Plant height is also a common metric used in tomato and other vegetable research, particularly in cases where destructive biomass sampling is not practical and in crops with an upright growth form (Dong et al., 2005; Fontes & de Araujo, 2006; Gianquinto et al., 2006; Muñoz-Carpena et al., 2008). Fresh-market tomatoes in particular are often staked or string trained, making their growth more upright and increasing the relevance of plant height as a proxy for biomass. Leaf chlorophyll concentration is another commonly used metric which can be sampled non-destructively. This parameter can be estimated rapidly and at relatively low cost at the leaf level using handheld meters such as the SPAD-502 chlorophyll meter, which delivers unitless readings proportional to plant chlorophyll content (Ling et al., 2011). SPAD readings have been used to characterize differences in plant N status and assess requirements for supplemental N fertilization in a wide variety of crops, including tomato (Chua et al., 2003; Gianquinto et al., 2006; Hussain et al., 2000; Sandoval-Villa et al., 2000; Varvel et al., 1997). Sampling leaf tissue (or petiole sap) to analyze foliar N concentration is a more direct measure of plant N status, but 3 the additional labor and costs for laboratory analysis can limit its use and sampling frequency. Indirect, non-destructive sampling methods therefore have advantages in practice. Research Needs and Objectives Compared to field crops like maize or wheat, little drone-based remote sensing research has been done in vegetables, and the broad range of biophysical variety across the vegetable crops means much of the work in cereals may not be directly transferable (Croft et al., 2014). The overall goal of this thesis was to assess the utility of drone-based remote sensing data in different vegetable cropping systems with their own specific research and production needs. This information can guide the adoption of this technology in these crops as well as further elucidate the possible applications of drone technology in vegetable cropping systems in general. In both chapters, comparisons of remote sensing to other data collection methods are included, each with their unique advantages and disadvantages. Chapter II presents results from a two-year trial conducted in active commercial processing carrot fields. The first objective of this trial was to determine the impacts of different topdress strategies on processing carrot yield, shoot biomass, and the risk of N loss, addressed through treatments representing a range of N rates, topdress application numbers, and split topdress timings. The second objective was to assess the possible utility of drone-based remote sensing as an effective topdress N management tool for processing carrots. Petiole sap nitrate testing is sometimes used for this purpose but is labor intensive and spatially incomplete. The N management treatments generated a range of plant N statuses and yield responses, allowing the comparison of petiole sap nitrate testing and other manual sampling regimes to remote sensing derived sufficiency indices as topdress decision support tools. 4 Chapter III presents results from a two-year fresh-market tomato trial conducted at the MSU Horticulture Teaching and Research Center where we evaluated the ability of drone-based remote sensing data to distinguish differences in tomato growth and N status. Using variability in fresh-market tomato biophysical characteristics within a two-year cover crop by N fertilizer rate experiment, the objective of this research was to compare drone-based remote sensing measurements (crop surface model [CSM]-derived height estimates and canopy spectral reflectance-derived vegetation indices [VIs]) to common manual measurements (plant height, leaf tissue N, and leaf chlorophyll meter [SPAD] readings). Manual sampling for physical parameters like plant height and foliar N is labor intensive and can involve additional laboratory costs, and SPAD readings can be highly variable depending on the crop or user sampling technique (Muñoz-Huerta et al., 2013). The range of plant biophysical characteristics generated through the different treatments, sampled via both manual data collection and drone-based remote sensing, created the opportunity to assess the different tools’ sensitivities to treatment differences and ability to estimate ground truth metrics. 5 REFERENCES 6 REFERENCES Batra, B. R., & Kalloo. (1990). Effect of different levels of irrigation and fertilization on growth and yield of carrot (Daucus carota L.) for root production. Vegetable Science, 17(2), 127– 139. Bendig, J., Bolten, A., & Bareth, G. (2013). UAV-based Imaging for Multi-Temporal, very high Resolution Crop Surface Models to monitor Crop Growth Variability. Photogrammetrie, Fernerkundung, Geoinformation, 2013(6), 551–562. Bendig, J., Yu, K., Aasen, H., Bolten, A., Bennertz, S., Broscheit, J., Gnyp, M. L., & Bareth, G. (2015). Combining UAV-based plant height from crop surface models, visible, and near infrared vegetation indices for biomass monitoring in barley. International Journal of Applied Earth Observation and Geoinformation, 39, 79–87. Brocks, S., & Bareth, G. (2018). Estimating barley biomass with crop surface models from oblique RGB imagery. Remote Sensing, 10(2). Chua, T. T., Bronson, K. F., Booker, J. D., Keeling, J. W., Mosier, A. R., Bordovsky, J. P., Lascano, R. J., Green, C. J., & Segarra, E. (2003). In-Season Nitrogen Status Sensing in Irrigated Cotton. Soil Science Society of America Journal, 67(5), 1428–1438. Croft, H., Chen, J. M., & Zhang, Y. (2014). The applicability of empirical vegetation indices for determining leaf chlorophyll content over different leaf and canopy structures. Ecological Complexity, 17(1), 119–130. De Souza, C. H. W., Lamparelli, R. A. C., Rocha, J. V., & Magalhães, P. S. G. (2017). Height estimation of sugarcane using an unmanned aerial system (UAS) based on structure from motion (SfM) point clouds. International Journal of Remote Sensing, 38(8–10), 2218–2230. Dong, J., Wu, F., & Zhang, G. (2005). Effect of cadmium on growth and photosynthesis of tomato seedlings. Journal of Zhejiang University SCIENCE, 6B(10), 974–980. Fontes, P. C. R., & de Araujo, C. (2006). Use of a chlorophyll meter and plant visual aspect for nitrogen management in tomato fertigation. Journal of Applied Horticulture, 8(1), 8–11. Gianquinto, G., Sambo, P., & Borsato, D. (2006). Determination of SPAD threshold values for the optimisation of nitrogen supply in processing tomato. Acta Horticulturae, 700, 159–166. Holman, F. H., Riche, A. B., Michalski, A., Castle, M., Wooster, M. J., & Hawkesford, M. J. (2016). High throughput field phenotyping of wheat plant height and growth rate in field plot trials using UAV based remote sensing. Remote Sensing, 8(12). Hussain, F., Bronson, K. F., Yadvinder-Singh, Bijay-Singh, & Peng, S. (2000). Use of chlorophyll meter sufficiency indices for nitrogen management of irrigated rice in Asia. 7 Agronomy Journal, 92(5), 875–879. Jackson, R. D. (1986). Remote Sensing of Biotic and Abiotic Plant Stress. Annual Review of Phytopathology, 24, 265–287. Jimenez-Berni, J. A., Deery, D. M., Rozas-Larraondo, P., Condon, A. T. G., Rebetzke, G. J., James, R. A., Bovill, W. D., Furbank, R. T., & Sirault, X. R. R. (2018). High throughput determination of plant height, ground cover, and above-ground biomass in wheat with LiDAR. Frontiers in Plant Science, 9(February), 1–18. Ling, Q., Huang, W., & Jarvis, P. (2011). Use of a SPAD-502 meter to measure leaf chlorophyll concentration in Arabidopsis thaliana. Photosynthesis Research, 107(2), 209–214. Mulla, D. J. (2013). Twenty five years of remote sensing in precision agriculture: Key advances and remaining knowledge gaps. Biosystems Engineering, 114(4), 358–371. Muñoz-Carpena, R., Icerman, J., Dukes, M. D., Zotarelli, L., & Scholberg, J. M. (2008). Tomato yield, biomass accumulation, root distribution and irrigation water use efficiency on a sandy soil, as affected by nitrogen rate and irrigation scheduling. Agricultural Water Management, 96(1), 23–34. Muñoz-Huerta, R. F., Guevara-Gonzalez, R. G., Contreras-Medina, L. M., Torres-Pacheco, I., Prado-Olivarez, J., & Ocampo-Velazquez, R. V. (2013). A review of methods for sensing the nitrogen status in plants: Advantages, disadvantages and recent advances. Sensors (Switzerland), 13(8), 10823–10843. Pérez-Harguindeguy, N., Díaz, S., Garnier, E., Lavorel, S., Poorter, H., Jaureguiberry, P., Bret- Harte, M. S., Cornwell, W. K., Craine, J. M., Gurvich, D. E., Urcelay, C., Veneklaas, E. J., Reich, P. B., Poorter, L., Wright, I. J., Ray, P., Enrico, L., Pausas, J. G., De Vos, A. C., … Cornelissen, J. H. C. (2013). New handbook for standardised measurement of plant functional traits worldwide. Australian Journal of Botany, 61(3), 167–234. Pinter, P. J., Hatfield, J. L., Schepers, J. S., Barnes, E. M., Moran, M. S., Daughtry, C. S. T., & Upchurch, D. R. (2003). Remote Sensing for Crop Management. Photogrammetric Engineering and Remote Sensing, 69(6), 647–664. Rouse, J. W., Haas, R. H., Schell, J. A., & Deering, D. W. (1973). Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation. Salas Fernandez, M. G., Becraft, P. W., Yin, Y., & Lübberstedt, T. (2009). From dwarves to giants? Plant height manipulation for biomass yield. Trends in Plant Science, 14(8), 454– 461. Sandoval-Villa, M., Guertal, E. A., & Wood, C. W. (2000). Tomato leaf chlorophyll meter readings as affected by variety, nitrogen form, and nighttime nutrient solution strength. Journal of Plant Nutrition, 23(5), 649–661. 8 Selk, G. E., LeValley, R. C., Highfill, G. A., & Buchanan, D. S. (2004). Comparing the Accuracy of the Cardy Portable Nitrate Meter with Laboratory Analysis of Nitrate Concentrations in Summer Annual Forages. 3, 1–2. Varvel, G. E., Schepers, J. S., & Francis, D. D. (1997). Ability for In-Season Correction of Nitrogen Deficiency in Corn Using Chlorophyll Meters. Soil Science Society of America Journal, 61(4), 1233–1239. Warncke, D., Dahl, J., & Zandstra, B. (2004). Nutrient Recommendations for Vegetable Crops in Michigan (E2904). In Michigan State University Extension Bulletin (Issue E2934). Westerveld, S. M., McKeown, A. W., McDonald, M. R., & Scott-Dupree, C. D. (2003). Chlorophyll and nitrate meters as nitrogen monitoring tools for selected vegetables in southern Ontario. Acta Horticulturae, 627, 259–266. Westerveld, Sean M., McDonald, M. R., & McKeown, A. W. (2007). Establishment of critical sap and soil nitrate concentrations using a Cardy nitrate meter for two carrot cultivars grown on organic and mineral soil. Communications in Soil Science and Plant Analysis, 38(13– 14), 1911–1925. 9 CHAPTER II: Topdress Strategies and Remote Sensing for Nitrogen Management in Processing Carrots ABSTRACT Improved nitrogen (N) fertilizer management strategies are required to maximize yield, quality, and harvest efficiency in processing carrots (Daucus carota L.) while minimizing risks of N loss. A 2-year field study was conducted in active commercial carrot fields near Hesperia and Pentwater, MI, U.S. to investigate the effects of four N fertilizer rates (29, 67, 135, and 202 kg ha-1 N), two N topdress strategies (single or 3 split), and timing of split applications on carrot production and N utilization. In addition, the potential for remote sensing-based vegetation indices (VIs) to guide in-season N topdress timing decisions was evaluated relative to conventional sampling methods, including petiole sap nitrate testing. Carrot root yield and shoot biomass increased with greater N rates without reaching a clear yield plateau. Multiple split N topdress applications did not affect yield relative to a single N application but did increase shoot biomass and N uptake in 2019. Both earlier and later than typical start timings for split N application exhibited the potential to increase N loss. While VI-based sufficiency indices explained at most 66% and 29% of the variation in carrot root yield in 2019 and 2020, respectively, the indices consistently explained greater variation as compared to in-season measures of petiole sap nitrate (6%), shoot N concentration (8%), and carrot root (10%) and shoot (14%) weights. Hypothetical N topdress decisions made using VI-based sufficiency indices were generally more conservative than petiole sap nitrate thresholds and recommended N 26% less often. Despite labor and technological tradeoffs, remote sensing may increase accuracy and scalability of N topdress decision support in processing carrot production. 10 INTRODUCTION Effective nitrogen (N) management is integral to optimizing carrot (Daucus carota L.) yield and quality while minimizing N losses to the environment. Processing carrot varieties in particular have high N demand, a long growing season of approximately six months, and are grown in sandy, irrigated soils creating a greater risk for leaching nitrate losses (Warncke et al., 2004). To satisfy carrot N demands, multiple topdress N applications are often applied throughout the season in addition to at-plant starter N. There are three components to this N management that carrot growers must decide upon each season. The first is the total N rate for the season, which informs the amount of topdress N applied regardless of strategy. The second is the number of split N applications with more applications costing additional labor and money. The third is the timing of individual topdress applications in order to synchronize N availability with plant uptake. Growers can address N risks and uncertainty in part by utilizing the petiole sap nitrate test. Petiole sap nitrate testing compares sampled nitrate values to a threshold value, prompting topdress N applications when sampled values fall below the threshold. Due to differences in weather conditions and carrot varieties, however, critical threshold values are not well- established and may only provide data concerning smaller sample areas within a larger field (Selk et al., 2004; Westerveld et al., 2003; Westerveld et al., 2007). Improving the spatial completeness of the data as well as reducing the labor required for sampling can improve the ability of farmers to make in-season N topdressing decisions. Drone-based remote sensing has been widely studied and used for crop N management, allowing for rapid, on-demand collection of high-resolution, spatially continuous data (Pinter et al., 2003; Usha & Singh, 2013). Vegetation indices (VIs) have been developed utilizing spectral 11 reflectance captured in aerial imagery to assess the health of crops and predict yields, but data inconsistencies between growing seasons and cultivars as well as challenges in identifying sources of variability have been observed (Pinter et al., 2003). The sufficiency index (SI) approach to N management was developed to address these issues by normalizing measured values, either VIs or other N status metrics like SPAD-502 chlorophyll meter readings, using values obtained from well-fertilized or non-N limiting reference strips. A SI value of 1.0 would indicate very similar readings/values to the reference area, while lower SI values would indicate lower readings/values than the reference area. By dividing by reference strip values, the SI approach allows the assessment of relative differences in these metrics with the assumption that differences are primarily caused by discrepancies in N status, with a minimum threshold value such as 0.95 often used to indicate N deficiencies (Blackmer & Schepers, 1995). SI values have been researched for use in a variety of crops and using a variety of sensors to determine N fertilizer needs across the broader field, increasing N use efficiency and mitigating consistency issues due to differences in soil, cultivar, and growth stage (Blackmer & Schepers, 1995; Mulla, 2013; Oliveira et al., 2013; Raun et al., 2002; Stamatiadis et al., 2020). Unlike the petiole sap nitrate test, however, SIs based on remote sensing metrics do not directly measure actual plant N concentrations, instead relying on plant coloration and size to assess N status. This can complicate the interpretation of SI values, as other nutrient deficiencies or diseases which are not uniformly distributed could impact SI calculations and violate the assumption that SI differences only reflect differences in plant N status. One objective of this research was to investigate the effects of total N rate, number of topdress applications (frontloaded or split), and topdress application timing on carrot yield, quality, and N uptake and removal. The second objective was to assess the suitability of drone- 12 based remote sensing utilizing VI-based SIs calculated from drone imagery as a topdress timing decision-making tool and to compare results against petiole sap nitrate testing. We hypothesized that the recommended N rate, split topdresses, and typical topdress timing would optimize yield. We further hypothesized that splitting topdresses would increase N uptake and removal compared to frontloaded topdresses. With respect to remote sensing, we hypothesized that SI values would be more strongly correlated with final yields than petiole sap nitrate concentrations at each sampling point during the season. 13 MATERIALS AND METHODS Field Trials From June to October 2019 an extensive field campaign was initiated in an active commercial agricultural field during carrot production in Hesperia, Michigan, USA (43.5552° N, 86.0958° W; “Site 1”). This region has a history of growing both field crops and vegetables in rotation. Site 1 soils included a Spinks-Okee complex and a Coloma-Toogood complex. The experiment was repeated at another commercial agricultural field in 2020, located northwest of Site 1 near Pentwater, MI (43.8267° N, 86.3827° W; “Site 2”). The soil at Site 2 was primarily Fern fine sand. A summary of initial soil characteristics is shown in Table 1.1. Weather data for both sites was sourced from weather stations located near Hart, Michigan (Enviroweather station 43.7366° N, -86.3594° W; NOAA station 43.6747° N, 86.4238° W), between 10 and 30 km from the sites, and may not perfectly represent actual weather conditions at the sites. Table 1.1. Summary of initial soil chemical characteristics at experimental sites. Site 1 was cropped to field corn in 2018; no small grain nurse cover crop was used to reduce windblown sediment in 2019 due to the presence of field corn residue. Site 2 was cropped in the same pattern during the 2019 and 2020 seasons. Both sites were prepared for carrot planting using a custom strip tiller consisting of a shank, bed forming disks, and rotary cultivator. Processing carrots of the cultivar Cupar were seeded in 3-row beds on 20 April 2019 and 24 April 2020 at Sites 1 and 2, respectively. Carrot rows were spaced 0.46 m apart within each bed with 1.63 m between bed centers. A total of 29 kg ha-1 N was applied as liquid starter fertilizer in both years. In 2019 the mixture contained 99 kg ha-1 UAN (28-0-0) and 13 kg ha-1 ATS (12-0-0- 14 YearSitepHCECOrganic MatterInorganic NP (Bray-1)KMgCacmol kg-1%201917.14.41.86.714514473692202026.53.71.18.97010786538----------------------------------------------- mg kg-1 ----------------------------------------------- 26S), and in 2020 it contained 62 kg ha-1 UAN, 10 kg ha-1 ammonium polyphosphate (10-34-0), and 13 kg ha-1 ATS. Experimental plots 4.88 m (three beds) across and 12.19 m deep were arranged in a randomized complete block design with four replications. Eight ground control points were placed in the field at both sites to allow precise georeferencing and photogrammetry of remote sensing imagery. A low-N control treatment used starter fertilizer only. A true zero-N check could not be included due to the commercial nature of the experimental sites where growers utilize starter N to ensure rapid early carrot growth; this is a key component of N management particularly in northern climates. Other treatments included a ramp of season-total N rates (67, 135, and 202 kg ha-1 N) corresponding to 50, 100, and 150% of the recommended season-total N rate for processing carrots on mineral soils in Michigan (Warncke et al., 2004). These rates were applied as starter fertilizer plus topdressed urea. Due to projections of a hot, dry summer in 2020, topdresses in that year were applied as urea mixed with Agrotain at 2.09 L t-1 to mitigate volatilization losses. Topdresses occurred as either a single frontloaded application or three equally sized split applications. Frontloaded applications occurred entirely on the dates corresponding to the first split applications, with split applications continuing at roughly 4 wk intervals. To further investigate the role of topdress application timing, two additional treatments were included that follow Michigan State University’s recommended grower practice (135 kg ha- 1 N applied as starter plus split topdress applications) with the exception that split topdress applications were initiated approximately 2 wk earlier and 2 wk later than the other split treatments. The 2019 and 2020 topdressing schedules for these different groups of treatments are shown in Table 1.2. 15 Table 1.2. Timeline of topdress events for different treatment groups. Soil and Plant Data Baseline inorganic N (sum of NO3 --N and NH4 + -N) levels were established by soil sampling to 20 cm depth concurrent with plot establishment on 5 June 2019 and 18 June 2020. Eight soil cores were collected from each plot for this purpose in 2020, whereas 12 soil cores were collected per replicate in 2019 due to labor constraints. Eight 20 cm cores were collected from each plot again immediately prior to harvest on 29 October 2019 and 13 October 2020 to measure residual soil N. After plot establishment, site visits occurred at approximately 2 wk intervals to accommodate the fertilizer application schedule and collect plant and soil data. In addition to the fertilizer application dates, data was also collected on 5 June, 25 September, 9 October, and 29 October 2019 as well as 23 September and 13 October 2020. During each site visit, five carrot plants were collected per plot on each date with the exception of 5 June 2019, 20 June 2019, and 1 July 2020 when 10 plants were sampled per plot due to their small size. Only plots corresponding to treatments 1 and 8 were sampled on 1 July 2020 due to labor limitations arising from the COVID-19 pandemic. The youngest fully-expanded leaf was removed from sampled plants, and the petiole sap was expressed and analyzed for nitrate concentration. Technical issues with available handheld nitrate meters led to the use of a Lachat QuikChem flow injection analysis system to test for 16 20192020Topdress Events20 June18 JuneEarly topdress #12 July1 JulyFrontloaded topdresses, split topdresses #116 July20 JulyEarly topdress #2, late topdress #11 Aug29 JulySplit topdresses #219 Aug12 AugEarly topdress #3, late topdress #229 Aug26 AugSplit topdresses #311 Sept10 SeptLate topdress #3 nitrate-N concentrations which were then converted to nitrate concentrations. The remaining shoots were separated from the roots, dried at 60˚C for 10-14 days, and weighed dry. Roots were weighed fresh and diameters measured at the widest point. To estimate final yield, carrots were harvested shortly before commercial harvest (29 October 2019, 13 October 2020) from the middle 6.10 m section of the center row of the center bed in each plot. Roots and shoots were separated and the shoots weighed fresh. Before being weighed, harvested carrots were graded into seven categories based on market specifications, including marketable (no major defects), forked (forked end), split (split along the skin), small (less than 3.5 cm diameter), nub (length and diameter approximately equal with rounded end), animal damaged, and rotten categories. Five additional carrot plants were sampled per plot on harvest dates. Roots and shoots were separated, and the roots were weighed fresh while the shoots were dried at 60˚C for 10-14 days before being weighed. Root tissue subsamples were taken from the center inch of each root, diced, weighed, and dried before being weighed again. The dried root and shoot tissue samples were then analyzed for total N percentage. The percent root moisture, derived from the fresh and dry weights, was used to convert the total yield to dry weight for root N calculations. Shoot samples were not weighed fresh prior to drying, so a 10% moisture level was assumed and used to convert shoot biomass at harvest to dry weight. Root and shoot harvest dry weights then were multiplied by their respective total N percentages to assess plant N uptake. N removal was assumed to be equivalent to root N content at harvest. Remote Sensing Data Collection A DJI Phantom 4 Pro quadcopter equipped with a 20-megapixel natural color camera was used to collect red-green-blue (RGB) composite imagery of the experimental sites. The Phantom 4 RedEdge Mount Kit (Sky Flight Robotics, Inc.) was also affixed to the quadcopter to allow the 17 use of a MicaSense RedEdge-MX camera and the collection of multispectral image data at the 475 nm (blue), 560 nm (green), 668 nm (red), 717 nm (red-edge), and 842 nm (near-infrared) spectral bands. Remote sensing data was collected on the same dates as field data at an altitude of 15 m AGL with 85% image front overlap and 75% side overlap for the RGB camera; grid maps defining the parameters for these flights were produced using the DJI GS Pro app run on an Apple iPad Mini. Prior to each flight, images of a MicaSense Calibrated Reflectance Panel were taken using the RedEdge-MX camera. Multispectral calibration images were not taken on 2 July 2019 due to an error, and locational data issues with the RGB imagery taken on 9 October 2019 rendered the images unusable. The RedEdge-MX camera also malfunctioned on 10 September 2020 resulting in total loss of data for this date. Remote sensing data was thereafter collected by the Michigan State University RS&GIS office on 22 September and 9 October 2020. On these dates, RGB data was collected at 15 m AGL using a DJI Phantom 4 Pro while RedEdge-MX imagery was collected using a DJI M210 V2 quadcopter at 25 m AGL. As normal, images of a MicaSense Calibrated Reflectance Panel were taken prior to each flight. Image Processing Orthomosaics were generated from both the RGB camera images and the RedEdge-MX images using the photogrammetry software Pix4Dmapper. For each flight date, a composite orthomosaic was produced from the RGB camera imagery and five single-band orthomosaics (blue, green, red, red-edge, and near-infrared) were produced from the RedEdge-MX imagery. The composite and single-band orthomosaics had resolutions of approximately 0.005 m and 0.01 m respectively. The UTM coordinates of the ground control points were derived from the first RGB composite orthomosaic at each site (5 June 2019, 18 June 2020) and used to georeference 18 following orthomosaics. Calibrated reflectance panel images as well as data from the downwelling light sensor of the RedEdge-MX camera were used with Pix4Dmapper to control for ambient light conditions in the multispectral images. The open-source geographic information system application QGIS (v3.2) was used to digitize experimental plots. For each site, the ‘Create Grid’ tool was used first to generate a grid of 4.88 m by 12.19 m rectangular polygons. This shapefile was manually aligned with the first RGB composite orthomosaic at each site using the ‘Advanced Digitizing’ toolbar, and the ‘Buffer’ tool was used to create a 0.3 m inward buffer to remove the plot edges. Further image processing and was performed using R statistical software (R Core Team, 2020). A buffer of 0.5 m was applied to the extent of the plot area shapefiles to avoid edge issues during resampling, and this extent was used to resample the orthomosaics from their native resolutions to standard 0.005 m (RGB composite) and 0.01 m (single-band) resolutions using the nearest neighbor method. The RGB composite was then separated into its constituent red, green, and blue bands, and each band was divided by the sum of the three bands to normalize the individual bands. Vegetation and Sufficiency Indices The normalized RGB and multispectral bands were used to calculate a selection of VIs associated with leaf chlorophyll concentration and canopy cover and, by extension, plant N status (Table 1.3). Three VIs calculated using only the RGB camera data (BGI, NPCI, GLI; “visible indices”) were included alongside three VIs calculated using multispectral camera data (NDVI, NDRE, GNDVI; “multispectral indices) to assess the relative performance of indices in these two groups given the differing equipment required. The visible indices BGI and NPCI have been found to be correlated with chlorophyll and carotenoid concentrations, while the GLI was 19 Table 1.3. Selected vegetation indices. R indicates measured reflectance with a subscript denoting the wavelength in nm. Subscripts r, g, and b represent the normalized red, green, and blue channels of the RGB camera, which does not have well-defined spectral bands. developed to separate vegetation in imagery based on chlorophyll’s absorption of red and blue light and reflectance of green (Louhaichi et al., 2001; Peñuelas et al., 1994; Zarco-Tejada et al., 2005). NDVI was first developed using satellite imagery to assess large-scale changes in green vegetation, but has since been used to monitor crop canopy development and biomass with a wide variety of sensors and platforms (Mulla, 2013; Rouse et al., 1973). NDRE and GNDVI were subsequently developed to improve estimation of leaf chlorophyll content, with NDRE using differences in red-edge reflectance between healthy and stressed vegetation and GNDVI using changes in chlorophyll’s green reflectance (Barnes et al., 2000; Buschmann & Nagel, 1993). Mean values of these VIs were calculated over the plot areas. In order to convert plot mean VI values to sufficiency indices, a reference value was first calculated for each sampling date by averaging the VI values of the four plots which received the 202 kg N ha-1, frontloaded topdress treatment. Of the treatments included in this study, this treatment (“the reference treatment”) was considered the best approximation of a non-N limited reference area due to the 20 Index NameAcronymFormulaSourceNormalized PigmentChlorophyll IndexNormalized DifferenceVegetation IndexNormalized DifferenceRed-Edge IndexGreen Leaf IndexBlue-Green IndexGreen NDVIZarco-Tejada et al., 2005Peñuelas et al., 1994Louhaichi et al., 2001Rouse et al., 1973Barnes et al., 2000Buschmann & Nagel, 1993GNDVI(R842 – R560) / (R842 + R560)(R842 – R717) / (R842 + R717)(R842 – R668) / (R842 + R668)(2Rg – Rr – Rb) / (2Rg + Rr + Rb)(Rr – Rb) / (Rr + Rb)BGINPCIGLINDVINDRERb / Rg high N rate applied at a single early topdress; however, subsequent analysis of yield response suggests this treatment was not truly non-N limiting under the conditions of our experiment (Blackmer & Schepers, 1994). Plot mean VI values were then divided by the reference value for the corresponding date, yielding SI values for each experimental plot on each sampling date from 20 June through 9 October 2019 and 18 June through 8 October 2019. For comparison with direct physical sampling, SIs were also calculated based on manual plant sampling throughout the season using root fresh weight, shoot dry weight, and shoot N percentage. Statistical Analysis Statistical analysis was performed using R statistical software. The fixed effects of total season N rate and topdress division on defect-free and total yield, shoot fresh weight, and root:shoot ratio were evaluated for each year, along with the percentage of harvest by weight of the forked, split, small, and nub defect categories. The data evaluated using this model excluded the starter-only treatment as well as the early and late topdress treatments to maintain balanced data. Dunnett’s test was then used to compare the individual treatments within this data subset to the low N starter-only treatment (“the control treatment”) to further develop the interpretation of topdress rate and division effects. The fixed effects of split topdress timing on the same response variables were evaluated using a subset of data which included the early and late topdress treatments as well as the split 135 kg ha-1 N rate treatment on the intermediate timing. Replicate was included as a random blocking factor in both models. The effects of season N rate, number of topdresses (frontloaded or split), and split topdress timing on variables related to soil and plant N were also analyzed using the same models and data subsets including the use of Dunnett’s test. These variables included total plant N uptake, net N uptake (applied N minus total plant N uptake), net N removal (applied N minus root N at harvest), and residual soil N at harvest. 21 In order to improve assumptions of normality, the percentages of harvest by weight of the forked, split, and nub defect categories required square-root transformation. In cases where heteroskedasticity was detected for one or more factors, weighted regression was used with different variance weights for each level of the relevant factors. Type 3 ANOVAs were performed to determine the significance of model terms, and estimated marginal means were computed for factor levels. Pairwise comparison of these means was performed using Tukey’s HSD (Hothorn et al., 2008; Lenth, 2020). In combination with this relatively conservative multiple comparison procedure, a less conservative significance level of 0.10 was used in order to allow for discussion of marginally significant results. In the case of the net N uptake and net N removal, two-tailed t tests were also used to assess the difference of marginal means from zero. For each sampling date, plot mean SI values were regressed against total yield. Additional regressions between total yield and petiole sap nitrate concentration in addition to SIs based on root fresh weight, shoot dry weight, and shoot N percentage were also calculated. Regression analysis demonstrated the relationships were best described by linear equations. The best-performing visible and multispectral VI-based SIs were selected on the basis of mean R2 values over all sampling dates within each group. In order to make hypothetical topdressing decisions in the same manner as the petiole sap nitrate test, i.e. by comparing sampled values to a minimum threshold, minimum SI thresholds for triggering topdress applications were calculated for the selected SIs. Using the regression equations between total yield and the SIs for each date, SI thresholds were calculated as the SIs corresponding to the mean total yield of the reference treatment in each year. Only highly significant (p < 0.01) regression equations were used. The average of these calculated thresholds for each year was then compared to plot SI values, with SI values below the threshold 22 corresponding to a hypothetical topdress event. The percentage of plots below the threshold in each treatment using SI methods as well as the petiole sap nitrate test were compared. Hypothetical topdressing decisions made using petiole nitrate concentration or SI metrics and their corresponding thresholds were also compared. For each plot at each sampling date, two metrics were considered to agree if they both showed that the plot required topdress N (sampled value < threshold; “yes”) or both showed that it did not (sampled value > threshold, “no”). Metrics were considered to disagree if one metric showed that a topdress was required while the other did not. Results were presented using a contingency table to summarize the decisions made using each metric as well as the agreement or disagreement between metrics. 23 RESULTS Weather Rainfall and temperature conditions differed between 2019 and 2020 with several large departures from the 30-year means (Table 1.4). May through August 2020 were 0.3-2.3°C warmer than the same months in 2019, with June 2020 in particular significantly hotter than both June 2019 and the 30-year June mean. Cumulative rainfall in April through October 2019 was particularly high, with 274 mm more rainfall than the same time period in 2020. In addition to being dryer than 2019, 2020 received 97 mm less rainfall than the 30-year average. Center pivot irrigation was managed by the grower at each location using irrigation schedules representative for the region. Table 1.4. Mean temperature and monthly rainfall for Hart, MI in 2019 and 2020, with NOAA 30-year (1981-2010) monthly averages for the same area. Carrot Yield and Quality The 202 kg N ha-1 rate demonstrated significantly higher defect-free yield, total yield, and shoot fresh weight along with significantly lower root:shoot ratio when compared to the lower N rates with the exception of the 2020 defect-free yield (Table 1.5). In that case, no significant defect-free yield differences were found between any N rate. Similarly, in 2020 only the 202 kg ha-1 N rate demonstrated significantly higher total yield than the control. In all other 24 2019202030-Yr2019202030-YrApr6.85.16.7686978May11.912.212.6956696June16.819.117.91118183July21.622.320.3608477Aug18.920.619.61199184Sept17.215.115.31044091Oct8.77.38.92187189RainfallMean Temperature--------- (°C) ------------------ (mm) --------- Table 1.5. Effects of N rate and topdress division on carrot root and shoot weight at harvest. † One treatment within group significantly different (p < 0.10) from low N control treatment. ‡ Two treatments within group significantly different (p < 0.10) from low N control treatment. Within each column, different letters indicate significant differences at p < 0.10 (Tukey HSD). cases, both 135 and 202 kg N ha-1 showed significantly higher defect-free and total yields and shoot fresh weights than the control as well as lower root:shoot ratios. There was no significant effect of topdress division on either defect-free or total yield, but treatments with split N applications consistently demonstrated significantly higher shoot fresh weight and lower root:shoot ratios than treatments with frontloaded N applications. The shapes of the yield and shoot fresh weight responses to season total N rate with the highest coefficient of determination varied between linear and logarithmic (Figure 1.1); second- order polynomials were also tested but no quadratic effects were found. The best fit over the range of data for total yield was logarithmic in 2019 but linear in 2020. Due to split topdresses increasing shoot biomass, treatments with split topdress applications were fitted separately from treatments with single frontloaded applications for shoot fresh weight. The best fit for shoot fresh 25 Low N treatment29 kg ha-169.066.178.772.08.611.39.46.4N Rate main effect67 kg ha-186.4A68.791.9A†75.3A10.4A12.3A9.0C6.2C135 kg ha-185.0A†69.892.8A†77.5A12.8A†15.0B†7.6B†5.2B‡202 kg ha-196.4B‡73.2103.3B‡83.2B‡15.8B‡17.7C‡6.6A‡4.8A‡Divison main effectFront87.7†72.395.9‡79.4†11.7A†14.4A†8.3B†5.6B†Split90.9‡68.896.1‡77.9†14.2B‡15.6B‡7.1A‡5.1A‡N Rate (N)Division (D)N × D0.044NSNSNSNSNSNSNS<0.001<0.001NSNSNSNS0.0170.074<0.0010.0060.026NS0.0310.008<0.001<0.001202020192020--------------------------------- (t ha-1) --------------------------------------- unitless ------------------------------------------------ Significance ------------------------------------------Defect-Free YieldTotal YieldShoot Fresh WtRoot:Shoot Ratio20192020201920202019 Figure 1.1. Lines of best fit for total yield and shoot fresh weight at harvest. Different fits were used for subsets of the data (front + none, split + none) where topdress division effects were found. Fit equations and coefficients of determination (R2) values are shown for corresponding regression fit lines. weight was consistently linear for treatments with the split topdress applications as well as frontloaded applications in both years. Split topdress timing showed no significant effects on defect-free or total yield and was only marginally significant (0.05 < p < 0.10) for shoot fresh weight and root:shoot ratio (Table 1.6). Late timing was associated with increased shoot biomass and decreased root:shoot ratio at Table 1.6. Effects of split topdress timing on carrot root and shoot weight at harvest. Within each column, different letters indicate significant differences at p < 0.10 (Tukey HSD). 26 Timing effectEarly77.670.785.276.610.7A14.3A7.9B5.4AMiddle87.969.493.378.915.0AB16.5AB6.5A4.8ALate87.976.395.882.815.6B17.3B6.3A4.8ATiming (T)0.075--------------------------------- (t ha-1) --------------------------------------- unitless ------------------------------------------------ Significance ------------------------------------------NSNSNSNS0.0850.0870.055201920202019202020192020Defect-Free YieldTotal YieldShoot Fresh WtRoot:Shoot Ratio20192020 harvest relative to the early timing, with the middle timing showing no significant difference from either group in shoot fresh weight. Effects of N rate, number of topdresses (frontloaded or split), and split topdress timing on carrot root cull categories were limited (Table 1.7). A low topdress N rate was associated with a higher proportion of small carrots in 2019, and split applications were associated with more split roots in 2020 compared to frontloaded applications. The highest season total N rate of 202 kg ha-1 also showed significantly higher proportions of split carrots in 2019 when compared to the low-N control treatment. No significant effects of split topdress timing were detected for any cull categories in either year (data not shown). Table 1.7. Effects of N rate and topdress division on carrot root cull categories at harvest. ‡ Two treatments within group significantly different (p < 0.10) from low-N control treatment. Within each column, different letters indicate significant differences at p < 0.10 (Tukey HSD). Nitrogen Uptake and Removal Higher season total N rates were associated with greater N uptake in both years (Table 1.8). Both total plant N uptake and N removal via roots exceeded applied N at 29 and 67 kg N ha-1 with the reverse occurring at 135 and 202 kg N ha-1. These mean differences (net N uptake, 27 Low N treatment29 kg ha-19.32.00.01.41.32.81.60.4N Rate main effect67 kg ha-12.43.10.72.01.1B2.40.50.8135 kg ha-13.73.20.92.80.6A2.20.31.0202 kg ha-12.34.41.8‡3.30.7A2.50.91.6Divison main effectFront3.63.60.61.8A0.82.30.51.0Split2.03.51.73.7B0.72.40.61.3N Rate (N)Division (D)N × DNSNSNS0.020NSNSNSNSNS----------------------------------------- Significance -----------------------------------------NSNSNSNSNSNS0.013NSNSNSNSNSNSNS0.017ForkSplitSmallNub20192020201920202019202020192020-------------------------------------- % harvest by weight -------------------------------------- Table 1.8. Effects of N rate and division on plant N uptake, net N uptake, net N removal, and residual soil N. * Residual soil N (nitrate + ammonium) calculated for a hectare-furrow-slice of 20 cm depth. † One treatment within group significantly different (p < 0.10) from low-N control treatment. ‡ Two treatments within group significantly different (p < 0.10) from low-N control treatment. Within each column, different letters indicate significant differences at p < 0.10 (Tukey HSD). net N removal) were also significantly different from zero in all cases except for net N uptake in 2019 for the 135 kg ha-1 N rate. Like plant N uptake, net N uptake and net N removal increased at higher N rates, with up to 80 kg N ha-1 difference between applied N and N removal. Split applications were associated with lower net N uptake and removal compared to frontloaded applications in 2019 but not 2020. There were also significant interaction effects between the season N rate and topdress division factors for total N uptake, net N uptake, and net N removal in 2019 (Figure 1.2). At 135 and 202 kg ha-1 N, split applications demonstrated higher plant N uptake, lower net N uptake, and lower net N removal compared to frontloaded treatments in 2019. No significant interaction effects were found for 2020 data. Some effects of split topdress timing on 28 Low N treatment29 kg ha-17386-44-57-34-444.35.9N Rate main effect67 kg ha-194A†95A-27A†-27A‡-15A†-12A‡4.66.6135 kg ha-1121B‡113B‡14B‡21B‡31B‡43B‡4.66.9202 kg ha-1166C‡151C‡36C‡50C‡59C‡80C‡5.413.3†Divison main effectFront109A‡118‡26B§16§40B§36§4.98.7Split145B§121‡-10A‡13§10A‡38§4.89.2N Rate (N)Division (D)N × DN Uptake- N Uptake- N RemovalSoil N *----------------------------------------- (kg N ha-1) -----------------------------------------NSNS0.005NS0.005NS0.006NSNSNS<0.001NS<0.001NS<0.001NSNSNS<0.001<0.001<0.001<0.001<0.001<0.001202020192020------------------------------------------ Significance ------------------------------------------Total PlantApplied NApplied NResidual20192020201920202019 Figure 1.2. 2019 effects of interaction effects between N rate and topdress division on N uptake and net N uptake and removal calculations. ** and *** indicate significant differences in marginal means of p < 0.01 and p < 0.001, respectively, at a given N rate. uptake and residual soil N metrics were detected, although most were marginally significant and only present in a single year (Table 1.9). Both net N uptake and net N removal were highest for the early timing in 2019 only, while in 2020 the late timing was associated with greater residual soil N (nitrate plus ammonium) at harvest. Table 1.9. Split topdress timing effects on plant N uptake, net N, and residual soil N. * Residual soil N (nitrate + ammonium) calculated for a hectare-furrow-slice of 20 cm depth. Within each column, different letters indicate significant differences at p < 0.10 (Tukey HSD). 29 Timing effectEarly107A11027B2541B444.56.5AMiddle137B112-3A2318A475.07.1ALate126AB1358AB-128AB284.310.5BTiming (T)0.037------------------------------------------ Significance ------------------------------------------0.062NS0.062NS0.081NSNSTotal PlantApplied NApplied NResidual201920202019----------------------------------------- (kg N ha-1) -----------------------------------------N Uptake- N Uptake- N RemovalSoil N *20202020201920202019 Sufficiency Indices Of the three visible indices tested, GLI showed the highest correlation to yield on average over both the 2019 and 2020 seasons with particularly strong correlations (up to R2 = 0.64) in the latter half of the 2019 season (Figure 1.3). All three multispectral indices showed similar correlations to final yield throughout the first half of the 2019 season and all of the 2020 season, but NDVI outperformed NDRE and GNDVI in the latter half of the 2019 season. GLI and NDVI were the visible and multispectral indices most highly correlated with total yield on average. Overall, VIs showed stronger correlations with final yield in 2019 than in 2020. By contrast, correlations between yield and root fresh weight, shoot dry weight, shoot N concentration, and petiole sap nitrate concentration showed R2 < 0.30 with the exception of sampled shoot dry Figure 1.3. Coefficients of determination between total yield and visible VIs, multispectral VIs, and manual sampling data (petiole nitrate, shoot N percentage, root fresh weight, shoot dry weight) over the 2019 and 2020 seasons. 30 weight, which achieved R2 = 0.34 on 25 September 2019. The highest R2 value between petiole sap nitrate concentration and final yield in either year was 0.15, and significant correlations (p < 0.10) were only detected on three dates: 2 July 2019, 26 August 2020, and 23 September 2020. As the indices most highly correlated with yield on average, the GLI- and NDVI-based SI (GLI-SI, NDVI-SI) values corresponding to the mean 2019 and 2020 reference treatment yields are shown in Table 1.10. GLI-SI values ranged from 0.981 to 1.019 while NDVI-SI values ranged from 0.986 to 1.003, with overall mean values of 1.001 and 0.995, respectively. Within 2019 only, the mean GLI-SI and NDVI-SI values were 1.007 and 0.997 while the 2020 means were 0.989 and 0.990. These within-year means were used as minimum threshold values for the following topdress application decision analysis. Table 1.10. Predicted SI values corresponding to the mean total yield of the reference treatment in 2019 and 2020. NS Regression model not significant at p < 0.01. n.d. No data. Topdress Application Decisions Comparisons of GLI-SI and NDVI-SI values to the corresponding minimum thresholds are shown in Figure 1.4. Overall, the patterns displayed by both the GLI-SI and the NDVI-SI are relatively similar. SI values for the 29 and 67 kg N ha-1 rates were typically close to or above the threshold in July and early August but fall below the threshold in late August and September, 31 GLINDVIGLINDVI06-20NSNS06-18NSNS07-02NSNS07-01NSNS07-16NS0.98607-20NSNS08-011.0190.99707-29NSNS08-191.0181.00008-12NSNS08-291.0160.99508-26NS0.98809-110.9980.99709-100.991n.d.09-250.9851.00109-220.9960.98910-09n.d.1.00310-080.9810.993--- SI (unitless) ---2020--- SI (unitless) ---2019 Figure 1.4. GLI- and NDVI-based sufficiency indices of treatments over the late July through early October sampling period, 2019 and 2020. Error bars show standard error of the treatment mean for each sampling date. Mean minimum thresholds for each year are shown as solid lines. while the higher N treatments remain closer to the threshold for the full season. With a few exceptions, frontloaded topdress treatments start with higher SI values than the split topdress treatments but drop lower over time, with the split topdress treatments remaining higher and closer to the minimum threshold through parts of August and September. Similarly, the late 32 topdress treatments tended to have lower SI values compared to the early topdress treatments in July and early August with the tendency reversing after that. Figure 1.5 shows similar comparisons between the petiole sap nitrate concentrations for different treatments to the minimum acceptable thresholds published by MSU Extension (2011). The majority of samples show concentrations below the recommended minimum, especially late in the season. In some instances, frontloaded topdress treatments at the 135 and 202 kg N ha-1 rates show higher petiole nitrate concentrations than split treatments with the pattern reversing before the end of the season. The early topdress treatments also tended to have higher petiole sap nitrate concentrations than the late treatments when the carrots were below 2-4 cm in root diameter. Figure 1.5. Petiole sap nitrate concentrations for different treatments. The minimum thresholds published by Michigan State University are shown as a solid line in both years for all treatments. In most cases, GLI-SI, NDVI-SI, and the petiole sap nitrate test showed deficiencies more often than N was applied on dates where all three data types were collected (Table 1.11). This held true even at the highest tested N rates. The petiole sap nitrate test indicated the need for 33 Table 1.11. Percentages of plots receiving topdress N during experiment or indicated to receive topdress N through decision support tools. Data aggregated across sampling dates. Only dates with all three data types included: 2 and 16 July, 1, 19, and 29 August, 11 and 25 September 2019; 20 and 29 July, 12 and 26 August, 22/23 September 2020. topdress N more often on average than the GLI-SI or NDVI-SI methods, with GLI-SI and NDVI- SI showing 79% agreement in recommending topdress application (Figure 1.6). The petiole sap nitrate test recommendation matched the GLI-SI and NDVI-SI recommendations in 57% and 64% of cases, respectively, with cases where the petiole sap nitrate test recommended a topdress but GLI-SI and NDVI-SI did not constituting the majority of the disagreement (34% and 28% of cases). Figure 1.6. Contingency table summarizing hypothetical topdress N decisions (Y “yes”; N “no”) indicated by three testing methods (petiole sap nitrate test, GLI-based SI, and NDVI-based SI) for all plots over sampling dates with all three data types (n=432). 34 SeasonTopdressN RateStrategy2019202020192020201920202019202029 kg ha-1None00865596100100100Front1404655718510090Split434071657580100100Front140755579458980Split4340544564558680Front140504543306450Split4340397054556160Early2940617046557185Late43605755717096100GLI-SI67 kg ha-1Actual Applied% plots receiving N----------------- % plots below minimum -----------------135 kg ha-1202 kg ha-1135 kg ha-1NDVI-SIPetiole Sap TestYNYNPetioleY50%34%57%28%Sap TestN9%7%9%7%Y52%14%N7%27%NDVI-SIGLI-SINDVI-SI DISCUSSION Tradeoffs in N Topdress Strategies The responses of shoot biomass to N rate were linear over the data range. The 2020 total yield N response was also linear, indicating that a yield plateau and associated optimal N rate were not found. The yield response was logarithmic in 2019 which indicated decreasing marginal yield returns as N rate increased, but a true yield plateau was not found; the recommended N rate was therefore not found to optimize yield under our study conditions, as had been hypothesized. Despite similarities in these trends, the magnitude of the effects of topdress N rate on shoot biomass and the root:shoot ratio at harvest were more pronounced than on total yield and especially defect-free yield, which showed modest differences associated with the different topdress N rates in only one year. This disparity agrees with the findings of Makries & Warncke (2013) and Westerveld et al. (2006), who generally concluded that carrot shoots are more responsive to N than carrot roots. The weaker yield responses may also be related to the ability of carrots to effectively scavenge N from deep in the soil (Westerveld et al., 2006). The lack of evidence for a yield plateau or optimal N rate in this trial is related to the experimental N rates used. 135 kg ha-1 is the optimal N rate recommended by MSU, so the 202 kg ha-1 was planned as the non-limiting N rate. Such recommendations are essentially an average, however, with evidence for higher optimal N rates of 150-180 kg ha-1 under different field conditions (Hochmuth et al., 1999; Warncke et al., 2004). Higher N rates should be included in future work to elucidate overall carrot plant response at these rates under different field conditions and ensure the inclusion of a non-N limited treatment. This trial also lacked a true zero-N check treatment due to the use of active commercial fields where starter N is used to increase early carrot growth, which is especially key for a high value crop in northern growing 35 regions like Michigan. The lack of a zero-N check and a non-N limited treatment present some limitations to interpretation of carrot N response and decision support in this trial. Compared to N rate, topdress application strategy (number of applications, timing) demonstrated fewer statistically significant effects at harvest. Neither the number of topdresses nor the split topdress timing showed significant effects on yield, and even the increases in the shoot fresh weight associated with later timing and split applications were generally only marginally significant. This again was contrary to the hypotheses, as neither split topdresses nor the typical topdress timing could be said to optimize yield if there were no significant differences between those and alternative strategies. Using less topdress N early and increasing the topdress rate over the season was associated with both increased shoot biomass and increased yield in a study by Colombari et al. (2018). While not a direct comparison, the application of topdress N later in the season showed similar effects to the later split topdress start time observed in this trial, i.e. an increase in mean yield with later topdress start times, although the differences were not statistically significant. Similar to this trial, however, other trials have found that the yields of carrots and other vegetable crops were unaffected by the number of topdress applications used (Kelling et al., 2015; Qu et al., 2019; Sanderson & Ivany, 1997; Smoleń & Sady, 2008). No consistent treatment differences were found among the proportions of the forked, split, small, and nub cull categories. Lower N rates were associated with higher proportions of small roots in 2019 but not 2020; this could be due to the higher 2019 rainfall exacerbating N deficiencies at the low rates. More split carrots were found in treatments using split applications compared to frontloaded applications in 2020, which may indicate that the more consistent supply of N lets the plants reach maturity faster and increases their risk of splitting. Higher N rates also led to average increases in the proportion of split carrot roots across both trial years. 36 While these differences were not statistically significant, the pattern suggests a similar influence on time to maturity and risk of splitting. N rate effects on N uptake, net N uptake, and net N removal were consistent across both years, significantly increasing plant N uptake and leading to positive values for both net N uptake and net N removal. Increases in mean carrot plant N uptake at higher N rates, although non-significant, were also found by Makries & Warncke (2013). However, topdress splitting and timing both showed impacts on N uptake, net N uptake, and net N removal in 2019 only. Frontloaded applications and early split applications in particular led to significantly lower plant N uptake and higher positive net N uptake/removal. This supports the hypothesis that splitting topdresses would increase N uptake and removal and may be partly attributable to leaching losses caused by the wet conditions in that year which were mitigated by the use of split applications and later topdress timings. Frontloading and split topdressing could therefore be equally viable strategies for achieving a given yield at a given N rate in some years, with a single topdress saving time, labor, and financial resources relative to split topdress but also increasing the risk of N loss in years with greater rainfall. In terms of timing, earlier split applications were associated with larger positive differences between applied N and uptake/removal in 2019 while in 2020 the later split applications were associated with significantly more residual soil N at harvest. Both of these findings present risks for N loss through in- or post-season leaching; the middle timing may mitigate these risks and provide an optimal balance by avoiding N applications before plants can uptake it while still providing sufficient time for uptake of split applications prior to harvest. Although some yield differences were not detected between treatments or treatment groups at a statistically significant level, there were average differences which could be 37 economically relevant. Processing carrots are a high-value specialty crop, so small yield increases can lead to economically significant returns and potentially justify increased logistical and management costs. For example, there was an average total yield increase of 8.1 t ha-1 at the 202 kg ha-1 season N rate compared to 135 kg N ha-1. At the average 2019 U.S. processing carrot price of $151 t-1 and a conservative urea cost of $0.50 kg-1, this constitutes an average return of $1223 ha-1 for an increased fertilizer cost of only $73 ha-1 (United States Department of Agriculture, 2020). Split topdresses with the late start time were also associated with an average yield increase of 8.4 t ha-1 compared to the early start time, a comparable but statistically insignificant yield advantage which could nevertheless possibly provide economic benefit. Adaptive N management, either through the established petiole sap nitrate test or another tool such as VI-based sufficiency indices, has the potential to realize some of these economic benefits by allowing for delivery of N closer to the optimal time points. Topdress Decision Making Tools The strength of the relationship between a given metric of plant N status and final yield is a crucial aspect of that metric’s relevance to N management. Comparing the VIs derived from remote sensing data and the manual plant sampling data in their respective correlations with total yield, the VIs generally showed significantly better correlations than shoot dry weight, root fresh weight, shoot N concentration, or petiole nitrate across much of the two trial years. This supports the hypothesis that SI values would correlate more strongly with yield than petiole sap nitrate, although this was mainly evidenced from mid-July through early October 2019. Prior to these dates, treatment differences had likely not developed yet as the plants were small, starter N had been applied at planting, and topdressing for most treatments only began on July 2 of that year. While the correlations with yield for most VIs also rose during those months in 2020, they were 38 weaker than at similar time points in 2019. This is likely related to the fact that yield was less N- responsive in 2020 compared to 2019, and an associated smaller variation between plots may have reduced the variability detectable by the SIs. Differences in soils, irrigation practices, and weather likely contributed to this disparity in yield response. In any case, the relationships found between VI-based SIs and total yield were not useful for identifying N deficiencies earlier than July or August, after the reference treatment had received its topdress N resulting in greater N variability across the trial. This disagrees somewhat with the findings of Makries (2005), who found that carrot shoot reflectance across a variety of wavelengths was highly correlated with available soil N over a broader time period. Still, the timing of this utility may vary depending on the amount of starter fertilizer used and when topdress N is introduced. Starter N is key to avoiding yield losses in many crops grown in Michigan, and the use of active commercial carrot fields for this trial required starter N, but it may also limit the utility of SIs by masking soil N deficiencies until later in the season. Both the depletion of starter N and the application of topdress N to the reference area would heavily influence the development of N status differences in the field, so careful analysis is required for interpreting SI differences as they develop. The use of VI-based metrics for identifying and correcting for N deficiencies in root vegetable crops like processing carrots is also complicated by imbalances in N response between roots and shoots. The results of this and other carrot trials suggest that shoots are more responsive to N than roots; as the aboveground component visible to inspection, differences in the size and appearance of plant shoots are essentially what VI differences show (Makries, 2005). The statistically significant relationships found between these VIs and yield demonstrate that there is some relationship between shoot size and root size, but care must be taken in 39 interpreting differences in shoot appearance (measured directly or indirectly) as potential differences in root yield and basing further topdress N decisions on that interpretation. This relates to a limitation of this study, which was the lack of a true non-N limited reference area and associated yield plateau. The manner in which the roots and shoots develop at these higher N rates, along with commensurate changes in SI values, is impossible to determine from the data collected in this trial. It is therefore not known at what N rates SIs may stop demonstrating useful correlations to yield and/or give misleading recommendations based on shoot changes not reflected by the roots. In order to use a SI to make topdress decisions, threshold values that trigger an application must be determined. In one early study utilizing a handheld chlorophyll meter in maize, the SI threshold for N application was based on a relative yield of 0.95 and defined accordingly as 0.95 (Blackmer & Schepers, 1994). This 0.95 threshold continued to be used for other metrics including remote sensing VIs and in other field crops such as cotton (Blackmer & Schepers, 1995; Chua et al., 2003; Varvel et al., 1997). Alternative thresholds were also used in other crops, such as a 0.90 threshold supported for improved resolution in rice (Hussain et al., 2000). The goal in developing the SI thresholds used in this research was to maximize the hypothetical yield without extending past the bounds of the linear relationships established between the SIs and total yield. The average yields of the high-N, frontloaded reference treatment in each year were selected as the bases for the SI thresholds in order to aim for a high yield while accounting for the varying yields of the two sites, although based on yield results this treatment was not non-N limiting. Though the ideal threshold could theoretically vary over time, particularly with potential luxury consumption and increased shoot growth without commensurate root growth, there was 40 limited variability in the calculated thresholds in this trial. Therefore, the mean of the calculated thresholds in each year was chosen to simplify the use of this overall methodology. Using the reference treatment to normalize both the VI values and the yield goal in this way theoretically results in a minimum threshold value of 1; this was borne out in this trial’s data as the calculated threshold values for GLI-SI and NDVI-SI fluctuated around 1 and their yearly averages were within 1% of that value. Processing carrots also have a higher monetary value relative to N inputs than many field crops in which lower thresholds were developed, so reductions in yield due to insufficient N supplies are an especially negative economic consequence to be avoided in this and other high-value crop systems. A higher minimum threshold works to minimize the risk of yield losses from N stress, and different thresholds may be justified with additional research. By contrast, the petiole sap nitrate concentration was not highly correlated to yield in this trial. This may be related to unknown variations in the time between irrigation events and sampling events, as concentrations assessed from samples taken soon after irrigation could be influenced by higher plant water content or uptake of N from irrigation water compared to sampling prior to irrigation (Qiu et al., 2014). The fact that the sampling schedule was not consistently synchronized with the irrigation schedule is a weakness of this study, as well as a practicality issue for the petiole sap nitrate test more broadly as it may not always be feasible to ensure consistent sampling practices with respect to irrigation. The weaker correlations could also be due to higher variability from luxury consumption or rapid changes related to the timing of topdress fertilization. Still, for the purposes of topdress decision making, the determining factor is simply whether or not the concentration meets the published minimum thresholds. There need not necessarily be a direct correlation between petiole nitrate concentration and yield in order to make useful topdress decisions using this tool. 41 While no treatments in this trial utilized the petiole sap nitrate test (or SI values) to actually trigger topdress applications, an analysis of hypothetical topdressing decisions made using petiole nitrate, GLI-SI, and NDVI-SI was conducted in order to evaluate those decisions relative to one another. The agreement between GLI-SI and NDVI-SI on whether or not to topdress was much greater than the agreement between the petiole sap nitrate test and either metric. Most of the disagreement between the petiole sap nitrate test and the other metrics was in the petiole nitrate recommending a topdress where the other metrics did not. In fact, the petiole sap nitrate test indicated the need for topdress N more often than the GLI-SI or NDVI-SI across all treatments. This was especially true in late September when the petiole sap nitrate test recommended topdressed N across all treatments. Typically, topdressed N is not applied close to harvest in order to minimize the risks of excess N and high root nitrate content (Warncke, 1996). Overall, the GLI-SI and NDVI-SI were more conservative in recommending topdress N applications compared to the petiole sap nitrate test using the calculated SI thresholds. Different SI thresholds would produce different results; at the popular SI threshold of 0.95 the recommendations from the GLI-SI and NDVI-SI would be even more conservative. The management of the reference area could also influence these recommendations. In this trial, large drops in the SI values of some treatments were recorded approximately two weeks after the initial topdress, likely due to the sudden growth of the plants in the high-N reference treatment relative to the other treatments. Growers may also use petiole sap nitrate thresholds different from those published by MSU Extension; using lower thresholds in particular would change the hypothetical decisions made in this trial and bring the recommendations more in line with those based on the drone-based SI metrics. Based on the very high proportion of decision points in 42 which the petiole nitrate concentration indicated a topdress was necessary, especially relative to the more conservative SI metrics, lower petiole nitrate thresholds may indeed be justified. In terms of sampling logistics, VI-based SIs offer several advantages over manual sampling regimes like the petiole sap nitrate test. Remote sensing, and drone-based remote sensing in particular, can improve data collection speed and spatial completeness of data compared to petiole sap nitrate testing or SIs based on manual sampling data like chlorophyll meter readings. The benefit of spatially complete data was seen even in this trial, where small sample sizes and high variability in manually sampled carrot plant data likely led to the poor correlations between that data and total yield compared to the SIs. Concomitant with these advantages is an improvement in the potential scalability of the technology for more advanced precision agriculture applications. Although adoption of variable rate N application technology has generally been slow, especially in specialty vegetable crops, future growers could take advantage of VI-based SIs for use with this technology much better than “lower-resolution” decision making tools like the petiole sap nitrate test. Measuring the concentrations of N in tissue or nitrates in petiole sap also only tell part of the story; VIs are influenced by plant size as well as color and thus indicative of plant growth in addition to N status. In this trial, this integrative component of VIs clearly demonstrated an advantage in predicting yield. However, measuring N status through foliage color and size rather than a direct physical assessment could also introduce bias from diseases or other nutrient deficiencies which have little to do with plant N status. There are also a number of practical issues that constrain the utility of VI-based SIs, particularly in processing carrot production. For example, it has been found in a variety of field crops that using SIs to trigger topdresses can decrease N fertilizer inputs while maintaining yields, effectively improving N use efficiency and reducing input costs (Chua et al., 2003; 43 Hussain et al., 2000; Pinter et al., 2003; Varvel et al., 1997). In processing carrots and other vegetable crops, however, high selling prices can make the risk of yield reduction from underapplying N outweigh the economic and environmental costs of overapplying. The technical demands of data collection using drones and image processing are also substantial, and the initial financial outlay for remote sensing equipment and software may be as well. For example, the drone used in this research cost approximately $1,500. The MicaSense RedEdge-MX multispectral camera and mount cost an additional $6,000 together, and the Pix4Dmapper Professor license cost approximately $2,000. These financial costs could be reduced through the use of open-source software such as OpenDroneMap, QGIS, and R for data processing and analysis, although particularly for photogrammetry, additional expertise may be required and results may not compare well to commercial products. The largest expense was the multispectral camera, however, so the largest cost reduction could come through the exclusive use of visible rather than multispectral indices. In this trial GLI-SI and NDVI-SI were similarly correlated to yield and the hypothetical topdressing recommendations derived from these metrics agreed 79% of the time. It may therefore be feasible to utilize only the drone’s integrated digital camera, although more research would be needed to assess the actual performance of these indices in N management. Even if equipment and software costs are minimized, specialized training, licensing, and skills are required to collect, process, and analyze remote sensing imagery. In research labs it can be feasible to perform these activities internally due to the high availability of relatively cheap skilled labor, but particularly in commercial operations the most practical implementation of this type of management may often be through a dedicated service. 44 CONCLUSIONS In this study, the effects of N rate, topdress splitting, and split topdress start time on a processing carrot crop were evaluated. Similar to other studies, the carrot shoots showed a more significant response to N management strategies than the roots, although at a rate 150% that of the MSU-recommended N rate there was a statistically and economically significant yield advantage. No yield plateau was reached across the N rates included in the study. Topdress splitting and the split topdress start time had no significant impact on yield, but frontloaded applications and late start times increased the risk of N loss to the environment. The possible utility of drone-based monitoring for topdress N decision support was also demonstrated via the correlations between total yield and sufficiency indices based on GLI and NDVI. Although a SI threshold value of approximately 1 was found to be reasonable in this trial, the lack of a true non-N limited reference area necessitates additional research to assess the broader applicability of that finding. This technology requires more equipment and expertise than the petiole sap nitrate test and may be more vulnerable to mis-identifying other nutrient or disease issues as N status issues, but it is also more scalable for variable rate N or other precision agriculture applications. Future research should ensure the presence of true non-N limited reference areas through inclusion of a larger range of N rates, which could also help clarify the nature of carrot yield response plateaus in high-yielding environments. Such trials should also include treatments with N management based on petiole sap nitrate testing and SIs, along with different thresholds for triggering topdress applications to examine the impacts and tradeoffs of these decision support tools more rigorously. 45 REFERENCES 46 REFERENCES Barnes, E. M., Clarke, T. R., Richards, S. E., Colaizzi, P. D., Haberland, J., Kostrzewski, M., Waller, P., Choi, C., Riley, E., Thompson, T., Lascano, R. J., Li, H., & Moran, M. S. (2000). Coincident Detection of Crop Water Stress, Nitrogen Status and Canopy Density Using Ground-Based Multispectral Data. Proceedings of the Fifth International Conference on Precision Agriculture. Blackmer, T. M., & Schepers, J. S. (1994). Techniques for monitoring crop nitrogen status in corn. Communications in Soil Science and Plant Analysis, 25(9–10), 1791–1800. Blackmer, T. M., & Schepers, J. S. (1995). Use of a chlorophyll meter to monitor nitrogen status and schedule fertigation for corn. Journal of Production Agriculture, 8(1), 56–60. Buschmann, C., & Nagel, E. (1993). In vivo spectroscopy and internal optics of leaves as basis for remote sensing of vegetation. International Journal of Remote Sensing, 14(4), 711–722. Chua, T. T., Bronson, K. F., Booker, J. D., Keeling, J. W., Mosier, A. R., Bordovsky, J. P., Lascano, R. J., Green, C. J., & Segarra, E. (2003). In-Season Nitrogen Status Sensing in Irrigated Cotton. Soil Science Society of America Journal, 67(5), 1428–1438. Colombari, L. F., Lanna, N. B. L., Guimarães, L. R. P., & Cardoso, A. I. I. (2018). Production and quality of carrot in function of split application of nitrogen doses in top dressing. Horticultura Brasileira, 36(3), 306–312. Hochmuth, G. J., Brecht, J. K., & Bassett, M. J. (1999). Nitrogen fertilization to maximize carrot yield and quality on a sandy soil. HortScience, 34(4), 641–645. Hothorn, T., Bretz, F., & Westfall, P. (2008). Simultaneous inference in general parametric models. Biometrical Journal, 50(3), 346–363. Hussain, F., Bronson, K. F., Yadvinder-Singh, Bijay-Singh, & Peng, S. (2000). Use of chlorophyll meter sufficiency indices for nitrogen management of irrigated rice in Asia. Agronomy Journal, 92(5), 875–879. Kelling, K. A., Arriaga, F. J., Lowery, B., Jordan, M. O., & Speth, P. E. (2015). Use of Hill Shape with Various Nitrogen Timing Splits to Improve Fertilizer Use Efficiency. American Journal of Potato Research, 92(1), 71–78. Lenth, R. (2020). emmeans: Estimated Marginal Means, aka Least-Squares Means. R Package Version 1.5.2-1. https://cran.r-project.org/package=emmeans Louhaichi, M., Borman, M. M., & Johnson, D. E. (2001). Spatially located platform and aerial photography for documentation of grazing impacts on wheat. Geocarto International, 16(1), 47 65–70. Makries, J. L. (2005). Managing Nitrogen For Quality Carrot Tops Using Remote Sensing. Michigan State University. Makries, J. L., & Warncke, D. D. (2013). Timing Nitrogen Applications for Quality Tops and Healthy Root Production in Carrot. Communications in Soil Science and Plant Analysis, 44(19), 2860–2874. MSU Extension. (2011). Petiole sap nitrate guidelines - MSU Extension. https://www.canr.msu.edu/news/petiole_sap_nitrate_guidelines Mulla, D. J. (2013). Twenty five years of remote sensing in precision agriculture: Key advances and remaining knowledge gaps. Biosystems Engineering, 114(4), 358–371. Oliveira, L. F., Scharf, P. C., Vories, E. D., Drummond, S. T., Dunn, D., Stevens, W. G., Bronson, K. F., Benson, N. R., Hubbard, V. C., & Jones, A. S. (2013). Calibrating Canopy Reflectance Sensors to Predict Optimal Mid-Season Nitrogen Rate for Cotton. Soil Science Society of America Journal, 77(1), 173–183. Peñuelas, J., Gamon, J. A., Fredeen, A. L., Merino, J., & Field, C. B. (1994). Reflectance Indices Associated With Physiological Changes in Nitrogen- and Water-Limited Sunflower Leaves (pp. 135–146). Pinter, P. J., Hatfield, J. L., Schepers, J. S., Barnes, E. M., Moran, M. S., Daughtry, C. S. T., & Upchurch, D. R. (2003). Remote Sensing for Crop Management. Photogrammetric Engineering and Remote Sensing, 69(6), 647–664. Qiu, W., Wang, Z., Huang, C., Chen, B., & Yang, R. (2014). Nitrate accumulation in leafy vegetables and its relationship with water. Journal of Soil Science and Plant Nutrition, 14(4), 761–768. Qu, Z. ming, Qi, X. chao, Wang, J., Chen, Q., & Li, C. liang. (2019). Effects of nitrogen application rate and topdressing times on yield and quality of Chinese cabbage and soil nitrogen dynamics. Environmental Pollutants and Bioavailability, 31(1), 1–8. R Core Team. (2020). R: A Language and Environment for Statistical Computing (3.5.1). R Foundation for Statistical Computing. https://www.r-project.org/ Raun, W. R., Solie, J. B., Johnson, G. V, Stone, M. L., Mullen, R. W., Freeman, K. W., Thomason, W. E., & Lukina, E. V. (2002). Improving Nitrogen Use Efficiency in Cereal Grain Production with Optical Sensing and Variable Rate Application. Agronomy Journal, 94, 815–820. Rouse, J. W., Haas, R. H., Schell, J. A., & Deering, D. W. (1973). Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation. 48 Sanderson, K. R., & Ivany, J. A. (1997). Carrot yield response to nitrogen rate. Journal of Production Agriculture, 10(2), 336–339. Selk, G. E., LeValley, R. C., Highfill, G. A., & Buchanan, D. S. (2004). Comparing the Accuracy of the Cardy Portable Nitrate Meter with Laboratory Analysis of Nitrate Concentrations in Summer Annual Forages (pp. 1–2). Spectrum Technologies, Inc. Smoleń, S., & Sady, W. (2008). Effect of various nitrogen fertilisation and foliar nutrition regimes on carrot (Daucus carota L.) yield. Journal of Horticultural Science and Biotechnology, 83(4), 427–434. Stamatiadis, S., Schepers, J. S., Evangelou, E., Glampedakis, A., Glampedakis, M., Dercas, N., Tsadilas, C., Tserlikakis, N., & Tsadila, E. (2020). Variable-rate application of high spatial resolution can improve cotton N-use efficiency and profitability. Precision Agriculture, 21(3), 695–712. United States Department of Agriculture. (2020). Vegetables 2019 summary. https://downloads.usda.library.cornell.edu/usda-esmis/files/02870v86p/0r967m63g/ sn00bf58x/vegean20.pdf Usha, K., & Singh, B. (2013). Potential applications of remote sensing in horticulture-A review. Scientia Horticulturae, 153, 71–83. Varvel, G. E., Schepers, J. S., & Francis, D. D. (1997). Ability for In-Season Correction of Nitrogen Deficiency in Corn Using Chlorophyll Meters. Soil Science Society of America Journal, 61(4), 1233–1239. Warncke, D. (1996). Soil and plant tissue testing for nitrogen management in carrots. Communications in Soil Science and Plant Analysis, 27(3–4), 597–605. Warncke, D., Dahl, J., & Zandstra, B. (2004). Nutrient Recommendations for Vegetable Crops in Michigan (E2904). In Michigan State University Extension Bulletin (Issue E2934). Westerveld, S. M., McKeown, A. W., McDonald, M. R., & Scott-Dupree, C. D. (2003). Chlorophyll and nitrate meters as nitrogen monitoring tools for selected vegetables in southern Ontario. Acta Horticulturae, 627, 259–266. Westerveld, S. M., McDonald, M. R., & McKeown, A. W. (2007). Establishment of critical sap and soil nitrate concentrations using a Cardy nitrate meter for two carrot cultivars grown on organic and mineral soil. Communications in Soil Science and Plant Analysis, 38(13–14), 1911–1925. Westerveld, S. M., McKeown, A. W., & McDonald, M. R. (2006). Distribution of nitrogen uptake , fibrous roots and nitrogen in the soil profile for fresh-market and processing carrot cultivars. Canadian Journal of Plant Science, 86, 1227–1237. 49 Zarco-Tejada, P. J., Berjón, A., López-Lozano, R., Miller, J. R., Martín, P., Cachorro, V., González, M. R., & De Frutos, A. (2005). Assessing vineyard condition with hyperspectral indices: Leaf and canopy reflectance simulation in a row-structured discontinuous canopy. Remote Sensing of Environment, 99(3), 271–287. 50 CHAPTER III: In-Season Remote Sensing of Growth and N Status in Fresh-Market Tomato ABSTRACT Non-destructive in-season measurements of tomato (Solanum lycopersicum) growth and N status have value in grower decision support. Utilizing variability in fresh-market tomato biophysical characteristics within a two-year cover crop by N fertilizer rate experiment in Holt, Michigan, U.S., the objective of this research was to compare drone-based remote sensing measurements (crop surface model [CSM]-derived height estimates and canopy spectral reflectance-derived vegetation indices [VIs]) to common manual measurements (plant height, leaf tissue N, and SPAD chlorophyll meter readings) as indicators for yield and treatment differences. CSM- derived heights were good estimators of measured heights (R2 = 0.89-0.96) overall, with equal or greater capacity to resolve significant treatment differences and equal or greater correlations with final tomato yields relative to measured heights. Neither SPAD nor VIs were capable of directly estimating foliar N concentration across sampling dates. However, VIs were significantly correlated with foliar N within dates, and generally exhibited equal or greater correlations with yield relative to foliar N >50 days after transplanting in 2019 (data not available for 2018). Foliar N concentration identified more statistically significant differences between treatments than SPAD or VIs. Drone-based remote sensing has the potential to detect relevant in-season tomato plant responses to varying N fertility with less in-field labor and more potential for scalability compared to manual measurements, with possible disadvantages in higher upfront costs and wages for data collection and processing. 51 INTRODUCTION Indirect measurements of plant biomass and N status are often used in agricultural research and management to assess plant growth and nutrient status rapidly, inexpensively, and/or non-destructively. For example, plant height can be used as a non-destructive proxy or as a direct estimator for plant biomass (Anthony et al., 2014; Bendig et al., 2015; Yue et al., 2017). Leaf chlorophyll content relates to plant N status and can be estimated non-destructively with sensors which measure the light reflected from or transmitted through plant leaves (Hunt et al., 2011; Penuelas et al., 1995; Westerveld et al., 2003). The SPAD-502 chlorophyll meter is a portable, relatively inexpensive tool that rapidly generates readings correlated with leaf chlorophyll content (Ling et al., 2011; Sandoval-Villa et al., 2000). Repeated sampling using these and similar metrics has been used to measure plant growth and health over time without destructive sampling of research plants, including use in fresh-market tomato (Solanum lycopersicum) (Dong et al., 2005; Gianquinto et al., 2006; Muñoz-Carpena et al., 2008). SPAD chlorophyll measurements have also been investigated for use in production systems with the goal of assessing plant health or N status for improved decision making at reduced costs (Fontes & de Araujo, 2006; Hussain et al., 2000; Varvel et al., 1997). The labor involved in collecting height and SPAD data can limit the amount of sampling which is practical in a given period or at a given interval, and in research can constrain the number of treatments which can be included in an experimental design. The SPAD meter also has practical limitations related to variability introduced by measurement technique as well as differences in plant cultivar, growth stage, and the presence of other nutrient or water deficiencies (Muñoz-Huerta et al., 2013). High-resolution remote sensing using drones allows for rapid collection of similar data over larger areas, potentially reducing issues related to variability 52 in small sample sizes (Pinter et al., 2003). Structure-from-motion photogrammetry using drone imagery has been studied for use in estimating plant height and biomass in crops such as barley and sugarcane (Bendig et al., 2015; De Souza et al., 2017). Remote sensing of plant reflectance in specific spectral bands can also be used to glean information about plant N status similar to SPAD meters, particularly though the calculation of vegetation indices (VIs) (Mulla, 2013). In measuring reflectance at the canopy level rather than the leaf level, high-resolution VIs such as those calculated from drone imagery integrate characteristics like canopy cover in their measurements in addition to chlorophyll content (Muñoz-Huerta et al., 2013). Using remote sensing techniques to estimate plant biophysical characteristics at the experimental plot level is typically done for large-scale phenotyping research or for precision agriculture applications in field crops (Bendig et al., 2014; Brocks & Bareth, 2018; Holman et al., 2016; Li et al., 2016). Still, some trials in vegetables including tomato have investigated estimating plant height from remote sensing data, either directly using high-resolution digital surface and terrain models or indirectly through VIs (Enciso et al., 2019; Kaplan et al., 2021). These trials found strong correlations between measured plant heights and estimated heights or VI values (R2 > 0.80) which indicates high potential for remote sensing in this area. These trials and others also investigated estimating attributes such as SPAD readings and leaf area index, although results were mixed and platforms other than drones were used for remote sensing (Wu et al., 2014). Analysis of treatment separability or relevance of these attributes to yield is less common, however. Using variability in growth and N status imposed by cover crop and N fertilizer treatments in an ongoing fresh market tomato experiment, one objective of this study was to assess the robustness of structure-from-motion estimates of plant height and remote sensing VIs 53 and as estimators for measured plant height and foliar N concentration, respectively. SPAD measurements were also assessed as possible estimators for foliar N. The second objective was to assess the sensitivity of measured and estimated plant heights, SPAD measurements, VIs, and foliar N concentration to yield and differences between treatment groups. We hypothesized that structure-from-motion estimates of plant height would accurately predict measured heights, and that both SPAD measurements and VIs would correlate strongly with foliar N. We further hypothesized that VIs would resolve more significant treatment differences than SPAD measurements. 54 MATERIALS AND METHODS Field Trial A field trial was initiated from August 2017 to September 2019 in a single research field located at Michigan State University’s Horticulture Teaching and Research Center in Holt, Michigan, USA (42.6714˚ N, 84.4845˚ W). The soil in the study area was a Conover loam. Initial soil chemical characteristics in 2017 included an average pH of 5.4 and 2.1% organic matter with 123 mg kg-1 P (Bray-1 equivalent), 216 mg kg-1 K, 92 mg kg-1 Mg, and 905 mg kg-1 Ca. Treatments were arranged in a split plot randomized complete block design with four replicates. Fall cover crop was the whole-plot factor and spring cash crop N rate was the split- plot factor. In 2017, the field was prepared for fall cover crop planting with three passes of a disc cultivator and cover crops were direct seeded into 7.62 x 7.62 m plots on 10 August. Oilseed radish (Raphanus sativus L. ‘Defender’), rapeseed (Brassica napus L. ‘Dwarf Essex’), and hairy vetch (Vicia villosa Roth) were seeded at 11, 9, and 42 kg ha-1 seed density, respectively, with 13 cm row spacing using a JP series push seeder (Jang Automation Company, Korea). A control treatment received no cover crop. All cover crop treatments were terminated the following spring on 17 May 2018. Each plot was divided into four subplots and a ramp of N rates (0, 56, 112, 168 kg ha-1 N) was applied as urea in a 1.22 m wide band in each subplot on 29 May 2018. Raised beds were formed following these bands, incorporating the urea, and drip irrigation and black plastic mulch were installed. Ten tomato plants (Solanum lycopersicum L. ‘Mountain Fresh’) were manually transplanted from trays into each subplot at 0.46 m spacing on 1 June. One additional guard plant (2018 cultivar Mariana, 2019 cultivar Roma VF) was planted on each end of the subplots. To emulate common management practices in the fresh market tomato industry, stakes were 55 installed with string around every three plants to encourage vertical growth and ease of harvest. This procedure was repeated in an adjacent section of the same field, with cover crops seeded on 16-17 August 2018 and terminated the following spring on 23 May 2019. Fertilizer was applied and beds formed on 31 May and tomatoes were transplanted on 3 June 2019. Data Collection Tomato plant heights were measured using a meter stick throughout the 2018 season on 26 June, 6 July, 13 July, 23 July, and 16 August. Heights were also measured approximately every two weeks during the 2019 season on 14 June, 26 June, 9 July, 23 July, 8 August, 21 August, and 9 September. On 13 July 2018, the youngest mature leaflets from each experimental plant were collected, dried for 7 days at 60˚C, weighed, and ground for tissue nutrient analysis. On 23 July 2018, a Konica Minolta SPAD-502 meter was used to collect SPAD data from the youngest mature leaflet on each experimental plant for a total of ten readings per subplot. In 2019, SPAD readings were taken on 9 July, 23 July, and 8 August and the youngest mature leaflets collected immediately afterward for nutrient analysis. Tomatoes were manually harvested five times in 2018 (24 and 31 August; 7, 17, and 26 September) and six times in 2019 (20 and 28 August; 3, 12, 18, and 26 September). Following each harvest, marketable tomatoes were graded based on the USDA grades U.S. No. 1, U.S. No. 2, U.S. No. 3; unmarketable tomatoes were grouped into a general cull category. A DJI Phantom 4 Pro quadcopter equipped with a 20-megapixel natural color camera was used to collect red-green-blue (RGB) composite aerial imagery of the experimental sites. In 2019, the Phantom 4 RedEdge Mount Kit (Sky Flight Robotics, Inc.) was also affixed to the quadcopter to allow the use of a MicaSense RedEdge-MX camera and the collection of multispectral image data at the 475 nm (blue), 560 nm (green), 668 nm (red), 717 nm (red-edge), 56 and 842 nm (near-infrared) spectral bands. RGB image data was collected weekly throughout both tomato growing seasons, i.e. 8 June to 16 October 2018 and 3 June to 26 October 2019, while multispectral data was only collected during the 2019 growing season. Remote sensing dates coincided with plant sampling dates with the exception of 6 July and 23 July 2018 which were instead flown on 5 July and 25 July 2018, as well as 14 June and 9 September 2019, flown 12 June and 6 September 2019. Remote sensing data was collected at an altitude of 14 m AGL with 85% image front overlap and 75% side overlap for the RGB camera; grid maps defining the parameters for these flights were produced using the DJI GS Pro app. Prior to each flight in which the RedEdge-MX multispectral camera was used, images of a MicaSense Calibrated Reflectance Panel were taken using that camera. Image Processing Orthomosaics were generated from both the RGB camera images and the RedEdge-MX images using the photogrammetry software Pix4Dmapper. For each flight date, an RGB composite orthomosaic was produced from the RGB camera imagery, and during the 2019 growing season five single-band orthomosaics (blue, green, red, red-edge, and near-infrared) were produced from the RedEdge-MX imagery. The composite and single-band orthomosaics had resolutions of approximately 0.005 m and 0.01 m respectively. A digital surface model (DSM) was also generated from the RGB imagery on each flight date with a resolution of approximately 0.005 m. Six ground control points (GCPs) were installed in the first year’s trial area on 8 June 2018 and another six GCPs were installed at the second year’s trial area on 17 August 2018. The UTM coordinates of the GCPs were derived from the RGB composite orthomosaics generated for each site using images taken on the installation dates and used to georeference following orthomosaics. Calibrated reflectance panel images as well as data from 57 the downwelling light sensor of the RedEdge-MX camera were used with Pix4Dmapper to control for ambient light conditions. The open-source geographic information system application QGIS (v3.2) was used to digitize experimental plots. For each experimental area, the ‘Create Grid’ tool was used first to generate a grid of 0.70 m by 4.57 m rectangular polygons; these polygons were of the same width as the raised plastic beds and extended over all experimental plants in each subplot. This shapefile was manually aligned with the first RGB composite orthomosaic at each site using the ‘Advanced Digitizing’ toolbar. Further image processing was performed using R statistical software (R Core Team, 2020). A buffer of 0.5 m was applied to the overall extent of the plot area shapefiles to avoid edge issues during resampling, and this extent was used to resample orthomosaics and DSMs from their native resolutions to standard 0.005 m (RGB composite, DSM) and 0.01 m (single-band) resolutions using the nearest neighbor method. The RGB composite was separated into its constituent red, green, and blue bands, and each band was divided by the sum of the three bands to normalize the individual bands. The normalized RGB and multispectral bands were used to calculate a selection of VIs associated with leaf chlorophyll concentration and canopy cover and, by extension, plant N status (Table 2.1). Three VIs calculated using only the RGB camera data (BGI, NPCI, GLI; “visible indices”) were included alongside three VIs calculated using multispectral camera data (NDVI, NDRE, GNDVI; “multispectral indices) to assess the relative performance of indices in these two groups given the differing equipment required. The visible indices BGI and NPCI have been found to be correlated with chlorophyll and carotenoid concentrations, while the GLI was developed to separate vegetation in imagery based on chlorophyll’s absorption of red and blue light and reflectance of green (Louhaichi et al., 2001; Peñuelas et al., 1994; Zarco-Tejada et al., 2005). 58 NDVI was first developed using satellite imagery to assess large-scale changes in green vegetation, but has since been used to monitor crop canopy development and biomass with a wide variety of sensors and platforms (Mulla, 2013; Rouse et al., 1973). NDRE and GNDVI were subsequently developed to improve estimation of leaf chlorophyll content, with NDRE using differences in red-edge reflectance between healthy and stressed vegetation and GNDVI using changes in chlorophyll’s green reflectance (Barnes et al., 2000; Buschmann & Nagel, 1993). Mean values of these six VIs were calculated over the plot areas defined by the experimental plot shapefiles. Table 2.1. Selected vegetation indices related to plant N status. R indicates measured reflectance with a subscript denoting the wavelength in nm. Subscripts r, g, and b represent the normalized channels of the RGB camera, which does not have well-defined spectral bands. Crop Surface Models Estimates of plant height were derived from crop surface models (CSMs) representing the elevation of the crop surface. DSMs obtained from photogrammetry of the RGB imagery were processed into CSMs by subtracting the ground elevation, isolating the height of the vegetation present in the imagery. Raster maps representing only the ground elevation, also known as digital 59 Index NameAcronymFormulaSourceNormalized PigmentChlorophyll IndexNormalized DifferenceVegetation IndexNormalized DifferenceRed-Edge IndexGreen Leaf IndexBlue-Green IndexGreen NDVIZarco-Tejada et al., 2005Peñuelas et al., 1994Louhaichi et al., 2001Rouse et al., 1973Barnes et al., 2000Buschmann & Nagel, 1993GNDVI(R842 – R560) / (R842 + R560)(R842 – R717) / (R842 + R717)(R842 – R668) / (R842 + R668)(2Rg – Rr – Rb) / (2Rg + Rr + Rb)(Rr – Rb) / (Rr + Rb)BGINPCIGLINDVINDRERb / Rg terrain models (DTMs), were generated using two distinct methods for comparison of the resultant CSMs and plant height estimates. Method 1 is derived from a method utilized by Bendig et al. (2014) in which a DSM is generated from imagery collected immediately prior to planting, i.e. when no vegetation is present, and used as the DTM for the entire season. Although no pre-transplant imagery was collected in either year, an approximation of a pre-transplant DTM was created from imagery collected immediately after transplanting on 8 June 2018 and 3 June 2019, respectively. Pixels representing plants were first identified by using the normalized RGB bands to calculate the Excess Green Minus Excess Red Index (ExGR = 3Rg – 2.4Rr – Rb), a VI used to classify imagery with ExGR > 0 indicating the presence of green vegetation (Meyer & Neto, 2008). For each date, this thresholding rule was applied to the ExGR raster to create a binary raster. In order to reduce excessive processing time, both the binary rasters and DSM rasters were resampled to a resolution of 0.100 m using the nearest neighbor and bilinear interpolation methods, respectively. The binary rasters were then used to mask vegetation pixels from the DSMs, creating spatially incomplete DTMs containing only pixels representing the terrain. In order to characterize the spatial covariance of the data in the incomplete DTMs, a linear variogram was generated for each incomplete DTM from a random sample of 10% of the pixels in each DTM. These variograms were created in R using the automap::autofitVariogram() function and used as arguments in the gstat::krige() function to interpolate the missing pixels via ordinary kriging (Gräler et al., 2016; Hiemstra et al., 2009). The “nmax” argument of the gstat:krige() function was set to 300, allowing only the nearest 300 observations to be used for prediction and substantially reducing processing time. Imagery collected immediately post- 60 transplant was used for this purpose due to the minimal plant size and absence of weeds, requiring minimal pixel removal and interpolation. Method 2 utilized the same kriging procedure on each individual flight date to form a new DTM for that date, with no flight required to image the field at or near transplanting. This method was used to evaluate the applicability of CSM methodologies to singular time points during the growing season as opposed to a sustained data collection campaign, similar to the work done by Holman et al. (2016). Subtracting the DTMs produced by these methods from the DSMs created for each flight date yielded two distinct CSMs. CSM1 was generated using the whole-season DTM created using Method 1 and CSM2 was generated using the individual-date DTMs created using Method 2. The 95th percentile values of CSM1 and CSM2 were calculated for each plot on each sampling date using the digitized plot shapefiles and the raster::extract() function in R (Hijmans, 2020). These values were taken as estimates of tomato plant height based on the work of Anthony et al. (2014), which demonstrated high accuracy in estimating maize plant heights using the 95th percentile elevation values of point clouds similar to those used to generate the DSMs in this trial. Statistical Analysis All statistical analysis was performed using R statistical software (R Core Team, 2020). Due to the detection of significant effects of year in ANOVA tests, the CSM-derived estimates of plant height were regressed against mean measured plant heights separately for each year. Data from all sampling dates and treatments within a given year were included in these regressions to capture the full range of observed plant heights and assess the overall predictive capacity of the CSM-derived plant height estimates. A more granular analysis was also performed, with simple linear regressions between measured and estimated plant heights utilizing only data from within 61 individual sampling dates to assess changes in predictive capacity over time. Similarly, the capacity for the SPAD readings and VI values to predict foliar N concentration was evaluated using simple linear regression with data from individual sampling dates as well as data pooled within year. Visual inspection of the data was used to assess the possibility of higher order effects. In order to evaluate the capacity of each metric (measured plant heights, CSM1/CSM2 estimated plant heights, SPAD readings, VI values, foliar N concentration) to distinguish treatment differences, marginal means were first computed from mixed-effects models utilizing each metric as a response variable (Hothorn et al., 2008). Cover crop and the cash crop N rate were considered as fixed factors in these models, as was the interaction of the two factors. Replicate and cover crop were included as the random main plot and sub-plot factors, respectively. Models were evaluated using data from each sampling date separately. When heteroskedasticity was detected for one or more factors, weighted mixed-effects modelling was used with unequal variance weights for each level of the relevant factors. Using the marginal means computed for each fixed factor level in the models as well as for each individual treatment, pairwise comparisons between these means was performed with Tukey’s HSD (Lenth, 2020). Pairwise comparisons with p < 0.1 were considered to show statistically significant differences between factor levels or treatments. The relative abilities of the response variables to detect significant treatment differences was evaluated by comparing the significant or non-significant determinations made for each pairwise comparison. Two metrics were considered to agree if they both showed a given pairwise comparison as significant or non- significant. They were considered to disagree if one showed a significant difference and the other did not. For plant height, the capacity for measured height to detect significant differences was 62 compared in this way against the capacity of plant height estimates derived from CSM1 and CSM2. For plant N status, the foliar N concentration was compared to the SPAD readings and VI values. Due to a lack of foliar N concentration data, pairwise comparisons for metrics of N status were not made for the 23 July 2018 sampling date. SPAD data and multispectral VI data (NDVI, NDRE, GNDVI) was also not collected on 13 July 2018; pairwise comparisons for these response variables were made using 2019 data only. Two in-field issues also led to the complete exclusion of some observations as outliers. In one plot, three transplants died shortly after planting and were replaced on 11 June 2018. The new transplants remained significantly smaller than the surrounding plants for several weeks, including the 26 June, 6 July, and 13 July 2018 sampling dates. Lodging caused by stake breakage in another plot on 16 August 2018 also artificially biased measurements taken in that plot on that date. The observations from these plots were excluded from the analysis as outliers on the relevant dates. 63 RESULTS Plant Height Estimates of plant height derived from both CSM1 (single initial DTM) and CSM2 (unique DTM per drone flight) were highly correlated (R2 = 0.89-0.96) with mean measured plant heights for all sampling dates within a given year combined (Figure 2.1). Height estimates derived from CSM1 underestimated measured plant height by greater than 12 cm early in the 2018 season, while CSM2 underestimated measured heights by ≤ 5 cm below 50 cm heights and overestimated by ≤ 8 cm above 50 cm. The fitted 2019 equation for CSM1 estimates showed the slope closest to 1 at 1.0 and the intercept closest to 0 at -3.6 making this fit the closest to a 1:1 estimation. Fitted slopes for both estimators were closest to 1 in 2019, with higher slopes of 1.2 for CSM1 and 1.4 for CSM2 in 2018. Within each year, CSM1 slopes were closer to 1 than CSM2 slopes with similar or lower-magnitude intercepts. Figure 2.1. Mean measured plant heights compared to CSM1 (single initial DTM) and CSM2 (unique DTM per drone flight) estimates of plant height with all data from each year combined (2018 n=316; 2019 n=448). The dotted line represents a 1:1 ratio, and the solid line represents the ordinary least squares fit of the data. 64 All correlations between measured and estimated plant heights by sampling date were statistically significant (p < 0.05), although the associated R2 values tended to be less than 0.40 in June and early July with particularly low R2 values for CSM2 estimates (Figure 2.2). The exception was the CSM1 estimates on 6 July 2018 with R2 = 0.53, and stronger correlations later in the season reached R2 values of up to 0.80. Correlations between measured heights and CSM1 estimates showed higher R2 values than between measured heights and CSM2 estimates on all sampling dates with the exception of 23 July 2018. Figure 2.2. R2 values for correlations between measured plant heights and estimates of plant height derived from CSM1 (single initial DTM) and CSM2 (unique DTM per drone flight) for each individual sampling date in 2018 and 2019. CSM1 and CSM2 showed similar performance relative to the measured plant heights in identifying significant differences between factor levels and between all treatments in pairwise comparisons (Figure 2.3). Using measured plant heights as the response variable, 40% of pairwise comparisons between levels of cash crop N rate found statistically significant (p < 0.1) differences while 14% of comparisons between cover crop levels found significant differences. For pairwise comparisons between all treatments, 4% showed significant differences. The same sets of pairwise comparisons (between cover crops, N rates, or all treatments) made using CSM1 estimates of plant height as the response variable agreed with the results of measured heights, i.e. both metrics found any given comparison to be significant or non-significant, in 85-93% of 65 Figure 2.3. Contingency tables summarizing the significance of pairwise comparisons of marginal means calculated from mixed-effects models with different response variables (measured plant heights vs. CSM1, CSM2 height estimates). Pairwise comparisons were performed between marginal means computed on each date for levels of cover crop (n=72), levels of cash crop N rate (n=72), and all individual treatments (n=1440). cases. Similarly, the significance or non-significance of pairwise comparisons using measured heights agreed with those using CSM2 height estimates in 79-96% of cases. Discrepancies between the results using measured and estimated plant heights more often took the form of significant differences found using estimated plant heights which were not significant using measured heights. Plant N Status Foliar N data were sampled once in 2018, sharing a sampling date (13 July) with a drone flight to collect RGB imagery and the associated visible indices data. In simple linear regressions for this date, the NPCI correlated best with foliar N concentration (R2 = 0.32) while BGI and GLI both demonstrated less significant correlations (R2 < 0.15). No data were available to assess the correlations between foliar N and other data types (SPAD, multispectral indices) in 2018. Foliar N, SPAD, and both visible and multispectral index data were collected three times in 2019 (9 July, 23 July, 8 August). When 2019 data were combined, NPCI was best correlated with foliar N (R2 = 0.14) while all other indices and SPAD were not well-correlated (R2 < 0.04). No 66 NSp < 0.1NSp < 0.1NSp < 0.1NS75.0%6.9%47.2%2.8%89.9%0.6%p < 0.111.1%6.9%12.5%37.5%6.5%3.1%NS75.0%9.7%47.2%8.3%93.0%0.8%p < 0.111.1%4.2%12.5%31.9%3.3%2.9%Cover CropCash Crop NAll TreatmentsCSM2MeasuredMeasuredMeasuredCSM1 evidence of higher order effects was found when 2019 data were pooled in this way. Foliar N concentrations ranged between 2.18% and 4.38% in 2018 with a similar range between 2.10- 4.92% in 2019. When linear regressions used only data within the three individual 2019 sampling dates, SPAD readings showed weak correlations to foliar N concentration (R2 < 0.10) while both the visible and multispectral indices showed stronger correlations (Figure 2.4). Correlations between foliar N concentration and VIs showed the highest R2 values on 9 July and 23 July 2019 with weaker correlations on 8 August with the exception of NPCI. NPCI was the visible index with the strongest linear correlation to foliar N (mean R2 = 0.44), with NDRE showing the strongest correlation of the multispectral indices (mean R2 = 0.52). Some evidence for a quadratic fit was found for the relationship between foliar N and GNDVI on 9 July 2019. Adding a quadratic term to the simple linear model on that sampling date resulted in R2 = 0.60 (Foliar N = -31 + 100*GNDVI – 71*GNDVI2) compared to the linear fit R2 = 0.58 (Foliar N = -5 + 15*GNDVI). Figure 2.4. R2 values for linear correlations between foliar N concentration and SPAD readings, visible indices, and multispectral indices for individual sampling dates. Both foliar N data and any other data type (in this case, visible indices) were collected on only one date in 2018; three dates with shared data are available and shown for all data types in 2019. 67 In detecting significant differences in pairwise comparisons between factor levels or all treatments, and particularly in the proportion of agreement with the results of pairwise comparisons using foliar N as the response variable, the performance of SPAD and VIs was varied. Using foliar N, significant differences between levels of cash crop N rate were found in 88% of cases overall (Figure 2.5). 38% of pairwise comparisons between cover crop factor levels were statistically significant, and 36% of pairwise comparisons between all treatments were significant. The results of pairwise comparisons between cover crop levels using other response variables agreed with those made using foliar N in 56% of cases using SPAD, 68% of cases on average using the tested visible indices, and 69% of cases on average using the tested multispectral indices. There was a larger discrepancy between the performance of SPAD and VIs in pairwise comparisons made between levels of the cash crop N rate factor. The significance or non-significance of SPAD comparisons between cash crop N levels agreed with the results of the foliar N comparisons in 28% of cases, while 42-63% of comparisons using visible indices and 50-61% of comparisons using multispectral indices agreed with determinations of significance or non-significance made using foliar N. Of the three groups of comparisons (between cover crops, between N rate factor levels, or between all treatments), agreement was highest at 65-73% for pairwise comparisons of all treatments for all metrics except BGI and GLI. A mean of 83% of the agreement between metrics for these pairwise comparisons was in finding no significant differences between treatments. Consistently across all metrics and groups of comparisons, the largest source of disagreement was significant treatment differences found using foliar N concentration as the response variable which were non-significant using other response variables. In the case of the comparisons using SPAD as well as cover crop comparisons using BGI and NPCI, no significant 68 Figure 2.5. Contingency tables summarizing the significance of Tukey-adjusted pairwise comparisons of marginal means calculated from mixed-effects models with different response variables (foliar N concentration vs. SPAD, visible indices, and multispectral indices). Pairwise comparisons were performed over levels of the cover crop and cash crop N rate factors (SPAD n=18; visible indices n=24; multispectral indices n=18) as well as the interactions (SPAD n=360; visible indices n=480; multispectral indices n=360). 69 NSp < 0.1NSp < 0.1NSp < 0.1NS55.6%44.4%11.1%72.2%65.3%34.7%p < 0.10.0%0.0%0.0%16.7%0.0%0.0%NS62.5%33.3%4.2%50.0%63.5%35.8%p < 0.10.0%4.2%8.3%37.5%0.2%0.4%NS62.5%33.3%4.2%29.2%61.7%27.5%p < 0.10.0%4.2%8.3%58.3%2.1%8.8%NS54.2%20.8%8.3%50.0%60.8%30.4%p < 0.18.3%16.7%4.2%37.5%2.9%5.8%NS44.4%22.2%0.0%38.9%60.8%22.8%p < 0.111.1%22.2%11.1%50.0%4.4%11.9%NS38.9%11.1%0.0%27.8%50.6%12.2%p < 0.116.7%33.3%11.1%61.1%14.7%22.5%NS38.9%16.7%0.0%27.8%56.7%19.4%p < 0.116.7%27.8%11.1%61.1%8.6%15.3%GNDVIVisible IndicesMutispectral IndicesGLINDVINDREManualSPADBGINPCICover CropCash Crop NAll TreatmentsFoliar N %Foliar N %Foliar N % differences were found at all which were not found using foliar N. This apparent under-detection of significant differences relative to foliar N could be large; on the median, significant differences found using only foliar N outnumbered differences found using only other metrics by a factor of 3 (excluding division by zero). The metric which showed the highest average agreement with foliar N (69%), NDRE, was also the only metric of those tested for which this pattern was reversed for pairwise comparisons of all treatments. Relationships to Total Yield Both measured and estimated plant heights were significantly correlated with total yield throughout 2018 and 2019 with one exception on 14 June 2019 when none were correlated with yield (Figure 2.6). Other exceptions occurred with measured plant heights on 26 June 2018 as well as with CSM2-based height estimates on 13 July 2018. Outside of these instances, however, R2 values for statistically significant correlations between yield and plant heights varied between 0.08 and 0.60 over the course of the season. These R2 values were highest in late July and August, with peaks in August for all three metrics in both years. The weakest correlations occurred in June and September. Within each date, CSM1 height estimates demonstrated a stronger correlation to total yield compared to measured plant heights across both years. Figure 2.6. R2 values for linear correlations between total tomato yield and plant heights (measured, CSM1 estimated, CSM2 estimated) in 2018 and 2019. 70 Foliar N concentration and both visible and multispectral indices showed significant correlations with total yield (Figure 2.7). The BGI was not statistically significantly correlated with total yield on either of the 2018 sampling dates. Foliar N concentration was most correlated with yield (R2 = 0.65) on 9 July 2019, when the R2 values for the visible and multispectral index correlations were 0.13-0.25 and 0.30-0.41, respectively. The R2 value for the correlation between foliar N and yield dropped over the two subsequent sampling dates which was opposite the trend observed for the VIs. By 8 August 2019, GLI showed the highest R2 value of the visible indices at 0.56 while NDRE outperformed the other multispectral indices with R2 = 0.74. SPAD showed the weakest average relationship (R2 < 0.12) with final yield across all five sampling dates. Figure 2.7. Coefficients of determination between total tomato yield and foliar N % (measured leaf N content), SPAD measurements (based on leaf spectral transmittance) and visible and multispectral indices from remote sensing (canopy-level averages) for individual sampling dates. Multispectral data available in 2019 only. Significant quadratic terms were found in the relationships between final yield, measured in tonnes per hectare, and BGI as well as GLI on 9 July 2019. Adding a quadratic term to the simple linear model for BGI on that sampling date gave R2 = 0.20 (Yield = -2846 + 6586*BGI – 71 3748*BGI2) compared to the linear fit with R2 = 0.13 (Yield = -202 + 175*BGI). Adding a quadratic term to the GLI model gave R2 = 0.31 (Yield = -148 + 2928*GLI – 10939*GLI2) compared to the linear fit with R2 = 0.25 (Yield = -7 + 430*GLI). A significant quadratic term was also found for the relationship between SPAD and final yield on 23 July 2019, when adding a quadratic term increased the R2 value to 0.21 (Yield = -228 + 7*SPAD – 0.04*SPAD2) from the linear R2 value of 0.11 (Yield = 11 + 0.4*SPAD). 72 DISCUSSION Accuracy of Tomato Plant Height Estimations Height estimates derived from CSM1 (single initial DTM) and CSM2 (unique DTM per flight date) were both highly correlated (R2 ≥ 0.89) with mean measured plant heights with regression slopes close to 1. There was also a high degree of agreement (80-95%) in the significant and non-significant differences for measured and estimated heights in pairwise comparisons for cash crop N and cover crop levels as well as between all treatments. These factors support the hypothesis that height estimates derived from both CSM1 and CSM2 are effective estimators for mean measured plant height. The estimated heights also demonstrated slightly higher R2 values in correlations with total yield across most sampling dates compared with the measured heights, and the main source of disagreement in the pairwise comparisons took the form of significant differences detected using estimated heights which were not significant using measured heights. The height estimates, particularly those derived from CSM1, may therefore have been more informative and sensitive to treatment differences compared to the manually measured heights. In any case, although the direct estimation of measured heights was imperfect, treatment differences were approximately as detectable using the estimates as using the manually collected data. The results of the linear regressions between estimated and measured plant heights are in line with those of other studies. Using methods similar to CSM1, i.e. subtracting a single pre- or post-season bare ground DTM from a DSM generated via structure-from-motion photogrammetry, height estimates of barley derived from CSMs were found to correlate well with measured plant heights with R2 = 0.92 (Bendig et al., 2014). Similar work in wheat yielded R2 = 0.94-0.95, while using the elevation of bare buffer zones to generate new DTMs for each 73 flight date (similar to the CSM2 method) gave wheat height estimates which correlated with measured heights at R2 = 0.93 (Holman et al., 2016). This study also found significant differences in plant height between treatments with different levels of applied N fertilizer, supporting the patterns found in this tomato trial. Other studies utilizing terrestrial laser scanning or other techniques to generate CSMs tended to have similar results, with R2 values for correlations between estimated and measured plant heights between 0.7 and 0.99 (Bareth et al., 2016; Bendig et al., 2013, 2015; Li et al., 2016). Additionally, using a method very similar to CSM2 (interpolating a DTM using non-vegetative area determined using the ExGR index), Geipel et al. (2014) investigated CSM height estimates for predicting yield in maize. Similar to the results of this work, they found that the relationship of the height estimates to final yield was significant in most cases, reached R2 values of up to 0.74, and varied with the growth stage of the crop. While there are clear differences in morphology and production practices between these grain crops and a specialty vegetable crop like fresh-market tomato, the similarity of these results indicates that the utility of remote sensing to meaningfully estimate plant heights is relatively consistent between these cropping systems. The differences between the CSM1 and CSM2 methodologies may account for the differences in their estimations of measured plant heights. A consistent underestimation of measured heights was observed for CSM1 with a particularly large magnitude in 2018, possibly indicating that using a higher percentile value would more precisely align with the manual measurements. The overestimation of measured height late in the season associated with CSM2 estimates may have been due to the use of raised beds and plant growth over time. While the height of these beds is accounted for in the single DTM used to derive CSM1, a new DTM is created using interpolation on every sampling date for CSM2. As the plants grew large enough to 74 cover the raised beds, the interpolation was likely unable to account for the beds and biased downward by the surrounding ground, increasing estimated height. Using a single DTM (the method associated with CSM1) is therefore likely preferable as it avoids this late season bias. For instances where this is not feasible, however, the outlined interpolation technique is still sufficient to deliver estimations with a similar capacity to detect treatment differences. Comparison of Methods for Assessing Plant N Status Foliar N concentration is the most direct measure of plant N status in this study and found the largest proportion of significant differences in pairwise comparisons between levels of the cover crop and cash crop N rate. The comparison between foliar N and SPAD is the most direct as both obtain measurements at the leaf level whereas the VIs measure reflectance at the canopy level and integrate additional information including size and spectral characteristics (i.e. color) of the plants (Muñoz-Huerta et al., 2013). In addition, foliar N is likely to be an earlier indicator of plant N status because characteristics measurable by proximal or remote sensing generally develop in response to plant N levels. Given the high performance in detecting treatment differences, foliar N was considered as a benchmark against which SPAD and VI metrics were assessed. Contrary to the hypothesis, neither SPAD nor any of the tested VIs were suitable for estimating foliar N concentration across more than a single date. These metrics may be more appropriate for assessing relative differences in foliar N concentration at a given point in time rather than for directly estimating the quantity. However, this finding represents data largely from only one year. All six tested VIs showed statistically significant relationships to foliar N on all dates with shared data with the exception of BGI on 8 August 2019 with R2 < 0.70. In other work investigating correlations between VIs and foliar N or chlorophyll concentration, the 75 correlations for well-performing indices commonly demonstrated R2 values in the range of 0.50- 0.80 (Croft et al., 2014; Hunt et al., 2011, 2012; Main et al., 2011). This range generally agrees with the best R2 values shown in this research, though the hypothesis that VIs would correlate strongly with foliar N is not supportable by the low R2 values found for the BGI and NPCI visible indices nor for the 8 August 2019 sampling date when R2 < 0.30 for most indices. Drone- based VIs may therefore represent a more rapid and repeatable method for monitoring relative differences in foliar N concentration and overall N status but could not replace direct foliar N assessments. SPAD was not shown to be appropriate even for assessing relative differences in this trial. Across the three dates with shared data, SPAD showed the weakest relationship to foliar N with only one correlation where p < 0.10 (R2 = 0.07). This overall lack of relationship between foliar N and SPAD readings conflicts somewhat with the literature; Ulissi et al. (2011) did find that visible-near infrared spectral reflectance data similar to VIs showed higher R2 values in correlations with foliar N compared to SPAD readings in tomato, but the lower R2 associated with the SPAD readings was 0.56 which is still substantially higher than those found in this trial. Also, while the strength of correlations between SPAD and foliar N varied between vegetable crops and growth stage, Westerveld et al. (2003) found some statistically significant correlations with R2 values typically in the 0.25-0.95 range for onion, carrot, and cabbage crops. Still, the results of this trial do not support the hypothesis that SPAD measurements would correlate strongly with foliar N. It is possible that chlorophyll saturation occurred in many plots, particularly those with higher N rates or N-rich cover crops, which could have reduced the variation in SPAD measurements such that correlations with foliar N were not significant (Muñoz-Huerta et al., 76 2013). This was evidenced by relatively high foliar N concentrations; 89% and 80% of samples were above the level considered sufficient for tomato leaf N in 2018 and 2019, respectively (University of Florida IFAS Extension, 2015). As a leaf-level measurement, SPAD may also saturate more easily than plot mean VIs, which are influenced by growth responses in addition to canopy spectral changes. In detecting treatment differences via pairwise comparison, both the selected VIs and the SPAD readings showed large departures from the significant differences found using foliar N concentration. For all pairwise comparisons, multispectral VIs typically had the highest agreement with foliar N (average 73% agreement) compared to visible VIs (average 67% agreement) or SPAD (65% agreement). In general, this disagreement came from significant differences in foliar N concentration which were not found with the other metrics. Given that foliar N was more correlated to yield than other metrics earlier in the season, it is likely that foliar N was also more sensitive to early treatment differences as changes in foliar N would precede later changes in spectral signature and biomass growth to which VIs would be more sensitive. In other studies across a range of crops including tomato, SPAD and VIs have both shown significant differences between N rate treatments, particularly later in the season, which also lends support to this possible pattern (Bohman et al., 2019; Bullock & Anderson, 1998; Hunt et al., 2018; Ihuoma & Madramootoo, 2020; Li et al., 2014; Minotti et al., 1994; Sandoval- Villa et al., 2000; Wu et al., 2007; Yang et al., 2014). VIs may therefore be more useful than SPAD for detecting significant differences, supporting that hypothesis, but foliar N concentrations still showed more significant differences between treatments than the VIs or SPAD. Drone based sampling could then be useful in rapidly and frequently collecting data related to N status given its lower labor requirements and per- 77 sampling-date costs in addition to sampling for foliar N at key time points, or as a possible replacement for SPAD sampling. Due to the stronger correlations between foliar N and yield early in the season relative to VIs, this would be most appropriate later in the season as foliar N is likely more useful for determining early-season differences. Research actually utilizing these tools for management would be required to assess whether potential delays in detecting N status differences using VIs relative to foliar N sampling would reduce the ability of growers to correct N deficiencies without impacting yield. Utility of Remote Sensing Metrics in Research and Management Both plant heights and measures of N status showed high correlations with final yield at some sampled time points, which supports the use of these metrics as useful indicators of plant health and future yield. However, there were some large discrepancies between sample data types in the strength of the relationships with total yield. Between foliar N concentration, SPAD readings, and plot mean VIs, SPAD consistently showed the weakest relationship with final yield (R2 < 0.12). Other work in tomatoes supports a stronger relationship between SPAD and yield; Gianquinto et al. (2006) found significant correlations between SPAD readings and yield (R2 = 0.30-0.75) on 19 sampling dates spanning 2 months in a variety of field conditions. Ultimately, SPAD readings indirectly estimate chlorophyll concentration, a metric associated with plant N status, while foliar N concentration is a more direct physical quantity and plot mean VI values integrate both the color and size of the plants for a more holistic assessment of plant growth and nutrient status. These advantages likely contributed to the stronger relationship between the foliar N concentration and the plot mean VI values compared to mean SPAD readings, along with the better correlations between these metrics and yield. Integrative 78 measures like the VIs may therefore be more valuable for assessing treatment outcomes in research environments compared to leaf-level spectral metrics like SPAD. The relationships between all plant height and N status metrics and total yield were much weaker in 2018 than 2019, likely due to poor field conditions in that year. There were more broken stakes and lodging which needed to be remedied in that year, adding a factor which may have impacted yield in a way not captured by sampling on limited dates. High between-row weed pressure in 2018 may also have impacted yield as well as remote sensing metrics. While the plot polygons were defined using the width of the raised beds, large weeds may have extended from the between-row soil to a position above the beds and biased the height estimates or mean VI values. Field maintenance is therefore key to acquiring high-quality remote sensing data free of biases. This is especially crucial in research environments where experimental fields and sample sizes are relatively small; in production, large fields with many sample plants would likely reduce the biases from isolated instances of lodging or small patches of weeds. The contrast in labor requirements for the differing sampling regimes is stark, although the comparison is complicated by data processing and disparities in skill level. Drone flights took one person approximately 15 minutes per flight to complete on average, compared with three people working for 90 minutes for manual height data collection. Collection of SPAD and foliar N data together also took two people approximately 90 minutes per sampling event in this trial. In total, 0.25 in-field labor-hours were required for a drone flight to collect height, biomass, and relative N status data compared to 7.50 in-field labor-hours to collect data manually. No accounting was taken for the time required to process each of these data types, but the expertise required for image processing and analysis is significantly higher than for drying and grinding biomass samples or data entry. The associated difference in wages may balance out the differing 79 labor times, although the importance of this difference depends upon whether time or money is a more limiting factor. More frequent sampling could also be possible using drone-based remote sensing compared to manual data collection depending on labor constraints and when data is needed. Equipment costs also contribute to the differences in resource requirements between manual sampling and remote sensing using drones. The initial investment for sampling via drone imagery, including the costs of a drone, photogrammetry software, and a multispectral camera, could be upwards of $10,000 in total. In this trial, the RedEdge-MX multispectral camera was the greatest of these expenses at $5,500. This camera is unnecessary for creating CSMs due to its lower resolution compared to the Phantom 4 Pro’s camera, so its utility in this trial was generating multispectral indices for assessing relative N status. These indices had a slight advantage over visible indices, showing more significant treatment differences and demonstrating higher agreement with differences found using foliar N on average. This advantage must be weighed against the camera’s high cost, however; relying only on the drone RGB camera for CSM and relevant index measurements and using a subscription service or open-source software for photogrammetry analysis could reduce up-front costs to below $2,000. For comparison, a SPAD (or similar) meter may cost between $2,000 and $3,000. While no special equipment is required outright for sampling leaf tissue for foliar N analysis, lab analysis fees may range from $10 to $20 per sample. This trial required processing of 64 samples per sampling date, translating to approximately $1,000 in lab analysis costs per sampling date. In general, sampling using multispectral drone imagery likely has higher initial costs and expertise required compared to manual sampling. However, significantly higher in-field labor-time requirements and ongoing processing costs associated with manual sampling regimes could 80 make the investment worthwhile, particularly if availability, cost, and ease of imagery collection and analysis decrease further in the future. 81 CONCLUSIONS Estimates of tomato plant height derived from crop surface models were effective estimators of measured plant height while the ability of SPAD or selected vegetation indices to predict foliar N concentration varied by sampling date. Within a sampling date, VIs were significantly correlated with foliar N concentration thus making them appropriate for use in assessing relative plant N status despite variation between dates reducing utility for directly predicting foliar N concentration. Both plant heights and the tested N status significantly correlated with final yield across multiple sampling dates reinforcing potential usage in research applications as indicators of plant growth, N status, and future yield. Treatment differences assessed using pairwise comparison of plant heights showed that plant heights estimated from CSMs and manually measured plant heights were similarly sensitive to these differences, with up to 96% agreement on the statistical significance of pairwise comparisons’ difference from zero. Foliar N concentration showed many more statistically significant differences via pairwise comparison compared to SPAD or the VIs, with SPAD in particular picking up few differences compared to foliar N. VIs may therefore represent an improved non-destructive sampling method compared to SPAD, but collecting foliar N data when practical remains an important method of gauging plant N status. The long-term costs of data collection could be comparable or even reduced using drone- based remote sensing compared to the manual data collection regimes used in this trial, providing a possible advantage to researchers along with reduced in-field labor. Another key research advantage of remote sensing in general is the creation of a rich stored record of plant data; other feature extractions and analyses (additional VIs, fruit counts, etc.) can be performed or tested using stored remote sensing data in parallel with planned goals. The scalability of 82 drone-based remote sensing may also support larger trials or higher-resolution decision support (i.e. precision agriculture applications) relative to manual sampling. The treatments involved in this study were the cover crop preceding the tomato crop and the N fertilizer rate applied to the tomato crop, but future work could involve integrating this type of analysis into research with different factors to assess the sensitivity of the metrics tested in this study to other factors. These methods could also be tested in other vegetable cropping systems, particularly those in which plant height can serve as a useful proxy for biomass. 83 REFERENCES 84 REFERENCES Anthony, D., Elbaum, S., Lorenz, A., & Detweiler, C. (2014). On crop height estimation with UAVs. IEEE International Conference on Intelligent Robots and Systems, Iros, 4805–4812. Bareth, G., Bendig, J., Tilly, N., Hoffmeister, D., Aasen, H., & Bolten, A. (2016). A comparison of UAV- and TLS-derived plant height for crop monitoring: Using polygon grids for the analysis of crop surface models (CSMs). Photogrammetrie, Fernerkundung, Geoinformation, 2016(2), 85–94. Barnes, E. M., Clarke, T. R., Richards, S. E., Colaizzi, P. D., Haberland, J., Kostrzewski, M., Waller, P., Choi, C., Riley, E., Thompson, T., Lascano, R. J., Li, H., & Moran, M. S. (2000). Coincident Detection of Crop Water Stress, Nitrogen Status and Canopy Density Using Ground-Based Multispectral Data. Proceedings of the Fifth International Conference on Precision Agriculture. Bendig, J., Bolten, A., & Bareth, G. (2013). UAV-based Imaging for Multi-Temporal, very high Resolution Crop Surface Models to monitor Crop Growth Variability. Photogrammetrie, Fernerkundung, Geoinformation, 2013(6), 551–562. Bendig, J., Bolten, A., Bennertz, S., Broscheit, J., Eichfuss, S., & Bareth, G. (2014). Estimating biomass of barley using crop surface models (CSMs) derived from UAV-based RGB imaging. Remote Sensing. Bendig, J., Yu, K., Aasen, H., Bolten, A., Bennertz, S., Broscheit, J., Gnyp, M. L., & Bareth, G. (2015). Combining UAV-based plant height from crop surface models, visible, and near infrared vegetation indices for biomass monitoring in barley. International Journal of Applied Earth Observation and Geoinformation, 39, 79–87. Bohman, B. J., Rosen, C. J., & Mulla, D. J. (2019). Evaluation of Variable Rate Nitrogen and Reduced Irrigation Management for Potato Production. Agronomy Journal, 111(4), 2005– 2017. Brocks, S., & Bareth, G. (2018). Estimating barley biomass with crop surface models from oblique RGB imagery. Remote Sensing, 10(2). Bullock, D. G., & Anderson, D. S. (1998). Evaluation of the Minolta SPAD‐502 chlorophyll meter for nitrogen management in corn. Journal of Plant Nutrition, 21(4), 741–755. Buschmann, C., & Nagel, E. (1993). In vivo spectroscopy and internal optics of leaves as basis for remote sensing of vegetation. International Journal of Remote Sensing, 14(4), 711–722. Croft, H., Chen, J. M., & Zhang, Y. (2014). The applicability of empirical vegetation indices for determining leaf chlorophyll content over different leaf and canopy structures. Ecological 85 Complexity, 17(1), 119–130. De Souza, C. H. W., Lamparelli, R. A. C., Rocha, J. V., & Magalhães, P. S. G. (2017). Height estimation of sugarcane using an unmanned aerial system (UAS) based on structure from motion (SfM) point clouds. International Journal of Remote Sensing, 38(8–10), 2218–2230. Dong, J., Wu, F., & Zhang, G. (2005). Effect of cadmium on growth and photosynthesis of tomato seedlings. Journal of Zhejiang University SCIENCE, 6B(10), 974–980. Enciso, J., Avila, C. A., Jung, J., Elsayed-Farag, S., Chang, A., Yeom, J., Landivar, J., Maeda, M., & Chavez, J. C. (2019). Validation of agronomic UAV and field measurements for tomato varieties. Computers and Electronics in Agriculture, 158(February), 278–283. Fontes, P. C. R., & de Araujo, C. (2006). Use of a chlorophyll meter and plant visual aspect for nitrogen management in tomato fertigation. Journal of Applied Horticulture, 8(1), 8–11. Geipel, J., Link, J., & Claupein, W. (2014). Combined spectral and spatial modeling of corn yield based on aerial images and crop surface models acquired with an unmanned aircraft system. Remote Sensing, 6(11), 10335–10355. Gianquinto, G., Sambo, P., & Borsato, D. (2006). Determination of SPAD threshold values for the optimisation of nitrogen supply in processing tomato. Acta Horticulturae, 700, 159–166. Gräler, B., Pebesma, E., & Heuvelink, G. (2016). Spatio-temporal interpolation using gstat. R Journal, 8(1), 204–218. Hiemstra, P. H., Pebesma, E. J., Twenhöfel, C. J. W., & Heuvelink, G. B. M. (2009). Real-time automatic interpolation of ambient gamma dose rates from the Dutch radioactivity monitoring network. Computers and Geosciences, 35(8), 1711–1721. Hijmans, R. J. (2020). raster: Geographic Data Analysis and Modeling. R Package Version 3.4- 5. https://cran.r-project.org/package=raster Holman, F. H., Riche, A. B., Michalski, A., Castle, M., Wooster, M. J., & Hawkesford, M. J. (2016). High throughput field phenotyping of wheat plant height and growth rate in field plot trials using UAV based remote sensing. Remote Sensing, 8(12). Hothorn, T., Bretz, F., & Westfall, P. (2008). Simultaneous inference in general parametric models. Biometrical Journal, 50(3), 346–363. Hunt, E. R., Daughtry, C. S. T., Eitel, J. U. H., & Long, D. S. (2011). Remote sensing leaf chlorophyll content using a visible band index. Agronomy Journal, 103(4), 1090–1099. Hunt, E. R., Donald, J., Spinelli, C. B., Turner, R. W., Bruce, A. E., Gadler, D. J., Brungardt, J. J., & Hamm, P. B. (2018). Monitoring nitrogen status of potatoes using small unmanned aerial vehicles. Precision Agriculture, 19(2), 314–333. 86 Hunt, E. R., Doraiswamy, P. C., McMurtrey, J. E., Daughtry, C. S. T., Perry, E. M., & Akhmedov, B. (2012). A visible band index for remote sensing leaf chlorophyll content at the Canopy scale. International Journal of Applied Earth Observation and Geoinformation, 21(1), 103–112. Hussain, F., Bronson, K. F., Yadvinder-Singh, Bijay-Singh, & Peng, S. (2000). Use of chlorophyll meter sufficiency indices for nitrogen management of irrigated rice in Asia. Agronomy Journal, 92(5), 875–879. Ihuoma, S. O., & Madramootoo, C. A. (2020). ScienceDirect Narrow-band reflectance indices for mapping the combined effects of water and nitrogen stress in field grown tomato crops. Biosystems Engineering, 192, 133–143. Kaplan, G., Fine, L., Lukyanov, V., Manivasagam, V. S., Malachy, N., & Tanny, J. (2021). Estimating Processing Tomato Water Consumption , Leaf Area Index , and Height Using Sentinel-2 and VENµS Imagery. Remote Sensing, 13, 1046. Lenth, R. (2020). emmeans: Estimated Marginal Means, aka Least-Squares Means. R Package Version 1.5.2-1. https://cran.r-project.org/package=emmeans Li, F., Miao, Y., Feng, G., Yuan, F., Yue, S., Gao, X., Liu, Y., Liu, B., Ustin, S. L., & Chen, X. (2014). Field Crops Research Improving estimation of summer maize nitrogen status with red edge-based spectral vegetation indices. Field Crops Research, 157, 111–123. Li, W., Niu, Z., Chen, H., Li, D., Wu, M., & Zhao, W. (2016). Remote estimation of canopy height and aboveground biomass of maize using high-resolution stereo images from a low- cost unmanned aerial vehicle system. Ecological Indicators, 67, 637–648. Ling, Q., Huang, W., & Jarvis, P. (2011). Use of a SPAD-502 meter to measure leaf chlorophyll concentration in Arabidopsis thaliana. Photosynthesis Research, 107(2), 209–214. Louhaichi, M., Borman, M. M., & Johnson, D. E. (2001). Spatially located platform and aerial photography for documentation of grazing impacts on wheat. Geocarto International, 16(1), 65–70. Main, R., Cho, M. A., Mathieu, R., O’Kennedy, M. M., Ramoelo, A., & Koch, S. (2011). An investigation into robust spectral indices for leaf chlorophyll estimation. ISPRS Journal of Photogrammetry and Remote Sensing, 66(6), 751–761. Meyer, G. E., & Neto, J. C. (2008). Verification of color vegetation indices for automated crop imaging applications. Computers and Electronics in Agriculture, 63(2), 282–293. Minotti, P. L., Halseth, D. E., & Sieczka, J. B. (1994). Field Chlorophyll Measurements to Assess the Nitrogen Status of Potato Varieties. HortScience, 29(12), 1497–1500. Mulla, D. J. (2013). Twenty five years of remote sensing in precision agriculture: Key advances 87 and remaining knowledge gaps. Biosystems Engineering, 114(4), 358–371. Muñoz-Carpena, R., Icerman, J., Dukes, M. D., Zotarelli, L., & Scholberg, J. M. (2008). Tomato yield, biomass accumulation, root distribution and irrigation water use efficiency on a sandy soil, as affected by nitrogen rate and irrigation scheduling. Agricultural Water Management, 96(1), 23–34. Muñoz-Huerta, R. F., Guevara-Gonzalez, R. G., Contreras-Medina, L. M., Torres-Pacheco, I., Prado-Olivarez, J., & Ocampo-Velazquez, R. V. (2013). A review of methods for sensing the nitrogen status in plants: Advantages, disadvantages and recent advances. Sensors (Switzerland), 13(8), 10823–10843. Penuelas, J., Baret, F., & Filella, I. (1995). Semi-empirical indices to assess carotenoids/chlorophyll a ratio from leaf spectral reflectance. Photosynthetica, 31(2), 221– 230. Peñuelas, J., Gamon, J. A., Fredeen, A. L., Merino, J., & Field, C. B. (1994). Reflectance Indices Associated With Physiological Changes in Nitrogen- and Water-Limited Sunflower Leaves (pp. 135–146). Pinter, P. J., Hatfield, J. L., Schepers, J. S., Barnes, E. M., Moran, M. S., Daughtry, C. S. T., & Upchurch, D. R. (2003). Remote Sensing for Crop Management. Photogrammetric Engineering and Remote Sensing, 69(6), 647–664. R Core Team. (2020). R: A Language and Environment for Statistical Computing (3.5.1). R Foundation for Statistical Computing. https://www.r-project.org/ Rouse, J. W., Haas, R. H., Schell, J. A., & Deering, D. W. (1973). Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation. Sandoval-Villa, M., Guertal, E. A., & Wood, C. W. (2000). Tomato leaf chlorophyll meter readings as affected by variety, nitrogen form, and nighttime nutrient solution strength. Journal of Plant Nutrition, 23(5), 649–661. Ulissi, V., Antonucci, F., Benincasa, P., Farneselli, M., Tosti, G., Guiducci, M., Tei, F., Costa, C., Pallottino, F., Pari, L., & Menesatti, P. (2011). Nitrogen concentration estimation in tomato leaves by VIS-NIR non-destructive spectroscopy. Sensors, 11(6), 6411–6424. University of Florida IFAS Extension. (2015). Nutrient Management of Vegetable and Row Crops Handbook. Varvel, G. E., Schepers, J. S., & Francis, D. D. (1997). Ability for In-Season Correction of Nitrogen Deficiency in Corn Using Chlorophyll Meters. Soil Science Society of America Journal, 61(4), 1233–1239. Westerveld, S. M., McKeown, A. W., McDonald, M. R., & Scott-Dupree, C. D. (2003). 88 Chlorophyll and nitrate meters as nitrogen monitoring tools for selected vegetables in southern Ontario. Acta Horticulturae, 627, 259–266. Wu, J., Wang, D., Rosen, C. J., & Bauer, M. E. (2007). Comparison of petiole nitrate concentrations , SPAD chlorophyll readings , and QuickBird satellite imagery in detecting nitrogen status of potato canopies. Field Crops Research, 101, 96–103. Wu, Q., Sun, H., Li, M. Z., & Yang, W. (2014). Development and application of crop monitoring system for detecting chlorophyll content of tomato seedlings. International Journal of Agricultural and Biological Engineering, 7(2), 138–145. Yang, H., Li, J., Yang, J., Wang, H., Zou, J., & He, J. (2014). Effects of Nitrogen Application Rate and Leaf Age on the Distribution Pattern of Leaf SPAD Readings in the Rice Canopy. PLoS ONE, 9(2), 1–11. Yue, J., Yang, G., Li, C., Li, Z., Wang, Y., Feng, H., & Xu, B. (2017). Estimation of winter wheat above-ground biomass using unmanned aerial vehicle-based snapshot hyperspectral sensor and crop height improved models. Remote Sensing, 9(7). Zarco-Tejada, P. J., Berjón, A., López-Lozano, R., Miller, J. R., Martín, P., Cachorro, V., González, M. R., & De Frutos, A. (2005). Assessing vineyard condition with hyperspectral indices: Leaf and canopy reflectance simulation in a row-structured discontinuous canopy. Remote Sensing of Environment, 99(3), 271–287. 89