The study site is associated with the River Pang adjacent to the Kennet Catchment. Site selection criteria are mainly linked with a project partner, e.g., Landwise NFM Project. They held the license for data acquisition for the catchment. This research was a part of modeling work packages for exploring NFM strategies for local farm production. Moreover, the site is near the University of Reading, Berkshire, United Kingdom. The University of Reading Weather Observatory provided the weather data used in this study. Access to resources such as data acquisition remained the key consideration for catchment selection for this research. The Pang River flows through the Berkshire Downs to join the River Thames at Pangbourne. The ground water table is the main fed water, and a little water from tributaries drains the cover of tertiary clay and sand in the lower catchment. The Pang has some challenges due to physical modifications, water quality, low flow, flooding and invasive species. Physical modifications are not limited to building structures such as culverts, bridges, tracks, infrastructure, etc., next to watercourses and water bodies. They lead to the loss of natural processes hampering the water environment (Pang River, 2020).

The Environment Agency monitors water bodies and catchment areas for their ecological, physio-chemical quality elements, hydro morphological supporting elements, chemical, priority substances and other pollutants. They classified the ecological status of the Pang as moderate and the water body type as a river, and the Pang water body administered that. The length of the river is 36.57 km^2, and the catchment area is 170.531 km² (Pang Water Body, 2024). There is no declared protected status. However, small sites of specific scientific interest are in the wider catchment, including Sulham, Tidmarsh Woods, Meadows, Coombe Wood, and Frilsham. The Pang Valley has a history of fluvial and groundwater flooding. The risk of flooding from surface water extent is about 0.1 percent annually (About the Pang, 2019). Local farming revolves around farming objectives and their interactions with businesses and lifestyles. The main land cover plus crops include grasses, wheat, oilseed rape, barley, maize, potato, beans, potatoes, etc. Winter wheat is a major cereal crop widely cultivated in arable farming systems in the catchment (Pearce, 1952; Perkin and Rehman, 1994).

This study utilized the Decision Support System for Agrotechnology Transfer (DSSAT), a software application that delivers dynamic crop growth simulation models for different crops. This tool evolved with multiple functions, incorporating weather, soil, crop management, and sample datasets for several crop models. The crop simulation models within DSSAT simulate crop development, growth, and yield based on the interaction of soil–plant-atmosphere dynamics (Jones et al., 2003). Employing three core modules, namely the Sbuild-module (for soil data file input), Weatherman-module (for climate data file input), and the Xbuild-module (for crop management data input), this tool facilitated wide-ranging crop growth modeling research (Decision Support System for Agrotechnology Transfer, n.d.). Hence, this research explored the DSSAT tool, found the most appropriate prototype experimental model for winter wheat and incorporated the available datasets (such as soils, winter wheat crop management, and climatic weather information) from the catchment area.

This research applied sample data features from the local soilscape, which describe the spatial variation and diversity of soils in a landscape. This data includes a range of soil types, textural composition, physical properties, and key features in a particular area. Various soil mapping and classification techniques delineate soilscapes and characterize spatial distribution and variability across landscapes (Farewell et al., 2011). Supplementary Table S1 shows the soilscapes explored for the Pang catchment. Soilscapes range from medium soils, deep clays/ silty, chalk and limestone soils with varying organic carbon contents. Important catchment features are highlighted in the map, shown in Supplementary Figure S1. This study used key features of data samples, which include soil features from soilscape, weather indices from climate, and winter wheat crop phenology from a prototype experiment. Data acquisition for three soils from local soilscape and weather indices of the Pang catchment area were accomplished through the relevant sources (NETMAP soilscape data) by the NFM Landwise Project and the University of Reading, United Kingdom Atmospheric Observatory, respectively. Three soilscape PDPG6, PDPG7 and PDPG8, containing key soil features, characteristics and spatial distribution, are summarized in Supplementary Table S1. These are not the entire soilscapes of the catchment but represent a considerable land in the Pang catchment area. This modeling research used detailed variability of real soils for given zones depth/ layers, soil texture, composition, SOC, bulk density and saturated hydraulic conductivity, lower limit, drained upper limit and saturated water contents are highlighted in the brown colored section of the Supplementary Table S2 for PDPG6, Supplementary Table S3 for PDPG7, and Supplementary Table S4 PDPG8, respectively. The “Weather indices” data file is available in the Supplementary material; the soilscape information is described in Supplementary Table S1, with experimental details in Supplementary Table S5 (Jamagne and King, 2002; Rameshwaram et al., 2023; Crop Cover Data Source, 2024; Geology Data Source, 2024; Land Cover Data Source, 2024; Soilscapes Data Source, 2024; Weather Data Source, 2024).

