The effect of non-uniform microscale distribution of sorption sites on solute diffusion in soil

,


Introduction
Soil is heterogeneous at all scales, but the degree of heterogeneity increases the finer the scale considered. For instance, at the scale of a root hair or fungal hypha, the distribution of mineral and organic matter surfaces on which nutrient and pollutant solutes are sorbed is patchy, and the patches have varying sizes and distributions. With advances in non-invasive imaging techniques, it is increasingly possible to quantify such heterogeneities and, in principle, to allow for them in models of soil processes (Leitner et al., 2010;Blunt et al., 2013;Keyes et al., 2013;Wildenschild & Sheppard, 2013). The research in this paper assesses the circumstances in which it is worth doing this.
with image-based modelling (Keyes et al., 2013;Daly et al., 2015;Tracy et al., 2015). Image-based modelling is the use of two-or three-dimensional images to derive true geometries on which to base mathematical models. Here we use three-dimensional images, generated with X-ray computed tomography (X-ray CT), to derive the true geometric impedance to solute transport in a particular experimental soil.
By combining image-based modelling with equations for time-dependent sorption reactions on the soil surfaces, we can allow for patchy distributions of the sorption sites. The resulting models enable us to test the assumption of microscale averaging used in conventional models.

Theory
The nomenclature is explained in Table 1. Table S1 in Supplementary Information provides further explanation about the mathematical symbols.

Dimensional model
We consider the transport of a solute through a soil volume across which there is a concentration gradient of the solute. The solute is distributed between the soil solution, where it is mobile, and the surfaces of soil particles, where it is sorbed and immobile. Unlike conventional models, the spatial distribution of sorption sites is considered to be patchy (i.e. there are areas of the soil surface where sorption does occur and areas where it does not). The sorption reactions we consider are either fast or slow in comparison with transport through the soil solution. For simplicity we consider a soil in which the pores are completely filled with solution and transport is solely by diffusion. Note, however, that it would be straightforward to model partially saturated soil by including air spaces in the pores.
To model the effects of a patchy distribution of sorption sites, we use the method of homogenization. The key requirement of this method is that the macroscale (observable) properties are, on average, independent of the microscale properties. Therefore, at the macroscale we do not need to know the exact paths taken by the solute. Rather, we need consider only the average solute flux through a microscale volume of soil. We use the microscale geometrical details to parameterize a set of macroscale equations that allow for the microscale detail through a representative 'cell problem' (Hornung, 1997;Pavliotis & Stuart, 2008). If we consider that the macro-and microscale spatial properties are independent of each other, we can define the macro-and microscale coordinate systems as x and y, respectively. This is illustrated in Figure 1. The characteristic macroscopic length-scale, L, is the length of the macroscopic unit in Figure 1(a) and the microscopic length-scale, l, is the length of the microscopic unit cell in Figure 1(b). Formally, homogenization theory allows us to treat x and y independently if the scaling ratio = l/L is sufficiently small (i.e. ≪ 1) (Pavliotis & Stuart, 2008).
We now consider the behaviour of the microscale domain in detail. Our aim is to derive a set of equations that can be applied Vector to represent the effect of random phase distribution for slow sorption reaction (dimensionless) -Γ Domain of total soil particle surface -Γ s Domain of sorption sites on soil surfaces -Γ ns Domain of non-sorbing sites on soil surfaces -Ω Domain of soil solution -∇ Partial differential operator in 3D space to the macroscale and a set of cell problems that capture the geometrical effects of the microscale domain. We define the soil solution domain as Ω and the total soil surface domain as Γ; Γ is subdivided into a domain of solute-sorbing sites, Γ s , and a domain of non-sorbing sites, Γ ns . Then, the diffusive transport of a solute in the soil solution is given by: where C l is the concentration of the solute in the soil solution and D l is its diffusion coefficient in free solution.
We define the amount of solute sorbed per unit of soil surface area in fast and slow reactions as C af and C as , respectively. The boundary conditions at the interface between soil surfaces and the solution for sorbing surfaces are: and for non-sorbing soil surfaces where n is a unit vector that is normal to the soil surface and points into the soil solution. Mass conservation is accommodated with linear first-order sorption reactions as follows: where k af , k df , k as and k ds are the spatially-dependent rate constants for fast and slow adsorption and desorption reactions, respectively.

