Determination of macro-scale soil properties from pore scale structures: image-based modelling of poroelastic structures

We show how a combination of X-ray computed tomography (X-CT) and image-based modelling can be used to calculate the effect of moisture content and compaction on the macroscopic structural properties of soil. Our method is based on the equations derived in Daly & Roose (2018 Proc. R. Soc. A 474, 20170141. (doi:10.1098/rspa.2017.0141)), which we have extended so they can be directly applied to the segmented images obtained from X-CT. We assume that the soils are composed of air-filled pore space, solid mineral grains and a mixed phase composed of both clay particles and water. We considered three different initial soil treatments, composed of two different compaction levels and two different moisture contents. We found that the effective properties of the soils were unaffected by compaction over the range tested in this paper. However, changing the moisture content significantly altered the hydraulic and mechanical properties of the soils. A key strength of this method is that it enables the optimization or even design of soils composed from different constituents, with specific mechanical and hydraulic properties.


1.
We want to obtain a model which is valid on the macroscopic length scale Lx and has displacement of typical size Ly. The final equations we derive are valid up to O( 0 ), i.e., they will have error of size associated with them. Hence, we non-dimensionalise with respect to these length scales, where G is the shear modulus of the mixed phase, µ w is the viscosity of water, and k is the permeability of the mixed phase. We note that in the general case derived in [5], the interaction between the mixed phase and the adjacent phases was described by a generalised Navier slip condition, parameterised by a pair of slip lengthsλ s andλ a on the solid and air boundaries respectively. After homogenisation it was found that the parameterλ a does not play a role in determining the macroscopic properties of the soil [5]. Hence, we do not consider it in this paper. The parameterλ s is the slip length of the mixed phase acting on the solid particles. Physically, this is a function of the porosity, particle size, particle density, and the viscosity of the fluid. For a low permeability mixed phase, there will only be a thin water film between the solid component of the mixed phase and the larger solid grains. Hence, we would expectλ s 1. For a higher permeability mixed phase, there would be a larger water layer resulting in a largerλ s . This situation is somewhat analogous to the movement of water through a pipe with porous walls; in which case the appropriate boundary condition would be the Beavers and Joseph condition [6], where the slip length associated with a permeable medium isλ s = α √ k and α is a constant of order 1. Assuming a similar relationship holds here, the slip length will only become important ifλ s /Lx 1. Typically, we are interested in considering soils on the length scale Lx ∼ 10 m and mixed phases with permeability in the range k 10 −10 m 2 . Hence, we chose to setλ s to zero and not consider it further in this analysis, further details on the validity range for these equations can be found in [5]. This scaling results in the following set of non-dimensional equations for the x ∈ Bm, (S1.2a) where the effective stress is given by Hence, the total stress in the mixed phase is given by σ = σ m − Ip m . The solid phase is described through rigid body translations, Bs,j 1 dxê 3 , (S1.2e) and the air phase is described as a Stokes incompressible fluid, for r = {m, s, a}. Here e(u m ) = (∇u m ) + (∇u m ) T is the strain and I is the identity matrix, ν is the drained Poisson ratio, φ is the porosity of the mixed phase,τ i are vectors tangent to the solid particle surface, and σ s,j is the stress tensor on the j-th solid particle. We write the boundary conditions on the mixed-solid interface aŝ n ms · ∇p m + g wê 3 = 0, x ∈ ∂Bms, (S1.2j) n ms · σ s =n ms · σ m −n ms p m , x ∈ ∂Bms, (S1.2k) On the air-mixed phase boundary we write equations aŝ Finally on the air-solid interface we write the boundary conditions v a = v s , x ∈ ∂Bas, (S1.2p) n as · σ a −n as p a =n as · σ s , x ∈ ∂Bas. (S1.2q) Here the non-dimensional constants are where ρ w is the water density, ρ m is the density of the solid part of the mixed phase, g is the acceleration due to gravity, κ is the mean curvature of the air-water interface, and γ is the surface tension. In order to confirm our assumption, that the soils are approximatly capillary static, we can estimate the capillary number (on the micro-scale) as Ca = (kG)/(Lyγ) ≈ 1.38 × 10 3 , assuming k = 10 −10 m 2 , G = 10 8 Pa, Ly = 10 −3 m, and γ = 0.72 × 10 −3 Pa m. Hence, we note that capillary forces dominate over viscous forces and our assumption that the soil is capillary static is reasonable. We note this estimate of Ca is likely to be an overestimate as the length scale used is that of the representative volume rather than the actual pore size (which will be smaller).

