Channelization in Porous Media driven by Erosion and Deposition

Channelization in Porous Media driven by Erosion and Deposition

R. Jäger ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials, Wolfgang-Pauli-Strasse 27, HIT, CH-8093 Zürich (Switzerland)    M. Mendoza ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials, Wolfgang-Pauli-Strasse 27, HIT, CH-8093 Zürich (Switzerland)    H. J. Herrmann ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials, Wolfgang-Pauli-Strasse 27, HIT, CH-8093 Zürich (Switzerland)

We develop and validate a new model to study simultaneous erosion and deposition in three-dimensional porous media. We study the changes of the porous structure induced by the deposition and erosion of matter on the solid surface and find that when both processes are active, channelization in the porous structure always occurs. The channels can be stable or only temporary depending mainly on the driving mechanism. Whereas a fluid driven by a constant pressure drop in general does not form steady channels, imposing a constant flux always produces stable channels within the porous structure. Furthermore we investigate how changes of the local deposition and erosion properties affect the final state of the porous structure, finding that the larger the range of wall shear stress for which there is neither erosion nor deposition, the more steady channels are formed in the structure.

I Introduction

Fluid flow through a porous medium can erode and/or deposit material and thereby change the shape of the solid boundaries, which in turn alter the path of the flow. These processes are responsible for shaping a variety of landscapes, some examples are meandering rivers middleton2005encyclopedia (), coastal erosion masselink2014introduction () and seepage cedergren1977seepage (). In industrial applications these processes also play a crucial role and are sometimes desired and sometimes disruptive, e.g. mechanical filters take advantage of deposition whereas in oil wells erosion produces sand that can damage the pumps, requiring expensive counter-measures. While these processes are easily observed in surface flows, it is much more difficult to do so when they occur within. Albeit there is an abundance of applications, the evolution of structures exposed to erosion and deposition is in general not easily predictable and usually requires the use of numerical techniques.

Erosion is the mechanical wearing of solid material by the shear force exerted by the fluid. The erosion rate is assumed to be proportional to the magnitude of the wall shear force FLM:57479 (); Bonelli2006555 (). For the erosion and deposition two thresholds for the shear force are introduced. Above a critical threshold the shear stress is high enough to erode matter from the solid surface whereas below a lower threshold the adhesive force of the suspended particles is dominant and deposition occurs. In the range between the thresholds there is locally no net change of mass at the surface. While this is an heuristic approach, it is well-founded, as it is consistent with the Hjulström curve Grabowski2011101 () which describes a gap between the erosion and the sedimentation flow velocity for all grain sizes. To accurately describe the erosion and deposition induced by the wall shear force and adhesive forces within a porous media, the Navier-Stokes equations must be solved and the trajectories of suspended particles have to be calculated. A commonly used approach is to describe the suspended particles as a solute SALLES19932839 (); Bouddour1996 (); yamamoto2009fluid (), whose evolution is given by the convection-diffusion equation. Even though these equations can be written down easily, solving them with realistic boundaries defined by the porous structure cannot be done analytically. Hence we develop a numerical model based on the lattice Boltzmann method (LBM), which is a very powerful tool to resolve fluid behavior in flows through porous media talon2012assessment (), that allows to describe the erosion induced by wall shear stress and also the deposition of fluid entrained matter.

Using our model we find that when both erosion and deposition are active, channels within the porous structure form. Mahadevan et al. mahadevan2012flow () and Kudrolli et al. PhysRevLett.117.028001 () already studied channelization in porous media using experiments and two-dimensional effective models for deposition and erosion. In both cases the fluid is modeled using Darcy’s law, which corresponds to a scale much larger than the pore size, and therefore, prohibiting the use of local shear stress forces at the solid boundaries to characterize the erosion and deposition processes. In particular, in the model of Mahadevan et al. the erosion depends on the pressure gradient of the fluid and Kudrolli et al. consider erosion and deposition depending on fluid velocity. In our model however, the Navier-Stokes equations are solved to resolve the fluid flow at pore scale and erosion can be considered to depend on the local shear stress exerted by the fluid onto the solid surface. Our goal is to study the microscopic physics of erosion and deposition and the resulting macroscopic changes of the three-dimensional porous structure. Since we want to study the interplay between erosion and deposition we only consider consolidated porous media; and hydraulic pressures that are high enough to erode material, but do not lead to a decompaction or fluidization of the porous medium PhysRevE.78.051302 ().

This paper is organised as follows: in Sec. II we present the set of differential equations needed to understand the erosion and deposition processes in porous media; in Sec. III, we introduce our model and perform some validations; and Sec. IV shows the results of our study on erosion and deposition in porous media including the channelization process. Finally, in Sec. V we summarize our results and conclusions.

Ii Theory

ii.1 Fluid and Particle Dynamics

There are several mechanisms that are responsible for the mechanical erosion and deposition and have to be considered to get the full description of evolution of the porous structure. First, the fluid flowing through a porous medium is described by the incompressible Navier-Stokes equations, which consist of the continuity equation,


and momentum conservation equation


where is the velocity field of the fluid, is the fluid density, the kinematic viscosity, the hydrostatic pressure, and an acceleration due to external forces. In our approach, the fluid contains suspended particles that are entrained by the fluid and can be deposited onto the solid surface. The suspended particles that we consider are small, get transported by the fluid, diffuse within it and are neutrally buoyant. Therefore they are described by the convection-diffusion equation, which yields the mass concentration of solid particles in the fluid. This equation reads:

