Multiple Scale Homogenisation of Nutrient Movement and Crop Growth in Partially Saturated Soil

In this paper, we use multiple scale homogenisation to derive a set of averaged macroscale equations that describe the movement of nutrients in partially saturated soil that contains growing potato tubers. The soil is modelled as a poroelastic material, which is deformed by the growth of the tubers, where the growth of each tuber is dependent on the uptake of nutrients via a sink term within the soil representing root nutrient uptake. Special attention is paid to the reduction in void space, resulting change in local water content and the impact on nutrient diffusion within the soil as the tubers increase in size. To validate the multiple scale homogenisation procedure, we compare the system of homogenised equations to the original set of equations and find that the solutions between the two models differ by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim 2 \%$$\end{document}≲2%. However, we find that the computation time between the two sets of equations differs by several orders of magnitude. This is due to the combined effects of the complex three-dimensional geometry and the implementation of a moving boundary condition to capture tuber growth.


Introduction
Application of solutes such as fertilisers and pesticides is important in modern agricultural practices (Godfray et al. 2010). However, more efficient solute application is needed in order to mitigate growing costs of fertilisers and environmental pollution, i.e. fertiliser and pesticide buffering, leaching and run off (Godfray et al. 2010 understanding water and solute movement in soil is vital for determining sustainable crop production for long-term food security (Comas et al. 2013). To aid with this goal, mathematical modelling of soil systems has been studied increasingly in recent years (Vereecken et al. 2016), since this offers one method to investigate plant-soil interactions while reducing time and resources compared to standard experimental practices. Combining mathematical modelling with traditional experiments allows us to efficiently improve our understanding of plant-soil interactions Daly et al. 2017). This can lead to further improvement of agricultural techniques for greater crop yield while minimising waste of resources.
Mathematical modelling of soil systems covers a wide range of spatial scales, including pore, plant and field scales (Darrah et al. 2006;Hopmans et al. 2002). As such, when studying transport of water and solutes in soil, complex geometries are often required to capture the intrinsic details contained in the microscopic structure of the scale that is considered. This typically requires vast amounts of computation time and resources (Daly and Roose 2018). Hence, it is often favourable to construct an averaged macroscopic geometry so that the macroscale transport properties can be attained directly from the microscale information (Bruna and Chapman 2015). One technique that is frequently used to obtain macroscale movement of fluids and solutes in soil or other porous media is multiple scale homogenisation (Hornung 2012). This mathematical technique is a method of devising a system of averaged macroscopic equations that are parameterised by associated cell problems, which are derived from the inherent microscopic structure of the domain (Pavliotis and Stuart 2008).
Multiple scale homogenisation has been successfully used in a wide range of porous media and soil applications, including modelling saturated fluid flow (Keller 1980), two-phase fluid flow (Daly and Roose 2015), wave propagation in poroelastic materials (Sharma 2007) and single-phase fluid flow in double porosity systems (Arbogast et al. 1990). One application that has been increasingly studied in recent years is homogenisation of moving interfaces for first-and second-order partial differential equations (Cardaliaguet et al. 2009;Lions and Souganidis 2005). Although there has been extensive research on the mathematical theory for the homogenisation of moving interfaces, few applications have been explored.
In this study, we demonstrate the utility of homogenisation by modelling the growth of potato tubers in soil, in which the growth is dependent on the quantity of nutrients the plant is able to draw up from the soil. We model the soil as a poroelastic material, such that any growth from a single crop will influence the water content adjacent to the plant and therefore the movement of nutrients in the vicinity. We use a combination of poroelastic theory and the diffusion equation in porous media to model the movement of nutrients in a deforming soil environment. We develop a series of approximate equations to describe nutrient movement, growth in tuber size and global nutrient uptake in soil.
There has been previous research which studied the effect of diffusion with spatially varying objects in porous media (Bruna and Chapman 2015), in which Rayleigh's multipole method was used to determine a spatially dependent effective diffusion coefficient based on the size of the sphere within the microscopic periodic geometry (Rayleigh 1892). Here, we extend this idea to model both spatially and temporally Fig. 1 Schematic of a dimensional poroelastic domain, where˜ is the total domain, Soil is the deformable poroelastic soil domain,˜ p j are the potato tubers,˜ Soil j are the poroelastic soil subdomains adjacent to each tuber and˜ j are the boundaries between˜ p j and˜ Soil . In addition, l x is the macroscale and l y is the microscale (Color figure online) varying objects in poroelastic media, which are coupled to the diffusion of the species within the material itself.
For simplicity, we choose to model the tubers as spherical objects in soil; however, this can be extended to any 3D geometry, including, but not limited to, ovoids, capsules and cylinders. To validate the homogenisation procedure, we compare the solution of the homogenised equations against the full system for a series of case studies. This shows the homogenised equations successfully capture the growth of each tuber and the change in nutrient diffusion from the reduction of volume within the domain.