The study incorporated soils as independent variables, utilizing real soil samples with original organic carbon content from various soilscapes. Classification and characterization details are presented in Supplementary Table S1; Supplementary Figure S1. The study applied the “RB209” soil classification approach for this project. Other independent variables include weather data of seasonal rainfall (mm) and temperature (C^°) from the Reading Atmospheric Observatory (Weather Data Source, 2024), as this research compared no-tillage versus conventional tillage practices in winter wheat crop production with tillage preference as treatment under a frequentist approach (Supplementary material S2) which represented as the main study input variable (independent variable) under the BBN modeling approach.

Four model soils were synthesized based on the three real soils mentioned above, with varying organic carbon contents in topsoil (1, 3, 5, 7%) (Table 1). The variations in organic carbon content within the modeled soils were adjusted in proportion to their existing levels in the real soil relative to the mineral particles. These variations were implemented using the Sbuild-module of the Decision Support System for Agrotechnology Transfer (DSSAT) version 4.8, with each layer adjusted accordingly (Decision Support System for Agrotechnology Transfer, n.d.). Fifteen soil profiles, representing each real soil from a soilscape along with its four modeled counterparts, were generated using the Sbuild-module of DSSAT. Refer to Supplementary Table S1 and Supplementary Tables S2–S4 for more information. Detailed mechanisms for synthesizing modeled soils are described as “Keynotes” under Supplementary Table S4. The Sbuild-module of the DSSAT tool applied the below regression equation on input data of soil texture to estimate the carbon in the stable organic matter based on the relationship developed by po Adiku (Porter et al., 2010; Prout et al., 2021). The Sbuild shell of the DSSAT tool calculated the bulk density of each layer of modeled soil.

The Cropping System Model (CSM) was employed within DSSAT v4.8 and contains various crop models to have simulated data for longer periods (Jones et al., 2003; Li et al., 2015). The crop models use the CERES models for wheat, and various controls & management scenarios are employed within the shell to simulate crop growth. Hence, 15 soils each for No-till and Conventional till (CT) executed in the CSM of DSSAT for having 1,440 simulations over 48 years from 1974 to 2022 and a summary of experiments is described in Table 2. The rationale for conducting simulations over an extended period is not only linked with the evaluation of historical data for improved forecasting of the variable of interest but also encompasses the provision for accessing the model with training and testing datasets. This method helped to find the responses of soils demonstrating local soilscapes containing varied carbon contents to deliver influence on wheat yield, surface runoff and GHG-CO₂ emissions through related variables over 48 years in No-till (NT) and conventional till (CT) management applications using the Cropping System Model (CSM) of the DSSAT (Zhao et al., 2024).

The BBN model structure was synthesized by retrieving evidence reported in the published literature and practical applications available and employed in comparable software. Summary of modules for variables interactions applied in the DSSAT and the list of variables with their relationships and evidence from the literature for BBN model structure are highlighted in Supplementary Tables S6, S7. The model structure for variable relationships was also confirmed through the source manuals available with the software, e.g., the DSSAT and AquaCrop. Key considerations revolved around ascertaining and validating the model assumptions as defined and identified in the relevant sections of this manuscript (Waffa and Benoit, 2015).

This study used the Netica software interface to build a graphical network structure and represented variables as nodes and their relationships as a link or arcs representing their connection. The identified conditional (ideally causal) relationships evident through the published literature determined the direction of the link or arc (Supplementary Table S7). A single relationship between two variables helped develop a simple BBN model structure foundation. This rationale is followed for every single relationship to build a Directed Acyclic Graph (DAG), which delivers a model structure having all interrelated variables connected in a network, as shown in Figure 1 (Ni et al., 2011). Additionally, the study proved the structural morphology of the BBN model, analyzed by the Netica software, which demonstrated key interactions and is shown in Supplementary Table S9. These interactions depicted their probabilistic estimates and performance accuracy, delivering the desired model outcomes. The model was evaluated through key performance indicators outlined for model outputs in the following sections.