where is the diffusion coefficient and is the particle concentration. Note that the concentration inside the porous medium will also change due to deposition and erosion as we will describe in the following section. Furthermore the concentration is assumed to be low, such that the change in rheology and the momentum exchange between fluid and suspended particles is negligible. By solving the Navier-Stokes and the convection-diffusion equations, we describe the evolution of the fluid and the suspended particles through the porous media.

ii.2 Erosion and Deposition

The wall shear stress (WSS) exerted by the fluid on a solid surface is defined by:


where is the dynamic viscosity, the fluid velocity parallel to the surface, and the distance from the surface. When the WSS overcomes the cohesive forces of the solid, solid matter is eroded and entrained by the fluid. On the other hand, when the adhesive forces, e.g. Van der Waals forces, are stronger than the erosive and other repulsive forces, particles within the fluid are deposited on the solid surface. Through these processes the solid surface and thus the boundary is altered and the flow changes. We consider slow erosion and deposition, where the change of the surface is much slower than the fluid flow velocity. Hence, the boundary can be considered quasi-static (or evolving slowly), i.e. non-moving, and only the magnitude of the deviatoric shear stress and not the hydraulic pressure has to be taken into account for the evolution of the surface. The shear force is calculated from the deviatoric shear tensor given by:


multiplied with the surface normal vector , yielding the shear force . The WSS is now the tangential part of the shear force which is obtained with the following expression:


the details for calculating the WSS are shown in section III.5. In the case of erosion, we assume that this WSS works as a wearing force producing sand or silt out of solid matter which is entrained with the flow. The erosion is assumed to be linearly proportional to the wall shear stress, this erosion law is widely used and experimentally verified Bonelli2006555 (); FLM:57479 (). The erosion is proportional to the stress gap, i.e. the WSS excess above some critical threshold and depends on the density and toughness of the solid matter. This dependency is incorporated in an erosion coefficient and a critical shear stress below which no erosion occurs. The porous medium is considered consolidated and incompressible such that the normal component of the stress tensor can be neglected and other erosive processes, such as fluidization PhysRevE.78.051302 () are not considered. For the deposition law we also assume a threshold of the WSS below which deposition occurs. The deposition of matter is also considered linearly proportional to the WSS below the threshold, hence at the critical deposition shear stress no deposition occurs. The deposition is naturally dependent on the concentration of solid matter in the fluid. This gives us an erosion and deposition law, which written in differential form, looks as follows:


where is the change in mass per unit area and is the fully saturated concentration. In the following the relative concentration will be simply referred to as . While the linear relation of erosion to wall shear stress is an empirical law, it is well established and reasonable for a variety of materials such as cohesive soils, clay and materials of similar properties FLM:57479 (); Bonelli2006555 (); bonelli2013erosion (). Furthermore, as we are considering wall shear stresses which are in the vicinity of the critical shear stress, higher order terms would only contribute minimally to erosion. The deposition of particles is usually assumed to be proportional to the concentration SALLES19932839 (); Bouddour1996 (); yamamoto2009fluid (), when the solubility is zero, however recent experimental findings Alem2015 () show that deposition decreases with increasing flow speed, thus we chose the simplest form that takes this into account, but converges to for low flow speeds. The erosion coefficient and the deposition coefficient depend on the properties of the solid matter, e.g. the cohesive forces Grabowski2011101 ().

Iii Numerical Implementation

Our numerical model is based on the lattice Boltzmann method (LBM) wolf2000lattice () with an improved two-relaxation-time (TRT) collision operator PhysRevE.83.056710 (); talon2012assessment (). In order to compute accurately solid boundaries and shear stress forces we also implement and analyze different interpolation and extrapolation schemes Filippova1998219 (); bouzidi (); Mei2000680 (). In particular, this avoids the use of the stair-case approximation for curved boundaries. In this work, parameters are written in numerical units, the conversion from numerical to physical units is illustrated in the appendix A.6.

iii.1 Lattice Boltzmann Method

To resolve the fluid flow at pore scale we use the D3Q19 lattice Boltzmann method (LBM) wolf2000lattice (), which means the model is based on a three-dimensional lattice with nineteen discrete velocity vectors . In contrast to other computational fluid methods that solve the Navier-Stokes equations directly, the LBM solves the underlying discrete Boltzmann equation FLM:409537 ():


where are distribution functions associated with the discrete velocity vectors , and a generic collision operator . Considering the Bhatnagar-Gross-Krook (BGK) approximation for the collision operator the lattice Boltzmann equation gets the following form FLM:409537 ():


where is a relaxation time, that is related with the kinematic viscosity of the fluid by . The equilibrium distribution function is given by wolf2000lattice ():


where and are the local density and velocity of the fluid, are weights (associated with the discrete velocity vectors) and the speed of sound is given by . The values of the weights and discrete velocities can be found in Ref. wolf2000lattice ().

The fluid density is computed from the zeroth moment of the distribution function, i.e. , and the fluid velocity from the first moment, . Using the Chapman-Enskog expansion it can be shown that Eq. (8) with Eq. (9) converges to the Navier-Stokes equations Eqs. (1) and (2) in the continuum limit wolf2000lattice ().