(b) Cell problems
In order to derive an averaged model for the poro-elastic behaviour of soils the method of homogenisation was applied [5]. The approach involved deriving a set of cell problems which capture the effect of the underlying geometry on the mechanical and fluid properties of the soil.
Here our approach differs slightly from the one described in [5]. The key assumption used in [5] was that the geometry is periodic, i.e., it is composed of regularly repeating units. In reality soils will not have a perfectly periodic structure. Hence, we follow the approach of enforcing periodicity via reflection of the imaged geometry [7]. The aim, for each cell problem, is to obtain a reflected cell problem based on the reflections y k → −y k for k = {1, 2, 3}. The reflected cell problems are identical to the original cell problems under an appropriate translation of the dependent variable, (u → −u for example). By finding these translations we can infer that if the system is invariant under the translation y k → −y k and u → −u, then u must be odd about the plane of reflection. Hence, an appropriate boundary condition would be u = 0 on this plane. By repeating this approach for all directions we can infer a set of boundary conditions to impose on the edges of our domain which respect the symmetry of the underlying cell problems.
This symmetry reduction allows us to overcome the lack of periodicity in the imaged soils [7,8]. However, it comes at the price of forcing us to make assumptions regarding the structure and anisotropy of the imaged geometry. Clearly the soil image, Figure 1, is not composed of regularly repeating units. Hence, enforcing periodicity changes the soil structure. However, assuming that the soil structure is sufficiently homogeneous, enforcing periodicity by reflection will only create a significantly altered soil structure on the boundaries of the geometry. As the volume considered is increased the total volume fraction that is altered by the reflection is reduced. Hence, by increasing the volume of soil considered we would expect the effective properties measured to converge to those of the soil sample.
From a physical perspective, a periodic geometry composed of arbitrary repeating units could be isotropic, anisotropic or even chiral. However, if we assume that the representative unit has reflective symmetry in the three principal Cartesian directions then we cannot construct a chiral material. In addition the only anisotropic materials which respect these symmetries are those which are already aligned along the principal coordinate axes about which we make our reflections. In other words, we are only able to detect the diagonal components of any anisotropic material. These assumptions are necessary in order to obtain an appropriate representation of the soil on which the cell problems can be solved. In addition we have no reason to believe that these soil samples will have a chiral structure or that they will exhibit any off-axes anisotropy. It is possible that, as gravity acts on the soil in a single direction (aligned with the coordinate axes), there will be some variation between the effective parameters in the y 1 , y 2 and y 3 directions. The symmetry reduction we have considered will allow us to detect this difference.
The cell problems derived in [5] describe the Darcy flow of air through the porespace, the impeded Darcy flow of water in the mixed phase, the local stiffness, and the pressure induced  stress. In order to illustrate the symmetry reduction we start with the Darcy flow of air through the porespace, illustrated schematically in Figure S1.1. The cell problem derived in [5] comes from the first order expansion of Stokes' equations as , (S1.4c) (S1.4d) combined with the observation that p a 0 ∼ p a 0 (x), i.e., p a 0 is independent of the microscale variable, see [5] for full details. Hence, by writing where ζ k and ω a k are the local velocity and pressure coefficients driven by a large scale pressure gradient in the air phase. These can be determined by solving the cell problem where periodic boundary conditions were assumed. For the imaged geometries we assume that the imaged geometry comprises 1/8 th of the representative volume of soil. The total geometry can be reconstructed by reflecting the imaged data in y 1 = 0, y 2 = 0 and y 3 = 0. By imposing the translations ya → −ya, for a = {1, 2, 3}, we derive the following symmetry boundary conditions ω a k = 0, ∂ ∂y k (ê k · ζ k ) = 0,êp · ζ k = 0, p = k, y ∈ ∂y k , (S1.6d) ∂ ∂y j ω a k = 0, ∂ ∂yp (êp · ζ k ) = 0,ê j · ζ k = 0, p = j, j = k, y ∈ ∂y j , (S1.6e) where Ω conditions (S1.6d) and (S1.6e) enforce a normal flow condition in the direction of forcing and a slip condition in the non-forced directions. The symmetry reduction is illustrated in Figure  S1.1. In order to increase the clarity of the symmetry reduction we have illustrated this for an idealized sample in two dimensions. The boundary conditions (S1.6) are also applicable to a three dimensional soil image.
The symmetry boundary conditions for the cell problems, arising from the compressibility of the mixed phase, are more complex. We consider first the effect of strain on the soil displacement. From a macroscopic point of view, applying a strain to the soil sample will cause deformation of all components in the soil. However, as we have assumed that the soil minerals move as rigid bodies, they cannot be deformed. Hence, any macro-scale deformation must be coupled with a micro-scale variation in soil deformation which is equal and opposite to the macroscale deformation on the solid grains (see [5] for a more formal derivation). This deformation is captured by imposing that on the solid particles the response to deformation can be calculated by solving where κ u pq is the displacement coefficient driven by a large scale strain gradient and α u pq is the corresponding stress coefficient. The boundary conditions on the boundary ∂y k can be written in the compact form e k · [êpêq +êqêp +êpêp +êqêq − I] · κ u pq = 0, y ∈ ∂y k , (S1.8d) e k · [êpêp +êqêq] α u pq [êpêp +êqêq] ·ê k = 0, p = q, y ∈ ∂y k , (S1.8e) orê j · α u pq ·ê j = 0, j = k, p = q, y ∈ ∂y k . (S1.8f ) These boundary conditions are illustrated for the case p = q in Figure 3. The boundary conditions on the solid grains are more complex and we consider the cases p = q and p = q separately. For the case p = q we observe that, for any solid grain that intersects a boundary, y k = 0 or y k = 0.5 for k = {1, 2, 3}, the net displacement of that grain in the directionê k must be zero. In addition we find that the stress condition projected into the directionê k is automatically satisfied. Therefore, for a particle that intersects the boundary y k = 0 or y k = 0.5 we have (S1.8g) Alternatively, if the particle does not intersect the boundary y k = 0 or y k = 0.5 then we apply the conditions whereκpp is an arbitrary constant. Similarly for the case p = q, if the solid grains intersect with y k = 0 or y k = 0.5 we apply ms · α u pq ·ê k dy = 0, (S1.8k) and κ u pq + ypêq + yqêp 2 ·ê j = 0, j = k, y ∈ ∂Γ 1/8 ms , (S1.8l) whereκpq is an arbitrary constant. It is possible that a single particle can intersect multiple coordinate boundaries and could have zero motion in the y 1 , y 2 and y 3 directions. Alternatively it may have zero motion in one of these directions but not others. Hence, considerable care must be taken in defining the location of the grain boundaries in relation to the domain boundaries, Figure 3. The final cell problem describes the local deformation due to changes in pressure that move the soil away from the capillary static equilibrium, i.e., p a 0 − p c − p m 0 = 0. We write this as ∇y · α p = 0, y ∈ Ω 1/8 m , (S1.9a) where κ p is the displacement coefficient driven by a the absolute pressure. The symmetry boundary conditions arê e k · κ p = 0, y ∈ ∂y k , (S1.9d) e j · α pê j = 0, j = k, y ∈ ∂y k . (S1.9e) For a particle that intersects the boundary y k = 0 or y k = 0.5 we have (S1.9f ) Alternatively, if the particle does not intersect the boundary y k = 0 or y k = 0.5 then we apply the conditions κ p −κ p,j ·ê k = 0, y ∈ Γ 1/8 ms , (S1.9g) Γ 1/8 ms,jn ms · α p ·ê k dy = 0, (S1.9h) whereκ p is a constant vector. Equations (S1.6), (S1.7), (S1.8) and (S1.9) are the cell problems presented in the main body of the paper, i.e., equations (3.2), (3.3), (3.4) and (3.5).

