Upscaling Gross Primary Production in Corn-Soybean Rotation Systems in the Midwest

The Midwestern US is dominated by corn (Zea mays L.) and soybean (Glycine max [L.] Merr.) production, and the carbon dynamics of this region are dominated by these production systems. An accurate regional estimate of gross primary production (GPP) is imperative and requires upscaling approaches. The aim of this study was to upscale corn and soybean GPP (referred to as GPPcalc) in four counties in Central Iowa in the 2016 growing season (DOY 145–269). Eight eddy-covariance (EC) stations recorded carbon dioxide fluxes of corn (n = 4) and soybean (n = 4), and net ecosystem production (NEP) was partitioned into GPP and ecosystem respiration (RE). Additional field-measured NDVI was used to calculate radiation use efficiency (RUEmax). GPPcalc was calculated using 16 MODIS satellite images, ground-based RUEmax and meteorological data, and improved land use maps. Seasonal NEP, GPP, and RE ( x ¯ ± SE) were 678 ± 63, 1483 ± 100, and −805 ± 40 g C m−2 for corn, and 263 ± 40, 811 ± 53, and −548 ± 14 g C m−2 for soybean, respectively. Field-measured NDVI aligned well with MODIS fPAR (R2 = 0.99), and the calculated RUEmax was 3.24 and 1.90 g C MJ−1 for corn and soybean, respectively. The GPPcalc vs. EC-derived GPP had a RMSE of 2.24 and 2.81 g C m−2 d−1, for corn and soybean, respectively, which is an improvement to the GPPMODIS product (2.44 and 3.30 g C m−2 d−1, respectively). Corn yield, calculated from GPPcalc (12.82 ± 0.65 Mg ha−1), corresponded well to official yield data (13.09 ± 0.09 Mg ha−1), while soybean yield was overestimated (6.73 ± 0.27 vs. 4.03 ± 0.04 Mg ha−1). The approach presented has the potential to increase the accuracy of regional corn and soybean GPP and grain yield estimates by integrating field-based flux estimates with remote sensing reflectance observations and high-resolution land use maps.


Introduction
Midwestern US agriculture is dominated by corn (Zea mays L.) and soybean (Glycine max [L.] Merr.) production and both crops are important export commodities. In 2017, corn and soybean production in the United States made up 33% and 34% of total world production, and 56% and 46% of US production came from the Corn Belt states (Iowa, Nebraska, Indiana, Illinois, and Kansas), respectively [1,2]. Vast agricultural areas are used in the production of these crops which can impact gross primary production (GPP) and the global carbon cycle. The Midwestern US region was found to have one of the highest GPP rates in the northern hemisphere, owing to the large-scale production of corn, a C4 crop [3]. In addition, conventional corn and soybean production with bare soil in the off-season and frequent tillage events can impact GPP compared to reduced till systems [4]. Despite recent improvements in

Study Sites
The study was conducted as part of the Soil Moisture Active Passive Validation Experiment in 2016 (SMAPVEX16) during the growing season (i.e., day of year (DOY) 145-269; total: 125 days) in four neighboring counties in North Central Iowa: Story, Franklin, Hamilton, and Hardin County [19,20]. The overall experiment was dedicated to exploring the land surface soil moisture conditions across a homogeneous landscape to improve vegetation characterization and soil moisture remote sensing algorithms via the Soil Moisture Active Passive Mission. Within each county, two field sites under a corn-soybean rotation, one each planted in corn and the other in soybeans (Figure 1). The sites are located on the Des Moines Lobe and soils belong to the Clarion-Nicollet-Webster association [21]. The site characteristics are summarized in Table 1. The sites reflect well the dominating crop and soil management conditions of the four counties, with 98% of the harvested cropland used for corn and soybean production, and only 13% of the cropped area under no-till and 3% planted with cover crops [5]. In terms of climatic parameters, the 30-year (1986-2016) minimum and maximum temperature between June and September averaged (± SE) over four nearby weather stations was 15.38 ± 0.16 and 26.94 ± 0.08 • C. All agricultural sites are rain-fed systems with a seasonal 30-year precipitation of 508 ± 17 mm [22].   [18]: Story County (lower), Hamilton County (center left), Hardin County (center right), and Franklin County (upper). Corn fields are marked dark green, and soybean fields are light green (each land use area > 55%). The red triangles and yellow circles show the location of the eddy-covariance (EC) stations in soybean and corn fields, respectively.