Three-Phase Poroelastic Soils
Let˜ ⊂ R 3 be an open bounded subset representing a soil system ( Fig. 1) that contains N potato tubers. We define˜ =˜ Soil ∪ N j=1˜ p j , where˜ Soil is the deformable poroelastic soil domain that is composed of water, air and solid components, and˜ p j are the j = 1, . . . , N potato tubers each with a boundary˜ j .
To describe the deformable poroelastic soil domain˜ Soil , we impose a system of equations that describe a three-phase poroelastic domain. To derive the system of equations, we use the conservation laws for mass and momentum. The conservation of mass equations for the three phases of air, water and soil solid is where φ a is the volumetric air content, φ w is the volumetric water content, φ s is the volumetric soil solid content,ṽ a is the air velocity,ṽ w is the water velocity,ṽ s is the velocity of the soil solid component andp w is the soil water pore pressure. Water uptake in our simulations is assumed to be dominated by transport through symplastic pathways, thus passively taken up by pressure gradients in the root xylem (Roose and Fowler 2004b). The ratio between the cortex and the xylem hydraulic conductivities along with the root surface area density is characterised by λ c , and the root xylem pressure is expressed as p r . Roots are assumed to be uniformly distributed throughout the soil domain. We note that we neglect the impact that tuber growth has on the root system. The expression − λ c (p w − p r ) represents water uptake by plant roots. Furthermore, Darcy's law for the relative phase velocity of air and water is written as wherep a is the soil air pore pressure, κ a and κ w are the air and water permeabilities, respectively, and μ a and μ w are the viscosities of air and water, respectively. The air and water pressuresp a andp w , and the air and water volume fractions φ a and φ w are related via the van Genuchten saturation expression (Van Genuchten 1980) where S w = φ w /(φ w + φ a ) is the relative water saturation, p c is the characteristic suction pressure and m is the van Genuchten parameter. The conservation of momentum equation is (Wang 2000) where G is the stress tensor,ũ s is the displacement of the solid soil matrix, G is the shear modulus of the soil solid, ν is the Poisson ratio, S a = φ a /(φ w + φ a ) is the relative air saturation and T is the identity tensor. The displacementũ s is related toṽ s by the relationshipṽ The system of Eqs.

Diffusion of Nutrients in Soil
Solutes such as nutrients typically exist in one of two states in soil, either sorbed to the soil solid surfaces or dissolved in the pore water (Roose et al. 2001). We state that the nutrient concentration in the sorbed state follows a reversible linear binding reaction such that, wherec s is the sorbed nutrient concentration and d s is the net transfer rate to the sorbed phase from the pore water phase. From the conservation of mass, the rate of change of the nutrient concentration in the pore water phase is wherec is the nutrient concentration in pore water, d l is the net transfer rate to the pore water phase from the sorbed phase, D is the diffusion coefficient and g is the nutrient uptake rate by plant roots. Adding (11) and (12) yields We assume there is a direct jump between the nutrients in the two states with no intermediate phase, such that d s + d l = 0. Furthermore, we define d s where k a is the adsorption rate of the nutrient in solution and k d is the desorption rate. We assume k d is sufficiently large such that d s /k d = ∂tc s /k d 1 and k a ∼ k d , theñ where b = k a /k d is the buffer power of the nutrient (Nye and Tinker 1977;Barber 1995;Roose and Fowler 2004a). This leads to the governing equation for nutrient movement in terms ofc only, i.e.

