Interactions between soil structure dynamics, hydrological processes, and organic matter cycling: A new soil‐crop model

The structure of soil is critical for the ecosystem services it provides since it regulates many key soil processes, including water, air and solute movement, root growth and the activity of soil biota. Soil structure is dynamic, driven by external factors such as land management and climate and mediated by a wide range of biological agents and physical processes operating at strongly contrasting time‐scales, from seconds (e.g., tillage) to many decades (e.g., faunal activity and soil aggregation). In this respect, positive feedbacks in the soil–plant system may lead in the longer term to soil physical degradation or to the recovery of structurally poor soils. As far as we are aware, no existing soil‐crop model can account for such processes. In this paper, we describe a new soil‐crop model (USSF, Uppsala model of Soil Structure and Function) that accounts for the effects of soil structure dynamics on water and organic matter cycling at the soil profile scale. Soil structure dynamics are expressed as time‐varying physical (bulk density, porosity) and hydraulic properties (water retention, hydraulic conductivity) responding to the activity of biological agents (i.e., earthworms, plant roots) and physical processes (i.e., tillage, soil swell‐shrink) at seasonal to decadal time‐scales. In this first application of the model, we present the results of 30‐year scenario simulations that illustrate the potential role and importance of soil structure dynamics for the soil water balance, carbon storage in soil, root growth, and winter wheat yields on two soils (loam and clay) in the climate of central Sweden. A sensitivity analysis was also performed for these two scenarios using the Morris method of elementary effects, which revealed that the most sensitive parameters controlling soil structure dynamics in the USSF model are those determining aggregation induced by organic matter turnover and swell/shrink. We suggest that the USSF model is a promising new tool to investigate a wide range of processes and phenomena triggered by land use and climate change. Results from this study show that feedback in the soil‐crop system mediated by the dynamics of soil physical and hydraulic properties are potentially of central importance for long‐term predictions of soil water balance, crop production, and carbon sequestration under global change.

Soil structure regulates critical processes such as water, air and solute movement, biological activity, and root growth, which directly or indirectly affect the ecosystem services provided by soil, such as crop production, the retention/removal of pollutants, water regulation and carbon sequestration, and climate regulation.The structure of agricultural soils is also dynamic at time scales ranging from seconds (e.g., tillage and compaction) to seasons and decades (e.g., root growth and the activity of soil biota; Meurer, Barron, et al., 2020).Soil structure dynamics triggered by changes in soil and crop management practices or climate (i.e., the critical external drivers) lead to feedbacks in the soil-plant system, which may result in the long term in either soil physical degradation or restoration (Henryson et al., 2018;Blanchy et al., 2023; Figure 1).
Soil-crop models that simulate water, matter, and energy flow in the soil-plant-atmosphere continuum are widely used as tools to evaluate the effects of management and climate on crop production and the environmental impacts of agriculture such as soil carbon sequestration and nutrient leaching (e.g., Constantin et al., 2019;Jones et al., 2017;Stöckle & Kemanian, 2020).Simple descriptions of within-season variations of soil physical and hydraulic properties induced by tillage and subsequent consolidation have been employed in some soil-crop models (e.g., Chandrasekhar et al., 2018;König et al., 2023;Maharjan et al., 2018).However, we are not aware of any existing model that can account for the effects of the broader range of physical and biological processes that determine soil structure dynamics at the much longer time scales relevant for agricultural systems under global change (i.e., from years to decades).Thus, in principle, no existing soil-crop model can forecast the effects of changes in land use and management, and climate on soil hydrological processes and crop production that are mediated by long-term changes in soil structure and soil physical and hydraulic properties (e.g., Vereecken et al., 2016;Vogel et al., 2018).This has hampered progress towards a sound understanding of the long-term impacts of changes in soil and crop management and climate on soil physical and hydraulic properties and the ecosystem services provided by soil.
The USSF model (Uppsala model of Soil Structure and Function) is designed to fill this gap.It describes the interactions between soil structure dynamics and soil hydrological processes, crop production, and organic matter cycling in the soil-plant-atmosphere system at the soil profile scale, for time scales ranging from seasons to centuries.In this paper, we provide a full description of the complete model and also present the results of 30-year scenario simulations that illustrate the potential role and importance of soil structure dynamics for the soil water balance, carbon storage in soil, root growth, and winter wheat yields on two soils (loam and clay) in the climate of central Sweden.Based on these scenario simulations, sensitivity analyses were performed to identify the most influential parameters regulating soil structure dynamics with respect to grain yields, carbon sequestration and the components of the water balance.

| MODEL SCOPE AND OVERVIEW
USSF describes water and organic matter cycling in the soil-crop system in a one-dimensional soil profile divided into four soil horizons, with the horizons sub-divided into

