BODIUM—A systemic approach to model the dynamics of soil functions

The increasing demand for biomass for food, animal feed, fibre and bioenergy requires optimization of soil productivity, while at the same time, protecting other soil functions such as nutrient cycling and buffering, carbon storage, habitat for biological activity and water filter and storage. Therefore, one of the main challenges for sustainable agriculture is to produce high yields while maintaining all the other soil functions. Mechanistic simulation models are an essential tool to fully understand and predict the complex interactions between physical, biological and chemical processes of soils that generate those functions. We developed a soil model to simulate the impact of various agricultural management options and climate change on soil functions by integrating the relevant processes mechanistically and in a systemic way. As a special feature, we include the dynamics of soil structure induced by tillage and biological activity, which is especially relevant in arable soils. The model operates on a 1D soil profile consisting of a number of discrete layers with dynamic thickness. We demonstrate the model performance by simulating crop growth, root growth, nutrient and water uptake, nitrogen cycling, soil organic matter turnover, microbial activity, water distribution and soil structure dynamics in a long‐term field experiment including different crops and different types and levels of fertilization. The model is able to capture essential features that are measured regularly including crop yield, soil organic carbon, and soil nitrogen. In this way, the plausibility of the implemented processes and their interactions is confirmed. Furthermore, we present the results of explorative simulations comparing scenarios with and without tillage events to analyse the effect of soil structure on soil functions. Since the model is process‐based, we are confident that the model can also be used to predict quantities that have not been measured or to estimate the effect of management measures and climate states not yet been observed. The model thus has the potential to predict the site‐specific impact of management decisions on soil functions, which is of great importance for the development of a sustainable agriculture that is currently also on the agenda of the ‘Green Deal’ at the European level.

over, microbial activity, water distribution and soil structure dynamics in a long-term field experiment including different crops and different types and levels of fertilization.The model is able to capture essential features that are measured regularly including crop yield, soil organic carbon, and soil nitrogen.In this way, the plausibility of the implemented processes and their interactions is confirmed.Furthermore, we present the results of explorative simulations comparing scenarios with and without tillage events to analyse the effect of soil structure on soil functions.Since the model is process-based, we are confident that the model can also be used to predict quantities that have not been measured or to estimate the effect of management measures and climate states not yet been observed.The model thus has the potential to predict the site-specific impact of management decisions on soil functions, which is of great importance for the development of a sustainable agriculture that is currently also on the agenda of the 'Green Deal' at the European level.

K E Y W O R D S
agriculture, computational model, simulation, soil microbiology, soil structure, sustainable soil