Eddy Covariance Datasets
Turbulent fluxes of carbon dioxide (F c ) were computed at each site (n = 8 EC stations) using the EC method [23]. One EC station was deployed in each field, so that the EC footprint covered the crop of interest in a representative portion of the field and in the main wind direction. The 90% footprint of each EC station was calculated [24], and mean (± SD) flux footprint area was 4.4 ± 2.3 ha ( Figure A1). The 90% footprint area at the corn-site in Franklin County surpassed the field, and the footprint was cropped to field boundaries. Note that most of the information recorded by the sensor is from nearby the EC station and not from the edges of the footprint. Also, the footprint area is an average of the entire season including day-and nighttime data, with higher fetch during night than during the day. An open-path infrared gas analyzer (IRGA) measured carbon dioxide and water density. Instantaneous high frequency sonic temperature and wind speed velocity components were measured with a 3-D sonic anemometer. The measurement height of IRGAs and sonic anemometers ranged between 200 and 500 cm, depending on crop height during the season. The F c were computed as the covariance between the instantaneous vertical wind velocity and carbon dioxide density. The F c was coordinate rotated [25] and corrected for temperature and air density variations [26]. The EC sampling rate (high frequency) was 20 Hz and was averaged to a 15-min time interval. Sonic and IRGA models and manufacturers are presented in Table 1.
A net radiometer measured incoming shortwave radiation (Table 1). An air temperature and humidity probe measured ambient air temperature (Ta) and relative humidity (rH) (Vaisala HMP45c, Vantaa, Finland). Precipitation was recorded using tipping bucket rain gauges (Model TR525 USW, Texas Electronics, Dallas, TX, USA). Soil moisture content was measured at 10 cm depth with soil moisture probes (Hydra Probe, Stevens Water Monitoring Systems, Portland, OR, USA). Thermocouples were buried at 2 and 6 cm to measure soil temperature, and heat flux plates (HFT-1, Radiation Energy Balance, Seattle, WA, USA) measured soil heat flux at 6 cm soil depth. An infrared temperature sensor (IRR, Apogee Instruments Inc., Logan, UT, USA) measured surface temperature and mounted at 45 • from perpendicular. The ancillary instrumentation was sampled at 0.1 Hz and averaged on a 15-min interval. The data were recorded on data loggers (CR5000, Campbell Scientific, Logan, UT). Further specifications of the EC systems are listed in Table 1.

Outlier Detection and Interpolation
The EC data were screened for statistical outliers and sensor failures and thereafter interpolated to fill data gaps [27]. Data was excluded under unfavorable weather conditions, such as during rain events, low wind turbulence and wind directions opposite to the sonic sensor direction [28], and time periods with high rH [29]. The internal sonic and IRGA warning flags and automated gain control were used to identify low-quality data. Data were screened for outliers using empirical set limits and the interquartile range (IQR) procedure, where outliers were defined as 1st/3rd Quartile ± 3.5 * IQR. Then the data were interpolated with an inverse weighted time average procedure [30]. The F c was separated for day-and nighttime (F c day , F c night ) and F c night (i.e., nighttime respiration) was interpolated using the relationship between air temperature and CO 2 fluxes [31], while F c day was gap filled with daily light response curves [32]. Approximately 31%-55% of the 15-min F c data was excluded, and F c night was more affected than F c day . Approximately 29%-42% of the F c data was interpolated. Daily values were calculated as the daily average multiplied by the number of samples per day.
Daily net ecosystem production (NEP) was calculated as the sum of daily F c day and F c night . Here the arithmetic sign convention is positive for F c from the atmosphere to the surface. The NEP was partitioned into GPP and RE. The GPP was calculated as the difference between NEP and RE. Daily RE was calculated as the sum of F c night and daytime respiration (Rd). Rd was derived from the non-linear Remote Sens. 2019, 11, 1688 6 of 23 relationship between incoming photosynthetic active radiation (IPAR) and F c day using light response curves, such as the rectangular hyperbola function [32]: where: F c day = daytime CO 2 flux (µmol m −2 s −1 ); AQY = apparent quantum yield (µmol CO 2 µmol PAR −1 ), P g max = maximum gross photosynthesis rate (µmol CO 2 m −2 s −1 ); IPAR = photosynthetic active radiation (µmol m −2 s −1 ); Rd = mean daytime respiration (µmol CO 2 m −2 s −1 ). IPAR was assumed to be 45% of incoming shortwave radiation [5] and was converted to photon flux density with 1 W m −2 = 4.6 µmol m −2 s −1 [33].
Daily NEP, RE, GPP and Rd values were screened again for outliers and interpolated using the inverse weighted time average procedure. Seasonal values were calculated as seasonal averages multiplied with number of days.