Boundary Conditions
Here, we define a series of boundary conditions on the interfaces˜ j , i.e. between the deformable poroelastic soil domain˜ Soil and the potato tubers˜ p j . To describe nutrient interaction on˜ j , we impose a zero flux condition, as the potato tubers take up nutrients through their rooting systems and not through the tuber surfaces: wheren is the unit normal vector pointing out of the geometry. Furthermore, on˜ j we assume the soil solid is displaced normally to the direction of the growing tuber, hence where ξ j is the displacement of the j th tuber given by, where r * is the initial radius of the tubers andr j is the radius of the jth tuber, which is related to the total amount of nutrients taken up by the roots. The growth of each tuber is expressed as whereṼ j is the tuber volume, α is the ratio between the rate of growth and uptake and˜ Soil j is the volume of soil adjacent to each potato tuber j (see Fig. 1). Here, we model the early-stage development of potato tubers (diameter 5-7 cm); hence, we approximate the tubers shape to be spherical. Therefore, Eq. (20) can be written in terms of the radiusr j only, i.e.
We state the water and air components of˜ Soil do not penetrate the tubers˜ p j ; thus, we require the Darcy velocities normal to the interface to be zero, i.e.
Finally, on˜ j we assume the air and water velocities are equal to the growth of the tubers, hence 2n ⊗n − T ·ṽ a =n∂tr j ,x ∈˜ j .

Non-dimensionalisation
To simplify the model and understand the magnitude of influence of each parameter, we non-dimensionalise the system of equations described above. We are interested in the macroscopic properties of the system of equations while retaining the influence of the microscopic structure. Hence, we identify two different length scales, the 'microscopic' length scale l y associated with the inner domain tuber geometry, and the macroscopic length scale l x associated with the full domain transport of water and nutrients. Under these scales, l y /l x = ε 1. We choose to non-dimensionalise using the scalingx where c max is the maximum concentration of the nutrient applied to˜ Soil and i = {w, a}. In (26) we use the macroscopic length scale l x as the spatial scaling to observe macroscale properties, the diffusion timescale l 2 x D for the time scaling and the shear modulus G for the pressure scaling. Shown in Fig. 2 are the non-dimensionalised macroscopic domain and microscopic domain . It follows that the air, water and solid phase continuity equations become: with the constitutive poroelastic mechanical law represented as: where the force balances and relative movement of the air and water in the mixture domain are represented as: The relationship between water and air is linked based on the Van Genuchten water retention relationship: The nutrients in the system follow the convection-diffusion equation: where we have no flux through the tuber surface: The soil solid phase displacement is equal to the increase in the tuber radius: Neither water nor air is assumed to flow through the tuber surface: The water and air velocities normal to the tuber surface also follow the rate of the tuber growth: Finally, the tuber growth rate is based on the rate of nutrient uptake out of the system: Here, the system was non-dimensionalised as follows:

