Pore-scale simulations of solute transport in urban soil substrates hosting roots

Abstract

The transport of nutrients in urban soil substrates, such as green roofs, is largely dictated by the need of the plants themselves, via adsorption by their root networks. This process pulls nutrients such as fertilizer from the soil substrate where biochar can act as a reservoir, holding the fertilizer for longer periods. Biochar is sometimes added to soils to store excess fertilizer resulting in reduced loss of fertilizer and leaching of contaminants into the groundwater. The biochar then desorbs the fertilizer (solute) into the surrounding soil. Due to limited experimental insights on the adsorption-desorption system of biochar and plants roots; we outline a complete methodology by which root structures are generated and incorporated into a fully 3D geometry. The lattice Boltzmann method is used for pore-scale simulations to compare the tap and fibrous (grass) roots’ behavior. The results show a difference in the mean contaminant breakthrough as well as in homogeneity, with the tap root having a lower mean but higher variation of concentration at the outlet. This is due to its lower surface area, which is directly linked to its adsorption capacity and the consequent quantity of biochar needed in the substrate. The lower biochar presence is quantified by a large $ℓ_{ξ}$ , the average distance between desorbing particles. This distance dictates the mixing strength within the system. At higher values, we see a lower mean contaminant breakthrough and a higher statistical variation within the breakthrough, engendering higher contaminant peak values during rainfall events. The opposite is true for lower values of $ℓ_{ξ}$ , represented by the fibrous case, which has a higher mean breakthrough but less variation. This methodology allows for testing varying plant and fertilizer deployment methods in order to create some best practices for minimizing contaminant breakthrough without the need for in situ tuning.

Keywords

advection-diffusion-reaction lattice Boltzmann biochar plant root green roofs

Introduction

Urban flooding is an important consideration given the limited catchment within most urban environments and its visible impact during extreme rainfall events. These events have become more frequent and intense in northern latitudes (Dore, 2005; Zolina et al., 2010). Green roofs (vegetated roofing) are an attractive option for sustainable urban development and are used to reduce urban flooding, improve protection and thermal performance of structures, and aid in removing contaminants (Johansson et al., 2024; Twohig et al., 2022). Although relatively “thin” in comparison to naturally occurring vegetated areas, they can be designed to maximize efficiency in water detention via understanding the dominant physiochemical mechanisms involved. Fertilizers are typically added to the soil substrate to nourish the vegetation. However, as green roofs generally lack dedicated reservoirs, excessive fertilizer application can lead to leaching and contamination of runoff with undesirable chemicals, creating an environmental burden. As a mitigation measure, biochar can be incorporated into the soil substrate to retain excess fertilizers. However, its application is typically limited (a few percent by volume) due to high cost and observed competition for adsorption sites with plant roots.

Biochar is essentially any organic matter that has been carbonized via high-temperature pyrolysis. It can adsorb both organic and inorganic contaminants such as heavy metals, fertilizers, pharmaceuticals, and dyes, among other substances (Abbas et al., 2018; Qiu et al., 2022). It does this through chemical interactions, electrostatic interaction and ion exchange, and sorption processes such as partitioning, surface sorption, pore filling, and complexation. This broad range of processes makes its application attractive in many situations requiring the remediation of a variety of contaminants. Biochar can also play a role in obstructing pesticide uptake in plants from contaminated soils, serving as an alternative method to phytoremediation by fixing the pesticide residue on the biochar itself and protecting the plants (Yu et al., 2009).

The amendment of biochar to soil substrates can make nutrients more available to the plants, but its quantity must be optimized, which is a delicate task, as different plants have varying adsorption capacities (Hou et al., 2022). Experimental and numerical studies on storm water management through green roofs also explore the role biochar has on water retention, although without the inclusion of vegetation (Gan et al., 2021, 2025). A review performed in Australia summarized globally diverse studies of reported contamination levels using indicators such as total phosphorus (Alim et al., 2022), which demonstrates the importance of understanding the relationship between water retention in green roofs with biochar and the addition of vegetation. Most experimental studies, in the field or in laboratory experiments, are often limited to a macroscopic scale in regard to the plant uptake itself, removing any insight into the microscopic (or pore) scale mechanisms which determine the process evolution. In city-scale numerical studies, the water retention capacity of green roofs is often simplified, without accounting for the saturation level of the soil substrate (Wei et al., 2023). Pore-scale numerical simulations allow for a more detailed analysis of how the local heterogeneities of the flow, reactivity, and concentration field determine the overall system dynamics; allowing finding the correlations between micro-scale phenomena with macro-scale performance and generalizing the findings.

The general numerical model used to capture adsorption-desorption within porous media is the advection-diffusion-reaction (ADR) equation. As the name suggests, a scalar quantity such as a contaminant (or solute) is transported via advection and diffusion as well as accumulated via a source or sink term coupled to a reaction rate. This equation can be treated as a continuum (Li et al., 2023; Oliveira et al., 2020; Shafabakhsh et al., 2024), a single-field model (Girault et al., 2024), or using the discrete lattice Boltzmann method (LBM), which solves the Boltzmann equation of particle motion (Zhang et al., 2024), linking microscopic particle motion with macroscopic fluid behavior. A comprehensive review of the current state-of-the-art on applied techniques can be found in Baqer and Chen (2022). At the macroscopic level, the porous medium is defined by the scalar parameters porosity $ϵ$ and permeability $κ$ , with additional parameters for multiple fluids. Permeability, a measure of the porous matrix’s resistance to the infiltrating liquid, is difficult to determine spatially, as most porous media is not isotropic. Moreover, most continuum formulations of flow through porous media simplify the liquid infiltration by removing any inertial component, as in Darcy’s law. In principle, solute transport through porous media can be viewed as a combination of two dominant processes. The advection-dominant regime is where the transport occurs at the beginning of the infiltration process or when the flow domain is shorter than the characteristic time, and the solute distribution is highly heterogeneous. In this regime, the solute forms plumes (or lamellae), which infiltrate the less resistant parts of the pore network. When the flow domain is longer than the characteristic time, the diffusion-dominant regime acts to redistribute the solute more evenly throughout the pore network, leading to a more homogeneous distribution throughout the network (Bonazzi et al., 2023; Le Borgne et al., 2013; Lei and Luo, 2021; Puyguiraud et al., 2020). If one considers the inclusion of spatially varying reactivity such as in cases where plant roots and biochar simultaneously affect the flow and adsorption of fertilizers, these phenomena become more unpredictable still. Any approach which assumes spatial homogeneity of the interacting physics will by definition be unable to capture the local variations of velocity, pore topology, solute accumulation, and molecular diffusion necessary to accurately capture the interactions.

As mentioned above, pore-scale modeling aims to resolve the assumptions employed in a macroscopic description of the system by solving it directly at the pore scale. This approach brings its own challenges; such as if one chooses to utilize a continuum approach, the complex geometry alone poses a significant challenge when defining appropriate mesh resolution. One common approach is pore network models (PNMs) which define the porous medium as a series of tubes connected to reservoirs, defined using a variety of parameters, and which combined aim to mimic the material permeability (Joekar-Niasar and Hassanizadeh, 2012; Zhao et al., 2024). Alternatively, the LBM has been utilized extensively in solving flows in porous media as well as the ADR equation. Ju et al. utilized the color gradient LBM to study the impact of pore morphology on immiscible fluid infiltration of an inhomogeneous porous domain (Ju et al., 2022). Su Yan has proposed a modification to LBM by applying nonlinear drag models relating to the mesoscale porosities and velocities to more accurately capture the porous structure evolution (Su, 2024). Yu et al. introduced a hybrid immersed-boundary LBM which utilizes a finite difference scheme for calculation of reactivity within the solid matrix itself (Yu et al., 2019). Even though there are many studies utilizing or expanding LBM for handling ADR there are very limited numerical investigations on the adsorption dynamics of plants. One such promising approach was taken by Koch et al. (2018), wherein adsorption across the plant membrane as well as internal transport within the plant itself is calculated, though the model remains biologically simplistic.

