## Abstract

Most water and nutrients essential for plant growth travel across a thin zone of soil at the interface between roots and soil, termed the rhizosphere. Chemicals exuded by plant roots can alter the fluid properties, such as viscosity, of the water phase, potentially with impacts on plant productivity and stress tolerance. In this paper, we study the effects of plant exudates on the macroscale properties of water movement in soil. Our starting point is a microscale description of two fluid flow and exudate diffusion in a periodic geometry composed from a regular repetition of a unit cell. Using multiscale homogenization theory, we derive a coupled set of equations that describe the movement of air and water, and the diffusion of plant exudates on the macroscale. These equations are parametrized by a set of cell problems that capture the flow behaviour. The mathematical steps are validated by comparing the resulting homogenized equations to the original pore scale equations, and we show that the difference between the two models is ≲7% for eight cells. The resulting equations provide a computationally efficient method to study plant–soil interactions. This will increase our ability to predict how contrasting root exudation patterns may influence crop uptake of water and nutrients.

## 1. Introduction

Jethro Tull's [1] 1762 observation that ‘roots are but as guts inverted… that spew out what is superfluous’ recognized the capacity of plants to exude chemicals from their roots to capture nutrients from the soil. More recent research has observed that these surface active chemicals alter fluid properties at the root–soil interface considerably [2]. This has the potential to affect storage and transport of water and nutrients. Plant breeding may be able to manipulate the chemistry and quantity of these exudates to improve resource capture and stress tolerance from droughts and floods, potentially addressing food security by improving yield. Certainly between species of plants, root exudate properties vary considerably. Understanding of root exudates and root–soil interactions can lead to advances in agricultural techniques to improve food production, particularly in extreme conditions [3]. Mathematical modelling provides one route through which the complexities of water and nutrient movement in soils can be understood [4]. Hence, developing new mathematical models, which use the vast computational resources that are now available, will lead to a significant improvement in root–soil interaction models. This in turn will further improve, and potentially optimize, crop yield.

Richards' equation is widely used to model the movement of water through partially saturated porous media, including soil, at large scales [5]. Traditionally, Richards' equation is derived by combining conservation of mass with Darcy's Law and parametrized by equilibrium measurements of the soil water characteristic curve (SWCC) and hydraulic conductivity [6,7]. More recently, Daly & Roose [8] used the Cahn–Hilliard–Stokes equations, which have been used to model two fluid flow in porous media [9–11], to show that these hydraulic properties of a partially saturated soil can be evaluated from the underlying geometry.

The equations derived by Daly & Roose [8] are appropriate for modelling bulk soil. However, they might not be directly applicable to the region of soil close to the roots over which the plants have influence, known as the rhizosphere [12]. The rhizosphere can have different structural, chemical, biological and hydraulic properties to the bulk soil [13–15]. This can be partially due to the presence of root exudates. Root exudates mix with soil pore water, creating a diffusion gradient away from the surface of the root. When exudates mix with soil pore water they can decrease the surface tension and soil–water content, increase viscosity and affect the contact angle between soil particles and pore water [2,16–18]. These impacts, however, vary between species of plants and it may be possible to breed future crops that produce exudates with desired physical impacts on soil pore water.

The impact of the altered hydraulic properties on the movement of water at large scales has been investigated, but is not well understood [13,15]. Raynaud [19] used Fickian diffusion to model exudates from a root in a simple cylindrical geometry. The water content was assumed constant and the water movement was controlled by the rate of uptake by the root, a sink term was used to model exudate decay and a source term was present at the root surface to model exudation from the root. The gradient of the adsorption isotherm, decay rate of the exudate, and soil water content primarily determined the time for the concentration to reach steady state. The effects of exudates on the macroscale have been considered in conditions where the viscosity dominates over surface tension in regions of high root exudate content [15,20]. They found that an increase in hydraulic connectivity of the rhizosphere due to the formation of liquid bridges. Daly *et al.* [13] used the model derived in Daly & Roose [8] to study the effect of increased contact angle, surface tension and viscosity. However, they did not explicitly model exudate diffusion or how this would affect the derivation of Richards' equation.