Parameter Estimation
Here, we estimate the parameters contained in Eqs. (27)-(41) in order to determine the magnitude of influence each parameter has on the system of equations. Since this model is motivated by the growth of potato tubers in soil, we assess the parameter values for silt soils as potatoes are frequently grown in this soil type (Shock et al. 1998).
Potato plants are typically grown in ridge and furrow type systems and are contained in the plough layer of soil, which is the top 30 cm of soil (Lesczynski and Tanner 1976). Hence, we choose the macroscopic length scale to be l x ≈ 0.3 m. Similarly, we assume that the tubers have an inter-tuber distance that is substantially less than the total length of the plough layer. We choose an inter-tuber distance of approximately l y ≈ 0.05 m, resulting in the ratio of the two length scales to be ε ≈ 0.1. We also assume an initial tuber radius of r * = O(0.05) m < l y .
One of the key nutrients responsible for plant growth and development is nitrogen (Nye and Tinker 1977). We choose to model this nutrient since plant growth is closely linked to abundance of nitrogen in soil. Nitrogen has a diffusion coefficient in soil water of D ≈ 2.5 × 10 −10 m 2 s −1 (Barber 1995). Furthermore, for the potato plant Solanum tuberosum L, the uptake rate of the nutrient nitrogen is g ≈ 1×10 −9 s −1 (Sattelmacher et al. 1990;Asfary et al. 1983). This was found to be in nitrogen concentrations in soil of c max ≈ 10 −1 kg m −3 (Asfary et al. 1983).
In early-stage growth of Solanum tuberosum L plants, the tuber radius growth rate is approximately 1 × 10 −9 m s −1 (Xu et al. 1998). If we assume that the quantity of nitrogen that is taken up by the plant is proportional to the growth of the tuber, then we can estimate the ratio between the rate of growth and the uptake, i.e. α ≈ 1 × 10 1 kg −1 m −1 (Sattelmacher et al. 1990;Asfary et al. 1983).
Using the values above, we find that the parameters κ a and κ w in equations (31) and (32) are κ a = O(10 9 ) and κ w = O(10 7 ). This is significantly larger than the other terms in the equations. Hence, we rewrite Eqs. (31) and (32) so that, which have the solutions p a = constant and p w = constant, i.e. the consolidation of the soil is substantially faster than the diffusion of solutes. Since p w = constant, we find that the sink term in Eq.
(2) representing root uptake is constant, i.e. λ c (p w − p r ) = F, where F is the water uptake rate by plant roots. The uptake rate of water from Solanum tuberosum L roots is F ≈ 1 × 10 −8 s −1 (Parker et al. 1989). From Equation (33), the solutions p a = constant and p w = constant result in S w = constant, and since S w + S a = 1, this leads to S a = constant. Although S w is constant, φ w will still change as a function of the changing domain geometry. Substituting 44 into 32 renders the domain water content to become dependent on the solid phase displacements: Thus, we reduce the system of Eqs. (27)-(41) to where F = Fl 2 x /D. Using the values discussed above, we find that the parameters contained in (46)-(51) have the approximate values For the remainder of this study, Eqs. (46)-(51) will be referred to as the 'full set' of equations to describe solute movement and tuber growth.

Homogenisation
In this section, we use multiple scale homogenisation to develop a set of averaged macroscale equations that describe the movement of nutrients and tuber growth in soil. From Equation (46), we observe that φ w is affected by two mechanisms: firstly by soil compression due to the growth of the tuber, i.e. ε∇·(φ w v s ), and secondly by root water uptake, i.e. F. From the non-dimensionalisation, we observe that the maximum displacement is bounded such that u s F. This leads to the results ∂ t u s F. If we consider a scenario when the tubers are not taking up water, Eq. 45 suggests that the primary change in water content is based on ε∇·(φ w ∂ t u s ). Since ε∇·(φ w ∂ t u s ) ∂ t u s and ∂ t u s F, then it follows that ε∇ · (φ w ∂ t u s ) F. Therefore, we find that the root water uptake term dominates the change in water content. Hence, for the homogenisation procedure, we neglect the term regarding soil compression, and the system of equations we homogenise reduces to To validate this assumption, we compare the full set of Eqs. (46)-(51) to the homogenised system of equations derived from (53)-(57) in the following section. We highlight that the horizontal boundaries for the full model preserve the periodicity presented here.
We observe there are two different length scales present in the geometry˜ , the macroscale l x and the microscale l y . Any change of O(1) on the length scale l x will result in a O(ε) change on the length scale l y . We can formalise this by assuming that the dependent variables φ w , c and r are functions of a small scale y and a large-scale x.
We denote the unit cell representing the microscale domain y ∈ ≡ [−1/2, 1/2] 3 . Using the two length scales and chain rule, the gradient operator is written as Furthermore, we expand φ w , c and r such that, The first step of the homogenisation procedure is to determine the most dominant terms in the system of Eqs. (53)-(57). To do this, we substitute Eqs. (59)-(61) into (53)-(57), collecting the largest terms O(ε −2 ). This results in the system of equations Theorem Equations (62)-(64) have the solution c 0 = c 0 (x, t), i.e. c 0 has large-scale dependence only.
From the theorem above, we observe that c 0 has large-scale dependence only and is independent of the small scale y; however, we receive no other information regarding the solution of c 0 .
To proceed with the homogenisation methodology, we collect the next most dominant terms in the system of equations. This is achieved by collecting terms O(ε −1 ) and using the results ∇ y c 0 = 0, i.e.
To ensure (70)-(72) form a well-posed problem, i.e. the equations have a solution that agrees with the boundary conditions, we check the solvability of the system. We can show the system is well-posed by applying the divergence theorem to equation (70) and use the boundary condition (71) such that, Next we choose to rescale c 1 such that, wherec 1 (x) is the large-scale component of c 1 (x, y). Substituting (74) into (70)-(72) yields the cell problem for χ k whereê k is the unit vector. We note that the tubers grow in the soil domain; hence, the cell problem solution χ k is dependent on the radius of the tuber. Since the cell problem is a representation of the impedance of nutrient movement due to tuber obstruction, and as the tuber grows, the impact on nutrient transport will change; therefore, we have the relationship χ k = χ k (r ), i.e. the cell problem solution is dependent on the radius of the tuber.
The last step of the homogenisation procedure is to collect terms O(ε 0 ), i.e.
To check (78)-(82) provide a well-posed problem, we check the solvability of the system of equations. To do this, we apply the divergence theorem to (79) and using boundary condition (80) yields We define to be the volume integral of the cell problem, which is dependent on the radius of the tuber. It follows that (84) can be written as This results in the approximate equations for φ w 0 , c 0 and r 0 where for k = (1, . . . , 3).
Here, the averaged terms || Soil (r 0 )|| and D e (r 0 ) are parameterised from the cell problem (75)

