Cropland Net Ecosystem Exchange Estimation for the Inland Pampas (Argentina) Using EVI, Land Cover Maps, and Eddy Covariance Fluxes

Estimations of Net Ecosystem Exchange (NEE) are crucial to assess the carbon sequestration/carbon source capacity of agricultural systems. Although several global models have been built to describe carbon flux patterns based on flux tower data, South American ecosystems (and croplands in particular) are underrepresented in the databases used to calibrate these models, leading to large uncertainties in regional and global NEE estimation. Despite the fact that almost half of the land surface is used worldwide for agricultural activities, these models still do not include variables related to cropland management. Using enhanced vegetation index (EVI) derived from MODIS imagery (250 m) and monthly CO2 exchange from a 9-year record of an eddy covariance (EC) flux tower in a crop field in the Inland Pampas region, we developed regression models to predict monthly NEE. We tested whether including a term for crop identity/land cover as a categorical variable (maize, soybean, wheat, and fallow) could improve model capability in capturing monthly NEE dynamics. NEE measured at the flux tower site was scaled to croplands across the Inland Pampa using crop-type maps, from which annual NEE maps were generated for the 2018–2019, 2019–2020, and 2020–2021 agricultural campaigns. The model based solely on EVI showed to be a good predictor of monthly NEE for the study region (r2 = 0.78), but model adjustment was improved by including a term for crop identity (r2 = 0.83). A second set of maps was generated taking into account carbon exports during harvest to estimate Net Biome Productivity (NBP) at the county level. Crops across the region as a whole acted as a carbon sink during the three studied campaigns, although with highly heterogeneous spatial and temporal patterns. Between 60% and 80% of the carbon sequestered was exported during harvest, a large decrease from the carbon sequestration capacity estimated using just NEE, which further decreased if fossil carbon emissions from agricultural supplies are taken into account. Estimates presented in this study are a first step towards upscaling carbon fluxes at the regional scale in a South American cropland area, and could help to improve regional to global estimations of carbon fluxes and refine national greenhouse gas (GHG) inventories.


