Predicting Soil Organic Carbon Dynamics of Integrated Crop-Livestock System in Brazil Using the CQESTR Model

Land degradation and reduction in productivity have resulted in losses of soil organic carbon (SOC) in agricultural areas in Brazil. Our objectives were to 1) evaluate the predictive performance of CQESTR model for a tropical savannah; and 2) examine the effect of integrated management systems, including Integrated Crop-Livestock System (ICLS) scenarios on SOC stocks. Two long-term paddocks, under similar edaphic and climate conditions were used in this study. In Paddock 4 (P4) the rotation was corn (Zea mays L.) and 3.5/4.5 years pasture (Urochloa ruziziensis), while rotations in Paddock 5 (P5) included 2.5 years of soybean (Glycine max L.), dryland rice (Oryza sativa L.), and corn followed by 2.5/3.5 years pasture (U. brizantha). Measured and CQESTR simulated values were significantly (0.0001) correlated (r = 0.94) with a mean square deviation (MSD) of 7.55, indicating that the model captured spatial-temporal dynamics of SOC. Predicted SOC increased by 18.0 and 12.04 Mg ha−1 at the rate of 0.90 and 0.60 Mg ha−1 year−1 under current ICLS management for P4 and P5, respectively, by 2039. ICLS increased soil C sequestration compared to simple grain cropping systems under both NT and CT due to high biomass input into the production system.


