Characterizing Livestock Production in Portuguese Sown Rainfed Grasslands: Applying the Inverse Approach to a Process-Based Model

: Grasslands are a crucial resource that supports animal grazing and provides other ecosystem services. We estimated the main properties of Portuguese sown biodiverse permanent pastures rich in legumes (SBP) starting from measured data for soil organic carbon (SOC) and using the Rothamsted Carbon Model. Starting from a dataset of SOC, aboveground production (AGP) and stocking rates (SR) in SBP, we used an inverse approach to estimate root to shoot (RS) ratios, livestock dung (LD), livestock intake (LI) and the ratio between easily decomposable and resistant plant material. Results for the best ﬁt show that AGP and belowground productivity is approximately the same (RS is equal to 0.96). Animals consume 61% of the AGP, which is within the acceptable range of protein and energy intake. Carbon inputs from dung are also within the range found in the literature (1.53 t C/livestock unit). Inputs from litter are equally distributed between decomposable and resistant material. We applied these parameters in RothC for a dataset from different sites that only comprises SOC to calculate AGP and SR. AGP and SR were consistently lower in this case, because these pastures did not receive adequate technical support. These results highlight the mechanisms for carbon sequestration in SBP.


Introduction
Livestock provide an important share of human nutrition, contributing to feeding the world population [1]. Grasslands are a critical resource for animal grazing. Grasslands also provide other ecosystem services, such as carbon sequestration (through soil organic carbon-SOC-accumulation), and support for fauna and flora biodiversity [2][3][4]. Carbon sequestration provides climate regulation and mitigation services [5]. These are a particularly relevant contribution for the sustainability of livestock production systems, where they are a countereffect to the large greenhouse gas emissions from enteric fermentation, manure management of production of feeds [6,7]. However, the pace and saturation of carbon storage in soils depends on the floristic composition and soil properties of grasslands, and their quality differs significantly considering grass and grazing management practices, climatic and site conditions, as well as presence or absence of tree cover [8,9].
Knowledge of grassland/pastures properties is fundamental to improving our understanding of this production system and managing these ecosystems towards the provision of multiple ecosystem

Materials and Methods
In this work, we used the RothC model to estimate SOC dynamics in Portuguese SBP. RothC was developed to model the carbon turnover in arable soils in the United Kingdom, but it has been expanded and successfully applied to model soil carbon dynamics also in grassland [35]. It takes into account the effects of temperature, moisture content and soil type. The model uses a monthly temporal step and is divided into five compartments: easily decomposable plant material (DPM), resistant plant material (RPM), microbial biomass, humified organic matter and inert organic matter. The inert pool is resistant to decomposition and does not change over time [36]. To initialize the model, the model distributes initial SOC through the five compartments according to Weihermüller et al. [37], which depends on SOC and clay fraction. The full list of parameters required is: initial SOC content (t C/ha) and clay content of the soil (%), soil cover period, monthly input of plant residues (t C/ha), monthly farmyard manure carbon inputs (t C/ha) and monthly irrigation (mm), average monthly mean air temperature ( • C), monthly rainfall (mm) and monthly open pan evaporation (mm). RothC uses only one parameter that depends strictly on the land use type, i.e., ratio between DPM and RPM (DPM/RPM ratio). This parameter is key in determining soil respiration, which is calculated differently for each SOC pool. It can vary between 0.25 (to forest land uses) and 1.44 (to annual croplands) [29]. For pastures, there is a wide range of variation depending on vegetation and management. Grassland DPM/RPM ratio can range between 0.67 for unimproved grassland and scrub and 1.44 for improved grassland [29].
We built an optimization problem (Section 2.1) to: (a) obtain the RS ratio, animal dung per livestock unit (LU) and animal intake (LI); and (b) use the parameters obtained previously to estimate AGP and SR in the cases where only SOM data were measured (Section 2.2).

Data Characterization
To calibrate the model, we used data from the AGRO 87 [38] and AGRO 71 [22] projects. All data used in this paper is available in Supporting Information File S1.
The AGRO 87 project was the longest to date dedicated to the SBP production system (4 production years-from autumn 2001 to autumn 2005), covered the most SBP farms (six in total located in two Portuguese NUTS II regions-Alentejo and Beira Interior), and measured the most relevant data. Data collected in this project included SOM content (%), SR (LU/ha) and AGP (t dry matter (DM)/ha), as indicated in Table 1. The six farms surveyed were farm 1 (near Vaiamonte); farm 2 (near Elvas), farm 3 (near Pavia); farm 4 (near Santiago do Cacém), farm 5 (near Coruche) and farm 6 (near Covilhã). AGRO 71 also covered SBP, but only in two farms and three production years (from autumn 2001 to autumn 2004). In this project, only SOM concentration was measured. The two farms are farm 7 (near Castro Verde) and farm 8 (near Mértola). All farms are in central and southern Portugal in the hot-summer Mediterranean climate region, according to the Köppen climate classification system [39,40]. Higher temperature and lower precipitation were found in farm 1 (the southernmost farm), and lower temperature and higher precipitation were observed in farm 4 (the northernmost farm). The location of the farms and the variables measured are presented in Section 2.2. Table 1. Summary information for each farm in Figure 2, including the project, period covered, measured data, temperature and precipitation.  SOM was determined through soil sampling analysis, at 10 cm depth. For each parcel of the project, one composite sample was obtained (mix of a variable number of sub-samples collected throughout each parcel). Samples were dried overnight at 35-37 • C and crumbled mechanically and passed through a 2 mm stainless steel sieve before being measure the organic matter content [22,38]. AGP was measured three or four times per season. Quadrants of 1m × 1m were used in each parcel. All plants inside were harvested, oven-dried at 65 • C for 48 h, and then weighted [41]. The number of animals (cows and/or sheep) per day in each of the parcels was registered by the farmers, and then averaged in a year (the days without grazing were also registered) [22,38]. AGP and SR on average for all farms and years, are approximately 7 t DM/ha and 1 LU/ha, respectively. Organic matter increments were, in the period between the autumn of 2001 (used as starting point) and the autumn of 2005, approximately 0.2% per year.