NDVI Ground-Based Measurements
Reflectance data on the ground were measured with a radiometer (CropScan, CROPSCAN, Inc., Rochester, MN, USA) in the corn and soybean field in Story County during the 2016 growing season. A total of 33 and 31 measurements were made on soybean and corn, respectively, in observation intervals of 2-14 days between DOY 124 and 308. Measurements were constrained to clear-sky conditions or with little cirrus clouds and with sun angles > 45 • (i.e., 10.00 and 14.00 CST). The radiometer measured reflected and incoming radiation at 16 wavebands between 460 and 1700 nm. Each measurement consisted of 30 individual scans, and an average is calculated by the unit. Measurements were taken above the crop at five locations north of the EC stations, and the measurement height was adjusted to yield a surface pixel size of 1.5 m. The radiometer was calibrated with a reflectance panel and a dark measurement prior to field observations. The Normalized Difference Vegetation Index (NDVI) was calculated as: where: ρ 660 and ρ 760 were the surface reflectance of the red (660 nm) and near-infrared (760 nm) wavelength bands, respectively. The individual observations at each location were averaged, and the averaged values interpolated with an inverse weighted time average procedure to yield a daily NDVI dataset for the 2016 growing season. Further information on CropScan measurements can be found in literature [34,35]. In addition to NDVI measurement, plant development and plant height was recorded (corn n = 21, soybean n = 40). The development stages are categorized in vegetative and reproductive growth as V and R Stages, respectively. The V-Stages in corn categorize plant emergence (VE), number of fully developed leaves (e.g., V1 for first leaf) and tasseling (VT), while the six R-Stages categorize corn flowering (R1), and ear development (R2-5) until maturation and senescence (R6) [36]. The V-stages in soybean categorize plant emergence (VE), first unifoliate leaf (VC), and number of fully developed trifoliate leaves on the main stem. In contrast to corn, reproductive stages occur during plant growth (V-stages), and categorize flowering (R1-2), pod development (R3-6), and maturation and senescence (R7-8) [37]. Plant height was measured from the ground to the top of the canopy, and daily height was calculated by fitting a three-parameter sigmoidal model. Due to resource restrictions, NDVI and intensive crop measurements were limited to the site in Story County.

MODIS fPAR Products, and Land Use Maps
The fraction of absorbed IPAR (400-700 nm wavelength), fPAR, can be derived from MODIS satellite images, and MODIS fPAR products (MCD15A2H) were provided by the U.S. Geological Survey [38,39]. The MODIS fPAR product had a resolution of 500 × 500 m and images are available at an 8-day interval. In total, 16 fPAR scenes between DOY 145 and 265 for the North Central Iowa Remote Sens. 2019, 11, 1688 7 of 23 counties were obtained. The quality of the fPAR scenes was checked with the Q/C layers provided by USGS, and areas of low quality were excluded.
The fPAR scenes were cropped to the 90% footprint area (i.e., two to four data points) of each EC station, as well as for the corn and soybean fields in each county. Corn and soybean fields were identified using the 2016 county land use maps ( Figure 1) provided by the Agricultural Conservation Planning Framework (ACPF) [18]. Only fields with a land use area >55% designated for corn and soybean production were used.