Non-dimensionalization
To proceed with the analysis, we express Equations (1)-(5) in non-dimensional forms so that there are fewer parameters to deal with and it is easier to determine which processes are most important for any given parameter regime. We do this by scaling space with the macroscale length, x = Lx*, and time with the macroscale diffusion time, t = L 2 D l t * , where the asterisks indicate dimensionless variables. We scale the solution concentrations with unit concentration as C l = [C]C l *, where C is a representative solute concentration in solution. Likewise, the sorbed concentrations are scaled as C af = [C]C * af and C as = ( k as k ds ) [C]C * as to balance the adsorption and desorption terms in Equations (4) and (5).
Next, we non-dimensionalize the equations to reduce the number of parameters in the model, and then re-scale to estimate the importance of each of these parameters, as follows (without the asterisks): n · C l = 0, x ∈ Γ ns , ( 3 a ) We have scaled Equations (1a)-(5a) such that the parameters in the equations are all of size one. The only exception is the small parameter that we use in the forthcoming derivation. The parameters af , df , as and ds are all dimensionless because of the difference in dimensions in the adsorption and desorption constants; their values are given in Table S1 in the supporting information.

Homogenized models
We apply the method of homogenization to Equations (1a)-(5a) to obtain effective macroscale models. Homogenization is a multiscale method based on a standard perturbation expansion in which both the dependent and independent variables (in this case the concentrations and the spatial coordinates) are expanded. This is essentially an averaging procedure that allows the effect of the micro-structure to be determined as a series of effective parameters. A full description of the method and mathematical details are included in the Supporting Information S1. In short, we consider the behaviour of Equations (1a)-(5a) on two scales, the micro-then the macroscale, and consider gradients on the macroscale as a perturbation to the microscale equations. There are three key steps. First, we show that, as a first approximation, the nutrient concentrations are independent of the microscale geometry. Physically, this means that diffusion at the microscale is so fast that the concentration is effectively constant. Second, we consider that although the concentration may be considered constant at the microscale, it may vary at the macroscale. If the concentration changes linearly by one unit over unit distance, it must also change by a quantity of size over a distance . Therefore, we treat any variation in the macroscale concentration as a perturbation on the microscale. In addition, we have perturbations that arise from reactions at the soil particle surfaces. The result is that we must solve a set of representative problems at the microscale that we refer to as cell problems. Typically, these problems describe how the microscale geometry affects diffusion at the macroscale and they must be solved numerically for all but the simplest geometries. Finally, by averaging over the microscale we can eliminate all effects of the microstructure from the problem and are left with a series of macroscale equations. These are parameterized by a set of effective parameters, which are calculated by solving the cell problems on the underlying microstructure.
We present two examples; each requires a separate homogenized model. In Case 1, the fast sorption reaction is effectively instantaneous and the slow reaction is slow in comparison with diffusion through the soil solution (discussed further in Parameter values). In Case 2, the slow reaction is 100-fold slower (i.e. k as is 100-fold smaller than in Case 1). We summarize the models for the two examples in this section and provide full details in the Supporting Information, SI.
To derive the macroscale equations, we use the following asymptotic expansion of the concentration variables with respect to as follows: where C 0 l (x), C 0 af (x) and C 0 as (x) are the leading order macroscale solute concentrations in solution and on the fast and slow sorption sites, respectively, and O( 2 ) indicates remainder terms. These leading order concentrations represent the effect of phase distribution for uniform distribution of sorbing phases. The concentrations for the non-uniform sorbing phase distributions on the soil surface are found as a perturbation to the average concentration and are given by C 1 l (x) , C 1 af (x) and C 1 as (x). If Equation (6) is substituted into Equations (1a)-(5a) we obtain, after the mathematical derivation given in Supporting Information B in Appendix S1, a homogenized macroscale model for Case 1: where D is a second-rank tensor that accounts for the geometric tortuosity of the solution diffusion pathway. Its value is obtained by solving the cell equations for the microscale pore space as shown in the supporting information. The reaction time constant ( ) in Equation (7) is: where ||Γ s || and ||Ω|| are the magnitudes of the non-dimensional sorbing phase area and soil solution volume, respectively. By considering terms of higher powers in , the macroscale homogenized model that allows for non-uniform spatial distribution of sorption phases is: where E is a third-rank tensor that accounts for the distribution of fast sorption sites and f is a vector that accounts for the distribution of slow sorption sites in the microscale domain (detail is provided in Supporting Information A and B). In addition: We now present the effective macroscale models for Case 2 where the slow sorption reaction is 100 times slower than in Case 1. The effective macroscale model with uniform sorption-site distribution is the conventional model for macro-scale diffusive transport: and the homogenized model at order O( ) is The diffusive behaviour of the solutes at the macroscale is the same as for Case 1. In this example however, the slow reaction is present, to leading order, only through the effective time constant , Equation (11). The effect of the slow reaction is described by Equation (12) and can be seen to contribute to C 1 l . In Case 1 the effect of the non-uniform sorption sites was described by the vector f, see Equation (9). In this example these terms do not occur until the next order. In other words, as the reaction is slowed down it produces a weaker effect on the overall solute concentration.
In summary, the homogenized macroscale equations suggest that the maximal effect of the spatial distribution of sorption sites occurs with the parameter regime in Case 1. At this point it is natural to ask whether there is a reaction rate for which the non-uniformity of sorption sites has a larger effect on solute transport. If we slow the reactions down further, then this effect is simply pushed to the higher orders (i.e. it becomes smaller). Alternatively, if we increase the speed of the slow reaction then it becomes comparable to, or faster than, diffusion. If this is the case then we see only the averaged effects. In all cases, the effect of non-uniform distribution of sorption sites appears to be at least an order later than that of averaged or uniform distributions, which indicates that the non-uniform distribution has a small effect only on the macroscale diffusion. Physically, this indicates that as the size of the macroscale domain (i.e. L) increases, the effects of the non-uniform distribution decrease. We illustrate this point further below.
We rescale the dimensionless terms in the macroscale homogenized equations with the corresponding scales given in the non-dimensionalization section above. The dimensional forms of Equations (7) and (9) are: ( and ( where D eff = D l D, E eff = LD l E and f eff = fD l /L are the effective homogenized model parameters with units cm 2 s −1 , cm 3 s − 1 and cm s − 1 , respectively, is the soil volumetric water content, and K f and K s are the equilibrium distributions of solute between the surface sites and the soil solution, given by: and The dimensional form of Equation (10) is: Similarly, we can obtain the dimensional forms of the homogenized macroscale Equations (11) and (12) for Case 2; these have been omitted for conciseness.

Parameter values
We parameterize the model for the sandy loam soil shown in Figure 1 and a strongly-sorbed solute such as phosphate or many heavy metal and radionuclide contaminants. The soil is a Eutric Cambisol (USDA, 2014) from the experimental farm of Bangor University (properties in Lucas & Jones, 2006, Soil B). It was sieved to <1.6 mm and the solid matter of the sieved soil was approximately 0.7 g g −1 quartz, 0.3 g g −1 clay minerals (mainly illite and feldspars) and <0.03 g g −1 organic matter.
Geometry. We consider an experimental system where the typical length of a soil column L = 1 cm. We define this as the macroscale because it is the largest scale length considered in this research. In addition, we have a scale length that is associated with the soil particle size distribution. We define this as the microscale l (Figure 1). We chose l = 0.01 cm to cover a range of particle sizes from fine clay to coarse sand in the microscale domain. This satisfies the requirement that the scaling ratio is = l/L ≪ 1, so that our two space-scale approximations are valid. We consider the effect of varying in 'Results and discussion' below. We obtained the value of the tensor D for the geometric impedance to diffusion in the soil pore space from the X-ray CT image in Figure 1; the value is D = 0.035 × I, where I is the identity matrix (Supporting Information B). For D l = 10 −5 cm 2 s −1 , which is typical of a simple inorganic solute, the effective diffusion coefficient in the soil pore network following Equation (13) is D eff = 3.5 × 10 −7 cm 2 s −1 . The mathematical procedure to obtain the effective macroscale parameters (i.e. tensor E and vector f (for slow reaction)) from the microscale domain is given in Supporting Information B. From the simulation we find the components of E eff take a maximum value of 8.25 × 10 −9 cm 3 s −1 and each of the components of f eff has value 1.1 × 10 −7 cm s −1 .
We calculated the surface area of clay minerals active in solute sorption as follows. From the X-ray CT image (Figure 1), the volume of soil solid per unit soil volume is 0.51 cm 3 cm −3 and the surface area of the soil solid per unit soil volume is 85.25 cm 2 cm −3 . Therefore, if the particle density is 2.65 g cm −3 , the mass of quartz, which is largely non-sorbing, per unit soil volume = 0.7 × 2.65 × 0.51 = 0.95 g cm −3 , and the mass of clay minerals per unit soil volume = 0.3 × 2.65 × 0.51 = 0.41 g cm −3 . Therefore, taking the quartz particles to be spherical with mean diameter 1.6 mm (the soil sieve size), the quartz specific surface is 14.15 cm 2 g −1 ; therefore, the quartz surface area per unit soil volume = 0.95 × 14.15 = 13.43 cm 2 cm −3 . The clay mineral surface area equals the total surface area minus the quartz surface area. This gives clay mineral surface area per unit volume = 85.25 − 13.43 = 71.82 cm 2 cm −3 (i.e. the clay minerals account for 84% of the total surface).
After estimating the surface area of the sorbing phase, we distributed sorption sites of equal total area over the microscale soil domain microscale soil domain. This was done by selecting portions of the soil CT-image surfaces as patches that collectively occupied 84% of the total surface (see above). We applied the sorption reactions on these spatially distributed soil surface patches following Equations (4) and (5). The remaining 16% of the total soil surface is non-sorbing; therefore, we solved Equation (3) on those surface areas. We give the mathematical detail of the distribution patterns (both uniform and non-uniform) of the sorption sites in Supporting Information A and B.
Sorption. We define the fraction of the total sorption surface area occupied by slow-reacting sites as , such that A s = A total and A f = (1 − )A total where A s and A f are the areas of slow and fast reacting sites, respectively, per unit volume of soil, and A total = A s + A f . To test the effects of varying (i) the ratio of A s to A f and (ii) the rate of the slow sorption reaction, we derived the relations at equilibrium from Equations (4) to (5) by: and k as C l − k ds C as = 0.
Combining Equations (16)- (19) gives: and We tested the effect of varying k as and (i.e. A s /A f ) while keeping the other variables constant (see 'Results and discussion'). For the fast reaction, we used k af = 1.78 × 10 −5 cm s −1 and k df = 2.57 × 10 −4 s −1 based on the data of Keyes et al. (2013) for the short-term (<2 hours) kinetics of phosphate sorption on the soil in Figure 1 and A f = 71.82 cm 2 cm −3 (i.e. K f = 5). For the slow reaction in Case 1 we used k as = 6.67 × 10 −7 cm s −1 and k ds = 7.05 × 10 −8 s −1 based on values in Ptashnyk et al. (2010) and A s = 71.82 cm 2 cm −3 (i.e. K s = 680). For the slow reaction in Case 2, k as = 6.67 × 10 −9 cm s −1 and k ds = 7.05 × 10 −10 s −1 .
We compare the time-scales of the fast and slow reactions with that of diffusion as follows. The half-time for diffusion into a soil column of length L maintained at a constant surface concentration is given by Crank (1979), equation 3.15). Therefore, for D = D 1 /K f = 2 × 10 −6 cm 2 s −1 and L = 1 cm, the diffusion half-time is t 1/2 = 27.3 hours. From the integral of Equation (4) for constant C l , the half-time for the fast sorption reaction is t 1/2 = ln 2/k df = 2.7 × 10 3 s or 0.75 hour, which is rapid compared with diffusion. Likewise, from Equation (5) the half-time for the slow sorption reaction is t 1/2 = ln 2/k ds = 9.8 × 10 6 s or 114 days for Case 1, which is slow compared with diffusion. Similarly for Case 2, the half-time for the slow reaction is 1.14 × 10 5 days. Therefore, the representations of sorption times in the two Cases are consistent with the model assumptions.
We ran the model for a constant solute concentration in solution (representative solute concentration, C) at one end of a soil column of 10 −4 μmol cm −3 , with zero transfer of solute across the other end of the column (L = 1 cm) and zero concentration initially throughout the column. The results below are for model runs of 25 days.

Model implementation
We used the open-source computational fluid dynamics software package OpenFOAM (Jasak et al., 2013) for the numerical simulations with the two homogenized models, and to assess their sensitivities to the spatial distributions of sorption sites. Figure 2(a,b) shows the effects of allowing for patchy microscale distributions of sorption sites on solute concentration profiles in the soil solution and soil surface when all the sorption sites react more rapidly than rates of diffusion. It appears that microscale patchiness has little effect on macroscale transport under these conditions. If sorption is rapid compared with diffusion (i.e. half-times for equilibration <∼1 hour), it is acceptable to consider the volume-averaged distribution of sorption sites without knowing their distribution patterns. Figure 2(c,d) shows the effect of both fast and slow sorption reactions, with the slow sorption being either moderately slower (Case 1) or much slower (Case 2) than diffusion. The effect of a patchy distribution of sites is evident for Case 1, although it is only slight, but is not for Case 2. Evidently there is an interaction between the reaction rate and the patchiness effect. In Case 2, the slow reaction is so slow that it does not contribute markedly to sorption on the time-scale of the model runs; the results are therefore similar to those in Figure 2(a,b), but with a smaller degree of sorption (smaller A f ) and faster macroscale diffusion. If the slow reaction is fast enough to affect diffusion, as in Case 1, a small effect of patchiness is evident in Figure 2(c,d). Figure 2(e,f) shows the results of when all the sorption sites react slowly. For Case 1, there is again a small effect of patchiness only. For Case 2, Figure 2(e,f) shows that there is now effectively no sorption and macroscale diffusion is correspondingly faster. Figure 3 shows the effect of a larger sorption surface area with a relatively fast slow reaction (i.e. Case 1). The concentration profiles are similar to those in Figure 2(e) (Case 1), which show no marked effect.

Results and discussion
The results show that the effect of non-uniform distribution of sorption sites on macroscale transport is detectable only for particular combinations of reaction rate and diffusion rate, as in Case 1, and then the effect is only slight and probably too small to be detectable experimentally. Therefore, the averaging assumptions used in conventional models with equilibrium reactions are justified.
We now consider the effect of increasing the micro-to macroscale ratio, . For Case 1 with a slow reaction rate, a tenfold increase in results in a large increase in the effect of patchily distributed sorption sites (Figure 4). The distance at which C l falls to half the value at the surface of the soil column is 0.014 cm for a uniform distribution of sites and 0.045 cm for a patchy distribution. This indicates the parameter range for which the effect of patchiness is important and correctly predicted. The effect depends on both and the sorption rate relative to the diffusion rate: must be small enough to satisfy the multiple space-scale perturbation requirements, but large enough to resolve the effect of patchiness with the right sorption rate. The reaction half-time is 114 days (Parameter values). The model suggests that for = 0.1, reaction rates 100-fold faster or slower than this will mean the reaction is either effectively instantaneous or insignificant compared with diffusion, and the effect of non-uniform distribution of sorption sites then becomes negligible. Runs with the model confirmed this (data not shown).
From a mathematical point of view, when is infinitely small the model predictions are exact. With an increase in , the error of prediction will also increase. The error is of the order 2 × 100%; it is just 1% for = 0.1. From a physical point of view, if we consider a soil column packed with soil particles that are small compared with the size of the column, the effect of microscale heterogeneities such as the distribution of sorption sites is insignificant. Likewise, at the pedon or field scale processes that take place on the micro-aggregate scale are irrelevant. Clearly, however, when the soil particle or aggregate size is an order of magnitude different only from the macroscale dimension that is being considered, the effects of particle or aggregate scale heterogeneities will be magnified and Figure 2 The effect of maintaining a constant concentration of a solute in solution at the surface of a soil column from which it was initially absent: (a, c and e) are solute concentrations in the soil solution (C l ), (b, d and f) are adsorbed solute concentrations in the soil solid (C af , C as for fast and slow reacting sites, respectively). (a), (b) All sorption sites were fast reacting (i.e. = 0). (c), (d) Half the sorption sites were fast reacting and half were slow (i.e. = 0.5), with slow reaction slow (Case 1) or very slow (Case 2) compared with diffusion. (e), (f) All the sorption sites were slow reacting (i.e. = 1) with slow reaction slow (Case 1) or very slow (Case 2) compared with diffusion. Solid lines are for a uniform distribution of sorption sites, points are for a non-uniform distribution and t = 25 days.
can be important. As the scaling ratio, , increases, the microscopic grain becomes closer to macroscopic soil and therefore imposes an increasing effect on model predictions.
The above conclusions are consistent with equivalent research on solute transport in groundwater systems Bosma et al., 1993;Espinoza & Valocchi, 1997;Dentz et al., 2011). Meile & Tuncay (2006) showed, based on pore-scale simulations, that in diffusion-driven transport the effect of non-homogeneity in reaction rates has only a small effect on the macroscale reactions.
The above conclusions will not necessarily hold if the solute reaction rates or other sources or sinks of the solutes themselves depend on microscale heterogeneities. The microscale distribution of Figure 3 Concentration-distance profiles in the soil solution as in Figure 2 for a tenfold larger sorption surface area (i.e. 1.1 × A total ). All sorption sites were slow reacting (i.e. = 1), with reaction rates fast enough to affect diffusion (i.e. Case 1). Solid lines are for a uniform distribution of sorption sites, points are for a non-uniform distribution and t = 25 days.

