Cotton Yield Estimate Using Sentinel-2 Data and an Ecosystem Model over the Southern US

Abstract

High-resolution data with nearly global coverage from Sentinel-2 mission open a new era for crop growth monitoring and yield estimation from remote sensing. The objective of this study is to demonstrate the potential of using Sentinel-2 biophysical data combined with an ecosystem modeling approach for estimation of cotton yield in the southern United States (US). The Boreal Ecosystems Productivity Simulator (BEPS) ecosystem model was used to simulate the cotton gross primary production (GPP) over three Sentinel-2 tiles located in Mississippi, Georgia, and Texas in 2017. Leaf area index (LAI) derived from Sentinel-2 measurements and hourly meteorological data from Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) reanalysis were used to drive the ecosystem model. The simulated GPP values at 20-m grid spacing were aggregated to the county level (17 counties in total) and compared to the cotton lint yield estimates at the county level which are available from National Agricultural Statistics Service in the United States Department of Agriculture. The results of the comparison show that the BEPS-simulated cotton GPP explains 85% of variation in cotton yield. Our study suggests that the integration of Sentinel-2 LAI time series into the ecosystem model results in reliable estimates of cotton yield.

1. Introduction

There are two major cotton types in the United States (US), one is the Pima cotton, located in the counties of Pinal and Graham in Arizona, and another type is the upland cotton (Gossypium hirsutum), a species that is native to Central America, Mexico, and the Caribbean and Gulf Coast regions, and is mainly located in southern states below 36°30′ north (N), since cotton grows well at higher temperature and in sunlight. Upland cotton is the most common type of cotton in the US, making up 95% of the planted cotton area in the US, and it is the focus of this study. Although cotton is naturally a perennial shrub, it is primarily grown as an annual crop to help pest control. In the US, the cotton is planted in spring, varying from the beginning of February to the beginning of June, and harvested in the fall. Its growing season is approximately 150 to 180 days, making it the longest of annually planted crops in the US. Since cotton is a thermophilic crop, the growing season length is important to the cotton yield. Although cotton is somewhat salt- and drought-tolerant and is an attractive crop for arid and semiarid regions, the cotton yield heavily relies on irrigation.

The cotton yield is determined by a few key factors, including species and genetics [1], soil [2] and nutrition [3,4] conditions, weather [5], irrigation [6,7], and management events [8]. Specifically, cotton is a perennial warm-season crop that requires a long growing period. Therefore, the length of growing season for cotton has a significant effect on cotton yield. For example, lint yield in double-cropping systems is significantly lower than that in monocultures [8,9].

Reliable estimation of cotton yield relies on the available information during the cotton growing season. Crop growth models, such as the Decision Support System for Agrotechnology Transfer (DSSAT) [10], are very mature in yield simulations at the plot level; however, they heavily rely on accurate inputs, such as numerous cotton traits, which are usually unavailable for large-scale applications. In contrast, remote sensing was proven as a practical approach for yield estimation at the regional and global scales [11,12,13,14,15,16,17,18]. However, the accuracy of cotton yield estimation based on remote sensing still relies on how much information is used to capture the important physiological stages during cotton development [19,20,21,22].

Using a single vegetation index [23], such as normalized difference vegetation index (NDVI), from a specific date (e.g., the maximum NDVI) can only provide limited information about cotton yield. For example, directly relating NDVI to cotton yield often gives low correlation (R² = 0.47) [24]. Zhao et al. [25] suggested that lint yield best correlates to the canopy reflectance indices measured at the early flower stage of cotton growth (R² of 0.56–0.89). Yang et al. [26] found that cotton yield variation (R² = 0.52) is mostly explained by airborne hyperspectral imagery collected 9.9 weeks after planting. Inclusion of multitemporal remote sensing to capture the phenology of cotton improves the accuracy of cotton yield estimation [27,28,29,30].