iii.2 Improved Collision Operator

The algorithm for the generic lattice Boltzmann equation (7) can be split into collision part

and streaming part

For the BGK collision operator (see Eq. (8)) the discretization error is viscosity dependent Pan2006898 (), however this problem can be mitigated in two different ways; one is to simply increase the resolution which will decrease the error; and the other is to use a multi-relaxation-time (MRT) collision operator Pan2006898 (). The two-relaxation-time (TRT) collision operator is by far the most efficient and its accuracy for porous media is as good as any other MRT method PhysRevE.83.056710 (); talon2012assessment (). The TRT method splits the distribution functions into symmetric and asymmetric parts , where is the index of the velocity vector pointing in the opposite direction of , i.e. . The projected equilibrium distribution functions read .

For the TRT method the full collision operator reads:


where the second relaxation time is a free parameter relaxing the asymmetric part of the distribution functions. To minimize the dependency of the viscosity the parameter is chosen to be constant () and therefore the second relaxation time is determined as porescalelbm (). The collision for the distribution functions and can be calculated simultaneously via Eq. (10), which allows the implementation of the TRT to be almost as fast as the one of the single-relaxation-time (BGK) approach.

iii.3 Solving the Convection-Diffusion Equation

To model the convection-diffusion with the lattice Boltzmann method a second set of distribution functions can be used to describe the mass concentration of the solute Yamamoto2584 () within the solvent:


Here, the equilibrium distribution is similar to the one for the fluid (see Eq. (9)) with the concentration instead of the density:

The diffusion coefficient is related to the relaxation time , . For the TRT method the evolution of the solution is calculated analogously to Eq. 10 with the symmetric and asymmetric splitting and the second relaxation parameter . Analogously to the calculation of the fluid density, the zeroth moment of the concentration distribution function is the concentration of the solute, , and via a Chapman-Enskog expansion, one can show that the lattice Boltzmann equation (11) converges to the convection-diffusion equation. Note that the positivity of the kinematic viscosity and diffusivity requires , and the model becomes unstable for . For the simulations shown in this paper we only used the relaxation times 1 and 0.6, however we found our model stable for values down to 0.503.

iii.4 Solid mass field to identify fluid and solid domains

To distinguish fluid and solid domains we introduce a new variable, a scalar mass field that is defined on the lattice, for fluid nodes the mass is zero, for solid nodes the mass is one and for interface cells the mass is between zero and one:

If the mass is one, the cell is still an interface node if there is a neighboring fluid node. Interface nodes have both solid and fluid matter and therefore also non-zero distribution functions . When an interface node becomes a solid node the distributions are set to zero . When a solid node becomes interface the distribution functions are linearly extrapolated so as not to introduce an artificial pressure or velocity gradient as was done by Yin et al. yin ():

where the sums run over all neighbors that are fluid nodes.

When the solute is locally deposited onto the solid surface or solid matter is eroded and added to the fluid the concentration field is changed locally by the mass difference . For this mass in-/ejection not to introduce any artificial momentum into the fluid, the mass change is equally distributed according to the weights: . This guarantees that the local momentum of the exchanged mass is zero as we have .

iii.5 Calculation of the Wall Shear Stress

The tangential wall shear stress acting on a surface is calculated from the shear stress tensor given by the fluid. There are two ways to compute the shear stress tensor. For the LBM the shear stress tensor can be calculated from the non-equilibrium part of the distribution functions:


where and denotes the component of the discrete velocity vector . Additionally, Thampi et al. Thampi20131 () propose a very accurate way to calculate derivatives of physical properties in LBM. From this we can calculate the shear tensor as:


where the prefactor depends only on the weights and discrete velocity vectors and needs to be computed just once:

We found that Eq. (13) is more accurate for calculating the shear tensor, but it requires that all neighboring nodes are fluid nodes. Therefore we use Eq. (13) wherever possible and Eq. (12) elsewhere. The wall shear force can be calculated as the product of the normal vector on the boundary surface times the deviatoric shear stress . Since the boundary in general does not lie on a fluid node we linearly extrapolate the shear stress onto the boundary as follows:

where is the shear force at the interface cell, is the relative distance from the node to the wall and the sum runs over all neighbors () that are fluid nodes. For the erosion and deposition law we only take into account the absolute value of the tangential component of the wall shear force (see Eq. (5)). The unitary normal vector to the surface can be calculated by the color gradient of the mass field (see appendix A.1). The normalized surface vector must be multiplied by the surface area to find the total WSS acting within an interface cell onto the surface, the total surface vector is then . For each interface node we need to estimate the surface area that is exposed towards the fluid, this area is estimated as the area of the plane defined by the surface vector , cut by the cells boundary, hence the surface is piecewise flat. Therefore the area is unity if the surface is parallel to the lattice planes and is maximum  , e.g. when . The total mass exchange caused by the erosion and deposition law within an interface cell is then .

iii.6 Boundary Interpolation Scheme