The work presented here is motivated by the effect of root exudates on soil, however, the theory can also be applied to areas such as geological waste disposal, oil production or oil-spill clean-up problems. Numerical methods have been used to investigate two fluid flow with mass transfer on the pore scale for applications in chemical engineering, such as determining the rate of CO_{2} capture [21,22]. Yang *et al.* [21] and Haroun *et al.* [22] implemented the Navier–Stokes equations, using the one fluid formulation with a characteristic function to define the interface, and are coupled to the mass transfer equation through the local velocity. In these studies, the solute concentration does not affect the behaviour of the fluid flow and the solute is able to diffuse across the interface between the two fluids, which have different diffusion coefficients [21,23,24]. Davidson & Rudman [23] considered a solute within a spherical drop of one fluid containing a solute within a second fluid. They validated the numerical calculations by comparison with the analytical solutions and considered mass transfer of a solute from a drop rising in a fluid column. Haroun *et al.* [24] examined the effect of a periodic corrugated geometry on mass transfer and found that recirculation zones, which held up the movement of the liquid, affected the mass transfer because it changed the shape of the fluid–fluid interface. Yang *et al.* [21] created a model of a microscale segmented flow microreactor in OpenFOAM, which shows the gas transfer between a gas and liquid phase, and could be used to optimize this type of system.

In this paper, we derive macroscale models for water movement in soil that take account of changes to fluid properties due to the presence of root exudates and the underlying pore scale geometry. To do this, we have extended the derivation of Daly & Roose [8] by developing a pore scale description of exudate diffusion, which we have coupled to a two fluid model for water movement. By including coupling terms to link the fluid properties to the exudate diffusion we were able to capture the effect of exudates on hydraulic properties. We have applied homogenization theory [25] to upscale the model from the pore scale to the macroscale, e.g. pot or field scale, and have obtained a set of coupled equations for water movement and the diffusion of exudates. The upscaling procedure used to develop the macroscale model has been validated against the underlying pore scale equations using an idealized geometry. The upscaled equations agree with the underlying pore scale equations within less than 7% error.

## 2. Derivation of equations

In this section, we describe the derivation of the macroscale coupled flow and diffusion equations. Our aim is to start with a set of equations on the pore scale and to use these to derive a set of macroscale equations. We will start with the Cahn–Hilliard two fluid model and couple this to a phase-dependent diffusion equation, which describes the movement of root exudates through the pore water.

### (a) The pore scale model

We consider a macroscale soil domain, *Ω*. On the pore scale, this is composed of repeating periodic units. The periodic units contain a fluid domain, *B*, and the soil particle surface, ∂*B*, as illustrated in figure 1. We start with the dimensional Cahn–Hilliard–Stokes equations, which we write following the notation used in Daly & Roose [8]:
*a*
*b*
*c*Here, the notation *ϕ* is the dimensionless fluid phase field, which takes the value *ϕ* = 1 in water and *ϕ* = 0 in air, *ϕ* = 0 and takes the value of the density of water when *ϕ* = 1, *f*(*ϕ*) = *ϕ*^{2}(1 − *ϕ*)^{2} is the energy of the two fluid system. *f*(*ϕ*) is chosen to have minima at *ϕ* = 0 and *ϕ* = 1. Hence, equation (2.2) has the solution *ϕ* varies rapidly.

Equations (2.1) and (2.2) are the Cahn–Hilliard–Stokes equations, where equation (2.1a) describes the movement of the air and water phases driven by the velocity of the combined velocity,

In order to be able to homogenize the equations, it is required that the component of velocity normal to the soil surface is zero, i.e. a zero penetration condition must be used. This can be combined with a no-slip, or generalized Navier slip condition. Here, we demonstrate the method using a no-slip condition on the soil boundary, ∂*B*, i.e.
*a*A flux condition that defines the behaviour between the water–air interface and the soil particle surface using the soil–liquid contact angle *θ*, as in Daly *et al.* [13], is given by
*b*and a zero flux condition for the capillary pressure is given by
*c*The initial saturation *B*∥ is the volume of the pore scale domain. The initial velocity is defined as * u* = 0.

Assuming the functions ^{−2}. This implies that the greater the viscosity of the water phase, the greater the drag between the water and air phases.

The viscosities of the water and air are combined into one function using *ϕ*. The viscosity, *ϕ* = 0, and *ϕ* = 1. Here,

