Intense Pasture Management in Brazil in an Integrated Crop-Livestock System Simulated by the DayCent Model

Process-based models (PBM) are important tools for understanding the benefits of Integrated Crop-Livestock Systems (ICLS), such as increasing land productivity and improving environmental conditions. PBM can provide insights into the contribution of agricultural production to climate change and help identify potential greenhouse gas (GHG) mitigation and carbon sequestration options. Rehabilitation of degraded lands is a key strategy for achieving food security goals and can reduce the need for new agricultural land. This study focused on the calibration and validation of the DayCent PBM for a typical ICLS adopted in Brazil from 2018 to 2020. We also present the DayCent parametrization for two forage species (ruzigrass and millet) grown simultaneously, bringing some innovation in the modeling challenges. We used aboveground biomass to calibrate the model, randomly selecting data from 70% of the paddocks in the study area. The calibration obtained a coefficient of determination (R2) of 0.69 and a relative RMSE of 37.0%. During the validation, we used other variables (CO2 flux, grain biomass, and soil water content) measured in the ICLS and performed a double validation for plant growth to evaluate the robustness of the model in terms of generalization. R2 validations ranged from 0.61 to 0.73, and relative RMSE from 11.3 to 48.3%. Despite the complexity and diversity of ICLS results show that DayCent can be used to model ICLS, which is an important step for future regional analyses and large-scale evaluations of the impacts


Introduction
Integrated Crop-Livestock Systems (ICLS) are defined by the diversification, rotation, consortium, and succession of agriculture and livestock inside the same area, allowing benefits between these two activities and the economic use of soil over the year [1]. The adoption of ICLS comprises 10-15 million ha (5% of total) of the area under agriculture and livestock in Brazil [2], which represents a significant opportunity to improve the sustainability of agricultural systems [3], especially in degraded areas and over sandy soils.
With the integration of soybean and mixed-pasture, ICLS impacts plant growth as the nitrogen (N) supplied by legumes (for instance, soybean) is absorbed more efficiently

Field Data
The data were obtained from a commercial area of 200 ha divided into four fields of 50 ha, located at the municipality of Caiuá, in the western region of São Paulo state, Brazil (21°38′ S and 51°54′ W; Figure 1). The climate is classified as tropical with a dry winter season (Aw), according to Köppen's classification [28]. The accumulated rainfall, mean minimum and maximum daily air temperatures during the monitored period were, respectively, 2021.8 mm, 20.4, and 32.4 °C, from October 2018 to April 2020 ( Figure 2). The meteorological data were obtained by prediction of worldwide energy resources (hereafter called NASA/POWER; https://power.larc.nasa.gov/data-access-viewer/, accessed on 2 July 2021), which has shown to be a good source for model simulations under Brazilian conditions [29][30][31].
The soil of the area is formed from rocks of Bauru Group sandstone of the Cretaceous period [32] and has the dominance of kaolinite in the clay fraction [33]. The soil of the region is classified as Oxisol with a sandy loam texture [34]. The slope between 3 to 7% is classified as "gently sloping" according to USDA classification [35]. Elevation range is between 310 m and 370 m above mean sea level.