The commonly used no-slip boundary condition in LBM is imposed by the bounce-back method, which puts the boundary either on lattice nodes (on-site bounce-back) or in the middle of two lattice nodes (half-way bounce-back). In order to allow for boundaries that are more generic and do not lie on points restricted by the lattice one has to modify the bounce-back method. There are several interpolation schemes to account for generic boundaries. Two of the most widely used schemes are from Bouzidi et al. bouzidi () and Filippova et al. Filippova1998219 (). The interpolation scheme from Filippova was improved by Mei et al. Mei2000680 () to allow for lower viscosity. We implemented both methods and compared them simulating the Poiseuille flow through a pipe of radius . The analytical solution for the velocity profile of the Poiseuille flow is as follows:

where is the distance to the center of the pipe, is the radius of the pipe, with a fluid driven by a pressure drop over length . For a force driven fluid the pressure drop can be substituted by the driving force, . For simulating the Poiseuille flow, we set the mass field to zero inside and outside to one . For interface cells the mass depends on the distance to the center . Here we switch off the erosion and deposition of solid material. The fluid is driven by a constant force in the direction of the pipe and we employ periodic boundary conditions at the inlet and outlet, hence only a one layer cross-section of the pipe has to be stored. The two-norm of the error of the stationary flow field , as defined in appendix A.2, was compared for the two interpolation schemes and the half-way bounce back (see figure 1). It was found that both interpolation schemes show very similar accuracy (second order) and are indeed superior to the classical half-way bounce back scheme (first order).

Figure 1: Poiseuille flow through a pipe with a given driving force and (in lattice units). The two-norm of the error (see appendix A.2) of the flow field for the Bouzidi and the Mei interpolation schemes and the half-way bounce back, combined with the BGK or the TRT (denoted by TRT) collision operator, are compared. The interpolation schemes show similar accuracy and scale better than the simple bounce back.

These interpolation schemes can also be used together with the TRT collision operator and we observe that there is no considerable difference in accuracy when compared with BGK (see figure 1). For stability reasons we settled on using the interpolation by Mei for our model.

iii.7 Validation

To verify that our model for erosion and deposition is working properly a setup is used where the flow field, the WSS, and hence also the surface evolution can be calculated analytically. For this purpose we consider pipe flow where the fluid is described by Poiseuille flow. The setup used for the simulations is the same as already described in section III.6, except that now the interface cells and hence also the surface evolve according to the erosion and deposition laws as described in section II.2.

iii.7.1 Validating Erosion: Erosion of a Pipe

To measure the erodibility of cohesive soils the hole erosion is commonly used (see e.g. Ref. wahl2010comparison ()). In this test water is flowed through a soil pipe under constant pressure drop and the soil is eroded and thus the pipe grows in width. Bonelli et al. Bonelli2006555 () developed an analytical theory to predict the evolution of the pipe given specific properties of the soil and fluid. They found a scaling law for the evolution of the radius that depends on the pressure driving the fluid, the critical shear stress and the erosion coefficient :


the erosion time is with the density of the erodible soil. We use the same setup except that instead of a pressure difference driving the fluid we introduce a constant driving force in the direction of the pipe. The force is implemented for the BGK as derived by Guo et al. PhysRevE.65.046308 () and for the TRT as described by Seta et al. seta2014implicit (). The pipe erosion law derived by Bonelli shows good compliance with hole erosion experiments and is analytically exact for slow erosion. The wall shear stress depends only on the driving force and pipe radius: . For low erosion coefficients our simulations show excellent agreement with the erosion scaling law as can be seen in figure 2. The flow is quasi-steady and for all times close to the analytical Poiseuille profile. However as the erosion rate increases with larger pipe radius the condition of slow erosion breaks down at some point. We can observe in figure 3 that for low erosion coefficients () the simulated erosion of the pipe is close to the theoretical scaling but deviates for higher erosion coefficients.

Figure 2: Pipe erosion for different erosion coefficients. Initial pipe radius is (in lattice units), the erosion was started after a steady flux was reached (). The time is rescaled by the erosion coefficient . A relaxation time was used, the erosion threshold is zero, , and the driving force is (in numerical units). The theoretical curve (blue line) is given by equation (14).
Figure 3: Convergence of the erosion time for small erosion coefficients. The plot shows the difference between the measured erosion time and the theoretical one . The setup is the same as in the simulations shown in figure 2, only the erosion coefficients are different. We see that the boundary interpolation scheme improves the accuracy of the erosion time . The theoretical erosion time is , denoted by the dashed line. The erosion time is measured by fitting the measured pipe radii (see figure 2) to the analytic equation (14).

iii.7.2 Validating Deposition: Clogging of a Pipe

To verify that our model also works properly for deposition, we start with an initial steady Poiseuille flow through a pipe, set a constant particle concentration and set the deposition threshold above the wall shear stress. For Poiseuille flow we know that the WSS is and the deposition law is given by . Thus we can derive a differential equation for the radius:

which is solved by the exponential function:


where we assume a constant concentration over time.

Figure 4: A pipe with initial radius 60 (lattice units) is clogged by the matter suspended within the fluid, the deposition threshold is set to , the relaxation time is , the driving force (in numerical units) and the concentration to . The smaller the radius becomes, the lower is the WSS and the deposition will attain a linear behavior . The theoretical expression is given by equation 15.

In figure 4 we see that the simulations are in good agreement with the analytic prediction, however this prediction is only valid for slow deposition where the flow field is quasi-steady and therefore close to the Poisseuille profile at all times.

Iv Channelization in Porous Media

