Modeling root growth in granular soils: effects of root stiffness and packing fraction

We use molecular dynamics simulations to investigate the effects of root bending stiffness and packing fraction on the path followed by a growing root in 2D packings of grains representing a soil. The root is modeled as a chain of elements that can grow in length and change their direction depending on the forces exerted by soil grains. We show that the root shape is mainly controlled by the bending stiffness of its apex. At low stiffness, the root randomly explores the pore space whereas at sufficiently high stiffness, of the order of soil hardness multiplied by mean grain size, the root follows a straight path across the soil. Between these two limits, the root shape can be characterized by the standard deviation of its re-directions at the scale of soil grains. We find that this shape parameter varies as a power-law function of the normalized bending stiffness.


Introduction
Besides many other factors such as developmental instability on the cellular level [1], the mechanical interactions between a growing root and its surrounding soil can have major impact on the growth variability and plant biomass production [2][3][4][5][6]. In particular, the inhomogeneous structure of soil at the particle scale and broad distributions of contact forces can significantly affect root growth trajectories [7][8][9][10][11]. Continuum representation of soils and roots has been used for describing soil stability [12][13][14][15]. In such models, the effect of soil is incorporated mainly through its average resistance to the root ingrowth and the interactions between soil grains and roots are not considered. There are, on the other hand, functional-structural plant models [16] that describe plant growth and incorporate biological and environmental factors but do not account for the soil strength. Many experiments and observations at the grain scale have shown that grain configuration and local cracks or obstacles have crucial effect on root growth [17][18][19][20][21][22][23].
In this paper, we report on discrete modeling of soilroot interactions and its application to a single root growing in a granular soil. Our model accounts for the elastic behavior of the elongation zone of the root and for stiffening at higher zones where root maturation occurs. We consider granular packings prepared by means of different procedures in order to evaluate the effects of packing fraction. We present here only the results regarding root trajectories. We first briefly describe the numerical model. Then, we analyze root trajectories and the effects of root and soil parameters. We conclude with a short summary of the results and scopes of this work.

Root-soil model
The root is modeled as a chain of connected elongated elements; see Fig. 1. The first root element is fixed on the top of soil sample. Starting from the bottom end of this element, a new circular element is added and elongated by moving the circle along a direction at a given growth rate until a length r is reached. Then, a new element is added by linking the top end of the second element to the bottom end point of the first element, and the same process is iteratively applied to each new element. Two successive elements are linked together via linear normal and angular springs defined by their stiffnesses. For angular stiffness, the reference angle α 0 is 180 degrees so that when no forces or torques are exerted on the root elements, they simply align themselves along the direction of the first element. All angles α are allowed and the torque is given by where K b is the bending stiffness. The grain-grain and grain-root contacts are governed by elastic forces acting along normal and tangential directions, as well as by dry friction law with coefficients μ pp and μ pr , respectively. Viscous terms are also added to control energy restitution. We use molecular dynamics method with velocity-Verlet time-stepping scheme to integrate the equations of motion of both soil grains and root elements [24,25]. The simulation box is partitioned into cells to keep track of the neighborhood lists of grains and root elements. Soil samples were prepared by gravitational deposit of grains in the simulation box. The grain diameter distribution is uniform in volume fractions between a minimum d min and a maximum d max (these are outer grain diameters; see below) [26].
An important feature of roots is that the bending stiffness is not uniform along the root. The cap and elongation zones are soft, allowing the root to explore the available pore space. Above elongation zone, the bending stiffness increases with time as a result of the accretion of new material, change in internal material properties and radial growth. We do not take radial growth into account but we attribute a high bending stiffness 10 3 K b to the root elements that do not grow anymore, where K b is the bending stiffness at the joint between the growing root element and the one to which it is attached. In other words, the root cap can more easily change its direction when it touches soil grains as compared to other root elements.
Another feature is the reference angle α 0 , which is initially set to 180 degrees, but every time the bending stiffness at a joint is increased, the angle α 0 is changed to its current value at that joint. In this way, the root can deform only elastically and its shape reflects the path followed by its growing cap element, which keeps its value of α 0 = 180 degrees as long as it continues to grow. It is frozen to its current value when it stops growing and a new root element is created. We see that this model combines two aspects that are important for a growing root: 1) The root growth is an outcome of both a constant growth rate and interactions with soil grains, and 2) The root shape fully reflects ('memorizes') the growth dynamics.
One issue in 2D modeling is that the pore space is not connected. In 3D, the pore space is connected throughout the soil and the roots can spread inside the pore space. To fix this problem in 2D, we use two different grain diameters: 1) An outer diameter 2R p used for grain-grain contacts, and 2) An inner diameter 2R p for grain-root contacts, as illustrated in Fig. 1. Hence, the root elements can cross the gap, which has a thickness d g = 2(R p − R p ) when two grains touch. A crucial parameter for root growth is thus the ratio β = 2R r /d g of root diameter 2R r to gap thickness. If β < 1, the root can cross the gap without interacting with the two grains whereas for β > 1 the root will be unable to cross the gap without pushing the two grains apart.