In order to couple the Cahn–Hilliard–Stokes equations with the diffusion of root exudates through water, we introduce a phase-dependent diffusion equation. We assume that the exudate can diffuse in water with a diffusion constant, *ϕ*. From a mathematical point of view, this term ensures that the concentration in the water phase is conserved and decays with the phase field at the air–water interface. Hence, the phase-dependent concentration equation is
*a*and
*b*Equations (2.6) are similar to the equation used by Haroun *et al.* [24,31] without the expression of Henry's Law which allows the chemical species to diffuse across the interface.

### (b) Dimensionless equations

Equations (2.1)–(2.6) describe the fluid flow in soil coupled with root exudate diffusion. Equations (2.1)–(2.6) are non-dimensionalized using *Y* = (0, 1)^{3}, is defined with fluid part *B* and solid boundary ∂*B* [8]. The macroscopic length scale is ^{−1}, and *γ*(*c*) is the dimensionless exudate concentration-specific surface tension. We choose *a*with
*b*where *M* = *ϕ*^{2}(1 − *ϕ*)^{2}. We define the dimensionless drag/interface parameter *a*
*b*
*c*
*d*with boundary conditions
*e*and
*f*and the transport equation
*g*with boundary condition
*h*

where *D* is assumed to be constant. Here, to simplify the analysis, we have neglected the influence of gravity on the air phase, i.e. *ε* and λ. The scaling results in a unit change in *μ* driving a fluid velocity of order *ε*^{−1} that corresponds to the fastest time scale, which is defined in the next section, and implies that the capillary forces dominate as we have assumed. Therefore, the velocity and gravity contributions first appear at order 1, corresponding to the intermediate timescale.

### (c) Homogenization

The full set of equations (2.8) is fourth order, stiff and nonlinear. Hence, it is time consuming and computationally expensive to solve them numerically on real soil geometries. To overcome this issue, we derive a set of equations that describe how the pore scale dynamics affect the macroscale behaviour using the method of multiple scale asymptotic homogenization, an averaging procedure for periodic structures [25]. Homogenization requires two assumptions: firstly, that the macro and micro length scales,

Homogenization is based on a linear expansion of the dimensionless equations in terms of *ε*. This leads to a cascade of numerical problems, in which the details of the pore scale geometry are captured by solving representative cell problems on one period of the domain [25]. The results from the cell problems then parametrize averaged macroscale equations that approximate the solution of the full set of equations. Importantly, the homogenized equations are solved on an averaged geometry, i.e. they only depend on the pore scale geometry through the averaged parameters. This makes the macroscale equations much more efficient to solve than the original full set of equations.

In the homogenization presented here, there are some differences to the standard procedure. Firstly, in addition to the macro and micro length scales, we also consider three different time scales. The fastest time scale, *τ*_{−1}, is associated with the leading-order pore scale dynamics, i.e. the equilibration of the air–water interface. The second time scale, *τ*_{0}, is the intermediate timescale associated with fluid flow driven by flux imbalances on the pore scale. Finally, the slowest time scale, *τ*_{1}, is associated with the macroscale behaviour, i.e. the slow variation in saturation due to macroscale pressure gradients. In addition, equations (2.8) are nonlinear and therefore the accuracy of the final macroscale approximation will depend on how well the equations can be approximated by linear expansions.