Further improvement of the estimation of cotton yield requires the integration of a crop growth simulation model, remote sensing, and climate data via a data assimilation approach [31,32,33]. Essentially, remote sensing provides the biophysical parameters, e.g., leaf area index (LAI), to the crop growth model for the simulation of gross primary production (GPP), net primary production (NPP), and cotton biomass, where GPP and NPP are the indicators of crop productivity [34]. GPP is the rate at which plants capture and store atmospheric carbon dioxide to generate oxygen and energy as biomass. NPP is the difference between GPP and the energy lost during plant autotrophic respiration. NPP, thus, represents the daily accumulation of standing biomass. Since yield is related to the cumulative production of biomass during a growing season, integration of multi-source datasets achieves high accuracy in cotton yield estimation, e.g., R² = 89% in the study of Hebbar et al. [35]. By incorporating geospatial data from remote sensing, such as LAI, into a crop growth model, Thorp et al. [36] reported that the root-mean-squared errors between measured and simulated seed cotton yield were 5% or less. In contrast to the use of complicated crop growth models, the light use efficiency (LUE) model is simple for implementation in cotton yield estimation; the LUE model can incorporate some meteorological parameters, such as temperature and vapor pressure deficit, into the yield estimation process [37], while it often lacks the capability to simulate soil water stress and the mechanism to explain the effects of CO₂ fertilization (i.e., the increased rate of photosynthesis in plants that results from increased levels of carbon dioxide in the atmosphere) [38].

In addition to the improvements in modeling, the satellite images for cotton yield estimation evolved from low to high spatial resolution, e.g., Advanced Very-High-Resolution Radiometer [39], Moderate Resolution Imaging Spectroradiometer (MODIS) [40], Landsat [41], and Sentinel-2 [42]. In particular, the Sentinel-2 mission is a constellation that consists of twin satellites (Sentinel-2A and Sentinel-2B) and is operated by the European Union (EU) Copernicus Program. With these two satellites, Sentinel-2 is able to map land surface with 13 bands in the visible, near-infrared, and shortwave-infrared spectrum, a five-day revisiting time, and spatial resolutions of up to 10 m. Given its powerful capability of data collection and free, full, and open data policy, Sentinel-2 opens a new era for operational agricultural monitoring.

With only Sentinel-2A data and a 10-day revisit interval, cloud contamination is an open issue for consistently tracking the growing conditions of cotton; this challenge is similar to using only Landsat 8 for crop growth monitoring. To mitigate the cloud contamination, the available approach is to blend MODIS and Landsat for practical mapping of crop progress [43,44]. Since 2017, both Sentinel-2A and B data are available; this gives an opportunity to rely on a single data source for operational agricultural applications, as the probability of cloud contamination is reduced.

So far, the utilities of Sentinel-2A and 2B data for cotton yield mapping are still limited [42,45,46]. The objective of this study was to evaluate the potential of Sentinel-2 data as the key input to simulate cotton GPP for the purpose of yield estimation in an ecosystem model. For this aim, we collected Sentinel-2 reflectance data to derive the time series of LAI as the input of an ecosystem model; the ecosystem model integrated LAI, meteorological data, and soil data to simulate cotton GPP; the cotton GPP–yield relationship was then calibrated using yield data collected from the National Agricultural Statistics Service (NASS) in the United States Department of Agriculture (Section 2 and Section 3). The uncertainty of the GPP simulation and conclusions are summarized in Section 4 and Section 5.

2. Materials and Methods

2.1. Ecosystem Model for GPP Simulation

The Boreal Ecosystems Productivity Simulator (BEPS), hourly version, is a process-based ecosystem model to simulate water, energy, and carbon budgets [47,48,49,50,51,52]; it was used in this study to simulate cotton GPP. BEPS is in development for the last 20 years [53]. Although BEPS was initially developed for boreal ecosystems, it was expanded and used for temperate and tropical ecosystems in regional and global scales [48,51,52,54,55,56,57,58,59] and for the estimation of winter wheat yield in China [18], as C3 plants share the same photosynthesis theory at the leaf level, i.e., Farquhar’s leaf-level biochemical model [60], and BEPS has a “two-leaf” approach to upscale the ecosystem GPP to canopy level for various ecosystems including crops [49,61].

In BEPS, the canopy-level photosynthesis (GPP; $A_{canopy}$) is simulated as the sum of the total photosynthesis of sunlit and shaded leaf groups [49]:

$$A_{canopy}=A_{sun}(g_{sc\_sun})L_{sun}+A_{sh}(g_{sc\_sh})L_{sh}$$

(1)