Highlights
• A new soil-crop model is described that accounts for the effects of soil structure dynamics • 30-year scenario simulations are performed to illustrate soil physical restoration due to manuring • A sensitivity analysis showed the importance of parameters regulating aggregation and swell/shrink • The model could be used to simulate degradation or recovery caused by land use or climate change numerical layers of variable thickness.The model, which is designed primarily for use as an explorative tool to investigate processes of soil physical degradation and restoration in agricultural systems, includes the most important interactions and feedbacks between soil structure, soil hydrology, and organic matter cycling in the soilplant system.These interactions determine the dynamics of soil structure and physical properties and their impacts on important ecosystem processes and services (e.g., crop growth and carbon sequestration) at temporal scales ranging from seasons to centuries (see Figure 1).Given that these are the spatial and temporal scales of interest, soil processes in the USSF model are represented at the continuum macroscopic scale rather than at the pore scale.However, to the extent that is possible, process representations in USSF have been informed by knowledge gained from pore-scale studies.In developing USSF, mechanistic process descriptions have been adopted when possible, with the aim to satisfy the principle of model parsimony and avoid all unnecessary empiricism (Jarvis, Larsbo, et al., 2022).Phenomenological approaches have been employed only where the understanding of the relevant underlying processes was considered insufficient to derive a more mechanistic model or where the parameter requirements of such an approach would be too demanding.In particular, some of the biological processes that alter soil structure (e.g., soil faunal activity; Young et al., 1998;Lucas et al., 2019) are only considered implicitly.This paper describes a first version of the USSF model and, as such, does not account for all processes affecting the dynamics of soil structure.In particular, nutrient cycling in the soil-crop system is currently neglected which means that crop growth can only be reliably predicted for conditions when nutrients are not limiting (e.g., fully fertilized crops).Furthermore, the effects of traffic compaction are not yet included in USSF.Nevertheless, this first version of the USSF model should be well suited to investigate the rates of recovery or regeneration by biological processes of the structure of previously compacted soil.The preferential flow of water in soil macropores (Jarvis, 2007;Nimmo, 2021) is another process that is not considered in USSF, which at first sight may seem surprising given that soil structure is a central feature of the model.Surface-connected macropores are known to determine the partitioning of incoming precipitation between infiltration and surface runoff (e.g., Beven & Germann, 1982;Fatichi et al., 2020;Or, 2020), and this critical aspect of the soil water balance is therefore captured by the model.In contrast, preferential water flow through macropores does not have such a strong impact on the partitioning of the infiltrated water between storage, evapotranspiration, and drainage from the soil profile (Jarvis, 1998), although it can dramatically affect the transport of solutes, particularly those that are strongly adsorbed or quickly degraded (Jarvis, 2007).The USSF model currently only simulates water flow in the soil-plant system and not the transport of solutes.If solute transport processes were included in future versions of USSF, the water flow model could, for example, be extended to a dual-permeability formulation to account for the impacts of preferential (i.e., non-equilibrium) flow (e.g., Šimu ˚nek et al., 2003).
It is not straightforward to validate a model such as USSF.This is because the model predicts structure dynamics and their effects at vastly contrasting time F I G U R E 1 Schematic illustration of positive feedbacks between crop production, carbon sequestration, soil structure and physical conditions.Note that the outer loop can also go into reverse, resulting in poorer crop growth, a loss of organic carbon and soil physical degradation.
scales from single seasons (e.g., runoff due to soil sealing) to many decades (e.g., earthworm effects on macroporosity, organic matter dynamics, and soil aggregation).Suitable field and experimental data will almost always be lacking for a test of all the model components at one site.In particular, soil organic matter contents and crop yields are often measured in long-term field trials set up to study the effects of tillage or crop management systems, but with only a few exceptions (e.g., Riley et al., 2008), soil physical and hydraulic properties are seldom measured, most often only at the beginning of the experiment.As testing the entire USSF model against data from one single field experiment seems impossible, we advocate testing one or more of the individual components of the model separately for different sites and experiments.This is the approach we have adopted in the development of USSF.Thus, the model describing interactions between soil physical and hydraulic properties and soil organic matter contents has been tested against time series data on soil organic carbon contents and occasional measurements of bulk density, porosity, soil surface elevation, and soil water retention curves available for a long-term field experiment at Ultuna in central Sweden (Meurer, Chenu, et al., 2020).The USSF module describing the growth of an arable crop has also recently been tested using comprehensive data for winter wheat obtained for one growing season at the same site (Ultuna) on the same soil type, although not on the same plot (in preparation).Earlier versions of the routines describing the effects of soil fauna and plant roots on soil structure were tested by Meurer, Barron, et al. (2020) using four years of data from a compaction recovery experiment in Switzerland (Keller et al., 2017), while the hydrological model has been tested against comprehensive data obtained in weighing lysimeters during six years at two sites with contrasting climates in north-west Germany (Jarvis, Groh, et al., 2022).Any test of one or more of the algorithms in a soil-crop model like USSF invariably involves some calibration.In building the model, we attempted to minimize the need for calibration by making use of 'hard-wired' pedotransfer functions as a way to estimate some of the soil parameters.Other parameters in the model, for example, some of those related to crop growth and development, can be estimated from literature data.

| Model description
Soil structure is considered in terms of the dynamics of bulk density, total porosity, and the soil water retention and hydraulic conductivity functions employed in the Richardson-Richards' equation used to calculate soil water flow, both of which reflect the underlying soil pore size distribution.The connectivity of macropores is another potentially important aspect of the pore space structure that is also treated in the USSF model, albeit in a rather simple fashion.
In this section, we focus primarily on how the USSF model describes the interactions between soil structure, expressed in terms of dynamic soil physical properties and hydraulic functions, and important soil processes such as root growth and soil organic matter turnover.Additional model algorithms dealing with abiotic processes in the soil (i.e., soil water balance and flows, soil temperatures), evapotranspiration, and the growth of an arable crop are described in detail in the supplementary information.

| Soil structure dynamics
Conceptual model of soil solids and pore space The soil solids in USSF comprise mineral matter and organic matter.Both are characterized by a fixed density and, thus, a constant relationship between volume and mass.The mass of mineral matter is assumed to be constant in any given soil horizon, while the mass of organic matter is dynamic.The pore space in the USSF model is classified according to both origin, with textural pore space distinguished from structural pore space, and pore diameter (pore size classes, see Figure 2).The total porosity is divided into three size (diameter) classes, namely macroporosity, mesoporosity, and microporosity, with the latter two classes comprising the matrix porosity.The structural pore space in USSF is found in all three of these size classes, while textural pore space is confined to the matrix.Macroporosity comprises three types of dynamic macropores generated by both physical (i.e., tillage, soil shrinkage cracks; Bronswijk, 1991) and biological processes (i.e., biopores created by plant roots and soil fauna; Meurer, Barron, et al., 2020) as well as a permanent (i.e., static) macroporosity.Structural porosity in the soil matrix consists of aggregation pore space generated by microbial activity and the turnover of soil organic matter.A fundamental assumption in the USSF model is that the volume of this aggregation pore space in the mesopore and micropore pore size ranges can be expressed as a linear function of the volume of soil organic matter stored in each pore region (e.g., Boivin et al., 2009;Emerson & McGarry, 2003;Meurer, Chenu, et al., 2020).In addition to opening and closing soil cracks, soil shrinkage and swelling also change the volume of the matrix pore space (e.g., Boivin et al., 2009;Bronswijk, 1991).

Phase relations
As a consequence of the processes generating soil structure dynamics, the thickness of a soil layer in the soil profile in USSF is variable.The thickness of a soil layer Δz (m) is given by: where A xs is a nominal cross-sectional area (= 1 m 2 ), and the total soil volume V t (m 3 ) is: where V p is the volume of pores (m 3 ), V s is the volume of solids (m 3 ), and V mat , V crack , V bio , V till and V sta (m 3 ) are the volumes of soil matrix pores, cracks, biopores, tillage voids, and static macropores, respectively.The volume of solids in Equation ( 2) is given by: where Δz ref and ϕ ref are reference values for the thickness (m) and matrix porosity (m 3 m À3 ) of a fully swollen layer of mineral soil without organic matter and aggregation, M som is the mass of organic matter in the layer (kg m À2 ), and γ som is the density of organic matter (kg m À3 ).