Field Data
The data were obtained from a commercial area of 200 ha divided into four fields of 50 ha, located at the municipality of Caiuá, in the western region of São Paulo state, Brazil (21°38′ S and 51°54′ W; Figure 1). The climate is classified as tropical with a dry winter season (Aw), according to Köppen's classification [28]. The accumulated rainfall, mean minimum and maximum daily air temperatures during the monitored period were, respectively, 2021.8 mm, 20.4, and 32.4 °C, from October 2018 to April 2020 ( Figure 2). The meteorological data were obtained by prediction of worldwide energy resources (hereafter called NASA/POWER; https://power.larc.nasa.gov/data-access-viewer/, accessed on 2 July 2021), which has shown to be a good source for model simulations under Brazilian conditions [29][30][31].
The soil of the area is formed from rocks of Bauru Group sandstone of the Cretaceous period [32] and has the dominance of kaolinite in the clay fraction [33]. The soil of the region is classified as Oxisol with a sandy loam texture [34]. The slope between 3 to 7% is classified as "gently sloping" according to USDA classification [35]. Elevation range is between 310 m and 370 m above mean sea level.   The soil of the area is formed from rocks of Bauru Group sandstone of the Cretaceous period [32] and has the dominance of kaolinite in the clay fraction [33]. The soil of the region is classified as Oxisol with a sandy loam texture [34]. The slope between 3 to 7% is classified as "gently sloping" according to USDA classification [35]. Elevation range is between 310 m and 370 m above mean sea level.
The pre-experiment started in the area defined by the Atlantic Forest biome (the vegetation is constituted of a dense forest of medium and large trees with a large amount of plant residue on the soil surface), which had extensive pasture without proper management (estimated beginning of extensive pasture in 1984 by farm data), either with fertilizers or conventional management practices, such as recovery/renewal pastures or animal management. In August 2007, 0.5 Mg ha −1 of dolomitic limestone was applied, followed by conventional practices of ploughing and harrowing ( Figure 3). Furthermore, 400 and 200 kg ha −1 were applied of reactive rock phosphate (35% P 2 O 5 total and 12.5% P 2 O 5 in citric acid solution) and soluble fertilizer NPK (4-21-15), respectively.
Single forage (Urochloa brizantha cv. Marandu) was sown in October 2007. Fertilization was carried out annually, consisting of the application of 300 kg ha −1 of the NPK formulation 20-10-10, and this fertilization was occasionally applied only to maintain the grazing.
Heifers aged 2 to 3 years (450 kg animal −1 ) were kept in extensive pasture in the experimental area throughout the period (2007-2018), with an average stocking rate around 1.6 to 1.8 animal unit (AU) ha −1 . All these land-use data were used to schedule the pre-experiment for the study area.
This pre-experiment simulation covers 1984 to 2017 with extensive pasture and management, shown in Figure 3. Soil cover reported by the farm was made of crops in previous periods through local information and confirmation via the SATVEG platform [36]. We modified the fertilization, crop, and grazing defaults for the pre-experiment to represent the land-use change and management. At the end of the pre-experiment simulation, total SOC (in the soil layer from 0-0.2 m) simulated was compared with those observed in the field from the initial soil analyses made before implementing the ICLS. Soil samples were taken in October 2018 for the chemical and physical characteriza tion of each paddock (Table 1 and Figure 1). Sampling field collections were performed t obtain soil physical and chemical characteristics, aboveground biomass and grain at har vest, and volumetric soil water content (VWC) ( Table A2). The initial measurements o the experimental area were performed before the ICLS implementation by stratified sys tematic unaligned sampling within each paddock and with a variation of 0-1.0 m for so and surface for plant data at the end of the extensive pasture period (Tables 1 and A2).  During the ICLS-experiment period from 2018 to 2020, the ICLS consisted of two phases: summer season (soybean) and winter season (mixed-pasture) ( Figure 3, Table A1). Preparation for the ICLS began in mid-August 2018, with the application of 1.5 Mg ha −1 of gypsum plus limestone. Then, approximately 5 Mg ha −1 of organic compost was applied to the terrace channel of the area in early September. The herbicide application was carried out to enter the new crop, following common practices in the other farms of the group.
The area was divided into four fields of about 50 ha each (Figure 1), in which soybean was sown under no-tillage. Different varieties were sown in each field (Table A1), with dates ranging from 18 to 22 November 2018, for the first soybean cycle. Monoammonium phosphate (MAP) and potassium chloride (KCl) were applied in 250 and 150 kg ha −1 in the soybean planting furrow under no-tillage. The soybean harvest took place from 30 March to 5 April 2019.
For the implantation of the pasture, the fields were divided into paddocks each to obtain a rotating grazing system, totaling 13 paddocks ( Figure 1). Winter season cultivation and pasture intercropping began simultaneously after the soybean harvest (2 to 6 April 2019). The pasture was composed of a consortium of millet (Pennisetum glaucum) and ruzigrass (Urochloa ruziziensis) with a seeding rate in the proportion of 15 kg ha −1 and 5 kg ha −1 of seeds, respectively. The rotational system consisted of two or three grazing events in each paddock depending on each field of pasture development, the first of which was of greater intensity (paddock average ≈ 156 cattle paddock −1 ), the second of lesser intensity (paddock average ≈ 84 cattle paddock −1 ), and the third of medium intensity (paddock average ≈ 109 cattle paddock −1 ), which was not performed for field 4 (Table A1). In November, the herbicide was applied in the total area, and the preparation for planting the second soybean in the area was carried out.
The second soybean cycle started between 11 and 15 December 2019, with the harvest being carried out from 4 to 9 April 2020. The same soybean varieties from the last cultivation were sown, but varied among the fields (Table A1). Like the first soybean cycle, the second was fertilized during the sowing with the previously mentioned fertilizer method.
Soil samples were taken in October 2018 for the chemical and physical characterization of each paddock (Table 1 and Figure 1). Sampling field collections were performed to obtain soil physical and chemical characteristics, aboveground biomass and grain at harvest, and volumetric soil water content (VWC) ( Table A2). The initial measurements of the experimental area were performed before the ICLS implementation by stratified systematic unaligned sampling within each paddock and with a variation of 0-1.0 m for soil and surface for plant data at the end of the extensive pasture period (Tables 1 and A2). Subsequently, physical (clay, sand, and bulk density) and chemical (total C and pH) analyses were performed following protocols from the Manual of soil analysis methods [38]. The bulk density (volumetric ring method) was determined by differentiating the sampling by depth with three repetitions for 0-0.20 m spaced by 0.05 m depth in each sample, one of 0.25-0.30 m, one sample of 0.35-0.40 m, a 0.65-0.70 m sample, and a 0.95-1.00 m sample. For the total C content, samples were collected at the same points and depth as the bulk density of the soil. These were air-dried, sieved with a 2.0 mm diameter sieve, then manually ground in a mortar and pistil, and a new 0.177 mm sieving [38]. Then, the samples were sent to the laboratory, which performed the analysis by dry oxidation by the elemental analyzer CHNS-O. This method is based on the oxidation of samples at high temperatures (approximately 1000 • C). The samples were placed in a tin capsule, which does not contain C. After total combustion, the gases containing each element were separated and the concentrations measured by different infrared detectors.
Regular collections of the aboveground biomass were performed after the ICLS phase started, varying from two to three repetitions for each paddock. For soybeans, one linear meter was considered, according to the sowing row, and for mixed-pasture, an area of 1 m 2 , in both cases cutting the biomass close to the ground (≈0.01 m) ( Figure 1, Table A2). The aboveground biomass material was dried in an oven with forced air circulation (65-70 • C) until reaching constant mass. Biomass C content was calculated considering the factor of 0.475 [39], with all the results expressed in g C m −2 .
Furthermore, VWC was determined by gravimetric water content measurements multiplied by soil density (0-0.05 m depth, drying in forced air circulation oven) using the same sampling scheme in which biomass measurements were taken. CO 2 flux was determined using a gas flow quantification chamber (LI 8100A, LI-COR Environmental, Lincoln, NE, USA), with nine samples per paddock on three different dates different from the aboveground biomass collection dates for the mixed-pasture, four for the second soybean, and two for the transition between seasons (Table A2). The measurement was carried out with a duration of 1.5 min for each collection point and carried out around 7 a.m. to avoid soil temperature instability. For the standardization of the CO 2 flux values, it was assumed that for the DayCent output, the CO 2 flux is constant throughout the day, that is ((CO 2 ) total = 24 (CO 2 ) hours) to compare with that measured and have a simplification of the soil respiration processes.