| INTRODUCTION
There is an increasing awareness within the political and public arena that soil plays a critical role for the functioning of the Earth system (Keesstra et al., 2016).At the European level, the new Green Deal formulates ambitious objectives with respect to more sustainable soil management (Montanarella & Panagos, 2021).
Agriculture is our essential basis for the production of food and many raw materials.At the same time, agriculture is the activity where man intervenes most strongly in the natural functions of soil.This is done, for example, by selecting crop rotations, tillage, and supply of organic and mineral fertilizers, as well as pesticides.In the first place, the aim is to achieve higher and more secure yields.However, in addition to the production function, there are other soil functions that we know are essential for the functioning of our ecosystems: buffering water flows, being a filter for clean groundwater, recycling nutrients, storing carbon and being the habitat for a myriad of organisms.For this reason, we need a sustainable or ecological intensification of agriculture that takes into account the protection of the other soil functions.Addressing this multifunctionality requires the development of appropriate policy frameworks and practical decision-support tools.A prerequisite for this development is the ability to quantify the impact of agricultural management on all these functions.Providing the necessary tools for this is a crucial task for soil scientists (Bouma, 2020).
The various soil functions can be considered to emerge from interactions of biological, physical and chemical soil processes which are closely linked.The typical approach to evaluate the state of soil functions today is to find suitable indicators that can be observed and that are clearly related to these functions (Bünemann et al., 2018).In this way, the actual state of soil functions can be evaluated and, based on time series, also their temporal dynamics.However, this approach is unsatisfactory in the sense that it is not possible to directly predict the effect of individual management measures on the various soil functions and to unravel the processes that are responsible for the change in soil functions.This can be achieved by a modelling approach representing the relevant soil processes and their interactions in a realistic way for given site conditions so that the soil functions emerge from these processes.Hence, to be site-specific and predictive the required model tools need to be process-based and mechanistic, and should represent soil as a complex system (H.-J.Vogel et al., 2018).There is a number of models following such a systemic approach but focusing on different aspects of the soil system such as dynamics of soil organic matter (Coleman & Jenkinson, 1996;Franko et al., 1995;Parton, 1996), nitrogen (Engel & Priesack, 1993;Godwin & Jones, 1991;Haas et al., 2013), water (Belmans et al., 1983;T. Vogel et al., 1996), solute transport (Vanclooster et al., 1994), biological activity (Deckmyn et al., 2020;Komarov et al., 2017) or integrating different existing model approaches within one framework (Holzworth et al., 2014).According to the respective focus, these models are missing some processes, which are, however, essential for the representation of the overall system.For arable soils, a typical example of such a missing component is soil structure formation by tillage.Some model approaches already incorporate several direct effects of tillage operations on different aspects of the soil system, such as turnover rates of soil organic matter (Jordon & Smith, 2022), mixing of soil organic matter and nutrients (Holzworth et al., 2014), surface residues, or a change in bulk density affecting hydrological processes (Andales et al., 2000).However, they do not account for the dynamic long-term effects following the explicit change These missing components can have implications for many other processes such as water (Fatichi et al., 2020) or carbon dynamics (H.-J.Vogel et al., 2022).Some other model concepts already account for soil structure dynamics induced by biological activity following different approaches on describing soil structure and with a focus on more specific soil functions such as carbon storage (Meurer et al., 2020) while BODIUM targets all soil functions and their long-term development in response to agricultural management.A mechanistic modelling approach that operates at a 'whole system' level, that is, including all necessary processes and providing predictive power with respect to soil functions has been proposed by H.-J. Vogel et al. (2018).The BODIUM model presented in this paper is motivated by this idea.
The intention of BODIUM is not to predict the yield or nutrient requirements for the current crop as is the case for typical soil-plant models in agriculture.The main intention is to evaluate different farming systems, characterized by crop rotations, tillage and fertilization strategies, in their effect on the different soil functions in the medium and long term.This is to be achieved by correctly modelling the relevant processes and their interactions according to our actual state of knowledge.In this way, we want to make sure that the differences between management measures in terms of their effect on soil functions are qualitatively represented correctly by the model.Hence, the aim is not to predict, for example, the yield or the nitrogen losses in quantitative terms (i.e., in kg/ha) but to evaluate the differences between the management option (i.e., is there a difference and if so, in what direction and how big is it approximately).This leads to another special feature of the BODIUM model, namely that the model parameters represent our current understanding of the underlying processes and are not calibrated based on any measured observations.To our knowledge, BODIUM is the first systemic model that follows this concept and can thus become a valuable decision-support tool in agriculture.BODIUM is intended to be a reliable process-based model to evaluate the long-term effect of future scenarios for which empirical data are missing and for which no experimental approach is available due to the long-term nature of the problem.The model is and will be in development as our understanding of the soil processes continues to improve.
The structure of this article is as follows: In the first section, we introduce the general model concept in terms of spatial and temporal scales, the type of processes included and the level of detail of their descriptions.This is followed by a short description of the technical implementation and the programming framework.Then we describe the long-term field experiment which is used to demonstrate the model performance.The model's plausibility is confirmed by comparing simulation results with experimental data for biomass production and the dynamics of carbon, nitrogen and microbial biomass, and analyse how new features of our model affect the outcome.Finally, the impact of tillage on soil functions is simulated as an example of explorative modelling, which allows the simulation of any management or climate scenario to predict future dynamics.The details of the implemented processes and their parametrization are provided in an extended supplement.
The BODIUM model is designed to explicitly represent all relevant components of the complex soil system.These components are interconnected by a large variety of processes.Many of these processes are well known and described by appropriate equations, and some of them are rather hypotheses that can be tested using a simulation model.The ambition of the model is to connect all known and presumed processes in a structurally correct way according to the current state of knowledge.The expectation is that this mechanistic approach will be in a position to reproduce and predict the effect of external drivers such as soil management and climate on the various soil functions.

| Spatial scale
Soil functions are highly sensitive to the type of soil and the local site conditions.Especially soil texture is directly related to soil structure, water capacity and the buffer capacity for nutrients.From this, it is obvious that different soil types have to be distinguished.Hence, the characteristic scale of our model is the field scale which is also the typical scale for soil management measures.The parametrization of many processes depends on inherent soil properties such as texture, mineralogy, stone content, depth of the soil profile and the more general hydrologic situation (e.g., depth of groundwater table or waterlogging).This can be accounted for at the field scale such that the model application is highly site-specific and accounts for the fact that agricultural soil management may have very different impacts depending on the local site conditions.
The field scale implies that the soil properties do not vary substantially in the horizontal plane.At this scale, flow and transport processes are mainly driven by vertical gradients in water potential or concentrations of substances.Lateral movement in the range of metres and beyond is limited to situations close to water saturation (H.-J.Vogel, 2019).Hence, the model operates, in fact, at the pedon scale, and flow and transport processes are modelled in 1D along the vertical soil profile with a flexible number of vertically connected sub-volumes (in the following referred to as 'nodes') and fluxes of water, solutes, energy and particulate matter in between.The spatial resolution in the vertical direction is high enough to represent the soil profile (i.e., horizons).While water and matter fluxes can be modelled for unsaturated soils, lateral fluxes close to water saturation and erosion are not represented in the present 1D model version.Simulating a field or watershed with substantial underlying heterogeneity is possible by modelling several profiles and summing up the resulting functions and fluxes according to the soils' proportions.