Matrix porosity
The matrix porosity comprises both textural pore space and aggregation porosity in two size classes (i.e., micropores and mesopores), and depends on soil texture, organic matter content, and swelling and shrinkage (Figure 2).For the sake of simplicity, we assume a linear (or proportional) soil shrinkage characteristic of slope p, where p = 0 holds for rigid soils and p = 1 is equivalent to so-called "normal" (or basic) shrinkage (i.e., the volume shrinkage of the matrix equals the volume of water lost).Clearly, under dry conditions, this linear model will not hold, as volume shrinkage usually becomes increasingly less than the volume of water lost (e.g., Groenevelt & Grant, 2001;Stewart et al., 2016).Neglecting this "residual" shrinkage phase should not unduly affect simulations of the soil water balance, since it is normally only strongly expressed when the soil is very dry and water flow is therefore very slow.For our simple linear model, the matrix pore volume V mat in Equation ( 2) can be expressed as: Schematic illustration of soil pore space structure and types of soil pores in the USSF model (Cajsa Lithell, Redcap Design).
where V w is the volume of water in the soil layer (m 3 ) and V mat(ref ) (m 3 ) is the matrix pore volume at the reference fully swollen state, which is given by: where f agg is the so-called aggregation factor (m 3 pores m À3 organic matter), which determines the volume of aggregation pore space created by soil organic matter (Meurer, Chenu, et al., 2020).
Combining Equations ( 4) and ( 5) gives the matrix porosity ϕ mat (= V mat /(A xs Δz)) as: We further assume that the size distribution of the textural pore space is unaffected by soil shrinkage, which means that the microporosity ϕ mic (m 3 m À3 ) can be expressed as: and the mesoporosity ϕ mes (m 3 m À3 ) as: where M som(mic) is the store of soil organic matter (kg m À2 ) in the micropore region and f mic(text) is the fraction of the textural pore space comprising micropores.In principle, f mic(text) could be estimated using one of several models that predict the soil pore size distribution from the particle size distribution (e.g., Arya & Heitman, 2015;Chang et al., 2019;Mohammadi & Vanclooster, 2011).However, USSF employs a simpler empirical approach to estimate f mic(text) which makes use of the observation that the water content at wilting point is usually strongly correlated with clay content but hardly at all with organic matter content (e.g., Kätterer et al., 2006;Ostovari et al., 2015;Pollacco, 2008).For the Brooks-Corey water retention model that is employed in USSF (see Equation ( 38) in "Soil hydraulic properties"), we can write: where ψ mes/mac and ψ mic/mes are pressure heads (m) defining the largest diameter mesopore and micropore and λ mat(t) is the Brooks-Corey pore size distribution index for a soil consisting only of textural pore space in the matrix (i.e., without organic matter), estimated from: where ψ w (m) is the wilting point pressure head (= À150 m) and the corresponding water content θ w (m 3 m À3 ) is estimated from the pedotransfer function derived by Ostovari et al. (2015) for European and North American soil data in the HYPRES and UNSODA databases: where f clay (kg kg À1 ) is the soil clay content.

Crack porosity
A simple approximate expression for the volume of cracks as a function of the extent of matrix shrinkage can be derived assuming that the change of the soil layer thickness is small in relation to the initial thickness (e.g., Aitchison & Holmes, 1953;Bronswijk, 1991): where r s is a dimensionless geometry factor (1 ≤ r s ≤3) such that r s = 3 for three-dimensional isotropic shrinkage and r s = 1 for unidimensional shrinkage (i.e., subsidence).
We assume here that r s is a constant, even though it may vary with soil water content (te Brake et al., 2013).

Bioporosity
A phenomenological model is used in USSF to represent the generation of biopores by root decay and faunal activity (Meurer, Barron, et al., 2020) and their partial destruction by soil tillage.Note that this simple heuristic model describing the effects of macro-fauna is only appropriate for temperate climate regions where earthworms are dominant "ecosystem engineers" (Blouin et al., 2013).
The change in soil bioporosity ϕ bio (m 3 m À3 ) (= V bio / (A xs Δz)) in USSF is given by: where I r (kg m À2 d À1 ) is the root decay rate (see equation 25), γ r is the specific density of roots (kg m À3 ), τ e and τ a (m 3 m À3 d À1 ) are the rates of bioporosity formation due to endogeic and anecic earthworms, k bio(loss) is a lumped first-order rate constant (d À1 ) to account for the loss of bioporosity as a consequence of soil in-filling by various biological and physical processes (Meurer, Barron, et al., 2020), Γ d (À) is a binary variable to indicate whether a tillage event occurs (Γ d = 1 on the day of tillage, zero otherwise), z* (0 ≤ z* ≤ 1) is the proportion of the soil layer included within the cultivated layer, the depth of which can be set to different values for harrowing and ploughing, and f d (À) is the proportion of bioporosity lost during a tillage event.
The formation rate constants τ e and τ a will depend on earthworm biomass and activity, neither of which is likely to remain constant under changing agroenvironmental conditions.However, rather than attempting to explicitly model earthworm populations, it seems preferable to use simpler indirect approaches to estimate their effects on macroporosity.As a more comprehensive pedotransfer scheme that could account for a wider range of the relevant factors affecting earthworm populations (e.g., soil organic matter, soil texture, bulk density, pH) is not available, we make use of a simple approach based on only one proxy variable.Thus, USSF currently assumes that the effects of endogeic earthworms vary with the soil organic matter content f som (kg kg À1 ) (Eriksen-Hamel et al., 2009;Fonte et al., 2009;Hendrix et al., 1992;Kladivko et al., 1997;Krück et al., 2006) expressed in the form of a threshold function: where f som(c,e) is a critical value of soil organic matter content and τ e,max is a maximum turnover rate (d À1 ) attained at an organic matter content of f som(m,e) (kg kg À1 ).The bioporosity generated by anecic earthworms is assumed to be determined by the inputs of above-ground crop residues and amendments I a (kg m À2 d À1 ) (Eriksen-Hamel et al., 2009;Ouellet et al., 2008): where z max is the maximum depth in the soil (m) exploited by burrowing anecic earthworms and f bio (m 3 kg À1 ) is the volume of biopores generated per mass of organic residue supplied.It can be noted that residues and amendments are input to the soil at harvest and ploughing, so that Equation 15gives an instantaneous impact on bioporosity.This is, of course, somewhat unrealistic, as such effects would occur over periods of weeks or even months.However, such timing errors should not be serious when the main focus of USSF is on longer-term dynamics.

Tillage porosity
We make use of a simple algorithm to simulate how tillage, subsequent consolidation, and soil surface sealing affect the soil macroporosity in cultivated horizons.The change in the porosity of macrovoids originating from soil tillage ϕ till (m 3 m À3 ) (= V till /(A xs Δz)) is given by: where τ till (m 3 m À3 day À1 ) represents the gain of porosity due to a tillage event, k con is a rate constant (d À1 ) to account for the loss of tillage porosity due to subsequent consolidation, R is the rainfall rate (cm d À1 ) f int is the fraction of the rain intercepted by the canopy, which is approximated as a function of canopy leaf area index using Beer's law (see equation S16 in the supplementary information), Γ 0 is a binary variable indicating the uppermost soil layer (for which Γ 0 = 1, zero otherwise), k seal is a rate constant for soil surface sealing (cm À1 ), and C (À) is a factor accounting for the difference in kinetic energy between drip and throughfall due to drop diameter and velocity (Moss & Green, 1987;Nearing & Bradford, 1987): where h (m) is the average height from which the canopy drip is generated, which is assumed to be equal to half the crop height, h min (m) is a threshold above which the sealing potential of the drip throughfall becomes significant, h max (m) is the height from which the velocity of the drip throughfall approaches the terminal velocity and C max (À) is a weighting factor to account for larger kinetic energy at terminal velocity of canopy drip compared to throughfall due to a larger average drop diameter.The sealing rate constant, k seal is assumed to depend on the organic matter content of the soil in the uppermost numerical soil layer (Grønsten & Børresen, 2009;Le Bissonais & Arrouays, 1997): where k s(max) (cm À1 ) is a maximum value attained when the organic matter content is less than a critical value f som(c,s) and f som(n,s) is an organic matter content above which sealing is negligible.