The homogenization procedure involves a large number of mathematical steps, which are somewhat analogous to the steps taken by Daly & Roose [8]. Hence, the full procedure has been detailed in the electronic supplementary material S1 and we have included here only the main steps and key results. We seek a perturbation expansion solution to equations (2.8), i.e. we use *η*_{1} = *c*_{1}(*δη*/*δc*)|_{c=c0}, *M*_{0} = *ϕ*^{2}_{0}(1 − *ϕ*_{0})^{2}, *M*_{1} = *ϕ*_{1}(*δM*/*δϕ*)|_{ϕ=ϕ0} and *δ*/*δv* is the functional derivative with respect to a function *v*. We substitute these expansions into equations (2.8) and solve the problems for increasing orders in *ε*.

#### (i) O ( ε − 1 )

The first step in the homogenization procedure is to collect the dominant terms from equations (2.8) using expansions (2.9). In this case, the largest terms are of size *ε*^{−1}. By collecting order *ε*^{−1} terms, we find
*a*
*b*
*c*
*d*with boundary conditions,
*e*
*f*
*g*

where *p*_{0}, *μ*_{0}, *ϕ*_{0} and *c*_{0} are periodic with period 1. The initial condition for the phase is defined, as per equation (2.4), for a chosen macroscale saturation, *S*, as *τ*_{−1}, i.e. longer than the time it takes for the air–water capillary interfaces to equilibrate on the pore scale. Hence, we are only interested in the steady-state behaviour of equations (2.10). However, as equation (2.10c) is nonlinear, the solution to (2.10a) is dependent on the initial condition. Therefore, to find the steady-state solutions of *p*_{0}, *μ*_{0}, *ϕ*_{0} and *c*_{0} from equations (2.10), it is necessary to include the time derivative in equation (2.10a) and numerically integrate to the steady state.

While we have to solve parts of equations (2.10) numerically, we are able to determine certain properties of the steady-state solution by analysing these equations. These properties will enable us to progress through the homogenization procedure without knowing the precise solution to equations (2.10). The numerical solutions can then be calculated to parametrize the resulting equations. We note that, at steady state, equation (2.10a) is satisfied for any *μ*_{0}∼*μ*_{0}(* x*), i.e.

*μ*

_{0}is constant in

*, and equation (2.10b) is satisfied for any*

**y***p*

_{0}∼

*p*

_{0}(

*). Following Daly & Roose [8], we note that at steady state*

**x***ϕ*

_{0}is a function of

*and*

**x***and can be written as*

**y***ϕ*

_{0}=

*S*(

*,*

**x***τ*

_{0},

*τ*

_{1}) +

*ϕ*

^{(m)}

_{0}(

*S*(

*,*

**x***τ*

_{0},

*τ*

_{1}),

*,*

**y***τ*

_{−1}), where

*ϕ*

^{(m)}

_{0}(

*S*(

*,*

**x***τ*

_{0},

*τ*

_{1}),

*,*

**y***τ*

_{−1}) is the modulated part of

*ϕ*

_{0}with zero average, and the only

*dependence in*

**x***ϕ*

_{0}comes through the saturation

*S*(

*). Also, at steady state, equations (2.10d,*

**x***g*) have the solution

*c*

_{0}.

As the value of *μ*_{0} is determined by equations (2.10a,*c*,*e*), it depends on the initial condition, *γ*(*c*_{0}) and the contact angle *θ*(*c*_{0}). We assume that, *c*. Rather, they are only dependent on the macroscale concentration of exudate and not the position of the air–water interface. As a result, we can scale the surface tension out of equations (2.10), using the scalings *a*
*b*
*c*
*d*
*e*

The result is that, for fixed values of *S* and *θ*, we can calculate *τ*_{−1} we write
*S* and *S* and

At this point, we return to the expression for the viscosity, equation (2.7b). In order to proceed, we need to address two issues with this equation. Firstly, since the Cahn–Hilliard equation allows the phase to be of order λ outside the defined range of (0, 1), inputting *ϕ*_{0} directly into the viscosity can result in negative viscosities, which are not physical. Secondly, as equation (2.7b) is not in the form *η*(* x*,

*) =*

**y***η*

^{y}(

*)*

**y***η*

^{x}(

*) it is impossible to separate the*

**x***and*

**x***dependence. Hence, we approximate the leading-order viscosity as*

**y***ϕ*

_{0}has been scaled using the maximum and minimum values to force the value to be strictly in the interval (0,1). The effect of having the macroscale viscosity outside the phase terms is that both water and air are dependent on the macroscale concentration, rather than just the water. As

*η*

^{(a)}≪

*η*

^{(w)}, the effect of this approximation should be small. We will quantify the effect of this approximation in §3.

#### (ii) O ( ε 0 )

We now collect terms of *τ*_{−1} that is much quicker than the one associated with bulk fluid movement. In addition, we show that the *τ*_{0}, electronic supplementary material S1. Hence, we can write
*a*
*b*
*c*
*d*where the boundary conditions are
*e*and
*f*and the exudate transport equation is
*g*with boundary condition
*h*

We note that the advection term from equation (2.8g) is not present as transport is dominated by the diffusion term when *a*
*b*
*c*
*d*
*e*

We substitute the solutions in a separable form, equations (2.17), into equations (2.16) and gather terms dependent on *μ*_{0}, the macroscale capillary pressure, to get,
*a*
*b*
*c*
*d*where *e*and
*f*

Equations (2.18) determine the fluid velocity driven by a large-scale variation in capillary pressure. Physically, this corresponds to the difference in pressure between the two phases. Hence, in the limit λ → 0, only the water phase is directly driven by the capillary pressure. The air phase is not directly driven by the capillary pressure, but can be set in motion by the water velocity at the air–water boundary.

Next, we gather terms dependent on the macroscale combined pressure, *p*_{0},
*a*
*b*
*c*
*d*where *e*and
*f*

Equations (2.19) determine the fluid velocity due to a unit pressure gradient, in this case both the air and water phases are driven by the combined pressure, *p*_{0}.

Next, we gather terms dependent on gravity, *g*,
*a*
*b*
*c*
*d*where *e*and
*f*

Equations (2.20) determine the fluid velocity due to gravity. As we chose to neglect the effect of gravity on the much less dense air phase, only the water phase is directly driven by gravity. Any induced movement of the air phase comes from the effect of water movement at the air–water interface.

Finally, we gather terms dependent on *c*_{0},
*a*with boundary condition
*b*

Equations (2.21) determine the local geometry-dependent diffusion offered by the soil pore structure and the position of the air–water interface. Physically, this will be combined with the unimpeded diffusion coefficient to calculate the effective diffusion coefficient in the water phase as a function of saturation.

Cell problems (2.18)–(2.21), for known *ϕ*_{0} and *ε*^{1}, which will relate these pore scale behaviours to the macroscale flow.

#### (iii) O ( ε 1 )

The final step in the homogenization procedure is to expand equations (2.8) to *B*, we find a series of macroscale equations that describe the averaged movement of air, water and exudate through the soil. Finally, we take the limit λ → 0 to obtain
*a*
*b*
*c*where
*d*and
*e*and we recall the SWCC is given by
*f*The derivation of these equations involves a large number of steps and the details are included in the electronic supplementary material S1. Equation (2.22a) is equivalent to the macroscale equation from Daly & Roose [8] for *a*
*b*
*c*

Equations (2.23a–*c*) describe the velocity of the water phase averaged over the pore scale domain driven by capillary pressure, combined pressure and gravity, respectively. Equation (2.22b) ensures the conservation of mass for the saturation equation (2.22a), where the concentration-dependent viscosity and surface tension are present, and are related to the pore scale behaviour by,
*a*
*b*
*c*

Equations (2.24a–*c*) describe the velocity of both the air and water phases averaged over the pore scale domain driven by capillary pressure, combined pressure and gravity, respectively. Finally, equation (2.22c) describes the movement of exudate on the macroscale, driven by diffusion and fluid movement, and is related to the pore scale behaviour by,

This is an advancement from the work of Daly & Roose [8] as we have combined the equations for fluid flow with the equations for exudate diffusion. If we neglect diffusion and root exudates then equations (2.22) reduce to those derived by Daly & Roose [8]. Alternatively, if we wanted to consider diffusion of solutes, which did not directly influence the properties of water, then setting

## 3. Validation of homogenization procedure

The homogenization procedure involved a large number of mathematical steps and assumptions. In order to test the accuracy of these assumptions and to transparently demonstrate the application of the method, we show that the macroscale model, equations (2.22) presented above, provides an accurate approximation of the pore scale equations, (2.8). We do this using an idealized geometry for two cases; firstly, we compare the homogenized equations to the pore scale equations with the original viscosity, equation (2.7b). Secondly, we compare the homogenized model to the original equations with the assumed viscosity, equation (2.15). Our overall aim is to test the accuracy of the mathematical steps, not to compare the predictions of this model with experimental results, which is beyond the scope of this paper.

### (a) Test geometry

We consider a soil column in which the saturation is perturbed by subtracting a linear function from the initial phase configuration. This drives a re-equilibration of both the fluid phase and the local concentration of exudate. All equations are solved numerically using Comsol Multiphysics (v5.2a). In order to simplify the analysis, we consider the case *p*_{0} is constant and the effects of gravity can be neglected.

For the purposes of validation, we focus on the effect of viscosity on water movement and exudate diffusion. Hence, we choose *ϕ*_{0}. The effect of changing the contact angle in this way has been investigated by Daly *et al.* [13]. Hence, here we neglect these changes and focus only on viscosity.

We consider the geometry shown in figure 2*a*,*b*. This idealized soil physical structure is composed of soil particles of two different sizes. The soil geometry is built from repetition of a unit cell. At the centre of the cell, we have placed a spherical soil particle of radius 0.49 [dimensionless] and at the corners of the cell we have placed a spherical soil particle of radius 0.35 [dimensionless]. These two sizes have been chosen to force the geometry to exhibit wetting–drying hysteresis, i.e. the SWCC will exhibit different values depending on whether the water content is increasing or decreasing. This hysteresis requires pores of different sizes, something which does not occur for a single soil particle size for a perfectly packed structure. The full geometry is formed from eight repetitions of the unit cell (which we label unit 1 to unit 8). However, to reduce the computational load we have used symmetry to quarter the domain, see figure 2*b*. The cell problem is calculated on one-eighth of the unit cell and it is of size 0.5^{3} (dimensionless) (figure 2*c*). It was meshed with tetrahedral elements, maximum size 0.0335 and minimum size 0.002, and one boundary layer on the surfaces of the soil particles. This mesh was sufficient to resolve the interface with width 0.02 (dimensionless) using quadratic elements. A mesh refinement study was conducted to ensure this was a suitable mesh. The same mesh was used for the full model. The homogenized model used a one-dimensional mesh with 96 evenly spaced elements. All the models were solved using the MUMPS solver with Newton's method, with a constant damping of 1 so that it has quadratic convergence, within the fully coupled solver in Comsol Multiphysics. The solution was declared converged when an absolute tolerance of 10^{−4} was reached.

The homogenized equations (2.22), are applied to a unit cuboid (figure 2*d*). In order to ensure that the homogenized equations (2.22) converge, we require that there are enough repetitions of the unit cell such that *L*_{y}/*L*_{x} = *ε*≪1. As the full equations are computationally expensive, we also require that we have a small enough domain that the resulting equations are computationally tractable. By gradually increasing the number of repetitions of the unit cell, we found that eight repetitions was sufficient for the homogenized equations (2.22) to converge. However, good convergence between the models can also be found for as few as four repetitions.

### (b) Numerical results

Before we solved the homogenized equations, the cell problem, equations (2.12), (2.18) and (2.21), were evaluated for a range of dimensionless capillary pressures between −7 and −4 [dimensionless] with λ = 2 × 10^{−3} [dimensionless] and *Υ* = 1 [dimensionless]. As we have neglected the effects of combined pressure and gravity, it is not necessary to solve cell problems (2.19) and (2.20). We solved equations (2.12), (2.18) and (2.21) for both increasing and decreasing capillary pressures in order to evaluate the wetting and drying curves (figure 3*a*). The parameters *F*(*S*, 0) and, hence, the effective parameters *D*_{eff}(*S*) exhibit wetting drying hysteresis and Haines' jumps [32], where the topology of the air–water interfaces changes resulting in a large change in saturation for a small change in capillary pressure. These parameters feed into equations (2.22) and are used to calculate the solution to the homogenized equations.

To establish the initial configuration for the phase of the full model, equation (2.8), the dimensionless capillary pressure was set to −4.25 and the model was run to steady state, resulting in an average saturation of 0.6 over each unit in the column. This was outside the hysteresis loop of the SWCC at the wet end (figure 3*a*). The saturation was perturbed in the full model by subtracting a linear function from the initial configuration of the phase, i.e. *A* is the magnitude of the perturbation and *h* is the vertical position of the column. For *A* > 0, this resulted in a greater saturation at the top of the column and lower saturation at the bottom. The equivalent condition was applied to the homogenized model using the expression *S*(*t* = 0) = 0.6 − *A*(1 − (*h*/8)). Two perturbations were used, a small perturbation with *A* = 0.1 [dimensionless] and a large perturbation with *A* = 0.25 [dimensionless] that allows for comparison between the models where Haines' jumps and hysteresis are present. The initial concentration of the root exudate was set to 0.2 [dimensionless] for the small perturbation and 0.5 [dimensionless] for the large perturbation. The expression for the concentration-dependent viscosity was chosen to be
*q* is 1059. This value comes from fitting the function to the concentration versus viscosity relation measured for barley root exudate from Naveed *et al.* [2] (M Naveed 2017, unpublished data, see electronic supplementary material, figure S1.1).