This study selected the output (dependent) variables based on the study objectives for assessing tillage preferences toward winter wheat production, as tillage practices are practiced to prepare land for crop establishment through farm management to achieve crop yield. Hence, crop yield suits a major output (dependent) variable (Jug et al., 2011). Surface runoff is another vital variable that can measure flood risk reduction dependent upon tilling choices to facilitate or resist the surface water flow during cropping seasons. For instance, decreased runoff facilitates flood alleviation while tilling wheat if no-till or reduced tilling preference is opted. This approach will cover the land through vegetation and limit the interface between soil cover and farm machinery. Hence, surface runoff serves as another vital output (dependent) variable. Wheat tilling on bare land or cover (such as crop residue or mulches on land) can create a difference in resisting surface runoff during wet winters (DeLaune and Sij, 2012). Last but not least, greenhouse gas emissions for carbon dioxide (GHG-CO₂) are another major output (dependent) variable because tillage preferences are linked with the disturbance of soil aggregates, impacting the decomposition rate of soil organic matter and carbon sequestration, etc. Hence, GHG-CO₂ emissions were considered a vital output (dependent) variable based on tillage preferences (Galic et al., 2019).

To quantify output (dependent) variables in response to input (independent) variables of rainfall, temperature, soils from soilscapes, and tillage systems, this research explored the output variables response accessible through the DSSAT cropping system model (CSM). This study included the most applicable variables with conditional or causal influence among conditional interactions. For instance, the abovementioned parameters are focused as input variables in the DSSAT shell. The study investigated output variables of biomass, yield, precipitation interception, available water space, runoff, carbon-soil organic matter, and GHG-CO₂ emissions. The investigations also explored the literature to confirm such key relationships. For instance, crop biomass accumulation is sensitive to weather indices during reproductive, flowering and grain-filling stages. A research study validated this relationship using the AquaCrop Model, which underwent Bayesian calibration involving meteorological, soil, crop, and management parameters (Zhang et al., 2018, 2019). Hence, these parameters (biomass and weather indices) were shortlisted and categorized in the BBN structure (Table 3).

The BBN model was developed based on a few assumptions, such as weather indices with SOC influencing wheat crop yield (Droste et al., 2020). Similarly, different tillage practices influence available water space and change in SOC contents over time (Sharma et al., 2016). Temperature intensity not only affects crop development & growth rate for biomass and influences GHG-CO₂ emissions (Tu and Li, 2017). Rainfall affects biomass production through water stress in drought conditions and can generate surface runoff where extreme rainfall occurs. Higher organic carbon content levels in soils increase soil fertility to generate higher crop yields, but lower OC levels affect the soil’s water-holding capacity. Soil textural compositions with higher rocks limit available water space in soils to impact surface runoff. Higher biomass intercepts precipitations to reduce surface runoff and deliver higher crop yield (Refer to details in Supplementary Tables S6, S7). Hence, the most pertinent variables were included in the BBN.

The BBN model included 11 nodes (variables), 17 links (relationships), and 2,800 conditional probabilities in the final BBN model structure, as shown in Figure 1. The structure of the BBN model is informed by evidence from published literature regarding variable relationships and their practical applications, as reported and utilized in other software. For example, a few software utilize crop growth simulation methods in hierarchical models, expressing these relationships through mathematical equations in the respective tools (Ahmed et al., 2020). These software functions on their fixed models for executing crop simulations. The DSSAT software uses the Crop Environment Resource Synthesis (CERES-Wheat) model to simulate wheat crop growth to climate, soil, and management, etc. and uses the default module, namely CERES-Godwin for soil organic matter and CENTURY model, to simulate SOC dynamics (Ahmed et al., 2023). A few modules and sub-modules for Actual CO₂, Mauna Loa, Hawaii (Keeling curve) are used for daily records of global atmospheric carbon dioxide concentration and others for soil conservation service for infiltration, modified soil profile for soil layer distribution, etc. (Keeling, 1986). Tillage practices impact soil carbon dynamics by altering soil structure, aeration, and microbial actions, which induce CO₂ emissions. The variable interface for the primary and sub-module in DSSAT is summarized in Supplementary Table S6. Then, the model executes to generate simulations after making an experimental set-up with all modules fed with the required input datasets and files.