INTRODUCTION
Estimations of Net Ecosystem Exchange (NEE) are crucial to assess the carbon sequestration or source capacity of agricultural systems (1). NEE of terrestrial ecosystems is quantified by the difference between CO 2 uptake via gross photosynthetic assimilation (GPP) and release of CO 2via ecosystem-wide respiration processes (R eco ) (2). GPP and net primary production (NPP) often used by ecologists do not fully describe the feedback of the ecosystems-atmosphere carbon exchange as NEE does, because they do not explicitly include soil-derived fluxes and other fluxes associated with heterotrophic respiration (3). In managed agricultural ecosystems, estimates of carbon balance have to be complemented NEE with "lateral" fluxes such as crop harvest, deriving what is known as Net Biome Productivity (NBP) (4), also called Net Ecosystem Carbon Budget (NECB) (5). To identify whether a managed ecosystem is a net source or a destination of carbon, it is necessary to calculate its NBP, since many times NEE suggests that the ecosystem is sequestering carbon when in fact it is acting as a net source (6,7).
During the last three decades, measurements of NEE between ecosystems and the atmosphere via direct quantification of net carbon flux have been carried out using the eddy covariance technique (EC; 8). Currently, there are at least 212 unevenly distributed sites around the world with EC micrometeorological towers included in the FLUXNET global network, and its FLUXNET2015 database, providing measurements of exchange of CO 2 , water vapor, and energy, plus associated biological and meteorological variables (9). Although the effective area of the EC measurements, its "footprint", is relatively small, often less than 1 km 2 , the variability of NEE at these sites is representative of variability at larger spatial scales (10). The combined use of the FLUXNET2015 database, spatially explicit climate information, and satellite products such as vegetation indices has enabled the extrapolation of carbon fluxes from the local scale to the regional, continental, and global scales (10)(11)(12)(13).
Although several global models have been built to describe carbon flux patterns based on EC data and vegetation indices (10,(13)(14)(15)(16), South American ecosystems (and croplands in particular) are still underrepresented in the databases used to calibrate these models, leading to large uncertainties in the regional and global estimatesof NEE. Global models show less than 50% efficiency in their ability to predict NEE (14). This may be attributable, at least in part, to approximately 85% of the FLUXNET sites used to calibrate these models being located in the northern hemisphere, mainly between 30°C and 55°C of latitude (17). Croplands cover about 12%-13% of the Earth's icefree land surface (18,19), and even when ca. 10% of FLUXNET sites are croplands, none of them is located in the southern hemisphere (9).
Human interventions are one of the key non-climatic controls of the carbon balance (13,20). Despite the fact that almost half of the land surface is already used worldwide for agricultural activities (livestock and crop production, 18), the inclusion of management-related variables in models remains a challenge (14), specially for croplands (15). In recent studies, variables related to leaf-area temporal dynamics were included as a proxy for the seasonality of different plant functional types (PFT) (16), but this is done using land cover maps derived from MODIS products (i.e., MCD12Q1), whose maximum level of detail is the "crop" category, thus masking the large heterogeneity of these managed ecosystems (11). To overcome this constraint for agricultural landscapes, accurate and detailed land cover maps for different seasons along the year are needed to improve the estimates of regional fluxes (21). Recently, several crop-type maps have been published for our study region as part of the MapBiomas project that can be used to include management in a regional model (22,23).
In this study, we aimed to develop a regression model to predict monthly NEE, combining enhanced vegetation index (EVI) derived from MODIS satellite imagery (250 m spatial resolution) with monthly CO 2 exchange from a 9-year record of an eddy covariance (EC) flux tower in a crop field in the Inland Pampas region of Argentina. We tested whether including a term for crop identity/land cover (maize, soybean, wheat, and fallow) as a categorical variable could improve the model ability for capturing monthly NEE dynamics. NEE measured at the flux tower site was scaled to croplands across the Inland Pampa (a subregion of the Pampas with similar soil and climate conditions) using the regression model, and annual NEE maps were generated for the 2018-2019, 2019-2020, and 2020-2021 agricultural campaigns (June to May is the usual 12-month period used to center the crop-year around December to February summer months). Spatially explicit crop yield data from national statistics at the county level was combined with NEE maps to estimate regional and county NBP.

Site Description
The flat Inland Pampa is part of the Pampas region ( Figure 1, 24), one of the most agriculturally productive areas worldwide and responsible for most of Argentina's total grain exports. The study site had a 3-year continuous crop rotation including soybean, wheat-soybean double cropping, and maize, representative of the region (22,25). This EC site is located in a 32-ha field in Carlos Casares, Buenos Aires Province, Argentina (35°37'15,52" S, 61°19'5,26" W, Figure 1). The climate is temperate sub-humid, with an average annual rainfall of 1,022 mm during the last 25 years, and mean temperatures ranging from 7.2°C in July to 23.8°C in January (26). Livestock and farming activities are carried out in different soils (27). Agriculture is concentrated in typic Hapludols, well drained and rich in organic matter (28). In deep sandy-loam soils, it is possible to find agriculture or a mesophyte pseudo-steppe with abundant grasses. In soils of flat and low areas close to lentic environments, there are halophyte steppes and livestock activities (24).