Our model is now applied to simulate erosion and deposition in a soft porous medium, by soft we mean that all solid matter within the simulation is erodible. For these simulations we use two basic setups, the first is to impose a constant pressure drop between the inlet and outlet, and the second enforces a constant flux at the inlet. For both setups the fluid driving forces are assumed to be much stronger than gravitational forces and therefore the latter is neglected. The flux is imposed using inlet boundary conditions as described by Kuzmin et al. Kuzmin2219 (). Due to computational power restrictions it is challenging to simulate a porous medium of the same size as occurring in nature or in experiments. To minimize the finite size effects we employ periodic boundary conditions at the walls. The porous medium consists of randomly placed spheres (of radius 7.1 grid points), if a sphere is cut by a wall, it will continue on the opposite side of the wall. Effectively this setup is equivalent to a porous medium with recurring pattern and infinite width. Spheres can also overlap and are placed until a porosity of 0.5 is reached. The inlet and outlet regions ( and cells, respectively) are kept free of solid matter. The dimensions of the simulation box are (, where is the flow direction). The particle concentration is fixed to at the inlet, however it changes within the porous medium due to deposition and erosion. Since we are interested in laminar flow; and higher relaxation time parameters achieve faster convergence to stationary states; the relaxation times are set to for all following simulations. For both setups we first use a sharp transition between deposition and erosion, i.e. .

iv.0.1 Constant Pressure Drop

For a constant pressure drop we find that there is a critical value for the erosion/deposition threshold (see figure 5). When the threshold is set below this critical value, the porous medium will completely erode and the flux will diverge. On the other hand, if it is set above the critical value the porous medium clogs completely and the flux goes to zero. We also see that the total time, until a final state is reached, increases when we go from eroded to clogged final states (see figure 5). The reason behind this is simple, clogging will reduce the flow and slow down the process, whereas erosion increases the flux and thus is a diverging process. The total mass of a clogged medium increases for low thresholds, but only slightly (see figure 5), this is due to a slower clogging process, which allows to fill up more pores with solid matter. The time it takes for the system to reach its final state grows the closer the threshold is to the critical one. However, before the porous medium is completely eroded or clogged, we see the formation of channels within the porous structure (see figure 6). The size of the largest channel formed will decide whether the medium will end up eroded or clogged, as already seen in the verification examples (section III.7), a pipe and analogously also a channel with a constant pressure drop will always continue to clog or erode (see appendix A.3). The critical threshold is however not generic, it depends not only on external parameters but also on the initial porosity and on the specific inner structure of the porous medium. For example, it increases for higher pressure drop or larger pores, because the occurring shear also increases. As will be seen in the next application, there is in theory a critical pipe radius for which there is neither erosion nor deposition, however such a configuration is very unstable and even small fluctuations in the flux or shape will lead to an outcome of complete erosion or clogging (see appendix A.3). Therefore we can conclude that when considering a sharp transition between erosion and deposition, there are no configurations for constant pressure drops that lead to steady non-trivial structures.

Figure 5: A porous medium with initial porosity is evolved by a fluid that is driven by a relative pressure drop of . The inlet solute concentration is , the coefficients , and . We measure the time until there is no more deposition or erosion and the total mass of erodible matter contained in the final structure. The critical threshold (dashed vertical line) is at . Below the critical threshold, the medium ends up completely eroded, there is no more solid mass in the simulation; above, it ends up clogged. The size of the simulation box is , where the inlet and outlet (each a length of cells) are constrained from deposition to avoid changing the inlet boundary conditions of the fluid. While specific numbers for time and total mass change with the system size, the two regimes of clogging and complete erosion remain.
Figure 6: Snapshots of the evolution of the porous medium for constant pressure drop and the shear threshold just below the critical value (see figure 5). The solid domain is shown in solid gray, the fluid in transparent color; low velocity fluid is shown in light blue, medium velocity in green, and high velocity in red. In this case the medium ends up completely eroded. Snapshot I () shows the initial stage, II () shows when discrete channels are formed, and at stage III () we see that smaller channels vanish while the biggest one grows.

iv.0.2 Constant Flux

Contrary to the constant pressure drop, a constant flux leads to a final state of the porous medium with stable channels. For the case of a sharp transition from deposition to erosion () the medium ends up having only one channel. As can be expected due to symmetry reasons it closely resembles a pipe, as the pipe is the only channel form that has a uniform constant WSS for steady flow. The flow then also closely resembles Poiseuille flow, the difference to a perfect pipe stems from the fact that the inlet and outlet are not a pipe but rectangular by construction and so disturb the form slightly, see figure 7 as an example.

Figure 7: Steady one channel state after erosion and deposition for a threshold . The flow direction is annotated by the (red) arrow. At the inlet the shape is more squared and wider, further inside the channel narrows and its shape more closely resembles a circular pipe.

We have found that the pipe is a bit wider at the inlet and then narrows slightly towards the outlet.

The pipe radius was measured in the middle of the simulation box (), which is of the same size as for the constant pressure drop simulations (). For Poiseuille flow we know that the relation between radius and flux is given by and the WSS by , from which we can find the critical radius