Validation of the Homogenisation Procedure
We validate the mathematical steps used in the homogenisation procedure by comparing the homogenised set of Eqs. (87)-(90) to the full set of Eqs. (46)-(51). We consider multiple comparisons by varying parameters for the buffer power b, root uptake rate F and initial volumetric water content φ w | t=0 to examine the accuracy of the averaging procedure.
We generate two geometries, one for the full set of Eqs. (46)-(51) containing potato tubers and the other uniform geometry for the homogenised Eqs. (87)-(90). We choose the domain length of each geometry to be composed of eight periodic cells. Due to the homogenisation procedure, the approximate Eqs. (87)-(90) do not require any tubers as the influence of the microscale geometry is contained in the parameterised terms || Soil (r 0 )|| and D e (r 0 ). Shown in Fig. 3 are the geometries used to validate the homogenisation procedure.
Lastly, we numerically illustrate that the full solution tends towards the homogenised solution as ε → 0. Since ε = l y l x , it suffices to show that the solutions become closer as l x → ∞, as this implies that ε → 0. We begin with by setting l x = 0.3, where the domain of the full solution can only consist of 3 tubers. We incrementally increase the domain size up to l x = 0.8, where the domain consists of 8 tubers. We take the percent difference between the homogenised solution for the concentration profile as: where c (ε) is the solute concentration profile in the full solution as a function of a given , c (ε) 0 is the solute concentration profile based on the homogenised solution, c (ε) − c (ε) 0 ∞ is the largest difference between the two solutions for all time points in the whole domain for a fixed ε value, and c (ε) ∞ is the supremum concentration value for all time in the full domain for the full solution for a fixed ε value.
To solve the systems of equations, we use the finite element package COMSOL Multiphysics ® 5.3 (www.comsol.com). We run our full model with a mesh consisting of 21729 tetrahedral elements and 1405 for the homogenised model. Simulations were run using the MUMPS (Multifrontal Massively Parallel Sparse) direct solver for a fully coupled physical system. In this section, we describe the implementation of each set of equations and show a comparison between them.