Total macroporosity and percolating macroporosity
The total macroporosity ϕ mac (m 3 m À3 ) results from all processes described above and is given as: where ϕ crack is the crack porosity (m 3 m À3 ) and ϕ sta is the permanent (static) macroporosity (m 3 m À3 ).Only a certain fraction of the total soil macroporosity will be longrange connected (Hunt et al., 2014;Jarvis et al., 2017).This so-called "percolating" macroporosity ϕ mac( p) (m 3 m À3 ) influences macropore hydraulic conductivity and root growth in USSF and is given by: where f con (À) is the long-range connected fraction, which is assumed to increase linearly once a percolation threshold ϕ mac(t) (m 3 m À3 ) is exceeded (Jarvis et al., 2017), reaching a maximum value f con(max) at a user-defined upper limit value of ϕ mac(u) : Layer thicknesses and soil bulk density Combining Equations ( 1) to (5) with Equation ( 12) gives the soil layer thickness as: while the bulk density γ b (kg m À3 ) is calculated as: where γ min (kg m À3 ) is the density of soil mineral matter.

| Impacts of soil structure dynamics on soil processes
In USSF, the soil structure dynamics outlined above affect root growth, soil organic matter turnover, and soil hydraulic properties.These relationships are described in the following.

Root growth and turnover
In contrast with many existing crop models, the crop growth module in USSF places as much emphasis on below-ground production as on above-ground growth and crop yields, recognizing the importance of crop rooting for soil structure (e.g., Uteau et al., 2013;Yunusa & Newton, 2003) and organic matter turnover and persistence (Kätterer et al., 2011;Rasse et al., 2005).In USSF, the root growth of an annual crop is affected by soil structure.In particular, it is assumed that connected macroporosity sustains the growth of roots downwards in the soil profile by allowing them to avoid strong soil layers (Gao, Hodgkinson, et al., 2016;Hatano et al., 1988;Jin et al., 2013;Kautz et al., 2013;White & Kirkegaard, 2010).Soil strength in the USSF model is calculated using a pedotransfer function for penetrometer resistance (Gao, Whalley, et al., 2016).
We employ a diffusion model to mimic root growth and the downward penetration of roots in the soil profile (e.g., Acock & Pachepsky, 1996;Dathe et al., 2014;Wang et al., 2021).With this model, changes in the root biomass in any soil layer B root (kg m À2 ) are given by: where z is depth (m), D is the root diffusion coefficient (m 2 d À1 ), and I r is the root decay rate (kg m À2 d À1 ), which is assumed to be negligible when the plants are young: where Γ h (À) is a binary variable to indicate root death at harvest (Γ h = 1 on the day of harvest, zero otherwise), k d is a rate constant for root decay during the growing season (d À1 ), and Γ d1 indicates that the crop development stage is earlier than S d1 (Γ d1 = 1, zero otherwise; equations S24 to S26 in supplementary information).The lower boundary condition to Equation ( 24) is given as a zero flux, while the upper boundary condition is the supply rate of assimilates from above ground that is available to be utilized for root growth A roots (kg m À2 d À1 ), given by: where A (kg m À2 d À1 ) is the assimilation rate of dry matter by the crop, f bg (À) is the fraction that is allocated below ground, and f ex (À) is the fraction of this belowground production which is directly input to the soil as root exudates, and therefore, does not contribute to root growth.Total below-ground production (= A f bg ) is calculated as the minimum of a demand for assimilates, A bg(d) (kg m À2 d À1 ), and the potential supply of assimilates to the roots (kg m À2 d À1 ): where f bg(max) is the maximum fraction of assimilates allocated below-ground which depends on the crop development stage (see Supplementary information).The below-ground demand for assimilates is assumed to be proportional to the sum of the current root biomass in each layer i modified by factors representing the limiting effects of aeration f a (0 ≤ f a ≤1, equation S28 in supplementary information), soil strength f s (0 ≤ f s ≤1), and temperature (0 ≤ f t,r ≤ 1): where R bg(pot) is a first-order rate coefficient (d À1 ).The soil strength factor f s is calculated from soil water potential, soil depth, and bulk density using a pedotransfer function for soil penetration resistance developed by Gao, Whalley, et al. (2016), while f t,r is estimated empirically as a piece-wise linear function of simulated soil temperature (S31 in Supplementary information).
The diffusion coefficient D for root penetration in the profile in Equation ( 24) is assumed to depend on crop development stage (with zero elongation at later stages) as well as soil temperature and soil strength.Similar to the approach proposed by Gaiser et al. (2013), we modify f s to reflect the ability of roots to avoid mechanically strong soil layers by growing downwards through soil macropores.For the sake of simplicity, we assume that this ability of roots depends on the connected macropore fraction f con and that soil strength does not limit downward root penetration in the soil profile when macropore connectivity is at a maximum: where Γ d3 is a binary variable to indicate that the development stage is earlier than S d3 (Γ d3 = 1, zero otherwise), and D max (m 2 d À1 ) is the maximum value of the root diffusion coefficient.Equations ( 28) and ( 29) show that although macropores may allow roots to penetrate unimpeded by mechanical resistance downwards in the profile, it is currently assumed in the USSF model that they do not alleviate reductions in the below-ground demand for assimilates, and thus, the potential effects on root growth in the surrounding strong soil matrix (Colombi et al., 2018).