GPP Upscaling and Radiation Use Efficiency
The 8-day GPP MODIS product from MODIS satellite imagery (MOD17A2H) is calculated as [6]: where: GPP MODIS = MODIS gross primary production product (g C m −2 8-day −1 ), RUE max = maximum radiation use efficiency (g C MJ −1 ), Tmin s = a temperature scalar (no units), VPD s = a vapor pressure deficit scalar (no units), fPAR = fraction of IPAR absorbed by the crop (no units), IPAR = incoming photosynthetic active radiation (MJ m −2 d −1 ). Note that the GPP MODIS is calculated as 8-day sums, and daily values can only be calculated as daily means (Appendix B). The scalars Tmin s and VPD s are proxies for RUE max limitations of temperature and vapor pressure deficit (VPD), respectively. The scalars Tmin s and VPD s are calculated as [6]: where: Tmin s = scalar for daily Tmin (no units), Tmin = daily minimum Ta ( • C), Tlim min = lower Tmin limit ( • C), Tlim max = upper Tmin limit ( • C), VPDs = scalar for daytime VPD (no units), VPD = average daytime vapor pressure deficit (kPa), VPD min = lower VPD limit (kPa), VPD max = upper VPD limit (kPa). The Tlim min and Tlim max are set to −8 and 12.02 • C, while VPD min and VPD max are set to 0.65 and 4.3 kPa, respectively [6].
The GPP was upscaled from field to county following the same concept in Equation (3) (referred to as GPP calc ). First, RUE max for corn and soybean was calculated using the EC station data and NDVI field measurements from the Story County site. Equation (3) was solved for RUE max and fPAR substituted with field-measured corn and soybean NDVI, assuming that fPAR ≈ NDVI [6]. The VPD was calculated as the difference between saturated (e s ) and actual vapor pressure (e a ), of which the former was calculated as a function of Ta [40], and the latter as e s multiplied with rH. Then VPD s and Tlim s were calculated. RUE max was calculated as the slope of a linear regression (p < 0.05) with the product of Tlim s , VPD s , NDVI and IPAR as independent variable and EC-derived GPP as dependent variable [7]. At last, GPP calc was upscaled for the whole county using Equation (3), with MODIS fPAR products cropped for county-wide corn and soybean fields, respectively, and crop-specific RUE max . GPP calc was set to zero when fPAR < 0.4 (see Section 4).

Data Preparation, and Statistical Analysis
The variation of corn and soybean GPP among EC stations was quantified by calculating the 95% confidence intervals (CI) using the bootstrap percentile method with two thousand replications. Bootstrapping was run with the boot package in R [41]. Relationships between NEP, GPP, RE and Remote Sens. 2019, 11, 1688 8 of 23 climate parameters were evaluated with Pearson correlation and the regression coefficient (r) was calculated. The ACPF land use maps were validated with the 2016 county-wide corn and soybean production area [1] using linear regression analysis (p < 0.05) and cross-checked with the MODIS land cover type product (MCD12Q1) from 2016 [38,42]. The median EC-footprint and county-wide fPAR, GPP calc , and GPP MODIS were calculated for further analysis. The footprint fPAR was validated against field-measured NDVI with linear regression analysis (p < 0.05). GPP calc was validated against EC-derived GPP using a linear regression analysis (p < 0.05) and was also compared to the GPP MODIS product [5,38]. The median county-wide GPP calc was interpolated with an inverse weighted time average procedure to produce a daily dataset, and grain yield was calculated [7,12] as: where: Yield GPP = calculated yield data (Mg dry weight ha −1 ), GPP calc = seasonal sum of daily median county GPP calc (g C m −2 ), HI = Harvest index, RS = root-shoot ratio, M = % moisture content, OC = % grain carbon content. The HI was set to 0.50 and 0.57 [43,44], the RS was 0.18 and 0.15 [45], M was 15.5% and 13%, and OC was 44.7% and 54% [46], for corn and soybean, respectively. The Yield GPP was compared to the 2016 county yield data [1]. All figures were created with the R programming language [41].