where the subscripts “sun” and “sh” denote the sunlit and shaded components of the photosynthesis (A) and LAI (L), whereas $g_{sc}$ is the stomatal resistance for carbon molecules. The sunlit and shaded LAI are separated by Equation (2) [49,61].

$$\begin{array}{l}{L}_{sun}=2\cos \theta \left(1-e^{-0.5\Omega L/\cos \theta }\right)\\ L_{sh}=L-L_{sun}\end{array}$$

(2)

where $\theta$ is the solar zenith angle, and $\Omega$ is the clumping index. In this study, the Sentinel-2 data were integrated into the GPP simulation through the LAI product as described in Section 2.4.

In BEPS, the vegetation biophysical and biochemical parameters are defined for plant function types including crops. For crop GPP simulation (cotton in this study), the clumping index was set to as 0.8. A detailed description of BEPS is provided in the Supplementary Materials.

2.2. Cotton-Specific Photosynthetic Parameters

Cotton leaf is characterized by a very high light saturation point (LSP). In well-watered conditions, the cotton LSP can reach 2304 µmol m⁻²∙s⁻¹, and the maximum CO₂ assimilation rate can reach 43.1 µmol∙m⁻²∙s⁻¹ according to the experiments by Zhang et al. [62]. In photosynthesis model, V_cmax and J_max are two important parameters to describe the photosynthetic rate at the leaf level. V_cmax is the maximum rate of carboxylation for Rubisco-limited (or light-saturated) conditions, and J_max is the maximum rate of electron transport for the light-limited condition [60]. According to the measurements of by Harley et al. [63], cotton has a temperature optimum near 40 °C for V_cmax and 35 °C for J_max. The V_cmax and J_max at the optimum temperatures are four and two times higher than their respective values at 25 °C [64]; therefore, the optimum temperatures for maximum carboxylation and electron transport were set to 313.15 K and 308.15K, respectively, in BEPS. The V_cmax at 25 °C varies with leaf nitrogen and phosphorus [65,66]. Since such nutrition information is not reliably available from Sentinel-2, the cotton V_cmax at 25 °C was set to 120 µmol∙m⁻²∙s⁻¹ to be consistent with our previous modeling for crops [52]. For BEPS simulation, the CO₂ concentration data were obtained from https://www.esrl.noaa.gov/gmd/ccgg/trends/global.html.

2.3. Yield Estimation and Prediction

Given the total cotton GPP (g∙C∙m⁻²∙year⁻¹) in the growing season simulated by BEPS, the cotton lint yield (Y, in units of kg∙ha⁻¹) was estimated as follows:

Y = GPP × CUE × b × HI × l × 10,

(3)

where CUE is the carbon use efficiency or the NPP-to-GPP ratio, ranging from 0 to 1 for different species [67] and often set to 0.5 as an average of global ecosystems [68]; b is the biomass to carbon ratio; since Schlesinger [69] noted that the carbon content of biomass is almost always found to be between 45 and 50% (by oven-dry mass), b is close to 2:1; HI is the harvest index as the ratio of harvested product (grain for food crops, lint and seed for cotton) to total dry biomass (sometimes refer to the total dry aboveground biomass); this is largely determined by the percentage of bolls in the total dry biomass of cotton; l is the lint percentage of cotton (lint weight to seed cotton weight), and 10 is the ratio from g∙m⁻² to kg∙ha⁻¹.

The cotton HI varies with species and is affected by human management. The HI (seed cotton yield to total aboveground biomass) shown in Hussein et al. [70] was 0.34, and the lint percentage was 41%. Maheswarappa et al. [71] reported an HI of 0.3; however, the cotton HI was less than 0.1 according to the report by Pettigrew et al. [72], with a lint percentage of 41.4% on average under the rotation conditions. Huang [73] reported HI values (seed cotton to aboveground biomass) from 0.28 to 0.33. According to Dowd et al. [74], the seed-to-fiber ratio was 1.41 ± 0.11, corresponding to a lint percentage of 41% (1/2.41). Pettigrew, Bruns, and Reddy [72] reported similar lint percentages (41%). If the HI includes both seed cotton and boll trash, then l is the gin turnout ratio, which is less than the lint percentage.