Storage and turnover of soil organic matter
Soil structure is assumed to influence organic matter storage and turnover in the USSF model in two distinct ways.Firstly, the mineralization rate of organic matter in microporous regions is slower due to physical protection (e.g.Dungait et al., 2012;Ekschmitt et al., 2008;Kravchenko & Guber, 2017).Secondly, in an attempt to mimic in a simple way the spatial distribution of root proliferation in soil, the partitioning of organic matter derived from both root decay/die-off and root exudation between micropore and mesopore regions of the soil is determined by their relative volumes.
Dual-porosity model of soil organic matter storage and turnover.The dual-porosity version of the ICBM model (Andrén & Kätterer, 1997;Wutzler & Reichstein, 2013) proposed by Meurer, Chenu, et al. (2020) and further developed here to account for the effects of soil temperature, soil moisture, and microbial biomass on organic matter decomposition rates (see supplementary information) is used to track four pools of soil organic matter (kg SOM m À2 ) representing two contrasting qualities of SOM stored in the mesopore and micropore regions of the matrix.This model considers two pools of young undecomposed organic matter, one stored in parts of the soil in close contact with mesopores and the other stored in microporous regions (M Y(mes) and M Y(mic) respectively).Likewise, the model also accounts for pools of older microbially processed organic matter stored in mesoporous and microporous regions of soil (M O(mes) and M O(mic) ).The total storage in a soil layer (kg SOM m À2 ) and changes in the mass of SOM in the four pools are given by: where f r(mic) is the fraction of the root-derived OM supplied to the micropore region (see supplementary information), f r (À) is the fraction of the total root biomass in the soil layer in question, k Y and k O (d À1 ) are first-order rate constants for the decomposition of younger and older organic matter, k t and k w (À) are response functions for the effects of soil temperature and moisture on decomposition rates, k u is an uptake limitation factor (À) varying from zero to unity that implicitly accounts for the dynamics of microbial biomass (Wutzler & Reichstein, 2013; see supplementary information), F prot (À) is a response factor varying from zero to unity that reduces decomposition rates in micropores to reflect physical protection, ε (À) is the microbial efficiency that varies from zero to unity, and T Y and T O (kg m À2 d À1 ) are source-sink terms for the mixing of organic matter between pore regions by both tillage and earthworm bioturbation (for a derivation, see Appendix S1 in the supplementary information): where τ d (m 3 m À3 d À1 ) is the soil mixing rate due to tillage.The USSF model also assumes that the organic matter stored in each of the four pools becomes mixed by tillage vertically within the depth of tillage on the day of cultivation.The efficiency of this vertical mixing by tillage is specified by the model user.Different values can be specified for harrowing and ploughing.
Supply of organic matter from crop residues and organic amendments.The supply of organic matter from incorporated above-ground harvest residues and organic amendments (Equations ( 15) and ( 31)), I a , (kg m À2 d À1 ) is given by: where I m is the amount of manure applied (kg m À2 ), Γ m is a binary variable indicating whether manuring occurs or not (one or zero), B ag (kg m À2 ) is the above-ground harvest residues, and f inc is the proportion of the supply of organic matter from above-ground sources that enters the soil layer in question.This depth distribution in the soil is determined by a user-specified proportion added to the surface layer, whilst assuming that the remainder is uniformly distributed in the soil to a maximum depth of incorporation.In contrast, root exudates are supplied to the soil proportionally with the distribution of root biomass with depth, as a fixed proportion of assimilates allocated below ground (see Equations ( 31) and ( 33)).It is further assumed in the USSF model that organic amendments and above-ground crop residues are added solely to the young organic matter pool in the mesopore region (Equations 31 and 33).

Soil hydraulic functions
In USSF, soil water flow is calculated with the Richardson-Richards' equation (see equation S1 in the supplementary information) with soil hydraulic properties that vary in time as a consequence of the changes in porosity and pore size distribution outlined above.The soil water retention and hydraulic conductivity functions required to solve this equation are obtained by combining the Mualem/Brooks-Corey model (Brooks & Corey, 1964;Mualem, 1976) for the soil matrix (i.e., micropores and mesopores) with the conceptually compatible model for the macropore region proposed by Jarvis (2008).Assuming that the residual water content is negligible, the composite soil water retention function in USSF is given by: where ψ max is the pressure head (cm) equivalent to the largest diameter macropore and λ mat and λ mac (À) are indices reflecting pore size distributions in the matrix and macropore regions.
Assuming that tortuosity is independent of pore size in the macropore region (Jarvis, 2008), hydraulic conductivity K (m d À1 ) is given by: where τ (À) is a parameter that reflects the tortuosity and connectivity of the matrix pore space and K s(mat) and K s (mac) are the saturated hydraulic conductivities of the soil matrix and macropores, respectively (m d À1 ), given by (Jarvis, 2008;Jarvis et al., 1999;Nasta et al., 2013): where C mat and C mac (m 3 day À1 ) are composite constants that can be derived in principle from the underlying physical theory, which is based on a conceptual model of the soil as a bundle of capillaries.However, pore networks in soil are more complex than this model allows (Hunt et al., 2013), so in practice C mat and C mac should be treated as empirical coefficients.
For the sake of simplicity, the pressure head equivalent to the largest macropore in soil ψ max is assumed to be constant, while the pore size distribution index in the macropores, λ mac is estimated from the soil clay content and macroporosity using the pedotransfer function derived by Jarvis (2008) from measurements of nearsaturated hydraulic properties: Temporal variations in the matrix porosity and its distribution among the two pore size classes will alter the pore size distribution index in the Brooks-Corey equation and, therefore, also the hydraulic conductivity function (Equations ( 40) and ( 42)).From Equation (38), we have: pore size distribution.In this example, a reduction in the organic matter content of a degraded soil has resulted in a loss of aggregation porosity in the matrix of 0.15 m 3 m À3 (0.1 and 0.05 m 3 m À3 in the mesopore and micropore regions respectively; Table 1), which consequently decreases λ mat (from 0.116 to 0.056), as the proportion of the matrix pore space comprising mesoporosity is smaller (see Equation 45).The hydraulic conductivity at matrix saturation is nearly three times smaller in the degraded soil as a result of these changes in the matrix porosity and its pore size distribution (see Equation 42; Figure 3b).The degraded soil is also characterized by a complete loss of macroporosity (Table 1).As a consequence, the total saturated hydraulic conductivity is nearly two orders of magnitude smaller than in the structured soil (see Figure 3b).