Full Equations
Implementation of the full set of Eqs. (46)-(51) requires the implementation of a complex moving boundary problem. This accounts for the uptake of nutrients by each tuber p j , the subsequent growth of p j and the reduction in volumetric water content φ w . The geometry we impose the full set of equations on can be seen in Fig. 3a. However, we require two versions of this geometry, an undeformed geometry that is constant in time, and a deforming geometry that is dependent on tuber growth, since different components of the system (46)-(51) are solved on either an undeformed or deforming frame of reference. There are three main components that are required to be implemented in order to solve (46)-(51); these are the poroelastic equations, the compaction and deformation of soil, and the nutrient movement equations.
To implement the poroelastic equations (46)- (47) and (49) for the local displacement u s and reduction in φ w is straightforward, since these equations are solved on the undeformed geometry regardless of tuber size. Using this solution at each time step, we can prescribe a deformation (for the deforming geometry) within the soil domain to correspond with the increase in tuber size.
The nutrient equations (48) and (50)-(51) are solved on the deforming geometry to correspond with the growth of the tubers. However, these equations use the poroelastic solution from the undeformed geometry. Hence, we implement a reference frame change such that poroelastic solution can be mapped from the undeformed geometry to the deformed geometry. This allows us to solve the nutrient equations on the deformed geometry corresponding with the prescribed tuber deformation.
Since the nutrient equations are solved on a deforming geometry, we are required to ensure that c is conserved. This is achieved by making two alterations to (48) and (50). Firstly, we note Reynolds transport theorem where dV and dA are volume and surface elements, respectively, ω is the velocity of the surface element,n is the normal vector pointing out of the geometry, F is any function of x and t and θ(t) is the domain. Reynolds transport theorem states that the change in nutrient concentration in a domain is equal to the change in concentration within the domain plus the rate at which nutrient is entering the domain. Applying equation (92) to the full set of equations we have leads to d dt where ω mesh is the velocity of the boundaries p j . This requires us to adapt equation (50) so that,n Equation (94) then satisfies the conservation law for moving boundaries. Secondly, as p j grows and S is deformed, this causes an advective movement effect on c within S . This can be interpreted as the boundaries of the tubers and j physically pushing the nutrients. Hence, we are required to add a conservative advection term to Eq. (48) accounting for the individual elements within the mesh moving, i.e.
This modified system of equations can then be successfully implemented to model coupled nutrient movement and poroelastic deformation from growing tubers.

Homogenised Equations
The geometry used to simulate the homogenised set of equations can be seen in Fig. 3a. However, to solve the set of homogenised Eqs. (87)-(90), we are required to solve a series of cell problems, i.e. Eqs. (75)-(77), to calculate the terms || Soil (r 0 )|| and D e (r 0 ) that parameterise Eqs. (88) and (89). Since the geometric properties of the domain are contained in || Soil (r 0 )|| and D e (r 0 ), we solve the cell problem for a series of different tuber radii to correspond with different levels of growth/displacement from the original tuber size. Using the results from the cell problems, we can construct interpolated functions to describe || Soil (r 0 )|| and D e (r 0 ) as functions of the homogenised radius r 0 .