INTRODUCTION
Land-use changes or agricultural management practices can change soil organic carbon (SOC) stocks. Furthermore, the soil is an important carbon (C) pool and can be a source or sink of atmospheric CO 2 depending on the agriculture management adopted (Carvalho et al., 2010). In Brazil, unsustainable soil management has resulted in land degradation, loss of productive capacity of agricultural areas, loss of SOC, resulting in CO 2 emissions. It is estimated that 60-70% of Brazilian pastures show signs of degradation (Soares et al., 2020;Lapig, 2021), and up to 80% of the degraded pastures are in Central Brazil (Balbino et al., 2012). In general, pastures are considered degraded when they support very low stocking rates (<0.5 AU), show low plant cover, are invaded by nonpalatable plant species, and are often densely populated with termite mounds . Cardoso et al. (2016) report that the area of pasture required to produce 1 kg of a carcass (dead weight of animal) on a degraded pasture is approximately 320 m 2 but this falls to 71 m 2 on managed pastures, like mixed grass/legume. Integrated Crop-Livestock Systems (ICLS) have emerged as an economically feasible management to restore degraded pastures (Cortner et al., 2019;Reis et al., 2019). Reis et al. (2019) report that the payback period of ICLS is shorter than continuous livestock or continuous cropping, and Cortner et al. (2019) showed that ranchers perceived ICLS as a necessity to maintain their livelihood amidst declining profits. In Brazil, integrated systems or mixed farming systems have received considerable attention from farmers and policy makers. They were part of the Low Carbon Emission Agriculture Plan (ABC Plan, known in Brazil as "Plano ABC", Brasil, 2012), which is the contribution of the agriculture sector to the National Climate Change Program, also launched as part of Brazil's Nationally Appropriate Mitigation Actions (NAMAs) in 2009, at the 15th UNFCCC-Conference of the Parties (COP 15). The ABC Plan consisted of six selected agricultural practices, one of which was integrated systems, supported by specific credit lines at low-interest rates for the period 2010-2020 (Mozzer, 2011;Brasil, 2012). Despite difficulties in measuring the effectiveness of positive incentive measures in agriculture such as the ABC Plan (Carauta et al., 2021), following the Paris Agreement during the 21st UNFCCC-Conference of the Parties, Brazil has presented a new edition of the ABC Plan for 2020-2030 (Plano ABC+, Brasil, 2021a), part of which is also included in the country's Nationally Determined Contributions (NDC, Brasil, 2020;Brasil, 2021b). The new edition of the ABC Plan aims at the additional adoption of 10 million hectares under ICLS and other integrated systems by farmers.
Well-designed ICLS both spatially and temporally allow better regulation of biogeochemical cycles and showing synergies between crop and livestock systems in system-wide evaluations of production and environmental quality Lemaire et al., 2014). Long-term ICLS enables constant and efficient nutrient cycling because the animal, pasture, and crop residues release nutrients at different rates (Assmann et al., 2017). When grasses are intercropped, especially in the form of consortium (e.g., palisade grass planted during the same cycle as a grain crop), there is a higher biomass input and consequently higher nutrient recycling (i.e., N, P, K, and Mg) in the production system (Pariz et al., 2017). Measured nitrous oxide emissions are reduced in both ICLS (Sato et al., 2017) and ICLF (integrated crop-livestock-forestry), the most complex version of integrated systems with tree components (Franzluebbers et al., 2016;Carvalho et al., 2017). Integrated crop-livestock systems are also shown to improve nutrient cycling by re-coupling nitrogen (N) and C cycles (Ryschawy et al., 2017). Sant-Anna et al. (2017) investigated changes in SOC of pastures, crop production systems, and ICLS in the Brazilian Cerrado and found the highest SOC stocks under ICLS; however, not all ICLS increased SOC even under zero tillage. Oliveira et al. (2018) reported that soil N deficiency negatively affected SOC accumulation in ICLF and concluded that N management is a key to increasing SOC accumulation in integrated production systems. Damian et al. (2021) evaluated SOC changes in a poorly managed pasture with more intensively fertilized and diversified pasture systems (FP) and ICLS. They concluded that fertilization every year (FP) and the implementation of a cropping phase alternating with pasture (ICLS) resulted in the highest SOC stocks. When the pasture is recovered by inserting a crop (usually a grain crop), the residual fertilizer after cropping enhances pasture productivity. In this case, the cost of fertilizer application is covered by the income received for the crop (Kluthcouski and Yokoyama, 2003;Macedo, 2009). Meanwhile, perennial grasses such as brachiaria grass (Urochloa spp.) with large above-ground and root biomass contribute to SOC and increase in grain yields more than in a production system specialized for crop production only, particularly in years with poor rainfall distribution (Salton et al., 2014;Bieluczyk et al., 2020). Moraes Sá et al. (2014) reported a close correlation between SOC stocks and grain yield of soybean and corn from a long-term experiment in an Oxisol in southern Brazil under a subtropical humid climate. The added plant residues in ICLS differ in quantity and quality, being more recalcitrant compared with continuous pasture. Therefore, ICLS have been suggested as promising agricultural management contributing to the decarbonization of Brazilian agriculture (Tadini et al., 2021).
Assessing SOC dynamics in agroecosystems is challenging (Tornquist et al., 2009) mainly because of the complex interactions among components such as soil, vegetation, grazing animals, and humans (Godde et al., 2020). These changes often occur gradually and are difficult to detect in the short-term against the larger background (McGill et al., 1986;Ghani et al., 1996;Bolinder et al., 1999). Since the ICLS consists of a mixture of different grain crops often in rotation with several forage species (e.g., grass or legume), with varying grazing intensity and livestock densities, it can be difficult to determine whether the agricultural system is a C source or sink. Detecting small management-induced changes in SOC during short periods of time is difficult because of large differences in spatial and temporal SOC stocks (Kravchenko and Robertson, 2011). Long-term experiments, with historical datasets, have been used to assess impacts of past agricultural management practices on SOC dynamics. Process-based C models are potential research tools to predict the SOC changes; however, they need validation for the edaphoclimatic conditions of each region or country. Additionally, process-based C models are useful tools to analyze soil management options and to compare impacts of different management scenarios on SOC stocks. Thus, these models can be used to supplement field experiments to study SOC dynamics and to estimate the distribution of C in soil (Al-Adamat et al., 2007;Gollany et al., 2012). Furthermore, C models are useful in predicting the effects of potential management changes on the soil C stocks (Gollany et al., 2021). This makes it possible to test different scenarios and seek strategies to mitigate the negative impacts of such changes.
Carbon models have been used extensively for ecosystems and soil types under temperate conditions, while their evaluation under tropical conditions is less common (Kamoni et al., 2007;Tornquist et al., 2009;Damian et al., 2021). In these regions soil organic matter cycling is very different from that observed in temperate regions, because of the predominance of highly weathered acidic soils with low cation exchange capacity, high temperatures and high annual precipitation, where processes, including SOC decomposition rate, can be tenfold faster (Moreira and Siqueira, 2006). Given the importance of Brazil for global agriculture, it is necessary to validate C models and use them to investigate management options to increase SOC stocks in complex and diversified agriculture systems like integrated crop-livestock production.
Because of its relative simplicity, the readily available model inputs, and the possibility to compute SOC in a soil up to 5 layers, the CQESTR model was selected to study how different agricultural management practices affect soil C dynamics under tropical climate over time in two agroecosystems under ICLS. The CQESTR model was used previously to predict soil C dynamics in two tillage systems under tropical soils (Ultisol and Oxisol) in southeastern and northeastern Brazil (Leite et al., 2009); however, never used to predict SOC stocks under ICLS in Brazilian Cerrado. Due to the paucity of measured long-term SOC stocks data for the intensified and diversified pastures in Brazil, we hypothesize that process-based models can be an efficient and cost-effective tool to predict soil C stocks under several management practices and to project SOC stocks change for several ICLS scenarios. Therefore, the objectives of this study were to: 1) validate the CQESTR model for a tropical savanna (Cerrado) and predict the effect of several agricultural management systems and practices, including ICLS and no-tillage (NT); and 2) simulate the effect of conventional tillage and NT production scenarios on SOC dynamics in diversified ICLS.

