The effect of root exudates on rhizosphere water dynamics

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.

TSG, 0000-0003-3231-2159; TR, 0000-0001-8710-1063 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 plantsoil interactions. This will increase our ability to predict how contrasting root exudation patterns may influence crop uptake of water and nutrients. 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][10][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][14][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][17][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  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.

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]: Here, the notation· denotes a dimensional value, φ is the dimensionless fluid phase field, which takes the value φ = 1 in water and φ = 0 in air,ũ u u = φũ u u (w) + (1 − φ)ũ u u (a) is the combined velocity, whereũ u u (w) andũ u u (a) are the water and air velocities, respectively, andσ σ σ = (Ṽũ u u) + (Ṽũ u u) T is the viscous stress tensor.ρ(φ) is the phase-dependent density, which takes the value of the density of air when φ = 0 and takes the value of the density of water when φ = 1,η(φ) is the phasedependent viscosity,g is the acceleration due to gravity,t is time,ζ is the drag coefficient between the water and air, andp is the combined pressure that enforces incompressibility of both the water and air phases. The mobility is generally free to choose, up to some structural requirements for tensors, and we have usedM = φ 2 (1 − φ) 2 /ζ for consistency with the homogenization of two fluids literature [8]. We note that in the case of Daly & Roose [8], it was derived directly from a free energy based on the idea that the two fluids exert a drag on each other with drag coefficientζ . The capillary pressure,μ, is given bỹ  whereγ is the surface tension between air and water, α = 6 √ 2,λ is the air-water interface thickness, which is small compared with the geometry on which the model is applied, and 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 φ = 0 + O(λ) and φ = 1 + O(λ) everywhere except for a region of widthλ, where φ 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,ũ u u, and the capillary pressure,μ. Equation (2.1b) is Stokes equation with capillary pressure and gravitational effects. Equation (2.1c) ensures that the mass of both the water and air is conserved. The Cahn-Hilliard model has been used here as it has previously been homogenized allowing us to build on existing theory [8,[26][27][28]. Here we extend from the previous application to a fuel cell in [26,27] to consider fluid flow in soil. Using a phase field variable within the Cahn-Hilliard framework allows the model to be more general than other models such as Korteweg theory, which is restricted to using density as the order parameter [29].
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. u u u = 0,x ∈ ∂B. (2.3a) 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 bŷ n n n ·Ṽφ = |Ṽφ| cos(θ(c)),x ∈ ∂B, (2.3b) and a zero flux condition for the capillary pressure is given bŷ The initial saturationS is chosen and used to establish the initial condition for the phase,φ, i.e. φ(x, 0) =φ(S), whereS the fluid model are dependent on the concentration of root exudate. In keeping with the literature, we assume that the root exudate concentrationc affects the air-water contact angle at the surface of the porous material, θ (c), surface tension,γ (c), viscosity,η(φ, c) [13,14,30] and the fluid-fluid drag coefficientζ (c). The functions θ(c),γ (c) andη(φ, c) can be parametrized using experimental measurements [2]. The concentration-dependent functions are expected to depend nonlinearly on the concentration of root exudate. Fitting to the experimental data from [2] shows that this dependence is quadratic. The fluid-fluid drag coefficient is particularly difficult to parametrize, and it is assumed to be linearly dependent on the viscosityζ =ςη(φ, c), whereς is a constant with units m −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, η(c, φ), takes the value,η (a) when φ = 0, andη (w) (c) when φ = 1. Here,η (a) is the viscosity of air andη (w) (c) is the viscosity of water, which depends on the concentration of root exudate. The viscosity function is defined asη (2.5) 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,D = 0, and that in the air-filled pore space the exudate is immobile. As the interface is mobile, an additional term needs to be included to ensure the exudate is advected with the interface as it moves. This is achieved by adding an advection term to the diffusion model, which accounts for the movement of the phase, φ. 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 V =L yṼ , so the unit cell, Y = (0, 1) 3 , is defined with fluid part B and solid boundary ∂B [8]. The macroscopic length scale isL x andL y /L x = ε 1. We assume the non-dimensional surface tension to beγ =γ γ (c), whereγ is a baseline surface tension, e.g. the surface tension of water in air at 20 • C, 0.072 N m −1 , and γ (c) is the dimensionless exudate concentration-specific surface tension. We choose λ =λ/L y . The non-dimensional variables, denoted without·, are given by Vμ, y y y ∈ B, with boundary conditions u u u = 0,n n n · Vφ = |Vφ| cos(θ(c)), y y y ∈ ∂B and the transport equation where σ σ σ = Vu u u + (Vu u u) T , and D is assumed to be constant. Here, to simplify the analysis, we have neglected the influence of gravity on the air phase, i.e.ρ (a) /ρ (w) O(ε) [8]. The scalings used here have been chosen so that Υ ∼ O(1) so that the only small parameters are ε 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,L x andL y , can be treated independently; and secondly, that the underlying geometry is periodic on the pore scale.
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 Using Taylor series expansions, we also write, 1/εη 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 ε.
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 and with boundary conditions,n n n · V y φ 0 = |V y φ 0 | cos(θ(c 0 )), y y y ∈ ∂B, 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 φ 0 (x x x, y y y, 0) =φ(S(x x x), y y y), where We are interested in the behaviour of the fluids and exudates on a timescale that is much longer than τ −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 x x), i.e. μ 0 is constant in y y y, and equation (2.10b) is satisfied for any p 0 ∼ p 0 (x x x). Following Daly & Roose [8], we note that at steady state φ 0 is a function of x x x and y y y and can be written as is the modulated part of φ 0 with zero average, and the only x x x dependence in φ 0 comes through the saturation S(x x x). Also, at steady state, equations (2.10d,g) have the solution As the value of μ 0 is determined by equations (2.10a,c,e), it depends on the initial condition,φ, the surface tension, γ (c 0 ) and the contact angle θ (c 0 ). We assume that, γ i.e. the surface tension and contact angle do not depend on higher order terms in 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τ −1 = τ −1 γ (C) andμ 0 = μ 0 /γ (C), to obtain a set of equations which depend only on the geometry and contact angle, (2.12c) n n n · Vφ 0 = |V y φ 0 | cos(θ(C)), y y y ∈ ∂B The result is that, for fixed values of S and θ , we can calculateμ 0 by solving equations (2.12). By repeating this for a range of different saturation values while keeping the contact angle fixed, we obtain where F[S, θ(C)] is the SWCC for contact angle θ (C). Hence, at steady state in τ −1 we write (2.14) , is the numerical analogue to fitting the well-known van-Genuchten or Brooks and Corey models to an experimentally measured SWCC [6,7]. In this case, however, the SWCC will be different for different contact angles. Hence, if θ (C) is known then the SWCC can be calculated as a function of both saturation and concentration. In reality, it will not be possible to solve equations (2.12) for all values of S andC. However, the curve can be numerically generated by solving for a fixed set of S andC and using interpolation to obtain the function at all saturation and concentration values. 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 x x, y y y) = η y (y y y)η x (x x x) it is impossible to separate the x x x and y y y dependence. Hence, we approximate the leading-order viscosity as where φ 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 O(ε 0 ) from equations (2.8). Our aim is to find an approximation for the effect of large-scale pressure and concentration gradients on the local phase, velocity, pressure and concentration. The steps required to derive these equations are included in the electronic supplementary material S1. We are interested in the behaviour of the porous medium on the long timescale, i.e. we assume that the equilibration of the air-water interface occurs on a timescale τ −1 that is much quicker than the one associated with bulk fluid movement. In addition, we show that the O( 0 ) fluxes balance and the resulting equations are independent of τ 0 , electronic supplementary material S1. Hence, we can write and where the boundary conditions are u u u 0 = 0,n n n · V x φ 0 +n n n · V y φ 1 = |V y φ 1 + V x φ 0 | cos(θ(C)), y y y ∈ ∂B, with boundary conditionn n n · (φ 0 V y c 1 − c 1 V y φ 0 ) = −n n n · φ 2 0 V xC , y y y ∈ ∂B. (2.16h) We note that the advection term from equation (2.8g) is not present as transport is dominated by the diffusion term when λ → 0, as shown in electronic supplementary material S1. Equations (2.16) describe the phase, velocity and pressure of the air and water phases and the exudate concentration of the water on the short space scale. We require that these equations are satisfied , (2.17a) and 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, 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 , 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 .
(2.20f ) 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.
(2.21b) 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 θ (C), provide a complete description of how the pore scale geometry and physical processes are dependent on large-scale variations in capillary pressure, combined pressure, acceleration due to gravity and the concentration of root exudates. In the next section, we will derive equations for terms of order ε 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 O(ε 1 ). Then, by applying a solvability condition and taking a volume average over the domain 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