Results
To validate the homogenisation procedure, we compare the homogenised Eqs. (87)-(90) against the original set of Eqs. (46)-(51). We choose to run a series of case studies by varying the parameters b, F and φ w | t=0 . For the buffer power b, we choose the values b ∈ {0.5, 5} since this covers a range of buffer powers for the nutrients nitrogen, boron, magnesium, zinc and molybdenum (Barber 1995). From the nondimensionalisation and parameter estimation, we observe the value for root water uptake is F = O(1). However, to test the homogenisation procedure, we select the values F ∈ {0.1, 10} for low and high levels of water uptake, respectively. Finally, for the initial water content φ w | t=0 we assign the values φ w | t=0 ∈ {0.4, 0.6} as these are approximate upper and lower bounds for silty soils (Das 2013).
In each of the simulations we impose a Dirichlet condition of c = c 0 = 1 on the top of each of the geometries shown in Fig. 3a. Additionally, we choose the initial nondimensionalised tuber radius to be r * = 0.025 and choose the remaining parameters to be g = α = 1. We also impose a stop condition on each of the simulations so that when the non-dimensionalised volume of a tuber has doubled in magnitude, the simulation is terminated. Finally, in order to construct interpolated functions to describe || Soil (r 0 )|| and D e (r 0 ) in Eqs. (88) and (89), we solve a series of 6 cell problems with varying sphere radii.
Shown in Fig. 4 are the nutrient profiles for c and c 0 down the length of the geometries shown in Fig. 3a. We observe for all buffer powers, root uptake values and initial porosities, that the homogenised nutrient profile for c 0 is qualitatively identical to the full nutrient concentration c. We find there to be a maximum error of 2% between the solutions across all scenarios.
Additionally, shown in Fig. 5 are the individual tuber radii r j for the full set of equations against the effect radius r 0 from the homogenised equations. Similar to the results from Fig. 4, we find that the effective radius r 0 successfully captures the growth of each tuber within the full domain shown in Fig. 3a. We find there to be a maximum error of 2% between the actual and effective tuber radius.
To highlight the accuracy of the homogenised set of equations, shown in Fig. 6 are detailed results for the simulation using the parameters F = 0.1, b = 0.5 and φ w | t=0 = 0.4. From Fig. 6a we observe that the effective radius r 0 is able to mimic the growth of the tubers in the full geometry. Figure 6b illustrates that as we have a system with more tubers, the full model converges to the homogenised model. For the 8 tuber scenario, the growing tubers can be seen in Fig. 6c, in which the tubers at the top of the full equation domain at the time point t = end have grown substantially larger than those at the base of the domain. Furthermore, we find that the solute concentration profiles exhibit identical trains between the full and homogenised domains.
As a final confirmation of our model accuracy, we increased our mesh density from 21,729 tetrahedral elements to 49,218 elements in order to insure that our results are sufficiently accurate. Percent difference between the refined vs course solutions was on the order of 0.1%. This grants us confidence that our mesh is sufficiently resolved given the current problem. It is worth noting that our current simulations do not consider any automatic mesh refinement, as the deformations modelled do not result in any large aspect ratios. Future considerations could include more general geometries or greater deformations (Dehghani et al. 2018). This would be more significant when considering large deformations, which are more common in soil materials (Yu 2000).
From Figs. 4 and 5, we observe that the homogenised equations successfully capture the nutrient movement and tuber growth in soil. However, the computation time between the two systems of equations differs by several orders of magnitude. We find that the full set of equations in three dimensions requires ≈ 5 min (300 s) to solve one simulation for eight periodic cells. Conversely, solving the homogenised equations requires ≈ 10 s to solve an analogous 3D simulation. Furthermore, the homogenised set of equations can be reduced to a 1D problem which will achieve the same results as the 3D problem due to the homogenisation procedure. We find that the computation time to solve the 1D problem is 1 s, which is substantially faster than the full set of equations. However, a set of 3D cell problems is required to parameterise the homogenised set of equations for the terms || Soil || and D e . In this case study, we chose to conduct six cell problems for varying sphere radii. Each of the cell problems requires ≈ 10 s to solve. However, these cell problems are only required to be solved once for each set of parameters. Hence, we find that the homogenised sets of equations can reduce the computation time substantially while retaining a high level of accuracy. Furthermore, we can highlight the influence that the tubers' radii have on the effective homogenised diffusion coefficient (Fig. 7). Under more dramatic growth scenarios where tubers increase their radii by a factor of 5, the effective diffusion in the system could be reduced by as much as 30%.