In DSSAT, crop models can simulate single cropping and subsequent cropping systems. The soil water balance is based on a formulation by Ritchie methods, applying the concept of a drained upper limit (DUL) and drained lower limit (LL) for available soil water. The water accounting procedure for each layer in the soil profile applies to the DSSAT by Porter and Ritchie methods. Water movement from an upper layer cascades to lower layers, resembling a series of linear reservoirs. The difference between rainfall and runoff calculates the infiltration. The drainage process between layers occurs if soil water in the layer exceeds its water-holding capacity. Root water uptake drives the upward flow due to transpiration and soil evaporation. Available soil water determines the potential root water extraction and the root length density of each layer in the soil (Ritchie et al., 1984; Parton et al., 1994; Ritchie et al., 1998; Suleiman and Ritchie, 2004; Boote et al., 2008; Porter et al., 2010; Asseng et al., 2013; Howe, 2015; Khosravi et al., 2022). The study included a parameter, “available water space,” and calculated that using figures of total available water in soil profile minus field capacity in entire simulations for all seasons (Cassel and Nielsen, 1986). This study also explored the DSSAT and Aqua Crop modeling tools manuals, confirming the underlying mechanisms and the variables’ conditional relations in delivering simulations and model outputs. Some models exhibit identical relationships between variables to give respective model outputs. For instance, the simulation of crop yield relies on biomass. The conditional relationships were identified and confirmed through literature and recorded in Supplementary Table S7. The next step was to execute the DSSAT model.

Data synthesis was achieved through the DSSAT tool that generated numerous simulations executed for 30 individual soils, each for 48 years, and a summary is described in Table 2. The DSSAT tool generated multiple variable responses based on its predefined algorithms for modular and sub-modular functions. However, this research explored only specified variables identified for their potential responses in comparing tillage preference choice in relation to wheat yield, surface runoff, and GHG-CO₂ emissions. The main user interface of the DSSAT tool was used to run the adapted model for simulations through window command selection. An information alert, “Simulations are completed,” popped up shortly with an option to click the next tab to access the analysis if the process is completed successfully. Alternatively, error messages and alerts may appear to navigate the inaccuracy for model calibrations at the appropriate modules, such as Xbuild, Sbuild, Weatherman, etc., and run the program again to complete the process successfully. A flowchart illustrating the process of generating simulated datasets in the DSSAT main interface, resulting in the creation of DSSAT OUTPUT files, is included as an annexure to Supplementary Table S6. This DSSAT data served as input data for the parametrization of the BBN (Jones et al., 2011).

Hence, the DSSAT served as a data simulator only, which provided a range of variable responses to soil attributes, weather, and crop management practices depicted in the Pang catchment area. Only pertinent variables and their responses were accessed following the DSSAT simulation compilations. This phase retrieved the system-generated output files (GHG.OUT, PlantGr2.OUT, PlantGro.OUT, SoilOrg.OUT, SoilWat.OUT, Weather.OUT) using the analysis tab and GBuild module (a plotting tool for data visualization). These output files contained the simulated datasets for all variables in the BBN model. This approach resulted in retrieving a total of 1,440 simulations. The acquired datasets were then formatted into a specified case file to parametrize the BBN.

Moreover, this study utilized the Netica software and its user interface to construct and parametrize the BBN model. The structural construction started from one variable (independent) as a parent node to the second (dependent) as a child and connected them based on their conditional relationship. All parent nodes are connected to their child nodes through an arrow link directing a conditional relationship from the parent to the child nodes by avoiding having a cycle in the network. This structure developed into a directed acyclic graph (DAG), a BBN model. All 11 nodes with 17 links were connected in a BBN and were now ready to parametrize using simulated datasets. The specified case file containing (the DSSAT outputs for specified variables) is used to parametrize the BBN model. The study applied the “Holdout Method” for model training and testing for validation. Therefore, seventy-five (75%) data were used for training the model by parametrization, as exemplified in the Supplementary Table S11. The assessment used this to populate the BBN model against unseen holdout datasets for the model, and the remaining twenty-five (25%) were then used to validate the BBN model by testing. The case files for the Holdout method’s application through 25% of datasets are available as specimens in Supplementary Table S12. The study also compared model validation results against error rates among “K-Fold Cross Validation” and “Train Test” Methods (Yadav and Shukla, 2016). The summary is highlighted in Table 4.