22e)
and we recall the SWCC is given by The derivation of these equations involves a large number of steps and the details are included in the electronic supplementary material S1. Equation   [8]. Alternatively, if we wanted to consider diffusion of solutes, which did not directly influence the properties of water, then setting γ (C) = η (w) (C) = 1 and θ (C) = θ (0) would provide a partially coupled set of equations describing the movement of air and water, from which the diffusion of solutes, which do not bind to the soil particle surfaces, could be calculated.

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Ū U U = 0, corresponding to the case in which the air viscosity is much smaller than the water viscosity. In this case, we also assume 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 γ (C) = 1 and note that, in the caseŪ U U = 0, the only effect of the surface tension changing is a variation in the timescales which appear in equation (2.22). We also make the simplification θ(C) = 0, corresponding to a fully wetted soil. Mathematically, the contact angle features in all the cell problems either directly, or indirectly through the function φ 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 2a,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 2b. The cell problem is calculated on one-eighth of the unit cell and it is of size 0.5 3 (dimensionless) (figure 2c). 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 2d). 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 3a). The parameters F(S, 0) and, hence, the effective parameters K[S, θ (C)] and 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 (  The initial exudate concentration also exhibits dynamics, due to the saturation perturbation (figure 4b). 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.
where 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 and are shown in table 2. These values are small (less than or equal to 7%) and similar to the variability between different experimental methods [33]. For the small perturbation all the models exhibit the same behaviour, resulting in a smaller maximum errors (less than or equal to 1.1%) and the results can be seen in the electronic supplementary material, figure S1.2.

(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. 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 V x S is effectively infinite at the Haines jump. In order to analyse the behaviour of the model at this point, future research could use matched asymptotics to attempt to capture the behaviour around a Haines' jump. This would allow the homogenized solution to be defined by appropriate jump conditions. 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.

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].