| Scenario simulations
We present here the results of simulations illustrating the potential impacts of soil structure dynamics on soil water balance, organic carbon sequestration in the soil, root growth and crop yields.The simulations were run using historical weather data (daily rainfall, maximum and minimum air temperatures, solar radiation) recorded at Ultuna, near Uppsala in east-central Sweden for the 30-year period of 1985 to 2014.The annual precipitation recorded during this period at Ultuna varied between 412 (in 1989) and 683 mm (in 2012), with a 30-year average of 541 mm.Daily rainfall totals were assumed to fall in a single storm at an intensity estimated using the rainfall disaggregation method described by Olsson (1998) with the parameters obtained for a site in south Sweden by Güntner et al. (2001).Potential evapotranspiration was calculated internally in USSF using the Hargreaves equation (Hargreaves & Samani, 1982).A soil profile was simulated which was initially 112 cm thick, using a unit hydraulic gradient as the bottom boundary condition in Richardson-Richards' equation.We simulated continuous winter wheat sown on 6th September and harvested on 18th August each year, with conventional tillage in autumn (ploughing and harrowing to 30 and 6 cm depth respectively).
We considered two contrasting soils.The first soil (hereafter termed "clay") is based on field measurements and calibration of the USSF model to two different experiments carried out at Ultuna in a clay soil (clay contents of ca.53 and 62% in topsoil and subsoil).We used crop parameters for winter wheat derived from a calibration of the USSF model against detailed field measurements of soil water contents and above-and below-ground biomass made during one growing season for a winter wheat crop (in preparation).The parameters determining organic matter turnover in the model were taken from Meurer, Chenu, et al. (2020), who calibrated the dualporosity SOM model against data from a long-term plot experiment with manuring and bare fallow treatments located elsewhere at Ultuna, but on the same type of clay soil.In the model applications shown here, we switched off the response functions for temperature, moisture and microbial limitation in the SOM model (i.e., k t = k w = k u = 1 in Equations 31-34) in order to maintain compatibility with these previous model calibrations.Model parameters C mat and C mac regulating the dynamics of soil hydraulic conductivity were set based on data available for the Ultuna clay soil (Messing & Jarvis, 1990, 1993).Based on site data reported by Lagerlöf et al. (2012), we assumed that anecic earthworms were absent in the Ultuna clay soil (i.e., τ a in Equation ( 15) was set to zero), while the rate coefficient for endogeic earthworms τ e,max in Equation ( 14) was estimated at 7 Â 10 À5 d À1 from a measured earthworm biomass of 16 g m À2 (Lagerlöf et al., 2012), known soil bulk densities in the Ultuna clay soil, and an assumed soil ingestion rate of 1 g g À1 d À1 (Curry & Schmidt, 2007).The slope of the shrinkage characteristic was estimated from measurements made on soil aggregates sampled from Ultuna clay at a depth of 40 to 50 cm (Messing & Jarvis, 1990).We simulated an initially degraded soil with an organic matter content of 1.8% in the ploughed horizon, 1% at 30-70 cm depth, and 0.5% at depths below 70 cm and no bioporosity throughout the soil profile.
We also performed scenario simulations for a hypothetical loam soil with a clay content of 15% throughout the soil profile.We assumed a lack of swelling and shrinking in this soil (p = 0).The reference porosity, which was set to 0.54 m 3 m À3 in Ultuna clay to represent measured bulk densities attained at a fully swollen state, was fixed at 0.35 m 3 m À3 in the loam.We also introduced bioporosity generation by anecic earthworms in the loam (Equations ( 13) and ( 15)), with initial values of bioporosity set to zero in the topsoil and 0.02 m 3 m À3 in the subsoil to reflect their presence.All other parameters were set to the same values as in the clay.For both soil types, we ran simulations with and without considering soil structure dynamics.In all four scenarios, manure incorporation in the ploughed horizon was simulated at a rate of 1 kg m À2 year À1 .In the case without soil structure dynamics, the soil physical and hydraulic properties remained constant during the simulation at their initial values and were thus unaffected by both the biological (root decay, earthworms, organic matter turnover) and physical (tillage, swell-shrink) structure-forming processes included in the USSF model.
The complete parameterization of the two soil scenarios can be found in the supplementary information in the documentation file "USSF model parameterization.xlsx".The values of sixteen parameters regulating soil structure dynamics included in a sensitivity analysis described in the following section are also given in Table 2 (termed "nominal values").

| Sensitivity analysis
A full sensitivity analysis including all the parameters in the USSF model would clearly be a formidable task and is well beyond the scope of this study.Instead, based on the two soil scenarios described above, we performed a sensitivity analysis to gain some insights into the relative importance of sixteen parameters in USSF that regulate soil structure dynamics through both physical and biological processes (see Table 2).As USSF is comparatively slow to run, we employed the Morris (or "elementary effects") method to quantify the sensitivity of various target output variables to these model parameters.This method is less computationally demanding than a complete Monte Carlo approach, yet avoids some of the limitations of the simple "one-at-a-time" method (Morris, 1991;Saltelli & Annoni, 2010).
In the Morris method, each parameter X i , i = 1, …., k, is normalized and assumed to be uniformly distributed between minimum and maximum values, with values set at j selected levels or steps, i.e., {0, 1/( j-1), 2/( j-1), …, 1}.The nominal and limiting values for the parameters shown in Table 2 were determined, where possible, from prior experience of the expected range of variation.For new parameters where such experience was lacking, the nominal values and ranges were set to achieve reasonable results in preliminary simulations.In the Morris method, T A B L E 2 Values used in the sensitivity analyses for parameters related to soil structure dynamics.