Eddy-Flux Station Data
Eight EC stations measured carbon fluxes and meteorological data in fields with corn-soybean rotation during the growing season 2016. Mean daily IPAR decreased from 8.98 MJ m −2 on DOY 159 to 1.02 MJ m −2 on DOY 253 (Figure 2A). Mean Ta ranged from 16.4 and 28.8 • C ( Figure 2B) and mean (± SE) seasonal precipitation was 515 ± 51 mm, which is warmer and with near normal rainfall compared to the 30-year seasonal average (see Section 2.1). The average daytime VPD ranged between 0.13-2.40 kPa ( Figure 2C).

MODIS fPAR vs NDVI
The NDVI was measured in one corn and soybean field in Story County and compared to MODIS fPAR products. Field-measured daily corn and soybean NDVI ranged from 0.15-0.88 and from 0.16-0.91 ( Figure 3A) ------NEP = net ecosystem production, GPP = gross primary production, RE = ecosystem respiration, YieldGPP = yield as calculated with Equation (6).
Soybean NDVI reached values > 0.75 later than corn at DOY 180, corresponding to plant development stage V8-10/R1-2 and an average height of approximately 0.40 m. There was a significant linear relationship (p < 0.05) between field-measured NDVI and MODIS fPAR within the 90% footprint area of the EC station, explaining 99% of the fPAR variation with a slope of 1.02 and 0.99 for corn and soybean, respectively ( Figure 3B).
The RUEmax was then calculated with Equation (3), substituting fPAR with NDVI and using meteorological data and GPP from the EC stations in the studied corn and soybean field in Story County. The RUEmax was the slope of a linear regression with GPP as dependent and the product of Tmin, VPDmin, IPAR, and NDVI as independent variables. The derived daily corn and soybean RUEmax was 3.24 and 1.90 g C MJ -1 . The linear relationship was significant at p < 0.05, explaining 93% and 89% of the GPP variation for corn and soybean, respectively ( Figure 3C).  (B) linear regression analysis (p < 0.05) with measured NDVI as independent and MODIS fPAR in the 90% footprint area as dependent variable; (C) linear regression analysis (p < 0.05) with the product of the Tmin s , VPD s , NDVI, and IPAR (MJ m −2 d −1 ) as independent and GPP (g C m −2 d −1 ) as dependent variable with the slope being RUE max . R 2 : coefficient of determination, Tmin s : temperature scalar, VPD s : vapor pressure deficit scalar, IPAR: incoming photosynthetic active radiation, GPP: gross primary production, RUE max : maximum radiation use efficiency.
Soybean NDVI reached values > 0.75 later than corn at DOY 180, corresponding to plant development stage V8-10/R1-2 and an average height of approximately 0.40 m. There was a significant linear relationship (p < 0.05) between field-measured NDVI and MODIS fPAR within the 90% footprint area of the EC station, explaining 99% of the fPAR variation with a slope of 1.02 and 0.99 for corn and soybean, respectively ( Figure 3B).
The RUE max was then calculated with Equation (3), substituting fPAR with NDVI and using meteorological data and GPP from the EC stations in the studied corn and soybean field in Story County. The RUE max was the slope of a linear regression with GPP as dependent and the product of Tmin, VPDmin, IPAR, and NDVI as independent variables. The derived daily corn and soybean RUE max was 3.24 and 1.90 g C MJ −1 . The linear relationship was significant at p < 0.05, explaining 93% and 89% of the GPP variation for corn and soybean, respectively ( Figure 3C).