Figure 4
Concentration-distance profiles in the soil solution as in Figure  2 for micro-to macroscale ratio ( ) = 0.1 (i.e. tenfold smaller than in Figure 2). All sorption sites were slow reacting (i.e. = 1), with reaction rates fast enough to affect diffusion (i.e. Case 1). Solid lines are for a uniform distribution of sorption sites, points are for a non-uniform distribution and t = 25 days. sorption sites might then be important. An example is the uptake of solutes by root hairs, which have lengths comparable to the soil geometric microscale (Keyes et al., 2013). Soil processes at such scales are increasingly measurable with non-invasive imaging techniques. Another example is the microbially-mediated turnover of soil organic matter, which Falconer et al. (2015) show, with a model that allows for non-linear microbial growth kinetics, depends strongly on the microscale distribution of both nutrients and microbes and associated local transport limitations. This implies that macroscale models of soil carbon turnover should account for microscale heterogeneities explicitly. This could be achieved with our modelling approach by incorporating the microbial growth kinetics in the microscale equations (the equivalents to Equations 2 and 3).
For simplicity, we considered a water-saturated soil in our simulations. Our approach could be applied easily to non-saturated soil, however. The geometry we have used is taken from X-ray CT images in which air and water can be clearly distinguished. By treating the air-water interface within a pore as a non-sorbing site, we could model partially saturated soil in the same way as saturated soil. The same macroscopic equations will apply, but the model parameters derived from the pore geometry will be different.
We have not considered advection. In general, for strongly-sorbed solutes rates of advection are small and diffusion is the dominant transport process (Roose & Kirk, 2009). Our approach could also be used, however, for cases where advection is important, such as for weakly sorbed solutes.

Conclusions
We have found that, in general, the distribution of sorption sites at the microscale has little effect on diffusion at the macroscale. This indicates that, in most cases, the assumption implicit in most conventional solute transport models that microscale heterogeneities are unimportant at the macroscale, is justified. Microscale heterogeneities, however, do become important if (i) the ratio of the microto macro-length scales, , is large and (ii) solute sorption on soil surfaces is slow relative to diffusion, but not so slow that it does not affect diffusion. Nevertheless, if we consider a small sample of soil, such as is often used in synchrotron computed-tomography experiments, then the effect of microscale heterogeneities at sorption sites will be significant. The modelling error will also be amplified, varying as 100 × 2 %. It will, however, remain smaller than the effects of the distribution of sorption sites.

Supporting Information
The following supporting information is available in the online version of this article: Supplementary Information, SI: Further modelling details.