The parametric details of the BBN model depend upon the relationships between parent and child nodes and their states defined for each variable in the network and their corresponding datasets (Supplementary Table S8). This aspect depends upon the nature and quality of datasets used to parametrize the BBN model. This study used synthetic datasets generated by the DSSAT tool, which delivered variable output responses following specified functions of certain phenomena, such as the pedotransfer function for soil water dynamics and CERES-Godwin and CENTURY-Parton SOC simulation. The probabilistic estimates for quantified variables in the BBN model were well-aligned with the known and established responses for the assessed parameters. The results highlight the adaptation of conservational tillage can lead to sustainable crop yield, reduced runoff, and GHG-CO₂ emissions by offering precipitation interception and available water space in soils (Rahman et al., 2021).

The output variables of the BBN model are wheat crop yield, surface runoff and total GHG-CO₂ emissions. The response of output variables is the reflection of model input variables and their vital interactions with associated variables based on their identified conditional phenomena in a hierarchical order, as shown in Figure 1. Results found model output variables have exhibited seasonal variations in their responses through long-term simulations. Weather indices have an imminent influence over crop yield, surface runoff, and total GHG- CO₂ emissions. Tillage preference for no-till or conventional methods has a comparable influence over SOC levels in the soils cultivated over extended periods. Tillage practices have profound effects on carbon sequestration and preserving soil structure. This investigation evaluated the model performance and validated that using test cases. The study also analyzed posterior probability distributions of interacting variables to assess their comparative impacts on each other.

Various metrics are applied to evaluate the performance of the BBN model to measure uncertainty (Table 5). These include model complexity (e.g., the number of nodes, links, etc.) (e.g., scenario evaluations) and sensitivity analysis (e.g., sensitivity to findings, case file simulation, etc.) using Netica software (Marcot, 2012; Namdari and Li, 2019).

The BBN model is validated using simulated datasets due to the scarcity of real data availability. Access to real-world data sources, resource constraints, and complications linked to data protection are common challenges in most modeling works. Establishing a long-term tillage-based experiment in competing time frames with challenges that arose due to COVID-19 pandemic implications has prompted this shift toward assessing data through simulation methods. This approach brings some limitations while using the simulated dataset for model validation. These may not capture the adaptation in real-world situations, face result validation concerns, have limited capacity in capturing factual variability, and carry bias embedded in algorithms defined in tools used for data simulations. However, simulated datasets offer a great opportunity to explore complex scientific systems where access to real-world data is a valid concern. This approach presents probabilistic solutions and diverse scenarios for comparing consequences, which could support decision-making for the right choices.

The assessment employed 25% of the synthetic dataset obtained from the DSSAT-generated simulations to test the BBN model for performing model validation in this research (Morgan and Granger and Max Henrion, 1990; Netica, 2024). The test cases representing datasets are used for response variables for 2012–22, the specimen listed in Supplementary Figure S7. This approach is also known as the holdout method for model validation (Jones et al., 2011). Confusion matrices GHG-CO₂ Emissions, Runoff, and Crop yield are tabulated in 6, 7, and 8, respectively. Error rate illustrated two types of error by the BBN model, e.g., type 1 error for false positive, which means rejecting a true hypothesis, and type 2 error means failing to reject a false hypothesis. The sum of these is reflected in the confusion matrix.

The study utilized Netica software to evaluate the model’s performance, which applied scoring rules based on its system-generated computations. The scoring rules are the model evaluation metrics and assess the model’s performance of probabilistic estimates or predictive distributions. These rules measure the predicted distribution, the observed outcome, and the model’s accuracy and calibration. Netica software retrieved the report with a confusion matrix carrying the scoring rule outcomes of the individual output variable. Scoring rule outcomes highlighted the actual beliefs of the states that followed the values in the case file. Results found three types of scoring rules (Logarithmic loss, Quadratic loss, and Spherical loss) performed by Netica software (Pearl, 1978; Netica, 2023).