Upscaling GPP to County Level
The total corn and soybean area from the ACPF land use maps corresponded well with the official production area in 2016 [1] (y = 1.01x; R 2 = 0.99, p < 0.05, n = 8). The land use maps also corresponded well with the MODIS land cover in 2016, and >99% of the ACPF map area used was classified as cropland in the MODIS land cover map. The average (± SE) ACPF land use area of corn and soybean in 2016 made up 55% ± 2% and 26% ± 1% of the total land area in four counties, respectively.
The average within-county fPAR range (difference between maximum and minimum fPAR) was 0.65 at DOY 145, decreased to 0.32 at DOY 201 and increased again to 0.65 at DOY 265 ( Figure 4, Video S1). The county median fPAR ranged from 0.25-0.9 with little variation among crops and counties during the 2016 growing season ( Figure 5A). Calculated daily corn and soybean GPP followed a seasonal pattern with highest values at 22.1 and 11.7 g C m −2 d −1 , respectively ( Figure 5B). A linear regression between EC-derived GPP and GPP calc within the 90% footprint was significant at p < 0.05, with a coefficient of determination (R 2 ) of 0.96 and 0.86 and a root mean square error (RMSE) of 2.24 and 2.81 g C m −2 d −1 , for corn and soybean, respectively ( Figure 5C). A linear regression between EC-derived GPP and daily GPP MODIS within the 90% footprint was significant at p < 0.05, with a R 2 of 0.89 and 0.21 and a RMSE of 2.44 and 3.30 g C m −2 d −1 for corn and soybean, respectively. GPP MODIS was not substantially different between crops, overestimated GPP at the onset and end of the season, and underestimated corn GPP (Figures A2 and A3).
The average within-county corn and soybean GPP calc range (difference between maximum and minimum GPP calc ) was 13.2 and 6.0 g C m −2 d −1 at DOY 145, decreased to 1.5 and 0.9 g C m −2 d −1 at DOY 201 and increased again to 5.1 and 3.6 g C m −2 d −1 at DOY 265 ( Figure 6). The Yield GPP was 12.82 ± 0.65 and 6.73 ± 0.27 Mg ha −1 , for corn and soybean, respectively ( Table 2). The official county-wide corn and soybean yield in 2016 was 13.09 ± 0.09 and 4.03 ± 0.04 Mg ha −1 , respectively [1].
The total corn and soybean area from the ACPF land use maps corresponded well with the official production area in 2016 [1] (y = 1.01x; R 2 = 0.99, p < 0.05, n = 8). The land use maps also corresponded well with the MODIS land cover in 2016, and >99% of the ACPF map area used was classified as cropland in the MODIS land cover map. The average (± SE) ACPF land use area of corn and soybean in 2016 made up 55% ± 2% and 26% ± 1% of the total land area in four counties, respectively. regression between EC-derived GPP and GPPcalc within the 90% footprint was significant at p < 0.05, with a coefficient of determination (R 2 ) of 0.96 and 0.86 and a root mean square error (RMSE) of 2.24 and 2.81 g C m -2 d -1 , for corn and soybean, respectively ( Figure 5C). A linear regression between ECderived GPP and daily GPPMODIS within the 90% footprint was significant at p < 0.05, with a R 2 of 0.89 and 0.21 and a RMSE of 2.44 and 3.30 g C m -2 d -1 for corn and soybean, respectively. GPPMODIS was not substantially different between crops, overestimated GPP at the onset and end of the season, and underestimated corn GPP (Figures A2 and A3).  . GPP calc as calculated with MODIS fPAR product for corn (black circles) and soybean (red triangles) for 2016: (A) median county fPAR, averaged (± SE) over county (n = 4), (B) daily median county GPP calc (g C M −2 d −1 ) averaged over county, and (C) linear regression analysis (p < 0.05) with daily EC-derived GPP as independent and 90% footprint GPP calc (g C m −2 d −1 ) as dependent variable, R 2 : coefficient of determination, GPP: gross primary production, fPAR: fraction of absorbed photosynthetic active radiation.

Discussion
In this study, the 2016 corn and soybean GPPcalc in four counties were upscaled using ground measurements of GPP, NDVI, and meteorological data as well as MODIS fPAR products. In contrast to previous studies, seasonal GPP, NEP, and RE from several EC stations concentrated on a relatively small land area were retrieved. It was therefore possible to estimate the spatial variation of GPP under similar climate conditions and for a relatively small land area as compared to previous studies [12][13][14][15]. The calculated 95% CI indicate a substantial variation in GPP, NEP, and RE among sites, despite