We compare the results from the full model and the homogenized model as a function of time. All equations were integrated to steady state and the dynamics can be seen in figure 4. For the large perturbation, evidence of a Haines' jump can be seen in figure 4*a* as not all the unit saturations converge to the same steady-state value. For the full model with assumed viscosity, equation (2.15), units 1 and 2 converge to a saturation of approximately 0.4, unit 3 converges to approximately 0.46, whereas the other units converge to approximately 0.51. The homogenized model exhibits the same behaviour as the full model with the assumed viscosity, and captures the effect of the Haines' jump.

The initial exudate concentration also exhibits dynamics, due to the saturation perturbation (figure 4*b*). It can be seen that unit 8, which has the greatest initial saturation, has the smallest exudate concentration due to the dilution effect. The opposite can be observed for unit 1. The homogenized model exhibits behaviour at time 10^{1} [dimensionless], where the exudate concentration of each unit overshoots the steady-state value of 0.63 for the exudate concentration, before all the units converge to 0.63 concentration. For the full model in which the assumed viscosity (2.15) is used, the concentration tends to the same value as the homogenized model, with values of 0.62–0.64 at steady state. The full model in which the original viscosity (2.7b) is used, displays different behaviour to the homogenized model for low saturation, units 1, 2 and 3. At time 10^{−2} [dimensionless], units 1 and 2 in the original viscosity model have converged to the same saturation, 0.38. They also converge to the same saturation at steady state as the homogenized and assumed viscosity model. Unit 3, however, converges to 0.49 at steady state, compared to the saturation of 0.46 reached by the homogenized model. This is due to a topological difference in the final phase configuration as can be seen by comparing figure 5*b*,*c*.