Model Inputs Calculation
SOM content was converted to SOC content by multiplying it by the average carbon fraction in organic matter (0.58 [42]), soil density (g/m 3 ) and soil depth (0.1 m). Soil density, as well as clay content (needed for RothC calibration, as described later), were obtained from the LUCAS-topsoil database [43]. LUCAS is a European-wide project where topsoil field measurements were taken from approximately 20,000 sites in European Union countries.
Additional parameters needed for RothC calibration are the soil cover period, monthly input of plant residues (t C/ha), monthly farmyard manure carbon inputs (t C/ha) and monthly irrigation (mm). Farmyard manure was not applied, and all parcels were rainfed pastures. Soil cover period is a binary monthly variable, where 1 means that the soil was covered with vegetation during that month and 0 means that the soil was uncovered. The soil is covered between September and June and it is uncovered in July and August, as it is good agronomic practice to fully graze the pasture before Summer to facilitate sward reestablishment through self-reseeding in the fall. Carbon inputs from plant residue (I plant ) and animal dung (I animal ) were calculated according to where RS is the root to shoot ratio, LD is the livestock dung excreted (t C/LU), LI is proportional livestock intake (kg DM eaten/kg DM pasture yield), AGP is the aboveground productivity (kg DM/ha), SR is the stocking rate (LU/ha), and CF is the carbon fraction of the legumes and grasses (0.4 t C/t DM) obtained from the Intergovernmental Panel on Climate Change (IPCC) [14].
Monthly air temperature and precipitation were respectively obtained from the "Global Precipitation Climatology Project (GPCP)" [44] and the Land Processes Distributed Active Archive Center (LP DAAC) project [45]. Monthly open pan evaporation was calculated as 75% of the potential evaporation, an assumption used by previous studies using RothC [46,47]. Potential evaporation was calculated using the Thornthwaite equation [48], which is a function of monthly average temperature and annual average temperature, number of days and average day length per month.

Optimization Procedure
Using the measured and third-party data described thus far, we implemented an optimization procedure to estimate the remaining parameters needed for RothC that were not measured and that cannot be estimated indirectly from any additional data sources. Those parameters are: RS ratio (belowground productivity (BGP)/AGP); LD (t C/LU); LI (livestock intake-kg DM/kg DM); and DPM/RPM. The estimated parameters are the same for all farms and production years. The determination of these parameters was the ultimate goal of the optimization procedure. Figure 1 depicts the three-step optimization procedure, which includes a data imputation procedure to fill in data gaps: (1) Estimate parameters using only data points without missing data, i.e., estimate RS, LI, LD and DPM/RPM for the AGRO 87 and 71 farms and years without data gaps, i.e., where AGP, SR and SOM were all measured. (2) Use the SOM data and parameters obtained to estimate the required AGP and SR for the SOM increments in the remaining farms and years where AGP and SR from AGRO 87 and 71 data were missing. (3) Estimate RS, LI, LD and DPM/RPM ratio using the complete data set (i.e., using field-measured AGP and SR used in the first step and AGP and SR estimated in the second step). With the new set of parameters, we estimated AGP and SR again.
Steps 2 and 3 are repeated until the values for all the parameters stabilize. The procedure (steps 1 to 3) was repeated in 100 iterative runs with different random initial conditions. Each run ended when the stop condition of the loss function was achieved. The set of parameters in the last iteration before the end of the run was recorded. The problem was formulated as where N is the total number of data points, i.e., the number of production years and farms. The algorithm minimizes the difference between field measured SOC stock (SOC measured ) and estimated SOC stock (SOC estimated ), in relative terms (divided by the field SOC stock), subject to constraints on acceptable values of the parameters gathered through literature revision. In each run, initialization parameters could range freely in the space defined by the problem constraints. RS ranged between 0.5 and 8 [13,14], LI between 0 and 1 kg DM/kg DM, indicating no grazing to total AGP grazing by animals. LD varied between 0.8 t C/ha.LU and 2.5 t C/haLU, depending on the animal size and C/N ratio of the dung [14], and the DPM/RPM ratio for grasslands ranged between 0.6 and 1.44.
subject to: where N is the total number of data points, i.e., the number of production years and farms. The algorithm minimizes the difference between field measured SOC stock (SOCmeasured) and estimated SOC stock (SOCestimated), in relative terms (divided by the field SOC stock), subject to constraints on acceptable values of the parameters gathered through literature revision. In each run, initialization parameters could range freely in the space defined by the problem constraints. RS ranged between 0.5 and 8 [13,14], LI between 0 and 1 kg DM/kg DM, indicating no grazing to total AGP grazing by animals. LD varied between 0.8 t C/ha.LU and 2.5 t C/haLU, depending on the animal size and C/N ratio of the dung [14], and the DPM/RPM ratio for grasslands ranged between 0.6 and 1.44.  Equation (3) is a nonlinear programming (NLP) problem. The objective function is a nonlinear function, and it is subject to linear equality and linear inequality constraints. Thus, the feasible region is a convex polytope. We used Matlab vR2017a to perform the optimization process, using the "fmincon" function. This function searches the feasible space for the minimum function value using an interior-point optimization algorithm (default), which solves NLP problems. It has been used in large, sparse problems, as well as small dense problems. The algorithm satisfies the boundary conditions for all iterations and can recover from non-admissible results or infinite results. We increased by 7 times the maximum function evaluations to 20,000 (default value is 3000), which allowed us to increase the search space. The stop condition was 10 −6 , which gave a very small search step, since the absolute value of the loss function ranges between 0.5 and 1. For validation of these choices, we ran the optimization procedure five times to ensure that the final optimum parameter set was always the same, i.e., better coverage of the variable space, indicating that the minimum of the loss function is global rather than local.
After arriving at the optimum set of parameters, we analyzed the relationship between the randomly selected initial parameter set and the final set of parameters after optimization. The goal was to ensure that the model was adjusting all parameters rather than fixate one or more parameters around the initial selection. We also performed a correlation analysis between the parameters of both Sustainability 2018, 10, 4437 7 of 21 systems, using Spearman's ρ [49] method. Spearman's ρ is a nonparametric measure of monotonic statistical dependence between two variables. This method does not make any assumptions about the distribution of the variables. We used IBM SPSS 25 to perform correlation analysis.