Discussion
In this study, the 2016 corn and soybean GPP calc in four counties were upscaled using ground measurements of GPP, NDVI, and meteorological data as well as MODIS fPAR products. In contrast to previous studies, seasonal GPP, NEP, and RE from several EC stations concentrated on a relatively small land area were retrieved. It was therefore possible to estimate the spatial variation of GPP under similar climate conditions and for a relatively small land area as compared to previous studies [12][13][14][15]. The calculated 95% CI indicate a substantial variation in GPP, NEP, and RE among sites, despite that the studied fields are well represented in the four counties [5] and the relatively high density of EC stations used. Individual spatial variation in GPP is probably connected to climate, soil and crop management ( Table 1, Table 2), such as precipitation, cloud cover, crop growth, crop varieties, fertilizer application, tillage, and amount of crop residues among others. Analyzing for these impacting factors was beyond the scope of this study and would require more information of each field but may be important to further increase the accuracy of GPP calc . In a previous study, substantial variation in corn and soybean F c and precipitation using 12 EC stations in a 5000-ha Midwestern US watershed was found, and this was accounted to difference in ground cover, plant growth and biomass accumulation among fields [47]. Differences in soil tillage has also impacted GPP [4], and soil respiration varies with nitrogen fertilizer application, irrigation, and tillage [48]. This indicates that a higher density of EC stations than typically used in field studies is needed to cover the spatial variation of GPP calc . The higher NEP values in corn are connected to higher photosynthetic rates of corn as a C4 crop, compared to soybean as a C3 crop. In addition, high amounts of corn residues from the previous season remain in the soybean fields, which further decreases NEP due to higher residue mineralization [27].
The GPP calc was more accurate than GPP MODIS , as evidenced by lower RMSE in comparison to EC-derived GPP. Field measurements of RUE max , improved land use maps, and a higher resolution of IPAR, Ta, and VPD measurements increased the accuracy of GPP calc compared to GPP MODIS .
One important tool to increase the accuracy of GPP calc compared to GPP MODIS was the ACPF land use map [18], which enabled distinction between corn and soybean fields within each county. The ACPF land use maps are based on an ensemble of aerial photographs, satellite imagery, and other land use map products [18]. Earlier studies on corn and soybean production in the Midwestern US were constraint to the Cropland Data Layer (CDL) and an area-weighted, crop-specific value was used to estimate GPP contribution by crop [12,13]. To the best of our knowledge, this is the first time that ACPF land use maps have been used to upscale GPP calc . The field boundaries provided by ACPF land use maps were the basis to create GPP calc maps as separated by crop. Approximately, 3 4 of the total land area in the studied counties were used for corn and soybean production, which demonstrates the economic importance and its impact on the carbon cycle in the region. There were more fields used for corn than for soybean production in the studied counties in 2016, and corn GPP (as a C4 crop) was higher than soybean GPP (as a C3 crop) ( Figure 2D). Consequently, the tempo-spatial variation in corn and soybean fields determines the amount of county-wide, seasonal cropland GPP calc . Accurate and annually updated land use maps are therefore imperative to upscale GPP calc and is important to improve GPP MODIS .
Another significant parameter to improve GPP MODIS is a crop-specific RUE max , which determines the maximum amount of GPP per unit of absorbed IPAR (i.e., fPAR × IPAR). Similar corn and soybean RUE max of 3.17 and 1.91 g C MJ −1 had previously been reported using EC data and MODIS fPAR in corn and soybean fields in the Midwestern US [12]. For the GPP MODIS product the RUE max is determined according to land cover categories (from MODIS MCD12Q1 land cover types), and a cropland RUE max of 1.044 g C MJ −1 was proposed [6]. This is lower than the corn and soybean RUE max in this study and does not distinguish between both crops. This consequently leads to the substantial underestimation of GPP MODIS , especially for corn as a C4 crop ( Figures A2 and A3, Video S2, Video S3).
The MODIS fPAR product within the EC footprint corresponded well with NDVI measurements on the ground for both crops, which is interesting considering the difference in resolution (i.e., 500 × 500 m satellite image vs. 1.5 m × 1.5 m pixel size of field measurement). It also indicates that MODIS fPAR displays the NDVI of both crops. It was also noted that Equation (3) overestimates GPP calc at the beginning and the end of the growing season, where the EC-measured GPP was zero ( Figure A2). This problem was somewhat resolved by setting GPP calc to zero, when fPAR < 0.4 (i.e., early and late in plant development). A field-measured NDVI value of 0.15 was found for bare soil in corn and soybean, similar to this study [28], and GPP is zero. It is likely that NDVI and MODIS fPAR do not represent GPP calc well in the early and late stages of plant development, despite that NDVI was found to be a good indicator of intercepted PAR [35]. A field-based fPAR as well as other MODIS imagery-derived VIs further improved GPP calc in corn and soybeans [14,15].
The importance of RUE max , fPAR, and land use maps should not diminish the importance of accurate meteorological data, as it determines the amount of energy intercepted by the crop (IPAR) and to what extent RUE is limited by temperature and VPD constraints (i.e., Tmin s and VPD s scalars). This eventually determines the amount of GPP calc . For example, on DOY 201, the median-county GPP calc was low, because IPAR was low at all eight EC stations (Figures 2A and 5B). Although the variation of IPAR, Ta, and VPD was low among EC stations (Figure 2), these variables impacted among-county daily GPP calc variation. For example, corn GPP calc on DOY 193 was higher in Hardin county owing to slightly higher VPD s and IPAR (Figure 6, Video S3). Political boundary delimitations (i.e., counties) are not connected to these climatic changes, and further methods are needed to model the spatial variation of IPAR, Ta, and VPD with higher resolution. Within-county variation of corn and soybean GPP calc is therefore only driven by the MODIS fPAR products in this study. In this study, IPAR, Ta, and VPD from ground measurements were used because of the relatively high resolution among stations. However, this may not always be the case. For example, the global GPP MODIS product is calculated with meteorological data provided from the NCEP-DOE reanalysis II datasets, which have a coarser resolution than the density of EC stations used in this study [49].
The seasonal sum of county-median GPP calc was used to estimate Yield GPP , which was compared to actual county yield data. While corn Yield GPP corresponded well with county corn yield in 2016, soybean Yield GPP was overestimated ( Table 2). There has been good agreement between Yield GPP and actual yield values in corn and soybean [13], and wheat [7], however, both studies found limitations of the approach. Yield GPP as shown in Equation (6) depends on more factors (HI, M, OC, and RS may differ among fields) than GPP calc and this estimation only represents a general trend. The overestimation of soybean Yield GPP could be connected to the mineralization of the 2015 corn residues in soybean-years in these corn-soybean rotation systems. The mineralization of corn residues from the previous growing season can have a substantial impact on current-year soybean NEP, GPP, and RE [4]; thus, the GPP calc as combined value of corn residue respiration and soybean assimilation may underrepresent soybean yield.

Conclusions
It was shown in this study, that GPP calc can be upscaled with higher accuracy (i.e., lower RMSE) than GPP MODIS . The GPP MODIS product does not distinguish among crops, which leads to a substantial underestimation of C4 crops such as corn. The fPAR MODIS product agreed with ground-measured NDVI. The meteorological data showed low variation among EC stations, indicating that a coarse-resolution meteorological dataset may be sufficient to upscale regional bulk GPP in this region of Iowa. However, variation of soil and crop management as well as micro-climate may require a higher density of EC and weather stations to estimate GPP calc on a field level. A substantial improvement of upscaling regional GPP can be achieved using crop-specific RUE max and improved land use maps. A more detailed categorization of MODIS land use maps (e.g., cropland C4 crops vs. cropland C3 crops), a crop-specific RUE max , and a different VI for the beginning and end of the growing season can increase the accuracy of GPP MODIS without the need to modify the algorithm or substantial re-processing of imagery.

Conflicts of Interest:
The authors declare no conflict of interest.