Significant experimental and numerical work has been done on the transport and uptake of nutrients by root systems of plants. It is worth noting that nutrients are typically added as granules to soil substrates, and plant roots can only absorb them once they have dissolved in water. Many approaches have been taken to model fertilizer release into infiltrating water, such as the various mathematical models for controlled-release fertilizers CRFs (Irfan et al., 2018). These models utilize a variety of mass conservation equations for normal and extended release fertilizers. Empirical models for a variety of nutrients are also used, as the individual behavior of each unique compound can vary significantly from a chemical reaction perspective. Once the fertilizer is released, many transport models have also been applied. Continuum models are usually coupled to plant uptake models, taking the form of intermittent or constant flux boundary conditions (Rein et al., 2011). The fertilizer source in these uptake simulations consists of bulk sources within the soil itself, which is appropriate for larger scale models aimed at determining the effect on entire catchment systems. However, these models do not explore any localized effect of fertilizer placement within the soil and the interplay between the root structures and this heterogeneous solute distribution (Elasbah et al., 2019). Genetic algorithms have also been applied to model bulk plant uptake and have been extended to take into account vegetative competition as well (Cropper and Comerford, 2005). The transport models may also be coupled to plant growth models, such as in Rein et al. (2011), which included growth, transpiration, translocation of the nutrients within the plant itself, and a variety of biological processes such as aging, defoliation, and leaching. A strong overview of the processes which require consideration and the simplifications in use can be found in Rengel (1993). From the purely plant root morphological side, several studies such as Barber and Silberbush (1984) and Cruz et al. (2004) showed that root morphology does have an impact on plant nutrient uptake, and that plants will grow their root systems in such a way as to maximize nutrient uptake, particularly when nutrients are scarce (heterogeneous nutrient allocation in the soil). To the knowledge of the authors, there are no models that look at the micro-scale desorption and adsorption process using a biochar release mechanism and “realistic” root system, respectively, and how proximity and root morphology affect plant nutrient uptake, which is relevant to reducing breakthrough from urban green roofing systems.

To address these complex processes, this work gives a framework by which pore-scale simulations can be undertaken to analyze mass transport within an urban soil substrate such as a green roof. In the cases included here we look at the interaction between randomly dispersed particle desorbers resembling biochar and an adsorbing root-like structure under extreme rainfall conditions (Pe $\approx 50$ , 180 mm/h). Such a high intensity rainfall is chosen because of their propensity to generate the most fertilizer runoff. Different root morphologies are created, a taproot and a thinner fibrous root system. Solute is desorbed from the biochar particles and adsorbed by the root structure at a constant rate under a steady flow profile. We fix the desorption/adsorption rate in order to isolate the dominant physical mechanisms within this system. The resultant concentration field is evaluated alongside breakthrough and a statistical analysis of the desorbing biochar particles’ effect on the system is undertaken. A parameter $ℓ_{ξ}$ is proposed to quantify the relationship between biochar particle distance on plume formation and homogeneity of the breakthrough as measured by the concentration moments. It is shown that when biochar particles are closer to each other, they promote a more spatially even distribution of fertilizer runoff, given the flow conditions and equal reactivity rates. This is the first time a link between the deployment of biochar in vegetated green roofs and breakthrough patterns has been undertaken, to the authors’ knowledge. Moreover, by testing different root configurations we show that one can design a green roof for differing breakthrough patterns as a result of plant choice and biochar distribution. These findings can be used to more accurately prescribe fertilizer dosage given vegetation type and biochar quantity and reduce contaminant breakthrough.

Materials and methods

Lattice Boltzmann method

In this work we solve the advection-diffusion-reaction (ADR) equation

\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c) - \nabla \cdot (u c) + R,

(1)

where $c$ is a scalar (in our case fertilizer concentration), $t$ is time, $D$ is molecular diffusion, $u$ is the advecting fluid velocity, and $R$ is a sink/source term related to reactive processes. At the pore scale this equation can be represented using the discrete Boltzmann equation for particle motion using LBM.

LBM is an ideal choice for solving our root-biochar transport problem due to the complex geometry involved and it allows for full field solutions of the dependent variables to be extracted at the pore scale. The system is solved on a regular lattice structure wherein each element consists of a centroid and nodes placed on a cubic convex shell. The fictive particles travel along the lattice nodes governed by probabilities appropriate for the chosen lattice geometry such that the macroscopic properties of the fluid are preserved (Succi, 2001). A 3D regular cubic lattice with 19 degrees of freedom (D3Q19) is used and the solved equation is of the form

f_{r} (x + c_{r} Δ t, t + Δ t) - f_{r} (x, t) = - τ^{- 1} (f_{r} (x, t) - f_{r}^{eq} (x, t)) + F_{r}

(2)

where $f_{r} (x, t)$ is the distribution function at position $x$ and time $t$ along the $r$ -th direction; $c_{r}$ is the so-called discrete velocity vector along the $r$ -th direction over time interval $Δ t$ ; $f_{r}^{eq}$ is the equilibrium distribution function; and $τ$ is the mean collision time and is related to kinematic viscosity by $ν = c_{s}^{2} (τ - 0.5 Δ t)$ . The equilibrium distribution function $f_{r}^{eq} (x, t)$ takes the form

f_{r}^{eq} = w_{r} ρ (1 + \frac{c_{r} \cdot u}{c_{s}^{2}} + \frac{{(c_{r} \cdot u)}^{2}}{2 c_{s}^{4}} - \frac{u \cdot u}{2 c_{s}^{2}}), r = 1 - 19

(3)

where $c_{r} = 0$ when $r = 1$ . The constants $w_{r}$ are scalar weighting parameters for the D3Q19 lattice specifically valued to recover the Navier-Stokes system at the macroscopic scale. Additionally, $ρ$ is the density; $c_{s}$ is the speed of sound; and $u$ is the velocity used for defining the equilibrium distribution functions. This velocity can differ from the fluid hydrodynamic velocity, on the basis of the specific forcing scheme used. Here, the macroscopic flow quantities density and velocity, $(ρ, u)$ are related to the hydrodynamic moments by:

ρ = \sum_{r} f_{r},

(4)

ρ u = \sum_{r} c_{r} f_{r} + \frac{Δ t}{2} (\frac{Δ P}{L}) .

(5)

The body force $F_{r}$ as formulated by Guo et al. (2002), mimics the flow rate intensity during a rain event in our case, and is applied to the fluid.

F_{r} = (1 - \frac{1}{2 τ}) w_{r} (\frac{c_{r} - u}{c_{s}^{2}} + \frac{c_{r} \cdot u}{c_{s}^{4}} c_{r}) (\frac{Δ P}{L}),

(6)

where $Δ P / L$ is the pressure gradient. It is fitted to induce the desired rainfall intensity which has Peclet number, Pe $\approx 50$ , where domain length $L$ and molecular diffusion $D$ are known. We convert Pe to horizontal rainfall intensity by

Pe D / d = 50 \times 10^{- 9} / 1 \times 10^{- 3} = 5 \times 10^{- 5} [m / s],

(7)

which can be converted easily to (mm/h), resulting in 180 mm/h, meaning the rainfall event is considered extreme or that channeling of runoff is taking place, a very possible scenario. For reference, the global record for rainfall intensity is currently 305 mm/h (WMO, 2025). The resultant mean flow velocity is averaged over the porous domain, given by $U = (\int_{V_{por}} u_{z} d V_{por}) / V_{por}$ . We are most interested in the limit of extreme rainfall events, where advection is comparable with adsorption, a relevant scenario for contaminants

Adsorption-desorption implementation