Shape parameter
The main parameters of our system are the root stiffness parameters, growth rate, maximum length r of root elements, root diameter R r , outer and inner grain diameters (defining the gap thickness), limit outer grain diameters d min and d max (the distribution being uniform in grain volume fractions), mean grain diameter d , friction coefficients μ pp and μ pr , packing fraction and time step. For the roots, the normal stiffness is set to high values. They are less important than bending stiffness K b that closely controls the root shape. We set μ pp = 0.4 and μ pr = 0.4. These can change depending on the nature of the roots and grains. The time step and growth rate are small compared to the grain relaxation time under gravity. We performed extensive simulations for different values of the parameters. The packings are composed of 1500 grains and prepared by gravitational deposit inside a rectangular box. Here, we consider only the effects of β and packing fraction Φ.  > 2 Nm), the root grows downwards by dislodging all grains on its way and is only slightly deviated from the vertical due to elastic deformation. In this limit, the root behaves as a penetrometer pushed vertically downwards into the granular material. For low values of K b , horizontal and vertical excursions are observed in the trajectories of the root, and the root becomes more and more "noisy" at increasingly lower values of K b . The root cap, due to its lower bending stiffness, follows a rectilinear motion inside a pore. It can change its direction only if it meets the surface of a grain. Therefore, the length over which the root re-directions are significant is the mean grain diameter. We will therefore quantify θ Figure 3. Definition of re-direction angle θ of a root. these re-directions θ of the root by the total change of direction of root elements over a distance equal to d ; see  Figure 4 shows the distributions of re-direction angles P(θ) for different values of K b . All distributions are peaked on zero but their width increases as K b declines. This means that the shape of the root may simply be characterized by the standard deviation S of the set {θ i } of redirection angles: to which we will refer as shape parameter. Figure 5 shows the shape parameter S as a function of K b for two different values of β. We distinguish two regimes separated by a characteristic value K * b of bending stiffness. For K b < K * b , S decreases as a power law as K b increases, i.e. the root becomes less and less noisy. For K b > K * b , S continues to decrease even faster. This latter range corresponds actually to very small elastic deviations of the root as a whole whereas the range K b < K * b reflects the structure of the root. The power-law behavior allows to characterize the shape by an exponent:

Effect of root bending stiffness
The exponent is γ 0. 16  the fact that in this case the root diameter is thicker than the gaps between grains, leading thus to a faster decrease of S . for three different values of the packing fraction Φ and for β = 1.54. S declines as K b increases and its values are bigger for higher values of the packing fraction. But for these higher values, the resistance to penetration is also higher. The resistance to penetration is the force H or the penetration stress H/(4π d 2 ) required to push a rigid rod into the granular material. This force H multiplied by the average grain size d is homogeneous to a torque that can be compared to K b . In other words, we expect that when the bending stiffness is normalized by H d , then the data points collapse on a single curve. This is what we observe in Fig. 6(b). We also see that the transition point is K * b H d , and we have γ 0.13. This scaling with K * b works also for other soil parameters such as particle size and cohesion. The variations of the exponent γ are generally small and, within the natural variability of soil configurations, γ varies between 0.12 and 0.2.

Conclusions
In this paper, we introduced a model for root growth inside granular materials. This model is characterized by a constant root growth and accounts for the mechanical interactions between the root cap and grains. The root bending stiffness appears as a key parameter for the root shape. We introduced a 'shape parameter' that accounts for the degree of noisiness of the root shape, and showed that it varies as a power law with bending stiffness. The data points for different values of the packing fraction can be collapsed by normalizing the bending stiffness by a characteristic torque defined from the resistance to penetration and mean grain size. For bending stiffness above this characteristic value, the root grows by dislodging grains so that the shape is nearly rectilinear and the shape parameter is negligibly small. Below the characteristic value, the root can not dislodge the grains but follows the pore space. We also investigated the forces exerted by soil grains on the growing root and their relation with root shape. Those results will be published elsewhere.