Crop-Livestock Management
Two areas under ICLS and a nearby area under native vegetation (reference site) were studied. The areas under ICLS are part of a Technological and Research Reference Unit (TRRU) which is part of the ILPF Network Association (Rede ILPF, 2021) that aims to accelerate the widespread adoption of crop-livestock-forestry integration technologies by rural producers as part of an effort aimed at the sustainable intensification of Brazilian agriculture. This TRRU is composed of six areas with sizes ranging from 5.1 to 9.2 ha.
The study area was covered by natural vegetation until 1933 when selective logging began and lasted until 1950 ( Figure 1). Deforestation was completed in the 1970s when agriculture started. Common bean, dryland rice and corn were planted as main crops, including the areas currently under ICLS. In 1983 rice cultivation stopped and crop rotations changed to common bean and corn. In 1994, soybean was cultivated for the first time in the area, and 5 years later, soybean entered the rotation. Finally, in 1995, the area was divided into 6 paddocks as ICLS was gradually implemented between 2000-2005. Few variations in the rotation have occurred since 2000. The introduction of the pasture to the system occurred by planting corn in consortium with brachiaria grass. Two of the six long-term paddocks ( Figure 2A) were selected for this simulation study, P4 (7.5 ha) and P5 (8.1 ha), because these two paddocks are under a long-term observation study including greenhouse gas (CO 2 , water vapor and CH 4 ), soil, weather, and radiation sensors; these paddocks also have long-term soil data, including SOC. Initially, corn was the common crop used in the ICLS rotation, alternated with several years of pasture. The crop phase of the rotation became more diversified by introducing soybean, aerobic rice and common bean to the crop phase starting in 2007/2008, first in P5, then in P4 during 2013/2014. From 2014, there is a more pronounced presence of pasture in P5 compared to P4. The current rotation is based on the alternation of crops with a pasture phase ( Table 1). The rotation in P4 includes 3 years of summer crops (rainy season) with pasture during the dry season and 2.5 years of pure pasture. In P5 the rotation consists of 3 years of summer crops with pasture in the dry season, and 3.5 years of continuous pasture phase. The whole system is conducted under no-tillage and pasture is always reintroduced by sowing corn in consortium with brachiaria grass. From 2020 on, another option for reintroducing pasture after the crop phase is a dryland ricebrachiaria grass consortium. The grass remaining after corn or rice harvest initiates the pasture phase. In 2020/2021, the summer crop was substituted by pasture in P4 because of the COVID-19 pandemic. Soil preparation until 1995 included tillage with a conventional disc plow (~10-cm depth) and leveling harrow (first operation at~7-cm and a second at~3-cm depth). Since then, exclusively direct seeding (no-tillage, NT) was used. Fertilizers were applied according to soil fertility analysis and crop need at seeding or as top-dressing. Nitrogen fertilizer application in P4 varied from 0 to 288 kg ha −1 (averaged 49.6 kg N ha −1 ): phosphorus 38-116 kg ha −1 (average 73.7 kg P 2 O 5 ha −1 ); and potassium ranged from 38 to 159 kg ha −1 (average 60.4 kg K 2 O ha −1 ); while application rates in P5 varied from 0 to 120 kg ha −1 (average 28.6 kg N ha −1 ); phosphorous 26-183 kg ha −1 (average 86.2 kg P 2 O 5 ha −1 ); and potassium 31-129 kg ha −1 (average 65.2 kg K 2 O ha −1 ). Fungicides and insecticides were used when needed according to manufacturer's recommendation and control level. Glyphosate was used to control weeds and brachiaria grass growing to allow seeding in NT.