Given the small variations in CUE, b, HI, and l for cotton, it is desirable to directly link GPP to yield by one coefficient C for practical application as follows:

$$Y=GPP\cdot C$$

(4)

Equation (4) assumes a linear relationship between GPP and yield; the correlation coefficient is examined below. Such a linear relationship is backed up by the strong correlation (R² = 0.93) between biomass and seed cotton yield as measured by Huang [73]; in contrast, the LAI can only explain 51% of variance in seed cotton yield as shown in the same study. Similarly, Lambert, Traoré, Blaes, Baret, and Defourny [42] reported that the aboveground biomass explains more variation in crop yield than the peak LAI.

Once the coefficient C is derived, the yield for each pixel can be estimated with cotton GPP at the end of the growing season, and the cotton production (P) for a given region can be estimated by

$$P={\displaystyle \sum _{i=1}^n{Y}_i\cdot A}$$

(5)

where Y_i is the cotton yield estimation for each pixel, A is the pixel area, and n is the number of pixels.

2.4. Sources of Model Input

The Copernicus Sentinel-2 mission comprises twin polar-orbiting satellites in the same sun-synchronous orbit, phased at 180°, aiming at monitoring land surface conditions at five-day intervals under cloud-free conditions. Each satellite carries a multispectral imager (MSI) with a swath of 290 km. The MSI provides a set of 13 spectral bands spanning from the visible and near-infrared to the shortwave-infrared, featuring four spectral bands at 10-m, six bands at 20-m, and three bands at 60-m spatial resolution. Sentinel-2A was launched on 23 June 2015 and Sentinel-2B was launched on 7 March 2017.

In this study, Sentinel-2 Level-1C data, top-of-atmosphere (TOA) reflectances, were used as input to the Sen2Cor atmospheric correction processor in order to convert the TOA reflectance images to the top-of-canopy (TOC) reflectance [75]. The Level-1C product results from using a digital elevation model (DEM) to project the images in cartographic geometry and is organized into tiles having an area of about 100 km² in the Universal Transverse Mercator projection and World Geodetic System 1984 ellipsoid; the Level-1C product is available from the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home). Three tiles (14SKB, 15SYT, and 17RKQ) in 2017 were chosen in order to cover a range of high to low cotton yields in southern US as summarized in

Table 1. Summary of Sentinel-2 images used in this study. The three tiles cover 17 United States (US) counties in total. ID—identifier; TX—Texas; MS—Mississippi; GA—Georgia.

Tile ID	Number of Images	Summary
14SKB	25	TX, mainly rainfed
15SYT	16	MS, largely irrigated
17RKQ	20	GA, mainly rainfed

. We chose 2017 as this is the first year when both Sentinel-2A and Sentinel-2B data were available in order to minimize the contamination of cloud for cotton yield mapping; furthermore, cotton yield data for validation in 2018 were unavailable at the time when the study was conducted. The tile in Georgia (GA; 17RKQ) represents the cotton along the coastal states and a cotton yield in the middle ranges, the tile in Mississippi (MS; 15SYT) represents the cotton with better irrigation conditions and a higher yield, and the tile in Texas (TX; 14SKB) represents drier conditions and a lower cotton yield. The locations of these three tiles are shown in Figure 1.

LAI is defined as one-half of the total green leaf area per unit ground surface area (m²/m²) for broadleaf canopies [76]. It is an important structural property of vegetation. Because leaf surfaces are the common boundary for atmosphere-biosphere energy and mass exchange, important processes such as canopy interception, evapotranspiration, and gross photosynthesis are directly proportional to LAI. In BEPS, LAI is one of the most important parameters to determine the total photosynthetic rate at the canopy level. The LAI data are derived using TOC reflectance as the input to the biophysical processor available from the Sentinel Application Platform (SNAP). In SNAP, the LAI retrieval is based on an artificial neural network that relates TOC reflectances in eight bands and sun-view geometries to LAI [77]. The eight bands are B3 (560 nm), B4 (665 nm), B5 (705 nm), B6 (740 nm), B7 (783 nm), B8a (865 nm), B11 (1610 nm), and B12 (2190 nm). For the neural network training, TOC reflectances were simulated with the PROSPECT + SAIL (Scattering by Arbitrarily Inclined Leaves) model [78] corresponding to the eight Sentinel-2 bands and considering variations of primary biophysical variables related to leaf optical properties and to the canopy structure. According to Verrelst, Muñoz, Alonso, Delegido, Rivera, Camps-Valls, and Moreno [77], the uncertainty of LAI retrievals from SNAP is 0.51 m²/m². The LAI products are derived in 20-m resolution.