Logarithmic loss value was calculated based on a natural log and ranged from zero to infinity, where zero represents the best performance at an acceptable level in all the above cases. For instance, the Logarithmic loss value is depicted as 0.9555 from the confusion matrix (for GHG-CO₂ emissions) under Table 6. This value is closer to zero than infinity and declares an acceptable performance level. Logarithmic loss penalizes incorrect predictions more seriously when they are confidently wrong. Hence, the Logarithmic loss value was calculated based on a natural log and ranges from zero to infinity, where zero represents the best performance and having the value of 0.9555 is an acceptable level of model performance in the given example.

Quadratic loss provides another metric to evaluate the performance of a model and is also known as the Brier score and ranges from zero to 2, where zero is the best. Quadratic loss penalizes larger errors more severely than smaller errors. Therefore, this value becomes sensitive to outliers, as larger errors contribute disproportionately to the overall loss. The quadratic loss value was depicted as 0.4859 from the confusion matrix (for GHG-CO₂ emissions) under Table 6. This value is closer to zero and sets an acceptable level for the model performance. Spherical payoff is another matrix to measure the performance of the BBN model by considering the accuracy and calibration of predictions. This scoring value is linked to the balanced risk–reward profile of a sphere. The spherical payoff scoring rule aims to reward the model for achieving high accuracy through making correct predictions following a well-calibrated model for accurate estimates of the likelihood of those predictions. This matrix measures the model performance using an index ranging [0,1], where 1 is the highest value for illustrating better model performance. In the above example, the Spherical loss value was depicted as 0.7122, which is closer to one and considered an acceptable level for the model performance. Similarly, this BBN model performed better than average for all output variables in dealing with the nuances of probability values under scoring rule results as depicted through the Confusion matrix (for Runoff) and Confusion matrix (for Crop yield) and shown under Tables 7, 8, respectively (Marcot, 2012; Netica, 2023).

However, the model was also validated using the K-Fold cross-validation method and tested through the train split approach. The comparative results are summarized in Table 4. The details are provided in the Supplementary material S1. Moreover, the multiple regression model was employed in this analysis and findings are detailed as Supplementary material S2. Such comparisons evaluate the modeling performance among different models, such as multiple regression and network models and promote the use of the business value approach (Bansal et al., 1993).

The BBN model was developed to integrate prior knowledge using simulated datasets. It highlights analyzing model inferences through probabilistic estimates represented by the modeled parameters. Salient features are described below. The BBN model captured uncertainties linked to seasonal variations by providing posterior probability data distributions reflecting probabilistic relationships across all interacting variables in the model and offering decision choices for competing tillage preferences. This study analyzed the posterior probability distribution of intermediary and output variables in the BBN model. This evaluation technique revealed the performance of the best and worst-performing variables individually and on average (Solheim, 2021).

This BBN model assessed tilling systems under climate change and sustainability challenges for winter wheat crop production. The assessment compared 48 cropping seasons of simulations for each included variable response and assessed their long-term impacts on wheat yield, surface runoff, and GHG-CO2 emissions as output variables (Nielsen and Vigil, 2010; Krauss et al., 2017; Nath et al., 2017). This model gained insight into those parameter responses where empirical data is not readily available except for soil, weather, wheat management and yield data. However, the synthetic datasets were generated using established machine-learning models in DSSAT. Then, the BBN model utilized a crop growth hierarchal model DSSAT output data from its prototype model experiment. The BBN model offers integration of identified variables using standardized responses to capture uncertainties among model parameters. These responses reflect conditional dependencies and are consistent with the established knowledge in the domain.

It is important to understand BBN, which provides a graphical probabilistic model representing a set of variables and their probabilistic dependencies using a directed acyclic graph. It modeled uncertain relationships between variables and made probabilistic inferences based on observed evidence. Contrarily, an emulator is a simplified mathematical model that approximates the behavior of a more complex, computationally expensive model. Emulators are typically used to speed up the simulation process by quickly estimating model outputs based on a subset of input values and mimicking the source (Aakula et al., 2024). However, the DSSAT is a mechanistic crop growth hierarchical model and serves as a simulator in this study. This tool operated based on its embedded machine learning capabilities and provided the DSSAT output data used as an input in this BBN. This BBN model presents probabilistic inferences and reasoning, serves different purposes, operates on different assumptions. However, this approach may have the potential to emulate functions not evaluated in this study (Liu and West, 2009). BBN models could have other challenges in building and maintaining large complex networks, and computational expansive, especially data-driven networks are highly dependent on the amount and quality of data. The BBN models can not represent feedback loops or cyclic dependencies between variables which could challenge applicability to certain domains. In such cases, alternative probabilistic models such as hidden Markov models or Markov random fields might be a better substitute. Nonetheless, this BBN model did not come across the challenges above (Wang, 2004).