| Process descriptions
We do not intend to 'reinvent the wheel'.BODIUM builds upon process descriptions that have been successfully applied in other modelling approaches that were typically designed for more specific applications such as crop growth (Groot, 1987;Penning de Vries, 1989;Priesack, 2006;Rahn et al., 2010;Wang & Engel, 1998, 2000;Yin & van Laar, 2005), root water uptake (Couvreur et al., 2012(Couvreur et al., , 2014)), or nitrogen leaching.Those are integrated with other processes where we developed new concepts, for example, for water dynamics, carbon turnover and microbial dynamics.This available knowledge is implemented into a homogenized programming framework that allows for a seamless coupling of biological, physical and chemical processes to reproduce the required system behaviour.
One of the major challenges is to represent all relevant processes at the scale defined by the discrete nodes of the 1D model which is in the range of centimetres to decimetres.For example, the degradation of organic molecules by organisms acting at the micro-scale needs to be described by variables that are relevant at the centimetre scale such as quality of organic matter, water content, temperature or soil structural attributes.This upscaling is a common challenge to all models operating at the scale of pedons since it is simply not possible to model the processes in full 3D at the pore scale.The main problem in the first place is not the lack of computing power, but the impossibility to describe the state of the heterogeneous soil system with the required resolution.Our general approach for upscaling is to start from the available process knowledge at the small scale, where the processes are actually happening, and, based on that, to develop a simplified effective description at the larger scale.An example is the switch from nitrification to denitrification.At the local spot this depends on the local availability of oxygen and organic substrate together with the local microbial activity.This has a sharp transition with the saturation of oxygen.At the large scale of soil pedons, due to small-scale soil heterogeneity, we typically observe a gradual increase in N 2 O fluxes as the soil moisture increases towards saturation.For making the scale change and to gain such a transition for BODIUM, the dynamic air-filled porosity and the derived distance distribution from any point within the soil matrix to the air-filled pore volume can be estimated as an upscaled attribute of soil structure.In this way, the anaerobic soil volume, that is, the fraction of denitrification, is calculated as a function of water content.The spatial distribution of organic matter is assumed to be homogeneous and the overall microbial activity depends on additional boundary conditions such as water content and temperature.In this way, the observed phenomena at the large scale is reproduced by translating the small-scale understanding (i.e., the microbial activity affected by boundary conditions) to an effective description based on available, measurable, macroscopic system variables (i.e., the mentioned boundary conditions driven by soil structure).Soil structure images for such a transformation are available via the soil structure library (Weller et al., 2022).

| Parametrization
A large number of represented processes leads to an even larger number of parameters for their description.Thus, it is not an option to calibrate parameters based on a limited number of observations in a highly under-determined system.In contrast, the required parameters are estimated based on a systematic assessment of the experimental evidence published in the relevant literature.Based on that, the model is operated in the forward mode (i.e., model parameters are not calibrated using measured data for each respective site) to evaluate the impact of various soil management strategies on the different soil functions in a sitespecific way.For this purpose, the interactions of the processes should be represented following the actual state of knowledge and/or well-defined hypotheses and the related parameters should be in a realistic range.The model is not primarily designed to predict a specific soil state variable of a given field.This would require also the management history of this specific field which is often unknown.
The performance of the model is validated based on data collected in a long-term field experiment, where different management strategies are known for a long time span of at least 20 years.The reliability of the model predictions can be assessed by comparison with the measured data observed in these experiments.This serves also as a plausibility check for assumptions that are not yet tested in experiments or values which cannot easily be measured, such as soil structure dynamics.

| Model implementation
In this section, the implementation of components and processes is described at a more general level.For a full description of the processes including equations and a full list of parameters, see Data S1 (Tables S1-S3).

| Model components
The 1D model operates at the scale of a soil profile represented by a vertical sequence of spatially discrete layers, that is, nodes.Each node contains all essential entities depicted as components: roots, mineral nutrients, water and air, fauna (only earthworms are represented explicitly, no fauna was simulated in the present study), microorganisms, soil structural attributes and soil organic matter of different quality (SOM) as illustrated in Figure 1.All these components co-exist and interact within the volume of each node without a spatial explicit location.Additionally, the model includes a stand-alone crop growth component connected to the uppermost soil node and interacting with the deeper nodes via the root component (Figure 1).The number and thickness of the nodes are flexible and can be chosen based on the structure of the soil profile with its characteristic soil layers.This way, both can be adapted to represent flow and transport processes along the soil profile in a realistic way.This is especially true close to the soil surface where the vertical thickness of a node should not exceed a few centimetres while further below a maximum node thickness of about 20 cm is adequate in most cases.There is no maximum depth, however, it is typically assumed that the profile discretization covers at least the maximum rooting depth.The volume and mass of each node may change during simulations to allow for local changes in bulk density as it is typical through tillage or swelling and shrinking.
F I G U R E 1 Scheme of the BODIUM model depicting components (coloured boxes), processes and interactions (arrows) and external factors (grey circles) within each soil node.On the right, the soil profile consisting of various soil nodes is presented.Note that this illustrates the current BODIUM version described in this publication, some of the missing arrows or components will be part of future versions (e.g., other nutrients).