Discussion
In this study, we developed a physical model for potato tuber growth that couples water and nutrient uptake with mechanical growth of potatoes in soil. The explicit consideration of the potato growth in the soil domain creates a physical impedance to nutrient transport through the soil. The geometry and the surface sinks due to the presence of potatoes impede the effective transport of nutrients through the soil domain (Fig. 6). If impedance to diffusion caused by the potato tubers was not considered, we would incur an error between 175 and 300% in the effective diffusion of the solute (where the effective diffusion is D eff = φ D φ+b in the case with no potatoes and (φ(r 0 )+b)|| Soil(r 0 ) || when impedance to diffusion caused by growing tubers is homogenised). Furthermore, an error of up to 62.5% in effective diffusion could occur if the tubers were modelled as a sphere with constant half of the final time-dependent radius. These errors in effective diffusion would greatly impact the predicted solute leaching or plant solute uptake of models, ignoring geometric impedance to diffusion.  2015), our biological agents grow proportional to the rate of nutrient uptake. However, our potato tubers also take up water, which impacts the advective fluxes associated with the nutrient transport in the unsaturated soil domain. As the focus of our system is to obtain a geometrically simplified model through our homogenisation procedure, the final equations that arise are convection-diffusion equations. By choosing a different re-scaling approach, it may be possible to obtain a similar Darcy type expression as demonstrated by O'Dea et al. (2015) and Penta et al. (2014); however, this was not within the scope of this study.
Previous studies have coupled fluid and solid mechanical systems to infer not only the impact that a solid inclusion would have on the fluid flow, but also the mechanical deformations that fluid flow would induce on the solid inclusion (Royer et al. 2019;Chen et al. 2019). Authors have found that the homogenised system parameters are impacted by the distribution of inclusions in the domain. While our modelling scheme does not explicitly account for the mechanical response of the inclusions to externally applied stresses, the distribution of our potato tubers impacts the flow and transport coefficients in a similar manner as demonstrated in previous studies (Royer et al. 2019). It is worth noting that the growth behaviour of our modelled tubers implies a compensation for external stresses. Plant roots are known to respond to mechanical stresses by increasing their radii and reducing their length (Abdalla et al. 1969). This behaviour does not readily lend itself to a simple coupling between mechanical stresses and growth responses, and future work should be conducted to better quantify these contrasting effects.
The full system of equations in this paper required the implementation of a complex moving boundary problem. This required the use of multiple domains to solve different components of the equations, and subsequent mappings of solutions across domains.
Not only does this system require considerable computational power to solve, the time required to correctly implement this system is substantial. This is due to ensuring conservation of mass and consistent mappings of solutions across domains. Using mathematical homogenisation, many of the more cumbersome modelling aspects were simplified into an effective media, where the tuber surfaces are treated as domain sinks, and the tuber geometries are accounted for in the diffusivity term shown in Figs. 4 and 5. Applying a similar method to a root system would facilitate a more rigorous quantification of bulk scale rhizosphere transport dynamics for both water and nutrients, generating better tools to disentangle the plant influence (rhizosphere soil) from the soil physical properties (bulk soil) .
Although the explicit model couples the poroelastic mechanical model to the transport equations for water and nutrients, a more specific mechanical coupling might be more appropriate to define the tubers expanding in partially saturated soil. Partially saturated soils are not subject to consolidation (Yu 2000); thus, considering the soil as an elasto-viscoplastic media may be important in this situation (Ghezzehei and Or 2000). Furthermore, the mechanical stresses likely exceed typical yield stress values found in soil under field saturation conditions (Ghezzehei and Or 2001). Previous models have utilised strictly linear-elastic parameters to quantify the mechanics of cavity expansion in unsaturated soil (Aravena et al. 2014); however, future work should attempt to remedy this by considering soil plasticity.
This study was motivated by the growth of tubers in soil; however, the system of equations is not limited to this particular problem. Other biological processes could also be modelled, including, but not limited too, clusters of lymph nodes swelling under an inflammatory response from a disease or virus moving through a biological tissue (Yang et al. 2014), the growth of roots in response to water and nutrients (Aravena et al. 2014;Drew and Saker 1978), or to model the effect of tumour growth on nutrient flow during angiogenesis (Alarcón et al. 2003).
Technical analysis regarding the homogenisation procedure showed encouraging results. Comparing the results from the homogenised sets of equations to the full set yielded less than about a 2% difference between nutrient concentrations at different depths, as shown in Fig. 4. Despite similarities in the results, the homogenised set of equations could be solved 1000 times faster than the full set of equations. Furthermore, the homogenised model could be physically scaled up with minimal increases to computational time, while increasing the domain size for the explicit geometry will substantially increase the computational time. This is important if we were to do combinatorial simulations spanning large numbers of soil and climate parameters to predict how potato crops grow. Thus, the averaged model will computationally allow extensive explorations of soil management and crop breeding strategies to be investigated in silico. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.