In figure 5, the result of the homogenized model is compared to sections from both of the full models. It can be seen that where the homogenized model has a smaller saturation, caused by the hysteresis and Haines' jumps present in the cell problem results, there is a different phase configuration in the full model. This configuration is different in the original viscosity and assumed viscosity models. The relative differences between the three models are calculated by
*N* is the total number of time steps taken by the solver for the homogenized model, *x*^{h}_{n} is the result of the homogenized model at time step *n* and *x*^{f}_{n} is the result of the full model (either original or assumed viscosity model) interpolated onto the time steps taken by the homogenized model solver. The relative differences are shown in table 1. The maximum errors between individual points of the models are calculated by

### (c) Discussion

The homogenized equations provide a good estimate for the saturation and concentration when compared with the full equations. This is particularly true when considering the assumed viscosity model (2.15). Overall, it is not surprising that the homogenized equations provide a better approximation to the assumed viscosity, rather than the original viscosity equations, as the assumed viscosity is the one used in the homogenization scheme. The key difference between these two viscosities is the location at which the Haines' jump is seen to occur.

A significant advantage to the homogenized model is that it is less computationally expensive than the full model. The cell problems for the homogenized model took 20 min and 3 GB ram for each capillary pressure evaluated, in this case 24. This resulted in a total computation time of 8 h on a laptop. The resulting homogenized model took just 2.5 min and 5.8 GB of ram on a standard desktop computer. The full model took up to 26 h and 35 GB ram using a segregated solver on the Iridis 4 supercomputer at the University of Southampton.