| Processes
The various model components within each node are linked by physical, chemical and biological processes illustrated through arrows in Figure 1.In the following, these processes are described only briefly.
Soil structure is considered to change not only after tillage operations but also during the consolidation thereafter.Thereby, soil structure is quantified through the volume fraction of different pore-size classes including biopores (i.e., channels formed by roots and earthworms) as a special class of pores with increased vertical connectivity.
Water flow is simulated based on the actual water saturation of the different pore size classes.These saturations define the local water potential and hydraulic conductivity.As a special feature during infiltration, water can invade larger pores prior to smaller ones to represent the effects of hydraulic non-equilibrium, hysteresis and preferential flow (Hassanizadeh & Gray, 1979).The transport of dissolved nutrients is directly coupled to water flow.The movement of gas components, especially oxygen and heat, is described by a diffusion process based on local gradients (Moldrup et al., 2004;Peng et al., 2016).
The turnover of organic matter is based on the separation of different pools, namely fresh organic matter, degraded organic matter (assimilable and non-assimilable) and physically protected organic matter bound to mineral surfaces or occluded in the mineral soil matrix.This is similar to most of the currently used carbon models (Abrahamsen & Hansen, 2000;Komarov et al., 2017;Parton, 1996).In addition to that, we consider the active microbial biomass as an individual component that changes its mass depending on the quality and quantity of available resources.Provided the availability of water and oxygen, the interaction between living microbes and available organic matter is mainly governed by stoichiometric relations in terms of carbon and nitrogen (Kirschbaum, 1995;König et al., 2017;Monod, 1949;Sinsabaugh et al., 2013;Stolpovsky et al., 2011;Svendsen et al., 1995;Wang & Post, 2012).
The simulation of plant growth is based on simplified modelling of photosynthesis and assimilate distribution including the plant-specific separation between aboveground biomass production and root growth.Water and nutrient uptake by plants is directly coupled to the actual photosynthesis.
The typical time step of the BODIUM model is 1 day, but can be adapted by the user.If external factors affect processes in a way requiring an increase in temporal resolution, the time step is automatically reduced for specific processes (i.e., water distribution after heavy rainfall).The duration of the simulation is also defined by the user, but should cover at least one vegetation period.

| Management options
With the generic plant component, BODIUM can simulate monocultures as well as crop rotations depending on available data.Parametrizations can be directly selected for winter wheat, spring wheat, winter barley, spring barley, sugar beet, potato, silage maize and corn maize (Hunt & Loomis, 1979;Jones et al., 1986;Lenz, 2007;Yang et al., 2004).Mineral and organic fertilizer can be applied at a defined depth and with a specific C:N ratio.Tillage operations can be simulated for different depths as well.This is reflected by a changing pore size distribution and the mixing of solid components along the depth of the tilled layer.Moreover, the oxygen concentration is set to the atmospheric value within this layer.The decrease in bulk density by tillage is followed by re-settling during the following time period depending on soil depth.It is further possible to set negative effects of soil management on earthworms and soil microorganisms in terms of a reduction, which is however not yet parameterized.