A desorption rate is assigned to a fraction of randomly chosen particles such that their surfaces $S$ follow the first-order kinetics given by:

{- \frac{\partial c^{*}}{\partial λ_{s}^{*}} |}_{S} = Da (c^{*} - c_{eq}) |_{S},

(8)

where Da is the Damköhler number, $c^{*} = c (x, t) / c_{0}$ is the dimensionless concentration at position $x = (x, y, z)$ , $c_{eq}$ is locally defined as 0 (adsorbing) or 1 (desorbing), and $λ_{S}^{*} = λ_{S} / d$ is the dimensionless direction pointing inward to the particle surface. The Damköhler number, a ratio of the reaction and diffusion rate, is given by

Da = k_{vol} d^{2} / D,

(9)

where $k_{vol} = kS / V_{root}$ is the measured adsorption rate on the particle surface, $d$ is the particle diameter, $D$ is the molecular diffusion rate, and $S, V_{root}$ are the root surface area and volume, respectively. We use the term “measured” for $k_{vol}$ because the rates used in this study are directly comparable to the values observed in biochar adsorption experiments (Pettersson et al., 2025). The previous study used only spherical reactive particles, thus a fixed specific surface. In this work the root specific surface varies by root type and the values of $k$ used in this work are chosen such that Da $\approx 1$ .

The concentration field $c (x, t)$ is calculated using a secondary population $g_{r}$ , which is advected by the fluid velocity $u$ , the steady-state solution of a prior simulation. The coupling is one-way, where the equilibrium function of the secondary population uses the macroscopic velocity $u$ , and the body force still acts. We employ a simple one-way coupling as the advected solute is assumed to not disturb the flow. The full system is

g_{r} (x + c_{r} Δ t, t + Δ t) - g_{r} (x, t) = - τ_{c}^{- 1} (g_{r} (x, t) - g_{r}^{eq} (x, t)) + G_{r},

(10)

g_{r}^{eq} = w_{r} c (1 + \frac{c_{r} \cdot u}{c_{s}^{2}} + \frac{{(c_{r} \cdot u)}^{2}}{2 c_{s}^{4}} - \frac{u \cdot u}{2 c_{s}^{2}}), r = 1 - 19,

(11)

G_{r} = (1 - \frac{1}{2 τ_{c}}) w_{r} (\frac{c_{r} - u}{c_{s}^{2}} + \frac{c_{r} \cdot u}{c_{s}^{4}} c_{r}) (\frac{Δ P}{L}),

(12)

where the diffusion is related to the scalar collision time $τ_{c}$ by $D = c_{s}^{2} (τ_{c} - 0.5 Δ t)$ . This second lattice population gives the local concentration field by

c (x, t) = \sum_{r} g_{r} (x, t) .

(13)

We impose a Neumann boundary condition at the adsorbing surfaces for the scalar $g_{r}$ and the distribution function at a fluid node $x$ in proximity of the adsorbing surface placed at $(x - c_{r})$ is corrected along the wall-normal direction $r$ as

g_{r} (x, t + 1) = \frac{- A_{1} + A_{2}}{A_{1} + A_{2}} g_{r} (x, t) + \frac{2 w_{r} A_{3}}{A_{1} + A_{2}},

(14)

where $A_{1} = k$ , $A_{2} = D$ , $A_{3} = 0$ , and $w_{r}$ are the same lattice weighting parameters used in equation (3). This condition was developed by Huang et al. (Huang et al., 2016; Huang and Yong, 2015).

As mentioned above, we fix the Damköhler number, $Da \approx 1$ , to ensure the balance of flow and reaction rate. This is based upon experiments and findings in our previous work (Pettersson et al., 2025), where the Damköhler number is calculated to be near unity for methylene blue adsorption by biochar particles. By combining both Da and Pe, we assume a realistic system in which the biochar desorbs into the advecting flow and is subsequently adsorbed by the plant root structures.

Root generation algorithm

Figure 1 displays the characteristics and geometries of the generated root systems used in the is work. In Figure 1(a) the averaged cross-sectional orientations are plotted, divided into $n_{i}$ where $i = z$ is the longitudinal facing (full line) and $i = x$ is the transverse facing (dashed line). The values are calculated as a sum of nodes with outward normal in the aforementioned orientation divided by the total nodes on the interface. This is calculated at depths of the domain measured by particle diameters and the sum of both orientations is unity. It is clear that as the root goes deeper into the medium the tap root maintains an even distribution of surface orientations, however the fibrous root displays a heavy bias toward transversal growth.

Figure 1.

(a) Cross-sectional averaged surface node orientations divided by total cross-sectional surface area, as a function of depth by particle diameters. Red, black lines denote fibrous, tap roots, respectively. Solid: longitudinal orientation, dash: transverse. (b) Fibrous root. (c) Tap root.

The root generation utilizes the space colonization algorithm (SCA) as per Runions et al. in 2007 (Runions et al., 2007) following their 2005 publication, which introduced a biologically-motivated algorithm for generating leaf patterns in 2D (Runions et al., 2005). The later paper includes an extended algorithm for 3D and treats competition for space as the main factor for the branching structure of trees. We briefly explain how the algorithm functions for the reader’s benefit.

A singular tree node was placed at the center of the upper porous boundary. It is from the tree node that the initial branches start growing. A number of attraction points $n$ are then randomly generated within a volume whichs depends on the shape being imitated, in this case the shape is cylindrical in nature. The cylindrical shape is defined by an inner and outer radius, $r_{i}, r_{o}$ , respectively, and a maximum length $z$ . The attraction points signal the algorithm that space is available for the root to grow in the point’s direction. Attraction points that are within a distance $d_{r}$ of the tree node, called the radius of influence, will influence its growth. Suppose we have a few such attraction points inside the radius of influence of the tree node; then a new tree node is created at an input distance $D_{i}$ from the initial tree node. The direction of the new tree node is determined by the average of the normalized vectors of the attraction points within the radius of influence. The line that’s drawn from the initial tree node to the new tree node becomes a branch.

One can alter this direction of growth by putting directional bias on it, referred to as tropisms. We have implemented gravitropism $g$ , which forces the branches vertically downwards, and a “horizontal” tropism $h$ , which forces the branches laterally to encourage a more spread-out growth. We set the strength of the tropism itself as well as a rate of change with root depth ( $k_{g}$ and $k_{h}$ ). One can see in Figure 1 that the fibrous root has more lateral branching closer to the surface, and more vertical branching deeper into the porous medium. The higher the value of $k_{g, h}$ , the faster the tropism influence disappears. The tropisms were necessary for generating the fibrous root and ensuring a nice transition from lateral to vertical growth.

This process continues iteratively until at least one tree node is closer to an attraction point than a user-defined kill distance $d_{k}$ , at which point the attraction point is removed. The algorithm terminates when a chosen number of iterations is reached or when all attraction points are removed.

In the taproot case the primary root is created separately by assuming a linear change in primary root thickness with depth, and then some random perturbations are added laterally. The “strength” of the lateral perturbations are set by $δ x, δ y$ , and determine the maximum possible lateral perturbations, and the change of root thickness with depth is determined using the rate $s$ . We then superimpose the root colonization algorithm on this primary root by sampling attraction points from a cylindrical sample space around the primary structure, not allowing overlapping structures. The thickness of each branch is determined by the root structure and the number of child branches each tree node has. An extended form of Bresenham’s line algorithm is used to select which lattice points in the voxelized domain should be selected to approximate a straight line between the two nodes of the root structure. We also impose a threshold on the minimum and maximum radius of each branch, $r_{\min}, r_{\max}$ , respectively. Table 1 gives the algorithm input parameters for each root geometry.

Table 1.

