^{1}

^{2}

^{1}

^{1}

^{*}

^{1}

^{2}

Edited by: Philippe C. Baveye, AgroParisTech Institut des Sciences et Industries du Vivant et de L'environnement, France

Reviewed by: Wulong Hu, Wuhan University of Technology, China; Albert J. Valocchi, University of Illinois at Urbana-Champaign, United States

This article was submitted to Soil Processes, a section of the journal Frontiers in Environmental Science

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Structure formation and self organization in soils determine soil functions and regulate soil processes. Mathematically based modeling can facilitate the understanding of organizing mechanisms at different scales, provided that the major driving forces are taken into account. In this research we present an extension of the mechanistic model for transport, biomass development and solid restructuring that was proposed in a former publication of the authors. Three main extensions are implemented. First, arbitrary shapes for the building units (e.g., spherical, needle-like, or platy particles), and also their compositions are incorporated into the model. Second, a gas phase is included in addition to solid, biofilm, and fluid phases. Interaction rules within and between the phases are prescribed using a cellular automaton method (CAM) and a system of partial differential equations (PDEs). These result in a structural self organization of the respective phases which define the time-dependent composition of the computational domain. Within the non-solid phases, chemical species may diffuse and react. In particular a kinetic Langmuir isotherm for heterogeneous surface reactions and a Henry condition for the transfer from/into the gas phase are applied. As third important model extension charges and charge conservation laws are included into the model for both the solid phase and ions in solution, as electrostatic attraction is a major driving force for aggregation. The ions move obeying the Nernst-Planck equations. A fully implicit local discontinuous Galerkin (LDG) method is applied to solve the resulting equation systems. The operational, comprehensive model allows to study structure formation as a function of the size and shape of the solid particles. Moreover, the effect of attraction and repulsion by charges is thoroughly discussed. The presented model is a first step to capture various aspects of structure formation and self organization in soils, it is a process-based tool to study the interplay of relevant mechanisms