| Programming tools
The BODIUM model is coded following the approach of object-oriented programming in the programming language C++ with the cross-platform development environment Qt (Version 6.2.3, https://www.qt.io/).It uses json files as input and is linked to an R script for directly analysing and plotting simulation results using the packages 'ggplot2', 'viridis', 'tidyft' and 'dplyr' (R Core Team, 2021).The source code including all input files used for the present study is available in a git repository (https://git.ufz.de/bodium/bodium_v1-0).Input by the user is needed for soil profile initialization, management and weather (Tables S4 and S5).Management comprises crop planting, nitrogen fertilization and tillage operations.

| Model validation
For model validation, we used data from long-term field experiments with a well-known management history, available weather data and a description of soil parameters.With that, all necessary boundary conditions are available and the model can be initialized properly.Depending on the specific motivation of the long-term field experiments, data sets including time series on different state variables are available.Typically, available data are plant yield, soil carbon, or nitrogen content.They always comprise only a part of the variables considered in the model but can be used for its validation and for a check of plausibility.This is also true for quantities that are measured only sporadically and are not monitored continuously.
In the following, data from the Static Fertilization Experiment Bad Lauchstädt were used (Merbach & Körschens, 2002;Merbach & Schulz, 2013).More information on the field experiment can be found on the experiments' website https://www.ufz.de/index.php?en=39220.The long-term field experiment is operated since 1902 and has an area of 4 ha split into eight strips, in which all crops of the rotation are grown simultaneously every year.The crop rotation is sugar beet-spring barley-potatoes-winter wheat until 2014, and silage maize-spring barley-silage maize-winter wheat thereafter.In total, 18 different fertilizer treatments are applied with three levels of organic fertilization (without, 20tha À1 farmyard manure (FYM) every 2 years, 30tha À1 FYM every 2 years) and six mineral fertilization treatments (without, NPK, NP, NK, N, PK) where the application rate depends on crop and FYM treatment.
We simulated four strips of the eight strips where the management was not essentially changed during the whole experiment for a duration of 40 years (1980-2019) so that we can assure no effects of management transition for the start of our simulations.As input, we used weather data from the weather station of the field site and documented management data on the date of seeding; date, type and amount of fertilization; and date and depth of tillage.As currently only nitrogen is implemented in BODIUM, we cannot simulate any limiting effects of phosphorous and potassium.In consequence, we have to assume an optimal supply of phosphorous and potassium within our simulations, reducing the mineral fertilizer treatments we are able to handle to two (NPK, PK).
Except for sugar beet, straw remains on the field after harvest and is incorporated during the next tillage event.Further, we include atmospheric nitrogen deposition with an amount of 50kgha À1 as a simple input rate distributed over the vegetation period (Merbach & Schulz, 2013).
The loamy soil used for the simulations is classified as Haplic Chernozem and the soil profile is initialized using the data provided in the related publication (Altermann et al., 2005;Leuther et al., 2022).Initial conditions were the same for all simulations except for nitrogen and carbon content.Here, we used the measured values of the year before the simulation period started.Initial partitioning of carbon pools was based on previous estimates for the stable pool (0.85) (Cécillon et al., 2021), and adjustment following 1 year of pre-simulation for the labile pools.
For comparison of our simulation results, we used the measured yield in terms of biomass of storage organs, yearly soil analyses of total organic carbon and total nitrogen in the depths of 0-20 cm, and of microbial biomass in the year 2012.See Merbach and Schulz (2013) for details on sampling and analysis.As an indicator for the model accuracy, we used the normalized root mean square error (nRMSE) with N the number of data, y s the simulated values, y o the observed values and y o,max and y o,min the maximum and minimum of the observed values respectively.Further, the model efficiency EF is calculated with y o the mean of the observed values.Here, an EF of 1 would indicate no difference between simulated and observed values, while a negative EF means that the model simulation outcome is worse than using the mean of the observed values.Additionally, we calculated the coefficient of determination r 2 , where appropriate.
To demonstrate the impact of the newly introduced concept of non-equilibrium water flow we compared the simulated yields with and without this feature (Ippisch et al., 2006).

| Explorative scenarios
In order to demonstrate the application of the BODIUM model for evaluating different management scenarios, two contrasting tillage schemes were simulated in an explorative way using the validation scenario of the static fertilization experiment as a basis.For the tillage scenario, the soil was loosened to the working depth between 5 and 28 cm, according to the documented management events.All components in the nodes down to the tillage depth were mixed by taking the average of the present values and the soil was aerated.The dates and depths were taken from the experimental station.
In a second no-till scenario all tillage practices were omitted and the soil structure was kept constant without any mixing.Due to the fact that we expect more earthworms to be active under no-till management, macropores were added to the soil structure.We assumed approximately 50 macropores per m 2 with a saturated conductivity of 10 À4 ms À1 (Chan, 2001).

| Evaluation of soil functions
For the comparison of the contrasting scenarios, we compare some quantitative measures that can be derived from the state variables of the model and provide substantial information for a direct evaluation of soil functions: Production: The yield is described as the biomass in storage organs at harvest for all simulated crops except for silage maize, where the overall aboveground biomass is used.Carbon storage: The total soil organic carbon stock summed up over the first 20 cm gives an expression for the carbon storage function.Nutrient cycling: As a measure for nitrogen buffering, the nitrogen use efficiency is calculated as the fraction of nitrogen that was uptaken by the plant from the total sum lost via plant uptake and nitrogen that has leached to the groundwater.Water storage and filtering: The available water capacity summed up over the first 30 cm gives the storage potential difference of the various scenarios.The difference was averaged over the whole time span.The filter function of soil is related to the adsorption of pollutants (which we do not cover here) and metabolization through biological activity.For the latter, we use the residence time of solutes within the soil profile weighted by the microbial activity, integrated over the first metre of the soil profile and over the vegetation period.

| Model validation
The simulated yield of plant biomass in storage organs (hereafter referred to as yield) showed good agreement with the measured yield for all treatments with an overall model efficiency (EF) of 0.71 and an nRMSE of 0.09.Further, the model is able to reproduce the general temporal dynamics of the crop rotation (Figure 2, see also Figures S3-S6).
For soil organic carbon, the overall trend between the simulated mean yearly outcomes for the different treatments fitted well with the observed data, with a model efficiency EF of 0.51 and a nRMSE of 0.12.However, when distinguishing between treatments and strips F I G U R E 3 Simulated soil organic carbon (lines) compared to observed soil organic carbon (dots) in topsoil (0-20 cm) for the six different treatments (rows) indicated by colour and the four strips (columns, SH, strip).
For the treatments without manure application, we observe a depletion of organic carbon over time which is faster than the observed decline.
The total nitrogen content shows a similar picture (Figure 4).
Although model accuracy for the overall simulations is reasonable (nRMSE: 0.15, EF: 0.31), model efficiency for the different treatments is between À13 and 0.08 and nRMSE ranges between 0.18 and 0.78.Again, we observe a faster depletion of total nitrogen for the treatments without manure application compared to the observed values.
Regarding microbial biomass, only data for one point of time in October 2012 are available for the different strips simulated here.As microbial biomass is a quite dynamic value, we compare the observed data with the daily simulated result of the whole simulation period (Figure 5).
In general, the model appears to overestimate the total microbial biomass.
To further understand the limiting processes for plant growth emerging from the model process interactions, we analysed the accumulated drought stress for each plant by summing up the total drought stress (according to S1.8.3) during the whole vegetation period (Figure 6).
When related to the simulated yield, the maximum yield limited by water becomes obvious as it decreases with increasing drought stress for all plants.Here, also an increase in nutrients in terms of higher fertilizer amount cannot further increase the yield.We observe a differentiation among the different crops as well as fertilizer treatments.The reason for the difference might be related to non-equilibrium water distribution along the profile.
To test the actual effect of the hydraulic nonequilibrium on yields, we repeated the simulations without non-equilibrium and compared the simulated yields (Figure 7).Overall, the yields are slightly higher and the prediction of the actual yield is slightly better (r 2 increased from 0.67 to 0.76, model efficiency from 0.59 to 0.73) for the simulations with non-equilibrium, however, the effect can be both positive or negative for specific years.The influence of modelled non-equilibrium differs depending on the water content profile and root development at times when heavy rainfalls occur.Two examples with different impacts on yields of sugar beet and maize are shown in Figures 8 and 9 for the years 2009 and 2017 respectively.In 2009, non-equilibrated deeper percolation occurs during the substantial rainfall events in June when plant roots are still distributed in the shallow soil and the water cannot be stored in the lower profile.This water is then missing later in the year when plants try to compensate the water deficacy by increased root growth around 50 cm depth.If non-equilibrium flow is not considered, the water remains in the shallower part of the profile and can support plant growth.In 2017, there was less rain in spring and some of the water that bypasses the topsoil can be stored in the lower, already desiccated profile.In the later season, the deeper developed plant roots can profit from this water.Regardless of nonequilibrium, the shallower water is used up until July.The following precipitation is distributed over a larger soil depth during non-equilibrium and does not remain near the surface, where it is partially lost through evaporation.The plants try to compensate for the deeper water shortage with root growth in the topsoil, but cannot balance the water loss.

| Explorative scenario simulations for tillage
As already stated earlier, the aim of BODIUM is to predict the effect of agricultural management on the identified five soil functions and their interrelated dynamics.In terms of productivity, no differences were observed between the tillage and no-till scenarios (Figure 10).As expected, higher yields are observed in the scenarios with nitrogen fertilization.
However, for the soil function 'carbon storage', the scenario without tillage shows a clear increase in topsoil organic carbon compared to the tilled scenarios for the fully fertilized treatment.This is true for all crops (Figure 11).Almost no differences are observed for the unfertilized treatment, where organic carbon input is mainly due to plant residues and roots.
As a measure for nitrogen buffering, the nitrogen use efficiency is used by calculating the fraction of nitrogen taken up by the plant from the total sum of lost nitrogen including uptaken nitrogen and nitrogen leached to the groundwater (Figure 12).
The nitrogen use efficiency is generally much higher in all N-unfertilized treatments-due to the nitrogen limitation in these simulations the plants take up almost all nitrogen.In the fully fertilized treatment, the nitrogen use efficiency is lower for all crops in the scenario with tillage (green).
In terms of water storage, the total available water capacity (AWC) within the top 30 cm (ploughing depth) was higher in the tillage scenario than in the no-till scenario, with an average increase of 5.7% across all simulations.However, the actual effect on the drought stress of plants is species-dependent because of different root systems.For silage maize with a deeper root system, the difference in drought stress between tillage and no-till is only 15% while for sugar beets with a more shallow root system, this amounts to 40%.The indicator for water filter was calculated based on the water residence time and the microbial activity, and was analysed for the months in which pesticides are usually applied (March-June).Higher values translate into a better filter capability due to a longer contact time of pesticide and degrading microbes and higher microbial activity in general.Overall, the filtering ability increases throughout the vegetation period with higher values in June compared to March due to a higher microbial activity-resulting from warmer temperatures and increasing energy sources, for example, organic matter and root exudates (Figure 13).
Although we observe significant differences between the tillage scenarios for some months and for both fertilization treatments, the filter function of both treatments is nearly identical.This is due to the high complexity of this function, as it can be either increased by an enhanced microbial activity or by a higher water residence time, which both are highly affected by soil structure dynamics induced by tillage.Our simulation results suggest that we successfully demonstrated the plausibility of the BODIUM model with observed data from the described field experiment.The validation results for yield fit reasonably well considering that the parameters are not calibrated and plant growth is modelled generically.The yields of the five crops differ in the simulation accuracy, indicating that optimization of plant parametrization would enhance model fit (see also Figure S2).The simulation of crop rotations is a challenging task also for dedicated plant growth models including multi-model ensembles (Kersebaum, 2022;Kollas et al., 2015), which is mainly improved by intensive calibration where comprehensive data sets are required (Kersebaum et al., 2015).However, the aim of BODIUM is to simulate dynamic trends in reaction to climate and management, but the focus is not on simulating accurate yields.Of course, model error of specific crops can accumulate in the following years and also affect other soil functions, and thus simulating the dynamic response is precious.We show that BODIUM is able to capture the trend in yield dynamics for most years, including extreme events such as drought.One example is the simulation of sugar beet in the years 1987 and 1991 (Figure 2).Although we have the highest deviation from observed results with sugar beet in general, we capture the high yield in the wet year of 1987 as well as the rather low yield following a drought period in 1991.However, there are also years where the simulations are far away from the observations-because our model cannot account for events like diseases, human failure, or some extreme weather events such as strong or very hot winds.The latter might explain the high differences between the simulation and observation of silage maize yields in 2019, where observed yields were quite low, especially in the manure treatments.
We also show that BODIUM simulates the relative differences in response to management, here the fertilization regime, compared to the observed differences (Figure 14).
The described differences in the model efficiency between the fertilizer treatments may indicate an underestimation of organic matter input, for example, resulting from crop growth and/or nitrogen mineralization due to biological processes.Indeed, the former is corroborated by the fact that the simulated soil organic carbon content for the treatments without organic fertilization is lower compared to the observed values.The soil organic carbon dynamics is also highly sensitive to the rates describing the stabilization and destabilization, which are still lumped values and not mechanistic in this model version.The maximum rates are similar for all treatments, and the actual processes are only indirectly affected by biological activity, as they depend on the current pool size.This is a pretty good example of how important a real mechanistic description is to simulate the complex response to management and environmental fluctuations.This is further strengthened by the differences in the various strips: we know from field measurements that there are differences in soil parameters such as pH between the different strips (Dierke & Werban, 2013;Werban et al., 2009), which we currently cannot simulate with our model approach, but which would likely have an influence on the dynamics of the whole system.Also, the fraction of stable carbon might differ between the strips and treatments, which would affect the simulated carbon dynamics, as we currently consider the same initial partitioning for all strips and treatments.
The underestimation of soil nitrogen indicates that the model still misses some processes related to nitrogen input like biological nitrogen fixation or nitrogen excretion by soil fauna (Lang & Russell, 2022).Another possibility is that the C:N ratio of the organic matter is not correct, for example, because of inadequate distribution of nitrogen in the different plant organs.This could lead to a higher C:N ratio of incoming plant material (roots, straw).Of course, an increased organic carbon content would also increase the total nitrogen content.
The overestimation of microbial carbon in our simulation results is in accordance with other microbial soil models of croplands (Horrigue et al., 2016;K. Wang et al., 2017).As mentioned, the microbial biomass is a quite dynamic value and therefore it is already promising that the model output is in the right range.The microbial biomass can fluctuate due to incoming plant material, root exudation, manure application, tillage events, interactions with soil fauna, pesticide application and local changes in temperature, soil water content and aeration.The two scenarios without organic fertilization show the best fit between observed and simulated data, likely because for these treatments the local conditions were less dynamic.However, together with the comparison of organic carbon and nitrogen content-where these two scenarios show the most differences-this may indicate that the stoichiometric considerations are not well balanced.Here, the simulated microbial biomass fits real observations, while organic carbon and nitrogen contents are underestimated.This may be due to a low carbon use efficiency resulting in a higher carbon loss via respiration and, at the same time, an imbalance of mineralization and immobilization.This can be caused either by a high C:N ratio of the soil organic matter assimilated by the microbes or a low internal C:N ratio of the microbes.In the current BODIUM version, the internal microbial C:N ratio is a static value, although the homeostatic flexibility of microbes, especially fungi, is discussed within the scientific community (Camenzind et al., 2021;Strickland & Rousk, 2010).In future model extensions, it is planned to also account for a more flexible C:N ratio as is already implemented by other modelling approaches (Kyker-Snowman et al., 2020;Manzoni et al., 2021).
In terms of drought stress, the differentiation in the effects between different crops can be explained by their traits: less drought stress of rather deep-rooting crops (silage maize) with increased drought resistance and high drought stress of crops with a higher overlap of the vegetation period and dry summer days (sugar beet).However, in the case of high nutrient deficiency, the limitation due to nutrients is more important than water limitation as the treatments with no nitrogen fertilization mostly show similar yields independent of the drought stress (yellow dots in Figure 6).
Non-equilibrium bypass flow has an impact on water and plant root development, it also increases nitrogen losses in the modelled scenarios.The annual yield pattern indicates that considering these dynamic effects can improve the prediction of the model.We are currently conducting a study on the validity and magnitude of the process using water content profiles from lysimeters and field observations.Earlier studies suggest also a considerable impact of this process (Hannes et al., 2016).

| Explorative scenario simulations for tillage
In summary, managing this specific site without tillage would cause an increase in soil carbon content and nitrogen use efficiency.Both are caused by the increased mineralization of organic matter in the wellaerated tilled topsoil and the related mobilization of nitrate.Moreover, this increased mineralization is also located in the lower part of the tilled horizon where root uptake is somewhat reduced and more nitrogen might be leached.At the same time, the available water capacity of the topsoil is reduced which might be critical in dry years.This example demonstrates the potential of the BODIUM model to evaluate the various soil functions simultaneously.The evaluation of agricultural management practices is typically done by analysing yields which are interpreted as an integral measure of soil quality and which is easily accessible also for farmers.Interestingly, the yield is actually not affected by the different tillage treatments while some of the other soil functions are.Thus, the model simulations allow for a much more detailed analysis of the soil system as a whole and, herewith, a more differentiated evaluation of the impact of management practices.

| OUTLOOK AND CONCLUSION
The BODIUM model presented here is a systemic soil model, developed in a modular way.This enables an easy integration of further components or processes of the soil system as well as a continuous improvement of the implemented modules according to newly developed process knowledge.First, model simulations for long-term field experiments demonstrated the plausibility of the implemented process interactions for soil as a complex system.To do so, the model was applied in the forward mode, that is, no model parameters are calibrated to measured data.
In principle, there are no restrictions for possible future extensions.We will focus here on a few points which are currently under development and will be addressed in the next model version.
A more comprehensive description of the nutrient cycle is crucial for plant growth and the dynamics of soil microbial activity based on stoichiometric principles (Achat et al., 2016).In particular, the dynamics of phosphorous will be included which will be built upon concepts that have been developed recently (Das et al., 2019;Nakhavali et al., 2022;Yu et al., 2020).Further, the effective rate parameters currently used to model the turnover of soil organic matter need to be replaced by mechanistic process descriptions.This is especially true for the stabilization of organic matter in interaction with the mineral soil matrix.We will consider the diffusion of dissolved organic carbon as an important process for mixing organic and mineral soil components.Additionally, bioturbation will be included which affects the depth distribution of surface litter and the mixing of mineral and organic components within individual soil layers (Meurer et al., 2020).This process is also influenced by agricultural management, thus, the effects of tillage and nitrogen fertilization on soil fauna will be implemented in the future based on information derived from different metaanalyses (Betancur-Corredor et al., 2022a, 2022b;Briones & Schmidt, 2017).
An unresolved critical question is in what detail small-scale microbial ecology in soils should be represented in models aiming at larger-scale soil functions (Smercina et al., 2021).Currently, this is lumped into an effective microbial pool.In a follow-up version, we will distinguish bacteria and fungi since they have some fundamentally different characteristics in terms of optimal living conditions, their stoichiometry and their way to explore the soil volume.This will be an important step towards a better characterization of the functional diversity of the soil microbiome.
The temporal changes in soil structure will be extended for compaction due to heavy machinery and for the impact of root growth, including cover crops and the activity of burrowing soil organisms.While erosion threatens all soil functions, it cannot be represented in our 1D model.However, we will address erodibility as a function of soil structure, texture and organic matter so that an interface to erosion models can be established.Our vision for the future development of BODIUM is to provide the source code and the related programming infrastructure for the wider scientific community.The process descriptions should be easy to edit even for users without in-depth programming knowledge.This can be accomplished using an appropriate metalanguage.In this way, BODIUM can be constantly updated according to the development of the current state of knowledge.Moreover, alternative hypotheses on soil processes can be investigated in terms of their implications for the context of the entire soil system.We believe that such a freely available community model will be a valuable tool to foster an improved understanding of soils and their response to external drivers, which is urgently needed under the current challenges of rapidly changing conditions in terms of land use and climate worldwide.
David Russell, Karin Schmelmer, Bastian Stößel and Andrey Zaytsev for valuable discussions and expert input during the conceptual model development.Data sets from the Static Fertilization Experiment Bad Lauchstädt were provided by the Helmholtz Centre for Environmental Research-UFZ.We thank the colleagues of the experimental field site for sampling and provision of input and validation data, especially Ines Merbach.Open Access funding enabled and organized by Projekt DEAL.

F
I G U R E 2 Yearly simulated yield (filled dots) compared to observed yield (unfilled dots) for four selected treatments.Crops are indicated by colour and difference between simulated and observed by line, with dashed lines indicating an overestimation of the simulation.Yearly precipitation is shown with light blue bars.

F
I G U R E 4 Simulated soil nitrogen content (lines) compared to observed soil nitrogen content (dots) in topsoil (0-20 cm) for the six different treatments (rows) indicated by colour and the four strips (columns, SH, strip).

F
I G U R E 5 Daily values of simulated microbial biomass in topsoil (0-20 cm) for the six different treatments indicated by colour (bars).Observed microbial biomass in topsoil (0-20 cm) of one sampling from October 2012 (points) for the six different treatments indicated by colour.F I G U R E 6 Accumulated drought stress related to simulated plant biomass at harvest for the six different treatments indicated by colour and for the five crops indicated with the different symbols.

F
I G U R E 7 Comparison of biomass production with and without dynamic non-equilibrium for all crops.

F
I G U R E 8 Development of profiles for nitrogen, roots and water content indicated by line type of simulations with (blue) and without (green) dynamic non-equilibrium of selected days in 2009, and the corresponding plant biomass and drought stress for 2009.

F
I G U R E 9 Development of profiles for nitrogen, roots and water content indicated by line type of simulations with (blue) and without (green) dynamic non-equilibrium of selected days in 2017, and the corresponding plant biomass and drought stress for 2017.F I G U R E 1 0 Plant yield for the two tillage scenarios indicated by colour and the two fertilization treatments for the five different crops.F I G U R E 1 1 Total organic soil carbon in topsoil (0-20 cm) for the two tillage scenarios indicated by colour and the two fertilization treatments for the five different crops.

F
I G U R E 1 2 Nitrogen use efficiency for the two tillage scenarios indicated by colour and the two fertilization treatments for the five different crops.F I G U R E 1 3 Daily values of the filter indicator for the two tillage scenarios indicated by colour and for the two fertilization treatments for the months of March, April, May and June.

F
I G U R E 1 4 Comparison of the difference in yield of observation (blue) and simulation (green) for different fertilization regimes.