Soil Sampling and Analysis
In 2010, 2015 and 2020 samples were collected from the paddocks (P4 and P5) to 100 cm depth at 10-cm increments. Four main soil profiles were used in 4 quadrants of each studied site, and the samples were taken to 30-cm depth from the soil profile, and from 4 satellite points at a 5-m distance from each profile. A composite sample for each profile and depth was prepared from these subsamples. In the other years, 30 individual samples were collected at 10 or 30 cm (0-10 cm in 1999 and 2013; 0-10, 10-30 cm in 2007) in a zigzag pattern. In this study, SOC stocks at 0-10 and 10-30 cm depths were used for each sampling year. In the case of quadrant sampling, each composite sample of 5 sub-samples represented 1.87 ha in P4 and 2.02 ha in P5; however, if considering sub-samples as individuals, the sampling density was 2.66 and 2.47 samples per ha for P4 and P5, respectively.
In the case of zigzag sampling, each individual sample represented 0.25 and 0.27 ha, at the sampling density of 4.00 and 3.70 samples per ha, respectively. Machado et al. (2009) studied the spatial variability of SOC for 0-5, 5-10, and 10-20 cm in a 12.5 ha soybean field under NT in a clay (555 g kg −1 ) Rhodic Ferralsol (Latosssolo Vermelho distroférrico, in the Brazilian Soil Classification System), comparable to our study area. They found that the recommended sampling density was a minimum of 0.64 samples per ha considering sampling depth between 0 and 20-cm depth. Before soil chemical analysis, plant tissues and other nonsoil material was removed by hand and sieving. Air-dried soil samples were sieved in a 2 mm sieve. Soil samples were collected for SOC analysis using hand auger. The SOC was determined by the Walkley-Black method (Nelson and Sommers, 1996) without external heating, using sulfuric acid to generate internal heat for the reaction, and a correction factor of 1.3 to calculate SOC. The Walkley-Black method was used because the Dumas method (elemental analysis) was not available before 2010. We opted to use the same method for the sake of comparability of the SOC measurements.
The hydrometer method was used to determine clay, sand, and silt contents, with a standard hydrometer with Bouyoucos scale (Gee et al., 2002). Soil bulk density was determined with the soil core method using Kopecky rings (Grossman et al., 2002). Soil pH was determined in water (Thomas, 1996). Exchangeable calcium, magnesium and extractable potassium were extracted as described by Kuo (1996), then Ca and Mg were determined by atomic absorption spectroscopy and K by flame emission spectrometry (Wright and Stuczynski, 1996). Potential acidity was determined according to Silva (2009). Aluminum was extracted according to Bertsch and Bloom (1996) as modified by Silva (2009). Selected soil properties of the soil surface (0-10 cm depth) are shown in Table 2.

The CQESTR Model Data Inputs
Most data inputs required by CQESTR were provided by Embrapa and other required information was obtained from literature. Weather information, such as average daily air temperature and monthly rainfall, was provided by the weather station at Embrapa. Paddock 4 received 9.5 Animal Unit (AU) ha −1 and P5 received 8.8 AU ha −1 . The grazing rotations were 10 grazing days and 56 rest days per cycle. The amount of manure each paddock received was estimated at~13 and 12 Mg ha −1 year −1 for P4 and P5, respectively. Embrapa provided information on above-ground biomass for Urochloa grass production, while below-ground biomass was estimated based on similar grass species root:shoot ratios (Bolinder et al., 2002). These grass species root:shoot estimates were successfully used (r = 0.987) in previous CQESTR simulations of reed canary grass and switchgrass pastures (Dell et al., 2018). For crop yield, we used an extensive dataset from the National Food Supply Company ["Companhia Nacional de Abastecimento"] (Conab, 2021). This information and harvest index were used to calculate the amount of above-ground crop biomass after harvest which was used to develop crop-specific yield files for each vegetation type. Average yields were calculated for different time spans where factors, such as weather and variety improvements, created obvious delineations in the yields for each crop individually. Corn yields of 2,225; 4,481; 4,939 and 7,620 kg ha −1 , dryland rice yields of 965 and 2,069 kg ha −1 , soybean yields of 1,845 and 2,801 kg ha −1 , and common bean yields of 522; 1,718 and 2,738 kg ha −1 were used in these simulations (Conab, 2021). Weeds and volunteer crop growth during fallow periods was estimated at 1,800 kg ha −1 (Pacheco et al., 2011). A 4,266 kg ha −1 residue yield was used for any millet crops; and yields for Urochloa brizantha and U. ruziziensis ranging from 4,850-12,736 and 3,100-15,095 kg ha −1 , respectively, were used for the simulation by duration of plant growth ranging from 0.5 to 4 years, based on data provided by Embrapa. A single pass no-till disc planter was used for seeding at 5-cm depth, in all years except summers of 2003 and 2007 for P4 and P5, respectively ( Table 1). Prior to seeding of soybean in those years, two passes were made with harrows at 3-and 9cm depth. Other information required by CQESTR: below ground biomass, nitrogen content of residue at decomposition initiation, fraction of pre-tillage residue weight remaining on the soil surface after each tillage, were based on literature (Santos et al., 2007;Leite et al., 2009;Pacheco et al., 2011;Mauad et al., 2012). Number and thickness of soil depths, SOC content and soil bulk density of each soil depth were provided by Embrapa ( Table 3). The CQESTR input values used for the initial SOC in the 0-10 and 10-30 cm depths were 15 and 25 g kg −1 , and corresponding soil bulk densities were 1.34 and 1.35 g cm −3 , respectively, which were used for the model spin-up period dating back to native vegetation conversion to agricultural management in the 1970s as a starting point. Concentration of SOC (g C kg −1 ) was converted to Mg C ha −1 using bulk density measured for each depth (Table 3) to determine soil mass per depth by area.