Micrometeorological Data
Fluxes of CO 2 were measured by the EC method between November 2012 and December 2020. The EC setup comprised a 3-D sonic anemometer (CSAT3; © Campbell Scientific, Logan, UT, USA) and an open path infrared gas analyzer (IRGA; LI-7500; © Li-Cor Inc., Lincoln, NE, USA), for measuring wind speed and sonic temperature and CO 2 concentrations, respectively, 2.5 m above the canopy. From high-frequency (10 Hz) EC data, fluxes were computed in 30-min average blocks with standard procedures (29) as the covariance between the concentration of CO 2 and the vertical component of wind speed. Invalid data (e.g., nighttime fluxes under non-turbulent conditions) were removed and gap filling was carried out following Reichstein et al. (2005) (30). A more detailed description of the methodology is found in Posse et al. (2018) (31). NEE was aggregated on a monthly basis (32). A footprint model (33) shows that measures cover about ca. 86% of the area of the crop field. Footprint area was defined as encompassing the 70% of accumulated total footprint area, covering ca. 30 ha (310 m average radius around the tower). By convention, positive flux values represent mass transfer from the surface to the atmosphere and negative values indicate the opposite.

Remote Sensing Data
Two regression models to estimate NEE were derived from EVI from NASA-MODIS. This index is more accurate than NDVI in areas with a dense canopy, is less sensitive to canopy background signal and atmosphere influences (34), has a better response to structural changes of deciduous canopies (35) and tends to show better correlation with croplands at different time scales (1,36). For the entire Flat Inland Pampa, EVI (MOD13Q1 MODIS product, 250 m spatial resolution) was obtained for each 16-day interval over the period November 2012-December 2020 from the Google Earth Engine repository (37). For the EC site, three pixels matching the footprint and centered at tower coordinates were chosen, and the monthly mean EVI for each of them was obtained and averaged [arithmetic mean as, e.g (38)., and (39) did]. Averaging over many pixels reduces gridding errors, the effects of geo-location error, and random noise (40). For regional analysis purposes, monthly mean EVI values of each pixel within the entire Flat Inland Pampa were obtained (32,41).

Model Development
Two linear regression analyses were performed to investigate the relationship between monthly averaged site EVI and monthly accumulated NEE. Within the general regression, we found statistically significant different slopes between different crops; thus, a second regression model was built including terms for land cover as categorical dummy variables. In a first attempt, categories were maize, soybean, wheat, and fallow. Fallow periods often exhibit carbon losses, and its inclusion in carbon balances is recommended for comprehensive quantification (42). Including fallow cover was problematic because these periods have low EVI values and mainly positive NEE values (around 0.2 and 0 to +100 gC m −2 , respectively; Figure 2). Thus, in order to include these during model development, we made two assumptions: (a) Fallow months were assigned to the preceding crop cover, the one responsible for the decomposing stubble; and (b) the situation during fallow is very similar regardless of the crop (low EVI and NEE are mainly attributed to respiration). Thus, we assumed that all crop covers had the same Y-intercept value. Therefore, mutually exclusive categories included in the land cover model were reduced to three: maize, soybean, and wheat. A contrast matrix ( Table 1) was built with 0's and 1's according to the absence or presence of each crop cover, respectively. For example, when categorical variables "soybean" and "wheat" have a value of 0, meaning the absence of these crop covers, the regression equation corresponds to maize cover ( Table 1). Months with no changes in land cover were assigned to the corresponding category. Sowing months were classified as fallow and harvest months were classified as the crop harvested, unless the change (sowing or harvest) involved more than 70% of the length of the month. The latter occurred in less than 4% of the months.
For both models, a 10-fold cross-validation was performed using the caret R package version 6.0-90 to assess the quality of the estimators. The fit criterion was characterized using the coefficient of determination (r 2 ) and the root mean square error (RMSE). The model assumptions were evaluated for normality and autocorrelation of the residuals and homoscedasticity. The general, aggregated-crop model slightly failed them but the land cover model improved this aspect (Shapiro-Wilk normality test p > 0.05; no residuals autocorrelation).