The homogenized model captures the fluid behaviour near a Haines jump with maximum error less than 7%. We only consider eight repetitions of the unit cell, which corresponds to *ε* = *L*_{y}/*L*_{x} = 0.125. Formally, we would expect an error of size *ε* to be present in the model. With this in mind, we see that the model provides a much more accurate description of the saturation and concentration than would be expected, i.e. the error is less than 12%. It is interesting to note that the homogenized model is able to capture the hysteresis and Haines jumps despite these violating one of the key assumptions used in the derivation. Formally, we require that *S* varies only on the long spatial scale *L*_{x}. However, by definition

In §2c, it is assumed that the surface tension and wetting angle are functions of the spatially averaged concentration in order to be able to separate the two length scales of the problem. This assumption is not valid in the case where significant concentrations of the root exudate are absorbed on the air–water interface or soil particle surfaces, and therefore some important physical processes may be ignored. However, in the present model, the concentration is dominated by the macroscale part of the concentration, *C*_{0}(*x*), and therefore we would expect that the effects of these missing physics would be small compared with the effect of the spatially averaged concentration. Furthermore, in soil systems, adsorption of organic compounds to mineral surfaces is complex and the physical characteristics of root exudates have only recently been explored [2,14]. Greater experimental characterization would be required before accounting for particle adsorption and air–water interface impacts from surfactants in root exudates, but the model could be adapted to account for these processes in the future.

## 4. Conclusion

In this paper, the Cahn–Hilliard–Stokes equations were combined with a phase-dependent diffusion equation for the root exudates and this full set of equations was homogenized. The homogenized equations were shown to reduce to a exudate concentration-dependent Richards' equation, (2.22a), and a saturation-dependent exudate diffusion equation, (2.22c). This homogenization procedure relied on two key assumptions. Firstly, the fluid–fluid drag coefficient is linearly dependent on the viscosity; and secondly, the viscosity of the air and water are sufficiently different that the error induced by assuming the air viscosity was dependent on exudate concentration was negligible.

In §3, the method was validated by comparing the homogenized, equations (2.22), and assumed viscosity model, equations (2.8) with (2.15), to the original viscosity model, equations (2.8) with (2.7b). The relative errors over the whole simulation for the homogenized and assumed viscosity solutions compared to the original viscosity solution were less than 1%. The maximum errors for a particular time point were less than 7%. This shows that the homogenized model is an appropriate approximation to the full set of equations.

The homogenized model is a much more attractive option for calculating the local macroscale flow and diffusion properties of soil than the pore scale description. This is particularly evident if the pore scale geometry is taken from three-dimensional images, a technique that is becoming more frequently used [13,34]. The model derived in this paper captures the effects of root exudates and could therefore be applied to rhizosphere soil. The homogenized model can also be used to investigate the impact of altered hydraulic properties on the movement of water at large scales and be used to improve our understanding of these effects. Detailed models of this kind have the distinct advantage that they can be used to relate specific porescale soil geometries to the observed macroscopic soil properties. The simulations presented here conserved the root exudate within the domain during wetting and drying. In a natural soil system, root exudates may leach and their physical properties may alter due to microbial decomposition over time [2]. Moreover, adsorption of the exudates onto soil minerals will also affect behaviour, potentially altering properties such as the contact angle and surface tension [2,14]. Future modelling could incorporate these impacts to assess the time variability and longevity of root exudate impacts on flow and retention in soil. This would be particularly useful in exploring how properties change with age along a growing root. Hence, this method could be used to gain greater understanding of soil hydraulic properties in and around the rhizosphere and even lead to the possibility of theoretically simulating the quantity of root exudates required to improve the chances of a plant thriving in particular soil environments. Ultimately this reveals the theoretical potential for plant breeders to manipulate root exudation in the development of crops with more effective root systems. The theory can be further extended to other applications such as geological waste disposal, oil production or oil-spill clean-up problems.

## Data accessibility

Data supporting this study are available on request from the University of Southampton repository at https://doi.org/10.5258/SOTON/D0609 [35].

## Authors' contributions

L.J.C. and K.R.D. co-developed the Comsol implementation, conducted the numerical simulations and drafted the manuscript. N.K. was consulted on the physical properties of soil and interpretation of the results. P.D.H., T.S.G. and T.R. designed the study and were consulted on the physical properties and of soil and interpretation of the results. All authors contributed to writing the manuscript and gave final approval for publication.

## Competing interests

The authors have no conflict of interest to declare.

## Funding

L.J.C. and N.K. are funded by BBSRC SARISA BB/L025620/1, L.J.C. is also funded by EPSRC EP/P020887/1. K.R.D. is funded by ERC 646809DIMR. P.D.H. and T.S.G. are funded by BBSRC BB/J00868/1. The James Hutton Institute receives funding from the Scottish Government. T.R. is funded by BBSRC SARISA BB/L025620/1, EPSRC EP/M020355/1, ERC 646809DIMR, BBSRC SARIC BB/P004180/1 and NERC NE/L00237/1.

## Acknowledgements

The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, and the *μ*-vis centre at the University of Southampton for the provision of necessary software and high performance computer support in the completion of the work. The authors would like to thank the ‘Rooty Team’ at the University of Southampton, Glyn Bengough at School of Science and Engineering, University of Dundee and The James Hutton Institute, Invergowrie, Dundee, and Muhammad Naveed at School of Biological Sciences, University of Aberdeen, for helpful discussions related to this work.

## Footnotes

Electronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.4203062.

- Received March 2, 2018.
- Accepted August 3, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.