Root algorithm parameters, where $r_{i}, r_{o}$ have units (mm), and $r_{\min}, r_{\min}, δ x / δ y$ have [l.u.], which are lattice units.

	$n$	$r_{i}$	$r_{o}$	$d_{r}$	$d_{k}$	$D_{i}$	$r_{\min}$	$r_{\max}$	$g$	$h$	$k_{g}$	$k_{h}$	$δ x / δ y$	$s$
Taproot	250	0	10	100	1	2	5	15	0.4	0.8	2	4	–	–
Fibrous	40	5	10	100	1	2	5	15	0.4	0.4	0	0	15.15	0.08

Packed bed generation and merging

A randomly packed bed was generated using Blender. An array of monodisperse spheres with a particle diameter of $1.0$ mm were initialized above a large cylindrical receptacle and allowed to drop with applied gravity and rigid body collision. Once the simulation was completed, the coordinates and sphere radii were exported.

The packed bed and root domains were merged using a different algorithm, in which the user first defines a rectangular subdomain within the cylindrical sphere packing, removing any wall effects on the packed bed. The subdomain sizes are chosen large enough to host root lateral spreading and allow for the investigation of contaminant transport beyond the root vertical elongation. The root geometry was then set to replace any spheres that conflict with the root geometry. An additional buffer zone with porosity 1 is added to the upstream end of the domain. Table 2 gives the domain geometry detailsdetails and the domain is shown in Figure 2.

Figure 2.

Domain schematic, where $L$ is the total domain length. Buffer volume prevents problems with periodic streamwise boundary condition. Flow Inlet and Outlet – Periodic boundary. Flow Lateral boundaries – Free slip. Only root structure and desorbing particles pictured. Re $= 0.08$ , Pe $= 50$ , Da $\approx 1.0$ .

Table 2.

Domain parameters. Re $= 0.08$ , Pe $= 50$ , Da $\approx 1$ . $k = 3.3 \times 10^{- 4} [s^{- 1}]$ .

	$NX$	$NY$	$NZ (L)$	$N_{TOTAL}$	$Δ X, Δ Y, Δ Z$	$d$
Lattice	256	256	768(1024)	67,108,864	1	22
Physical (mm)	11.6	11.6	35.0/46.6	–	0.045	1.0

Table 3 displays the root morphological differences between the cases. Note that while volume is similar in both cases, the adsorbing surface area varies significantly. The number of desorbing particles $n_{p, des}$ is chosen such that the total reactive surface area of the particle desorbers matches that of the root adsorber in each case, to keep Da $\approx 1$ . The difference in adsorbing and desorbing surface area is less than 4% for both cases, thus we consider them of equal influence within the system.

Table 3.

Parameters of the roots and adsorbing biochar, Damköhler number Da $\approx 1$ , $δ_{S}$ is the difference in surface area between adsorbers and desorbers.

	$S$	$V_{root}$	$S / V_{root}$	$δ_{S}$	$n_{p, des}$
Taproot	141,896	1,527,495	0.09	+2.26%	99
Fibrous root	234,134	1,525,486	0.15	+3.90%	166

Operating conditions

An initial concentration $c_{0} = 0$ is applied everywhere. The roots and particles are treated as solids with infinite capacity. Thus, desorbers release the solute at a constant rate into the domain from the particle surface for the duration of the simulation. The reverse occurs with the adsorbing roots, the reaction occurs at the root surface and the solute is removed from the domain at a constant rate. No flow or reaction is calculated in the root or particle interiors. This is not reflective of the real world, wherein nutrients absorbed by the plants’ roots are transported away from the root system. This process is difficult to emulate with any accuracy as the plant may remove the nutrients from the soil as needed and may not even store them locally in the roots, further complicating this model. The desorbing particles would also have a maximum release capacity however since we are interested in how the reactive element proximity and location within the packing affect the breakthrough we consider our simplification a logical starting point.

The solute distribution evolves in a steady flow field with a Re $= 0.08$ . The velocity streamwise boundaries are assigned a periodic condition, such that any solute flowing into the outlet, is returned to the inlet. We impose a buffer zone at the inlet boundary to prevent any effect on the transport within the porous zone, our area of interest. The lateral boundaries are assigned a free-slip condition. The streamwise scalar boundary conditions are $c = 0$ at the inlet and $\nabla c \cdot n = 0$ at the outlet. The lateral scalar boundaries are assigned the symmetry condition.

Model validation

Figure 3 displays the validation for flow between an array of periodic cubic centered spheres and includes a grid convergence study. For a complete validation of the passive scalar reactivity model, we refer the reader to the previous work by Maggiolo et al. (2023) which includes a full validation in the Appendix.

Figure 3.

Flow validation: (a) $F_{d}^{*}$ is dimensionless drag as the ratio of the drag over Stokes drag. The line shows the solution of Sangani and Acrivos (1982), with $ϵ$ as the volume fraction. Inset: Grid convergence, shows minimum number of nodes per particle diameter to be eight, as per Ladd (1994) and (b) Solute flow rates for validation: void marks are equation (16), filled are equation (15). Upper panel – Fibrous root. Lower panel – Tap root. $q^{*}$ is the flux rate, $k_{vol} t$ is dimensionless time.

In Figure 3(a) we plot $F_{d}^{*}$ , the dimensionless drag as a function of the volume fraction. The line shows the solution of Sangani and Acrivos (1982) for slow flow past an array of spheres. We see a good agreement as the volume fraction $ϵ$ is increased from 0.56 (the minimum possible with the sphere arrangement) to 1.0, thus we are satisfied that the numerical model can capture flow within a porous environment such as a packed bed. The inset in Figure 3(a) displays a grid convergence study undertaken on the system. The results show a required minimum number of nodes per particle diameter to be around 8–10, as per Ladd (1994); our particle size is fixed at 22 nodes per particle. Figure 3(b) plots the mass balance of the solute given by

q_{reac} = k (n_{p} π d^{2} (1 - c_{1})) - S_{root} c_{2}),

(15)

the desorption/adsorption flux difference, where $c_{1}, c_{2}$ are the mean concentration at the desorbing/adsorbing surfaces, respectively. The outlet mass flux is

q_{out} = \int_{A_{out}} (u_{out} c) dA,

(16)

where $A_{out}$ is the outlet cross-sectional area and $u$ the local flow velocity. The upper panel displays these rates for the fibrous root, and the tap root (lower panel). The fluxes $q^{*}$ are plotted as a function of dimensionless time $k_{v o l} t$ based upon the adsorption rate. The small difference at long times between equations (15) and (16) is likely due to the estimation of the surface concentration, where an underestimation of $c_{1}$ gives a slightly larger $q_{reac}$ value. We approach a steady state wherein both rates match, indicating mass is conserved correctly.

Results and discussion

A qualitative representation of the dynamics within our system is shown in Figure 4; where the fibrous root is shown on the left (a), and the tap root on the right (b). The dimensionless concentration $c / c_{1}$ is shown, ranging from 0 (green) to 1 (red/magenta), indicating the mean desorber surface concentration. One can see the particles as transparent outlines, along with the structure of the root surface which is represented by the solid light green color.

Figure 4.

Root geometries: (a) fibrous root and (b) taproot. Dimensionless concentration $c / c_{1}$ is shown where $c_{1}$ is the mean concentration at the desorbing surface. The root volume is shown in light green. Particle volumes are transparent.

A more abrupt change (seen as a magenta/white border rather than a gradient from magenta to white, e.g.) signifies a high concentration gradient. In the tap root case (right panel) the magenta has a more gradual gradient and is more spread in terms of area, particularly in the longitudinal direction. We also note that longitudinal stretches of low concentration (denoted by the green plumes) are more visible for the tap root in this particular slice, implying the potential of more efficient solute removal in this case. These observations may be due to the location of some desorbing particles, which are likely responsible for localized pockets of higher concentration below the root structures. Another cause of these pockets is the presence of concentration plumes extending from the desorbers downstream. It is also the case that these plumes can coalesce, leading to larger areas of concentration at the outlet visible in both cases. A deeper analysis is necessary to determine what is occurring within the concentration field, particularly in relation to the location of the random particle desorbers and the plumes extending therefrom.