NEE Upscaling and NBP Maps
Land cover was determined using crop-type maps for the 2018-2019, 2019-2020, and 2020-2021 agricultural campaigns (43-45). These maps describe the presence of the main single and double crops from June to May (i.e., southern hemisphere winter to fall), based on field surveys and supervised classifications with a 30-m resolution and accuracy higher than 85%. For the 2019-2020 and 2020-2021 campaigns, winter and summer crop maps have been published separately, so we combined them into a single map for each campaign. We used the land cover categories included in the land cover model (maize, soybean, and wheat) masking the rest. Thus, our study is not an extrapolation of the carbon balance to the whole region, because we are not making inference on pixels not included in the NEE calculation, some of which correspond to livestock activities. Crop maps did not discriminate between winter cereals [wheat, barley, rye, oats, and winter fodder grasses (44)]. However, as most of the winter cereal sown area in this region is wheat [National Crop Statistics (46)], we assumed that all pixels included in that category behaved as wheat (6).
By overlaying crop-type maps with the monthly mean EVI maps, we selected MODIS pixels that contained only one of our three crops included in the model ("pure" pixels). On these selected MODIS pixels, the land cover model was used to estimate NEE. Pure pixels belonging to maize, wheat, and soybean represented together 21%, 28%, and 33% of the agricultural area mapped in the 2018-2019, 2019-2020, and 2020-2021 crop maps, respectively. The remaining cropland area corresponds to "impure" pixels that contain other crops, more than one crop cover, or combinations of crops with non-cultivated land. Fallow NEE was estimated as stated above with maize and soybean cover from the previous summer. Since we lacked these data for the 2018-2019 campaign, we assumed that if soybean was grown in the 2018-2019 campaign, maize had previously been grown in 2017-2018 and vice versa, which occurs in more than 80% of the Rolling Pampas cropland area, next to the Inland Pampa (25). Crop yield data from national statistics at the county level [National Crop Statistics (46)] and an average carbon content of 0.45 (47) were used to calculate carbon exported per crop and NBP.  Carbon footprint for each agricultural campaign, weighted by the area of each crop, was estimated combining a region-wide average NBP in MgC/ha and greenhouse gas (GHG) emissions per Mg of grain harvested from carbon inputs in fertilizers, fuels, machinery, and pesticides. For soybean and maize, GHG emissions of 0.150 MgCO 2 -eq/Mg grain and 0.172 MgCO 2 -eq/ Mg grain were obtained from (48), as the mean GHG emissions for the Inland Pampas region. For winter wheat, we used a mean value of 0.45 MgCO 2 -eq/Mg (49,50).
NEE gC m -2 month -1 À Á ¼ 104:83 -390:57 Ã EVI (Eq:1) The land cover model (Eq. 2) relates monthly mean site EVI values to monthly NEE distinguishing between maize, soybean, and wheat slopes and assuming the same Y-intercept (see Materials and Methods) (n = 95, p < 0.001, r 2 = 0.83, RMSE: 39.3, Figure 2B and Table 2). Parameter estimation using K-fold cross-validation (k = 10) gives the same results and a higher r 2 (0.85). The slope for wheat does not differ from that of maize ( Besides the land cover model that includes all three cover types, monthly NEE for maize, soybean, and wheat can be individually predicted through their monthly mean EVI. In isolation, all three showed that they have significant and different slopes and Y-intercept value ( Figure 3). Soybean reaches higher EVI values than maize and wheat. Maize reaches NEE values equal to those of soybeans at lower EVI values.
Both models mirrored the temporal dynamics of observed data ( Figure 4). They were able to establish, with an efficiency of 85%, whether the field behaved as a carbon source or sink. For values close to zero, between ±40 gC m −2 month −1 (in line with model RMSE), such efficiency drops to 76%. Both models underestimated the peaks of carbon loss above 80 gC m −2 month −1 that mostly occur after the harvest of summer crops, maize, and soybean. The land cover model provided an overall slightly better description of NEE monthly dynamics. Yearly, cumulative agricultural campaigns NEE (June to May, Table 3) from observed monthly data were predicted with reasonable accuracy by both models (RMSE: land cover model: 120 gC m −2 year −1 ; general model: 160 gC m −2 year −1 ). They seemed to compensate month-to-month residuals leading to a close matching between 9-year NEE and the observed data.

NEE Upscaling and Inland Pampa Cropland NBP
Estimates made with this approach show that croplands of the region as a whole acted as a carbon sink during the three agricultural campaigns (Table 4) At the regional scale, NBP was 80%, 62%, and 78% higher than NEE for the 2018-2019, 2019-2020, and 2020-2021 campaigns, respectively, due to carbon export during harvest ( Table 4). Among counties, NBP showed a large spatial and temporal heterogeneity. Figures 5D-F show that the strong net carbon uptake in croplands for all seasons was mainly distributed from the center towards the southeast of the study area in all three agricultural campaigns. This may have been explained, at least in part, by the high percentage of area with double cropping (Figures 5G-I). Throughout the three campaigns, counties seemed to show trajectories of NBP relatively independent of each other. Carbon footprint analysis shows that the actual C balance (NBP minus input fossil emissions) decreases between 20% and 47% when inputs are taken into account ( Table 5).
NEE spatial patterns seem to exhibit scale dependency, as net fluxes can change from broader to fine scales. Figures 6A, B show the detailed heterogeneous spatial distribution of NEE in   Figure 5E), but includes sites with an NEE of more than −900 gC m −2 year −1 , equivalent to 9 MgC ha −1 year −1 . Temporal dynamics of the monthly averaged NEE of the Inland Pampa were similar for the three campaigns (Figure 7). Carbon sequestration seems to be concentrated in the summer months and carbon losses in the fall and winter. During the spring, the study region seems to function on average as a weak sink.