The BBN model is appraised per established metrics and summarized in Table 5. No contradictory responses were observed from model performance except for a higher error rate in predictive accuracy highlighted in the confusion matrix. The error rate is higher for GHG, runoff, and yield and is reported as 29.39, 33.64, and 43.94%, respectively. This outcome could be due to a greater variation in the organic carbon contents defined with varying fluctuations among 30 soils (Tibshirani and Tibshirani, 2009). Three real soils with built in soil organic level contents have been used to simulate each into four varying OC levels. Hence, the variation in GHG, runoff, and yield must appear in the results and are also subject to a higher error rate. However, the model’s predictive accuracy is still over 56% among all output variables of interest using the holdout method of model validation. The model was trained using data from 1975 to 2011 and tested against 2012–2022. However, the model was also validated using the K-Fold Cross Validation (CV) method using K-fold splits into training and test datasets. Then, the model was trained using data from 1975 to 2014 (K1 + K2 + K3 + K4) and tested using data from 2015 to 2022 (K5) to yield consistent and reliable results, proved (Table 4; Supplementary material S1). This challenge commonly appears among models in several agroecological studies, and an example is presented in Figure 5 by comparing published crop yield datasets (2016–2021) from research covering the same catchment area (Hazarika et al., 2009; Karimi et al., 2021; Mahmood et al., 2023). Model performance could be improved by employing parameter tuning, varying (scaling) parametric categories, changing (normalization) datasets, data homogeneousness, etc (Singh and Singh, 2022). Further exploration of the DSSAT tool for its interacting functionalities could be another area of interest for comparing model parameters’ performance against those of other competing tools.

The BBN model illustrated that the No-till (NT) choice of cultivation outclassed the Conventional till (CT) by reducing GHG-CO2 emissions probability estimation in the highest band from “7.34 to 7.31%.” Probabilistic estimates are expressed in percentages (e.g., 0 to 100%), representing the likelihood or probability of different outcomes based on the available evidence and model assumptions. They can be presented as ranges or uncertainty intervals, with one bound (lower or upper) specified (Hohle and Teigen, 2018). NT also exhibited a probability of a rising trend in higher wheat yield bands, increasing from “7.46 to 7.56%.” Cumulative runoff probability decreased in higher bands, with the “high band” diminishing from “8.52 to 8.50%.”

In contrast, the CT option had a snowballing trend in GHG-CO₂ emissions in the highest band from “7.34 to 7.36%.” CT resulted in probabilistic estimates of lessening in wheat yield probability in the highest band from “7.46 to 7.35%,” while cumulative runoff amplified in higher bands, with a probability of the “high band” rising from “8.52 to 8.55%” (Radford, 2007; Lu et al., 2015). The impact of tillage methods on GHG-CO₂ emission is linked with the potential phenomenon of creating soil disturbance in adapting tilling frequency and intensity, which can trigger disturbing soil aggregates, impacting water retention and moisture availability for plant uptake and influencing root development. These methods also influence soil microbial communities and their activities in decomposing soil organic matter. Nutrient cycling is another aspect affected by tilling choices, which can favor or resist underlying mechanisms. Tilling preference influences SOC contents. Higher levels of OC contents favor maintaining soil health conditions for building better soil structure and more resilient soils to help achieve better crop production.

The BBN model effectively apprehended uncertainties linked to seasonal variations by providing posterior probability data distributions reflecting probabilistic relationships across all variables and offered decision choice for NT favoring soil carbon stocks (highest among soil “NT.OC-7%PDPG8,” e.g., 286,634 kg/ha) in winter wheat. Furthermore, the model favors reduced GHG-CO₂ emissions through NT practices and naturally restraining surface runoff generations for flood mitigation (Maraseni and Cockfield, 2011). For additional improvement of the BBN model, examination of early and late seasonal planting, supply rates of nutrients, and the use of various crop varieties under diverse farming systems are to be explored. A few constraints are not limited to access to empirical and synthetic datasets used in this model. Simulated datasets are used to calibrate and validate the BBN model (Aliferis and Cooper, 1994; Sebastiani and Ramoni, 2001; Newlands et al., 2014; Ali, 2023), so results variation could be expected where entire real/ empirical datasets will be applied.