(c) The averaged model
Using the results from the cell problems the final stage of the homogenisation procedure constitutes a volume average [5]. This approach was used to obtain a set of macroscopic equations for the solid displacement u m , the air pressure p a and the fluid pressure in the mixed phase p m , parametrised by the underlying geometry. In the original paper [5], this final stage of averaging resulted in 20 effective tensor quantities. Due to the symmetry of the cell problems in this study, several of these parameters become identically zero. Hence, we present the final dimensionless equations for the symmetric case. For the elastic displacement of the solid the averaged force balance is (S1.10) where C u ijkl is the stiffness tensor and C p ij is the pressure induced stress coefficient respectively, these are given by ∂κ u klq ∂y j + ∂κ u klj ∂yq + δ qk δ jl + ν 1 − 2ν δ qj δ kl + ∇y · κ u kl δ qj dy, (S1.11a) ∇y · κ p δ qj dy, (S1.11b) and The equation for the fluid pressure in the mixed phase is ∂ ∂t (S1.12) where K w ij is the effective water mobility, A u ij and A p are the increment of water content due to strain and pressure changes respectively. They are given by (S1.14) where K a ij is the effective air mobility, B u ij and B p are the increment of air content due to strain and pressure changes respectively. They are given by ∂Bas,jn am · γ u ij dy, (S1.15b) ∂Bas,jn am · γ p dy. (S1.15c) These equations provide a complete description of how the micro-scale geometry and mechanical properties influences the macro-scale mechanical properties of the poro-elastic composite. We note that if we consider the capillary static limit with an incompressible mixed phase then A u ij = A u δ ij and B u ij = B u δ ij , and equations (S1.10) and (S1.14) reduce to Darcy's law for the water and air phases respectively.
At this point we consider the effect of solving equations (S1.6) on a 1/8-th geometry. We notice that for the k-th cell problem, asê k · ζ k is even on all boundaries, we would need to multiply the result by a factor of 8 to account for the 7/8ths of the geometry not considered. However, as the effective parameters are all normalized by the cell volume, this factor of 8 will cancel out. In addition, as each of the remaining velocity components is odd about two boundaries, we find Ωaê j · ζ k dy = 0, j = k. (S1.16) Hence, K a ij is zero for p = q. By a similar analysis as for (S1.6) we find that K w ij , C p ij , A u ij and B u ij are zero if i = j, and C u ijkl is zero if i = k and j = l. In summary, the effective equations describe a three phase poroelastic material and are parameterised by 9 effective parameters. These are; g ef f the effective gravitational force; C u ijkl , the macro-scale or effective stiffness tensor; C p ij the stress coefficient for pressure variations about the capillary static point; K w ij and K a ij , the water and air relative mobilities; and A u ij , A p , B u ij and B p , the strain and pressure induced variations in water and air content.