In order to quantify the existence of plumes, we examine the concentration differences in pores between the desorbing particles and the global mean value. By examining these fluctuation correlations we can determine if a pore upstream has an effect on those downstream, essentially the formation of a longitudinal plume. Figure 5 gives the average fluctuation correlations between the concentration in any given pair of pores within a prescribed downstream search angle $α$ , as shown in Figure 5(a). We define fluctuation as the difference between the mean concentration in the $i$ th pore and the average concentration among all pores that belong to the same volumetric slice of thickness $d$ , such that $c_{i}^{'} = c_{i} - 〈 c_{i} 〉$ . The correlations $c_{i}^{'} c_{j}^{'}$ are then the ensemble average among all pairs of pores $i, j$ that satisfy the criterion of the search angle $α$ , separated by longitudinal distance $z_{i} - z_{j}$ , normalized by the average distance between desorbers $ℓ_{ξ}$ . To identify the pore space, Delauney tesellation, is used. It is a method of partitioning volumes between packed spheres using tetrahedrons. A complete description of this technique can be found in Ben-Noah et al. (2025).

Figure 5.

(a) Visualization of Delauney tesselation used to calculate the pores between desorbing particles. $i, j$ are pores, $α$ is the search angle, (b) pore concentration $c / c_{1}$ , the spheres are the pores not the solid particles in this representation, and (c) averaged correlation of pore concentrations as a function of normalized distance between pores; normalized by $ℓ_{ξ}$ , the average distance between desorbing particles. The correlation is normalized by the variance. Search angles $α$ range from $5^{°}$ to $45^{°}$ . Inset: Correlations with distance $z_{i} - z_{j} = 5 ℓ_{ξ}$ .

In Figure 5(a) we see the partitioning of pores $i$ and $j$ which are defined by a tetrahedron. We have visualized the pores in Figure 5(b) as spheres, colored by their deviation from the mean solute concentration at the desorbing surface $c_{i} / c_{1}$ . We stress that these spheres do not correspond to the solid particles, but rather the pore space between them. It is clear that the downstream effect does exist, with clear vertical structures of concentration away from the mean solute concentration.

Figure 5(c) quantifies the mean correlation of pore concentration, normalized by the variance, as a function of longitudinal distance between any two pores, as well as normalized by the average distance between desorbers. Several search angles are used, ranging from $α = 5 - 45^{°}$ . The upper panel displays the correlations for the fibrous case with the lower for the tap root. The insets show corrrelation values for a fixed distance $z_{i} - z_{j} = 5 ℓ_{ξ}$ . We see a high correlation between the concentrations, even at distances greater than the average between two desorbers, which follows since the Peclet number is much higher than the Damköhler number. This is true for all search angles, with an increasing correlation as the search angle is restricted; meaning a dominant longitudinal correlation, or a greater advective effect. If one examines the inserts, we see that when the distance between slices is fixed, a faster decorrelation is apparent for search angles $α > 30^{°}$ . This is quantitative confirmation that plume formation is occurring in both root cases, with the tap root configuration exhibiting slightly higher correlation, but qualitatively both curves show a very similar trend. This is in line with the structure of the low concentration plume formations visible in Figure 4. Once the distances increases, we see a decrease in correlation until it decorrelates at six average desorber distances. This suggests there is a maximum depth after which plumes cease to have any relevance and molecular diffusion distributes the solute more evenly within the porous volume. The thinner the green roof (or less biochar applied), the more the plume formation will affect the breakthrough distribution. The choice of plant type can dictate the depth of the substrate or vice versa, given the application of biochar as a consequence of either choice.

Figure 6 displays the probability distribution functions (PDFs) for three different locations within the domain. Figure 6(a) is calculated at the outlet over a small volume with depth of one particle diameter $d$ . Figure 6(b) and (c) are calculated on the surface of the adsorbing root structures and the desorbing biochar particles, respectively. All concentrations are normalized by the maximum concentration $c_{\max} = 1$ .

Figure 6.

PDFs of normalized concentration at a small outlet volume (a), on the adsorber surface (b), and on the desorber surface (c). The tighter fibrous root curve in (a) implies a more homogeneous breakthrough front.

The surface behaviors for both cases are similar in Figure 6(b) and (c), with similar distribution curves. The concentration magnitudes also span the same range, though the peaks differ, with maxima at 0.01 and 0.1 for adsorbers and desorbers, respectively. Figure 6(a), however, shows a clear difference in the outlet concentration profile between the fibrous and tap roots. The fibrous root has a tighter curve, implying a more homogeneous outlet concentration front, or breakthrough. This would indicate plumes with similar concentration. Since Pe $>>$ Da, meaning the advective forces dominate those of the reaction, we expect to see plume formation. This could take the form of isolated plumes with dissimilar concentrations or many plumes in close proximity with more equal concentrations. Additionally, the Pe and Da values are fixed for both root geometries, thus any difference in breakthrough concentration will not be due to this imbalance. The concentration magnitudes themselves are not dissimilar, with an almost identical maximum value around 0.03.

In order to look at the plume spatial distribution and its contribution to the concentration distribution, the evolution of the normalized concentration variance as a function of depth in the porous medium is plotted in Figure 7.

\frac{σ_{c}^{2} (z)}{μ_{c}^{2}} = \frac{1}{μ_{c}^{2}} \int_{A_{0} d} {(c - μ_{c} (z))}^{2} d A_{0} d,

(17)

where $σ_{c}^{2} (z)$ is calculated over a transverse slice $A_{0}$ with longitudinal thickness $d$ , $μ_{c} (z)$ is the volumetric slice mean concentration, $μ_{c}$ is the global mean concentration.

Figure 7.

Normalized variance as a function of normalized depth. Blue line represents analytical relationship between the quantities.

We apply the random aggregation model (Duplat and Villermaux, 2008) to show how the overlap of plumes gives the relationship in Figure 7. This model states that you assume each concentration field in every plume to be identically distributed and that overlaps occur randomly in space. This allows for the calculation of the global variance as the variance obtained from the random sum of identically distributed concentrations; and as a consequence of the random plume overlaps.

We start with defining that each plume diffuses in time and space such that its transverse spread $s_{p}$ is treated as a 1D Fickian diffuse system with the standard proportionality characterization

s_{p} = \sqrt{Dt} \propto t^{1 / 2},

(18)

and in time $t$ it travels the longitudinal distance $Δ z = Ut$ from its point of origin (a desorbing particle), where $Δ z = z - z_{ξ}$ , with $z_{ξ}$ the vertical position of the desorbing particle. Substituting this into equation (18) gives

s_{p} (z) = \sqrt{D Δ z / U} .

(19)

We can then define the ensemble average of the transversal spread of all plumes at a given depth $z$ as

〈 s_{p} 〉_{p} (z) = \frac{1}{n_{p} (z)} \sum_{j = 1}^{n_{p} (z)} s_{p, j} (z),

(20)

s_{p, j} (z) = \sqrt{D (z - z_{ξ, j}) / U},

(21)

where $n_{p} (z) = z / ℓ_{ξ}$ . $n_{p} (z)$ is the number of plumes present at depth $z$ in a cross-section volume $ℓ_{ξ}^{2}$ and $ℓ_{ξ}$ is the average distance between desorbers. $s_{p, j}$ is the spread of the $j^{th}$ plume. We convert this discrete formulation to a continuum formulation by defining a volume with length $L$ and cross-sectional area $ℓ_{ξ}^{2}$ within which we wish to calculate the same ensemble average as in equation (20).

〈 s_{p} 〉_{p} = \frac{1}{n_{p}} \int_{0}^{n_{p} = z / ℓ_{ξ}} P_{ξ} s_{p} (z) d (\frac{z}{ℓ_{ξ}}),

(22)

where we introduce the probability $P_{ξ}$ of a plume existing at depth $z$ . This probability is defined by the average distance between desorbers such that $P_{ξ} = 1 / ℓ_{ξ}$ . We then rewrite equation (19) in terms of $z / ℓ_{ξ}$

s_{p} (z / ℓ_{ξ}) = {(\frac{D}{U})}^{1 / 2} {(\frac{Δ z}{ℓ_{ξ}})}^{1 / 2} \frac{1}{ℓ_{ξ}^{1 / 2}}

(23)

Utilizing our definition $Δ z = z - z_{ξ}$ we employ a change of variable to the integral in equation (22) such that

d (\frac{z_{ξ}}{ℓ_{ξ}}) = - \frac{d (Δ z)}{ℓ_{ξ}},

(24)

z_{ξ} = 0 \Rightarrow Δ z = z,

(25)

\frac{z_{ξ}}{ℓ_{ξ}} = \frac{z}{ℓ_{ξ}} \Rightarrow Δ z = 0 .

(26)

Equation (22) is rewritten as

〈 s_{p} 〉_{p} = \frac{1}{n_{p}} \int_{z / ℓ_{ξ}}^{0} P_{ξ} s_{p} (\frac{Δ z}{ℓ_{ξ}}) [- d \frac{Δ z}{ℓ_{ξ}}] .

(27)

Finally, by combining equation (23) and rewriting the integral in equation (27) we can write the plume transversal spread ensemble average as

〈 s_{p} 〉_{p} = \frac{1}{n_{p}} \int_{0}^{z / ℓ_{ξ}} P_{ξ} {(\frac{D}{U})}^{1 / 2} {(\frac{Δ z}{ℓ_{ξ}})}^{1 / 2} \frac{1}{ℓ_{ξ}^{1 / 2}} d (\frac{Δ z}{ℓ_{ξ}}),

(28)

which is integrated over the number of plumes at a given depth $z$ , not over the depth itself. This integral can be solved directly, giving

〈 s_{p} 〉_{p} = \frac{1}{n_{p}} \frac{1}{ℓ_{ξ}^{1 / 2}} {(\frac{D}{U})}^{1 / 2} \frac{z^{3 / 2}}{ℓ_{ξ}^{3 / 2}} .

(29)

We model the blending of the plumes at a certain depth $z$ as the random aggregation of plumes with maximum concentration

〈 c_{\max} 〉_{p} \propto 〈 s_{p}^{- 2} 〉_{p},

(30)

and with a variance

σ_{p}^{2} \propto {〈 s_{p}^{- 4} 〉}_{p} .

(31)

We expect the variance of a single isolated plume to approximately scale inversely with the square of the maximum concentration (Heyman et al., 2021), giving the relationship in equation (31). The mean transversal spacing between the plumes $L_{p}^{2} (z)$ , at depth $z$ , can be estimated by

L_{p}^{2} \propto ℓ_{ξ}^{2} / n_{p},

(32)

The number of overlaps $N$ within this area can be estimated by the spread of any plume divided by the spacing between plumes such that

N \propto s_{p}^{2} / L_{p}^{2} \propto s_{p}^{2} n_{p}^{2} .

(33)

As mentioned above, we assume that each plume’s concentration is described by a small, identical, Gaussian distribution when compared to the total physical volume in which the plume itself is located (Heyman et al., 2021). This is a very strong assumption, however it is extremely difficult to calculate a PDF for the plume structures we are attempting to isolate, thus we are forced to use a model that attempts to capture the PDF via aggregation of smaller, individual plumes. The random aggregation model states that the global variance is obtained from the random sum of identically distributed concentrations of randomly overlapping plumes, that is $σ_{c}^{2} = N σ_{p}^{2}$ . The global mean concentration increases with $n_{p}$ such that $μ_{c} \propto n_{p}$ , irrespective of the random aggregation scenario. One can write the relation between the concentration variance and the depth by inputting equations (29), (31), and (33) such that the global variance is proportional to

\frac{σ_{c}^{2}}{μ_{c}^{2}} = \frac{σ_{p}^{2} N}{μ_{c}^{2}} \propto \frac{〈 s_{p}^{- 4} 〉_{p} 〈 s_{p}^{2} 〉_{p}}{L_{p}^{2} n_{p}^{2}} = \frac{1}{ℓ_{ξ}^{2}} \frac{ℓ_{ξ}^{3}}{z^{3}} \frac{ℓ_{ξ}}{1} \frac{z^{2}}{ℓ_{ξ}^{2}} \frac{ℓ_{ξ}}{z} = \frac{ℓ_{ξ}}{z^{2}} .

(34)

The characterization above is valid when diffusion is the main mechanism for plume spread and when inertial effects are negligible such as with low Re. High values of Da can also invalidate this equation. A very low Da can lead to the possible occurrence of advective mixing, particularly for low concentration in high velocity pores, whereas a higher Da will serve to homogenize the concentration by suppressing mixing effects.

The blue line represents the analytical relationship between the moments of the PDF and the evolution of $z / ℓ_{ξ}$ as given in equation (34).

In Figure 7 all cases scale very similarly to the analytical relationship between the moments of plumes as a function of depth when described as an aggregate of plumes with a given spread and overlap. The discrepancy at low $z$ values can be attributed to the lack of sufficient root surface area and desorbing particles near the soil surface. We see that for all cases the total variance decreases in depth $z$ at the predicted rate, meaning that the mixing is dictated by the average distance between desorbers. That distance is a consequence of the random distribution of desorbing particles, set to match the surface area of the roots. We stress that the concentration variance is directly proportional to $ℓ_{ξ}$ , whereas the mean concentration is inversely proportional. This means that with increasing distance between desorbers we expect to see increasingly inhomogeneous concentration, meaning plumes within the substrate. Similarly, the mean concentration within the domain will decrease at larger desorbing particle distances, at any fixed depth $z$ . It is important to mention here that matching the surface areas was imposed by the authors. The choice is motivated by a few factors, the largest of which is that the adsorption and desorption rates are identical, thus this is the simplest method to balance these forces within the full domain. It is difficult to determine the adsorptivity rates of plants as it is a multi-step process and difficult to measure in isolation. We have therefore chosen the reactivity of biochar as per our earlier investigations (Pettersson et al., 2025), a physically relevant benchmark to begin such an analysis. This allows an analysis of the dynamics without significant interference generated from large differences in surface area/reactivity stemming from purely the root geometries. In practice, it may be that the desorbing surface area is much smaller than that of the adsorbing due to the requirements for minimal solute runoff, or that there exists a maximum limit wherein additional matter is no longer desirable. Then, given the choice of equalized surface area, a direct link can be seen between the root surfaces and the mixing of the plumes. The fibrous root has a larger surface area, meaning that the desorbing particles are more tightly spaced and act to homogenize the concentration at the outlet more effectively than the tap root, as seen in Figure 6(a) by the narrower breakthrough distribution curve than that of the tap root, which spans a larger range of values.

We have also varied $ℓ_{ξ}$ by doubling the number of desorbers in order to check the sensitivity to $ℓ_{ξ}$ . We have maintained Da $\approx 1$ by halving the reaction rate, to preserve the system dynamics for comparison. It is clear that there is almost no effect on the result with this modification. The cases with doubled desorber quantity are given by $ℓ_{ξ} = 3.3, 2.8$ in Figure 7, while the base cases are $ℓ_{ξ} = 4.2, 3.5$ , respectively.

We can say that by prescribing the surface area of randomly distributed desorbing particles we effectively assign an average vertical packing distance $ℓ_{ξ}$ . This in turn modifies the distribution of the breakthrough. Given that we have enforced equal surface area in our system, the choice of plant dictates the quantity of biochar added to the substrate if a more homogeneous breakthrough is desired. In a sense, limits are set on the distance between biochar required by the adsorbing root surface area. Anything higher will result in lower breakthrough but a less homogeneous one, and vice versa if too much biochar is added.

By increasing the number of desorbing particles we expect to have a more homogeneous mixing, with the concentration variance $σ_{c}^{2} / μ_{c}^{2}$ scaling with depth by $(ℓ_{ξ} / z^{2})$ , allowing us to predict the outlet concentration moments from those closer to the upstream surface. These measures can allow for tuning of the desorber placement and frequency such that vegetation can receive the nutrients they require, without causing excessive breakthrough of the solute.

Conclusions

This work proposes a methodology by which the interaction between different roots and the biochar desorbers within a soil substrate can be quantified. We have analyzed the adsorptive efficiency of two different root configurations under steady flow conditions corresponding to an extreme rainfall event by solving the advection-diffusion-reaction equation using the lattice Boltzmann method. An analysis of the fertilizer runoff is undertaken, with each root geometry varying the homogeneity of the breakthrough. It is shown that fibrous root structures act to homogenize the concentration breakthrough profiles to a greater degree than the tap root. An investigation of plume statistics and concentration fluctuations near the reactive particles shows that we can adequately characterize the evolution of the concentration profile as an overlap of plumes and a consequence of the root surface area. This relationship is quantified using a measure of the average distance between desorbing particles $ℓ_{ξ}$ which dictates the mixing strength within the porous medium. In essence, with the assumption of a similar reaction rate, one can manipulate the distribution and magnitude of the breakthrough by characterizing the surface area of a given plant type. This can be used to minimize unwanted contaminants from reaching the groundwater and lessen the costs of biochar amendment in green roof substrates by more accurately predicting the amount required.

If one applies this to a green roof design, there are several relevant lessons to be taken. Choice of plant type will have an impact on the breakthrough under heavy rainfall conditions. This is exacerbated if the green roof is thin, on the order of a few centimeters, such as extensive green roofs, which are typically around 5 cm thick and more common in practice due to their lower installation costs. If the designer chooses more fibrous root plants such as sedum or grasses they should expect to see more homogeneous breakthrough (a constant quantity of fertilizer being washed out) whereas if a taproot plant such as dandelion is used, the mean breakthrough will be lower but there may be events which have significantly higher peaks, potentially problematic given the lack of bespoke drainage systems attach to typical extensive green roof systems.

Alternatively. greater breakthrough homogeneity is achieved by adding more biochar, thus reducing the distance between each active element, however this will increase overall breakthrough magnitude and may even act to poison or damage the vegetation if too much fertilizer is adsorbed and subsequently released in a short period. Thus, there exists a maximum amount of biochar that should be added to any green roof. Our work included 10% biochar by volume, around the maximum distribution that remains safe for the plants, though this can be exceeded in certain circumstances, particularly if bespoke drainage systems are in place and/or the chosen plant type is accustomed to higher levels of fertilizer in the soil.

In general, it is more beneficial to have the biochar dispersed in closer proximity to the plant root structure, though of course the roots will grow to the nutrient source in time regardless. The plant adsorption will induce a greater solute gradient in proximity to the biochar desorbers, thus increasing the efficiency of the transfer irrespective of the flow velocity. An interesting design consideration can be taken from this; if one is concerned about root growth, one can incentivize them to grow in a particular direction by providing the nutrient source in a particular location, which can also act as a reservoir for any excess solute as well, preventing unnecessary contamination of the rainwater flowing through the roof.

This study focused on the impact of a singular root on the distribution of concentration within our domain. If one expands the domain or includes additional roots in close proximity, one would expect the total breakthrough to decrease, but not linearly, since the roots will compete for the desorbed “nutrients”. The authors looked into the effect of proximity between the adsorbing roots and desorbing particles, where it is known that plants will grow their roots to reduce the distance between the source of the nutrients and themselves. The concept of competition is certainly interesting, and to a degree, addressed in this work in the form of the fibrous root, where more root structures exist within close proximity to each other. As for neighboring root structures, it is more difficult. It has been pointed out that even within a singular root tendril there is heterogeneous root uptake, driven by entirely biological processes that are extremely difficult to predict. From a fundamentally mechanical/reactive perspective we expect the presence of additional root tendrils, be they from the same plant or a neighbor, will increase desorption by removing more of the solute from the system. If the roots are more spatially distributed, we should see a decrease in plume formation, as the solute will be removed faster and the roots systems will not grow in overlapping patterns, thus forcing the particles into closer proximity to the biochar, creating a more efficient transfer of the solute.

Footnotes

ORCID iD

Kaj Pettersson

Ethical considerations

No ethical approval was required in this study.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Swedish National Infrastructure for Computing (SNIC) and the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (FORMAS) [2019-01261, 2018-05973].

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.*

References

Abbas

Ali

Rizwan

(2018) A critical review of mechanisms involved in the adsorption of organic and inorganic contaminants through biochar. Arabian Journal of Geoscience 11(16). https://doi.org/10.1007/s12517-018-3790-1

Alim

Rahman

Tao

, et al. (2022) Green roof as an effective tool for sustainable urban development: An Australian perspective in relation to stormwater and building energy management. Journal of Cleaner Production 362: 132561. https://doi.org/10.1016/j.jclepro.2022.132561

Baqer

Chen

(2022) A review on reactive transport model and porosity evolution in the porous media. Environmental Science and Pollution Research 29(32): 47873–47901. https://doi.org/10.1007/s11356-022-20466-w

Barber

Silberbush

(1984) Plant Root Morphology and Nutrient Uptake, Chapter 4. John Wiley & Sons, Ltd., pp.65–87.

Ben-Noah

Hidalgo

Dentz

(2025) Pore-space partitioning in geological porous media using the curvature of the distance map. Transport in Porous Media 152(1). https://doi.org/10.1007/s11242-024-02142-4

Bonazzi

Dentz

de Barros

FPJ

(2023) Mixing in multidimensional porous media: A numerical study of the effects of source configuration and heterogeneity. Transport in Porous Media 146(1–2): 369–393. https://doi.org/10.1007/s11242-022-01822-3

Cropper

Comerford

(2005) Optimizing simulated fertilizer additions using a genetic algorithm with a nutrient uptake model. Ecological Modelling 185(2–4): 271–281. https://doi.org/10.1016/j.ecolmodel.2004.12.010

Cruz

Green

Watson

, et al. (2004) Functional aspects of root architecture and mycorrhizal inoculation with respect to nutrient uptake capacity. Mycorrhiza 14(3): 177–184. https://doi.org/10.1007/s00572-003-0254-5

Dore

(2005) Climate change and changes in global precipitation patterns: What do we know? Environment International 31(8): 1167–1181. https://doi.org/10.1016/j.envint.2005.03.004

10.

Duplat

Villermaux

(2008) Mixing by random stirring in confined mixtures. Journal of Fluid Mechanics 617: 51–86. https://doi.org/10.1017/S0022112008003789

11.

Elasbah

Selim

Mirdan

, et al. (2019) Modeling of fertilizer transport for various fertigation scenarios under drip irrigation. Water 11(5): 893. https://doi.org/10.3390/w11050893

12.

Gan

Garg

Huang

, et al. (2025) Experimental and numerical investigation on rainwater management of dual substrate layer green roofs using biochar-amended soil. Biomass Conversion and Biorefinery 15: 27387–27396. https://doi.org/10.1007/s13399-022-02754-0

13.

Gan

Garg

Wang

, et al. (2021) Influence of biochar amendment on stormwater management in green roofs: Experiment with numerical investigation. Acta Geophysica 69(6): 2417–2426. https://doi.org/10.1007/s11600-021-00685-4

14.

Girault

Chadil

Masi

, et al. (2024) Two-field and single-field representations of gas–solid reactive flow with surface reactions. International Journal of Multiphase Flow 175: 104796. https://doi.org/10.1016/j.ijmultiphaseflow.2024.104796

15.

Guo

Zheng

Shi

(2002) Discrete lattice effects on the forcing term in the lattice Boltzmann method. Physical Review E 65: 046308. https://doi.org/10.1103/physreve.65.046308

16.

Heyman

Lester

Le Borgne

(2021) Scalar signatures of chaotic mixing in porous media. Physical Review Letters 126: 034505. https://doi.org/10.1103/PhysRevLett.126.034505

17.

Hou

Pugazhendhi

Sindhu

, et al. (2022) An assessment of biochar as a potential amendment to enhance plant nutrient uptake. Environmental Research 214: 113909. https://doi.org/10.1016/j.envres.2022.1139092022.

18.

Huang

Yong

(2016) Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations. Journal of Computational Physics 310: 26–44. https://doi.org/10.1016/j.jcp.2016.01.008

19.

Huang

Yong

(2015) Boundary conditions of the lattice Boltzmann method for convection–diffusion equations. Journal of Computational Physics 300: 70–91. https://doi.org/10.1016/j.jcp.2015.07.045

20.

Irfan

Razali

KuShaari

, et al. (2018) A review of mathematical modeling and simulation of controlled-release fertilizers. Journal of Controlled Release 271: 45–54. https://doi.org/10.1016/j.jconrel.2017.12.017

21.

Joekar-Niasar

Hassanizadeh

(2012) Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Critical Reviews in Environmental Science and Technology 42(18): 1895–1976.

22.

Johansson

Fedje

Modin

, et al. (2024) Removal and release of microplastics and other environmental pollutants during the start-up of bioretention filters treating stormwater. Journal of Hazardous Materials 468: 133532. https://doi.org/10.1016/j.jhazmat.2024.133532

23.

Gong

Zheng

(2022) Effects of pore topology on immiscible fluid displacement: Pore-scale lattice Boltzmann modelling and experiments using transparent 3D printed models. International Journal of Multiphase Flow 152: 104085. https://doi.org/10.1016/j.ijmultiphaseflow.2022.104085

24.

Koch

Heck

Schröder

, et al. (2018) A new simulation framework for soil–root interaction, evaporation, root growth, and solute transport. Vadose Zone Journal 17(1): 1–21. https://doi.org/10.2136/vzj2017.12.0210

25.

Ladd

AJC

(1994) Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results. Journal of Fluid Mechanics 271: 311–339. https://doi.org/10.1017/S0022112094001783

26.

Le Borgne

Dentz

Villermaux

(2013) Stretching, coalescence, and mixing in porous media. Physical Review Letters 110: 204501. https://doi.org/10.1103/PhysRevLett.110. 204501

27.

Lei

Luo

(2021) Pore-scale simulation of miscible viscous fingering with dissolution reaction in porous media. Physics of Fluids 33(3): 034134. https://doi.org/10.1063/5.0045051

28.

Wang

, et al. (2023) Phase-field modeling of coupled reactive transport and pore structure evolution due to mineral dissolution in porous media. Journal of Hydrology 619: 129363. https://doi.org/10.1016/j.jhydrol.2023.129363

29.

Maggiolo

Modin

Kalagasidis

(2023) Transition from diffusion to advection controlled contaminant adsorption in saturated chemically heterogeneous porous subsurfaces. Physical Review Fluids 8: 024502. https://doi.org/10.1103/physrevfluids.8.024502

30.

Oliveira

TDS

Blunt

Bijeljic

(2020) Multispecies reactive transport in a microporous rock: Impact of flow heterogeneity and reversibility of reaction. Water Resources Research 56(12): e2020WR027317. https://doi.org/10.1029/2020wr027317

31.

Pettersson

Nordlander

Sasic Kalagasidis

, et al. (2025) Dynamics of contaminant flow through porous media containing random adsorbers. Transport in Porous Media 152. https://doi.org/10.1007/s11242-025-02150-y

32.

Puyguiraud

Perez

Hidalgo

, et al. (2020) Effective dispersion coefficients for the upscaling of pore-scale mixing and reaction. Advances in Water Resources 146: 103782. https://doi.org/10.1016/j.advwatres.2020.103782

33.

Qiu

Liu

Ling

, et al. (2022) Biochar for the removal of contaminants from soil and water: A review. Biochar 4(1). https://doi.org/10.1007/s42773-022-00146-1

34.

Rein

Legind

Trapp

(2011) New concepts for dynamic plant uptake models. SAR and QSAR in Environmental Research 22(1–2): 191–215. https://doi.org/10.1080/1062936X.2010.548829

35.

Rengel

(1993) Mechanistic simulation models of nutrient uptake: A review. Plant and Soil 152(2): 161–173. https://doi.org/10.1007/BF00029086

36.

Runions

Fuhrer

Lane

, et al. (2005) Modeling and visualization of leaf venation patterns. ACM Transactions on Graphics 24: 702–711. https://doi.org/10.1145/1073204.1073251

37.

Runions

Lane

Prusinkiewicz

(2007) Modeling Trees with a space colonization algorithm. In: Ebert

Merillou

(eds) Eurographics Workshop on Natural Phenomena. The Eurographics Association.

38.

Sangani

Acrivos

(1982) Slow flow through a periodic array of spheres. International Journal of Multiphase Flow 8(4): 343–360. https://doi.org/10.1016/0301-9322(82)90047-7

39.

Shafabakhsh

Le Borgne

Renard

, et al. (2024) Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media. Advances in Water Resources 192: 104791. https://doi.org/10.1016/j.advwatres.2024.104791

40.

Succi

(2001) The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond. Oxford University Press.

41.

(2024) An improved weighted topology optimization lattice Boltzmann model for porous structures of advection–diffusion chemical reaction systems. Chemical Engineering Journal 495: 153267. https://doi.org/10.1016/j.cej.2024.153267

42.

Twohig

Casali

Aydin

(2022) Can green roofs help with stormwater floods? A geospatial planning approach. Urban Forestry & Urban Greening 76: 127724. https://doi.org/10.1016/j.ufug.2022.127724

43.

Wei

Cai

, et al. (2023) Investigating the effect of biochar application on raindrop-driven soil erosion under laboratory rainfall experiments. Geoderma 430: 116291. https://doi.org/10.1016/j.geoderma.2022.116291

44.

WMO (2025) World meteorological organization’s world weather & climate extremes archive. Available at: https://wmo.asu.edu/content/world-greatest-sixty-minute-one-hour-rainfall (accessed 26 June 2025).

45.

Regenauer-Lieb

Tian

(2019) A hybrid immersed boundary-lattice Boltzmann/finite difference method for coupled dynamics of fluid flow, advection, diffusion and adsorption in fractured and porous media. Computational Geosciences 128: 70–78. https://doi.org/10.1016/j.cageo.2019.04.005

46.

Ying

Kookana

(2009) Reduced plant uptake of pesticides with biochar additions to soil. Chemosphere 76(5): 665–671. https://doi.org/10.1016/j.chemosphere.2009.04.001

47.

Zhang

Mao

Hilty

, et al. (2024) Volumetric lattice Boltzmann method for pore-scale mass diffusion-advection process in geopolymer porous structures. Journal of Rock Mechanics and Geotechnical Engineering 16(6): 2126–2136. https://doi.org/10.1016/j.jrmge.2024.03.006

48.

Zhao

Zhang

, et al. (2024) Dynamic pore network modeling of imbibition in real porous media with corner film flow. Langmuir 40(14): 7364–7374. https://doi.org/10.1021/acs.langmuir.3c03534

49.

Zolina

Simmer

Gulev

, et al. (2010) Changing structure of European precipitation: Longer wet periods leading to more abundant rainfalls. Geophysical Research Letters 37(6). https://doi.org/10.1029/2010gl042468