CQESTR Validation
The exact rotation (as shown in Table 1) was used in the simulations for each of the paddocks. Corn, soybean, and dryland rice were grown in annual rotations from 1990 until 2003 except for an entire fallow crop year in 1994 for P4 and during 1994/1995 for P5. From 2003 to 2039, a 5-years ICLS cycle was used starting with a crop in 2004/2005 for P4 and a 6-years ICLS cycle was used starting with a crop in 2003/2004 for P5, respectively, typified by a first year of a corn intercropped with pasture followed by 3-or 4-years of pasture only. For 3 years following pasture, soybean, dryland rice and common bean were grown in different phases before restarting the cycle. The yield file used for each year during the simulation was crop-specific, as described above, and a soil operation file, created specifically for each management operation, as well as estimated animal manure, were used in the simulations. The SOC dynamics in the 0-10 and 10-30 cm depths were simulated for 50 years , and the simulation results were compared to observed SOC stocks in 1999SOC stocks in , 2007SOC stocks in , 2010SOC stocks in , 2013SOC stocks in , 2015 and 2020 for the 0-10 cm depth, and in 2007, 2010, 2015 (P5 only), and 2020 for the 10-30 cm depth.

CQESTR Simulation of Crop Rotation Change Scenarios
At the end of the validation period (1990-2020), SOC was simulated for additional management scenarios to represent the adoption of cropping systems and the transition to ICLS with different crop rotations for nearly two additional decades (2021-2039). Two crop-livestock rotation simulation scenarios were prepared for the two paddocks. Paddock 4 had a 5-years rotation, which included corn as a summer crop followed by 4.5 years pasture. The rotation in P5 was 6 years and included soybean followed by fallow in the first year, dryland rice and fallow in the second year and corn followed by 3.5 years pasture. Actual crop biomass yields estimated pasture residue inputs and animal manure were used in the simulations.
Additionally, four scenarios were prepared including soybean and corn in rotation. In the first scenario, the same crop sequence as P4/P5 was used, but tillage was changed from NT to conventional tillage (CT). The second scenario included soybean as a summer crop followed by corn in a second harvest in the same crop year under no-till (Soy-CS(NT)); in the third scenario, soybean and corn alternated each year as summer crops under no-till (Soy/Corn(NT)) followed by a fallow period; and finally in the fourth scenario, soybean and corn alternated each year as summer crops under conventional   tillage (Soy/Corn(CT)). Conventional tillage operations in the simulation consisted of tilling with a conventional disc plow (10cm depth) and leveling harrow (first operation at 7-cm and a second at 3-cm depth).

CQESTR Model Evaluation
Model performance was evaluated as described by Liang et al. (2009) using regression analysis and mean square deviation (MSD) statistics. The three components of MSD are squared bias (SB), nonunity slope (NU) and lack of correlation (LC), which are entirely independent and related to terms of the linear regression equation (Y = a + bX) and the regression coefficient (r 2 ) (Gauch et al., 2003).

Evaluation of CQESTR Model Performance
The linear fit of simulated vs. observed SOC explained 94% of the variation (Figure 3). Regression analysis of 18 pairs of observed and simulated SOC values for Paddock 4 (P4) and Paddock 5 (P5) was significantly (p < 0.0001) correlated (r = 0.94), with a slope of 0.84 not significantly different from 1.0 for soil sampling depths of 0-10 and 10-30 cm (Figure 3). The total mean square deviation (MSD) was 7.55. It indicates that CQESTR can accurately predict measured SOC stocks at the two soil sampling depths at both ICLS sites without calibration. The MSD was partitioned into its components: lack of correlation (LC), square bias (SB), and non-unity slope (NU), and were 7.55, 0.15, and 0.08 Mg SOC ha −1 , respectively. The lack of correlation (LC), the highest contributing component of MSD accounted for 97% of the total MSD, which indicated that prediction errors were associated with data scattering and high standard deviation of observed SOC data. The square bias (SB) accounted for 2%, and non-unity slope accounted for 1% of the total MSD.