This study focuses only on limited output variables, e.g., GHG-CO₂ emissions and does not cater to nitrous (N₂O) and methane (CH₄) emissions. Surface runoff measurement presents only standard soil-specific features in the S-module of the DSSAT and could consider slope and topographic features for future exploration. The study utilized simulated wheat yield and compared that with real field data from the catchment area for 2016–2021. Utilizing catchment-wise yield data could enhance yield estimations on a larger scale. While simulated datasets were generated solely using the DSSAT tool, exploring other competing tools could yield comparable results. These areas offer futuristic research on this model. Moreover, variation in model response could also be expected due to spatio-temporal factors and the extent of noisy datasets once employed. Hence, perfect datasets can deliver more reliable model outcomes (Merino et al., 2016).

This research employed simulated datasets, offering valuable research, experimentation, and decision-making prospects. However, using simulated datasets entails implications regarding predictive accuracy, cost-effectiveness, experimental flexibility, risk mitigation, and ethical considerations. Researchers must carefully consider these factors when utilizing simulated data in their studies. Seeking real-world datasets whenever possible allows for comparing predicted results with observed data, mitigating some of these concerns. This study compared wheat yield datasets for comparing against simulated datasets and found few soils exhibited comparable results.

Overall, valuable improvements were made to the BBN, which would not have been possible without data access. This was particularly relevant in this study, where dataset access for a longer period was nonexistent. This approach helps integrate DSSAT simulated data into a BBN, provides an opportunity to measure against actual wheat yield (Figure 5) and provides a good insight to compare the model and identify the higher-performing soils under competing tilling methods. Hence, the model can potentially improve the accuracy, realism, and robustness of the model’s predictions and enable more informed decision-making in agricultural systems.

The current iteration of the BBN is in its alpha stage and recommends testing in diverse real-world scenarios to progress into beta versions. Successful scalability may lead to the development of a commercial gamma-level BBN model, customizable for diverse applications. Furthermore, enrichments in decision support functionalities, informed by real-life data testing, could enable the introduction of individual predictive models as independent software tools for mobile apps and portal-based interfaces.

The following are evident advantages achieved by this BBN model.

The model can be used to inform decision-making and improve agricultural sustainability. For instance, adapting suitable tillage preferences is a crucial management choice with a commercial focus on obtaining crop yield. Conventional tillage choice could deliver instant yield gains at the cost of hampering long-term soil health and structural resilience. This approach will cater only to the interest of acquiring optimum crop yields at the cost of attracting several complicated issues of soil compaction erosion, consequently attracting higher surface runoff & greater flooding risk in the long run. For example, CT promotes frequent and intense tillage application, exerting mechanical pressure through interfacing soils. Heavy and frequent tillage can expose the soil organic matter to increased oxygen levels and microbial decomposition, affecting soil functioning for crop growth and triggering GHG-CO₂ emissions. Tillage methods can also influence soil nutrient cycling processes by disturbing soil physical, chemical, and biological properties. Conventional tillage practices could accelerate the mineralization of organic matter, releasing nutrients and raising nutrient availability for plant uptake with an increased risk of loss of nutrients. These changes in nutrient availability affect crop growth, yield and nutrient use efficiency, which could influence soil carbon dynamics and GHG-CO₂ emissions.

However, the reduced-to-no-till choice of cultivation can bring multifaceted benefits along with attaining sustainable crop yields over time, limiting the loss of SOC, preserving soil health conditions, limiting surface runoff due to land cover vegetation preference favoring flood alleviation, and reducing the risk of soil compaction and erosion (Farahani et al., 2022). At the helm of climate change and sustainability challenges, the importance of carbon sequestration has become manifold. This model supports reduced to no-till (NT) as the right choice of tillage practices through which crop residues and re can be retained on the land to get incorporated into soils, supporting land cover and increasing SOC contributions. This aligns with sustainable crop farming practices.