SOC Reference Data from the Literature
To start the equilibrium modeling, it was necessary to use literature data on SOC from native vegetation close to the study area in Caiuá, São Paulo (SP). Published data from Rigolin et al. [40], compared the C stocks (0-0.2 m) and the different managements in an area in Presidente Prudente, SP (22 • 07 S and 51 • 27 W), and used data from native vegetation, without human intervention for 50 years; C stocks were compared with the one simulated by DayCent because it was from the same biome and with similar climatic and soil conditions ( Figure A1 detail d). For the extensive pasture, we used data collected on 28 October 2018, which collected soil bulk density and total soil C content up to 1 m depth to estimate C stocks and aboveground biomass for all paddock points in the study area in Caiuá, SP.
To increase confidence that the model is applicable for different areas of ICLS after the calibration, we employed data of the aboveground biomass from the literature called other regions (OR) using [41,42] datasets, both with an ICLS with mixed-pastures and soybean. The same calibrated parameters were used without any change, changing only the input data and the management adopted in the model. Souza's [41] study was carried out in Rondonópolis, Mato Grosso state (MT), Brazil (16 • 27 S and 54 • 34 W), located in the Cerrado biome.
The simulations for equilibrium were set to represent such a biome. The experiment was carried out in a succession of the soybean crop and had a mixture of ruzigrass with millet. The mixed-pasture aboveground biomass was collected during the experiment, and the soybean yields. The experimental design was a randomized block design, with seven treatments with three replications per treatment (can be found on the GitHub platform (https://github.com/yanefsilva/DayCent_ICLS.git, 21 February 2022)).
For the experiments published by Machado [42], ICLSs were evaluated in two locations, São Gabriel do Oeste (19 • 24 S and 54 • 34 W) and Dourados (22 • 14 S and 54 • 49 W), both in Mato Grosso do Sul state (MS), Brazil. Both sites are in the Cerrado biome (can be found on the GitHub platform (https://github.com/yanefsilva/DayCent_ICLS.git, 21 February 2022)). Thus, the same equilibrium file built for [41] was used, changing only the soil and climate inputs. In this study, the mixed-pasture in rotation with the soybean crop were also grown in the ICLS. The study used an experimental design which was a 9 × 2 × 2 factorial, in which the A factor is forages and their mixtures, factor B is the two locations (São Gabriel do Oeste and Dourados), and factor C is the two times of evaluation (2007 and 2008), with four replications; the aboveground biomass and soybean grain yield data were measured and used in the modeling study.
The soil and management data presented in the above-mentioned studies were adopted, and the climate data were acquired from NASA/POWER ( Figure A1 details a, b and c).

Model Description and Simulation Settings
The software version of the DayCent model (DD15centEVI.exe) (Colorado State University, Natural Resources Ecology Laboratory) [14,43] was used to perform the simulations. Briefly, the DayCent model includes sub-models to estimate plant growth, soil temperature, VWC, and the biogeochemical cycle of nutrients, including C and N fluxes between soils and plants [23]. DayCent is the daily version of the Century model [44,45], differing from its predecessor. It provides results on a daily time scale and dynamically represents the emission processes of GHGs, as the release of CO 2 and N 2 O in the soil [14,23].
The plant growth sub-model simulates the potential biomass production and yield using solar radiation and plant photosynthetic efficiency determined by species, variety, or clone. Potential yield is then corrected by multipliers related to air temperature, soil water, nutrient availability (N, P, and S), shading, phenological stage, and atmospheric CO 2 concentration [23]. The achievable yield is less than or equal to the potential yield depending on the factors mentioned above and their interaction with the management practices applied in the area, which include different crops (no-till, conventional tillage, and mowing), mineral and organic fertilization, irrigation, pasture, and fire [23]. N 2 O and CH 4 were not evaluated for the study because of the lack of data to generate results.
In each paddock, simulations for ICLS were carried out at two to three points, and the points were averaged to represent each paddock studied, simulating the impacts of the ICLS implementation in the initial integration period between mixed-pasture and soybean ( Figure 4). The grazing intensity was considered to impact the aboveground biomass production, evaluated in the model. All files used in equilibrium, pre-experiment and experiment can be found on the GitHub platform (https://github.com/yanefsilva/ DayCent_ICLS.git, 21 February 2022).  To simulate the status of the C pools for the ICLS fully detailed in Section 2.2, an equilibrium simulation with the historical land cover and land-use changes is necessary. Usually, it is considered a first run with the original natural vegetation (Atlantic Forest in the present case) for ~6000 years, time to reach the C equilibrium between the different C pools solved by the model. During this equilibrium simulation, the C is gained through litter decomposition and senesced roots and lost by heterotrophic respiration. The C fluxes between the active, slow, and passive SOC pools tend to reach equilibrium [46]. To generate the model file corresponding to the Atlantic Forest equilibrium, we chose to modify the default values, such as the PRDX (2) parameter (Table A3), inside the tree.100 file, which indicates the potential C production by the aerial part of the vegetation, to approximate the average values corresponding to the literature's reported values for the Atlantic Forest ( Figure 4). To simulate the status of the C pools for the ICLS fully detailed in Section 2.2, an equilibrium simulation with the historical land cover and land-use changes is necessary. Usually, it is considered a first run with the original natural vegetation (Atlantic Forest in the present case) for~6000 years, time to reach the C equilibrium between the different C pools solved by the model. During this equilibrium simulation, the C is gained through litter decomposition and senesced roots and lost by heterotrophic respiration. The C fluxes between the active, slow, and passive SOC pools tend to reach equilibrium [46]. To generate the model file corresponding to the Atlantic Forest equilibrium, we chose to modify the default values, such as the PRDX (2) parameter (Table A3), inside the tree.100 file, which indicates the potential C production by the aerial part of the vegetation, to approximate the average values corresponding to the literature's reported values for the Atlantic Forest ( Figure 4).
Afterwards, we ran the model considering this equilibrium, in the total SOC from the Atlantic Forest simulation, as the initial condition for an extended simulation covering the period from the deforestation to the present, considering the land cover changes and typical crop management practices in the period, i.e., the history of land cover before the beginning of the ICLS.

Calibration and Validation Procedures
The main input variables of DayCent are minimum and maximum air temperature ( • C), precipitation (cm), soil texture, pH, bulk density (g cm −3 ), and management practices. These variables were collected from fields, and the missing information was estimated from the use of DayCent in the previous literature calibrations, considering the soil physical, water flow and soil temperature, crop growth, and soil organic matter sub-models. In addition, historical information on soil (pedological map from São Paulo state), climate (NASA/POWER), and general information on the biome and the history of land use (<mapbiomas.org/>) were collected for use as input data in DayCent. Soil profile water retention characteristics were derived from field texture samples and the pedotransfer functions of Saxton (<https://hrsl.ba.ars.usda.gov/soilwater/Index.htm>, 21 February 2022 [47]).
As the dataset collected was only in one farm, we separated the dataset into two parts, for calibration and validation, following Garrison et al. [48]. For the calibration, we used 70% of the data sets collected in the field (9 paddocks from the study area, randomly chosen), and 30% of the data were used for model validation (4 paddocks, one from each field) and also general model evaluation.
Calibrations were made through systematic variations in key model parameters until the best fit (details of statistical metrics provided later) was obtained between the Day-Cent output and the corresponding observed data, within the acceptable range for each parameter [11,23]. We performed the model calibration by adjusting parameters (details in the Results section, Table 2), as proposed by Del Grosso et al. [22], which considers the submodel of plant growth, with aboveground biomass C (Aglivc) as the main output parameter to be considered in the calibration, as it is a complex integrated system, with cattle in and out, thus making it impossible to collect different data. As we have a satisfactory amount of Aglivc data, we focus on Aglivc as the main parameter for calibration. As the pasture was submitted to a variable cattle stocking rate, rotational management system, we also changed the grazing file (FLGREM and FDGREM, Table A3) to adjust the number of animals in the paddocks.
The parameters calibrated for the harvest were the fraction of aboveground live which will not be affected by harvest operations (AGLREM) and the fraction of belowground live which will not be affected by harvest operations (BGLREM) (Table A3); as the use of clearing was not standardized to reduce the height of the canopy and thus achieve leveling of the pasture, we used pasture cut-off values similar to clearing done for new management in the area [49,50] (Tables 2 and A1). For fertilization, the calibrated took into account the MAP and the fertilizers used during the pre-experiment, and the amount of N left by the fertilizer [51][52][53] (Tables 2 and A1). For the crop model, we added the default for millet with ruzigrass to represent different growth rates for the mixed-pastures according to the literature [54][55][56][57]. As in the experiment ICLS, for each field, there were stocking rates of different animals between the paddocks, we made the necessary changes for grazing, considering the low, medium, and intermediate grazing (Tables 2 and A1). The calibrated parameters considered available information from the literature on ICLS and data collected in the field to estimate the parameters as best as possible.  We performed a double validation for the ICLSs described in Section 2.1. First, we used 30% of the study area paddocks (4 paddocks) for analysis, with the same input data for each paddock variable and keeping the same climate file as the calibration simulations, as well as changing the soils and soil management practices of each of the four paddocks; all calibrated parameters were kept unchanged for the validation runs. Grain (for soybean) and aboveground biomass for the paddocks not used in the calibration were used for a set of validation. At this stage, we also compared model simulations with measured CO 2 fluxes and VWC for all paddocks to test the model capability to predict other agroecosystem variables besides the plant ones.
Additionally, we tested the model with biometric data from three different sites described in Section 2.1, from now on called 'OR', to evaluate the possibility of applying the model in different environments in which farmers adopt similar ICLS. In this second validation, the same parameters were used without any change, changing only the input data in the model, which covers data from two studies published by Souza [41] and Machado [42]. Details about the equilibrium simulations for these other ICLS can be found on the GitHub platform https://github.com/yanefsilva/DayCent_ICLS.git, 21 February 2022.
We used the Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), Coefficient of determination (R 2 ), and Nash-Sutcliffe model Efficiency (NSE) as statistical metrics for evaluating the model calibration, as we did not have repetitions of observations of the data collected in each paddock. It was found that the RMSE could not be used for intercomparisons for many state variables with different units [58], such as aboveground biomass, water content, and CO 2 flux. In these cases, the relative RMSE (in percentage) was used as a relative measure for intercomparisons of different variables or different models [59].
For the optimal calibration adjustment, we considered the lowest mean absolute error (MAE) and NSE values close to 1. NSE equal to 1 represents a perfect agreement, and NSE above 0 indicates that the simulated values describe a trend for the observed data better than the average of the observations as the predictor. The performance of the model was based on NSE [60], for which a performance rating of 'very good' is equivalent to NSE 0.75-1.00, 'good' NSE 0.65-0.75, 'satisfactory' NSE 0.50-0.65, and the NSE is 'unsatisfactory' if <0.50. Although these classifications were proposed by Moriasi et al. [60] for hydrological studies, we assumed they are also applicable to agricultural or forestry ones. Statistical analyses were carried out in the free R environment [61], using caret package [62].

Equilibrium Modeling
The total simulated SOM was around 39.0 Mg C ha −1 , where the major fraction was stored in Slow SOM, followed by the Passive SOM and a relatively small fraction in the Active SOM ( Figure 5). The model suggested that SOM may take millennials to reach an equilibrium between the C addition and loss for the Atlantic Forest biome. The C stocks found for a native forest near the farm was reported as 36.4 Mg C ha −1 , a difference of 6.7% compared to that simulated by DayCent ( Figure 6).
For the simulations with the extensive pasture, in all paddocks in the study area in Caiuá, SP, the model calibration achieved good results because the values of total SOC measured in the literature compared their averages with a difference of 17.3% by DayCent ( Figure 6).

Model Calibration
We considered the PRDX parameter in the mixed-pasture and soybean crop to achieve a refined calibration based on the observed values ( Figure 7 and Table 2). For the first soybean cycle, the model responded well, while for the second soybean cycle and mixedpasture, there were some paddocks (namely F3P4, F4P1, and F4P3) in which the model had a reasonable to poor response (Figure 7). equilibrium between the C addition and loss for the Atlantic Forest biome. The C stocks found for a native forest near the farm was reported as 36.4 Mg C ha −1 , a difference of 6.7% compared to that simulated by DayCent ( Figure 6).
For the simulations with the extensive pasture, in all paddocks in the study area in Caiuá, SP, the model calibration achieved good results because the values of total SOC measured in the literature compared their averages with a difference of 17.3% by DayCent ( Figure 6).

Model Calibration
We considered the PRDX parameter in the mixed-pasture and soybean crop to achieve a refined calibration based on the observed values ( Figure 7 and Table 2). For the first soybean cycle, the model responded well, while for the second soybean cycle and mixed-pasture, there were some paddocks (namely F3P4, F4P1, and F4P3) in which the model had a reasonable to poor response (Figure 7).
The parameters modified during equilibrium calibration, pre-experiment, and ICLS were considered using the range suggested by the DayCent manual (Tables 2 and A3). The values of the parameters were estimated according to the literature; in the case of the C/N ratio, we used the changes in the PRAMX and PRAMN parameters (Table 2), following [63]. Plants with high N levels in their biomass, such as legumes, provide residues with a low C/N ratio, which means faster decomposition and, consequently, a higher N mineralization rate than other vegetables [63].
In the model calibration, the R 2 for Aglivc was 0.69, which was statistically strong considering the observed data, with relative RMSE of 37.0%, MAE of 24.2%, and an NSE of 0.66 classified as satisfactory according to Moriasi [60] (Figure 8).

Model Calibration
We considered the PRDX parameter in the mixed-pasture and soybean crop to achieve a refined calibration based on the observed values ( Figure 7 and Table 2). For the first soybean cycle, the model responded well, while for the second soybean cycle and mixed-pasture, there were some paddocks (namely F3P4, F4P1, and F4P3) in which the model had a reasonable to poor response (Figure 7).
The parameters modified during equilibrium calibration, pre-experiment, and ICLS were considered using the range suggested by the DayCent manual (Tables 2 and A3). The values of the parameters were estimated according to the literature; in the case of the C/N ratio, we used the changes in the PRAMX and PRAMN parameters (Table 2), following [63]. Plants with high N levels in their biomass, such as legumes, provide residues with a low C/N ratio, which means faster decomposition and, consequently, a higher N mineralization rate than other vegetables [63].
In the model calibration, the R 2 for Aglivc was 0.69, which was statistically strong considering the observed data, with relative RMSE of 37.0%, MAE of 24.2%, and an NSE of 0.66 classified as satisfactory according to Moriasi [60] (Figure 8).  The parameters modified during equilibrium calibration, pre-experiment, and ICLS were considered using the range suggested by the DayCent manual (Tables 2 and A3). The values of the parameters were estimated according to the literature; in the case of the C/N ratio, we used the changes in the PRAMX and PRAMN parameters (Table 2), following [63]. Plants with high N levels in their biomass, such as legumes, provide residues with a low C/N ratio, which means faster decomposition and, consequently, a higher N mineralization rate than other vegetables [63].
In the model calibration, the R 2 for Aglivc was 0.69, which was statistically strong considering the observed data, with relative RMSE of 37.0%, MAE of 24.2%, and an NSE of 0.66 classified as satisfactory according to Moriasi [60] (Figure 8).

Model Validation
The model presented a consistent growth pattern along the seasons for the extensive pasture, before installing the ICLS, to the second soybean cycle and when the mixed-pasture was established (Figure 9). The mixed-pasture with values below 150 g C m −2 , the vegetative peak of the first soybean cycle values above 300 g C m −2 , and the second soybean, perhaps due to climatic factors, Aglivc was lower than the values of the first planting in the four paddocks. The model captured the monthly variations in VWC and the variations between paddocks related to the soil sandy fraction variability and the crop management (cattle management between paddocks) after calibration, shown in Figure 10. Except for the F4P2 paddock (Figure 10d), all the others had similar behavior between simulated and observed during the study period yielding an R 2 of 0.73, relative RMSE of 37.0%, MAE of 27.1%, and NSE of 0.51, indicating a good agreement (Figure 10b).
Regarding soil CO2 flux, the model overestimated three of the four paddocks used to validate such an agroecosystem variable. However, their average data were within the range of observed data during the ICLS (Figure 11). The lower NSE concerning the other analyzed parameters was 0.44, and relative RMSE was considered high at 48.3%, MAE of 31.9% and R 2 of 0.61. By Tukey's test, two of the four paddocks were considered to have equal means between simulated and observed ( Figure 11).

Model Validation
The model presented a consistent growth pattern along the seasons for the extensive pasture, before installing the ICLS, to the second soybean cycle and when the mixedpasture was established (Figure 9). The mixed-pasture with values below 150 g C m −2 , the vegetative peak of the first soybean cycle values above 300 g C m −2 , and the second soybean, perhaps due to climatic factors, Aglivc was lower than the values of the first planting in the four paddocks.

Model Validation
The model presented a consistent growth pattern along the seasons for the extensive pasture, before installing the ICLS, to the second soybean cycle and when the mixed-pasture was established (Figure 9). The mixed-pasture with values below 150 g C m −2 , the vegetative peak of the first soybean cycle values above 300 g C m −2 , and the second soybean, perhaps due to climatic factors, Aglivc was lower than the values of the first planting in the four paddocks. The model captured the monthly variations in VWC and the variations between paddocks related to the soil sandy fraction variability and the crop management (cattle management between paddocks) after calibration, shown in Figure 10. Except for the F4P2 paddock (Figure 10d), all the others had similar behavior between simulated and observed during the study period yielding an R 2 of 0.73, relative RMSE of 37.0%, MAE of 27.1%, and NSE of 0.51, indicating a good agreement (Figure 10b).
Regarding soil CO2 flux, the model overestimated three of the four paddocks used to validate such an agroecosystem variable. However, their average data were within the range of observed data during the ICLS (Figure 11). The lower NSE concerning the other analyzed parameters was 0.44, and relative RMSE was considered high at 48.3%, MAE of 31.9% and R 2 of 0.61. By Tukey's test, two of the four paddocks were considered to have equal means between simulated and observed ( Figure 11). The model captured the monthly variations in VWC and the variations between paddocks related to the soil sandy fraction variability and the crop management (cattle management between paddocks) after calibration, shown in Figure 10. Except for the F4P2 paddock (Figure 10d), all the others had similar behavior between simulated and observed during the study period yielding an R 2 of 0.73, relative RMSE of 37.0%, MAE of 27.1%, and NSE of 0.51, indicating a good agreement (Figure 10b).  In Figure 12, we brought together all Aglivc values from the validated paddocks and the two other regions (three ICLSs) to analyze how well the model responded to the observed values. For this generalized validation, model performance for Aglivc yielded R 2 of 0.67, RMSE of 41.6%, and NSE of 0.63, indicating a good model fit.
We are considering a system with an output of two soybean cycles, one mixed-pasture cycle and the other data from farms that also used ICLS as a basis for mixed-pasture ( Figure 12). The relative RMSE at 41.6%, MAE with 28.5%, and NSE with 0.63, we can consider as indicative of a low bias in the model, which was reinforced by visual analyzing of the convergence between the confidence interval line (gold) (95%) and the standard deviation (gray) with the regression line (observed vs. simulated) (Figure 12), with the majority of the simulated points falling within this confidence interval. Regarding soil CO 2 flux, the model overestimated three of the four paddocks used to validate such an agroecosystem variable. However, their average data were within the range of observed data during the ICLS (Figure 11). The lower NSE concerning the other analyzed parameters was 0.44, and relative RMSE was considered high at 48.3%, MAE of 31.9% and R 2 of 0.61. By Tukey's test, two of the four paddocks were considered to have equal means between simulated and observed ( Figure 11).  In Figure 12, we brought together all Aglivc values from the validated paddocks and the two other regions (three ICLSs) to analyze how well the model responded to the observed values. For this generalized validation, model performance for Aglivc yielded R 2 of 0.67, RMSE of 41.6%, and NSE of 0.63, indicating a good model fit.
We are considering a system with an output of two soybean cycles, one mixed-pasture cycle and the other data from farms that also used ICLS as a basis for mixed-pasture ( Figure 12). The relative RMSE at 41.6%, MAE with 28.5%, and NSE with 0.63, we can consider as indicative of a low bias in the model, which was reinforced by visual analyzing of the convergence between the confidence interval line (gold) (95%) and the standard deviation (gray) with the regression line (observed vs. simulated) (Figure 12), with the majority of the simulated points falling within this confidence interval. In Figure 12, we brought together all Aglivc values from the validated paddocks and the two other regions (three ICLSs) to analyze how well the model responded to the observed values. For this generalized validation, model performance for Aglivc yielded R 2 of 0.67, RMSE of 41.6%, and NSE of 0.63, indicating a good model fit.
We are considering a system with an output of two soybean cycles, one mixed-pasture cycle and the other data from farms that also used ICLS as a basis for mixed-pasture ( Figure 12). The relative RMSE at 41.6%, MAE with 28.5%, and NSE with 0.63, we can consider as indicative of a low bias in the model, which was reinforced by visual analyzing of the convergence between the confidence interval line (gold) (95%) and the standard deviation (gray) with the regression line (observed vs. simulated) (Figure 12), with the majority of the simulated points falling within this confidence interval. In the validation set, simulations of grain biomass (cgrain) presented similar average and quantiles' values ( Figure 13). The model was able to capture the variability among paddocks. Moreover, for the OR that the modeling was performed to prove the robustness of the model in ICLS, the model slightly overestimated the cgrain but represented well the variability. Tukey's test indicated equal means for the four validated paddocks and OR, showing no significant difference between simulated and observed. Figure 13. Boxplots of grain biomass for the four validation fields and for the OR (other regions) that had the same ICLSs with the observed data versus the predicted data. Tukey's test (p < 0.05) was performed to compare the average between the observed and the simulated. The letters indicate the similarity between the validation sites and the simulated by DayCent.

Equilibrium
The replacement of natural ecosystems by agroecosystems with crops usually provides the decline in soil C content due to the reduction in the microbial C and the increase of SOM decomposition [64]. With the removal of the native vegetation and the beginning of the extensive pasture, the C stocks dropped in the soil (Figure 6). This may have been caused by the soil preparation (disturbing), breaking the aggregates offering protection to SOM and stimulating the biological activity [65]. In the validation set, simulations of grain biomass (cgrain) presented similar average and quantiles' values ( Figure 13). The model was able to capture the variability among paddocks. Moreover, for the OR that the modeling was performed to prove the robustness of the model in ICLS, the model slightly overestimated the cgrain but represented well the variability. Tukey's test indicated equal means for the four validated paddocks and OR, showing no significant difference between simulated and observed. In the validation set, simulations of grain biomass (cgrain) presented similar average and quantiles' values ( Figure 13). The model was able to capture the variability among paddocks. Moreover, for the OR that the modeling was performed to prove the robustness of the model in ICLS, the model slightly overestimated the cgrain but represented well the variability. Tukey's test indicated equal means for the four validated paddocks and OR, showing no significant difference between simulated and observed. Figure 13. Boxplots of grain biomass for the four validation fields and for the OR (other regions) that had the same ICLSs with the observed data versus the predicted data. Tukey's test (p < 0.05) was performed to compare the average between the observed and the simulated. The letters indicate the similarity between the validation sites and the simulated by DayCent.

Equilibrium
The replacement of natural ecosystems by agroecosystems with crops usually provides the decline in soil C content due to the reduction in the microbial C and the increase of SOM decomposition [64]. With the removal of the native vegetation and the beginning of the extensive pasture, the C stocks dropped in the soil (Figure 6). This may have been caused by the soil preparation (disturbing), breaking the aggregates offering protection to Figure 13. Boxplots of grain biomass for the four validation fields and for the OR (other regions) that had the same ICLSs with the observed data versus the predicted data. Tukey's test (p < 0.05) was performed to compare the average between the observed and the simulated. The letters indicate the similarity between the validation sites and the simulated by DayCent.

Equilibrium
The replacement of natural ecosystems by agroecosystems with crops usually provides the decline in soil C content due to the reduction in the microbial C and the increase of SOM decomposition [64]. With the removal of the native vegetation and the beginning of the extensive pasture, the C stocks dropped in the soil ( Figure 6). This may have been caused by the soil preparation (disturbing), breaking the aggregates offering protection to SOM and stimulating the biological activity [65].
In tropical ecosystems, the losses of C observed after the land clearing and cultivation of soils are more accelerated than in temperate regions [66]. Thus, the decrease of nutrients in an extensive pasture in native vegetation is due to soil revolving, even if minimal, since there is a negative effect on SOM [67][68][69]. Climatic conditions, especially temperature, also affect the decomposition rate of plant residues and the stabilization of SOM [44]. The occurrence of higher temperatures tends to accelerate the decomposition rate of SOM [70]. One alternative to increase the total C stock in the agriculture systems is by applying sustainable management practices and adopting crop systems that increase primary production and reduce C removal and soil degradation. Therefore, improved management practices affect the factors that regulate the synthesis and decomposition of SOM, contributing to the sustainable productivity of the ecosystem [71].

Calibrated Model
The DayCent model adequately simulated aboveground biomass (Figure 7), mainly for the beginning of the soybean crop and the mixed-pasture growth phases, periods that had no direct grazing effect on the simulations; even though it was an area with intensive pasture management, there were periods just after the mixed-pasture implementation where the animals were not present in some paddocks. All nine paddocks had the same altered parameters, and the behavior throughout the experiment was similar. Thus, minor model failures in predicting aboveground biomass did not influence each paddock's overall behavior. Due to the study area being a commercial farm, minor problems with pests, diseases, and weeds may also have influenced the biomass production to some extent, all of which are not currently captured by the DayCent model unless we force the model to reproduce the reductions of potential plant growth.
Del Grosso et al. [22] reported that most of the errors in DayCent model outputs were associated with imperfections in model algorithms and parameters instead of uncertainty in model drivers. Therefore, efforts to improve the model should compare model outputs using numerous observations for various C and N components from field experiments to identify weaknesses and rectify model shortcomings [72].
There was no under-or overestimation tendency of the aboveground biomass simulations ( Figure 8). The model explained 0.69 of the mixed-pasture, soybean, and extensive pasture measured data (Figure 8), which many authors found a 0.25-0.69 variation over aboveground biomass during the calibration stage [73].
In ICLS, there is a differentiated contribution of plant residues regarding conventional grain production systems, both in the surface and soil subsurface [74]. In intensive ICLS, the often succession of crop and pasture planting contributes to the root growth of both plants and potentially increases the SOM in the deeper soil layers [75]. SOM can be investigated by detailing its behavior over time in the future. Not only the inclusion of crops in pasture areas can contribute to improving the system, but also planting pastures in crop areas can increase crop yield by improving edaphic properties, the presence of straw and pasture roots, increasing C levels and significantly ameliorating the conditions of aeration and water infiltration capacity in the soil [76]. In future, we need to investigate the GHG mitigation which ICLS can potentially offer, which farmers can use to understand their carbon footprint, encouraging management strategies to improve agricultural sustainability [77].

Validated Model
In the validation stage, the DayCent model performed well in representing the ICLS, with the model results reflecting the dynamic phases of such a system (e.g., soybean replaced by mixed-pasture with a rotational and intensive grazing system) and management operations (e.g., sowing, harvesting, mowing, fertilization, grazing events) (Figures 7 and 9). Even though three different soybean varieties were sowing in the area in both seasons, the model presented consistent aboveground biomass and yield simulation results. Perhaps the difference between the rainfall distribution and the difference in the average air temperature among cycles have not influenced the differentiation of aboveground biomass production.
For the four paddocks independently evaluated, the behavior of aboveground biomass shows that even with some differences, it was still possible to represent the ICLS (Figure 9). The model did not overestimate or underestimate aboveground biomass for the study period; previous studies have reported a DayCent tendency to underestimate grain yield and aboveground biomass [11,78,79]. Generally, significant errors may occur because the model does not consider several processes that depend on crop management, such as those associated with biotic stresses from weeds, diseases, or insects [80]. Additionally, improvements in the solving processes, such as replacing the generic phenology and partitioning rules by implementing traditional crop models such as the CROPGRO [81], could potentially overcome some of these challenges in predicting soybean development and growth and yields.
The DayCent model adequately simulated VWC ( Figure 10) after adjustments during the calibration simulations. We highlight that the DayCent presents consistent results for these sand soils, similar to reported performance for other soil types in Brazil [82,83]. We had a median accuracy not only for the rainy season (December and January) but also for the driest months (July and August), unlike some studies that reported DayCent difficulties in simulating VWC in very dry periods [84], based on the reference described in Section 2.4.
The VWC and its dynamics are of utmost importance since it affects plant growth, SOM mineralization, and, therefore, the C fluxes in the model. Moreover, in crops such as soybean, water stress impacts not only direct the plant growth but also triggers morphophysiological responses mainly in the reproductive phase, demonstrated by premature senescence of leaves and flowers, death of pods and, consequently, grain yield reduction [85]. Water stress effects on aboveground biomass (Figures 7 and 9) and grain biomass ( Figure 13) were noticed at some stages within the paddocks. In the DayCent modeling approach, these different species responses may be corrected by inserting different soil root distributions as model inputs, for example, adjusting the soil.in file (model input) to better represent the growth and distribution and add a water response factor in the crop.100 file to better represent crop response sensitivity to drought events. However, it is important to emphasize that these properties must be kept relatively constant for the same species. Their alteration must be based on quantifiable aspects of the production environments.
DayCent captured the highest observed values of CO 2 emission linked to management events ( Figure 11). The soil CO 2 release occurs soon after the no-till practices. A few weeks later, when more C is incorporated into the soil, the broking of structures protecting C, make it available for the soil microorganisms. Consequently, more C is decomposed and released to the atmosphere, resulting, in the short term, in a reduction of the total SOC, and, if frequent, in the SOC content in the medium and longterm [86,87]. The flow of CO 2 in the soil is very complex and incorporates both autotrophic and heterotrophic respiration and the passive diffusion of CO 2 , encompassing biological and physicochemical processes that control the CO 2 fluxes into and out the soil, and between different C compartments of the soil system [88].
DayCent has received little assessment of CO 2 predicted in tropical and ICLS conditions. As mentioned, this model was previously tested in Brazil in its capabilities to predict SOC changes [89], crop yield for common succession/rotation systems, N 2 O [15], and flows of CH 4 [79], but not yet for CO 2 fluxes, despite its importance in accurately simulating the conditions of a more intensive and complex system. DayCent overestimates the GHG over different environments [11,79,90], including our study. With an R 2 of 0.61, and relative RMSE of 48.3%, when comparing the means, we only had differences in the F2P2 paddock, which can be explained by the plant residues having a higher C/N ratio and having a slower decomposition, thus staying longer in the ground, and the CO 2 flux becomes lower than the others [91]. The results presented here show that CO 2 flux can vary with the implementation of different ICLS, and depends on the management practices. Therefore, results support the importance of measuring and validating ecosystem models, such as DayCent, to predict the SOC and CO 2 fluxes.
The validation for applications in OR showed consistent results ( Figure 13). The simulated values were similar to those observed considering the two soybean crops. With a minimum precision error about the other parameters (relative RMSE = 11.3 %, MAE = 8.0 % and NSE = 0.52), we had a good result for DayCent within the ICLS for the soybean crop considering the adjusted model also for different sites. Study [92] found R 2 = 0.32 for soybean, and in our validation, it was found at 0.61. Despite the dispersion, DayCent can estimate the productivity of various crops in different places in the world. This result suggests the model potential in estimating plant growth for monoculture systems, such as soybeans, and more complex systems, such as ICLS. The DayCent model considers abiotic influences on plant growth and development, such as solar incidence, temperature, and precipitation [86]. Thus, we expected that perfect coupling of the grain production estimated by the model would not occur ( Figure 13).
A long-term ICLS implemented on sandy soil in a farm of western São Paulo state (SP) in Brazil, using a scheme of two years of soybean in the rainy season followed by pasture (dry season) and two years of only pasture, showed soybean yields ranging from 2.9 to 4.3 Mg ha −1 , with the lowest values obtained under severe dry spells conditions [90]. These were considerably higher soybean yields than other farms in the same region in the conventional system (average of 1.8 Mg ha −1 ) [93], indicating that the increment in agricultural outputs may benefit from ICLSs.
We found similar performance for the ICLS in our area and for the other three regions ( Figure 12) as previously reported for other pasture species (R 2 = 0.67) [94]. Although the results suggest that the model can predict an ICLS, additional validations are needed due to the diversity of ICLS and its adoption over different production environments, such as cultivated species, animals raised, soil textures and classes, climates, and management practises. This is a fundamental step for any model application beyond the environmental conditions that the model was evaluated, to check how stable the calibration is and to what extent the model can be extrapolated.
The study area benefited from the adoption of the ICLS by deciding to replace a degraded pasture and adopt this system. The system creates a synergy between soybeanpasture-soybean, being beneficial for both monoculture and mixed-pasture, and helps increase soil nutrients, improving land-use diversification based on the spatial and temporal integration of crop and livestock components. The result can reduce pressure to open new areas for production, proving to be an excellent option for sustainable food production. However, gaps in knowledge about the long-term effects of ICLS can make it difficult for farmers to evaluate the benefits that ICLS may provide to their farm production and impose a barrier to payments for environmental services, such as carbon credits [95]. By properly modeling ICLS, the long-term effect can be accounted for contributing with positive credits in the C trading or even participating in traceability regarding the sustainability within the beef cattle market.

Conclusions
In this study, we calibrated and validated the DayCent model in ICLS farms in tropical Brazil, particularly over sandy soils. With the simulation of two pastures grown simultaneously, it was possible to fill in a gap that had not yet been addressed with DayCent, especially when using a system as complex and intense as the ICLS which we analyzed in our study. By using the aboveground biomass data to calibrate DayCent parameters, it was demonstrated that the model could be simplified according to the amount of data available, and it was possible to simulate plant growth, grain biomass, soil VWC, total SOC and CO 2 fluxes, considering spatio-temporal evaluation and precision and accuracy metrics. Consistent simulations of grain biomass and plant growth data for other ICLSs were also found. Our study is an important step towards understanding the benefits of ICLSs and future regional analyses of the large-scale impacts of new agricultural practices, particularly over tropical environments.

Acknowledgments:
The authors would like to thank the owner, manager, and staff of the Campina farm (CV Nelore Mocho) for their support and assistance. We are also grateful to the undergraduate and graduate students, postdoctoral researchers, and also to the technicians for helping with the field data collection and preparation. NIPE and FEAGRI/UNICAMP are acknowledged for the infrastructure and support provided for this project development. Meteorological data were obtained from the NASA Langley Research Center (LaRC) POWER Project funded through the NASA Earth Science/Applied Science Program.

Conflicts of Interest:
The authors declare no conflict of interest.   4  Table A3. The description about the modified parameters used in the DayCent model.

AGLREM
Fraction of aboveground live which will not be affected by harvest operations

BGLREM
Fraction of belowground live which will not be affected by harvest operations BIOMAX Aboveground biomass level above which the minimum and maximum C/E ratios of new shoot increments equal pramn(*,2) and pramx(*,2), respectively CLAYPG Number of soil layers that crop roots can occupy. The value used as CLAYPG for annual plants will vary from 1 on the day that plant growth starts to CLAYPG as read from the CROP option on day FRTC(3) of plant growth

EFRGRN(1)
Fraction of aboveground N which goes to grain FDGREM Fraction of standing dead (stdedc) removed by a grazing event over a one-month period. The daily removal rate is approximately fdgrem/30

FERAMT(1)
Amount of N to be added

FERAMT(2)
Amount of P to be added

FLGREM
Fraction of live shoots (aglivc) removed by a grazing event over a one-month period. The daily removal rate is approximately flgrem/30

FULCAN
Value of aboveground live C (aglivc) at full canopy cover, above which potential production is not reduced (above which there is no restriction on seedling growth)

PPDF(3)
Left curve shape for parameterization of a Poisson Density Function curve to simulate temperature effect on growth

PPDF(4)
Right curve shape for parameterization of a Poisson Density Function curve to simulate temperature effect on growth

PRDX(1)
Coefficient for calculating total monthly potential production as a function of solar radiation outside the atmosphere. It functions as a radiation use efficiency scalar on potential production. It reflects the relative genetic potential of the plant; larger PRDX(1) values indicate greater growth potential

PRDX(2)
Coefficient for calculating total monthly potential production as a function of solar radiation outside the atmosphere. It functions as a radiation use efficiency scalar on potential production. It reflects the relative genetic potential of the plant; larger PRDX (2)

PRDX(2)
Coefficient for calculating total monthly potential production as a function of solar radiation outside the atmosphere. It functions as a radiation use efficiency scalar on potential production.
It reflects the relative genetic potential of the plant; larger PRDX(2) values indicate greater growth potential

SNFXMX
Maximum symbiotic N fixation for forest (actual symbiotic N fixation will be less if available mineral N is sufficient for growth)