Although Sentinel-2 has a high revisit frequency, the availability of cloud-free and cloud shadow-free images is limited depending on the tile location. The LAI retrievals from SNAP are often contaminated by cloud and cloud shadow; therefore, abrupt decreased LAI values are shown in the time series. To minimize the effects of cloud contamination, a locally adjusted cubic-spline capping (LACC) method was used to screen affected LAI data points and to replace them through temporal interpolation. [79]. The interpolated LAI data at the daily step were used as input to BEPS for GPP simulation.

A global reanalysis dataset, Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), was chosen to drive BEPS simulation [80]. Six key meteorological elements including relative humidity, wind speed, and air temperature at 2 m above the surface, surface atmosphere pressure and incoming solar shortwave flux, and total precipitation at the surface level were resampled into 20-m resolution as input.

In BEPS, the soil is stratified into five layers (up to 2 m depth) for the simulation of soil water and heat fluxes. For BEPS input, four soil attributes, including texture, bulk density, organic content, and reference depth, were extracted from Soil Grids (https://soilgrids.org) [81], which provides global soil attributes at 250-m resolution.

After the cotton is harvested from the field, the raw cotton (or seed cotton) goes through the ginning processing, in order to separate the raw fiber (cotton lint) from the seeds and other debris. The raw fiber is compressed into bales. According to the Joint Cotton Industry Bale Packaging Committee (JCIBPC) in the US established in 2001, a standard cotton bale weighs 500 pounds. However, specific cotton bale weight varies between 480 and 500 pounds in different states; for example, the average bale weight in 2017 was 490 pounds (https://www.nass.usda.gov/Publications/Todays_Reports/reports/ctgnan18.pdf). For statistical purposes, NASS often reports cotton yield based on 480-pound net weight bales. In this study, we used the reported planted area (acres) and quantity of bales to calculate yield, and the bale weight was taken as 490 pounds according to the NASS statistics. The cotton yield data at the county level in 2017 from NASS were used in this study. The NASS cotton yield data are based on surveys and are the only publicly available data that cover statistics of cotton yield for all US counties.

The USDA Cropland Data Layer (CDL) product in 2017 [82,83] was used to mask out non-cotton areas for GPP simulation in the same year. The simulated cotton GPP values at 20-m resolution for each US county were averaged according to their Federal Information Processing Standard (FIPS) county code; the county-average GPP for cotton was then compared to the county-level yield data from NASS.

3. Results

3.1. Spatial Distribution of Cotton GPP in Three Sentinel-2 Tiles

One challenge to accurately simulate cotton GPP using BEPS is to determine the irrigation information. We developed an irrigation module in the GPP simulation; when this module is turned on, the soil moisture in each layer is increased to 0.42 m³/m³ (volumetric soil water content) if it is below this level every few days (e.g., 14 days). We chose 0.42 m³/m³ because this is the maximum field capacity for all soil textures [84]. In this study, there is no ideal way to identify which and when a pixel is irrigated. According to the NASS statistics, the ratios of irrigated cotton area are 32%, 34%, and 51% for Texas (TX), Georgia (GA), and Mississippi (MS), respectively; therefore, the irrigation module was turned on for tile 15SYT (MS) only.

Another challenge is to determine the start and end of growing season for cotton. In MS, the cotton is usually planted in later April or May, and the harvest begins in October and November. In the southern part of the Great Plains (covering 14SKB), the plant date can be in June. Since both cotton plots in tiles 14SKB and 15SYT are single-cropping systems, the LAI out of the growing season is low; therefore, the annual GPP is a good proxy for cotton GPP in the whole growing season. This is not true for the tile 17RKQ in GA. In the coastal plain, especially in lower south Georgia where the growing season is long, double cropping is popular. GA is characterized by a wheat–cotton cropping system, and the cotton is often planted in June; therefore, the cotton growing season for BEPS simulation was limited from June 1 to October 31.

Mainly rainfed, the tile 14SKB in TX is characterized by low GPP values (726 ± 74 g∙C∙m⁻²∙year⁻¹) from the BEPS simulation, although some high GPP values (~1500 g∙C∙m⁻²∙year⁻¹) appear to be associated with the center-pivot irrigation (Figure 2). The GPP in the cotton growing season for tile 17RKQ is 1482 ± 74 g∙C∙m⁻²∙year⁻¹ (Figure 3). Owing to the long growing season, the annual GPP for tile 15SYT is 2005 ± 177 g∙C∙m⁻²∙year⁻¹ (Figure 4). The statistics of average GPP values for each county are listed in

Table 2. The statistics of cotton area, production, and yield from the National Agricultural Statistics Service (NASS) and the Boreal Ecosystems Productivity Simulator (BEPS) in 2017 for the three Sentinel-2 tiles. FIPS—Federal Information Processing Standard; GPP—gross primary production.

County	FIPS	Area	Production	Lint Yield	GPP (BEPS)
		(ha)	(kg)	(kg/ha)	(g∙C∙m−2∙year−1)
TX (14SKB)
Lynn	48,305	120,522	79,922,356	663	715
Garza	48,169	16,593	12,322,983	743	752
Dawson	48,115	105,198	69,978,487	665	805
Borden	48,033	23,329	12,449,007	534	630
MS (15SYT)
Coahoma	28,027	36,881	42,187,930	1144	1722
Quitman	28,119	10,192	11,590,623	1137	1985
Tallahatchie	28,135	17,839	23,602,659	1323	2116
Sunflower	28,133	8379	10,650,446	1271	1891
Leflore	28,083	17,748	21,326,898	1202	2200
Carroll	28,015	8859	10,161,688	1147	2118
GA (17RKQ)
Worth	13,321	22,541	21,292,891	945	1532
Tift	13,277	9227	8,268,221	896	1409
Colquitt	13,071	19,223	20,181,571	1050	1536
Cook	13,075	7244	5,267,657	727	1364
Berrien	13,019	11,088	9,001,692	812	1451
Thomas	13,275	12,505	10,779,804	862	1552
Brooks	13,027	15,459	14,291,575	924	1529

.

3.2. The GPP–Lint Yield Relationship for Upland Cotton

As shown in Figure 5, BEPS-simulated GPP strongly correlated with cotton lint yield from the NASS statistics for the 17 counties (R² = 0.85). This suggests that cotton GPP from simulation can be used as a reliable proxy for cotton yield. It is noteworthy to mention that Equations (3) and (4) have no intercepts, while the regression equation in Figure 5 has an intercept. If the intercept is forced to zero in Figure 5, the slope is 0.61 and the R² is 0.67. A similar intercept was also observed when cotton biomass was used to explain the variation in seed cotton as shown in Huang [73]. This small intercept suggests that there may be systematic uncertainty in the GPP simulation or that the correlation between GPP and cotton yield is nonlinear when both GPP and cotton yield are small. For example, the cotton V_cmax was fixed to 120 µmol∙m⁻²∙s⁻¹ in our simulations as, currently, there is no such information available for each pixel from the satellite.

As shown in Table 2, the cotton lint yield in tile 17RKA (GA) was lower than the full-season cotton yield in MS. Cotton needs enough heat to grow. With respect to cotton grown in a double-crop system following wheat, the length of fruiting period and adequate temperature accumulation for boll maturation are reduced, thereby lowering the lint yield potential. Dodds et al. [85] found that, for cotton after wheat in a double-cropping system, the lint yield is generally low depending on the plant dates (lint yield range: 750–944 kg/ha). As shown in the 2012 Georgia cotton production guide (http://www.ugacotton.com/vault/productionguide/2012CottonProductionGuide.pdf), the cotton yield in a double-cropping system can have a yield reduction of up to 30% due to the limited growing season, when compared to full-season cotton planted in early May. In Arkansas, the yield potential decreases by 2% every day if cotton is planted past May 20 (https://www.uaex.edu/publications/PDF/FSA-2163.pdf). This suggests the importance of capturing planting date for cotton yield estimation, while our GPP simulation explained such variation in cotton growing season.

4. Discussions

4.1. Crop Growth Modeling

One of the core tasks for yield mapping is the modeling of crop photosynthesis. The crop growth models share in common the use of either a light use efficiency model or an enzyme kinetics-based model. For example, the modeling of cotton photosynthesis in Agricultural Production Systems Simulator (APSIM) [86] is to relate photosynthesis to radiation using an empirical equation [87]. In Crop Environment REsource Synthesis (CERES)-Wheat, the daily photosynthesis is modeled using an LUE approach [88]; the WOFOST (World Food Studies) [89] also uses Farquhar’s model [60]; the Simple Algorithm for Yield Estimation model (SAFY) also uses an LUE model for photosynthesis [90]. In common, BEPS is a process-based two-leaf enzyme kinetic ecosystem model, and the photosynthesis module of BEPS is sufficient for accurate simulation of crop photosynthesis rate. Although BEPS was initially developed for boreal ecosystems, it was expanded and used for estimations of winter wheat yield in China [18,91], for yield estimation of rapeseed in Hubei province [92], for estimation of maize yield in China [93,94], and for crop GPP modeling in Europe [95] and the US [96].

The abovementioned crop growth models have the strengths of crop trait modeling and handling of soil nutrition and water management. However, such an approach necessitates availability of these crop traits at large scales. Operational cotton yield mapping may require a compromise between large scale and high spatial resolution, and simplicity of the model structure. Similar to this work, two recent studies used Sentinel-2-derived LAI for biomass and yield mapping. Punalekar et al. [97] reported that the R² between pasture biomass and LAI ranged from 0.16–0.73 to 0.22–0.76 for the NDVI-derived LAI and radiative transfer model-based LAI, respectively. Lambert, Traoré, Blaes, Baret, and Defourny [42] demonstrated that the R² values between LAI and yield were 0.68, 0.62, 0.8, and 0.48 for cotton, maize, millet, and sorghum, respectively. Here, we show that the R² value between GPP and lint yield was 0.85. Despite a single factor C in Equation (4) being used, our study demonstrates the advantage of integrating both biotic and abiotic factors for cotton yield mapping at very high spatial resolutions.

4.2. Yield Variations in Southern US Tiles

As expected, the modeled GPP did not show much variations for these tiles (14SKB and 15SYT), showing no obvious variations in cotton yield among counties either. Within tile 17RKQ, the modeled GPP explained 50% of variations in cotton yield within a narrow range. For tile 14SKB with variations in irrigation conditions (Figure 2), there were large variations in modeled GPP ranging from 500 g∙C∙m⁻²∙year⁻¹ to 2500 g∙C∙m⁻²∙year⁻¹, although these variations were averaged at the county level; the area with center-pivot irrigation (circles) was characterized by high GPP values up to 2500 g∙C∙m⁻²∙year⁻¹. Although our study showed the important climatological effect on cotton yield, the region with highest yield was not located in GA where the condition is most ideal because of the short cotton growing season due to crop rotations; instead, the highest yield in our study was shown in MS due to their longer growing season and irrigation. This suggests that our model captured the role of crop management, in addition to the climatological effect.

5. Conclusions

In this study, we evaluated the potential of Sentinel-2A and 2B data for mapping cotton lint yield in an ecosystem model. Three Sentinel-2 tiles in southern US were chosen, covering regions with low to high cotton yields. BEPS-simulated GPP values at the pixel level were averaged into the county level and compared to the cotton lint yield data collected from NASS. Our study suggests that BEPS-simulated GPP explains 85% of variation in cotton lint yield and can be used as a reliable indicator for estimation of cotton lint yield. We also found that there may be systemic bias in our GPP simulation due to limited knowledge about cotton V_cmax, which is related to leaf nitrogen content. Accurate GPP simulation for cotton in a double-cropping system relies on the knowledge of the start and end dates of the cotton growing season, and this posed a challenge in our simulation due to certain cloud contamination in Sentinel-2 images in southern Georgia; therefore, we suggest the fusion of multi-satellite data, such as from Landsat 7 and 8, Sentinel-2, or synthetic aperture radar data from Sentinel-1 at similar spatial resolutions in order to track the key phenology stages of cotton development for a more reliable estimation of cotton yield.