Estimation of Aboveground Productivity and Stocking Rate for an Additional SOM Data Set
After estimating RS, LD, LI and DPM/RPM of SBP, we applied them in the distinct case of the EDP project [25]. This project was sponsored by Energias de Portugal, EDP (the main Portuguese electrical utility) with the aim of demonstrating how land use activities could help Portugal comply with its Kyoto Protocol target. SBP installation and management was included in the land uses covered (alongside afforestation, forest management, and no-tillage in croplands). Soil samples were collected during the project, but only SOM was measured. SBP was the only pasture system covered, in seven production years (autumn 2006 to autumn 2012) for six farms, for a total of 77 parcels. The farms were: farm 6 (also present in AGRO 87 project), farm 9 (near Samora Correia), farm 10 (near Coruche), farm 11 (near Coruche); farm 12 (near Grândola) and farm 13 (near Montemor-o-Novo). Here, AGP and SR were unknown, only SOM was measured (332 soil samples). We used the same optimization problem in Equation (3) to minimize the difference between estimated and measured SOC stock. In this case, the optimization variables were not the pasture-specific parameters, but rather AGP and SR. Estimated AGP and SR, in each farm and production year (i.e., 48 parameters in total-24 AGP and SRs) were iteratively adjusted to minimize the loss function value, using the parameters previously obtained in Section 2.1. We also initialized the optimization process for 100 iterations to find the parameter set which leads to the lower function loss value in each iteration and to try to avoid local minima. We compared results with the AGRO project. We also performed a correlation analysis using Spearman's ρ method between AGP and SR. Table 1 summarizes measured data and main climatic variables, i.e., mean annual temperature and annual precipitation, for all farms covered in this paper. Here, AGP and SR were unknown, only SOM was measured (332 soil samples). We used the same optimization problem in Equation (3) to minimize the difference between estimated and measured SOC stock. In this case, the optimization variables were not the pasture-specific parameters, but rather AGP and SR. Estimated AGP and SR, in each farm and production year (i.e., 48 parameters in total-24 AGP and SRs) were iteratively adjusted to minimize the loss function value, using the parameters previously obtained in Section 2.1. We also initialized the optimization process for 100 iterations to find the parameter set which leads to the lower function loss value in each iteration and to try to avoid local minima. We compared results with the AGRO project. We also performed a correlation analysis using Spearman's ρ method between AGP and SR. Table 1 summarizes measured data and main climatic variables, i.e., mean annual temperature and annual precipitation, for all farms covered in this paper. Figure 3 depicts the score (i.e., loss function value) as a function of the optimum in each iteration RS ratio, LD, LI and DPM/RPM ratio, where each iteration corresponds to a different initial set of parameters given as input for the optimization process to run. Parameters in the "best" iteration, i.e., the iteration with the lowest loss function value, which was equal to 0.61 (represented in red in The final values of RS, LI and DPM/RPM for all runs are always ±5% with respect to the global optimum. The variation is higher for LD: 11%. The variability of the final set of parameters in each iteration does not influence predicted SOC dynamics significantly. As an example, we ran the model during a ten-year period starting from an arbitrary point with initial SOC equal to 25 t C/ha and considering an AGP of 7 t DM/ha and a SR of 1.0 LU/ha, which were the average of the AGROs project, as well as average climate conditions. The average difference of the estimated average annual C accumulation using any two sets of final parameters is 10 −10 t C/ha.yr. For LD, we originally provided the model with rather broad constraints on the parameters based on IPCC data. That resulted in LD ranging between 1.30 and 1.56 t C/LU (Figure 3c). Deposition usually ranges between 80 and 120 kg N/LU in the most common breeds of cattle used in Portugal [50].

Estimation of Production Parameters for the Pasture Systems
Considering that the C/N ratio of cattle dung in Portugal ranges between 15 and 25 [51], and using this C/N ratio to convert from mass of N into mass of C, this leads to LD in the range 1.2 to 3 t C/LU. The repetition of the optimization procedure produced differences between the best runs that were never higher than 2% in terms of the minimum loss function value. The parameters never deviated more than approximately 4%. Therefore, there is no evidence that increasing the number of iterations would change results, and it can be expected that the minimum reported here is the global minimum of the loss function. Figure 4 depicts the final parameters after optimization per iteration as a function of the respective randomly selected initial value. Through visual inspection, it is clear that there is no relationship between the two variables, for all parameters (e.g., the initial LI varied between 0 and 1, The RS ratio of 0.98 is close to values reported in the literature for "temperate arid shrubland" (value of 1.063). Another trend shown by the results in this work, as well as previous studies [13], is that when AGP is lower, the RS ratio tends to increase. Higher LIs are observed during production years and in farms with the highest AGP, but there are exceptions. For example, AGP in the last production year in farm 2 is lower than AGP in the second production year in farm 5, but the LI is higher (farm 2: 0.63 t DM/ha; farm 5: 0.62 t DM/ha). LI varies between 1.4 t DM/ha in the second year in farm 2 and 6.1 t DM/ha in the first year in farm 1. Farm 1 is where animals were fed the most on pasture (5.87 t DM/ha, on average), whereas farm 5 is the one where animals were fed the least on pasture (2.48 t DM/ha).
For LD, we originally provided the model with rather broad constraints on the parameters based on IPCC data. That resulted in LD ranging between 1.30 and 1.56 t C/LU (Figure 3c). Deposition usually ranges between 80 and 120 kg N/LU in the most common breeds of cattle used in Portugal [50].
Considering that the C/N ratio of cattle dung in Portugal ranges between 15 and 25 [51], and using this C/N ratio to convert from mass of N into mass of C, this leads to LD in the range 1.2 to 3 t C/LU.
The repetition of the optimization procedure produced differences between the best runs that were never higher than 2% in terms of the minimum loss function value. The parameters never deviated more than approximately 4%. Therefore, there is no evidence that increasing the number of iterations would change results, and it can be expected that the minimum reported here is the global minimum of the loss function. Figure 4 depicts the final parameters after optimization per iteration as a function of the respective randomly selected initial value. Through visual inspection, it is clear that there is no relationship between the two variables, for all parameters (e.g., the initial LI varied between 0 and 1, but the final value only ranged between 0.60 and 0.63- Figure 4b). The outliers mentioned above are also unrelated to the initial parameters Additionally, the best iteration never corresponds to a corner solution within the range of the variable.  Regarding the correlation analysis, all parameters are significantly correlated (at 5% level) except for the DPM/RPM ratio. DPM/RPM is the least correlated parameter (average: 0.120), e.g., the highest DPM/RPM correlation is with LD, 0.306. The highest correlation, in absolute value, is between the RS and LD (−0.814). The negative sign means that there is a tradeoff between LD and RS. SBP have high carbon inputs due to high AGP, and the model only requires BGP or animal inputs to arrive at the measured level of OM. The second highest correlation is between LI and LD (0.343).
When we used the optimum parameters to simulate SOC dynamics in the original farms, we found that the model overestimates measured SOC by 1 ± 1.5 t C/ha (minimum: 0; maximum: 17). Figure 5 shows the fit for each farm. Disregarding the outlier farm (farm 4) and the last production year of farm 2 (SOM increase of 2.7%-average: 0.4%), the simulated SOC underestimates measured SOC, by only 0.1 ± 0.7 t C/ha. The quality of the fit is a function of the farm and of the year. Examples in which the Regarding the correlation analysis, all parameters are significantly correlated (at 5% level) except for the DPM/RPM ratio. DPM/RPM is the least correlated parameter (average: 0.120), e.g., the highest DPM/RPM correlation is with LD, 0.306. The highest correlation, in absolute value, is between the RS and LD (−0.814). The negative sign means that there is a tradeoff between LD and RS. SBP have high carbon inputs due to high AGP, and the model only requires BGP or animal inputs to arrive at the measured level of OM. The second highest correlation is between LI and LD (0.343).
When we used the optimum parameters to simulate SOC dynamics in the original farms, we found that the model overestimates measured SOC by 1 ± 1.5 t C/ha (minimum: 0; maximum: 17). Figure 5 shows the fit for each farm. Disregarding the outlier farm (farm 4) and the last production year of farm 2 (SOM increase of 2.7%-average: 0.4%), the simulated SOC underestimates measured SOC, by only 0.1 ± 0.7 t C/ha. The quality of the fit is a function of the farm and of the year. Examples in which the results are a function of the farm are farms 2 and 6 during the second production year. In farm 2, the estimated SOC stock is only 1 t C/ha lower than measured SOC in farm 1 (Figure 5b), but 6 t C/ha higher (Figure 5f) in farm 6. An example of the influence of the characteristics of the year is farm 6 (Figure 5f). In the first year, the difference between the estimated and measured SOC stock is about 2 t C/ha, but in the second year the difference is about 6 t C/ha. The difference is higher in the last year on farm 2 (17 t C/ha). Regarding farm 1 (Figure 5a), SOC stock is always underestimated on average by 4 t C/ha in the second year due to the fact that measured SOC stocks in this farm are above the range of stocks used to calibrate the model, and as such they are outliers. Nevertheless, in the first year the difference between field measured and estimated SOC stock is close to zero. Farm 1 and farm 2 had similar measured SRs, but annual SOM increments were twice as high in farm 1 (0.54% vs. 0.23%).   Figure 6a presents the estimated SOC stock as a function of field-measured SOC stock, after application of RothC using the previously estimated parameters estimated to the farms involved in the EDP project. Figure 6a also presents the linear regression fits, where the intercept is significantly different from zero (p-value equal to 2.1 × 10 −6 ), and the slope is significantly different of 1 (p-value equal to 6.9 × 10 −10 ), which means that our model has a non-negligible bias. In the non-biased regression, the intercept would be zero and the slope would be one. In this case, the intercept is about 4 t C/ha and the slope approximately 0.8, meaning that the model overestimates the lower SOC stocks Pasture productivity in farm 1 was estimated, rather than field-measured as in farm 2, but the estimation provided an AGP close to the maximum observed in the AGRO 87 project (average estimated productivity in farm 1 was 9.5 t DM/ha; the maximum for all years and farms was 10.6 t DM/ha). Regarding farm 6, a reduction of about 2 t C/ha SOC was measured in the third year, but the model estimates an increase by the same amount. This was because, given the AGP and SR measured, an increase in SOC was expected. During the third year, measured AGP and SR in this farm are higher than AGP and SR during the first year on farm 5 (AGP: 5.8 t DM/ha versus 3.9 t DM/ha; SR: 1.29 versus 1.18), but in farm 5 SOC increases and in farm 6 it decreases. Weather data cannot justify this difference, because average monthly temperature is similar in both cases (temperature: 17 • C versus 18 • C) and annual accumulated precipitation is lower in the case of farm 5 (precipitation: 475 mm versus 590 mm). Thus, in the second year of farm 6, climatic conditions were more favorable for SOC accumulation than in the second year of farm 5.

Estimating AGP and SR for an Unrelated SOC Database
Multiple factors could justify these outliers. First, it could happen due to an experimental error during field measurements or in their laboratorial analysis. Second, SBP in Alentejo region are complex agroforestry systems in "Montado" areas of low-density cork or holm oak forest, where trees can affect SOC changes due to extra deposited litter (from leaves and dead branches) on the field or the influence of root systems [52,53]. Data points were selected as much as possible away from the influence of tree cover, but it cannot be fully guaranteed that trees did not affect measurements of AGP and SOM data. Even if trees did not influence results, animal grazing could have played a role. Only average SR for entire plots were recorded, but grazing and consequent dung accumulation is uneven throughout the plot due to gradients of plant establishment that are selectively grazed. Sample points could have been taken from a non-representative site where grazing was particularly more or less intense. Additionally, the nature of the organic materials in these plots and the dynamics of their assimilation into soils could have made our assumption that C input into soils and their conversion to SOC are simultaneous. We assumed that litter is transformed to SOC in the same month that it is produced. This assumption worked well in general, but in those outlier cases, it is possible that litter took longer to be incorporated in the soil [54]. Factors such as litter composition affect the incorporation rate. Finally, RothC is a relatively light model in terms of the number of parameters and processes depicted. The processes of SOC accumulation in these cases could possibly require more complex and detailed models, such as the Century model [55], which takes into account more variables. Figure 6a presents the estimated SOC stock as a function of field-measured SOC stock, after application of RothC using the previously estimated parameters estimated to the farms involved in the EDP project. Figure 6a also presents the linear regression fits, where the intercept is significantly different from zero (p-value equal to 2.1 × 10 −6 ), and the slope is significantly different of 1 (p-value equal to 6.9 × 10 −10 ), which means that our model has a non-negligible bias. In the non-biased regression, the intercept would be zero and the slope would be one. In this case, the intercept is about 4 t C/ha and the slope approximately 0.8, meaning that the model overestimates the lower SOC stocks and underestimates the higher SOC stocks. There is an overestimation of the SOC stock of 1.5 t C/ha in the 0-20 t C/ha range, while for higher stocks, there is an overestimation of 2 t C/ha. This is probably due to the fact that about 85% of the measured SOC stocks used for the optimization procedure were in the 10-30 t C/ha range. This is the region where the two linear fits in Figure 6a are closer and the model overestimates SOC stocks by 0.04 t C/ha only.

Estimating AGP and SR for an Unrelated SOC Database
Two additional parameters originally not measured, AGP and SR, were co-determined. No constraints were imposed on either. However, if the model is valid, AGP and SR should increase simultaneously. This is due to the empirical observation that, in general, if pastures are more productive, farmers tend to use them for animal feed to a larger extent [56]. There are some cases, namely when the size of herd is small, when SR does not influences the AGP [57,58]. In the case of SBP, however, the SR, farm and herd size are typically high (e.g., in the AGRO87 project, typically 0.5-1 LU/ha in 14 ha plots), and as such we would expect to observe on average a direct relationship between SR and AGP.
consequently carbon inputs into the soil were low. The only effect of precipitation in these farms was the increase of soil respiration by affecting SOM decomposition rates, particularly in the DPM pool, resulting in low measured SOC accumulation. The most likely explanation for the lower AGP, SR and SOC in EDP farms is that pastures during the AGRO project were accompanied by technical services ensuring best management practices but not before or during the EDP project. For this reason, these SBP areas showed more signs of degradation, likely because of inadequate grazing management, lack of fertilization and/or reseeding when necessary. A clear marker of the degradation is the deficient annual re-establishment of the pasture from the seed bank, as highlighted by the decrease in legume species (which are almost all introduced).    Figure 6b presents the SR estimated with the model as a function of the estimated AGP, for the parameter set with lower loss function value. The list of estimated AGP and SR that minimize the loss function for all farms in the EDP project is in Table 2. As expected, farms and years with higher AGP also have higher SR as a general rule, i.e., when there is more feed available there are more grazing animals. The r 2 is close to one (approximately 0.943) and the RMSE is 1.589. Additionally, there are also strong and significant (at 1% level) Spearman's correlations between estimated SOC (average per farm and year) with AGP and SR, at 0.524 and 0.547, respectively.
The inverse of the regression slope can be interpreted as the animal pasture animal intake per LU, 3.77 t DM/LU. This intake can be converted into the crude protein (CP) and gross energy intake (GE) provided to cattle per head, and compared with the necessities of the animals.
We assumed, for comparison purposes, that the animals on the farms can be either cows or calves. We considered one cow to be 1 LU, and one calf to be 0.4 LU, and used the pasture CP and GE content obtained by Morais et al. [24] for SBP. Allocating pasture feed according to the proportions in LUs, for an adult cow, the pasture provides 525 kg CP/cow.yr and 55 × 10 3 MJ GE/cow.yr, and for the calf 261 kg CP/calf.yr and 27 × 10 3 MJ GE/calf.yr.
Regarding the necessities of the animals, CP and GE intake are highly dependent on the animal liveweight, breed, type of diet, month since calving, peak milk (for adult cows) and other variables. We therefore compare the CP and GE intakes estimated above with plausible variation ranges in the literature [59][60][61][62]. In these calculations, we considered live weights for adult cows and calves of 677 kg and 320 kg, respectively, for the Alentejana breed [24].
For an adult cow, the need for CP ranges between 215 kg CP/cow.year for a cow with 600 kg and 7 months after calving [60] and 690 kg CP/cow.yr with 720 kg and 2 months after calving [61]. Our results of 525 kg CP/cow.yr are within this range.
For the calf, CP intake required ranges between 213 kg CP/calf.yr for a calf with 300 kg [59] and 530 kg CP/calf.yr for a calf with 350 kg [60]. Our result of 261 kg CP/calf.yr is also within this range. Regarding GE intake for an adult cow, it ranges between 49 × 10 3 MJ GE/cow.yr for a low energy intake diet [62] and 60 × 10 3 MJ GE/cow.yr for a high energy intake diet [63], and we obtained a mid-range result of 55 × 10 3 MJ GE/cow.yr. Table 2. Aboveground productivity and stocking rate obtained in the model application for EDP farms. DM-Dry matter, LU-Livestock unit. The year 2007-2008 showed generalized higher SOC stock increases for all farms. According to data measurements, the average annual increase was 3 t C/ha. The model estimates the increase at 1 t C/ha. Nevertheless, this is not fully reflected in the AGP and SR estimated. The year with highest AGP and SR is 2010-2011, although the AGP difference is only 0.1 t DM/ha between 2007-2008 and 2010-2011, and the difference of SRs is only 0.04 LU/ha. Table 2 presents the AGP and SR per farm and production year for farms surveyed in the EDP project. Average AGP and SR are each about one third of the AGP and SR measured in the AGRO projects (EDP: 2.1 t DM/ha and 0.53 LU/ha; AGROs: 6.8 t DM/ha and 0.91 LU/ha). This could happen for several reasons. It could be due to the fact that the pasture plots in these farms are older, as hinted by the lower SOC stocks accumulation (only 0.3 t C/ha while in the AGROs farms it is 2 t C/ha) [64]. Thus, the model would use less carbon inputs from plants and animals. However, older SBP also benefit from the positive effects of accumulated higher SOM, which should improve N availability from mineralization. The reason for the lower AGP could also be the fact that the composition of the pastures in both projects was also significantly different. In AGRO farms, legumes were on average about 40% [38], while EDP pastures were dominated by grasses, and legumes were less than 25% [64]. For example, in farm 10, there was a predominance of grasses more commonly found in natural pastures [64]. Nevertheless, there was also considerable variation of species and varieties among SBP sown during the AGRO project alongside stable yields. Climatic data also influence the results obtained. Monthly average temperature is similar for both datasets (EDP: 18 • C; AGRO: 19 • C), but annual precipitation is significantly higher in the period between 2006 and 2012 (EDP: 700 mm; AGRO: 500 mm). Higher precipitation did not produce increased yields in EDP project farms, and consequently carbon inputs into the soil were low. The only effect of precipitation in these farms was the increase of soil respiration by affecting SOM decomposition rates, particularly in the DPM pool, resulting in low measured SOC accumulation. The most likely explanation for the lower AGP, SR and SOC in EDP farms is that pastures during the AGRO project were accompanied by technical services ensuring best management practices but not before or during the EDP project. For this reason, these SBP areas showed more signs of degradation, likely because of inadequate grazing management, lack of fertilization and/or reseeding when necessary. A clear marker of the degradation is the deficient annual re-establishment of the pasture from the seed bank, as d by the decrease in legume species (which are almost all introduced).
As most prior studies on SBP have focused on particular aspects of the system (such as SOM increases), other examples of co-determination of AGP, SR and SOC is difficult to find. As shown before in this paper, these variables are highly related with the farm, due to different climatic conditions, management and other aspects, such as soil characteristics. Due to the difference between SBP and other pasture systems, especially in terms of AGP and SR [38], the results obtained in this paper for SBP are not directly comparable with other pastures systems. Consequently, we found only one study in the literature [64] that (a) focused on SBP, and (b) included one of the same farms, within the period covered in this paper. In 2010-2011 in farm 11, we estimated that the AGP was 2.2 t DM/ha, which is similar to 2.8 t DM/ha obtained by Simionesei et al. [64] for the same year and farm. Their work used a physically based soil-water dynamics model, i.e., the MOHID model [65], which was calibrated with AGP field measurements.

Main Limitations and Future Perspectives
SBP have been the subject of special attention in the recent past due to their potential for carbon sequestration. For this reason, most monitoring initiatives have focused on SOM only. Measuring SOM helps quantify carbon sequestration potentials, but it is insufficient to explain why or how SOC increases in SBP, and consequently to extrapolate for larger regions or to estimate what would happen in different areas (important for the expansion of this system). Our paper introduces an idea that, to our knowledge, is innovative or at least rarely applied-that we can use computer optimization techniques and process-based modelling to go back to the incomplete datasets and estimate missing information that provides a deeper understanding of what happened during that period. As demonstrated in Section 3.2, the approach presented here is useful to suppress those data gaps. With the model parameterization performed in this work, it is possible to estimate AGP and SR knowing only SOC stock changes. AGP is a crucial variable for livestock production, e.g., if AGP increases, in principle, grazing can be increased and animals can consume less concentrated feeds. AGP is usually measured through destructive methods of biomass estimation and is limited to small area due to the destructive nature, time, expense and labor involved [66][67][68]. Other non-destructive methods have been proposed, e.g., using remote sensing imaging [69,70]. Nevertheless, the influence of tree cover in "Montado" agri-forestry regions is a relevant limitation to the use of remote sensing methods as well as our approach, as high tree cover density also influences AGP [71]. This is a relevant issue in SBP systems, particularly in this sample. For example, in farm 11, the tree density is about 170 trees/ha of Quercus suber [64]. Thus, there is no single indirect estimation method for AGP, and ideally a combination of approaches could be preferable, depending on the specific case.
The presence of legumes in the pastures is particularly relevant for SOC (which is the case of SBP). In nutrient-poor regions, multifunctional pastures tend to produce more in comparison with monoculture pastures, which translates into more available feed [16,72]. This is more relevant because SBP are a mix of grasses and legumes, while other pastures systems usually are mostly composed by grasses. Nevertheless, our work focused only on the pastures systems and not on the species that are in the pasture, i.e., the effect of presence of legumes was not assessed. For example, RS is one parameter highly related with the species, e.g., grasses tend to have higher RS than legumes [73,74]. Thus, the RS of the pasture is highly dependent on the share of grasses/legumes (in AGROs SBP parcels the fraction of legumes varied between 10 and 73%). The fraction of legumes also influences nitrogen intake and consequently nitrogen excretion. Nitrogen content of dung (and therefore the C/N ratio) is one of the most relevant features influencing dung production and decomposition [75][76][77]. Animal excretion usually is addressed in terms of nitrogen content, while our work focuses on carbon, which introduces the issue of the carbon content or C/N ratio of dung. C/N ratios in dung have a considerable large degree of variation (between 15 and 25 [51]) depending on multiple variables, e.g., temperature and moisture [78]. In our work, none of these effects of legumes could be evaluated.
Carbon and nitrogen fluxes in livestock production are one of the major contributors to GHG emissions [79][80][81] and nutrient losses [82,83]. A complete assessment of GHG balance requires a more comprehensive analysis of the production system, including animal dung left on the field and sensitivity of emissions to feed, plant litter formation and also include nitrogen in soil sub-systems [84]. Low nutrient retention in soils due to nutrient leaching is a particular important point in semi-arid regions [85]. Future work should connect nutrient leaching with other aspects of the carbon and nitrogen balance, such as the fraction of legumes or organic matter content [23,72]. Furthermore, so far, mass balances of livestock production analysis have used one hectare as the unit of analysis. A mass unit of product can also be used, as it is the choice of most studies targeting meat production, as in life cycle assessment studies [86][87][88][89], with some exceptions. This mass unit should, however, always be adjusted for nutritional value to encompass the complex nutritional functions provided by food products (such as some units assessed by [90]) and connected them to global impacts of different human diets [91].
As future work, the model should also be calibrated and tested with more recent data and include more sites and longer periods. SOC saturation is believed to be attained after about 10 years [92], but few measurements have been conducted in SBP that were more than 10 years old. Barriers to obtaining longer SOC time series are the difficulty and cost of obtaining data. Merging our work with indirect data collection of these variables-for example, through remote sensing, such as, for example, aerial visible and near-infrared photographs and/or satellite data-can avoid the burden of data collection.
The outcomes of this paper have already been used in the life cycle assessment of the SBP system [24,93]. They were also used to understand C cycling in SBP ecosystems [84]. Another potential application of this work is the calculation of C sequestration factors in SBP. With regionalized AGP and SR, the results of this paper could be used to update the Portuguese C sequestration factor for SBP (6.6 t CO 2 e/ha.yr), obtained by Teixeira et al. [22] (based on the SOM data from the AGRO 71 and 87 datasets also used here), which was used in the Portuguese National Inventory Report [50]. The Teixeira et al. model [21] is a simple SOC dynamic model, considering a single carbon pool with time-invariant carbon input and mineralization rate. The carbon input is a linear function of SOM 0 , the initial SOM level, taken as a proxy for local features (mostly edaphoclimatic) that characterize the ecological quality of each site. This proxy was introduced to capture the pattern in the data that locations with higher initial SOM typically had higher increases in SOM. This model [22] avoids overfitting issues that can occur in optimization approaches, but it was obtained for specific conditions of the farms where data was collected without biogeochemical process base. In our case, site-dependency in the model is provided by AGP and SR, which determine C inputs to soil, and by temperature and precipitation, which are important to determine mineralization.
Therefore, C sequestration according to the model by Teixeira et al. [22] is pre-determined by the initial SOM concentration. C sequestration according to the model calibrated in this paper takes into account more site and management-specific parameters. These parameters would enable more accurate regionalized estimations if regionalized AGP and SR were available, but for the moment it is at least possible to compare the results of the present model with the C sequestration factor previously established under the same conditions. For comparison purposes only, we ran the model using the average initial SOM concentration assumed in Teixeira et al. [22] (0.87%), a soil bulk density of 1.48 g/cm 3 , an AGP of 7 t DM/ha, SR of 1 LU/ha, and the average temperature and precipitation for the region of Alentejo (16 • C and 75 mm). Our model provides a slightly higher estimation of the potential carbon sequestration under these average conditions, 7.2 t CO 2 e/ha.yr. Note, however, that Teixeira et al. [22] considered an initial SOC stock of 7.4 t C/ha (equivalent to the SOM concentration of 0.87% in 10 cm topsoil), while, according to the LUCAS-topsoil database [43], the average SOM concentration in pastures and agroforestry areas in Alentejo region is about 1.7% (almost double of SOM concentration considered in Teixeira et al. [22], also at 10 cm depth, assuming uniform distribution in the 0-20 cm layer). Furthermore, Teixeira et al. [22] implicitly obtained a mineralization rate of 19%, while using this paper results in a significantly lower mineralization rate (only 7%). In the formulation of Teixeira et al. [22], mineralization rate and annual carbon inputs are arrived at statistically using information on SOM only. This interdependency is the justification for the higher mineralization rate, as the model arrived at an adjustment with faster carbon turnover than estimated here.
The results used in this paper therefore have two main advantages compared to Teixeira et al. [22]: (a) they reflect more accurately the dynamics of C accumulation in SBP; and (b) they enable a regionalized calculation of C sequestration based on more site-specific parameters.

Conclusions
This paper is the first application of a process-based soil model to SBP. The procedure established here can successfully reduce the need for data collection. Combining theoretical knowledge (models) with machine learning techniques can successfully help infer additional information from incomplete data sets. The main outcomes of this work are (a) the parameterization of SBP in Portugal, using a soil process-based model (RothC model), and (b) the demonstration of the application potential of the calibration results to estimate AGP and SR in different SOC datasets.
SBP is a unique pasture system that increases SOC accumulation, which is particularly relevant in semi-arid regions strongly depleted of organic soil materials [22,94]. We provide a better characterization of these livestock production systems, looking specifically at the dynamics of SOC as affected by (1) the plant sub-system (RS), which feeds livestock and provides carbon inputs to soils; (2) the soil sub-system (SOC stock changes), which accumulates SOC and emits carbon through mineralization of organic materials; and (3) the animal sub-system (LD and LI), which takes up plant materials and re-deposits them during grazing as animal carbon inputs to soil. The optimization procedure points to a likely RS ratio of 0.98, an input from animal excreta of 1.53 t C/LU, and a proportional animal intake of 0.61. This estimation of animal intake was independently validated for the crude protein and energy content of dry matter. Regarding DPM/RPM, the high value obtained here is characteristic of improved grasslands, where organic matter turnover is quicker. The same model was applied to a different dataset consisting of SOM measurements in older sown pastures, where we estimated AGP (average: 2.0 t DM/ha) and SR (average: 0.52 LU/ha). AGP and SR in this sample are significantly lower than in the dataset used for calibration, likely because these plots were not duly accompanied by technical support systems and therefore had no specific maintenance management operations. In the future, the model and approach established here can be used to improve the carbon sequestration factor of the SBP system in Portugal, as well as to evaluate the likely effects on yield and carbon sequestration of SBP installation in other farms within the Mediterranean region.
Author Contributions: T.G.M. performed the analysis and was the primary writer of the manuscript. N.R.R. collected data and assisted in the analysis of results. R.F.M.T. and T.D. conceived and supervised the study, and assisted in the writing process.