Group
Parameter Only included in the analysis for the clay scenario (for the loam scenario p is fixed at zero).b Only included in the loam scenario (f bio is set to zero in the clay scenario).
a number of so-called "trajectories", n, are run, each one comprising k + 1 simulations for which parameter values are randomly selected to vary one at a time from randomized starting values, while all the others remain fixed.The total number of simulations to be performed therefore equals n(k + 1).The elementary effects E i (X) are defined by: where Y is the target output variable of interest and Δ is a normalized step size for changes in parameter values.
Parameter sensitivity is usually quantified using the means and standard deviations of the elementary effects.
For monotonic responses, a large mean value of the F I G U R E 4 Simulated temporal variations in hydraulic properties of the clay soil (available water capacity is defined here as the difference in water contents between pressure heads of À0.5 and À 150 m).
elementary effects means that the target variable is sensitive to the parameter in question.However, calculated mean elementary effects can be affected by cancellation effects if the response of the target variable to the parameter is non-monotonic.We therefore make use of the mean absolute elementary effect to quantify sensitivity.Elementary effects do not take into account parameter correlation directly, but a large standard deviation indicates that the parameter interacts with other parameters or that its effects are non-linear.
For the parameters j and Δ it is advantageous to set j to an even value and Δ equal to j/(2(j-1)).Although this does not guarantee equal-probability sampling of the inputs, it does improve the likelihood of this (Morris, 1991).We used j = 4 so that Δ ¼ 0:6 _ 6. Saltelli and Annoni (2010) suggest that reliable qualitative F I G U R E 5 Simulated temporal variations in hydraulic properties of the loam soil (available water capacity is defined here as the difference in water contents between pressure heads of À0.5 and À150 m).results for screening purposes can be obtained with the Morris method for analyses supported by between four and ten trajectories.The number of trajectories should also balance the number of steps or levels in order to obtain a reasonably representative exploratory sample.We used eight trajectories giving 136 simulations.
Parameter sensitivity was assessed for six target output variables: four terms of the water balance (accumulated recharge, transpiration, evapotranspiration, and surface runoff plus interflow for the 30-year simulation period), harvest (grain) yields and the change of organic matter stocks during the simulation.Note that water lost by interflow (see supplementary information) was a small component of the water balance in all simulations, which is why it was assessed together with surface runoff.These output variables were normalized by their mean values to enable us to compare their relative sensitivity to the input parameters.The analyses were carried out using the Sensitivity Analysis Library in Python (SALib).and 5 show some selected physical and hydraulic properties for the simulations considering soil structure dynamics in the clay and loam soils, respectively.Seasonal variations in macroporosity are simulated in both soils related to tillage and subsequent consolidation in the ploughed topsoil as well as biopore generation due to root decay at harvest in the upper subsoil.In the clay soil, seasonal variations in the crack porosity also occur due to soil swelling and shrinkage (Figure 4).The hydraulic properties simulated for the soil matrix also vary on seasonal time-scales, as a consequence of the aggregation porosity formed due to the input of crop F I G U R E 6 Simulated soil organic matter stocks in the topsoil and total soil profile, with soil structure dynamics modelled using the nominal parameter values shown in Table 2.
F I G U R E 7 Water balances simulated for the clay and loam scenarios, accounting for soil structure dynamics using the nominal parameter values shown in Table 2.
residue OM to the soil and, in the clay soil, swelling and shrinkage of the soil matrix.Some significant longer-term trends are superimposed on these seasonal fluctuations.In the ploughed topsoil, total porosity steadily increases throughout the 30-year period in both soils as a consequence of increases in the stock of organic matter (see Figure 6) and, therefore, aggregation porosity in the soil matrix.Smaller increases in matrix (aggregation) porosity and SOM stocks were also simulated in the loam subsoil (Figure 5), whereas these effects were less apparent in the clay subsoil (Figures 4 and 6) probably due to the absence of anecic earthworms.The organic amendments applied in the model also influenced the shape of the water retention function, both in the topsoil of the clay and throughout the loam soil profile, with the Brooks-Corey pore size distribution index λ mat increasing according to Equation ( 45) because the mesoporous region gained comparatively more organic matter and, therefore, aggregation pore space (not shown).Figures 4 and 5 show that these changes in the shapes of the soil water function resulted in only modest increases in the available water capacity in the loam and clay topsoil.The saturated matrix hydraulic conductivity steadily increased in the plouged horizons of both soils and in the upper subsoil of the loam, as a result of the increases in matrix porosity and pore size distribution index (see Equation ( 42)).The long-term trend in the total saturated hydraulic conductivity, K s (i.e., K s = K s(mac) + K s(mat) ) in the upper subsoil of the loam is even more dramatic, increasing by one order of magnitude during the 30-year period, mostly as a consequence of the enhanced generation of root and earthworm biopores due to improved crop growth and carbon inputs to the soil.In contrast, the clay soil showed significant seasonal variations due to tillage and soil shrinkage and cracking, whereas long-T A B L E 3 Dry matter balance (gains and losses in kg m À2 year À1 ) for the 30-year simulations.
Winter wheat grain yields simulated for the clay and loam scenarios, accounting for soil structure dynamics using the nominal parameter values shown in Table 2.
term trends in K s during the simulation were much less apparent.No long-term improvements in the physical and hydraulic properties were simulated in the subsoil (see Figures 4 and 5) at depths beyond those influenced by crop roots and earthworms.The trends in the lower subsoil shown in Figures 4 and 5 were generated by a partial loss of the macroporosity assumed to be initially present and small changes in organic matter stocks.We have no data with which to validate the simulations for the loam soil, since this is entirely a virtual scenario.However, a qualitative "reality check" of the simulations of dynamic hydraulic properties for the clay soil is feasible, since it is based on the Ultuna field site for which occasional measurements of soil hydraulic properties have been made.Messing and Jarvis (1990) reported seasonal variations of both macroporosity and saturated hydraulic conductivity K s in the topsoil and upper subsoil in Ultuna clay of a similar range and magnitude as those shown in Figure 4, with K s varying from ca. 100 cm day À1 in spring and autumn to 1600 cm day À1 in summer, corresponding to increases in ϕ mac from ca. 0.06 to ca. 0.2 m 3 m À3 .
Figure 7 shows predicted water balances for both soil scenarios for the 30-year simulations that accounted for soil structure dynamics.The water balances show similar interannual variations in both soils, except that transpiration is ca.25% larger in the clay, while recharge to groundwater and losses by surface runoff and interflow are accordingly smaller (Figure 7).These differences are presumably a consequence of soil shrinkage and cracking in the clay, which reduced surface runoff and allowed roots to penetrate deeper into the soil profile.Figure 8 shows the yields of winter wheat, again for the simulations accounting for structure dynamics.Yields are, on average, almost 40% larger in the clay soil as a consequence of the larger transpiration rates (mean yields of 0.72 kg m À2 in the clay vs. 0.52 kg m À2 in the loam; Table 3).Winter wheat grain yields in this range are rather typical for the Uppsala region (the average yield for the county of Uppsala in this period was 0.53 kg m À2 , SCB, Statistics Sweden).F I G U R E 1 0 Differences in annual recharge and surface runoff for simulations with and without soil structure dynamics.

| Influence of soil structure dynamics
No clear long-term trends are apparent in the simulated water balance terms and crop yields for both soils shown in Figures 7 and 8, despite large increases in the simulated stocks of soil OM due to the organic amendments and the consequent improvements in soil hydraulic properties (Figures 4 and 5).This can be largely attributed to a series of drier summers during the later years of the simulation period (not shown).However, clear long-term effects of the recovery of soil structure are revealed by a comparison of the results of the simulations run with and without considering structure dynamics.This is especially true for the loam soil, for which small reductions in surface runoff and ca.10% increases in both transpiration and grain yields are simulated during the 30-year period in the simulation with structure dynamics, compared with the case for a static soil structure (Figures 9 and 10).In the first two years, grain yields on the loam were larger in the simulation with static soil structure, because the initial soil physical properties, which were maintained throughout the simulation, were more favourable for crop growth.In the clay soil, crop transpiration and grain yields for the 30-year period are ca.26% and 40% larger in the simulation accounting for structure dynamics (Figure 9).These differences are predominantly a consequence of the within-season dynamics of soil structure, with the development of shrinkage cracks which reduce surface runoff and enhance root penetration in the soil profile.In contrast to the loam, long-term trends in transpiration and grain yields were hardly discernable in the clay soil, in line with the more modest changes in hydraulic properties during the simulation, especially in the subsoil (Figures 4 and 5).
Figure 11 contrasts root distributions simulated in the soil profile at mid-summer of the final year for simulations with and without a consideration of soil structure dynamics.For both soils, the root biomass was larger in the simulations that account for soil structure dynamics.The difference between these scenarios is especially clear in the deeper subsoil of the loam, where improved soil physical conditions (i.e., smaller bulk density and larger macroporosity due to the presence of anecic earthworms) led to a greater proliferation of roots (Figure 11).
For both soils, accounting for soil structure dynamics in the model also led to increases in simulated OM stocks in the soil profile that are roughly proportional to the additional amounts of OM supplied to the soil (ca.16% and 2% in the clay and loam soils, respectively; Table 3).