Understanding the structural formation in soils with respect to space and time is highly demanding at any scale. Pedogenesis induces an aggregate structure in most soils (see Totsche et al.,

However, to our knowledge no currently available model is able to simulate a fully dynamic and mechanistic evolution of the soil structure. A first step towards that direction applying cellular automaton methods (CAMs) has been made in the work of Crawford et al. (

CAMs provide a flexible way to describe the restructuring of soil particles and pore-filling phases. The variability in soil structure as a consequence of the self organization of the soil-microbe system has been investigated in Crawford et al. (

CAMs have also been applied in the context of biofilm models. In combination with experiments, Tang and coworkers prescribed biomass spreading rules in Tang and Valocchi (

In Ray et al. (

In this research, we present the related model extensions and study their impact using simulation scenarios. The CAM rules are adapted for the effects of solid surface charges and ion transport, and their impact is considered also in the solute. Rotations of solid particles are for instance included into the model to facilitate the aggregation of solid particles with opposite face charges.

To study the formation of microaggregates from building units of diverse types new prototypes have been implemented to represent different geometric structures. Moreover, a gas phase has been included to consider moisture variations. Thus situations can be taken into account where the aggregates are not fully saturated but, e.g., coated only by a thin film of water. A restructuring of the solid phase may then induce the necessity of a reorganization of the gas phase to maintain the non-wetting properties (section 2.3.2). Exchange between and transport within the different phases may become prominent since mobility of species varies within water, gas, or biofilm. As electric forces are an important driving force for aggregation—as highlighted in Totsche et al. (

Several simulation scenarios are investigated in two dimensions to demonstrate the effects of the novel mechanisms. The underlying model equations are discretized by means of a local discontinuous Galerkin (LDG) method as presented and analyzed in Rupp and Knabner (

The paper is structured as follows: In section 2, we establish the underlying mathematical model. We shortly review the model introduced in Ray et al. (

Within our model, we essentially consider the following prototypical time- and space-dependent model parts:

a solid phase (

a fluid phase (

a bio phase (

a gas phase (

several types of mobile, potentially charged chemical species that may participate in heterogeneous or homogeneous reactions. Species that are present in the gas and the fluid/bio phase follow Henry's Law as transfer condition. Examples for relevant chemical species are for instance oxygen, or nitrate as nutrients for bacteria and microorganisms creating extracellular polymeric substances (EPS) and biomass.

a gluing agent (e.g., EPS) which is produced as a consequence of local biological activity (c.f. Crawford et al.,

The cellular automaton acts on a regular quadrilateral domain ^{2} squares ^{i} with ^{i}. At first, one of the following cell states is (randomly, but in desired proportions) assigned to each of the squares: “bio” (

Domain

The cells ^{i} in the cellular automaton correspond to the smallest physical units in the model. In our new approach different inseparable building units with various shapes may be defined that are composed of cells ^{i} (c.f. Figure

Possible shapes of building units (size not to scale); from left to right: Spherical, platy, or needle-like, and composites. Black and white to indicate equal charge for the composites.

The union _{s} with boundary Γ: = ∂_{s}. Likewise, the union _{b}, the union _{f}, and the union _{g} (c.f. Figure _{LG}: = (∂(_{f} ∪ _{b})) ∩ ∂ _{g} of the gas with the liquid phase. In each

Within the fluid, bio, and gas phases, the continuum parts of the model come into play. Here, (possibly coupled, partial) differential equations are solved for the transported, potentially charged chemical species, and also the immobile biomass. Likewise, an ordinary differential equation is considered for the

The cellular automaton rules are based on the movement of the building units (section 2.2). For single cells this was already described in Ray et al. (

Stencil of size 0 (

The application of the cellular automaton method was already discussed in Ray et al. (

The biomass spreading rule is based on the CAM described in Tang and Valocchi (

The CAM rules for uncharged solid cells are inspired by Crawford et al. (

For the new CAM, “movable cells”

Here, _{α, îk} [1/L] denotes the average concentration of the gluing agent at the common face of the neighboring solid cells ^{2}] balances the impact of gluing and electric forces. If (1) holds true for all neighbors of

The attraction or affinity _{i} [−] of potential target fluid or gas cells ^{i} is then a weighted superposition of the different attracting forces related to gluing agent, type of neighbors, and electric forces:

with proportional constants ^{2}], and

Effects of charges in solid restructuring. Initial configuration with potential movements to target cells for a stencil of 2

As in the case without charges, the target cell with the largest _{i} is selected. If several target cells are identified as equally attractive, one is randomly chosen and the restructuring is carried out. A conflict may occur if the same target cell is selected by different mobile solid cells. To resolve such conflicts one solid cell is randomly chosen to move for each of the conflicts.

In Figure

In the following, we briefly discuss the extensions of the CAM due to the occurrence of a gas phase.

In natural soils the

Reorganization of gas cells being tangent to solid cells.

In Figure

CAM including the reorganization of the gas phase illustrating its non-wetting property. Random initial configuration

The nonlinearly coupled Nernst-Planck-Poisson Equations (2) describe the movement of a potentially charged chemical species in the liquid, i.e., the time-dependent domain _{f} ∪ _{b} whose structural arrangement is defined in each time step by the CAM (see section 2.3). In contrast to this, the Poisson Equation (4) is defined on the whole periodic domain ^{2}] denotes the volume concentrations of the _{s}.

Here, the advective flux is proportional to the ^{3})] with discontinuous, phase-dependent mobility _{r} [IT^{2}/M]. Likewise, the diffusivities _{r} [L^{2}/T] are discontinuous. The _{r} [1/(L^{2} T)] and obey the mass action law while the

Since the mobility _{r} vanishes for uncharged species

with the inverse

The electric field E and the ^{2}/(T^{3}I)] are computed in terms of a

where ^{2}] and ^{-1}]. Finally, _{r} [−] is the _{0} [I^{2}T^{4}/(ML^{2})] denotes the _{r} [−] is the discontinuous, phase depending

The consistency condition

complements the model meaning that sources and sinks of the electric potential cancel out globally, i.e., the whole domain must be electrically neutral. Hence, _{r} cannot be chosen randomly but has to ensure the conservation of charge, while

The equations of the continuum model part (2), (3), (4) are discretized using the local discontinuous Galerkin method as described in Rupp and Knabner (

The overall algorithm including the continuum model part and the discrete CAM is depicted in Table _{1} [−], ϵ_{2} [−] are constants. The frequency of updates of the geometric structure (defined by

Within the discrete model part—which is solved for all global time steps—the main factors related to structural changes, apart from the generation of biomass from bacteria, are evaluated, namely the biomass spreading, solid restructuring and reorganization of the gas phase. The implementation is written in C++ and based on M++ (Wieners,

To demonstrate the impact of the implemented mechanisms on the formation of structures in soils, we present several simulation scenarios. The first scenarios focus on the demonstration of single effects and are chosen in such a way that the self organization according to CAM rules is illustrated. To this end, various model components such as charges, biomass, gluing agent, or solutes in the liquid etc. are switched off. Thereafter, combined effects are investigated to illustrate the overall capability of the model. Finally, the interplay of the discrete and continuum model component is shown.

First, we investigate the influence of the range of attraction of particles on each other for structure formation which is represented by different stencil widths. Since this scenario focuses on the demonstration of a single effect, no charges, biomass, gluing agent, or solutes in the liquid etc. are taken into account. Thus, attraction of the cells to each other is uniform and can be interpreted as the sum of attracting forces as, e.g., van der Waals forces, that lead also to homoaggregation. It is thus only determined by the number of neighbors [see Equation (2) with

Self organization depending on range of attraction:

Since we study the formation of structures and do not consider their disaggregation here, the simulations run into a quasi-stationary state. The resulting aggregated structures are depicted in the middle of Figure ^{-1}], stencil of size 1: 29, 720/32, 768 ≈ 0.91 [L^{-1}], stencil of size 3: 9, 208/32, 768 ≈ 0.28 [L^{-1}] after 500 CAM steps; Figure

The shape and size of the inseparable building units strongly influence the formation of structural patterns. For an illustration of this effect, we consider domains of sizes 256 × 256 cells where 50% are filled by uncharged solid particles (either single cells or needles of length 5). In Figure ^{-1}]), and for the needles by a factor of 0.84 (from 46, 096/32, 768 ≈ 1.41 to 38, 950/32, 768 ≈ 1.19 [L^{-1}]).

Self organization depending on shape and size constraints:

Electrostatic forces are a major driving force for particle aggregation (c.f. Totsche et al.,

As initial configuration, we consider 20% solid particles and 64 × 64 cells in total for each case. For the second scenario, we randomly distribute constant charge numbers between −4 and +4 on all solid cell faces ^{−1}], for uncharged particles 444/819 ≈ 0.54 [L^{−1}], and for charged particles 1, 328/819 ≈ 1.62 [L^{−1}]; c.f. Figure

Effects of charges in solid restructuring: Initial configuration for both simulations, porosity

In this section, we combine the restructuring of phases according to the CAM with the PDE model to illustrate the capability of our comprehensive model: A three phase system (as in Figure ^{2}], respectively, is present initially. It is degraded with a first order rate and rate constant −0.375 [1/T] in the fluid phase representing, e.g., the consumption of a nutrient, such as oxygen by aerobic organisms in the fluid phase. The diffusion is faster in the gas phase (^{−4} [L^{2}/T]) compared to the fluid phase (^{−8} [L^{2}/T]) and the transfer between the phases is determined via Henry's law with solubility constant

Aerobic bacteria in fluid phase combined with solid restructuring and gas reorganization: The scale depicts the concentration of a nutrient. Thus, solid cells are red, fluid cells are violet to dark blue, and gas cells are pale blue to white in the right picture.

In this section, we combine the restructuring of phases according to the CAM with the PDE model for movement of potentially charged species in solution. Thus, we show the interplay of ions in solution that adsorb to charged particles, i.e., we combine the CAM for charged solids with the PDEs for a charged fluid phase. The sorption of ions to the surfaces defined by

changes the attraction of particles and the resulting structures. The left picture in Figure

Interplay of ions in solution with charged solids: Initial configuration (first image) of solid (red), quasi stationary, final configuration of charged solid in neutral solution (second), and final configuration of charged solid in ionic solution (third image) when heterogeneous reactions alter the total charges of solid cells' edges. The zoom highlights charges on solid edges, the rainbow scale corresponds to the surface concentrations.

We compare the resulting structures under the influence of an uncharged chemical species and positively/negatively charged species. The uniform initial concentration of the uncharged species is 10 [1/L^{2}] and homogeneous Neumann boundary is applied at the solid's surface Γ. Likewise, the uniform initial concentration of the positively charged chemical species is set to 10 [1/L^{2}]. To ensure electroneutrality (also balancing the total charge of the solid), a negatively charged species with an initially homogeneous concentration of 13.27 [1/L^{2}] is necessary (not plotted in Figure ^{2}] as depicted in the middle picture of Figure

In this research, we presented a comprehensive mathematical model for structure formation. A novel discrete–continuum approach was taken, combining a model of partial differential equations for charged, reactive multicomponent transport with a cellular automaton method for the interactive self-organization of solid, bio, fluid and gas phases. A thorough illustration by means of numerical simulations was performed. This versatile approach has the potential to study the interplay of different aggregation mechanisms

Although our results have already contributed towards enhancing our understanding of structure formation and self organization in soils, there are several aspects that may be added to the model.

First, more research is needed to investigate unsaturated fluid flow. This was for instance done for non-evolving angular pore networks representing soil aggregates in Ebrahimi and Or (

Furthermore other process mechanisms could be identified to broaden the applicability of our model. Those have to be determined and validated with the help of experimental studies which itself is a challenge. One step into that direction would be a model extension with respect to disaggregation of the resulting structures. Different aggregation and disaggregation mechanisms and their respective ratios need to be investigated. Moreover, building units and their composites (i.e., smaller units of building units) naturally undergo some random movement (representing self diffusion) independently of the CAM rules. Such a movement depending on the respective size of the particles needs to be included into the model. A related research question is how time scales are related to diffusion and how interaction processes can be balanced in a reasonable way.

Finally, a quantitative evaluation of the resulting structures would be possible by means of Minkowski functionals and further geometric measures characterizing among others the connectivity and compactness of a structure.

AR discretized and implemented the numerical model, and conducted the simulations. KT contributed to developing the modeling concepts and process mechanisms. AP contributed to the model development, the design of the simulation studies and the interpretation of the results. NR developed the model concepts and formulations, and the simulation scenarios.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This research was kindly supported by the DFG RU 2179 “MAD Soil—Microaggregates: Formation and turnover of the structural building blocks of soils”. Discussion within the RU significantly promoted the present research. We particularly appreciate the input of Thomas Ritschel (Friedrich-Schiller Universität Jena) and Stefan Dultz (Leibniz Universität Hannover), and thank our students Andreas Meier and Markus Musch.