The critical WSS is the same as the erosion and deposition threshold . The formula for the critical radius shows that the radius depends on the imposed flux , the dynamic viscosity ( and the shear threshold and we can show that a pipe with this radius is indeed stable (see appendix A.4). We compare the measured radius of the channel with the expected radius and find that they are in very good agreement (see figure 8). Theoretically, having two or more pipes with the same radius could also be a steady configuration, it would however only be meta-stable, small fluctuations would result in clogging one pipe and widening the other (see appendix A.5).

Figure 8: At the inlet a constant flux with constant velocity is imposed, the inlet solute concentration is set to 0.1 and the relaxation parameters to one, . The erosion and deposition coefficients are the same as in the previous setup (fig. 5). The radius of the channel is measured in the middle of the simulation box (), the expected radii are calculated from equation 16. The system has to be larger than the diameter of the pipe, otherwise additional boundary effects will occur.

iv.0.3 Constant Flux and Gap

After finding the steady state solutions for a sharp transition we now study a system with a gap between deposition and erosion. This means that there is a range of WSS for which matter is neither removed nor attached. For small gaps we found that the final configuration still features one channel (see figure 9), which however deviates more and more from a perfect pipe, the wider the gap, its form depending on the original porous structure. For even larger gaps we found that there are steady states with more than one channel, the number of channels depends not only on the gap but also on the flux and the size of the simulation box. The smaller the pipe formed by erosion/deposition without a gap, the more channels can be expected for a final state with gap. We can also see in figure 9 that the permeability monotonously decreases for increasing gap until there is a minimum, which is also where many distinct channels are formed. The permeability increases again when going to even larger gaps which leads to final structures that do not show distinguishable channels but depend strongly on the original porous structure. For very large gaps there is neither erosion nor deposition and the porous structure is static, in the setup of figure 9 this is the case for gaps larger than . Note that corresponds to the solution shown in figure 8, while changing or the system size changes the permeability and the number of channels; it also limits the range of ; the qualitative behavior is the same as long as is such that both erosion and deposition are active. The final state depends strongly on the local deposition and erosion, the stronger the flow changes the porous structure, the less the final state depends on the original structure. Although the gap in this case cannot be adjusted, the final structure can be influenced by adjusting the flow velocity and hence the shear within the porous medium. Mahadevan et al. mahadevan2012flow () concluded that channelization can arise due to the preferential erosion induced by the flow. As channelization cannot progress further when there is only one channel left we showed that this effect is stronger the smaller the difference between deposition and erosion thresholds.

Figure 9: Starting from the simulation shown in figure 8 with the threshold we introduce and increase a gap between the deposition and the erosion threshold . In addition to measuring the final permeability, the number of channels was counted at the inlet and the outlet. For each gap, we made simulations for four random porous media made of solid spheres (with radius 7.1 grid points) with different seeds for the random number generator (for the placement of the spheres). The first range (I) includes final states with only one channel, in the second range (II) there are multiple channels and bifurcations. In the third range (III) channels are not distinguishable anymore, but the medium is a porous structure that depends strongly on the original structure.
Figure 10: Snapshots of the final state of the solid structure after erosion and deposition for different gaps between erosion and deposition thresholds. These snapshots show examples for the three different regions painted in figure 9. The first snapshot (I) corresponds to a gap of , the second (II) to and the third (III) to .

V Conclusion

We have found that erosion and deposition within porous media always lead to channelization. However stable channels for any material can only be expected if the fluid features a constant flux. For a fluid driven by constant pressure drop, the final configuration is either completely clogged or completely eroded. On the other hand a constant flux will always leave at least one stable channel. This statement is also supported by the experimental findings of Alem et al. Alem2015 (), they tested deposition in a filter, and found only completely clogged final states for constant pressure drop, but steady state configurations for constant flux. We also showed that the final state of a porous medium altered by erosion and deposition crucially depends on the balance between the two effects and a circular pipe is formed when the transition is sharp. While specific values such as permeability or number of channels change with the system size; and the threshold for erosion and deposition; the qualitative behavior remains the same. Our numerical findings may help to better understand and mitigate internal erosion in embankment dams erosiondams () and also in naturally occurring flows through porous media. In our model gravity is not considered, if gravitational forces are not negligible, buoyancy and sedimentation of suspended particles play a role, hydrostatic pressure increases towards the lower part of the fluid and also the rate of deposition changes PhysRevLett.97.138001 (). To study these effects in the future, it is possible to incorporate gravitational forces into our model by adding an external force to the fluid dynamics equations and by making the deposition dependent on particle density.

While our model allows to simulate erosion and deposition for a variety of fluids and materials, it cannot capture other interesting phenomena in flows through porous media such as decompaction and fluidization PhysRevE.78.051302 (). However, extensions of the model to study these effects can be a subject of future research.

We acknowledge financial support from the European Research Council (ERC) Advanced Grant 319968-FlowCCS.

Appendix A Appendix

a.1 Color Gradient

The solid and fluid domains are defined by the mass scalar field (see section III.4). To get the normal vector of the surface we define the color gradient in the following way:

This equation defines the orientation of the surface and allows to compute the normal and tangential part of fluid properties, e.g. the shear force.

a.2 Error measurement

To measure the difference between an analytic velocity field and a numerical computed field a measure for the error must be defined. It is common practice to use the two-norm of the difference:

where the sum runs over all nodes in the fluid domain. This measure allows to evaluate the accuracy of a fluid dynamics computation (see figure 1).

a.3 Constant pressure drop channel instability

For Poiseuille flow through a pipe the wall shear stress is given by:

The critical radius for which there is neither erosion nor deposition would be . However, for a sharp transition between erosion and deposition () a small fluctuation in radius or pressure difference will lead to complete erosion or complete clogging:


These equations illustrate why a channel is unstable under constant pressure drop (see section IV.0.1).

a.4 Constant flux channel stability

For a Poiseuille flow through a pipe with constant flux, the flux is determined by , and the shear by . Thus, we can derive an equation for the shear depending on the flux:

From this it is easy to see that the pipe with critical radius is stable, for a sharp transition from erosion to deposition () :


These equations show that for the same flux a larger pipe will result in lower flow velocity and thus in a lower shear, which results in deposition; and vice versa for pipes with radii smaller than the critical one. This is the reason why a constant flux will always form one stable channel (see figure 8).

a.5 Stability analysis for a two channel state

A final state with two (or more) pipes is theoretically possible for sharp transitions () and would yield a critical radius of

where we assume the pipes lie far apart and hence assume a Poiseuille profile for both flows. This is stable for perturbations that are symmetrically occurring for both pipes. However a perturbation of just one pipe will widen one and clog the other pipe. Assume there is a perturbation of flux and the flux in one channel increases by , since the total flux is constant, the flux in the other channel decreases by the same amount. The WSS in the first channel now increases:

thus the first channel will start to erode and the second will clog until it vanishes. Even though this is only a sketch it shows that configurations with two (or more) channels are unstable.

a.6 Conversion from numerical to physical units

For comparing numerical with experimental findings, numerical units have to be converted to physical units, a very instructive recipe how to do this is written by J. Latt latt2008choice (). Here we provide an example for converting the units of the simulation shown in figure 5, for this purpose we match three parameters, namely the fluid density , the viscosity and the length scale . The parameters are matched using the assumption that the diameter of the spheres building the porous medium is 1.4 mm; and the fluid is a heavy crude oil Hart2014 () with density kg/m and viscosity m/s. Hence all other parameters can be calculated using the unit conversion scheme shown in latt2008choice (), e.g. for the time scale:


where is the viscosity in numerical units, and in physical units.

parameter numerical units physical units
1 kg/m
1/6 10 m/s
1 0.1 mm
1 s
1 0.017 s/m
2.8 9689 Pa
Table 1: This table shows an example conversion from numerical to physical units. The first three parameters () were set to match the fluid properties of a heavy crude oil and a porous structure of pore size 1.4mm, the other parameters () can be calculated from the first three as shown in Eq. (21).

Note that and are difficult to match with experimental data, since they depend highly on the material and the formation of the porous medium. For instance, while the erosion coefficient is close to values reported by Bonelli et al. Bonelli2006555 (), the erosion threshold is almost two orders of magnitude higher.


  • (1) V. Middleton, M. Church, M. Coniglio, L. Hardie, and F. Longstaffe, Encyclopedia of Sediments and Sedimentary Rocks. Encyclopedia of Sediments and Sedimentary Rocks, Springer Netherlands, 2005.
  • (2) G. Masselink, M. Hughes, and J. Knight, Introduction to coastal processes and geomorphology. Routledge, 2014.
  • (3) H. Cedergren, Seepage, drainage, and flow nets. Wiley Classics in Ecology and Environmental Science Series, Wiley, 1977.
  • (4) G. Parker and N. Izumi, “Purely erosional cyclic and solitary steps created by flow over a cohesive bed,” Journal of Fluid Mechanics, vol. 419, pp. 203–238, 9 2000.
  • (5) S. Bonelli, O. Brivois, R. Borghi, and N. Benahmed, “On the modelling of piping erosion,” Comptes Rendus Mécanique, vol. 334, no. 8–9, pp. 555 – 559, 2006. Observation, analysis and modelling in complex fluid mediaSpecial issue for the 60th birthday of Professor Roland Borghi.
  • (6) R. C. Grabowski, I. G. Droppo, and G. Wharton, “Erodibility of cohesive sediment: The importance of sediment properties,” Earth-Science Reviews, vol. 105, no. 3–4, pp. 101 – 120, 2011.
  • (7) J. Sallès, J. Thovert, and P. Adler, “Deposition in porous media and clogging,” Chemical Engineering Science, vol. 48, no. 16, pp. 2839 – 2858, 1993.
  • (8) A. Bouddour, J. L. Auriault, and M. Mhamdi-Alaoui, “Erosion and deposition of solid particles in porous media: Homogenization analysis of a formation damage,” Transport in Porous Media, vol. 25, no. 2, pp. 121–146, 1996.
  • (9) K. Yamamoto, S. Satake, H. Yamashita, N. Takada, and M. Misawa, “Fluid simulation and x-ray ct images for soot deposition in a diesel filter,” The European physical journal special topics, vol. 171, no. 1, pp. 205–212, 2009.
  • (10) L. Talon, D. Bauer, N. Gland, S. Youssef, H. Auradou, and I. Ginzburg, “Assessment of the two relaxation time lattice-boltzmann scheme to simulate stokes flow in porous media,” Water Resources Research, vol. 48, no. 4, 2012.
  • (11) A. Mahadevan, A. Orpe, A. Kudrolli, and L. Mahadevan, “Flow-induced channelization in a porous medium,” EPL (Europhysics Letters), vol. 98, no. 5, p. 58003, 2012.
  • (12) A. Kudrolli and X. Clotet, “Evolution of porosity and channelization of an erosive medium driven by fluid flow,” Phys. Rev. Lett., vol. 117, p. 028001, Jul 2016.
  • (13) O. Johnsen, C. Chevalier, A. Lindner, R. Toussaint, E. Clément, K. J. Måløy, E. G. Flekkøy, and J. Schmittbuhl, “Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas,” Phys. Rev. E, vol. 78, p. 051302, Nov 2008.
  • (14) S. Bonelli, Erosion in Geomechanics Applied to Dams and Levees. Civil engineering and geomechanics series, Wiley, 2013.
  • (15) A. Alem, N.-D. Ahfir, A. Elkawafi, and H. Wang, “Hydraulic operating conditions and particle concentration effects on physical clogging of a porous medium,” Transport in Porous Media, vol. 106, no. 2, pp. 303–321, 2015.
  • (16) D. A. Wolf-Gladrow, Lattice-gas cellular automata and lattice Boltzmann models: an introduction. Springer Science & Business Media, 2000.
  • (17) L.-S. Luo, W. Liao, X. Chen, Y. Peng, and W. Zhang, “Numerics of the lattice boltzmann method: Effects of collision models on the lattice boltzmann simulations,” Phys. Rev. E, vol. 83, p. 056710, May 2011.
  • (18) O. Filippova and D. Hänel, “Grid refinement for lattice-bgk models,” Journal of Computational Physics, vol. 147, no. 1, pp. 219 – 228, 1998.
  • (19) M. Bouzidi, M. Firdaouss, and P. Lallemand, “Momentum transfer of a boltzmann-lattice fluid with boundaries,” Physics of Fluids, vol. 13, no. 11, pp. 3452–3459, 2001.
  • (20) R. Mei, W. Shyy, D. Yu, and L.-S. Luo, “Lattice boltzmann method for 3-d flows with curved boundary,” Journal of Computational Physics, vol. 161, no. 2, pp. 680 – 699, 2000.
  • (21) X. Shan, X.-F. Yuan, and H. Chen, “Kinetic theory representation of hydrodynamics: a way beyond the navier–stokes equation,” Journal of Fluid Mechanics, vol. 550, pp. 413–441, 3 2006.
  • (22) C. Pan, L.-S. Luo, and C. T. Miller, “An evaluation of lattice boltzmann schemes for porous medium flow simulation,” Computers & Fluids, vol. 35, no. 8–9, pp. 898 – 909, 2006. Proceedings of the First International Conference for Mesoscopic Methods in Engineering and Science.
  • (23) E. Fattahi, C. Waluga, B. I. Wohlmuth, U. Rüde, M. Manhart, and R. Helmig, “Pore-scale lattice boltzmann simulation of laminar and turbulent flow through a sphere pack,” CoRR, vol. abs/1508.02960, 2015.
  • (24) K. Yamamoto, K. Yamauchi, N. Takada, M. Misawa, H. Furutani, and O. Shinozaki, “Lattice boltzmann simulation on continuously regenerating diesel filter,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 369, no. 1945, pp. 2584–2591, 2011.
  • (25) X. Yin, G. Le, and J. Zhang, “Mass and momentum transfer across solid-fluid boundaries in the lattice-boltzmann method,” Phys. Rev. E, vol. 86, p. 026701, Aug 2012.
  • (26) S. P. Thampi, S. Ansumali, R. Adhikari, and S. Succi, “Isotropic discrete laplacian operators from lattice hydrodynamics,” Journal of Computational Physics, vol. 234, pp. 1 – 7, 2013.
  • (27) T. Wahl, “A comparison of the hole erosion test and jet erosion test,” in Proceedings of the Joint Federal Interagency Conference on Sedimentation and Hydrologic Modeling, pp. 1–11, 2010.
  • (28) Z. Guo, C. Zheng, and B. Shi, “Discrete lattice effects on the forcing term in the lattice boltzmann method,” Phys. Rev. E, vol. 65, p. 046308, Apr 2002.
  • (29) T. Seta, R. Rojas, K. Hayashi, and A. Tomiyama, “Implicit-correction-based immersed boundary–lattice boltzmann method with two relaxation times,” Physical Review E, vol. 89, no. 2, p. 023307, 2014.
  • (30) A. Kuzmin, Z. L. Guo, and A. A. Mohamad, “Simultaneous incorporation of mass and force terms in the multi-relaxation-time framework for lattice boltzmann schemes,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 369, no. 1944, pp. 2219–2227, 2011.
  • (31) R. Fell, C. F. Wan, J. Cyganiewicz, and M. Foster, “Time for development of internal erosion and piping in embankment dams,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 129, no. 4, pp. 307–314, 2003.
  • (32) A. D. Araújo, J. S. Andrade, and H. J. Herrmann, “Critical role of gravity in filters,” Phys. Rev. Lett., vol. 97, p. 138001, Sep 2006.
  • (33) J. Latt, “Choice of units in lattice boltzmann simulations,” Freely available online at http://lbmethod. org/_media/howtos: lbunits. pdf, 2008.
  • (34) A. Hart, “A review of technologies for transporting heavy crude oil and bitumen via pipelines,” Journal of Petroleum Exploration and Production Technology, vol. 4, no. 3, pp. 327–336, 2014.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description