DISCUSSION
Both models presented in this study were able to capture the temporal dynamics of measured NEE through crop leaf-area phenology. The land cover model provided a better description of the timing of carbon sequestration, perhaps because it captures more variability by differentiating the slopes of each crop. This improvement in NEE estimation by explicitly including crop identity was also observed in more complex models with shorter temporal resolution and coarser spatial resolution (51). Both pixel-based models presented in this paper allow direct estimation of NEE from satellite data and are thus easy to implement with regularly available EVI data, as they do not need external inputs such as meteorology, light use efficiency (LUE) or FPAR (1,52). Use of coarser resolution weather, biophysical, and ecophysiological data to estimate GPP may introduce significant errors (41), although it remains unclear how the EVI-NEE relationship varies with time scales (32) and fine spatial resolutions (41).
Underestimation of carbon losses during the post-harvest months in both models may be explained by the fact that fallow months have a large variability in their NEE over a narrow range of EVI values. EVI cannot differentiate between different fallow situations because it is only sensitive to photosynthetic activity. A few studies (41,53) have proposed using surface-temperature satellite products (such as LST) to overcome this issue, since soil temperature is closely related to R eco . Periods without photosynthetic activity generate difficulties in estimating FIGURE 4 | Temporal dynamics of observed (gray line) and predicted NEE with both models. The yellow line represents general model prediction and the green line shows land cover model prediction. Gray dashed lines around zero enclose the area between ±40 gC m −2 month −1 where efficiency to determine whether the field behaved as a carbon source or sink drops to 76%. carbon fluxes with a remote sensing approach; thus, removing winter months tends to improve EVI-derived GPP estimates (36,53). That may work for northern hemisphere sites affected by snow or cold temperatures. However, in the milder southern hemisphere temperate croplands, fall and winter months without photosynthetic activity can have a strong impact on annual carbon balance (31,54) and not including them may lead to think that the ecosystem behaves like a carbon destination, when in fact it is acting as a net source (6,7).
The spatial pattern we observed for NBP is partly explained by the percentage of area under double cropping within each county. In turn, it is possible that it also depends on the yields of the different crops (Supplementary Material S1). Wheat yields are higher in the southern counties than in the northern ones, where soybean and maize yields are higher than in the south. Soybean and maize early single-cropping, sown in November, usually have higher yields than when it is preceded by wheat or another winter crop (46). These yield differences would encourage, at least in part, wheat and soybean double-cropping in the south and early maize and soybean single-cropping further north. Other factors that may have influence but were not taken into account in this study are farmer decisions (as to which species to plant and how much to fertilize them) based on changes in international prices, land tenure regime, agricultural policies and weather conditions (25).
Estimates of carbon fluxes are highly dependent on both spatial and temporal scales, with spatial prediction being usually more difficult (55). According to our estimations, the Inland Pampa croplands seem to be acting as a carbon sink as a whole, even though the spatial distribution of NEE and NBP at the county and pixel-site scales is highly heterogeneous ( Figure 5, Figure 6). Despite the fact that NEE and NBP estimations of heterogeneous regional ecosystems are relevant for assessing ecosystem services in land-use planning (56,57), studies comprising upscaling at the regional scale are scarce. Temporal variations in croplands are highly variable and partially dependent on management (58), but understanding and predicting its monthly mean dynamics seem to be less challenging than spatial patterns. We propose that estimates made for a full rotation or at least a full year are more realistic and comparable than those made for a single growing season (6, 7).
Our regional estimates of NEE are roughly in agreement with those for other cropland regions of the world, although somewhat higher (11,52,59). Overestimation of carbon fluxes from croplands have been reported in upscaling models, particularly the ones including maize and soybean (11,51,60). This study is based on pure pixels of the most productive crops (46) in a high-yielding agricultural region (61) and coupled with a possible underestimation of carbon losses during fallow months as stated above. Therefore, it is possible that we have overestimated their carbon sink capacity. Furthermore, using a single average crop yield value per county (best information available) in sites with NEE values much lower than the average may have resulted in an overestimation of the NBP.
Between 60% and 80% of carbon sequestered was exported during harvest on a regional scale, generating a large decrease in the carbon sequestration capacity estimate expressed using just NEE. This reinforces the idea that in order to evaluate the complete carbon balance of a managed ecosystem, it is necessary to estimate its NBP (7). NEE indicates carbon sequestration capacity but the carbon absorbed by crops and allocated to grains will be released back into the atmosphere due to grain and sub-product consumption by humans and cattle.  Ecosystems with high carbon exchange rates such as the one presented in this study might play an important role in ecosystem services supply related to food provisioning (57,62,63). Thus, estimating their NBP seems a more appropriate way to assess the potential carbon remaining in their soils.
If GHG emissions from carbon inputs in fertilizers, fuels, machinery, and pesticides are taken into account, the average carbon balance for croplands comprised in this study becomes closer to neutrality. However, a small carbon sink still remains for the three agricultural campaigns. Estimates presented here  suggested that GHG emissions from fossil carbon inputs can represent up to nearly half of NBP in croplands. These results are not in agreement with studies in other cropland regions of the world (58,64,65) showing that croplands act as a net source of GHGs. Carbon footprint estimates presented here are a coarse regional average and more detailed estimations that take into account spatial internal heterogeneity are clearly needed. Nevertheless, two recent studies in the region suggest that Pampean cropland soils might have a considerable carbon sequestration potential (66,67), which opens the way for devising effective mitigation. GHG emissions from agriculture (crops and livestock) are estimated to be growing by 1% annually and surpassing land use change emissions (68). In addition, management practices in croplands can have a significant impact on climate (69)(70)(71)(72)(73). Estimates presented in this study are a first attempt made in upscaling carbon fluxes at the regional scale in a South American cropland area, showing the carbon sequestration capacity of croplands in the Inland Pampa. Local/regional models and upscaling processes such as the one presented in this study could become an important tool both for local environmental assessment and food-security policies, and globally as a benchmark for global models, given their uncertainties. Consideration of data from the models here presented could help to improve regional to global estimations of carbon fluxes and refine national GHG inventories for our country and perhaps carbon sequestration capacity assessments for similar regions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.