CQESTR Simulated Rotation Scenarios
CQESTR predicted the measured SOC values for P4 at both depths very well (r = 0.94), except underestimating the SOC values by 3.30 Mg SOC ha −1 in the top 10-cm depth in 2020 (Figures 3, 4). CQESTR underestimated the measured SOC values for P5 by 4.14 and 4.33 Mg SOC ha −1 in the top 10-cm depth during 2010 and 2013 sampling periods and overestimated measured SOC values by 4.94 Mg SOC ha −1 in 10-30 cm during the 2015 sampling period. Measured and simulated SOC stocks increased in the topsoil for P4 more with time than for P5 ( Figure 3). The measured SOC values for P4 ( Figure 4) had a similar pattern to that in P5 ( Figure 5), although the measured and simulated SOC stocks increased faster in P4 than in P5.

CQESTR Projected Soil Organic Carbon Stock and Rotation Scenarios
A comparison of simulated SOC stocks changes over 40 years for Paddock 4 (P4) and Paddock 5 (P5) indicates more SOC increase in P4 than in P5 (

CQESTR Model Performance
The regression results illustrate that CQESTR accurately simulates long-term SOC dynamics under ICLS management for the tropical savannah conditions (Figure 3). The linear fit of simulated vs. observed SOC explained 94% of the variation with the MSD of 7.55, indicating that the model captured spatialtemporal dynamics of SOC in the topsoil well despite limited SOC data. Regression analysis for the long-term data from across North America (Liang et al., 2009), however, resulted in a higher r (0.98), feasibly due to a large number of pairs of predicted and measured data analyzed (306) and for the much wider range of SOM examined (7.3-57.9 g kg −1 ). The relatively small (7.32 Mg ha −1 ) lack of correlation (LC) between the observed and simulated values indicated that CQESTR prediction errors were mainly associated with data scattering. The factors such as the variability of SOC in the field, especially under paddock conditions, and sensitivity limits of the precision of SOC measurement related to sample collection, processing, and analysis. Leite et al. (2009) reported that square bias (SB) was the highest component of the MSD when simulating SOC under several tillage management systems under tropical conditions. CQESTR simulated SOC dynamics under different management systems in some soils better than others under tropical conditions (Leite et al., 2009).
Measured and simulated SOC stocks in the topsoil increased faster in P4 (Figure 4) than in P5. This increase in SOC stocks is most likely because Urochloa spp. is cultivated in P4 as a cover crop in the rotation for a longer duration than in P5. The measured SOC values in P4 ( Figure 4) were similar to that for P5 ( Figure 5). The rapid increase of SOC stocks in P4 and P5 is most likely due to the high capacity of these soils to retain SOC because of their high clay content ( Table 2) and the presence of iron and aluminum oxides, which results in C stabilization (Leite et al., 2009;Bayer et al., 2011). Also, it could be due to the insufficient removal of plant and dung residues from the samples during soil preparation as evident from the large standard deviation of the samples (Gollany et al., 2013). However, the large standard deviation could also occur because the studies areas were production sites and not small plot experiments with replication.
Higher SOC values in the top 10 cm for P4 in 2010 and 2015 were observed (Figure 3), although during the ICLS, P5 had the more complex crop rotation cycle in that period. Both areas were under similar edaphic conditions; however, P4 had more time in the pasture phase than P5 (Table 1) (4.5 vs. 3.5 years pasture). Another reason could be the soil disturbance at seeding events, even under NT, which occurred more often in P5 at the beginning of the ICLS. Crop intensification started in P4 after 2013/14. In 2019/20, common bean and rice + Urochloa were planted during the summer season, which added more crop residues to P4. The rice + Urochloa consortium is a modification in the rotation of P4 to increase biomass and C input into the system.
Additionally, P4 is smaller (7.5 ha) than P5 (8.1 ha) but has received the same total amount of manure. Therefore, it is expected to reach the pre-agricultural SOC stocks earlier than P5 (Figures 3, 4). Also, the influence of grazing cattle stimulates root growth and exudate production, which can modify the ratio of root and above-ground biomass and the quality of the C added to the soil (Bayer et al., 2011) and consequently influence C stocks and soil organic matter decomposition in the soil profile.
Paddock 5 received more legumes, typically soybean, than P4 and this could be another reason for the difference in SOC stocks of the two paddocks. Soybean with high biological N fixation efficiency is usually cultivated without N fertilizer, using only inoculation with bacteria of the genus Rhizobium spp. Most of this N is exported in grains, and the negative or null net N balance in cropping systems reduces biomass yields and C accumulation  (Sisti et al., 2004;Urquiaga et al., 2006;Bayer et al., 2011). These authors recommend the introduction of legumes as green manure in the rotation to promote soil C sequestration. According to Sisti et al. (2004) and Souza et al. (2009), N from legumes that biologically fix N can promote more C accumulation than N from mineral sources. These studies, however, were carried out in long-term experiments of grain crop systems. Our understanding of the C and N cycles in complex systems such as ICLS is still limited (Soussana and Lemaire, 2014).

CQESTR Simulated Scenarios
A comparison among the simulated scenarios by 2039 indicate differences in SOC stocks at the 0-10 cm depth (  (Table 5), relative to P4 under NT. This is likely because of less soil disturbance and residue decomposition, consequently reduced mineralization rate and less C losses under NT compared to CT. The rate of SOC stocks changes in the 0-30 cm soil depth for ICLS(NT) and ICLS(CT) in P4 were 0.90 and 0.82 Mg ha −1 year −1 , respectively ( Table 4). This is consistent with other soil tillage studies. No-till promotes slower crop residue turnover and mineralization than CT (Sherrod et al., 2003), and less SOC loss (West and Post, 2002). In P4 under ICLS, especially under NT, most of the simulated SOC increase occurred in the 0-10 cm soil depth, while little changes in SOC were predicted for the 10-30 cm soil depth SOC, leading to more SOC in the top 10 cm than in the underlying 20-cm layer by 2039. The changes in SOC were due to the difference in the SOC accumulation rates of 0.82 vs. 0.08 Mg ha −1 year −1 for the 0-10 cm and 10-30 cm depths, respectively. Whereas in P5, larger SOC increases predicted at 10-30 cm depth in all the simulations. One possible explanation can be the deposition of large amounts of stubble in the ICLS, which under NT stays on the surface and is slowly incorporated into the soil. The stratification of SOC in NT is a widely reported characteristic of this practice (Mrabet, 2002;Sá and Lal, 2009). The proportionally longer participation of the pasture phase in P5 compared to P4 in the 2012-2020 period, that allowed better grass root development, is likely the reason for the higher SOC accumulation rate at the 10-30 cm depth in P5 (0.41 Mg ha −1 year −1 ), compared to P4 (0.08 Mg ha −1 year −1 ). The higher SOC accumulation rate in the lower soil layer under CT can be explained by the incorporation of the surface litter into the lower soil layer. The role of roots in SOC accumulation or decomposition is still not completely understood (Dijkstra et al., 2020). On one hand, rhizodeposition is a great contributor to SOC stabilization, mainly through stimulating microorganisms. Mineral associated organic matter (MAOM) is protected from further microbial decomposition, and according to Cotrufo et al. (2013) MAOM is predominantly formed from microbial products. On the other hand, plant roots are suggested to be responsible for the destabilization of SOC in a process called rhizosphere priming effect, that is, rhizodeposition is used as substrate by a group of microbes, enhancing SOC decomposition (Huo et al., 2017). Nitrogen uptake by plants can increase competition with microbes that can further stimulate SOC decomposition by microorganisms. In fact, the contribution of the root system to SOC is the result of a balance between its SOC stabilizing and reactivating/ priming effect.
Organic matter formation and persistence in soils depend on environmental conditions, soil microbiota, the quality of soil minerals, and soil chemical and physical properties, especially soil structure (Hunt et al., 2020). Soil C stabilization and storage can occur by chemical stabilization, biochemical resistance, and physical protection. The adoption of NT increases water retention and lowers soil temperature, which favors microbial activity. Furthermore, fewer soil perturbations, favor build-up of larger soil aggregates (Madari et al., 2005), conferring more physical protection (Barreto et al., 2009). The amount of crop residue is of key importance in accumulating soil C, particularly in the Cerrado biome where high temperatures and humidity during the rainy season do not favor crop residue maintenance over the soil surface. Soil organic C stocks were predicted to increase only by 2.42 and 1.89 Mg ha −1 at 0-30 cm between 2019 and 2039 for Paddock 4 and 5, respectively, under Soy/ Corn(NT), which has alternating years of soybean and corn crops ( Table 4). Introducing corn as a second harvest in the same year (Soy-CS (NT)) which resulted in somewhat higher annual biomass input due to corn residues (13.3 Mg ha −1 year −1 vs. 9.5 in Soy-CS and Soy/Corn, respectively), did not increase SOC stocks substantially (4.25 and 2.25 Mg ha −1 in P4 and P5, respectively), especially compared to ICLS (Tables 4 and 5). This indicates that providing N through grain legumes without large increase in crop biomass will not result in substantial soil C accrual. Also, as mentioned before, in the case of grain legumes, like soybean, most of the N is removed with the grain, which results in no positive N balance in the system (Sisti et al., 2004). It is well known that positive N balance in the soil is necessary to achieve C accumulation due to a narrow range (10-14) of soil C:N ratio in most soils.
Therefore, even under NT, low disturbance systems such as Soy/ CS(NT) with relatively low residue inputs to the soil are less likely to improve SOC as in ICLS with high residue and manure inputs. Soil organic C stocks loss are predicted for scenarios that have annual alternating corn and soybean crop with fallow under conventional tillage (Soy/Corn(CT)). Decreases in SOC stocks of 2.78 and 3.34 Mg ha −1 at rates of −0.14 and −0.17 Mg ha −1 year −1 in the 0-30 cm soil depth under the Soy/Corn(CT) scenario; respectively, for P4 and P5 are predicted. Stockmann et al. (2013) reported that SOC negatively correlated with tillage. The above discussion shows that single promising conservative management practice adoption will not necessarily result in soil C sequestration; therefore, interactions of all the components need to be considered when managing ICLS (Valkama et al., 2020).
In summary, the CQESTR model predicted an increase in SOC stocks for both the NT and CT scenario under ICLS. Predicted SOC increased by 18.0 (28%) and 12.04 Mg ha −1 (19%) at the rate of 0.90 and 0.60 Mg ha −1 year −1 under current ICLS management for Paddock 4 and Paddock 5, respectively, by 2039. In single crop rotations under NT (i.e., Soy-CS and Soy/Corn) SOC accumulation at 0-30 cm was still predicted (between 6.67 and 2.92%), but in Soy/Corn under NT, at the 0-10 cm depth, SOC loss was predicted. Clearly, SOC accumulation in ICLS was favored by the pasture phase, and by introducing brachiaria grass in the rotation.

CONCLUSION
The CQESTR model was validated for the edaphoclimatic conditions of the Cerrado biome in Brazil, and successfully predicted the effect of several agricultural management practices in Integrated Crop-Livestock Systems (ICLS). The model captured spatial-temporal dynamics of SOC very well. The CQESTR predicted SOC increases by 18.0 (28%) and 12.04 Mg ha −1 (19%) for Paddock 4 (with the long pasture phase of Urochloa spp.), and Paddock 5 (with the shorter pasture phase of Urochloa spp.), respectively, by 2039. The use of the extended pasture phase without the crop phase was found to be the best management to increase carbon stocks and could assist Brazilian national initiatives aimed at restoring degraded pasture areas (i.e., "ABC Plan"), as well as reducing carbon dioxide emission from the soil.
Because of the limited measured long-term data for the ICLS under tropical climate, process-based models can be a costeffective tool to predict soil C stocks change under several ICLS scenarios, and to analyze soil management options and compare impacts of each scenario on SOC stocks.
The grass root biomass and root distribution under tropical savanna conditions may diverge widely from the characteristics of temperate grasses. CQESTR's estimation of SOC stocks could be improved with site-specific below-ground biomass and root distribution data. More long-term studies, SOC data, and root biomass for tropical grasses from diverse tropical biomes are needed to improve SOC stocks prediction. Furthermore, in comparing different management with complex rotations such as ICLS, soil samples must be taken during the same phases of the ICLS rotation. Careful sampling and sample cleaning can reduce plant and animal residue and reduce particulate organic matter or light fraction C incorrectly quantified as part of the stable SOC pool of the sample and improve long-term SOC stocks prediction.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JO: sampling collection and analysis, data collection and curation, investigating, methodology, validation, visualization, formal analysis, writing, review and editing. HG: conceptualization, data curation, CQESTR software, formal analysis, funding, investigating, methodology, project administration, resources, supervision, validation, visualization, formal analysis, writing-original draft, writing-review and editing. RP, data curation, visualization, writing-review and editing. LL, writing-review and editing. BM: data curation, formal analysis, funding, investigating, methodology, project administration, resources, supervision, writing-review and editing. MC: data curation, formal analysis, funding, methodology, project administration, resources, supervision, writing-review and editing. PM: data curation, formal analysis, funding, investigating, methodology, project administration, resources, supervision, writing-review and editing.