| Sensitivity analysis
Figure 12 shows the calculated mean absolute elementary effects for the parameters included in the Morris sensitivity analysis for the clay and loam soils, respectively, while the means and standard deviations of the elementary effects are presented in Figures S2-S5 in the supplementary information.These results suggest that the inputoutput relationships in USSF are mostly monotonic, while none of the sixteen parameters considered seems markedly more affected by non-linear and/or interaction effects than any of the others.
Figure 12 shows that the parameter sensitivities are broadly similar on both soils.For example, three parameters related to the effects of the crop on soil sealing (h min , h max , C max ) are by far the least sensitive for all six target variables on both soils.The aggregation factor, f agg , is clearly the most sensitive parameter for all target output variables, apart from surface runoff on the clay soil, for which it is the second most sensitive parameter.Figures S2 and S3 show that, as expected, improved soil aggregation reduces surface runoff and recharge, whilst increasing evapotranspiration, grain yields and the sequestration of organic matter in soil.However, some differences between the soils are also apparent.For example, in the loam, for which swell/shrink is neglected, the parameters regulating bioporosity generation and loss are comparatively more sensitive, especially for transpiration, yields, and organic matter stocks.As might be expected, all three of these target variables show similar responses to the sixteen parameters, since along with radiation and temperature, water stress controls assimilation, while changes in soil organic matter stocks are regulated by the supply of crop residues.
Of the six output variables considered, surface runoff is clearly the most sensitive to the parameters controlling structure dynamics in both soils.With the exception of h min , h max and C max , all of the parameters showed some sensitivity towards this target variable.As expected, surface runoff is moderately sensitive to the parameters controlling loosening by tillage and subsequent consolidation and sealing (i.e., τ till , k con , k s(max) , f som(c,s) , f som(n,s) ), but it is also equally sensitive or even more sensitive to parameters regulating aggregation and earthworm activity.
The analysis revealed that the most sensitive parameters controlling soil structure dynamics are those determining aggregation and swell/shrink.Fortunately, these parameters are also relatively straightforward to estimate from field measurements.The aggregation factor, f agg , can be estimated from the relationship between soil organic matter content and bulk density (e.g., Meurer, Chenu, et al., 2020), while the slope of the shrinkage characteristic, p, can be directly measured, although this is not a routine analysis in most soil laboratories.Much of the soil parameterization in the USSF model is already based on simple pedotransfer functions, so that the development of an additional pedotransfer function for p would be useful.We believe that sufficient data to support this should be available in the literature.
F I G U R E 1 2 Mean absolute elementary effects and 95% confidence intervals calculated for sixteen parameters regulating soil structure dynamics in the two soil scenarios (see Table 2 for explanations of parameter names).Some moderately sensitive parameters in USSF are also highly uncertain, with very little information currently available concerning typical values.This is especially the case for the parameters controlling the generation and loss of bioporosity, which is described in USSF with a phenomenological approach using simple proxy variables for faunal biomass and activity.The capability of USSF to reliably predict these processes would likely be enhanced by new and improved estimation schemes for the relevant parameters combined with additional testing of the model against suitable datasets (e.g., Keller et al., 2017;Lucas et al., 2019).Clearly, it would also be worthwhile to test other novel components of the model, in particular the description of root growth and how root penetration is affected by soil strength and structure, using existing data sets that are suitable for this purpose (e.g.Colombi et al., 2018).

| CONCLUSIONS AND OUTLOOK
Here, we described a new soil-crop model (USSF) that accounts for soil structure dynamics mediated by both biological and physical processes.The model is designed primarily as an explorative tool to better understand the key controlling factors and relevant time-scales of soil degradation and restoration driven by changes in soil and crop management and climate.
In this first study, we performed scenario simulations and sensitivity analyses to illustrate the effects of soil structure dynamics on soil physical and hydraulic properties, soil water balance, crop growth and stocks of soil organic matter, using as a case study the application of manure to restore physically degraded soils that are depleted of organic matter.As a consequence of the build-up of soil organic matter stocks, the model simulated what appear to be reasonably realistic long-term trends in bulk density and porosity, hydraulic conductivity, water retention and root growth, all of which are superimposed on seasonal patterns induced by tillage and, in the clay soil, swelling and shrinking.Each of these individual simulated effects was relatively small, but taken together they steadily reduced surface runoff and enhanced transpiration rates and grain yields during the 30-year simulation period, especially in the loam soil.These changes in the water balance induced by structure dynamics would also be important for predictive modelling of the environmental impacts of agriculture, for example, of erosion or pollution of groundwater and surface water bodies by nutrients and agro-chemicals.
Although some of the key components of the model have been successfully tested in previous studies, there is clearly a need for further model testing for a wider range of soil types and climates.Furthermore, this paper describes a first version of the USSF model that neglects some potentially important processes affecting soil structure, such as the effects of traffic compaction and soil freezing and thawing, while others are treated in a very simplistic manner, such as soil consolidation after tillage and the dynamics of soil bioporosity.In a wider perspective, nutrient cycling in the soil-plant system has so far also been neglected in USSF, which means that low-input farming systems cannot yet be considered.Clearly, there is some scope for further model development.Nevertheless, we conclude that the USSF model is a promising new tool to investigate a wide range of processes and phenomena triggered by land use and climate change that other soil-crop models cannot simulate.Results from this study show that feedbacks in the soil-crop system mediated by the dynamics of soil physical and hydraulic properties, are potentially of central importance for longterm predictions of soil water balance, crop production and carbon sequestration under global change.

Figure
Figure 3a,b presents an illustrative example of how the soil water retention and hydraulic conductivity functions in USSF would respond to changes in porosity and T A B L E 1 Porosities used to calculate the water retention and hydraulic conductivity functions shown in Figure3, with common parameter values of ψ max = À1 cm, ψ mes/mac = À10 cm, ψ mic/ mes = À600 cm, λ mac = À0.5, τ = 1 and C mat and C mac = 0.01 m 3 d À1 .Soil m 3 m À3 Figures4 and 5show some selected physical and hydraulic properties for the simulations considering soil structure dynamics in the clay and loam soils, respectively.Seasonal variations in macroporosity are simulated in both soils related to tillage and subsequent consolidation in the ploughed topsoil as well as biopore generation due to root decay at harvest in the upper subsoil.In the clay soil, seasonal variations in the crack porosity also occur due to soil swelling and shrinkage (Figure4).The hydraulic properties simulated for the soil matrix also vary on seasonal time-scales, as a consequence of the aggregation porosity formed due to the input of crop

F
I G U R E 9 Differences in annual transpiration and grain yields for simulations with and without soil structure dynamics (% difference = 100*(with-without)/without).

F
I G U R E 1 1 Example of winter wheat root distributions for simulations with and without structure dynamics (root distributions simulated at mid-summer in the final year of the 30-year simulation).